Introducción

La máquina de cafés de la escuela de informática falla con cierta regularidad. Para intentar estimar la probabilidad con la que cada café es derramado, pedimos 100 cafés y anotamos el número de cafés derramados. El resultado dio lugar a un total de 7 cafés derramados. ¿Podríamos concluir que la probabilidad de error de la máquina es de un 7%?¿Como podemos cuantificar la incertidumbre de nuestra estimación?

Modelando el proceso de servir un café.

Los cafés de Bernouilli

Vamos a hacer una serie de suoposiciones sobre la máquina de cafés:

  • La máquina tiene una probabilidad p de cometer un error al servir un café. Este valor p es desconocido, y a través de experimentos queremos estimar su verdadero valor.

  • La probabilidad p de cometer un error permanece constante en el tiempo (al menos en el tiempo que llevaban haciendo el doctorado)

  • Cuando la máquina sirve un café bien o mal, no influye en la probabilidad de que el siguiente café esté bien o mal. Los cafés son independientes.

Con lo cual, los resultados del experimento serán: {bien servido, mal servido} a lo que asignaremos los valores \(\{0,1\}\). \[Prob(X = x: x \in \{0, 1\}) = \left\{ \begin{array}{lcc} p & si & x = 1 \\ \\ 1 - p & si & x = 0 \end{array}\right.\]

Lo que también se puede escribir como:

\[Prob(X = x: x \in \{0, 1\}) = p^x(1 -p)^{1 - x}\]

Todo esto lo podremos resumir más diciendo que el resultado de servir un café con éxito o no, es aleatorio y su probabilidad queda distribuida por un proceso de Bernoulli, lo cual queda matemáticamente expresado como:

\[X \sim Bern(p)\]

Los cafés binomiales

Al pedir 100 cafés, hemos realizado el experimento de Bernouilli 100 veces, y esto se traduce en que tenemos una “colección” de resultados.

\[X = \left\{ X_1 = x_1, X_2 = x_2, X_3 = x_3...X_{100} = x_100\right\} \;\;\;\;\;\; \forall X_i \sim Bern(p)\]

Estamos interesados en sumar el número de cafés mal servidos. EL número de cafés mál servidos será otra variable aleatoria calculada a partir de los 100 experimentos independientes de Bernoulli. Llamemos a esa nueva variable aleatoria \(Y\) tal que:

\[Y = y: y = \sum_i^{100}{x_i}: \;\; y \in \mathbf\{0, 1 , 2 ...\}\].

Como vemos \(Y\) es una nueva variable aleatoria y queda soportada por el conjunto de los números naturales con el cero incluido.

Es fácil comprobar que la probabilidad de una configuración determinada de tamaño n de experimentos de Bernouilli es:

\[prob\left\{X_1 = x_1, X_2 = x_2, ..., X_n = x_n\right\} = p^{y}(1 - p)^{n - y} = p^{\bar{x}n}(1 - p)^{n - \bar{x}n}\] Nos interesa el número de cafés derramados o no derramados, independientemente del orden en el que estos hayan ido apareciendo. Para calcular esta probabilidad, solo tendríamos que sumar todas las configuraciones en las que resulte que \(Y = y\) dado un tamaño de \(n\). Por ejemplo, si resulta que hay \(C^n_y\) configuraciones distintas en las que \(Y = y\) dado el tamaño \(n\) entonces simplemente multiplicaríamos por \(C_y^n\) (que equivaldría a sumar \(C_y^n\) veces) la probabilidad de la configuración, algo así:

\[prob\{Y = y\} = C^n_y*p^{y}(1 - p)^{n - y}\] El número de configuraciones en el que resulta que \(Y = y\) dado \(n\) es el coeficiente binomial, y se calcula:

\[C_y^n = \binom{n}{y} = \frac{n!}{y!(n-y)!}\]

Visto todo esto, podemos generalizar el fenómeno aleatorio que gobierna nuestra experimento sobre la máquina de cafés de la siguiente forma: \[Prob\left\{Y = y\right\} = \binom{n}{y} p^{y}(1 - p)^{n - y}\]

De la misma forma que hicimos con el proceso de Bernouilli, una variable aleatoria \(Y\) cuya probabilidad se distribuye de forma binomial de parámetros p y n se expresará:

\[Y \sim Binom(n, p)\]

Conclusión y reflexión sobre la construcción del modelo

Para construir un modelo binomial, solo nos ha hecho falta unas cuantas hipótesis:

  1. Cada experimento es independiente
  2. Existe un valor p que gobierna la probabilidad de cada experimento
  3. Cada experimento solo tiene dos posibles resultados complementarios: éxito o fracaso

No necesitamos asumir nada más para construir este modelo.

Inferencia sobre el parámetro

Justificación de la inferencia bayesiana

Plausibilidad relativa de un conjunto de hipótesis

Una vez que ya tenemos desarrollado nuestro modelo binomial podemos utilizarlo para simular nuestra máquina de cafés.

Supongamos que asumimos que la verdadera probabilidad de error p = 0.07. Entonces:

prob_mle <- dbinom(x = 7, size = 100, prob = 0.07)
print(prob_mle)
## [1] 0.154499

Supongamos ahora que nos apetece asumir como verdadero otro valor, uno muy cercano al que maximiza la similitud a la verdad que traen los dato, por ejemplo \(p = 0.08\), entonces:

prob_08 <- dbinom(x = 7, size = 100, prob = 0.08)
print(prob_08)
## [1] 0.1439538

¿cuántas veces más probable sería la hipótesis de máxima verosimilitud \(p = 0.07\) frente a una hipótesis muy cercana a ésta\(p = 0.08\)?

odds_mle_08 <- prob_mle/prob_08
print(odds_mle_08)
## [1] 1.073254

Como vemos, la hipótesis de la máxima verosimilitud es aproximadamente un \(7\%\) más probable(en términos relativos) que si asumieramos que en realidad la verdadera probabilidad \(p = 0.08\). Esto significa que sería muy verosimil que el verdadero parámetro fuera \(p = 0.08\) y no \(0.07\), aunque ciertamente es más verosimil el segundo, pero eso no lo convierte en el verdadero. ¿por qué deberíamos escoger el valor de máxima verosimilud y descartar el otro, pudiendo escoger ambos? Es más, la estadística bayesiana nos ofrece un marco lógico para poder escoger todos los valores a la vez en proporción a su plausibilidad relativa.

Siguiendo con la línea de razonamiento que hemos adquirido, supongamos que enfoco mis sospechas en estos 5 valores: \(p \in {0.05, 0.07, 0.09, 0.11, 0.13}\), teniendo en cuenta que nuestro experimento con el que venimos tratando (y = 7, n = 100) podemos simular el experimento tantas veces como queramos:

nsim = 1e4
param_space <- list("0.05" = 0.05, "0.07" = 0.07, "0.09" = 0.09, "0.11" = 0.11, "0.13" = 0.13)
df_sim_each_param <- map_df(param_space, ~rbinom(n = nsim,  size = 100, prob = .x))
df_sim_each_param
## # A tibble: 10,000 x 5
##    `0.05` `0.07` `0.09` `0.11` `0.13`
##     <int>  <int>  <int>  <int>  <int>
##  1      1      3      6     13     17
##  2      3      7      9      8     18
##  3      6     13      9     14     12
##  4      6      8      4     14     18
##  5      4      9      6     10     14
##  6      9      9      9     15     13
##  7      3      5     14     13     14
##  8      7     10      8      9     10
##  9      6     10      7      8     11
## 10      4      6     10      7     13
## # ... with 9,990 more rows

En primer lugar vamos a ver la distribución del experimento para \(p = 0.05\). Como es natural, este experimento encuentra su máxima verosimilitud en \(5\) (es decir, bajo estas condiciones, lo más verosímil sería que el número de cafés derramados de 100 fueran 5). Sin embargo observamos que bajo estas condiciones, derramar 7 cafés de 100 también es probable. De hecho, resultó en 7/100 en aproximadamente 1000 simulaciones de 10000 (téngase en cuenta que el 5 salió unas 1700 sobre 10000, lo cual sugiere que el mejor estimador no llegaría a doblar la probabilidad).

df_sim_each_param %>% 
  ggplot(aes(x = `0.05`)) +
  geom_bar(position = "dodge", alpha = 0.6, fill = "red") +
  geom_vline(xintercept = 7)

Lo que nos preguntamos ahora es: ¿para cada uno de los parámetros que hemos usado \((0.05, 0.07, ..., 0.13)\), cuantas veces el experimento resultó en 7 fallos sobre 100 en cada uno de ellos? Vamos a verlo.

df_summary <- df_sim_each_param %>% 
  gather() %>%
  group_by(key, value) %>%
  summarise(n = n()) %>%
  mutate(prob = n/sum(n)) %>% 
  mutate(alp = if_else(value == 7, 1, 0.9))

plot1 <- df_summary  %>%
  ggplot(aes(x = value, y = n, fill = key, alpha = alp )) +
  geom_bar(stat = "identity", position = "dodge") +
  geom_vline(xintercept = 7)

plot2 <- df_summary %>%
  ggplot(aes(x = value, y = n, fill = key, alpha = alp )) +
  geom_bar(stat = "identity") +
  geom_vline(xintercept = 7)

grid.arrange(plot1, plot2)

Como vemos, cada parámetro queda ponderado por el número de veces que “la evidencia (los datos)” es obtenida al simular la máquina frente a distintas hipótesis. Es como si las simulaciones salieran a la defensa de las hipótesis, y luego observáramos cuantas razones tiene tal parámetro para defender que él es el verdadero causante de los datos, siendo las simulaciones sus argumentos. Vamos ahora a calcularlo:

df_freq <- df_sim_each_param %>%
  gather() %>%
  filter(value == 7) %>%
  group_by(key) %>%
  summarise(freq_abs_seven = n()) %>%
  mutate(freq_rel_seven = freq_abs_seven / sum(freq_abs_seven) ) 

print(df_freq)
## # A tibble: 5 x 3
##   key   freq_abs_seven freq_rel_seven
##   <chr>          <int>          <dbl>
## 1 0.05            1072         0.230 
## 2 0.07            1550         0.333 
## 3 0.09            1151         0.247 
## 4 0.11             626         0.134 
## 5 0.13             261         0.0560
df_freq %>%
  ggplot(aes(x = key, y = freq_abs_seven, fill = key)) +
  geom_bar(stat = "identity", alpha = 0.7) +
  geom_text(aes(label = freq_abs_seven))

Si tengo en cuenta todos las formas en las que resultado es \(y = 7\) entonces tendré que sumar 1035 + 1479 + 1157 + 620 + 241 = 4532. Es decir, si tengo en cuenta los 5 modelos, \(y = 7\) se obtuvo un total de 4532 veces. Teniendo en cuenta esto podríamos decir que de 4532 formas en las que se consiguió que \(y = 7\). Vamos a resumir lo que tenemos:

  • De 50000 posibilidades(5 hipotesis con 10000 simulaciones cada una), tengo 4542 formas de probar como se puede conseguir unos datos tal que \(y = 7\)
  • De estas 4542 formas, tengo 1035 razones para pensar que \(p = 0.05\)
  • De estas 4542 formas, tengo 1479 razones para pensar que \(p = 0.07\)
  • De estas 4542 formas, tengo 1157 razones para pensar que \(p = 0.09\)
  • De estas 4542 formas, tengo 620 razones para pensar que \(p = 0.11\)
  • De estas 4542 formas, tengo 241 razones para pensar que \(p = 0.13\)
df_freq %>%
  ggplot(aes(x = key, y = freq_rel_seven, fill = key)) +
  geom_bar(stat = "identity", alpha = 0.7) +
  geom_text(aes(label = round(freq_rel_seven, 2)))

Cada una de estas barras estaría mostrándonos la “plausibilidad de relativa de cada una de las hipótesis” o “dadas todas las formas posibles de conseguir los datos, en que proporción cada hipótesis se haría responsable de dichos datos”. ¿No es esto una medida de probabilidad?. Como pueden observar, justo aquí es donde nos metemos en el gran debate sobre “¿qué es exactamente la probabilidad?”. Históricamente, el punto de vista frecuentista, no aceptaría esta “plausubilidad relativa” como una verdadera medida de probabilidad, y desde el punto de vista bayesiano, no hay nada que impida interpretar esta medida como probabilidad. Al fin y al cabo, por un lado cumple con la definción axiomática de kolmogorov y por otro lado va en linea con nuestra intuición interna del significado de probabilidad.

Desmitificando el teorema de bayes

El teorema de Bayes suele estar definido de la siguiente forma:

\(P(p|datos) = \frac{P(datos|p)P(p)}{p(datos)} = \frac{P(datos|p)P(p)}{\sum p(datos|p)}\)

Por ahora vamos a modificar esta fórmula, tal que: \(P(p|datos) = \frac{P(datos|p)}{p(datos)} = \frac{P(datos|p)}{\sum p(datos|p)*1}\)

Volviendo la vista arriba, vamos a realizar el siguiente cálculo, pero ahora con probabilidades. En primer lugar la suma de todas las formas posibles en las que obtengo los datos, no es más que la frecuencia relativa con la que se obtienen los datos, o dicho de otra forma la PROBABILIDAD DE LOS DATOS, la $ + + + + = $ es decir: \(0.1035 + 0.1479 + 0.1157 + 0.0620 + 0.0241 = 0.4532\). Entonces

Plausibilidad relativa para cada valor del espacio parámetrico

Lo que vamos a suponer ahora es que el parámetro p puede ser cualquier valor situado entre cero y uno. Es decir, el espacio parámetrico (en el que nuestro parámetro puede existir) está situado en intervalo continuo y cerrado delimitado por cero y uno. Obviamente ya no tiene sentido asignar una cantidad de plausibilidad o probabilidad, sino una densidad de plausibilidad o probabilidad.

\(p \in [0,1]\)

Para simular esta situación vamos a hacer lo siguiente:

nsim = 1e6
# Dividimos el intervalo 0 y 1 en 10000000 (asi modelaremos el continuo)
p_grid <- seq(from = 0, to = 1, length.out = nsim)
# simulamos con este grid
rbin_sim <- rbinom(nsim, prob = p_grid, size = 100)
# Lo ponemos en un data_frame
df_sim <- as_tibble(data.frame(p_grid, rbin_sim))
print(df_sim)
## # A tibble: 1,000,000 x 2
##        p_grid rbin_sim
##         <dbl>    <int>
##  1 0                 0
##  2 0.00000100        0
##  3 0.00000200        0
##  4 0.00000300        0
##  5 0.00000400        0
##  6 0.00000500        0
##  7 0.00000600        0
##  8 0.00000700        0
##  9 0.00000800        0
## 10 0.00000900        0
## # ... with 999,990 more rows

Vamos ahora a dibujar un histograma como con los parámetros responsables de haber producido un resultado de \(y = 7\). En primer lugar dibujamos un histograma y después aproximamos la su función de densidad. Además añadimos una linea vertical en donde se encuentra la máxima versosimilitud \(p = 0.07\),

plot1 <- df_sim %>% 
  filter(rbin_sim == 7) %>%
  ggplot(aes(x = p_grid)) +
  geom_histogram(aes(y = stat(count / sum(count))), fill = "red", alpha = 0.6) +
  geom_vline(xintercept = 0.07)

plot2 <- df_sim %>% 
  filter(rbin_sim == 7) %>%
  ggplot(aes(x = p_grid)) +
  geom_density(fill = "red", alpha = 0.6) +
  geom_vline(xintercept = 0.07)

grid.arrange(plot1, plot2)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Solo por curiosidad, vamos a ver como quedaría si usaramos \(y \in \{0, 5, 50, 70, 100\}\)

df_sim %>% 
  filter(rbin_sim %in% c(0, 20, 40, 60, 80, 100)) %>%
  ggplot(aes(x = p_grid, y = rbin_sim, group = rbin_sim)) +
  geom_density_ridges(fill = "red")
## Picking joint bandwidth of 0.00458

Aquí podemos observar como se verían las distintas distribuciones de probabilidad sobre el espacio parámetrico dados distintos valores de \(y\).

Inferencia Bayesiana

Hemos comenzado mostrando que efectivamente podemos cuantificar la incertidumbre que existe sobre las causas que generan las evidencias, es decir, los datos. Vamos ahora a dar un paso adelante introduciendo el concepto de la probabilidad a priori.

Probabilidad a priori+

Volviendo al ejemplo en el que teníamos un conjunto de hipótesis sobre cual podría ser el parámetro, teníamos que:

  • 1035 se justifican con \(p = 0.05\)
  • 1479 se justifican con \(p = 0.07\)
  • 1157 se justifican con \(p = 0.09\)
  • 620 se justifican con \(p = 0.11\)
  • 241 se justifican con \(p = 0.13\)

Supongamos que hablamos con una persona experta en estas máquinas de cafés, alguien que lleva el control de calidad de éstas. Esta persona nos dice lo siguiente: “Dada mi experiencia como ingeniero fabricando estas máquinas, sé a ciencia cierta que en el mercado hay circulando máquinas de café tal que: el 70% de estas máquinas tienen una tasa de fallo de 0.05, el 20% de 0.07 y el restante de 0.09. Además, es físicamente imposible que la máquina pueda tener una tasa de fallo superior a 0.09”

¿Como podríamos utilizar esta información para inferir el verdadero parámetro?