# Gaussian process regression functions. # CREATE THE COVARIANCE MATRIX FOR TRAINING RESPONSES. The arguments are # the matrix of training inputs, the covariance function (taking an input # vector and a matrix as arguments), and the vector of hyperparameters. # The first element in the hyperparameter vector is the residual # variance, which is added to the diagonal of the covariance matrix, # but which is ignored by the covariance function. The result is the # matrix of covariances for responses in training cases. gp_cov_matrix <- function (x_train, covf, hypers) { n <- nrow(x_train) C <- matrix(NA,n,n) for (i in 1:n) { C[i,] <- covf(x_train[i,],x_train,hypers) } diag(C) <- diag(C) + hypers[1]^2 C } # PREDICT THE RESPONSE FOR TEST CASES. The arguments are the matrix of # inputs for training cases, the vector of responses for training cases, # the matrix of inputs for test cases, and the covariance function and # hyperparameter vector as described for gp_cov_matrix. The result is # a vector of predictions of the responses for test cases (the predictive # mean). gp_predict <- function (x_train, y_train, x_test, covf, hypers) { n <- nrow(x_train) C <- gp_cov_matrix(x_train,covf,hypers) b <- solve(C,y_train) r <- numeric(nrow(x_test)) for (i in 1:nrow(x_test)) { k <- covf(x_test[i,],x_train,hypers) r[i] <- k %*% b } r } # FIND THE LOG LIKELIHOOD BASED ON RESPONSES IN THE TRAINING SET. # The arguments are the matrix of inputs for training cases, the # vector of responses for training cases, and the covariance function # and hyperparameter vector as described for gp_cov_matrix. The # result is the log of the probability density for the training # responses (with all normalizing factors included). gp_log_likelihood <- function (x_train, y_train, covf, hypers) { n <- nrow(x_train) C <- gp_cov_matrix(x_train,covf,hypers) b <- solve(C,y_train) - (n/2)*log(2*pi) - determinant(C)$modulus/2 - y_train%*%b/2 } # FIND HYPERPARAMETER VALUES MAXIMIZING LOG PROBABILITY. Arguments # are the matrix of training inputs, the vector of training responses, # the covariance function (as for gp_cov_matrix), and a vector of # initial values for the hyperparameters. (Any additional arguments are # passed on to the nlm function.) The result is a vector of values # for the hyperparameters that should usually at least a locally maximize # the likelihood. # # The maximization is done in terms of the log likelihood. It is assumed # that hyperparameter values are squared before use, so that negative # values are equivalent to the corresponding positive value. gp_find_hypers <- function (x_train, y_train, covf, hypers0, ...) { m <- nlm ( function (h) - gp_log_likelihood(x_train,y_train,covf, h), hypers0, ... ) cat ("Maximum log likelihood:", -m$minimum, "\n") abs(m$estimate) }