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