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