##rm(list=ls())
##
## Modelo binomial: inferência para o parâmetro de proporção.
##
## Este tutorial mostra aproximação computacional por simulação
## para resultados analíticos da análise Bayesiana
##
## Associe os comandos a seguir com os resultados teóricos para este modelo.
## Faça interpretações no contexto do exemplo das aulas (canhotos) ou similar 

##
## PARTE 1: O modelo (mecanismo gerador de dados)
##
## Assume-se que há incerteza sobre o valor do parâmetro.
## A incerteza é expressa por uma distribuição de probabilidades para o parâmetro.

a <- 3; b <- 15
curve(dbeta(x, a, b),
      xlab = expression(theta),
      ylab = expression(group("[", theta, "]"))
      )

## Assuma-se que, para um valor do parâmetro, as observações vem de um modelo binomial,
## ou seja, com probabilidade para cada possível valor dada pelo modelo Binomial.
(th <- rbeta(1, a, b))
n <- 27
plot(0:27, dbinom(0:27, size = n, prob = th), type = "h",
     xlab = "y", ylab = "[Y]")
## E na prática uma amostra fornece um dos possíveis valores
(yobs <- rbinom(1, size = n, prob = th))
points(yobs, 0, cex = 2, col = 2, pch = 19)

## repetir algumas vezes o comando acima. Interprete.

## mas na realidade não sabemos th !
rm(th)

##
## PARTE 2: O suposto mecanismo gerador de dados: aproximação por simulação
##
Nsim <- 10E06                
th <- rbeta(Nsim, a, b)                ## simulações para valor de theta (priori)
y <- rbinom(Nsim, size = n, prob = th) ## respectivos valores simulados para y (preditiva priori)

## Gráficos das priori e preditiva priori: simuladas e teóricas
par(mfrow = c(1,2), mar = c(3.5,3.5,0,0), mgp = c(2,1,0))
hist(th, nclass = 30, freq = FALSE, col = "transparent", main = "",
     xlab = expression(theta), ylab = expression(group("[", theta, "]")))
curve(dbeta(x, a , b), add = TRUE, col = 4)
legend("topright", c("Simulada","Teorica"), lty = 1, col = c(1,4))

(pred.priori <- prop.table(table(factor(y, levels = 0:n))))
plot(0:n, pred.priori, type = "h",
     xlab = "y", ylab = expression(group("[", Y, "]")))
dbetabin <- function(x, size, shape1, shape2, log = FALSE){
    res <- lchoose(size, x) + lbeta(x+shape1, size-x+shape2) - lbeta(shape1,shape2)
    if(!log) res <- exp(res)
    res
}
dy.teo <- dbetabin(0:n, n, a, b)
lines((0:n)+0.1, dy.teo, type = "h", col = 4)
legend("topright", c("Simulada","Teorica"), lty = 1, col = c(1,4))

## Conclusão: podemos aproximar (bem) as distribuições teóricas
##            por simulações 

## 
## PARTE 3: análise com dado obtido de uma amostra
##
## - Posteriori
## O resultado obtido em uma amostra de tamanho "n":
yobs = 5

##  [\theta|y]: condicionando no dado observado (obtenção da posteriori)
head(cbind(th, y), n = 15)
th.y <- th[y == yobs]
length(th)         
length(th.y)         
length(th.y)/Nsim

hist(th.y, nc = 30, freq = FALSE, col = "transparent", xlim = c(0,1))

## comparando [\theta] e [\theta|y]
par(mfrow=c(1,2))
hist(th, nc = 30, freq = FALSE, col = "transparent", xlim = c(0,1))
hist(th.y, nc = 30, freq = FALSE, col = "transparent", xlim = c(0,1))

##
## Posteriori: Aproximação por simulação vs resultado analítico 
par(mfrow = c(1,1))
hist(th.y, nc = 30, freq = FALSE, col = "transparent")
curve(dbeta(x, shape1 = yobs+a, shape2 = n-yobs+b), col = 4, lwd = 2, add = TRUE)
## forçando xlim em (0,1)
hist(th.y, nc = 30, freq = FALSE, col = "transparent", xlim = c(0,1))
curve(dbeta(x, shape1 = yobs+a, shape2 = n-yobs+b), col = 4, lwd = 2, add = TRUE)

## A estimativa de MV 
(th.hat <- yobs/n)

## Preditiva (a posteriori, ou marginal) [Y_p|y]
np <- 40
y.y <- rbinom(length(th.y), size = np, prob = th.y)
round(prop.table(table(y.y)), dig = 5)
round(prop.table(table(factor(y.y, levels = 0:np))), dig = 5)
dy.y <- prop.table(table(factor(y.y, levels = 0:np)))
plot(0:np, dy.y, type = "h")

## Preditiva "plug-in" (não Bayesiana) 
y.ythHat <- dbinom(0:np, size = np, prob = th.hat)
plot(0:np, y.ythHat, type = "h")
##
## comparando
plot((0:np)+0.1, y.ythHat, type = "h", col = "darkorange")
lines((0:np)-0.1, dy.y, type = "h")
legend("topright", c("Bayesiana","plugin"), lty = 1,
       col = c("black","darkorange"))

## Teórica
dbetabin <- function(x, size, shape1, shape2, log = FALSE){
    res <- lchoose(size, x) + lbeta(x+shape1, size-x+shape2) - lbeta(shape1,shape2)
    if(!log) res <- exp(res)
    res
}
dy.y.teo <- dbetabin(0:np, np, yobs+a, n-yobs+b)
#dy.y.teo1 <- VGAM::dbetabinom.ab(0:np, np, yobs+a, n-yobs+b)
#all.equal(dy.y.teo, dy.y.teo1) 
plot((0:np)-0.1, dy.y, type = "h")
lines((0:np)+0.1, dy.y.teo, type = "h", col = 4)
legend("topright", c("Simulada","Teorica"), lty = 1,
       col = c(1,4))

## Preditivas para outro tamanho de amostra
np <- 120
y120.y <- rbinom(length(th.y), size = np, prob = th.y)
dy120.y <- prop.table(table(factor(y120.y, levels = 0:np)))
plot(0:np, dy120.y, type = "h")
y120.ythHat <- dbinom(0:np, size = np, prob = th.hat)
plot(0:np, y120.ythHat, type = "h")
##
plot((0:np)+0.1, y120.ythHat, type = "h", col = 2)
lines((0:np)-0.1, dy120.y, type = "h")
##
## 
dy120.y.teo <- dbetabin(0:np, np, a+yobs, n-yobs+b)
plot((0:np)-0.1, dy120.y, type = "h")
lines((0:np)+0.1, dy120.y.teo, type = "h", col = 4)

plot((0:np)-0.1, dy120.y, type = "h", xlim = c(0,60))
lines((0:np)+0.1, dy120.y.teo, type = "h", col = 4, xlim = c(0,60))


## Resumos, probabilidades,  etc
## Resumos da posteriori
## aproximação por simulação e
## exatos (analítico) por propriedades da distribuição Beta
(th.y.media <- mean(th.y))
(yobs+a)/(yobs+a + n-yobs+b)

(th.y.mediana <- median(th.y))
(a+yobs - 1/3)/(yobs+a + n-yobs+b - 2/3)

dth.y <- density(th.y)
plot(dth.y)
(dth.y.moda <- dth.y$x[which.max(dth.y$y)])
(yobs+a-1)/(yobs+a-1 + n-yobs+b-1)

hist(th.y, nc = 30, freq = FALSE, col = "transparent")
curve(dbeta(x, shape1 = yobs+a, shape2 = n-yobs+b), col = 4, lwd = 2, add = TRUE)
abline(v = c(th.y.media, th.y.mediana, dth.y.moda), lty = 1:3, col = "blue") 
legend("topright", c("moda","mediana","média"), lty = 1:3)

mean(th.y < 0.08)
mean(th.y > 0.15)
mean(th.y > 0.10 & th.y < 0.25)
median(th.y)
quantile(th.y)
quantile(th.y, probs = c(0.025, 0.975))
##
pbeta(0.08, a+yobs, b+n-yobs)
pbeta(0.15, a+yobs, b+n-yobs, lower = FALSE)
diff(pbeta(c(0.10, 0.25), a+yobs, b+n-yobs))
qbeta(0.50, a+yobs, b+n-yobs)
qbeta(c(0,0.25, 0.50, 0.75,1), a+yobs, b+n-yobs)
qbeta(c(0.025, 0.975), a+yobs, b+n-yobs)
##
th.y2 <- th.y^2
mean(th.y2 < 0.05)
mean(th.y2 > 0.01)

th.ylg <- log(th.y/(1-th.y))
mean(th.ylg < -2)

par(mfrow=c(1,3), mar=c(3,3,1,0), mgp = c(2,1,0))
hist(th.y)
hist(th.y2)
hist(th.yl)

##
##
tb <- table(y120.y)
as.numeric(names(tb)[which.max(tb)])

median(y120.y)
quantile(y120.y, prob = c(0.025, 0.975))

median(y120.y)
quantile(y120.y, prob = c(0.025, 0.975))

mean(y120.y)
120 * (yobs+a)/(yobs+a+n-yobs+b)

mean(y120.y <= 10)
mean(y120.y >= 40)
mean(y120.y > 20 & y120.y <= 30)

VGAM::pbetabinom.ab(10, 120, yobs+a, n-yobs+b)
1 - VGAM::pbetabinom.ab(39, 120, yobs+a, n-yobs+b)
diff(VGAM::pbetabinom.ab(c(20, 30), 120, yobs+a, n-yobs+b))