# How could I improve this Central Limit Theorem demonstration speed- and content-wise?

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:

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 ?

-
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

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]


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]


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 *)

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