#
# Inferencia para el parámetro p de una distribución binomial(n,theta).
#
rm(list=ls())
#
# Generar datos
#
tiradas <- 20
theta <- 0.3
cruces <- rbinom(1,tiradas,theta)
caras <- tiradas-cruces
#
# INFERENCIA CLÁSICA
# a) EMV
#
emvtheta <- cruces/tiradas
emvtheta 
#
# b) Intervalo de confianza
#
ci <- 0.95 
c <- qnorm((1+ci)/2,0,sqrt(emvtheta*(1-emvtheta)/tiradas))
intervalo <- c(emvtheta-c,emvtheta+c)
intervalo
#
# c) Contraste y p-valor
#
alpha <- 0.05
theta0 <- 0.5
H1 <- "menor"
if (H1=="menor"){
  pvalor <- pnorm(emvtheta,theta0,sqrt(theta0*(1-theta0)/tiradas))
}else if (H1=="mayor"){
  pvalor <- 1-pnorm(emvtheta,theta0,sqrt(theta0*(1-theta0)/tiradas))
}
pvalor
#
# d) Predicción
#
nuevastiradas <- 10
pxpredclas <- dbinom(c(0:nuevastiradas),nuevastiradas,emvtheta)
#
# INFERENCIA BAYESIANA
# a) A prioris
#
a <- 1
b <- 1
#
# b) A posterioris
#
apost <- a+cruces
bpost <- b+caras
#
# c) Media, moda y mediana a posteriori
#
thetamedia <- apost/(apost+bpost)
thetamoda <- (apost-1)/(apost+bpost-2)
thetamediana <- qbeta(0.5,apost,bpost)
c(thetamedia,thetamoda,thetamediana)
#
# d) Intervalo de credibilidad
#
intervalobayes <- c(qbeta(alpha/2,apost,bpost),qbeta(1-alpha/2,apost,bpost))
intervalobayes
#
# e) Probabilidad a posteriori de H0.
#
if (H1=="menor"){
  pH0 <- 1-pbeta(theta0,apost,bpost)
}else if (H1=="mayor"){
  pH0 <- pbeta(theta0,apost,bpost)
}
#
# f) Predicción: el paquete "extraDistr"" contiene la distribución beta-binomial.
#
if(!require("extraDistr")) install.packages("extraDistr")
library(extraDistr)
pxpredpost <- dbbinom(x=c(0:nuevastiradas),size=nuevastiradas,alpha=apost,beta=bpost,log=FALSE)
pxpredpri <- dbbinom(x=c(0:nuevastiradas),size=nuevastiradas,alpha=a,beta=b,log=FALSE)
#
# Gráfico de la distribución a priori (rojo), la verosimilitud escalada (azul) y la a posteriori (verde).
#
thetagrid <- seq(0,1,0.001)
vero <- dbeta(thetagrid,cruces+1,caras+1)
priori <- dbeta(thetagrid,a,b)
posteriori <- dbeta(thetagrid,apost,bpost)
fmax <- max(vero,priori,posteriori)
plot(thetagrid,vero,type='l',lwd=2,col='blue',ylim=c(0,fmax),xlab=expression(theta),ylab='f')
lines(thetagrid,priori,type='l',lwd=2,col="red")
lines(thetagrid,posteriori,type='l',lwd=2,col='green')
#
# Gráfico comparando las distribuciones predictivas: frecuentista (azul) bayesiana a priori (rojo) y a posteriori (verde)
# 
pxmax <- max(pxpredpri,pxpredpost,pxpredclas)
plot(c(0:nuevastiradas),pxpredpost,type='h',lwd=2,col='green',xlab='x',ylab='p',ylim=c(0,pxmax))
lines(c(0:nuevastiradas)-0.1,pxpredpri,type='h',lwd=2,col="red")
lines(c(0:nuevastiradas)+0.1,pxpredclas,type='h',lwd=2,col="blue")