# How to calculate the variance of a partition of variables

I'm running an experiment where I'm gathering (independent) samples in parallel, I compute the variance of each group of samples and now I want to combine then all to find the total variance of all the samples.

I'm having a hard time finding a derivation for this as I'm not sure of terminology. I think of it as a partition of one RV.

So I want to find $Var(X)$ from $Var(X_1)$, $Var(X_2)$, ..., and $Var(X_n)$, where $X$ = $[X_1, X_2, \dots, X_n]$.

EDIT: The partitions are not the same size/cardinality, but the sum of the partition sizes equal the number of samples in the overall sample set.

EDIT 2: There is a formula for a parallel computation here, but it only covers the case of a partition into two sets, not $n$ sets.

• Is this the same as my question here: mathoverflow.net/questions/64120/…
– user4500
Commented May 6, 2011 at 20:29
• What does that last bracket mean? And what do you mean by "total variance"? Is it anything other than the variance of the combined dataset?
– whuber
Commented May 6, 2011 at 20:38
• @whuber which last bracket? "total variance" means the variance of the total dataset. Commented May 9, 2011 at 12:37
• The expression $[X_1, X_2, \dots, X_n]$ could mean many things (although conventionally it would be a vector): I was looking for a clarification.
– whuber
Commented May 9, 2011 at 13:59

The formula is fairly straightforward if all the sub-sample have the same sample size. If you had $g$ sub-samples of size $k$ (for a total of $gk$ samples), then the variance of the combined sample depends on the mean $E_j$ and variance $V_j$ of each sub-sample: $$Var(X_1,\ldots,X_{gk}) = \frac{k-1}{gk-1}(\sum_{j=1}^g V_j + \frac{k(g-1)}{k-1} Var(E_j)),$$ where by $Var(E_j)$ means the variance of the sample means.

A demonstration in R:

> x <- rnorm(100)
> g <- gl(10,10)
> mns <- tapply(x, g, mean)
> vs <- tapply(x, g, var)
> 9/99*(sum(vs) + 10*var(mns))
[1] 1.033749
> var(x)
[1] 1.033749


If the sample sizes are not equal, the formula is not so nice.

EDIT: formula for unequal sample sizes

If there are $g$ sub-samples, each with $k_j, j=1,\ldots,g$ elements for a total of $n=\sum{k_j}$ values, then $$Var(X_1,\ldots,X_{n}) = \frac{1}{n-1}\left(\sum_{j=1}^g (k_j-1) V_j + \sum_{j=1}^g k_j (\bar{X}_j - \bar{X})^2\right),$$ where $\bar{X} = (\sum_{j=1}^gk_j\bar{X}_j)/n$ is the weighted average of all the means (and equals to the mean of all values).

Again, a demonstration:

> k <- rpois(10, lambda=10)
> n <- sum(k)
> g <- factor(rep(1:10, k))
> x <- rnorm(n)
> mns <- tapply(x, g, mean)
> vs <- tapply(x, g, var)
> 1/(n-1)*(sum((k-1)*vs) + sum(k*(mns-weighted.mean(mns,k))^2))
[1] 1.108966
> var(x)
[1] 1.108966


By the way, these formulas are easy to derive by writing the desired variance as the scaled sum of $(X_{ji}-\bar{X})^2$, then introducing $\bar{X}_j$: $[(X_{ji}-\bar{X}_j)-(\bar{X}_j-\bar{X})]^2$, using the square of difference formula, and simplifying.

• thanks. Unfortunately, I can't guarantee that my partitions are all the same size. I'm running a massively parallel process where I need to calculate the variances of each partition in parallel then combine in the end, but the results/samples from each parallel process are not equal (it's a Monte Carlo simulation of received photons). Commented May 9, 2011 at 12:43
• I can't +1 this enough, super helpful formula for parallel computation in a data warehouse environment Commented Mar 8, 2013 at 17:31

This is simply an add-on to the answer of aniko with a rough sketch of the derivation and some python code, so all credits go to aniko.

# derivation

Let $$X_j \in X = \{X_1, X_2, \ldots, X_g\}$$ be one of $$g$$ parts of the data where the number of elements in each part is $$k_j = |X_j|$$. We define the mean and the variance of each part to be \begin{align*} E_j & = \mathrm{E}\left[X_j\right] = \frac{1}{k_j} \sum_{i=1}^{k_j} X_{ji}\\ V_j & = \mathrm{Var}\left[X_j\right] = \frac{1}{k_j-1} \sum_{i=1}^{k_j} (X_{ji} - E_j)^2 \end{align*} respectively. If we set $$n = \sum_{j=1}^g k_j$$, the variance of the total dataset is given by: \begin{align*} \mathrm{Var}\left[X\right] & = \frac{1}{n-1} \sum_{j=1}^{g} \sum_{i=1}^{k_j} (X_{ji} - \mathrm{E}\left[X\right])^2 \\ & = \frac{1}{n-1} \sum_{j=1}^{g} \sum_{i=1}^{k_j} \big((X_{ji} - E_j) - (\mathrm{E}\left[X\right] - E_j)\big)^2 \\ & = \frac{1}{n-1} \sum_{j=1}^{g} \sum_{i=1}^{k_j} (X_{ji} - E_j)^2 - 2(X_{ji} - E_j)(\mathrm{E}\left[X\right] - E_j) + (\mathrm{E}\left[X\right] - E_j)^2 \\ & = \frac{1}{n-1} \sum_{j=1}^{g} (k_j - 1) V_j + k_j (\mathrm{E}\left[X\right] - E_j)^2. \end{align*} If we have the same size $$k$$ for each part, i.e. $$\forall j: k_j = k$$, above formula simplifies to \begin{align*} \mathrm{Var}\left[X\right] & = \frac{1}{n-1} \sum_{j=1}^g (k-1) V_j + k(g-1) \mathrm{Var}\left[E_j\right] \\ & = \frac{k-1}{n-1} \sum_{j=1}^g V_j + \frac{k(g-1)}{k-1} \mathrm{Var}\left[E_j\right] \end{align*}

# python code

The following python function works for arrays that have been splitted along the first dimension and implements the "more complex" formula for differently sized parts.

import numpy as np

def combine(averages, variances, counts, size=None):
"""
Combine averages and variances to one single average and variance.

# Arguments
averages: List of averages for each part.
variances: List of variances for each part.
counts: List of number of elements in each part.
size: Total number of elements in all of the parts.
# Returns
average: Average over all parts.
variance: Variance over all parts.
"""
average = np.average(averages, weights=counts)

# necessary for correct variance in case of multidimensional arrays
if size is not None:
counts = counts * size // np.sum(counts, dtype='int')

squares = (counts - 1) * variances + counts * (averages - average)**2
return average, np.sum(squares) / (size - 1)


It can be used as follows:

# sizes k_j and n
ks = np.random.poisson(10, 10)
n = np.sum(ks)

# create data
x = np.random.randn(n, 20)
parts = np.split(x, np.cumsum(ks[:-1]))

# compute statistics on parts
ms = [np.mean(p) for p in parts]
vs = [np.var(p, ddof=1) for p in parts]

# combine and compare
combined = combine(ms, vs, ks, x.size)
numpied = np.mean(x), np.var(x, ddof=1)
distance = np.abs(np.array(combined) - np.array(numpied))
print('combined --- mean:{: .9f} - var:{: .9f}'.format(*combined))
print('numpied  --- mean:{: .9f} - var:{: .9f}'.format(*numpied))
print('distance --- mean:{: .5e} - var:{: .5e}'.format(*distance))

• I believe you need a pair of brackets in your final equation (same size k) after you extract a factor of $k-1$. Commented Sep 29, 2020 at 12:08
• @DzamoNorton I am afraid I do not understand where you want additional brackets Commented Sep 30, 2020 at 5:43
• I've illustrated it here mathb.in/45578. You can also see said brackets in the same-size-k formula in the other answer to this question. Commented Sep 30, 2020 at 9:19
• But I also see now I'd lost track of the scope of the summand. Now that I see the summand must run to the end of the line, the brackets are indeed unnecessary. Apologies! Commented Sep 30, 2020 at 9:27

Many thanks to those who answered and commented! Your answers helped me with a problem where I wanted long-term retention of descriptive statistics, but was not able to retain the original data.

My solution was to create rollup records that would have one row for each combination of desired dimensions of the data (one of which might be time). Depending on the cardinality of the dimensions, this can be a huge space savings. In this example, we have 100 rollup records for 1 million measurements, since the cardinality of each dimension is 10.

I thought I'd demonstrate this method in SQL (PostgreSQL), and perhaps others might find this useful. This is my first answer, so do be kind. I'm not sampling the data, so I use the population formulas that use the divisor of N.

First, we create some measurement records, then roll them up into summary rollup records. All we store is sum, count, and population variance for each combination.

stats=> \timing
Timing is on.
stats=> SELECT setseed(.5);
setseed
---------

(1 row)

Time: 0.580 ms
stats=> drop table if exists measurements ;
NOTICE:  table "measurements" does not exist, skipping
DROP TABLE
Time: 0.450 ms
stats=> SELECT setseed(.5);
setseed
---------

(1 row)

Time: 0.497 ms
stats=> create table measurements as
SELECT generate_series(1,1000000) AS id,
floor(random() *10) AS dimension1,
floor(random() *10) AS dimension2,
random() * 10 AS value;
SELECT 1000000
Time: 776.785 ms
stats=> drop table if exists rollups;
NOTICE:  table "rollups" does not exist, skipping
DROP TABLE
Time: 0.839 ms
stats=> create table rollups as
select dimension1, dimension2, count(1) as measurements,
sum(value) as value_sum, var_pop(value) as value_var_pop
from measurements
group by dimension1, dimension2;
SELECT 100
Time: 880.328 ms


Then we calculate the mean and population standard deviation, first from the rollup records, then the complete dataset. This uses a "Subquery" to calculate the population mean used in the calculation. Note that the query times from summary records are considerably faster and should be acceptable for real-time display.

stats=> SELECT
SUM (value_sum) / SUM (measurements) mean,
SQRT (
(
SUM (
ru.value_var_pop * ru.measurements
) + SUM (
ru.measurements * (
(ru.value_sum / ru.measurements) /* row avg */
- (                              /* grand avg subquery */
SELECT
SUM (value_sum) / SUM (measurements) grand_avg
FROM
rollups ga
)
) ^ 2
)
) / SUM (ru.measurements)
) as stddev_pop
FROM
rollups ru;
mean       |    stddev_pop
------------------+------------------
4.99954057566761 | 2.88618225050068
(1 row)

Time: 1.255 ms

stats=> select avg(value), stddev_pop(value) from measurements;
avg        |    stddev_pop
------------------+------------------
4.99954057566761 | 2.88618225050069
(1 row)

Time: 71.719 ms


As a bonus, we calculate mean and standard deviations by "dimension1," first from the rollup records, then the complete dataset. This uses a "Correlated Subquery" to calculate the means for each dimension1 value.

stats=> SELECT dimension1,
SUM (value_sum) / SUM (measurements) mean,
SQRT (
(
SUM (
ru.value_var_pop * ru.measurements
) + SUM (
ru.measurements * (
(ru.value_sum / ru.measurements) /* row avg */
- (                              /* grand avg subquery */
SELECT
SUM (value_sum) / SUM (measurements) grand_avg
FROM
rollups ga
WHERE ga.dimension1 = ru.dimension1   /* same (correlated) group as outer query  */
)
) ^ 2
)
) / SUM (ru.measurements)
) as stddev_pop
FROM
rollups ru
GROUP BY
dimension1
order by 1 ;
dimension1 |       mean       |    stddev_pop
------------+------------------+------------------
0 | 5.00667407841466 | 2.89025062220051
1 | 4.98845700180817 | 2.88379096643363
2 | 4.99899669999969 | 2.88825163824503
3 | 4.98553153843827 | 2.88365611510306
4 | 4.99844263698516 |  2.8860230080168
5 |  5.0179497778027 | 2.88749365173117
6 | 5.00517112927829 | 2.88346174527641
7 | 5.00532821601654 | 2.89266358278946
8 | 4.98688468436205 | 2.88383465545076
9 | 5.00219626955073 | 2.88222690226276
(10 rows)

Time: 3.828 ms
stats=> select dimension1, avg(value), stddev_pop(value) from measurements group by dimension1 order by 1;
dimension1 |       avg        |    stddev_pop
------------+------------------+------------------
0 | 5.00667407841466 | 2.89025062220053
1 | 4.98845700180817 | 2.88379096643364
2 | 4.99899669999969 | 2.88825163824505
3 | 4.98553153843827 | 2.88365611510306
4 | 4.99844263698516 | 2.88602300801678
5 |  5.0179497778027 | 2.88749365173116
6 | 5.00517112927829 |  2.8834617452764
7 | 5.00532821601654 | 2.89266358278946
8 | 4.98688468436205 | 2.88383465545078
9 | 5.00219626955073 | 2.88222690226277
(10 rows)

Time: 113.342 ms