Check out the Pastebin Gadgets Shop. We have thousands of fun, geeky & affordable gadgets on sale :-)Want more features on Pastebin? Sign Up, it's FREE!

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

By: arcanosam on Aug 15th, 2012  |  syntax: Python  |  size: 3.50 KB  |  views: 268  |  expires: Never
Text below is selected. Please press Ctrl+C to copy to your clipboard. (⌘+C on Mac)
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'))
clone this paste RAW Paste Data
Top