]>
git.proxmox.com Git - ceph.git/blob - ceph/src/boost/libs/math/example/find_root_example.cpp
1 // find_root_example.cpp
3 // Copyright Paul A. Bristow 2007, 2010.
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)
10 // Example of using root finding.
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!
17 First we need some includes to access the normal distribution
18 (and some std output of course).
21 #include <boost/math/tools/roots.hpp> // root finding.
23 #include <boost/math/distributions/normal.hpp> // for normal_distribution
24 using boost::math::normal
; // typedef provides default type is double.
27 using std::cout
; using std::endl
; using std::left
; using std::showpoint
; using std::noshowpoint
;
29 using std::setw
; using std::setprecision
;
31 using std::numeric_limits
;
39 cout
<< "Example: Normal distribution, root finding.";
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.
52 double mean
= 3.; // kg
53 double standard_deviation
= 0.1; // kg
54 normal
packs(mean
, standard_deviation
);
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)
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?
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
;
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
;
84 Setting the packer to 3.06449 will mean that fraction of packs >= 2.9 is 0.95.
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.
89 Alternatively, we could invest in a better (more precise) packer with a lower standard deviation.
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
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
; //
98 Quantile of 0.05 = 2.83551, mean = 3, sd = 0.1
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.
104 Let's start by guessing that it (now 0.1) needs to be halved, to a standard deviation of 0.05
106 normal
pack05(mean
, 0.05);
107 cout
<< "Quantile of " << p
<< " = " << quantile(pack05
, p
)
108 << ", mean = " << pack05
.mean() << ", sd = " << pack05
.standard_deviation() << endl
;
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
;
115 Fraction of packs >= 2.9 with a mean of 3 and standard deviation of 0.05 is 0.9772
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
122 normal
pack06(mean
, 0.06);
123 cout
<< "Quantile of " << p
<< " = " << quantile(pack06
, p
)
124 << ", mean = " << pack06
.mean() << ", sd = " << pack06
.standard_deviation() << endl
;
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
;
130 Fraction of packs >= 2.9 with a mean of 3 and standard deviation of 0.06 is 0.9522
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.
136 But in this normal distribution case, we could be even smarter and make a direct calculation.
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.
147 "\n""Message from thrown exception was:\n " << e
.what() << std::endl
;
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
168 //] [/root_find_output]