I want to calculate the distribution that comes from summing up random numbers between -0.5 and 0.5.

A computationally intensive version of doing this is to simply do the following:

   Sum[RandomReal[] - 0.5, {i, 1, 
     numberOfRandomNumbersAdded}], {1000000}], \
{numberOfRandomNumbersAdded, 1, 5}]]

enter image description here

This won't work however if the number of numbers to be summed is very large, like a 1000000 for example.

So I tried doing it by first calculating the distribution of one random number, then calculating the distribution of two random numbers from the first distribution. The first distribution is simply:


The next distribution can be calculated by integrating the previous distribution from x-0.5 to x+0.5:

Integrate[UnitBox[a], {a, x - 1/2, x + 1/2}]

I can't figure out how to iterate this to, for example, draw a plot from the result.

I also tried the following:

    Integrate[# /. x -> a, {a, x - 1/2, x + 1/2}, 
     Assumptions -> x \[Element] Reals], {x, -2, 2}], x] &, 
 UnitBox[x], 5]
Plot[%, {x, -2, 2}]

enter image description here

  • $\begingroup$ Any reason not to use built-ins? $\endgroup$
    – ciao
    Apr 19, 2015 at 3:26
  • $\begingroup$ No. What built-in would help? $\endgroup$
    – user
    Apr 19, 2015 at 3:27
  • 2
    $\begingroup$ maybe UniformSumDistribution. E.g., dis[n_] := UniformSumDistribution[n, {-.5, .5}]; Plot[Evaluate[PDF[dis@#, x] & /@ Range[10]], {x, -2, 2}]? $\endgroup$
    – kglr
    Apr 19, 2015 at 3:30
  • $\begingroup$ It works! Thanks a lot. $\endgroup$
    – user
    Apr 19, 2015 at 3:32
  • $\begingroup$ @MathematicaUser39386, posted the comment as an answer. $\endgroup$
    – kglr
    Apr 19, 2015 at 4:03

2 Answers 2

  • summing random numbers between -0.5 and 0.5.

  • the number of numbers to be summed is very large, like a 1000000 for example.

The sum of $n$ identical Uniform random variables is known as a generalised Irwin-Hall distribution, implemented in Mathematica as the UniformSumDistribution [ see kguler's answer]. The latter takes an $n$-part piecewise form, so, for example, if $n = 10$, the pdf is:

PDF[UniformSumDistribution[10, {-1/2, 1/2}], x]

enter image description here

The above is exact and works beautifully for small $n$. However, the OP poses the problem of $n$ being very large -- such as $n$ = 1 million -- which involves creating a 1 million part piecewise structure which is certainly going to fail. In fact, for $n = 1000$, simply making a plot of the pdf takes about 400 seconds on my Mac Pro:

(AA = Plot[ 
 PDF[UniformSumDistribution[1000, {-1/2, 1/2}], x], {x, -50, 50}, 
 PlotRange -> All]) // AbsoluteTiming 

390 seconds

  • For $n=10,000$, the above plot is simply unworkable ... which leaves the OP's question: how to proceed for large $n$??

Large $n$: apply Central Limit Theorem

For large $n$, I would suggest applying the (Lindeberg-Lévy) Central Limit Theorem:

  • If the random variables $(X_1, X_2, \dots)$ are independent and identically distributed, each with finite mean $\mu$ and finite variance $\sigma^2$, then the sample sum:

$$S_n \overset{a}{\sim} N\left(n \mu,n \sigma ^2\right)$$

In our case, $X \sim \text{Uniform}(-\frac12,\frac12)$ so $\mu = 0$ and $\sigma^2 = \frac{1}{12}$, so the asymptotic distribution of the sample sum is:

$$S_n \overset{a}{\sim} N\left(0, \frac{n}{12} \right)$$

with pdf $f(x)$:

       f = Exp[-6 (x^2/n)] / Sqrt[n Pi /6];
domain[f] = {x, -Infinity, Infinity} && {n > 0};

Here is a plot of pdf $f(x)$ when $n = 1000$:

BB = Plot[f /. n -> 1000, {x, -50, 50}, PlotRange -> All, PlotStyle -> Red]

This is, of course, instantaneous and works for arbitrarily large $n$. The following plot compares the exact solution AA to the asymptotic solution BB when $n = 1000$:

Show[BB, AA]

enter image description here

There is no discernible visual difference between the two plots here. Whereas the exact AA solution fails to evaluate for very large $n$, the asymptotic BB solution will always evaluate immediately, and with ever improving accuracy as $n$ increases.


The built-in function UniformSumDistribution may be useful:

usd[n_] := UniformSumDistribution[n, {-.5, .5}];
Plot[Evaluate[PDF[usd@#, x] & /@ Range[20]], {x, -2, 2}, 
 PlotStyle -> (ColorData[{"Rainbow", {1, 20}}] /@ Range[20]),
 Exclusions -> None, PlotRange -> All, ImageSize -> 500, 
 PlotLegends -> ("usd (" <> ToString[#] <> ")" & /@ Range[20])]

enter image description here

skd[n_] := SmoothKernelDistribution[RandomVariate[usd[n], 5000]];
Plot[Evaluate[PDF[#, x] & /@ {skd[5], usd[5], skd[3], usd[3]}], {x, -2, 2}, 
 ImageSize -> 500, Exclusions -> None, 
 Filling -> {1 -> {{2}, {Orange, Yellow}}, 3 -> {{4}, {Blue, Green}}}]

enter image description here


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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