gpr.R

# algorithm 2.1 in Rasmussen's book Gaussian Processes for Machine Learning
# input: X, y, k, s2, xs

# read rawdata generated by gpml-matlab-v3.6-2015-07-07/doc/index.html
data <- read.csv("gprdata.csv", header = TRUE, sep = ",")
X <- data$x
y <- data$y
n <- 20

# fit a simple gpr with default cov function and hyper parameters
# hyper for cov
l <- 1
sigmaf <- 1
# hyper for likli
s2 <- 0.1^2
# test data
xss <- seq(-1.9, 1.9, len = 101)

# cov functions square exponential
cov <- function(X, Y) {
  SE <- function(xi, xj) sigmaf ^ 2 * exp(-0.5 * (xi - xj) ^ 2 / l ^ 2)
  outer(X, Y, SE)
}

# line 2 in algorithm 2.1
# in R, chol calcs a higher triangular instead of a lower in the book
L <- t(chol(cov(X, X) + s2 * diag(n)))
# line 3 in algorithm 2.1
alpha <- solve(t(L), solve(L, y))

# line 4 to 6 in algorithm 2.1, repeat for each xs

CalcMean <- function(xs) crossprod(cov(X, xs), alpha)
CalcV <- function(xs) cov(xs, xs) - crossprod(solve(L, cov(X, xs)))
fsmean <- sapply(xss, FUN = CalcMean)
fsv <- sapply(xss, FUN = CalcV)

# line 7, compute negativce log margical likelihood
logmarg <- -0.5 * crossprod(y, alpha) - sum(log(diag(L)))- n / 2 * log(2 * pi)

# output
ysmean <- fsmean
ysv <- fsv + s2
logmarg

# plotting
predlb <- ysmean - 2 * sqrt(ysv)
predub <- ysmean + 2 * sqrt(ysv)

# a plot that resembles one produced by Matlab gpml package
windows()
plot(xss, ysmean, xlim = c(-2, 2.5), ylim = c(-1, 5))
polygon(c(xss, rev(xss)), c(predub, rev(predlb)), col = "gray85")
lines(xss, ysmean, col = "red")
points(X, y, pch = "+", col = "tan1")