Defining a set of random numbers

Question

I am attempting to use the sympy.stats module to help compute the expectation value of an expression. In the simplest form my problem is the following:

I know how to define two normally distributed random variables and compute the expectation

from sympy.stats import Normal, E
import sympy as sp
sigma = sp.symbols('sigma', positive=True)

x1 = Normal("x1", 0, sigma)
x2 = Normal("x2", 0, sigma)
print E(x1 * x2)

Out: 0

However, I would like to be able to define say X which is a vector (or list) of N normally distributed random variable. The trouble is I need N to remain arbitrary so I can't create a list or anything.

Does anyone have any suggestions on how to proceed?

The reason I need to do this is that I will be taking summation over the components of X, some of them will be correlated and some won't.

Answer 1

I take it you're looking for something like the following?

mu = MatrixSymbol('mu', n, 1)
Sigma = MatrixSymbol('Sigma', n, n)

X = Normal('X', mu, sigma)  # A multivariate normal distribution

If so then you're currently out of luck. SymPy stats doesn't currently support anything more than the combintation of several univariate distributions, we don't handle multivariate distributions in generality.

If you're willing to go numeric then you may find PyMC3 of interest.