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
  • 1
    $\begingroup$ 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, $\endgroup$
    – user64494
    Nov 4, 2019 at 19:43
  • $\begingroup$ 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. $\endgroup$
    – JimB
    Nov 6, 2019 at 17:25
  • $\begingroup$ 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 $\endgroup$ Nov 6, 2019 at 18:11
  • $\begingroup$ The answer to computing equation 3 may be somewhere here sciencedirect.com/science/article/pii/S0047259X13000298 $\endgroup$ Nov 6, 2019 at 18:17
  • $\begingroup$ 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. $\endgroup$
    – JimB
    Nov 6, 2019 at 20:11

3 Answers 3


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;

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;

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

Dimension 2 equation

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

Dimension 3 equation

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]]]
(* 0 *)

d3 = TransformedDistribution[(x1 x2/(x1^2 + x2^2 + x3^2)), 
  {x1, x2, x3} \[Distributed] MultinormalDistribution[IdentityMatrix[3]]]
(* 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]}];
(* 1/2 *)

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

d4 = TransformedDistribution[y1/(y1 + y2plusy3plusy4), {y1 \[Distributed] ChiSquareDistribution[1], 
  y2plusy3plusy4 \[Distributed] ChiSquareDistribution[3]}];
(* 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.

  • $\begingroup$ 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)}$ $\endgroup$ Nov 7, 2019 at 18:50
  • $\begingroup$ 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]$ $\endgroup$ Nov 7, 2019 at 18:53
  • $\begingroup$ 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. $\endgroup$
    – JimB
    Nov 7, 2019 at 20:58
  • $\begingroup$ 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. $\endgroup$ 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 *)
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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.