Fortran Programming Language :: Well cited Fortran Gaussian PRNG?


	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.hiroshima-u.ac.jp/~m-mat/MT/VERSIONS/FORTRAN/TAKA... I know it is possible to transform a uniform distribution to a Gaussian one through the Box-Muller 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.hiroshima-u.ac.jp/~m-mat/MT/VERSIONS/FORTRAN/TAKA... > I know it is possible to transform a uniform distribution to a > Gaussian one through the Box-Muller 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 Box-Muller 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/Box-MullerTransformation.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 2007-04-13 11:09:58 -0300, glen herrmannsfeldt <g@ugcs.caltech.edu> said: > 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.hiroshima-u.ac.jp/~m-mat/MT/VERSIONS/FORTRAN/TAKA... I - Hide quoted text - >> know it is possible to transform a uniform distribution to a >> Gaussian one through the Box-Muller 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 Box-Muller 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/Box-MullerTransformation.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 Box-Muller or Ziggurat techniques... Start with the LAPACK {SDCZ}LARNV routine. Distribution option #3 uses Box-Muller. 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 one-half. In radians, 30 degrees is pi/6. It's problematical whether the calculated value of SIN(pi/6) will round to exactly one-half or not. In cycles, 30 degrees is 1/12. Again, it's problematical whether your cycle-based sine function will give exactly one-half 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 floating-point 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 floating-point - and, with proper care, even in decimal floating-point). Does computing the tangent of pi/4 radians exactly yield one? Maybe, for some floating-point 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 floating-point 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 2007-04-13 13:03:08 -0300, "jamesgi@att.net" <jamesgi@att.net> said: - Hide quoted text - > 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 one-half. In radians, 30 degrees > is pi/6. It's problematical whether the calculated > value of SIN(pi/6) will round to exactly one-half or not. > In cycles, 30 degrees is 1/12. Again, it's problematical > whether your cycle-based sine function will give exactly > one-half 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 > floating-point 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 floating-point - > and, with proper care, even in decimal floating-point). > Does computing the tangent of pi/4 radians exactly yield > one? Maybe, for some floating-point 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 floating-point > 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: - Hide quoted text - >On 2007-04-13 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 one-half. In radians, 30 degrees >> is pi/6. It's problematical whether the calculated >> value of SIN(pi/6) will round to exactly one-half or not. >> In cycles, 30 degrees is 1/12. Again, it's problematical >> whether your cycle-based sine function will give exactly >> one-half 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 >> floating-point 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 floating-point - >> and, with proper care, even in decimal floating-point). >> Does computing the tangent of pi/4 radians exactly yield >> one? Maybe, for some floating-point 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 floating-point >> 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 do-loop before, that does warm-up 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 Large-Scale 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 do-loop before, that does warm-up 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.uk\|cern.ch] Room 40-1-B12, 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: - Hide quoted text - > 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 do-loop before, that does warm-up 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.uk\|cern.ch] Room 40-1-B12, 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 Box-Muller (Marsaglia-Bray) 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 central-limit-theorem hand- waving are still valid given a "good" uniform PRNG. So the OP is still using a form of Box-Muller (or Marsaglia-Bray, 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. Win-win is better that win-lose. -- 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 2007-04-13 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. Win-win is better that win-lose. 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.