Variance of product of a matrix and vector using Python 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)
&gt;&gt; 1815

Vb = np.dot(np.identity(b.shape[0]), statistics.variance(b))
np.dot(A.dot(Vb), np.transpose(A))
&gt;&gt; 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!
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.