quantile.hpp
Go to the documentation of this file.
1 
12 #ifndef MLPACK_CORE_MATH_QUANTILE_HPP
13 #define MLPACK_CORE_MATH_QUANTILE_HPP
14 
15 #include <mlpack/prereqs.hpp>
16 
17 namespace mlpack {
18 namespace math {
19 
28 inline double ErfInverse(double x)
29 {
30  double w, p;
31 
32  w = -log((1.0 - x) * (1.0 + x));
33 
34  if (w < 6.250000)
35  {
36  w = w - 3.125000;
37  p = -3.6444120640178196996e-21;
38  p = -1.685059138182016589e-19 + p * w;
39  p = 1.2858480715256400167e-18 + p * w;
40  p = 1.115787767802518096e-17 + p * w;
41  p = -1.333171662854620906e-16 + p * w;
42  p = 2.0972767875968561637e-17 + p * w;
43  p = 6.6376381343583238325e-15 + p * w;
44  p = -4.0545662729752068639e-14 + p * w;
45  p = -8.1519341976054721522e-14 + p * w;
46  p = 2.6335093153082322977e-12 + p * w;
47  p = -1.2975133253453532498e-11 + p * w;
48  p = -5.4154120542946279317e-11 + p * w;
49  p = 1.051212273321532285e-09 + p * w;
50  p = -4.1126339803469836976e-09 + p * w;
51  p = -2.9070369957882005086e-08 + p * w;
52  p = 4.2347877827932403518e-07 + p * w;
53  p = -1.3654692000834678645e-06 + p * w;
54  p = -1.3882523362786468719e-05 + p * w;
55  p = 0.0001867342080340571352 + p * w;
56  p = -0.00074070253416626697512 + p * w;
57  p = -0.0060336708714301490533 + p * w;
58  p = 0.24015818242558961693 + p * w;
59  p = 1.6536545626831027356 + p * w;
60  }
61  else if (w < 16.000000)
62  {
63  w = sqrt(w) - 3.250000;
64  p = 2.2137376921775787049e-09;
65  p = 9.0756561938885390979e-08 + p * w;
66  p = -2.7517406297064545428e-07 + p * w;
67  p = 1.8239629214389227755e-08 + p * w;
68  p = 1.5027403968909827627e-06 + p * w;
69  p = -4.013867526981545969e-06 + p * w;
70  p = 2.9234449089955446044e-06 + p * w;
71  p = 1.2475304481671778723e-05 + p * w;
72  p = -4.7318229009055733981e-05 + p * w;
73  p = 6.8284851459573175448e-05 + p * w;
74  p = 2.4031110387097893999e-05 + p * w;
75  p = -0.0003550375203628474796 + p * w;
76  p = 0.00095328937973738049703 + p * w;
77  p = -0.0016882755560235047313 + p * w;
78  p = 0.0024914420961078508066 + p * w;
79  p = -0.0037512085075692412107 + p * w;
80  p = 0.005370914553590063617 + p * w;
81  p = 1.0052589676941592334 + p * w;
82  p = 3.0838856104922207635 + p * w;
83  }
84  else
85  {
86  w = sqrt(w) - 5.000000;
87  p = -2.7109920616438573243e-11;
88  p = -2.5556418169965252055e-10 + p * w;
89  p = 1.5076572693500548083e-09 + p * w;
90  p = -3.7894654401267369937e-09 + p * w;
91  p = 7.6157012080783393804e-09 + p * w;
92  p = -1.4960026627149240478e-08 + p * w;
93  p = 2.9147953450901080826e-08 + p * w;
94  p = -6.7711997758452339498e-08 + p * w;
95  p = 2.2900482228026654717e-07 + p * w;
96  p = -9.9298272942317002539e-07 + p * w;
97  p = 4.5260625972231537039e-06 + p * w;
98  p = -1.9681778105531670567e-05 + p * w;
99  p = 7.5995277030017761139e-05 + p * w;
100  p = -0.00021503011930044477347 + p * w;
101  p = -0.00013871931833623122026 + p * w;
102  p = 1.0103004648645343977 + p * w;
103  p = 4.8499064014085844221 + p * w;
104  }
105  return p * x;
106 }
107 
115 inline double Quantile(double p, double mu = 0.0, double sigma = 1.0)
116 {
117  return mu + sigma * std::sqrt(2.0) * ErfInverse(2 * p - 1);
118 }
119 
120 } // namespace math
121 } // namespace mlpack
122 
123 #endif
double Quantile(double p, double mu=0.0, double sigma=1.0)
Computes the quantile function of Guassian distribution at given probability.
Definition: quantile.hpp:115
Linear algebra utility functions, generally performed on matrices or vectors.
The core includes that mlpack expects; standard C++ includes, Armadillo, cereal, and a few basic mlpa...
double ErfInverse(double x)
Computes the inverse erf function using the rational approximation from Numerical Recipes...
Definition: quantile.hpp:28