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: 197  |  expires: Never
download  |  raw  |  embed  |  report abuse  |  print
Text below is selected. Please press Ctrl+C to copy to your clipboard. (⌘+C on Mac)
  1. # -*- coding: utf-8 -*-
  2. """ this code refer to this article: http://www.codeproject.com/Articles/432194/How-to-Calculate-the-Chi-Squared-P-Value
  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