From: rusin@olympus.math.niu.edu (Dave Rusin)
Newsgroups: sci.math.num-analysis
Subject: Re: random generator from user distribution
Date: 24 Aug 1995 14:34:54 GMT

In article <41dqgf$kbe@rover.ucs.ualberta.ca>,
Steve Dew <sdew@ee.ualberta.ca> wrote:
>I'm looking for references/programs which can accurately but efficiently
>generate a sequence of random numbers from an arbitrary user-supplied
>distribution (specified in the form of a sequence of x,y pairs).

I don't know if the plotted  (x,y) pairs are supposed to be the
probability density function or the cumulative probability (that is,
as  x--> oo, does  y-->0 or  y --> 1 ?), but I'll suppose the latter.
(If not, you'll need to pick an interpolation between the points and
integrate. If you pick linear interpolation, this is a trivial task).
Then invert this function, that is, let  g  be the function
g: [0,1] --> R, a table of whose values is the set of pairs  (y,x).

In that case, you need only generate values  y  in [0,1]  at random
and compute  x=g(y); the  x  values will have the desired distribution.
Of course you still run into the interpolation question, but I would
think a linear scheme would be sufficient for a random number 
generator if you have enough data points. (To be consistent with a
linear interpolation of a probability density function, I guess you'd
have to use a quadratic approximation to  g  with continuous first
derivative.)

Here's the underlying theory: if  f(t)  is the probability density function,
the integral C(x)=\int_0^x f(t)dt is the cumulative distribution function.
I am asking you to define  g  as the inverse to  C and then generate the
random numbers by taking  x = g(y)  where the inputs  y  are randomly
generated in [0,1] with uniform density. Indeed, given an interval
(a,b),  the probability that the number  x  so generated lies in the
interval is the probability that the chosen  y  lies in the interval
(g^-1(a), g^-1(b) ), i.e., (C(a), C(b)). Since  y  is generated randomly
with a uniform distribution, this probability is just the length
(=C(b)-C(a) ) of the last interval. For  b  close to  a, this is roughly
C'(a).(b-a), i.e., f(a).(b-a). But that's what it means to say that
f  is the probability density function.

dave

==============================================================================