# Moments of Multinormal in symbolic form

Gaussian expectations often simplify (20.5 of Seber's "A Matrix Handbook for Statisticians" gives many examples). I'm wondering if there's a way to use Mathematica to produce the formulas semi-automatically. I've recently needed to get the following formulas:

$$E[x'xx'x]=\text{tr}(\Sigma)^2+2\text{tr}(\Sigma^2)$$ $$E[(x'x)xx']=\text{tr}(\Sigma)\Sigma+2\Sigma^2$$ $$E\left[\frac{x'\Sigma x}{x'x}\right]=?$$ My approach has been to take a couple of diagonal $$\Sigma$$, run Expectation and then try to guess the formula from the result. Is there a smarter approach?

Clear[x, s];
d = 4;
vec = Array[x, {d}];
sigma = DiagonalMatrix[Table[k, {k, 1, d}]];
mean = Array[0 &, {d}];
result = Expectation[(vec.vec) Outer[Times, vec, vec],
vec \[Distributed] MultinormalDistribution[mean, sigma]]
result == Tr[sigma] sigma + 2 sigma^2

• How about \ [CapitalSigma] = {{Subscript[\ [Sigma], 1]^2, \ [Rho] Subscript[\ [Sigma], 1] Subscript[\ [Sigma], 2]}, {\ [Rho] Subscript[\ [Sigma], 1] Subscript[\ [Sigma], 2], Subscript[\ [Sigma], 2]^2}}; Moment[ MultinormalDistribution[{Subscript[\ [Mu], 1], Subscript[\ [Mu], 2]}, \ [CapitalSigma]], {3, 2}]? The result is too long for this comment, Nov 4, 2019 at 19:43
• I wonder if all one can do is write up a bunch of replacement rules to mimic what's in Seber's book. My feeble attempts are solving the "?" equation does not result in a symbolic form even for a bivariate normal.
– JimB
Nov 6, 2019 at 17:25
• For bivariate normal and ??? equation I get Tr[sigma]-Sqrt[Det[sigma]] using the method above, however, couldn't figure out a way to extend it to 3 variables, both Expectation and NExpectation take too long to generate guesses Nov 6, 2019 at 18:11
• The answer to computing equation 3 may be somewhere here sciencedirect.com/science/article/pii/S0047259X13000298 Nov 6, 2019 at 18:17
• I was using a non-diagonal $\Sigma$ and integrating directly using Integrate but wasn't having any success. I'll take a closer look at the reference you provided.
– JimB
Nov 6, 2019 at 20:11

Following @user64494 suggestion, looks like Moment works better than Expectation here.

The first equation can be rewritten as a sum of moments

xs = {x1, x2, x3};
momentSpecs = First /@ CoefficientRules[Total[xs*xs]*Total[xs*xs]];
moments =
Moment[MultinormalDistribution[{0, 0, 0},
DiagonalMatrix[{s1, s2, s3}]], #] & /@ momentSpecs;
Total[moments]


For second equation, since result is diagonal, symmetry means it's sufficient to get expression for the first entry

xs = {x1, x2, x3};
momentSpecs = First /@ CoefficientRules[x1^2*Total[xs*xs]];
moments =
Moment[MultinormalDistribution[{0, 0, 0},
DiagonalMatrix[{s1, s2, s3}]], #] & /@ momentSpecs;
Total[moments]


For third equation, the method above doesn't quite work. The answer may lay somewhere in https://www.sciencedirect.com/science/article/pii/S0047259X13000298 , but for the case of bivariate normal, the following seems to hold

$$E\left[\frac{x'\Sigma x}{x'x}\right]=\text{tr}(\Sigma)-\sqrt{\text{Det}(\Sigma)}$$

Here is what might be considered a partial answer.

Suppose that $$x=(x_1, x_2,\ldots,x_n)'$$ has a multivariate normal distribution with all means equal to $$0$$, variances equal to $$1$$ and all covariances equal to zero. We want to know $$E[x' \Sigma x/x'x]$$ where $$\Sigma$$ is a symmetric square matrix.

The expectation is $$\text{tr}(\Sigma)/n$$. Below is a handwaving argument.

Define two functions that give $$\Sigma$$ and $$x$$:

Σ[n_] := Table[ToExpression["σ" <> ToString[Min[i, j]] <> ToString[Max[i, j]]], {i, n}, {j, n}]
x[n_] := Table[ToExpression["x" <> ToString[i]], {i, n}]


We see the following for $$x'\Sigma x/(x'x)$$:

x[2].Σ[2].x[2]/(x[2].x[2]) // Expand


x[3].Σ[3].x[3]/(x[3].x[3]) // Expand


There are two kinds of terms. One with $$x_i^2$$ in the numerator and one with $$x_i x_j$$ in the numerator. The expectation of the terms with $$x_i x_j$$ in the numerator is zero.

d2 = TransformedDistribution[(x1 x2/(x1^2 + x2^2)),
{x1, x2} \[Distributed] MultinormalDistribution[IdentityMatrix[2]]]
Mean[d2]
(* 0 *)

d3 = TransformedDistribution[(x1 x2/(x1^2 + x2^2 + x3^2)),
{x1, x2, x3} \[Distributed] MultinormalDistribution[IdentityMatrix[3]]]
Mean[d3]
(* 0 *)


For the terms with $$x_i^2$$ in the numerator, all of the random variables in those terms are squared which have independent $$\chi^2_1$$ distributions. The means can be determined in the following manner:

d2 = TransformedDistribution[y1/(y1 + y2), {y1 \[Distributed] ChiSquareDistribution[1],
y2 \[Distributed] ChiSquareDistribution[1]}];
Mean[d2]
(* 1/2 *)

d3 = TransformedDistribution[y1/(y1 + y2plusy3), {y1 \[Distributed] ChiSquareDistribution[1],
y2plusy3 \[Distributed] ChiSquareDistribution[2]}];
Mean[d3]
(* 1/3 *)

d4 = TransformedDistribution[y1/(y1 + y2plusy3plusy4), {y1 \[Distributed] ChiSquareDistribution[1],
y2plusy3plusy4 \[Distributed] ChiSquareDistribution[3]}];
Mean[d4]
(* 1/4 *)


Note that $$y_2+y_3 \sim \chi^2_2$$ and $$y_2+y_3+y_4 \sim \chi^2_3$$.

So putting all of this together the expectation is $$\text{tr}(\Sigma)/n$$.

It appears that your original question has $$\Sigma$$ (the symmetric square matrix) identical to the covariance matrix of the multivariate normal. But the reference you gave starts out the $$x_i$$ random variables having independent unit normals. This is why I've labeled this as a "partial" answer.

• Hm, this doesn't seem to add up for bivariate normal, where I verified the following formula to give correct result $\text{Tr}(\Sigma)-\sqrt{\text{Det}(\Sigma)}$ Nov 7, 2019 at 18:50
• BTW, you can assume diagonal sigma because of symmetry. The paper I gave further reduces it to the case of standard normal by applying a change of coordinates. It's sufficient to solve the following for standard normal and diagonal sigma: $E[x\Sigma^2 x]/E[x \Sigma x]$ Nov 7, 2019 at 18:53
• I think you need to put in your post the complete set of assumptions you're making. What I have explicitly assumes that the covariance matrix for the random vector $x$ is the identity matrix. It appears that you are using $\Sigma$ for both the covariance matrix and the matrix $A$ in $x'Ax/(x'x)$. But making that assumption explicit would be helpful to others to know exactly what you need.
– JimB
Nov 7, 2019 at 20:58
• Sorry I assumed it would be more clear from the code supplied, $\Sigma$ is the covariance matrix of my Gaussian, there are no other matrices involved. Nov 7, 2019 at 21:29

For multivariate normal distributions with mean vector $$\mu$$ and covariance matrix $$\Sigma$$ I've not had any success in getting a closed form for $$E[x'\Sigma x/(x'x)]$$ other than for a bivariate normal distribution when the mean vector is all zeros.

For whatever it's worth here's a numerical solution using the generating function in the reference you supplied:

(* From https://www.sciencedirect.com/science/article/pii/S0047259X13000298 *)

(* Generating function *)
ϕ[t1_, t2_, μ_, A_, B_] := Det[IdentityMatrix[Length[μ]] - 2 t1 A - 2 t2 B]^(-1/2) *
Exp[μ.(Inverse[IdentityMatrix[Length[μ]] - 2 t1 A - 2 t2 B] - IdentityMatrix[Length[μ]]).μ/2]

(* Expectation of (x.A.x)^p/(x.B.x)^q where p and q are integers > 0 *)
(* given x is a multivariate normal with mean μ and covariance matrix Σ. *)
expect[μ_, Σ_, A_, B_, p_, q_] := Module[{a, b, d, g, gInv},
g = MatrixPower[Σ, 0.5];
gInv = MatrixPower[Σ, -0.5];
a = g.A.g;
b = g.B.g;
μ0 = gInv.μ;
d = t^(q - 1) ((D[ϕ[t1, t2, μ0, a, b], {t1, p}]) /. t1 -> 0 /. t2 -> -t);
NIntegrate[d, {t, 0, ∞}]/Gamma[q]
]

(* Example *)
μ = {35, 5, 2};
Σ = {{1, 3/4, 1/2}, {3/4, 3, 1/8}, {1/2, 1/8, 7}};
{p0, q0} = {3, 2};
A = Σ;
B = IdentityMatrix[Length[μ]];
e = expect[μ, Σ, A, B, p0, q0]
(* 3331.9498311917214 *)

(* Check using random sample *)
SeedRandom[1234];
z = RandomVariate[MultinormalDistribution[μ, Σ], 1000000];
Mean[(#.A.#)^p0/(#.B.#)^q0 & /@ z]
(* 3331.9016258162883 *)


There appears to be a closed form for a trivariate normal when there is a diagonal covariance matrix and 2 of the 3 variances are equal and the mean is $$(0,0,0)$$. Other possibilities could be determined by changing NIntegrate to Integrate and 0.5 in the MatrixPower functions to 1/2.