then each can be expressed in terms of the equivalent raw moments: where is generally taken to be the duration of the time-history, or the number of points if is constant. The benefit of expressing the statistical moments in terms of is that the sets can be combined by addition, and there is no upper limit on the value of . where the subscript represents the concatenated time-history or combined . These combined values of can then be inversely transformed into raw moments representing the complete concatenated time-history Known relationships between the raw moments () and the central moments () are then used to compute the central moments of the concatenated time-history. Finally, the statistical moments of the concatenated history are computed from the central moments: Covariance Very similar algorithms can be used to compute the covariance. Naïve algorithm The naïve algorithm is For the algorithm above, one could use the following Python code: def naive_covariance(data1, data2): n = len(data1) sum1 = sum(data1) sum2 = sum(data2) sum12 = sum([i1 * i2 for i1, i2 in zip(data1, data2)]) covariance = (sum12 - sum1 * sum2 / n) / n return covariance With estimate of the mean As for the variance, the covariance of two random variables is also shift-invariant, so given any two constant values and it can be written: and again choosing a value inside the range of values will stabilize the formula against catastrophic cancellation as well as make it more robust against big sums. Taking the first value of each data set, the algorithm can be written as: def shifted_data_covariance(data_x, data_y): n = len(data_x) if n < 2: return 0 kx = data_x[0] ky = data_y[0] Ex = Ey = Exy = 0 for ix, iy in zip(data_x, data_y): Ex += ix - kx Ey += iy - ky Exy += (ix - kx) * (iy - ky) return (Exy - Ex * Ey / n) / n Two-pass The two-pass algorithm first computes the sample means, and then the covariance: The two-pass algorithm may be written as: def two_pass_covariance(data1, data2): n = len(data1) mean1 = sum(data1) / n mean2 = sum(data2) / n covariance = 0 for i1, i2 in zip(data1, data2): a = i1 - mean1 b = i2 - mean2 covariance += a * b / n return covariance A slightly more accurate compensated version performs the full naive algorithm on the residuals. The final