Goodness-of-fit for discrete distribution - tails issues

Question

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?

Answer 1

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))

Answer 2

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.