Introducción

En esta sesión, ilustramos como hacer inferencia clasica y inferencia bayesiana para \(\theta = P(\text{cruz})\) en problemas asociadas con tiradas de monedas. Empezamos generando unos datos. Como alternative podríamos haber generado datos de una distribución binomial negativa o una distribución geométrica. Lo importante es tener definido el número de tiradas, el número de caras y el número de cruces.

rm(list=ls())
tiradas <- 20
theta <- 0.25
cruces <- rbinom(1,tiradas,theta)
caras <- tiradas-cruces

Inferencia clásica

El primer paso en la inferencia clásica es calcular el EMV, \(\hat{\theta} = \frac{\text{cruces}}{\text{tiradas}} \).

emvtheta <- cruces/tiradas
emvtheta
## [1] 0.15

En segundo lugar, podemos calcular un intervalo de confianza, \(\hat{\theta} \pm z_{1-\alpha/2} \sqrt{\frac{\hat{\theta}(1-\hat{\theta})}{n}} \). ¿Sería correcto usar este intervalo si hubieramos generado los datos de una binomial negativa?

ci <- 0.95 
c <- qnorm((1+ci)/2,0,sqrt(emvtheta*(1-emvtheta)/tiradas))
intervalo <- c(emvtheta-c,emvtheta+c)
intervalo
## [1] -0.006490575  0.306490575

Luego, podemos contrastar, por ejemplo \(H_0: \theta=\theta_0\) frente a \(H_1:\theta < \theta_0\). Recordamos que bajo \(H_0\), \(\hat{\theta} \sim N\left(\theta_0, \frac{\theta_0(1-\theta_0)}{n} \right) \). En este caso, fijamos \(\theta_0=0.5\) para contrastar si la moneda es sesgada a favor de caras.

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
## [1] 0.0008725593

¿El contraste sería correcto con datos de la binomial negativa?

Finalmente, quizás queremos hacer predicción. Supongamos que queremos predecir el número de cruces, \(X\), en 10 tiradas más de la moneda. Luego, la predicción clásica es \(X \sim\) Binomial\(\left(10,\hat{\theta}\right)\).

nuevastiradas <- 10
pxpredclas <- dbinom(c(0:nuevastiradas),nuevastiradas,emvtheta)

Inferencia bayesiana

El primer paso para implementar la inferencia bayesiana es definir una distribución a priori apropiada para \(\theta\). Podemos seleccionar por ejemplo una distribución uniforme, \(f(\theta) = 1\) para \(0 < \theta < 1\).

¿Es una distribución beta?

a <- 1
b <- 1

Después, calculamos la distribución a posteriori. Dada una distribución a priori Beta\((a,b)\), la distribución a posteriori es otra distribución beta, Beta\((a+\text{cruces},b+\text{caras})\).

apost <- a+cruces
bpost <- b+caras

La media, mediana o moda a posteriori pueden ser buenas estimadores puntuales de \(\theta\).

thetamedia <- apost/(apost+bpost)
thetamoda <- (apost-1)/(apost+bpost-2)
thetamediana <- qbeta(0.5,apost,bpost)
c(thetamedia,thetamoda,thetamediana)
## [1] 0.1818182 0.1500000 0.1720895

También es fácil calcular un intervalo de credibilidad.

intervalobayes <- c(qbeta(alpha/2,apost,bpost),qbeta(1-alpha/2,apost,bpost))
intervalobayes
## [1] 0.05446357 0.36342399

Podemos calcular la probabilidad a posteriori de \(H_0\).

if (H1=="menor"){
  pH0 <- 1-pbeta(theta0,apost,bpost)
}else if (H1=="mayor"){
  pH0 <- pbeta(theta0,apost,bpost)
}
pH0
## [1] 0.0007448196

¿Es similar al p-valor clásico? ¿La interpretación es la misma?

Recordamos que la distribución predictiva en este problema es la distribución beta-binomial. Como no es una distribución estandár en R, tenemos que instalar un paquete de distribuciones.

require(extraDistr)
## Loading required package: extraDistr
## Warning: package 'extraDistr' was built under R version 3.5.3
library(extraDistr)

Luego, podemos hacer la predicción basada en la distribución a priori \[ f(x) = \int_0^1 f(x|\theta) f(\theta) \, d\theta \] y la predicción basada en la a posteriori \[ f(x|\text{datos}) = \int_0^1 f(x|\theta) f(\theta|\text{datos}) \, d\theta . \]

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)

Ahora podemos comparar la inferencia bayesiana con la clásica de manera visual. Primero graficamos la verosimilitud escalada (en azul), la distribución a priori (en rojo) y la distribución a posteriori (en 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')

¿Dónde está la verosimilitud escalada?

Por último, podemos comparar la distribución predictiva clásica (en azul) con la predicción basada en la a priori (en rojo) y la predicción bayesiana basada en la a posteriori (en 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")