suggest_N {bssm}

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,
  theta,
  candidates = seq(10, 100, by = 10),
  replications = 100,
  seed = sample(.Machine$integer.max, size = 1)
)

Arguments

model

Model of class nongaussian or ssm_nlg.

theta

A vector of theta corresponding to the model, at which point the standard deviation of the log-likelihood is computed. Typically MAP estimate from the (approximate) MCMC run. Can also be an output from run_mcmc which is then used to compute the MAP estimate of theta.

candidates

A vector of positive integers containing the candidate number of particles to test. Default is seq(10, 100, by = 10).

replications

Positive integer, how many replications should be used for computing the standard deviations? Default is 100.

seed

Seed for the C++ RNG (positive integer).

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

Value

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

References

Doucet, A, Pitt, MK, Deligiannidis, G, Kohn, R (2015). Efficient implementation of Markov chain Monte Carlo when using an unbiased likelihood estimator, Biometrika, 102(2) p. 295-313, https://doi.org/10.1093/biomet/asu075

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

Examples

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 <- 1.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_prior(1, 2, 0.001),
  mu = normal(0, 0, 10),
  xreg = cbind(x1,x2), beta = normal(c(0, 0), 0, 5),
  u = u)

# theta from earlier approximate MCMC run
# out_approx <- run_mcmc(model, mcmc_type = "approx", 
#   iter = 5000) 
# theta <- out_approx$theta[which.max(out_approx$posterior), ]

theta <- c(rho = 0.64, sigma = 1.16, mu = 1.1, x1 = 0.56, x2 = 1.28)

estN <- suggest_N(model, theta, candidates = seq(10, 50, by = 10),
  replications = 50, seed = 1)
plot(x = estN$results$N, y = estN$results$sd, type = "b")
estN$N