sns.run: Drawing multiple samples using Stochastic Newton Sampler

In sns: Stochastic Newton Sampler (SNS)

Drawing multiple samples using Stochastic Newton Sampler

Description

This is a wrapper around sns, allowing one to draw multiple samples from a distribution while collecting diagnostic information.

Usage

sns.run(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()
  , ...)

Arguments

init	Initial value for the MCMC chain.
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.
niter	Number of iterations to perform (in ‘nr’ and ‘mcmc’ mode combined).
nnr	Number of initial iterations to spend in ‘nr’ mode.
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.
print.level	If greater than 0, print sampling progress report.
report.progress	Number of sampling iterations to wait before printing progress reports.
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 parameters to be passed to fghEval.

Value

sns.run returns an object of class sns with elements:

samplesMat	A matrix object with nsample rows and K cols.
acceptance	Metropolis proposal percentage acceptance.
burn.iters	Number of burn-in ierations.
sample.time	Time in seconds spent in sampling.
burnin.time	Time in seconds spent in burn-in.

Note

1. sns.run cannot be used if SNS is being run as part of a Gibbs cycle, such that the conditional distribution being sampled by SNS changes from one iteration to next. In such cases, sns must be used instead, inside an explicit Gibbs-cycle for loop.

2. 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

See Also

sns, summary.sns, plot.sns, predict.sns

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.init <- rep(0.0, K)

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

# sampling of likelihood
beta.smp <- sns.run(init = beta.init
  , fghEval = loglike.poisson, niter = 1000
  , nnr = 20, X = X, y = y)
smp.summ <- summary(beta.smp)

# compare mean of samples against ML estimate (from glm)
cbind(beta.glm, smp.summ$smp$mean)

# 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.smp <- sns.run(init = beta.init
  , fghEval = loglike.poisson.fonly, niter = 1000, nnr = 20
  , X = X, y = y, numderiv = 2)
smp.summ <- summary(beta.smp)
cbind(beta.glm, smp.summ$smp$mean)


## End(Not run)

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