Python

Point and interval estimation in Python

Let us consider a random sample $x_1,x_2,\dots,x_N$ drawn from a normal distribution with unknown mean $\mu$ and unknown variance $\sigma^2$. The best point estimates for the mean and the variance are given by the sample mean\begin{equation}
\hat\mu = \frac{1}{N}\sum_{i=1}^{N} x_i
\end{equation}
and the sample variance
\begin{equation}
\hat\sigma^2 = \frac{1}{N-1}\sum_{i=1}^{N} (x_i-\hat\mu)^2\,,
\end{equation}
respectively. The sample mean and the sample variance and interval estimates for the mean and the variance may be computed easily using Python.

#!/usr/bin/env python
# -*- coding: utf-8 -*-

from pylab import *
import scipy.stats 

N=1000.      # sample size
gamma=0.95  # confidence level

mu=10.      # true mean
sigma=2.    # true standard diviation
x=randn(N)*sigma+mu  # surrogate data

mu_hat=mean(x)           # sample mean
sigma_hat=std(x, ddof=1) # sample standard deviation

print('sample mean mu_hat                  : %f' % mu_hat)
print('sample standard deviation sigma_hat : %f' % sigma_hat)
l=scipy.stats.t.ppf( (1-gamma)/2, N-1)    # lower percentile
u=scipy.stats.t.ppf( 1-(1-gamma)/2, N-1)  # upper percentile
print('confidence interval mu_hat          : (%f, %f)' % 
      (mu_hat+l*sigma_hat/sqrt(N), mu_hat+u*sigma_hat/sqrt(N)))
l=scipy.stats.chi2.ppf( (1-gamma)/2, N-1)    # lower percentile
u=scipy.stats.chi2.ppf( 1-(1-gamma)/2, N-1)  # upper percentile
print('confidence interval sigma_hat       : (%f, %f)' % 
      ( sqrt((N-1)/u)*sigma_hat, sqrt((N-1)/l)*sigma_hat))

Leave a Reply

Your email address will not be published. Required fields are marked *