bssm (version 1.1.5)

suggest_N: Suggest Number of Particles for \(\psi\)-APF Post-correction

Description

Function estimate_N estimates suitable number particles needed for accurate post-correction of approximate MCMC

Usage

suggest_N(
  model,
  mcmc_output,
  candidates = seq(10, 100, by = 10),
  replications = 100,
  seed = sample(.Machine$integer.max, size = 1)
)

Arguments

model

Model of class nongaussian or ssm_nlg.

mcmc_output

An output from run_mcmc used to compute the MAP estimate of theta. While the intended use assumes this is from approximate MCMC, it is not actually checked, i.e., it is also possible to input previous (asymptotically) exact output.

candidates

Vector containing the candidate number of particles to test. Default is seq(10, 100, by = 10).

replications

How many replications should be used for computing the standard deviations? Default is 100.

seed

Seed for the random number generator.

Value

List with suggested number of particles N and matrix containing estimated standard deviations of the log-weights and corresponding number of particles.

Details

Function suggest_N estimates the standard deviation of the logarithm of the post-correction weights at approximate MAP of theta, using various particle sizes and suggest smallest number of particles which still leads standard deviation less than 1. Similar approach was suggested in the context of pseudo-marginal MCMC by Doucet et al. (2015), but see also Section 10.3 in Vihola et al (2020).

References

A. Doucet, M. K. Pitt, G. Deligiannidis, R. Kohn, Efficient implementation of Markov chain Monte Carlo when using an unbiased likelihood estimator, Biometrika, Volume 102, Issue 2, 2015, Pages 295<U+2013>313, https://doi.org/10.1093/biomet/asu075 Vihola, M, Helske, J, Franks, J. Importance sampling type estimators based on approximate marginal Markov chain Monte Carlo. Scand J Statist. 2020; 1<U+2013> 38. https://doi.org/10.1111/sjos.12492

Examples

# NOT RUN {
set.seed(1)
n <- 300
x1 <- sin((2 * pi / 12) * 1:n)
x2 <- cos((2 * pi / 12) * 1:n)
alpha <- numeric(n)
alpha[1] <- 0
rho <- 0.7
sigma <- 2
mu <- 1
for(i in 2:n) {
  alpha[i] <- rnorm(1, mu * (1 - rho) + rho * alpha[i-1], sigma)
}
u <- rpois(n, 50)
y <- rbinom(n, size = u, plogis(0.5 * x1 + x2 + alpha))

ts.plot(y / u)

model <- ar1_ng(y, distribution = "binomial", 
  rho = uniform(0.5, -1, 1), sigma = gamma(1, 2, 0.001),
  mu = normal(0, 0, 10),
  xreg = cbind(x1,x2), beta = normal(c(0, 0), 0, 5),
  u = u)

out_approx <- run_mcmc(model, mcmc_type = "approx", 
 iter = 5000)

estN <- suggest_N(model, out_approx, candidates = seq(10, 50, by = 10))
plot(x = estN$results$N, y = estN$results$sd, type = "b")
estN$N
# }