Package {FastGP}

Posterior Sampling Using Elliptical Slice Sampling

Description

The ess function uses elliptical slice sampling to draw from the posterior distribution within a Bayesian model in which the prior is a Gaussian process/multivariate normal. If the covariance matrix corresponding to the prior is also Toeplitz, use ess_Toep for a faster version.

Elliptical slice sampling is a valid Markov chain Monte Carlo (MCMC) algorithm for posterior sampling when the prior is a Gaussian process/multivariate normal. The prior is assumed to have 0 mean, but a non-zero mean can be used by shfiting the location in log.lik, and shifting the location of the returned samples likewise.

Murray, I., Adams, R., & MacKay, D. (2010, March). Elliptical slice sampling. In Proceedings of the thirteenth international conference on artificial intelligence and statistics (pp. 541-548). JMLR Workshop and Conference Proceedings.

Usage

ess(log.lik,Y, Sig, N_mcmc,burn_in,N,flag)
ess_Toep(log.lik,Y, Sig, N_mcmc,burn_in,N)

Arguments

Author(s)

Giri Gopalan gopalan88@gmail.com

Examples

# See demo/FastGPdemo.r.
set.seed(123)
library(MASS)

N <- 100
sigma <- 20 #variance parameter in the covariance function
phi <- 100 #decay parameter for the squared-exponential kernel
x_seq <- seq(N)
Sig <- as.matrix(sigma*exp(-as.matrix(dist(x_seq))^2/phi)) + .001*diag(N) #test covariance function

A <- 1 #amplitude of the sin function used for our signal
T <- 5 #period of the sin function used as our signal

#A function that creates a signal, which can be toggled for multiple shapes. 
#In the present iteration we use a sin function with variable period.
signal <- function(t) {
  return(A*sin(t/T))
}
#Generate a copy of the signal on the time scale of 1...N
t <- mvrnorm(n = 1, mu = rep(0,N), Sigma = Sig , tol = 1e-6)
s <- signal(seq(1:N)+t)+rnorm(N,sd=sqrt(.001))

#log-lik function
log_lik <- function(w,s){
  return (rcpp_log_dmvnorm(S = .001*diag(1,N),mu=signal(seq(1:N)+w),s,istoep=TRUE))
}

mcmc_samples <- ess_Toep(log_lik, s, Sig, 1000, 250, N)
par(mfrow=c(1,1))
plot(colMeans(mcmc_samples), type="l")
points(t, type="l", lwd=3, col="red")
points(colMeans(mcmc_samples) + 2* apply(mcmc_samples, 2, sd), type="l", lty="dashed")
points(colMeans(mcmc_samples) - 2* apply(mcmc_samples, 2, sd), type="l", lty="dashed")

mcmc_samples <- ess(log_lik, s, Sig, 1000, 250, N, flag = TRUE)
par(mfrow=c(1,1))
plot(colMeans(mcmc_samples), type="l")
points(t, type="l", lwd=3, col="red")
points(colMeans(mcmc_samples) + 2* apply(mcmc_samples, 2, sd), type="l", lty="dashed")
points(colMeans(mcmc_samples) - 2* apply(mcmc_samples, 2, sd), type="l", lty="dashed")

mcmc_samples <- ess(log_lik, s, Sig, 1000, 250, N, flag = FALSE)
par(mfrow=c(1,1))
plot(colMeans(mcmc_samples), type="l")
points(t, type="l", lwd=3, col="red")
points(colMeans(mcmc_samples) + 2* apply(mcmc_samples, 2, sd), type="l", lty="dashed")
points(colMeans(mcmc_samples) - 2* apply(mcmc_samples, 2, sd), type="l", lty="dashed")

Matrix Operations Using Rcpp and RcppEigen

Description

Performs matrix operations using Rcpp and RcppEigen, and some custom C++ to handle symmetric and positive-definite Toeplitz matrices.

Usage

rcppeigen_invert_matrix(A)
rcppeigen_get_det(A)
rcppeigen_get_chol(A)
tinv(A)
tdet(A)
rcpp_distance(A, nrow, ncol)

Arguments

Details

Functions with "rcppeigen" directly call RcppEigen implementations of the associated functions, specifically for inverse, determinant, and lower triangular Cholesky factor. The function tinv computes an inverse and tdet returns the log-determinant of a symmetric and positive-definite Toeplitz matrix using methods from Trench and Durbin from "Matrix Computations" by Golub and Van Loan. Rcpp_distance computes the pairwise distances between rows of input matrix.

Author(s)

gopalan88@gmail.com

Examples

# See demo/FastGPdemo.R

N <- 5
sigma <- 1 #variance parameter in the covariance function
phi <- 5 #decay parameter for the squared-exponential kernel
x_seq <- seq(N)
Sig <- as.matrix(sigma*exp(-as.matrix(dist(x_seq))^2/phi)) + .001*diag(N) #test covariance function

rcppeigen_invert_matrix(Sig)
rcppeigen_get_det(Sig)
rcppeigen_get_chol(Sig)
tinv(Sig)
#NOTE: tdet returns the log of determinant
tdet(Sig)
rcpp_distance(matrix(seq(1,N^2),nrow=N), N, N)

Multivariate Normal Sampling and Log-Density Evaluation

Description

These functions allow for the sampling of and evaluation of the log-density of a multivariate normal vector.

Usage

rcpp_log_dmvnorm(S,mu,x, istoep)
rcpp_rmvnorm(n,S,mu)
rcpp_rmvnorm_Toep(n,S,mu)
rcpp_rmvnorm_stable(n,S,mu)

Arguments

Author(s)

Giri Gopalan gopalan88@gmail.com

Examples

#See demo/FastGPdemo.R, which includes timing comparisons
set.seed(123)

N <- 5
sigma <- 1 #variance parameter in the covariance function
phi <- 5 #decay parameter for the squared-exponential kernel
x_seq <- seq(N)
Sig <- as.matrix(sigma*exp(-as.matrix(dist(x_seq))^2/phi)) + .001*diag(N) #test covariance function
mu <- rep(0, N)

rcpp_log_dmvnorm(Sig,mu,rep(.1,N), istoep = TRUE)
rcpp_rmvnorm(n=1,Sig,mu)
rcpp_rmvnorm_Toep(n=1,Sig,mu)
rcpp_rmvnorm_stable(n=1,Sig,mu)