View source: R/post_correction.R
post_correct  R Documentation 
\psi
APFFunction post_correct
updates previously obtained approximate MCMC
output with postcorrection weights leading to asymptotically exact
weighted posterior, and returns updated MCMC output where components
weights
, posterior
, alpha
, alphahat
, and
Vt
are updated (depending on the original output type).
post_correct(
model,
mcmc_output,
particles,
threads = 1L,
is_type = "is2",
seed = sample(.Machine$integer.max, size = 1)
)
model 
Model of class 
mcmc_output 
An output from 
particles 
Number of particles for 
threads 
Number of parallel threads (positive integer, default is 1). 
is_type 
Type of IScorrection. Possible choices are

seed 
Seed for the C++ RNG (positive integer). 
The original object of class mcmc_output
with updated
weights, logposterior values and state samples or summaries (depending on
the mcmc_output$mcmc_type
).
Doucet A, Pitt M K, Deligiannidis G, Kohn R (2018). Efficient implementation of Markov chain Monte Carlo when using an unbiased likelihood estimator. Biometrika, 102, 2, 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 < 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)
out_approx < run_mcmc(model, mcmc_type = "approx",
local_approx = FALSE, iter = 50000)
out_is2 < post_correct(model, out_approx, particles = 30,
threads = 2)
out_is2$time
summary(out_approx, return_se = TRUE)
summary(out_is2, return_se = TRUE)
# latent state
library("dplyr")
library("ggplot2")
state_approx < as.data.frame(out_approx, variable = "states") >
group_by(time) >
summarise(mean = mean(value))
state_exact < as.data.frame(out_is2, variable = "states") >
group_by(time) >
summarise(mean = weighted.mean(value, weight))
dplyr::bind_rows(approx = state_approx,
exact = state_exact, .id = "method") >
filter(time > 200) >
ggplot(aes(time, mean, colour = method)) +
geom_line() +
theme_bw()
# posterior means
p_approx < predict(out_approx, model, type = "mean",
nsim = 1000, future = FALSE) >
group_by(time) >
summarise(mean = mean(value))
p_exact < predict(out_is2, model, type = "mean",
nsim = 1000, future = FALSE) >
group_by(time) >
summarise(mean = mean(value))
dplyr::bind_rows(approx = p_approx,
exact = p_exact, .id = "method") >
filter(time > 200) >
ggplot(aes(time, mean, colour = method)) +
geom_line() +
theme_bw()
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.