# Estimating angle between two random vectors

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]

• I haven't checked in detail, but almost certainly all the delayed assignments are horribly slowing down your function. – LLlAMnYP Mar 8 '17 at 19:45
• Are you specifically interested in the answer in 4D? – Thies Heidecke Mar 8 '17 at 20:12
• I think if I can show it for 4D with high numeric confidence, I'll have good confidence this formula works for arbitrary dimensions – Yaroslav Bulatov Mar 9 '17 at 0:26
• Have you already tried transforming the integral in Carl's answer to hyperspherical coordinates? – J. M. is away Mar 22 '17 at 4:36

TracePrint[

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 -> StatisticsNExpectationDump$NExpectationAccuracyGoal, Compiled -> StatisticsNExpectationDump$NExpectationCompiled, PrecisionGoal -> StatisticsNExpectationDump$NExpectationPrecisionGoal, WorkingPrecision -> StatisticsNExpectationDump$NExpectationWorkingPrecision, MinRecursion -> StatisticsNExpectationDumpminrec$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?