Negative probability?

Bug introduced in 9.0.1 and fixed in 10.0.2

I am trying to get the sum of the squares of seven random variables, all uniformly distributed. This is what I tried.

Probability[
a^2 + b^2 + c^2 + d^2 + e^2 + f^2 + g^2 < 1,
{a, b, c, d, e, f, g} \[Distributed]
UniformDistribution[{{0, 1}, {0, 1}, {0, 1}, {0, 1}, {0, 1}, {0, 1}, {0, 1}}]]


The output was

-π^3/840


How can I get Mathematica to get the right answer?

Warning: the code takes about seven minutes to run on my machine.

• What version of Mathematica are you using? Version 8.0.4 returns Pi^3/840, which is correct. Version 9.0.1 has been working for much longer than 8.0.4 and finally returned -Pi^3/840, which looks to be a bug (possibly in the symbolic intergration code). – Szabolcs May 30 '14 at 18:42
• Just for fun: "Negative energies and probabilities should not be considered as nonsense. They are well-defined concepts mathematically, like a negative of money." - Paul Dirac – eldo May 30 '14 at 18:47
• That's weird: I'm using MMA 9.0.1, Windows, 64bit, and I get \[Pi]^3/840. – Niki Estner May 30 '14 at 18:47
• @nikie I ran this on OS X. Integrate[Boole[a^2 + b^2 + c^2 + d^2 + e^2 + f^2 + g^2 < 1], {a, -1, 1}, {b, -1, 1}, {c, -1, 1}, {d, -1, 1}, {e, -1, 1}, {f, -1, 1}, {g, -1, 1}] also gives a wrong result: it has an extra minus sign. (Somehow I though that extending the integration domain would help...) – Szabolcs May 30 '14 at 18:53
• 9.0.1.0 on Windows (64-bit) reproduces the issue. Weirdly the platform and version seem to be the same as @nikie. – ziyuang May 30 '14 at 21:24

This is not an answer (yet). Rather it explores the question in more depth.

n = 8;
parameters = ConstantArray[{0, 1}, n];
variables =
Symbol /@ CharacterRange["a", FromCharacterCode[ToCharacterCode["a"] + n - 1]];


The following takes a long time to evaluate, but the results it produces reveal give us a better view of the problem with Probability.

Table[
With[{params = parameters[[;; i]], vars = variables[[;; i]]},
Probability[vars.vars < 1, vars \[Distributed] UniformDistribution[params]]],
{i, n}]

{1, π/4, π/6, π^2/32, π^2/60, π^3/384, -(π^3/840), π^4/6144}


Only in the penultimate case is the result incorrect (its magnitude is correct but its sign is bogus). Why is 7 a (black) magical number?

I am going to send a query about this to WRI tech support. I will update this post, quoting their response, after I receive it.

Update

I have received an answer to the query I sent to WRI tech support. I quote the relevant part:

The function Probability does behave inappropriately in Mathematica 9 on Mac for i = 7 in the given notebook. However, this issue has been fixed and the performance is also improved in the prerelease version of Mathematica. Thank you for bringing it to our attention.

It seems that WRI was already aware of the problem and has fixed it in Mathematica 10.

This is also not an answer but a brief study in

\$Version

"8.0 for Microsoft Windows (64-bit) (October 7, 2011)"


which might be of interest. It considers three methods of calculating the requested probability.

It shows that in this version there is no negative probability but there is still a "critical number" which amounts to 8. For n = 8 the calculation for two methods is not completely performed but the result is returned in the form of commands still to be done.

Method 1

use Probability[] as in the original formulation of the problem

p[n_] := Probability[Sum[a[i]^2, {i, 1, n}] < 1,
Table[a[i], {i, 1, n}] \[Distributed]
UniformDistribution[Array[{0, 1} &, n]]]


Calculate some values and measure time required

{#, Timing[p[#]] // Reverse} & /@ Range

{{1, {1, 0.016}}, {2, {\[Pi]/4, 0.14}}, {3, {\[Pi]/6, 0.265}}, {4, {\[Pi]^2/
32, 1.155}}, {5, {\[Pi]^2/60, 4.789}}, {6, {\[Pi]^3/384, 15.007}},
{7, {\[Pi]^3/840, 34.336}}, {8, {Probability[
a^2 + a^2 + a^2 + a^2 + a^2 + a^2 + a^2 + a^2 <
1, {a, a, a, a, a, a, a, a} \[Distributed]
UniformDistribution[{{0, 1}, {0, 1}, {0, 1}, {0, 1}, {0, 1}, {0, 1}, {0,
1}, {0, 1}}]], 836.774}}}


Method 2

use Integrate with Boole

Notice that f(n) = 2^n p(n)

f[n_] := Integrate[Boole[Sum[a[i]^2, {i, 1, n}] < 1], Sequence @@ Table[{a[i], -1, 1}, {i, 1, n}]]


Calculate some values and measure time required

{#, Timing[f[#]] // Reverse} & /@ Range

{{1, {2, 0.015}}, {2, {\[Pi], 0.203}}, {3, {(4 \[Pi])/3,
0.328}}, {4, {\[Pi]^2/2, 1.716}}, {5, {(8 \[Pi]^2)/15,
8.018}}, {6, {\[Pi]^3/6, 32.682}}, {7, {(16 \[Pi]^3)/105,
145.877}}, {8, {(-(4/3))*Pi*Integrate[
Sqrt[-(1/Sqrt[1 - a^2 - a^2 - a^2 - a^2 - a^2])]*
(-Sqrt[1 - a^2 - a^2 - a^2 - a^2 - a^2])^(3/2)*
(1 - a^2 - a^2 - a^2 - a^2 - a^2), {a, -1, 1},
{a, -Sqrt[1 - a^2], Sqrt[1 - a^2]},
{a, -Sqrt[1 - a^2 - a^2], Sqrt[1 - a^2 - a^2]},
{a, -Sqrt[1 - a^2 - a^2 - a^2],
Sqrt[1 - a^2 - a^2 - a^2]},
{a, -Sqrt[1 - a^2 - a^2 - a^2 - a^2],
Sqrt[1 - a^2 - a^2 - a^2 - a^2]}], 2290.718}}}


Method 3

volume of n-dimensional unit sphere

Exact values up to n=10 are swiftly calculated by this method for reference. Using polar ccordinates (see e.g. http://de.wikipedia.org/wiki/Polarkoordinaten#n-dimensionale_Polarkoordinaten) we have for the volume

v[n_] := Integrate[
r^(n - 1) Product[Sin[u[j]]^j, {j, 1, n - 2}], {r, 0, 1}, {phi, 0, 2 \[Pi]},
Sequence @@ Table[{u[j], 0, \[Pi]}, {j, 1, n - 2}]]


Notice that v(n) = 2^n p(n).

Calculate some values and measure time required

{#, Timing[v[#]] // Reverse} & /@ Range

{{1, {2 \[Pi], 0.}}, {2, {\[Pi], 0.}}, {3, {(4 \[Pi])/3,
0.063}}, {4, {\[Pi]^2/2, 0.14}}, {5, {(8 \[Pi]^2)/15,
0.234}}, {6, {\[Pi]^3/6, 0.375}}, {7, {(16 \[Pi]^3)/105,
0.483}}, {8, {\[Pi]^4/24, 0.702}}, {9, {(32 \[Pi]^4)/945,
2.621}}, {10, {\[Pi]^5/120, 4.103}}}


I have checked that there were no memory shortages on my PC. Hence there must be an intrisic reason in MMA for magic number 8.

Regards, Wolfgang

Fixed in 10.0.2 Probability[a^2 + b^2 + c^2 + d^2 + e^2 + f^2 + g^2 < 1, {a, b, c, d, e, f, g}
\[Distributed]  UniformDistribution[{{0, 1}, {0, 1}, {0, 1}, {0, 1}, {0, 1},
{0, 1}, {0, 1}}]] Just another way to calculate volume of hypersphere and then relevant probability using recursion:

vs[n_] :=
Most@Nest[{#[]/(#[] + 1), 2 Pi #[], #[] + 1} &, {1, 2,
0}, n]
v[n_] := vs[n][]


The probabilities:

Grid[Prepend[{#1, #2, N@#2} & @@@ ({#, v[#]/2^#} & /@ Range),
Style[#, Bold] & /@ {"n",
"P[\!$$\*UnderoverscriptBox[\(\[Sum]$$, $$j = 1$$, \
$$n$$]\)\!$$\*SuperscriptBox[SubscriptBox[\(x$$, $$j$$], $$2$$]\)<1|x\
\[Distributed]Unif(0,1\!$$\*SuperscriptBox[\()$$, $$n$$]\)]", ""}]] 