# Simulated chi-square distribution doesn't match theoretical

Can someone explain why the distribution of Chi-square values I'm getting (using Pearson goodness-of-fit test) doesn't match the expected Chi-sqaure distribution?

The test seems in this case to be massively conservative...my type 1 error rate never gets close to the .05 alpha level.

ncell=5
npercell=100

Chi.values=P.values=c()
for(i in 1:10000){
observed=rbinom(ncell,npercell,.5) #all cells are sampled with the same proportion
expected=rep(sum(observed)/ncell,ncell) #expected values are mean value of observed
Chi.values[i]=sum(((observed-expected)^2)/expected)
P.values[i]=1-pchisq(Chi.values[i],ncell-1)
}

hist(Chi.values,freq = F,breaks=seq(0,100,.1),xlim=c(0,max(Chi.values)))
#my observed Chi Squares

#the theoretical distribution of Chi squares
mean(P.values<.05) # proportion of type 1 errors


You are applying a test for "all cells have same proportion" however your sampling scheme is not appropriate.

Here is an appropriate scheme with a constant total size $N = 100$ ; I use a multinomial draw:

ncell <- 5
size <- 100

N <- 1e4
chi.values <- p.values <- numeric(N)
for(i in 1:N) {
expe <- rep(size/ncell, ncell)
obse <- rmultinom(1, size, prob=rep(1/ncell, ncell))
chi.values[i] <- chi2 <- sum( (obse-expe)**2/expe )
p.values[i] <- pchisq(chi2, df = ncell-1, lower.tail=FALSE)
}


The qq-plot of the $p$-values is ok :

plot(sort(p.values), pch=".") You can also do that with a varying size:

ncell <- 5; size <- 100; N <- 1e4
chi.values <- p.values <- numeric(N)
for(i in 1:N) {
size <- rpois(1, lambda=100)
expe <- rep(size/ncell, ncell)
obse <- rmultinom(1, size, prob=rep(1/ncell, ncell))
chi.values[i] <- chi2 <- sum( (obse-expe)**2/expe )
p.values[i] <- pchisq(chi2, df = ncell-1, lower.tail=FALSE)
}
plot(sort(p.values), pch=".") In your sampling scheme, you're drawing five independent values from a binomial $\text{Bin}(n = 100, p = 0.5)$. This is not the same as choosing a total size $N$ and them randomly distributing $N$ elements in five bins, which gives five non-independent values.

Post Scriptum You are computing $$T = {1\over \overline X} \sum_{i=1}^5 \left(X_i - \overline X\right)^2$$ with $X_i \sim \text{Bin}(n=100, p=0.5)$. Using a normal approximation for the binomial, we get easily that $\sum_{i=1}^5 \left(X_i - \overline X\right)^2 \sim 25\times \chi^2(4)$ and $\overline X \sim \mathcal N(50, 5)$, and these two are independent. Thus the law of the statistic you compute is the quotient of a scale chi-square by a normal. You can check it

> qqplot(25*rchisq(10000, df=4)/rnorm(10000, mean=50, sd=sqrt(5)), Chi.values)
> abline(0,1,col="red") • +1 The explanation in the post script makes this thread valuable beyond the R code of the question. – Glen_b Mar 7 '15 at 0:10
• Thank you @Glen_b ;) I like this question, I would be happy if my students were able to spot the mistake in the original code and to give the post-scriptum as well... – Elvis Mar 7 '15 at 8:54