Home     |     .Net Programming    |     cSharp Home    |     Sql Server Home    |     Javascript / Client Side Development     |     Ajax Programming

Ruby on Rails Development     |     Perl Programming     |     C Programming Language     |     C++ Programming     |     IT Jobs

Python Programming Language     |     Laptop Suggestions?    |     TCL Scripting     |     Fortran Programming     |     Scheme Programming Language


 
 
Cervo Technologies
The Right Source to Outsource

MS Dynamics CRM 3.0

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.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

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:

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:

  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:

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.

Add to del.icio.us | Digg this | Stumble it | Powered by Megasolutions Inc