I'm interested in angles between random vectors on unit ball. Expected value of cosine is 0 by symmetry, so, expectation of cosine squared is the interesting bit.

The following NExpectation takes several minutes and gives error margins that are 50% of the answer ...is this a fundamentally hard problem for numerical integration, or are there tricks to make it easier?

x := {x1, x2, x3, x4};
y := {y1, y2, y3, y4};
normal := MultinormalDistribution[{0, 0, 0, 0}, IdentityMatrix[4]];
vars := {x \[Distributed] normal, 
  y \[Distributed] 
NExpectation[(x.y/(Norm[x] Norm[y]))^2, vars]


NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small.

PS, the real motivation is to verify that the following formula (obtained through guessing) works for vectors sampled from more general distribution --

val1[a_, b_, c_, 
  d_] := (a + b + c + d)/(Sqrt[a] + Sqrt[b] + Sqrt[c] + Sqrt[d])^2
val2[a_, b_, c_, d_] := Block[{},
  x := {x1, x2, x3, x4};
  y := {y1, y2, y3, y4};
  normal := 
   MultinormalDistribution[{0, 0, 0, 
     0}, {{a, 0, 0, 0}, {0, b, 0, 0}, {0, 0, c, 0}, {0, 0, 0, d}}];
  vars := {x \[Distributed] normal, y \[Distributed] normal};
  NExpectation[(x.y/(Norm[x] Norm[y]))^2, vars]]
val1[0.5, 1.5, 2.5, 3.5] - val2[0.5, 1.5, 2.5, 3.5]
  • $\begingroup$ I haven't checked in detail, but almost certainly all the delayed assignments are horribly slowing down your function. $\endgroup$
    – LLlAMnYP
    Mar 8, 2017 at 19:45
  • $\begingroup$ Are you specifically interested in the answer in 4D? $\endgroup$ Mar 8, 2017 at 20:12
  • $\begingroup$ I think if I can show it for 4D with high numeric confidence, I'll have good confidence this formula works for arbitrary dimensions $\endgroup$ Mar 9, 2017 at 0:26
  • $\begingroup$ Have you already tried transforming the integral in Carl's answer to hyperspherical coordinates? $\endgroup$ Mar 22, 2017 at 4:36

1 Answer 1


We can peek at the integral Mathematica is performing under the hood:

    NExpectation[(x.y/(Norm[x] Norm[y]))^2, vars], 

HoldForm[NIntegrate[(E^((-x1^2 - x2^2 - x3^2 - x4^2)/2 + (-y1^2 - y2^2 - y3^2 - y4^2)/2)*(x1*y1 + x2*y2 + x3*y3 + x4*y4)^2)/(16*Pi^4*(Abs[x1]^2 + Abs[x2]^2 + Abs[x3]^2 + Abs[x4]^2)*(Abs[y1]^2 + Abs[y2]^2 + Abs[y3]^2 + Abs[y4]^2)), {x1, -Infinity, Infinity}, {x2, -Infinity, Infinity}, {x3, -Infinity, Infinity}, {x4, -Infinity, Infinity}, {y1, -Infinity, Infinity}, {y2, -Infinity, Infinity}, {y3, -Infinity, Infinity}, {y4, -Infinity, Infinity}, AccuracyGoal -> Statistics`NExpectationDump`$NExpectationAccuracyGoal, Compiled -> Statistics`NExpectationDump`$NExpectationCompiled, PrecisionGoal -> Statistics`NExpectationDump`$NExpectationPrecisionGoal, WorkingPrecision -> Statistics`NExpectationDump`$NExpectationWorkingPrecision, MinRecursion -> Statistics`NExpectationDump`minrec$24078, Sequence[]]]


As you can see, it is an 8-dimensional integral, so I'm not surprised that Mathematica is having difficulties. Perhaps you can work with this integral directly instead, and play with Method options, or alternate parametrizations?


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.