gaussian process regression

gp.clj

(ns demo.gp
  (:use (incanter core stats)))
(set! *warn-on-reflection* true)


(def X (matrix [-4 -3 -1 0 2]))
(def Y (matrix [-2 0 1 2 -1]))
(def X* (matrix  (range -5 5 0.2040816)))


(defn k
  "Calculates the covariance matrix using simplified squared exponential function."
  [X1 X2]
  (let [Σ (for [i (range (count X1))
                j (range (count X2))]
            (exp (* -0.5 (pow (abs (- (sel X1 i 0)
                                      (sel X2 j 0)))
                              2))))]
    ;; build matrix representation
    (matrix (partition (count X2) Σ))))


;; covariance matrices
(def Kxx   (k X X))
(def Kxx*  (k X X*))
(def Kx*x  (k X* X))
(def Kx*x* (k X* X*))

;; correspond to equation (2.19) from book
(def Kxx-1 (solve Kxx))
(def f*   (mmult Kx*x Kxx-1 Y))
(def cov* (minus Kx*x* (mmult Kx*x Kxx-1 Kxx*)))



(defn mvrnorm
  "NOTE: incanter/sample-mvn uses cholesky, which requires positive semi def. So we use eigenvalue decomposition instead."
  [size & {:keys [mean sigma]}]
  (let [p (count mean)
        eS (decomp-eigenvalue sigma)
        X (matrix (sample-normal (* size p)) p)
        values (map (fn [x] (if (< x 0) 0 (sqrt x))) (:values eS))]

    (plus mean (mmult (:vectors eS)
                      (diag values)
                      (trans X)))))


(defn plot-samples
  "Sample form normal distributing with infered mean and covariance."
  []
  (let [plot (xy-plot X* f*)]
    (dotimes [i 5]
      (add-lines plot X* (mvrnorm 1 :mean f* :sigma cov*)))
    (add-points plot X Y)

    (view plot)))

based on the R code from http://www.jameskeirstead.ca/r/gaussian-process-regression-with-r/