1. Trabajando con distribuciones de probabilidad dadas

En un esfuerzo por describir de manera simplificada el comportanmiento de aquellas variables de la naturaleza que están sujetas a variación no explicada o a incertidumbre, los estadísticos han diseñado una serie de modelos matemáticos que buscan simular dicho comportamiento. En este informe abordaremos algunos cálculos implicados en el uso de esos modelos para resolver problemas de probabilidad. Las distribuciones de probabilidad consideradas son:

  • Distribución normal
  • Distribución binomial
  • Distribución de Poisson

1.1 Distribución normal

Ejemplo:

La vida de un láser de semiconductores con una alimentación de energía constante tiene una distribución normal con una media de 7000 horas y una desviación estándar de 600 horas.

  1. ¿Cuál es la probabilidad de que un láser falle antes de 5000 horas?
  2. ¿Cuál es la vida en horas que exceden el 95% de los láseres?
  3. Si se usan tres láseres en un producto y se supone que fallan de manera independiente, ¿cuál es la probabilidad de que los tres sigan funcionando después de 7000 horas?
1.1.1 Solución item a

Dado que se pregunta por la probabilidad del evento x < 5000 donde x representa las horas de funcionamiento del láser, tenemos que:

\[P(x < 5000) = P(z < z_{0})\] Con \(z_{0} = (5000 - \mu) / \sigma\), y \(z_{0}\) distribuida normal estándar.

Es decir que necesitamos una función en R que nos entregue la probabilidad acumulada para una variable aleatoria normal estándar cuando se suministra el valor \(z_{0}\) hasta el cual deseamos acumular la probabilidad, a esto se le llama el cálculo de una probabilidad acumulada desde \(-\infty\) hasta \(z_0\), donde \(z_0\) es un cuantil(valor) dado o calculado a partir de estandarización. A continuación un código en R que brinda ese cálculo.

mu = 7000
sigma = 600
x = 5000

z_0 = (x-mu) / sigma # Acá calculamos z_0 a partir de la estandarización del valor x dado en el enunciado
p = pnorm(z_0, 0, 1)

# Otra alternativa es no estandarizar y pasar la media y la desviación estándar
# original, al igual que el valor de x tal y como se define en el evento

p2 = pnorm(x, mu, sigma)
print(paste("Esta es la probabilidad estandarizando:", p, "y esta es la probabilidad sin estandarizar", p2))
## [1] "Esta es la probabilidad estandarizando: 0.000429060333196837 y esta es la probabilidad sin estandarizar 0.000429060333196837"

Acá es importante aclarar que la invocación a la función pnorm significa lo siguiente:

\(pnorm(z_0, 0, 1) = P(Z <= z_0)\) con \(Z\) distribuida normal estándar

Y de forma aún más general:

\(pnorm(x_0, \mu, \sigma) = P(X <= x_0)\) con \(X\) distribuida normal con media \(\mu\) y desviación estándar \(\sigma\)

1.1.2 Solución item b

Acá se nos ofrece la pregunta opuesta: dada una probabilidad asociada a un evento, ¿cuál es el cuantil(valor) que la genera?. Simbolicamente se trata de:

\(P(Z > z_0) = 0.95\)

Con lo cual:

\(1 - P(Z <= z_0) = 0.95\), asi que:

\(P(Z <= Z_0) = 0.05\) con Z normal estándar, o lo que es equivalente: \(P(X <= x_0) = 0.05\) con \(X\) normal con media \(\mu\) y desviación estándar \(\sigma\).

Hacemos las anteriores modificaciones con el fin de llegar a un evento de tipo “menor que”, pues tal como se explicó en el apartado anterior, esta es justamente la pregunta que nos responde la función pnorm. Sin embargo invocar pnorm implica conocer \(z_0\) o conocer \(x_0\) pero esto es justamente lo que se nos pregunta, ¿Qué podemos hacer?. Un enfoque sería probar con diferentes valores de \(z_0\) o de \(x_0\) hasta obtener la probabilidad deseada. Veamos:

seq_z = seq(-5, 5, 0.1)
print(seq_z)
##   [1] -5.0 -4.9 -4.8 -4.7 -4.6 -4.5 -4.4 -4.3 -4.2 -4.1 -4.0 -3.9 -3.8 -3.7 -3.6
##  [16] -3.5 -3.4 -3.3 -3.2 -3.1 -3.0 -2.9 -2.8 -2.7 -2.6 -2.5 -2.4 -2.3 -2.2 -2.1
##  [31] -2.0 -1.9 -1.8 -1.7 -1.6 -1.5 -1.4 -1.3 -1.2 -1.1 -1.0 -0.9 -0.8 -0.7 -0.6
##  [46] -0.5 -0.4 -0.3 -0.2 -0.1  0.0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9
##  [61]  1.0  1.1  1.2  1.3  1.4  1.5  1.6  1.7  1.8  1.9  2.0  2.1  2.2  2.3  2.4
##  [76]  2.5  2.6  2.7  2.8  2.9  3.0  3.1  3.2  3.3  3.4  3.5  3.6  3.7  3.8  3.9
##  [91]  4.0  4.1  4.2  4.3  4.4  4.5  4.6  4.7  4.8  4.9  5.0
seq_x = seq(5000, 10000, 100)
print(seq_x)
##  [1]  5000  5100  5200  5300  5400  5500  5600  5700  5800  5900  6000  6100
## [13]  6200  6300  6400  6500  6600  6700  6800  6900  7000  7100  7200  7300
## [25]  7400  7500  7600  7700  7800  7900  8000  8100  8200  8300  8400  8500
## [37]  8600  8700  8800  8900  9000  9100  9200  9300  9400  9500  9600  9700
## [49]  9800  9900 10000
resultados_z <- data.frame(
  "z_0" = seq_z,
  "x_asociado" = seq_z * sigma + mu,
  "p" = pnorm(seq_z, 0, 1),
  "error" = abs(0.05 - pnorm(seq_z, 0, 1))
  )

resultados_x <- data.frame(
  "x_0" = seq_x,
  "p" = pnorm(seq_x, 7000, 600),
  "error" = abs(0.05 - pnorm(seq_x, 7000, 600))
  )

# Calculando el mejor valor para z0
print(head(resultados_z))
##    z_0 x_asociado               p      error
## 1 -5.0       4000 0.0000002866516 0.04999971
## 2 -4.9       4060 0.0000004791833 0.04999952
## 3 -4.8       4120 0.0000007933282 0.04999921
## 4 -4.7       4180 0.0000013008075 0.04999870
## 5 -4.6       4240 0.0000021124547 0.04999789
## 6 -4.5       4300 0.0000033976731 0.04999660
index_min_z0 = which.min(resultados_z$error)
min_error_z0 = resultados_z[index_min_z0, ]
print("Este es el mejor valor usando z0")
## [1] "Este es el mejor valor usando z0"
print(min_error_z0)
##     z_0 x_asociado          p       error
## 35 -1.6       6040 0.05479929 0.004799292
# Calculando el mejor valor para x0
print(head(resultados_x))
##    x_0            p      error
## 1 5000 0.0004290603 0.04957094
## 2 5100 0.0007709848 0.04922902
## 3 5200 0.0013498980 0.04865010
## 4 5300 0.0023032661 0.04769673
## 5 5400 0.0038303806 0.04616962
## 6 5500 0.0062096653 0.04379033
index_min_x0 = which.min(resultados_x$error)
min_error_x0 = resultados_x[index_min_x0, ]
print("Este es el mejor valor usando x0")
## [1] "Este es el mejor valor usando x0"
print(min_error_x0)
##     x_0          p       error
## 11 6000 0.04779035 0.002209648

El proceso anterior arroja una buena estimación del cuantil(valor) deseado, con un error pequeño. Sin embargo, ¿hay una función nativa en R que me haga este tipo de búsquedas?. La respuesta es si:

# Usando estandarización:
z0_optimo <- qnorm(0.05, 0, 1)
x0_estandarizado <- z0_optimo * sigma + mu
print(paste("z0_optimo:", z0_optimo))
## [1] "z0_optimo: -1.64485362695147"
print(paste("x0_optimo basado en estandarizacion:", x0_estandarizado))
## [1] "x0_optimo basado en estandarizacion: 6013.08782382912"
# Sin usar estandarización:
x0_optimo <- qnorm(0.05, 7000, 600)
print(paste("x0_optimo:", x0_optimo))
## [1] "x0_optimo: 6013.08782382912"

Conclusión:

Para calcular un cuantil dada una probabilidad, asegúrese de expresar su evento en la forma “menor que”, y a continuación invoque la función qnorm, bien sea para una normal estándar:

\[z_0 = qnorm(p, 0, 1)\]

En cuyo caso deberá revertir el valor estandarizado al valor de su variable original usando:

\[x_0 = z_0 * \sigma + \mu\]

O bien para una normal con media \(\mu\) y desviación estándadar \(\sigma\):

\[x_0 = qnorm(p, \mu, \sigma)\] En cuyo caso la función le entregará directamente el valor \(x_0\) deseado.

1.1.3 Solución item c

Sean A, B, C los eventos de que el primer, segundo y tercer láser funcionen después de 7000 horas respectivamente. Entonces se solicita

\(P(A \cap B \cap C) = P(A) * P(B) * P(C)\)

La descomposición de la probabilidad del evento intersección como productos de los eventos individuales es posible por que se asume la independencia de los eventos. Como cada evento es en esencia el mismo, ya que cada duración de cada láser se conduce bajo la misma distribución de probabilidad, tenemos:

\(P(A \cap B \cap C) = [P(A)]^3 = [P(X > 7000)]^3 = [1-P(X<=7000)]^3 = [1-pnorm(7000,7000,600)]^3\)

p = (1 -pnorm(7000, 7000, 600))**3
print(paste("La probabilidad es ", p))
## [1] "La probabilidad es  0.125"
1.1.4 Simulación Teorema Central del Límite

El enunciado del teorema establece:

Dado un conjunto de observaciones independientes e idénticamente distribuidas (i.i.d.) \(X_1, X_2, X_3, ..., X_n\) con una media \(\mu\) y una desviación estándar \(\sigma\), la distribución de la media muestral \(\bar{x}\) de muestras de tamaño \(n\) se aproxima a una distribución normal con una media \(\mu\) y una desviación estándar \(\sigma/\sqrt{n}\) a medida que n tiende a infinito.

En otras palabras, a medida que el tamaño de la muestra aumenta, la distribución de las medias muestrales se vuelve más centrada en la media poblacional y menos dispersa, siguiendo una distribución normal independientemente de la distribución subyacente de los datos originales

Este teorema es de gran importancia en inferencia estadística y tiene aplicaciones en diversos campos, ya que permite realizar suposiciones sobre las propiedades de las medias muestrales y calcular intervalos de confianza y pruebas de hipótesis con mayor precisión.

Para el caso particular de una población normal, extraer muestras de tamaño \(n\) garantiza gratuitamente la normalidad de la media muestral para cualquier tamaño de muestra \(n\), es decir, no es requerido que el tamaño de muestra sea especialmente grande. Si por el contrario, la población de la cual se extraen las muestras no es normal, entonces el mínimo tamaño de muestra requerido para que la media muestral se conduzca como una distribución normal con media \(\mu\) y desviación estándar \(\sigma/\sqrt{n}\) dependerá mucho de la forma de la distribución poblacional pero en casi todos los casos n = 30 será suficiente para alcanzar un comportamineto decentemente normal.

1.1.4.1 Preparar diferentes poblaciones de orígen
library(ggplot2)

# Definir posibles poblaciones
poblacion_normal <- rnorm(1000, mean = 50, sd = 10)
poblacion_exponencial <- rexp(1000, rate = 5)
poblacion_binomial <- rbinom(n=1000, p=0.95, size=20)

poblaciones <- list(poblacion_normal, poblacion_exponencial, poblacion_binomial)
names_pob <-c("normal", "exponencial", "binomial")

# Dibujar histogramas para observar las diferencias entre las poblaciones

graph_hist <- function(poblacion, name_pob) {
  data_pob <- data.frame("valor" = poblacion)
  histograma_pob <- ggplot(data = data_pob , aes(x = valor)) +
    geom_histogram(fill = "lightblue", color = "black") +
    labs(title = paste('Histograma para una población {name_pob}'),
         x = "Valores",
         y = "Frecuencia")
  return(histograma_pob)
}

for (i in 1:length(poblaciones)) {
  graph <- graph_hist(poblaciones[[i]], names_pob[i])
  print(graph)
}

# Construyo una función que calcula varias medias muestrales de una población estadística dada
# Cada media es calculada sobre la base de un muestra de tamaño n.
# La cantidad de medias muestrales que calcula dependerá del parámetro num_simulaciones. 

calcular_medias_muestrales <- function(n, poblacion, num_simulaciones) {
  medias_muestrales <- replicate(num_simulaciones, mean(sample(poblacion, n)))
  return(medias_muestrales)
}

# Defino una función para generar todas las medias muestrales de todos los tamaños que deseamos investigar
# Los tamaños a investigar se pasan en el parámetro "tamanos"

get_data <- function(poblacion, tamanos, num_simulaciones){
 
  # Crear un dataframe para almacenar los resultados
  resultados <- data.frame()
  # Simular y almacenar resultados para diferentes tamaños de muestra
  for (n in tamanos) {
    medias_muestrales <- calcular_medias_muestrales(n, poblacion, num_simulaciones)
    resultados <- rbind(resultados, data.frame(TamanoMuestra = n, MediasMuestrales = medias_muestrales))
  }
  return(resultados)
}

# Defino una función que grafica una serie de histogramas para cada tamaño de muestra considerado
plot_histograms <- function(resultados, binwidth, xlim, ylim, name_pob) {
  p <- ggplot(resultados, aes(x = MediasMuestrales)) +
  geom_histogram(binwidth = binwidth, fill = "lightblue", color = "black") +
  facet_wrap(~TamanoMuestra, scales = "free_x", ncol = 5) +
  labs(title = paste("Aproximación del Teorema Central del Límite para poblacion:", name_pob),
       x = "Media Muestral",
       y = "Frecuencia") +
  xlim(xlim) +
  ylim(ylim) +
  theme_minimal()
  print(p)
}

# Defino una función que crea una animación en formato gift con los cambios graduales
# que sufre la distribución de medias muestrales conforme el tamaño de muestra crece
create_animation <- function(resultados, binwidth, xlim, ylim, name_pob, mu, sigma) {
  
  get_media <- function(closest_state) {
  # Filtrar los datos correspondientes a un tamaño de muestra particular
  data_frame_actual <- resultados[resultados$TamanoMuestra == closest_state, ]
  
  # Realizar el cálculo de la media de medias muestrales
  media_actual <- mean(data_frame_actual$MediasMuestrales)
  return(media_actual)
  }

  sd_teorica <- function(closest_state) {
    # Obtener la desviación estándar teórica para un tamaño de muestra pasado en el parámetro
    # "closest_state"
    sd  <- sigma / sqrt(as.double(closest_state))
    return(sd)
  }
  
  get_sd <- function(closest_state) {
    # Filtrar los datos correspondientes al un tamaño de muestra particular
    data_frame_actual <- resultados[resultados$TamanoMuestra == closest_state, ]
    
    # Realizar los cálculos de la desviación estándar de las medias muestrales
    sd_actual <- sd(data_frame_actual$MediasMuestrales) 
    return(sd_actual)
  }
  
  title <- paste('Animación para población', name_pob, '\nTamaño de muestra: {closest_state}')
  x_label <- paste('Media de medias: {format(get_media(closest_state), digits=4)}')
  x_label <- paste(x_label, ',  Media esperada:', mu)
  x_label <- paste(x_label, ', Error: {format(100 * (get_media(closest_state) - mu) / mu, digits=4)}%')
  x_label <- paste0(x_label, '\n')
  x_label <- paste0(x_label, 'sd de medias: {format(get_sd(closest_state), digits=4)}')
  x_label <- paste0(x_label, ', sd de medias teórico: {format(sd_teorica(closest_state), digits=4)}')
  x_label <- paste0(x_label, 
                    ', Error: {format(100 * (get_sd(closest_state) - sd_teorica(closest_state)) 
                    / sd_teorica(closest_state), digits=4)}%')

  myPlot <- ggplot(resultados, aes(x= MediasMuestrales)) +
  geom_histogram(binwidth = binwidth, fill = "lightblue", color = "black") +
  transition_states(TamanoMuestra) +
  labs(title = title,
       x=x_label, y= "Frecuencia") +
  xlim(xlim) +
  ylim(ylim) +
  ease_aes('linear')

  animate(myPlot, duration = 10, fps = 20, renderer = gifski_renderer())
  file = paste0(name_pob, ".gif")
  anim_save(file)
}
1.1.4.2 Realizar comprobación de TCL para población normal
max_pob = max(poblacion_normal)
min_pob = min(poblacion_normal)
binwidth = (max_pob - min_pob) / 20
print(binwidth)
## [1] 2.973772
tamano_muestra_inicial <- 5
tamano_muestra_final <- 100
num_simulaciones <- 1000  # Número de simulaciones para cada tamaño de muestra
tamanos <- seq(tamano_muestra_inicial, tamano_muestra_final, by = 5)

resultados_normal <- get_data(poblacion_normal, tamanos, num_simulaciones)
# Crear el gráfico de distribución de medias muestrales

plot_histograms(resultados_normal, binwidth=0.5, c(40, 60), c(0, 240), "Normal")

Ahora presentemos esto en una animación:

# install.packages("gganimate")
# install.packages("gapminder")
# install.packages("gifski")
library(ggplot2)
library(gganimate)
library(gifski)

create_animation(resultados_normal, binwidth=0.5, c(40, 60), c(0, 240), "Normal", 50, 10)

Primer ejercicio:

Realizar la comprobación de TCL para las otras dos poblaciones. Pista: Tendrá que ajustar de forma conveniente los diferentes parámetros de la función plot_histograms y create_animation, elgiendo el ancho de la barra más apropiado y los límites para los ejes x y y.

max_pob = max(poblacion_exponencial)
min_pob = min(poblacion_exponencial)
binwidth = (max_pob - min_pob) / 100
print(binwidth)
## [1] 0.01638112
tamano_muestra_inicial <- 5
tamano_muestra_final <- 100
num_simulaciones <- 1000  # Número de simulaciones para cada tamaño de muestra
tamanos <- seq(tamano_muestra_inicial, tamano_muestra_final, by = 5)

resultados_exponencial <- get_data(poblacion_exponencial, tamanos, num_simulaciones)
# Crear el gráfico de distribución de medias muestrales

plot_histograms(resultados_exponencial, binwidth=binwidth, c(0, 0.5), c(0, 320), "Exponencial")

Ahora animación para la exponencial:

create_animation(resultados_exponencial, binwidth=binwidth, c(0, 0.5), c(0, 320), "Exponencial", 1/5, 1/5)