Services & ResourcesWolfram ForumsMathGroup Archive

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

Re: defining random variable with piecewise-continuous

• To: mathgroup at smc.vnet.net
• Subject: [mg79522] Re: [mg79462] defining random variable with piecewise-continuous
• From: Darren Glosemeyer <darreng at wolfram.com>
• Date: Sat, 28 Jul 2007 05:38:15 -0400 (EDT)
• References: <200707270953.FAA02959@smc.vnet.net>

sdw wrote:
> any ideas on how to do this?
>
> basically seeking a computationally efficient way to define the PDF, which 
> may look like a set of continuous or discontinuous rectangle functions and 
> have the Mathematica random number generator produce values from this PDF
>
> thanks 
>
>   

Here is function that will generate numbers from a piecewise rectangular 
density, which sounds like what you are looking for. The function takes 
as arguments a list of lower bounds for the intervals, a list of upper 
bounds, a list of density values for the intervals, and the number of 
random numbers desired.


In[1]:= rand = Compile[{{lower, _Real, 1}, {upper, _Real, 1}, {height, 
_Real,
            1}, {n, _Integer}},
          Module[{rc, mult = (upper - lower), rr = RandomReal[1., n]},
           rc = RandomChoice[height*mult -> Range[Length[height]], n];
           lower[[rc]] + rr*mult[[rc]]]];

The code works by randomly selecting regions based on the probabilities 
of following in each region and then scaling random uniforms into the 
selected regions.

So for instance, if the pdf is the following

In[2]:= pdf[x_] :=
         Piecewise[{{.2, 0 <= x <= 2}, {.1, 4 <= x < 6}, {.4, 6 <= x <= 7}}]


ten numbers from this distribution can be generated with this input.

In[3]:= rand[{0., 4., 6.}, {2., 6., 7.}, {.2, .1, .4}, 10]

Out[3]= {6.4957, 1.16766, 5.80021, 1.27402, 6.50937, 6.70427, 0.728181, 
1.25066, 6.49188,
 
 >    6.92062}


Histograms can be used to compare the random numbers with the 
theoretical density.

In[4]:= <<Histograms`

In[5]:= Show[Histogram[rand[{0., 4., 6.}, {2., 6., 7.}, {.2, .1, .4}, 
10^4],
  HistogramScale -> 1], Plot[pdf[x], {x, 0, 7}, PlotStyle -> {Thick}]]

In[6]:= Show[Histogram[rand[{0., 4., 6.}, {2., 6., 7.}, {.2, .1, .4}, 
10^4],
          HistogramScale -> 1, HistogramCategories -> {0, 2, 4, 6, 7}],
         Plot[pdf[x], {x, 0, 7}, PlotStyle -> {Thick}]]


Note that for more general piecewise densities efficient methods may 
vary greatly depending on the shapes of the pieces.

Darren Glosemeyer
Wolfram Research

• References:
  • defining random variable with piecewise-continuous PDF Mathematica 6
    • From: "sdw" <warwick@jps.net>