# Estimating population variance through simulation in R

I want to estimate the variance of the exponential distribution with a rate of $\lambda=0.2$.

I'm drawing a sample of 5 exponentials 1000 times, and know that the theoretical variance of my 'population' should be:
$$\sigma^2=\frac{1}{\lambda^2}=\frac{1}{0.2^2} = 25$$ From my understanding, dividing the deviations from the mean by the sample size will give a biased estimate of the population variance, even if the calculation will technically provide the variance of the individual sample. Instead, I should be able to get an unbiased estimate of the variance by accounting for the n-1 degrees of freedom:
$$s^2 = \sum\limits_{i=1}^n \frac{(x_i-\bar{x})^2}{N-1}$$ I now want to explore this bias in R by running two simulations, one accounting for the reduced degrees of freedom and another in which the sum is only divided by $N$.

variance <- NULL
for (i in 1:1000) {
variance <- c(variance, sum((rexp(5, 0.2) - mean(rexp(5, 0.2)))^2)/(5-1))
}
mean(variance)


However, my supposedly unbiased estimator, accounting for the fewer degrees of freedom, ends up being consistently more biased than the uncorrected version, which is expected to be biased.

For example with set.seed(200), I'm getting the following estimates of the population variance:

Dividing by n-1: 35.64
Dividing by n: 28.51

By comparison, I get much better results using R's built-in var() function. However, I would like to write out as many of the calculations by hand as possible, to help improve my understanding.

Am I misunderstanding the theory or am I misunderstanding the functions in my R code (or both)?

• Is $5$ too small sample size for a decent simulation? Also, you look confuse the sample size and the number of simulation times. – Zhanxiong Jul 24 '15 at 16:48
• I tried increasing the sample size, but it only seems to make the difference smaller - the bias is still there. Could you elaborate on how I may be confusing the sample size and number of simulations? – ageil Jul 24 '15 at 16:53
• Oh, I am sorry, you got that right (in your for loop). – Zhanxiong Jul 24 '15 at 17:03
• I am not quite sure what you meant by "bias". If you mean the bias between the simulation result and the true population variance --- it is always there, which is called "intrinsic uncertainty" and can be made closer to the true population value when increasing sample size. If you mean the bias between your own coding and using the built-in function 'var', this is weird, as I did this by myself and the result agree (I see your possible problem, you have two "rexp(5, 0.2)" in your code, they should be the same, try to make an assignment first in your code.) – Zhanxiong Jul 24 '15 at 17:03
• Just an R comment. Use variance[i] <- ... rather then variance <- c(variance, ..., in many cases using c, cbind etc. makes things slow in R so better to avoid them. variance in your example is already a vector so use it as so! – Tim Jul 24 '15 at 17:41

Your code has some errors. I increased the number of iterations in your simulation (to 10k) to get a better approximation of where the simulated distribution is centered. The biggest problem is that in your code, you generate two different sets of data. One set is used for the data, and the other set is used to calculate the mean. To further deepen your understanding, you may also want to try using the known population mean and not using Bessel's correction to estimate the variance. Here is some code:

set.seed(200)
B          <- 10000
varianceN  <- vector(length=B)
varianceN1 <- vector(length=B)
varianceP  <- vector(length=B)
for (i in 1:B) {
data          <- rexp(5, 0.2)
varianceN[i]  <- sum((data-mean(data))^2) / (5)
varianceN1[i] <- sum((data-mean(data))^2) / (5-1)
varianceP[i]  <- sum((data-5)^2)          / (5)
}
mean(varianceN)   # [1] 19.85737
mean(varianceN1)  # [1] 24.82172
mean(varianceP)   # [1] 24.85525