Created
September 24, 2012 01:08
-
-
Save pedroteixeira/3773682 to your computer and use it in GitHub Desktop.
gaussian process regression
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
(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))) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
based on the R code from http://www.jameskeirstead.ca/r/gaussian-process-regression-with-r/