## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(echo = TRUE, eval = FALSE)

## ----eval=FALSE---------------------------------------------------------------
# library(BayesRTMB)
# 
# code <- rtmb_code(
#   setup = {
#     # Data preprocessing and constant definitions
#   },
#   parameters = {
#     # Declaration of parameters to be estimated
#   },
#   transform = {
#     # Calculation of derived variables using parameters and data
#   },
#   model = {
#     # Definition of priors and likelihoods
#   },
#   generate = {
#     # Calculation of posterior predictions and generated quantities
#   }
# )

## ----eval=FALSE---------------------------------------------------------------
# parameters = {
#   # Scalar parameter (with lower bound of 0)
#   sigma = Dim(lower = 0)
# 
#   # Vector parameter (length N)
#   mu = Dim(N)
# 
#   # Matrix parameter (N rows, M columns)
#   X = Dim(c(N, M))
# 
#   # Array parameter (I x J x K)
#   A = Dim(c(I, J, K))
# }

## ----eval=FALSE---------------------------------------------------------------
# model = {
#   # Priors
#   mu ~ normal(0, 10)
#   sigma ~ exponential(1)
# 
#   # Likelihood
#   Y ~ normal(mu, sigma)
# }

## ----eval=FALSE---------------------------------------------------------------
# model = {
#   lp <- lp + normal_lpdf(mu, 0, 10)
#   lp <- lp + exponential_lpdf(sigma, 1)
#   lp <- lp + normal_lpdf(Y, mu, sigma)
# }

## ----eval=FALSE---------------------------------------------------------------
# setup = {
#   N <- length(Y)  # Get the sample size of the observed data
#   P <- ncol(X)    # Get the number of columns in the design matrix (number of predictors)
# }

## -----------------------------------------------------------------------------
# library(BayesRTMB)
# 
# simulate_exgaussian <- function(N, mu, sigma, lambda) {
#   rnorm(N, mean = mu, sd = sigma) + rexp(N, rate = lambda)
# }
# 
# set.seed(1234)
# 
# true <- c(mu = 0.45, sigma = 0.06, lambda = 8)
# rt <- simulate_exgaussian(
#   N = 500,
#   mu = true["mu"],
#   sigma = true["sigma"],
#   lambda = true["lambda"]
# )
# 
# hist(rt)

## -----------------------------------------------------------------------------
# code_exgaussian <- rtmb_code(
#   setup = {
#     Y <- rt
# 
#     exgaussian_lpdf <- function(x, mu, sigma, lambda, sum = TRUE) {
#       z <- (x - mu) / sigma - lambda * sigma
# 
#       res <-
#         log(lambda) +
#         lambda * (mu - x) +
#         0.5 * (lambda * sigma)^2 +
#         RTMB::pnorm(z, log.p = TRUE)
# 
#       if (sum) sum(res) else res
#     }
#   },
# 
#   parameters = {
#     mu     <- Dim(lower = 0)
#     sigma  <- Dim(lower = 0)
#     lambda <- Dim(lower = 0)
#   },
# 
#   model = {
#     Y ~ exgaussian(mu, sigma, lambda)
# 
#     mu     ~ exponential(1)
#     sigma  ~ exponential(1)
#     lambda ~ exponential(1 / 10)
#   },
# 
#   generate = {
#     log_lik <- exgaussian_lpdf(Y, mu, sigma, lambda, sum = FALSE)
#     report(log_lik)
#   }
# )

## -----------------------------------------------------------------------------
# init_exgaussian <- list(
#   mu = median(rt) - sd(rt) / 2,
#   sigma = sd(rt) / 2,
#   lambda = 1 / (sd(rt) / 2)
# )
# 
# mdl_exgaussian <- rtmb_model(
#   data = list(rt = rt),
#   code = code_exgaussian,
#   init = init_exgaussian
# )
# 
# fit_map <- mdl_exgaussian$optimize(num_estimate = 5)
# 
# est_map <- unlist(fit_map$estimate(pars = "parameters"))[names(true)]
# data.frame(
#   parameter = names(true),
#   true = as.numeric(true),
#   estimate = as.numeric(est_map),
#   error = as.numeric(est_map - true)
# )

## -----------------------------------------------------------------------------
# fit_mcmc <- mdl_exgaussian$sample(parallel = TRUE)
# fit_mcmc
# fit_mcmc$WAIC()
# 
# fit_vb <- mdl_exgaussian$variational()
# fit_vb

## ----eval=FALSE---------------------------------------------------------------
# model = {
#   eta <- rtmb_vector(0, N)
#   for (i in seq_len(N)) {
#     eta[i] <- alpha + X[i, ] %*% beta
#   }
#   Y ~ normal(eta, sigma)
# }

## ----eval=FALSE---------------------------------------------------------------
# model = {
#   logit_x <- rtmb_array(0, dim = c(N_time, C, D))
#   for (t in seq_len(N_time)) {
#     for (c in seq_len(C)) {
#       for (d in seq_len(D)) {
#         logit_x[t, c, d] <- alpha[d] + beta[d] * X[t, c]
#       }
#     }
#   }
# }

## ----eval=FALSE---------------------------------------------------------------
# fit2 <- upgrade_fit(fit)