##
## Pacotes para amostragem MCMC
##
## Pacotes considerados aqui:
##  - MCMCpack
##  - mcmc
##  - JAGS
##  - Stan
##
## Pacotes auxiliares para diagnóstico: coda (usado por MCMCpack, rjags),
##   posterior e bayesplot (ecossistema Stan)
##
## O modelo:
## Y|theta ~ B(n, theta)        (verossimilhança)
## theta ~ U[0,1] \equiv B(1,1) (priori)
## theta|y ~ Be(y+a, (n-y) + b) (mas vamos supor que não temos esta forma fechada)

##
## Fingimos não conhecer a posteriori, mas a usamos para conferir os resultados.
##

## rm(list=ls()))

## No exemplo:
a <- 1; b <-  1
n <- 80; y <- 19 

set.seed(2026)

## Resumos "exatos" da posteriori, para comparação
post.exata <- c(media = (y + a) / (n + a + b),
                dp    = sqrt((y + a) * (n - y + b) / ((n + a + b)^2 * (n + a + b + 1))),
                qbeta(c(0.025, 0.5, 0.975), y + a, n - y + b))
names(post.exata)[3:5] <- c("2.5%", "50%", "97.5%")
post.exata

## Função auxiliar: resume amostra de theta e compara com o exato
resumo <- function(th) {
    th <- as.numeric(th)
    c(media = mean(th), dp = sd(th), quantile(th, c(0.025, 0.5, 0.975)))
}

## OBS:
## Transformações logito e logito inverso:
## qlogis(p) = log(p/(1-p))   e   plogis(x) = exp(x)/(1+exp(x))


##
## 1. Pacote MCMCpack
##

library(MCMCpack)  ## tb carrega "coda"
## Implementa algoritmo de Metrópolis
##
## Usando a função MCMCmetrop1R() que permite "user-defined continuous distribution"
## argumento "fun" : unnormalised (log)density
## theta|y \propto theta^(y+a-1) (1-theta)^(n-y+b-1)
## log(theta|y) \propto (y+a-1) * log(\theta) + (n-y+b-1) * log(1-theta)
##
## Na escala eta = logito(theta) entra o jacobiano |d theta/d eta| = theta(1-theta):
## eta|y \propto theta^(y+a) (1-theta)^(n-y+b),  com theta = plogis(eta)
## (o jacobiano cancela os "-1"; daí apost = y+a e bpost = n-y+b em "fun")
## Ver Jacobiano.txt

## Documentação:
##  - ANTES do MCMC a função faz uma otimização (optim) e usa a inversa da
##    hessiana no máximo (V) como covariância da proposta; "tune" multiplica
##    o desvio padrão da proposta (valores ~ 2.4/sqrt(dim) costumam ser bons)

## usando reparametrização eta = log(theta/(1-theta)) para parâmetro na escala do logito 
fun <- function(eta, apost = y+a, bpost = n-y+b){
    ## o calculo aqui seria:
    #p <- exp(eta)/(1+exp(eta))
    #return(apost * log(p) + bpost * log(1-p))
    ## mas a forma equivalente a seguir é mais estável numericamente
    logp <- ifelse(eta < 0, eta - log1p(exp(eta)), - log1p(exp(- eta)))
    logq <- ifelse(eta < 0, - log1p(exp(eta)), - eta - log1p(exp(- eta)))
    ## logl <- sum(logp[y == 1]) + sum(logq[y == 0])
    apost * logp + bpost * logq
}

ams1 <- MCMCmetrop1R(fun, theta.ini = 0, 
                    ## optim.method = "L-BFGS-B", optim.lower = 0, optim.upper = 1,
                    mcmc = 10000, burnin = 500, thin = 10, tune = 2.4)
str(ams1)
summary(ams1)
plot(ams1)
plot(acf(ams1))
autocorr.plot(ams1)
## A taxa de aceitação correta é a impressa pelo MCMCmetrop1R ao final.
## Com thin = 10, rejectionRate() compara amostras guardadas a cada 10
## iterações (quase sempre diferentes) e superestima a aceitação (~0.99)
1-rejectionRate(ams1)  
effectiveSize(ams1)  
##geweke.diag(ams1)

## Transformando para escala da proporção
ams1.th <- plogis(ams1)
## comparando com resultado teórico
plot(density(ams1.th), main = "MCMCpack")
curve(dbeta(x, y + a, n - y + b), add = TRUE, col = 2, lty = 2)
rbind(MCMCpack = resumo(ams1.th), exato = post.exata)
## 
## O pacote possui funções "pré prontas" para certos modelos
## Nestes casos não é necessário definir "fun" 
##
detach("package:MCMCpack")

##
## 2. Pacote mcmc
##
## Usa algoritmo de Metrópolis
## Assim como em MCMCpack::MCMCmetrop1r() requer "log unnormalized density"

## mcmc::metrop():
#metrop(obj, initial, nbatch, blen = 1, nspac = 1, scale = 1, outfun, debug = FALSE, ...)
require(mcmc)

ams2 <- metrop(fun, initial = 0, nbatch = 1000)
str(ams2)
ams2$accept

# Recomeça onde parou anterior (e desta forma, emula "burnin")
ams2 <- mcmc::metrop(ams2, nbatch = 10000, scale=0.6)
str(ams2)
names(ams2)
ams2$accept
ams2$nbatch
ams2$blen
str(ams2$batch)

with(ams2, plot(ts(batch)))
with(ams2, acf(batch))
with(ams2, plot(density(batch)))

## quantidade/parametro de interesse
## comparando com resultado teórico
with(ams2, plot(density(plogis(batch))))
curve(dbeta(x,y+a,n-y+b), add=T)

## alternativa transformando direto na função
ams2 <- mcmc::metrop(fun, initial = 0, nbatch = 10000, scale=0.6,
                    outfun = \(x){exp(x)/(1+exp(x))})
str(ams2)
names(ams2)
ams2$accept
ams2$nbatch
ams2$blen
str(ams2$batch)

with(ams2, plot(ts(batch)))
with(ams2, acf(batch))
with(ams2, plot(density(batch)))
curve(dbeta(x,y+a,n-y+b), add=T)

ams2 <- mcmc::metrop(fun, initial = 0, nbatch = 1000, scale=0.6,
                    outfun = \(x){c(x, exp(x)/(1+exp(x)))})
ams2 <- mcmc::metrop(ams2, nbatch = 10000, scale=0.6, nspac = 10,
                    outfun = \(x){c(x, exp(x)/(1+exp(x)))})
str(ams2)
str(ams2$batch)
head(ams2$batch)
with(ams2, plot(ts(batch)))
with(ams2, acf(batch))
with(ams2, plot(density(batch[,1])))
with(ams2, plot(density(batch[,2])))
curve(dbeta(x,y+a,n-y+b), add=T)
rbind(mcmc = resumo(ams2$batch[, 2]), exato = post.exata)

## "blen" armazena as MÉDIAS dos lotes ("batch means").
## Servem para estimar o erro de Monte Carlo da média a posteriori.
## OBS: argumentos omitidos (scale, nspac, outfun) são herdados do objeto ams2
ams2.bm <- metrop(ams2, nbatch = 200, blen = 100)
dim(ams2.bm$batch)
colMeans(ams2.bm$batch)                              # estimativas (eta, theta)
apply(ams2.bm$batch, 2, sd) / sqrt(ams2.bm$nbatch)   # EP de Monte Carlo

detach("package:mcmc")

##
## 3. Jags (instalado "stand alone" no sistema)
##    interface com R via via pacotes rjags ou runjags
##
require(rjags)

## Notação BUGS: dbin(p, n) (ordem de dbin != dbinom do R)
modelo.jags <- "model {
  ## verossimilhança [Y|theta]
  y ~ dbin(theta, n)
  ## priori [theta]
  theta ~ dbeta(a, b)
  ## quantidade derivada (não interfere na amostragem)
  phi <- logit(theta)
}"

## 3 cadeias começando em locais diferentes
inis <- list(list(theta = 0.1),
             list(theta = 0.5),
             list(theta = 0.9)
             )

jags <- jags.model(textConnection(modelo.jags),
                   data = list(y = y, n = n, a = a, b = b),
                   n.chains = 3, inits = inis,
                   n.adapt = 1000)   # fase de adaptação dos amostradores

update(jags, n.iter = 1000)          # burn-in (descartado)
ams3 <- coda.samples(jags, c("theta", "phi"), n.iter = 10000)

str(ams3)
plot(ams3)
plot(ams3[[1]])
plot(ams3[[2]])
plot(ams3[[3]])

par(mfrow=c(2,2))
plot(ams3)
str(ams3)
plot(acf(ams3[[1]]))
plot(acf(ams3[[2]]))
plot(acf(ams3[[3]]))
par(mfrow=c(1,1))
autocorr.plot(ams3[[1]])
effectiveSize(ams3)
summary(ams3)
HPDinterval(ams3)
rbind(JAGS = resumo(as.matrix(ams3)[, "theta"]), exato = post.exata)

gelman.diag(ams3)   # Rhat (fator de redução de escala potencial): ~ 1
gelman.plot(ams3)
##
detach("package:rjags")

## pacote/interface alternativa
## require(runjags)
## datalist <- dump.format(list(x = y, n = n))
## params <- c("theta")
## inicial <- dump.format(list(theta = 0.5))
## ## Modelo
## mod <- "model{
## x ~ dbin(theta, n)
## theta ~ dbeta(1, 1)    # a = b = 1
## }"
## ## Ajuste
##m.jags <- run.jags(model = mod,
##                   monitor = params,
##                   data = datalist,
##                   inits = c(inicial, inicial),
##                   n.chains = 2,
##                   burnin = 5000,
##                   thin = 5,
##                   sample = 10000
##)

## Resultados
## m.jags
## qbeta(c(0.025, 0.5, 0.975), a + y, b + n - y)
## plot(m.jags)

##
## 4. Stan  (Hamiltonian Monte Carlo, algoritmo NUTS)
##    interfaces: rstan (dentro do R) ou cmdstanr (via CmdStan, mais atualizado)
##
require(rstan)

options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)    # não recompila o modelo a cada execução

modelo.stan <- "// exemplo binomial
data {
  int<lower=0> n;
  int<lower=0, upper=n> y;
  real<lower=0> a;
  real<lower=0> b;
}
parameters {
  real<lower=0, upper=1> theta;  // Stan transforma internamente para logito
}                                // e inclui o jacobiano automaticamente
model {
  theta ~ beta(a, b);
  y ~ binomial(n, theta);
}
generated quantities {
  real phi = logit(theta);
}"

dados <- list(n = n, y = y, a = a, b = b)

## O modelo é traduzido para C++ e compilado (demora na 1a vez).
## Com NUTS em geral não é necessário "thin"; o default (iter = 2000,
## warmup = iter/2, chains = 4) costuma ser suficiente.
fit <- stan(model_code = modelo.stan, data = dados,
            chains = 4, iter = 2000, seed = 2026, refresh = 0)

fit                         # inclui n_eff e Rhat
traceplot(fit, pars = c("theta", "phi"))
plot(fit, pars = "theta")
check_hmc_diagnostics(fit)  # divergências, profundidade da árvore, E-BFMI

ams4 <- extract(fit)$theta
plot(density(ams4), main = "Stan")
curve(dbeta(x, y + a, n - y + b), add = TRUE, col = 2, lty = 2)
rbind(Stan = resumo(ams4), exato = post.exata)

## Gráficos e diagnósticos do ecossistema Stan (valem para qualquer MCMC)
library(bayesplot)
mcmc_trace(as.array(fit), pars = c("theta", "phi"))
mcmc_acf(as.array(fit), pars = "theta")
## posterior::summarise_draws(as.array(fit))   # Rhat e ESS "bulk"/"tail"

## Exploração interativa (abre o navegador)
if (interactive()) {
    library(shinystan)
    launch_shinystan(as.shinystan(fit))
}