# GP.R - Gaussian process regression functions.
#
# R. M. Neal, April 2006


# CREATE THE COVARIANCE MATRIX FOR TRAINING RESPONSES.  The arguments
# are the matrix of training inputs, the covariance funcition (taking
# two input vectors 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)
  { for (j in i:N)
    { C[i,j] <- covf(x.train[i,],x.train[j,],hypers)
      C[j,i] <- C[i,j]
    }
  }

  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 the the response 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)

  k <- numeric(N)
  r <- numeric(nrow(x.test))

  for (i in 1:nrow(x.test))
  { for (j in 1:N)
    { k[j] <- covf(x.test[i,],x.train[j,],hypers)
    }
    r[i] <- k %*% b
  }

  r
}


# FIND THE LOG PROBABILITY DENSITY FOR 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.prob <- 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.  The result is a vector of
# values for the hyperparameters that should usually at least a
# locally maximize the probability density of the training responses.
#
# The maximization is done in terms of the log probability density,
# with hyperparameter h transformed to sqrt(h-0.01) in order to avoid
# problems with hyperparameters that are negative or near zero (less
# than 0.01).

gp.find.hypers <- function (x.train, y.train, covf, hypers0)
{
  m <- nlm (
    function (h) - gp.log.prob(x.train,y.train,covf, h^2+0.01), 
    sqrt(hypers0-0.01) )

  m$estimate^2 + 0.01
}