## S3 method for class 'nongaussian'
run_mcmc(
  model,
  iter,
  particles,
  output_type = "full",
  mcmc_type = "is2",
  sampling_method = "psi",
  burnin = floor(iter/2),
  thin = 1,
  gamma = 2/3,
  target_acceptance = 0.234,
  S,
  end_adaptive_phase = FALSE,
  local_approx = TRUE,
  threads = 1,
  seed = sample(.Machine$integer.max, size = 1),
  max_iter = 100,
  conv_tol = 1e-08,
  ...
)

set.seed(1)
n <- 50 
slope <- cumsum(c(0, rnorm(n - 1, sd = 0.001)))
level <- cumsum(slope + c(0, rnorm(n - 1, sd = 0.2)))
y <- rpois(n, exp(level))
poisson_model <- bsm_ng(y, 
  sd_level = halfnormal(0.01, 1), 
  sd_slope = halfnormal(0.01, 0.1), 
  P1 = diag(c(10, 0.1)), distribution = "poisson")
  
# Note small number of iterations for CRAN checks
mcmc_out <- run_mcmc(poisson_model, iter = 1000, particles = 10, 
  mcmc_type = "da")
summary(mcmc_out, what = "theta", return_se = TRUE)

set.seed(123)
n <- 50
sd_level <- 0.1
drift <- 0.01
beta <- -0.9
phi <- 5

level <- cumsum(c(5, drift + rnorm(n - 1, sd = sd_level)))
x <- 3 + (1:n) * drift + sin(1:n + runif(n, -1, 1))
y <- rnbinom(n, size = phi, mu = exp(beta * x + level))

model <- bsm_ng(y, xreg = x,
  beta = normal(0, 0, 10),
  phi = halfnormal(1, 10),
  sd_level = halfnormal(0.1, 1), 
  sd_slope = halfnormal(0.01, 0.1),
  a1 = c(0, 0), P1 = diag(c(10, 0.1)^2), 
  distribution = "negative binomial")

# run IS-MCMC
# Note small number of iterations for CRAN checks
fit <- run_mcmc(model, iter = 5000,
  particles = 10, mcmc_type = "is2", seed = 1)

# extract states   
d_states <- as.data.frame(fit, variable = "states", time = 1:n)

library("dplyr")
library("ggplot2")

 # compute summary statistics
level_sumr <- d_states %>% 
  filter(variable == "level") %>%
  group_by(time) %>%
  summarise(mean = Hmisc::wtd.mean(value, weight, normwt = TRUE), 
    lwr = Hmisc::wtd.quantile(value, weight, 
      0.025, normwt = TRUE), 
    upr = Hmisc::wtd.quantile(value, weight, 
      0.975, normwt = TRUE))

# visualize
level_sumr %>% ggplot(aes(x = time, y = mean)) + 
  geom_line() +
  geom_line(aes(y = lwr), linetype = "dashed", na.rm = TRUE) +
  geom_line(aes(y = upr), linetype = "dashed", na.rm = TRUE) +
  theme_bw() + 
  theme(legend.title = element_blank()) + 
  xlab("Time") + ylab("Level")

# theta
d_theta <- as.data.frame(fit, variable = "theta")
ggplot(d_theta, aes(x = value)) + 
 geom_density(aes(weight = weight), adjust = 2, fill = "#92f0a8") + 
 facet_wrap(~ variable, scales = "free") + 
 theme_bw()
 
 
# Bivariate Poisson model:

set.seed(1)
x <- cumsum(c(3, rnorm(19, sd = 0.5)))
y <- cbind(
  rpois(20, exp(x)), 
  rpois(20, exp(x)))

prior_fn <- function(theta) {
  # half-normal prior using transformation
  dnorm(exp(theta), 0, 1, log = TRUE) + theta # plus jacobian term
}

update_fn <- function(theta) {
  list(R = array(exp(theta), c(1, 1, 1)))
}

model <- ssm_mng(y = y, Z = matrix(1,2,1), T = 1, 
  R = 0.1, P1 = 1, distribution = "poisson",
  init_theta = log(0.1), 
  prior_fn = prior_fn, update_fn = update_fn)
  
# Note small number of iterations for CRAN checks
out <- run_mcmc(model, iter = 5000, mcmc_type = "approx")

sumr <- as.data.frame(out, variable = "states") %>% 
  group_by(time) %>% mutate(value = exp(value)) %>%
  summarise(mean = mean(value), 
    ymin = quantile(value, 0.05), ymax = quantile(value, 0.95))
ggplot(sumr, aes(time, mean)) + 
geom_ribbon(aes(ymin = ymin, ymax = ymax),alpha = 0.25) + 
geom_line() + 
geom_line(data = data.frame(mean = y[, 1], time = 1:20), colour = "tomato") + 
geom_line(data = data.frame(mean = y[, 2], time = 1:20), colour = "tomato") +
theme_bw()