# How to calculate the confidence interval of the mean of means? (2)

Note: this question is different than this other one as follows:

• That one asks for the CI mean of means.
• This one asks for the CI mean of means, only from the perspective of the variations within the small observations from each system, and not from the overall variations between the evaluated systems.

For example, say that we got these sample arrays:

• $a = [56,0,57]$.
• $b = [227,265,263]$.
• $c = [448,406,488]$.

Where:

• Array $a$ has 3 samples from process A.
• Array $b$ has 3 samples from process B.
• Array $c$ has 3 samples from process C.

And their sample statistics (mean $\bar x$, variance $\sigma_x^2$, stdev $\sigma_x$, standard error $\text{SE}_{\bar x}$) as I've calculated from the arrays above are:

• $\bar a=37.667$, $\sigma^2_a=1064.333$, $\sigma_a=32.6241$, $\text{SE}_{\bar a}=\frac{\sigma_a}{\sqrt{3}} = 18.8355$.
• $\bar b=251.667$, $\sigma^2_b=457.333$, $\sigma_b=21.3854$, $\text{SE}_{\bar b}=\frac{\sigma_b}{\sqrt{3}} = 12.3468$.
• $\bar c=447.333$, $\sigma^2_c=1681.333$, $\sigma_c=41.004$, $\text{SE}_{\bar c}=\frac{\sigma_c}{\sqrt{3}} = 23.6737$.

Now, I calculate the $0.95$ confidence intervals of each array as follows:

$$\text{CI}_{\bar a} = 4.303 \times \text{SE}_{\bar a} = \pm 81.0428$$

$$\text{CI}_{\bar b} = 4.303 \times \text{SE}_{\bar b} = \pm 53.1242$$

$$\text{CI}_{\bar c} = 4.303 \times \text{SE}_{\bar c} = \pm 101.8597$$

So far all seem straightforward.

## The question:

Now suppose that I calculate the mean of the means. Let's call it $\overline {abc} = \frac{\bar a + \bar b + \bar c}{3} = 245.556$

• What are the confidence intervals of the mean of means $\overline{abc}$? Let's call it $\text{CI}_{\overline{abc}}$?

I don't want $\text{CI}_{\overline{abc}}$ to account for variations between systems A, B, and C. I only want it to account for variations within them.

So, for example, the means of the systems A, B and C could vary drastically (because they are very different systems). However, the samples within each mean might not vary drastically (thus low $\text{CI}_{\bar x}$ per each system). I therefore expect $\text{CI}_{\overline{abc}}$ to be small because the $\text{CI}_{\bar x}$ is small (for any $x$).

## My attempt:

Python code:

import scipy
import scipy.stats

# data
a = [56, 0, 57]
b = [227, 265, 263]
c = [448, 406, 488]

# how many SE needed for a CI of .95
h = scipy.stats.t._ppf((1+0.95)/2.0, len(a)-1)

# compute many means
means = []
for i in range(len(a)):
mean = scipy.mean(a[i] + b[i] + c[i])
means.append(mean)

# find the CI of those many means
means_se = scipy.stats.sem(means)
means_ci = h * means_se
print "mean of means CI is: +/- " + str(means_ci)


Output:

mean of means CI is: +/- 170.599562802


Note that this mean of means CI that only accounts for intra-process variation, is different then the mean of means CI that only accounts for inter-process variation. Here is the other CI (note that it's much bigger):

import scipy
import scipy.stats

# data
a = [56, 0, 57]
b = [227, 265, 263]
c = [448, 406, 488]

# how many SE needed for a CI of .95
h = scipy.stats.t._ppf((1+0.95)/2.0, 2)

mean_a = scipy.mean(a)
mean_b = scipy.mean(b)
mean_c = scipy.mean(c)

means_se = scipy.stats.sem([mean_a,mean_b,mean_c])
means_ci = h * means_se
print "other mean of means CI is: +/- " + str(means_ci)


Output:

other mean of means CI is: +/- 509.004022253


I know two approaches. One is uncertainty propagation. For python linear uncertainty propagation is implemented in the uncertainties package:

import uncertainties

ma = uncertainties.ufloat(scipy.mean(a), scipy.stats.sem(a))
mb = uncertainties.ufloat(scipy.mean(b), scipy.stats.sem(b))
mc = uncertainties.ufloat(scipy.mean(c), scipy.stats.sem(c))

print scipy.mean([ma, mb, mc])


Which gives 246 +- 11.

The other approach is bootstrapping which is maybe similar to your first approach. A case resampling scheme could for instance be implemented like this:

import scipy
import scipy.stats

# data
a = [56, 0, 57]
b = [227, 265, 263]
c = [448, 406, 488]

N = 1000
means = []
for i in range(N):
mean = scipy.mean(scipy.mean(scipy.random.choice(a, len(a))) +\
scipy.mean(scipy.random.choice(b, len(b))) +\
scipy.mean(scipy.random.choice(c, len(c)))
)
means.append(mean)
print scipy.mean(means), scipy.std(means)


Which gives something like 245 +-9.