# Implementing chi-square in python and testing on scipy's poisson and norm variates

I need to implement Pearson's chi-squared test to test random variates. But I get very different results using different sequence length,degrees of freedom or even seeds. Only a few times sequences pass the test (for demonstration purposes tests are on scipy's generators, not custom). Is my implementation wrong or is it the behaviour of chi-squared test? It looks quite strange that I often fail tests for scipy's functions (and custom too), and the effect by changing seed suprises me too.

Some examples of my implementation output:

For normal distribution generated by scipy with mean 10 and scale Chi2 statistic is 1.48554564247, when critical value is 1.14547622606 (significance alpha=0.05, degrees of freedom = 5, sample size = 200000)

For poisson distribution generated by scipy with mean 10 Chi2 statistic is 3.08213050263, while critical value is 1.14547622606 (significance alpha=0.05, degrees of freedom = 5, sample size = 100000)

And I also tried with a much smaller sample size = 2000:

For poisson distribution generated by scipy with mean 10 Chi2 statistic is 24.5618663076 while critical value is 11.5913052088 (significance alpha=0.05, degrees of freedom = 21)

For n degrees of freedom bins are selected in form of 2 bins (-Infinity;min observed value), (max observed value;+Infinity) and n-2 equal intervals beetween previous two intervals.

UPDATE: Know I think that I make decision to reject using wrong values. F.e. if I have a significance value alpha = 0.05, should I use chi2 percent point function with arguments of alpha and dof or 1-alpha and dof? I think I was mislead by russian Wikipedia page notation

import numpy as np
from scipy.stats import chi2
from scipy.stats import chisquare

def getBins(xmin,xmax,n_bins):
r = np.linspace(xmin,xmax,num=n_bins+1,endpoint = True)
r = r+10**(-10) # including rightmost
r=r-2*10**(-10) # excluding xmin from (-Inf;xmin] bin
return np.concatenate((np.array([float('-inf')]), r, np.array([float('inf')])))

# Calculates probabilities for each bin (a,b] within given cumulative distribution function
def piCalcDecoratorNew(bins, *args):
def real_piCalcDecorator(cdfFunc):
def piCalc(*args):
piA = np.zeros(len(bins)-1)
if len(args)==1:
args = args
piA =cdfFunc(bins,args)
piA[-1] =1-cdfFunc(bins[-2],args)
for i in range(1,len(bins)-2):
piA[i]=cdfFunc(bins[i+1],args)-cdfFunc(bins[i],args)
else: #number of params >1
piA =cdfFunc(bins,*args)
piA[-1] =1-cdfFunc(bins[-2],*args)
for i in range(1,len(bins)-2):
piA[i]=cdfFunc(bins[i+1],*args)-cdfFunc(bins[i],*args)
return piA
return piCalc
return real_piCalcDecorator

# similar to scipy's chisquare()
def chi2statistic(obs = np.array([16, 18, 16, 14, 12, 12], dtype='float'), exp = np.array([16, 16, 16, 16, 16, 8],dtype='float')):
temp = np.square(obs-exp,dtype='float')
with np.errstate(divide='ignore',invalid='ignore'):
temp = temp / exp
temp[exp == 0] = 0
return sum(temp)
# like return chisquare(obs,exp)

def chi2test(df,x, alpha,cdfFunc,*args):
N = len(x)
xmin = min(x)
xmax = max(x)
bins = getBins(xmin,xmax,df-2)
print "Bins for histogram are "
print bins
piCalc = piCalcDecoratorNew(bins,*args)(cdfFunc)
piks = piCalc(*args)
print "Expected probability to be in a bin"
print piks
a = piks*float(N)
b = np.histogram(x,bins)
print "Observed probabilities for bins"
print b/float(N)
print "Chi2 statistic is {0}".format(chi2statistic(b,a))
print "Critical value is {0}".format(chi2.ppf(alpha,df))
return (chi2statistic(b,a),chi2.ppf(alpha,df))


And here's some of my tests and outputs:

# Testing on scipy's norm
from scipy.stats import norm
alpha = 0.05
test_sequence = norm.rvs(loc=10.0, scale=2.0, size=100000, random_state=42)
print chi2test(5,test_sequence,alpha,norm.cdf,10.0,2.0)

Bins for histogram are
[        -inf   1.06879227   7.03191768  12.99504309  18.9581685
inf]
Expected probability to be in a bin
[  3.99216127e-06   6.88950087e-02   8.63972197e-01   6.71250536e-02
3.74819718e-06]
Observed probabilities for bins
[ 0.       0.06836  0.86407  0.06757  0.     ]
Chi2 statistic is 1.48554564247
Critical value is 1.14547622606
(1.4855456424733853, 1.1454762260617695)

# Testing on scipy's poisson
from scipy.stats import poisson
alpha = 0.05
test_sequence = poisson.rvs(mu=10, size=100000, random_state=42)
print chi2test(5,test_sequence,alpha,poisson.cdf,10.0)

Bins for histogram are
[            -inf  -1.00000000e-10   9.00000000e+00   1.80000000e+01
2.70000000e+01              inf]
Expected probability to be in a bin
[  0.00000000e+00   4.57929714e-01   5.34883781e-01   7.18425107e-03
2.25353405e-06]
Observed probabilities for bins
[ 0.       0.45658  0.53583  0.00759  0.     ]
Chi2 statistic is 3.08213050263
Critical value is 1.14547622606
(3.0821305026304886, 1.1454762260617695)

# Testing on scipy's poisson with other parameters
from scipy.stats import poisson
alpha = 0.05
test_sequence = poisson.rvs(mu=10, size=300000, random_state=42)

print chi2test(max(test_sequence),test_sequence,alpha,poisson.cdf,10.0)

Bins for histogram are
[            -inf  -1.00000000e-10   1.08000000e+00   2.16000000e+00
3.24000000e+00   4.32000000e+00   5.40000000e+00   6.48000000e+00
7.56000000e+00   8.64000000e+00   9.72000000e+00   1.08000000e+01
1.18800000e+01   1.29600000e+01   1.40400000e+01   1.51200000e+01
1.62000000e+01   1.72800000e+01   1.83600000e+01   1.94400000e+01
2.05200000e+01   2.16000000e+01   2.26800000e+01   2.37600000e+01
2.48400000e+01   2.59200000e+01   2.70000000e+01              inf]
Expected probability to be in a bin
[  0.00000000e+00   4.99399227e-04   2.26999649e-03   7.56665496e-03
1.89166374e-02   3.78332748e-02   6.30554580e-02   9.00792257e-02
1.12599032e-01   1.25110036e-01   1.25110036e-01   1.13736396e-01
9.47803301e-02   1.24985051e-01   3.47180696e-02   2.16987935e-02
1.27639962e-02   7.09110899e-03   3.73216263e-03   1.86608131e-03
8.88610150e-04   4.03913704e-04   1.75614654e-04   7.31727725e-05
2.92691090e-05   1.54267384e-05   2.25353405e-06]
Observed probabilities for bins
[  0.00000000e+00   5.03333333e-04   2.24000000e-03   7.73666667e-03
1.89433333e-02   3.75066667e-02   6.30300000e-02   9.04766667e-02
1.13150000e-01   1.24163333e-01   1.25440000e-01   1.13123333e-01
9.52500000e-02   1.25460000e-01   3.40800000e-02   2.17266667e-02
1.27300000e-02   7.00333333e-03   3.82000000e-03   1.97000000e-03
9.30000000e-04   4.46666667e-04   1.56666667e-04   6.00000000e-05
4.33333333e-05   1.00000000e-05   0.00000000e+00]
Chi2 statistic is 20.8855527231
Critical value is 16.1513958497
(20.885552723131507, 16.151395849664102)

• Can you clarify this for people who don't read Python, & make it more software-neutral? If you are asking for code check, that would be off topic here (it would belong on the Code Review SE site). Questions about how the chi-squared test works are on topic here, but you may want to bring those front and center. – gung - Reinstate Monica Mar 20 '16 at 12:49
• @zero I clarified outputs, put clear question in bold. What else I need to do? – Pavel Akapian Mar 20 '16 at 13:37
• Good. Please also review the Zen of Python, because that is very Unpythonic code. – user46925 Mar 20 '16 at 13:41
• my guess is that df=max(test_sequence) is wrong. It should be number of bins - 1 if I remember correctly. – Josef Mar 20 '16 at 18:04
• @user333700 I found recommendations to take number of bins as number of unique discrete values. But it was for finite random variables. In my case random variable is infinite, so i think that the best df is max(test_sequence)+1+1 (for all unique values observed starting from zero + 1 for tail) – Pavel Akapian Mar 20 '16 at 19:09

• ppf is the inverse cdf, so it should be 1-alpha for the upper tail of the chisquare distribution. However, you can also use the corresponding survival function, isf for the critical value and sf for the pvalue, for those it should be alpha (0.05) directly. – Josef Mar 20 '16 at 22:22