INTRODUCCIÓN

El informe presenta el desarrollo de la simulación estocástica para evaluar el comportamiento de variables aleatorias discretas y continuas bajo diferentes tamaños de muestra (10, 50, 100, 1000, 10000). El objetivo es analizar visual y numéricamente la convergencia de los promedios muestrales hacia sus valores esperados teóricos.

DISTRIBUCIÓN DE POISSON

Ejercicio de agronomía, pulgones en tomate.

Un centro de investigación agrícola monitorea la presencia de pulgones en plantas de tomate de tres invernaderos (I1, I2, I3). El número de pulgones por planta se modela con una distribución de Poisson. Estudios previos en zonas vecinas sufieren que el promedio verdadero es de lambda=4 pulgones por planta. El centro planea inspeccionar n=150 plantas por invernadero.

Datos del ejercicio

#Definición de variables

n <- 150 # Número de plantas inspeccionadas 
lambda <- 4 # Promedio real de pulgones por planta
z_critico <- qnorm(1 - 0.05/2)

Simulación con 10 iteraciones

B <- 10 # Escenarios simulados
set.seed(4321)

media_estimada <- replicate(B, mean(rpois(n, lambda = lambda)))

media_simulada <- mean(media_estimada)
ic_inferior <- media_simulada - z_critico * sqrt(media_simulada / n)
ic_superior <- media_simulada + z_critico * sqrt(media_simulada / n)

print(paste("B =", B, "| Media simulada =", round(media_simulada, 4),
            "| IC 95% = (", round(ic_inferior, 4), ",", round(ic_superior, 4), ")"))
## [1] "B = 10 | Media simulada = 4.0933 | IC 95% = ( 3.7696 , 4.4171 )"
par(mfrow = c(1, 2))

hist(media_estimada, breaks = 30, prob = T,
     main = paste("Histograma, B =", B), xlab = "Media muestral")
abline(v = c(ic_inferior, ic_superior), lty = 2, col = "red")
curve(dnorm(x, mean = lambda, sd = sqrt(lambda / n)), add = T, col = "blue")

prom_acumulado <- cumsum(media_estimada) / (1:B)
plot(prom_acumulado, type = "l", log = "x",
     main = paste("Convergencia, B =", B),
     xlab = "Iteraciones", ylab = "Promedio acumulado")
abline(h = lambda, lty = 2, col = "red")

Simulación con 50 iteraciones

#B=50
B <- 50
set.seed(4321)

media_estimada <- replicate(B, mean(rpois(n, lambda = lambda)))

media_simulada <- mean(media_estimada)
ic_inferior <- media_simulada - z_critico * sqrt(media_simulada / n)
ic_superior <- media_simulada + z_critico * sqrt(media_simulada / n)

print(paste("B =", B, "| Media simulada =", round(media_simulada, 4),
            "| IC 95% = (", round(ic_inferior, 4), ",", round(ic_superior, 4), ")"))
## [1] "B = 50 | Media simulada = 4.0329 | IC 95% = ( 3.7116 , 4.3543 )"
par(mfrow = c(1, 2))

hist(media_estimada, breaks = 30, prob = T,
     main = paste("Histograma, B =", B), xlab = "Media muestral")
abline(v = c(ic_inferior, ic_superior), lty = 2, col = "red")
curve(dnorm(x, mean = lambda, sd = sqrt(lambda / n)), add = T, col = "blue")

prom_acumulado <- cumsum(media_estimada) / (1:B)
plot(prom_acumulado, type = "l", log = "x",
     main = paste("Convergencia, B =", B),
     xlab = "Iteraciones", ylab = "Promedio acumulado")
abline(h = lambda, lty = 2, col = "red")

Simulación con 100 iteraciones

#B= 100
B <- 100
set.seed(4321)

media_estimada <- replicate(B, mean(rpois(n, lambda = lambda)))

media_simulada <- mean(media_estimada)
ic_inferior <- media_simulada - z_critico * sqrt(media_simulada / n)
ic_superior <- media_simulada + z_critico * sqrt(media_simulada / n)

print(paste("B =", B, "| Media simulada =", round(media_simulada, 4),
            "| IC 95% = (", round(ic_inferior, 4), ",", round(ic_superior, 4), ")"))
## [1] "B = 100 | Media simulada = 4.0114 | IC 95% = ( 3.6909 , 4.3319 )"
par(mfrow = c(1, 2))

hist(media_estimada, breaks = 30, prob = T,
     main = paste("Histograma, B =", B), xlab = "Media muestral")
abline(v = c(ic_inferior, ic_superior), lty = 2, col = "red")
curve(dnorm(x, mean = lambda, sd = sqrt(lambda / n)), add = T, col = "blue")

prom_acumulado <- cumsum(media_estimada) / (1:B)
plot(prom_acumulado, type = "l", log = "x",
     main = paste("Convergencia, B =", B),
     xlab = "Iteraciones", ylab = "Promedio acumulado")
abline(h = lambda, lty = 2, col = "red")

Simulación con 1000 iteraciones

#B=1000
B <- 1000
set.seed(4321)

media_estimada <- replicate(B, mean(rpois(n, lambda = lambda)))

media_simulada <- mean(media_estimada)
ic_inferior <- media_simulada - z_critico * sqrt(media_simulada / n)
ic_superior <- media_simulada + z_critico * sqrt(media_simulada / n)

print(paste("B =", B, "| Media simulada =", round(media_simulada, 4),
            "| IC 95% = (", round(ic_inferior, 4), ",", round(ic_superior, 4), ")"))
## [1] "B = 1000 | Media simulada = 3.9906 | IC 95% = ( 3.6709 , 4.3103 )"
par(mfrow = c(1, 2))

hist(media_estimada, breaks = 30, prob = T,
     main = paste("Histograma, B =", B), xlab = "Media muestral")
abline(v = c(ic_inferior, ic_superior), lty = 2, col = "red")
curve(dnorm(x, mean = lambda, sd = sqrt(lambda / n)), add = T, col = "blue")

prom_acumulado <- cumsum(media_estimada) / (1:B)
plot(prom_acumulado, type = "l", log = "x",
     main = paste("Convergencia, B =", B),
     xlab = "Iteraciones", ylab = "Promedio acumulado")
abline(h = lambda, lty = 2, col = "red")

Simulación con 10000 iteraciones

#B=10000
B <- 10000
set.seed(4321)

media_estimada <- replicate(B, mean(rpois(n, lambda = lambda)))

media_simulada <- mean(media_estimada)
ic_inferior <- media_simulada - z_critico * sqrt(media_simulada / n)
ic_superior <- media_simulada + z_critico * sqrt(media_simulada / n)

print(paste("B =", B, "| Media simulada =", round(media_simulada, 4),
            "| IC 95% = (", round(ic_inferior, 4), ",", round(ic_superior, 4), ")"))
## [1] "B = 10000 | Media simulada = 4.0017 | IC 95% = ( 3.6816 , 4.3218 )"
par(mfrow = c(1, 2))

hist(media_estimada, breaks = 30, prob = T,
     main = paste("Histograma, B =", B), xlab = "Media muestral")
abline(v = c(ic_inferior, ic_superior), lty = 2, col = "red")
curve(dnorm(x, mean = lambda, sd = sqrt(lambda / n)), add = T, col = "blue")

prom_acumulado <- cumsum(media_estimada) / (1:B)
plot(prom_acumulado, type = "l", log = "x",
     main = paste("Convergencia, B =", B),
     xlab = "Iteraciones", ylab = "Promedio acumulado")
abline(h = lambda, lty = 2, col = "red")

Análisis

En la producción agraria es indispensable la optimización de recursos, la simulación nos permite realizar el muestreo miles de veces para tener una visión de la presencia de plagas sin la necesidad de salir al campo. Al tomar una muestra de n = 150 plantas hay límites: el costo de inspeccionar una por una, la necesidad de que la muestra sea realmente aleatoria para no sesgar el resultado, y el supuesto de que los pulgones se distribuyen sin agruparse (si hay focos de plaga, el modelo Poisson se queda corto).

En los gráficos, con 50 y 100 repeticiones el promedio todavía se mueve bastante; con 1000 ya se pega a lambda = 4 y se queda ahí, y con 10000 solo lo confirma. La estabilización clara está en n = 1000, y en la práctica significa que la media de pulgones por planta es confiable y se puede comparar con el umbral de intervención para decidir si se aplica control químico.

En todo esto actúan dos ideas claves: La Ley de los Grandes Números, es lo que vemos en el gráfico de convergencia: mientras más repeticiones hacemos, más se acerca el promedio al valor real, lambda = 4. El Teorema del Límite Central, se ve en el histograma: aunque los conteos de pulgones por planta no son simétricos, las medias de 150 plantas se distribuyen de forma casi normal alrededor de 4, y por eso la curva azul calza tan bien y podemos usar z = 1.96 para el intervalo de confianza. Dicho simple: la primera nos dice hacia dónde vamos, y el segundo, cómo se reparten los resultados alrededor de ese punto.

DISTRIBUCIÓN EXPONENCIAL

Ejercicio en manufactura, pérdidas por falla de maquinaria.

Una empresa manufacturera estudia las pérdidas económicas X, expresadas en miles de dólares, ocasionadas por fallas inesperadas de maquinaria. Para su análisis, las pérdidas por evento se modelan mediante una distribución exponencial con tasa lambda=0.10. La pérdida media teórica por evento es de 10 mil dólares.

Datos del ejercicio

# 1) Definición de variables 

n<- 100 # fallas inesperadas de maquinaria 
lambda <- 0.10 # Parámetro de riesgo exponencial
media_teorica <- 1/lambda # Perdida media estimada (en miles de dólares)
media_teorica
## [1] 10
z_critico <- qnorm(1-0.05/2)
z_critico
## [1] 1.959964

Simulación con 10 iteraciones

#B=10
B <- 10
set.seed(1223)

perdida_media <- replicate(B, mean(rexp(n, rate = lambda)))
media_simulada <- mean(perdida_media); media_simulada
## [1] 9.971302
se_simulado <- sd(perdida_media)
se_simulado
## [1] 1.050609
ic_inferior <- media_simulada - z_critico*se_simulado 
ic_inferior
## [1] 7.912146
ic_superior <- media_simulada + z_critico*se_simulado 
ic_superior
## [1] 12.03046
print(paste("B =", B, "| Media simulada =", round(media_simulada, 4),
            "| IC 95% = (", round(ic_inferior, 4), ",", round(ic_superior, 4), ")"))
## [1] "B = 10 | Media simulada = 9.9713 | IC 95% = ( 7.9121 , 12.0305 )"
par(mfrow = c(1, 2))

hist(perdida_media, breaks = 5, prob = T, xlim = c(6, 14), ylim = c(0,0.5), col = "gold",
     main = paste("Histograma, B =", B), xlab = "Media muestral", ylab = "Probabilidad")
abline(v = c(ic_inferior, ic_superior), lty = 2, col = "blue")
curve(dnorm(x, mean = media_teorica, sd = sd(perdida_media)), add = T, col = "darkred", lwd = 2)

promedio_acumulado <- cumsum(perdida_media)/ (1:B)

plot(promedio_acumulado, type = "l", col = "darkgreen", lwd = 2,
     main = "Convergencia B=10", xlab = "# Escenarios",
     ylab = "Perdida media acumulada (miles de $)",
     log = "x")
abline(h = media_teorica, col = "red", lty = 2)

Simulación con 100 iteraciones

#B=100
B <- 100
set.seed(1223)

perdida_media <- replicate(B, mean(rexp(n, rate = lambda)))
media_simulada <- mean(perdida_media); media_simulada
## [1] 9.949976
se_simulado <- sd(perdida_media)
se_simulado
## [1] 1.053545
ic_inferior <- media_simulada - z_critico*se_simulado 
ic_inferior
## [1] 7.885066
ic_superior <- media_simulada + z_critico*se_simulado 
ic_superior
## [1] 12.01489
print(paste("B =", B, "| Media simulada =", round(media_simulada, 4),
            "| IC 95% = (", round(ic_inferior, 4), ",", round(ic_superior, 4), ")"))
## [1] "B = 100 | Media simulada = 9.95 | IC 95% = ( 7.8851 , 12.0149 )"
par(mfrow = c(1, 2))

hist(perdida_media, breaks = 15, prob = T, xlim = c(6, 14), ylim = c(0,0.5), col = "gold",
     main = paste("Histograma, B =", B), xlab = "Media muestral", ylab = "Probabilidad")
abline(v = c(ic_inferior, ic_superior), lty = 2, col = "blue")
curve(dnorm(x, mean = media_teorica, sd = sd(perdida_media)), add = T, col = "darkred", lwd = 2)

promedio_acumulado <- cumsum(perdida_media)/ (1:B)

plot(promedio_acumulado, type = "l", col = "darkgreen", lwd = 2,
     main = "Convergencia con B = 100", xlab = "# Escenarios",
     ylab = "Perdida media acumulada (miles de $)",
     log = "x")
abline(h = media_teorica, col = "red", lty = 2)

Simulación con 1000 iteraciones

#B=1000
B <- 1000
set.seed(1223)

perdida_media <- replicate(B, mean(rexp(n, rate = lambda)))
media_simulada <- mean(perdida_media); media_simulada
## [1] 9.981046
se_simulado <- sd(perdida_media)
se_simulado
## [1] 1.010299
ic_inferior <- media_simulada - z_critico*se_simulado 
ic_inferior
## [1] 8.000896
ic_superior <- media_simulada + z_critico*se_simulado 
ic_superior
## [1] 11.9612
print(paste("B =", B, "| Media simulada =", round(media_simulada, 4),
            "| IC 95% = (", round(ic_inferior, 4), ",", round(ic_superior, 4), ")"))
## [1] "B = 1000 | Media simulada = 9.981 | IC 95% = ( 8.0009 , 11.9612 )"
par(mfrow = c(1, 2))

hist(perdida_media, breaks = 15, prob = T, xlim = c(6, 14), ylim = c(0,0.4), col = "gold",
     main = paste("Histograma, B =", B), xlab = "Media muestral", ylab = "Probabilidad")
abline(v = c(ic_inferior, ic_superior), lty = 2, col = "blue")
curve(dnorm(x, mean = media_teorica, sd = sd(perdida_media)), add = T, col = "darkred", lwd = 2)

promedio_acumulado <- cumsum(perdida_media)/ (1:B)

plot(promedio_acumulado, type = "l", col = "darkgreen", lwd = 2,
     main = "Convergencia con B = 1000", xlab = "# Escenarios",
     ylab = "Perdida media acumulada (miles de $)",
     log = "x")
abline(h = media_teorica, col = "red", lty = 2)

Simulación con 10000 iteraciones

#B=10000
B <- 10000
set.seed(1223)

perdida_media <- replicate(B, mean(rexp(n, rate = lambda)))
media_simulada <- mean(perdida_media); media_simulada
## [1] 9.99699
se_simulado <- sd(perdida_media)
se_simulado
## [1] 0.9974146
ic_inferior <- media_simulada - z_critico*se_simulado 
ic_inferior
## [1] 8.042093
ic_superior <- media_simulada + z_critico*se_simulado 
ic_superior
## [1] 11.95189
print(paste("B =", B, "| Media simulada =", round(media_simulada, 4),
            "| IC 95% = (", round(ic_inferior, 4), ",", round(ic_superior, 4), ")"))
## [1] "B = 10000 | Media simulada = 9.997 | IC 95% = ( 8.0421 , 11.9519 )"
par(mfrow = c(1, 2))

hist(perdida_media, breaks = 15, prob = T, xlim = c(6, 14), ylim = c(0,0.4), col = "gold",
     main = paste("Histograma, B =", B), xlab = "Media muestral", ylab = "Probabilidad")
abline(v = c(ic_inferior, ic_superior), lty = 2, col = "blue")
curve(dnorm(x, mean = media_teorica, sd = sd(perdida_media)), add = T, col = "darkred", lwd = 2)

promedio_acumulado <- cumsum(perdida_media)/ (1:B)

plot(promedio_acumulado, type = "l", col = "darkgreen", lwd = 2,
     main = "Convergencia con B = 10000", xlab = "# Escenarios",
     ylab = "Perdida media acumulada (miles de $)",
     log = "x")
abline(h = media_teorica, col = "red", lty = 2)

Análisis

Es de importancia obtener los valores a los que se estabilzan los valores promedio sin exponer a la empresa a grandes pérdidas económicas, además de pérdidas de tiempo en realizar repetidamente este proceso o esperar a que naturalmente sucedan estas fallas para registrarlas, por el contrario, la simulación nos facilita los datos necesarios en minutos.

Igualmente, están presentes: El Teorema del Límite Central, al realizar simulaciones cambiando de número de escenarios o iteraciones es fácil darse cuenta el cumplimiento del Teorema del límite Central, ya que al realizar 10 iteraciones se ven varios espacios en blaco que están dentro de la curva normal, además de algunas barras del histograma que salen de dicha curva normal, lo que sugiere que las medias no dan una indicación clara de distribución normal, sin embargo, al ir incrementando las iteraciones, B=100, B=1000, hasta llegar a B=10000, es fácil darse cuenta que el histograma se apega claramente a la curva normal, dejando pocos espacios en blanco dentro de la curva y pocos espacios en amarillo fuera de la curva normal.

La Ley de los grandes Números, de manera similar al observar el incremento de las iteraciones, partiendo de 10 iteraciones, se puede observar que a menos iteraciones mayor fluctuación tiene las medias simuladas (más alejadas del valor teórico), sin embargo al trabajar con 1000 iteraciones ya se ve que las fluctuaciones se reducen considerablemente, dando a entender que desde este punto se puede considerar que hemos alcanzado o nos hemos acercado lo suficiente al valor teórico de la media (media poblacional) y con mínimas fluctuaciones.