Probabilidades y Distribuciones para Variables Aleatorias

Objetivos:

  • Familiarizarse con distintas distribuciones de probabilidades y con la manera de usarlas en R.
  • Simular procesos estocásticos (aleatorios).
  • Usar un bestiario de distribuciones. Con este documento se pueden familiarizar con distintas distribuciones de probabilidades (e.g. Poisson, Binomial, Normal, etc.)

Las distribuciones de probabilidad en general son ecuaciones (pueden ser tablas también) que definen la probabilidad de ocurrencia para cada posible resultado de un experimento o procedimiento de muestreo. R tiene funciones para todas las distribuciones de probabilidad estándares, y para cada una de estas distribuciones tenemos funciones para:

  • generar valores,
  • calcular probabilidades,
  • probabilidades acumuladas y
  • cuantiles

En R estas funciones comienzan con las letras r, d ,p y q respectivamente. Por ejemplo, para la distribución de Possion tenemos: rpois, dpois, ppois, y qpois. Vamos a ver qué hace cada una de estas funciones.

En muchos casos es útil poder generar valores aleatorios (muestras) de una distribución en particular. Al decir “de una distribución particular”, nos referimos no solo al tipo de distribución (e.g. Poisson) sino también a determinados valores de parámetros (e.g. \(\lambda = 1\)). Asumimos que estas muestras generadas por la computadora son una “muestra aleatoria” pero en realidad provienen de un generador de números aleatorios, por lo que es más correcto decir que son “pseudo-aleatorios”. Un aspecto importante, sobre todo pensando en la reproducibilidad de nuestro trabajo, es que si iniciamos al generador de números aleatorios con un valor determinado, la secuencia de números pseudo-aleatorios se va a repetir, y por lo tanto podemos reproducir exactamente una simulación estocástica. Existen muchos algoritmos para generar números pseudo-aleatorios, pero en general no nos metemos demasiado en estos detalles y confiamos en que R sabe lo que hace.

En R usamos la función set.seed para inicializar el generador de números aleatorios en un valor determinado. Por ejemplo set.seed(1234). El número que le pasamos a set.seed es arbitrario (\(12345\) en este caso) pero debe ser un número entero.

Distribuciones discretas

Simulamos una muestra de \(10\) valores de una distribución de Poisson con media \(1.2\). Imaginemos que son datos de conteo de semillas en \(10\) trampas de \(1 \text{m}^2\).

set.seed(12345)
rpois(n = 10, lambda = 1.2)
##  [1] 2 2 2 3 1 0 1 1 2 4

Nuestros scripts pueden ser más fáciles de leer y modificar si definimos variables por fuera de las funciones. Por ejemplo:

n <- 10
lambda <- 1.2
y <- rpois(n = n, lambda = lambda)

Veamos la frecuencia del número de semillas por trampa simuladas en un histograma

op <- par(cex.main = 1.5, cex.lab = 1.5, font.lab = 2, cex.axis = 1.3, bty = "n", 
    las = 1)  #Estos argumentos en la función par son para definir características de nuestro plot. Si quieren ver en detalle de que se trata y que se puede hacer escriban en la consola de R ?par

hist(y, xlab = "Valores", ylab = "frecuencia", main = "")

par(op)

Una vez que tenemos “datos” como estos, podemos ver las frecuencias relativas y compararlas con la distribución teórica usando la función dpois. Para hacer una comparación visual, graficamos las frecuencias relativas “observadas” y le agregamos la probabilidad teórica de observar cada uno de esos valores. En lugar de usar un histograma, vamos a usar la función table para estimar la proporción de la muestra que cae en cada valor de \(y\) para luego graficar esos valores.

obsprobs <- table(y)/n  # frecuencias relativas para cada valor en los datos

op <- par(cex.main = 1.5, cex.lab = 1.5, font.lab = 2, cex.axis = 1.3, bty = "n", 
    las = 1)
plot(obsprobs, xlab = "Valores", ylab = "Probabilidad")
tprobs <- dpois(x = 0:max(y), lambda = lambda)  # acá calculamos las probabilidades teóricas
points(0:max(y), tprobs, pch = 16, col = 2)

par(op)

Otra función útil es la Función de Distribución, que nos da la probabilidad de encontrar un valor menor o igual a un valor dado. Por ejemplo, la probabilidad de obtener \(y<=3\) con \(\lambda=\) 1.2 es:

ppois(3, lambda = lambda)
## [1] 0.966231

Si queremos saber la probabilidad de obtener un valor mayor a \(3\) hacemos

1 - ppois(3, lambda = lambda)
## [1] 0.03376897

Para ver la probabilidad de un valor en particular, por ejemplo \(y=3\), podemos usar ppois ó dpois

ppois(3, lambda = lambda) - ppois(2, lambda = lambda)
## [1] 0.08674393
dpois(3, lambda = lambda)
## [1] 0.08674393

Por último, la función cuantil qpois es la inversa de la función de distribución, y para un valor de probabilidad acumulada nos devuelve el valor de la variable

qpois(0.95, lambda = lambda)  # esto nos daria el valor del cuantil 0.95
## [1] 3

La función qpois sirve también para calcular intervalos que contienen un porcentaje de los valores de la distribución. Por ejemplo, el \(95\)% de los valores de una distribución están entre el cuantil \(0.025\) y el cuantil \(0.975\). El equivalente empírico (i.e. de los datos) es la función quantile.

qpois(c(0.025, 0.975), lambda)
## [1] 0 4
quantile(y, probs = c(0.025, 0.975))
##  2.5% 97.5% 
## 0.000 2.775

Ejercicios:

  1. Asumiendo una distribución de Poisson:
    • ¿Cuál es la probabilidad teórica de obtener \(y=0\) para \(\lambda = 1\)?
    • ¿Cuál es la probabilidad teórica de \(y > 2\) para \(\lambda = 1\)?
    • Comparar la probabilidad de \(y=0\) para \(\lambda = 1\) de la distribución teórica con la obtenida empíricamente con una muestra de \(10\), \(100\), \(1000\) y \(10000\) observaciones generadas usando rpois.
    • Ahora con las mismas muestras, calcular los intervalos de \(95\)% teóricos y empíricos.
  2. ¿Cómo cambia la forma de la distribución de Poisson a medida que cambia \(\lambda\)?
  3. Otras distribuciones discretas. Ver otras distribuciones del bestiario. ¿Qué distribución les parece más adecuada para los siguientes tipos de datos? (Justificar la respuesta!):
    • Las hembras de ñu pueden tener una sola cría por año, pero no todas las hembras se reproducen todos los años. Se registró en una muestra aleatoria de \(100\) hembras de una población si habían tenido o no una cría ese año.
    • Se está muestreando la lluvia de semillas de coihue en una zona que se incendió recientemente. El coihue se dispersa por el viento y queremos ver cuántas semillas llegan a la zona incendiada para estimar el potencial de recuperación del bosque luego del incendio. Para ello se diseñó un muestreo en el que se colocaron trampas de semillas de \(0.25 \text{m}^2\) dentro de una grilla y se registró el número de semillas que cayó en cada una de las trampas durante la temporada de dispersión.
    • En la manzana, las flores tienen \(20\) óvulos en total. En una muestra aleatoria compuesta de \(1000\) flores de manzana seleccionadas dentro de un cultivo, se determinó la cantidad de óvulos fecundados por flor.
    • En un cultivo de frambuesa se seleccionaron al azar \(100\) flores, en las cuales se registró el número de polinizadores que visitaron cada flor en un período de \(5\) minutos.

Distribuciones continuas

Podemos hacer algo parecido a lo que hicimos con la distribución de Poisson para distribuciones continuas. Una diferencia importante entre las distribuciones discretas y contiunas es que para variables continuas tenemos infinitos valores posibles. Entonces ya no hablamos de la probabilidad de obtener un valor en particular sino de la densidad de probabilidad alrededor de un valor determinado. Esta densidad de probabilidad puede ser mayor que \(1\). En general, no prestamos demasiada atención a los valores particulares de densidad de probabilidad sino a comparaciones relativas.

Veamos por ejemplo la distribución Normal. Simulamos cambios de un año para otro en el área de uso (home range) de una especie de ratón. Esto quiere decir que el ratón puede: (1) aumentar el tamaño de su área de uso, y la diferencia observada será positiva; (2) reducir el tamaño del área y la diferencia será negativa; o bien (3) mantener estable su área de uso, por lo que la diferencia será cercana a 0. Suponiendo que observamos \(60\) individuos:

set.seed(1234)
n <- 60
mu <- 0
sigma <- 1
y <- rnorm(n, mean = mu, sd = sigma)

Ahora graficamos los datos simulados con un histograma y le agregamos una línea (en rojo) con la estimación de densidades calculadas a partir de los datos simulados. Esto es análogo al cálculo de las frecuencias relativas en las distribuciones discretas. Luego agregamos la densidad de probabilidades teórica (en negro).

op <- par(cex.lab = 1.5, font.lab = 2, cex.axis = 1.3, las = 1, bty = "n")
hist(y, breaks = 20, freq = FALSE, main = "")
lines(density(y), lwd = 3, lty = 3, col = "red")
curve(dnorm(x, mean = mu, sd = sigma), lwd = 3, add = TRUE, xlim = c(-4, 4))
Histograma de datos simulados, densidad empírica (en negro) y teórica (en rojo) de una distribución Normal con $\mu=$ 0 y $\sigma=$ 1

Histograma de datos simulados, densidad empírica (en negro) y teórica (en rojo) de una distribución Normal con \(\mu=\) 0 y \(\sigma=\) 1 

par(op)
  • Visualmente, ¿Qué tan parecida es la distribución empírica con la teórica?
  • ¿Cómo cambia esa similitud si duplico el esfuerzo de muestreo?
  • ¿Y si logramos tener \(800\) observaciones?

Veamos ahora cómo obtener distintos valores de probabilidades a partir de la distribución teórica. Por ejemplo, la densidad de probabilidad de encontrar un cero:

dnorm(0, mean = mu, sd = sigma)
## [1] 0.3989423

¿Cuál es la probabilidad de \(y > 1.96\)?

pnorm(1.96, mean = mu, sd = sigma, lower.tail = FALSE)  # pnorm por defecto calcula la probabilidad de obtener un valor que sea mayor o igual. En este caso queremos la probabilidad de que sea menor por eso ponemos 'lower.tail = FALSE'  
## [1] 0.0249979

Ahora calculamos el intervalo del \(95\)%

qnorm(c(0.025, 0.975), mean = mu, sd = sigma)
## [1] -1.959964  1.959964

Ejercicios:

  1. ¿Cómo cambia la forma de la distribución cuando cambia \(\sigma\)?
  2. Otras distribuciones continuas. Ver otras distribuciones del bestiario para familiarizarse con las formas que éstas pueden tener y el tipo de datos que pueden llegar a representar. Qué distribución les parece más adecuada para los siguientes tipos de datos:
    • Luego de un incendio, se registró la proporción de copa de árboles que queda sana.
    • Para ver el efecto que tiene un determinado tratamiento de manejo ganadero sobre el peso del ganado, se registró el peso antes y después de aplicado dicho tratamiento en una muestra de 300 individuos.
    • En condiciones de laboratorio, se quiere registrar el tiempo que tarda un frugívoro en defecar las semillas de los frutos que ingiere.
    • Se colocaron collares GPS en ciervos colorados, y se quiere ver la distancia que recorren por día los individuos.

Extendiendo distribuciones

En muchos casos, un aspecto clave en el análisis de datos consiste en modelar cambios en los parámetros de una distribución en relación con alguna variable explicadora (covariable). Por ejemplo, veamos datos de peso (en gramos) de individuos de monito del monte.

Para cargar los datos, podemos leerlos desde un repositorio online (Opción \(1\)) o desde un archivo guardado en nuestro directorio de trabajo (Opción \(2\)). Para la opción \(2\), tenemos que setear el directorio de trabajo usando la función setwd(). Hay diferentes funciones que pueden utilizarse para cargar datos a R según el formato en el que tengamos nuestra base de datos. En este caso, como es un archivo de texto (.txt) vamos a usar la función read.table, y como nuestro archivo tiene los nombres de las variables en el encabezado, especificamos header = TRUE.

# Opción 1
datos <- read.table("https://sites.google.com/site/modelosydatos/peso.txt", 
    header = TRUE)

# Opción 2 datos <- read.table('peso.txt', header = TRUE)

op <- par(cex.lab = 1.5, font.lab = 2, cex.axis = 1.3, las = 1, bty = "n")
hist(datos$peso, 20, freq = FALSE, xlab = "Peso (gramos)", ylab = "Densidad", 
    main = "")
lines(density(datos$peso), lwd = 2)

par(op)

Vemos que los pesos de estos animales tienen una distribución sesgada. Pero quizás esto se debe a que estamos mezclando individuos hembras y machos. Si separamos por sexo, las distribuciones siguen siendo aparentemente sesgadas pero se puede ver que parte del sesgo original se debe a que hay cierto dimorfismo sexual. Es decir que la distribución que vemos de pesos en realidad es una mezcla de dos distribuciones, una para los pesos de los machos y otra para los pesos de las hembras.

op <- par(cex.lab = 1.5, font.lab = 2, cex.axis = 1.3, las = 1, bty = "n")
hist(datos$peso, 20, freq = FALSE, xlab = "Peso (gramos)", ylab = "Densidad", 
    main = "")
lines(density(datos$peso[datos$sexo == 1]), lwd = 2)
lines(density(datos$peso[datos$sexo == 0]), lwd = 2)

par(op)

Otra forma común de extender distribuciónes es agregando un “modelo” a algún parámetro de una distribución. En este caso, además del peso de los individuos, tenemos datos de largo total.

op <- par(cex.lab = 1.5, font.lab = 2, cex.axis = 1.3, las = 1, bty = "n")
plot(datos$largo, datos$peso, pch = 21, bg = "gray", xlab = "Largo (cm.)", ylab = "Peso (gr.)")

par(op)

Para simular datos “parecidos” a estos podemos hacer que la media de una distribución normal sea una función lineal del largo total del monito del monte. Probemos por ejemplo con \(\mu = -60 + \text{largo} \times 4.2\) y asumiendo que el desvío estándar es de 6:

n <- length(datos$largo)
y <- rnorm(n, mean = -60 + datos$largo * 4.2, sd = 6)

op <- par(cex.lab = 1.5, font.lab = 2, cex.axis = 1.3, las = 1, bty = "n")
plot(datos$largo, y, pch = 21, bg = "gray", xlab = "Largo (cm.)", ylab = "Peso (gr.)")

par(op)

¿Parecido? Cómo saberlo… Durante el curso veremos formas de estimar la relación entre parámetros de distribuciones y variables predictoras (y otras cosas) y cómo ver si los modelos que ajustamos son apropiados para nuestros datos.

Ejercicios:

  1. Genere un set de datos de diferencias en el caudal de un río según las precipitaciones. Asumimos que la relacion entre las precipitaciones (\(pp\)) y la diferencia en el caudal del rio (delta) es la siguiente: \[\delta = -20 + 0.2 \cdot{pp} \] Además asumimos que la variabilidad en el caudal del río es de \(\sigma =5 \text{mm}\)

  2. Simulemos ahora el número de crías de renacuajos que sobrevivieron del total de crías que nacieron según como fue de frío el invierno. Sabemos que la relacion entre la probabilidad de supervivencia de las crias (psup) y la temperatura media del invierno (temp) es la siguiente: \[ \delta = -20 \cdot \left( 1 - e^{0.02\cdot \text{temp}} \right) \]

Juan Manuel Morales . 6 de Septiembre del 2015. Última actualización: 2019-02-19