Each Answer to this Q is separated by one/two green lines.

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?

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.

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.

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.*