0
$\begingroup$

I am trying to calculate the variance of the product of matrix A and vector b. As it was mathematically discussed on another post https://math.stackexchange.com/questions/2365166/what-is-the-variance-of-a-constant-matrix-times-a-random-vector, I have tried to do the calculation in python.

To recap, we have \begin{equation} Y = \bf{Ab} \end{equation} where \begin{equation} \bf{b} \sim N(\mu_b, \bf{I}\sigma_b^2) \end{equation}

Therefore \begin{equation} Var(Y) = \bf{A}\sigma_b^2\bf{A}^T \end{equation} I have tried this code with

import numpy as np
import statistics

A = np.arange(12).reshape(4, 3)   # some 4 by 3 array
b = np.array([2,4,5])              # some 3 by 1 vector

Y = np.dot(A,b)
statistics.variance(Y)
>> 1815

Vb = np.dot(np.identity(b.shape[0]), statistics.variance(b))
np.dot(A.dot(Vb), np.transpose(A))
>> array([[ 10.,  28.,  46.,  64.],
       [ 28., 100., 172., 244.],
       [ 46., 172., 298., 424.],
       [ 64., 244., 424., 604.]])

However, the statistic.variance(Y) outputs a different result compared to np.dot(A.dot(Vb), np.transpose(A)).

I wonder how should I do it properly in python?

Your help is really appreciated!

$\endgroup$
0
$\begingroup$

To test formulas like these numerically, you need to compute sample statistics, like this:

    import numpy as np
    A = np.arange(12).reshape(4, 3)  # Note: A is a fixed matrix
    b_avg = np.array([2, 4, 5])  # b is random, with a mean vector
    b_var = 1e0 * np.eye(b_avg.size)  # and a (co)variance matrix

    # generate a random sample of b vectors    
    n_sample = 10**4
    b_sample = np.random.multivariate_normal(
        mean=b_avg, cov=b_var, size=n_sample).T

    # compute sample moments (each column is a sample)
    b_avg_smp = b_sample.mean(axis=1)
    b_var_smp = np.cov(b_sample)

    # compare sample moments to population moments
    print('b vector: statistics')
    print('Averages (population vs. sample):',
          b_avg.tolist(), b_avg_smp.round(1).tolist(), sep='\n')
    print('Variances (population vs. sample):',
          b_var.tolist(), b_var_smp.round(1).tolist(), sep='\n')

This will give something like this:

b vector: statistics
Averages (population vs. sample):
[2, 4, 5]
[2.0, 4.0, 5.0]
Variances (population vs. sample):
[[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]
[[1.0, -0.0, 0.0], [-0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]

i.e. the sample statistics are consistent with the specified $b$ distribution.

Then we can transform the $b$ sample into a $y$ sample, and do the analogous comparison, using your analytically computed $y$ moments, i.e.

    # transform b moments into y moments (@ is matrix multiply)
    y_avg = A @ b_avg
    y_var = A @ b_var @ A.T

    # transform b sample into y sample, and compute sample moments
    y_sample = A @ b_sample
    y_avg_smp = y_sample.mean(axis=1)
    y_var_smp = np.cov(y_sample)
    
    # compare population vs. sample moments
    print('y vector: statistics')
    print('Averages (population vs. sample):',
          y_avg.tolist(), y_avg_smp.round(1).tolist(), sep='\n')
    print('Variances (population vs. sample):',
          y_var.tolist(), y_var_smp.round(1).tolist(), sep='\n')

This will give something like this:

y vector: statistics
Averages (population vs. sample):
[14, 47, 80, 113]
[13.9, 46.8, 79.7, 112.7]
Variances (population vs. sample):
[[5.0, 14.0, 23.0, 32.0], [14.0, 50.0, 86.0, 122.0], [23.0, 86.0, 149.0, 212.0], [32.0, 122.0, 212.0, 302.0]]
[[5.0, 14.0, 23.1, 32.1], [14.0, 50.1, 86.1, 122.1], [23.1, 86.1, 149.0, 212.0], [32.1, 122.1, 212.0, 302.0]]

i.e. the sample $y$ statistics are consistent with the analytically computed $y$ distribution.

(Alternatively, the uncertainties package will automatically do linearized error propagation for Gaussian random variables.)

$\endgroup$
1
  • $\begingroup$ Thank you very much, @GeoMatt22, I got the idea $\endgroup$ – user2842390 Jan 9 at 2:29
0
$\begingroup$

Your are mixing some concepts. It's hard to explain but i will try.

statistics.variance(b) gives you the variance of n samples in variable b in 1 dimension. The variance matrix for b is a 1x1 matrix.

Vb is the variance matrix for 3 dimensional sample with several samples and zero covarinace between the dimensions. But it is not the variance matrix for b.

The definition of A for b is wrong. It works because np.dot uses automatically b.T, but conceptually is wrong.

I hope this help to understand it.

$\endgroup$
1
  • $\begingroup$ Thanks. Would you suggest any simulations as I showed in the code above to demonstrate how to calculate the variance of Y as in the scenario with Y=Ab where b ~ $N(\mu,\bf{I\sigma_b^2})$ ? Thank you! $\endgroup$ – user2842390 Jan 8 at 16:38

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.