3
$\begingroup$

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)))
$\endgroup$
2
  • 9
    $\begingroup$ Neither of these estimators is unbiased. Comparing them to one another doesn't seem to have anything at all to do with bias, anyway. $\endgroup$
    – whuber
    Oct 8, 2017 at 22:00
  • $\begingroup$ 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. $\endgroup$ Oct 8, 2017 at 22:15

1 Answer 1

2
$\begingroup$

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.

Motivation for my answer

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.

Histograms from simulation

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

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.