#
# Inferencia para theta = P(cruz) con una distribución a priori mixtura:
#   f(theta) = w_1 f(theta|a_1,b_1) + ... + w_k f(theta|a_k,b_k)
# donde w_i > 0 y w_1 + ... + w_k = 1 son pesos y f(theta|a,b) es una densidad beta.
#
rm(list=ls())
#
# Distribución a priori
#
k <- 2
w <- c(0.5,0.5)
a <- c(5,6)
b <- c(5,3)
#
# Generación de datos
#
tiradas <- 20
theta <- 0.3
cruces <- rbinom(1,tiradas,theta)
caras <- tiradas-cruces
#
# Distribución a posteriori
#
apost <- a+cruces
bpost <- b+caras
wpost <- w*beta(apost,bpost)/beta(a,b)
wpost <- wpost/sum(wpost)
#
# Gráfico de la distribución a priori (rojo), verosimilitud escalada (azul) y a posteriori (verde).  
#
thetagrid <- c(0:1000)/1000
verosimilitud <- dbeta(thetagrid,cruces+1,caras+1)
fthetapri <- rep(0,1001)
fthetapost <- rep(0,1001)
for (i in 1:k){
  fthetapri <- fthetapri+w[i]*dbeta(thetagrid,a[i],b[i])
  fthetapost <- fthetapost+wpost[i]*dbeta(thetagrid,apost[i],bpost[i])
}
maxf <- max(c(verosimilitud,fthetapri,fthetapost))
plot(thetagrid,fthetapri,type='l',lwd=2,col='red',xlab=expression(theta),ylab='f',ylim=c(0,maxf))
lines(thetagrid,fthetapost,type='l',lwd=2,col='green')
lines(thetagrid,verosimilitud,type='l',lwd=2,col='blue')
#
# Predicción
#
if (!require("extraDistr")) install.packages("extraDistr")
library(extraDistr)
nuevastiradas <- 10
pxpred <- rep(0,(nuevastiradas+1))
for (i in 1:k){
  pxpred <- pxpred+wpost[i]*dbbinom(x=c(0:nuevastiradas),size=nuevastiradas,alpha=apost[i],beta=bpost[i],log=FALSE)
}
plot(c(0:nuevastiradas),pxpred,type='h',lwd=2,col="green",xlab="x",ylab="P")