>>>>> "Ray" == Raymond Toy <toy@rtp.ericsson.se> writes:
>>>>> "Tuukka" == Tuukka Toivonen <tuukkat@ees2.oulu.fi> writes:
Tuukka> On 4 Feb 2002, Raymond Toy wrote:
>>> Maxima has a Gaussian random number generator (look in bessel.lisp!),
>>> but it's rather odd. It sums up 12 numbers generated from a uniform
Tuukka> Just if someone decides to code a better generator, maybe this helps:
Tuukka> http://www.mathworks.com/company/newsletter/clevescorner/spring01_cleve.shtml
Tuukka> Somewhat interesting article that describes algorithms used in the
Tuukka> Matlab RNG.
Ray> As a matter of fact, I was going to replace it with the Ziggurat
Ray> method because I already have a Lisp implemenation and on my tests it
Ray> is far faster than any other method I've implemented.
Probably won't win on speed but I've always thought the Box-Muller
algorithm was quite cool.
(defun gaussian-random (&optional (mean 0.0) (standard-deviation 1.0))
"Generate Gaussian (normal distribution) random numbers two at a time.
This is the Box-Muller algorithm, or polar method.
See Knuth, Seminumerical Algorithms, Algorithm P, Section 3.4.1."
(multiple-value-bind
(v1 v2 sumsq)
(loop for v1 = (- (random 2.0) 1) ; number between -1 and 1 (x)
for v2 = (- (random 2.0) 1) ; number between -1 and 1 (y)
for sumsq = (+ (* v1 v1) (* v2 v2)) ; sum of the squares, r^2 on polar
until (< 0.0 sumsq 1.0) ; pick only those within unit disk
finally (return (values v1 v2 sumsq)))
(let* ((fac (sqrt (/ (* -2 (log sumsq)) sumsq))))
(values
(+ mean (* standard-deviation fac v1))
(+ mean (* standard-deviation fac v2))))))
--
Liam Healy
Liam.Healy@nrl.navy.mil