So I'm going to stop pretending like I'm an expert and actually post
some Clojure code. Be constructively critical 'cause I'm a n00b in
that regard ;-) This is a pseudorandom number generator for the
Gaussian (0,1) distribution.
(defn next-gaussrand-state [current-
state]
^{:doc "Given the current state (current-state) of the Gaussian-
(0,1)
distribution pseudorandom number generator, return the next state.
This
operation is nondestructive as long as current-state's underlying
uniform
pseudorandom number generator is nondestructive. next-gaussrand-state
is
a low-level routine out of which one can build a seq of
Gaussian-
distribution pseudorandom numbers. See notes in the text
below."}
(if (= (:phase current-state)
0)
(loop [state current-
state]
(let [prng (:uniform-prng-seq
state)
U1 (first
prng)
U2 (frest
prng)
prng-rrest (rrest
prng)
V1 (- (* 2.0 U1)
1.0)
V2 (- (* 2.0 U2)
1.0)
S (+ (* V1 V1) (* V2
V2))
X (* V1 (. Math sqrt (/ (* -2.0 (. Math log S))
S)))]
(if (or (>= S 1) (= S
0))
;; Reject the out-of-range sample. Recall that U1 and U2
have
;; been shifted from (0,1] to (-1, 1]. It's entirely
possible
;; for S to be > 1 -- no more than 2 in exact arithmetic,
and
I
;; suspect V1 and V2 are computed exactly for the extreme
case
;; of U1=1 and
U2=1.
(recur {:uniform-prng-seq prng-rrest :X X :V1 V1 :V2 V2 :S
S :phase (:phase
state)})
{:uniform-prng-seq prng-rrest :X X :V1 V1 :V2 V2 :S S :phase
(- 1 (:phase state))})))
(let [V1 (:V1 current-
state)
V2 (:V2 current-
state)
S (:S current-
state)]
{:uniform-prng-seq (:uniform-prng-seq current-
state)
:X (* V2 (. Math sqrt (/ (* -2.0 (. Math log S))
S)))
:V1 V1 :V2 V2 :S S :phase (- 1.0 (:phase current-
state))})))
defn first-gaussrand-state [uniform-prng-
seq]
^{:doc "Given a seq of uniformly distributed pseudorandom numbers on
(0,1),
return the initial state of a seq of Gaussian-(0,1) pseudorandom
numbers.
"}
(next-gaussrand-state {:X 0.0 :V1 0.0 :V2 0.0 :S 0.0 :phase
0 :uniform-prng-seq uniform-prng-
seq}))
(defn gaussrand-value [gaussrand-
state]
(:X gaussrand-
state))
(defn new-gaussrand-seq [uniform-prng-
seq]
^{:doc "Return a seq of Gaussian-(0,1) (mean 1 and standard
deviation
0)
pseudorandom numbers, given a seq of uniformly distributed (on
(0,1))
pseudorandom floating-point
numbers."}
(map gaussrand-value (iterate next-gaussrand-state (first-gaussrand-
state uniform-prng-
seq))))
(defn test-gaussrand
[n]
^{:doc "Simple test for new-gaussrand-seq. Should return a small
floating-
point
number that gets smaller as n (a positive integer) gets
bigger."}
(let [PRNG (new
java.util.Random)
;; Somewhat broken because java.util.Random is not stateless,
and
;; the nextDouble() method changes its state. However,
this
;; particular java.util.Random object is not shared outside of
this
;; routine, so we can guarantee that uniform-prng-seq is a
true
seq
;; in the context of this
method.
uniform-prng-seq (iterate (fn [_] (. PRNG nextDouble)) (. PRNG
nextDouble))
gaussrand-seq (new-gaussrand-seq uniform-prng-
seq)]
(/ (reduce + 0 (take n gaussrand-seq)) n)))
Notes on next-gaussrand-state:
This implementation was ported almost naively from the C FAQ
online:
http://c-faq.com/lib/gaussian.html
Knuth discusses this method in his classic tome ("The Art of Computer
Programming, Volume 2: Seminumerical Algorithms", section 3.4.1,
subsection C, algorithm P), in which he cites the following:
G. Marsaglia and T. A. Bray, "A convenient method for generating
normal variables," SIAM Review, volume 6, issue 3, pp. 260-4, July
1964.
After porting the C FAQ code, I discovered that the same method is
also used by java.util.Random.nextDouble(), at least as of Java
version 1.4.2: see
http://java.sun.com/j2se/1.4.2/docs/api/java/util/Random.html#nextGaussian()
There it is called "the polar method of G. E. P. Box, M. E. Muller,
and G. Marsaglia." (Where is this G. E. P. Box in the list of
authors in the SIAM Review article? I don't own a copy of Knuth so
I'll have to dig a little to find out.) Note also that Marsaglia and
others have written further works on the subject, even in the last few
years, so the problem is by no means dead. I haven't looked more
closely at this algorithm so I include it here merely as a useful
exercise for myself and not because I intend others to go and copy it
without questioning it. Anyway, the C FAQ makes a mistake multiple
times in the aforementioned article: it tries to get uniformly
distributed pseudorandom doubles on (0,1) by scaling the output of rand
(). rand() only gives you 32 (or 31?) bits of entropy but you need 53
bits for a double-precision floating-point number. That's why Java
calls nextDouble() (which according to Sun's Java documentation does
the Right Thing). So take note if you roll your own PRNG not to do
that ;-)
mfh
--~--~---------~--~----~------------~-------~--~----~
You received this message because you are subscribed to the Google Groups
"Clojure" group.
To post to this group, send email to clojure@googlegroups.com
To unsubscribe from this group, send email to
clojure+unsubscr...@googlegroups.com
For more options, visit this group at
http://groups.google.com/group/clojure?hl=en
-~----------~----~----~----~------~----~------~--~---
