probabilities: revision, new submodules random-numbers and monte-carlo

8 files changed, 376 insertions, 8 deletions

diff --git a/build.xml b/build.xml
index ce29688b..32b05b2a 100644
--- a/build.xml
+++ b/build.xml
@@ -122,7 +122,9 @@
       <arg value="clojure.contrib.pprint.PrettyWriter"/>
       <arg value="clojure.contrib.pprint"/>
       <arg value="clojure.contrib.pprint.utilities"/>
-      <arg value="clojure.contrib.probabilities.dist"/>
+      <arg value="clojure.contrib.probabilities.finite-distributions"/>
+      <arg value="clojure.contrib.probabilities.monte-carlo"/>
+      <arg value="clojure.contrib.probabilities.random-numbers"/>
       <arg value="clojure.contrib.prxml"/>
       <arg value="clojure.contrib.repl-ln"/>
       <arg value="clojure.contrib.repl-utils"/>
diff --git a/src/clojure/contrib/gen_html_docs.clj b/src/clojure/contrib/gen_html_docs.clj
index 7a6d2729..28db5b32 100644
--- a/src/clojure/contrib/gen_html_docs.clj
+++ b/src/clojure/contrib/gen_html_docs.clj
@@ -523,8 +523,9 @@ emits the generated HTML to the path named by path."
      'clojure.contrib.json.write
      'clojure.contrib.lazy-xml.with-pull
      'clojure.contrib.miglayout.internal
-     'clojure.contrib.probabilities.dist
-     'clojure.contrib.probabilities.dist.examples
+     'clojure.contrib.probabilities.finite-distributions
+     'clojure.contrib.probabilities.monte-carlo
+     'clojure.contrib.probabilities.random-numbers
      'clojure.contrib.sql.internal
      'clojure.contrib.test-clojure.evaluation
      'clojure.contrib.test-clojure.for
diff --git a/src/clojure/contrib/load_all.clj b/src/clojure/contrib/load_all.clj
index eab7f060..fa0e9a27 100644
--- a/src/clojure/contrib/load_all.clj
+++ b/src/clojure/contrib/load_all.clj
@@ -63,7 +63,9 @@ math
 mmap
 monads
 ns-utils
-probabilities.dist
+probabilities.finite-distributions
+probabilities.monte-carlo
+probabilities.random-numbers
 prxml
 repl-ln
 repl-utils
diff --git a/src/clojure/contrib/probabilities/dist/examples.clj b/src/clojure/contrib/probabilities/examples_finite_distributions.clj
index 17bf401e..a7ae7618 100644
--- a/src/clojure/contrib/probabilities/dist/examples.clj
+++ b/src/clojure/contrib/probabilities/examples_finite_distributions.clj
@@ -6,8 +6,8 @@
 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 
-(ns clojure.contrib.probabilities.dist.examples
- (:use [clojure.contrib.probabilities.dist
+(ns clojure.contrib.probabilities.examples-finite-distributions
+ (:use [clojure.contrib.probabilities.finite-distributions
 	:only (uniform prob cond-prob join-with dist-m choose
 	      normalize certainly cond-dist-m normalize-cond)])
  (:use [clojure.contrib.monads
diff --git a/src/clojure/contrib/probabilities/examples_monte_carlo.clj b/src/clojure/contrib/probabilities/examples_monte_carlo.clj
new file mode 100644
index 00000000..13796e2f
--- /dev/null
+++ b/src/clojure/contrib/probabilities/examples_monte_carlo.clj
@@ -0,0 +1,69 @@
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+;;
+;; Monte-Carlo application examples
+;;
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+
+(ns clojure.contrib.probabilities.random.examples-monte-carlo
+  (:require [clojure.contrib.generic.collection :as gc])
+  (:use [clojure.contrib.probabilities.random-numbers
+	 :only (lcg rand-stream)])
+  (:use [clojure.contrib.probabilities.finite-distributions
+	 :only (uniform)])
+  (:use [clojure.contrib.probabilities.monte-carlo
+	 :only (random-stream discrete interval normal lognormal exponential
+		n-sphere
+		sample sample-sum sample-mean sample-mean-variance)]
+	:reload)
+  (:use [clojure.contrib.monads
+	:only (domonad state-m)]))
+
+; Create a linear congruential generator
+(def urng (lcg 259200 7141 54773 1))
+
+;; Use Clojure's built-in random number generator
+;(def urng rand-stream)
+
+; Sample transformed distributions
+(defn sample-distribution
+  [n rt]
+  (take n (gc/seq (random-stream rt urng))))
+
+; Interval [-2, 2)
+(sample-distribution 10 (interval -2 2))
+; Compare with a direct transformation
+(= (sample-distribution 10 (interval -2 2))
+   (map (fn [x] (- (* 4 x) 2)) (take 10 (gc/seq urng))))
+
+; Normal distribution
+(sample-distribution 10 (normal 0 1))
+
+; Log-Normal distribution
+(sample-distribution 10 (lognormal 0 1))
+
+; Exponential distribution
+(sample-distribution 10 (exponential 1))
+
+; n-sphere distribution
+(sample-distribution 10 (n-sphere 2 1))
+
+; Discrete distribution
+(sample-distribution 10 (discrete (uniform (range 1 7))))
+
+; Compose distributions in the state monad
+(def sum-two-dists
+  (domonad state-m
+    [r1 (interval -2 2)
+     r2 (normal 0 1)]
+    (+ r1 r2)))
+
+(sample-distribution 10 sum-two-dists)
+
+; Distribution transformations
+(sample-distribution  5 (sample 2 (interval -2 2)))
+(sample-distribution 10 (sample-sum 10 (interval -2 2)))
+(sample-distribution 10 (sample-mean 10 (interval -2 2)))
+(sample-distribution 10 (sample-mean-variance 10 (interval -2 2)))
+
diff --git a/src/clojure/contrib/probabilities/dist.clj b/src/clojure/contrib/probabilities/finite_distributions.clj
index 21008fc7..7dcc280f 100644
--- a/src/clojure/contrib/probabilities/dist.clj
+++ b/src/clojure/contrib/probabilities/finite_distributions.clj
@@ -1,7 +1,7 @@
 ;; Finite probability distributions
 
 ;; by Konrad Hinsen
-;; last updated March 2, 2009
+;; last updated April 16, 2009
 
 ;; Copyright (c) Konrad Hinsen, 2009. All rights reserved.  The use
 ;; and distribution terms for this software are covered by the Eclipse
@@ -11,7 +11,7 @@
 ;; agreeing to be bound by the terms of this license.  You must not
 ;; remove this notice, or any other, from this software.
 
-(ns clojure.contrib.probabilities.dist
+(ns clojure.contrib.probabilities.finite-distributions
   (:use [clojure.contrib.monads
 	 :only (defmonad domonad with-monad maybe-t m-lift m-chain)]
 	 [clojure.contrib.def :only (defvar)]))
diff --git a/src/clojure/contrib/probabilities/monte_carlo.clj b/src/clojure/contrib/probabilities/monte_carlo.clj
new file mode 100644
index 00000000..83bda5bf
--- /dev/null
+++ b/src/clojure/contrib/probabilities/monte_carlo.clj
@@ -0,0 +1,237 @@
+;; Monte-Carlo algorithms
+
+;; by Konrad Hinsen
+;; last updated April 21, 2009
+
+;; Copyright (c) Konrad Hinsen, 2009. All rights reserved.  The use
+;; and distribution terms for this software are covered by the Eclipse
+;; Public License 1.0 (http://opensource.org/licenses/eclipse-1.0.php)
+;; which can be found in the file epl-v10.html at the root of this
+;; distribution.  By using this software in any fashion, you are
+;; agreeing to be bound by the terms of this license.  You must not
+;; remove this notice, or any other, from this software.
+
+(ns clojure.contrib.probabilities.monte-carlo
+  "Monte-Carlo method support
+
+   Monte-Carlo methods transform an input random number stream
+   (usually having a continuous uniform distribution in the interval [0, 1))
+   into a random number stream whose distribution satisfies certain
+   conditions (usually the expectation value is equal to some desired
+   quantity). They are thus transformations from one probability distribution
+   to another one.
+
+   This library represents a Monte-Carlo method by a function that takes
+   as input the state of a random number stream with uniform distribution
+   (see clojure.contrib.probabilities.random-numbers) and returns a
+   vector containing one sample value of the desired output distribution 
+   and the final state of the input random number stream. Such functions
+   are state monad values and can be composed using operations defined
+   in clojure.contrib.monads.
+   "
+  (:use [clojure.contrib.macros :only (const)])
+  (:use [clojure.contrib.types :only (deftype)])
+  (:use [clojure.contrib.stream-utils :only (defstream stream-next)])
+  (:use [clojure.contrib.monads
+	 :only (with-monad state-m m-lift m-seq m-fmap)])
+  (:require [clojure.contrib.generic.arithmetic :as ga])
+  (:require [clojure.contrib.accumulators :as acc]))
+
+;; Random number transformers and random streams
+;;
+;; A random number transformer is a function that takes a random stream
+;; state as input and returns the next value from the transformed stream
+;; plus the new state of the input stream. Random number transformers
+;; are thus state monad values.
+;;
+;; Distributions are implemented as random number transformers that
+;; transform a uniform distribution in the interval [0, 1) to the
+;; desired distribution. Composition of such distributions allows
+;; the realization of any kind of Monte-Carlo algorithm. The result
+;; of such a composition is always again a distribution.
+;;
+;; Random streams are defined by a random number transformer and an
+;; input random number stream. If the randon number transformer represents
+;; a distribution, the input stream must have a uniform distribution
+;; in the interval [0, 1).
+
+; Random stream definition
+(deftype ::random-stream random-stream
+  "Define a random stream by a distribution and the state of a
+   random number stream with uniform distribution in [0, 1)."
+  {:arglists '([distribution random-stream-state])}
+  (fn [d rs] (list d rs)))
+
+(defstream ::random-stream
+  [[d rs]]
+  (let [[r nrs] (d rs)]
+    [r (random-stream d nrs)]))
+
+; Rejection of values is used in the construction of distributions
+(defn reject
+  "Return the distribution that results from rejecting the values from
+   dist that do not satisfy predicate p."
+  [p dist]
+  (fn [rs]
+    (let [[r nrs] (dist rs)]
+      (if (p r)
+	(recur nrs)
+	[r nrs]))))
+
+; Draw a value from a discrete distribution given as a map from
+; values to probabilities.
+; (see clojure.contrib.probabilities.finite-distributions)
+(with-monad state-m
+  (defn discrete
+    "A discrete distribution, defined by a map dist mapping values
+     to probabilities. The sum of probabilities must be one."
+    [dist]
+    (letfn [(pick-at-level [l dist-items]
+	      (let [[[x p] & rest-dist] dist-items]
+		(if (> p l)
+		  x
+		  (recur (- l p) rest-dist))))]
+      (m-fmap #(pick-at-level % (seq dist)) stream-next))))
+
+; Uniform distribution in an finite half-open interval
+(with-monad state-m
+  (defn interval
+    [a b]
+    "Transform a sequence of uniform random numbers in the interval [0, 1)
+     into a sequence of uniform random numbers in the interval [a, b)."
+    (let [d (- b a)
+	  f (if (zero? a)
+	      (if (= d 1)
+		identity
+		(fn [r] (* d r)))
+	      (if (= d 1)
+		(fn [r] (+ a r))
+		(fn [r] (+ a (* d r)))))]
+      (m-fmap f stream-next))))
+
+; Normal (Gaussian) distribution
+(defn normal
+  "Transform a sequence urs of uniform random number in the interval [0, 1)
+   into a sequence of normal random numbers with mean mu and standard
+   deviation sigma."
+  [mu sigma]
+  ; This function implements the Kinderman-Monahan ratio method:
+  ;  A.J. Kinderman & J.F. Monahan
+  ;  Computer Generation of Random Variables Using the Ratio of Uniform Deviates
+  ;  ACM Transactions on Mathematical Software 3(3) 257-260, 1977
+  (fn [rs]
+    (let [[u1  rs] (stream-next rs)
+	  [u2* rs] (stream-next rs)
+	  u2 (- 1. u2*)
+	  s (const (* 4 (/ (. Math exp (- 0.5)) (. Math sqrt 2.))))
+	  z (* s (/ (- u1 0.5) u2))
+	  zz (+ (* 0.25 z z) (. Math log u2))]
+      (if (> zz 0)
+	(recur rs)
+	[(+ mu (* sigma z)) rs]))))
+
+; Lognormal distribution
+(with-monad state-m
+  (defn lognormal
+    "Transform a sequence of uniform random numbesr in the interval [0, 1)
+     into a sequence of lognormal random numbers with mean mu and standard
+     deviation sigma."
+    [mu sigma]
+    (m-fmap #(. Math exp %) (normal mu sigma))))
+
+; Exponential distribution
+(with-monad state-m
+  (defn exponential
+    "Transform a sequence of uniform random numbers in the interval [0, 1)
+     into a sequence of exponential random numbers with parameter lambda."
+    [lambda]
+    (when (<= lambda 0)
+      (throw (IllegalArgumentException.
+  	    "exponential distribution requires a positive argument")))
+    (let [neg-inv-lambda (- (/ lambda))
+	  ; remove very small numbers to prevent log from returning -Infinity
+	  not-too-small  (reject #(< % 1e-323) stream-next)]
+      (m-fmap #(* (. Math log %) neg-inv-lambda) not-too-small))))
+
+; Another implementation of the normal distribution. It uses the
+; Box-Muller transform, but discards one of the two result values
+; at each cycle because the random number transformer interface cannot
+; handle two outputs at the same time.
+(defn normal-box-muller
+  "Transform a sequence of uniform random numbers in the interval [0, 1)
+   into a sequence of normal random numbers with mean mu and standard
+   deviation sigma."
+  [mu sigma]
+  (fn [rs]
+    (let [[u1 rs] (stream-next rs)
+	  [u2 rs] (stream-next rs)
+	   v1 (- (* 2.0 u1) 1.0)
+	   v2 (- (* 2.0 u2) 1.0)
+	   s  (+ (* v1 v1) (* v2 v2))
+	   ls (. Math sqrt (/ (* -2.0 (. Math log s)) s))
+	   x1 (* v1 ls)
+	   x2 (* v2 ls)]
+	  (if (or (>= s 1) (= s 0))
+	    (recur rs)
+	    [x1 rs]))))
+
+; Finite samples from a distribution
+(with-monad state-m
+
+    (defn sample
+      "Return the distribution of samples of length n from the
+       distribution dist"
+      [n dist]
+      (m-seq (replicate n dist)))
+
+    (defn sample-reduce
+      "Returns the distribution of the reduction of f over n samples from the
+       distribution dist."
+      ([f n dist]
+	 (if (zero? n)
+	   (m-result (f))
+	   (let [m-f    (m-lift 2 f)
+		 sample (replicate n dist)]
+	     (reduce m-f sample))))
+      ([f val n dist]
+	 (let [m-f    (m-lift 2 f)
+	       m-val  (m-result val)
+	       sample (replicate n dist)]
+	   (reduce m-f m-val sample))))
+
+    (defn sample-sum
+      "Return the distribution of the sum over n samples from the
+       distribution dist."
+      [n dist]
+      (sample-reduce ga/+ n dist))
+
+    (defn sample-mean
+      "Return the distribution of the mean over n samples from the
+       distribution dist"
+      [n dist]
+      (let [div-by-n (m-lift 1 #(ga/* % (/ n)))]
+	(div-by-n (sample-sum n dist))))
+
+    (defn sample-mean-variance
+      "Return the distribution of the mean-and-variance (a vector containing
+       the mean and the variance) over n samples from the distribution dist"
+      [n dist]
+      (let [extract (m-lift 1 (fn [mv] [(:mean mv) (:variance mv)]))]
+	(extract (sample-reduce acc/add acc/empty-mean-variance n dist))))
+
+)
+
+; Uniform distribution inside an n-sphere
+(with-monad state-m
+  (defn n-sphere
+    "Return a uniform distribution of n-dimensional vectors inside an
+     n-sphere of radius r."
+    [n r]
+    (let [box-dist    (sample n (interval (- r) r))
+	  sq          #(* % %)
+	  r-sq        (sq r)
+	  vec-sq      #(apply + (map sq %))
+	  sphere-dist (reject #(> (vec-sq %) r-sq) box-dist)
+	  as-vectors  (m-lift 1 vec)]
+      (as-vectors sphere-dist))))
+
diff --git a/src/clojure/contrib/probabilities/random_numbers.clj b/src/clojure/contrib/probabilities/random_numbers.clj
new file mode 100644
index 00000000..451a10c5
--- /dev/null
+++ b/src/clojure/contrib/probabilities/random_numbers.clj
@@ -0,0 +1,57 @@
+;; Random number generators
+
+;; by Konrad Hinsen
+;; last updated April 16, 2009
+
+;; Copyright (c) Konrad Hinsen, 2009. All rights reserved.  The use
+;; and distribution terms for this software are covered by the Eclipse
+;; Public License 1.0 (http://opensource.org/licenses/eclipse-1.0.php)
+;; which can be found in the file epl-v10.html at the root of this
+;; distribution.  By using this software in any fashion, you are
+;; agreeing to be bound by the terms of this license.  You must not
+;; remove this notice, or any other, from this software.
+
+(ns clojure.contrib.probabilities.random-numbers
+  "Random number streams
+
+   This library provides various random number generators with a common
+   stream interface. They all produce pseudo-random numbers that are uniformly
+   distributed in the interval [0, 1), i.e. 0 is a possible value but 1 isn't.
+   For transformations to other distributions, see
+   clojure.contrib.probabilities.monte-carlo."
+  (:use [clojure.contrib.types :only (deftype)])
+  (:use [clojure.contrib.stream-utils :only (defstream)])
+  (:use [clojure.contrib.def :only (defvar)]))
+
+;; Linear congruential generator
+;; http://en.wikipedia.org/wiki/Linear_congruential_generator
+
+(deftype ::lcg lcg
+  "Create a linear congruential generator"
+  {:arglists '([modulus multiplier increment seed])}
+  (fn [modulus multiplier increment seed]
+    {:m modulus :a multiplier :c increment :seed seed})
+  (fn [s] (map s (list :m :a :c :seed))))
+
+(defstream ::lcg
+  [lcg-state]
+  (let [{m :m a :a c :c seed :seed} lcg-state
+	value (/ (float seed) (float m))
+	new-seed (rem (+ c (* a seed)) m)]
+    [value (assoc lcg-state :seed new-seed)]))
+
+;; A generator based on Clojure's built-in rand function
+;; (and thus random from java.lang.Math)
+;; Note that this generator uses an internal mutable state.
+;;
+;; The state is *not* stored in the stream object and can thus
+;; *not* be restored!
+
+(defvar rand-stream (with-meta 'rand {:type ::rand-stream})
+  "A random number stream based on clojure.core/rand. Note that this
+   generator uses an internal mutable state. The state is thus not stored
+   in the stream object and cannot be restored.")
+
+(defstream ::rand-stream
+  [dummy-state]
+  [(rand) dummy-state])