• API
• FAQ
• Tools
• Trends
• Archive
daily pastebin goal
69%
SHARE
TWEET

# How to Calculate the Chi-Squared P-Value Python version

arcanosam Aug 15th, 2012 442 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. # -*- coding: utf-8 -*-
3.
4. Points of interest
5.
6. - Well, I cannot see how to be made the original code in more Pythonic way... or improve legibility, if anyone could help...
7.
8. It was need to use only Decimal, because float overflow in loop on lines 69 to 74. Not throw an error, but cause NaN and Infinite in some these variables (recomend debug to see for yourself)
9.
10. - So Decimal seems to providing much more accuracy in calculations as far I could see... more than C++ Builder on my runs
11.
12. -getcontext function is used only because I could use power method instead pow function. I see an advantage because I don't need to cast float to Decimal when using pow.
13.
14. - same thing to Decimal exp method, only Decimal data used
15.
16. - Bonus: here you found a C/C++ version tested on C++ Builder XE: http://pastebin.com/Jv85Eu7Y
17.
18. - it's almost equal the original C source, exception to IsNan and IsInfinite wich from VCL Math.
19.
20. """
21.
22. from decimal import Decimal, getcontext
23. #from math import gamma
24.
25. A = Decimal('15')
26.
27. def mygamma(z):
28.
29.     """
30.        The constant SQRT2PI is defined as sqrt(2.0 * PI);
31.        For speed the constant is already defined in decimal
32.        form.  However, if you wish to ensure that you achieve
33.        maximum precision on your own machine, you can calculate
34.        it yourself using (sqrt(atan(1.0) * 8.0))
35.
36.    """
37.
38.     #const long double SQRT2PI = sqrtl(atanl(1.0) * 8.0);
39.     SQRT2PI = Decimal('2.5066282746310005024157652848110452530069867406099383')
40.
41.     f = Decimal('1')
42.     sum_v = SQRT2PI
43.
44.     sc = getcontext().power(z+A,z+Decimal('0.5'))
45.
46.     sc *= Decimal(Decimal('-1') * (z+A)).exp()
47.
48.     sc /= z
49.
50.     for k in range(1,15):
51.         z+=Decimal('1')
52.         ck = getcontext().power(A - Decimal(k) , Decimal(k) - Decimal('0.5'))
53.         ck *= Decimal(A -Decimal(k)).exp()
54.         ck /= f
55.
56.         sum_v += (ck / z)
57.
58.         f *= (Decimal('-1') * k)
59.
60.     return sum_v * sc
61.
62. def approx_gamma(z):
63.
64.     RECIP_E = Decimal('0.36787944117144232159552377016147')
65.     TWOPI = Decimal('6.283185307179586476925286766559')
66.
67.     d = Decimal('1.0') / (Decimal('10.0') * z)
68.     d = Decimal('1.0') / ((Decimal('12') * z) - d)
69.     d = (d + z) * RECIP_E
70.     d = getcontext().power(d,z)
71.     d *= Decimal.sqrt(TWOPI / z)
72.
73.     return d
74.
75. def igf(s, z):
76.
77.     if z < Decimal('0'):
78.         return Decimal('0')
79.
80.     sc = Decimal('1') / s
81.     sc *= getcontext().power(z,s)
82.     sc *= Decimal(-z).exp()
83.
84.     sum_v = Decimal('1')
85.     nom = Decimal('1')
86.     denom = Decimal('1')
87.
88.     for i in range(0,200):
89.         nom *= z
90.         s+=Decimal('1')
91.         denom *= s
92.
93.         sum_v += (nom / denom)
94.
95.     return sum_v * sc
96.
97.
98. def chisqr(dof, cv):
99.
100.     if cv < Decimal('0') or dof < Decimal('1'):
101.         return Decimal('0')
102.
103.     k = dof * Decimal('0.5')
104.     x = cv * Decimal('0.5')
105.
106.     if dof == Decimal('2'):
107.         print Decimal(Decimal('-1') * x).exp()
108.         return
109.
110.     pvalue = igf(k,x)
111.
112.     if pvalue.is_nan() or pvalue.is_infinite() or pvalue <= 1e-8:
113.         return 1e-14
114.
115.     pvalue /= mygamma(k)       # res: 0.0636423441336337067152539023
116.     #pvalue /= approx_gamma(k)  # res: 0.0636423441307368208316220413
117.     #pvalue = Decimal(gamma(k)) # res: -2.670673575075121885604374585E+212
118.
119.
120.     print 'Dof: %s Cv: %s = %s' % (dof, cv, Decimal('1') - pvalue)
121.
122. if __name__ == '__main__':
123.
124.     chisqr(Decimal('255'), Decimal('290.285192'))
RAW Paste Data
Top