1- /* statistics accelerator C extensor : _statistics module. */
1+ /* statistics accelerator C extension : _statistics module. */
22
33#include "Python.h"
44#include "structmember.h"
@@ -10,11 +10,13 @@ module _statistics
1010[clinic start generated code]*/
1111/*[clinic end generated code: output=da39a3ee5e6b4b0d input=864a6f59b76123b2]*/
1212
13-
14- static PyMethodDef speedups_methods [] = {
15- _STATISTICS__NORMAL_DIST_INV_CDF_METHODDEF
16- {NULL , NULL , 0 , NULL }
17- };
13+ /*
14+ * There is no closed-form solution to the inverse CDF for the normal
15+ * distribution, so we use a rational approximation instead:
16+ * Wichura, M.J. (1988). "Algorithm AS241: The Percentage Points of the
17+ * Normal Distribution". Applied Statistics. Blackwell Publishing. 37
18+ * (3): 477–484. doi:10.2307/2347330. JSTOR 2347330.
19+ */
1820
1921/*[clinic input]
2022_statistics._normal_dist_inv_cdf -> double
@@ -34,7 +36,7 @@ _statistics__normal_dist_inv_cdf_impl(PyObject *module, double p, double mu,
3436 // Algorithm AS 241: The Percentage Points of the Normal Distribution
3537 if (fabs (q ) <= 0.425 ) {
3638 r = 0.180625 - q * q ;
37- // Hash sum AB: 55.88319 28806 14901 4439
39+ // Hash sum- 55.8831928806149014439
3840 num = (((((((2.5090809287301226727e+3 * r +
3941 3.3430575583588128105e+4 ) * r +
4042 6.7265770927008700853e+4 ) * r +
@@ -54,11 +56,11 @@ _statistics__normal_dist_inv_cdf_impl(PyObject *module, double p, double mu,
5456 x = num / den ;
5557 return mu + (x * sigma );
5658 }
57- r = q <= 0.0 ? p : 1.0 - p ;
59+ r = ( q <= 0.0 ) ? p : ( 1.0 - p ) ;
5860 r = sqrt (- log (r ));
5961 if (r <= 5.0 ) {
6062 r = r - 1.6 ;
61- // Hash sum CD: 49.33206 50330 16102 89036
63+ // Hash sum- 49.33206503301610289036
6264 num = (((((((7.74545014278341407640e-4 * r +
6365 2.27238449892691845833e-2 ) * r +
6466 2.41780725177450611770e-1 ) * r +
@@ -77,7 +79,7 @@ _statistics__normal_dist_inv_cdf_impl(PyObject *module, double p, double mu,
7779 1.0 );
7880 } else {
7981 r -= 5.0 ;
80- // Hash sum EF: 47.52583 31754 92896 71629
82+ // Hash sum- 47.52583317549289671629
8183 num = (((((((2.01033439929228813265e-7 * r +
8284 2.71155556874348757815e-5 ) * r +
8385 1.24266094738807843860e-3 ) * r +
@@ -96,23 +98,30 @@ _statistics__normal_dist_inv_cdf_impl(PyObject *module, double p, double mu,
9698 1.0 );
9799 }
98100 x = num / den ;
99- if (q < 0.0 ) x = - x ;
101+ if (q < 0.0 ) {
102+ x = - x ;
103+ }
100104 return mu + (x * sigma );
101105}
102106
107+
108+ static PyMethodDef statistics_methods [] = {
109+ _STATISTICS__NORMAL_DIST_INV_CDF_METHODDEF
110+ {NULL , NULL , 0 , NULL }
111+ };
112+
103113static struct PyModuleDef statisticsmodule = {
104114 PyModuleDef_HEAD_INIT ,
105115 "_statistics" ,
106116 _statistics__normal_dist_inv_cdf__doc__ ,
107117 -1 ,
108- speedups_methods ,
118+ statistics_methods ,
109119 NULL ,
110120 NULL ,
111121 NULL ,
112122 NULL
113123};
114124
115-
116125PyMODINIT_FUNC
117126PyInit__statistics (void )
118127{
0 commit comments