<jupyter_start><jupyter_text>Incremental Update of The Mean and Covariance of Running Samples Problem StatementWe have samples of $\mathbf{X}=\{\mathbf{x}_i\}_{i=1}^M$ with $M$ degree of freedom (DoF) with $M$ being a finite positive integer. For a sample in a given time, the state vector is $X_n = [x_{1,n}, x_{2,n}, ..., x_{M,n}]$ with $n$ being the sequence number of the sample. The mean of $N$ such samples is $\overline{X}(N)$ (an M-DoF vector):$$\overline{X}(N) = \frac{1}{N} \Sigma_{n=1}^{N} X_n$$or, to write it in another way:$$\overline{X}(N) = [\bar{x}_1(N), \bar{x}_2(N), ..., \bar{x}_M(N)] = [\frac{1}{N}\Sigma_{n=1}^N x_{1,n}, \frac{1}{N}\Sigma_{n=1}^N x_{2,n}, ..., \frac{1}{N}\Sigma_{n=1}^N x_{M,n}]$$And the corresponding unbiased covariance $\sigma_{i,j}(N)$ is:$$\sigma_{i,j}(N) = \frac{1}{N-1} \Sigma_{n=1}^{N} [x_{i,n} - \bar{x}_i(N)] [x_{j,n} - \bar{x}_j(N)]$$For a long running sampling process in which $N$ can grow infinitely, we may not be able to keep all samples of $\mathbf{X}$ to calculate the current mean and covariance. But we can calculate them incrementally with O(1) space complexity. DerivationHere we derive the formulae to calculate $\overline{X}(N+1)$ and $\sigma_{i,j}(N+1)$ given the current sample $X_{N+1}$ and previous mean $\overline{X}(N)$ and covariance $\sigma_{i,j}(N)$The new mean is pretty straightforward:$$\begin{align*}\overline{X}(N+1) &= \frac{1}{N+1} \Sigma_{n=1}^{N+1} X_n \\&= \frac{1}{N+1} (\Sigma_{n=1}^{N} X_n + X_{N+1}) \\&= \frac{1}{N+1} (N\overline{X}(N) + X_{N+1})\end{align*}$$The new covariance is:$$\begin{align*}\sigma_{i,j}(N+1) &= \frac{1}{N} \Sigma_{n=1}^{N+1} [x_{i,n} - \bar{x}_i(N+1)] [x_{j,n} - \bar{x}_j(N+1)] \\&= \frac{1}{N} (\Sigma_{n=1}^{N} [x_{i,n} - \bar{x}_i(N+1)] [x_{j,n} - \bar{x}_j(N+1)] + [x_{i,N+1} - \bar{x}_i(N+1)] [x_{j,N+1} - \bar{x}_j(N+1)])\end{align*}$$Since new term near the end of the equation above $[x_{i,N+1} - \bar{x}_i(N+1)] [x_{j,N+1} - \bar{x}_j(N+1)]$ can be easily calculated at this point, we will focus on getting the first term in the brackets:Let's define the mean update (which can also be easily calculated at this point):$$\Delta\bar{x}_i (N) = \bar{x}_i (N+1) - \bar{x}_i (N)$$Then:$$\begin{align*}&\Sigma_{n=1}^{N} [x_{i,n} - \bar{x}_i(N+1)] [x_{j,n} - \bar{x}_j(N+1)] \\&= \Sigma_{n=1}^{N} [x_{i,n} - (\bar{x}_i(N) + \Delta\bar{x}_i(N))] [x_{j,n} - (\bar{x}_j(N) + \Delta\bar{x}_j(N))] \\&= \Sigma_{n=1}^{N} [x_{i,n} - \bar{x}_i(N)][x_{j,n} - \bar{x}_j(N)] + \Sigma_{n=1}^{N} [x_{i,n} - \bar{x}_i (N)] \Delta\bar{x}_j(N) + \Sigma_{n=1}^{N} [x_{j,n} - \bar{x}_j (N)] \Delta\bar{x}_i(N)  + \Sigma_{n=1}^{N} \Delta\bar{x}_i(N) \Delta\bar{x}_j(N) \\&= \Sigma_{n=1}^{N} [x_{i,n} - \bar{x}_i(N)][x_{j,n} - \bar{x}_j(N)] + \Sigma_{n=1}^{N} \Delta\bar{x}_i(N) \Delta\bar{x}_j(N)  \\&= (N-1) \sigma_{i,j}(N) + N \Delta\bar{x}_i(N) \Delta\bar{x}_j(N)\end{align*}$$Therefore, the new covariance is:$$\sigma_{i,j}(N+1) = \frac{N-1}{N} \sigma_{i,j}(N) + [\bar{x}_i(N+1) - \bar{x}_i(N)] [(\bar{x}_j(N+1) - \bar{x}_j(N)] + \frac{1}{N} [x_{i,N+1} - \bar{x}_i(N+1)] [x_{j,N+1} - \bar{x}_j(N+1)]$$ Example Code<jupyter_code>import numpy as np

# Generate random samples
N_DoF = 9    # degree of freedom of states
N_sample = 10000
samples = np.random.default_rng().random((N_DoF, N_sample))

# Method1: use numpy built-in functions
np_mean = np.sum(samples, axis=1) / samples.shape[1]
np_cov = np.cov(samples)
print(f'np_mean: {np_mean}')
print(f'np_cov:\n{np_cov}')

# Method 2: use incremental update above
curr_mean = np.zeros(N_DoF)
curr_cov = np.zeros((N_DoF,N_DoF))
curr_N = 0    # number of samples till the time of calculation
for sample in samples.transpose():
    if curr_N == 0:
        curr_mean += sample
    else:
        next_mean = (curr_N * curr_mean + sample) / (curr_N + 1)
        delta_mean = next_mean - curr_mean
        for i in range(9):
            for j in range(9):
                curr_cov[i][j] = \
                        (curr_N - 1)/curr_N * curr_cov[i][j]    \
                        + delta_mean[i] * delta_mean[j] \
                        + 1/ curr_N * (sample[i] - next_mean[i]) * (sample[j] - next_mean[j])
        curr_mean = next_mean
    curr_N += 1
print(f'i_mean: {curr_mean}')
print(f'i_cov:\n{curr_cov}')<jupyter_output>np_mean: [0.49263978 0.4983813  0.50171179 0.49758825 0.49862414 0.49854833
 0.50209534 0.50408652 0.49486627]
np_cov:
[[ 8.47182198e-02  6.70486401e-04  2.34839021e-04 -1.49486481e-03
   1.76914667e-04  7.15548746e-04  3.49897398e-04  2.28240254e-04
  -8.87970718e-04]
 [ 6.70486401e-04  8.22198306e-02 -2.02858671e-04  3.23969968e-04
  -1.36810583e-03  1.01783075e-03  8.24766713e-04  3.03025339e-04
  -6.15513653e-04]
 [ 2.34839021e-04 -2.02858671e-04  8.34315486e-02 -3.96289016e-04
   1.24640315e-03 -1.14439229e-04  6.95774527e-04 -5.14504935e-05
   6.61015953e-04]
 [-1.49486481e-03  3.23969968e-04 -3.96289016e-04  8.32135561e-02
   8.04660923e-04 -8.17886412e-04  7.44674522e-04  9.17769318e-04
  -1.84914370e-03]
 [ 1.76914667e-04 -1.36810583e-03  1.24640315e-03  8.04660923e-04
   8.26705046e-02 -1.29864126e-03  1.53125393e-04 -1.24475335e-03
  -7.49379362e-04]
 [ 7.15548746e-04  1.01783075e-03 -1.14439229e-04 -8.17886412e-04
  -1.29864126e-03  8.36152597e-02 -6.24623760e-05 -2.1975491[...]
{'fraction_non_alphanumeric': 0.2299090219863533, 'fraction_numerical': 0.14234268385140258, 'mean_word_length': 3.788566243194192, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '357944', 'license': 'MIT', 'n_tokens_mistral': 2841, 'n_tokens_neox': 2435, 'n_words': 456, 'path': 'notes/incremental-update-mean-cov.ipynb', 'repo_name': 'lijuno/math-notes', 'size': '5276'}