Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

share|improve this question
2  
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 at 18:42
2  
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 at 18:47
2  
That's weird: I'm using MMA 9.0.1, Windows, 64bit, and I get \[Pi]^3/840. –  nikie May 30 at 18:47
1  
@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 at 18:53
2  
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 at 21:24

2 Answers 2

up vote 13 down vote accepted

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.

share|improve this answer

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

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.