predict.mcmc_output  R Documentation 
Draw samples from the posterior predictive distribution for future
time points given the posterior draws of hyperparameters \theta
and
latent state alpha_{n+1}
returned by run_mcmc
.
Function can also be used to draw samples from the posterior predictive
distribution p(\tilde y_1, \ldots, \tilde y_n  y_1,\ldots, y_n)
.
## S3 method for class 'mcmc_output'
predict(
object,
model,
nsim,
type = "response",
future = TRUE,
seed = sample(.Machine$integer.max, size = 1),
...
)
object 
Results object of class 
model 
A 
nsim 
Positive integer defining number of samples to draw. Should be
less than or equal to 
type 
Type of predictions. Possible choices are

future 
Default is 
seed 
Seed for the C++ RNG (positive integer). Note that this affects
only the C++ side, and 
... 
Ignored. 
A data.frame consisting of samples from the predictive posterior distribution.
fitted
for insample predictions.
library("graphics")
y < log10(JohnsonJohnson)
prior < uniform(0.01, 0, 1)
model < bsm_lg(window(y, end = c(1974, 4)), sd_y = prior,
sd_level = prior, sd_slope = prior, sd_seasonal = prior)
mcmc_results < run_mcmc(model, iter = 5000)
future_model < model
future_model$y < ts(rep(NA, 25),
start = tsp(model$y)[2] + 2 * deltat(model$y),
frequency = frequency(model$y))
# use "state" for illustrative purposes, we could use type = "mean" directly
pred < predict(mcmc_results, model = future_model, type = "state",
nsim = 1000)
library("dplyr")
sumr_fit < as.data.frame(mcmc_results, variable = "states") >
group_by(time, iter) >
mutate(signal =
value[variable == "level"] +
value[variable == "seasonal_1"]) >
group_by(time) >
summarise(mean = mean(signal),
lwr = quantile(signal, 0.025),
upr = quantile(signal, 0.975))
sumr_pred < pred >
group_by(time, sample) >
mutate(signal =
value[variable == "level"] +
value[variable == "seasonal_1"]) >
group_by(time) >
summarise(mean = mean(signal),
lwr = quantile(signal, 0.025),
upr = quantile(signal, 0.975))
# If we used type = "mean", we could do
# sumr_pred < pred >
# group_by(time) >
# summarise(mean = mean(value),
# lwr = quantile(value, 0.025),
# upr = quantile(value, 0.975))
library("ggplot2")
rbind(sumr_fit, sumr_pred) >
ggplot(aes(x = time, y = mean)) +
geom_ribbon(aes(ymin = lwr, ymax = upr),
fill = "#92f0a8", alpha = 0.25) +
geom_line(colour = "#92f0a8") +
theme_bw() +
geom_point(data = data.frame(
mean = log10(JohnsonJohnson),
time = time(JohnsonJohnson)))
# Posterior predictions for past observations:
yrep < predict(mcmc_results, model = model, type = "response",
future = FALSE, nsim = 1000)
meanrep < predict(mcmc_results, model = model, type = "mean",
future = FALSE, nsim = 1000)
sumr_yrep < yrep >
group_by(time) >
summarise(earnings = mean(value),
lwr = quantile(value, 0.025),
upr = quantile(value, 0.975)) >
mutate(interval = "Observations")
sumr_meanrep < meanrep >
group_by(time) >
summarise(earnings = mean(value),
lwr = quantile(value, 0.025),
upr = quantile(value, 0.975)) >
mutate(interval = "Mean")
rbind(sumr_meanrep, sumr_yrep) >
mutate(interval =
factor(interval, levels = c("Observations", "Mean"))) >
ggplot(aes(x = time, y = earnings)) +
geom_ribbon(aes(ymin = lwr, ymax = upr, fill = interval),
alpha = 0.75) +
theme_bw() +
geom_point(data = data.frame(
earnings = model$y,
time = time(model$y)))
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.