2
$\begingroup$

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)
$\endgroup$
0

1 Answer 1

5
$\begingroup$

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.

You asked 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)

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, 
      from=0, to=6*sigmasq, col="red", add=TRUE) 

enter image description here

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
$\endgroup$
5
  • $\begingroup$ Thanks a lot for your detailed answer and also your explanations regarding my wrong mathematical notations! That already helped a lot! $\endgroup$
    – Mexx
    Commented Jun 13 at 20:17
  • $\begingroup$ Second comment (I was too slow with editing the first one) My questions might sound strange, but I'm actually looking for the sample variance distribution. If I get confidence intervals and a maximum likelihood estimator (MLE) for s², I should also get a distribution for a sample variance. I thought that the MLE of s² for a sample is equal to the expected value of that sample variance, right? So, again instead of the population variance distribution, I'm looking for the sample variance distribution. $\endgroup$
    – Mexx
    Commented Jun 13 at 20:37
  • 1
    $\begingroup$ @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. $\endgroup$
    – Henry
    Commented Jun 13 at 22:46
  • $\begingroup$ Sorry, I was not precise enough in my formulation and thanks for pointing that out. What I wanted to say is, in your simulation your repeatedly sample from the population with known µ and σ² and calculate s². I want to get the probability distribution of s² for one sample, when the true µ and σ² are not known. $\endgroup$
    – Mexx
    Commented Jun 14 at 8:40
  • $\begingroup$ Another thought. Was my initial simulation with the reciprocal-chi-squared distribution actually the distribution I am looking for? The simulated confidence intervals match the calculated confidence intervals for s² almost perfectly (see my result table). I just got unsure, because I couldn't find my initial point estimator for s² in this distribution anymore by averaging over all simulated values. However, as you pointed out, the expected value from a reciprocal-chi-squared distribution is 1/(n-3) and this would allow me to calculate from the distribution back to s² again. $\endgroup$
    – Mexx
    Commented Jun 14 at 12:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.