# Unbiased estimator of standard deviation. Simulation results

As an exercise I decided to check if unbiased estimator of standard deviation of sample is giving better results than biased estimator . So far it looks that only in aprox 55% cases. Am I doing something wrong or this is normal?

My methodology:

• generate sample of 100 numbers from range 1 to 1000
• 10 000 times choose 10 numbers from above
• for each 10 numbers calculate biased and unbiased estimator of standard deviation
• check how many times unbiased estimator was closer to standard deviation in population (in comparison to biased estimator).

My code:

import numpy as np
import random
from math import sqrt

# Lets generate set of 100 random numbers
myarray = np.random.randint(1,1000,100)
myarray_sd = np.std(myarray)

right = 0
wrong= 0

for k in range(10000):
# lets choose sample of numbers
sample= random.sample(set(myarray),10)

# calculate mean
sample_mean=np.mean(sample)

# calculate sd for n and n-1
sample_sd_n = sqrt(sum((sample-sample_mean)**2)/len(s))
sample_sd_n1 = sqrt(sum((sample-sample_mean)**2)/(len(s)-1))

# callculate diffrences between both sd and sd in population
res_n= abs(sample_sd_n - myarray_sd)
res_n1=abs(sample_sd_n1 - myarray_sd)

# check if std calculated using n-1 is more accurate
if res_n1<res_n:
right +=1
else:
wrong +=1

print ('The theory is correct in: %f cases' % (round(right/(right+wrong),2)))

• Neither of these estimators is unbiased. Comparing them to one another doesn't seem to have anything at all to do with bias, anyway.
– whuber
Oct 8, 2017 at 22:00
• What is your population distribution? Is it uniform on 1-1000 or something else such as uniform over the selected 10? How are you computing the sample variance (i.e. what is N? What is the population variance? Unbiasedness means that the sample variance averages in large samples the population variance. It say nothing about how often it is closer to the population variance than some other estimate. Oct 8, 2017 at 22:15

I thought adding an answer might be useful for anyone else stumbling across this question.

As is mentioned in the comments, both the estimators for the standard deviation you have mentioned in the question are biased. See this question for great explanations of this.

In the original question a simulation is mentioned. Simulations can be a good way to check if we are unsure about things.

## Simulation

I take $$10,000$$ samples of $$10$$ points from a normal distribution with mean $$=0$$ and standard deviation $$=1$$. Then for each of these $$10,000$$ samples apply the estimators $$s$$ and $$s_N$$ (defined in the original question, $$s$$ has $$N-1$$ in the denominator $$s_N$$ has $$N$$ in the denominator) to estimate the standard deviation.

The results of both of these are plotted in the histograms above. The black lines show the mean of the estimates for the standard deviation. The orange lines shows us the true value of the standard deviation.

We see that in both cases the means of the estimates are less than the true value. Based on this simulation we might conclude that the expected value for both $$s$$ and $$s_N$$ is strictly less than the true value of $$\sigma = 1$$.

## Exercises/Things to explore

It might be fun/useful to experiment with the code below. For example use the biased and unbiased estimators for the variance instead and see how the plots above change.

## Python Code

import numpy as np
import matplotlib.pyplot as plt
mean=0
standard_deviation=1
#create a 10000 samples containign 10 points
data_samples = np.random.normal(loc=mean, scale=standard_deviation, size=(10000,10))
#For each of our 1000 samples, apply each of our estimators.
sn_results = np.std(data_samples, axis=1)
s_results = np.std(data_samples, axis=1, ddof=1)
#mean of each estimator
mean_sn = np.mean(sn_results)
mean_s = np.mean(s_results)

#plotting
fig, axs = plt.subplots(2, sharex=True, figsize=(15,7))
for a in fig.axes:
a.tick_params(
axis='x',           # changes apply to the x-axis
which='both',       # both major and minor ticks are affected
bottom=True,
top=False,
labelbottom=True)
fig.suptitle('Comparison of the two estimators for the standard deviation (true value σ=1)')
axs[0].hist(sn_results)
axs[0].axvline(mean_sn, c="black", label="Mean of estimator s_n")
axs[0].axvline(1, c="orange", label="Actual")
axs[0].legend(loc='upper left')
axs[1].hist(s_results)
axs[1].axvline(mean_s, c="black", label="Mean of estimator s")
axs[1].axvline(1, c="orange", label="Actual")
axs[1].legend(loc='upper left')