Teorema de Bayes

El teorema de Bayes, en la teoría de la probabilidad, es una proposición planteada por el matemático inglés Thomas Bayes (1702-1761)1 y publicada póstumamente en 1763, que expresa:

la probabilidad condicional de un evento aleatorio A dado B en términos de la distribución de probabilidad condicional del evento B dado A y la distribución de probabilidad marginal de solo A.

Bayes

En términos más generales y menos matemáticos, el teorema de Bayes es de enorme relevancia puesto que vincula la probabilidad de A dado B con la probabilidad de B dado A. Es decir, por ejemplo, que sabiendo la probabilidad de tener un dolor de cabeza dado que se tiene gripe, se podría saber (si se tiene algún dato más), la probabilidad de tener gripe si se tiene un dolor de cabeza. Muestra este sencillo ejemplo la alta relevancia del teorema en cuestión para la ciencia en todas sus ramas, puesto que tiene vinculación íntima con la comprensión de la probabilidad de aspectos causales dados los efectos observados.

Fórmula de Bayes

Con base en la definición de probabilidad condicionada se obtiene la Fórmula de Bayes, también conocida como Regla de Bayes:

Fórmula del Teorema de Bayes

Visualización del teorema de bayes

Teorema de Bayes

Análisis Bayesianos

Los análisis Bayesianos son similares a los que vimos en en el sentido que dependen explícitamente de modelos probabilísticos para los datos.

Es decir, tenemos que definir un modelo de datos. La gran diferencia es que con Bayes, podemos obtener distribuciones de probabilidades para todas las cantidades no observadas, incluyendo parámetros, valores perdidos o nuevas observaciones que todavía no hemos hecho. De esta manera, los análisis Bayesianos nos permiten cuantificar incertidumbre y armar modelos realistas que tienen en cuenta por ejemplo observaciones imperfectas.

Como vimos en la teórica, la regla de Bayes planteada en términos de datos y parámetros es:

\[ p(\boldsymbol{\theta} \lvert \boldsymbol{y}) = \frac{p(\boldsymbol{y} \lvert \boldsymbol{\theta}) p(\boldsymbol{\theta)}}{\int p(\boldsymbol{y} \lvert \boldsymbol{\theta)} p(\boldsymbol{\theta)} d \boldsymbol{\theta} } \]

Es decir que la probabilidad posterior de los parámetros θ dado que observamos los datos y es igual al likelihood multiplicado por las previas y dividido por la probabilidad total de los datos. La función de likelihood nos da la probabilidad de observar los datos condicional al valor de los parámetros p(y|θ). La previa de los parámetros p(θ)

refleja los posibles valores de los parámetros de acuerdo con nuestras “creencias” previas, o los resultados de estudios anteriores, o lo que nos parece que tiene sentido para el sistema de estudio (en definitiva, en base a información previa). Finalmente, la probabilidad total de los datos se obtiene integrando la función de lilkelihood sobre los posibles valores de los parámetros que define la previa. Como veremos más adelante, los análisis Bayesianos combinados con métodos numéricos permiten analizar modelos con muchos parámetros, niveles de variabilidad y variables “ocultas”, pero primero vamos a empezar por casos simples donde podemos calcular las posteriores directamente.

Para entender bien cómo es todo el proceso, vamos a simular los datos.

Imaginen que queremos estudiar la remoción de frutos en 30 plantas. En cada una de las plantas marcamos 20 frutos y contamos cuántos son removidos por dispersores luego de un tiempo fijo. Si suponemos que un buen modelo para este tipo de datos es una distribución Binomial con una probabilidad de éxito fija hacemos:

set.seed(1234)
nobs <- 30 #número de obserbaciones (plantas)
frutos <- rep(20,nobs)
p_rem <- 0.2 #Probabilidad de remoción por fruto Previa (a priori)
removidos <- rbinom(nobs, size = frutos, prob=p_rem)
removidos

##  [1] 2 4 4 4 6 5 0 3 5 4 5 4 3 7 3 6 3 3 2 3 3 3 2 1 3 6 4 7 6 1

El modelo de datos (cuántos frutos son removidos) es una Binomial con número de pruebas (la cantidad de frutos disponibles) conocido. Para hacer un análisis Bayesiano de estos datos, tenemos que definir una previa para la probabilidad de éxito (p_rem). Esa previa tiene que tomar valores continuos entre 0 y 1. Una opción sería una distribución uniforme con esos límites, pero si usamos una distribución Beta, es posible obtener un resultado analítico para la posterior. En este caso, la posterior es otra distribución Beta pero con sus parámetros actualizados en base a las observaciones. Se dice entonces que la distribución Beta es la conjugada de la Binomial. Si la previa de la tasa de remoción por fruto es una distribución Beta con parámetros α y β, actualizamos los valores de α y β en base a la cantidad de éxitos y fracasos obervados.

La posterior de la tasa de remoción por fruto es entonces una Beta con α=∑y, β=∑(n−y) donde y representa a los frutos removidos de los n disponibles. Veamos como hacer esto en R.

#previas
alpha <- 1
beta <- 1
alpha_p <- alpha + sum(removidos)
beta_p <- alpha + sum(frutos-removidos)

Para obtener el valor esperado de una distribución Beta hacermos

alpha_p/(alpha_p + beta_p)

## [1] 0.1877076

Eso nos da un estimador puntual de la probabilidad de remoción por fruto p_rem. Para tener una medida de incertidumbre alrededor de este valor, podemos ver los cuantiles de la posterior

qbeta (c(0.025, 0.975), alpha_p, beta_p)

## [1] 0.1575462 0.2198340

También podemos graficar la distribución posterior y compararla con la previa para ver cuánto aprendimos haciendo el análisis.

op <- par(cex.lab = 1.5, font.lab = 2, cex.axis = 1.3, las = 1, bty = "n")
curve(dbeta(x, alpha + sum(removidos), beta + sum(frutos - removidos)), lwd = 2, 
    ylab = "Densidad de probabilidad", xlab = "Probabilidad de remoción")
curve(dbeta(x, 1, 1), lwd = 2, col = "gray", add = TRUE)
text(0.6, 2.5, "previa")
text(0.35, 12, "posterior")

¿Cómo sabemos si estas estimaciones tienen sentido para nuestros datos? En este caso, la pregunta es trivial porque conocemos cómo se generaron los datos, pero cuando trabajamos con datos de verdad, el modelo de datos es un supuesto y tenemos que ver si ese supuesto tiene sentido.

Una opción para contestar esa pregunta es hacer simulaciones a partir de la posterior y compararlas con los datos.

nreps <- 10000
vals <- 0:20  # posibles valores de remoción
res <- matrix(NA, nreps, length(vals) - 1)  # matriz para resultados
p_sim <- rbeta(nreps, alpha_p, beta_p)  # muestra aleatoria de la posterior

for (i in 1:nreps) {
    tmp <- rbinom(nobs, frutos, p_sim[i])
    res[i, ] <- hist(tmp, right = FALSE, breaks = vals, plot = FALSE)$density
}

plot(table(removidos)/nobs, xlim = c(0, 10), ylim = c(0, 0.6), ylab = "frecuencia", 
    type = "p", pch = 19)
library(coda)
ci <- HPDinterval(as.mcmc(res))
lines(0:19, ci[, 2])
lines(0:19, ci[, 1])

Ahora, según lo aprendido hasta ahora y lo que puedan investigar.

Haga 2 ejercicios del teorema de bayes aplicado a su área de estuido.

1.- Una resistencia para una placa madre de una marca de televisores es de 100 Ohms, estas resistencias tienen un error del +- 1% lo cual ocasiona que se caliente un disipador cuando está por debajo de ese margen de error y que no llegue suficiente voltaje a un circuito siguiente si está por encima del margen. Además el fabricante asegura que el 0.1% de estas resistencias salen defectuosas desde su manufacturación. Esto quiere decir que una de cada mil son defectuosas desde su fabricación y de las funcionales a la hora de manufacturarse una de cada cien trabaja fuera del margen nominal resistivo.

Código de colores de resistencias

set.seed(1234)
nobs <- 100 #número de observaciones (resistencias)
r <- rep(50,nobs) #De 50 resistencias que agarremos
p_rem <- 0.01 #es la probabilidad de que salga una defectuosa por error o de fábrica luego de que salga a ensamblaje de las placas.
defectuosas <- rbinom(nobs, size = r, prob=p_rem)
defectuosas

##   [1] 0 1 1 1 1 1 0 0 1 0 1 0 0 2 0 1 0 0 0 0 0 0 0 0 0 1 0 2 1 0 0 0 0 0 0 1 0
##  [38] 0 3 1 0 1 0 1 0 0 1 0 0 1 0 0 1 0 0 0 0 1 0 1 1 0 0 0 0 1 0 0 0 0 0 1 0 1
##  [75] 0 0 0 0 0 1 2 0 0 0 0 1 0 0 0 1 0 1 0 0 0 0 0 0 0 1

#previas
alpha <- 1
beta <- 1
alpha_p <- alpha + sum(defectuosas)
beta_p <- alpha + sum(r - defectuosas)

Para obtener el valor esperado de una distribución Beta hacermos

alpha_p/(alpha_p + beta_p)

## [1] 0.007596961

Eso nos da un estimador puntual de la probabilidad de remoción por fruto p_rem. Para tener una medida de incertidumbre alrededor de este valor, podemos ver los cuantiles de la posterior.

qbeta (c(0.025, 0.975), alpha_p, beta_p)

## [1] 0.005382586 0.010183600

También podemos graficar la distribución posterior y compararla con la previa para ver cuánto aprendimos haciendo el análisis.

op <- par(cex.lab = 1.5, font.lab = 2, cex.axis = 1.3, las = 1, bty = "n")
curve(dbeta(x, alpha + sum(defectuosas), beta + sum(r - defectuosas)), lwd = 2, 
    ylab = "Densidad de probabilidad", xlab = "Probabilidad de remoción")
curve(dbeta(x, 1, 1), lwd = 2, col = "gray", add = TRUE)
text(0.1, 2.5, "previa")
text(0.1, 15, "posterior")

Estos datos nos dicen que es muy, pero muy difícil e improbable que falle una resistencia en una televisión nueva, sin uso.

Una opción para contestar esa pregunta es hacer simulaciones a partir de la posterior y compararlas con los datos.

nreps <- 10000
vals <- 0:20  # posibles valores de remoción
res <- matrix(NA, nreps, length(vals) - 1)  # matriz para resultados
p_sim <- rbeta(nreps, alpha_p, beta_p)  # muestra aleatoria de la posterior

for (i in 1:nreps) {
    tmp <- rbinom(nobs, r, p_sim[i])
    res[i, ] <- hist(tmp, right = FALSE, breaks = vals, plot = FALSE)$density
}

plot(table(defectuosas)/nobs, xlim = c(0, 10), ylim = c(0, 0.6), ylab = "frecuencia", 
    type = "p", pch = 19)
library(coda)
ci <- HPDinterval(as.mcmc(res))
lines(0:19, ci[, 2])
lines(0:19, ci[, 1])

2.- Una pantalla OLED de un teléfono celular según el fabricante comienza a fallar después de 20 meses, asegurando que de cada 1000 terminales el 25% presentaría efecto fantasma por pantalla quemada, y el 1% presenta un vector de pixeles (vertical u horizontal) completo que deja de funcionar.

Pánel de pantalla AMOLED quemada

Líneas de pixeles muertos

set.seed(1234)
nobs <- 1000 #número de observaciones (celulares)
cel <- rep(100,nobs) #De 100 celulares después de 20 meses que analicemos
p_rem <- 0.03 #es la probabilidad de que salga una defectuosa de un vector de pixeles pero no de pantalla fantasma.
defecto <- rbinom(nobs, size = r, prob=p_rem)
defecto

##    [1] 0 2 2 2 3 2 0 1 2 1 2 1 1 3 1 3 1 1 0 1 1 1 0 0 1 2 1 3 3 0 1 1 1 1 0 2 0
##   [38] 1 5 2 1 2 1 2 1 1 2 1 1 2 0 1 2 1 0 1 1 2 0 3 3 0 1 0 1 2 1 1 0 2 0 3 0 2
##   [75] 0 1 1 0 1 2 3 1 0 1 0 3 1 1 0 3 0 3 0 0 0 1 1 0 1 2 0 2 1 0 0 1 0 0 1 0 2
##  [112] 0 4 0 1 3 4 1 0 2 2 3 5 4 1 1 1 1 1 1 4 2 0 1 3 1 3 2 2 3 1 5 1 1 1 2 2 2
##  [149] 4 2 1 1 0 3 1 5 2 6 1 1 1 2 1 1 0 2 1 0 2 1 2 2 2 1 1 2 1 2 0 1 1 1 2 0 4
##  [186] 0 3 2 1 2 2 5 0 3 2 3 3 2 4 2 2 1 1 2 1 2 1 1 0 5 2 1 0 2 2 3 2 2 1 3 1 0
##  [223] 1 1 0 1 2 1 1 1 2 2 0 2 2 0 2 1 1 5 1 1 0 2 4 1 2 2 4 2 1 2 1 1 0 2 1 0 2
##  [260] 1 3 0 2 1 0 0 2 0 0 2 1 2 1 2 0 1 3 1 1 1 1 1 3 2 2 1 1 1 2 3 2 1 4 1 2 2
##  [297] 2 0 2 1 1 4 1 1 4 1 0 5 1 2 3 2 1 1 0 0 0 1 1 3 0 3 1 0 1 1 1 0 2 3 2 3 4
##  [334] 0 1 1 1 0 3 1 4 1 1 4 1 1 0 2 3 0 0 2 0 2 6 4 2 2 1 2 0 2 1 1 3 4 2 1 3 2
##  [371] 4 1 3 1 1 1 2 0 1 2 1 0 0 0 0 0 2 0 5 4 2 2 2 2 0 3 2 1 0 1 2 1 2 2 2 0 2
##  [408] 0 2 1 2 2 4 0 3 1 1 2 0 1 2 3 1 0 3 0 0 1 1 1 2 2 1 3 1 2 1 3 2 2 3 3 2 1
##  [445] 5 1 3 2 2 0 1 2 6 4 2 2 1 2 2 0 1 3 1 1 0 0 1 1 2 2 0 1 0 4 1 4 1 1 2 0 1
##  [482] 2 2 1 0 0 2 2 2 0 2 2 3 1 2 1 2 0 2 2 2 2 3 0 2 3 2 3 0 2 1 4 3 2 0 2 2 0
##  [519] 2 1 3 1 1 3 2 2 0 0 2 1 1 2 1 2 0 5 2 2 2 1 1 2 4 0 2 1 3 3 4 3 1 1 2 0 3
##  [556] 4 1 0 2 3 0 5 2 3 1 3 2 2 0 0 2 2 0 2 1 1 0 1 1 3 3 0 2 1 0 1 1 4 4 2 4 3
##  [593] 0 2 1 2 2 2 2 0 1 0 0 2 1 1 3 0 1 1 1 2 1 0 4 1 1 2 2 1 3 0 2 3 1 3 1 4 1
##  [630] 2 3 2 1 2 4 3 1 3 1 1 2 0 1 3 1 1 0 3 1 2 1 2 3 1 0 1 3 0 1 3 2 4 3 0 1 1
##  [667] 2 1 1 2 0 3 1 1 4 3 1 3 0 4 1 0 2 1 3 1 3 0 1 1 3 1 1 1 0 0 1 1 1 2 0 2 2
##  [704] 1 0 0 1 2 0 2 2 4 4 0 0 0 1 0 4 0 3 0 4 0 1 5 1 1 1 5 3 0 0 1 2 2 1 2 5 3
##  [741] 1 0 0 0 2 0 1 2 0 2 1 2 1 3 1 3 2 3 1 2 0 2 0 2 1 2 3 2 2 1 3 3 3 0 2 1 0
##  [778] 1 2 1 4 0 0 1 0 2 3 0 2 0 2 0 2 1 2 1 1 3 0 3 0 3 1 1 1 2 0 1 2 1 2 0 3 2
##  [815] 0 1 1 2 1 1 0 1 1 1 1 2 0 1 1 1 1 2 3 4 1 2 2 0 2 2 0 2 2 0 1 3 2 4 2 2 1
##  [852] 4 4 0 1 2 0 0 2 4 0 1 0 1 0 3 1 4 1 1 4 3 4 1 2 0 0 0 1 2 2 1 1 1 2 3 0 0
##  [889] 0 0 3 2 2 1 1 1 2 2 4 3 2 1 3 1 3 0 1 1 2 1 1 1 0 1 0 0 1 4 4 2 1 1 0 1 2
##  [926] 1 2 1 0 0 0 5 0 2 0 1 0 1 1 2 3 2 3 2 2 1 1 0 1 1 2 4 1 3 3 0 3 1 0 2 1 1
##  [963] 3 2 0 0 2 1 2 1 7 2 0 2 0 2 2 1 4 2 2 3 0 2 5 1 0 5 0 0 1 2 5 0 2 0 2 1 4
## [1000] 0

#previas
alpha <- 1
beta <- 1
alpha_p <- alpha + sum(defecto)
beta_p <- alpha + sum(cel - defecto)

Para obtener el valor esperado de una distribución Beta hacermos

alpha_p/(alpha_p + beta_p)

## [1] 0.01547969

Eso nos da un estimador puntual de la probabilidad de remoción por fruto p_rem. Para tener una medida de incertidumbre alrededor de este valor, podemos ver los cuantiles de la posterior.

qbeta (c(0.025, 0.975), alpha_p, beta_p)

## [1] 0.01472378 0.01625395

También podemos graficar la distribución posterior y compararla con la previa para ver cuánto aprendimos haciendo el análisis.

op <- par(cex.lab = 1.5, font.lab = 2, cex.axis = 1.3, las = 1, bty = "n")
curve(dbeta(x, alpha + sum(defecto), beta + sum(cel - defecto)), lwd = 2, 
    ylab = "Densidad de probabilidad", xlab = "Probabilidad de remoción")
curve(dbeta(x, 1, 1), lwd = 2, col = "gray", add = TRUE)
text(0.1, 30, "previa")
text(0.1, 390, "posterior")

Estos datos nos dicen que es muy, pero muy difícil e improbable que un celular con pantalla OLED se descomponga en un vector de pixeles sin que la pantalla esté quemada provocando el efecto fantasma, pues la tecnología OLED trabaja pixel a pixel, y la probabilidad tan alta de que la pantalla haga el efecto fantasma hace muy difícil esta situación particular.

Una opción para contestar esa pregunta es hacer simulaciones a partir de la posterior y compararlas con los datos.

nreps <- 10000
vals <- 0:20  # posibles valores de remoción
res <- matrix(NA, nreps, length(vals) - 1)  # matriz para resultados
p_sim <- rbeta(nreps, alpha_p, beta_p)  # muestra aleatoria de la posterior

for (i in 1:nreps) {
    tmp <- rbinom(nobs, r, p_sim[i])
    res[i, ] <- hist(tmp, right = FALSE, breaks = vals, plot = FALSE)$density
}

plot(table(defecto)/nobs, xlim = c(0, 10), ylim = c(0, 0.6), ylab = "frecuencia", 
    type = "p", pch = 19)
library(coda)
ci <- HPDinterval(as.mcmc(res))
lines(0:19, ci[, 2])
lines(0:19, ci[, 1])

Redacción

Este análisis nos hace predecir que es probable que pase con una muestra futura, pero no con simplme consideración de la probabilidad simple, sino también considerando lo que ya se conoce aparte de los datos, o sease, aunque tenga una probabilidad de que suceda algo, no sucederá así siempre, y este análisis se aproxima a esa hipótesis realista, y se confirmó con los dos ejemplos de la carrera, los cuales no fueron exactamente lo que la probabilidad a priori suponía, sino que hubo una aproximación a esta, que es más realista que la misma probabilidad simple.

U2A12

Jorge Valenzuela Parra

11/11/2020

Teorema de Bayes

Fórmula de Bayes

Análisis Bayesianos

Redacción