MathGroup Archive 1998

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

Search the Archive

Re: random number generation


  • To: mathgroup@smc.vnet.net
  • Subject: [mg11205] Re: random number generation
  • From: Daniel Lichtblau <danl@wolfram.com>
  • Date: Mon, 2 Mar 1998 23:10:58 -0500
  • Organization: Wolfram Research, Inc.
  • References: <6cnabv$bha@smc.vnet.net>

Christopher Moss wrote:
> 
> Hello,
> I'm a college student and have been given a problem where I must
> generate random numbers using both the lognormal distribution and the
> gumbel (also known as extreme value type 1) distribution.  The problem
> is that the old school professor won't let the class use any built-in
> functions in software packages to do this.  I know Mathematica has this
> and it would be very easy to use :-(.  Does anyone out there have
> "programs" written in Mathematica to generate random numbers.  I can't
> imagine why anyone would go through the exercise of doing this but I
> thought I would take a shot.
> 
> Thankyou,
> ___________________________________________ Christopher Moss
> "The DenKeeper"
> http://home.earthlink.net/~mosser/denkeeper.html
> ___________________________________________


I really do not know much about the subject; apologies to the experts if
I mangle the terminology. A nice reference that will get the
terminology straight and explain the relevant idea is Ch. 5 (sections 2
and 4) of "Introduction to Probability Theory" by Hoel, Port, and
Stone. There are no doubt dozens of more recent references, this
happens to be the one I have around.

There is a standard way to generate elements in a given continuous
probability distribution, given access to a uniform distribution random
generator. I'll illustrate with the normal distribution
(1/Sqrt[Pi])*Exp[-x^2]. First note that in Mathematica the function
call Random[] gives a pseudo-random double precision value uniformly
distributed in the range from zero to one.

The idea can be pictured graphically as follows. You have a continuous
distribution function f, so the range is zero to one. Graph its inverse
(call this function g) whose domain is now zero to one. A value p
between zero and one corresponds to g(p). This is because g(p) is the
value for which there is cumulative probability (that is, area between
graph and x-axis) p to the left in the graph of the distribution
function f.

How to go about coding this? The key is to notice that if we generate
.5, say, from calling Random[], then this should correspond to the
value in the normal distribution that gives probability .5, in other
words, 0. More generally, if we have a between zero and one, we need to
find the value t such that the probability of getting a in the normal
distribution is t. The equation to set up is
1/Sqrt[Pi]*Integrate[Exp[-x^2],{x,-Infinity,t}]==a.

In[23]:= 1/ 1/Sqrt[Pi]*Integrate[Exp[-x^2],{x,-Infinity,t}]//InputForm

Out[23]= If[t < 0, (Sqrt[Pi]*(1 + Erf[t]))/2, Integrate[E^(-x^2), {x,
-Infinity, t}]]/ Sqrt[Pi]

This is just undue pessimism; the integral is also correct for t>=0. 

In[25]:=  1/Sqrt[Pi]*Integrate[Exp[-x^2],
	{x,-Infinity,t},Assumptions->t>0]
         1 + Erf[t]
Out[25]= ----------
             2

In[28]:= ee = %25;


Say we had obtained a=.7 from calling Random[]. In the normal
distribution this would correspond to the value

In[29]:= FindRoot[ee==.7, {t,0}]
Out[29]= {t -> 0.370807}

Let's check this result.

In[30]:= ee /. %
Out[30]= 0.7


Since this is homework I'll stop here. One thing I will add is that it
may be possible to gain a bit of speed by replacing the distribution
function of interest with an InterpolatingFunction approximation.


Daniel Lichtblau
Wolfram Research



  • Prev by Date: Re: lattice definition: help
  • Next by Date: RE: Efficient use of coefficient--Efficient simplification
  • Prev by thread: Re: Importing .pix files .
  • Next by thread: init.m on IBM