123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169 |
- // find_root_example.cpp
- // Copyright Paul A. Bristow 2007, 2010.
- // Use, modification and distribution are subject to the
- // Boost Software License, Version 1.0.
- // (See accompanying file LICENSE_1_0.txt
- // or copy at http://www.boost.org/LICENSE_1_0.txt)
- // Example of using root finding.
- // Note that this file contains Quickbook mark-up as well as code
- // and comments, don't change any of the special comment mark-ups!
- //[root_find1
- /*`
- First we need some includes to access the normal distribution
- (and some std output of course).
- */
- #include <boost/math/tools/roots.hpp> // root finding.
- #include <boost/math/distributions/normal.hpp> // for normal_distribution
- using boost::math::normal; // typedef provides default type is double.
- #include <iostream>
- using std::cout; using std::endl; using std::left; using std::showpoint; using std::noshowpoint;
- #include <iomanip>
- using std::setw; using std::setprecision;
- #include <limits>
- using std::numeric_limits;
- #include <stdexcept>
-
- //] //[/root_find1]
- int main()
- {
- cout << "Example: Normal distribution, root finding.";
- try
- {
- //[root_find2
- /*`A machine is set to pack 3 kg of ground beef per pack.
- Over a long period of time it is found that the average packed was 3 kg
- with a standard deviation of 0.1 kg.
- Assuming the packing is normally distributed,
- we can find the fraction (or %) of packages that weigh more than 3.1 kg.
- */
- double mean = 3.; // kg
- double standard_deviation = 0.1; // kg
- normal packs(mean, standard_deviation);
- double max_weight = 3.1; // kg
- cout << "Percentage of packs > " << max_weight << " is "
- << cdf(complement(packs, max_weight)) << endl; // P(X > 3.1)
- double under_weight = 2.9;
- cout <<"fraction of packs <= " << under_weight << " with a mean of " << mean
- << " is " << cdf(complement(packs, under_weight)) << endl;
- // fraction of packs <= 2.9 with a mean of 3 is 0.841345
- // This is 0.84 - more than the target 0.95
- // Want 95% to be over this weight, so what should we set the mean weight to be?
- // KK StatCalc says:
- double over_mean = 3.0664;
- normal xpacks(over_mean, standard_deviation);
- cout << "fraction of packs >= " << under_weight
- << " with a mean of " << xpacks.mean()
- << " is " << cdf(complement(xpacks, under_weight)) << endl;
- // fraction of packs >= 2.9 with a mean of 3.06449 is 0.950005
- double under_fraction = 0.05; // so 95% are above the minimum weight mean - sd = 2.9
- double low_limit = standard_deviation;
- double offset = mean - low_limit - quantile(packs, under_fraction);
- double nominal_mean = mean + offset;
- normal nominal_packs(nominal_mean, standard_deviation);
- cout << "Setting the packer to " << nominal_mean << " will mean that "
- << "fraction of packs >= " << under_weight
- << " is " << cdf(complement(nominal_packs, under_weight)) << endl;
- /*`
- Setting the packer to 3.06449 will mean that fraction of packs >= 2.9 is 0.95.
- Setting the packer to 3.13263 will mean that fraction of packs >= 2.9 is 0.99,
- but will more than double the mean loss from 0.0644 to 0.133.
- Alternatively, we could invest in a better (more precise) packer with a lower standard deviation.
- To estimate how much better (how much smaller standard deviation) it would have to be,
- we need to get the 5% quantile to be located at the under_weight limit, 2.9
- */
- double p = 0.05; // wanted p th quantile.
- cout << "Quantile of " << p << " = " << quantile(packs, p)
- << ", mean = " << packs.mean() << ", sd = " << packs.standard_deviation() << endl; //
- /*`
- Quantile of 0.05 = 2.83551, mean = 3, sd = 0.1
- With the current packer (mean = 3, sd = 0.1), the 5% quantile is at 2.8551 kg,
- a little below our target of 2.9 kg.
- So we know that the standard deviation is going to have to be smaller.
- Let's start by guessing that it (now 0.1) needs to be halved, to a standard deviation of 0.05
- */
- normal pack05(mean, 0.05);
- cout << "Quantile of " << p << " = " << quantile(pack05, p)
- << ", mean = " << pack05.mean() << ", sd = " << pack05.standard_deviation() << endl;
- cout <<"Fraction of packs >= " << under_weight << " with a mean of " << mean
- << " and standard deviation of " << pack05.standard_deviation()
- << " is " << cdf(complement(pack05, under_weight)) << endl;
- //
- /*`
- Fraction of packs >= 2.9 with a mean of 3 and standard deviation of 0.05 is 0.9772
- So 0.05 was quite a good guess, but we are a little over the 2.9 target,
- so the standard deviation could be a tiny bit more. So we could do some
- more guessing to get closer, say by increasing to 0.06
- */
- normal pack06(mean, 0.06);
- cout << "Quantile of " << p << " = " << quantile(pack06, p)
- << ", mean = " << pack06.mean() << ", sd = " << pack06.standard_deviation() << endl;
- cout <<"Fraction of packs >= " << under_weight << " with a mean of " << mean
- << " and standard deviation of " << pack06.standard_deviation()
- << " is " << cdf(complement(pack06, under_weight)) << endl;
- /*`
- Fraction of packs >= 2.9 with a mean of 3 and standard deviation of 0.06 is 0.9522
- Now we are getting really close, but to do the job properly,
- we could use root finding method, for example the tools provided, and used elsewhere,
- in the Math Toolkit, see __root_finding_without_derivatives.
- But in this normal distribution case, we could be even smarter and make a direct calculation.
- */
- //] [/root_find2]
- }
- catch(const std::exception& e)
- { // Always useful to include try & catch blocks because default policies
- // are to throw exceptions on arguments that cause errors like underflow, overflow.
- // Lacking try & catch blocks, the program will abort without a message below,
- // which may give some helpful clues as to the cause of the exception.
- std::cout <<
- "\n""Message from thrown exception was:\n " << e.what() << std::endl;
- }
- return 0;
- } // int main()
- /*
- Output is:
- //[root_find_output
- Autorun "i:\boost-06-05-03-1300\libs\math\test\Math_test\debug\find_root_example.exe"
- Example: Normal distribution, root finding.Percentage of packs > 3.1 is 0.158655
- fraction of packs <= 2.9 with a mean of 3 is 0.841345
- fraction of packs >= 2.9 with a mean of 3.0664 is 0.951944
- Setting the packer to 3.06449 will mean that fraction of packs >= 2.9 is 0.95
- Quantile of 0.05 = 2.83551, mean = 3, sd = 0.1
- Quantile of 0.05 = 2.91776, mean = 3, sd = 0.05
- Fraction of packs >= 2.9 with a mean of 3 and standard deviation of 0.05 is 0.97725
- Quantile of 0.05 = 2.90131, mean = 3, sd = 0.06
- Fraction of packs >= 2.9 with a mean of 3 and standard deviation of 0.06 is 0.95221
- //] [/root_find_output]
- */
|