######## BENCHMARKING FUNCTIONS ######## library("sns") library("MfUSampler") library("SamplerCompare") library("RegressionFactory") # collection of functions used for in performance benchmarking section of # manuscript benchmark.gendata <- function(N = 1000, K = 10, family = c("binomial", "poisson", "exponential"), off.diag = 0) { family <- match.arg(family) my.sigma <- matrix(off.diag, nrow = K - 1, ncol = K - 1) diag(my.sigma) <- 1 X <- rmvnorm(N, sigma = my.sigma) X <- cbind(1, X) beta <- runif(K, min = -0.5, max = +0.5) Xbeta <- X %*% beta if (family == "binomial") { y <- 1 * (runif(N) < 1/(1 + exp(-Xbeta))) } else if (family == "poisson") { y <- rpois(n = N, lambda = exp(Xbeta)) } else if (family == "exponential") { y <- rexp(n = N, rate = exp(Xbeta)) } ret <- list(X = X, y = y, beta = beta, family = family, N = N, K = K) return(ret) } benchmark.sample <- function(X, y, family = c("binomial", "poisson", "exponential"), n.mcmc = 100, n.init = 10, sampler = c("SNS", "slice", "ARS", "AM")) { family <- match.arg(family) sampler <- match.arg(sampler) N <- nrow(X) K <- ncol(X) n.iter <- n.init + n.mcmc if (family == "binomial") { my.fbase <- fbase1.binomial.logit } else if (family == "poisson") { my.fbase <- fbase1.poisson.log } else if (family == "exponential") { my.fbase <- fbase1.exponential.log } loglike.sns <- function(beta, X, y) { regfac.expand.1par(beta, X, y, fbase1 = my.fbase, fgh = 2) } loglike.slice <- function(beta, X, y) { regfac.expand.1par(beta, X, y, fbase1 = my.fbase, fgh = 0) } loglike.ars <- function(beta, X, y, grad) { if (grad) { return(regfac.expand.1par(beta, X, y, fbase1 = my.fbase, fgh = 1)$g) } else { return(regfac.expand.1par(beta, X, y, fbase1 = my.fbase, fgh = 0)) } } loglike.am <- make.dist(ndim = K, name = "glm", log.density = function(u) { regfac.expand.1par(u, X, y, my.fbase, fgh = 0) }) t <- NA beta.init <- rep(0, K) for (i in 1:10) beta.init <- sns(x = beta.init, fghEval = loglike.sns, rnd = F, X = X, y = y) if (sampler == "SNS") { t <- proc.time()[3] ret <- sns.run(init = beta.init, fghEval = loglike.sns, niter = n.mcmc, nnr = 0, X = X, y = y) t <- proc.time()[3] - t } else if (sampler == "slice") { lp <- rep(NA, n.iter) ret <- array(NA, dim = c(n.mcmc, K)) beta <- beta.init t <- proc.time()[3] for (n in 1:n.mcmc) { beta <- MfU.Sample(x = beta, f = loglike.slice, uni.sampler = "slice", X = X, y = y) ret[n, ] <- beta # lp[n] <- loglike.slice(beta, X, y) } t <- proc.time()[3] - t # attr(ret, 'lp') <- lp } else if (sampler == "ARS") { lp <- rep(NA, n.iter) ret <- array(NA, dim = c(n.mcmc, K)) beta <- beta.init t <- proc.time()[3] for (n in 1:n.mcmc) { beta <- MfU.Sample(x = beta, f = loglike.ars, uni.sampler = "ars", X = X, y = y) ret[n, ] <- beta # lp[n] <- loglike.slice(beta, X, y) } t <- proc.time()[3] - t # attr(ret, 'lp') <- lp } else if (sampler == "AM") { t <- proc.time()[3] am.out <- adaptive.metropolis.sample(target.dist = loglike.am, x0 = beta.init, sample.size = n.iter, burn.in = n.init/n.iter) ret <- am.out$X t <- proc.time()[3] - t } attr(ret, "t") <- t return(ret) } benchmark.driver <- function(seed = 0, op = mean, nsmp.num = 10000, config.nk = data.frame(N = c(1000), K = c(10)), family.vec = c("binomial", "poisson", "exponential"), samplers = c("SNS", "slice", "ARS", "AM"), off.diag = 0) { set.seed(seed) config <- merge(merge(merge(config.nk, data.frame(family = family.vec)), data.frame(sampler = samplers)), data.frame(off.diag = off.diag)) config$nsmp <- nsmp.num/config$K data.config <- unique(config[, c("N", "K", "family", "off.diag")]) data.list <- lapply(1:nrow(data.config), function(i) { benchmark.gendata(N = data.config$N[i], K = data.config$K[i], family = as.character(data.config$family[i]), off.diag = data.config$off.diag[i]) }) nconfig <- nrow(config) config$t <- NA config$ess <- NA for (i in 1:nconfig) { data.index <- which(data.config$N == config$N[i] & data.config$K == config$K[i] & data.config$family == config$family[i] & data.config$off.diag == config$off.diag[i]) dat <- data.list[[data.index]] mcmc <- benchmark.sample(X = dat$X, y = dat$y, family = as.character(config$family[i]), n.mcmc = config$nsmp[i], n.init = round(config$nsmp[i]/10), sampler = as.character(config$sampler[i])) config$t[i] <- attr(mcmc, "t") config$ess[i] <- op(ess(mcmc, method = "ise")) cat("finished run", i, "of", nconfig, "\n") } config$iss <- config$ess/config$t return(config) } benchmark.plot <- function(rslt) { unique.family.data <- unique(rslt[, c("family", "N", "K")]) unique.family <- unique(rslt[, c("family")]) if (length(unique.family.data) > length(unique.family)) stop("invalid result set for plotting") nfamily <- length(unique.family) pch.vec <- 1:nfamily ylim <- range(rslt$iss) for (i in 1:nfamily) { rslt.extract <- rslt[which(rslt$family == unique.family[i]), ] if (i == 1) { plot(rslt.extract$iss, xlab = "sampler", ylab = "independent samples per second", type = "b", ylim = ylim, xaxt = "n", pch = pch.vec[i]) axis(side = 1, at = 1:nrow(rslt.extract), labels = rslt.extract$sampler) } else { lines(rslt.extract$iss, pch = pch.vec[i], type = "b") } } legend("topright", legend = unique.family, pch = pch.vec) } plot.bench.N <- function(df, samplers = c("SNS", "slice")) { nsamplers <- length(samplers) ycols <- paste("iss", samplers, sep = ".") ylim <- range(df[, ycols]) pch.vec <- 1:nsamplers plot(df$N, df[, ycols[1]], type = "b", ylim = ylim, log = "x", xlab = "N", ylab = "independent samples / sec", pch = pch.vec[1]) if (nsamplers > 1) { for (i in 2:nsamplers) { lines(df$N, df[, ycols[i]], type = "b", pch = pch.vec[i]) } } legend("topleft", legend = samplers, pch = pch.vec, lty = rep(1, nsamplers)) } plot.bench.corr <- function(df, samplers = c("SNS", "slice")) { nsamplers <- length(samplers) ycols <- paste("iss", samplers, sep = ".") ylim <- range(df[, ycols]) pch.vec <- 1:nsamplers plot(df$off.diag, df[, ycols[1]], type = "b", ylim = ylim #, log = 'x' , xlab = "correlation coefficient", ylab = "independent samples / sec", pch = pch.vec[1]) if (nsamplers > 1) { for (i in 2:nsamplers) { lines(df$off.diag, df[, ycols[i]], type = "b", pch = pch.vec[i]) } } legend("right", legend = samplers, pch = pch.vec, lty = rep(1, nsamplers)) } ######## SECTION 4 USING SNS ########## ################# library("sns") library("mvtnorm") my.seed <- 0 ################# logdensity.mvg <- function(x, mu, isigsq) { f <- dmvnorm(x = as.numeric(x), mean = mu, sigma = solve(isigsq), log = TRUE) g <- -isigsq %*% (x - mu) h <- -isigsq return(list(f = f, g = g, h = h)) } ################# set.seed(my.seed) K <- 3 mu <- runif(K, min = -0.5, max = +0.5) isigsq <- matrix(runif(K * K, min = 0.1, max = 0.2), ncol = K) isigsq <- 0.5 * (isigsq + t(isigsq)) diag(isigsq) <- rep(0.5, K) x.init <- rep(0, K) x.smp <- sns.run(x.init, logdensity.mvg, niter = 500, mh.diag = TRUE, mu = mu, isigsq = isigsq) summary(x.smp) ################# library("RegressionFactory") loglike.poisson <- function(beta, X, y) { regfac.expand.1par(beta, X = X, y = y, fbase1 = fbase1.poisson.log) } ################# set.seed(my.seed) K <- 5 N <- 1000 X <- matrix(runif(N * K, -0.5, +0.5), ncol = K) beta <- runif(K, -0.5, +0.5) y <- rpois(N, exp(X %*% beta)) ################# beta.init <- rep(0, K) beta.glm <- glm(y ~ X - 1, family = "poisson", start = beta.init)$coefficients ################# beta.sns <- sns.run(beta.init, fghEval = loglike.poisson, niter = 20, nnr = 20, X = X, y = y) beta.nr <- beta.sns[20, ] cbind(beta.glm, beta.nr) ################# beta.smp <- sns.run(beta.init, loglike.poisson, niter = 200, nnr = 20, mh.diag = TRUE, X = X, y = y) plot(beta.smp, select = 1) summary(beta.smp) ################# beta.smp <- sns.run(beta.init, loglike.poisson, niter = 1000, nnr = 20, mh.diag = TRUE, X = X, y = y) predmean.poisson <- function(beta, Xnew) exp(Xnew %*% beta) ymean.new <- predict(beta.smp, predmean.poisson, nburnin = 100, Xnew = X) predsmp.poisson <- function(beta, Xnew) rpois(nrow(Xnew), exp(Xnew %*% beta)) ysmp.new <- predict(beta.smp, predsmp.poisson, nburnin = 100, Xnew = X) summary(ymean.new) summary(ysmp.new) ################# set.seed(my.seed) K <- 100 X <- matrix(runif(N * K, -0.5, +0.5), ncol = K) beta <- runif(K, -0.5, +0.5) y <- rpois(N, exp(X %*% beta)) beta.init <- glm(y ~ X - 1, family = "poisson")$coefficients beta.smp <- sns.run(beta.init, loglike.poisson, niter = 100, nnr = 10, mh.diag = TRUE, X = X, y = y) summary(beta.smp) ################# beta.smp.part <- sns.run(beta.init, loglike.poisson, niter = 100, nnr = 10, mh.diag = TRUE, part = sns.make.part(K, 10), X = X, y = y) summary(beta.smp.part) ################# par(mfrow = c(1, 2)) plot(beta.smp, select = 1) plot(beta.smp.part, select = 1) ################# loglike.linreg.het <- function(coeff, X, Z, y) { K1 <- ncol(X) K2 <- ncol(Z) beta <- coeff[1:K1] gamma <- coeff[K1 + 1:K2] mu <- X %*% beta sigma <- sqrt(exp(Z %*% gamma)) f <- sum(dnorm(y, mu, sigma, log = TRUE)) return(f) } ################# set.seed(my.seed) K1 <- 5 K2 <- 5 N <- 1000 X <- matrix(runif(N * K1, -0.5, +0.5), ncol = K1) Z <- matrix(runif(N * K2, -0.5, +0.5), ncol = K2) beta <- runif(K1, -0.5, +0.5) gamma <- runif(K1, -0.5, +0.5) mu <- X %*% beta var <- exp(Z %*% gamma) y <- rnorm(N, X %*% beta, sd = sqrt(var)) ################# coeff.init <- rep(0, K1 + K2) check.logdensity <- sns.check.logdensity(coeff.init, loglike.linreg.het, X = X, Z = Z, y = y, dx = 1, nevals = 10, blocks = list(1:(K1 + K2), 1:K1, K1 + 1:K2)) check.logdensity ################# loglike.linreg.het.beta <- function(beta, gamma, X, Z, y) { K1 <- length(beta) ret <- regfac.expand.2par(c(beta, gamma), X, Z, y, fbase2 = fbase2.gaussian.identity.log) return(list(f = ret$f, g = ret$g[1:K1], h = ret$h[1:K1, 1:K1])) } loglike.linreg.het.gamma <- function(gamma, beta, X, Z, y) { K1 <- length(beta) K2 <- length(gamma) ret <- regfac.expand.2par(c(beta, gamma), X, Z, y, fbase2 = fbase2.gaussian.identity.log) return(list(f = ret$f, g = ret$g[K1 + 1:K2], h = ret$h[K1 + 1:K2, K1 + 1:K2])) } ################# nsmp <- 100 beta.iter <- rep(0, K1) gamma.iter <- rep(0, K2) beta.smp <- array(NA, dim = c(nsmp, K1)) gamma.smp <- array(NA, dim = c(nsmp, K1)) for (n in 1:nsmp) { beta.iter <- sns(beta.iter, loglike.linreg.het.beta, gamma = gamma.iter, X = X, Z = Z, y = y, rnd = nsmp > 10) gamma.iter <- sns(gamma.iter, loglike.linreg.het.gamma, beta = beta.iter, X = X, Z = Z, y = y, rnd = nsmp > 10) beta.smp[n, ] <- beta.iter gamma.smp[n, ] <- gamma.iter } beta.est <- colMeans(beta.smp[(nsmp/2 + 1):nsmp, ]) gamma.est <- colMeans(gamma.smp[(nsmp/2 + 1):nsmp, ]) cbind(beta, beta.est) cbind(gamma, gamma.est) ################# library("MfUSampler") loglike.linreg.het.gamma.fonly <- function(gamma, beta, X, Z, y) { return(regfac.expand.2par(c(beta, gamma), X, Z, y, fbase2 = fbase2.gaussian.identity.log, fgh = 0)) } beta.iter <- rep(0, K1) gamma.iter <- rep(0, K2) for (n in 1:nsmp) { beta.iter <- sns(beta.iter, loglike.linreg.het.beta, gamma = gamma.iter, X = X, Z = Z, y = y, rnd = nsmp > 10) gamma.iter <- MfU.Sample(gamma.iter, loglike.linreg.het.gamma.fonly, beta = beta.iter, X = X, Z = Z, y = y, uni.sampler = "slice") beta.smp[n, ] <- beta.iter gamma.smp[n, ] <- gamma.iter } beta.est.mix <- colMeans(beta.smp[(nsmp/2 + 1):nsmp, ]) gamma.est.mix <- colMeans(gamma.smp[(nsmp/2 + 1):nsmp, ]) cbind(beta, beta.est.mix) cbind(gamma, gamma.est.mix) ################# loglike.linreg.het.gamma.numaug <- sns.fghEval.numaug(loglike.linreg.het.gamma.fonly, numderiv = 2) beta.iter <- rep(0, K1) gamma.iter <- rep(0, K2) for (n in 1:nsmp) { beta.iter <- sns(beta.iter, loglike.linreg.het.beta, gamma = gamma.iter, X = X, Z = Z, y = y, rnd = nsmp > 10) gamma.iter <- sns(gamma.iter, loglike.linreg.het.gamma, beta = beta.iter, X = X, Z = Z, y = y, rnd = nsmp > 10) beta.smp[n, ] <- beta.iter gamma.smp[n, ] <- gamma.iter } beta.est.num <- colMeans(beta.smp[(nsmp/2 + 1):nsmp, ]) gamma.est.num <- colMeans(gamma.smp[(nsmp/2 + 1):nsmp, ]) cbind(beta, beta.est.num) cbind(gamma, gamma.est.num) ################# sessionInfo() ######## SECTION 5 BENCHMARKING ####### # The results were generated under R version 3.3.1 with sns 1.1.2, # MfUSampler 1.0.2, RegressionFactory 0.7.2, and mvtnorm 1.0-5. # NOTE: CPU times are bound to vary from run to run, even on the same machine, # resulting in different ISS (independent samples per second) numbers (ESS # numbers should remain the same as long as same random seed is used) # on a Windows laptop with Intel Core i7 and 16GB installed RAM, the total # execution time for this script was about 1.5 hours t <- proc.time()[3] # Table 1: performance comparison of SNS and several other samplers for N=1000, # K=10, off.diag=0.0 using mutiple distributions bench <- benchmark.driver(seed = 0, nsmp.num = 1e+05) save(bench, file = "bench.df") # calculating speedup numbers binomial which.binomial.other <- which(bench$family == "binomial" & bench$sampler != "SNS") which.binomial.sns <- which(bench$family == "binomial" & bench$sampler == "SNS") speedup.binomial <- bench$iss[which.binomial.sns]/max(bench$iss[which.binomial.other]) cat("SNS speedup - binomial:", speedup.binomial, "\n") # poisson which.poisson.other <- which(bench$family == "poisson" & bench$sampler != "SNS") which.poisson.sns <- which(bench$family == "poisson" & bench$sampler == "SNS") speedup.poisson <- bench$iss[which.poisson.sns]/max(bench$iss[which.poisson.other]) cat("SNS speedup - poisson:", speedup.poisson, "\n") # exponential which.exponential.other <- which(bench$family == "exponential" & bench$sampler != "SNS") which.exponential.sns <- which(bench$family == "exponential" & bench$sampler == "SNS") speedup.exponential <- bench$iss[which.exponential.sns]/max(bench$iss[which.exponential.other]) cat("SNS speedup - exponential:", speedup.exponential, "\n") # Figure 3: impact of correlation structure in X (and hence posterior # distribution for beta's) on relative performance of SNS, slicer, and adaptive # Meteropolis samplers.corr <- c("SNS", "slice", "AM") bench.corr <- benchmark.driver(seed = 0, nsmp.num = 1e+05, samplers = samplers.corr, off.diag = c(0, 0.25, 0.5, 0.75, 0.9)) save(bench.corr, file = "bench.corr.df") bench.corr.wide.binomial <- reshape(bench.corr[which(bench.corr$family == "binomial"), ], v.names = c("t", "ess", "iss"), idvar = c("off.diag", "family"), timevar = "sampler", direction = "wide") bench.corr.wide.poisson <- reshape(bench.corr[which(bench.corr$family == "poisson"), ], v.names = c("t", "ess", "iss"), idvar = c("off.diag", "family"), timevar = "sampler", direction = "wide") bench.corr.wide.exponential <- reshape(bench.corr[which(bench.corr$family == "exponential"), ], v.names = c("t", "ess", "iss"), idvar = c("off.diag", "family"), timevar = "sampler", direction = "wide") pdf(file = "fig_bench_corr_binomial.pdf") plot.bench.corr(bench.corr.wide.binomial, samplers = samplers.corr) title(main = "binomial") dev.off() pdf(file = "fig_bench_corr_poisson.pdf") plot.bench.corr(bench.corr.wide.poisson, samplers = samplers.corr) title(main = "poisson") dev.off() pdf(file = "fig_bench_corr_exponential.pdf") plot.bench.corr(bench.corr.wide.exponential, samplers = samplers.corr) title(main = "exponential") dev.off() # impact of data size (N) on relative performance of samplers samplers.N <- c("SNS", "slice", "AM") bench.N <- benchmark.driver(seed = 0, nsmp.num = 1e+05, config.nk = data.frame(N = c(50, 100, 250, 1000, 5000), K = 10), samplers = samplers.N, off.diag = 0) save(bench.N, file = "bench.N.df") bench.N.wide.binomial <- reshape(bench.N[which(bench.N$family == "binomial"), ], v.names = c("t", "ess", "iss"), idvar = c("N", "K", "family"), timevar = "sampler", direction = "wide") bench.N.wide.poisson <- reshape(bench.N[which(bench.N$family == "poisson"), ], v.names = c("t", "ess", "iss"), idvar = c("N", "K", "family"), timevar = "sampler", direction = "wide") bench.N.wide.exponential <- reshape(bench.N[which(bench.N$family == "exponential"), ], v.names = c("t", "ess", "iss"), idvar = c("N", "K", "family"), timevar = "sampler", direction = "wide") pdf(file = "fig_bench_N_binomial.pdf") plot.bench.N(bench.N.wide.binomial, samplers = samplers.N) title(main = "binomial") dev.off() pdf(file = "fig_bench_N_poisson.pdf") plot.bench.N(bench.N.wide.poisson, samplers = samplers.N) title(main = "poisson") dev.off() pdf(file = "fig_bench_N_exponential.pdf") plot.bench.N(bench.N.wide.exponential, samplers = samplers.N) title(main = "exponential") dev.off() t <- proc.time()[3] - t cat("total execution time:", t/60, "min\n")