How to deal with complicated gaussian integrals in Mathematica?

Question

As we know, for most Gaussian integrals, we can get the analytical result. Now I have many Gaussian integrals to treat, which have the following general form,

Integrate[
   x1^n1 x2^n2 x3^n3 ... xd^
     nd Exp[-{x1, x2, x3 ... xd}.matrix.{x1, x2, x3 ... xd} + c1 x1 + 
        c2 x2 ... cd xd], {x1, -Infinity, Infinity},
     {x2, -Infinity, Infinity} ... {xd, -Infinity, Infinity}]. 

The matrix is a symmetrical positive definite matrix.

We can get the analytical integral result for this general form. However, it is a difficult task for Mathematica to deal with the similiar expressions, which takes so much time. I have written a package for treating this problem.

But, it is still not sufficient because I need to deal with lots of these kinds of integral expressions. So the computation time is important.

So have you ever seen the package for doing this kind of Gaussian integral? Or can you give me some suggestions about this problem? Thank you very much!

Thanks all of you. For the concept MultinormalDistribution, it is very helpful for me.

By the way, How can I transform the following integral expression to the standard form of MultinormalDistribution in a simple way?

vec1 := {r1x, r1y, r1z};
vec2 := {r2x, r2y, r2z};
vect := {r1x, r1y, r1z, r2x, r2y, r2z};
vcof := {c1, c2, c3, c4, c5, c6};
Integrate[Exp[-vect.mat.vect + vcof.vect] (vec1.vec2)^n (vec1.vec1)^
l (vec2.vec2)^m, {r1x, -Infinity, Infinity}, {r1y, -Infinity, 
Infinity}, {r1z, -Infinity, Infinity}, {r2x, -Infinity, 
Infinity}, {r2y, -Infinity, Infinity}, {r2z, -Infinity, Infinity}];

The mat is a symmetrical positive definite sparse matrix. And n, l, m are integers.
Thanks!

Answer 1

What you have is a MultinormalDistribution. The quadratic and linear forms in the exponential can be rewritten in terms of $\frac12(\vec{x}-\vec{\mu})^\top\Sigma^{-1}(\vec{x}-\vec{\mu})$ where $\vec{\mu}$ represents the mean and $\Sigma$ the covariance matrix, see the documentation.

With this, you can do integrals of the type given in the question by invoking Expectation, as in this example:

Expectation[
 x^2 y^3, {x, y} \[Distributed] 
  MultinormalDistribution[{μ1, μ2}, 
  {{σ1^2, ρ σ1 σ2}, 
   {ρ σ1 σ2, σ2^2}}]]

The result is:

$\text{$\mu $1}^2 \text{$\mu $2}^3+\text{$\mu $2}^3 \text{$\sigma $1}^2+6 \text{$\mu $1} \text{$\mu $2}^2 \rho \text{$\sigma $1} \text{$\sigma $2}+3 \text{$\mu $1}^2 \text{$\mu $2} \text{$\sigma $2}^2+3 \text{$\mu $2} \text{$\sigma $1}^2 \text{$\sigma $2}^2+6 \text{$\mu $2} \rho ^2 \text{$\sigma $1}^2 \text{$\sigma $2}^2+6 \text{$\mu $1} \rho \text{$\sigma $1} \text{$\sigma $2}^3$

Edit

Regarding the normalization prefactor mentioned in Sjoerd's comment, you can use the fact that for any dimension $n$

$\iint\exp(-\frac{1}{2}\vec{z}^\top \Sigma^{-1}\vec{z})\mathrm dz^n = (2\pi)^{n/2}\sqrt{\det(\Sigma)}$

Hopefully these hints will be enough for you to fill in the missing linear-algebra steps to make the connection to your given matrix matrix.

Edit 2

In response to the comment by chris, for polynomials as prefactors one can also use the slightly simpler but equivalent form

Moment[
 MultinormalDistribution[{μ1, μ2}, 
  {{σ1^2, ρ σ1 σ2}, 
   {ρ σ1 σ2, σ2^2}}],
 {2, 3}]

This is the same example as above, with the powers of x and y appearing in the second argument. See the documentation for Moment.

The difference between Moment and Expectation is that Moment is restricted to the expectation values of polynomials.

Edit 3

Before going on with the symbolic manipulations that I assumed are desired here, let me also point out that you can do your integrals pretty straightforwardly if your integrand contains no symbolic parameters. Then you just need to do a numerical integral by replacing Integrate with NIntegrate.

But now back to the symbolic part:

A follow-up question arose how to complete the square in the exponential to get to the standard form of the multinormal distribution, starting from a form like this:

$$\exp(\,\vec{x}^\top A\vec{x}+\vec{v}^\top\vec{x})$$

The matrix $A$ in the exponential is symmetric, $A^\top=A$, and positive definite. Therefore $A$ is invertible, and the inverse is symmetric,

$$\left(A^{-1}\right)^\top=A^{-1}$$

With this, you can verify

$$\left(\vec{x}+\frac{1}{2}A^{-1}\vec{v}\right)^\top A\left(\vec{x}+\frac{1}{2}A^{-1}\vec{v}\right)=\vec{x}^\top A\vec{x}+\vec{v}^\top\vec{x}+\frac{1}{4}\vec{v}^\top A^{-1}\vec{v} $$

by directly multiplying out the factors on the left. Therefore,

$$\exp(\vec{x}^\top A\vec{x}+\vec{v}^\top \vec{x})=\exp(\left(\vec{x}-\vec{\mu}\right)^\top A\left(\vec{x}-\vec{\mu}\right)-\frac{1}{4}\vec{v}^\top A^{-1}\vec{v})$$

where

$$\vec{\mu}\equiv-\frac{1}{2}A^{-1}\vec{v} $$

Compare this to the standard form of the Gaussian integral, and you see that in the notation of Mathematica's documentation

$$A \equiv \Sigma^{-1}$$

and our integral differs from the standard Gaussian one by a factor

$$\exp(-\frac{1}{4}\vec{v}^\top A^{-1}\vec{v})$$

Now we have all the pieces that are needed, except that you still have to calculate the inverse matrix $A^{-1} \equiv \Sigma$, using

Inverse[mat]

if I go back to your original notation where the matrix $A$ is called mat.