Moved examples directory out of the src directory.

examples/bayes.clj

@@ -0,0 +1,156 @@
+
+;;; examples/bayes.clj -- Bayesian estimation library for Clojure
+
+;; by David Edgar Liebke http://incanter.org
+;; March 11, 2009
+
+;; Copyright (c) David Edgar Liebke, 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.
+
+;; CHANGE LOG
+;; March 11, 2009: First version
+
+
+
+(ns examples.bayes 
+  (:use (incanter core stats)))
+
+
+
+(defn bayes-regression-noref [N x y]
+"
+  This function implments the Gibbs sampling example using full conditional in OLS
+  from Scott Lynch book 'Introduction to Applied Bayesian Statistics in the Social
+  Sciences (page 171). This version is purely function with no immutability.
+"
+  (let [lm (linear-model y x :intercept false)
+        pars (trans (:coefs lm))
+        xtxi (solve (mmult (trans x) x))
+        nx (ncol x)
+        shape (/ (- (nrow x) (ncol x)) 2)
+       ]
+      (loop [
+             coefs (transient [[0 0 0 0 0 0 0 0 0]])
+             variances (transient [1]) 
+             i 0
+            ]
+        (if (= i N)
+          {:coef (matrix (persistent! coefs)) :var (persistent! variances)}
+          (let [
+                b (to-list (plus pars (mmult (trans (sample-normal nx)) 
+                                             (decomp-cholesky (mult xtxi (variances i))))))
+                resid (minus y (mmult x b))
+                s2 (/ 1 (sample-gamma 1 :shape shape :rate (mult (mmult (trans resid) resid) 0.5) ))
+               ]
+            (recur (conj! coefs b) (conj! variances s2) (inc i)))))))
+
+
+
+
+(defn bayes-regression-full [N x y]
+"
+  This function implments the Gibbs sampling example using full conditional in OLS
+  from Scott Lynch book 'Introduction to Applied Bayesian Statistics in the Social
+  Sciences (page 171). This version uses immutability (i.e. references)
+"
+  (let [lm (linear-model y x :intercept false)
+        pars (trans (:coefs lm))        
+        xtxi (solve (mmult (trans x) x))
+        nx (ncol x)
+        b (ref [[0 0 0 0 0 0 0 0 0]])
+        s2 (ref [1])
+        resid (ref 0)
+        shape (/ (- (nrow x) (ncol x)) 2)
+       ]
+    (do
+      (dotimes [i N]
+        (dosync
+          (alter b conj 
+            (to-list (plus pars (mmult (trans (sample-normal nx)) (decomp-cholesky (mult xtxi (@s2 i)))))))
+          (ref-set resid (minus y (mmult x (@b (inc i)))))
+          (alter s2 conj (/ 1 (sample-gamma 1 :shape shape :rate (mult (mmult (trans @resid) @resid) 0.5) )))))
+        ;; return a map with the estimated coefficients and variances
+       {:coef (matrix @b) :var @s2})))
+
+
+
+
+(defn bayes-regression [N x y]
+"
+  This function implments the Gibbs sampling example using the composition method
+  in OLS from Scott Lynch book 'Introduction to Applied Bayesian Statistics in the 
+  Social Sciences (page 173)
+"
+  (let [lm (linear-model y x :intercept false)
+        pars (:coefs lm)
+        xtxi (solve (mmult (trans x) x))
+        resid (:residuals lm)
+        shape (/ (- (nrow x) (ncol x)) 2)
+        rate (mult 1/2 (mmult (trans resid) resid))
+        s-sq (div 1 (sample-gamma N :shape shape :rate rate))]
+    ;; return a map with the estimated coefficients and variances
+    {:coef 
+      (matrix 
+        ;(pmap ;; run a parallel map over the values of s-sq
+        (map
+          (fn [s2] 
+            (to-list (plus (trans pars)
+                (mmult (trans (sample-normal (ncol x))) 
+                  (decomp-cholesky (mult s2 xtxi))))))
+          (to-list (trans s-sq)))) 
+     :var s-sq}))
+
+
+
+
+(defn bayes-regression-mh [N x y]
+"
+  This function implments the Gibbs sampling example using Metropolis Hastings
+  in OLS from Scott Lynch book 'Introduction to Applied Bayesian Statistics in the 
+  Social Sciences (page 168)
+"
+  (let [ lm (linear-model y x :intercept false)
+         b-scale (to-list (div (sqrt (:std-errors lm)) 2))
+         s2-scale (/ (sd (mult (:residuals lm) (div (dec (nrow x)) (- (nrow x) (ncol x))))) 2)
+         post-fn (fn [x y b s-sq] 
+                  (let [resid (minus y (mmult x b))]
+                    (plus (mult -1157.5 (log s-sq)) 
+                          (mult (div -0.5 s-sq) (mmult (trans resid) resid))))) 
+        reject? (fn [x y cand-b cand-s2 old-b old-s2] 
+                  (< (- (post-fn x y cand-b cand-s2) 
+                        (post-fn x y old-b old-s2))
+                     (log (rand))))
+        ncol-x (ncol x)
+       ]
+    (loop [
+           coefs (transient [[0 0 0 0 0 0 0 0 0]])
+           variances (transient [1])
+           i 0
+          ]
+      (if (= i N)
+          {:coef (matrix (persistent! coefs)) :var (persistent! variances)}
+          (let [
+                old-b (coefs i)
+                old-s2 (variances i)
+                cand-b (into [] (map #(+ (sample-normal 1 :mean 0 :sd %1) %2) b-scale old-b))
+                new-b (loop [b old-b j 0]
+                          (if (= j ncol-x)
+                            b
+                            (recur 
+                              (if (reject? x y (assoc b j (cand-b j)) old-s2 old-b old-s2) 
+                                b
+                                (assoc b j (cand-b j)))
+                              (inc j)))) 
+                cand-s2 (+ (sample-normal 1 :mean 0 :sd s2-scale) old-s2)
+                new-s2 (if (or (< cand-s2 0) (reject? x y new-b cand-s2 new-b old-s2)) 
+                         old-s2 
+                         cand-s2)
+               ]
+        (recur (conj! coefs new-b) (conj! variances new-s2) (inc i)))))))
+
+

examples/benfry/locations.tsv

@@ -0,0 +1,50 @@
+AL	439	270
+AK	94	325
+AZ	148	241
+AR	368	247
+CA	56	176
+CO	220	183
+CT	576	120
+DE	556	166
+FL	510	331
+GA	478	267
+HI	232	380
+ID	143	101
+IL	405	168
+IN	437	165
+IA	357	147
+KS	302	194
+KY	453	203
+LA	371	302
+ME	595	59
+MD	538	162
+MA	581	108
+MI	446	120
+MN	339	86
+MS	406	274
+MO	365	197
+MT	194	61
+NE	286	151
+NV	102	157
+NH	580	89
+NJ	561	143
+NM	208	245
+NY	541	107
+NC	519	221
+ND	283	65
+OH	472	160
+OK	309	239
+OR	74	86
+PA	523	144
+RI	589	117
+SC	506	251
+SD	286	109
+TN	441	229
+TX	291	299
+UT	154	171
+VT	567	86
+VA	529	189
+WA	92	38
+WV	496	178
+WI	392	103
+WY	207	125

examples/benfry/random.tsv

@@ -0,0 +1,50 @@
+AL	0.1
+AK	-5.3
+AZ	3
+AR	7
+CA	11
+CO	1.5
+CT	-6.7
+DE	-4
+FL	9
+GA	2
+HI	-3.3
+ID	6.6
+IL	7.2
+IN	7.1
+IA	6.9
+KS	6
+KY	1.8
+LA	7.5
+ME	-4
+MD	0.1
+MA	-6
+MI	1.7
+MN	-2
+MS	-4.4
+MO	-2
+MT	1.0
+NE	1.2
+NV	1.6
+NH	0.5
+NJ	0.2
+NM	8.8
+NY	1.4
+NC	9.7
+ND	5.4
+OH	3.2
+OK	6
+OR	-4
+PA	-7
+RI	-2
+SC	1
+SD	6
+TN	5
+TX	-3.4
+UT	2.3
+VT	4.8
+VA	3
+WA	2.2
+WV	5.4
+WI	3.1
+WY	-6

examples/blog/annotations.clj

@@ -0,0 +1,64 @@
+;; from the following blog post
+;; http://incanter-blog.org/2009/06/07/annotating-incanter-plots/
+
+(use '(incanter core charts))
+
+(def x (range (* -2 Math/PI) (* 2 Math/PI) 0.01))
+(def plot (xy-plot x (sin x)))
+(view plot)
+
+;; annotate the plot
+(doto plot
+  (add-pointer (- Math/PI) (sin (- Math/PI)) 
+               :text "(-pi, (sin -pi))")
+  (add-pointer Math/PI (sin Math/PI) 
+               :text "(pi, (sin pi))" :angle :ne)
+  (add-pointer (* 1/2 Math/PI) (sin (* 1/2 Math/PI)) 
+               :text "(pi/2, (sin pi/2))" :angle :south))
+
+
+;; try the different angle options
+(doto plot
+  (add-pointer 0 0 :text "north" :angle :north)
+  (add-pointer 0 0 :text "nw" :angle :nw)
+  (add-pointer 0 0 :text "ne" :angle :ne)
+  (add-pointer 0 0 :text "west" :angle :west)
+  (add-pointer 0 0 :text "east" :angle :east)
+  (add-pointer 0 0 :text "south" :angle :south)
+  (add-pointer 0 0 :text "sw" :angle :sw)
+  (add-pointer 0 0 :text "se" :angle :se))
+
+
+
+;; PCA chart example
+(use '(incanter core stats charts datasets))
+;; load the iris dataset
+(def iris (to-matrix (get-dataset :iris)))
+;; run the pca
+(def pca (principal-components (sel iris :cols (range 4))))
+;; extract the first two principal components
+(def pc1 (sel (:rotation pca) :cols 0))
+(def pc2 (sel (:rotation pca) :cols 1))
+
+;; project the first four dimension of the iris data onto the first
+;; two principal components
+(def x1 (mmult (sel iris :cols (range 4)) pc1)) 
+(def x2 (mmult (sel iris :cols (range 4)) pc2)) 
+
+;; now plot the transformed data, coloring each species a different color
+(def plot (scatter-plot x1 x2 
+            :group-by (sel iris :cols 4)
+            :x-label "PC1" :y-label "PC2" :title "Iris PCA"))
+
+(view plot)
+
+
+;; add some text annotations
+(doto plot
+  (add-text -2.5 -6.5 "Setosa")
+  (add-text -5 -5.5 "Versicolor") 
+  (add-text -8 -5.5 "Virginica"))
+
+;; put box around the first group
+(add-polygon plot [[-3.2 -6.3] [-2 -6.3] [-2 -3.78] [-3.2 -3.78]])
+

examples/blog/bayesian_multinomial.clj

@@ -0,0 +1,50 @@
+;; from the following blog post:
+;; http://incanter-blog.org/2009/07/01/bayes-multinomial/
+
+;; Bayesian inference of multinomial distribution parameters
+
+(use '(incanter core stats bayes charts))
+
+
+(def y [727 583 137])
+
+(div y 1447.) ;; (0.502 0.403 0.095)
+
+
+(def  theta (sample-multinomial-params 1000 y))
+(def theta1 (sel theta :cols 0))
+(def theta2 (sel theta :cols 1))
+(def theta3 (sel theta :cols 2))
+
+;; view means, 95% CI, and histograms of the proportion parameters
+(mean theta1) ;; 0.502
+(sd theta1) ;; 0.0129
+(quantile theta1 :probs [0.025 0.975]) ;; (0.476 0.528)
+(view (histogram theta1
+                 :title "Bush, Sr. Support"))
+
+(mean theta2) ;; 0.403
+(sd theta2) ;; 0.013
+(quantile theta2 :probs [0.025 0.975]) ;; (0.376 0.427)
+(view (histogram theta2
+                 :title "Dukakis Support"))
+
+(mean theta3) ;; 0.095
+(sd theta3) ;; 0.008
+(quantile theta3 :probs [0.025 0.975]) ;; (0.082 0.111)
+(view (histogram theta3
+                 :title "Other/Undecided"))
+
+;; view  a histogram of the difference in proportions between the first
+;; two candidates
+(view (histogram (minus theta1 theta2) 
+                 :title "Bush, Sr. and Dukakis Diff"))
+
+
+
+;; sample the proportions directly from a Dirichlet distribution
+(def alpha [1 1 1])
+(def  props (sample-dirichlet 1000 (plus y alpha)))
+
+
+