Introducción

“El proceso de simulación constituye una herramienta poderosa para la estadística que se pueden emplear para entender relaciones complejas y estimar valores difíciles de calcular directamente. Para entenderlo utilizaremos se plantean los siguientes problemas:”

En el presente documento se encuentra la solución propuesta para los problemas de simulación estadística correspondientes a la actividad número 2 del curso de simulación y métodos estadística de la maestría en Ciencia de Datos de la Pontificia Universidad Javeriana.

Problema 1

Estimación del valor de \(\pi\) \ La siguiente figura sugiere como estimar el valor de π con una simulación. En la figura, un circuito con un área igual a \(\pi/4\), está inscrito en un cuadrado cuya área es igual a 1. Se elige de forma aleatoria n puntos dentro del cuadrado . La probabilidad de que un punto esté dentro del círculo es igual a la fracción del área del cuadrado que abarca a este, la cual es \(\pi/4\) Por tanto, se puede estimar el valor de \(\pi/4\) al contar el número de puntos dentro del círculo, para obtener la estimación de \(\pi/4\) De este último resultado se encontrar una aproximación para el valor de \(\pi\)

Pasos sugeridos:\

  1. Genere \(n\) coordenadas \(x: X_{1}, . . . , X_{n}\). Utilice la distribución uniforme con valor mínimo de \(0\) y valor máximo de \(1\) La distribución uniforme genera variables aleatorias que tienen la misma probabilidad de venir de cualquier parte del intervalo \((0,1)\).
# Definir la función para determinar si un punto está dentro del círculo
esta_dentro_del_circulo = function(coordenada_x, coordenada_y) {
  distancia_desde_centro = (coordenada_x - 0.5)^2 + (coordenada_y - 0.5)^2
  return(distancia_desde_centro < 0.25)
}

Ejercicio para una simulación de n = 1.000

  1. Genere \(1000\) coordenadas \(y: Y_{1},...,Y_{n}\), utilizando nuevamente la distribución uniforme con valor mínimo de \(0\) y valor máximo de \(1\).

  2. Cada punto \((X_{i},Y_{i})\) se encuentra dentro del círculo si su distancia desde el centro \((0.5,0.5)\) es menor a \(0.5\). Para cada par \((X_{i},Y_{i})\) determine si la distancia desde el centro es menor a \(0.5\). Esto último se puede realizar al calcular el valor \((X_{i}−0.5)^{2}+(Y_{i}−0.5)^{2}\), que es el cuadrado de la distancia, y al determinar si es menor que \(0.25\).

  3. ¿Cuántos de los puntos están dentro del círculo? ¿Cuál es su estimación de \(\pi\)?

n = 1000
# Generar n coordenadas de X y Y
set.seed(2023)
coordenada_x = runif(n, min = 0, max = 1)
coordenada_y = runif(n, min = 0, max = 1)

# Determinar si los puntos están dentro del círculo
dentro_del_circulo = esta_dentro_del_circulo(coordenada_x, coordenada_y)

# Contar cuántos puntos están dentro del círculo
puntos_dentro = sum(dentro_del_circulo)

# Calcular la estimación de π
estimacion_pi = 4 * (puntos_dentro / n)

cat("Número de puntos dentro del círculo para n =", n, ":", puntos_dentro, "\n")
## Número de puntos dentro del círculo para n = 1000 : 801
cat("Estimación de π para n =", n, ":", estimacion_pi, "\n")
## Estimación de π para n = 1000 : 3.204
# Crear un gráfico n = 1000
plot(coordenada_x, coordenada_y, col = ifelse(dentro_del_circulo, "#63B8FF", "#EEE9E9"), pch = 20,
     main = "Simulación de puntos aleatorios n = 1.000",
     xlab = "Coordenada X", ylab = "Coordenada Y")

# Dibujar el círculo sin borrar los puntos existentes
theta = seq(0, 2 * pi, length.out = 100)
circ_x = 0.5 + 0.5 * cos(theta)
circ_y = 0.5 + 0.5 * sin(theta)
lines(circ_x, circ_y, col = "#0000CD", lwd = 2)

Ejercicio para una simulación de n = 10.000

n = 10000
set.seed(2030)
coordenada_x = runif(n, min = 0, max = 1)
coordenada_y = runif(n, min = 0, max = 1)

dentro_del_circulo = esta_dentro_del_circulo(coordenada_x, coordenada_y)
puntos_dentro = sum(dentro_del_circulo)
estimacion_pi = 4 * (puntos_dentro / n)

print(paste("Número de puntos dentro del círculo para n =", n, ":", puntos_dentro))
## [1] "Número de puntos dentro del círculo para n = 10000 : 7792"
print(paste("Estimación de π para n =", n, ":", estimacion_pi))
## [1] "Estimación de π para n = 10000 : 3.1168"
# Crear un gráfico n = 10.000
plot(coordenada_x, coordenada_y, col = ifelse(dentro_del_circulo, "#63B8FF", "#EEE9E9"), pch = 20,
     main = "Simulación de puntos aleatorios n = 10.000",
     xlab = "Coordenada X", ylab = "Coordenada Y")
# Dibujar el círculo
theta = seq(0, 2 * pi, length.out = 100)
circ_x = 0.5 + 0.5 * cos(theta)
circ_y = 0.5 + 0.5 * sin(theta)
lines(circ_x, circ_y, col = "#0000CD", lwd = 2)

Ejercicio para una simulación de n = 100.000

n = 100000
set.seed(2030)
coordenada_x = runif(n, min = 0, max = 1)
coordenada_y = runif(n, min = 0, max = 1)

dentro_del_circulo = esta_dentro_del_circulo(coordenada_x, coordenada_y)
puntos_dentro = sum(dentro_del_circulo)
estimacion_pi = 4 * (puntos_dentro / n)

print(paste("Número de puntos dentro del círculo para n =", n, ":", puntos_dentro))
## [1] "Número de puntos dentro del círculo para n = 1e+05 : 78624"
print(paste("Estimación de π para n =", n, ":", estimacion_pi))
## [1] "Estimación de π para n = 1e+05 : 3.14496"
# Crear un gráfico n = 100.000
plot(coordenada_x, coordenada_y, col = ifelse(dentro_del_circulo, "#63B8FF", "#EEE9E9"), pch = 20,
     main = "Simulación de puntos aleatorios n = 100.000",
     xlab = "Coordenada X", ylab = "Coordenada Y")
# Dibujar el círculo
theta = seq(0, 2 * pi, length.out = 100)
circ_x = 0.5 + 0.5 * cos(theta)
circ_y = 0.5 + 0.5 * sin(theta)
lines(circ_x, circ_y, col = "#0000CD", lwd = 2)

Conclusión: En los ejercicios de simulación al aumentar \(n\) bajo la elaboración de este ejercicio el número esperado es igual al número \(\pi\)

Ejercicio 2

Propiedades de los estimadores La simulación ayuda a entender y validad las propiedades de los estimadores estadísticos como son. insesgadez, eficiencia y la consistencia principalmente. El siguiente problema permite evidenciar las principales características de un grupo de estimadores propuestos para la estimación de un parámetro asociado a un modelo de probabilidad.

Sean \(X_{1}, X_{2}, X_{3}\) y \(X_{4}\), una muestra aleatoria de tamaño \(n=4\) cuya población la conforma una distribución exponencial con parámetro desconocido. Determine las características de cada uno de los siguientes estimadores propuestos:

  1. \(\hat{\theta_{1}} = \frac{X_{1}+X_{2}}{6} + \frac{X_{3}+X_{4}}{3}\)

  2. \(\hat{\theta_{2}} = \frac{(X_{1}+2X_{2}+3X_{3}+4X_{4})}{5}\)

  3. \(\hat{\theta_{3}} = \frac{X_{1}+X_{2}+X_{3}+X_{4}}{4}\)

  4. \(\hat{\theta_{4}} = \frac{min \{X_{1},X_{2},X_{3},X_{4} \} + max \{X_{1},X_{2},X_{3},X_{4} \}}{2}\)

  1. Genere una muestras de \(n=20\), \(50\), \(100\) y \(1000\) para cada uno de los estimadores planteados.
Funcion_estimadores = function(estimador, muestra, lambda) {
  # Crear matriz de estimadores
  mx = matrix(rexp(estimador * muestra, rate = lambda), nrow = estimador, ncol = muestra, byrow = TRUE) 
  
  # Seleccionar subconjunto de estimadores
  m_estimadores = mx[1:estimador, ]
  
  # Calcular estimadores para cada muestra
  t1 = (m_estimadores[, 1] + m_estimadores[, 2] / 6) + (m_estimadores[, 3] + m_estimadores[, 4] / 3)
  t2 = (m_estimadores[, 1] + (2 * m_estimadores[, 2]) + (3 * m_estimadores[, 3]) + (4 * m_estimadores[, 4])) / 5
  t3 = ((m_estimadores[, 1] + m_estimadores[, 2] + m_estimadores[, 3] + m_estimadores[, 4]) / 4)
  t4 = (min(m_estimadores[, 1], m_estimadores[, 2], m_estimadores[, 3], m_estimadores[, 4]) + max(m_estimadores[, 1], m_estimadores[, 2], m_estimadores[, 3], m_estimadores[, 4])) / 2
  
  # Graficar boxplot
  boxplot(m_estimadores, xlab = "Estimadores", ylab = "Valores")
  abline(h = lambda,  col = "red") 
  
  # Calcular propiedades de los estimadores
  m_medias = apply(m_estimadores, 2, mean)
  m_varianza = apply(m_estimadores, 2, var)
  sesgo = lambda - m_medias 
  
  # Imprimir propiedades de los estimadores
  valores = data.frame(Estimador = paste("Estimador", 1:4), Media = m_medias, Varianza = m_varianza, Sesgo = sesgo)
  print(valores)
}

  1. En cada caso evalue las propiedades de insesgadez, eficiencia y consistencia

  2. Suponga un valor para el parámetro \(\theta\)

# Llamar a la función con diferentes parámetros
lambda = 1
mue = 4                  
sim = 20     
set.seed(1015)
Funcion_estimadores(sim, mue, lambda)

##     Estimador    Media  Varianza       Sesgo
## 1 Estimador 1 1.099279 0.8542797 -0.09927873
## 2 Estimador 2 1.349316 2.8253037 -0.34931568
## 3 Estimador 3 1.322809 3.1608728 -0.32280896
## 4 Estimador 4 1.066961 0.6389195 -0.06696082
lambda = 1
mue = 4                  
sim = 50     
set.seed(1015)
Funcion_estimadores(sim, mue, lambda)

##     Estimador     Media  Varianza       Sesgo
## 1 Estimador 1 1.0965000 0.9183035 -0.09650002
## 2 Estimador 2 1.1346336 1.6719717 -0.13463364
## 3 Estimador 3 1.0892831 1.6662160 -0.08928305
## 4 Estimador 4 0.9491779 0.4712841  0.05082210
lambda = 1
mue = 4                  
sim = 100     
set.seed(1015)
Funcion_estimadores(sim, mue, lambda)

##     Estimador     Media  Varianza        Sesgo
## 1 Estimador 1 1.1214879 1.4072443 -0.121487883
## 2 Estimador 2 0.9679395 1.0553541  0.032060531
## 3 Estimador 3 1.0250448 1.2348286 -0.025044757
## 4 Estimador 4 0.9918490 0.7367481  0.008150997
lambda = 1
mue = 4                  
sim = 1000    
set.seed(1015)
Funcion_estimadores(sim, mue, lambda)

##     Estimador     Media  Varianza        Sesgo
## 1 Estimador 1 1.0019591 1.0430985 -0.001959072
## 2 Estimador 2 0.9653345 0.8652995  0.034665511
## 3 Estimador 3 0.9575673 0.9374961  0.042432724
## 4 Estimador 4 0.9875813 0.9374443  0.012418745

Problema 3

Teorema del Límite Central

El Teorema del Límite Central es uno de los más importantes en la inferencia estadística y habla sobre la convergencia de los estimadores como la proporción muestral a la distribución normal. Algunos autores afirman que esta aproximación es bastante buena a partir del umbral \(n>30\). A continuación se describen los siguientes pasos para su verificación:

  1. Realice una simulación en la cual genere una población de n=1000 (Lote), donde el porcentaje de individuos (supongamos plantas) enfermas sea del \(50\%\).
set.seed(1997)
LOTE = rbinom(1000 , 1 , .5)

  1. Genere una función que permita: Obtener una muestra aleatoria de la población y Calcule el estimador de la proporción muestral \(\hat{p}\) para un tamaño de muestra dado \(n\).

  2. Repita el escenario anterior (b) \(n=500\) veces y analice los resultados en cuanto al comportamiento de los \(500\) resultados del estimador \(\hat{p}\). ¿Qué tan simétricos o sesgados son los resultados obtenidos? y ¿qué se puede observar en cuanto a la variabilidad?. Realice en su informe un comentario sobre los resultados obtenidos.

muestra_aleatoria_aleatoria = function(vector) {
  tamano_muestra = sample(1:length(vector), 1)
  muestra = sample(vector, size = tamano_muestra, replace = FALSE)
  return(muestra)
}

set.seed(1997)
muestra_1 = muestra_aleatoria_aleatoria(LOTE) ## Se tomo el 17% de la muestra

set.seed(1928)
LOTE_500 = rbinom(500 , 1 , .5)

set.seed(1928)
muestra_2 = muestra_aleatoria_aleatoria(LOTE_500) ## Se tomo el 26% de la muestra
prop.table(table(muestra_2))
## muestra_2
##         0         1 
## 0.5419847 0.4580153
0.5402-0.458 ;rm(LOTE ,LOTE_500 , muestra_1 , muestra_2 )
## [1] 0.0822

Los resultados varian en la unidad de elección 1 Mientras en la primera simulación se obtiene un 54% de los datos hacia plantas enfermas dato la seleción del \(70\%\) de la muestra. en el segundo ejercicio este resultado es el \(45\%\) de los datos una diferencia de 0.08 puntos porcentuales.

  1. Repita los puntos b y c para tamaños de muestra \(n=5, 10, 15, 20, 30, 50, 60, 100, 200, 500\). Compare los resultados obtenidos para los diferentes tamaños de muestra en cuanto a la normalidad. Utilice pruebas de bondad y ajuste (shapiro wilks :shapiro.test()) y métodos gráficos (gráfico de normalidad: qqnorm()). Comente en su informe los resultados obtenidos
set.seed(1943)
LOTE_5 = rbinom(5 , 1 , .5)
LOTE_10 = rbinom(10 , 1 , .5)
LOTE_15 = rbinom(15 , 1 , .5)
LOTE_20 = rbinom(20 , 1 , .5)
LOTE_30 = rbinom(30 , 1 , .5)
LOTE_50 = rbinom(50 , 1 , .5)
LOTE_60 = rbinom(60 , 1 , .5)
LOTE_100 = rbinom(100 , 1 , .5)
LOTE_200 = rbinom(200 , 1 , .5)
LOTE_500 = rbinom(500 , 1 , .5)

# 
set.seed(1943)
muestra_1 = muestra_aleatoria_aleatoria(LOTE_5)
muestra_2 = muestra_aleatoria_aleatoria(LOTE_10)
muestra_3 = muestra_aleatoria_aleatoria(LOTE_15)
muestra_4 = muestra_aleatoria_aleatoria(LOTE_20)
muestra_5 = muestra_aleatoria_aleatoria(LOTE_30)
muestra_6 = muestra_aleatoria_aleatoria(LOTE_50)
muestra_7 = muestra_aleatoria_aleatoria(LOTE_60)
muestra_8 = muestra_aleatoria_aleatoria(LOTE_100)
muestra_9 = muestra_aleatoria_aleatoria(LOTE_200)
muestra_10 = muestra_aleatoria_aleatoria(LOTE_500)

# Realizar la prueba de Shapiro para cada muestra en la lista
muestras = list(muestra_2, muestra_3, muestra_4, muestra_5, muestra_6, muestra_7, muestra_8, muestra_9, muestra_10)

for (muestra in muestras) {
  print(shapiro.test(muestra))
  print(qqnorm(muestra))
}
## 
##  Shapiro-Wilk normality test
## 
## data:  muestra
## W = 0.75, p-value < 2.2e-16

## $x
## [1] -0.8694238  0.0000000  0.8694238
## 
## $y
## [1] 0 0 1
## 
## 
##  Shapiro-Wilk normality test
## 
## data:  muestra
## W = 0.6457, p-value = 0.0001673

## $x
##  [1]  0.0000000 -1.7688250  0.1940281 -1.1983797  0.3957253 -0.8694238
##  [7]  0.6151411 -0.6151411  0.8694238 -0.3957253  1.1983797 -0.1940281
## [13]  1.7688250
## 
## $y
##  [1] 1 0 1 0 1 0 1 0 1 0 1 0 1
## 
## 
##  Shapiro-Wilk normality test
## 
## data:  muestra
## W = 0.63748, p-value = 7.387e-06

## $x
##  [1]  0.18911843  0.31863936  0.45376219 -1.95996398  0.59776013 -1.43953147
##  [7]  0.75541503 -1.15034938 -0.93458929  0.93458929 -0.75541503  1.15034938
## [13] -0.59776013 -0.45376219 -0.31863936 -0.18911843  1.43953147 -0.06270678
## [19]  0.06270678  1.95996398
## 
## $y
##  [1] 1 1 1 0 1 0 1 0 0 1 0 1 0 0 0 0 1 0 0 1
## 
## 
##  Shapiro-Wilk normality test
## 
## data:  muestra
## W = 0.62622, p-value = 4.249e-07

## $x
##  [1] -2.08535557  0.28221615 -1.59321882  0.38032564  0.48224821 -1.32495769
##  [7] -1.12814365 -0.96742157  0.58945580  0.70392179 -0.82846465 -0.70392179
## [13] -0.58945580 -0.48224821  0.82846465 -0.38032564 -0.28221615 -0.18675612
## [19]  0.96742157  1.12814365  1.32495769 -0.09297185  0.00000000  0.09297185
## [25]  1.59321882  2.08535557  0.18675612
## 
## $y
##  [1] 0 1 0 1 1 0 0 0 1 1 0 0 0 0 1 0 0 0 1 1 1 0 0 0 1 1 0
## 
## 
##  Shapiro-Wilk normality test
## 
## data:  muestra
## W = 0.61691, p-value = 1.651e-07

## $x
##  [1]  0.35293399  0.44658820 -2.11438077 -1.62836141 -1.36448875  0.54434091
##  [7] -1.17154617  0.64760358 -1.01449875 -0.87916772  0.75829256 -0.75829256
## [13]  0.87916772 -0.64760358  1.01449875 -0.54434091 -0.44658820 -0.35293399
## [19] -0.26228277  1.17154617 -0.17374106  1.36448875 -0.08654337  0.00000000
## [25]  1.62836141  0.08654337  0.17374106  0.26228277  2.11438077
## 
## $y
##  [1] 1 1 0 0 0 1 0 1 0 0 1 0 1 0 1 0 0 0 0 1 0 1 0 0 1 0 0 0 1
## 
## 
##  Shapiro-Wilk normality test
## 
## data:  muestra
## W = 0.63682, p-value = 2.616e-09

## $x
##  [1] -2.28654795  0.05573169  0.11163715 -1.83391464 -1.59321882 -1.42017907
##  [7]  0.16789400 -1.28155157  0.22468772 -1.16394958 -1.06056224  0.28221615
## [13]  0.34069483 -0.96742157 -0.88199821  0.40036338  0.46149369 -0.80257189
## [19]  0.52440051 -0.72791329 -0.65710857  0.58945580  0.65710857 -0.58945580
## [25]  0.72791329  0.80257189  0.88199821  0.96742157 -0.52440051 -0.46149369
## [31]  1.06056224  1.16394958  1.28155157  1.42017907 -0.40036338 -0.34069483
## [37] -0.28221615  1.59321882 -0.22468772  1.83391464 -0.16789400 -0.11163715
## [43]  2.28654795 -0.05573169  0.00000000
## 
## $y
##  [1] 0 1 1 0 0 0 1 0 1 0 0 1 1 0 0 1 1 0 1 0 0 1 1 0 1 1 1 1 0 0 1 1 1 1 0 0 0 1
## [39] 0 1 0 0 1 0 0
## 
## 
##  Shapiro-Wilk normality test
## 
## data:  muestra
## W = 0.63647, p-value = 5.953e-08

## $x
##  [1] -2.17792307 -0.03687053  0.03687053 -1.70478090 -1.44999891  0.11081291
##  [7]  0.18536702 -1.26496924  0.26096740  0.33809166  0.41728414  0.49918784
## [13]  0.58458986 -1.11533736  0.67448975  0.77020808  0.87356913 -0.98723099
## [19]  0.98723099 -0.87356913 -0.77020808  1.11533736 -0.67448975  1.26496924
## [25]  1.44999891 -0.58458986  1.70478090 -0.49918784 -0.41728414  2.17792307
## [31] -0.33809166 -0.26096740 -0.18536702 -0.11081291
## 
## $y
##  [1] 0 1 1 0 0 1 1 0 1 1 1 1 1 0 1 1 1 0 1 0 0 1 0 1 1 0 1 0 0 1 0 0 0 0
## 
## 
##  Shapiro-Wilk normality test
## 
## data:  muestra
## W = 0.63638, p-value = 3.414e-13

## $x
##  [1] -2.51912447  0.00000000  0.02949402  0.05901372 -2.10496772 -1.88950996
##  [7]  0.08858496 -1.73719396  0.11823387  0.14798711  0.17787196  0.20791656
## [13]  0.23815006  0.26860288 -1.61698114 -1.51642562  0.29930691 -1.42921942
## [19] -1.35170224 -1.28155157  0.33029579  0.36160523 -1.21719663 -1.15752281
## [25]  0.39327329 -1.10170893  0.42534083  0.45785193  0.49085443 -1.04913140
## [31] -0.99930480 -0.95184328  0.52440051  0.55854749  0.59335862  0.62890422
## [37]  0.66526290 -0.90643457 -0.86282209 -0.82079209 -0.78016439 -0.74078547
## [43]  0.70252320  0.74078547 -0.70252320  0.78016439 -0.66526290 -0.62890422
## [49] -0.59335862 -0.55854749  0.82079209  0.86282209  0.90643457 -0.52440051
## [55] -0.49085443 -0.45785193 -0.42534083  0.95184328  0.99930480  1.04913140
## [61]  1.10170893 -0.39327329 -0.36160523 -0.33029579 -0.29930691  1.15752281
## [67] -0.26860288  1.21719663  1.28155157  1.35170224 -0.23815006  1.42921942
## [73] -0.20791656  1.51642562 -0.17787196 -0.14798711 -0.11823387  1.61698114
## [79]  1.73719396  1.88950996  2.10496772 -0.08858496 -0.05901372  2.51912447
## [85] -0.02949402
## 
## $y
##  [1] 0 1 1 1 0 0 1 0 1 1 1 1 1 1 0 0 1 0 0 0 1 1 0 0 1 0 1 1 1 0 0 0 1 1 1 1 1 0
## [39] 0 0 0 0 1 1 0 1 0 0 0 0 1 1 1 0 0 0 0 1 1 1 1 0 0 0 0 1 0 1 1 1 0 1 0 1 0 0
## [77] 0 1 1 1 1 0 0 1 0
## 
## 
##  Shapiro-Wilk normality test
## 
## data:  muestra
## W = 0.62378, p-value < 2.2e-16

## $x
##   [1] -0.228075394 -0.214833790 -2.797207657 -2.421327991 -0.201629750
##   [6] -0.188460772 -0.175324397 -2.229612250 -2.095996807 -1.991811422
##  [11] -1.905608404 -0.162218209 -1.831604880 -0.149139828 -1.766454744
##  [16] -1.708040295 -0.136086908 -0.123057134 -0.110048219 -0.097057902
##  [21] -1.654932729 -0.084083943 -0.071124123 -1.606120563 -0.058176239
##  [26] -0.045238100 -1.560860364 -0.032307530 -0.019382360 -1.518589009
##  [31] -0.006460427  0.006460427 -1.478869270  0.019382360 -1.441354576
##  [36]  0.032307530  0.045238100  0.058176239 -1.405765350 -1.371872628
##  [41] -1.339486424  0.071124123  0.084083943  0.097057902  0.110048219
##  [46] -1.308447267  0.123057134  0.136086908  0.149139828  0.162218209
##  [51] -1.278619931  0.175324397 -1.249888705 -1.222153765  0.188460772
##  [56]  0.201629750  0.214833790 -1.195328362 -1.169336603  0.228075394
##  [61]  0.241357111  0.254681544  0.268051349 -1.144111692  0.281469243
##  [66] -1.119594506  0.294938007  0.308460491 -1.095732453  0.322039619
##  [71]  0.335678392  0.349379897  0.363147310  0.376983904  0.390893054
##  [76] -1.072478525  0.404878245  0.418943080 -1.049790524 -1.027630415
##  [81] -1.005963788 -0.984759399  0.433091285  0.447326723  0.461653396
##  [86] -0.963988789 -0.943625954  0.476075461  0.490597238  0.505223222
##  [91]  0.519958094 -0.923647061 -0.904030212 -0.884755232  0.534806735
##  [96]  0.549774243  0.564865945  0.580087415 -0.865803484  0.595444493
## [101] -0.847157720  0.610943307 -0.828801932  0.626590290 -0.810721240
## [106]  0.642392210  0.658356192  0.674489750  0.690800818 -0.792901778
## [111]  0.707297783 -0.775330604  0.723989527  0.740885469 -0.757995614
## [116]  0.757995614  0.775330604  0.792901778  0.810721240  0.828801932
## [121] -0.740885469 -0.723989527 -0.707297783 -0.690800818 -0.674489750
## [126]  0.847157720 -0.658356192  0.865803484 -0.642392210  0.884755232
## [131] -0.626590290  0.904030212  0.923647061  0.943625954  0.963988789
## [136] -0.610943307 -0.595444493 -0.580087415  0.984759399 -0.564865945
## [141]  1.005963788  1.027630415 -0.549774243  1.049790524 -0.534806735
## [146] -0.519958094 -0.505223222 -0.490597238 -0.476075461  1.072478525
## [151] -0.461653396  1.095732453  1.119594506  1.144111692  1.169336603
## [156]  1.195328362 -0.447326723  1.222153765  1.249888705  1.278619931
## [161]  1.308447267  1.339486424 -0.433091285 -0.418943080 -0.404878245
## [166] -0.390893054 -0.376983904  1.371872628  1.405765350 -0.363147310
## [171]  1.441354576  1.478869270 -0.349379897  1.518589009  1.560860364
## [176] -0.335678392 -0.322039619  1.606120563  1.654932729  1.708040295
## [181]  1.766454744 -0.308460491 -0.294938007 -0.281469243  1.831604880
## [186]  1.905608404 -0.268051349  1.991811422  2.095996807  2.229612250
## [191] -0.254681544  2.421327991  2.797207657 -0.241357111
## 
## $y
##   [1] 1 1 0 0 1 1 1 0 0 0 0 1 0 1 0 0 1 1 1 1 0 1 1 0 1 1 0 1 1 0 1 1 0 1 0 1 1
##  [38] 1 0 0 0 1 1 1 1 0 1 1 1 1 0 1 0 0 1 1 1 0 0 1 1 1 1 0 1 0 1 1 0 1 1 1 1 1
##  [75] 1 0 1 1 0 0 0 0 1 1 1 0 0 1 1 1 1 0 0 0 1 1 1 1 0 1 0 1 0 1 0 1 1 1 1 0 1
## [112] 0 1 1 0 1 1 1 1 1 0 0 0 0 0 1 0 1 0 1 0 1 1 1 1 0 0 0 1 0 1 1 0 1 0 0 0 0
## [149] 0 1 0 1 1 1 1 1 0 1 1 1 1 1 0 0 0 0 0 1 1 0 1 1 0 1 1 0 0 1 1 1 1 0 0 0 1
## [186] 1 0 1 1 1 0 1 1 0
qqnorm(muestra_1)

  1. Repita toda la simulación (puntos a – d), pero ahora para lotes con \(10\%\) de plantas enfermas y de nuevo para lotes con un \(90\%\) de plantas enfermas. Concluya sobre los resultados del ejercicio.
# 10%
set.seed(1933)
LOTE_5 = rbinom(5 , 1 , .1)
muestra_1 = muestra_aleatoria_aleatoria(LOTE_5) ## Se tomo el 17% de la muestra
prop.table(table(muestra_1))
## muestra_1
## 0 
## 1
#shapiro.test(muestra_1)
qqnorm(muestra_1)

LOTE_10 = rbinom(10 , 1 , .1)
muestra_2 = muestra_aleatoria_aleatoria(LOTE_10) ## Se tomo el 17% de la muestra
prop.table(table(muestra_2))
## muestra_2
##   0   1 
## 0.9 0.1
#shapiro.test(muestra_2)
qqnorm(muestra_2)

LOTE_15 = rbinom(15 , 1 , .1)
muestra_3 = muestra_aleatoria_aleatoria(LOTE_15) ## Se tomo el 17% de la muestra
prop.table(table(muestra_3))
## muestra_3
##     0     1 
## 0.875 0.125
shapiro.test(muestra_3)
## 
##  Shapiro-Wilk normality test
## 
## data:  muestra_3
## W = 0.4184, p-value = 1.047e-06
qqnorm(muestra_3)

LOTE_20 = rbinom(20 , 1 , .1)
muestra_4 = muestra_aleatoria_aleatoria(LOTE_20) ## Se tomo el 17% de la muestra
prop.table(table(muestra_4))
## muestra_4
## 0 
## 1
#shapiro.test(muestra_4)
qqnorm(muestra_4)

LOTE_30 = rbinom(30 , 1 , .1)
muestra_5 = muestra_aleatoria_aleatoria(LOTE_30) ## Se tomo el 17% de la muestra
prop.table(table(muestra_5))
## muestra_5
##   0   1 
## 0.5 0.5
shapiro.test(muestra_5)
## 
##  Shapiro-Wilk normality test
## 
## data:  muestra_5
## W = 0.72863, p-value = 0.02386
qqnorm(muestra_5)

LOTE_50 = rbinom(50 , 1 , .1)
muestra_6 = muestra_aleatoria_aleatoria(LOTE_50) ## Se tomo el 17% de la muestra
prop.table(table(muestra_6))
## muestra_6
##         0         1 
## 0.9047619 0.0952381
shapiro.test(muestra_6)
## 
##  Shapiro-Wilk normality test
## 
## data:  muestra_6
## W = 0.34141, p-value = 9.399e-09
qqnorm(muestra_6)

LOTE_60 = rbinom(60 , 1 , .1)
muestra_7 = muestra_aleatoria_aleatoria(LOTE_60) ## Se tomo el 17% de la muestra
prop.table(table(muestra_7))
## muestra_7
##          0          1 
## 0.91666667 0.08333333
shapiro.test(muestra_7)
## 
##  Shapiro-Wilk normality test
## 
## data:  muestra_7
## W = 0.32693, p-value = 1.207e-06
qqnorm(muestra_7)

LOTE_100 = rbinom(100 , 1 , .1)
muestra_8 = muestra_aleatoria_aleatoria(LOTE_100) ## Se tomo el 17% de la muestra
prop.table(table(muestra_8))
## muestra_8
## 0 
## 1
#shapiro.test(muestra_8)
qqnorm(muestra_8)

LOTE_200 = rbinom(200 , 1 , .1)
muestra_9 = muestra_aleatoria_aleatoria(LOTE_200) ## Se tomo el 17% de la muestra
prop.table(table(muestra_9))
## muestra_9
##         0         1 
## 0.8571429 0.1428571
shapiro.test(muestra_9)
## 
##  Shapiro-Wilk normality test
## 
## data:  muestra_9
## W = 0.41588, p-value < 2.2e-16
qqnorm(muestra_9)

LOTE_500 = rbinom(500 , 1 , .1)
muestra_10 = muestra_aleatoria_aleatoria(LOTE_500) ## Se tomo el 17% de la muestra
prop.table(table(muestra_10))
## muestra_10
##         0         1 
## 0.8985507 0.1014493
shapiro.test(muestra_10)
## 
##  Shapiro-Wilk normality test
## 
## data:  muestra_10
## W = 0.34486, p-value < 2.2e-16
qqnorm(muestra_10)

## Con 90%
set.seed(1944)
LOTE_5 = rbinom(5 , 1 , .1)
muestra_1 = muestra_aleatoria_aleatoria(LOTE_5) ## Se tomo el 17% de la muestra
prop.table(table(muestra_1))
## muestra_1
## 0 
## 1
#shapiro.test(muestra_1)
qqnorm(muestra_1)

LOTE_10 = rbinom(10 , 1 , .1)
muestra_2 = muestra_aleatoria_aleatoria(LOTE_10) ## Se tomo el 17% de la muestra
prop.table(table(muestra_2))
## muestra_2
## 0 
## 1
#shapiro.test(muestra_2)
qqnorm(muestra_2)

LOTE_15 = rbinom(15 , 1 , .1)
muestra_3 = muestra_aleatoria_aleatoria(LOTE_15) ## Se tomo el 17% de la muestra
prop.table(table(muestra_3))
## muestra_3
## 0 
## 1
#shapiro.test(muestra_3)
qqnorm(muestra_3)

LOTE_20 = rbinom(20 , 1 , .1)
muestra_4 = muestra_aleatoria_aleatoria(LOTE_20) ## Se tomo el 17% de la muestra
prop.table(table(muestra_4))
## muestra_4
##         0         1 
## 0.8823529 0.1176471
shapiro.test(muestra_4)
## 
##  Shapiro-Wilk normality test
## 
## data:  muestra_4
## W = 0.38525, p-value = 1.586e-07
qqnorm(muestra_4)

LOTE_30 = rbinom(30 , 1 , .1)
muestra_5 = muestra_aleatoria_aleatoria(LOTE_30) ## Se tomo el 17% de la muestra
prop.table(table(muestra_5))
## muestra_5
##          0          1 
## 0.95454545 0.04545455
shapiro.test(muestra_5)
## 
##  Shapiro-Wilk normality test
## 
## data:  muestra_5
## W = 0.22147, p-value = 7.417e-10
qqnorm(muestra_5)

LOTE_50 = rbinom(50 , 1 , .1)
muestra_6 = muestra_aleatoria_aleatoria(LOTE_50) ## Se tomo el 17% de la muestra
prop.table(table(muestra_6))
## muestra_6
## 0 
## 1
#shapiro.test(muestra_6)
qqnorm(muestra_6)

LOTE_60 = rbinom(60 , 1 , .1)
muestra_7 = muestra_aleatoria_aleatoria(LOTE_60) ## Se tomo el 17% de la muestra
prop.table(table(muestra_7))
## muestra_7
##         0         1 
## 0.8571429 0.1428571
shapiro.test(muestra_7)
## 
##  Shapiro-Wilk normality test
## 
## data:  muestra_7
## W = 0.41806, p-value = 1.207e-10
qqnorm(muestra_7)

LOTE_100 = rbinom(100 , 1 , .1)
muestra_8 = muestra_aleatoria_aleatoria(LOTE_100) ## Se tomo el 17% de la muestra
prop.table(table(muestra_8))
## muestra_8
##         0         1 
## 0.8214286 0.1785714
#shapiro.test(muestra_8)
qqnorm(muestra_8)

LOTE_200 = rbinom(200 , 1 , .1)
muestra_9 = muestra_aleatoria_aleatoria(LOTE_200) ## Se tomo el 17% de la muestra
prop.table(table(muestra_9))
## muestra_9
##          0          1 
## 0.90140845 0.09859155
shapiro.test(muestra_9)
## 
##  Shapiro-Wilk normality test
## 
## data:  muestra_9
## W = 0.34016, p-value = 3.067e-16
qqnorm(muestra_9)

LOTE_500 = rbinom(500 , 1 , .1)
muestra_10 = muestra_aleatoria_aleatoria(LOTE_500) ## Se tomo el 17% de la muestra
prop.table(table(muestra_10))
## muestra_10
##         0         1 
## 0.8902954 0.1097046
shapiro.test(muestra_10)
## 
##  Shapiro-Wilk normality test
## 
## data:  muestra_10
## W = 0.36041, p-value < 2.2e-16
qqnorm(muestra_10)

Problema 4

Estimación boostrap

Cuando se extrae una muestra de una población que no es normal y se requiere estimar un intervalo de confianza se pueden utilizar los métodos de estimación bootstrap. Esta metodología supone que se puede reconstruir la población objeto de estudio mediante un muestreo con reemplazo de la muestra que se tiene. Existen varias versiones del método. Una presentación básica del método se describe a continuación:

El artículo de In-use Emissions from Heavy Duty Dissel Vehicles (J.Yanowitz, 2001) presenta las mediciones de eficiencia de combustible en millas/galón de una muestra de siete camiones. Los datos obtenidos son los siguientes: \(7.69, 4.97, 4.56, 6.49, 4.34, 6.24 y 4.45\). Se supone que es una muestra aleatoria de camiones y que se desea construir un intervalo de confianza del \(95\%\) para la media de la eficiencia de combustible de esta población. No se tiene información de la distribución de los datos. El método bootstrap permite construir intervalos de confianza del \(95\%\) - Para ilustrar el método suponga que coloca los valores de la muestra en una caja y extrae uno al azar. Este correspondería al primer valor de la muestra bootstrap \(X^{∗}_{1}\). Después de anotado el valor se regresa X∗1 a la caja y se extrae el valor \(X^{∗}_{2}\), regresandolo nuevamente. Este procedimiento se repite hasta completar una muestra de tamaño \(n\), \(X^{∗}_{1}\), \(X^{∗}_{2}\),\(X^{∗}_{3}\),\(X^{∗}_{n}\), conformando la muestra bootstrap.

Es necesario extraer un gran número de muestras (suponga \(k = 1000\)). Para cada una de las muestra bootstrap obtenidas se calcula la media {x^{*}_{i}} , obteniéndose un valor para cada muestra. El intervalo de confianza queda conformado por los percentiles \(P2.5\) y \(P97.5\). Existen dos métodos para estimarlo:

Métodos Formula
Método 1 \((P_{2.5};P_{97.5})\)
Método 2 \((2\bar{x}-P_{97.5};2\bar{x}-P_{2.5})\)

Construya el intervalo de confianza por los dos métodos y compare los resultados obtenidos. Comente los resultados. Confiaría en estas estimaciones?

# Estimar los intervalos de confianza para la media de una muestra de datos

x = c(7.69, 4.97, 6.49, 4.34, 6.24, 4.45) # Datos de la muestra
n = 1000

muest_bootstrap = matrix(NA, nrow = n, ncol = length(x))

set.seed(1015)
for (i in 1:n) {
  muestra_bootstrap = sample(x, replace = TRUE)
  muest_bootstrap[i, ] = muestra_bootstrap
}

medias_x_muestra = apply(muest_bootstrap, 1, mean)

metodo1 = quantile(medias_x_muestra, probs = c(0.025, 0.975)) #Se calcula el IC para el método 1
metodo1
##     2.5%    97.5% 
## 4.778167 6.678625
metodo2 = c(2 * mean(medias_x_muestra) - metodo1[2], 2 * mean(medias_x_muestra) - metodo1[1]) #Se calcula el IC para el método 2
metodo2
##    97.5%     2.5% 
## 4.715962 6.616420
# Calcular intervalo de confianza para la media.
media_metodos = c(mean(medias_x_muestra) - sd(medias_x_muestra), mean(medias_x_muestra) + sd(medias_x_muestra))

# Graficar el histograma con líneas para intervalos de confianza

hist(medias_x_muestra, 
     main = "Histograma con intervalos de confianza y medida", 
     xlab = "Media Bootstrap", 
     ylab = "Frecuencia", col = "grey")

abline(v = metodo1, col = "red", lwd = 2, lty = 2)
abline(v = metodo2, col = "blue", lwd = 2, lty = 2)
abline(v = mean(medias_x_muestra), col = "green", lwd = 2, lty = 2)

legend("topright", legend = c("Método 1", "Método 2", "Media"), col = c("red", "blue", "green"), lty = 2, lwd = 2)