FilesExpand file tree

 main

/

_statisticsmodule.c

Copy path

Blame

More file actions

Blame

More file actions

Latest commit

 

History

History

History

151 lines (137 loc) · 4.81 KB

 main

/

_statisticsmodule.c

Copy path

Top

File metadata and controls

• Code
• Blame

151 lines (137 loc) · 4.81 KB

Raw

Copy raw file

Download raw file

Open symbols panel

Edit and raw actions

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

54

55

56

57

58

59

60

61

62

63

64

65

66

67

68

69

70

71

72

73

74

75

76

77

78

79

80

81

82

83

84

85

86

87

88

89

90

91

92

93

94

95

96

97

98

99

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

/* statistics accelerator C extension: _statistics module. */

#include "Python.h"

#include "clinic/_statisticsmodule.c.h"

/*[clinic input]

module _statistics

[clinic start generated code]*/

/*[clinic end generated code: output=da39a3ee5e6b4b0d input=864a6f59b76123b2]*/

/*

* There is no closed-form solution to the inverse CDF for the normal

* distribution, so we use a rational approximation instead:

* Wichura, M.J. (1988). "Algorithm AS241: The Percentage Points of the

* Normal Distribution". Applied Statistics. Blackwell Publishing. 37

* (3): 477–484. doi:10.2307/2347330. JSTOR 2347330.

*/

/*[clinic input]

_statistics._normal_dist_inv_cdf -> double

p: double

mu: double

sigma: double

/

[clinic start generated code]*/

static double

_statistics__normal_dist_inv_cdf_impl(PyObject *module, double p, double mu,

double sigma)

/*[clinic end generated code: output=02fd19ddaab36602 input=24715a74be15296a]*/

{

double q, num, den, r, x;

if (p <= 0.0 || p >= 1.0) {

goto error;

}

q = p - 0.5;

if(fabs(q) <= 0.425) {

r = 0.180625 - q * q;

// Hash sum-55.8831928806149014439

num = (((((((2.5090809287301226727e+3 * r +

3.3430575583588128105e+4) * r +

6.7265770927008700853e+4) * r +

4.5921953931549871457e+4) * r +

1.3731693765509461125e+4) * r +

1.9715909503065514427e+3) * r +

1.3314166789178437745e+2) * r +

3.3871328727963666080e+0) * q;

den = (((((((5.2264952788528545610e+3 * r +

2.8729085735721942674e+4) * r +

3.9307895800092710610e+4) * r +

2.1213794301586595867e+4) * r +

5.3941960214247511077e+3) * r +

6.8718700749205790830e+2) * r +

4.2313330701600911252e+1) * r +

1.0);

if (den == 0.0) {

goto error;

}

x = num / den;

return mu + (x * sigma);

}

r = (q <= 0.0) ? p : (1.0 - p);

if (r <= 0.0 || r >= 1.0) {

goto error;

}

r = sqrt(-log(r));

if (r <= 5.0) {

r = r - 1.6;

// Hash sum-49.33206503301610289036

num = (((((((7.74545014278341407640e-4 * r +

2.27238449892691845833e-2) * r +

2.41780725177450611770e-1) * r +

1.27045825245236838258e+0) * r +

3.64784832476320460504e+0) * r +

5.76949722146069140550e+0) * r +

4.63033784615654529590e+0) * r +

1.42343711074968357734e+0);

den = (((((((1.05075007164441684324e-9 * r +

5.47593808499534494600e-4) * r +

1.51986665636164571966e-2) * r +

1.48103976427480074590e-1) * r +

6.89767334985100004550e-1) * r +

1.67638483018380384940e+0) * r +

2.05319162663775882187e+0) * r +

1.0);

} else {

r -= 5.0;

// Hash sum-47.52583317549289671629

num = (((((((2.01033439929228813265e-7 * r +

2.71155556874348757815e-5) * r +

1.24266094738807843860e-3) * r +

2.65321895265761230930e-2) * r +

2.96560571828504891230e-1) * r +

1.78482653991729133580e+0) * r +

5.46378491116411436990e+0) * r +

6.65790464350110377720e+0);

den = (((((((2.04426310338993978564e-15 * r +

1.42151175831644588870e-7) * r +

1.84631831751005468180e-5) * r +

7.86869131145613259100e-4) * r +

1.48753612908506148525e-2) * r +

1.36929880922735805310e-1) * r +

5.99832206555887937690e-1) * r +

1.0);

}

if (den == 0.0) {

goto error;

}

x = num / den;

if (q < 0.0) {

x = -x;

}

return mu + (x * sigma);

error:

PyErr_SetString(PyExc_ValueError, "inv_cdf undefined for these parameters");

return -1.0;

}

static PyMethodDef statistics_methods[] = {

_STATISTICS__NORMAL_DIST_INV_CDF_METHODDEF

{NULL, NULL, 0, NULL}

};

PyDoc_STRVAR(statistics_doc,

"Accelerators for the statistics module.\n");

static struct PyModuleDef_Slot _statisticsmodule_slots[] = {

{0, NULL}

};

static struct PyModuleDef statisticsmodule = {

PyModuleDef_HEAD_INIT,

"_statistics",

statistics_doc,

0,

statistics_methods,

_statisticsmodule_slots,

NULL,

NULL,

NULL

};

PyMODINIT_FUNC

PyInit__statistics(void)

{

return PyModuleDef_Init(&statisticsmodule);

}