Mercurial > lasercutter
comparison src/clojure/contrib/test_contrib/probabilities/examples_finite_distributions.clj @ 10:ef7dbbd6452c
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 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; | |
2 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; | |
3 ;; | |
4 ;; Probability distribution application examples | |
5 ;; | |
6 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; | |
7 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; | |
8 | |
9 (ns | |
10 #^{:author "Konrad Hinsen" | |
11 :skip-wiki true | |
12 :doc "Examples for finite probability distribution"} | |
13 clojure.contrib.probabilities.examples-finite-distributions | |
14 (:use [clojure.contrib.probabilities.finite-distributions | |
15 :only (uniform prob cond-prob join-with dist-m choose | |
16 normalize certainly cond-dist-m normalize-cond)]) | |
17 (:use [clojure.contrib.monads | |
18 :only (domonad with-monad m-seq m-chain m-lift)]) | |
19 (:require clojure.contrib.accumulators)) | |
20 | |
21 ;; Simple examples using dice | |
22 | |
23 ; A single die is represented by a uniform distribution over the | |
24 ; six possible outcomes. | |
25 (def die (uniform #{1 2 3 4 5 6})) | |
26 | |
27 ; The probability that the result is odd... | |
28 (prob odd? die) | |
29 ; ... or greater than four. | |
30 (prob #(> % 4) die) | |
31 | |
32 ; The sum of two dice | |
33 (def two-dice (join-with + die die)) | |
34 (prob #(> % 6) two-dice) | |
35 | |
36 ; The sum of two dice using a monad comprehension | |
37 (assert (= two-dice | |
38 (domonad dist-m | |
39 [d1 die | |
40 d2 die] | |
41 (+ d1 d2)))) | |
42 | |
43 ; The two values separately, but as an ordered pair | |
44 (domonad dist-m | |
45 [d1 die | |
46 d2 die] | |
47 (if (< d1 d2) (list d1 d2) (list d2 d1))) | |
48 | |
49 ; The conditional probability for two dice yielding X if X is odd: | |
50 (cond-prob odd? two-dice) | |
51 | |
52 ; A two-step experiment: throw a die, and then add 1 with probability 1/2 | |
53 (domonad dist-m | |
54 [d die | |
55 x (choose (/ 1 2) d | |
56 :else (inc d))] | |
57 x) | |
58 | |
59 ; The sum of n dice | |
60 (defn dice [n] | |
61 (domonad dist-m | |
62 [ds (m-seq (replicate n die))] | |
63 (apply + ds))) | |
64 | |
65 (assert (= two-dice (dice 2))) | |
66 | |
67 (dice 3) | |
68 | |
69 | |
70 ;; Construct an empirical distribution from counters | |
71 | |
72 ; Using an ordinary counter: | |
73 (def dist1 | |
74 (normalize | |
75 (clojure.contrib.accumulators/add-items | |
76 clojure.contrib.accumulators/empty-counter | |
77 (for [_ (range 1000)] (rand-int 5))))) | |
78 | |
79 ; Or, more efficiently, using a counter that already keeps track of its total: | |
80 (def dist2 | |
81 (normalize | |
82 (clojure.contrib.accumulators/add-items | |
83 clojure.contrib.accumulators/empty-counter-with-total | |
84 (for [_ (range 1000)] (rand-int 5))))) | |
85 | |
86 | |
87 ;; The Monty Hall game | |
88 ;; (see http://en.wikipedia.org/wiki/Monty_Hall_problem for a description) | |
89 | |
90 ; The set of doors. In the classical variant, there are three doors, | |
91 ; but the code can also work with more than three doors. | |
92 (def doors #{:A :B :C}) | |
93 | |
94 ; A simulation of the game, step by step: | |
95 (domonad dist-m | |
96 [; The prize is hidden behind one of the doors. | |
97 prize (uniform doors) | |
98 ; The player make his initial choice. | |
99 choice (uniform doors) | |
100 ; The host opens a door which is neither the prize door nor the | |
101 ; one chosen by the player. | |
102 opened (uniform (disj doors prize choice)) | |
103 ; If the player stays with his initial choice, the game ends and the | |
104 ; following line should be commented out. It describes the switch from | |
105 ; the initial choice to a door that is neither the opened one nor | |
106 ; his original choice. | |
107 choice (uniform (disj doors opened choice)) | |
108 ] | |
109 ; If the chosen door has the prize behind it, the player wins. | |
110 (if (= choice prize) :win :loose)) | |
111 | |
112 | |
113 ;; Tree growth simulation | |
114 ;; Adapted from the code in: | |
115 ;; Martin Erwig and Steve Kollmansberger, | |
116 ;; "Probabilistic Functional Programming in Haskell", | |
117 ;; Journal of Functional Programming, Vol. 16, No. 1, 21-34, 2006 | |
118 ;; http://web.engr.oregonstate.edu/~erwig/papers/abstracts.html#JFP06a | |
119 | |
120 ; A tree is represented by two attributes: its state (alive, hit, fallen), | |
121 ; and its height (an integer). A new tree starts out alive and with zero height. | |
122 (def new-tree {:state :alive, :height 0}) | |
123 | |
124 ; An evolution step in the simulation modifies alive trees only. They can | |
125 ; either grow by one (90% probability), be hit by lightning and then stop | |
126 ; growing (4% probability), or fall down (6% probability). | |
127 (defn evolve-1 [tree] | |
128 (let [{s :state h :height} tree] | |
129 (if (= s :alive) | |
130 (choose 0.9 (assoc tree :height (inc (:height tree))) | |
131 0.04 (assoc tree :state :hit) | |
132 :else {:state :fallen, :height 0}) | |
133 (certainly tree)))) | |
134 | |
135 ; Multiple evolution steps can be chained together with m-chain, | |
136 ; since each step's input is the output of the previous step. | |
137 (with-monad dist-m | |
138 (defn evolve [n tree] | |
139 ((m-chain (replicate n evolve-1)) tree))) | |
140 | |
141 ; Try it for zero, one, or two steps. | |
142 (evolve 0 new-tree) | |
143 (evolve 1 new-tree) | |
144 (evolve 2 new-tree) | |
145 | |
146 ; We can also get a distribution of the height only: | |
147 (with-monad dist-m | |
148 ((m-lift 1 :height) (evolve 2 new-tree))) | |
149 | |
150 | |
151 | |
152 ;; Bayesian inference | |
153 ;; | |
154 ;; Suppose someone has three dice, one with six faces, one with eight, and | |
155 ;; one with twelve. This person throws one die and gives us the number, | |
156 ;; but doesn't tell us which die it was. What are the Bayesian probabilities | |
157 ;; for each of the three dice, given the observation we have? | |
158 | |
159 ; A function that returns the distribution of a dice with n faces. | |
160 (defn die-n [n] (uniform (range 1 (inc n)))) | |
161 | |
162 ; The three dice in the game with their distributions. With this map, we | |
163 ; can easily calculate the probability for an observation under the | |
164 ; condition that a particular die was used. | |
165 (def dice {:six (die-n 6) | |
166 :eight (die-n 8) | |
167 :twelve (die-n 12)}) | |
168 | |
169 ; The only prior knowledge is that one of the three dice is used, so we | |
170 ; have no better than a uniform distribution to start with. | |
171 (def prior (uniform (keys dice))) | |
172 | |
173 ; Add a single observation to the information contained in the | |
174 ; distribution. Adding an observation consists of | |
175 ; 1) Draw a die from the prior distribution. | |
176 ; 2) Draw an observation from the distribution of that die. | |
177 ; 3) Eliminate (replace by nil) the trials that do not match the observation. | |
178 ; 4) Normalize the distribution for the non-nil values. | |
179 (defn add-observation [prior observation] | |
180 (normalize-cond | |
181 (domonad cond-dist-m | |
182 [die prior | |
183 number (get dice die) | |
184 :when (= number observation) ] | |
185 die))) | |
186 | |
187 ; Add one observation. | |
188 (add-observation prior 1) | |
189 | |
190 ; Add three consecutive observations. | |
191 (-> prior (add-observation 1) | |
192 (add-observation 3) | |
193 (add-observation 7)) | |
194 | |
195 ; We can also add multiple observations in a single trial, but this | |
196 ; is slower because more combinations have to be taken into account. | |
197 ; With Bayesian inference, it is most efficient to eliminate choices | |
198 ; as early as possible. | |
199 (defn add-observations [prior observations] | |
200 (with-monad cond-dist-m | |
201 (let [n-nums #(m-seq (replicate (count observations) (get dice %)))] | |
202 (normalize-cond | |
203 (domonad | |
204 [die prior | |
205 nums (n-nums die) | |
206 :when (= nums observations)] | |
207 die))))) | |
208 | |
209 (add-observations prior [1 3 7]) |