Correct way to obtain confidence interval with scipy

I have a 1-dimensional array of data:

a = np.array([1,2,3,4,4,4,5,5,5,5,4,4,4,6,7,8])

for which I want to obtain the 68% confidence interval (ie: the 1 sigma).

The first comment in this answer states that this can be achieved using scipy.stats.norm.interval from the scipy.stats.norm function, via:

from scipy import stats
import numpy as np
mean, sigma = np.mean(a), np.std(a)

conf_int = stats.norm.interval(0.68, loc=mean, 
    scale=sigma)

But a comment in this post states that the actual correct way of obtaining the confidence interval is:

conf_int = stats.norm.interval(0.68, loc=mean, 
    scale=sigma / np.sqrt(len(a)))

that is, sigma is divided by the square-root of the sample size: np.sqrt(len(a)).

The question is: which version is the correct one?

Solution 1:

The 68% confidence interval for a single draw from a normal distribution with mean mu and std deviation sigma is

stats.norm.interval(0.68, loc=mu, scale=sigma)

The 68% confidence interval for the mean of N draws from a normal distribution with mean mu and std deviation sigma is

stats.norm.interval(0.68, loc=mu, scale=sigma/sqrt(N))

Intuitively, these formulas make sense, since if you hold up a jar of jelly beans and ask a large number of people to guess the number of jelly beans, each individual may be off by a lot -- the same std deviation sigma -- but the average of the guesses will do a remarkably fine job of estimating the actual number and this is reflected by the standard deviation of the mean shrinking by a factor of 1/sqrt(N).

If a single draw has variance sigma**2, then by the Bienaymé formula, the sum of N uncorrelated draws has variance N*sigma**2.

The mean is equal to the sum divided by N. When you multiply a random variable (like the sum) by a constant, the variance is multiplied by the constant squared. That is

Var(cX) = c**2 * Var(X)

So the variance of the mean equals

(variance of the sum)/N**2 = N * sigma**2 / N**2 = sigma**2 / N

and so the standard deviation of the mean (which is the square root of the variance) equals

sigma/sqrt(N).

This is the origin of the sqrt(N) in the denominator.

Here is some example code, based on Tom's code, which demonstrates the claims made above:

import numpy as np
from scipy import stats

N = 10000
a = np.random.normal(0, 1, N)
mean, sigma = a.mean(), a.std(ddof=1)
conf_int_a = stats.norm.interval(0.68, loc=mean, scale=sigma)

print('{:0.2%} of the single draws are in conf_int_a'
      .format(((a >= conf_int_a[0]) & (a < conf_int_a[1])).sum() / float(N)))

M = 1000
b = np.random.normal(0, 1, (N, M)).mean(axis=1)
conf_int_b = stats.norm.interval(0.68, loc=0, scale=1 / np.sqrt(M))
print('{:0.2%} of the means are in conf_int_b'
      .format(((b >= conf_int_b[0]) & (b < conf_int_b[1])).sum() / float(N)))

prints

68.03% of the single draws are in conf_int_a
67.78% of the means are in conf_int_b

Beware that if you define conf_int_b with the estimates for mean and sigma based on the sample a, the mean may not fall in conf_int_b with the desired frequency.

If you take a sample from a distribution and compute the sample mean and std deviation,

mean, sigma = a.mean(), a.std()

be careful to note that there is no guarantee that these will equal the population mean and standard deviation and that we are assuming the population is normally distributed -- those are not automatic givens!

If you take a sample and want to estimate the population mean and standard deviation, you should use

mean, sigma = a.mean(), a.std(ddof=1)

since this value for sigma is the unbiased estimator for the population standard deviation.

Solution 2:

I just checked how R and GraphPad calculate confidence intervals, and they increase the interval in case of small sample size (n). E.g., more than 6-fold for n=2 compared to a large n. This code (based on shasan's answer) matches their confidence intervals:

import numpy as np, scipy.stats as st

# returns confidence interval of mean
def confIntMean(a, conf=0.95):
  mean, sem, m = np.mean(a), st.sem(a), st.t.ppf((1+conf)/2., len(a)-1)
  return mean - m*sem, mean + m*sem

For R, I checked against t.test(a). GraphPad's confidence interval of a mean page has "user level" info on the sample size dependency.

Here the output for Gabriel's example:

In [2]: a = np.array([1,2,3,4,4,4,5,5,5,5,4,4,4,6,7,8])

In [3]: confIntMean(a, 0.68)
Out[3]: (3.9974214366806184, 4.877578563319382)

In [4]: st.norm.interval(0.68, loc=np.mean(a), scale=st.sem(a))
Out[4]: (4.0120010966037407, 4.8629989033962593)

Note that the difference between the confIntMean() and st.norm.interval() intervals is relatively small here; len(a) == 16 is not too small.

Solution 3:

I tested out your methods using an array with a known confidence interval. numpy.random.normal(mu,std,size) returns an array centered on mu with a standard deviation of std (in the docs, this is defined as Standard deviation (spread or “width”) of the distribution.).

from scipy import stats
import numpy as np
from numpy import random
a = random.normal(0,1,10000)
mean, sigma = np.mean(a), np.std(a)
conf_int_a = stats.norm.interval(0.68, loc=mean, scale=sigma)
conf_int_b = stats.norm.interval(0.68, loc=mean, scale=sigma / np.sqrt(len(a)))


conf_int_a
(-1.0011149125527312, 1.0059797764202412)
conf_int_b
(-0.0076030415111100983, 0.012467905378619625)

As the sigma value should be -1 to 1, the / np.sqrt(len(a)) method appears to be incorrect.

Edit

Since I don't have the reputation to comment above I'll clarify how this answer ties into unutbu's thorough answer. If you populate a random array with a normal distribution, 68% of the total will fall within 1-σ of the mean. In the case above, if you check that you see

b = a[np.where((a>-1)&(a <1))]
len(a)
> 6781

or 68% of the population falls within 1σ. Well, about 68%. As you use a larger and larger array, you will approach 68% (In a trial of 10, 9 were between -1 and 1). That's because the 1-σ is the inherent distribution of the data, and the more data you have the better you can resolve it.

Basically, my interpretation of your question was If I have a sample of data I want to use to describe the distribution they were drawn from, what is the method to find the standard deviation of that data? whereas unutbu's interpretation appears to be more What is the interval to which I can place the mean with 68% confidence?. Which would mean, for jelly beans, I answered How are they guessing and unutbu answered What do their guesses tell us about the jelly beans.