Selecting random numbers from a specified distribution

The default setting for Random is a uniform distribution on the interval {0, 1}. There are a number of built-in ways to modify this, such as specifying the type of number you want or to specify a domain different from {0, 1}. You can also use the Statistics packages and specify any of the built-in distributions from which to select the random numbers. There are occasions when you need to select random numbers from a distribution which is not built in to Mathematica. This notebook illustrates how to do that.

A density function

Here is a density function to illustrate the process. It is based on a half period of the sine function with support interval {h, h+ π/B}.

f2[t_,B_,h_] :=  B UnitStep[(t-h)(h+ π/B-t)]Sin[B(t-h)]/2

A check using plot

Plot[f2[t,5,2],{t,1,4}]

A check using  Integrate

Integrate[f2[t,2.5,3],{t,1,6}]//Simplify

A cumulative distribution function

In order to select random numbers from this distribution we need first to find the distribution function. You can check that the following is the correct distribution function for f.

F[x_,B_,h_] :=(1/2) (-(-1+Cos[B (h-x)]) UnitStep[(h+(π/B)-x) (-h+x)]+2 UnitStep[-h-(π/B)+x])

A check using Plot

Plot[F[x,3,2],{x,1,4}]

The next step

The next step in the process is to find the inverse of the distribution function. If you have used UnitStep to define your density and distribution functions you will need to extract the relevant parts and use them in Solve.

sol = Solve[(1/2) (-(-1+Cos[B (h-x)]))==y,x]

We need the second solution here as we will see. This will always be defined on the interval {0, 1} thus the UnitStep as a factor.

Inv[y_,B_,h_] := ((B h+ArcCos[1-2 y])/B)UnitStep[y (1-y)]

A check of our inverse using Plot.

Plot[Inv[y,3,2],{y,0,1}]

Defining our SinRandom function

All that is left is to define our random function.

SinRandom[B_,h_] :=  Inv[Random[],B,h];

A check of our SinRandom function

To illustrate we have done what we claim we will need some functions from the Statistics`DataManipulation` package. So we load it now.

<<Statistics`

pts = Table[SinRandom[3,2],{t,100000}];

pts[[1]]

Next we get the frequencies of the values in our list.

Clear[freqs]

freqs = 18.72 BinCounts[pts,{2.,2.+ π/3,π/60}]/100000;

bins = Range[2.+ π/120,2.+ π/3-π/120,π/60];

A graphical check

We can see that this is a reasonable set of points from the given distribution by plotting the frequencies of the points against the density function after some appropriate scaling. There will always be some variation from the curve but if you run these last few cells a number of times you will see that the samples follow the distribution closely.

Plot[f2[x,3,2],{x,1,4},Epilog-> {Hue[.95],AbsolutePointSize[3], Point/@Transpose[{bins, freqs}]}]