-
Notifications
You must be signed in to change notification settings - Fork 14
/
post_correct.Rd
134 lines (118 loc) · 4.12 KB
/
post_correct.Rd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/post_correction.R
\name{post_correct}
\alias{post_correct}
\title{Run Post-correction for Approximate MCMC using \eqn{\psi}-APF}
\usage{
post_correct(
model,
mcmc_output,
particles,
threads = 1L,
is_type = "is2",
seed = sample(.Machine$integer.max, size = 1)
)
}
\arguments{
\item{model}{Model of class \code{nongaussian} or \code{ssm_nlg}.}
\item{mcmc_output}{An output from \code{run_mcmc} used to compute the MAP
estimate of theta.
While the intended use assumes this is from approximate MCMC, it is not
actually checked, i.e., it is also possible to input previous
(asymptotically) exact output.}
\item{particles}{Number of particles for \eqn{\psi}-APF (positive integer).
Suitable values depend on the model and the data, but often relatively
small value less than say 50 is enough. See also \code{suggest_N}}
\item{threads}{Number of parallel threads (positive integer, default is 1).}
\item{is_type}{Type of IS-correction. Possible choices are
\code{"is3"} for simple importance sampling (weight is computed for each
MCMC iteration independently),
\code{"is2"} for jump chain importance sampling type weighting (default), or
\code{"is1"} for importance sampling type weighting where the number of
particles used forweight computations is proportional to the length of the
jump chain block.}
\item{seed}{Seed for the C++ RNG (positive integer).}
}
\value{
The original object of class \code{mcmc_output} with updated
weights, log-posterior values and state samples or summaries (depending on
the \code{mcmc_output$mcmc_type}).
}
\description{
Function \code{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
\code{weights}, \code{posterior}, \code{alpha}, \code{alphahat}, and
\code{Vt} are updated (depending on the original output type).
}
\examples{
\donttest{
Sys.setenv("OMP_THREAD_LIMIT" = 2) # For CRAN
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_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()
}
}
\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
}