4
$\begingroup$

My purpose is to demonstrate the Central Limit theorem by assuming two custom distributions, taking some random samples, calculating their means and ploting the results. Ideally, the histograms should tend to a normal distribution.

The code I am using is the following:

\[ScriptCapitalD]1 = ProbabilityDistribution[
   1/5 Exp[-(1/5) Abs[2 x - 3]],
   {x, -∞, +∞}];

\[ScriptCapitalD]2 = ProbabilityDistribution[
    If[x >= -1/2 && x <= +1/2, 1, 0],
    {x, -∞, +∞}];

f[n_] := ParallelTable[
    Histogram[
       Mean /@ Table[RandomVariate[n, 10], {k}], PlotLabel -> "n=" <> ToString@k],
    {k, 10, 170, 40}]

cltHistPlots = GraphicsGrid[
    {f[\[ScriptCapitalD]1], f[\[ScriptCapitalD]2]},
    Spacings -> 0]

And the graph that is generated is this:

enter image description here

My questions are:

  1. How could I speed the process up ? The code is very slow! I think the bottleneck is in the part where RandomVariate[] calculates the samples.

  2. A bit offtopic, but, could you think of any other addition I could make so that the audience would better understand the theorem ?

$\endgroup$
2
  • $\begingroup$ On my machine, it is 3 to 4 times faster to write RandomVariate[dist, {10, k}] than Table[RandomVariate[dist,10],{k}]. $\endgroup$ Commented Oct 30, 2013 at 20:01
  • $\begingroup$ @MichaelE2 thanks, fixed. I was also passing the plot range, but I removed it for brevity. $\endgroup$
    – stathisk
    Commented Oct 30, 2013 at 20:05

1 Answer 1

6
$\begingroup$

You could use InverseCDF and map a uniform RandomReal distribution onto yours. Note: this assumes that the inverse cdf results in something easy to calculate. For demonstration purposes, this should be easily achieved.

dist1 = ProbabilityDistribution[
   1/5 Exp[-(1/5) Abs[2 x - 3]], {x, -\[Infinity], +\[Infinity]}];

dist2 = ProbabilityDistribution[
    If[x >= -1/2 && x <= +1/2, 1, 
     0], {x, -\[Infinity], +\[Infinity]}];;

invcdf1 = Function[{q}, Evaluate@InverseCDF[dist1, q], Listable]
invcdf2 = Function[{q}, Evaluate@InverseCDF[dist2, q], Listable]

Mathematica graphics

f[icdf_] := 
 ParallelTable[
  Histogram[Mean /@ Table[icdf[RandomReal[1, k]], {k}], 
   PlotLabel -> "n=" <> ToString@k], {k, 10, 170, 40}]

cltHistPlots = GraphicsGrid[{f[invcdf1], f[invcdf2]}, Spacings -> 0]

Mathematica graphics

I suppose you might overlay a plot of the theoretical normal distribution for comparison purposes.


If you're willing to be a little risky, you could use the following compiled version. It assumes that neither 0 or 1 will be fed to the inverse cdfs.

invcdf1 = 
 Compile[{q}, Evaluate@Simplify[InverseCDF[dist1, q], 0 < q < 1], 
  RuntimeAttributes -> {Listable}, Parallelization -> True];
invcdf2 = 
 Compile[{q}, Evaluate@Simplify[InverseCDF[dist2, q], 0 < q < 1], 
  RuntimeAttributes -> {Listable}, Parallelization -> True];

GraphicsGrid[{f[invcdf1], f[invcdf2]}, Spacings -> 0] // AbsoluteTiming // First
(* 0.083894 *)
$\endgroup$
3
  • $\begingroup$ @Zet I added a compiled version that is much faster (on my machine). $\endgroup$
    – Michael E2
    Commented Oct 30, 2013 at 20:11
  • $\begingroup$ Wow @MichaelE2, just wow! Even your first solution was fast enough for me, let alone the compiled version. I'm so happy, I will accept your answer right away :D How did you know BTW that using the inverse transform sampling would be so much faster than RandomVariate[] ? $\endgroup$
    – stathisk
    Commented Oct 30, 2013 at 20:16
  • $\begingroup$ @Zet I thought it worth a try to see if Mathematica would come up with a closed form for the inverse cdf. If it did, it would be fast. If not, then it would probably use some numerical algorithm that would not be so fast. In some sense, we were lucky. RandomVariate perhaps uses the numerical algorithm in all cases, but I'm not sure. $\endgroup$
    – Michael E2
    Commented Oct 30, 2013 at 21:51

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.