Sampling distribution for a statistic sometimes is hard to find. Simulation is a good substitute to provide you a direct taste of how the distribution looks like. Well, before starting, you need to make sure you have the followings:
Now let's see an example where I will demonstrate that when you have i.i.d. sample from Normal distribution, \( \frac{(n-1)S_n^2}{\sigma^2} \sim \chi_{n-1} \).
Let's first do some preperation.
mu <- 4
sigma2 <- 1
n <- 50 # sample size
n.rep <- 1000
record_stat <- rep(0, n.rep) # Initialize a vector to record the statistics.
The next is the loop. Let R do it!
#### BEGIN Loop
for (i in 1:n.rep) {
x <- rnorm(n, mu, sqrt(sigma2)) # sample size n and the population is N(mu,sigma2).
S2 <- sum((x - mean(x))^2)/(n - 1)
record_stat[i] <- (n - 1) * S2/sigma2
}
#### END Loop
Ok, now we have the statistics calculated and recorded in a vector named record_stat, let's plot it out and then see how it looks like.
hist(record_stat, freq = F, main = "histogram of the statistic") ## plot histogram
lines(density(record_stat)) ## add it's estimated density.
To compare with the chi-square(n-1), let's simulate some chi-square(n-1) random variable and add the density to the existing plot for comparision.
hist(record_stat, freq = F, main = "histogram of the statistic") ## plot histogram
lines(density(record_stat))
lines(density(rchisq(n.rep, (n - 1))), col = "red") ## add chisq density.
Hurray!! We're done! Now think about this, for the example we have, do we need the conditions that \( n \rightarrow \infty \) in order to get the \( \chi_{n-1} \)?