# $Id: cass.R 648 2010-03-12 02:46:56Z mthompson $ # Sample Code for Covariance-Adaptive Slice Sampling # (c) 2010 Madeleine Thompson # # This file is intended as a demonstration. You may use it however # you like, entirely at your own risk. # For information about covariance-adaptive slice sampling, see: # # Madeleine Thompson and Radford M. Neal, "Covariance-Adaptive Slice # Sampling," University of Toronto Department of Statistics Technical # Report 1002, 2010. # # This technical report can be found at: # # http://www.cs.utoronto.ca/~radford/ # # DCHUD from LINPACK is included with this file. The original comes from: # # http://www.netlib.org/linpack # Building chud.so # ---------------- # # covariance.match.sample() can use LINPACK's DCHUD function to perform # updates that would normally be O(p^3) in O(p^2) time. If you # wish to use it, you will need to compile dchud.f, then link it as # a shared object. # # On MacOS X (64-bit Intel) with gfortran (from gcc 4.2): # # gfortran -O3 -arch x86_64 -c dchud.f # R --vanilla CMD SHLIB -o chud.so dchud.o # # On Solaris x86 with g77: # # g77 -O3 -c dchud.f # R --vanilla CMD SHLIB -o chud.so dchud.o # # If you plan to put chud.so somewhere other than the current directory, # you can change the path used by dyn.load. If it can't find chud.so, # it will fall back to an O(p^3) R implementation. # Calling the Samplers # -------------------- # # This file contains the functions covariance.match.sample() and # shrinking.rank.sample(), each implementing a sampler from the # technical report referenced above. # # Here is a simple example of calling each sampler on a standard # Gaussian in two dimensions: # # source("cass.r") # logf <- function(x) -(x[1]^2+x[2]^2)/2 # dlogf <- function(x) c(-x[1], -x[2]) # X1 <- covariance.match.sample(logf, dlogf, x0=c(0,0), N=200, sigma.c=5) # X2 <- shrinking.rank.sample(logf, dlogf, x0=c(0,0), N=200, sigma.c=5) # # The first argument to each sampler function is logf, which should # be a function that takes a single vector argument and returns the # log density at that point. # # The second argument also takes a vector argument and returns the # gradient of logf at that point. # # The third argument is the starting state of the chain; the dimension # of the target distribution is inferred from its length. # # sigma.c (and theta for covariance.match.sample()) are optional tuning # parameters described in the technical report. # # shrinking.rank.sample() can optionally take the arguments downscale # and min.dimension. # # downscale is the factor by which the crumb standard deviation in # the reduced subspace is scaled each time a proposal is rejected. # # min.dimension is the limit on how small the dimension of the subspace # from which crumbs are drawn is allowed to get. To simulate a sampler # equivalent to the "Non-Adaptive Crumbs" sampler described in the # technical report, set min.dimension=length(x0). # # Both sampler functions return an N by length(x0) matrix of samples. ### Covariance-Matching Slice Sampler ### covariance.match.sample <- function(logf, dlogf, x0, N, sigma.c=1, theta=1) { p <- length(x0) X <- array(NA,c(N,p)) X[1,] <- x0 for (samp in 2:N) { x0 <- X[samp-1,] # current target coordinates M <- logf(x0) # current estimate of the log density mode logy0 <- M - rexp(1) # \tilde y_0, the log density at the slice R <- diag(1/sigma.c, nrow=p) # Cholesky factor of the proposal covar. F <- diag(1/sigma.c, nrow=p) # Cholesky factor of the crumb covar. cbar.star <- numeric(p) # \bar c_k^*, the un-normalized crumb mean repeat { crumb <- x0 + backsolve(F, rnorm(p)) cbar.star <- cbar.star + crossprod(F, F %*% crumb) cbar <- backsolve(R, backsolve(R, cbar.star, transpose=TRUE)) x <- cbar + backsolve(R, rnorm(p)) # draw proposal logy <- logf(x) if (logy>=logy0) # accept the proposal? break G <- as.numeric(dlogf(x)) stopifnot(all(is.finite(G))) G.norm <- sqrt(sum(G^2)) g <- G / G.norm delta <- sqrt(sum((x-crumb)^2)) u <- x + delta * g l.u <- logf(u) kappa <- -2/delta^2 * ( l.u - logy - G.norm * delta ) l.xu <- 0.5 * G.norm^2/kappa + logy # peak of parabola M <- max(M, l.xu) sigma.sq <- 2/3 * (M-logy0)/kappa # estimate conditional variance Rg <- R %*% g alpha <- max(0, 1/sigma.sq - (1+theta)*sum(Rg * Rg)) F <- chud(sqrt(theta)*R, sqrt(alpha)*g) R <- chud(sqrt(1+theta)*R, sqrt(alpha)*g) } X[samp,] <- x } return(X) } ### Shrinking-Rank Slice Sampler ### shrinking.rank.sample <- function(logf, dlogf, x0, N, sigma.c=1, downscale=0.9, min.dimension=1) { P <- function(M,a) a - M %*% crossprod(M,a) p <- length(x0) X <- array(NA,c(N,p)) X[1,] <- x0 for (samp in 2:N) { x0 <- X[samp-1,] logy0 <- logf(x0) - rexp(1) J <- array(NA, c(p,0)) k <- 0 cbar.star <- numeric(p) # un-normalized mean of crumbs Lambda <- 0 # proposal precision repeat { k <- k + 1 sigma.ck <- sigma.c * downscale^(k-1) # crumb standard deviation crumb <- sigma.ck * rnorm(p) Lambda <- Lambda + sigma.ck^(-2) cbar.star <- cbar.star + crumb / sigma.ck^2 x <- x0 + P(J, cbar.star/Lambda + rnorm(p)/sqrt(Lambda)) if (logf(x)>=logy0) break if (ncol(J) < p-min.dimension) { gx <- as.numeric(dlogf(x)) gstar <- P(J, gx) if (sum(gstar*gx)/sqrt(sum(gstar^2))/sqrt(sum(gx^2)) > 0.5) J <- cbind(J, P(J,gstar/sqrt(sum(gstar^2)))) } } X[samp,] <- x } return(X) } ### R Interface for LINPACK DCHUD ### # Computes Q such that Q^T Q = R^T R + x x^T. chud <- function(R, x) { p <- nrow(R) stopifnot(is.matrix(R) && p==ncol(R)) stopifnot(is.numeric(x) && length(x)==p) L <- .Fortran('dchud', R, p, p, x, 0, 0, 0, 0, 0, numeric(p), numeric(p)) return(L[[1]]) } tryCatch(dyn.load('chud.so', local=TRUE), error=function(e) { warning("Could not load chud.so. Falling back to O(p^3) chud.") chud <<- function(R, x) chol(crossprod(R)+tcrossprod(x)) })