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

   1 // distribution_construction.cpp
   2
   3 // Copyright Paul A. Bristow 2007, 2010, 2012.
   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 // Caution: this file contains Quickbook markup as well as code
  11 // and comments, don't change any of the special comment markups!
  12
  13 #ifdef _MSC_VER
  14 #  pragma warning (disable : 4996) // disable -D_SCL_SECURE_NO_WARNINGS C++ 'Checked Iterators'
  15 #endif
  16
  17 #include <iostream>
  18 #include <exception>
  19
  20 //[distribution_construction_1
  21
  22 /*`
  23 The structure of distributions is rather different from some other statistical libraries,
  24 for example, those written in less object-oriented language like FORTRAN and C:
  25 these provide a few arguments to each free function.
  26
  27 Boost.Math library provides each distribution as a template C++ class.
  28 A distribution is constructed with a few arguments, and then
  29 member and non-member functions are used to find values of the
  30 distribution, often a function of a random variate.
  31
  32 For this demonstration, first we need some includes to access the
  33 negative binomial distribution (and the binomial, beta and gamma distributions too).
  34
  35 To demonstrate the use with a high precision User-defined floating-point type
  36 `cpp_dec_float` we also need an include from Boost.Multiprecision.
  37 */
  38
  39 #include <boost/math/distributions/negative_binomial.hpp> // for negative_binomial_distribution
  40   using boost::math::negative_binomial_distribution; // default type is double.
  41   using boost::math::negative_binomial; // typedef provides default type is double.
  42 #include <boost/math/distributions/binomial.hpp> // for binomial_distribution.
  43 #include <boost/math/distributions/beta.hpp> // for beta_distribution.
  44 #include <boost/math/distributions/gamma.hpp> // for gamma_distribution.
  45 #include <boost/math/distributions/normal.hpp> // for normal_distribution.
  46
  47 #include <boost/multiprecision/cpp_dec_float.hpp> // for cpp_dec_float_100
  48 /*`
  49 Several examples of constructing distributions follow:
  50 */
  51 //] [/distribution_construction_1 end of Quickbook in C++ markup]
  52
  53 int main()
  54 {
  55   try
  56   {
  57 //[distribution_construction_2
  58 /*`
  59 First, a negative binomial distribution with 8 successes
  60 and a success fraction 0.25, 25% or 1 in 4, is constructed like this:
  61 */
  62   boost::math::negative_binomial_distribution<double> mydist0(8., 0.25);
  63   /*`
  64   But this is inconveniently long, so we might be tempted to write
  65   */
  66   using namespace boost::math;
  67   /*`
  68   but this might risk ambiguity with names in `std random` so
  69   [*much] better is explicit `using boost::math::` statements, for example:
  70   */
  71   using boost::math::negative_binomial_distribution;
  72   /*`
  73   and we can still reduce typing.
  74
  75   Since the vast majority of applications use will be using `double` precision,
  76   the template argument to the distribution (`RealType`) defaults
  77   to type `double`, so we can also write:
  78   */
  79
  80   negative_binomial_distribution<> mydist9(8., 0.25); // Uses default `RealType = double`.
  81
  82   /*`
  83   But the name `negative_binomial_distribution` is still inconveniently long,
  84   so, for most distributions, a convenience `typedef` is provided, for example:
  85
  86      typedef negative_binomial_distribution<double> negative_binomial; // Reserved name of type double.
  87
  88   [caution
  89   This convenience typedef is [*not provided] if a clash would occur
  90   with the name of a function: currently only `beta` and `gamma`
  91   fall into this category.
  92   ]
  93
  94   So, after a using statement,
  95   */
  96
  97   using boost::math::negative_binomial;
  98
  99   /*`
 100   we have a convenient typedef to `negative_binomial_distribution<double>`:
 101   */
 102   negative_binomial mydist(8., 0.25);
 103
 104   /*`
 105   Some more examples using the convenience typedef:
 106   */
 107   negative_binomial mydist10(5., 0.4); // Both arguments double.
 108   /*`
 109   And automatic conversion takes place, so you can use integers and floats:
 110   */
 111   negative_binomial mydist11(5, 0.4); // Using provided typedef double, int and double arguments.
 112   /*`
 113   This is probably the most common usage.
 114   */
 115   negative_binomial mydist12(5., 0.4F); // Double and float arguments.
 116   negative_binomial mydist13(5, 1); // Both arguments integer.
 117
 118   /*`
 119   Similarly for most other distributions like the binomial.
 120   */
 121   binomial mybinomial(1, 0.5); // is more concise than
 122   binomial_distribution<> mybinomd1(1, 0.5);
 123
 124   /*`
 125   For cases when the typdef distribution name would clash with a math special function
 126   (currently only beta and gamma)
 127   the typedef is deliberately not provided, and the longer version of the name
 128   must be used.  For example do not use:
 129
 130      using boost::math::beta;
 131      beta mybetad0(1, 0.5); // Error beta is a math FUNCTION!
 132
 133   Which produces the error messages:
 134
 135   [pre
 136   error C2146: syntax error : missing ';' before identifier 'mybetad0'
 137   warning C4551: function call missing argument list
 138   error C3861: 'mybetad0': identifier not found
 139   ]
 140
 141   Instead you should use:
 142   */
 143   using boost::math::beta_distribution;
 144   beta_distribution<> mybetad1(1, 0.5);
 145   /*`
 146   or for the gamma distribution:
 147   */
 148   gamma_distribution<> mygammad1(1, 0.5);
 149
 150   /*`
 151   We can, of course, still provide the type explicitly thus:
 152   */
 153
 154   // Explicit double precision:  both arguments are double:
 155   negative_binomial_distribution<double>        mydist1(8., 0.25);
 156
 157   // Explicit float precision, double arguments are truncated to float:
 158   negative_binomial_distribution<float>         mydist2(8., 0.25);
 159
 160   // Explicit float precision, integer & double arguments converted to float:
 161   negative_binomial_distribution<float>         mydist3(8, 0.25);
 162
 163   // Explicit float precision, float arguments, so no conversion:
 164   negative_binomial_distribution<float>         mydist4(8.F, 0.25F);
 165
 166   // Explicit float precision, integer arguments promoted to float.
 167   negative_binomial_distribution<float>         mydist5(8, 1);
 168
 169   // Explicit double precision:
 170   negative_binomial_distribution<double>        mydist6(8., 0.25);
 171
 172   // Explicit long double precision:
 173   negative_binomial_distribution<long double>   mydist7(8., 0.25);
 174
 175   /*`
 176   And you can use your own RealType,
 177   for example, `boost::math::cpp_dec_float_50` (an arbitrary 50 decimal digits precision type),
 178   then we can write:
 179   */
 180   using namespace boost::multiprecision;
 181
 182   negative_binomial_distribution<cpp_dec_float_50>  mydist8(8, 0.25);
 183   // `integer` arguments are promoted to your RealType exactly, but
 184   // `double` argument are converted to RealType,
 185   // possibly losing precision, so don't write:
 186
 187   negative_binomial_distribution<cpp_dec_float_50>  mydist20(8, 0.23456789012345678901234567890);
 188  // to avoid truncation of second parameter to `0.2345678901234567`.
 189
 190   negative_binomial_distribution<cpp_dec_float_50>  mydist21(8, cpp_dec_float_50("0.23456789012345678901234567890") );
 191
 192   // Ensure that all potentially significant digits are shown.
 193   std::cout.precision(std::numeric_limits<cpp_dec_float_50>::digits10);
 194   cpp_dec_float_50 x("1.23456789012345678901234567890");
 195   std::cout << pdf(mydist8, x) << std::endl;
 196 /*` showing  0.00012630010495970320103876754721976419438231705359935
 197
 198 [warning When using multiprecision, it is all too easy to get accidental truncation!]
 199
 200 For example, if you write
 201 */
 202   std::cout << pdf(mydist8, 1.23456789012345678901234567890) << std::endl;
 203 /*`
 204 showing  0.00012630010495970318465064569310967179576805651692929,
 205 which is wrong at about the 17th decimal digit!
 206
 207 This is because the value provided is truncated to a `double`, effectively
 208   `double x = 1.23456789012345678901234567890;`
 209
 210 Then the now `double x` is passed to function `pdf`,
 211 and this truncated `double` value is finally promoted to `cpp_dec_float_50`.
 212
 213 Another way of quietly getting the wrong answer is to write:
 214 */
 215   std::cout << pdf(mydist8, cpp_dec_float_50(1.23456789012345678901234567890)) << std::endl;
 216 /*`
 217 A correct way from a multi-digit string value is
 218 */
 219   std::cout << pdf(mydist8, cpp_dec_float_50("1.23456789012345678901234567890")) << std::endl;
 220 /*`
 221
 222 [tip Getting about 17 decimal digits followed by many zeros is often a sign of accidental truncation.]
 223 */
 224
 225 /*`
 226 [h4 Default arguments to distribution constructors.]
 227
 228 Note that default constructor arguments are only provided for some distributions.
 229 So if you wrongly assume a default argument, you will get an error message, for example:
 230
 231    negative_binomial_distribution<> mydist8;
 232
 233 [pre error C2512 no appropriate default constructor available.]
 234
 235 No default constructors are provided for the `negative binomial` distribution,
 236 because it is difficult to chose any sensible default values for this distribution.
 237
 238 For other distributions, like the normal distribution,
 239 it is obviously very useful to provide 'standard'
 240 defaults for the mean (zero) and standard deviation (unity) thus:
 241
 242     normal_distribution(RealType mean = 0, RealType sd = 1);
 243
 244 So in this case we can write:
 245 */
 246   using boost::math::normal;
 247
 248   normal norm1;       // Standard normal distribution.
 249   normal norm2(2);    // Mean = 2, std deviation = 1.
 250   normal norm3(2, 3); // Mean = 2, std deviation = 3.
 251
 252   }
 253   catch(std::exception &ex)
 254   {
 255     std::cout << ex.what() << std::endl;
 256   }
 257
 258   return 0;
 259 }  // int main()
 260
 261 /*`There is no useful output from this demonstration program, of course. */
 262
 263 //] [/end of distribution_construction_2]
 264
 265 /*
 266 //[distribution_construction_output
 267
 268   0.00012630010495970320103876754721976419438231705359935
 269   0.00012630010495970318465064569310967179576805651692929
 270   0.00012630010495970318465064569310967179576805651692929
 271   0.00012630010495970320103876754721976419438231705359935
 272
 273 //] [/distribution_construction_output]
 274
 275 */
 276
 277
 278