1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
|
<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=US-ASCII">
<title>Find Scale (Standard Deviation) Example</title>
<link rel="stylesheet" href="../../../../math.css" type="text/css">
<meta name="generator" content="DocBook XSL Stylesheets V1.78.1">
<link rel="home" href="../../../../index.html" title="Math Toolkit 2.1.0">
<link rel="up" href="../find_eg.html" title="Find Location and Scale Examples">
<link rel="prev" href="find_location_eg.html" title="Find Location (Mean) Example">
<link rel="next" href="find_mean_and_sd_eg.html" title="Find mean and standard deviation example">
</head>
<body bgcolor="white" text="black" link="#0000FF" vlink="#840084" alink="#0000FF">
<table cellpadding="2" width="100%"><tr>
<td valign="top"><img alt="Boost C++ Libraries" width="277" height="86" src="../../../../../../../../boost.png"></td>
<td align="center"><a href="../../../../../../../../index.html">Home</a></td>
<td align="center"><a href="../../../../../../../../libs/libraries.htm">Libraries</a></td>
<td align="center"><a href="http://www.boost.org/users/people.html">People</a></td>
<td align="center"><a href="http://www.boost.org/users/faq.html">FAQ</a></td>
<td align="center"><a href="../../../../../../../../more/index.htm">More</a></td>
</tr></table>
<hr>
<div class="spirit-nav">
<a accesskey="p" href="find_location_eg.html"><img src="../../../../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../find_eg.html"><img src="../../../../../../../../doc/src/images/up.png" alt="Up"></a><a accesskey="h" href="../../../../index.html"><img src="../../../../../../../../doc/src/images/home.png" alt="Home"></a><a accesskey="n" href="find_mean_and_sd_eg.html"><img src="../../../../../../../../doc/src/images/next.png" alt="Next"></a>
</div>
<div class="section">
<div class="titlepage"><div><div><h5 class="title">
<a name="math_toolkit.stat_tut.weg.find_eg.find_scale_eg"></a><a class="link" href="find_scale_eg.html" title="Find Scale (Standard Deviation) Example">Find
Scale (Standard Deviation) Example</a>
</h5></div></div></div>
<p>
First we need some includes to access the <a class="link" href="../../../dist_ref/dists/normal_dist.html" title="Normal (Gaussian) Distribution">Normal
Distribution</a>, the algorithms to find scale (and some std output
of course).
</p>
<pre class="programlisting"><span class="preprocessor">#include</span> <span class="special"><</span><span class="identifier">boost</span><span class="special">/</span><span class="identifier">math</span><span class="special">/</span><span class="identifier">distributions</span><span class="special">/</span><span class="identifier">normal</span><span class="special">.</span><span class="identifier">hpp</span><span class="special">></span> <span class="comment">// for normal_distribution</span>
<span class="keyword">using</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">normal</span><span class="special">;</span> <span class="comment">// typedef provides default type is double.</span>
<span class="preprocessor">#include</span> <span class="special"><</span><span class="identifier">boost</span><span class="special">/</span><span class="identifier">math</span><span class="special">/</span><span class="identifier">distributions</span><span class="special">/</span><span class="identifier">find_scale</span><span class="special">.</span><span class="identifier">hpp</span><span class="special">></span>
<span class="keyword">using</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">find_scale</span><span class="special">;</span>
<span class="keyword">using</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">complement</span><span class="special">;</span> <span class="comment">// Needed if you want to use the complement version.</span>
<span class="keyword">using</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">policies</span><span class="special">::</span><span class="identifier">policy</span><span class="special">;</span> <span class="comment">// Needed to specify the error handling policy.</span>
<span class="preprocessor">#include</span> <span class="special"><</span><span class="identifier">iostream</span><span class="special">></span>
<span class="keyword">using</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span><span class="special">;</span> <span class="keyword">using</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span>
<span class="preprocessor">#include</span> <span class="special"><</span><span class="identifier">iomanip</span><span class="special">></span>
<span class="keyword">using</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">setw</span><span class="special">;</span> <span class="keyword">using</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">setprecision</span><span class="special">;</span>
<span class="preprocessor">#include</span> <span class="special"><</span><span class="identifier">limits</span><span class="special">></span>
<span class="keyword">using</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special">;</span>
</pre>
<p>
For this example, we will use the standard <a class="link" href="../../../dist_ref/dists/normal_dist.html" title="Normal (Gaussian) Distribution">Normal
Distribution</a>, with location (mean) zero and standard deviation
(scale) unity. Conveniently, this is also the default for this implementation's
constructor.
</p>
<pre class="programlisting"><span class="identifier">normal</span> <span class="identifier">N01</span><span class="special">;</span> <span class="comment">// Default 'standard' normal distribution with zero mean</span>
<span class="keyword">double</span> <span class="identifier">sd</span> <span class="special">=</span> <span class="number">1.</span><span class="special">;</span> <span class="comment">// and standard deviation is 1.</span>
</pre>
<p>
Suppose we want to find a different normal distribution with standard
deviation so that only fraction p (here 0.001 or 0.1%) are below a certain
chosen limit (here -2. standard deviations).
</p>
<pre class="programlisting"><span class="keyword">double</span> <span class="identifier">z</span> <span class="special">=</span> <span class="special">-</span><span class="number">2.</span><span class="special">;</span> <span class="comment">// z to give prob p</span>
<span class="keyword">double</span> <span class="identifier">p</span> <span class="special">=</span> <span class="number">0.001</span><span class="special">;</span> <span class="comment">// only 0.1% below z = -2</span>
<span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"Normal distribution with mean = "</span> <span class="special"><<</span> <span class="identifier">N01</span><span class="special">.</span><span class="identifier">location</span><span class="special">()</span> <span class="comment">// aka N01.mean()</span>
<span class="special"><<</span> <span class="string">", standard deviation "</span> <span class="special"><<</span> <span class="identifier">N01</span><span class="special">.</span><span class="identifier">scale</span><span class="special">()</span> <span class="comment">// aka N01.standard_deviation()</span>
<span class="special"><<</span> <span class="string">", has "</span> <span class="special"><<</span> <span class="string">"fraction <= "</span> <span class="special"><<</span> <span class="identifier">z</span>
<span class="special"><<</span> <span class="string">", p = "</span> <span class="special"><<</span> <span class="identifier">cdf</span><span class="special">(</span><span class="identifier">N01</span><span class="special">,</span> <span class="identifier">z</span><span class="special">)</span> <span class="special"><<</span> <span class="identifier">endl</span><span class="special">;</span>
<span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"Normal distribution with mean = "</span> <span class="special"><<</span> <span class="identifier">N01</span><span class="special">.</span><span class="identifier">location</span><span class="special">()</span>
<span class="special"><<</span> <span class="string">", standard deviation "</span> <span class="special"><<</span> <span class="identifier">N01</span><span class="special">.</span><span class="identifier">scale</span><span class="special">()</span>
<span class="special"><<</span> <span class="string">", has "</span> <span class="special"><<</span> <span class="string">"fraction > "</span> <span class="special"><<</span> <span class="identifier">z</span>
<span class="special"><<</span> <span class="string">", p = "</span> <span class="special"><<</span> <span class="identifier">cdf</span><span class="special">(</span><span class="identifier">complement</span><span class="special">(</span><span class="identifier">N01</span><span class="special">,</span> <span class="identifier">z</span><span class="special">))</span> <span class="special"><<</span> <span class="identifier">endl</span><span class="special">;</span> <span class="comment">// Note: uses complement.</span>
</pre>
<pre class="programlisting">Normal distribution with mean = 0 has fraction <= -2, p = 0.0227501
Normal distribution with mean = 0 has fraction > -2, p = 0.97725
</pre>
<p>
Noting that p = 0.02 instead of our target of 0.001, we can now use
<code class="computeroutput"><span class="identifier">find_scale</span></code> to give a
new standard deviation.
</p>
<pre class="programlisting"><span class="keyword">double</span> <span class="identifier">l</span> <span class="special">=</span> <span class="identifier">N01</span><span class="special">.</span><span class="identifier">location</span><span class="special">();</span>
<span class="keyword">double</span> <span class="identifier">s</span> <span class="special">=</span> <span class="identifier">find_scale</span><span class="special"><</span><span class="identifier">normal</span><span class="special">>(</span><span class="identifier">z</span><span class="special">,</span> <span class="identifier">p</span><span class="special">,</span> <span class="identifier">l</span><span class="special">);</span>
<span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"scale (standard deviation) = "</span> <span class="special"><<</span> <span class="identifier">s</span> <span class="special"><<</span> <span class="identifier">endl</span><span class="special">;</span>
</pre>
<p>
that outputs:
</p>
<pre class="programlisting">scale (standard deviation) = 0.647201
</pre>
<p>
showing that we need to reduce the standard deviation from 1. to 0.65.
</p>
<p>
Then we can check that we have achieved our objective by constructing
a new distribution with the new standard deviation (but same zero mean):
</p>
<pre class="programlisting"><span class="identifier">normal</span> <span class="identifier">np001pc</span><span class="special">(</span><span class="identifier">N01</span><span class="special">.</span><span class="identifier">location</span><span class="special">(),</span> <span class="identifier">s</span><span class="special">);</span>
</pre>
<p>
And re-calculating the fraction below (and above) our chosen limit.
</p>
<pre class="programlisting"><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"Normal distribution with mean = "</span> <span class="special"><<</span> <span class="identifier">l</span>
<span class="special"><<</span> <span class="string">" has "</span> <span class="special"><<</span> <span class="string">"fraction <= "</span> <span class="special"><<</span> <span class="identifier">z</span>
<span class="special"><<</span> <span class="string">", p = "</span> <span class="special"><<</span> <span class="identifier">cdf</span><span class="special">(</span><span class="identifier">np001pc</span><span class="special">,</span> <span class="identifier">z</span><span class="special">)</span> <span class="special"><<</span> <span class="identifier">endl</span><span class="special">;</span>
<span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"Normal distribution with mean = "</span> <span class="special"><<</span> <span class="identifier">l</span>
<span class="special"><<</span> <span class="string">" has "</span> <span class="special"><<</span> <span class="string">"fraction > "</span> <span class="special"><<</span> <span class="identifier">z</span>
<span class="special"><<</span> <span class="string">", p = "</span> <span class="special"><<</span> <span class="identifier">cdf</span><span class="special">(</span><span class="identifier">complement</span><span class="special">(</span><span class="identifier">np001pc</span><span class="special">,</span> <span class="identifier">z</span><span class="special">))</span> <span class="special"><<</span> <span class="identifier">endl</span><span class="special">;</span>
</pre>
<pre class="programlisting">Normal distribution with mean = 0 has fraction <= -2, p = 0.001
Normal distribution with mean = 0 has fraction > -2, p = 0.999
</pre>
<h5>
<a name="math_toolkit.stat_tut.weg.find_eg.find_scale_eg.h0"></a>
<span class="phrase"><a name="math_toolkit.stat_tut.weg.find_eg.find_scale_eg.controlling_how_errors_from_find"></a></span><a class="link" href="find_scale_eg.html#math_toolkit.stat_tut.weg.find_eg.find_scale_eg.controlling_how_errors_from_find">Controlling
how Errors from find_scale are handled</a>
</h5>
<p>
We can also control the policy for handling various errors. For example,
we can define a new (possibly unwise) policy to ignore domain errors
('bad' arguments).
</p>
<p>
Unless we are using the boost::math namespace, we will need:
</p>
<pre class="programlisting"><span class="keyword">using</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">policies</span><span class="special">::</span><span class="identifier">policy</span><span class="special">;</span>
<span class="keyword">using</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">policies</span><span class="special">::</span><span class="identifier">domain_error</span><span class="special">;</span>
<span class="keyword">using</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">policies</span><span class="special">::</span><span class="identifier">ignore_error</span><span class="special">;</span>
</pre>
<p>
Using a typedef is convenient, especially if it is re-used, although
it is not required, as the various examples below show.
</p>
<pre class="programlisting"><span class="keyword">typedef</span> <span class="identifier">policy</span><span class="special"><</span><span class="identifier">domain_error</span><span class="special"><</span><span class="identifier">ignore_error</span><span class="special">></span> <span class="special">></span> <span class="identifier">ignore_domain_policy</span><span class="special">;</span>
<span class="comment">// find_scale with new policy, using typedef.</span>
<span class="identifier">l</span> <span class="special">=</span> <span class="identifier">find_scale</span><span class="special"><</span><span class="identifier">normal</span><span class="special">>(</span><span class="identifier">z</span><span class="special">,</span> <span class="identifier">p</span><span class="special">,</span> <span class="identifier">l</span><span class="special">,</span> <span class="identifier">ignore_domain_policy</span><span class="special">());</span>
<span class="comment">// Default policy policy<>, needs using boost::math::policies::policy;</span>
<span class="identifier">l</span> <span class="special">=</span> <span class="identifier">find_scale</span><span class="special"><</span><span class="identifier">normal</span><span class="special">>(</span><span class="identifier">z</span><span class="special">,</span> <span class="identifier">p</span><span class="special">,</span> <span class="identifier">l</span><span class="special">,</span> <span class="identifier">policy</span><span class="special"><>());</span>
<span class="comment">// Default policy, fully specified.</span>
<span class="identifier">l</span> <span class="special">=</span> <span class="identifier">find_scale</span><span class="special"><</span><span class="identifier">normal</span><span class="special">>(</span><span class="identifier">z</span><span class="special">,</span> <span class="identifier">p</span><span class="special">,</span> <span class="identifier">l</span><span class="special">,</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">policies</span><span class="special">::</span><span class="identifier">policy</span><span class="special"><>());</span>
<span class="comment">// New policy, without typedef.</span>
<span class="identifier">l</span> <span class="special">=</span> <span class="identifier">find_scale</span><span class="special"><</span><span class="identifier">normal</span><span class="special">>(</span><span class="identifier">z</span><span class="special">,</span> <span class="identifier">p</span><span class="special">,</span> <span class="identifier">l</span><span class="special">,</span> <span class="identifier">policy</span><span class="special"><</span><span class="identifier">domain_error</span><span class="special"><</span><span class="identifier">ignore_error</span><span class="special">></span> <span class="special">>());</span>
</pre>
<p>
If we want to express a probability, say 0.999, that is a complement,
<code class="computeroutput"><span class="number">1</span> <span class="special">-</span>
<span class="identifier">p</span></code> we should not even think
of writing <code class="computeroutput"><span class="identifier">find_scale</span><span class="special"><</span><span class="identifier">normal</span><span class="special">>(</span><span class="identifier">z</span><span class="special">,</span> <span class="number">1</span> <span class="special">-</span>
<span class="identifier">p</span><span class="special">,</span>
<span class="identifier">l</span><span class="special">)</span></code>,
but use the <a class="link" href="../../overview/complements.html" title="Complements are supported too - and when to use them">complements</a>
version (see <a class="link" href="../../overview/complements.html#why_complements">why complements?</a>).
</p>
<pre class="programlisting"><span class="identifier">z</span> <span class="special">=</span> <span class="special">-</span><span class="number">2.</span><span class="special">;</span>
<span class="keyword">double</span> <span class="identifier">q</span> <span class="special">=</span> <span class="number">0.999</span><span class="special">;</span> <span class="comment">// = 1 - p; // complement of 0.001.</span>
<span class="identifier">sd</span> <span class="special">=</span> <span class="identifier">find_scale</span><span class="special"><</span><span class="identifier">normal</span><span class="special">>(</span><span class="identifier">complement</span><span class="special">(</span><span class="identifier">z</span><span class="special">,</span> <span class="identifier">q</span><span class="special">,</span> <span class="identifier">l</span><span class="special">));</span>
<span class="identifier">normal</span> <span class="identifier">np95pc</span><span class="special">(</span><span class="identifier">l</span><span class="special">,</span> <span class="identifier">sd</span><span class="special">);</span> <span class="comment">// Same standard_deviation (scale) but with mean(scale) shifted</span>
<span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"Normal distribution with mean = "</span> <span class="special"><<</span> <span class="identifier">l</span> <span class="special"><<</span> <span class="string">" has "</span>
<span class="special"><<</span> <span class="string">"fraction <= "</span> <span class="special"><<</span> <span class="identifier">z</span> <span class="special"><<</span> <span class="string">" = "</span> <span class="special"><<</span> <span class="identifier">cdf</span><span class="special">(</span><span class="identifier">np95pc</span><span class="special">,</span> <span class="identifier">z</span><span class="special">)</span> <span class="special"><<</span> <span class="identifier">endl</span><span class="special">;</span>
<span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"Normal distribution with mean = "</span> <span class="special"><<</span> <span class="identifier">l</span> <span class="special"><<</span> <span class="string">" has "</span>
<span class="special"><<</span> <span class="string">"fraction > "</span> <span class="special"><<</span> <span class="identifier">z</span> <span class="special"><<</span> <span class="string">" = "</span> <span class="special"><<</span> <span class="identifier">cdf</span><span class="special">(</span><span class="identifier">complement</span><span class="special">(</span><span class="identifier">np95pc</span><span class="special">,</span> <span class="identifier">z</span><span class="special">))</span> <span class="special"><<</span> <span class="identifier">endl</span><span class="special">;</span>
</pre>
<p>
Sadly, it is all too easy to get probabilities the wrong way round, when
you may get a warning like this:
</p>
<pre class="programlisting">Message from thrown exception was:
Error in function boost::math::find_scale<Dist, Policy>(complement(double, double, double, Policy)):
Computed scale (-0.48043523852179076) is <= 0! Was the complement intended?
</pre>
<p>
The default error handling policy is to throw an exception with this
message, but if you chose a policy to ignore the error, the (impossible)
negative scale is quietly returned.
</p>
<p>
See <a href="../../../../../../example/find_scale_example.cpp" target="_top">find_scale_example.cpp</a>
for full source code: the program output looks like this:
</p>
<pre class="programlisting">Example: Find scale (standard deviation).
Normal distribution with mean = 0, standard deviation 1, has fraction <= -2, p = 0.0227501
Normal distribution with mean = 0, standard deviation 1, has fraction > -2, p = 0.97725
scale (standard deviation) = 0.647201
Normal distribution with mean = 0 has fraction <= -2, p = 0.001
Normal distribution with mean = 0 has fraction > -2, p = 0.999
Normal distribution with mean = 0.946339 has fraction <= -2 = 0.001
Normal distribution with mean = 0.946339 has fraction > -2 = 0.999
</pre>
</div>
<table xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision" width="100%"><tr>
<td align="left"></td>
<td align="right"><div class="copyright-footer">Copyright © 2006-2010, 2012-2014 Nikhar Agrawal,
Anton Bikineev, Paul A. Bristow, Marco Guazzone, Christopher Kormanyos, Hubert
Holin, Bruno Lalande, John Maddock, Johan Råde, Gautam Sewani, Benjamin Sobotta,
Thijs van den Berg, Daryle Walker and Xiaogang Zhang<p>
Distributed under the Boost Software License, Version 1.0. (See accompanying
file LICENSE_1_0.txt or copy at <a href="http://www.boost.org/LICENSE_1_0.txt" target="_top">http://www.boost.org/LICENSE_1_0.txt</a>)
</p>
</div></td>
</tr></table>
<hr>
<div class="spirit-nav">
<a accesskey="p" href="find_location_eg.html"><img src="../../../../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../find_eg.html"><img src="../../../../../../../../doc/src/images/up.png" alt="Up"></a><a accesskey="h" href="../../../../index.html"><img src="../../../../../../../../doc/src/images/home.png" alt="Home"></a><a accesskey="n" href="find_mean_and_sd_eg.html"><img src="../../../../../../../../doc/src/images/next.png" alt="Next"></a>
</div>
</body>
</html>
|
