Clojureを用いた確率的プログラミング(Probabilistic Programming)

> lein new app gr-test(任意のプロジェクト名)

(defproject gr-test "0.1.0-SNAPSHOT"
  :description "FIXME: write description"
  :url "http://example.com/FIXME"
  :license {:name "EPL-2.0 OR GPL-2.0-or-later WITH Classpath-exception-2.0"
            :url "https://www.eclipse.org/legal/epl-2.0/"}
  :dependencies [[org.clojure/clojure "1.11.1"]]
  :main ^:skip-aot gr-test.core
  :target-path "target/%s"
  :plugins [[org.clojars.benfb/lein-gorilla "0.6.0"]]
  :profiles {:uberjar {:aot :all
                       :jvm-opts ["-Dclojure.compiler.direct-linking=true"]}})

>lein gorilla
………………………………………
Gorilla-REPL: 0.6.0
Unable to reach update server.
Started nREPL server on port 49989
Running at http://127.0.0.1:49991/worksheet.html .
Ctrl+C to exit.

(defproject gr-anglican "0.1.0-SNAPSHOT"
  :description "FIXME: write description"
  :url "http://example.com/FIXME"
  :license {:name "EPL-2.0 OR GPL-2.0-or-later WITH Classpath-exception-2.0"
            :url "https://www.eclipse.org/legal/epl-2.0/"}
  :dependencies [[org.clojure/clojure "1.11.1"]
                 [anglican "1.1.0"]
                 [net.mikera/core.matrix "0.62.0"]]
  :main ^:skip-aot gr-anglican.core
  :target-path "target/%s"
  :plugins [[org.clojars.benfb/lein-gorilla "0.6.0"]]
  :profiles {:uberjar {:aot :all
                       :jvm-opts ["-Dclojure.compiler.direct-linking=true"]}})

(defdist gda [pi feature-dists]

  [class-dist (discrete pi)
   C (count pi)]

  (sample* [this]
           (let [my-class (sample* class-dist)
                 my-features (sample* (nth feature-dists my-class))]
             [my-class my-features]))

  (observe* [this [my-features my-class]]
            (if (not (nil? my-class))
              (+
                (observe* class-dist my-class)
                (observe* (nth feature-dists my-class) my-features))

              (reduce log-sum-exp
                      (map

                        #(+'
                           (observe* class-dist %)
                           (observe* (nth feature-dists %) my-features))

                        (range C))))))

(defn classification-log-probabilities [gda-dist feature]
  (map
    (fn [c] (observe* gda-dist [feature c]))
    (range (:C gda-dist))))

(defn classification-probabilities [gda-dist feature]
  (let [unnormalized (classification-log-probabilities gda-dist feature)
        sum (reduce log-sum-exp unnormalized)
        normalized (map #(- % sum) unnormalized)]
    (map #(Math/exp %) normalized)))



(defn index-of-maximum [lst]
  (second
    (apply max-key first
           (map
             (fn [x y] [x y])
             lst
             (range)))))

(defn classify-feature [gda-dist feature]
  (let [unnormalized (classification-log-probabilities gda-dist feature)]
    (index-of-maximum unnormalized)))

(def feature-dist-1 (mvn [0 0] [[4 3] [3 4]]))
(def feature-dist-2 (mvn [0 10] [[10 -3] [-3 5]]))
(def feature-dist-3 (mvn [10 10] [[6 0] [0 8]]))

(def pi [1 1 1])

(def my-gda (gda pi (list feature-dist-1 feature-dist-2 feature-dist-3))) ;;training data

(def x [5 5]) ;; test data

(def probs (map #(Math/round (* 100 %)) (classification-probabilities my-gda x)))

(plot/bar-chart '(0 1 2) probs)
(reduce str (map #(str (Math/round (* 100 %)) "% ") (classification-probabilities my-gda x)))
;; "43% 41% 15% "
(str "The feature " x " has classification " (classify-feature my-gda x))
;;The feature [5 5] has classification 0"

;; HYPERPARAMETERS (for priors):

(def nu [1 1 1])
(def alpha [1 1])


;; PRIORS:

(def mu-priors (mvn nu (mat/identity-matrix 3)))
(def sigma (mat/identity-matrix 3))
(def pi-prior (dirichlet alpha))


(defquery gda-trainer [data x]
  (declare :primitive mat/zero-matrix mat/to-vector mat/identity-matrix gda classification-probabilities)
  (let [N (count data)

        ;; sample parameters

        mu1 (sample mu-priors)
        mu2 (sample mu-priors)
        pi (sample pi-prior)

        ;; construct GDA distribution
        
        feature-dists (list (mvn mu1 sigma) (mvn mu2 sigma))
        gda-dist (gda pi feature-dists)]

    
    ;; condition on data
    (loop [n 0]
      (if (< n N)
        (do
          (observe gda-dist (nth data n))
          (recur (+ n 1)))))
    
    ;; now project the type of x
    (classification-probabilities gda-dist x)))

(def dataset [
               [[3 3 3] 1]
               [[1 1 1]]
               [[0 0 0]]
               [[5 5 5] 1]
               [[4 4 4] 1]
               [[-1 -1 -1]]
               [[-2 0 0]]

               ])

(def x [1.5 1 1])

(def S 100000)
(def particles 100)
(def nodes 10)

(def samples (take S (drop S 
                           (doquery :ipmcmc gda-trainer [dataset x]
                                    :number-of-particles particles :number-of-nodes nodes))))

(def weights (map :log-weight samples))
(def projections (map :result samples))

(take-last 10 projections)

;;((0.538848314823797 0.46115168517620253) (0.538848314823797 0.46115168517620253) 
   (0.538848314823797 0.46115168517620253) (0.538848314823797 0.46115168517620253) 
   (0.538848314823797 0.46115168517620253) (0.9837451961913958 0.016254803808604182) 
   (0.9837451961913958 0.016254803808604182) (0.9837451961913958 0.016254803808604182) 
   (0.9837451961913958 0.016254803808604182) (0.9837451961913958 0.016254803808604182))

(defn calculate-mean [series]
  (empirical-mean
    (partition 2
               (interleave series weights))))

;; this will work better
(def percentages (map #(Math/round (* 100 %)) (calculate-mean projections)))

(plot/bar-chart (range (count percentages)) percentages)

(reduce str (map #(str % "% ") percentages))
;;"89% 11%"