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
main()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
514