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.

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 ?

share|improve this question
    
On my machine, it is 3 to 4 times faster to write RandomVariate[dist, {10, k}] than Table[RandomVariate[dist,10],{k}]. –  KennyColnago Oct 30 '13 at 20:01
    
@MichaelE2 thanks, fixed. I was also passing the plot range, but I removed it for brevity. –  Zet Oct 30 '13 at 20:05

1 Answer 1

up vote 6 down vote accepted

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 *)
share|improve this answer
    
@Zet I added a compiled version that is much faster (on my machine). –  Michael E2 Oct 30 '13 at 20:11
    
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[] ? –  Zet Oct 30 '13 at 20:16
    
@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. –  Michael E2 Oct 30 '13 at 21:51

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.