Suggest Number of Particles for ψ-APF Post-correction
Function estimate_N estimates suitable number particles needed for accurate
post-correction of approximate MCMC
suggest_N( model, mcmc_output, candidates = seq(10, 100, by = 10), replications = 100, seed = sample(.Machine$integer.max, size = 1) )
model |
Model of class |
mcmc_output |
An output from |
candidates |
Vector containing the candidate number of particles to test. Default
is |
replications |
How many replications should be used for computing the standard deviations? Default is 100. |
seed |
Seed for the random number generator. |
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).
List with suggested number of particles N and matrix containing
estimated standard deviations of the log-weights and corresponding number of particles.
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–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– 38. https://doi.org/10.1111/sjos.12492
## 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
## End(Not run)Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.