]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
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 | ||
34 | ||
35 | //] //[/root_find1] | |
36 | ||
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 | */ |