views:

66

answers:

1

Hi, I am trying to write a chi square goodness-of-fit test for Beta distribution from scratch, without using any external functions. The code below reports '1' for a fit, even though kstest from scipy.stats returns a zero. Data is distributed normally, so my function should also return zero.

import numpy as np
from scipy.stats import chi2
from scipy.stats import beta
from scipy.stats import kstest
from scipy.stats import norm

preds = norm.rvs(5,2,size=200)
preds.sort()

bin_size = 30
bins = np.linspace(0,10,bin_size)
counts = np.digitize(preds, bins)
mean = 5
var = 2

sum = 0
for i in range(len(bins)-1):
    p = beta.cdf(bins[i+1], mean, var) - beta.cdf(bins[i], mean, var)  
    freq = len(counts[counts==i]) / float(len(counts))    
    sum = sum + ((freq - p)**2)/p

dof = len(counts)-2
pval = 1 - chi2.cdf(sum, dof)
print pval

In the code, I create bins, measure frequencies based on the bins, calculate expected frequency using Beta distribution CDF, and sum it up resulting in the X^2 test statistic.

The kstest call is

print kstest(preds, 'beta', [mean, var])

What am I doing wrong here?

Thanks,

+1  A: 

Problem was with the DOF definition:

dof = len(preds)-2

is the correct choice. Also, I had to reduce bin size to 15 in order to get consistent '0' result. It is known that Chi^2 tests are sensitive on bin size.