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:

- The population where your sample comes from. For example, \[ N(\mu=4,\sigma=1) \]
- The sample size \( n \). Hence you will have \[ X_1, \ldots, X_n \]
- The statistic you are intersted in. Of course!! For example, \[ \frac{(n-1)S_n^2}{\sigma^2} \]
- The number of replications \( n.rep \).

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