Simulation to find sampling distribution

How I make this file?

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.

plot of chunk plot

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.

plot of chunk chisq

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} \)?