# 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. Mar 8, 2017 at 19:45
• Are you specifically interested in the answer in 4D? Mar 8, 2017 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 Mar 9, 2017 at 0:26
• Have you already tried transforming the integral in Carl's answer to hyperspherical coordinates? Mar 22, 2017 at 4:36

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