asymptotic_var: Asymptotic Variance of IS-type Estimators

In bssm: Bayesian Inference of Non-Linear and Non-Gaussian State Space Models

Asymptotic Variance of IS-type Estimators

Description

The asymptotic variance MCMCSE^2 is based on Corollary 1 of Vihola et al. (2020) from weighted samples from IS-MCMC. The default method is based on the integrated autocorrelation time (IACT) by Sokal (1997) which seem to work well for reasonable problems, but it is also possible to use the Geyer's method as implemented in ess_mean of the posterior package.

Usage

asymptotic_var(x, w, method = "sokal")

Arguments

x	A numeric vector of samples.
w	A numeric vector of weights. If missing, set to 1 (i.e. no weighting is assumed).
method	Method for computing IACT. Default is "sokal", other option "geyer".

Value

A single numeric value of asymptotic variance estimate.

References

Vihola M, Helske J, Franks J. (2020). Importance sampling type estimators based on approximate marginal Markov chain Monte Carlo. Scand J Statist. 1-38. https://doi.org/10.1111/sjos.12492

Sokal A. (1997). Monte Carlo Methods in Statistical Mechanics: Foundations and New Algorithms. In: DeWitt-Morette C, Cartier P, Folacci A (eds) Functional Integration. NATO ASI Series (Series B: Physics), vol 361. Springer, Boston, MA. https://doi.org/10.1007/978-1-4899-0319-8_6

Gelman, A, Carlin J B, Stern H S, Dunson, D B, Vehtari A, Rubin D B. (2013). Bayesian Data Analysis, Third Edition. Chapman and Hall/CRC.

Vehtari A, Gelman A, Simpson D, Carpenter B, Bürkner P-C. (2021). Rank-normalization, folding, and localization: An improved Rhat for assessing convergence of MCMC. Bayesian analysis, 16(2):667-718. https://doi.org/10.1214/20-BA1221

Examples

set.seed(1)
n <- 1e4 
x <- numeric(n)
phi <- 0.7
for(t in 2:n) x[t] <- phi * x[t-1] + rnorm(1)
w <- rexp(n, 0.5 * exp(0.001 * x^2))
# different methods:
asymptotic_var(x, w, method = "sokal")
asymptotic_var(x, w, method = "geyer")

data("negbin_model")
# can be obtained directly with summary method
d <- suppressWarnings(as_draws(negbin_model))
sqrt(asymptotic_var(d$sd_level, d$weight))

bssm documentation built on Nov. 2, 2023, 6:25 p.m.