View source: R/post_correction.R
suggest_N  R Documentation 
Function estimate_N
estimates suitable number particles needed for
accurate postcorrection of approximate MCMC.
suggest_N( model, theta, candidates = seq(10, 100, by = 10), replications = 100, seed = sample(.Machine$integer.max, size = 1) )
model 
Model of class 
theta 
A vector of theta corresponding to the model, at which point
the standard deviation of the loglikelihood is computed. Typically MAP
estimate from the (approximate) MCMC run. Can also be an output from

candidates 
A vector of positive integers containing the candidate
number of particles to test. Default is 
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). 
Function suggest_N
estimates the standard deviation of the
logarithm of the postcorrection 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 pseudomarginal 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 logweights and
corresponding number of particles.
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. 295313, 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. 138. https://doi.org/10.1111/sjos.12492
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[i1], 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
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.