1 // find_scale.cpp
2
3 // Copyright Paul A. Bristow 2007, 2010.
4
5 // Use, modification and distribution are subject to the
6 // Boost Software License, Version 1.0.
7 // (See accompanying file LICENSE_1_0.txt
8 // or copy at http://www.boost.org/LICENSE_1_0.txt)
9
10 // Example of finding scale (standard deviation) for normal (Gaussian).
11
12 // Note that this file contains Quickbook mark-up as well as code
13 // and comments, don't change any of the special comment mark-ups!
14
15 //[find_scale1
16 /*`
17 First we need some includes to access the __normal_distrib,
18 the algorithms to find scale (and some std output of course).
19 */
20
21 #include <boost/math/distributions/normal.hpp> // for normal_distribution
22 using boost::math::normal; // typedef provides default type is double.
23 #include <boost/math/distributions/find_scale.hpp>
24 using boost::math::find_scale;
25 using boost::math::complement; // Needed if you want to use the complement version.
26 using boost::math::policies::policy; // Needed to specify the error handling policy.
27
28 #include <iostream>
29 using std::cout; using std::endl;
30 #include <iomanip>
31 using std::setw; using std::setprecision;
32 #include <limits>
33 using std::numeric_limits;
34 //] [/find_scale1]
35
main()36 int main()
37 {
38 cout << "Example: Find scale (standard deviation)." << endl;
39 try
40 {
41 //[find_scale2
42 /*`
43 For this example, we will use the standard __normal_distrib,
44 with location (mean) zero and standard deviation (scale) unity.
45 Conveniently, this is also the default for this implementation's constructor.
46 */
47 normal N01; // Default 'standard' normal distribution with zero mean
48 double sd = 1.; // and standard deviation is 1.
49 /*`Suppose we want to find a different normal distribution with standard deviation
50 so that only fraction p (here 0.001 or 0.1%) are below a certain chosen limit
51 (here -2. standard deviations).
52 */
53 double z = -2.; // z to give prob p
54 double p = 0.001; // only 0.1% below z = -2
55
56 cout << "Normal distribution with mean = " << N01.location() // aka N01.mean()
57 << ", standard deviation " << N01.scale() // aka N01.standard_deviation()
58 << ", has " << "fraction <= " << z
59 << ", p = " << cdf(N01, z) << endl;
60 cout << "Normal distribution with mean = " << N01.location()
61 << ", standard deviation " << N01.scale()
62 << ", has " << "fraction > " << z
63 << ", p = " << cdf(complement(N01, z)) << endl; // Note: uses complement.
64 /*`
65 [pre
66 Normal distribution with mean = 0 has fraction <= -2, p = 0.0227501
67 Normal distribution with mean = 0 has fraction > -2, p = 0.97725
68 ]
69 Noting that p = 0.02 instead of our target of 0.001,
70 we can now use `find_scale` to give a new standard deviation.
71 */
72 double l = N01.location();
73 double s = find_scale<normal>(z, p, l);
74 cout << "scale (standard deviation) = " << s << endl;
75 /*`
76 that outputs:
77 [pre
78 scale (standard deviation) = 0.647201
79 ]
80 showing that we need to reduce the standard deviation from 1. to 0.65.
81
82 Then we can check that we have achieved our objective
83 by constructing a new distribution
84 with the new standard deviation (but same zero mean):
85 */
86 normal np001pc(N01.location(), s);
87 /*`
88 And re-calculating the fraction below (and above) our chosen limit.
89 */
90 cout << "Normal distribution with mean = " << l
91 << " has " << "fraction <= " << z
92 << ", p = " << cdf(np001pc, z) << endl;
93 cout << "Normal distribution with mean = " << l
94 << " has " << "fraction > " << z
95 << ", p = " << cdf(complement(np001pc, z)) << endl;
96 /*`
97 [pre
98 Normal distribution with mean = 0 has fraction <= -2, p = 0.001
99 Normal distribution with mean = 0 has fraction > -2, p = 0.999
100 ]
101
102 [h4 Controlling how Errors from find_scale are handled]
103 We can also control the policy for handling various errors.
104 For example, we can define a new (possibly unwise)
105 policy to ignore domain errors ('bad' arguments).
106
107 Unless we are using the boost::math namespace, we will need:
108 */
109 using boost::math::policies::policy;
110 using boost::math::policies::domain_error;
111 using boost::math::policies::ignore_error;
112
113 /*`
114 Using a typedef is convenient, especially if it is re-used,
115 although it is not required, as the various examples below show.
116 */
117 typedef policy<domain_error<ignore_error> > ignore_domain_policy;
118 // find_scale with new policy, using typedef.
119 l = find_scale<normal>(z, p, l, ignore_domain_policy());
120 // Default policy policy<>, needs using boost::math::policies::policy;
121
122 l = find_scale<normal>(z, p, l, policy<>());
123 // Default policy, fully specified.
124 l = find_scale<normal>(z, p, l, boost::math::policies::policy<>());
125 // New policy, without typedef.
126 l = find_scale<normal>(z, p, l, policy<domain_error<ignore_error> >());
127 /*`
128 If we want to express a probability, say 0.999, that is a complement, `1 - p`
129 we should not even think of writing `find_scale<normal>(z, 1 - p, l)`,
130 but use the __complements version (see __why_complements).
131 */
132 z = -2.;
133 double q = 0.999; // = 1 - p; // complement of 0.001.
134 sd = find_scale<normal>(complement(z, q, l));
135
136 normal np95pc(l, sd); // Same standard_deviation (scale) but with mean(scale) shifted
137 cout << "Normal distribution with mean = " << l << " has "
138 << "fraction <= " << z << " = " << cdf(np95pc, z) << endl;
139 cout << "Normal distribution with mean = " << l << " has "
140 << "fraction > " << z << " = " << cdf(complement(np95pc, z)) << endl;
141
142 /*`
143 Sadly, it is all too easy to get probabilities the wrong way round,
144 when you may get a warning like this:
145 [pre
146 Message from thrown exception was:
147 Error in function boost::math::find_scale<Dist, Policy>(complement(double, double, double, Policy)):
148 Computed scale (-0.48043523852179076) is <= 0! Was the complement intended?
149 ]
150 The default error handling policy is to throw an exception with this message,
151 but if you chose a policy to ignore the error,
152 the (impossible) negative scale is quietly returned.
153 */
154 //] [/find_scale2]
155 }
156 catch(const std::exception& e)
157 { // Always useful to include try & catch blocks because default policies
158 // are to throw exceptions on arguments that cause errors like underflow, overflow.
159 // Lacking try & catch blocks, the program will abort without a message below,
160 // which may give some helpful clues as to the cause of the exception.
161 std::cout <<
162 "\n""Message from thrown exception was:\n " << e.what() << std::endl;
163 }
164 return 0;
165 } // int main()
166
167 //[find_scale_example_output
168 /*`
169 [pre
170 Example: Find scale (standard deviation).
171 Normal distribution with mean = 0, standard deviation 1, has fraction <= -2, p = 0.0227501
172 Normal distribution with mean = 0, standard deviation 1, has fraction > -2, p = 0.97725
173 scale (standard deviation) = 0.647201
174 Normal distribution with mean = 0 has fraction <= -2, p = 0.001
175 Normal distribution with mean = 0 has fraction > -2, p = 0.999
176 Normal distribution with mean = 0.946339 has fraction <= -2 = 0.001
177 Normal distribution with mean = 0.946339 has fraction > -2 = 0.999
178 ]
179 */
180 //] [/find_scale_example_output]
181