1
$\begingroup$

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] 
   normal};
NExpectation[(x.y/(Norm[x] Norm[y]))^2, vars]

Get

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]
Improve this question
$\endgroup$
4
  • $\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$
    – Thies Heidecke
    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$
    – Yaroslav Bulatov
    Mar 9, 2017 at 0:26
  • $\begingroup$ Have you already tried transforming the integral in Carl's answer to hyperspherical coordinates? $\endgroup$
    – J. M.'s persistent exhaustion ♦
    Mar 22, 2017 at 4:36

1 Answer 1

Reset to default
5
$\begingroup$

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

TracePrint[
    NExpectation[(x.y/(Norm[x] Norm[y]))^2, vars], 
    _NIntegrate,
    TraceAction->((Print[InputForm@#];Abort[])&)
]

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[]]]

$Aborted

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?

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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