ceph/src/boost/libs/math/example/find_scale_example.cpp

   1 // find_scale.cpp
   2
   3 // Copyright Paul A. Bristow 2007, 2010.
   4
   5 // Use, modification and distribution are subject to the
   6 // Boost Software License, Version 1.0.
   7 // (See accompanying file LICENSE_1_0.txt
   8 // or copy at http://www.boost.org/LICENSE_1_0.txt)
   9
  10 // Example of finding scale (standard deviation) for normal (Gaussian).
  11
  12 // Note that this file contains Quickbook mark-up as well as code
  13 // and comments, don't change any of the special comment mark-ups!
  14
  15 //[find_scale1
  16 /*`
  17 First we need some includes to access the __normal_distrib,
  18 the algorithms to find scale (and some std output of course).
  19 */
  20
  21 #include <boost/math/distributions/normal.hpp> // for normal_distribution
  22   using boost::math::normal; // typedef provides default type is double.
  23 #include <boost/math/distributions/find_scale.hpp>
  24   using boost::math::find_scale;
  25   using boost::math::complement; // Needed if you want to use the complement version.
  26   using boost::math::policies::policy; // Needed to specify the error handling policy.
  27
  28 #include <iostream>
  29   using std::cout; using std::endl;
  30 #include <iomanip>
  31   using std::setw; using std::setprecision;
  32 #include <limits>
  33   using std::numeric_limits;
  34 //] [/find_scale1]
  35
  36 int main()
  37 {
  38   cout << "Example: Find scale (standard deviation)." << endl;
  39   try
  40   {
  41 //[find_scale2
  42 /*`
  43 For this example, we will use the standard __normal_distrib,
  44 with location (mean) zero and standard deviation (scale) unity.
  45 Conveniently, this is also the default for this implementation's constructor.
  46 */
  47   normal N01;  // Default 'standard' normal distribution with zero mean
  48   double sd = 1.; // and standard deviation is 1.
  49 /*`Suppose we want to find a different normal distribution with standard deviation
  50 so that only fraction p (here 0.001 or 0.1%) are below a certain chosen limit
  51 (here -2. standard deviations).
  52 */
  53   double z = -2.; // z to give prob p
  54   double p = 0.001; // only 0.1% below z = -2
  55
  56   cout << "Normal distribution with mean = " << N01.location()  // aka N01.mean()
  57     << ", standard deviation " << N01.scale() // aka N01.standard_deviation()
  58     << ", has " << "fraction <= " << z
  59     << ", p = "  << cdf(N01, z) << endl;
  60   cout << "Normal distribution with mean = " << N01.location()
  61     << ", standard deviation " << N01.scale()
  62     << ", has " << "fraction > " << z
  63     << ", p = "  << cdf(complement(N01, z)) << endl; // Note: uses complement.
  64 /*`
  65 [pre
  66 Normal distribution with mean = 0 has fraction <= -2, p = 0.0227501
  67 Normal distribution with mean = 0 has fraction > -2, p = 0.97725
  68 ]
  69 Noting that p = 0.02 instead of our target of 0.001,
  70 we can now use `find_scale` to give a new standard deviation.
  71 */
  72    double l = N01.location();
  73    double s = find_scale<normal>(z, p, l);
  74    cout << "scale (standard deviation) = " << s << endl;
  75 /*`
  76 that outputs:
  77 [pre
  78 scale (standard deviation) = 0.647201
  79 ]
  80 showing that we need to reduce the standard deviation from 1. to 0.65.
  81
  82 Then we can check that we have achieved our objective
  83 by constructing a new distribution
  84 with the new standard deviation (but same zero mean):
  85 */
  86   normal np001pc(N01.location(), s);
  87 /*`
  88 And re-calculating the fraction below (and above) our chosen limit.
  89 */
  90   cout << "Normal distribution with mean = " << l
  91     << " has " << "fraction <= " << z
  92     << ", p = "  << cdf(np001pc, z) << endl;
  93   cout << "Normal distribution with mean = " << l
  94     << " has " << "fraction > " << z
  95     << ", p = "  << cdf(complement(np001pc, z)) << endl;
  96 /*`
  97 [pre
  98 Normal distribution with mean = 0 has fraction <= -2, p = 0.001
  99 Normal distribution with mean = 0 has fraction > -2, p = 0.999
 100 ]
 101
 102 [h4 Controlling how Errors from find_scale are handled]
 103 We can also control the policy for handling various errors.
 104 For example, we can define a new (possibly unwise)
 105 policy to ignore domain errors ('bad' arguments).
 106
 107 Unless we are using the boost::math namespace, we will need:
 108 */
 109   using boost::math::policies::policy;
 110   using boost::math::policies::domain_error;
 111   using boost::math::policies::ignore_error;
 112
 113 /*`
 114 Using a typedef is convenient, especially if it is re-used,
 115 although it is not required, as the various examples below show.
 116 */
 117   typedef policy<domain_error<ignore_error> > ignore_domain_policy;
 118   // find_scale with new policy, using typedef.
 119   l = find_scale<normal>(z, p, l, ignore_domain_policy());
 120   // Default policy policy<>, needs using boost::math::policies::policy;
 121
 122   l = find_scale<normal>(z, p, l, policy<>());
 123   // Default policy, fully specified.
 124   l = find_scale<normal>(z, p, l, boost::math::policies::policy<>());
 125   // New policy, without typedef.
 126   l = find_scale<normal>(z, p, l, policy<domain_error<ignore_error> >());
 127 /*`
 128 If we want to express a probability, say 0.999, that is a complement, `1 - p`
 129 we should not even think of writing `find_scale<normal>(z, 1 - p, l)`,
 130 but use the __complements version (see __why_complements).
 131 */
 132   z = -2.;
 133   double q = 0.999; // = 1 - p; // complement of 0.001.
 134   sd = find_scale<normal>(complement(z, q, l));
 135
 136   normal np95pc(l, sd); // Same standard_deviation (scale) but with mean(scale) shifted
 137   cout << "Normal distribution with mean = " << l << " has "
 138     << "fraction <= " << z << " = "  << cdf(np95pc, z) << endl;
 139   cout << "Normal distribution with mean = " << l << " has "
 140     << "fraction > " << z << " = "  << cdf(complement(np95pc, z)) << endl;
 141
 142 /*`
 143 Sadly, it is all too easy to get probabilities the wrong way round,
 144 when you may get a warning like this:
 145 [pre
 146 Message from thrown exception was:
 147    Error in function boost::math::find_scale<Dist, Policy>(complement(double, double, double, Policy)):
 148    Computed scale (-0.48043523852179076) is <= 0! Was the complement intended?
 149 ]
 150 The default error handling policy is to throw an exception with this message,
 151 but if you chose a policy to ignore the error,
 152 the (impossible) negative scale is quietly returned.
 153 */
 154 //] [/find_scale2]
 155   }
 156   catch(const std::exception& e)
 157   { // Always useful to include try & catch blocks because default policies
 158     // are to throw exceptions on arguments that cause errors like underflow, overflow.
 159     // Lacking try & catch blocks, the program will abort without a message below,
 160     // which may give some helpful clues as to the cause of the exception.
 161     std::cout <<
 162       "\n""Message from thrown exception was:\n   " << e.what() << std::endl;
 163   }
 164   return 0;
 165 }  // int main()
 166
 167 //[find_scale_example_output
 168 /*`
 169 [pre
 170 Example: Find scale (standard deviation).
 171 Normal distribution with mean = 0, standard deviation 1, has fraction <= -2, p = 0.0227501
 172 Normal distribution with mean = 0, standard deviation 1, has fraction > -2, p = 0.97725
 173 scale (standard deviation) = 0.647201
 174 Normal distribution with mean = 0 has fraction <= -2, p = 0.001
 175 Normal distribution with mean = 0 has fraction > -2, p = 0.999
 176 Normal distribution with mean = 0.946339 has fraction <= -2 = 0.001
 177 Normal distribution with mean = 0.946339 has fraction > -2 = 0.999
 178 ]
 179 */
 180 //] [/find_scale_example_output]