Multivariate Normal Random Variables are extraordinarily convenient.
The probability density of a multivariate normal random variable is proportional to the following:
Where is an n-dimensional state vector, is the mean of the distribution, and is an n by n covariance matrix. Pictorally a 2-D density might be represented like this:
With contour lines showing decreasing probability levels dropping off around the mean (blue x). This distribution is entirely defined by the two quantities, which gives the center of the distribution and which effectively gives the shape of the ellipses. That is, rather than carry around the functional form above, we can simply define X as and forget the rest.
Multivariate normals are convenient for three reasons
- They are easy to represent – we only need a mean and covariance matrix
- Linear transformations of normals are again normals
- All operations are represented with linear algebra
First off, multivariate normals are simple to represent. This ends up being a big deal for functions on very high dimensional spaces. Imagine writing down a general function on 1000 variables.
Second, linear functions of normals are again normals. This is huge. For example this means that we could project the image above to one of the coordinate axes (or any axis) and get out our old friend the bell curve. As we work on our random variables the three conveniences remain true.
Third, the computation to perform these linear transformations of random variables is done solely through linear algebra on the mean and covariance matrices. Fortunately, linear algebra is something about which we know quite a bit.
So, as long as we’re willing to say that our variables are normally distributed (which is often not far from the truth) we can efficiently represent and compute on huge spaces of interconnected variables.
Multivariate Normals (MVNs) have been a goal of mine for some time while working on this project. They’re where this project starts to intersect with my actual work. I do lots of manipulations on MVNs and would like to stop dealing with all the matrix algebra.
In order to build them correctly it was clear I would need a relatively powerful symbolic matrix expression system. I’ve been working on something over at this branch.
Now we can represent symbolic matrices and, using them, represent MVN Random Variables
# Lets make a Multivariate Normal Random Variable >>> mu = MatrixSymbol('mu', n, 1) # n by 1 mean vector >>> Sigma = MatrixSymbol('Sigma', n, n) # n by n covariance matrix >>> X = Normal(mu, Sigma, 'X') # a multivariate normal random variable # Density is represented just by the mean and covariance >>> Density(X) (μ, Σ) >>> H = MatrixSymbol('H', k, n) # A linear operator >>> Density(H*X) # What is the density of X after being transformed by H? (H⋅μ, H⋅Σ⋅H') # Lets make some measurement noise >>> zerok = ZeroMatrix(k, 1) # mean zero >>> R = MatrixSymbol('R', k, k) # symbolic covariance matrix >>> noise = Normal(zerok, R, 'eta') # Density after noise added in? >>> Density(H*X + noise) # This is a Block matrix ⎛[H I]⋅⎡μ⎤, [H I]⋅⎡Σ 0⎤⋅⎡H'⎤⎞ ⎜ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥⎟ ⎝ ⎣0⎦ ⎣0 R⎦ ⎣I ⎦⎠ # When we collapse the above expression it looks much nicer >>> block_collapse(Density(H*X + noise)) (H⋅μ, R + H⋅Σ⋅H') # Now lets imagine that we observe some value of HX+noise, # what does that tell us about X? How does our prior distribution change? >>> data = MatrixSymbol('data', k, 1) >>> Density(X , Eq(H*X+noise, data) ) # Density of X given HX+noise==data # I'm switching to the latex printer for this
# Again, this block matrix expression can be collapsed to the following >>> block_collapse(Density(X, Eq(H*X+noise, data) )) μ + Σ⋅H'⋅(R + H⋅Σ⋅H')^-1⋅(-H⋅μ + -data) , (I + -Σ⋅H'⋅(R + H⋅Σ⋅H')^-1⋅H)⋅Σ
This is the multivariate case of my previous post on data assimilation. Effectively all I’ve done here is baked in the logic behind the Kalman Filter and exposed it through my statistics operators Density, Given, etc… so that it has become more approachable.
1) This is all untested. Please let me know if something is wrong. Already I see an error with the latex printing.
2) For organizational reasons it seems unlikely that Matrix Expressions will make it into SymPy in their current form. As a result this code probably won’t make it into SymPy any time soon.
My active branch is over here:
with the multivariate normal code here:
The matrices live here: