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

   1 // wald_example.cpp or inverse_gaussian_example.cpp
   2
   3 // Copyright Paul A. Bristow 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 using the Inverse Gaussian (or Inverse Normal) distribution.
  11 // The Wald Distribution is
  12
  13
  14 // Note that this file contains Quickbook mark-up as well as code
  15 // and comments, don't change any of the special comment mark-ups!
  16
  17 //[inverse_gaussian_basic1
  18 /*`
  19 First we need some includes to access the normal distribution
  20 (and some std output of course).
  21 */
  22
  23 #ifdef _MSC_VER
  24 # pragma warning (disable : 4224)
  25 # pragma warning (disable : 4189)
  26 # pragma warning (disable : 4100)
  27 # pragma warning (disable : 4224)
  28 # pragma warning (disable : 4512)
  29 # pragma warning (disable : 4702)
  30 # pragma warning (disable : 4127)
  31 #endif
  32
  33 //#define BOOST_MATH_INSTRUMENT
  34
  35 #define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error
  36 #define BOOST_MATH_DOMAIN_ERROR_POLICY ignore_error
  37
  38 #include <boost/math/distributions/inverse_gaussian.hpp> // for inverse_gaussian_distribution
  39   using boost::math::inverse_gaussian; // typedef provides default type is double.
  40   using boost::math::inverse_gaussian_distribution; // for inverse gaussian distribution.
  41
  42 #include <boost/math/distributions/normal.hpp> // for normal_distribution
  43 using boost::math::normal; // typedef provides default type is double.
  44
  45 #include <boost/array.hpp>
  46 using boost::array;
  47
  48 #include <iostream>
  49   using std::cout; using std::endl; using std::left; using std::showpoint; using std::noshowpoint;
  50 #include <iomanip>
  51   using std::setw; using std::setprecision;
  52 #include <limits>
  53   using std::numeric_limits;
  54 #include <sstream>
  55   using std::string;
  56 #include <string>
  57   using std::stringstream;
  58
  59 // const double tol = 3 * numeric_limits<double>::epsilon();
  60
  61 int main()
  62 {
  63   cout << "Example: Inverse Gaussian Distribution."<< endl;
  64
  65  try
  66   {
  67
  68       double tolfeweps = numeric_limits<double>::epsilon();
  69       //cout << "Tolerance = " << tol << endl;
  70
  71       int precision = 17; // traditional tables are only computed to much lower precision.
  72       cout.precision(17); // std::numeric_limits<double>::max_digits10; for 64-bit doubles.
  73
  74      // Traditional tables and values.
  75      double step = 0.2; // in z
  76       double range = 4; // min and max z = -range to +range.
  77       // Construct a (standard) inverse gaussian distribution s
  78       inverse_gaussian w11(1, 1);
  79       // (default mean = units, and standard deviation = unity)
  80       cout << "(Standard) Inverse Gaussian distribution, mean = "<< w11.mean()
  81           << ", scale = " << w11.scale() << endl;
  82
  83 /*` First the probability distribution function (pdf).
  84  */
  85       cout << "Probability distribution function (pdf) values" << endl;
  86       cout << "  z " "      pdf " << endl;
  87       cout.precision(5);
  88       for (double z = (numeric_limits<double>::min)(); z < range + step; z += step)
  89       {
  90         cout << left << setprecision(3) << setw(6) << z << " "
  91           << setprecision(precision) << setw(12) << pdf(w11, z) << endl;
  92       }
  93       cout.precision(6); // default
  94  /*`And the area under the normal curve from -[infin] up to z,
  95       the cumulative distribution function (cdf).
  96 */
  97
  98       // For a (default) inverse gaussian distribution.
  99       cout << "Integral (area under the curve) from 0 up to z (cdf) " << endl;
 100       cout << "  z " "      cdf " << endl;
 101       for (double z = (numeric_limits<double>::min)(); z < range + step; z += step)
 102       {
 103         cout << left << setprecision(3) << setw(6) << z << " "
 104           << setprecision(precision) << setw(12) << cdf(w11, z) << endl;
 105       }
 106       /*`giving the following table:
 107 [pre
 108     z       pdf
 109   2.23e-308 -1.#IND
 110   0.2    0.90052111680384117
 111   0.4    1.0055127039453111
 112   0.6    0.75123750098955733
 113   0.8    0.54377310461643302
 114   1      0.3989422804014327
 115   1.2    0.29846949816803292
 116   1.4    0.2274579835638664
 117   1.6    0.17614566625628389
 118   1.8    0.13829083543591469
 119   2      0.10984782236693062
 120   2.2    0.088133964251182237
 121   2.4    0.071327382959107177
 122   2.6    0.058162562161661699
 123   2.8    0.047742223328567722
 124   3      0.039418357969819712
 125   3.2    0.032715223861241892
 126   3.4    0.027278388940958308
 127   3.6    0.022840312999395804
 128   3.8    0.019196657941016954
 129   4      0.016189699458236451
 130   Integral (area under the curve) from 0 up to z (cdf)
 131     z       cdf
 132   2.23e-308 0
 133   0.2    0.063753567519976254
 134   0.4    0.2706136704424541
 135   0.6    0.44638391340412931
 136   0.8    0.57472390962590925
 137   1      0.66810200122317065
 138   1.2    0.73724578422952536
 139   1.4    0.78944214237790356
 140   1.6    0.82953458108474554
 141   1.8    0.86079282968276671
 142   2      0.88547542598600626
 143   2.2    0.90517870624273966
 144   2.4    0.92105495653509362
 145   2.6    0.93395164268166786
 146   2.8    0.94450240360053817
 147   3      0.95318792074278835
 148   3.2    0.96037753019309191
 149   3.4    0.96635823989417369
 150   3.6    0.97135533107998406
 151   3.8    0.97554722413538364
 152   4      0.97907636417888622
 153 ]
 154
 155 /*`We can get the inverse, the quantile, percentile, percentage point, or critical value
 156 for a probability for a few probability from the above table, for z = 0.4, 1.0, 2.0:
 157 */
 158       cout << quantile(w11, 0.27061367044245421 ) << endl; // 0.4
 159       cout << quantile(w11, 0.66810200122317065) << endl; // 1.0
 160       cout << quantile(w11, 0.88547542598600615) << endl; // 2.0
 161 /*`turning the expect values apart from some 'computational noise' in the least significant bit or two.
 162
 163 [pre
 164   0.40000000000000008
 165   0.99999999999999967
 166   1.9999999999999973
 167 ]
 168
 169 */
 170
 171     //  cout << "pnorm01(-0.406053) " << pnorm01(-0.406053) << ", cdfn01(-0.406053) = " << cdf(n01, -0.406053) << endl;
 172    //cout << "pnorm01(0.5) = " << pnorm01(0.5) << endl; // R pnorm(0.5,0,1) = 0.6914625  == 0.69146246127401312
 173     // R qnorm(0.6914625,0,1) = 0.5
 174
 175     // formatC(SuppDists::qinvGauss(0.3649755481729598, 1, 1), digits=17)  [1] "0.50000000969034875"
 176
 177
 178
 179   // formatC(SuppDists::dinvGauss(0.01, 1, 1), digits=17) [1] "2.0811768202028392e-19"
 180   // formatC(SuppDists::pinvGauss(0.01, 1, 1), digits=17) [1] "4.122313403318778e-23"
 181
 182
 183
 184   //cout << " qinvgauss(0.3649755481729598, 1, 1) = " << qinvgauss(0.3649755481729598, 1, 1) << endl;  // 0.5
 185  // cout << quantile(s, 0.66810200122317065) << endl; // expect 1, get 0.50517388467190727
 186   //cout << " qinvgauss(0.62502320258649202, 1, 1) = " << qinvgauss(0.62502320258649202, 1, 1) << endl; // 0.9
 187   //cout << " qinvgauss(0.063753567519976254, 1, 1) = " << qinvgauss(0.063753567519976254, 1, 1) << endl; // 0.2
 188   //cout << " qinvgauss(0.0040761113207110162, 1, 1) = " << qinvgauss(0.0040761113207110162, 1, 1) << endl; // 0.1
 189
 190   //double x = 1.; // SuppDists::pinvGauss(0.4, 1,1) [1] 0.2706137
 191   //double c = pinvgauss(x, 1, 1); // 0.3649755481729598 ==   cdf(x, 1,1) 0.36497554817295974
 192   //cout << "  pinvgauss(x, 1, 1) = " << c << endl; //  pinvgauss(x, 1, 1) = 0.27061367044245421
 193   //double p = pdf(w11, x);
 194   //double c = cdf(w11, x); // cdf(1, 1, 1) = 0.66810200122317065
 195   //cout << "cdf(" << x << ", " << w11.mean() << ", "<< w11.scale() << ") = " << c << endl; // cdf(x, 1, 1) 0.27061367044245421
 196   //cout << "pdf(" << x << ", " << w11.mean() << ", "<< w11.scale() << ") = " << p << endl;
 197   //double q = quantile(w11, c);
 198   //cout << "quantile(w11, " << c <<  ") = " << q << endl;
 199
 200   //cout  << "quantile(w11, 4.122313403318778e-23) = "<< quantile(w11, 4.122313403318778e-23) << endl; // quantile
 201   //cout << "quantile(w11, 4.8791443010851493e-219) = " << quantile(w11, 4.8791443010851493e-219) << endl; // quantile
 202
 203   //double c1 = 1 - cdf(w11, x); //  1 - cdf(1, 1, 1) = 0.33189799877682935
 204   //cout << "1 - cdf(" << x << ", " << w11.mean() << ", " << w11.scale() << ") = " << c1 << endl; // cdf(x, 1, 1) 0.27061367044245421
 205   //double cc = cdf(complement(w11, x));
 206   //cout << "cdf(complement(" << x << ", " << w11.mean() << ", "<< w11.scale() << ")) = " << cc << endl; // cdf(x, 1, 1) 0.27061367044245421
 207   //// 1 - cdf(1000, 1, 1) = 0
 208   //// cdf(complement(1000, 1, 1)) = 4.8694344366900402e-222
 209
 210   //cout << "quantile(w11, " << c << ") = "<< quantile(w11, c) << endl; // quantile = 0.99999999999999978 == x = 1
 211   //cout << "quantile(w11, " << c << ") = "<< quantile(w11, 1 - c) << endl; // quantile complement. quantile(w11, 0.66810200122317065) = 0.46336593652340152
 212 //  cout << "quantile(complement(w11, " << c << ")) = " << quantile(complement(w11, c)) << endl; // quantile complement                = 0.46336593652340163
 213
 214   // cdf(1, 1, 1) = 0.66810200122317065
 215   // 1 - cdf(1, 1, 1) = 0.33189799877682935
 216   // cdf(complement(1, 1, 1)) = 0.33189799877682929
 217
 218   // quantile(w11, 0.66810200122317065) = 0.99999999999999978
 219   // 1 - quantile(w11, 0.66810200122317065) = 2.2204460492503131e-016
 220   // quantile(complement(w11, 0.33189799877682929)) = 0.99999999999999989
 221
 222
 223   // qinvgauss(c, 1, 1) = 0.3999999999999998
 224   // SuppDists::qinvGauss(0.270613670442454, 1, 1) [1] 0.4
 225
 226
 227   /*
 228   double qs = pinvgaussU(c, 1, 1); //
 229     cout << "qinvgaussU(c, 1, 1) = " << qs << endl; // qinvgaussU(c, 1, 1) = 0.86567442459240929
 230     // > z=q - exp(c) * p [1] 0.8656744 qs 0.86567442459240929 double
 231     // Is this the complement?
 232     cout << "qgamma(0.2, 0.5, 1) expect 0.0320923 = " << qgamma(0.2, 0.5, 1) << endl;
 233     // qgamma(0.2, 0.5, 1) expect 0.0320923 = 0.032092377333650807
 234
 235
 236   cout << "qinvgauss(pinvgauss(x, 1, 1) = " << q
 237   << ", diff = " << x - q << ", fraction = " << (x - q) /x << endl; // 0.5
 238
 239  */   // > SuppDists::pinvGauss(0.02, 1,1)  [1] 4.139176e-12
 240   // > SuppDists::qinvGauss(4.139176e-12, 1,1) [1] 0.02000000
 241
 242
 243     // pinvGauss(1,1,1) = 0.668102  C++  == 0.66810200122317065
 244   // qinvGauss(0.668102,1,1) = 1
 245
 246    //  SuppDists::pinvGauss(0.3,1,1) = 0.1657266
 247   // cout << "qinvGauss(0.0040761113207110162, 1, 1) = " << qinvgauss(0.0040761113207110162, 1, 1) << endl;
 248   //cout << "quantile(s, 0.1657266) = " << quantile(s, 0.1657266) << endl; // expect 1.
 249
 250   //wald s12(2, 1);
 251   //cout << "qinvGauss(0.3, 2, 1) = " << qinvgauss(0.3, 2, 1) << endl; // SuppDists::qinvGauss(0.3,2,1) == 0.58288065635052944
 252   //// but actually get qinvGauss(0.3, 2, 1) = 0.58288064777632187
 253   //cout  << "cdf(s12, 0.3) = " << cdf(s12, 0.3) << endl; //  cdf(s12, 0.3) = 0.10895339868447573
 254
 255  // using boost::math::wald;
 256   //cout.precision(6);
 257
 258  /*
 259  double m = 1;
 260   double l = 1;
 261   double x = 0.1;
 262   //c = cdf(w, x);
 263   double p = pinvgauss(x, m, l);
 264   cout << "x = " << x << ",  pinvgauss(x, m, l) = " << p << endl; // R 0.4 0.2706137
 265   double qg = qgamma(1.- p, 0.5, 1.0, true, false);
 266   cout << " qgamma(1.- p, 0.5, 1.0, true, false) = " << qg << endl; // 0.606817
 267   double g = guess_whitmore(p, m, l);
 268   cout << "m = " << m << ", l = " << l << ",   x = " << x << ", guess = " << g
 269     << ", diff = " << (x - g) << endl;
 270
 271   g = guess_wheeler(p, m, l);
 272    cout << "m = " << m << ", l = " << l << ",   x = " << x << ", guess = " << g
 273     << ", diff = " << (x - g) << endl;
 274
 275    g = guess_bagshaw(p, m, l);
 276    cout << "m = " << m << ", l = " << l << ",   x = " << x << ", guess = " << g
 277     << ", diff = " << (x - g) << endl;
 278
 279    // m = 1, l = 10,   x = 0.9, guess = 0.89792, diff = 0.00231075 so a better fit.
 280   // x = 0.9, guess = 0.887907
 281   // x = 0.5, guess = 0.474977
 282   // x = 0.4, guess = 0.369597
 283   // x = 0.2, guess = 0.155196
 284
 285   // m = 1, l = 2,   x = 0.9, guess = 1.0312, diff = -0.145778
 286   // m = 1, l = 2,   x = 0.1, guess = 0.122201, diff = -0.222013
 287   //  m = 1, l = 2,   x = 0.2, guess = 0.299326, diff = -0.49663
 288   //   m = 1, l = 2,   x = 0.5, guess = 1.00437, diff = -1.00875
 289   // m = 1, l = 2,   x = 0.7, guess = 1.01517, diff = -0.450247
 290
 291   double ls[7] = {0.1, 0.2, 0.5, 1., 2., 10, 100}; // scale values.
 292   double ms[10] = {0.001, 0.02, 0.1, 0.2, 0.5, 0.9, 1., 2., 10, 100};  // mean values.
 293    */
 294
 295     cout.precision(6); // Restore to default.
 296   } // try
 297   catch(const std::exception& e)
 298   { // Always useful to include try & catch blocks because default policies
 299     // are to throw exceptions on arguments that cause errors like underflow, overflow.
 300     // Lacking try & catch blocks, the program will abort without a message below,
 301     // which may give some helpful clues as to the cause of the exception.
 302     std::cout <<
 303       "\n""Message from thrown exception was:\n   " << e.what() << std::endl;
 304   }
 305   return 0;
 306 }  // int main()
 307
 308
 309 /*
 310
 311 Output is:
 312
 313 inverse_gaussian_example.cpp
 314   inverse_gaussian_example.vcxproj -> J:\Cpp\MathToolkit\test\Math_test\Debug\inverse_gaussian_example.exe
 315   Example: Inverse Gaussian Distribution.
 316   (Standard) Inverse Gaussian distribution, mean = 1, scale = 1
 317   Probability distribution function (pdf) values
 318     z       pdf
 319   2.23e-308 -1.#IND
 320   0.2    0.90052111680384117
 321   0.4    1.0055127039453111
 322   0.6    0.75123750098955733
 323   0.8    0.54377310461643302
 324   1      0.3989422804014327
 325   1.2    0.29846949816803292
 326   1.4    0.2274579835638664
 327   1.6    0.17614566625628389
 328   1.8    0.13829083543591469
 329   2      0.10984782236693062
 330   2.2    0.088133964251182237
 331   2.4    0.071327382959107177
 332   2.6    0.058162562161661699
 333   2.8    0.047742223328567722
 334   3      0.039418357969819712
 335   3.2    0.032715223861241892
 336   3.4    0.027278388940958308
 337   3.6    0.022840312999395804
 338   3.8    0.019196657941016954
 339   4      0.016189699458236451
 340   Integral (area under the curve) from 0 up to z (cdf)
 341     z       cdf
 342   2.23e-308 0
 343   0.2    0.063753567519976254
 344   0.4    0.2706136704424541
 345   0.6    0.44638391340412931
 346   0.8    0.57472390962590925
 347   1      0.66810200122317065
 348   1.2    0.73724578422952536
 349   1.4    0.78944214237790356
 350   1.6    0.82953458108474554
 351   1.8    0.86079282968276671
 352   2      0.88547542598600626
 353   2.2    0.90517870624273966
 354   2.4    0.92105495653509362
 355   2.6    0.93395164268166786
 356   2.8    0.94450240360053817
 357   3      0.95318792074278835
 358   3.2    0.96037753019309191
 359   3.4    0.96635823989417369
 360   3.6    0.97135533107998406
 361   3.8    0.97554722413538364
 362   4      0.97907636417888622
 363   0.40000000000000008
 364   0.99999999999999967
 365   1.9999999999999973
 366
 367
 368
 369 > SuppDists::dinvGauss(2, 1, 1) [1] 0.1098478
 370 > SuppDists::dinvGauss(0.4, 1, 1) [1] 1.005513
 371 > SuppDists::dinvGauss(0.5, 1, 1) [1] 0.8787826
 372 > SuppDists::dinvGauss(0.39, 1, 1) [1] 1.016559
 373 > SuppDists::dinvGauss(0.38, 1, 1) [1] 1.027006
 374 > SuppDists::dinvGauss(0.37, 1, 1) [1] 1.036748
 375 > SuppDists::dinvGauss(0.36, 1, 1) [1] 1.045661
 376 > SuppDists::dinvGauss(0.35, 1, 1) [1] 1.053611
 377 > SuppDists::dinvGauss(0.3, 1, 1) [1] 1.072888
 378 > SuppDists::dinvGauss(0.1, 1, 1) [1] 0.2197948
 379 > SuppDists::dinvGauss(0.2, 1, 1) [1] 0.9005211
 380 >
 381 x = 0.3 [1, 1] 1.0728879234594337  // R SuppDists::dinvGauss(0.3, 1, 1) [1] 1.072888
 382
 383 x = 1   [1, 1] 0.3989422804014327
 384
 385
 386  0 "                NA"
 387  1 "0.3989422804014327"
 388  2 "0.10984782236693059"
 389  3 "0.039418357969819733"
 390  4 "0.016189699458236468"
 391  5 "0.007204168934430732"
 392  6 "0.003379893528659049"
 393  7 "0.0016462878258114036"
 394  8 "0.00082460931140859956"
 395  9 "0.00042207355643694234"
 396 10 "0.00021979480031862676"
 397
 398
 399 [1] "                NA"     " 0.690988298942671"     "0.11539974210409144"
 400  [4] "0.01799698883772935"    "0.0029555399206496469"  "0.00050863023587406587"
 401  [7] "9.0761842931362914e-05" "1.6655279133132795e-05" "3.1243174913715429e-06"
 402 [10] "5.96530227727434e-07"   "1.1555606328299836e-07"
 403
 404
 405 matC(dinvGauss(0:10, 1, 3), digits=17)  df = 3
 406 [1] "                NA"     " 0.690988298942671"     "0.11539974210409144"
 407  [4] "0.01799698883772935"    "0.0029555399206496469"  "0.00050863023587406587"
 408  [7] "9.0761842931362914e-05" "1.6655279133132795e-05" "3.1243174913715429e-06"
 409 [10] "5.96530227727434e-07"   "1.1555606328299836e-07"
 410 $title
 411 [1] "Inverse Gaussian"
 412
 413 $nu
 414 [1] 1
 415
 416 $lambda
 417 [1] 3
 418
 419 $Mean
 420 [1] 1
 421
 422 $Median
 423 [1] 0.8596309
 424
 425 $Mode
 426 [1] 0.618034
 427
 428 $Variance
 429 [1] 0.3333333
 430
 431 $SD
 432 [1] 0.5773503
 433
 434 $ThirdCentralMoment
 435 [1] 0.3333333
 436
 437 $FourthCentralMoment
 438 [1] 0.8888889
 439
 440 $PearsonsSkewness...mean.minus.mode.div.SD
 441 [1] 0.6615845
 442
 443 $Skewness...sqrtB1
 444 [1] 1.732051
 445
 446 $Kurtosis...B2.minus.3
 447 [1] 5
 448
 449   Example: Wald distribution.
 450   (Standard) Wald distribution, mean = 1, scale = 1
 451   1 dx =      0.24890250442652451, x =      0.70924622051646713
 452   2 dx =    -0.038547954953794553, x =      0.46034371608994262
 453   3 dx =   -0.0011074090830291131, x =      0.49889167104373716
 454   4 dx = -9.1987259926368029e-007, x =      0.49999908012676625
 455   5 dx =  -6.346513344581067e-013, x =      0.49999999999936551
 456   dx = 6.3168242705156857e-017 at i = 6
 457    qinvgauss(0.3649755481729598, 1, 1) = 0.50000000000000011
 458   1 dx =       0.6719944578376621, x =       1.3735318786222666
 459   2 dx =     -0.16997432635769361, x =      0.70153742078460446
 460   3 dx =    -0.027865119206495724, x =      0.87151174714229807
 461   4 dx =  -0.00062283290009492603, x =      0.89937686634879377
 462   5 dx = -3.0075104108208687e-007, x =      0.89999969924888867
 463   6 dx = -7.0485322513588089e-014, x =      0.89999999999992975
 464   7 dx =   9.557331866250277e-016, x =      0.90000000000000024
 465   dx = 0 at i = 8
 466    qinvgauss(0.62502320258649202, 1, 1) = 0.89999999999999925
 467   1 dx =   -0.0052296256747447678, x =      0.19483508278446249
 468   2 dx =  6.4699046853900505e-005, x =      0.20006470845920726
 469   3 dx =  9.4123530465288137e-009, x =      0.20000000941235335
 470   4 dx =  2.7739513919147025e-016, x =      0.20000000000000032
 471   dx = 1.5410841066192808e-016 at i = 5
 472    qinvgauss(0.063753567519976254, 1, 1) = 0.20000000000000004
 473   1 dx =                       -1, x =     -0.46073286697416105
 474   2 dx =      0.47665501251497061, x =      0.53926713302583895
 475   3 dx =       -0.171105768635964, x =     0.062612120510868341
 476   4 dx =     0.091490360797512563, x =      0.23371788914683234
 477   5 dx =     0.029410311722649803, x =      0.14222752834931979
 478   6 dx =     0.010761845493592421, x =      0.11281721662666999
 479   7 dx =    0.0019864890597643035, x =      0.10205537113307757
 480   8 dx =  6.8800383732599561e-005, x =      0.10006888207331327
 481   9 dx =  8.1689466405590418e-008, x =      0.10000008168958067
 482   10 dx =   1.133634672475146e-013, x =      0.10000000000011428
 483   11 dx =  5.9588135045224526e-016, x =      0.10000000000000091
 484   12 dx =   3.433223674791152e-016, x =      0.10000000000000031
 485   dx = 9.0763384505974048e-017 at i = 13
 486    qinvgauss(0.0040761113207110162, 1, 1) = 0.099999999999999964
 487
 488
 489      wald_example.vcxproj -> J:\Cpp\MathToolkit\test\Math_test\Debug\wald_example.exe
 490   Example: Wald distribution.
 491   Tolerance = 6.66134e-016
 492   (Standard) Wald distribution, mean = 1, scale = 1
 493   cdf(x, 1,1) 4.1390252102096375e-012
 494   qinvgauss(pinvgauss(x, 1, 1) = 0.020116801973767886, diff = -0.00011680197376788548, fraction = -0.005840098688394274
 495   ____________________________________________________________
 496   wald 1, 1
 497   x =                     0.02, diff x - qinvgauss(cdf) = -0.00011680197376788548
 498   x =      0.10000000000000001, diff x - qinvgauss(cdf) = 8.7430063189231078e-016
 499   x =      0.20000000000000001, diff x - qinvgauss(cdf) = -1.1102230246251565e-016
 500   x =      0.29999999999999999, diff x - qinvgauss(cdf) = 0
 501   x =      0.40000000000000002, diff x - qinvgauss(cdf) = 2.2204460492503131e-016
 502   x =                      0.5, diff x - qinvgauss(cdf) = -1.1102230246251565e-016
 503   x =      0.59999999999999998, diff x - qinvgauss(cdf) = 1.1102230246251565e-016
 504   x =      0.80000000000000004, diff x - qinvgauss(cdf) = 1.1102230246251565e-016
 505   x =      0.90000000000000002, diff x - qinvgauss(cdf) = 0
 506   x =      0.98999999999999999, diff x - qinvgauss(cdf) = -1.1102230246251565e-016
 507   x =                    0.999, diff x - qinvgauss(cdf) = -1.1102230246251565e-016
 508
 509
 510 */
 511
 512
 513