sns.R
sns <- function(x, fghEval, rnd=TRUE, gfit=NULL, mh.diag = FALSE, part = NULL
, numderiv = 0, numderiv.method = c("Richardson", "simple"), numderiv.args = list()
, ...)
{
fitGaussian <- function(x, f, ...)
{
ret <- f(x,...) # Evaluate the function at 'x'
Sigma <- solve(-ret$h)
mu <- x + Sigma %*% ret$g
return (list(mu=as.vector(mu), # Newton method solution
Sigma=Sigma, # Inverse Hessian or Covariance matrix
iSigma=-ret$h, # Inverse covariance or Hessian
f=ret$f, # function value
g=ret$g)) # gradient
}
numderiv <- as.integer(numderiv)
if (numderiv > 0) {
numderiv.method <- match.arg(numderiv.method)
if (numderiv == 1) { # we need numeric hessian
fghEval.int <- function(x, ...) {
fg <- fghEval(x, ...)
h <- hessian(func = function(x, ...) fghEval(x, ...)$f, x = x, ..., method = numderiv.method, method.args = numderiv.args)
return (list(f = fg$f, g = fg$g, h = h))
}
} else { # we need numeric gradient and hessian
fghEval.int <- function(x, ...) {
f <- fghEval(x, ...)
g <- grad(func = fghEval, x = x, ..., method = numderiv.method, method.args = numderiv.args)
h <- hessian(func = fghEval, x = x, ..., method = numderiv.method, method.args = numderiv.args)
return (list(f = f, g = g, h = h))
}
}
} else {
fghEval.int <- fghEval
}
if (!is.null(part)) { # code for 'state space partitioning', recursive call
fghEval.part <- function(xsub, xfull, subset, ...) {
xfull[subset] <- xsub
ret <- fghEval.int(xfull, ...)
return (list(f = ret$f, g = ret$g[subset], h = ret$h[subset, subset]))
}
npart <- length(part)
accept <- logical(npart)
if (mh.diag) {
diag <- array(NA, dim = c(4, npart))
dimnames(diag) <- list(c("log.p", "log.p.prop", "log.q", "log.q.prop"), 1:npart)
}
for (n in 1:npart) {
subset <- part[[n]]
rettmp <- sns(x = x[subset], fghEval = fghEval.part, rnd = rnd, gfit = NULL, part = NULL, mh.diag = mh.diag, xfull = x, subset = subset, ...)
x[subset] <- rettmp
accept[n] <- attr(rettmp, "accept")
if (mh.diag) diag[, n] <- attr(rettmp, "mh.diag")
}
attr(x, "lp") <- attr(rettmp, "lp")
attr(x, "accept") <- accept
attr(x, "gfit") <- fitGaussian(x, fghEval.int, ...)
if (mh.diag) attr(x, "mh.diag") <- mh.diag
return (x)
}
# rnd: if FALSE, perform Newton's optimization (non-stochastic)
# Fit Gaussian at x
if (is.null(gfit)) gfit <- fitGaussian(x = x, f = fghEval.int, ...)
mu <- gfit$mu
Sigma <- gfit$Sigma # Covariance
iSigma <- gfit$iSigma # Inverse covariance
K <- length(x)
if (rnd) {
# Draw sample from proposal distribution (Gaussian fit at x)
x.prop <- as.vector(rmvnorm(n=1, mean=mu, sigma=Sigma))
} else {
# Run (non-stochastic) Newton optimization
rho <- 0.5; c <- 0.5;
alphak <- 1;
d <- mu - x; # use newton's direction as step
search_x <- as.vector(mu);
fk <- gfit$f; # Values at the current point
gk <- gfit$g;
fk1 <- fghEval.int(search_x, ...)$f; # Function value at searching point
ls_iter <- 1;
# Linesearch by backtracking from full Newton step
while (fk1 < fk + c*alphak*(t(gk)%*%d) && ls_iter < 20) {
alphak <- alphak*rho; # if so, then go half way
search_x <- x + alphak*d;
fk1 <- fghEval.int(search_x, ...)$f;
ls_iter <- ls_iter + 1;
}
x.prop <- as.vector(search_x);
}
log.q.prop <- dmvnorm(as.vector(x.prop), mu, Sigma, log=TRUE)
# fit Gaussian at x.prop
gfit.prop <- fitGaussian(x=x.prop,f=fghEval.int,...)
mu.prop <- gfit.prop$mu
Sigma.prop <- gfit.prop$Sigma
iSigma.prop <- gfit.prop$iSigma
# create MH acceptance ratio
log.q <- dmvnorm(as.vector(x), mu.prop, Sigma.prop, log=TRUE)
log.p <- gfit$f
log.p.prop <- gfit.prop$f
log.ratio <- (log.p.prop-log.p) + (log.q-log.q.prop)
ratio <- min(1, exp(log.ratio))
# perform acceptance test
if (!rnd || ratio==1 || runif(1)<ratio) {
gfit <- gfit.prop
x <- x.prop;
attr(x,"accept") <- TRUE
attr(x,"lp") <- log.p.prop
} else {
attr(x,"accept") <- FALSE
attr(x,"lp") <- log.p
}
attr(x,"gfit") <- gfit
if (mh.diag) {
diag <- c(log.p, log.p.prop, log.q, log.q.prop)
names(diag) <- c("log.p", "log.p.prop", "log.q", "log.q.prop")
attr(x, "mh.diag") <- diag
}
return (x)
}
sns.run <- function(init, fghEval, niter = 100, nnr = min(10, round(niter/4))
, mh.diag = FALSE, part = NULL, print.level = 0
, report.progress = ceiling(niter/10)
, numderiv = 0, numderiv.method = c("Richardson", "simple"), numderiv.args = list()
, ...)
{
fghEval.int <- sns.fghEval.numaug(fghEval, numderiv, numderiv.method, numderiv.args)
# checking arguments
stopifnot(niter >= 1)
if (report.progress <= 0) {
warning("invalid value specific for 'report.progress', using default.")
report.progress <- ceiling(niter/10)
}
if (missing(init)) stop("starting point for MCMC chain must be specified")
npart <- max(length(part),1)
# initialization
K <- length(init)
x <- init
accept <- matrix(logical(niter * npart), ncol = npart)
lp <- double(niter)
chain <- matrix( , nrow = niter, ncol = K)
if (mh.diag) {
diagnostic <- array(NA, dim = c(niter, 4, npart))
dimnames(diagnostic) <- list(1:niter, c("log.p", "log.p.prop", "log.q", "log.q.prop"), 1:npart)
if (niter > nnr && npart == 1) f.reldev <- rep(NA, niter - nnr)
}
# performing nr/mcmc iterations
for (i in 1:niter) {
x <- sns(x, fghEval.int, rnd = i > nnr, gfit = attr(x, "gfit"), mh.diag = mh.diag, part = part
, numderiv = 0#, numderiv.method = numderiv.method, numderiv.args = numderiv.args
, ...)
accept[i, ] <- attr(x, "accept")
lp[i] <- attr(x, "lp")
chain[i, ] <- x
if (mh.diag) {
diagnostic[i, , ] <- attr(x, "mh.diag")
if (npart == 1) {
if (i == nnr) {
f.ref <- attr(x, "gfit")$f
x.ref <- attr(x, "gfit")$mu
} else if (i > nnr) {
if (attr(x, "gfit")$f < f.ref) f.reldev[i - nnr] <- (attr(x, "gfit")$f - f.ref + 0.5 * t(x - x.ref) %*% attr(x, "gfit")$iSigma %*% (x - x.ref)) / (f.ref - attr(x, "gfit")$f)
}
}
}
if (print.level && (i %% report.progress == 0))
cat(paste0("finished iter ", i, " of ", niter, "\n"))
}
# assembling output
attr(chain, "init") <- init
attr(chain, "lp.init") <- fghEval.int(init, ...)$f
attr(chain, "accept") <- accept
attr(chain, "lp") <- lp
attr(chain, "nnr") <- nnr
attr(chain, "part") <- part
if (mh.diag) {
attr(chain, "mh.diag") <- drop(diagnostic)
if (niter > nnr && npart == 1) attr(chain, "reldev") <- f.reldev
}
class(chain) <- "sns"
return (chain)
}