sharing weights really mixes between models

src/probabilistic_clojure/embedded/fit_poly.clj

@@ -24,11 +24,11 @@
     ^{:author "Nils Bertschinger"
       :doc "Part of probabilistic-clojure.embedded. Demo of fitting polynomials."}
   probabilistic-clojure.embedded.fit-poly
-  (:use [probabilistic-clojure.embedded.api :only (det-cp gv trace-failure cond-data metropolis-hastings-sampling)]
+  (:use [probabilistic-clojure.embedded.api :only (det-cp gv trace-failure cond-data metropolis-hastings-sampling def-prob-cp)]
 	[probabilistic-clojure.utils.stuff :only (indexed)]
 	[probabilistic-clojure.utils.sampling :only (normalize sample-from)]
 	[probabilistic-clojure.embedded.choice-points
-	 :only (gaussian-cp discrete-cp log-pdf-discrete)]
+	 :only (gaussian-cp discrete-cp log-pdf-discrete *gaussian-proposal-sdev*)]
 
 	[incanter.core   :only (view)]
 	[incanter.charts :only (histogram add-lines xy-plot scatter-plot)]
@@ -43,17 +43,27 @@ f(x) = coeffs[0] + coeffs[1]*x + ... + coeffs[k]*x^(k-1)"
 
   (reduce + (map * coeffs (iterate (partial * x) 1))))
 
-(def test-poly [-1 1 0.5 -0.2])
+(def test-poly [-1 1 -0.7 0.3])
 
 (defn uni-rand [low high]
   (+ low (rand (- high low))))
 
 (def demo-data
-     (repeatedly 42
+     (repeatedly 17
 		 (fn []
 		   (let [x (uni-rand -5 5)
 			 y (poly x test-poly)]
-		     [x (sample-normal 1 :mean y :sd 0.5)]))))
+		     [x (sample-normal 1 :mean y :sd 2)]))))
+
+;; special gaussian-cp which initially returns mu
+(def-prob-cp gaussian0-cp [mu sdev]
+  :sampler [] mu
+  :calc-log-lik [x] (Math/log (pdf-normal x :mean mu :sd sdev))
+  :proposer [old-x] (let [proposal-sd *gaussian-proposal-sdev*
+			  new-x (sample-normal 1 :mean old-x :sd proposal-sd)]
+		      [new-x
+		       (Math/log (pdf-normal new-x :mean old-x :sd proposal-sd))
+		       (Math/log (pdf-normal old-x :mean new-x :sd proposal-sd))]))
 
 (defn poly-fit-fixed
   [rank data]
@@ -83,8 +93,8 @@ f(x) = coeffs[0] + coeffs[1]*x + ... + coeffs[k]*x^(k-1)"
 	   			(for [rank rank-range]
 	   			  [rank (doall
 	   				 (for [r (range (inc rank))]
-	   				   (gv (gaussian-cp [:coeff rank r]
-	   						    0.0 5.0))))])))
+	   				   (gv (gaussian0-cp [:coeff rank r]
+							     0.0 5.0))))])))
 	   obs-noise (det-cp :noise
 	   		     (into {} (for [rank rank-range]
 	   				[rank (Math/abs (gv (gaussian-cp [:noise rank]
@@ -105,4 +115,39 @@ f(x) = coeffs[0] + coeffs[1]*x + ... + coeffs[k]*x^(k-1)"
 	graph (scatter-plot (map first demo-data) (map second demo-data))]
     (doseq [r ranks]
       (add-lines graph xs (map #(poly % (get (second fitted) r)) xs)))
-    (view graph)))
+    (view graph)))
+
+(defn poly-fit-shared
+  ([data] (poly-fit-shared [1] data))
+  ([rank-range data]
+     (let [n (count rank-range)
+	   rank (if (= n 1)
+		  (det-cp :rank (first rank-range))
+		  (discrete-cp :rank (zipmap rank-range (repeat (/ 1 n)))))
+	   coeffs (det-cp :coeffs
+			  ;; one set of coeffs which are shared by all models
+			  ;; and removed if unused
+			  (doall
+			   (for [r (range (inc (gv rank)))]
+			     (gv (gaussian0-cp [:coeff r]
+					       0.0 5.0)))))
+	   obs-noise (det-cp :noise
+			     (Math/abs (gv (gaussian-cp [:noise rank]
+							0.0 5.0))))]
+       (doseq [[i [x y]] (indexed data)]
+	 (cond-data (gaussian-cp [:obs i]
+				 (poly x (gv coeffs))
+				 (gv obs-noise))
+		    y))
+       (det-cp :fit
+	       [(gv rank) (gv coeffs) (gv obs-noise)]))))
+
+(defn poly-demo-shared [ranks data]
+  (let [xs (range -5 5 0.05)
+
+	fitted (last (take 75000 (metropolis-hastings-sampling
+				 (fn [] (poly-fit-shared ranks demo-data)))))
+	graph (scatter-plot (map first demo-data) (map second demo-data))]
+    (doseq [r ranks]
+      (add-lines graph xs (map #(poly % (take r (second fitted))) xs)))
+    (view graph)))