[Boost-commit] svn:boost r39325 - sandbox/math_toolkit/libs/math/example

From: pbristow_at_[hidden]
Date: 2007-09-16 11:27:51

Author: pbristow
Date: 2007-09-16 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	2007-09-16 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 mark-up as well as code
+// and comments, don't change any of the special comment mark-ups!
+
+//[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 re-calculating 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 re-used,
+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	2007-09-16 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 mark-up as well as code
+// and comments, don't change any of the special comment mark-ups!
+
+//[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 1-58488-635-8, (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 Cauchy-Lorentz 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  Cauchy-Lorentz 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 over-fill to make sure we meet a minimum specification (or under-fill 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, X-ray 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	2007-09-16 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 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;
+//] //[/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	2007-09-16 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 mark-up as well as code
+// and comments, don't change any of the special comment mark-ups!
+
+//[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 re-calculating 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 re-used,
+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]