# Expected Value Chi Square distribution

I'm trying to simulate the distribution from the sample variance $$s^2$$ and compare it with the theoretical distribution. Therefore, I perform a fairly simple simulation (upfront, I'm not a mathematician and I tried my very best to write down all mathematical notations as good as possible, but I'm pretty sure that there are some mistakes, so please correct me, I'm happy to learn):

1. $$X$$ is a random variable, sampled from the standard normal distribution: $$X∼N(μ,σ^2)$$ with $$μ = 0$$ and $$σ^2 = 1$$
2. Sample $$n$$-times getting $$n$$ observations
3. Calculate the sample variance $$s^2$$
4. Calculate the theoretical 2.5% and 97.5% quantiles (var_lower_qt_theo and var_upper_qt_theo) using the following inequality: $$\frac{(n−1)s^2}{\chi^2_{0.025,n-1}}<σ^2<\frac{(n−1)s^2}{\chi^2_{0.975,n-1}}$$
5. In order to reproduce the results via a simulation and to generate a distribution for $$σ^2$$, I changed $$\frac{(n−1)s^2}{\sigma^2}∼\chi^2_{n-1}$$ to $$\frac{(n−1)s^2}{\chi^2_{n-1}}∼\sigma^2$$. Am I allowed to do this? I sampled $$m = 100,000$$ times to reproduce the results from point 4 by using $$\frac{(n−1)s^2}{Z}$$, where $$Z ∼ \chi^2_{n-1}$$
6. Calculating the quantiles (var_lower_qt_sim and var_upper_qt_sim) and mean (mean_variance_sim) from these observations
7. Repeat 1-6 for $$n$$ with $$\{2,3,4,...,20\}$$

The theoretical and simulated quantiles match together (except for some small deviations due to the sampling procedure)

However, calculating the mean from these simulated results does not give me back $$s^2$$, instead I get an estimator exactly $$\frac{n-1}{n-3}$$ bigger than $$s^2$$ for $$n>3$$. Can someone explain this systematic deviation? It seems that: $$E(\frac{(n−1)s^2}{Z}) \neq \frac{(n−1)s^2}{E(Z)} = E(s^2)$$

Obviously, I did something wrong. So, I thought maybe I need to change $$\frac{(n−1)s^2}{\sigma^2}∼\chi^2$$ to $$\frac{\chi^2 \sigma^2}{n−1}∼s^2$$ and repeat the steps 1-7 (However, here I would assume that the sample variance is the population variance, what is wrong again....). This time the mean of the simulated results (var_sim2 in the code) is identical to $$s^2$$, however the quantiles are wrong (not shown in the table). So it seems that I need to choose between the correct confidence intervals of sample variance or the correct expected value. I both cases, I wouldn't trust the simulated distribution.

Can someone explain these observations?

EDIT:

Question

To reduce all of this down to one specific topic. I'm looking for the distribution of the sample variance, which has the "correct" quantiles (e.g., contains in 95% of all cases the true population variance) and the "correct" expected value (e.g. expected value is equal to the sample variance)

Result Table:

n sample_variance mean variance_sim quotient var_lower_qt_sim var_lower_qt_theo var_upper_qt_sim var_upper_qt_theo
2 1.396 84600.341 0.000 0.283 0.278 1401.237 1421.894
3 0.307 4.258 0.072 0.083 0.083 11.860 12.115
4 1.249 3.772 0.331 0.401 0.401 17.036 17.358
5 0.993 1.970 0.504 0.357 0.357 8.089 8.201
6 0.286 0.476 0.600 0.111 0.111 1.720 1.720
7 1.115 1.673 0.667 0.464 0.463 5.347 5.408
8 1.396 1.953 0.715 0.613 0.610 5.800 5.782
9 0.955 1.274 0.749 0.434 0.436 3.502 3.504
10 0.382 0.493 0.776 0.181 0.181 1.280 1.274
11 1.626 2.035 0.799 0.795 0.794 5.002 5.007
12 1.301 1.594 0.816 0.653 0.653 3.754 3.751
13 0.731 0.877 0.834 0.374 0.376 1.972 1.991
14 1.285 1.516 0.847 0.677 0.675 3.339 3.335
15 0.744 0.868 0.857 0.400 0.399 1.858 1.850
16 1.536 1.771 0.867 0.839 0.838 3.673 3.679
17 1.334 1.526 0.874 0.742 0.740 3.095 3.089
18 1.668 1.892 0.881 0.939 0.939 3.747 3.748
19 0.854 0.960 0.889 0.487 0.487 1.876 1.867
20 0.955 1.068 0.895 0.551 0.552 2.042 2.037

set.seed(09062024)
sample_var <- NULL
theo_qt_upper <- NULL
theo_qt_lower <- NULL
var_qt_lower <- NULL
var_qt_upper <- NULL
var_mean <- NULL
theo_mean <- NULL
var_mean2 <- NULL
var_sim2 <- NULL

n <- seq(2,20,1)
sd <- 1

for(i in 1:length(n)){
x <- rnorm(n = n[i], 0,sd)
sample_var[i] <- var(x)
var_sim <- (n[i]-1)*sample_var[i]/rchisq(n = 100000, df = n[i]-1)
var_sim2 <- (rchisq(n = 100000, df = n[i]-1)*sample_var[i])/(n[i]-1)

var_qt_lower[i] <- quantile(var_sim, probs = 0.025)
var_qt_upper[i] <- quantile(var_sim, probs = 0.975)
theo_qt_upper[i] <- (n[i]-1)*sample_var[i]/qchisq(p = 0.025, df = n[i]-1)
theo_qt_lower[i] <- (n[i]-1)*sample_var[i]/qchisq(0.975, df = n[i]-1)
var_mean[i] <- mean(var_sim)
var_mean2[i] <- mean(var_sim2)
w <- dchisq(seq(1e-5,250,1e-3),df = n[i]-1)/sum(dchisq(seq(1e-5,250,1e-3),df = n[i]-1))
chi <- seq(1e-5,250,1e-3)

theo_mean[i] <- sum(((n[i]-1)*sample_var[i]/chi)*w)
}

df <- data.frame(n = n,
sample_variance = sample_var,
variance_sim = var_mean,
variance_sim2 = var_mean2,
variance_theo_distr = theo_mean,
quotient = sample_var/theo_mean,
var_lower_qt_sim = var_qt_lower,
var_lower_qt_theo = theo_qt_lower,
var_upper_qt_sim = var_qt_upper,
var_upper_qt_theo = theo_qt_upper)


Your (5) "I changed $$\frac{(n−1)s^2}{\sigma^2}∼\chi^2_{n-1}$$ to $$\frac{(n−1)s^2}{\chi^2_{n-1}}∼\sigma^2$$" does not make much sense as $$\sigma^2$$ does not have a distribution while your $$(n-1)s^2$$ and the associated value from $$\chi^2_{n-1}$$ should be dependent on each other.

Dividing $$(n−1)s^2$$ by an independent $$Z\sim \chi^2_{n-1}$$ will indeed give a different result. $$\frac 1Z$$ has a so-called inverse-chi-squared distribution (I would prefer reciprocal-chi-squared distribution) and $$E\left[\frac 1 Z\right]=\frac1{n-3}$$ when $$n>3$$, which when multiplied by $$(n-1)$$ leads to your $$\frac{n-1}{n-3}$$ observation.

the distribution of the sample variance, which has the "correct" quantiles (e.g., contains in 95% of all cases the true population variance) and the "correct" expected value (e.g. expected value is equal to the sample variance)

and again your final words are a little strange: it would make more sense to be looking for a distribution whose expected value is equal to the population variance if that is what is important to you.

You already have that at the start of your question when looking at $$\frac{(n−1)s^2}{\sigma^2}∼\chi^2_{n-1}$$ so $$s^2$$ has a scaled chi-square distribution with expectation $$\sigma^2$$, the same distribution as $$\frac{\sigma^2}{(n-1)}Z$$ if $$Z\sim \chi^2_{n-1}$$.

Your $$\frac{(n−1)s^2}{\chi^2_{0.025,n-1}}$$ and $$\frac{(n−1)s^2}{\chi^2_{0.975,n-1}}$$ are sensible endpoints for a $$95\%$$ confidence interval, though the former is the upper point and the latter the lower point while you have them the other way round - the reciprocal reverses the inequality.

Here is an R simulation illustrating this for $$n=5, \mu=0, \sigma^2=1$$ and $$95\%$$ confidence intervals, though all of those are easily adjusted if you wish:

samplevariance <- function(n, mu, sigmasq){
X <- rnorm(n, mu, sqrt(sigmasq))
Xbar <- sum(X) / n
ssq <- sum((X - Xbar)^2) / (n-1)
return(ssq)
}

set.seed(2024)
n <- 5
mu <- 0
sigmasq <- 1
confidence <- 0.95
cases <- 10^6
simssq <- replicate(cases, samplevariance(n, mu, sigmasq))
upperci <- simssq * (n-1) / qchisq((1-confidence)/2, n-1)
lowerci <- simssq * (n-1) / qchisq((1+confidence)/2, n-1)

mean(simssq)                                # should be close to sigmasq
# 1.000121
mean(upperci > sigmasq & lowerci < sigmasq) # should be close to confidence
# 0.949927


which is what you are asking for, allowing for simulation noise.

We can plot the distribution of $$s^2$$ (in black below) and compare it to the theoretical scaled chi-square distribution (in red) to see the match.

plot(density(simssq))
curve(dchisq(x*(n-1)/sigmasq, n-1) * (n-1)/sigmasq,


We can even see your $$\frac{n-1}{n-3}$$ result:

mean((n-1) * simssq / rchisq(cases, n-1))
# 1.995052
sigmasq * (n-1) / (n-3)
# 2

• @Mexx My graph is the distribution of $s^2=\frac1{n-1}\sum (X_i-\bar X)^2$ - in this particular case assuming $n=5$ and normal i.i.d. $X_i$ with $\sigma^2=1$. $s^2$ is an unbiased estimator of $\sigma^2$ but is not a MLE of $\sigma^2$, though the slightly smaller $\frac1{n}\sum (X_i-\bar X)^2$ is an MLE of $\sigma^2$. Your "the MLE of $s^2$ for a sample is equal to the expected value of that sample variance" again does not make much sense to me. Commented Jun 13 at 22:46