Introducción

En las sesiones previas, en problemas de estimar \(\theta = P(\textrm{cruz})\) en problemas de tiradas de monedas, siempre hemos utilizado una distribución a priori (beta) conjugada. Supongamos ahora que empleamos otra distribución a priori, por ejemplo una distribución normal truncada para tener soporte entre 0 y 1.

rm(list=ls())
mu <- 0.5
sigma2 <- 0.09
sigma <- sqrt(sigma2)

Ya, al observar una muestra de \(y\) cruces y \(n-y\) caras , la distribución a posteriori es \[ f(\theta|\textrm{datos}) \propto \exp\left(-\frac{1}{2\sigma^2}(\theta-\mu)^2 \right) \theta^y (1-\theta)^{n-y} \] para \(0 < \theta < 1\) que no es fácil de integrar de manera analítica.

tiradas <- 20
theta <- 0.25
cruces <- rbinom(1,tiradas,theta)
caras <- tiradas-cruces

Integración numérica

Como es un problema simple, podemos utilizar la integración numérica (por ejemplo la regla de Simpson) para calcular la constante de la distribución.

fprop <- function(theta,mu,sigma,cruces,caras){
  f <- dnorm(theta,mu,sigma)*(theta^cruces)*((1-theta)^caras)
  return(f)
}
c <- integrate(fprop,lower=0,upper=1,mu=mu,sigma=sigma,cruces=cruces,caras=caras)$value
c
## [1] 3.162158e-05

Dividiendo por la constante, podemos calcular y graficar la densidad (y la a priori y verosimilitud escalada)

fdens <- function(theta,mu,sigma,cruces,caras,c){
  f <- fprop(theta,mu,sigma,cruces,caras)/c
  return(f)
}
rejilla <- c(0:1000)/1000
ftheta <- fdens(theta=rejilla,mu,sigma,cruces,caras,c)
plot(rejilla,ftheta,type='l',lwd=2,col='green',xlab=expression(theta),ylab='f')
lines(rejilla,dbeta(rejilla,cruces+1,caras+1),type='l',lwd=2,col='blue')
lines(rejilla,dnorm(rejilla,mean=mu,sd=sigma)/(pnorm(1,mean=mu,sd=sigma)-pnorm(0,mean=mu,sd=sigma)),type='l',lwd=2,
      col='red')

No obstante, si queremos calcular la media, varianza, … cada vez tenemos que hacer otra integración numérica.

fmedia <- function(theta,mu,sigma,cruces,caras,c){
  f <- theta*fdens(theta,mu,sigma,cruces,caras,c)
  return(f)
}
media <- integrate(fmedia,lower=0,upper=1,mu=mu,sigma=sigma,cruces=cruces,caras=caras,c=c)$value
media
## [1] 0.2030858

Muestreo de importancia

Implementar el método de Monte Carlo básico tampoco es fácil. ¿Cómo podemos muestrear la distribución a posteriori sin saber ni su constante de integración?

Supogamos que sí podemos muestrear de otra distribución, \(g\).

Luego, observamos que \[ \int f(\theta|\textrm{datos}) \, d\theta = \int \frac{f(\theta|\textrm{datos})}{g(\theta)} g(\theta) \, d\theta = 1. \]

Esta fórmula nos proporciona un método para estimar la constante de la distribución. Sea \(f_{\propto}\) la densidad a posteriori sin constante. Luego, si tomamos una muestra de tamaño \(N\) de \(g\), podemos estimar \(c\) con \[ \hat{c} = \frac{1}{N} \sum_{i=1}^N w_i \] donde \(w_i = \frac{f_{\propto}(\theta_i)}{g(\theta_i)}\) y \(\theta_1,...,\theta_N\) son los valores muestreados.

En nuestro ejemplo, supongamos que \(g\) es una distribución uniforme.

N <- 100000
theta <- runif(N)
w <- fprop(theta,mu,sigma,cruces,caras)/1
cest <- mean(w)
c(c,cest)
## [1] 3.162158e-05 3.167237e-05

Luego, recordamos que la media es \[ E[\theta|\textrm{datos}] = \int \theta f(\theta|\textrm{datos}) \, d\theta = \int \frac{\theta f(\theta|\textrm{datos})}{g(\theta)} g(\theta) \, d\theta \] y podemos estimarla con \[ \hat{E}[\theta|\textrm{datos}] = \frac{1}{N\hat{c}} \sum_{i=1}^{n} w_i \theta_i. \] Es posible estimar los otros momentos de manera parecida.

El bootstrap bayesiano

La muestra de importancia es una muestra de \(g\) en lugar de la verdadera distribución a posteriori. Podemos generar una muestra de la a posteriori mediante el submuestreo: remuestreando los valores \(\theta_1,...,\theta_N\) con reemplazamiento y con probabilidades proporcionales a sus pesos, \(w_1,...,w_N\).

N2=N
thetapost <- sample(theta,size=N2,replace=TRUE,prob=w)
hist(thetapost,probability=TRUE,xlab=expression(theta),ylab='f',main='')
lines(rejilla,ftheta,type='l',lwd=2,col='green')

Problemas del muestreo de importancia

Es fundamental que la distribución \(g\) no está muy alejado y tiene colas más largas que \(f \). Si no es el caso, la mayoría de los pesos van a ser muy pequeños.

N <- 100000
theta <- rbeta(N,10,1)
w <- fprop(theta,mu,sigma,cruces,caras)/dbeta(theta,N,1)
cest <- mean(w)
c(c,cest)
## [1] 3.162158e-05          Inf

En nuestro ejemplo, \(g(\theta)\) es mucho más pequeño que \(f(\theta) \) en la cola izquierda de la distribución.

El método de rechazo

Una alternativa al muestreo de importancia es el algoritmo de rechazo. En este caso, se muestra de una distribución \(g\) con colas más largas que \(f\) como anteriormente pero aqúi se aceptan algunos de los valores muestrados y rechazan otros.

El algoritmo es el siguiente

  1. Definir una constante \(M \) tal que \(f(\theta|\text{datos}) \le M g(\theta) \, \forall \theta . \)
  2. Generar valores \(\theta_i\) de \(g\) y aceptar \(\theta_i\) con probabilidad \(p = \frac{f(\theta_i|\text{datos})}{Mg(\theta_i)}.\)

Observamos que a pesar de no saber la constante de la distribución, el algoritmo funcionará, ya que \[ f(\theta|\text{datos}) \le M g(\theta) \Rightarrow f_{\propto} (\theta|\text{datos}) \le c M g(\theta) \] y escribiendo \(M^* = cM\), aceptamos un valor generado con probabilidad \[ p = \frac{f_{\propto}(\theta_i|\textrm{datos})}{M^* g(\theta_i)}. \]

En nuestro ejemplo, sea \(g\) una densidad uniforme como anteriomente. Luego necesitamos elegir \(M^*\) tal que \[ f_{\propto}(\theta|\text{datos}) = \exp\left(-\frac{1}{2\sigma^2}(\theta-\mu)^2 \right) \theta^y (1-\theta)^{n-y} \le M^* \times 1.\]

Mstar <- max(fprop(theta=rejilla,mu,sigma,cruces,caras))
Mstar <- Mstar*1.01
N <- 100000
theta <- runif(N)
pac <- fprop(theta,mu,sigma,cruces,caras)/(Mstar*1)
u <- runif(N)
thetaac <- theta[u<pac]
Nac <- length(thetaac)
Nac
## [1] 20792

Observamos que el número de puntos aceptados es aproximadamente un cuarto de los puntos generados. El muestreador no es muy eficiente. No obstante, el histograma es bastante cercana a la densidad calculada mediante integración numérica y la media estimada también es muy cercana a la ``verdadera’’ media.

hist(thetaac,probability=TRUE,xlab=expression(theta),ylab='f',main='')
lines(rejilla,ftheta,type='l',lwd=2,col='green')

c(media,mean(thetaac))
## [1] 0.2030858 0.2031903