Modelo de Bernoulli para muestras pequeñas
Alfredo Zea Valenzuela
09/12/2022
1. Modelo de Bernoulli clásico
Un modelo Bernoulli clasico se puede definir de la siguiente manera: \[ y_{i} \sim Bernoulli(p)~~p\in [0,1]\]
Procedemos a simular una muestra pequeña de 25 observaciones con parámetro \(p\) conocido:
# Eliminar los objetos de Global Environment
rm(list = ls())
options("digits" = 3)
set.seed(123)
# ==============================================================================
# Simulando una muestra Bernoulli
# ==============================================================================
N = 25
p = 0.65
y = rbinom ( size = 1, prob = p , n = N )
hist(y,main = sprintf("Histograma distribución Bernoulli con p=%.2f",p) ,
sub = paste("Tamaño de muestra:", N) ,
ylab = "Frecuencia",
xlab = "y", prob = TRUE)
# ==============================================================================
# Estimacion del párametros por MV
# ==============================================================================
pmv = round(mean(y),3)
df_resumen =
data.frame( Parametro = c('p','α','β'),
Teorico = c(p,'',''),
MediaMV = c(pmv,'',''))
df_resumen## Parametro Teorico MediaMV
## 1 p 0.65 0.52
## 2 α
## 3 βComo se aprecia se puede estimar el parámetro \(p\) con el estimador clásico de máxima verosimilitud.
2. Modelo de Bernoulli bayesiano
En ocasiones se requiere estimar el parámetro del modelo Bernoulli mediante el enfoque bayesiano y considerarademás una muestra pequeña. En el enfoque bayesiano el parámetro de la distribución se considera una variable aletoria y por consiguiente tiene asociada una distribución.
Un modelo Bernoulli aplicando el enfoque bayesiano se puede extender forma jerárquica utilizando distribuciones a priori: \[ y_{i} \sim Bernoulli(p)~~con~~p\in [0,1]\] \[ p \sim Beta(\alpha,\beta)~~con~~\alpha>0~y~\beta>0\] \[\alpha \sim Gamma(a,b)~~con~~a>0~y~b>0\] \[\beta \sim Gamma(c,d)~~con~~c>0~y~d>0\] Siendo los parámetros de este modelo \(p\),\(\alpha\) y \(\beta\). Mientras que \(a\),\(b\),\(c\) y \(d\) son los denominados Hiperparámetros y son valores conocidos.
3. Distribución A posteriori y distribuciones condicionales
Aplicando el teorema de bayes se puede deducir que la distribución A Posteriori de los parámetros dado la muestra es:
\[\pi(\theta|Y) \propto (p^{\sum_{i=1}^{N}y_{i}}(1-p)^{N-\sum_{i=1}^{N}y_{i}}p^{\alpha-1}(1- p)^{\beta-1})(\alpha^{a-1}e^{-b\alpha})(\beta^{c-1}e^{-d\beta})\] A partir de esta distribución y debido que se han utilizado distribuciones a priori conjugadas en la extensión del modelo, se puede deducir la distribución condicional para cada uno de los parámetros.
Parámetro: \(p\) \[p|\alpha,\beta,Y \sim Beta(\sum_{i=1}^{N}y_{i} + \alpha, N - \sum_{i=1}^{N}y_{i}+ \beta)\] Parámetro: \(\alpha\) \[\alpha|p,\beta,Y \sim Gamma(a, b - ln(p))\] Parámetro: \(\beta\) \[\beta|p,\alpha,Y \sim Gamma(c, d - ln(1- p))\]
Con estas tres distribuciones condicionales aplicamos el algoritmo de Gibbs para obtención de muestra aleatoria de los parámetros.
4. Aplicación del algoritmo de Gibbs
# ------------------------------------------------------------------------------
# Algoritmo de Gibs
# ------------------------------------------------------------------------------
fMuestradorGibbs = function( y,m,k,Oi,Hi){
a = Hi[1]
b = Hi[2]
c = Hi[3]
d = Hi[4]
Os = matrix(0.00, nrow = m , ncol = length(Oi))
Os[1,] = Oi
suma_y = sum(y)
N = length(y)
for (i in 2:m){
Os[i,1] = rbeta(n = 1 ,suma_y + Os[(i-1),2],N - suma_y + Os[(i-1),3] )
Os[i,2] = rgamma(n = 1 ,a , b - log(Os[i,1]) )
Os[i,3] = rgamma(n = 1 ,c , d - log(1 - Os[i,1]) )
}
par(mfrow=c(3,1))
plot( Os[,1],type="l",xlab = sprintf("Iteraciones: %.0f",m), ylab ="p",
main = "Convergencia de parámetro: p")
plot( Os[,2],type="l",xlab = sprintf("Iteraciones: %.0f",m), ylab ="α",
main = "Convergencia de parámetro: α")
plot( Os[,3],type="l",xlab = sprintf("Iteraciones: %.0f",m), ylab ="β",
main = "Convergencia de parámetro: β")
par(mfrow=c(1,1))
return(Os[(m - k + 1):m,])
}Procedemos a ejecutar el algoritmo para obtener una muestra de tamaño k mediante una simulación de m iteraciones.
# ==============================================================================
# Simulaciones
# ==============================================================================
# Número de iteraciones
m = 100000
# Tamaño de la muestra
k = 1000
# Vector de parametros inicial
Oi = c(0.1,0.1,0.1)
# Vector de Hiperparámetros
Hi = c(0.1,0.1,0.1,0.1)
Os = fMuestradorGibbs(y,m,k,Oi,Hi)
tail(Os,50)## [,1] [,2] [,3]
## [951,] 0.462 1.12e-01 8.25e-07
## [952,] 0.500 1.25e-06 1.62e-09
## [953,] 0.514 2.20e-13 1.55e-02
## [954,] 0.537 1.62e-01 2.12e-01
## [955,] 0.486 1.00e-08 1.16e-03
## [956,] 0.487 3.13e-02 8.77e-07
## [957,] 0.410 4.36e-01 3.28e-02
## [958,] 0.559 1.85e+00 1.87e-02
## [959,] 0.679 7.55e-06 1.61e-21
## [960,] 0.323 1.79e-03 6.31e-09
## [961,] 0.498 3.51e-02 3.86e-01
## [962,] 0.381 4.83e-08 4.09e-13
## [963,] 0.510 3.47e-03 3.05e-01
## [964,] 0.605 2.26e-03 2.65e-12
## [965,] 0.589 9.12e-01 1.82e-05
## [966,] 0.611 4.85e-02 2.14e-01
## [967,] 0.563 3.07e-09 1.55e-12
## [968,] 0.464 4.63e-08 6.79e-04
## [969,] 0.626 3.11e-01 2.11e-08
## [970,] 0.516 2.83e-05 2.88e-16
## [971,] 0.650 8.53e-01 8.35e-06
## [972,] 0.717 7.30e-09 1.68e-03
## [973,] 0.569 8.43e-01 2.28e-06
## [974,] 0.408 5.82e-11 9.43e-01
## [975,] 0.396 9.03e-02 2.89e-14
## [976,] 0.370 1.67e-02 1.49e-01
## [977,] 0.400 8.84e-04 4.72e-03
## [978,] 0.626 1.06e-01 5.23e-05
## [979,] 0.673 5.05e-05 1.03e-04
## [980,] 0.491 6.37e-02 7.55e-02
## [981,] 0.754 1.49e-17 6.87e-03
## [982,] 0.371 3.69e-03 2.73e-01
## [983,] 0.577 5.67e-05 3.78e-01
## [984,] 0.478 2.80e-02 1.26e-03
## [985,] 0.529 7.12e-05 3.42e-01
## [986,] 0.550 2.14e-10 4.42e-02
## [987,] 0.485 2.93e-11 1.13e-01
## [988,] 0.543 1.70e-10 1.38e-08
## [989,] 0.683 1.69e+00 2.86e-03
## [990,] 0.630 1.95e+00 4.96e-08
## [991,] 0.497 4.57e-01 7.25e-03
## [992,] 0.553 2.08e-04 5.29e-02
## [993,] 0.325 2.53e-01 9.83e-01
## [994,] 0.395 6.79e-03 4.80e-02
## [995,] 0.513 1.41e-04 7.94e-07
## [996,] 0.439 1.33e-05 4.38e-04
## [997,] 0.502 7.01e-01 4.25e-04
## [998,] 0.449 7.00e-07 6.17e-07
## [999,] 0.556 1.24e-01 8.00e-06
## [1000,] 0.517 5.03e-04 3.52e-055. Estimación puntual de los parámetros
A continuacipon se muestran gráficos que permiten analizar la convergencia y distribución de los parámetros.
fGetMode <- function(v) {
uniqv <- unique(v)
uniqv[which.max(tabulate(match(v, uniqv)))]
}
# ------------------------------------------------------------------------------
# Analisis descriptivo de parámetro
# ------------------------------------------------------------------------------
fAnalisisDescriptivo = function( Os,s){
xmedia = mean (Os)
xmedian = median(Os)
xmoda = fGetMode(Os)
xvar = var(Os)
cat("Estadisticos de la muestra \n",
"media :" ,xmedia,"\n",
"mediana :" ,xmedian,"\n",
"moda :" ,xmoda,"\n",
"varianza:" ,xvar)
ls_title = "Distribucion de parámetro "
par(mfrow=c(2,2))
plot( Os,type="l",xlab = "Iteraciones" , ylab = s,
main = "Convergencia de parámetro")
acf(Os,main = "Autocorrelación de la muestra")
hist(Os,main = ls_title,
ylab = "Frecuencia",
xlab = s, prob = TRUE, breaks = 50)
abline(v=xmedia, col='red', lwd=3)
par(mfrow=c(1,1))
return(c(xmedia,xmedian,xmoda,xvar))
}Parámetro: \(p\)
parm1 = fAnalisisDescriptivo(Os[,1],"p")## Estadisticos de la muestra
## media : 0.515
## mediana : 0.513
## moda : 0.561
## varianza: 0.00974
Parámetro: \(\alpha\)
parm2 = fAnalisisDescriptivo(Os[,2],"α")## Estadisticos de la muestra
## media : 0.137
## mediana : 0.000669
## moda : 0.0343
## varianza: 0.156
Parámetro: \(\beta\)
parm3 = fAnalisisDescriptivo(Os[,3],"β")## Estadisticos de la muestra
## media : 0.133
## mediana : 0.000786
## moda : 0.968
## varianza: 0.214
6. Cuadro comparativo de resultados
Presentamos un cuadro comparativo para comparar la estimación de \(p\) por máxima verosimilitud y el enfoque bayesiano versus su valor teórico.
df_resumen[,"MediaBY"] = c(parm1[1],parm2[1],parm3[1])
df_resumen## Parametro Teorico MediaMV MediaBY
## 1 p 0.65 0.52 0.515
## 2 α 0.137
## 3 β 0.133