## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----srr-tags, eval = FALSE, echo = FALSE-------------------------------------
#  #' roxygen_block_name
#  #'
#  #' @srrstats {G1.3} Vignette presents consistent statistical terminology to be
#  #' used throughout ernest.
#  #' @srrstats {BS1.2, BS1.2b} Contains examples for specifying priors.
#  #' @srrstats {BS1.3, BS1.3a} Describes how to save and continue previous runs
#  #' with `generate()`.

## ----message=FALSE, echo=FALSE------------------------------------------------
library(ernest)
library(brms)
library(posterior)

## -----------------------------------------------------------------------------
data("epilepsy")
frame <- model.frame(count ~ zAge + zBase * Trt, epilepsy)
X <- model.matrix(count ~ zAge + zBase * Trt, epilepsy)
Y <- model.response(frame)

## -----------------------------------------------------------------------------
poisson_log_lik <- function(predictors, response, link = "log") {
  force(predictors)
  force(response)
  link <- make.link(link)
  
  function(theta) {
    eta <- predictors %*% theta
    mu <- link$linkinv(eta)
    sum(dpois(response, lambda = mu, log = TRUE))
  }
}

epilepsy_log_lik <- create_likelihood(poisson_log_lik(X, Y))
epilepsy_log_lik

epilepsy_log_lik(c(1.94, 0.15, 0.57, -0.20, 0.05))

## -----------------------------------------------------------------------------
epilepsy_log_lik(c(1.94, 0.15, 0.57, -0.20, Inf))

## -----------------------------------------------------------------------------
poisson_vec_lik <- function(predictors, response, link = "log") {
  force(predictors)
  force(response)
  link <- make.link(link)

  function(theta_mat) {
    eta_mat <- predictors %*% t(theta_mat)
    mu_mat <- link$linkinv(eta_mat)
    colSums(apply(mu_mat, 2, \(col) dpois(response, lambda = col, log = TRUE)))
  }
}

epilepsy_vec_lik <- create_likelihood(vectorized_fn = poisson_vec_lik(X, Y))
epilepsy_vec_lik

## -----------------------------------------------------------------------------
theta_mat <- matrix(
  c(1.94, 0.15, 0.57, -0.20, 0.05, 1.94, 0.0, 0.0, 0.0, 0.0),
  byrow = TRUE,
  nrow = 2
)
epilepsy_vec_lik(theta_mat)
epilepsy_log_lik(theta_mat)

## -----------------------------------------------------------------------------
coef_names <- c("Intercept", "zAge", "zBase", "Trt1", "zBase:Trt1")

norm_transform <- function(unit) {
  qnorm(unit, sd = 2.5)
}
custom_prior <- create_prior(norm_transform, names = coef_names)
custom_prior

## -----------------------------------------------------------------------------
model_prior <- create_normal_prior(names = coef_names, sd = 2.5)
model_prior

## -----------------------------------------------------------------------------
multi_ellipsoid()
multi_ellipsoid(enlarge = 1.5)

rwmh_cube()
rwmh_cube(steps = 30, target_acceptance = 0.4)

## -----------------------------------------------------------------------------
sampler <- ernest_sampler(
  epilepsy_log_lik,
  model_prior,
  sampler = rwmh_cube(),
  nlive = 300,
  seed = 42
)
sampler

## -----------------------------------------------------------------------------
run_1k <- generate(sampler, max_iterations = 1000, show_progress = FALSE)
run_1k

## -----------------------------------------------------------------------------
run_1k$samples$log_weight |> summary()

## -----------------------------------------------------------------------------
run <- generate(run_1k, show_progress = FALSE)
run

## -----------------------------------------------------------------------------
summary(run)
plot(run, which = c("weight", "likelihood"))

## -----------------------------------------------------------------------------
sim_run <- calculate(run, ndraws = 1000)
sim_run
plot(sim_run, which = c("weight", "likelihood"))

## -----------------------------------------------------------------------------
as_draws(run) |>
  resample_draws() |>
  summarise_draws()
visualize(run, .which = "trace")
visualize(run, -Intercept, .which = "density")