Simulating a Normal Distribution

Andy Long
Spring, 2020

Suppose we want to simulate draws from a normal distribution:

Let’s say with mean 2 and standard deviation 3, and we want to generate 1000 values from this distribution:

mu=2;
sigma=3;
samplesize=1000;
(* we'll look at stuff between plus and minus 4 standard deviations from the mean: *)
leftint=mu-4*sigma;
rightint=mu+4*sigma;

Just so you know, here are the Probability Density Function and Cumulative Distribution Function for this Normal (the cumulative is just the integral of the pdf, from -infinity to x):

PDF[NormalDistribution[mu,sigma],x]
Plot[%,{x,leftint,rightint},Filling->Axis]
CDF[NormalDistribution[mu,sigma],x]
Plot[%,{x,leftint,rightint},Filling->Axis]

Well, Mathematica makes that easy, with RandomVariate, so this is the thing that we’re going to do “the hard way”:

randomData=RandomVariate[NormalDistribution[mu,sigma],samplesize];
Show[
Histogram[randomData,20,"ProbabilityDensity"],
Plot[PDF[NormalDistribution[mu,sigma],x],{x,leftint,rightint},PlotStyle->Thick]]

We could do this ourselves, however, for any distribution.  I’m just going to illustrate it here with this normal.

We start with a uniform probability density function (pdf), and pull a samplesize of random values from it:

PDF[UniformDistribution[{0,1}],x]
Plot[%,{x,-1,2},Filling->Axis]
us=RandomVariate[UniformDistribution[{0,1}],samplesize];

Now we solve 
u=normalcdf(0,1,z) for z, 
and then we have our x from the z: 
x = mu + sigma*z

zs=Table[Solve[CDF[NormalDistribution[0,1],z]==us[[i]]],{i,1,samplesize}];
zs[[1]]
zs=Table[z/.zs[[i]][[1]],{i,1,samplesize}];
zs[[1]]
xs=mu + sigma * zs;
xs[[1]]

Show[Histogram[xs,20,"ProbabilityDensity"],
Plot[PDF[NormalDistribution[mu,sigma],x],{x,leftint,rightint},PlotStyle->Thick]]