Simulation to estimate the standard deviation of a normal distribution

Question

I am trying to estimate the mean and the standard deviation by simulation from a normal distribution. I have the following code in R:

mu=c()
sigma=c()

for(i in 1:100000){
  y=rnorm(20,0,1)
  sh=shapiro.test(y)
  if(sh$p.value>0.05){
    mu[i]=mean(y)
    sigma[i]=sd(y)
  } 
}
mean(mu,na.rm = T)
mean(sigma, na.rm = T)

The problem now is that the mean values of the standard deviation are always less than the true value entered, should it be sometimes greater than the true value because it is an unbiased estimator in case of normality?

Answer 1

This is because the sample standard deviation is a biased estimator of the standard deviation as explained here: https://en.wikipedia.org/wiki/Unbiased_estimation_of_standard_deviation.

Let $S^2$ denote the Bessel-corrected sample variance and $S$ denote its square root, i.e. the Bessel-corrected sample standard deviation. It is well known that the statistic $$ S^2=\frac{1}{n-1} \sum_{i=1}^{n} \bigg( X_i - \bar{X} \bigg)^2$$ is an unbiased estimator of the variance. So, why isn't its square root an unbiased estimator of the standard deviation? Because, by Jensen's inequality, $\mathbb{E}[S] \leq \sqrt{\mathbb{E}[S^2]}$. Actually, because $S$ is a non-degenerate random variable, Jensen's inequality is strict!

If you are interested in an unbiased estimator, use the following code inspired by the Results for the normal distribution section from the same article:

# n: sample size
# returns the correction factor for a normal sample of size n
corr.factor <- function(n){
  return(sqrt(2/(n-1))*gamma(n/2)/gamma((n-1)/2))
}

corr <- corr.factor(20)

mu=c()
sigma=c()

for(i in 1:100000){
  y=rnorm(20,0,1)
  sh=shapiro.test(y)
  if(sh$p.value>0.05){
    mu[i]=mean(y)
    sigma[i]=sd(y)/corr
  } 
}
mean(mu,na.rm = T)
mean(sigma, na.rm = T)

EDIT: The problem is not connected with the Shapiro-Wilk bit. It is due to the fact that the expectation operator does not commute with non-linear functions. If you use my code with the correction factor, there shouldn't be a problem anymore.

Of course, the more 'Gaussian' the data as per Shapiro-Wilk, the more accurate my proposal because the correction factor is meant to correct the bias for normal samples but this is not the root cause of the systematic underestimation. The bottom line is that the arithmetic mean of the $10^5$ sample standard deviations is not the same as the square root of the arithmetic mean of the sample variances…