R/rcpp_mvn.R

In FastGP: Efficiently Using Gaussian Processes with Rcpp and RcppEigen

#R functions to deal with multivariate normals efficiently using Rcpp and RcppEigen functions

#draw n samples from a multivariate normal distribution with covariance matrix S and mean mu
rcpp_rmvnorm <- function(n,S,mu)
{
  chol_S <- rcppeigen_get_chol(S)
  m <- dim(chol_S)[1]
  mat <- matrix(rnorm(n*m),nrow=m,ncol=n)
  return (t(chol_S%*%mat+mu))
}
#faster version using O(n^2) Cholesky when S is a Toeplitz covariance matrix
rcpp_rmvnorm_Toep <- function(n,S,mu)
{
  tchol_out <- toep_chol_inv(S[1,])
  A <- tchol_out[[1]]
  D <- tchol_out[[2]]
  m <- dim(A)[1]
  mat <- sqrt(D) * matrix(rnorm(n*m), nrow = m, ncol = n)
  V <- forwardsolve(A, mat)
  return(t(V + mu))
}
#draw n samples from a multivariate normal distribution with covariance matrix S
#and mean mu with a more stable version of Cholesky
rcpp_rmvnorm_stable <- function(n,S,mu)
{
  if(length(mu) == 1){
    return(rnorm(n,mean=mu,sd=sqrt(S)))
  }
  else{
    chol_S <- rcppeigen_get_chol_stable(S)
    diag_s <- rcppeigen_get_chol_diag(S)
    diag_s[which(diag_s < 0)] <- 0
    diag_s <- diag(sqrt(diag_s))
    m <- dim(chol_S)[1]
    mat <- matrix(rnorm(n*m),nrow=m,ncol=n)
    return (t(chol_S%*%diag_s%*%mat+mu))
  }
}
#invert a symmetric, positive-definite, Toeplitz matrix 
#can be handy for GPs where grid points are evenly spaced
tinv <- function(A){
  M <- A/(A[1,1])
  N <- dim(M)[1]
  r <- M[1,2:N]
  y <- durbin(r,N-1)
  return(trench(r,y,N)/A[1,1])
}
#calculate log-determinant of symmetric, positive-definite, Toeplitz matrix
tdet <- function(A){
  r <- A[1,]
  return(toep_logdet(r))
}
#evaluate the log density of a multivariate normal distribution with covariance
#S and mean mu at point x, with a flag for whether the covariance matrix is Toeplitz
rcpp_log_dmvnorm <- function(S,mu,x, istoep)
{
  n <- length(x)
  z <- x-mu
  if (istoep == TRUE)
  {
    tchol_out <- toep_chol_inv(S[1,])
    A <- tchol_out[[1]]
    D <- tchol_out[[2]]
    log_det <- sum(log(D))
    w <- drop(A %*% z)
    q <- as.numeric(crossprod(w, w / D))
  }
  else
  {
    diag <- rcppeigen_get_diag(S)
    log_det <- 2*sum(log(diag))
    chol_S <- rcppeigen_get_chol(S)
    v <- forwardsolve(chol_S, z)
    q <- sum(v^2)
  }
  return((-n/2)*log(2*pi)-(1/2)*log_det-(1/2)*q)
}