rm(list = ls())

set.seed(2026)  

# Simulação Priori, Verossimilhança => Posteriori

a <- 5     
b <- 10    
n <- 5     
y <- 21    

N_sim <- 2e5   


a_post <- a + y
b_post <- b + n



lambda_grid <- seq(0, 15, by = 0.05)

priori_grid <- dgamma(lambda_grid, shape = a, rate = b)
priori_grid <- priori_grid / sum(priori_grid)

vero_grid <- dpois(y, lambda = n * lambda_grid)
vero_grid <- vero_grid / sum(vero_grid)     

conj_grid <- dgamma(lambda_grid, a, b) * dpois(y, n * lambda_grid)
post_grid <- conj_grid / sum(conj_grid)

x11()

plot(
  lambda_grid, priori_grid,
  type = "h", lwd = 2, col = "red",
  xlim = c(0, 10),
  ylim = c(0, max(priori_grid, vero_grid, post_grid) * 1.1),
  xlab = expression(lambda), ylab = "Probabilidade",
  main = "Priori, Verossimilhanca e Posteriori (via grid)"
)
lines(lambda_grid + 0.015, vero_grid, type = "h", lwd = 2, col = "blue")
lines(lambda_grid + 0.03,  post_grid, type = "h", lwd = 2, col = "black")
legend(
  "topright",
  legend = c("Priori", "Verossimilhanca", "Posteriori"),
  col = c("red", "blue", "black"), lwd = 2, bty = "n"
)


# =====================================================================

# SImulando da posteriori e comparando com a curva analítica e com o 
#resultado da posteriori calculada anteriormente pelo teorema de Bayes

lambda_post_sim <- rgamma(N_sim, shape = a_post, rate = b_post)

x11()

hist(
  lambda_post_sim, probability = TRUE, breaks = 60,
  col = "gray85", border = "white",
  xlim = c(0, 10),
  main = "Posteriori de lambda: simulacao x formula fechada",
  xlab = expression(lambda),
)
curve(
  dgamma(x, shape = a_post, rate = b_post),
  add = TRUE, lwd = 2, col = "black"
)


lines(lambda_grid, post_grid / diff(lambda_grid)[1], type = "l", lwd = 2, col = "red", lty = 2)
legend(
  "topright",
  legend = c("Simulacao (Monte Carlo)", "Densidade Gama(a+y, b+n)", "Posteriori via grid"),
  col = c("gray70", "black", "red"), lwd = 2, lty = c(NA, 1, 2), pch = c(15, NA, NA),
  bty = "n"
)

#Média da posteriori:

#Analítica:
a_post / b_post

#Via Simulação:
mean(lambda_post_sim)

#Variância da posteriori:

#Analítica:
a_post / b_post^2

#Via Simulação: 
var(lambda_post_sim)

# =====================================================================

# Simulando lambda e utilizando este lambda na simulação de Y

lambda_priori_sim <- rgamma(N_sim, shape = a, rate = b)

y_priori_sim <- rpois(N_sim, lambda = n * lambda_priori_sim)

tab_priori_sim <- prop.table(table(y_priori_sim))

tab_priori_sim

valores_priori <- as.numeric(names(tab_priori_sim))

prob_priori_bn <- b / (b + n)     
dens_priori_bn <- dnbinom(valores_priori, size = a, prob = prob_priori_bn)

plot(
  valores_priori, tab_priori_sim,
  type = "h", lwd = 2, col = "darkorange",
  xlim = c(0, max(valores_priori)),
  ylim = c(0, max(tab_priori_sim, dens_priori_bn) * 1.1),
  xlab = "y", ylab = "P(Y = y)",
  main = "Preditiva a priori: simulacao x BinNeg(a, b/(b+n))"
)
lines(valores_priori + 0.3, dens_priori_bn, type = "h", lwd = 2, col = "steelblue")
legend(
  "topright",
  legend = c("Simulacao Monte Carlo", "BinNeg(size=a, prob=b/(b+n))"),
  col = c("darkorange", "steelblue"), lwd = 2, bty = "n"
)


# =====================================================================

# Simulando Y's utilizando lambdas da POSTERIORI de lambda

m <- 1   

ynovo_sim <- rpois(N_sim, lambda = lambda_post_sim)   

tab_post_sim <- prop.table(table(ynovo_sim))

tab_post_sim

valores_post <- as.numeric(names(tab_post_sim))

prob_post_bn <- b_post / (b_post + m)
dens_post_bn <- dnbinom(valores_post, size = a_post, prob = prob_post_bn)

plot(
  valores_post, tab_post_sim,
  type = "h", lwd = 2, col = "darkorange",
  xlim = c(0, max(valores_post)),
  ylim = c(0, max(tab_post_sim, dens_post_bn) * 1.1),
  xlab = expression(y[novo]), ylab = "P(Y_novo = y)",
  main = "Preditiva a posteriori: simulacao x BinNeg(a+y, (b+n)/(b+n+m))"
)
points(valores_post + 0.3, dens_post_bn, type = "h", lwd = 2, col = "steelblue")
legend(
  "topright",
  legend = c("Simulacao Monte Carlo", "BinNeg(size=a+y, prob=(b+n)/(b+n+m))"),
  col = c("darkorange", "steelblue"), lwd = 2, bty = "n"
)


# =====================================================================



x11()
par(mfrow = c(2, 2))

#  priori vs posteriori de lambda
curve(dgamma(x, a, b), from = 0, to = 10, lwd = 2, col = "red",
      ylim = c(0, max(dgamma(seq(0,10,0.01), a_post, b_post))*1.7),
      xlab = expression(lambda), ylab = "densidade", main = "lambda: priori x posteriori")
curve(dgamma(x, a_post, b_post), add = TRUE, lwd = 2, col = "black")
legend("topright", legend = c("Priori", "Posteriori"), col = c("red","black"), lwd = 2, bty = "n")

# Y observado (dado) marcado sobre a preditiva a priori
plot(valores_priori, tab_priori_sim, type = "h", lwd = 2, col = "darkorange",
     xlab = "y", ylab = "P(Y=y)", main = "Preditiva a priori")
#bline(v = y, col = "black", lty = 2, lwd = 2)
#legend("topright", legend = "y observado", lty = 2, lwd = 2, bty = "n")

# priori x posteriori preditivas juntas
plot(valores_priori, tab_priori_sim, type = "h", lwd = 2, col = "darkorange",
     xlim = c(0, max(valores_priori, valores_post)),
     ylim = c(0, max(tab_priori_sim, tab_post_sim)*1.1),
     xlab = "y", ylab = "probabilidade", main = "Preditivas: priori x posteriori")
lines(valores_post + 0.3, tab_post_sim, type = "h", lwd = 2, col = "steelblue")
legend("topright", legend = c("Preditiva a priori", "Preditiva a posteriori"),
       col = c("darkorange", "steelblue"), lwd = 2, bty = "n")

# preditiva a posteriori
plot(valores_post, tab_post_sim, type = "h", lwd = 2, col = "steelblue",
     xlab = expression(y[novo]), xlim = c(0,20), ylab = "P(Y_novo=y)", main = "Preditiva a posteriori")


par(mfrow = c(1, 1))


# =====================================================================

# Comparando as preditivas a priori por simulação e pelo resultado analítico

a2 <- 86
b2 <- 5
N2 <- 1e6   

lambda2 <- rgamma(N2, a2, b2)
y2 <- rpois(N2, lambda = lambda2)

vy2 <- prop.table(table(y2))
valores_y2 <- as.numeric(names(vy2))

p2 <- b2 / (1 + b2)
ybn2 <- rnbinom(N2, size = a2, prob = p2)
vybn2 <- prop.table(table(ybn2))

plot(
  valores_y2, vy2,
  type = "h", col = "darkorange", lwd = 2,
  xlab = "y", ylab = "P(Y=y)",
  main = "Mistura Poisson-Gama (simulada) x Binomial Negativa (rnbinom)",
  ylim = c(0, max(vy2, vybn2) * 1.1)
)
lines(
  as.numeric(names(vybn2)) + 0.3, vybn2,
  type = "h", col = "steelblue", lwd = 2
)
legend("topright", legend = c("Mistura Poisson-Gama simulada", "rnbinom(size=a, prob=b/(b+1))"),
       col = c("darkorange", "steelblue"), lwd = 2, bty = "n")