0
$\begingroup$

Please help me understanding basic features of chi-squared test for distribution comparison. As far as I understand, if we draw two samples from the same distribution, and treat one sample as "observed" and another as "expected" then most of the time the associated chi2 value will be low and the respective p-value to reject the null hypothesis that the two samples come from the same distribution will be high. Therefore, I expect that if I repeat this process many times and test this p-value against some predefined threshold (say 0.05), then only a in a small fraction of times (~5% in our case) the associated p-value will be below the critical threshold. However, in my simulations this null hypothesis is rejected in about 23% of sample pairs. What am I missing? Below is my Python code. In this code, i generate two samples of 100 random integers between 0 and 3. I then calculate the p-value using my implementation of chi2 tests (I verified the correctness of my chi2 test function against an online calculator) and compare this value to a threshold (ALPHA = 0.05). I repeat this process 1000 times and calculate the percentage of times where the p-value was below ALPHA

import numpy as np
import scipy.stats as stats
def chi2testTwoSamples(observed, expected):
    countO = defaultdict(int)
    countE = defaultdict(int)
    nPoints = len(observed)
    assert nPoints == len(expected)
    for o in observed: countO[o] += 1
    for e in expected: countE[e] += 1

    keysO = set(countO.keys())
    keysE = set(countE.keys())
    keys = list(keysO.union(keysE))
    keys.sort()
    countExpected = np.array([countE[k] for k in keys], dtype=float) 
    countObserved = np.array([countO[k] for k in keys], dtype=float) 
    sel = (countExpected != 0)
    countExpected = countExpected[sel]
    countObserved = countObserved[sel]
    chi2 = np.sum((countExpected - countObserved) ** 2 / countExpected)
    df = len(keys) - 1
    pVal = 1.0 - stats.distributions.chi2.cdf(chi2, df)
    return (chi2, pVal)



TIMES = 1000
SAMPLES = 100
ALPHA = 0.05
timesBelowAlpha = 0
np.random.seed(9999)
for i in range(TIMES):
    a = np.random.randint(0, 3, SAMPLES) #draw many integers between 0 and 3
    b = np.random.randint(0, 3, SAMPLES)
    p = chi2testTwoSamples(a, b)[1]
    if p < ALPHA:
        timesBelowAlpha += 1
percentage = 100.0 * timesBelowAlpha / TIMES
print '%.1f%%'%(percentage)
$\endgroup$
2
$\begingroup$

I hardly read Python but what I understand from the code of your chi2testTwoSamples function is that you consider the second sample as "expected values"; that’s not the way it works!

  • Either you compare the four observed values $n_0, n_1, n_2, n_3$ obtained from counts of a single random sample of $N$ elements in $\{0, 1, 2, 3\}$ with the expected values $e_0 = e_1 = e_2 = e_3 = {1\over 4} N$ ;

  • or you have two samples of total size $N$ (resp. $M$) with counts $n_0, n_1, n_2, n_3$ (resp. $m_0, m_1, m_2, m_3$), and assuming that both samples come from the same distribution you infer the frequencies $f_0, f_1, f_2, f_3$ of the four outcomes by $f_i = {1 \over N+M} (n_i + m_i)$ and expected values are $N f_i$ (resp. $M f_i$).

In both cases, the result should follow a $\chi^2$ distribution with $3$ degrees of freedom. The computation of the number of degrees of freedom is different in the two cases, it’s not every time as simple as number of categories - 1...

$\endgroup$
  • $\begingroup$ you understood my code correctly. Can you please clarify what are the N and M variables in your second bullet point? $\endgroup$ – Boris Gorelik Jan 3 '12 at 9:28
  • $\begingroup$ @bgbg The total size of the two samples, $N = n_0 + n_1 + n_2 + n_3$ and $M = m_0 + m_1 + m_2 + m_3$. In your case $N = M = 100$. $\endgroup$ – Elvis Jan 3 '12 at 9:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.