
BoostCommit : 
From: pbristow_at_[hidden]
Date: 20070916 11:27:51
Author: pbristow
Date: 20070916 11:27:50 EDT (Sun, 16 Sep 2007)
New Revision: 39325
URL: http://svn.boost.org/trac/boost/changeset/39325
Log:
1st commit
Added:
sandbox/math_toolkit/libs/math/example/find_location_example.cpp (contents, props changed)
sandbox/math_toolkit/libs/math/example/find_mean_and_sd_normal.cpp (contents, props changed)
sandbox/math_toolkit/libs/math/example/find_root_example.cpp (contents, props changed)
sandbox/math_toolkit/libs/math/example/find_scale_example.cpp (contents, props changed)
Added: sandbox/math_toolkit/libs/math/example/find_location_example.cpp
==============================================================================
 (empty file)
+++ sandbox/math_toolkit/libs/math/example/find_location_example.cpp 20070916 11:27:50 EDT (Sun, 16 Sep 2007)
@@ 0,0 +1,170 @@
+// find_location.cpp
+
+// Copyright Paul A. Bristow 2007.
+
+// 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 finding location (mean)
+// for normal (Gaussian) & Cauchy distribution.
+
+// Note that this file contains Quickbook markup as well as code
+// and comments, don't change any of the special comment markups!
+
+//[find_location1
+/*`
+First we need some includes to access the normal distribution,
+the algorithms to find location (and some std output of course).
+*/
+
+#include <boost/math/distributions/normal.hpp> // for normal_distribution
+ using boost::math::normal; // typedef provides default type is double.
+#include <boost/math/distributions/cauchy.hpp> // for cauchy_distribution
+ using boost::math::cauchy; // typedef provides default type is double.
+#include <boost/math/distributions/find_location.hpp>
+ using boost::math::find_location;
+ using boost::math::complement; // Needed if you want to use the complement version.
+ using boost::math::policies::policy;
+
+#include <iostream>
+ using std::cout; using std::endl;
+#include <iomanip>
+ using std::setw; using std::setprecision;
+#include <limits>
+ using std::numeric_limits;
+//] [/find_location1]
+
+int main()
+{
+ cout << "Example: Find location (mean)." << endl;
+ try
+ {
+//[find_location2
+/*`
+For this example, we will use the standard normal distribution,
+with mean (location) zero and standard deviation (scale) unity.
+This is also the default for this implementation.
+*/
+ normal N01; // Default 'standard' normal distribution with zero mean and
+ double sd = 1.; // normal default standard deviation is 1.
+/*`Suppose we want to find a different normal distribution whose mean is shifted
+so that only fraction p (here 0.001 or 0.1%) are below a certain chosen limit
+(here 2, two standard deviations).
+*/
+ double z = 2.; // z to give prob p
+ double p = 0.001; // only 0.1% below z
+
+ cout << "Normal distribution with mean = " << N01.location()
+ << ", standard deviation " << N01.scale()
+ << ", has " << "fraction <= " << z
+ << ", p = " << cdf(N01, z) << endl;
+ cout << "Normal distribution with mean = " << N01.location()
+ << ", standard deviation " << N01.scale()
+ << ", has " << "fraction > " << z
+ << ", p = " << cdf(complement(N01, z)) << endl; // Note: uses complement.
+/*`
+[pre
+Normal distribution with mean = 0, standard deviation 1, has fraction <= 2, p = 0.0227501
+Normal distribution with mean = 0, standard deviation 1, has fraction > 2, p = 0.97725
+]
+We can now use ''find_location'' to give a new offset mean.
+*/
+ double l = find_location<normal>(z, p, sd);
+ cout << "offset location (mean) = " << l << endl;
+/*`
+that outputs:
+[pre
+offset location (mean) = 1.09023
+]
+showing that we need to shift the mean just over one standard deviation from its previous value of zero.
+
+Then we can check that we have achieved our objective
+by constructing a new distribution
+with the offset mean (but same standard deviation):
+*/
+ normal np001pc(l, sd); // Same standard_deviation (scale) but with mean (location) shifted.
+/*`
+And recalculating the fraction below our chosen limit.
+*/
+cout << "Normal distribution with mean = " << l
+ << " has " << "fraction <= " << z
+ << ", p = " << cdf(np001pc, z) << endl;
+ cout << "Normal distribution with mean = " << l
+ << " has " << "fraction > " << z
+ << ", p = " << cdf(complement(np001pc, z)) << endl;
+/*`
+[pre
+Normal distribution with mean = 1.09023 has fraction <= 2, p = 0.001
+Normal distribution with mean = 1.09023 has fraction > 2, p = 0.999
+]
+
+[h4 Controlling Error Handling from find_location]
+We can also control the policy for handling various errors.
+For example, we can define a new (possibly unwise)
+policy to ignore domain errors ('bad' arguments).
+
+Unless we are using the boost::math namespace, we will need:
+*/
+ using boost::math::policies::policy;
+ using boost::math::policies::domain_error;
+ using boost::math::policies::ignore_error;
+
+/*`
+Using a typedef is often convenient, especially if it is reused,
+although it is not required, as the various examples below show.
+*/
+ typedef policy<domain_error<ignore_error> > ignore_domain_policy;
+ // find_location with new policy, using typedef.
+ l = find_location<normal>(z, p, sd, ignore_domain_policy());
+ // Default policy policy<>, needs "using boost::math::policies::policy;"
+ l = find_location<normal>(z, p, sd, policy<>());
+ // Default policy, fully specified.
+ l = find_location<normal>(z, p, sd, boost::math::policies::policy<>());
+ // A new policy, ignoring domain errors, without using a typedef.
+ l = find_location<normal>(z, p, sd, policy<domain_error<ignore_error> >());
+/*`
+If we want to use a probability that is the
+[link complements complement of our probability],
+we should not even think of writing `find_location<normal>(z, 1  p, sd)`,
+but, [link why_complements to avoid loss of accuracy], use the complement version.
+*/
+ z = 2.;
+ double q = 0.95; // = 1  p; // complement.
+ l = find_location<normal>(complement(z, q, sd));
+
+ normal np95pc(l, sd); // Same standard_deviation (scale) but with mean(location) shifted
+ cout << "Normal distribution with mean = " << l << " has "
+ << "fraction <= " << z << " = " << cdf(np95pc, z) << endl;
+ cout << "Normal distribution with mean = " << l << " has "
+ << "fraction > " << z << " = " << cdf(complement(np95pc, z)) << endl;
+ //] [/find_location2]
+ }
+ 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()
+
+
+
+//[find_location_example_output
+/*`
+[pre
+Example: Find location (mean).
+Normal distribution with mean = 0, standard deviation 1, has fraction <= 2, p = 0.0227501
+Normal distribution with mean = 0, standard deviation 1, has fraction > 2, p = 0.97725
+offset location (mean) = 1.09023
+Normal distribution with mean = 1.09023 has fraction <= 2, p = 0.001
+Normal distribution with mean = 1.09023 has fraction > 2, p = 0.999
+Normal distribution with mean = 0.355146 has fraction <= 2 = 0.95
+Normal distribution with mean = 0.355146 has fraction > 2 = 0.05
+]
+*/
+//] [/find_location_example_output]
Added: sandbox/math_toolkit/libs/math/example/find_mean_and_sd_normal.cpp
==============================================================================
 (empty file)
+++ sandbox/math_toolkit/libs/math/example/find_mean_and_sd_normal.cpp 20070916 11:27:50 EDT (Sun, 16 Sep 2007)
@@ 0,0 +1,409 @@
+// find_mean_and_sd_normal.cpp
+
+// Copyright Paul A. Bristow 2007.
+
+// 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 finding mean or sd for normal distribution.
+
+// Note that this file contains Quickbook markup as well as code
+// and comments, don't change any of the special comment markups!
+
+//[normal_std
+/*`
+First we need some includes to access the normal distribution,
+the algorithms to find location and scale
+(and some std output of course).
+*/
+
+#include <boost/math/distributions/normal.hpp> // for normal_distribution
+ using boost::math::normal; // typedef provides default type is double.
+#include <boost/math/distributions/cauchy.hpp> // for cauchy_distribution
+ using boost::math::cauchy; // typedef provides default type is double.
+#include <boost/math/distributions/find_location.hpp>
+ using boost::math::find_location;
+#include <boost/math/distributions/find_scale.hpp>
+ using boost::math::find_scale;
+ using boost::math::complement;
+ using boost::math::policies::policy;
+
+#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;
+//] [/normal_std Quickbook]
+
+int main()
+{
+ cout << "Find_location (mean) and find_scale (standard deviation) examples." << endl;
+ try
+ {
+
+//[normal_find_location_and_scale_eg
+
+/*`
+[h4 Using find_location and find_scale to meet dispensing and measurement specifications]
+
+Consider an example from K Krishnamoorthy,
+Handbook of Statistical Distributions with Applications,
+ISBN 1584886358, (2006) p 126, example 10.3.7.
+
+"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.
+Assume the packing is normally distributed."
+
+We start by constructing a normal distribution with the given parameters:
+*/
+
+double mean = 3.; // kg
+double standard_deviation = 0.1; // kg
+normal packs(mean, standard_deviation);
+/*`We can then find the fraction (or %) of packages that weigh more than 3.1 kg.
+*/
+
+double max_weight = 3.1; // kg
+cout << "Percentage of packs > " << max_weight << " is "
+<< cdf(complement(packs, max_weight)) * 100. << endl; // P(X > 3.1)
+
+/*`We might want to ensure that 95% of packs are over a minimum weight specification,
+then we want the value of the mean such that P(X < 2.9) = 0.05.
+
+Using the mean of 3 kg, we can estimate
+the fraction of packs that fail to meet the specification of 2.9 kg.
+*/
+
+double minimum_weight = 2.9;
+cout <<"Fraction of packs <= " << minimum_weight << " with a mean of " << mean
+ << " is " << cdf(complement(packs, minimum_weight)) << endl;
+// fraction of packs <= 2.9 with a mean of 3 is 0.841345
+
+/*`This is 0.84  more than the target fraction of 0.95.
+If we want 95% to be over the minimum weight,
+what should we set the mean weight to be?
+
+Using the KK StatCalc program supplied with the book and the method given on page 126 gives 3.06449.
+
+We can confirm this by constructing a new distribution which we call 'xpacks'
+with a safety margin mean of 3.06449 thus:
+*/
+double over_mean = 3.06449;
+normal xpacks(over_mean, standard_deviation);
+cout << "Fraction of packs >= " << minimum_weight
+<< " with a mean of " << xpacks.mean()
+ << " is " << cdf(complement(xpacks, minimum_weight)) << endl;
+// fraction of packs >= 2.9 with a mean of 3.06449 is 0.950005
+
+/*`Using this Math Toolkit, we can calculate the required mean directly thus:
+*/
+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;
+// mean + (mean  low_limit  quantile(packs, under_fraction));
+
+normal nominal_packs(nominal_mean, standard_deviation);
+cout << "Setting the packer to " << nominal_mean << " will mean that "
+ << "fraction of packs >= " << minimum_weight
+ << " is " << cdf(complement(nominal_packs, minimum_weight)) << endl;
+// Setting the packer to 3.06449 will mean that fraction of packs >= 2.9 is 0.95
+
+/*`This calculation is generalized as the free function called find_location.
+
+
+TODO link to find_location docs
+To use this we will need to
+*/
+
+#include <boost/math/distributions/find_location.hpp>
+ using boost::math::find_location;
+/*`and then use find_location function to find safe_mean,
+& construct a new normal distribution called 'goodpacks'.
+*/
+double safe_mean = find_location<normal>(minimum_weight, under_fraction, standard_deviation);
+normal good_packs(safe_mean, standard_deviation);
+/*`with the same confirmation as before:
+*/
+cout << "Setting the packer to " << nominal_mean << " will mean that "
+ << "fraction of packs >= " << minimum_weight
+ << " is " << cdf(complement(good_packs, minimum_weight)) << endl;
+// Setting the packer to 3.06449 will mean that fraction of packs >= 2.9 is 0.95
+
+/*`
+[h4 Using CauchyLorentz instead of normal distribution]
+
+After examining the weight distribution of a large number of packs, we might decide that,
+after all, the assumption of a normal distribution is not really justified.
+We might find that the fit is better to a
+[link __cauchy_distrib CauchyLorentz distribution].
+This distribution has wider 'wings', so that whereas most of the values
+are closer to the mean than the normal, there are also more values than 'normal'
+that lie further from the mean than the normal.
+
+This might happen because a larger than normal lump of meat is either included or excluded.
+
+We first create a __cauchy_distrib with the original mean and standard deviation,
+and estimate the fraction that lie below our minimum weight specification.
+*/
+
+cauchy cpacks(mean, standard_deviation);
+cout << "Cauchy Setting the packer to " << mean << " will mean that "
+ << "fraction of packs >= " << minimum_weight
+ << " is " << cdf(complement(cpacks, minimum_weight)) << endl;
+// Cauchy Setting the packer to 3 will mean that fraction of packs >= 2.9 is 0.75
+
+/*`Note that far fewer of the packs meet the specification, only 75% instead of 95%.
+Now we can repeat the find_location, using the cauchy distribution as template parameter,
+in place of the normal used above.
+*/
+
+double lc = find_location<cauchy>(minimum_weight, under_fraction, standard_deviation);
+cout << "find_location<cauchy>(minimum_weight, over fraction, standard_deviation); " << lc << endl;
+// find_location<cauchy>(minimum_weight, over fraction, packs.standard_deviation()); 3.53138
+/*`Note that the safe_mean setting needs to be much higher, 3.53138 instead of 3.06449,
+so we will make rather less profit.
+
+And again confirm that the fraction meeting specification is as expected.
+*/
+cauchy goodcpacks(lc, standard_deviation);
+cout << "Cauchy Setting the packer to " << lc << " will mean that "
+ << "fraction of packs >= " << minimum_weight
+ << " is " << cdf(complement(goodcpacks, minimum_weight)) << endl;
+// Cauchy Setting the packer to 3.53138 will mean that fraction of packs >= 2.9 is 0.95
+
+/*`Finally we could estimate the effect of a much tighter specification,
+that 99% of packs met the specification.
+*/
+
+cout << "Cauchy Setting the packer to "
+ << find_location<cauchy>(minimum_weight, 0.99, standard_deviation)
+ << " will mean that "
+ << "fraction of packs >= " << minimum_weight
+ << " is " << cdf(complement(goodcpacks, minimum_weight)) << endl;
+
+/*`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 kg per pack.
+
+Of course, this calculation is not limited to packs of meat, it applies to dispensing anything,
+and it also applies to a 'virtual' material like any measurement.
+
+The only caveat is that the calculation assumes that the standard deviation (scale) is known with
+a reasonably low uncertainty, something that is not so easy to ensure in practice.
+And that the distribution is well defined, __normal_distrib or __cauchy_distrib, or some other.
+
+If one is simply dispensing a very large number of packs,
+then it may be feasible to measure the weight of hundreds or thousands of packs.
+With a healthy 'degrees of freedom', the confidence intervals for the standard deviation
+are not too wide, typically about + and  10% for hundreds of observations.
+
+For other applications, where it is more difficult or expensive to make many observations,
+the confidence intervals are depressingly wide.
+
+See [link math_toolkit.dist.stat_tut.weg.cs_eg.chi_sq_intervals Confidence Intervals on the standard deviation]
+for a worked example
+[@../../example/chi_square_std_dev_test.cpp chi_square_std_dev_test.cpp]
+of estimating these intervals.
+
+
+[h4 Changing the scale or standard deviation]
+
+Alternatively, we could invest in a better (more precise) packer
+(or measuring device) with a lower standard deviation, or scale.
+
+This might cost more, but would reduce the amount we have to 'give away'
+in order to meet the specification.
+
+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 kg.
+*/
+normal pack05(mean, 0.05);
+cout << "Quantile of " << p << " = " << quantile(pack05, p)
+ << ", mean = " << pack05.mean() << ", sd = " << pack05.standard_deviation() << endl;
+// Quantile of 0.05 = 2.91776, mean = 3, sd = 0.05
+
+cout <<"Fraction of packs >= " << minimum_weight << " with a mean of " << mean
+ << " and standard deviation of " << pack05.standard_deviation()
+ << " is " << cdf(complement(pack05, minimum_weight)) << endl;
+// Fraction of packs >= 2.9 with a mean of 3 and standard deviation of 0.05 is 0.97725
+/*`
+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 standard deviation to 0.06 kg,
+constructing another new distribution called pack06.
+*/
+normal pack06(mean, 0.06);
+cout << "Quantile of " << p << " = " << quantile(pack06, p)
+ << ", mean = " << pack06.mean() << ", sd = " << pack06.standard_deviation() << endl;
+// Quantile of 0.05 = 2.90131, mean = 3, sd = 0.06
+
+cout <<"Fraction of packs >= " << minimum_weight << " with a mean of " << mean
+ << " and standard deviation of " << pack06.standard_deviation()
+ << " is " << cdf(complement(pack06, minimum_weight)) << endl;
+// Fraction of packs >= 2.9 with a mean of 3 and standard deviation of 0.06 is 0.95221
+/*`
+Now we are getting really close, but to do the job properly,
+we might need to use root finding method, for example the tools provided,
+and used elsewhere, in the Math Toolkit, see
+[link math_toolkit.toolkit.internals1.roots2 Root Finding Without Derivatives].
+
+But in this (normal) distribution case, we can and should be even smarter
+and make a direct calculation.
+*/
+
+/*`Our required limit is minimum_weight = 2.9 kg, often called the random variate z.
+For a standard normal distribution, then probability p = N((minimum_weight  mean) / sd).
+
+We want to find the standard deviation that would be required to meet this limit,
+so that the p th quantile is located at z (minimum_weight).
+In this case, the 0.05 (5%) quantile is at 2.9 kg pack weight, when the mean is 3 kg,
+ensuring that 0.95 (95%) of packs are above the minimum weight.
+
+Rearranging, we can directly calculate the required standard deviation:
+*/
+normal N01; // standard normal distribution with meamn zero and unit standard deviation.
+p = 0.05;
+double qp = quantile(N01, p);
+double sd95 = (minimum_weight  mean) / qp;
+
+cout << "For the "<< p << "th quantile to be located at "
+ << minimum_weight << ", would need a standard deviation of " << sd95 << endl;
+// For the 0.05th quantile to be located at 2.9, would need a standard deviation of 0.0607957
+
+/*`We can now construct a new (normal) distribution pack95 for the 'better' packer,
+and check that our distribution will meet the specification.
+*/
+
+normal pack95(mean, sd95);
+cout <<"Fraction of packs >= " << minimum_weight << " with a mean of " << mean
+ << " and standard deviation of " << pack95.standard_deviation()
+ << " is " << cdf(complement(pack95, minimum_weight)) << endl;
+// Fraction of packs >= 2.9 with a mean of 3 and standard deviation of 0.0607957 is 0.95
+
+/*`This calculation is generalized in the free function find_scale,
+as shown below, giving the same standard deviation.
+*/
+double ss = find_scale<normal>(minimum_weight, under_fraction, packs.mean());
+cout << "find_scale<normal>(minimum_weight, under_fraction, packs.mean()); " << ss << endl;
+// find_scale<normal>(minimum_weight, under_fraction, packs.mean()); 0.0607957
+
+/*`If we had defined an over_fraction, or percentage that must pass specification
+*/
+double over_fraction = 0.95;
+/*`And (wrongly) written
+
+ double sso = find_scale<normal>(minimum_weight, over_fraction, packs.mean());
+
+With the default policy, we would get a message like
+
+[pre
+Message from thrown exception was:
+ Error in function boost::math::find_scale<Dist, Policy>(double, double, double, Policy):
+ Computed scale (0.060795683191176959) is <= 0! Was the complement intended?
+]
+
+But this would return a *negative* standard deviation  obviously impossible.
+The probability should be 1  over_fraction, not over_fraction, thus:
+*/
+
+double ss1o = find_scale<normal>(minimum_weight, 1  over_fraction, packs.mean());
+cout << "find_scale<normal>(minimum_weight, under_fraction, packs.mean()); " << ss1o << endl;
+// find_scale<normal>(minimum_weight, under_fraction, packs.mean()); 0.0607957
+
+/*`But notice that using '1  over_fraction'  will lead to a
+[link why_complements loss of accuracy, especially if over_fraction was close to unity.]
+In this (very common) case, we should instead use the __complements,
+giving the most accurate result.
+*/
+
+double ssc = find_scale<normal>(complement(minimum_weight, over_fraction, packs.mean()));
+cout << "find_scale<normal>(complement(minimum_weight, over_fraction, packs.mean())); " << ssc << endl;
+// find_scale<normal>(complement(minimum_weight, over_fraction, packs.mean())); 0.0607957
+
+/*`Note that our guess of 0.06 was close to the accurate value of 0.060795683191176959.
+
+We can again confirm our prediction thus:
+*/
+
+normal pack95c(mean, ssc);
+cout <<"Fraction of packs >= " << minimum_weight << " with a mean of " << mean
+ << " and standard deviation of " << pack95c.standard_deviation()
+ << " is " << cdf(complement(pack95c, minimum_weight)) << endl;
+// Fraction of packs >= 2.9 with a mean of 3 and standard deviation of 0.0607957 is 0.95
+
+/*`Notice that these two deceptively simple questions:
+
+* Do we overfill to make sure we meet a minimum specification (or underfill to avoid an overdose)?
+
+and/or
+
+* Do we measure better?
+
+are actually extremely common.
+
+The weight of beef might be replaced by a measurement of more or less anything,
+from drug tablet content, Apollo landing rocket firing, Xray treatment doses...
+
+The scale can be variation in dispensing or uncertainty in measurement.
+*/
+//] [/normal_find_location_and_scale_eg Quickbook end]
+
+ }
+ 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:
+
+Find_location and find_scale examples.
+Percentage of packs > 3.1 is 15.8655
+Fraction of packs <= 2.9 with a mean of 3 is 0.841345
+Fraction of packs >= 2.9 with a mean of 3.06449 is 0.950005
+Setting the packer to 3.06449 will mean that fraction of packs >= 2.9 is 0.95
+Setting the packer to 3.06449 will mean that fraction of packs >= 2.9 is 0.95
+Cauchy Setting the packer to 3 will mean that fraction of packs >= 2.9 is 0.75
+find_location<cauchy>(minimum_weight, over fraction, standard_deviation); 3.53138
+Cauchy Setting the packer to 3.53138 will mean that fraction of packs >= 2.9 is 0.95
+Cauchy Setting the packer to 0.282052 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
+For the 0.05th quantile to be located at 2.9, would need a standard deviation of 0.0607957
+Fraction of packs >= 2.9 with a mean of 3 and standard deviation of 0.0607957 is 0.95
+find_scale<normal>(minimum_weight, under_fraction, packs.mean()); 0.0607957
+find_scale<normal>(minimum_weight, under_fraction, packs.mean()); 0.0607957
+find_scale<normal>(complement(minimum_weight, over_fraction, packs.mean())); 0.0607957
+Fraction of packs >= 2.9 with a mean of 3 and standard deviation of 0.0607957 is 0.95
+
+*/
+
+
+
Added: sandbox/math_toolkit/libs/math/example/find_root_example.cpp
==============================================================================
 (empty file)
+++ sandbox/math_toolkit/libs/math/example/find_root_example.cpp 20070916 11:27:50 EDT (Sun, 16 Sep 2007)
@@ 0,0 +1,239 @@
+// find_root_example.cpp
+
+// Copyright Paul A. Bristow 2007.
+
+// 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 markup as well as code
+// and comments, don't change any of the special comment markups!
+
+//[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;
+//] //[/root_find1]
+
+namespace boost{ namespace math { namespace tools
+{
+
+template <class F, class T, class Tol>
+inline std::pair<T, T> bracket_and_solve_root(F f, // functor
+ const T& guess,
+ const T& factor,
+ bool rising,
+ Tol tol, // binary functor specifying termination when tol(min, max) becomes true.
+ // eps_tolerance most suitable for this continuous function
+ boost::uintmax_t& max_iter); // explicit (rather than default) max iterations.
+// return interval as a pair containing result.
+
+namespace detail
+{
+
+ // Functor for finding standard deviation:
+ template <class RealType, class Policy>
+ struct standard_deviation_functor
+ {
+ standard_deviation_functor(RealType m, RealType s, RealType d)
+ : mean(m), standard_deviation(s)
+ {
+ }
+ RealType operator()(const RealType& sd)
+ { // Unary functor  the function whose root is to be found.
+ if(sd <= tools::min_value<RealType>())
+ { //
+ return 1;
+ }
+ normal_distribution<RealType, Policy> t(mean, sd);
+ RealType qa = quantile(complement(t, alpha));
+ RealType qb = quantile(complement(t, beta));
+ qa += qb;
+ qa *= qa;
+ qa *= ratio;
+ qa = (df + 1);
+ return qa;
+ } // operator()
+ RealType mean;
+ RealType standard_deviation;
+ }; // struct standard_deviation_functor
+} // namespace detail
+
+template <class RealType, class Policy>
+RealType normal_distribution<RealType, Policy>::estimate_standard_deviation(
+ RealType difference_from_mean,
+ RealType mean,
+ RealType sd,
+ RealType hint) // Best guess available  current sd if none better?
+{
+ static const char* function = "boost::math::normal_distribution<%1%>::estimate_standard_deviation";
+
+ // Check for domain errors:
+ RealType error_result;
+ if(false == detail::check_probability(
+ function, sd, &error_result, Policy())
+ )
+ return error_result;
+
+ if(hint <= 0)
+ { // standard deviation can never be negative.
+ hint = 1;
+ }
+
+ detail::standard_deviation_functor<RealType, Policy> f(mean, sd, difference_from_mean);
+ tools::eps_tolerance<RealType> tol(policies::digits<RealType, Policy>());
+ boost::uintmax_t max_iter = 100;
+ std::pair<RealType, RealType> r = tools::bracket_and_solve_root(f, hint, RealType(2), false, tol, max_iter, Policy());
+ RealType result = r.first + (r.second  r.first) / 2;
+ if(max_iter == 100)
+ {
+ policies::raise_evaluation_error<RealType>(function, "Unable to locate solution in a reasonable time:"
+ " either there is no answer to how many degrees of freedom are required"
+ " or the answer is infinite. Current best guess is %1%", result, Policy());
+ }
+ return result;
+} // estimate_standard_deviation
+} // namespace tools
+} // namespace math
+} // namespace boost
+
+
+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
+[link math_toolkit.toolkit.internals1.roots2 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:
+
+
+
+*/
Added: sandbox/math_toolkit/libs/math/example/find_scale_example.cpp
==============================================================================
 (empty file)
+++ sandbox/math_toolkit/libs/math/example/find_scale_example.cpp 20070916 11:27:50 EDT (Sun, 16 Sep 2007)
@@ 0,0 +1,180 @@
+// find_scale.cpp
+
+// Copyright Paul A. Bristow 2007.
+
+// 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 finding scale (standard deviation) for normal (Gaussian).
+
+// Note that this file contains Quickbook markup as well as code
+// and comments, don't change any of the special comment markups!
+
+//[find_scale1
+/*`
+First we need some includes to access the __normal_distrib,
+the algorithms to find scale (and some std output of course).
+*/
+
+#include <boost/math/distributions/normal.hpp> // for normal_distribution
+ using boost::math::normal; // typedef provides default type is double.
+#include <boost/math/distributions/find_scale.hpp>
+ using boost::math::find_scale;
+ using boost::math::complement; // Needed if you want to use the complement version.
+ using boost::math::policies::policy; // Needed to specify the error handling policy.
+
+#include <iostream>
+ using std::cout; using std::endl;
+#include <iomanip>
+ using std::setw; using std::setprecision;
+#include <limits>
+ using std::numeric_limits;
+//] [/find_scale1]
+
+int main()
+{
+ cout << "Example: Find scale (standard deviation)." << endl;
+ try
+ {
+//[find_scale2
+/*`
+For this example, we will use the standard __normal_distrib,
+with location (mean) zero and standard deviation (scale) unity.
+Conveniently, this is also the default for this implementation's constructor.
+*/
+ normal N01; // Default 'standard' normal distribution with zero mean
+ double sd = 1.; // and standard deviation is 1.
+/*`Suppose we want to find a different normal distribution with standard deviation
+so that only fraction p (here 0.001 or 0.1%) are below a certain chosen limit
+(here 2. standard deviations).
+*/
+ double z = 2.; // z to give prob p
+ double p = 0.001; // only 0.1% below z = 2
+
+ cout << "Normal distribution with mean = " << N01.location() // aka N01.mean()
+ << ", standard deviation " << N01.scale() // aka N01.standard_deviation()
+ << ", has " << "fraction <= " << z
+ << ", p = " << cdf(N01, z) << endl;
+ cout << "Normal distribution with mean = " << N01.location()
+ << ", standard deviation " << N01.scale()
+ << ", has " << "fraction > " << z
+ << ", p = " << cdf(complement(N01, z)) << endl; // Note: uses complement.
+/*`
+[pre
+Normal distribution with mean = 0 has fraction <= 2, p = 0.0227501
+Normal distribution with mean = 0 has fraction > 2, p = 0.97725
+]
+Noting that p = 0.02 instead of our target of 0.001,
+we can now use `find_scale` to give a new standard deviation.
+*/
+ double l = N01.location();
+ double s = find_scale<normal>(z, p, l);
+ cout << "scale (standard deviation) = " << s << endl;
+/*`
+that outputs:
+[pre
+scale (standard deviation) = 0.647201
+]
+showing that we need to reduce the standard deviation from 1. to 0.65.
+
+Then we can check that we have achieved our objective
+by constructing a new distribution
+with the new standard deviation (but same zero mean):
+*/
+ normal np001pc(N01.location(), s);
+/*`
+And recalculating the fraction below (and above) our chosen limit.
+*/
+ cout << "Normal distribution with mean = " << l
+ << " has " << "fraction <= " << z
+ << ", p = " << cdf(np001pc, z) << endl;
+ cout << "Normal distribution with mean = " << l
+ << " has " << "fraction > " << z
+ << ", p = " << cdf(complement(np001pc, z)) << endl;
+/*`
+[pre
+Normal distribution with mean = 0 has fraction <= 2, p = 0.001
+Normal distribution with mean = 0 has fraction > 2, p = 0.999
+]
+
+[h4 Controlling how Errors from find_scale are handled]
+We can also control the policy for handling various errors.
+For example, we can define a new (possibly unwise)
+policy to ignore domain errors ('bad' arguments).
+
+Unless we are using the boost::math namespace, we will need:
+*/
+ using boost::math::policies::policy;
+ using boost::math::policies::domain_error;
+ using boost::math::policies::ignore_error;
+
+/*`
+Using a typedef is convenient, especially if it is reused,
+although it is not required, as the various examples below show.
+*/
+ typedef policy<domain_error<ignore_error> > ignore_domain_policy;
+ // find_scale with new policy, using typedef.
+ l = find_scale<normal>(z, p, l, ignore_domain_policy());
+ // Default policy policy<>, needs using boost::math::policies::policy;
+
+ l = find_scale<normal>(z, p, l, policy<>());
+ // Default policy, fully specified.
+ l = find_scale<normal>(z, p, l, boost::math::policies::policy<>());
+ // New policy, without typedef.
+ l = find_scale<normal>(z, p, l, policy<domain_error<ignore_error> >());
+/*`
+If we want to express a probability, say 0.999, that is a complement, `1  p`
+we should not even think of writing `find_scale<normal>(z, 1  p, l)`,
+but [link why_complements instead], use the __complements version.
+*/
+ z = 2.;
+ double q = 0.999; // = 1  p; // complement of 0.001.
+ sd = find_scale<normal>(complement(z, q, l));
+
+ normal np95pc(l, sd); // Same standard_deviation (scale) but with mean(scale) shifted
+ cout << "Normal distribution with mean = " << l << " has "
+ << "fraction <= " << z << " = " << cdf(np95pc, z) << endl;
+ cout << "Normal distribution with mean = " << l << " has "
+ << "fraction > " << z << " = " << cdf(complement(np95pc, z)) << endl;
+
+/*`
+Sadly, it is all too easy to get probabilities the wrong way round,
+when you may get a warning like this:
+[pre
+Message from thrown exception was:
+ Error in function boost::math::find_scale<Dist, Policy>(complement(double, double, double, Policy)):
+ Computed scale (0.48043523852179076) is <= 0! Was the complement intended?
+]
+The default error handling policy is to throw an exception with this message,
+but if you chose a policy to ignore the error,
+the (impossible) negative scale is quietly returned.
+*/
+//] [/find_scale2]
+ }
+ 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()
+
+//[find_scale_example_output
+/*`
+[pre
+Example: Find scale (standard deviation).
+Normal distribution with mean = 0, standard deviation 1, has fraction <= 2, p = 0.0227501
+Normal distribution with mean = 0, standard deviation 1, has fraction > 2, p = 0.97725
+scale (standard deviation) = 0.647201
+Normal distribution with mean = 0 has fraction <= 2, p = 0.001
+Normal distribution with mean = 0 has fraction > 2, p = 0.999
+Normal distribution with mean = 0.946339 has fraction <= 2 = 0.001
+Normal distribution with mean = 0.946339 has fraction > 2 = 0.999
+]
+*/
+//] [/find_scale_example_output]
BoostCommit list run by bdawes at acm.org, david.abrahams at rcn.com, gregod at cs.rpi.edu, cpdaniel at pacbell.net, john at johnmaddock.co.uk