1 // find_root_example.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 using root finding.
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 //[root_find1
16 /*`
17 First we need some includes to access the normal distribution
18 (and some std output of course).
19 */
20
21 #include <boost/math/tools/roots.hpp> // root finding.
22
23 #include <boost/math/distributions/normal.hpp> // for normal_distribution
24 using boost::math::normal; // typedef provides default type is double.
25
26 #include <iostream>
27 using std::cout; using std::endl; using std::left; using std::showpoint; using std::noshowpoint;
28 #include <iomanip>
29 using std::setw; using std::setprecision;
30 #include <limits>
31 using std::numeric_limits;
32 #include <stdexcept>
33 using std::exception;
34
35 //] //[/root_find1]
36
main()37 int main()
38 {
39 cout << "Example: Normal distribution, root finding.";
40 try
41 {
42
43 //[root_find2
44
45 /*`A machine is set to pack 3 kg of ground beef per pack.
46 Over a long period of time it is found that the average packed was 3 kg
47 with a standard deviation of 0.1 kg.
48 Assuming the packing is normally distributed,
49 we can find the fraction (or %) of packages that weigh more than 3.1 kg.
50 */
51
52 double mean = 3.; // kg
53 double standard_deviation = 0.1; // kg
54 normal packs(mean, standard_deviation);
55
56 double max_weight = 3.1; // kg
57 cout << "Percentage of packs > " << max_weight << " is "
58 << cdf(complement(packs, max_weight)) << endl; // P(X > 3.1)
59
60 double under_weight = 2.9;
61 cout <<"fraction of packs <= " << under_weight << " with a mean of " << mean
62 << " is " << cdf(complement(packs, under_weight)) << endl;
63 // fraction of packs <= 2.9 with a mean of 3 is 0.841345
64 // This is 0.84 - more than the target 0.95
65 // Want 95% to be over this weight, so what should we set the mean weight to be?
66 // KK StatCalc says:
67 double over_mean = 3.0664;
68 normal xpacks(over_mean, standard_deviation);
69 cout << "fraction of packs >= " << under_weight
70 << " with a mean of " << xpacks.mean()
71 << " is " << cdf(complement(xpacks, under_weight)) << endl;
72 // fraction of packs >= 2.9 with a mean of 3.06449 is 0.950005
73 double under_fraction = 0.05; // so 95% are above the minimum weight mean - sd = 2.9
74 double low_limit = standard_deviation;
75 double offset = mean - low_limit - quantile(packs, under_fraction);
76 double nominal_mean = mean + offset;
77
78 normal nominal_packs(nominal_mean, standard_deviation);
79 cout << "Setting the packer to " << nominal_mean << " will mean that "
80 << "fraction of packs >= " << under_weight
81 << " is " << cdf(complement(nominal_packs, under_weight)) << endl;
82
83 /*`
84 Setting the packer to 3.06449 will mean that fraction of packs >= 2.9 is 0.95.
85
86 Setting the packer to 3.13263 will mean that fraction of packs >= 2.9 is 0.99,
87 but will more than double the mean loss from 0.0644 to 0.133.
88
89 Alternatively, we could invest in a better (more precise) packer with a lower standard deviation.
90
91 To estimate how much better (how much smaller standard deviation) it would have to be,
92 we need to get the 5% quantile to be located at the under_weight limit, 2.9
93 */
94 double p = 0.05; // wanted p th quantile.
95 cout << "Quantile of " << p << " = " << quantile(packs, p)
96 << ", mean = " << packs.mean() << ", sd = " << packs.standard_deviation() << endl; //
97 /*`
98 Quantile of 0.05 = 2.83551, mean = 3, sd = 0.1
99
100 With the current packer (mean = 3, sd = 0.1), the 5% quantile is at 2.8551 kg,
101 a little below our target of 2.9 kg.
102 So we know that the standard deviation is going to have to be smaller.
103
104 Let's start by guessing that it (now 0.1) needs to be halved, to a standard deviation of 0.05
105 */
106 normal pack05(mean, 0.05);
107 cout << "Quantile of " << p << " = " << quantile(pack05, p)
108 << ", mean = " << pack05.mean() << ", sd = " << pack05.standard_deviation() << endl;
109
110 cout <<"Fraction of packs >= " << under_weight << " with a mean of " << mean
111 << " and standard deviation of " << pack05.standard_deviation()
112 << " is " << cdf(complement(pack05, under_weight)) << endl;
113 //
114 /*`
115 Fraction of packs >= 2.9 with a mean of 3 and standard deviation of 0.05 is 0.9772
116
117 So 0.05 was quite a good guess, but we are a little over the 2.9 target,
118 so the standard deviation could be a tiny bit more. So we could do some
119 more guessing to get closer, say by increasing to 0.06
120 */
121
122 normal pack06(mean, 0.06);
123 cout << "Quantile of " << p << " = " << quantile(pack06, p)
124 << ", mean = " << pack06.mean() << ", sd = " << pack06.standard_deviation() << endl;
125
126 cout <<"Fraction of packs >= " << under_weight << " with a mean of " << mean
127 << " and standard deviation of " << pack06.standard_deviation()
128 << " is " << cdf(complement(pack06, under_weight)) << endl;
129 /*`
130 Fraction of packs >= 2.9 with a mean of 3 and standard deviation of 0.06 is 0.9522
131
132 Now we are getting really close, but to do the job properly,
133 we could use root finding method, for example the tools provided, and used elsewhere,
134 in the Math Toolkit, see __root_finding_without_derivatives.
135
136 But in this normal distribution case, we could be even smarter and make a direct calculation.
137 */
138 //] [/root_find2]
139
140 }
141 catch(const std::exception& e)
142 { // Always useful to include try & catch blocks because default policies
143 // are to throw exceptions on arguments that cause errors like underflow, overflow.
144 // Lacking try & catch blocks, the program will abort without a message below,
145 // which may give some helpful clues as to the cause of the exception.
146 std::cout <<
147 "\n""Message from thrown exception was:\n " << e.what() << std::endl;
148 }
149 return 0;
150 } // int main()
151
152 /*
153 Output is:
154
155 //[root_find_output
156
157 Autorun "i:\boost-06-05-03-1300\libs\math\test\Math_test\debug\find_root_example.exe"
158 Example: Normal distribution, root finding.Percentage of packs > 3.1 is 0.158655
159 fraction of packs <= 2.9 with a mean of 3 is 0.841345
160 fraction of packs >= 2.9 with a mean of 3.0664 is 0.951944
161 Setting the packer to 3.06449 will mean that fraction of packs >= 2.9 is 0.95
162 Quantile of 0.05 = 2.83551, mean = 3, sd = 0.1
163 Quantile of 0.05 = 2.91776, mean = 3, sd = 0.05
164 Fraction of packs >= 2.9 with a mean of 3 and standard deviation of 0.05 is 0.97725
165 Quantile of 0.05 = 2.90131, mean = 3, sd = 0.06
166 Fraction of packs >= 2.9 with a mean of 3 and standard deviation of 0.06 is 0.95221
167
168 //] [/root_find_output]
169 */
170