sns: Stochastic Newton Sampler (SNS)

In sns: Stochastic Newton Sampler (SNS)

Stochastic Newton Sampler (SNS)

Description

SNS is a Metropolis-Hastings MCMC sampler with a multivariate Gaussian proposal function resulting from a local, second-order Taylor series expansion of log-density. The mean of the Gaussian proposal is identical to the full Newton-Raphson step from the current point. During burn-in, Newton-Raphson optimization can be performed to get close to the mode of the pdf which is unique due to convexity, resulting in faster convergence. For high dimensional densities, state space partitioning can be used to improve mixing. Support for numerical differentiation is provided using numDeriv package. sns is the low-level function for drawing one sample from the distribution. For drawing multiple samples from a (fixed) distribution, consider using sns.run.

Usage

sns(x, fghEval, rnd = TRUE, gfit = NULL, mh.diag = FALSE
  , part = NULL, numderiv = 0
  , numderiv.method = c("Richardson", "simple")
  , numderiv.args = list(), ...)

Arguments

x	Current state vector.
fghEval	Log-density to be sampled from. A valid log-density can have one of 3 forms: 1) return log-density, but no gradient or Hessian, 2) return a list of f and g for log-density and its gradient vector, respectively, 3) return a list of f, g, and h for log-density, gradient vector, and Hessian matrix. Missing derivatives are computed numerically.
rnd	Runs 1 iteration of Newton-Raphson optimization method (non-stochastic or 'nr' mode) when FALSE. Runs Metropolis-Hastings (stochastic or 'mcmc' mode) for drawing a sample when TRUE.
gfit	Gaussian fit at point init. If NULL then sns will compute a Gaussian fit at x.
mh.diag	Boolean flag, indicating whether detailed MH diagnostics such as components of acceptance test must be returned or not.
part	List describing partitioning of state space into subsets. Each element of the list must be an integer vector containing a set of indexes (between 1 and length(x) or length(init)) indicating which subset of all dimensions to jointly sample. These integer vectors must be mutually exclusive and collectively exhaustive, i.e. cover the entire state space and have no duplicates, in order for the partitioning to represent a valid Gibbs sampling approach. See sns.make.part and sns.check.part.
numderiv	Integer with value from the set 0,1,2. If 0, no numerical differentiation is performed, and thus fghEval is expected to supply f, g and h. If 1, we expect fghEval to provide f amd g, and Hessian will be calculated numerically. If 2, fghEval only returns log-density, and numerical differentiation is needed to calculate gradient and Hessian.
numderiv.method	Method used for numeric differentiation. This is passed to the grad and hessian functions in numDeriv package. See the package documentation for details.
numderiv.args	Arguments to the numeric differentiation method chosen in numderiv.method, passed to grad and hessian functions in numDeriv. See package documentation for details.
...	Other arguments to be passed to fghEval.

Value

sns returns the sample drawn as a vector, with attributes:

accept	A boolean indicating whether the proposed point was accepted.
ll	Value of the log-density at the sampled point.
gfit	List containing Gaussian fit to pdf at the sampled point.

Note

1. Since SNS makes local Gaussian approximations to the density with the covariance matrix of the Gaussian proposal being the log-density Hessian, there is a strict requirement for the log-density to be concave.

2. Proving log-concavity for arbitrary probability distributions is non-trvial. However, distributions generated by replacing parameters of a concave distribution with linear expressions are known to be log-concave. This negative-definiteness invariance as well as expressions for full gradient and Hessian in terms of derivatives of low-dimensional base distributions are discussed in the vignette. The GLM expansion framework is available in the R package RegressionFactory.

3. See package vignette for more details on SNS theory, software, examples, and performance.

Author(s)

Alireza S. Mahani, Asad Hasan, Marshall Jiang, Mansour T.A. Sharabiani

References

Mahani A.S., Hasan A., Jiang M. & Sharabiani M.T.A. (2016). Stochastic Newton Sampler: The R Package sns. Journal of Statistical Software, Code Snippets, 74(2), 1-33. doi:10.18637/jss.v074.c02

Hastings, W. K. (1970). Monte Carlo sampling methods using Markov chains and their applications. Biometrika, 57(1), 97-109.

Qi, Y., & Minka, T. P. (2002). Hessian-based markov chain monte-carlo algorithms. 1st Cape Cod Workshop on Monte Carlo Methods.

See Also

sns.run, sns.fghEval.numaug

Examples

## Not run: 

# using RegressionFactory for generating log-likelihood and its derivatives
library(RegressionFactory)

loglike.poisson <- function(beta, X, y) {
  regfac.expand.1par(beta, X = X, y = y,
                     fbase1 = fbase1.poisson.log)
}

# simulating data
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.0, K)

# glm estimate, for reference
beta.glm <- glm(y ~ X - 1, family = "poisson",
                start = beta.init)$coefficients

# running SNS in non-stochastic mode
# this should produce results very close to glm
beta.sns <- beta.init
for (i in 1:20)
  beta.sns <- sns(beta.sns, loglike.poisson, X = X, y = y, rnd = F)

# comparison
all.equal(as.numeric(beta.glm), as.numeric(beta.sns))

# trying numerical differentiation
loglike.poisson.fonly <- function(beta, X, y) {
  regfac.expand.1par(beta, X = X, y = y, fgh = 0,
                     fbase1 = fbase1.poisson.log)
}

beta.sns.numderiv <- beta.init
for (i in 1:20)
  beta.sns.numderiv <- sns(beta.sns.numderiv, loglike.poisson.fonly
                  , X = X, y = y, rnd = F, numderiv = 2)
all.equal(as.numeric(beta.glm), as.numeric(beta.sns.numderiv))

# add numerical derivatives to fghEval outside sns
loglike.poisson.numaug <- sns.fghEval.numaug(loglike.poisson.fonly
  , numderiv = 2)

beta.sns.numaug <- beta.init
for (i in 1:20)
  # set numderiv to 0 to avoid repeating 
  # numerical augmentation inside sns
  beta.sns.numaug <- sns(beta.sns.numaug, loglike.poisson.numaug
                           , X = X, y = y, rnd = F, numderiv = 0)
all.equal(as.numeric(beta.glm), as.numeric(beta.sns.numaug))


## End(Not run)

sns documentation built on Nov. 2, 2022, 5:15 p.m.