7.1.1 Ejemplo: ruleta sesgada

Ahora supongamos que nuestra ruleta no está del todo balanceada. Por ejemplo, podría ser estuviera colgada en una pared, y al girar la flecha es un poco más probable que la flecha apunte hacia el piso en lugar de hacia el cielo.

En este caso, si la dirección hacia arriba es 90 grados y hacia abajo es 270 grados, quisiéramos por ejemplo que \[P(260 < X <280) < P(80 < X < 100)\] Y nótese que debe ser posible asignar probabilidades a cualquier sector de la ruleta con el nuevo modelo que propongamos. ¿Cómo podríamos modificar nuestra asignación de probabilidades?

Una de las maneras más fáciles es pensando que nuestra probabilidad se obtiene integrando una funcion constante:

• Si \([a,b]\) es un sector de la ruleta con \(a<b\), podríamos poner \[P(X\in [a,b]) = \int_a^b \frac{1}{360} \,dx = \frac{b-a}{360}\] De forma que si \(f(x)= 1/360\) para valores \(0 \leq x < 360\), nuestra probabilidad se escribe como la integral

\[P(X\in [a,b]) = \int_a^b f(x) \,dx \]

• En este caso, probabilidad es área bajo la curva \(f(x)=1/360\) que se calcula integrando sobre el intervalo de interés

Para generalizar la idea es la siguiente:

• Usamos la fórmula anterior, pero modificamos o perturbamos la función \(f(x) = 1/360\) para que \(f\) sea un poco más alta alrededor de 270 grados (abajo), y un poco más baja alrededor de 90 grados (arriba).
• Lo único que necesitamos es que \(f(x)\) no puede tomar valores negativos (por que si no obtendríamos probabilidades negativas en algunos sectores), y la integral sobre la ruleta completa debe ser uno:

\[P(X\in [0, 360]) = \int_0^{360} f(x)\,dx = 1\] Podríamos utilizar por ejemplo:

f_dens <- function(x){
    x_rad <- 2 * pi * x / 360
    (1/360) +  0.0002 * cos(x_rad - 3 * pi  / 2)
}
graf_1_tbl <- tibble(angulo = seq(0, 360, 1), tipo = "uniforme",
                   f = 1 / 360) 
graf_2_tbl <- tibble(angulo = seq(0, 360, 1), tipo = "colgada") %>% 
    mutate(f = f_dens(angulo))
graf_tbl <- bind_rows(graf_1_tbl, graf_2_tbl)
ggplot(graf_tbl, aes(x = angulo, y = f, colour = tipo)) +
    geom_line() +
    ylim(c(0, 0.003)) + facet_wrap(~tipo, nrow = 1)

El cálculo se hace ahora con área bajo la curva. Para calcular la probabilidad \[P(X\in [50, 130]),\] integramos la función \(f\) correspondiente, que corresponde a calcular área bajo la curva:

ggplot(graf_tbl, aes(x = angulo, y = f, colour = tipo)) +
    geom_line() +
    ylim(c(0, 0.003)) + facet_wrap(~tipo, nrow = 1) +
    geom_area(aes(x = ifelse(angulo > 50 & angulo < 130, angulo, 0)), 
              fill="salmon", alpha = 0.5)

Y ahora vemos que para la versión perturbada, más de la probabilidad se concentra alrededor de 270 grados que alrededor de 90. Por las propiedades de la integral, todas las propiedades usuales de probabilidad se cumplen.

7.5.1 Ejemplo: distribución uniforme

Simularemos 500 datos uniformes en \([0, 10]\):

x_sim_u <- runif(500, 0, 10)

Podríamos calcular algunos cuantiles empíricos:

quantile(x_sim_u, c(0.10, 0.50, 0.90))
#>      10%      50%      90% 
#> 1.223266 4.998312 8.946309

Por el ejercicio anterior sabemos cuáles son los cuantiles teóricos correspondientes a una uniforme en \([0,10]\). Podemos calcularlos también como sigue:

qunif(c(0.10, 0.5, 0.90), 0, 10)
#> [1] 1 5 9

Y vemos que son muy similares los cuantiles empíricos y teóricos. Una mejor manera de considerar esta similitud es graficando todos los cuantiles empíricos y comparándolos con los teóricos:

ggplot(tibble(x = x_sim_u), aes(sample = x)) +
    geom_abline(colour = "red") +
    geom_qq(distribution = stats::qunif, dparams = list(min = 0, max = 10)) +
    xlab("Cuantiles Teóricos U(0,10)") + ylab("Cuantiles de datos")

Y vemos que la forma de las dos distribuciones es muy similar: los cuantiles empíricos son muy similares a los teóricos. Existen algunas fluctuaciones debidas al muestreo aleatorio.

7.5.2 Ejemplo: distribución triangular

Repetimos para la distribución triangular. Los cuantiles que calculamos arriba son:

qtri(c(0.25, 0.75), a = 0, b = 4)
#> [1] 1.414214 2.585786

x_sim_tri <- rtri(500, 0, 4)
ggplot(tibble(x = x_sim_tri), aes(sample = x)) +
    geom_abline(colour = "red") +
    geom_qq(distribution = qtri, dparams = list(a = 0, b = 4)) +
    xlab("Cuantiles Teóricos triangular(0,4)") + ylab("Cuantiles de datos")

Nótese que otra vez, los cuantiles teóricos se alinean bien con los teóricos.

7.7.1 Media teórica y empírica

Sabemos que la media de un conjunto de datos \(x_1,x_2,\ldots,x_n\) está dada por

\[\bar{x} = \frac{x_1 + x_2 + \cdots x_n}{n}.\] Ahora consideremos qué pasa con esta cantidad cuando las \(x_i's\) son observaciones independientes de una distribución con densidad teórica \(f(x)\). Utilizaremos simulaciones de la distribución triangular que vimos arriba

datos <- tibble(x = rtri(5000, 0, 4))
ggplot(datos, aes(x = x)) + stat_bin(breaks = seq(0, 4, 0.25))

Podemos aproximar la media de estos datos promediando los valores iniciales de los intervalos de corte ponderado por el número de casos. Lo hacemos así:

particion <- seq(0, 4, 0.25)
iniciales <- head(particion, -1) # quitar último
agrupados_cubeta <- datos %>% 
  mutate(inicial = cut(x, breaks = particion, labels = iniciales)) %>%  
  mutate(inicial = as.numeric(as.character(inicial))) %>% 
  group_by(inicial) %>%
  count() %>% 
  ungroup() 
agrupados_cubeta
#> # A tibble: 16 x 2
#>    inicial     n
#>      <dbl> <int>
#>  1    0       44
#>  2    0.25   110
#>  3    0.5    199
#>  4    0.75   256
#>  5    1      330
#>  6    1.25   441
#>  7    1.5    494
#>  8    1.75   569
#>  9    2      567
#> 10    2.25   535
#> 11    2.5    450
#> 12    2.75   377
#> 13    3      282
#> 14    3.25   178
#> 15    3.5    130
#> 16    3.75    38

Y calculamos la aproximación a la media como sigue:

agrupados_cubeta %>% 
  summarise(media_aprox = sum(inicial * n) / sum(n))
#> # A tibble: 1 x 1
#>   media_aprox
#>         <dbl>
#> 1        1.89

Apliquemos esta idea cuando tenemos una densidad \(f(x)\). Dividimos el rango de la densidad en cubetas, y aproximamos la densidad por rangos, por ejemplo:

particion <- seq(0, 4, 0.25)
valor <- dtri(particion, 0, 4)
approx_tbl <- tibble(x = particion, densidad = valor)
densidad_tri <- tibble(x = seq(0, 4, 0.01)) %>% 
  mutate(densidad = dtri(x, 0, 4))
ggplot(densidad_tri) + 
  ylab('f(x)') + 
  geom_line(aes(x = x, y= densidad), alpha = 0.8) +
  geom_step(data = approx_tbl, aes(x = x, y = densidad), colour = "red") +
  theme_minimal() 

Y repetimos el mismo proceso: ponderamos los valores iniciales por la altura de la densidad:

approx_tbl %>% 
  summarise(media_approx = sum(x * densidad) / sum(densidad))
#> # A tibble: 1 x 1
#>   media_approx
#>          <dbl>
#> 1            2

Y esta es una aproximación a la media de esta distribución.

Nótese que la cantidad que estamos calculando es

\[\sum_i x_i f(x_i) \Delta\]

donde \(\Delta\) es igual a

approx_tbl %>% summarise(suma_densidad = 1 / sum(densidad))
#> # A tibble: 1 x 1
#>   suma_densidad
#>           <dbl>
#> 1          0.25

que es el ancho del intervalo de las particiones. Recordamos por cálculo que esta es una la aproximación a la siguiente integral:

\[\sum_i x_i f(x_i) \Delta \approx \int xf(x)\,dx\]

De modo que para pasar de media de los datos \(\bar{x}\) a la media de \(\mu_f\) de una distribución la equivalencia es:

\[ \bar{x} = \frac{1}{n}\sum x_i \longleftrightarrow \int xf(x) = \mu \] donde

• Usamos la densidad en lugar de frecuencias relativas para ponderar
• Usamos integral en lugar de suma

Es decir, cuando tenemos una densidad teórica continua, es necesario integrar en lugar de sumar para calcular su media.

7.7.2 Varianza y desviación estándar

Otra cantidad importante es la varianza de una muestra. Es una medida de dispersión, y se calcula como

\[ \hat{\sigma}^2 = \frac{1}{n}\sum_{i=1}^n (x_i - \bar{x})^2 \]

Nótese que cuando los datos están altamente concentrados alrededor del valor de la media, la varianza es chica, y cuando hay más dispersión alrededor de la media, la varianza es grande.

La desviación estándar es la raíz cuadrada de esta cantidad:

\[ \hat{\sigma} = \sqrt{\frac{1}{n}\sum_{i=1}^n (x_i - \bar{x})^2} \] que tiene la ventaja de que está en las mismas unidades que la variable original.

Siguiendo el mismo patrón que arriba, tenemos que integrar en lugar de sumar, y ponderar por la densidad en lugar de la frecuencia: El equivalente en una distribución teórica es también una integral, y la varianza está definida por

\[\sigma^2 = \int (x - \mu)^2 f(x) \,dx\]

La desviación estándar \(\sigma\) es la raíz cuadrada de esta cantidad.

\[\sigma = \sqrt{\int (x - \mu)^2 f(x)\,dx}\]

En resumen la densidad indica la frecuencia relativa de datos que esperamos observar alrededor de cada punto, y

• Resúmenes de una distribución teórica (como cuantiles, media, varianza, etc.) se calculan integrando ponderado por la densidad.
• Resúmenes de una distribución empírica se calculan sumando ponderando por la frecuencia relativa.

7.7.3 ¿Cuándo usar media y desviación estándar?

Algunos modelos probabilísticos son más fáciles de tratar analíticamente usando media y varianza. En esos casos, conviene usar estas medidas. Esto es especialmente cierto cuando las distribuciones son simétricas y no tienen colas muy largas.

En cuanto a datos observados, conviene usar media y varianza cuando pretendemos modelarlos con densidades como las del párrafo anterior. En este caso, la interpretación de estos valores se hace a través de la forma de la densidad. Es importante checar (por ejemplo usando gráficas de cuantiles) que esas densidades describen apropiadamente a los datos.

7.8.1 Promedio de variables

Consideremos que tiramos 40 dados justos de 6 lados, y consideramos su promedio \(\bar{X}\) como resultado de nuestro experimento aleatorio. ¿Cómo se ve la distribución de probabilidades de esta variable \(\bar{X}\)?

Comenzamos haciendo simulacion:

simular_bolsa <- function(num_dados = 40){
    dados <- sample(1:6, num_dados, replace = TRUE)
    media <- mean(dados)
    media
}
simular_bolsa()
#> [1] 3.85

Veamos cómo se ven los resultados si hacemos este experimento un número grande de veces:

set.seed(23)
sims_dados <- map_dbl(1:10000, ~ simular_bolsa())
ggplot(tibble(resultado = sims_dados), 
       aes(x = resultado)) +
    geom_histogram(binwidth = 0.1)

Y notamos una forma de “campana” interesante. Esto se explica porque típicamente tendremos tiros bajos y altos, de modo que muchos resultados de este experimento se concentran alrededor de la media de un dado (3.5). Además, existen fluctuaciones aleatorios, y a veces tenemos un poco más de tiros altos o de tiros bajos, de forma que existe dispersión alrededor de 3.5.

Sin embargo, estas desviaciones de 3.5 no pueden ser muy grandes: por ejemplo, para tener un promedio de 1, todas las tiradas de los 40 dados tendrían que dar 1, y eso es muy poco probable. Igualmente, para que el promedio fuera cercano a 6, la gran mayoría de los 40 dados deberían de dar 6, lo cual otra vez es muy poco probable. Esto explica al menos la forma general de la forma de las colas derecha e izquierda de esta distribución.

Los dados podrían ser diferentes (por ejemplo, un poco cargados a 1 o 6, o más cargados a valores centrales), y los argumentos de arriba también se cumplirían. Lo que es sorprendente es que, independientemente de cómo sean las particularidades de los dados, la forma analítica de esta distribución que acabamos de mostrar es la misma.

Esta forma está descrita por la densidad normal estándar:

\[f(x) = \frac{1}{\sqrt{2\pi}} e^{-\frac{x^2}{2}}\]

cuya gráfica presentamos a continuación:

tibble(x = seq(-3.5, 3.5, 0.01)) %>% 
    mutate(f = (1/sqrt(2*pi)) * exp(-x^2 / 2)) %>% 
ggplot(aes(x = x, y = f)) + geom_line()

A una variable \(Z\) que tiene esta densidad le llamamos una variable con distribución normal estándar. Comparemos cuantiles en nuestro ejemplo:

ggplot(tibble(resultado = sims_dados),
       aes(sample = resultado)) +
    geom_qq(distribution = stats::qnorm) +
    xlab("Normal estándar teórica") +
    ylab("Promedio de 40 dados")

Y notamos que los cuantiles no corresponden, pero el espaciamiento entre los cuantiles de los datos y los teóricos de la normal estándar es el mismo. Quiere decir que estas dos distribuciones tienen la misma forma, aunque estén escaladas y centradas en distintos valores.

Probemos con promedios de 20 observaciones triangulares en \((0,1)\) por ejemplo. El resultado es el mismo:

set.seed(23)
sims_tri <- map_dbl(1:10000, ~ mean(rtri(20, 0 ,1)))
ggplot(tibble(resultado = sims_tri), 
       aes(x = resultado)) +
    geom_histogram(binwidth = 0.01)

ggplot(tibble(resultado = sims_tri),
       aes(sample = resultado)) +
    geom_qq(distribution = stats::qnorm) +
    xlab("Normal estándar teórica") +
    ylab("Promedio de 20 triangulares (0,1)")

Otra vez, la forma general es la misma, aún cuando los datos están centrados y escalados de manera distinta.

7.12.1 Variables aleatorias exponenciales

Supongamos que estamos modelando tiempos a la ocurrencia \(X\) de un evento (por ejemplo en análisis de supervivencia en estudios clínicos). Esta es una variable que toma valores en los números positivos. ¿Cómo podría ser su forma?

Consideremos por ejemplo que el tiempo de espera no tienen memoria. Es decir: si hay esperamos un periodo de \(t\) días por ejemplo, la distribución del tiempo restante que tenemos que esperar no depende de \(s\). En términos de probabilidad, podríamos escribir: \[P(X > s + t | X > t) = P(X > s)\] Esto se lee: dado que el evento ocurre en más de \(t\) días, la probabilidad de que tarde al menos otros \(s\) días no dependen de \(t\). A la función \(S(t) = P(X>t)\) muchas veces se le llama la función de supervivencia.

Este es un modelo base útil, que después puede extenderse a procesos donde los eventos ocurren aceleradamente (envejecen), o donde los eventos dan evidencia de robustez (los que han sobrevivido hasta cierto tiempo se espera que duren más que lo que inicialmente), o quizá una combinación de los dos dependiendo de el valor de \(t\).

Con el supuesto de falta de memoria, la ecuación de arriba se cumple, y entonces (¿por qué?): \[P(X > s + t) = P(X > s) P(X > t),\] así que \[\frac{1}{s}(P(X > s+ t) - P(X>t)) = \frac{1}{s}(P(X>s) - 1)P(X>t).\] Que podemos reescribir usando la densidad \(f(x)\) como \[\frac{1}{s}\int_t^{t+s} f(x)\,dx = P(X>t)\frac{1}{s}\int_0^s f(x)\,dx\] Conforme \(s\) se hace más chica, el lado izquierdo converge a \(f(t)\). El lado derecho, por otra parte, converge a \(f(0)\), y obtenemos \[f(t) = f(0)\int_t^\infty f(x)\,dx,\] y ahora derivamos de ambos lados para obtener \[f'(t) = -f(0)f(t)\] La única función que satisface esta propiedad (su derivada es proporcional a ella misma) es \[f(t) = \lambda e^{-\lambda t}\] Puedes checar que efectivamente esta densidad cumple que la ecuación anterior. A esta densidad le llamamos la densidad exponencial con tasa \(\lambda\).

lambda <- 1
sim_exp <- rexp(1000, rate = lambda)
ggplot(tibble(x = sim_exp), aes(x = x)) +
    geom_histogram(boundary = 0)
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Ejercicio: demuestra que una variable aleatoria exponencial con tasa \(\lambda\) tiene media \(E(X) = 1/\lambda\) y $Var(X) = 1 / ^2

El parámetro \(\lambda\) se llama tasa por la interpretación de tiempo de espera que mostramos arriba. Supongamos que \(\lambda = 10\). Eso quiere decir que esperamos observar el evento en \(1/\lambda = 1/10\) minutos, por ejemplo, o lo que es lo mismo, a una tasa de \(\lambda = 10\) eventos por minuto.

Nótese finalmente que todas las variables exponenciales son escalamientos de una variable exponencial con \(\lambda = 1\). Esto es porque si \(Y=kX\), y \(X\) es exponencial con \(\lambda= 1\) entonces:

\[P(Y>t)=P(kX > t) = P(X > t/k) = \int_{t/k}^\infty e^{-x}\,dx = e^{-t/k}\] De modo que derivando, encontramos que la densidad de \(Y\) es \(\frac{1}{k} e^{-t/k}\), que es una exponencial con tasa \(\lambda = 1/k\).

7.12.2 Ejemplo: exponencial

Supongamos que un tipo de focos tienen tiempos de vida exponencial con una vida media de 10 años. ¿Cuál es la probabilidad de que un foco dure más de 15 años? Si tenemos un foco que ya vivió 10 años, cuál es la probabilidad de que dure otros 15 años.

Tenemos que la vida de un foco es una variable \(X\) exponencial con parámetro \(\lambda = 1/10\). La probabilidad de que dure más de 15 años es entonces \[P(X>15) = \int_{15}^\infty \frac{1}{10} e^{x/10}\,dx\]

Podemos calcular a mano, o usar rutinas usuales de R:

pexp(15, rate = 1/10, lower.tail = FALSE)
#> [1] 0.2231302

7.12.3 Variables aleatorias gamma

Esta es otra familia que extiende la familia de distribuciones exponenciales. La forma analítica de una densidad gamma con parámetro de forma \(k>0\) y tasa \(\lambda\)

\[f(x) = C x^{k-1} e^{-\lambda x}\] donde la constante \(C\) de normalización depende de \(k\) y \(\lambda\).

Abajo vemos datos simulados de densidades Gamma con distintas combinaciones de parámetros:

params_tbl <- crossing(k = c(1, 2, 5, 10), lambda = c(1/4, 1/2))
sims_tbl <- params_tbl %>% 
    mutate(sims = map2(k, lambda, 
                       ~ rgamma(10000, shape = .x, rate = .y))) %>% 
    unnest(cols = sims)

ggplot(sims_tbl, aes(x = sims)) + 
    geom_histogram(boundary = 0, bins = 50) +
    facet_grid(k~lambda)

Nótese que cada columna es un rescalamiento de la otra, pero las densidades de los renglones tienen distinta forma. Puedes ver aquí parámetros como esperanza, varianza de esta estas distribuciones, junto con otras propiedades y aplicaciones.