Mercurial > lasercutter

comparison src/clojure/contrib/probabilities/monte_carlo.clj @ 10:ef7dbbd6452c

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