# Simulation to estimate the standard deviation of a normal distribution

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? • What's the point of using the Shapiro-Wilks test ? – Sal Mangiafico Dec 23 '19 at 11:23 • I think the normality test is filtering out the high-variance distributions. When you get a large value of y, it convinces the Shapiro test of non-normality, and your code does not consider that standard deviation. Now I must wonder...why do the normality test? – Dave Dec 23 '19 at 12:10 • the normality test has been added because we found that 5% of the generated vectors fail to follow the normal distribution. we still have the same problem without using the normality test. – Hossam Dec 23 '19 at 12:35 • You know the population distribution. The population is normal. – Dave Dec 23 '19 at 16:09 ## 1 Answer 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…

• we still have the same problem without using the normality test – Hossam Dec 23 '19 at 12:38
• No. You should get much closer to $1$ with my code. You are confusing taking the standard deviation for one single sample of $20$ observations and taking the arithmetic mean of the standard deviation of $10^5$ samples of $20$ observations. It is a bit subtle but if you conduct a few simulations, you will be convinced I am quite right. I don't deserve the $-1$. LOL – Mickybo Yakari Dec 23 '19 at 14:03
• Yes you right. I was talking about using the normality test only in my code. It is my first time here and I don't know -1 meaning. – Hossam Dec 25 '19 at 8:22
• Welcome to Cross Validated. I meant that someone had downvoted my answer by clicking on the downward arrow next to my answer. – Mickybo Yakari Dec 25 '19 at 9:48