# Goodness-of-fit for discrete distribution - tails issues

In general, I want to see if a set of values comes from a certain discrete distribution (in this case integers from Poisson, or maybe gamma-Poisson). Chi-squared consistently comes up as the best GOF test, but occasionally the test behaves erratically, mostly due to some issues on distribution tails.

How would one handle tail values? They frequently have low expected counts which invalidates the chi-squared test. On the other hand, excluding them might be wrong if there are many observed values in those bins. Is there an exact version that sums the contributions of each bin same as the chi-square version?

Is there a method that avoids the problem of extreme values with low expected counts?

I ultimately opted for using Fisher at the boundaries, then converting each bin to equivalent chi-square.

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

def pearson_chi(x,y):
arr = (x-y)**2/y
return sum(arr)

def neyman_chi(x,y):
vars = np.maximum(x,1)
arr = (x-y)**2/vars
return sum(arr)

def mid_chi(x,y):
arr = y - x + x*np.log(x/y)
return 2*sum(arr)

def gof_chi(fobs, fexp, ddof, chi_fun=mid_chi, thresh=10):
fisher_candidates = fexp <= thresh
chi_candidates = ~fisher_candidates
if sum(chi_candidates) < 2: # chi-squared doesn't work on single-element arrays
fisher_candidates[:] = True
chi_total = 0
else:
chi_total = chi_fun(fobs[chi_candidates], fexp[chi_candidates])

# Handle Fisher candidates individually
fisher_total = 0
sum_obs = sum(fobs)
sum_exp = sum(fexp)
candidate_pairs = zip(fobs[fisher_candidates], fexp[fisher_candidates])
for cand in candidate_pairs:
table = np.array([[cand[0], sum_obs-cand[0]],[cand[1], sum_exp-cand[1]]])
_, pval = scipy.stats.fisher_exact(table)
fisher_total += chi2.ppf(1-pval,1)

# Find full statistic and pval
stat = chi_total + fisher_total
pval = 1 - chi2.cdf(stat, len(fobs)-1-ddof)

return stat, pval

def poisson_fit(samp, average_fun=np.mean, chi_fun=mid_chi, thresh=10):

N = len(samp)

# Find lambda of Poisson
lam = average_fun(samp)

# Find upper edge
observed_edge = round(max(samp)) + 1
theoretical_edge = round(lam + 3*(lam**0.5)) # mean + n*SD
while round(scipy.stats.poisson.pmf(theoretical_edge, lam) * N) > 0:
theoretical_edge += 5
upper_edge = max(observed_edge, theoretical_edge)
edges = np.arange(-0.5, upper_edge + 0.5)

# Perform histogram binning
fobs, _ = np.histogram(samp, bins=np.arange(-0.5, upper_edge + 0.5))
fexp = np.round(scipy.stats.poisson.pmf(np.arange(upper_edge), lam) * N)

# Select only nonzero pairs (to avoid increasing DF too much)
nonzero_bool = (fobs + fexp) > 0
nz_obs = fobs[nonzero_bool]
nz_exp = fexp[nonzero_bool]

return gof_chi(nz_obs, nz_exp, ddof=1, chi_fun=mid_chi)

# Generate data
for i in range(50):
sample = np.random.poisson(lam=4.1, size=30)
print(poisson_fit(sample))


Without looking closely at this, I can't very well judge what you've done, but my suggestion would just be to use Pearson's chi-square test with a simulated distribution instead of relying on an approximation to the chi-squared distribution.

I'm not sure how Fisher's exact test is implemented in SciPy, but it's often implemented via contingency tables with fixed margins. That's probably not a valid assumption for the data that you're looking at.