Ejemplo: simulación tres cartas

baraja <- tibble(numero = 1:13) %>% 
   crossing(tibble(figura = c("C", "D", "T", "P"))) %>% 
   mutate(carta = paste(numero, figura))
nrow(baraja)
#> [1] 52
baraja
#> # A tibble: 52 x 3
#>    numero figura carta
#>     <int> <chr>  <chr>
#>  1      1 C      1 C  
#>  2      1 D      1 D  
#>  3      1 P      1 P  
#>  4      1 T      1 T  
#>  5      2 C      2 C  
#>  6      2 D      2 D  
#>  7      2 P      2 P  
#>  8      2 T      2 T  
#>  9      3 C      3 C  
#> 10      3 D      3 D  
#> # … with 42 more rows

cartas <- baraja$carta
exp_3_cartas <- function(cartas){
   sample(cartas, 3) 
}
set.seed(132185)
exp_3_cartas(cartas)
#> [1] "7 C"  "8 P"  "13 T"

Simulamos el experimento

sims <- map(1:20000, ~ exp_3_cartas(cartas))
sims[1:5]
#> [[1]]
#> [1] "4 C"  "8 T"  "12 D"
#> 
#> [[2]]
#> [1] "11 D" "6 P"  "7 C" 
#> 
#> [[3]]
#> [1] "12 C" "1 D"  "5 C" 
#> 
#> [[4]]
#> [1] "12 D" "9 D"  "13 T"
#> 
#> [[5]]
#> [1] "7 T"  "13 D" "13 T"

Sin más información, la probabilidad de corazón en la primera extracción es:

sims %>% 
   map_lgl(~ str_detect(.x[1], "C")) %>% 
   mean()
#> [1] 0.2474

Que debe estar alrededor de 1/4. Ahora condicionamos a que las cartas 2 y 3 son corazones:

sims_F <- sims %>%
   keep(~ str_detect(.x[2], "C") & str_detect(.x[3], "C"))
length(sims_F)
#> [1] 1203
sims_F[1:5]
#> [[1]]
#> [1] "9 T"  "3 C"  "12 C"
#> 
#> [[2]]
#> [1] "5 D"  "2 C"  "13 C"
#> 
#> [[3]]
#> [1] "5 C"  "10 C" "1 C" 
#> 
#> [[4]]
#> [1] "4 P" "6 C" "1 C"
#> 
#> [[5]]
#> [1] "8 C"  "13 C" "11 C"

Y sobre estas simulaciones hacemos el mismo cálculo de arriba:

sims_F  %>% 
   map_lgl(~ str_detect(.x[1], "C")) %>% 
   mean()
#> [1] 0.2236076

Esto nos da una aproximación de \(P(E|F)\). Nótese que si \(F\) es un evento con probabilidad baja, entonces será necesario correr más veces el experimento, pues el número de veces que ocurre \(F\) es relativamente bajo.

Podemos definir también la probabilidad condicional en general, para cualquier probabilidad \(P\) no necesariamente resultante de un modelo equiprobable:

Ejercicio

Se divide al azar una baraja de 52 de cartas en 4 pilas iguales. Calcula la probabilidad de cada pila tenga exactamente un as.

Podemos hacer el evento \(E_1\) que el as de corazones y el de diamantes están en diferentes pilas, \(E_2\) el evento de que el as de corazones, el de diamantes, y el de tréboles están en diferentes pilas, y finalmente \(E_3\) el evento de que todos los ases están en diferentes pilas. Nótese que buscamos la probabilidad \(P(E_3)\), pero será más fácil si escribimos:

\[P(E_3) = P(E_3E_2E_1) = P(E_1)P(E_2|E_1)P(E_3|E_1E_2)\]

Primero, el as de corazones está en alguna de las pilas. La probabilidad de el as de diamantes no esté en esa pila es \(P(E_1) = 1 - 12/51 = 39/51\) (¿por qué?), pues la pila que tiene el as de corazones tiene otras 12 cartas de las 51 disponibles. La probabilidad de que el as de diamantes sea una de esas 12 cartas es entonces 12/51.

Si se cumple \(E_1\), entonces el as de corazones y de diamantes están en pilas distintas. Entonces la probabilidad de que el as de tréboles caiga en una de esas dos pilas es \(24/50\), así que \[P(E_2|E_1) = 1 -24/50 = 25/50\]

Finalmente, si el as de corazones, diamantes y tréboles están en distintas pilas, entonces la probabilidad de que el de espadas caiga en la una de esas pilas es \(36/49\), de modo que \[P(E_3|E_2E_1) = 1 - 36/49 = 13/49\].

Usando la reglal producto, \[P(E_1E_2E_3) = \frac{39(26)(13)}{51(50)(49)} \approx 0.105\]

Ejemplo

Sacamos dos cartas de una bajara de 52 cartas. Vimos un ejemplo donde queríamos calcular la probabilidad de \(N_2=\) la segunda carta que sacamos es negra. Argumentamos usando espacio equiprobables que \(P(N_2) = 0.5\).

Si \(N_1 =\) la primera carta es negra, entonces la ley de probabilidad total explica por qué pasa esto tomando en cuenta el resultado de la primera extracción:

• Tenemos que \(P(N_2|N_1) = 25/51\) y \(P(N_1) = 1/2\)
• Además, \(P(N_2|N_1^c) = 26/51\) y \(P(N_1^c) = 1/2,\)

de forma que

\[P(N_2) = (25/51)(1/2) + (26/51)(1/2) = \frac{(25 + 26)}{51(2)} = /2\]

La ley de probabilidades totales consiste de las reglas de ponderación usuales que conocemos.

Ejemplo

En un país hay 20% de adultos de 20 años o menos, y 80% de adultos de 21 años o más. Entre los adultos de 20 años o menos, el 90% solo usa plataformas digitales para ver televisión. Entre los adultos de 21 años o más, el 15% solo usa plataformas digitales para ver televisión. Si tomamos un adulto al azar de esta población, ¿cuál es la probabilidad de que solo use digital para ver TV?

La respuesta es

\[ 0.90(0.20) + 0.15(0.80) = 0.3\]

Ejemplo

Supongamos que tenemos una baraja de 8 cartas, 4 negras y 4 blancas. Extraemos sucesivamente dos cartas a azar. Sea \(N_1=\) la primera carta es negra y \(B_2=\) la segunda carta es blanca. ¿Cuál es la probabilidad condicional de que la primera carta haya sido negra si nos dicen que la segunda fue blanca?

Aunque ya resolvimos problemas como este, parece confuso en un principio. Sin embargo, podemos calcular: \[P(N_1|B_2) = \frac{P(B_2|N_1)P(N_1)}{P(B_2)}\] Sabemos que \(P(B_2|N_1) = 4/7\), y que \(P(N_1)=1/2\) y \(P(B_2) = 1/2\), de modo que

\[P(N_1|B_2) = 4/7 > 1/2\]

es decir, si sabemos que la segunda carta fue blanca, eso hace más probable que la primera carta haya sido negra.

Ejemplo (parte 1)

Supongamos que una aseguradora cree que hay dos tipos de personas: unos con más riesgo de tener accidentes y otros con menos riesgo. Los datos muestran que una persona con riesgo alto tendrá un accidente en algún momento del año con probabilidad 0.04, y esta probabilidad baja a 0.01 para una persona de riesgo bajo. Si 10% de la población tiene riesgo alto, ¿cuál es la probabilidad de que un asegurado nuevo tenga un accidente un año después de comprar su póliza? (Nota: no sabemos si la persona nueva es de riesgo alto o bajo).

Puedes resolver son simulación, o usar la ley de probabilidad total. Si \(A\) es el evento de tener un accidente, \(R_A\) es el evento de que la persona es de riesgo alto y \(R_B\) es el evento que la persona es de riesgo bajo, entonces

\[P(A) = P(A|R_A)P(R_A) + P(A|R_B)P(R_B)\]

pues \(R_A\) y \(R_B\) cubren todas las posibilidades. Entonces

\[P(A) = 0.04(0.10) + (0.01)(0.90) = 0.004 + 0.009 = 0.013\] Es decir, su probabilidad es de 1.3% de tener una accidente en el primer año.

Ejemplo (parte 2)

Ahora vemos que un cliente tuvo un accidente en su primer año. ¿Cuál es la probabilidad de que sea un cliente de riesgo alto?

La pregunta es de probabilidad condicional, porque ya tenemos información. Queremos calcular

\[P(R_A| A)\]

Si usamos la regla de Bayes obtenemos

\[P(R_A|A)= \frac{P(A|R_A)P(R_A)}{P(A)}\]

Sustituimos los datos que tenemos

\[P(R_A|A)= \frac{0.04(0.10)}{P(A)}\] y \(P(A)\) que calculamos en el ejercicio anterior:

\[P(R_A|A)= \frac{0.04(0.10)}{0.013} \approx 0.3077\]

De manera que al principio la probabilidad no condicionada de ser de alto riesgo era de 10%, y cuando tiene un accidente esta probabilidad se triplica.

Ejemplo

Supongamos que en un concurso de TV tenemos tres puertas y debemos escoger una. Atrás de una de ellas hay un premio, y no hay nada detrás de las otras dos. Escogemos una de las puertas.

Ahora el conductor (que sabe dónde está el premio), abre una puerta vacía, y nos pregunta si queremos cambiar o no de puerta. ¿Cuál es la mejor estrategia, cambiar o quedarnos con la que escogimos al principio?

Veamos la estrategia de quedarnos con la puerta que escogimos. Sin perdida de generalidad, suponemos que escogemos la puerta 1.

Ahora observamos que el conductor abre la puerta 2.

Sea \(E_1=\) el premio está en la puerta 1, y \(A_2=\) el conductor abre la puerta 2, donde no hay premio. Por la regla de bayes:

\[P(E_1 | A_2) = \frac{P(A_2 | E_1) P(E_1)}{P(A_2)}\]

Sabemos que \(P(E_1)= 1/3\), y que \(P(A_2|E_1)=1/2\) (pues el conductor pudo abrir la puerta 2 o 3). Adicionalmente

\[P(A_2) = P(A_2|E_1)P(E_1) + P(A_2|E_1^c)P(E_1^c) = (1/2)(1/3) + (1/2) (2/3) = 1/2\]

pues \(P(A_2|E_1^c) = 1/2\), es decir, si el premio no está en la puerta 1, la probabilidad de que abrir la puerta 2 es la probabilidad de que el premio esté en la puerta 2 dado que no está en la puerta 1

Sustituyendo,

\[P(E_1 | A_2) = \frac{(1/2)(1/3)}{1/2} = 1/3\]

Asi que si cambiamos, la probabilidad de ganar es de 2/3.

Simula para verificar tus resultados.

Ejemplos

1. Muestra que si sabemos que sacar un corazón de una baraja de 52 cartas es independiente de sacar un as
2. Muestra que en nuestro modelo equiprobable de dos volados, el resultado de un volado es independiente del resultado del otro.
3. El evento “la suma de la tirada de dos dados” es independiente del resultado del primer dado.

Ejercicio

Hacemos un número indefinido de pruebas independientes, y cada una de ellas puede resultar en éxito con probabilidad \(p\) y fracaso con probabilidad \(1-p\). Calcula la probabilidad de que 1) al menos un éxito suceda en la primeras 20 pruebas. 2) todas las pruebas sean éxito y 3) exactamente 5 de las 2o pruebas sean éxito.

Ejemplo: número de seises

Construye un modelo para el número de seises en dos tiradas de dados. Escribimos \(X\)= número de seises que obtenemos en dos tiradas de dado.

Los resultados posibles son 0, 1 y 2 seises. Para calcular la probablidad de tirar dos seises hacemos:

\[P(X=2) = P(S_2\, y \, S_1) = P(S_2 | S_1) P(S_1).\] ahora, si suponemos que el primer tiro no afecta de ninguna manera el resultado del segundo tiro, entonces \(P(S_2|S_1) = P(S_2)\), y la fórmula es

\[P(X=2) = P(S_2\, y \, S_1) = P(S_2) P(S_1) = (1/6)(1/6) = 1/36.\] Usando el mismo argumento podemos calcular de la probabilidad de obtener ningún seis es

\[P(X=0) = P(S_2^c \, y \, S_1^c ) = P(S_2^c)P(S_1^c) = (5/6)(5/6) = 25/36\]

La probabilidad restante se puede calcularse directamente, o notar que como 0, 1 y 2 son los únicos posibles resultados, entonces

\[P(X=1) = 1- P(X=0) - P(X=2) = 1 - 25/36 - 1/36 = 10/36\]

Checa tus resultados usando simulación.

Observación: en muchos casos, la independencia se construye como un supuesto para construir modelos más complejos, cuando este supuesto es adecuado. Un ejemplo es cuando tomamos muestras de una población: si tomamos cada muestra independientemente de las otras, analizar los resultados es mucho más fácil que cuando hay esquemas complejos de dependencias entre los datos que xtraemos.

Ejemplo 1

Supongamos que una máquina produce dos tipos de bolsas de 25 galletas: la mitad de las veces produce una bolsa con 5 galletas de avena y 20 de chispas de chocolate, y la otra mitad produce galletas con 23 galletas de avena y 2 de chispas de chocolate.

Tomamos una bolsa, y no sabemos qué tipo de bolsa es (parámetro desconocido). Extraemos al azar una de las galletas, y es de chispas de chocolate (observación).

Por máxima verosimilitud, inferimos que la bolsa que estamos considerando tiene 5 galletas de avena. Esto es porque es más probable observar una galleta de chispas en las bolsas que contienen 5 galletas de avena que en las bolsas que contienen 23 galletas de avena. Podemos cuantificar la probabilidad que “acertemos” en nuestra inferencia.

Cómo se aprecia en el ejemplo anterior, el esquema general es:

1. Existe un proceso del que podemos obtener observaciones de algún sistema o población real.
2. Tenemos un modelo probabilístico que dice cómo se producen esas observaciones a partir del sistema o población real.
3. Usualmente este modelo tiene algunas cantidades que no conocemos, que rigen el proceso y cómo se relaciona el proceso con las observaciones.

Nuestro propósito es:

4. Extraemos observaciones del proceso:
   \[x_1, x_2, \ldots, x_n\]
5. Queremos aprender de los parámetros desconocidos del proceso para calcular cantidades de interés acerca del sistema o población real

En principio, los modelos que consideramos pueden ser complicados y tener varias partes o parámetros. Veamos primero un ejemplo clásico con un solo parámetro, y cómo lo resolveríamos usando máxima verosimilitud.

Nota: Cuando decimos muestra en general nos referimos a observaciones independientes obtenidas del mismo proceso (ver la sección anterior para ver qué significa que sea independientes). Este esquema es un supuesto que simplifica mucho los cálculos, como discutimos antes. Muchas veces este supuesto sale del diseño de la muestra o del estudio, pero en todo caso es importante considerar si es razonable o no para nuestro problema particular.

Ejemplo

Supongamos que queremos saber qué proporción de registros de una base de datos tiene algún error menor de captura. No podemos revisar todos los registros, así que tomamos una muestra de 8 registros, escogiendo uno por uno al azar de manera independiente. Revisamos los 8 registros, y obtenemos los siguientes datos:

\[x_1 = 0, x_2 = 1, x_3 = 0, x_4 = 0, x_5 =1, x_6 =0, x_7 =0, x_8 =0\]

donde 1 indica un error menor. Encontramos dos errores menores. ¿Cómo estimamos el número de registros con errores leves en la base de datos?

Ya sabemos una respuesta razonable para nuestro estimador puntual, que sería \(\hat{p}=2/8=0.25\). Veamos cómo se obtendría por máxima verosimilitud.

Según el proceso con el que se construyó la muestra, debemos dar una probabilidad de observar los 2 errores en 8 registros. Supongamos que en realidad existe una proporción \(p\) de que un registro tenga un error. Entonces calculamos

Probabilidad de observar la muestra:

\[P(X_1 = 0, X_2 = 1, X_3 = 0, X_4 = 0, X_5 =1, X_6 =0, X_7 =0, X_8 =0)\]

es igual a

\[P(X_1 = 0)P(X_2 = 1)P(X_3 = 0)P( X_4 = 0)P(X_5 =1)P(X_6 =0)P(X_7 =0)P(X_8 =0)\] pues la probabilidad de que cada observación sea 0 o 1 no depende de las observaciones restantes (la muestra se extrajo de manera independiente).

Esta ultima cantidad tiene un parámetro que no conocemos: la proporcion \(p\) de registros con errores. Así que lo denotamos como una cantidad desconocida \(p\). Nótese entonces que \(P(X_2=1) = p\), \(P(X_3=0) = 1-p\) y así sucesivamente, así que la cantidad de arriba es igual a

\[p(1-p)p(1-p)(1-p)p(1-p)(1-p)(1-p) \]

que se simplifica a

\[ L(p) = p^2(1-p)^6\]

Ahora la idea es encontrar la p que maximiza la probabilidad de lo que observamos. En este caso se puede hacer con cálculo, pero vamos a ver una gráfica de esta función y cómo resolverla de manera numérica.

verosimilitud <- function(p){
  p^2 * (1-p)^6
}
dat_verosim <- tibble(x = seq(0,1, 0.001)) %>% mutate(prob = map_dbl(x, verosimilitud))
ggplot(dat_verosim, aes(x = x, y = prob)) + geom_line() +
  geom_vline(xintercept = 0.25, color = "red") +
  xlab("p")

Nótese que esta gráfica:

• Depende de los datos, que pensamos fijos.
• Cuando cambiamos la \(p\), la probabilidad de observar la muestra cambia. Nos interesa ver las regiones donde la probabilidad es relativamente alta.
• El máximo está en 0.25.
• Así que el estimador de máxima verosimilitud es \(\hat{p} = 0.25\)

Obsérvese que para hacer esto usamos:

• Un modelo teórico de cómo es la población con parámetros y
• Información de como se extrajo la muestra,

y resolvimos el problema de estimación convirtiéndolo en uno de optimización.

Probamos esta idea con un proceso más complejo:

Aspectos numéricos

Cuando calculamos la verosimilitud arriba, nótese que estamos multiplicando números que pueden ser muy chicos (por ejemplo \(p^6\), etc). Esto puede producir desbordes numéricos fácilmente. Por ejemplo para un tamaño de muestra de 1000, podríamos tener que calcular

p <- 0.1
proba <- (p ^ 800)*(1-p)^200
proba
#> [1] 0

En estos casos, es mejor hacer los cálculos en la escala logarítmica. El logaritmo convierte productos en sumas, y baja exponentes multiplicando. Si calculamos en escala logaritmica la cantidad de arriba, no tenemos problema:

log_proba <- 800 * log(p) + 200 * log(1-p)
log_proba
#> [1] -1863.14

Ahora notemos que

• Maximizar la verosimilitud es lo mismo que maximizar la log-verosimilitud, pues el logaritmo es una función creciente. Si \(x_{max}\) es el máximo de \(f\), tenemos que \(f(x_{max})>f(x)\) para cualquier \(x\), entonces tomando logaritmo, \[log(f(x_{max}))>log(f(x))\] para cualquier \(x\), pues el logaritmo respeta la desigualdad por ser creciente.
• Usualmente usamos la log verosimilitud para encontrar estimador de máxima verosimilitud
• Hay razónes teóricas y de interpretación por las que también es conveniente hacer esto.