Distribuciones de Probabilidad e Intervalos de Confianza usando RStudio

Cálculo de Probabilidad, Cuantiles y generación de números aleatorios

Una distribución de probabilidad describe las posibles ocurrencias de un evento y la probabilidad asociada a cada una de esas ocurrencias. En otras palabras, proporciona una función que asigna una probabilidad a cada resultado posible de un experimento aleatorio.

En RStudio, para trabajar con distribuciones de probabilidad, se utilizan funciones específicas para calcular tanto la función de densidad de probabilidad (PDF) como la función de probabilidad acumulada (CDF) de una variable aleatoria. Las funciones que calculan la densidad de probabilidad suelen comenzar con el prefijo d(density) mientras que las funciones que calculan la función de probabilidad acumulada suelen comenzar con el prefijo p(probability). Las funciones que calculan los cuantiles suelen comenzar con el prefijo q(quantile) y para generar números aleatorios el prefijio r(random).

El sufijo añadido se refiere al nombre de la distribución, el cual ayuda a identificar la versión específica de la distribución. Por ejemplo norm para la distribución normal, chisq para una Ji-Cuadrada (Chi-Squared), binom para una Binomial, f para una F, t para una t-Student, pois para una Poisson.

Ejemplo con la Distribución Normal (Caso Continuo)

Para la distribución Normal

• dnorm: Calcula el valor de la PDF en el punto x, f(x).

• pnorm: Calcula el valor de la CDF en el punto x, F(x).

• qnorm: Calcula el percentil \alpha, es decir, retorna el valor de x que acumula dicha probabilidad, X_{\alpha}.

• rnorm: Retorna n valores simulados de la distribución normal.

Los argumentos de la función varían ligeramente dependiendo de la distribución, específicamente, en sus parámetros.

Sea X una variable aleatoria tal que X \sim N(\mu, \sigma^{2}), su PDF esta dada por

f(x) = \frac{1}{\sigma\sqrt{2\pi}}exp\left\{-\frac{1}{2} \left(\frac{x-\mu}{\sigma}\right)^{2} \right\} \quad ;x \in \mathbb{R} del cual su CDF es

F(x) = \mathbb{P}(X < x) = \int_{-\infty}^{x} f(t)dt = \Phi(x) \quad ;x \in \mathbb{R} Para \mu = 2, \sigma^{2} = 4, la densidad evaluada en x = 0 es

\begin{aligned} f(0) &= \frac{1}{2\sqrt{2\pi}}exp\left\{-\frac{1}{2} \left(\frac{0-2}{2}\right)^{2} \right\} \\ & = \frac{1}{2\sqrt{2\pi}}exp\left\{-\frac{1}{2}\right\} \approx 0.1209 \end{aligned}

Empleando la función dnorm, ajustando los argumentos mean = 2 para la media, sd = 2, para la desviación estándar obtenemos:

(dnorm(x = 0, mean = 2, sd = 2))

[1] 0.1209854

Por otra parte, al calcular F(0) = \Phi(0)

\begin{aligned} F(0) = \mathbb{P}(X < 0) &= \int_{-\infty}^{0} f(t)dt \\ &= \int_{-\infty}^{0}\frac{1}{2\sqrt{2\pi}}exp\left\{-\frac{1}{2} \left(\frac{t-2}{2}\right)^{2} \right\}dt \\ &\approx 0.1586 \end{aligned}

Al emplear la función pnorm, la diferencia está en el argumento lower.tail, el cual es igual a TRUE por default y especifica que X \leq x, si es FALSEentonces calcula para el evento X > x.

(pnorm(0, mean = 2, sd = 2, lower.tail = TRUE))

[1] 0.1586553

(1 - pnorm(0, mean = 2, sd = 2, lower.tail = FALSE)) # Por complemento P(X < x) = 1 - P(X > x)

[1] 0.1586553

Con la función qnorm, podemos comprobar el resultado

qnorm(0.1586553, mean = 2, sd = 2, lower.tail = TRUE)

[1] 0

Luego para simular números n = 500 aleatorios de esta distribución es posible emplear la función rnorm

(x <- rnorm(n = 500, mean = 2, sd = 2))

 [1]  3.3132307  4.0900007  4.5601664  1.1921006  2.8920896 -1.1265188
 [7]  2.8049742  5.8920951 -0.6250543  3.9439397

(Se muestran solo las 10 primeras simulaciones)

Si se realiza un histograma de los valores simulados junto a la densidad teórica de X, tenemos

El procedimiento se puede extender a otras distribuciones, con la ligera variante en la definición de sus argumentos. Por ejemplo

• dchisq: Calcula para un x dado su evaluación en la PDF de una distribución Ji-Cuadrada (Chi-squared) con \nu grados de libertad (degrees of freedom) (X \sim\chi_{(\nu)}^{2}). Los grados de libertad se definen en el argumento df.

• pt: Calcula para un x dado su evaluación en la CDF de una distribución t-Student con \nu grados de libertad (X \sim t(\nu)). Los grados de libertad se definen en el argumento df.

• qf: Calcula el cuantil (quantile) para una probabilidad \alpha dada, el valor x que acumula dicha probabilidad para una distribución F de Fisher (X \sim F_{(m_{1}, m_{2})}), es decir X_{\alpha}. Los grados de libertad se definen en el argumento df1 y df2.

• rt: Retorna n valores simulados de la distribución t-Student con \nu grados de libertad(X \sim t(\nu)). Los grados de libertad se definen en el argumento df.

Ejemplo con Distribución Binomial (Caso Discreto)

El caso discreto de variables aleatorias está marcado por la sutileza de que, teóricamente, no hay como tal una función de densidad, sino una Función de Distribución de Probabilidad (PDF).

Si X es una variable aleatoria tal que X \sim Binomial(n, p), su PDF esta dada por

\mathbb{P}(X = x) = \binom{n}{x}p^{x}(1-p)^{n-x} \quad; x \in \{0,1,2,...,n\} mientras que su función de probabilidad acumulada (CDF) está determinada por

\mathbb{P}(X \leq x) = \sum_{k = 0}^{x} \mathbb{P}(X = k) = \sum_{k = 0}^{x} \binom{n}{k}p^{k}(1-p)^{n-k} \quad; x \in \{0,1,2,...,n\}

En este caso las funciones para la distribución Binomial

• dbinom: Calcula el valor de la PF en el punto x, \mathbb{P}(X = x).

• pbinom: Calcula el valor de la CDF en el punto x, \mathbb{P}(X \leq x).

• qbinom: Calcula el cuantil (quantile) para una probabilidad \alpha dada, esto es, retorna el valor de x que acumula dicha probabilidad, X_{\alpha}.

• rbinom: Retorna n valores simulados de la distribución binomial.

Por ejemplo, con una probabilidad de éxito p = 0.5 en n = 30 ensayos,

\begin{aligned} \mathbb{P}(X = 15) & = \binom{30}{15}(0.5)^{15}(0.5)^{30-15} \\ & = \binom{30}{15}(0.25)^{15} \\ & \approx 0.1444 \end{aligned} Usando la función dbinom, ajustando los argumentos size = 30, prob = 0.5

(dbinom(x = 15, size = 30, prob = 0.5))

[1] 0.1444644

Por otra parte si deseamos calcular \mathbb{P}(X\leq15)

\begin{aligned} \mathbb{P}(X\leq15) &= \sum_{k = 0}^{15} \binom{30}{k}0.5^{k}(1-0.5)^{30-k} \\ &= \binom{30}{1}0.5^{1}(1-0.5)^{30-1} + \binom{30}{2}0.5^{2}(1-0.5)^{30-2} + \dots +\binom{30}{15}0.5^{15}(1-0.5)^{30-15} \\ & \approx 0.5722 \end{aligned}

En la función pbinom al igual que en el ejemplo de la distribución Normal, el argumento lower.tail define la probabilidad que está siendo calculada.. Si lower.tail = TRUE, valor por default calcula la probabilidad para el evento X \leq x, si es FALSE entonces se calcula la probabilidad para el evento X > x.

pbinom(q = 15, size = 30, prob = 0.5, lower.tail = TRUE)

[1] 0.5722322

(1 - pbinom(q = 15, size = 30, prob = 0.5, lower.tail = FALSE)) # Por complemento P(X =< x) = 1 - P(X > x)

[1] 0.5722322

Comprobemos el resultado con qbinom

qbinom(0.5722322,size = 30,prob = 0.5,lower.tail = TRUE)

[1] 15

Realizamos una simulación de n = 500 números aleatorios de esta distribución, ahora con el prefijo de la distribución binomial queda rbinom y los parámetros correspondientes, size = 30, prob = 0.5.

(x <- rbinom(n = 500, size = 30, prob = 0.5))

 [1] 17 16 15 12 12 14 14 18 22 15

(Se muestran solo las primeras 10 simulaciones)

Graficamos las proporciones de los valores simulados junto a las proporciones teóricas de X

Las barras representan las probabilidades empíricas para cada valor X simulado, mientras que los puntos rojos representan la probabilidad teórica de una distribución Binomial \sim (30, 0.5).

El procedimiento se puede extender a otras distribuciones discretas, con la ligera variante en la definicion de sus argumentos (parámetros).

Intervalos de Confianza

Un estimador intervalar está definido por dos variables aleatorias \hat{\theta}_{L} = \hat{\theta}_L (X_1, \ldots, X_n) y \hat{\theta}_U = \hat{\theta}_U (X_1, \ldots, X_n), entre las cuales se dice que un parámetro de población \theta se encuentra con alta probabilidad. Es decir

{P} (\hat{\theta}_L \leq \theta \leq \hat {\theta} _U) = 1 - \alpha , en que

• \alpha es un valor pequeño llamado nivel de significancia.
• 1- \alpha es el coeficiente de confianza.
• \hat{\theta}_L se llama límite de confianza inferior.
• \hat{\theta}_U se llama límite de confianza superior.

IC para la Media \mu

\sigma^{2} Conocida

Sea X_1, \ldots, X_n una muestra aleatoria (secuencia de v.a. iid) de una población con media \mu y varianza \sigma^2 conocida. Entonces

Z = \frac{\bar{X} - \mu}{\sigma / \sqrt{n}} \sim N(0, 1)

Un intervalo de 100 (1 - \alpha) \% de confianza para \mu es dado por:

IC(\mu, 1 - \alpha) = \left[ \bar{x} - z_{1 - \alpha / 2}\frac{\sigma}{\sqrt{n}} ; \bar{x} + z_{1 - \alpha / 2} \frac{\sigma}{\sqrt{n}} \right] Adicionalmente, el margen de error, denotado por m, es dado por: m = z_{1-\alpha/2} \dfrac{\sigma}{\sqrt{n}}

Ejemplo:

Se mide el tiempo de respuesta de un sistema informático y se obtiene una muestra de 25 tiempos de respuesta los cuales tienen un tiempo promedio de 10 segundos. Se sabe que la desviación estándar poblacional de estos tiempos de respuesta es de 2.5 segundos. Con base en esta muestra, se desea calcular un intervalo de confianza del 95% para la media del tiempo de respuesta del sistema.

Tenemos la siguiente información:

• n = 25
• \sigma = 2.5
• \bar{x} = 10
• (1 - \alpha)\times 100\% = 95\%; \alpha = 0.05

n <- 25
sigma <- 2.5
xBar <- 10
alpha <- 0.05

Inicialmente calculamos el margen de error, para esto necesitamos z_{1-\alpha/2}

(z_alpha <- qnorm(1 - (alpha/2), mean = 0, sd = 1, lower.tail = TRUE))

[1] 1.959964

Luego calculamos m = z_{1-\alpha/2} \dfrac{\sigma}{\sqrt{n}}

m <- z_alpha*(sigma/sqrt(n))

El IC esta dado por

c(LI = xBar - m, LS = xBar + m)

       LI        LS 
 9.020018 10.979982 

Con un nivel de confianza del 95%, se estima que el verdadero valor medio del tiempo de respuesta del sistema informático está entre 9.02 y 10.98 segundos, es decir, si se repitiera el proceso de muestreo muchas veces y se calculara un intervalo de confianza para la media en cada ocasión, aproximadamente el 95% de estos intervalos contendrían el verdadero valor medio del tiempo de respuesta del sistema.

\sigma^{2} Desconocida

Sea X_1, \ldots, X_n una muestra aleatoria de la población con distribución normal con media \mu y varianza \sigma^2 desconocida. Usaremos la varianza muestral s^2 como una estimación de \sigma^2. Como consecuencia, obtenemos una distribución t-student con n - 1 grados de libertad:

T = \frac{\bar{X} - \mu}{s/\sqrt{n}} \sim t_{(n - 1)}

Un intervalo de 100(1 - \alpha) \% de confianza para \mu es dado por:

IC(\mu, 1 - \alpha) = \left[ \bar{x} - t_{(n-1, 1-\alpha/2)} \frac{s}{\sqrt{n}}; \bar{x} + t_{(n-1,1-\alpha/2)} \frac{s}{\sqrt{n}} \right] Adicionalmente, el margen de error, denotado por m, es dado por: m = t_{(n-1,1-\alpha/2)} \frac{s}{\sqrt{n}}

Ejemplo:

En un laboratorio de salud se está investigando el efecto de un nuevo tratamiento para reducir el peso en adultos con sobrepeso. Para evaluar la efectividad del tratamiento, se realiza un estudio donde se registran los pesos de una muestra de 30 participantes antes de iniciar el tratamiento. Los datos recolectados se almacenan en un conjunto de datos llamado pesos.

pesos <- c(68.5, 70.2, 72.1, 65.8, 63.4, 67.9, 71.3, 69.6, 68.7, 70.4,
           66.2, 64.9, 68.3, 70.1, 73.0, 72.5, 68.6, 69.3, 67.2, 71.9,
           70.8, 72.6, 66.7, 68.8, 70.5, 69.4, 71.2, 68.0, 67.3, 65.5)

El objetivo es calcular un intervalo de confianza del 90% para la media del peso de la población de adultos con sobrepeso.

Se cuenta con la siguiente información

• n = 30
• (1 - \alpha)\times 100\% = 90\%; \alpha = 0.1

n <- 30
alpha <- 0.1

Primero calculamos las estimaciones \bar{x} = \frac{1}{n} \sum_{i = 1}^{n}x_{i} y s = \sqrt{\frac{1}{n-1} \sum_{i=1}^n (x_i - \bar{x})^2}

(xBar <- mean(pesos))

[1] 69.02333

(s <- sd(pesos))

[1] 2.443949

Inicialmente calculamos el margen de error, para esto necesitamos t_{(n-1, 1-\alpha/2)}

(t_alpha <- qt(1 - (alpha/2), df = n - 1, lower.tail = TRUE))

[1] 1.699127

Luego calculamos el margen de error m = t_{(n-1, 1-\alpha/2)} \frac{s}{\sqrt{n}}

(m <- t_alpha*(s/sqrt(n)))

[1] 0.7581538

El IC esta dado por

c(LI = xBar - m, LS = xBar + m)

      LI       LS 
68.26518 69.78149 

Con un nivel de confianza del 90%, se estima que el verdadero valor medio del peso de la población de adultos con sobrepeso está entre 68.26 kg y 69.78 kg.

IC para la Varianza \sigma^{2}

Sabiendo que X_1, \ldots, X_n es una muestra aleatoria de X \sim N(\mu, \sigma^2), con \mu y \sigma^2 desconocidas. Es posible mostrar que:

\frac{(n-1) s^2}{\sigma^2} \sim \chi^2_{(n-1)} Un intervalo de 100 (1 - \alpha)\% de confianza para \sigma^2 es dado por:

IC(\sigma^2, 1 - \alpha) = \left[ \frac{(n-1) s^2}{q_2} ; \frac{(n - 1) s^2}{q_1} \right]

En que q_1 = \chi_{(n-1, \alpha/2)}^2 y q_2 = \chi_{(n-1, 1-\alpha/2)}^2

Ejemplo:

Una empresa de materiales compuestos está investigando la resistencia de un nuevo material para su uso en aplicaciones de ingeniería. Se selecciona una muestra aleatoria de 20 especímenes del material y se mide su resistencia a la tracción. Los datos de resistencia a la tracción de los especímenes son los siguientes (en MPa):

resistencia <- c(73.2, 71.5, 75.1, 69.8, 70.3, 72.6, 76.5, 74.2, 73.9, 70.7,
72.1, 71.3, 75.8, 76.0, 72.9, 70.4, 74.7, 72.8, 74.3, 73.1)

La empresa desea calcular un intervalo de confianza del 95% para la varianza de la resistencia a la tracción de este nuevo material compuesto.

Se tiene la siguiente información

• n = 20

• (1 - \alpha)\times 100\% = 95\%; \alpha = 0.5

n <- 20
alpha <- 0.05

Primero calculamos s^{2}

(s_2 <- var(resistencia))

[1] 4.002526

Luego q_1 y q_2

(q_1 <- qchisq(alpha/2, df = n - 1,lower.tail = TRUE))

[1] 8.906516

(q_2 <- qchisq(1 - (alpha/2), df = n - 1, lower.tail = TRUE))

[1] 32.85233

El IC esta dado por

c(LI = ((n - 1)*s_2)/q_2, LS = ((n - 1)*s_2)/q_1)

      LI       LS 
2.314844 8.538467 

Con un nivel de confianza del 95%, se estima que la verdadera varianza de la resistencia a la tracción del nuevo material compuesto está entre 2.3148 MPa^2 y 8.5384 MPa^2.

IC para la Proporción p

Por el Teorema del Límite Central: la distribución muestral de \widehat{p} se aproxima de la distribución normal cuando n es suficientemente grande:

\widehat{p} \stackrel{a}{\sim} N\left( p, \frac{p (1 - p)}{n} \right) Si usamos la estimación de \widehat{p} para también estimar el error estándar \sqrt{\dfrac{p (1 - p)}{n}}, podemos construir el siguiente IC de 100 ( 1 - \alpha) \%:

IC(p, 1 - \alpha) = \left[ \widehat{p} - z_{1-\alpha/2} \sqrt{\frac{\widehat{p} (1 - \widehat{p})}{n}}; \widehat{p} + z_{1-\alpha/2} \sqrt{\frac{\widehat{p} (1 - \widehat{p})}{n}} \right] Adicionalmente, el margen de error, denotado por m, es dado por: m = z_{1-\alpha/2} \sqrt{\frac{\widehat{p} (1 - \widehat{p})}{n}}

Ejemplo:

En una empresa de investigación de mercado, un analista de datos se encarga de determinar la proporción de clientes satisfechos con un nuevo producto recientemente lanzado. Para ello, se aplicó un formulario a 400 clientes seleccionados de forma aleatoria que han utilizado el producto en los últimos tres meses de los cuales 280 estan satisfechos con el nuevo producto. Esta empresa desea obtener un intervalo de confianza del 97% para la proporción de clientes satisfechos.

Se tiene la siguiente información:

• n = 400, suficientemente grande para ser soportado por TLC
• \hat{p} = \frac{\text{Número de clientes Satisfechos}}{\text{Número total de clientes}} = \frac{280}{400}
• (1 - \alpha)\times 100\% = 97\%; \alpha = 0.03

n <- 400
hat_p <- 180/400 
alpha <- 0.03

Para el calculo del margen de error calculamos inicialmente z_{1-\alpha/2}

z_alpha <- qnorm(1 - (alpha/2), mean = 0, sd = 1, lower.tail = TRUE)

Calculamos m = z_{1-\alpha/2} \sqrt{\frac{\widehat{p} (1 - \widehat{p})}{n}}

m <- z_alpha*sqrt((hat_p * (1 - hat_p))/n)

El IC esta dado por

c(LI = hat_p - m, LS = hat_p + m)

       LI        LS 
0.3960197 0.5039803 

Con un nivel de confianza del 95%, se puede afirmar que aproximadamente entre el 39.60% y el 50.03% de los clientes están satisfechos con el nuevo producto.