How to calculate the variance of a partition of variables

Question

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.

Answer 1

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.

Answer 2

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))

Answer 3

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