Article View: sci.math.stat
Article #563Re: normal random number generato
From: Naftaly Stramer
Date: Mon, 28 Sep 1992 23:23
Date: Mon, 28 Sep 1992 23:23
412 lines
14069 bytes
14069 bytes
In article <1992Sep28.183535.19303@cbfsb.cb.att.com>, rizzo@cbnewsf.cb.att.com (anthony.r.rizzo) writes:
>In article <1992Sep28.155713.6316@dazixco.ingr.com> naftaly@dazixco.ingr.com writes:
>>
>>Thank to all the neters who send me information on the subject. It really helped.
>>
>>Naftaly
>>
>
>How about SHARING all that information with the rest of us?
>
>I, for one, am quite interested in the subject too.
>
>Tony Rizzo
You are right. I am sorry I didn't do it early.
---
Naftaly Stramer (303)581-2370, FAX (303)581-9972 DAZIX, An Intergraph Company
Internet: naftaly@dazixco.ingr.com 6095 Spine Road
UUCP: ..uunet!ingr!dazixco!naftaly Boulder, CO 80301
Please excuse my English
From regeorge@magnus.acs.ohio-state.edu Mon Sep 21 14:33:21 1992
Return-Path: <regeorge@magnus.acs.ohio-state.edu>
Date: Mon, 21 Sep 92 16:33:54 -0400
From: Robert E George <regeorge@magnus.acs.ohio-state.edu>
To: naftaly@dazixco.ingr.com
Subject: Re: normal random number generator
Newsgroups: sci.math.stat
Organization: The Ohio State University
Cc:
Content-Length: 862
Status: RO
X-Lines: 23
I would suggest you use IMSL if it is available to you.
If not, try generating Unif(0,1) random variables and transforming
them into normal random variables -- almost any simulation book (Morgan or
Ripley for certain) can show you how.
But *be wary* of the uniform random number generator. If one comes included
with your system, be skeptical: find out how good it is before commiting
to it. (Morgan suggests some ways to do that.) If you need to build one from
scratch, that too can be done. Morgan will outline the process -- it
is based on modular arithmetic.
If you need further help, let me know . . .
Good luck,
Robert George
(speaking only for myself)
"The very essence of individual freedom is equal justice under a rule of
law, a law to which every man shall be subject and which no executive
can modify."
Senator Robert Taft , March 11, 1944
From T.Moore@massey.ac.nz Mon Sep 21 15:42:03 1992
Return-Path: <T.Moore@massey.ac.nz>
Date: Tue, 22 Sep 92 09:32:26 +1200
From: T.Moore@massey.ac.nz (Terry Moore)
To: naftaly@supergas.dazixco.ingr.com
Subject: Re: normal random number generator
Organization: Massey University
Content-Length: 914
Status: RO
X-Lines: 26
In article <1992Sep21.150254.19566@dazixco.ingr.com>, naftaly@supergas.dazixco.ingr.com (Naftaly Stramer) writes:
>
>
> Hi,
>
> I am looking for a routine (or the algorithem) that will generate random numbers,
> normally distribute.
>
There was a discussion of this a few weeks ago. Someone described
the polar method:
Let U1, U2 be uniform on (0, 1), let R = sqrt(-ln(U1)), theta = 2*pi*U2
then Z1 = R*sin(theta), and Z2 = R*cos(theta) are a pair of independent
random standard normals.
Someone else pointed out that this doesn't work as well as it might
if the uniform numbers are not as independent as they should be.
However, for a good uniform generator, I think the polar method
is quite good if you don't want too many.
See Knuth's "Art of Computer Programming" Vol?
for more information. He also describes a method due
to Marsaglia et al which is very fast for large sets of
numbers.
Terry Moore
From clint@leland.Stanford.EDU Mon Sep 21 16:43:51 1992
Return-Path: <clint@leland.Stanford.EDU>
Date: Mon, 21 Sep 92 15:43:57 PDT
From: clint@leland.stanford.edu
To: naftaly@dazixco.ingr.com
Subject: Re: normal random number generator
Newsgroups: sci.math.stat
Organization: DSG, Stanford University, CA 94305, USA
Cc:
Content-Length: 450
Status: RO
X-Lines: 8
That's indeed a FAQ. Usually you generate uniform deviates, using the
most common 2**31-1 multiplicative congruential algorithm. A good
reference is L'Ecuyer, Pierre, "Random Numbers for Simulation", Communications
of the ACM, October 1990, pp. 85-97. To generate normal deviates from
uniform, several approaches are possible. Some of the algorithms are
described in Abramovitz and Stegun, "Handbook of Mathematical Functions...".
Clint Cummins
From alan@mel.dms.csiro.au Mon Sep 21 17:41:05 1992
Return-Path: <alan@mel.dms.csiro.au>
Date: Tue, 22 Sep 92 09:41:20 EST
From: alan@mel.dms.csiro.au (Alan Miller)
To: naftaly@dazixco.ingr.com
Subject: Re: normal random number generator
Newsgroups: sci.math.stat
Organization: CSIRO, Div. of Mathematics & Statistics, Melbourne, Vic., AUSTRALIA
Cc:
Content-Length: 6473
X-Lines: 193
Status: RO
In article <1992Sep21.150254.19566@dazixco.ingr.com> you write:
>
>Hi,
>
>I am looking for a routine (or the algorithem) that will generate random numbers,
>normally distribute.
>
>thanks,
>
>Naftaly
>
>---
>Naftaly Stramer (303)581-2370, FAX (303)581-9972 DAZIX, An Intergraph Company
>Internet: naftaly@dazixco.ingr.com 5101 Spine Road
>UUCP: ..uunet!ingr!dazixco!naftaly Boulder, CO 80301
>Please excuse my English
Naftaly,
Here are two routines to generate random normals. The first is from ranlib,
which is available from statlib, and the second is my own.
Cheers
Alan Miller
REAL FUNCTION snorm()
C**********************************************************************C
C C
C C
C (STANDARD-) N O R M A L DISTRIBUTION C
C C
C C
C**********************************************************************C
C**********************************************************************C
C C
C FOR DETAILS SEE: C
C C
C AHRENS, J.H. AND DIETER, U. C
C EXTENSIONS OF FORSYTHE'S METHOD FOR RANDOM C
C SAMPLING FROM THE NORMAL DISTRIBUTION. C
C MATH. COMPUT., 27,124 (OCT. 1973), 927 - 937. C
C C
C ALL STATEMENT NUMBERS CORRESPOND TO THE STEPS OF ALGORITHM 'FL' C
C (M=5) IN THE ABOVE PAPER (SLIGHTLY MODIFIED IMPLEMENTATION) C
C C
C Modified by Barry W. Brown, Feb 3, 1988 to use RANF instead of C
C SUNIF. The argument IR thus goes away. C
C C
C**********************************************************************C
C
DIMENSION a(32),d(31),t(31),h(31)
C
C THE DEFINITIONS OF THE CONSTANTS A(K), D(K), T(K) AND
C H(K) ARE ACCORDING TO THE ABOVEMENTIONED ARTICLE
C
DATA a/0.0,.3917609E-1,.7841241E-1,.1177699,.1573107,.1970991,
+ .2372021,.2776904,.3186394,.3601299,.4022501,.4450965,
+ .4887764,.5334097,.5791322,.6260990,.6744898,.7245144,
+ .7764218,.8305109,.8871466,.9467818,1.009990,1.077516,
+ 1.150349,1.229859,1.318011,1.417797,1.534121,1.675940,
+ 1.862732,2.153875/
DATA d/5*0.0,.2636843,.2425085,.2255674,.2116342,.1999243,
+ .1899108,.1812252,.1736014,.1668419,.1607967,.1553497,
+ .1504094,.1459026,.1417700,.1379632,.1344418,.1311722,
+ .1281260,.1252791,.1226109,.1201036,.1177417,.1155119,
+ .1134023,.1114027,.1095039/
DATA t/.7673828E-3,.2306870E-2,.3860618E-2,.5438454E-2,
+ .7050699E-2,.8708396E-2,.1042357E-1,.1220953E-1,.1408125E-1,
+ .1605579E-1,.1815290E-1,.2039573E-1,.2281177E-1,.2543407E-1,
+ .2830296E-1,.3146822E-1,.3499233E-1,.3895483E-1,.4345878E-1,
+ .4864035E-1,.5468334E-1,.6184222E-1,.7047983E-1,.8113195E-1,
+ .9462444E-1,.1123001,.1364980,.1716886,.2276241,.3304980,
+ .5847031/
DATA h/.3920617E-1,.3932705E-1,.3950999E-1,.3975703E-1,
+ .4007093E-1,.4045533E-1,.4091481E-1,.4145507E-1,.4208311E-1,
+ .4280748E-1,.4363863E-1,.4458932E-1,.4567523E-1,.4691571E-1,
+ .4833487E-1,.4996298E-1,.5183859E-1,.5401138E-1,.5654656E-1,
+ .5953130E-1,.6308489E-1,.6737503E-1,.7264544E-1,.7926471E-1,
+ .8781922E-1,.9930398E-1,.1155599,.1404344,.1836142,.2790016,
+ .7010474/
C
10 u = ranf()
s = 0.0
IF (u.GT.0.5) s = 1.0
u = u + u - s
20 u = 32.0*u
i = int(u)
IF (i.EQ.32) i = 31
IF (i.EQ.0) GO TO 100
C
C START CENTER
C
30 ustar = u - float(i)
aa = a(i)
40 IF (ustar.LE.t(i)) GO TO 60
w = (ustar-t(i))*h(i)
C
C EXIT (BOTH CASES)
C
50 y = aa + w
snorm = y
IF (s.EQ.1.0) snorm = -y
RETURN
C
C CENTER CONTINUED
C
60 u = ranf()
w = u* (a(i+1)-aa)
tt = (0.5*w+aa)*w
GO TO 80
70 tt = u
ustar = ranf()
80 IF (ustar.GT.tt) GO TO 50
90 u = ranf()
IF (ustar.GE.u) GO TO 70
ustar = ranf()
GO TO 40
C
C START TAIL
C
100 i = 6
aa = a(32)
GO TO 120
110 aa = aa + d(i)
i = i + 1
120 u = u + u
IF (u.LT.1.0) GO TO 110
130 u = u - 1.0
140 w = u*d(i)
tt = (0.5*w+aa)*w
GO TO 160
150 tt = u
160 ustar = ranf()
IF (ustar.GT.tt) GO TO 50
170 u = ranf()
IF (ustar.GE.u) GO TO 150
u = ranf()
GO TO 140
END
c--------------------------------------------------------------------------
real function rnorm(second)
c
c Generate a random normal deviate using the polar method.
c Reference: Marsaglia,G. & Bray,T.A. 'A convenient method for
c generating normal variables', Siam Rev., vol.6,
c 260-264, 1964.
c
c N.B. Set second = .false. the first time that this function
c is called.
c
logical second
c
c Local variables
c
real one, u, v, sum, sln
real rand
save v, sln
data one /1.0/
c
if(second) go to 20
c
c First call; generate a pair of random normals
c
second = .true.
10 u = rand()
v = rand()
u = u + u - one
v = v + v - one
sum = u*u + v*v
if(sum .ge. one) go to 10
sln = - log(sum)
sln = sqrt((sln + sln)/sum)
rnorm = u*sln
return
c
c If second, use the second random number generated on last call
c
20 second = .false.
rnorm = v*sln
return
end
--
Alan Miller, Quality Improvement Project
CSIRO Division of Mathematics & Statistics, Melbourne, Australia
Phone: +61 3 542-2266 Fax: +61 3 542-2474 E-mail: alan@mel.dms.csiro.au
Mail: CSIRO DMS, Private Bag 10, Clayton, Vic. 3168, Australia
From mzf@castle.edinburgh.ac.uk Tue Sep 22 08:11:21 1992
Return-Path: <mzf@castle.edinburgh.ac.uk>
From: M Firat <mzf@castle.edinburgh.ac.uk>
Subject: Re: normal random number generator
To: naftaly@dazixco.ingr.com
Date: Tue, 22 Sep 92 15:11:07 WET DST
Content-Length: 677
X-Lines: 28
Status: RO
> Hi,
>
> I am looking for a routine (or the algorithem) that will generate random numbers,
> normally distribute.
>
> thanks,
>
> Naftaly
>
> ---
> Naftaly Stramer (303)581-2370, FAX (303)581-9972 DAZIX, An Intergraph Company
> Internet: naftaly@dazixco.ingr.com 5101 Spine Road
> UUCP: ..uunet!ingr!dazixco!naftaly Boulder, CO 80301
> Please excuse my English
Fortran 77 NAG library routine GO5CBF generate normally distributed
random variables. If you have any problems please do not hasitate to
contact me.
Good luck.
Mehmet
Dept. of Math. and stats.
Univ. of Edinburgh
JCMB
Edinburg EH9 3JZ
Scotland, UK
Tel : (031) 650 8575
From whart@cs.UCSD.EDU Tue Sep 22 09:11:59 1992
Return-Path: <whart@cs.UCSD.EDU>
Date: Tue, 22 Sep 92 08:12:10 PDT
From: whart@cs.UCSD.EDU (Bill Hart)
To: naftaly@dazixco.ingr.com
Subject: Re: normal random number generator
Newsgroups: sci.math.stat
Content-Length: 278
X-Lines: 12
Status: RO
In sci.math.stat you write:
>Hi,
>I am looking for a routine (or the algorithem) that will generate random numbers,
>normally distribute.
Check out the software in Numerical Recipies in C by Press et. al. There
is a routine for doing just that in this collection.
--Bill
From miller@cogsci.UCSD.EDU Tue Sep 22 13:44:02 1992
Return-Path: <miller@cogsci.UCSD.EDU>
Date: Tue, 22 Sep 92 12:44:25 PDT
From: miller@cogsci.UCSD.EDU (Jeff Miller)
To: naftaly@dazixco.ingr.com
Subject: Re: normal random number generator
Newsgroups: sci.math.stat
Organization: University of California, San Diego
Cc:
Content-Length: 74
Status: RO
X-Lines: 2
Try "Numerical Recipes: The art of scientific computing."
by Press et al.
Naftaly Stramer (303)581-2370, FAX (303)581-9972 DAZIX, An Intergraph Company
Internet: naftaly@dazixco.ingr.com 6095 Spine Road
UUCP: ..uunet!ingr!dazixco!naftaly Boulder, CO 80301
Please excuse my English
Technical Details
Message-ID:
<1992Sep28.232300.8559@dazixco.ingr.com>Path:
arc1-ca.pugleaf.net!sparky!uunet!infonode!ingr!dazixcon!naftalyReferences:
<1992Sep28.155713.6316@dazixco.ingr.com> <1992Sep28.183535.19303@cbfsb.cb.att.com>Show Full Headers
Complete Article Headers
Newsgroups: sci.math.stat Path: arc1-ca.pugleaf.net!sparky!uunet!infonode!ingr!dazixcon!naftaly From: naftaly@supergas.dazixco.ingr.com (Naftaly Stramer) Subject: Re: normal random number generato Message-ID: <1992Sep28.232300.8559@dazixco.ingr.com> Sender: naftaly@supergas (Naftaly Stramer) Nntp-Posting-Host: supergas Reply-To: naftaly@dazixco.ingr.com Organization: Dazix References: <1992Sep28.155713.6316@dazixco.ingr.com> <1992Sep28.183535.19303@cbfsb.cb.att.com> Date: Mon, 28 Sep 1992 23:23:00 GMT Lines: 412