Services & ResourcesWolfram ForumsMathGroup Archive

[Date Index] [Thread Index] [Author Index]

Re: Random number with custom distribution

• To: mathgroup at smc.vnet.net
• Subject: [mg107461] Re: [mg107445] Random number with custom distribution
• From: DrMajorBob <btreat1 at austin.rr.com>
• Date: Sat, 13 Feb 2010 05:23:19 -0500 (EST)
• References: <201002121008.FAA15765@smc.vnet.net>
• Reply-to: drmajorbob at yahoo.com

NOTE: there's a BUG indicated below, in RandomReal or Interpolation,  
probably the latter.

Here's an example PDF, first just a raw non-negative function:

fRaw[x_] =
   Abs@Piecewise[{{Cos@x, 0 <= x < Pi}, {Sin@x, Pi <= x <= 2 Pi}}];
Integrate[fRaw@x, {x, 0, 2 Pi}]

4

(Use NIntegrate, if Integrate doesn't work.)

To make it a PDF, scale it by dividing by the total integral (4 in this  
case):

pdf[x_] =
   Abs@Piecewise[{{Cos@x, 0 <= x < Pi}, {Sin@x, Pi <= x <= 2 Pi}}]/4;
pdfPlot =
  Plot[pdf@x, {x, 0, 2 Pi},
   PlotStyle -> Directive[Red, Thickness[.01]]]

Next, we need the CDF:

Clear[cdf]
cdf = cdf /.
    First@NDSolve[{cdf'[x] == pdf@x, cdf@0 == 0}, cdf, {x, 0, 2 Pi}];
cdfPlot = Plot[cdf@x, {x, 0, 2 Pi}]

In order to sample, we need the inverse CDF.

First, retrieve data points from the plot, reversing {x,y} in each  
datapoint:

data = First@Cases[cdfPlot, Line[x_] :> Reverse /@ x, Infinity];
ListPlot@data

That's a mirror image of the CDF plot. Interpolating gives an inverse to  
the CDF:

invCDF = Interpolation@data

InterpolatingFunction[{{3.20571*10^-8,1.}},<>]

That doesn't cover the whole range, unfortunately, so we should have made  
sure that "data" includes the points {0,0} and {2Pi,1}.

The current first and last points are

data[[{1, -1}]]

{{3.20571*10^-8, 1.28228*10^-7}, {1., 6.28319}}

We'll simply replace them with the end-points we want:

data[[{1, -1}]] = {{0, 0}, {1, 2 Pi}}

{{0, 0}, {1, 2 \[Pi]}}

and recompute the interpolation:

invCDF=Interpolation@data

InterpolatingFunction[{{0.,1.}},<>]

And now, random sampling is easy:

random[] := invCDF@RandomReal[]

Table[random[], {100}]

{0.173608, 0.135088, 6.11428, 4.90872, 4.30839, 5.06734, 2.80241, \
0.187845, 3.12592, 4.07651}

But there's a BUG lurking, since random[] is sometimes BADLY negative:

Select[Table[random[], {100000}], ! 0 <= # <= 2 Pi &]

{-0.0385181, -9.55779, -6.85895, -10.0228, -8.54523, -9.86493, \
-10.2442, -3.56191, -9.99411, -7.73028, -8.6684, -8.18295, -8.29155, \
-7.13628, -0.458468, -5.05335, -1.66142, -9.86214, -7.29216, \
-2.61625, -4.49397, -9.74907, -8.68753, -8.00137, -9.98978}

To fix this (and eliminate outputs above 2Pi, if any), redefine

Clear[random]
random[] :=
  Module[{try = invCDF@RandomReal[]},
   If[0 <= try <= 2 Pi, try, random[]]]

n = 10000;
Select[samples = Table[random[], {n}], ! 0 <= # <= 2 Pi &]
hist = Histogram[samples, {.01}, "ProbabilityDensity",
    ChartStyle -> Directive[Blue, EdgeForm[]]];

{}

Show[hist, pdfPlot]

A million samples makes the two plots match far better... but if your  
computer is slow, don't try it.

Bobby

On Fri, 12 Feb 2010 04:08:16 -0600, wiso <giurrero at gmail.com> wrote:

> Dear experts, is there a simple way to generate a random sample of a
> custom distribution? With Random I can generate random numbers
> according to the predefined distribution (Uniform, Normal, ...) but
> what about custom distribution (in my case a polynomial function)?
>


-- 
DrMajorBob at yahoo.com

• References:
  • Random number with custom distribution
    • From: wiso <giurrero@gmail.com>