Teorema de bayes

setwd("~/PYE1213")

Importar paquetes

library(pacman)
p_load('coda','prettydoc','readr')

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:

Bayes

  • Visualización del teorema de bayes

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} } \]

Likehood : verosimilitud

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.

  • Ejercicio 1: Análisis bayesiano

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 observaciones (Plantas)
frutos <- rep(20,nobs)
p_rem <- 0.2 # Probabilidad de remoción por fruto
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
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 hacemos

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")

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)
ci <- HPDinterval(as.mcmc(res))
lines(0:19, ci[, 2])
lines(0:19, ci[, 1])

url <- "https://sites.google.com/site/modelosydatos/quintral.txt"
quintral <- read.table(url,header = TRUE)

alpha <- 1
beta <- 1

alpha_p <- alpha + sum(quintral$Removidos)
beta_p <- alpha + sum(quintral$Frutos-quintral$Removidos)

op <- par(cex.lab = 1.5, font.lab = 2, cex.axis = 1.3, las = 1, bty = "n")
curve(dbeta(x, alpha + sum(quintral$Removidos), beta + sum(quintral$Frutos-quintral$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")

nreps <- 10000
nobs <- length(quintral$Removidos)

vals <- 0:60  # 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, quintral$Frutos, p_sim[i])
    res[i, ] <- hist(tmp, right = FALSE, breaks = vals, plot = FALSE)$density
}

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

Asignación:

Encuentre 2 ejercicios de aplicaciones del teorema de Bayes y ejecútelos, estos ejercicios pueden o no ser relacionados con su carrera.

Ejercicio 1:

El 20% de los empleados de una empresa son ingeniero y otro 20% son economistas. El 75% de los ingenieros ocupan un puesto directivo y el 50% de los economistas también, mientas que los no ingenieros y los no economistas solamente el 20% ocupa un puesto directivo. ¿Cuál es la probabilidad de que un empleado directivo elegido al azar sea un ingeniero?

#set.seed(1234)
# Probabilidades primarias
probIng <- c(0.2)
probEco <- c(0.2)
probOtro <- c(0.6)
# Ingreso de las probabilidades secundarias
probDirectalIng <- c(0.75)
probDirectalEco <- c(0.50)
probDirectalOtr <- c(0.20)
Probabilidad1 <- probDirectalIng*probIng/sum(probDirectalIng*probIng,probDirectalEco*probEco,probDirectalOtr*probOtro)
Probabilidad2 <- probDirectalEco*probEco/sum(probDirectalIng*probIng,probDirectalEco*probEco,probDirectalOtr*probOtro)
Probabilidad3 <- probDirectalOtr*probOtro/sum(probDirectalIng*probIng,probDirectalEco*probEco,probDirectalOtr*probOtro)
A <- rbind(Probabilidad1,Probabilidad2,Probabilidad3)
A
##                    [,1]
## Probabilidad1 0.4054054
## Probabilidad2 0.2702703
## Probabilidad3 0.3243243

Ejercicio 2:

Almacen de ropa

El 60% de las ventas en unos grandes almacenes corresponde a artículos con precios rebajados. Los clientes devuelven el 15% de los artículos que compran rebajados, porcentaje que disminuye al 8% si los artículos han sido adquiridos sin rebajas.

  1. Determina el porcentaje global de artículos devueltos.

  2. ¿Qué porcentaje de artículos devueltos fueron adquiridos con precios rebajados?

\(PR=60\%\) \(PDR=15\%\) \(PNDR=85\%\)

\(PNR=40\%\) \(PDNR=8\%\) \(PNDNR=92\%\)

  1. Determina el porcentaje global de artículos devueltos.
PR <- 0.6
PDR <- 0.15
PNDR <- 0.85
PNR <- 0.4
PDNR <- 0.08
PNDNR <- 0.92

PD <- (PR * PDR) + (PNR*PDNR)
cat("Probabilidad de que regresen un producto es de: ",PD*100,"%")
## Probabilidad de que regresen un producto es de:  12.2 %
  1. ¿Qué porcentaje de artículos devueltos fueron adquiridos con precios rebajados?
(PR*PDR)/PD
## [1] 0.7377049