1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/src/clojure/contrib/probabilities/monte_carlo.clj	Sat Aug 21 06:25:44 2010 -0400
     1.3 @@ -0,0 +1,240 @@
     1.4 +;; Monte-Carlo algorithms
     1.5 +
     1.6 +;; by Konrad Hinsen
     1.7 +;; last updated May 3, 2009
     1.8 +
     1.9 +;; Copyright (c) Konrad Hinsen, 2009. All rights reserved.  The use
    1.10 +;; and distribution terms for this software are covered by the Eclipse
    1.11 +;; Public License 1.0 (http://opensource.org/licenses/eclipse-1.0.php)
    1.12 +;; which can be found in the file epl-v10.html at the root of this
    1.13 +;; distribution.  By using this software in any fashion, you are
    1.14 +;; agreeing to be bound by the terms of this license.  You must not
    1.15 +;; remove this notice, or any other, from this software.
    1.16 +
    1.17 +(ns
    1.18 +  ^{:author "Konrad Hinsen"
    1.19 +     :doc "Monte-Carlo method support
    1.20 +
    1.21 +           Monte-Carlo methods transform an input random number stream
    1.22 +           (usually having a continuous uniform distribution in the
    1.23 +           interval [0, 1)) into a random number stream whose distribution
    1.24 +           satisfies certain conditions (usually the expectation value
    1.25 +           is equal to some desired quantity). They are thus
    1.26 +           transformations from one probability distribution to another one.
    1.27 +
    1.28 +           This library represents a Monte-Carlo method by a function that
    1.29 +           takes as input the state of a random number stream with
    1.30 +           uniform distribution (see
    1.31 +           clojure.contrib.probabilities.random-numbers) and returns a
    1.32 +           vector containing one sample value of the desired output
    1.33 +           distribution and the final state of the input random number
    1.34 +           stream. Such functions are state monad values and can be
    1.35 +           composed using operations defined in clojure.contrib.monads."}
    1.36 +  clojure.contrib.probabilities.monte-carlo
    1.37 +  (:refer-clojure :exclude (deftype))
    1.38 +  (:use [clojure.contrib.macros :only (const)])
    1.39 +  (:use [clojure.contrib.types :only (deftype)])
    1.40 +  (:use [clojure.contrib.stream-utils :only (defstream stream-next)])
    1.41 +  (:use [clojure.contrib.monads
    1.42 +	 :only (with-monad state-m m-lift m-seq m-fmap)])
    1.43 +  (:require [clojure.contrib.generic.arithmetic :as ga])
    1.44 +  (:require [clojure.contrib.accumulators :as acc]))
    1.45 +
    1.46 +;; Random number transformers and random streams
    1.47 +;;
    1.48 +;; A random number transformer is a function that takes a random stream
    1.49 +;; state as input and returns the next value from the transformed stream
    1.50 +;; plus the new state of the input stream. Random number transformers
    1.51 +;; are thus state monad values.
    1.52 +;;
    1.53 +;; Distributions are implemented as random number transformers that
    1.54 +;; transform a uniform distribution in the interval [0, 1) to the
    1.55 +;; desired distribution. Composition of such distributions allows
    1.56 +;; the realization of any kind of Monte-Carlo algorithm. The result
    1.57 +;; of such a composition is always again a distribution.
    1.58 +;;
    1.59 +;; Random streams are defined by a random number transformer and an
    1.60 +;; input random number stream. If the randon number transformer represents
    1.61 +;; a distribution, the input stream must have a uniform distribution
    1.62 +;; in the interval [0, 1).
    1.63 +
    1.64 +; Random stream definition
    1.65 +(deftype ::random-stream random-stream
    1.66 +  "Define a random stream by a distribution and the state of a
    1.67 +   random number stream with uniform distribution in [0, 1)."
    1.68 +  {:arglists '([distribution random-stream-state])}
    1.69 +  (fn [d rs] (list d rs)))
    1.70 +
    1.71 +(defstream ::random-stream
    1.72 +  [[d rs]]
    1.73 +  (let [[r nrs] (d rs)]
    1.74 +    [r (random-stream d nrs)]))
    1.75 +
    1.76 +; Rejection of values is used in the construction of distributions
    1.77 +(defn reject
    1.78 +  "Return the distribution that results from rejecting the values from
    1.79 +   dist that do not satisfy predicate p."
    1.80 +  [p dist]
    1.81 +  (fn [rs]
    1.82 +    (let [[r nrs] (dist rs)]
    1.83 +      (if (p r)
    1.84 +	(recur nrs)
    1.85 +	[r nrs]))))
    1.86 +
    1.87 +; Draw a value from a discrete distribution given as a map from
    1.88 +; values to probabilities.
    1.89 +; (see clojure.contrib.probabilities.finite-distributions)
    1.90 +(with-monad state-m
    1.91 +  (defn discrete
    1.92 +    "A discrete distribution, defined by a map dist mapping values
    1.93 +     to probabilities. The sum of probabilities must be one."
    1.94 +    [dist]
    1.95 +    (letfn [(pick-at-level [l dist-items]
    1.96 +	      (let [[[x p] & rest-dist] dist-items]
    1.97 +		(if (> p l)
    1.98 +		  x
    1.99 +		  (recur (- l p) rest-dist))))]
   1.100 +      (m-fmap #(pick-at-level % (seq dist)) stream-next))))
   1.101 +
   1.102 +; Uniform distribution in an finite half-open interval
   1.103 +(with-monad state-m
   1.104 +  (defn interval
   1.105 +    [a b]
   1.106 +    "Transform a sequence of uniform random numbers in the interval [0, 1)
   1.107 +     into a sequence of uniform random numbers in the interval [a, b)."
   1.108 +    (let [d (- b a)
   1.109 +	  f (if (zero? a)
   1.110 +	      (if (= d 1)
   1.111 +		identity
   1.112 +		(fn [r] (* d r)))
   1.113 +	      (if (= d 1)
   1.114 +		(fn [r] (+ a r))
   1.115 +		(fn [r] (+ a (* d r)))))]
   1.116 +      (m-fmap f stream-next))))
   1.117 +
   1.118 +; Normal (Gaussian) distribution
   1.119 +(defn normal
   1.120 +  "Transform a sequence urs of uniform random number in the interval [0, 1)
   1.121 +   into a sequence of normal random numbers with mean mu and standard
   1.122 +   deviation sigma."
   1.123 +  [mu sigma]
   1.124 +  ; This function implements the Kinderman-Monahan ratio method:
   1.125 +  ;  A.J. Kinderman & J.F. Monahan
   1.126 +  ;  Computer Generation of Random Variables Using the Ratio of Uniform Deviates
   1.127 +  ;  ACM Transactions on Mathematical Software 3(3) 257-260, 1977
   1.128 +  (fn [rs]
   1.129 +    (let [[u1  rs] (stream-next rs)
   1.130 +	  [u2* rs] (stream-next rs)
   1.131 +	  u2 (- 1. u2*)
   1.132 +	  s (const (* 4 (/ (. Math exp (- 0.5)) (. Math sqrt 2.))))
   1.133 +	  z (* s (/ (- u1 0.5) u2))
   1.134 +	  zz (+ (* 0.25 z z) (. Math log u2))]
   1.135 +      (if (> zz 0)
   1.136 +	(recur rs)
   1.137 +	[(+ mu (* sigma z)) rs]))))
   1.138 +
   1.139 +; Lognormal distribution
   1.140 +(with-monad state-m
   1.141 +  (defn lognormal
   1.142 +    "Transform a sequence of uniform random numbesr in the interval [0, 1)
   1.143 +     into a sequence of lognormal random numbers with mean mu and standard
   1.144 +     deviation sigma."
   1.145 +    [mu sigma]
   1.146 +    (m-fmap #(. Math exp %) (normal mu sigma))))
   1.147 +
   1.148 +; Exponential distribution
   1.149 +(with-monad state-m
   1.150 +  (defn exponential
   1.151 +    "Transform a sequence of uniform random numbers in the interval [0, 1)
   1.152 +     into a sequence of exponential random numbers with parameter lambda."
   1.153 +    [lambda]
   1.154 +    (when (<= lambda 0)
   1.155 +      (throw (IllegalArgumentException.
   1.156 +  	    "exponential distribution requires a positive argument")))
   1.157 +    (let [neg-inv-lambda (- (/ lambda))
   1.158 +	  ; remove very small numbers to prevent log from returning -Infinity
   1.159 +	  not-too-small  (reject #(< % 1e-323) stream-next)]
   1.160 +      (m-fmap #(* (. Math log %) neg-inv-lambda) not-too-small))))
   1.161 +
   1.162 +; Another implementation of the normal distribution. It uses the
   1.163 +; Box-Muller transform, but discards one of the two result values
   1.164 +; at each cycle because the random number transformer interface cannot
   1.165 +; handle two outputs at the same time.
   1.166 +(defn normal-box-muller
   1.167 +  "Transform a sequence of uniform random numbers in the interval [0, 1)
   1.168 +   into a sequence of normal random numbers with mean mu and standard
   1.169 +   deviation sigma."
   1.170 +  [mu sigma]
   1.171 +  (fn [rs]
   1.172 +    (let [[u1 rs] (stream-next rs)
   1.173 +	  [u2 rs] (stream-next rs)
   1.174 +	   v1 (- (* 2.0 u1) 1.0)
   1.175 +	   v2 (- (* 2.0 u2) 1.0)
   1.176 +	   s  (+ (* v1 v1) (* v2 v2))
   1.177 +	   ls (. Math sqrt (/ (* -2.0 (. Math log s)) s))
   1.178 +	   x1 (* v1 ls)
   1.179 +	   x2 (* v2 ls)]
   1.180 +	  (if (or (>= s 1) (= s 0))
   1.181 +	    (recur rs)
   1.182 +	    [x1 rs]))))
   1.183 +
   1.184 +; Finite samples from a distribution
   1.185 +(with-monad state-m
   1.186 +
   1.187 +    (defn sample
   1.188 +      "Return the distribution of samples of length n from the
   1.189 +       distribution dist"
   1.190 +      [n dist]
   1.191 +      (m-seq (replicate n dist)))
   1.192 +
   1.193 +    (defn sample-reduce
   1.194 +      "Returns the distribution of the reduction of f over n samples from the
   1.195 +       distribution dist."
   1.196 +      ([f n dist]
   1.197 +	 (if (zero? n)
   1.198 +	   (m-result (f))
   1.199 +	   (let [m-f    (m-lift 2 f)
   1.200 +		 sample (replicate n dist)]
   1.201 +	     (reduce m-f sample))))
   1.202 +      ([f val n dist]
   1.203 +	 (let [m-f    (m-lift 2 f)
   1.204 +	       m-val  (m-result val)
   1.205 +	       sample (replicate n dist)]
   1.206 +	   (reduce m-f m-val sample))))
   1.207 +
   1.208 +    (defn sample-sum
   1.209 +      "Return the distribution of the sum over n samples from the
   1.210 +       distribution dist."
   1.211 +      [n dist]
   1.212 +      (sample-reduce ga/+ n dist))
   1.213 +
   1.214 +    (defn sample-mean
   1.215 +      "Return the distribution of the mean over n samples from the
   1.216 +       distribution dist"
   1.217 +      [n dist]
   1.218 +      (let [div-by-n (m-lift 1 #(ga/* % (/ n)))]
   1.219 +	(div-by-n (sample-sum n dist))))
   1.220 +
   1.221 +    (defn sample-mean-variance
   1.222 +      "Return the distribution of the mean-and-variance (a vector containing
   1.223 +       the mean and the variance) over n samples from the distribution dist"
   1.224 +      [n dist]
   1.225 +      (let [extract (m-lift 1 (fn [mv] [(:mean mv) (:variance mv)]))]
   1.226 +	(extract (sample-reduce acc/add acc/empty-mean-variance n dist))))
   1.227 +
   1.228 +)
   1.229 +
   1.230 +; Uniform distribution inside an n-sphere
   1.231 +(with-monad state-m
   1.232 +  (defn n-sphere
   1.233 +    "Return a uniform distribution of n-dimensional vectors inside an
   1.234 +     n-sphere of radius r."
   1.235 +    [n r]
   1.236 +    (let [box-dist    (sample n (interval (- r) r))
   1.237 +	  sq          #(* % %)
   1.238 +	  r-sq        (sq r)
   1.239 +	  vec-sq      #(apply + (map sq %))
   1.240 +	  sphere-dist (reject #(> (vec-sq %) r-sq) box-dist)
   1.241 +	  as-vectors  (m-lift 1 vec)]
   1.242 +      (as-vectors sphere-dist))))
   1.243 +