




Fortran Programming Language









Well cited Fortran Gaussian PRNG?
Hi group, Does anyone have recommendations for a well cited Fortran Gaussian pseudorandom number generator? I currently use the following uniform PRNGs based on the Mersenne Twister Algorithm: http://www.math.sci.hiroshimau.ac.jp/~mmat/MT/VERSIONS/FORTRAN/TAKA... I know it is possible to transform a uniform distribution to a Gaussian one through the BoxMuller or Ziggurat techniques but it will need to be tested, and testing of PRNGs is not a trivial exercise, AIUI. Thanks everyone :) skate xx
On Apr 13, 5:57 am, "sk8terg1rl" <sk8terg1rl_2@yahoo.co.uk> wrote: > Hi group, > Does anyone have recommendations for a well cited Fortran Gaussian > pseudorandom number generator?
There are some at statistician Alan Miller's site http://users.bigpond.net.au/amiller/random.html . The codes cite the sources of the algorithms. If you want to ensure that your results are not distorted by a poor Gaussian RNG, I suggest putting the code for the RNG in a separate module, say random_normal_mod, so that you can later replace the module with one using another algorithm and see if the results are effectively the same.
sk8terg1rl wrote: > Does anyone have recommendations for a well cited Fortran Gaussian > pseudorandom number generator? > I currently use the following uniform PRNGs based on the Mersenne > Twister Algorithm: > http://www.math.sci.hiroshimau.ac.jp/~mmat/MT/VERSIONS/FORTRAN/TAKA... > I know it is possible to transform a uniform distribution to a > Gaussian one through the BoxMuller or Ziggurat techniques but it will > need to be tested, and testing of PRNGs is not a trivial exercise, It is, I believe, difficult to get a Gaussian RNG to have a proper tail. To have the tail make sense you really need a lot of samples. With BoxMuller it isn't too hard to figure out the tail distribution if you know the resolution of the uniform PRNG behind it. http://mathworld.wolfram.com/BoxMullerTransformation.html Also, this reminds me again of the usefulness of cosd() and sind(), that is, in degrees. That avoids the problem of using an approximation to pi different than the one used by sin().  glen
On 20070413 11:09:58 0300, glen herrmannsfeldt <g@ugcs.caltech.edu> said: I
>> know it is possible to transform a uniform distribution to a >> Gaussian one through the BoxMuller or Ziggurat techniques but it will >> need to be tested, and testing of PRNGs is not a trivial exercise, > It is, I believe, difficult to get a Gaussian RNG to have a proper > tail. To have the tail make sense you really need a lot of samples. > With BoxMuller it isn't too hard to figure out the tail distribution > if you know the resolution of the uniform PRNG behind it. > http://mathworld.wolfram.com/BoxMullerTransformation.html > Also, this reminds me again of the usefulness of cosd() and sind(), > that is, in degrees. That avoids the problem of using an approximation > to pi different than the one used by sin(). >  glen
At the risk of sounding like a broken record it would make a lot of sense to have alternate periodic functions with period EXACTLY one. In introductory calculus one defines sine to have a derivative at zero of one. Then with some amount of derivation one can calculate the period. The result is radian measure of angles. In a first graduate course in harmonic analysis one defines the periodic character function to have period of one. Then with some amount of derivation one can calculate the derivative at zero. The result is cyclic measure of angles. The magic value that shows up in both schemes is two pi. If one goes and buys a frequency meter the result will by displayed in cycles. All of the instances that arise in practice of the problems of angle reduction requiring an accurate value of two pi are symptoms of using radian measure when cyclic measure would be more appropriate. It is interesting that physical instruments act like they have taken the graduate course in harmonic analysis. Perhaps a few more "analysts" should learn to do likewise. It is really a rather trivial exercise in where one places the two pi. Sometimes the analysts even know that this is more sensible but they have never been given "permission" to do the sensible thing. They are reading instruments in cycles but trying to calculate in radians.
sk8terg1rl wrote: > Does anyone have recommendations for a well cited Fortran Gaussian > pseudorandom number generator? > ... > I know it is possible to transform a uniform distribution to a > Gaussian one through the BoxMuller or Ziggurat techniques...
Start with the LAPACK {SDCZ}LARNV routine. Distribution option #3 uses BoxMuller. Walter
On Apr 13, 8:27 am, Gordon Sande <g.sa@worldnet.att.net> wrote: ... > [...] All of the instances that arise in practice of the > problems of angle reduction requiring an accurate value > of two pi are symptoms of using radian measure when cyclic > measure would be more appropriate. [...]
Well, I guess that depends on your definition of "all". Suppose you have a value of 30 degrees. The sine of 30 degrees is exactly onehalf. In radians, 30 degrees is pi/6. It's problematical whether the calculated value of SIN(pi/6) will round to exactly onehalf or not. In cycles, 30 degrees is 1/12. Again, it's problematical whether your cyclebased sine function will give exactly onehalf or not. In both cases the problem is that the angle in question is not exactly representable in floating point using the angle measure chosen. For degrees, the angle is an integer and is exactly representable in all practical floatingpoint kinds. (Similarly for 60 degrees, and so on.) A similar problem exists for radians when the angle is 45 degrees (that's an eighth of a cycle, so your cycle based measure would work for binary floatingpoint  and, with proper care, even in decimal floatingpoint). Does computing the tangent of pi/4 radians exactly yield one? Maybe, for some floatingpoint KINDs. But can you always count on it? If all angles that were an integral multiple of 1/24 of a cycle were exactly representable in floatingpoint these identities could all be met. Degrees meets that requirement: 1/24 of a cycle is 15 degrees. And it has the advantage of being a measure of angles in common usage. I've always been surprised the engineering community didn't demand degree based trig functions.  J. Giles "I conclude that there are two ways of constructing a software design: One way is to make it so simple that there are obviously no deficiencies and the other way is to make it so complicated that there are no obvious deficiencies."  C. A. R. Hoare
On 20070413 13:03:08 0300, "jamesgi@att.net" <jamesgi@att.net> said:
> On Apr 13, 8:27 am, Gordon Sande <g.sa @worldnet.att.net> wrote: > ... >> [...] All of the instances that arise in practice of the >> problems of angle reduction requiring an accurate value >> of two pi are symptoms of using radian measure when cyclic >> measure would be more appropriate. [...] > Well, I guess that depends on your definition of "all". > Suppose you have a value of 30 degrees. The sine of > 30 degrees is exactly onehalf. In radians, 30 degrees > is pi/6. It's problematical whether the calculated > value of SIN(pi/6) will round to exactly onehalf or not. > In cycles, 30 degrees is 1/12. Again, it's problematical > whether your cyclebased sine function will give exactly > onehalf or not. In both cases the problem is that the > angle in question is not exactly representable in floating > point using the angle measure chosen. For degrees, the angle > is an integer and is exactly representable in all practical > floatingpoint kinds. (Similarly for 60 degrees, and so on.) > A similar problem exists for radians when the angle is > 45 degrees (that's an eighth of a cycle, so your cycle > based measure would work for binary floatingpoint  > and, with proper care, even in decimal floatingpoint). > Does computing the tangent of pi/4 radians exactly yield > one? Maybe, for some floatingpoint KINDs. But can you > always count on it? > If all angles that were an integral multiple of 1/24 of > a cycle were exactly representable in floatingpoint > these identities could all be met. Degrees meets that > requirement: 1/24 of a cycle is 15 degrees. And it has > the advantage of being a measure of angles in common > usage. I've always been surprised the engineering community > didn't demand degree based trig functions.
There are two problems. Mixing them together does not help. One is angle reduction which is what I mentioned and you have ignored. The other is representation in finite precision. You have done a marvelous job of arguing for a large composite base, like 60, and cyclic measure. Use of degrees is a pragmatic and commonly used approximation to that. Physics repeatedly gets into angle reduction problems. I have not seen any of physics where degrees arise naturally instead of cycles. The places where degrees are used do not seem to lead to much need of angle reduction which would be quite easy if it were required.
In a previous article, Gordon Sande <g.sa@worldnet.att.net> wrote:
>On 20070413 13:03:08 0300, "jamesgi @att.net" <jamesgi @att.net> said: >> On Apr 13, 8:27 am, Gordon Sande <g.sa@worldnet.att.net> wrote: >> ... >>> [...] All of the instances that arise in practice of the >>> problems of angle reduction requiring an accurate value >>> of two pi are symptoms of using radian measure when cyclic >>> measure would be more appropriate. [...] >> Well, I guess that depends on your definition of "all". >> Suppose you have a value of 30 degrees. The sine of >> 30 degrees is exactly onehalf. In radians, 30 degrees >> is pi/6. It's problematical whether the calculated >> value of SIN(pi/6) will round to exactly onehalf or not. >> In cycles, 30 degrees is 1/12. Again, it's problematical >> whether your cyclebased sine function will give exactly >> onehalf or not. In both cases the problem is that the >> angle in question is not exactly representable in floating >> point using the angle measure chosen. For degrees, the angle >> is an integer and is exactly representable in all practical >> floatingpoint kinds. (Similarly for 60 degrees, and so on.) >> A similar problem exists for radians when the angle is >> 45 degrees (that's an eighth of a cycle, so your cycle >> based measure would work for binary floatingpoint  >> and, with proper care, even in decimal floatingpoint). >> Does computing the tangent of pi/4 radians exactly yield >> one? Maybe, for some floatingpoint KINDs. But can you >> always count on it? >> If all angles that were an integral multiple of 1/24 of >> a cycle were exactly representable in floatingpoint >> these identities could all be met. Degrees meets that >> requirement: 1/24 of a cycle is 15 degrees. And it has >> the advantage of being a measure of angles in common >> usage. I've always been surprised the engineering community >> didn't demand degree based trig functions. >There are two problems. Mixing them together does not help. >One is angle reduction which is what I mentioned and you have ignored. >The other is representation in finite precision. You have done a marvelous >job of arguing for a large composite base, like 60, and cyclic measure. Use >of degrees is a pragmatic and commonly used approximation to that. >Physics repeatedly gets into angle reduction problems. I have not seen any >of physics where degrees arise naturally instead of cycles. The places where >degrees are used do not seem to lead to much need of angle reduction which >would be quite easy if it were required.
When someone comes up with a cyclic A/D, I'll think about it Chris
On Apr 13, 12:34 pm, "Beliavsky" <beliav@aol.com> wrote: > On Apr 13, 5:57 am, "sk8terg1rl" <sk8terg1rl_2 @yahoo.co.uk> wrote: > > Hi group, > > Does anyone have recommendations for a well cited Fortran Gaussian > > pseudorandom number generator? > There are some at statistician Alan Miller's sitehttp://users.bigpond.net.au/amiller/random.html. The codes cite the > sources of the algorithms. If you want to ensure that your results are > not distorted by a poor Gaussian RNG, I suggest putting the code for > the RNG in a separate module, say random_normal_mod, so that you can > later replace the module with one using another algorithm and see if > the results are effectively the same.
Thanks everyone. I ended up going for the small (function) program rnorm.f90. Very concise coding and makes a standard normal distribution with mean 0 and stdev 1. http://users.bigpond.net.au/amiller/random/rnorm.f90 Would have been nice to be able to control a random seed for the process though, as I need to repeat my own tests at least several times with different seeds to get averages. This one seems to use a Fortran intrinsic uniform RNG with its own assumed seed. No worries though, I'll just have a dummy doloop before, that does warmup runs on "rnorm" so I get different PRN sequences based on how long the warm up has been run for. Inefficient but workable. I didn't know the references for the PRNGs were available in the code, it was handy as I've also downloaded and skimmed through the cited paper by Marsaglia & Bray [1964] from Jstor. I wasn't able to get Muller's 1959 review paper although it looks handy. skate xx P.S. There appears to be a small error in the paper. Page 260 says "Muller [7], gives a good summary of the methods available in 1958." It should probably be 1959. 6. Mervin E. Muller, An Inverse Method for the Generation of Random Normal Deviates on LargeScale Computers (1958) 7. Mervin E. Muller, A Comparison of Methods for Generating Normal Deviates on Digital Computers (1959)
On 13 Apr 2007 11:33:54 0700, sk8terg1rl <sk8terg1rl_2@yahoo.co.uk> wrote in <1176489234.070503.237@q75g2000hsh.googlegroups.com>: > http://users.bigpond.net.au/amiller/random/rnorm.f90 > Would have been nice to be able to control a random seed for the > process though, as I need to repeat my own tests at least several > times with different seeds to get averages. This one seems to use a > Fortran intrinsic uniform RNG with its own assumed seed. No worries > though, I'll just have a dummy doloop before, that does warmup runs > on "rnorm" so I get different PRN sequences based on how long the warm > up has been run for. Inefficient but workable. Well, no, it uses the standard f90 intrinsic RANDOM_NUMBER() which has an associated subroutine RANDOM_SEED() that you can use to set the seed. Check your Language Reference. You'll need to use the SIZE= optional argument to find out how large an array of integer seeds to provide, unless the default initialisation from (usually) time and date is sufficient.  Ivan Reid, School of Engineering & Design, _____________ CMS Collaboration, Brunel University. Ivan.Reid@[brunel.ac.ukcern.ch] Room 401B12, CERN KotPT  "for stupidity above and beyond the call of duty".
On Apr 13, 2:53 pm, "Dr Ivan D. Reid" <Ivan.R@brunel.ac.uk> wrote:
> On 13 Apr 2007 11:33:54 0700, sk8terg1rl <sk8terg1rl_2 @yahoo.co.uk> > wrote in <1176489234.070503.237 @q75g2000hsh.googlegroups.com>: > >http://users.bigpond.net.au/amiller/random/rnorm.f90 > > Would have been nice to be able to control a random seed for the > > process though, as I need to repeat my own tests at least several > > times with different seeds to get averages. This one seems to use a > > Fortran intrinsic uniform RNG with its own assumed seed. No worries > > though, I'll just have a dummy doloop before, that does warmup runs > > on "rnorm" so I get different PRN sequences based on how long the warm > > up has been run for. Inefficient but workable. > Well, no, it uses the standard f90 intrinsic RANDOM_NUMBER() which > has an associated subroutine RANDOM_SEED() that you can use to set the seed. > Check your Language Reference. You'll need to use the SIZE= optional > argument to find out how large an array of integer seeds to provide, > unless the default initialisation from (usually) time and date is sufficient. >  > Ivan Reid, School of Engineering & Design, _____________ CMS Collaboration, > Brunel University. Ivan.Reid@[brunel.ac.ukcern.ch] Room 401B12, CERN > KotPT  "for stupidity above and beyond the call of duty".
1. The code cited has been adapted to use F90+ instinsics  obviously not available in 1958. 2. It uses an alternate form of the BoxMuller (MarsagliaBray) transform. Instead of choosing an angle between 0 and two_pi, then using sine and cosine to resolve the components, the program generates a value inside the square which encloses the unit circle. It then rejects points outside the unit circle. [The general rejection technique goes back to vonNeumann.] 3. This technique *may* have better performance than the method which uses trig functions directly. However, the paper usually cited  from 1973  by Naylor? IIRC  tested a uniform generator with a very small multiplier. Such uniform PRNGS are known to have a problem with serial correlation. I don't know if that author's theoretical calculations, which made some assumptions and did some centrallimittheorem hand waving are still valid given a "good" uniform PRNG. So the OP is still using a form of BoxMuller (or MarsagliaBray, if you prefer).  elliot
On Apr 13, 11:23 am, Gordon Sande <g.sa@worldnet.att.net> wrote: ... > [...] The places where degrees are used do not seem to lead > to much need of angle reduction which would be quite easy if > it were required.
Argument reduction doesn't help if you then have to convert to radians or cycles in order to compute the values of the trig functions themselves. Whereas, if you had the trig functions defined on degrees, you lose no information converting your cycles to degrees and using the degree based trig functions. Winwin is better that winlose.  J. Giles "I conclude that there are two ways of constructing a software design: One way is to make it so simple that there are obviously no deficiencies and the other way is to make it so complicated that there are no obvious deficiencies."  C. A. R. Hoare
On 20070413 18:31:56 0300, "jamesgi@att.net" <jamesgi@att.net> said: > On Apr 13, 11:23 am, Gordon Sande <g.sa @worldnet.att.net> wrote: > ... >> [...] The places where degrees are used do not seem to lead >> to much need of angle reduction which would be quite easy if >> it were required. > Argument reduction doesn't help if you then have to convert > to radians or cycles in order to compute the values of the > trig functions themselves. Whereas, if you had the trig > functions defined on degrees, you lose no information > converting your cycles to degrees and using the degree > based trig functions. Winwin is better that winlose.
There seems to be three natural situations where periodic functions arise: 1. cyclic phenomena where period 1 is natural and large cycle counts commonly arise. 2. calculus/trigonometry where the "conventional" sine etc arise. The value of two pi used by the functions internally for argument reduction is often not available to users as it has a higher precision. This curious situation is a byproduct of the technical definitions used by the implementors so is not an error. 3. rational fractions of a cycle are natural and degrees are conventional. It is easy to define user functions to convert any of these to the other. It is also a nuisance. The repeated questions that are posted seem to be by those who are solving cyclic problems and dealing with the problem of the awkward periodicity of the conventional sine. Their mathematics is often equally well expressed in either cycles or radians. All too often the usual mathematical training does not make it clear that cycles are as natural radian measure. Some vendors provide degree based periodic functions as an extension. Period 1 cyclic funtions would be natural for many users. But Bessel functions are not provided either so any real improvement is not likely to happen.





