# post_correct: Run Post-correction for Approximate MCMC using psi-APF In bssm: Bayesian Inference of Non-Linear and Non-Gaussian State Space Models

 post_correct R Documentation

## Run Post-correction for Approximate MCMC using ψ-APF

### Description

Function post_correct updates previously obtained approximate MCMC output with post-correction 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).

### Usage

post_correct(
model,
mcmc_output,
particles,
is_type = "is2",

### References

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, 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 <- 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_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,
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()



bssm documentation built on May 4, 2022, 1:06 a.m.