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]