Negative probability?

Question

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.

Answer 1

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.

Answer 2

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

{{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[1]^2 + a[2]^2 + a[3]^2 + a[4]^2 + a[5]^2 + a[6]^2 + a[7]^2 + a[8]^2 < 
     1, {a[1], a[2], a[3], a[4], a[5], a[6], a[7], a[8]} \[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[8]

{{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[1]^2 - a[2]^2 - a[3]^2 - a[4]^2 - a[5]^2])]*
         (-Sqrt[1 - a[1]^2 - a[2]^2 - a[3]^2 - a[4]^2 - a[5]^2])^(3/2)*
         (1 - a[1]^2 - a[2]^2 - a[3]^2 - a[4]^2 - a[5]^2), {a[1], -1, 1}, 
       {a[2], -Sqrt[1 - a[1]^2], Sqrt[1 - a[1]^2]}, 
       {a[3], -Sqrt[1 - a[1]^2 - a[2]^2], Sqrt[1 - a[1]^2 - a[2]^2]}, 
       {a[4], -Sqrt[1 - a[1]^2 - a[2]^2 - a[3]^2], 
         Sqrt[1 - a[1]^2 - a[2]^2 - a[3]^2]}, 
       {a[5], -Sqrt[1 - a[1]^2 - a[2]^2 - a[3]^2 - a[4]^2], 
         Sqrt[1 - a[1]^2 - a[2]^2 - a[3]^2 - a[4]^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[10]

{{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

Answer 3

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

Answer 4

Just another way to calculate volume of hypersphere and then relevant probability using recursion:

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

The probabilities:

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