rm(list=ls())
set.seed(2026)  
##
## Y|\theta ~ Poisson(theta)
## \theta ~ HN(\sigma = 2)
##

##
## 1. Aproximação discreta da posteriori
##

## quantis por interpolação linear da CDF discreta
qdiscreta <- function(probs, theta, p.theta) {
    pac.theta <- cumsum(p.theta)/sum(p.theta)
    approx(pac.theta, theta, xout = probs, rule = 2)$y
}

##
## 1.1 Grade de 5 pontos
##
th <- c(0.5, 1.5, 2.5, 3.5, 4.5)
th.prior <- dnorm(th, sd = 2) ## 2* é cte!
(th.prior <- th.prior/sum(th.prior))

y <- c(0,3,2,1,0,2)
Sy <- sum(y)
n <- length(y)
th.vero <- exp(Sy*log(th) - n*th)
(th.vero <- th.vero/sum(th.vero))

th.prod <- th.prior*th.vero
(th.post <- th.prod/sum(th.prod))

par(mfrow=c(1,2), mar=c(3,3,0.5,0.5), mgp=c(1.8,0.8,0))
plot(th-0.05, th.prior, type="h", col="red", xlim=c(0,5), ylim=c(0, 1),
     xlab=expression(theta), ylab="probabilidade (grade 5 pontos)")
lines(th, th.vero, type = "h", col="blue")
lines(th+0.05, th.post, type = "h", col="black")
legend("topright", legend=c("priori","verossimilhança","posteriori"),
       col=c("red","blue","black"), lty=1)

plot(th, cumsum(th.prior), type="s", col="red", xlim=c(0,5), ylim=c(0, 1),
     xlab=expression(theta), ylab="CDF (grade 5 pontos)")
points(th, cumsum(th.prior), pch=19, col="red")
points(c(th[1], th), cumsum(c(0, th.prior)), col="red")
lines(th, cumsum(th.post), type="s")
points(th, cumsum(th.post), pch=19)
points(c(th[1], th), cumsum(c(0, th.post)))
par(mfrow=c(1,1))

resumo.grade5 <- list(media = sum(th * th.post),
                      moda = th[which.max(th.post)],
                      P2 = sum(th.post[th > 2]),
                      IC = qdiscreta(c(0.025, 0.975), th, th.post))
unlist(resumo.grade5)

##
## 1.2 Grade fina
##
gr <- 0.1
thF <- seq(0.01, 6, by = gr)
thF.prior <- 2*dnorm(thF, sd = 2)
thF.prior <- thF.prior/sum(thF.prior)
thF.vero <- sapply(thF, function(th) prod(dpois(y, th)))
thF.vero <- thF.vero/sum(thF.vero)  
thF.post <- thF.prior*thF.vero/sum(thF.prior*thF.vero)

par(mfrow=c(1,2), mar=c(3,3,0.5,0.5), mgp=c(1.8,0.8,0))
plot(thF-0.05, thF.prior, type="h", col="red", xlim=c(0,5), ylim=c(0, max(thF.post)),
     xlab=expression(theta), ylab="probabilidade (grade fina)")
lines(thF, thF.vero, type = "h", col="blue")
lines(thF+0.05, thF.post, type = "h", col="black")
legend("topright", legend=c("priori","verossimilhança","posteriori"),
       col=c("red","blue","black"), lty=1)

plot(thF, cumsum(thF.prior), type="s", col="red", xlim=c(0,5), ylim=c(0, 1),
     xlab=expression(theta), ylab="CDF (grade fina)")
lines(thF, cumsum(thF.post), type="s")
par(mfrow=c(1,1))

resumo.gradefina <- list(
    media = sum(thF*thF.post),
    moda  = thF[which.max(thF.post)],
    P2    = sum(thF.post[thF>2]),
    IC    = qdiscreta(c(0.025, 0.975), thF, thF.post))
unlist(resumo.gradefina)

##
## Comparando resultados das duas grades
##
par(mfrow=c(1,2), mar=c(3,3,0.5,0.5), mgp=c(1.8,0.8,0))
plot(th-0.05, th.prior, type="h", col="red", xlim=c(0,5), ylim=c(0, 1),
     xlab=expression(theta), ylab="probabilidade (grade 5 pontos)")
lines(th+0.05, th.post, type = "h", col="black")
legend("topright", legend=c("priori","posteriori"), col=c("red","black"), lty=1)
##
plot(thF-0.01, thF.prior, type="h", col="red", xlim=c(0,6), ylim=c(0, max(thF.post)),
     xlab=expression(theta), ylab="probabilidade (grade fina)")
lines(thF+0.01, thF.post, type="h", col="black")
legend("topright", legend=c("priori","posteriori"), col=c("red","black"), lty=1)
par(mfrow=c(1,1))
##
## 2. Distribuição contínua posteriori
##
##
## Outras aproximações para a posteriori
##

## função auxiliar para plotar a posteriori contínua e a grade fina
plot.disc <- function(xlim = c(0, 6), ...) {
    ## divide probabilidades da grade fina pelo espaçamento para comparar com contínuas
    plot(thF-0.01, thF.post/gr, type="h", col="gray", xlim=xlim,
         ylim=c(0, 1.1*max(thF.post)/gr),
         xlab=expression(theta), ylab= expression(paste("P[", theta,"]")), ...)
    dHN <- function(x, sigma) 2*dnorm(x, sd = sigma)
    curve(dHN(x, sigma = 2), from = 0, to = 6, n = 601, col = "red", lwd = 2, add = TRUE)
}
##plot.disc()
## Numerador do Teorema de Bayes
## f(y|theta) * f(theta) 
num <- function(theta, Sy, n, sigma, log = FALSE) {
    logvero <- Sy*log(theta) - n*theta
    logprior <- dnorm(theta, sd = sigma, log = TRUE) # + log(2) 
    if(log) return(logprior+logvero)
    else return(exp(logprior+logvero))
}

##
## 2.1 Integração Numérica
##
den <- function(Sy, n, sigma) {
    integrate(num, lower = 0, upper = Inf, Sy = Sy, n = n, sigma = sigma)$value
}

dpost <- function(theta, Sy, n, sigma) {
    Num <- num(theta, Sy = Sy, n = n, sigma = sigma)
    Den <- den(Sy = Sy, n = n, sigma = sigma)
    Num/Den
}

## Verificando que a integral da posteriori é 1
integrate(dpost, lower = 0, upper = Inf, Sy = Sy, n = n, sigma = 2)$value

## CDF (por integração) e quantis (por inversão numérica da CDF)
ppost <- function(q, Sy, n, sigma) {
    integrate(dpost, lower = 0, upper = q, Sy = Sy, n = n, sigma = sigma)$value
}
qpost <- function(p, Sy, n, sigma) {
    uniroot(function(th) ppost(th, Sy = Sy, n = n, sigma = sigma) - p,
            lower = 1e-6, upper = 30)$root
}

plot.disc()
curve(dpost(x, Sy = Sy, n = n, sigma = 2), from = 0, to = 6, n = 601, add = TRUE, lwd = 2)

resumo.num <- list(
     media = integrate(function(th) th*dpost(th, Sy = Sy, n = n, sigma = 2), 0, Inf)$value,
     moda  = optimize(function(th) dpost(th, Sy = Sy, n = n, sigma = 2),
                       interval = c(0.001, 10), maximum = TRUE)$maximum,
     P2    = 1 - ppost(2, Sy = Sy, n = n, sigma = 2),
     IC    = c(qpost(0.025, Sy = Sy, n = n, sigma = 2),
               qpost(0.975, Sy = Sy, n = n, sigma = 2)))
unlist(resumo.num)

##
## 2.2 Aproximação normal
##
(approxNorm <- optim(0.5, fn = num, Sy = Sy, n = n, sigma = 2, log = TRUE,
                     lower = 0, upper = 10, method = "Brent",
                     control = list(fnscale = -1), hessian = TRUE))
(mu.hat <- approxNorm$par)
(sigma.hat <- drop(sqrt(-1/approxNorm$hessian)))

plot.disc()
curve(dnorm(x, mean = mu.hat, sd = sigma.hat), from = 0, to = 6, n = 601,
      col = "blue", lwd = 2, add = TRUE)

## IC por quantis da normal (pode incluir valores negativos...)
resumo.norm <- list(
     media = mu.hat,
     moda  = mu.hat,
     P2    = 1 - pnorm(2, mean = mu.hat, sd = sigma.hat),
     IC    = qnorm(c(0.025, 0.975), mean = mu.hat, sd = sigma.hat))
unlist(resumo.norm)

##
## 2.3 Variational Bayes
##

integra <- function(theta, VBfun, Sy, n, sigma, lam) {
    VBfun(theta, lam) * (num(theta, Sy = Sy, n = n, sigma = sigma, log = TRUE) -
                         VBfun(theta, lam, log = TRUE))
}

## Os limites de integração vêm dos quantis da própria q (qVBfun): a integral é
## feita onde q tem massa.  Inf causa erro em alguns valores de lam.
ELBO <- function(lam, VBfun, qVBfun, Sy, n, sigma, eps = 1e-10) {
    lim <- qVBfun(c(eps, 1 - eps), lam)
    integrate(integra, lower = lim[1], upper = lim[2],
              VBfun = VBfun, Sy = Sy, n = n, sigma = sigma, lam = lam)$value
}

## Pontos iniciais baseados nos dados: a verossimilhança é proporcional a uma
## Gama(Sy+1, n) e a posteriori é esse núcleo vezes o fator suave da priori. 
a0 <- Sy + 1
b0 <- n
mean0 <- a0/b0
sd0 <- sqrt(a0)/b0
sdlog0 <- sqrt(log(1 + 1/a0))   ## casamento de momentos da log-normal

## opções para optim() (maximização do ELBO)
ctrl <- list(fnscale = -1, reltol = 1e-12)

##
## 2.3.1 Família Variacional Normal Truncada (em zero):
##       theta|y ~~ N(muv, sigma_v^2) I(theta>0)
##
VB.Tnorm <- function(theta, lam, log = FALSE) {
    ## lam = c(muv, log(sigmav))
    muv <- lam[1]
    sigmav <- exp(lam[2])
    res <- -log(sigmav) + dnorm((theta - muv)/sigmav, log = TRUE) -
        pnorm((0 - muv)/sigmav, lower.tail = FALSE, log.p = TRUE)
    if(log) return(res)
    else return(exp(res))
}
qVB.Tnorm <- function(p, lam) {
    muv <- lam[1]
    sigmav <- exp(lam[2])
    lkappa <- pnorm(0, muv, sigmav, lower.tail = FALSE, log.p = TRUE)
    qnorm(lkappa + log1p(-p), muv, sigmav, lower.tail = FALSE, log.p = TRUE)
}

(fit.vb.norm <- optim(c(m = mean0, ls = log(sd0)), fn = ELBO,
                       VBfun = VB.Tnorm, qVBfun = qVB.Tnorm, Sy = Sy, n = n, sigma = 2,
                       control = ctrl))

(m.vb.norm <- unname(fit.vb.norm$par[1]))
(s.vb.norm <- unname(exp(fit.vb.norm$par[2])))

## Comparação com aprox. normal (2.2) centrada na moda, curvatura da log-posteriori na moda.
## VB minimiza KL(q||posteriori) em vez de usar a curvatura
rbind("aprox. normal" = c(m = mu.hat, s = sigma.hat),
      "VB normal"     = c(m = m.vb.norm, s = s.vb.norm))

plot.disc()
curve(dnorm(x, mean = mu.hat, sd = sigma.hat), from = 0, to = 6, n = 601,
      col = "blue", lwd = 2, add = TRUE)
curve(VB.Tnorm(x, fit.vb.norm$par), from = 0, to = 6, n = 601,
      col = "cyan3", lwd = 2, lty = 2, add = TRUE)

## média, moda, P(theta>2) e IC de theta sob a normal truncada
Ptr.vb.norm <- pnorm(0, mean = m.vb.norm, sd = s.vb.norm, lower.tail = FALSE)
(resumo.vb.norm <- list(
     media = m.vb.norm + s.vb.norm*dnorm(m.vb.norm/s.vb.norm)/Ptr.vb.norm,
     moda  = max(m.vb.norm, 0),
     P2    = pnorm(2, mean = m.vb.norm, sd = s.vb.norm, lower.tail = FALSE)/Ptr.vb.norm,
     IC    = qnorm(pnorm(0, mean = m.vb.norm, sd = s.vb.norm) + c(0.025, 0.975)*Ptr.vb.norm,
                   mean = m.vb.norm, sd = s.vb.norm)))
unlist(resumo.vb.norm)
##
## 2.3.2 Família Variacional Log-Normal:
##       theta|y ~~ LN(muv, sigma_v^2)
##
VB.lognorm <- function(theta, lam, log = FALSE) {
    ## lam = c(muv, log(sigmav))
    muv <- lam[1]
    sigmav <- exp(lam[2])
    dlnorm(theta, meanlog = muv, sdlog = sigmav, log = log)
}
qVB.lognorm <- function(p, lam) qlnorm(p, meanlog = lam[1], sdlog = exp(lam[2]))

(fit.vb <- optim(c(mu = log(mean0) - sdlog0^2/2, ls = log(sdlog0)), fn = ELBO,
                  VBfun = VB.lognorm, qVBfun = qVB.lognorm, Sy = Sy, n = n, sigma = 2,
                  control = ctrl))

(mu.vb <- fit.vb$par[1])
(s.vb  <- exp(fit.vb$par[2]))

plot.disc()
curve(dlnorm(x, meanlog = mu.vb, sdlog = s.vb), from = 0.001, to = 6, n = 601,
      col = "purple", lwd = 2, add = TRUE)

## média, moda, P(theta>2) e IC de theta sob a aproximação log-normal
(resumo.vb <- list(
     media = exp(mu.vb + s.vb^2/2),
     moda  = exp(mu.vb - s.vb^2),
     P2    = plnorm(2, meanlog = mu.vb, sdlog = s.vb, lower = FALSE),
     IC    = qlnorm(c(0.025, 0.975), meanlog = mu.vb, sdlog = s.vb)))
unlist(resumo.vb)
##
## 2.3.3 Família variacional: 
##       theta|y ~~ Gama(av, bv)
## 
VB.gama <- function(theta, lam, log = FALSE) {
    ## lam = c(log(shape), log(rate))
    shape <- exp(lam[1])
    rate <- exp(lam[2])
    dgamma(theta, shape = shape, rate = rate, log = log)
}
qVB.gama <- function(p, lam) qgamma(p, shape = exp(lam[1]), rate = exp(lam[2]))

(fit.vb.gama <- optim(c(la = log(a0), lb = log(b0)), fn = ELBO,
                       VBfun = VB.gama, qVBfun = qVB.gama, Sy = Sy, n = n, sigma = 2,
                       control = ctrl))

(a.vb <- unname(exp(fit.vb.gama$par[1])))
(b.vb <- unname(exp(fit.vb.gama$par[2])))

## optim devolve convergence = 0 quando converge
stopifnot(fit.vb.norm$convergence == 0,
          fit.vb$convergence == 0,
          fit.vb.gama$convergence == 0)

plot.disc()
curve(dgamma(x, shape = a.vb, rate = b.vb), from = 0.001, to = 6, n = 601,
      col = "darkgreen", lwd = 2, add = TRUE)

## média, moda, P(theta>2) e IC de theta sob a aproximação gama
(resumo.vb.gama <- list(
     media = a.vb/b.vb,
     moda  = max(a.vb - 1, 0)/b.vb,
     P2    = pgamma(2, shape = a.vb, rate = b.vb, lower.tail = FALSE),
     IC    = qgamma(c(0.025, 0.975), shape = a.vb, rate = b.vb)))
unlist(resumo.vb.gama)

##
## Comparando as famílias:
## KL(q||posteriori) = log p(y) - ELBO 
## log p(y) = log den(),
logpy <- log(den(Sy = Sy, n = n, sigma = 2))
## Nesse exemplo simples é (apenas 1 parâmetro) permite obter a KL exata de cada aproximação.
##
## Como referência, inclui-se a normal de 2.2 (truncada em theta > 0, como em 2.3.1)
elbo.laplace <- ELBO(c(mu.hat, log(sigma.hat)), VB.Tnorm, qVB.Tnorm,
                     Sy = Sy, n = n, sigma = 2)
ELBOs <- c("aprox. normal" = elbo.laplace, "VB normal" = fit.vb.norm$value,
           "VB log-normal" = fit.vb$value, "VB gama" = fit.vb.gama$value)
rbind(ELBO = ELBOs, KL = logpy - ELBOs)

## Comparando as posterioris das famílias variacionais
plot.disc()
curve(VB.Tnorm(x, fit.vb.norm$par), from = 0, to = 6, n = 601,
      col = "cyan3", lwd = 2, lty = 2, add = TRUE)
curve(dlnorm(x, meanlog = mu.vb, sdlog = s.vb), from = 0.001, to = 6, n = 601,
      col = "purple", lwd = 2, add = TRUE)
curve(dgamma(x, shape = a.vb, rate = b.vb), from = 0.001, to = 6, n = 601,
      col = "darkgreen", lwd = 2, lty = 2, add = TRUE)
legend("topright",
       legend=c("VB normal","VB log-normal","VB gama"),
       col=c("cyan3","purple","darkgreen"), lty=c(2,1,2), lwd=2)

##
## 2.4 MCMC básico (Metrópolis)
##
burnin <- 1000 ; thin <- 10 ; niter <- 101000
theta <- numeric(niter)
theta[1] <- 0.5
acc <- 0
for (i in 2:niter) {
    theta.p <- rnorm(1, mean = theta[i-1], sd = 1.5)
    if (theta.p > 0) {
        p1  <- num(theta.p, Sy = Sy, n = n, sigma = 2)
        p2  <- num(theta[i-1], Sy = Sy, n = n, sigma = 2)
        r <- p1/p2
        if (runif(1) < r) {
            theta[i] <- theta.p
            acc <- acc + 1
        } else {
            theta[i] <- theta[i-1]
        }
    } else {
        theta[i] <- theta[i-1]
    }
}
acc/niter
theta <- theta[-(1:burnin)]
theta <- theta[seq(1, length(theta), by = thin)]

par(mfrow = c(1, 2))
plot(theta, type = "l", xlab = "iteracao", ylab = expression(theta))
acf(theta, lag.max = 20, main = "ACF")
par(mfrow = c(1, 1))

plot.disc()
lines(density(theta), col = "tomato", lwd = 2)

theta.m <- theta  
(ess.m <- coda::effectiveSize(theta.m))

dd.m <- density(theta.m)
resumo.mcmc1 <- list(
     media = mean(theta.m),
     moda  = dd.m$x[which.max(dd.m$y)],
     P2    = mean(theta.m > 2),
     IC    = unname(quantile(theta.m, c(0.025, 0.975))))

unlist(resumo.mcmc1)

## Com reparametrização (I): passeio em log(theta), a mesma reparametrização da
## log-normal do VB (2.3.2)
burnin <- 1000 ; thin <- 10 ; niter <- 101000
theta <- numeric(niter)
theta[1] <- 0.5
acc <- 0
for (i in 2:niter) {
    theta.p <- exp(rnorm(1, mean = log(theta[i-1]), sd = 1.2))
    p1  <- num(theta.p, Sy = Sy, n = n, sigma = 2)
    p2  <- num(theta[i-1], Sy = Sy, n = n, sigma = 2)
    r <- (p1/p2) * (theta.p/theta[i-1])  ## jacobiano da reparametrização
    if (runif(1) < r) {
        theta[i] <- theta.p
        acc <- acc + 1
    } else {
        theta[i] <- theta[i-1]
    }
}
acc/niter
theta <- theta[-(1:burnin)]
theta <- theta[seq(1, length(theta), by = thin)]

par(mfrow = c(1, 2))
plot(theta, type = "l", xlab = "iteracao", ylab = expression(theta))
acf(theta, lag.max = 20, main = "ACF")
par(mfrow = c(1, 1))

plot.disc()
lines(density(theta.m), col = "tomato", lwd = 2)
lines(density(theta), col = "orange", lwd = 2, lty = 2)

(ess.m.log <- coda::effectiveSize(theta))

dd.log <- density(theta)
resumo.mcmc2 <- list(
     media = mean(theta),
     moda  = dd.log$x[which.max(dd.log$y)],
     P2    = mean(theta > 2),
     IC    = unname(quantile(theta, c(0.025, 0.975))))

unlist(resumo.mcmc2)

c(ESS.theta = ess.m, ESS.log.theta = ess.m.log)

## Todos métodos

par(mfrow=c(1,1))
plot.disc(xlim = c(0, 3.5))
curve(dpost(x, Sy = Sy, n = n, sigma = 2), from = 0, to = 3.5, n = 201, add = TRUE, lwd = 2)
curve(dnorm(x, mean = mu.hat, sd = sigma.hat), from = 0, to = 3.5, n = 201,
      col = "blue", lwd = 2, add = TRUE)
lines(density(theta.m), col = "tomato", lwd = 2)
lines(density(theta), col = "orange", lwd = 2, lty = 2)
curve(VB.Tnorm(x, fit.vb.norm$par), from = 0, to = 3.5, n = 201,
      col = "cyan3", lwd = 2, lty = 2, add = TRUE)
curve(dlnorm(x, meanlog = mu.vb, sdlog = s.vb), from = 0.001, to = 3.5, n = 201,
      col = "purple", lwd = 2, add = TRUE)
curve(dgamma(x, shape = a.vb, rate = b.vb), from = 0.001, to = 3.5, n = 201,
      col = "darkgreen", lwd = 2, lty = 2, add = TRUE)

legend("topright",
       legend=c("priori","grade fina","integração numérica","aprox. normal",
                "VB normal","VB log-normal","VB gama",
                "MCMC (passeio em theta)","MCMC (passeio em log theta)"),
       col=c("red","gray","black","blue","cyan3","purple","darkgreen","tomato","orange"),
       lty=c(1,1,1,1,2,1,2,1,2), lwd=2)

##

resumo <- rbind(
    "grade 5 pontos"      = unlist(resumo.grade5),
    "grade fina"          = unlist(resumo.gradefina),
    "integração numérica" = unlist(resumo.num),
    "aprox. normal"       = unlist(resumo.norm),
    "VB normal"           = unlist(resumo.vb.norm),
    "VB log-normal"       = unlist(resumo.vb),
    "VB gama"             = unlist(resumo.vb.gama),
    "MCMC (theta)"        = unlist(resumo.mcmc1),
    "MCMC (log theta)"    = unlist(resumo.mcmc2)
)
colnames(resumo) <- c("media", "moda", "P(theta>2)", "IC 2.5%", "IC 97.5%")
round(resumo, 3)

##
## 3. Distribuição preditiva a posteriori
##
## Qualquer aproximação da posteriori que permita simular $\theta$ gera a preditiva
## de uma nova observação.
## Para cada theta simulado, sorteia-se Y_p ~ Poisson(theta).
## Aqui o objetivo é apenas mostrar como obter a preditiva a partir de cada aproximação;
## a comparação entre os métodos é feita na posteriori.

## Como referência (sem erro de simulação) a preditiva exata é
##     P(Y_p = k | y) = \int dpois(k, theta) * dpost(theta) dtheta
nsim <- 10000   ## mesmo tamanho das cadeias MCMC (após burnin e thinning)

## normal truncada em theta > 0 por inversão da CDF (rpois() exige theta > 0; a normal
## de 2.2 tem P(theta < 0) > 0, e a de 2.3.1 já é truncada por construção)
rnormtr <- function(n, m, s) {
    qnorm(runif(n, min = pnorm(0, mean = m, sd = s), max = 1), mean = m, sd = s)
}

## integração numérica: inversão da CDF integrada numericamente numa grade fina
## (qdiscreta interpola a CDF acumulada; qpost com uniroot seria lento demais aqui)
thG <- seq(0.001, 12, length.out = 6000)
pG  <- dpost(thG, Sy = Sy, n = n, sigma = 2)
## descarta as pontas da grade com massa desprezível (< 1e-12 por ponto): a CDF acumulada
## satura em 1 já a partir de theta ~ 7.6 e os empates fariam approx() emitir um aviso
keep <- pG/sum(pG) > 1e-12
thG <- thG[keep] ; pG <- pG[keep]

## (na ordem das seções: VB antes de MCMC)
theta.sim <- list(
    "grade 5 pontos"      = sample(th, nsim, replace = TRUE, prob = th.post),
    "grade fina"          = sample(thF, nsim, replace = TRUE, prob = thF.post),
    "integração numérica" = qdiscreta(runif(nsim), thG, pG),
    "aprox. normal"       = rnormtr(nsim, mu.hat, sigma.hat),
    "VB normal"           = rnormtr(nsim, m.vb.norm, s.vb.norm),
    "VB log-normal"       = rlnorm(nsim, meanlog = mu.vb, sdlog = s.vb),
    "VB gama"             = rgamma(nsim, shape = a.vb, rate = b.vb),
    "MCMC (theta)"        = theta.m,
    "MCMC (log theta)"    = theta)
yp.sim <- lapply(theta.sim, function(t) rpois(length(t), t))

## probabilidades preditivas; valores >= ypmax são agrupados em uma única categoria
## (factor() é aplicado aos sorteios brutos, não às proporções)
ypmax <- 8
predtab <- function(yp) {
    yp |> pmin(ypmax) |> factor(levels = 0:ypmax) |> table() |> prop.table() |> unclass()
}
pred.exata <- sapply(0:(ypmax - 1), function(k) {
    integrate(function(t) dpois(k, t)*dpost(t, Sy = Sy, n = n, sigma = 2), 0, Inf)$value
})
pred <- rbind("exata (integração)" = c(pred.exata, 1 - sum(pred.exata)),
              t(sapply(yp.sim, predtab)))
colnames(pred) <- c(0:(ypmax - 1), paste0(ypmax, "+"))
round(pred, 3)

cores <- c("sienna", "gray40", "black", "blue", "cyan3", "purple", "darkgreen",
           "tomato", "orange")
par(mfrow = c(1, 1), mar = c(3, 3, 0.5, 0.5), mgp = c(1.8, 0.8, 0), lend = "butt")
plot(0:ypmax, pred[1, ], type = "h", lwd = 10, col = "gray85", ylim = c(0, max(pred)),
     xlab = expression(y[p]), ylab = "probabilidade preditiva")
matlines(0:ypmax, t(pred[-1, ]), type = "b", pch = 19, cex = 0.6, lty = 1, col = cores)
legend("topright", legend = c("exata", rownames(pred)[-1]), lty = 1,
       col = c("gray85", cores), lwd = c(10, rep(1, 9)), cex = 0.7)