# Generic Matrices

A multivariate Gaussian distribution for $k$ dimensions looks like this:

$$\frac{1}{\sqrt{(2\pi)^k|\mathbf \Sigma|}}\mathbf e^{ -\frac 1 2 \mathbf x^T \mathbf \Sigma^{-1}\mathbf x }$$

$\mathbf x$ is a $k$-dimensional vector, and $\mathbf \Sigma$ a $k$ by $k$ symmetric positive definite matrix. Now I want Mathematica to integrate this expression over all $\mathbf x$, the result should be 1. How can I do this? Nothing else is known, in particular we don't know the number of dimensions.

• If you can diagonalize sigma ... – Dr. belisarius Jul 2 '14 at 13:49
• @belisarius $\Sigma$ can be diagonalized, because it's symmetric. However, I don't have a specific matrix here, I'd like a general solution to this and similar problems. – pentadecagon Jul 2 '14 at 13:52
• Mathematica cannot do an integral where the number of dimensions is not specified. It is up to you to put the problem in a form where the integral can be carried out over known-dimensional spaces. (For example diagonalization decomposes it into 1D problems.) It is at that point only that Mathematica can start helping you. People might be able to show you how to do this if you explain your problem in a bit more detail. – Szabolcs Jul 2 '14 at 13:56
• Related quation: 3N-dimensional integral – Jens Jul 2 '14 at 16:04

There is indeed a generic expression for (the essential part of) your integral which leaves the dimension n open.

(We will use n istead of k and s = Sigma^(-1) for the matrix in the exponent).

Main idea

The main idea is to use the function Sequence (and, of course, delayed assingment).

Consider an informal expression of the type

Integrate[ f[x[1],x[2],...,x[n]], {x[1],a,b},{x[2],a,b}, ... , {x[n],a,b}]

How can we deal with the part {x[1],a,b},{x[2],a,b}, ... , {x[n],a,b} ?

The answer is provided by Sequence:

{x[1],a,b},{x[2],a,b}, ... , {x[n],a,b} = Sequence@@Table{x[i],a,b}, {i,1,n}]

Normalization factor

Before we formulate the integral we consider the normalization factor. The essential part is the determinant of the matrix s. Notice that we are using s here insted of sigma^(-1). It turns out that it is convenient to take into account that s is without loss of generality a symmetric matrix.

We therefore define first the general determinant as

In[1]:= det[n_] := Det[Table[s[i, j], {i, 1, n}, {j, 1, n}]]


and then the symmterized form

In[2]:= dets[n_] :=
Simplify[det[n],
And @@ Flatten[Table[s[i, j] == s[j, i], {i, 1, n}, {j, 1, n}], 1]]


The integral

Now we have all ingredients available to write down the first form of the integral in question

In[3]:= h1[n_] := Sqrt[dets[n]/(2 \[Pi])^n]
Integrate[Exp[- 1/2 Sum[x[i] s[i, j] x[j], {i, 1, n}, {j, 1, n} ]],
Sequence @@ Table[{x[i], -\[Infinity], \[Infinity]}, {i, 1, n}]]


When you call the function h1 with a specific n you get the result in the form of a ConditionalExpression where the condition turns out to be in part that of the positive definiteness of the matrix s.

It is also possible to include some obvious conditions to be imposed on the matrix s as Assumptions. The next form of the integral is thus:

In[4]:= h2[n_] := Sqrt[dets[n]/(2 \[Pi])^n]
Integrate[Exp[- 1/2 Sum[x[i] s[i, j] x[j], {i, 1, n}, {j, 1, n} ]],
Sequence @@ Table[{x[i], -\[Infinity], \[Infinity]}, {i, 1, n}],
Assumptions -> Table[s[i, j], {i, 1, n}, {j, 1, n}] \[Element] Reals]


h2 takes into account, that the Matrix s is supposed to be real.

In[5]:= h3[n_] := Sqrt[dets[n]/(2 \[Pi])^n]
Integrate[Exp[- 1/2 Sum[x[i] s[i, j] x[j], {i, 1, n}, {j, 1, n} ]],
Sequence @@ Table[{x[i], -\[Infinity], \[Infinity]}, {i, 1, n}],
Assumptions -> {And @@ Table[s[i, i] > 0, {i, 1, n}],
Flatten[Table[s[i, j], {i, 1, n}, {j, 1, n}], 1] \[Element] Reals}]


h3 additionally takes into account the simplest reqirement of positive definiteness, viz. all diagonal elements have to be >0.

Examples

Examples for n = 1

In[6]:= n = 1;

In[7]:= h1[n]

Out[7]= ConditionalExpression[1, Re[s[1, 1]] > 0]


Mapping Simplify to this expression simplifies the result under the condition given:

In[8]:= Simplify @@ %

Out[8]= 1


Next form h2

In[9]:= h2[n]
Simplify @@ %

Out[9]= ConditionalExpression[1, s[1, 1] > 0]

Out[10]= 1


Next form h3

In[11]:= h3[n]

Out[11]= 1


Example n = 2

For simplicity we will employ only h3

In[12]:= n = 2;

In[13]:= h3[n]

Out[13]= ConditionalExpression[(
2 Sqrt[-s[2, 1]^2 + s[1, 1] s[2, 2]] Sqrt[
s[2, 2]/(-(s[1, 2] + s[2, 1])^2 + 4 s[1, 1] s[2, 2])])/Sqrt[
s[2, 2]], (s[1, 2] + s[2, 1])^2 < 4 s[1, 1] s[2, 2]]

In[14]:= Simplify @@ %

Out[14]= 2 Sqrt[(-s[2, 1]^2 + s[1, 1] s[2, 2])/(-(s[1, 2] + s[2, 1])^2 +
4 s[1, 1] s[2, 2])]


Now apply the condition of symmetry of the matrix s:

In[15]:= Simplify[%,
And @@ Flatten[Table[s[i, j] == s[j, i], {i, 1, n}, {j, 1, n}], 1]]

Out[15]= 1


Example n = 3

In[16]:= n = 3;

In[17]:= Timing[h3[n]]

Out[17]= {148.747, ConditionalExpression[(2 Sqrt[
1/(-(s[2, 3] + s[3, 2])^2 + 4 s[2, 2] s[3, 3])] Sqrt[
2 s[2, 1] s[3, 1] s[3, 2] - s[1, 1] s[3, 2]^2 - s[2, 1]^2 s[3, 3] +
s[2, 2] (-s[3, 1]^2 + s[1, 1] s[3, 3])])/(\[Sqrt](1/(s[2, 3]^2 +
2 s[2, 3] s[3, 2] + s[3, 2]^2 -
4 s[2, 2] s[3, 3])(s[1, 3]^2 s[2, 2] - s[1, 2] s[2, 3] s[3, 1] -
s[2, 1] s[2, 3] s[3, 1] + s[2, 2] s[3, 1]^2 -
s[1, 2] s[3, 1] s[3, 2] - s[2, 1] s[3, 1] s[3, 2] -
s[1, 3] (-2 s[2, 2] s[3, 1] + s[1, 2] (s[2, 3] + s[3, 2]) +
s[2, 1] (s[2, 3] + s[3, 2])) + s[1, 2]^2 s[3, 3] +
2 s[1, 2] s[2, 1] s[3, 3] + s[2, 1]^2 s[3, 3] +
s[1, 1] (s[2, 3]^2 + 2 s[2, 3] s[3, 2] + s[3, 2]^2 -
4 s[2, 2] s[3, 3])))), (s[2, 3]^2 + 2 s[2, 3] s[3, 2] +
s[3, 2]^2 - 4 s[2, 2] s[3, 3]) (s[1, 3]^2 s[2, 2] -
s[1, 2] s[2, 3] s[3, 1] - s[2, 1] s[2, 3] s[3, 1] + s[2, 2] s[3, 1]^2 -
s[1, 2] s[3, 1] s[3, 2] - s[2, 1] s[3, 1] s[3, 2] -
s[1, 3] (-2 s[2, 2] s[3, 1] + s[1, 2] (s[2, 3] + s[3, 2]) +
s[2, 1] (s[2, 3] + s[3, 2])) + s[1, 2]^2 s[3, 3] +
2 s[1, 2] s[2, 1] s[3, 3] + s[2, 1]^2 s[3, 3] +
s[1, 1] (s[2, 3]^2 + 2 s[2, 3] s[3, 2] + s[3, 2]^2 -
4 s[2, 2] s[3, 3])) > 0]}

In[18]:= Simplify @@ %[[2]]

Out[18]= (2 Sqrt[1/(-(s[2, 3] + s[3, 2])^2 + 4 s[2, 2] s[3, 3])] Sqrt[
2 s[2, 1] s[3, 1] s[3, 2] - s[1, 1] s[3, 2]^2 - s[2, 1]^2 s[3, 3] +
s[2, 2] (-s[3, 1]^2 + s[1, 1] s[3, 3])])/(\[Sqrt](1/(s[2, 3]^2 +
2 s[2, 3] s[3, 2] + s[3, 2]^2 - 4 s[2, 2] s[3, 3])(s[1, 3]^2 s[2, 2] -
s[1, 2] s[2, 3] s[3, 1] - s[2, 1] s[2, 3] s[3, 1] + s[2, 2] s[3, 1]^2 -
s[1, 2] s[3, 1] s[3, 2] - s[2, 1] s[3, 1] s[3, 2] -
s[1, 3] (-2 s[2, 2] s[3, 1] + s[1, 2] (s[2, 3] + s[3, 2]) +
s[2, 1] (s[2, 3] + s[3, 2])) + s[1, 2]^2 s[3, 3] +
2 s[1, 2] s[2, 1] s[3, 3] + s[2, 1]^2 s[3, 3] +
s[1, 1] (s[2, 3]^2 + 2 s[2, 3] s[3, 2] + s[3, 2]^2 -
4 s[2, 2] s[3, 3]))))


Again applying the condition of symmetry of the matrix s:

In[19]:= Simplify[%,
And @@ Flatten[Table[s[i, j] == s[j, i], {i, 1, n}, {j, 1, n}], 1]]

Out[19]= (Sqrt[1/(-s[3, 2]^2 + s[2, 2] s[3, 3])] Sqrt[
2 s[2, 1] s[3, 1] s[3, 2] - s[1, 1] s[3, 2]^2 - s[2, 1]^2 s[3, 3] +
s[2, 2] (-s[3, 1]^2 + s[1, 1] s[3, 3])])/Sqrt[(-2 s[2, 1] s[3, 1] s[3, 2] +
s[1, 1] s[3, 2]^2 + s[2, 1]^2 s[3, 3] +
s[2, 2] (s[3, 1]^2 - s[1, 1] s[3, 3]))/(s[3, 2]^2 - s[2, 2] s[3, 3])]


This is still not suffient to reveal it to be equal to 1. Taking the square removes the square root and considering that the integral is a positive quantity finally does the job.

In[35]:= Simplify[%^2]

Out[35]= 1


Remark: although it is nice that we have a theoretical framework available, its practical usefulness for higher n must be checked. The speed would be much higher if the matrix s is given explicitly.

Thanks for your question and best regards,

Wolfgang

This may not be your intention, however, please note Mathematica 9.0 has built-in MultinormalDistribution, where you can enter covariance matrix $\mathbb\Sigma$. For example:

mnd = MultinormalDistribution[{1, 2,
3}, {{2, 1/2, -1/3}, {1/2, 1, 0}, {-1/3, 0, 2/3}}];
funs = {Mean, Variance, Skewness, Kurtosis,
Composition[MatrixForm, Correlation]};
Grid[Transpose[{funs, Through[funs[mnd]]}], Frame -> All]