Standard Deviation addition and subtraction

Question

Is there a formula to calculate the Standard deviation biased on the sum or subtraction of others data sets ?

Example:

Dataset1 (5 elements to count):
values: 5,10,15,20,25
mean: 15
Sum of Squared mean: 275 (5^2+10^2+...)/5
Population variance: 50
Population Standard deviation: 7,071067812
Population Max STD  22,07106781
Population Min STD  7,928932188

Dataset2 (5 elements to count):
values: 2,4,11,7,16
mean: 8
Sum of Squared mean: 89,2 (2^2+4^2+...)/5
Population variance: 25,2
Population Standard deviation: 5,019960159
Population Max STD  13,01996016
Population Min STD  2,980039841

Dataset3 (5 elements to count):
The elements are a sum of the previous dataset
values: 7,14,26,27,41
mean: 23 (<-- Ok, sum of the previous means)
Sum of Squared mean: 666,2
Population variance: 137,2
Population Standard deviation: 11,71324037
Population Max STD  34,71324037
Population Min STD  11,28675963

The mean of Data3 is easily computed as Mean of Data1 + Mean of data2

But,...how to I calculate the their values?

For example, knowing that the Squared Sum can be used to calculate the variance. Is there a way to calculate directly the Squared Sum of Data3 using a formula biased on the Data1 and Data2?

If not, is there a way to calculate variance of Data3, without using covariance ? (This is because, covariance will assume I´ll have to perform another calculation of sums). I was thinking ins a formula more directly, instead of calculating each element all over again.

Answer 1

The Standard deviation, mean, variance, etc. can be calculated by simply keeping a running tally of some calculations. The sum of the sets can be added by also keeping a running sum of the product of the 2 data points. Ref

The calculation of STD is sensitive to subtraction. Keeping sumxx, sumxy, sumyy to higher precision is recommended. The below uses long, and long long. A FP implementation could use double and long double. Details of higher precision issues are not expressed deeply here other to to use them with the running summations.

#include <stdio.h>
#include <math.h>

struct stat2 {
  long sumx;
  long sumy;
  long long sumxx;
  long long sumxy;
  long long sumyy;
  size_t count;
};

void stat2_add(struct stat2 *stat, long x, long y) {
  stat->sumx += x;
  stat->sumxx += 1LL * x * x;
  stat->sumy += y;
  stat->sumyy += 1LL * y * y;

  // This is the only extra reoccurring work needed to meet OP's goal
  stat->sumxy += 1LL * x * y; 

  stat->count++;
}

double stat2_avg(const struct stat2 *stat, int index) {
  switch (index) {
    case 'x':
      return 1.0 * stat->sumx / stat->count;
    case 'y':
      return 1.0 * stat->sumy / stat->count;
    default:
      return 1.0 * (stat->sumx + stat->sumy) / stat->count;
  }
}

double stat2_std(const struct stat2 *stat, int index) {
  double offset = 0.0;  // or 1.0 depending on STD model
  double var;
  switch (index) {
    case 'x':
      var = (stat->sumxx - 1.0 * stat->sumx * stat->sumx / stat->count)
          / (stat->count - offset);
      break;
    case 'y':
      var = (stat->sumyy - 1.0 * stat->sumy * stat->sumy / stat->count)
          / (stat->count - offset);
      break;
    default: {
      // SUM(x+y) = SUM(x) + SUM(y)
      double z = stat->sumx + stat->sumy;
      // SUM((x+y)*(x+y)) = SUM(x*x) + 2*SUM(x*y) + SUM(y*y)
      double zz = stat->sumxx + 2LL * stat->sumxy + stat->sumyy;
      var = (zz - 1.0 * z * z / stat->count) / (stat->count - offset);
    }
  }
  return sqrt(var);
}

void stat2_report(const struct stat2 *stat, const char *title) {
  printf("%s\n", title);
  printf("  x   Avg:%9f  STD:%f\n", stat2_avg(stat, 'x'), stat2_std(stat, 'x'));
  printf("  y   Avg:%9f  STD:%f\n", stat2_avg(stat, 'y'), stat2_std(stat, 'y'));
  printf("  x+y Avg:%9f  STD:%f\n", stat2_avg(stat, 'z'), stat2_std(stat, 'z'));
}

int main(void) {
  size_t i;
  struct stat2 A = { 0, 0, 0, 0, 0 };
  int dataA[] = { 5, 10, 15, 20, 25 };
  int dataB[] = { 2, 4, 11, 7, 16 };
  for (i = 0; i < 5; i++)
    stat2_add(&A, dataA[i], dataB[i]);
  stat2_report(&A, "A");
  return 0;
}

Output

A
  x   Avg:15.000000  STD:7.071068
  y   Avg: 8.000000  STD:5.019960
  x+y Avg:23.000000  STD:11.713240