I am looking to compute the covariance matrix of an inferred Gaussian process in R. Below I outline how I would do this manually, but I realize that the kernlab
package provides the function gausspr
, which would be preferable for many reasons. Unfortunately, while I see how to get the expected value (predict.gausspr
) using the package, I do not see how to get the associated covariance matrix.
To provide a minimal working example for comparison, I loosely follow Rasmussen & Williams (2006) Chapter 2:
Consider we have the observed x,y
points and x
values where we desire predicted y
values:
obs <- data.frame(x = c(-4, -3, -1, 0, 2),
y = c(-2, 0, 1, 2, -1))
x_predict <- seq(-5,5,len=50)
Use a radial basis kernel:
SE <- function(Xi,Xj, l=1) exp(-0.5 * (Xi - Xj) ^ 2 / l ^ 2)
cov <- function(X, Y) outer(X, Y, SE)
Calculate mean and covariance matrix:
sigma.n <- 0.3
cov_xx_inv <- solve(cov(obs$x, obs$x) + sigma.n^2 * diag(1, length(obs$x)))
Ef <- cov(x_predict, obs$x) %*% cov_xx_inv %*% obs$y
Cf <- cov(x_predict, x_predict) - cov(x_predict, obs$x) %*% cov_xx_inv %*% cov(obs$x, x_predict)
kernlab approach
I think I see how I get the equivalent expected values in kernlab:
library(kernlab)
gp <- gausspr(obs$x, obs$y, kernel="rbfdot", kpar=list(sigma=1/(2*l^2)), fit=FALSE, scaled=FALSE, var=0.8)
Ef_kernlab <- predict(gp, x_predict)
but I don't see how I get the associated covariances? Can the covariance matrix be reconstructed by knowing alpha
?
Ideally there is a solution that doesn't involve explicitly inverting a matrix again, since that inverse was already calculated in determining the alpha
values during the fit done by gausspr...
There are many cases where it would be nice to have access to the resulting Gaussian process, such as in generating plots as in Rasmussen and Williams:
require(ggplot2)
dat <- data.frame(x=x_predict, y=(Ef), ymin=(Ef-2*sqrt(diag(Cf))), ymax=(Ef+2*sqrt(diag(Cf))))
ggplot(dat) +
geom_ribbon(aes(x=x,y=y, ymin=ymin, ymax=ymax), fill="grey80") + # Var
geom_line(aes(x=x,y=y), size=1) + #MEAN
geom_point(data=obs,aes(x=x,y=y)) + #OBSERVED DATA
scale_y_continuous(lim=c(-3,3), name="output, f(x)") + xlab("input, x")