|
rlm@10
|
1 ;; Monte-Carlo algorithms
|
|
rlm@10
|
2
|
|
rlm@10
|
3 ;; by Konrad Hinsen
|
|
rlm@10
|
4 ;; last updated May 3, 2009
|
|
rlm@10
|
5
|
|
rlm@10
|
6 ;; Copyright (c) Konrad Hinsen, 2009. All rights reserved. The use
|
|
rlm@10
|
7 ;; and distribution terms for this software are covered by the Eclipse
|
|
rlm@10
|
8 ;; Public License 1.0 (http://opensource.org/licenses/eclipse-1.0.php)
|
|
rlm@10
|
9 ;; which can be found in the file epl-v10.html at the root of this
|
|
rlm@10
|
10 ;; distribution. By using this software in any fashion, you are
|
|
rlm@10
|
11 ;; agreeing to be bound by the terms of this license. You must not
|
|
rlm@10
|
12 ;; remove this notice, or any other, from this software.
|
|
rlm@10
|
13
|
|
rlm@10
|
14 (ns
|
|
rlm@10
|
15 ^{:author "Konrad Hinsen"
|
|
rlm@10
|
16 :doc "Monte-Carlo method support
|
|
rlm@10
|
17
|
|
rlm@10
|
18 Monte-Carlo methods transform an input random number stream
|
|
rlm@10
|
19 (usually having a continuous uniform distribution in the
|
|
rlm@10
|
20 interval [0, 1)) into a random number stream whose distribution
|
|
rlm@10
|
21 satisfies certain conditions (usually the expectation value
|
|
rlm@10
|
22 is equal to some desired quantity). They are thus
|
|
rlm@10
|
23 transformations from one probability distribution to another one.
|
|
rlm@10
|
24
|
|
rlm@10
|
25 This library represents a Monte-Carlo method by a function that
|
|
rlm@10
|
26 takes as input the state of a random number stream with
|
|
rlm@10
|
27 uniform distribution (see
|
|
rlm@10
|
28 clojure.contrib.probabilities.random-numbers) and returns a
|
|
rlm@10
|
29 vector containing one sample value of the desired output
|
|
rlm@10
|
30 distribution and the final state of the input random number
|
|
rlm@10
|
31 stream. Such functions are state monad values and can be
|
|
rlm@10
|
32 composed using operations defined in clojure.contrib.monads."}
|
|
rlm@10
|
33 clojure.contrib.probabilities.monte-carlo
|
|
rlm@10
|
34 (:refer-clojure :exclude (deftype))
|
|
rlm@10
|
35 (:use [clojure.contrib.macros :only (const)])
|
|
rlm@10
|
36 (:use [clojure.contrib.types :only (deftype)])
|
|
rlm@10
|
37 (:use [clojure.contrib.stream-utils :only (defstream stream-next)])
|
|
rlm@10
|
38 (:use [clojure.contrib.monads
|
|
rlm@10
|
39 :only (with-monad state-m m-lift m-seq m-fmap)])
|
|
rlm@10
|
40 (:require [clojure.contrib.generic.arithmetic :as ga])
|
|
rlm@10
|
41 (:require [clojure.contrib.accumulators :as acc]))
|
|
rlm@10
|
42
|
|
rlm@10
|
43 ;; Random number transformers and random streams
|
|
rlm@10
|
44 ;;
|
|
rlm@10
|
45 ;; A random number transformer is a function that takes a random stream
|
|
rlm@10
|
46 ;; state as input and returns the next value from the transformed stream
|
|
rlm@10
|
47 ;; plus the new state of the input stream. Random number transformers
|
|
rlm@10
|
48 ;; are thus state monad values.
|
|
rlm@10
|
49 ;;
|
|
rlm@10
|
50 ;; Distributions are implemented as random number transformers that
|
|
rlm@10
|
51 ;; transform a uniform distribution in the interval [0, 1) to the
|
|
rlm@10
|
52 ;; desired distribution. Composition of such distributions allows
|
|
rlm@10
|
53 ;; the realization of any kind of Monte-Carlo algorithm. The result
|
|
rlm@10
|
54 ;; of such a composition is always again a distribution.
|
|
rlm@10
|
55 ;;
|
|
rlm@10
|
56 ;; Random streams are defined by a random number transformer and an
|
|
rlm@10
|
57 ;; input random number stream. If the randon number transformer represents
|
|
rlm@10
|
58 ;; a distribution, the input stream must have a uniform distribution
|
|
rlm@10
|
59 ;; in the interval [0, 1).
|
|
rlm@10
|
60
|
|
rlm@10
|
61 ; Random stream definition
|
|
rlm@10
|
62 (deftype ::random-stream random-stream
|
|
rlm@10
|
63 "Define a random stream by a distribution and the state of a
|
|
rlm@10
|
64 random number stream with uniform distribution in [0, 1)."
|
|
rlm@10
|
65 {:arglists '([distribution random-stream-state])}
|
|
rlm@10
|
66 (fn [d rs] (list d rs)))
|
|
rlm@10
|
67
|
|
rlm@10
|
68 (defstream ::random-stream
|
|
rlm@10
|
69 [[d rs]]
|
|
rlm@10
|
70 (let [[r nrs] (d rs)]
|
|
rlm@10
|
71 [r (random-stream d nrs)]))
|
|
rlm@10
|
72
|
|
rlm@10
|
73 ; Rejection of values is used in the construction of distributions
|
|
rlm@10
|
74 (defn reject
|
|
rlm@10
|
75 "Return the distribution that results from rejecting the values from
|
|
rlm@10
|
76 dist that do not satisfy predicate p."
|
|
rlm@10
|
77 [p dist]
|
|
rlm@10
|
78 (fn [rs]
|
|
rlm@10
|
79 (let [[r nrs] (dist rs)]
|
|
rlm@10
|
80 (if (p r)
|
|
rlm@10
|
81 (recur nrs)
|
|
rlm@10
|
82 [r nrs]))))
|
|
rlm@10
|
83
|
|
rlm@10
|
84 ; Draw a value from a discrete distribution given as a map from
|
|
rlm@10
|
85 ; values to probabilities.
|
|
rlm@10
|
86 ; (see clojure.contrib.probabilities.finite-distributions)
|
|
rlm@10
|
87 (with-monad state-m
|
|
rlm@10
|
88 (defn discrete
|
|
rlm@10
|
89 "A discrete distribution, defined by a map dist mapping values
|
|
rlm@10
|
90 to probabilities. The sum of probabilities must be one."
|
|
rlm@10
|
91 [dist]
|
|
rlm@10
|
92 (letfn [(pick-at-level [l dist-items]
|
|
rlm@10
|
93 (let [[[x p] & rest-dist] dist-items]
|
|
rlm@10
|
94 (if (> p l)
|
|
rlm@10
|
95 x
|
|
rlm@10
|
96 (recur (- l p) rest-dist))))]
|
|
rlm@10
|
97 (m-fmap #(pick-at-level % (seq dist)) stream-next))))
|
|
rlm@10
|
98
|
|
rlm@10
|
99 ; Uniform distribution in an finite half-open interval
|
|
rlm@10
|
100 (with-monad state-m
|
|
rlm@10
|
101 (defn interval
|
|
rlm@10
|
102 [a b]
|
|
rlm@10
|
103 "Transform a sequence of uniform random numbers in the interval [0, 1)
|
|
rlm@10
|
104 into a sequence of uniform random numbers in the interval [a, b)."
|
|
rlm@10
|
105 (let [d (- b a)
|
|
rlm@10
|
106 f (if (zero? a)
|
|
rlm@10
|
107 (if (= d 1)
|
|
rlm@10
|
108 identity
|
|
rlm@10
|
109 (fn [r] (* d r)))
|
|
rlm@10
|
110 (if (= d 1)
|
|
rlm@10
|
111 (fn [r] (+ a r))
|
|
rlm@10
|
112 (fn [r] (+ a (* d r)))))]
|
|
rlm@10
|
113 (m-fmap f stream-next))))
|
|
rlm@10
|
114
|
|
rlm@10
|
115 ; Normal (Gaussian) distribution
|
|
rlm@10
|
116 (defn normal
|
|
rlm@10
|
117 "Transform a sequence urs of uniform random number in the interval [0, 1)
|
|
rlm@10
|
118 into a sequence of normal random numbers with mean mu and standard
|
|
rlm@10
|
119 deviation sigma."
|
|
rlm@10
|
120 [mu sigma]
|
|
rlm@10
|
121 ; This function implements the Kinderman-Monahan ratio method:
|
|
rlm@10
|
122 ; A.J. Kinderman & J.F. Monahan
|
|
rlm@10
|
123 ; Computer Generation of Random Variables Using the Ratio of Uniform Deviates
|
|
rlm@10
|
124 ; ACM Transactions on Mathematical Software 3(3) 257-260, 1977
|
|
rlm@10
|
125 (fn [rs]
|
|
rlm@10
|
126 (let [[u1 rs] (stream-next rs)
|
|
rlm@10
|
127 [u2* rs] (stream-next rs)
|
|
rlm@10
|
128 u2 (- 1. u2*)
|
|
rlm@10
|
129 s (const (* 4 (/ (. Math exp (- 0.5)) (. Math sqrt 2.))))
|
|
rlm@10
|
130 z (* s (/ (- u1 0.5) u2))
|
|
rlm@10
|
131 zz (+ (* 0.25 z z) (. Math log u2))]
|
|
rlm@10
|
132 (if (> zz 0)
|
|
rlm@10
|
133 (recur rs)
|
|
rlm@10
|
134 [(+ mu (* sigma z)) rs]))))
|
|
rlm@10
|
135
|
|
rlm@10
|
136 ; Lognormal distribution
|
|
rlm@10
|
137 (with-monad state-m
|
|
rlm@10
|
138 (defn lognormal
|
|
rlm@10
|
139 "Transform a sequence of uniform random numbesr in the interval [0, 1)
|
|
rlm@10
|
140 into a sequence of lognormal random numbers with mean mu and standard
|
|
rlm@10
|
141 deviation sigma."
|
|
rlm@10
|
142 [mu sigma]
|
|
rlm@10
|
143 (m-fmap #(. Math exp %) (normal mu sigma))))
|
|
rlm@10
|
144
|
|
rlm@10
|
145 ; Exponential distribution
|
|
rlm@10
|
146 (with-monad state-m
|
|
rlm@10
|
147 (defn exponential
|
|
rlm@10
|
148 "Transform a sequence of uniform random numbers in the interval [0, 1)
|
|
rlm@10
|
149 into a sequence of exponential random numbers with parameter lambda."
|
|
rlm@10
|
150 [lambda]
|
|
rlm@10
|
151 (when (<= lambda 0)
|
|
rlm@10
|
152 (throw (IllegalArgumentException.
|
|
rlm@10
|
153 "exponential distribution requires a positive argument")))
|
|
rlm@10
|
154 (let [neg-inv-lambda (- (/ lambda))
|
|
rlm@10
|
155 ; remove very small numbers to prevent log from returning -Infinity
|
|
rlm@10
|
156 not-too-small (reject #(< % 1e-323) stream-next)]
|
|
rlm@10
|
157 (m-fmap #(* (. Math log %) neg-inv-lambda) not-too-small))))
|
|
rlm@10
|
158
|
|
rlm@10
|
159 ; Another implementation of the normal distribution. It uses the
|
|
rlm@10
|
160 ; Box-Muller transform, but discards one of the two result values
|
|
rlm@10
|
161 ; at each cycle because the random number transformer interface cannot
|
|
rlm@10
|
162 ; handle two outputs at the same time.
|
|
rlm@10
|
163 (defn normal-box-muller
|
|
rlm@10
|
164 "Transform a sequence of uniform random numbers in the interval [0, 1)
|
|
rlm@10
|
165 into a sequence of normal random numbers with mean mu and standard
|
|
rlm@10
|
166 deviation sigma."
|
|
rlm@10
|
167 [mu sigma]
|
|
rlm@10
|
168 (fn [rs]
|
|
rlm@10
|
169 (let [[u1 rs] (stream-next rs)
|
|
rlm@10
|
170 [u2 rs] (stream-next rs)
|
|
rlm@10
|
171 v1 (- (* 2.0 u1) 1.0)
|
|
rlm@10
|
172 v2 (- (* 2.0 u2) 1.0)
|
|
rlm@10
|
173 s (+ (* v1 v1) (* v2 v2))
|
|
rlm@10
|
174 ls (. Math sqrt (/ (* -2.0 (. Math log s)) s))
|
|
rlm@10
|
175 x1 (* v1 ls)
|
|
rlm@10
|
176 x2 (* v2 ls)]
|
|
rlm@10
|
177 (if (or (>= s 1) (= s 0))
|
|
rlm@10
|
178 (recur rs)
|
|
rlm@10
|
179 [x1 rs]))))
|
|
rlm@10
|
180
|
|
rlm@10
|
181 ; Finite samples from a distribution
|
|
rlm@10
|
182 (with-monad state-m
|
|
rlm@10
|
183
|
|
rlm@10
|
184 (defn sample
|
|
rlm@10
|
185 "Return the distribution of samples of length n from the
|
|
rlm@10
|
186 distribution dist"
|
|
rlm@10
|
187 [n dist]
|
|
rlm@10
|
188 (m-seq (replicate n dist)))
|
|
rlm@10
|
189
|
|
rlm@10
|
190 (defn sample-reduce
|
|
rlm@10
|
191 "Returns the distribution of the reduction of f over n samples from the
|
|
rlm@10
|
192 distribution dist."
|
|
rlm@10
|
193 ([f n dist]
|
|
rlm@10
|
194 (if (zero? n)
|
|
rlm@10
|
195 (m-result (f))
|
|
rlm@10
|
196 (let [m-f (m-lift 2 f)
|
|
rlm@10
|
197 sample (replicate n dist)]
|
|
rlm@10
|
198 (reduce m-f sample))))
|
|
rlm@10
|
199 ([f val n dist]
|
|
rlm@10
|
200 (let [m-f (m-lift 2 f)
|
|
rlm@10
|
201 m-val (m-result val)
|
|
rlm@10
|
202 sample (replicate n dist)]
|
|
rlm@10
|
203 (reduce m-f m-val sample))))
|
|
rlm@10
|
204
|
|
rlm@10
|
205 (defn sample-sum
|
|
rlm@10
|
206 "Return the distribution of the sum over n samples from the
|
|
rlm@10
|
207 distribution dist."
|
|
rlm@10
|
208 [n dist]
|
|
rlm@10
|
209 (sample-reduce ga/+ n dist))
|
|
rlm@10
|
210
|
|
rlm@10
|
211 (defn sample-mean
|
|
rlm@10
|
212 "Return the distribution of the mean over n samples from the
|
|
rlm@10
|
213 distribution dist"
|
|
rlm@10
|
214 [n dist]
|
|
rlm@10
|
215 (let [div-by-n (m-lift 1 #(ga/* % (/ n)))]
|
|
rlm@10
|
216 (div-by-n (sample-sum n dist))))
|
|
rlm@10
|
217
|
|
rlm@10
|
218 (defn sample-mean-variance
|
|
rlm@10
|
219 "Return the distribution of the mean-and-variance (a vector containing
|
|
rlm@10
|
220 the mean and the variance) over n samples from the distribution dist"
|
|
rlm@10
|
221 [n dist]
|
|
rlm@10
|
222 (let [extract (m-lift 1 (fn [mv] [(:mean mv) (:variance mv)]))]
|
|
rlm@10
|
223 (extract (sample-reduce acc/add acc/empty-mean-variance n dist))))
|
|
rlm@10
|
224
|
|
rlm@10
|
225 )
|
|
rlm@10
|
226
|
|
rlm@10
|
227 ; Uniform distribution inside an n-sphere
|
|
rlm@10
|
228 (with-monad state-m
|
|
rlm@10
|
229 (defn n-sphere
|
|
rlm@10
|
230 "Return a uniform distribution of n-dimensional vectors inside an
|
|
rlm@10
|
231 n-sphere of radius r."
|
|
rlm@10
|
232 [n r]
|
|
rlm@10
|
233 (let [box-dist (sample n (interval (- r) r))
|
|
rlm@10
|
234 sq #(* % %)
|
|
rlm@10
|
235 r-sq (sq r)
|
|
rlm@10
|
236 vec-sq #(apply + (map sq %))
|
|
rlm@10
|
237 sphere-dist (reject #(> (vec-sq %) r-sq) box-dist)
|
|
rlm@10
|
238 as-vectors (m-lift 1 vec)]
|
|
rlm@10
|
239 (as-vectors sphere-dist))))
|
|
rlm@10
|
240
|
