Problema 1

Estimación del valor de π


La siguiente figura sugiere como estimar el valor de π con una simulación. En la figura, un circuito con un área igual a π/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 π/4. Por tanto, se puede estimar el valor de π/4 al contar el número de puntos dentro del círculo, para obtener la estimación de π/4. De este último resultado se encontrar una aproximación para el valor de π.

Pasos sugeridos:
  1. Genere n coordenadas x: X1,. . ., Xn. 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)
  2. Genere 1000 coordenadas y : Y1,…,Yn, utilizando nuevamente la distribución uniforme con valor mínimo de 0 y valor máximo de 1.
  3. Cada punto (Xi,Yi) se encuentra dentro del círculo si su distancia desde el centro (0.5,0.5) es menor a 0.5. Para cada par (Xi,Yi) determine si la distancia desde el centro es menor a 0.5. Esto último se puede realizar al calcular el valor (Xi−0.5)2 + (Yi−0.5)2, que es el cuadrado de la distancia, y al determinar si es menor que 0.25.
  4. ¿Cuántos de los puntos están dentro del círculo? ¿Cuál es su estimación de π?

Desarrollo de Problema 1



Para desarrollar este problema empezamos generando puntos usando una distribución uniforme, tanto para la coordenadas “X” como para las coordenadas “Y”, para evitar que cada vez que se realice el problema se genere una iteración diferente declaramos una semilla diferente (función “set.seed()”) para que los resultados sean reproducibles cada vez. 
## [1] "Mostrando los 10 primeros números de cada vector obtenido"
##  [1] 0.26550866 0.37212390 0.57285336 0.90820779 0.20168193 0.89838968
##  [7] 0.94467527 0.66079779 0.62911404 0.06178627
##  [1] 0.1848823 0.7023740 0.5733263 0.1680519 0.9438393 0.9434750 0.1291590
##  [8] 0.8334488 0.4680185 0.5499837
#Código usado para generar los valores:
set.seed(1)
x=runif(1000, min = 0, max = 1)
set.seed(2)
y=runif(1000, min = 0, max = 1)
Tras ello procedemos a calcular la distancia de cada punto (x,y) respecto del centro (0.5,0.5), con lo que obtenemos un vector “z” que representará la distancia de todos los puntos generados. 
## [1] "Mostrando los 10 primeros números del vector distancia"
##  [1] 0.3927918 0.2393899 0.1033652 0.5261398 0.5347776 0.5961412 0.5790157
##  [8] 0.3701946 0.1330160 0.4410552
#Código usado para generar los valores:
distancia <- function(x,y)sqrt(((x-0.5)^2)+((y-0.5)^2))
z= distancia(x,y)
Para hallar π tenemos que la proporción de puntos es equivalente al área del círculo inscrito en el cuadrado, que es π/4. Entonces para ello contamos los puntos cuya distancia es menor a 0.5 respecto del centro y lo dividimos por 1000 para hallar la proporción, finalmente despejamos π de la ecuación para hallar la estimación por este método. 
## [1] "La proporción de puntos dentro del cículo es: 0.769"
## [1] "La estimación de π es: 3.076"
## [1] "El error de estimación de π es: 0.0208787900986598"
#Contamos cuantos puntos estan dentro del circulo, es decir menores a 0.5
puntos<-sum(z<=0.5)

#Calulamos la probabilidad a partir de los resultados obtenidos
p<-puntos/1000
print(paste("La proporción de puntos dentro del cículo es:",p))

#Si π/4=p, entonces π es igual a:
pi_valor <- p*4
print(paste("La estimación de π es:",pi_valor))

#El error de este metodo de estimación fue igual a:
error_estimacion = (pi - pi_valor)/pi
print(paste("El error de estimación de π es:",error_estimacion))
Con ello tenemos que para una muestra de 1,000 puntos la proporción de puntos se ubicó en 76,9%, con ello el valor de π obtenido fue de 3.076…, y con ello con respecto al verdadero valor de π nuestra estimación tuvo un error de 2.09%. Como ya tenemos las funciones base podemos crear una función que abarque las anteriores y cuyo valor de ingreso sea el número de muestras para ver como mejora la estimación a medida que aumenta el número de muestras, tomamos un n de 10,000 y de 100,000 como referencia. Con un n de 10,000: 
## [1] "La proporción de puntos dentro del cículo es: 0.7811"
## [1] "La estimación de π es: 3.1244"
## [1] "El error de estimación de π es: 0.00547259160736437"
#Crear función que abarque todo el proceso de cálculo de proporción
proporcion_puntos <- function(n){
  set.seed(1)
  x=runif(n, min = 0, max = 1)
  set.seed(2)
  y=runif(n, min = 0, max = 1)
  z= distancia(x,y)
  return(sum(z<=0.5)/n)}

#Calculando la estimación para n=10,000
p<-proporcion_puntos(10000)
print(paste("La proporción de puntos dentro del cículo es:",p))
pi_valor <- p*4
print(paste("La estimación de π es:",pi_valor))
error_estimacion = (pi - pi_valor)/pi
print(paste("El error de estimación de π es:",error_estimacion))
Para un n=10,000 tenemos una proporción de 78.11%, con ello logramos una estimación de π de 3.124…, esto nos genera un error con respecto al valor verdadero de π de solo 0.54%. Ahora revisemos con valor de n=100,000. 
## [1] "La proporción de puntos dentro del cículo es: 0.78323"
## [1] "La estimación de π es: 3.13292"
## [1] "El error de estimación de π es: 0.00276059137707851"
#Calculando la estimación para n=100,000
p<-proporcion_puntos(100000)
print(paste("La proporción de puntos dentro del cículo es:",p))
pi_valor <- p*4
print(paste("La estimación de π es:",pi_valor))
error_estimacion = (pi - pi_valor)/pi
print(paste("El error de estimación de π es:",error_estimacion))
Finalmente para un n=100,000 tenemos una proporción de 78.32%, la estimación de π es ahora 3.133…, con ello el error respecto del valor verdadero de π se ha reducido a tan solo 0.27%.

Problema 2

Propiedades de los estimadores


La simulación ayuda a entender y validar 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 X1, X2, X3 y X4, una muestra aleatoria de tamaño 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:

Desarrollo de Problema 2


Para desarrollar este problema empezamos con el desarrollo del caso base para una muestra pequeña, y posteriormente con la formulación ya planteada generalizamos para casos con muestras más grandes, que es el objetivo de este problema, para ello asumiremos como parámetro (1/λ)=2. Primero empecemos con una muestra de tamaño 20. 
## [1] "Mostrando los 10 primeros números de cada vector obtenido"
##  [1] 0.37759092 0.59082139 0.07285336 0.06989763 0.21803431 1.44748427
##  [7] 0.61478103 0.26984142 0.47828375 0.07352300
##  [1] 0.93267622 0.20237404 0.07332633 0.86535486 0.04476309 0.33344882
##  [7] 0.53718343 0.75581465 0.65713795 0.07826514
##  [1] 0.865313286 0.307516399 0.616458293 0.502042225 0.102100675 0.104394054
##  [7] 1.141821671 0.005023914 0.034035353 0.057249436
##  [1] 0.08580031 2.15197247 0.43405281 0.40132351 0.31357421 0.36742913
##  [7] 0.22440589 0.02534994 0.25467503 0.41857483
#Generando vectores de 20 variables para cada coordenada
set.seed(1)
x1=rexp(20,2)
set.seed(2)
x2=rexp(20,2)
set.seed(3)
x3=rexp(20,2)
set.seed(4)
x4=rexp(20,2)

print("Mostrando los 10 primeros números de cada vector obtenido")
show(x1[1:10])
show(x2[1:10])
show(x3[1:10])
show(x4[1:10])
Creamos las correspondientes funciones para cada parámetro y calculamos el vector de 20 muestras de cada uno de los parámetros, luego procedemos a revisar las características de los parámetros para ver con el actual tamaño de muestra como es la insesgadez y la eficiencia de los parámetros.
Tabla 1. Comparación entre los parámetros para n=20
θ1 θ2 θ3 θ4
Promedio 0.4817846 0.9590842 0.5023997 0.5585915
Desviación 0.2759010 0.5337224 0.2695575 0.3088400 
#Creando funciones para los parámetros:
parametro1<- function(x1,x2,x3,x4){((x1+x2)/6)+((x3+x4)/3)}
parametro2<- function(x1,x2,x3,x4){(x1+2*(x2)+3*(x3)+4*(x4))/5}
parametro3<- function(x1,x2,x3,x4){(x1+x2+x3+x4)/4}
parametro4<- function(x1,x2,x3,x4){(pmin(x1,x2,x3,x4)+pmax(x1,x2,x3,x4))/2}

#Calculando parámetros:
θ1<-parametro1(x1,x2,x3,x4)
θ2<-parametro2(x1,x2,x3,x4)
θ3<-parametro3(x1,x2,x3,x4)
θ4<-parametro4(x1,x2,x3,x4)

#Mostrando el comportamiento de los cuatro parámetros
boxplot(θ1,θ2,θ3,θ4, names = c("θ1","θ2","θ3","θ4"), main="Gráfico 1: Comparación entre los parámetros para n=20")
abline(h=0.5)

#Creando tabla para mostrar las características de los parámetros
datos20<-t(data.frame("Promedio"=c(mean(θ1),mean(θ2),mean(θ3),mean(θ4)),"Desviación"=c(sd(θ1),sd(θ2),sd(θ3),sd(θ4))))
colnames(datos20) <- c("θ1","θ2","θ3","θ4")

library(kableExtra)
kbl(datos20, caption = "Tabla 1. Comparación entre los parámetros para n=20")%>%
  kable_classic(full_width = F)
Podemos ver los parámetros tienden a 0.5, solo el parámetro θ2 está fallando en la propiedad de insesgadez, ya que la media de este no está apuntando de la misma forma, respecto a la eficiencia θ1 y θ3 tienen las menores desviaciones, siendo θ3 el estimador de mejor comportamiento en ambos atributos.


Para probar la propiedad de consistencia, necesitamos hacer pruebas en tamaños de muestra mayores, y de paso validar que las propiedades anteriormente evaluadas se mantengan, creamos una función cuya variable sea n, y como ya tenemos las funciones anteriores solo es comprimir buena parte de la formulación anterior en una única función y ahora evaluamos con n=50, 100 y 1000. También guardamos los vectores de cada parámetro obtenido anteriormente para realizar comparaciones más adelante. 
#Almacenando los parámetros para n=20
θ_n20<-data.frame(θ1,θ2,θ3,θ4)

#Crear función que abarque todo el proceso de cálculo de proporción
Calc_parametros <- function(n){
  set.seed(1)
  x1=rexp(n,2)
  set.seed(2)
  x2=rexp(n,2)
  set.seed(3)
  x3=rexp(n,2)
  set.seed(4)
  x4=rexp(n,2)
  parametro1<- function(x1,x2,x3,x4){((x1+x2)/6)+((x3+x4)/3)}
  parametro2<- function(x1,x2,x3,x4){(x1+2*(x2)+3*(x3)+4*(x4))/5}
  parametro3<- function(x1,x2,x3,x4){(x1+x2+x3+x4)/4}
  parametro4<- function(x1,x2,x3,x4){(pmin(x1,x2,x3,x4)+pmax(x1,x2,x3,x4))/2}
  θ1<-parametro1(x1,x2,x3,x4)
  θ2<-parametro2(x1,x2,x3,x4)
  θ3<-parametro3(x1,x2,x3,x4)
  θ4<-parametro4(x1,x2,x3,x4)
  return(θ_n<-data.frame(θ1,θ2,θ3,θ4))}

#Calculando la estimación para n=50, 100 y 1000
θ_n50<-Calc_parametros(50)
θ_n100<-Calc_parametros(100)
θ_n1000<-Calc_parametros(1000)
Con la información ya calculada, ahora procedemos a visualizar la información para ver cómo se han comportado los parámetros en cada uno de los casos.
Tabla 2. Comparación entre los parámetros para n=20
θ1 θ2 θ3 θ4
Promedio 0.4817846 0.9590842 0.5023997 0.5585915
Desviación 0.2759010 0.5337224 0.2695575 0.3088400
Tabla 3. Comparación entre los parámetros para n=50
θ1 θ2 θ3 θ4
Promedio 0.5059073 1.0133409 0.5106041 0.5908491
Desviación 0.2843369 0.5917188 0.2675514 0.3374944
Tabla 4. Comparación entre los parámetros para n=100
θ1 θ2 θ3 θ4
Promedio 0.5014690 0.9959711 0.5034396 0.5803656
Desviación 0.2504838 0.5177877 0.2418094 0.2964896
Tabla 5. Comparación entre los parámetros para n=1000
θ1 θ2 θ3 θ4
Promedio 0.5078035 1.0172257 0.5100604 0.5913347
Desviación 0.2688620 0.5574844 0.2566013 0.3151266 


#Mostrando el comportamiento de los cuatro parámetros para n=20, 50, 100 y 1000
par(mfrow = c(1, 2))
boxplot(θ_n20, main="Gráfico 2: Comparación entre los parámetros para n=20", cex.main=0.65)
abline(h=0.5)
boxplot(θ_n50, main="Gráfico 3: Comparación entre los parámetros para n=50", cex.main=0.65)
abline(h=0.5)

par(mfrow = c(1, 2))
boxplot(θ_n100, main="Gráfico 4: Comparación entre los parámetros para n=100", cex.main=0.65)
abline(h=0.5)
boxplot(θ_n1000, main="Gráfico 5: Comparación entre los parámetros para n=1000", cex.main=0.65)
abline(h=0.5)

#Creando tabla para mostrar las características de los parámetros
datos50<-t(data.frame("Promedio"=c(mean(θ_n50$θ1),mean(θ_n50$θ2),mean(θ_n50$θ3),mean(θ_n50$θ4)),"Desviación"=c(sd(θ_n50$θ1),sd(θ_n50$θ2),sd(θ_n50$θ3),sd(θ_n50$θ4))))
colnames(datos50) <- c("θ1","θ2","θ3","θ4")
datos100<-t(data.frame("Promedio"=c(mean(θ_n100$θ1),mean(θ_n100$θ2),mean(θ_n100$θ3),mean(θ_n100$θ4)),"Desviación"=c(sd(θ_n100$θ1),sd(θ_n100$θ2),sd(θ_n100$θ3),sd(θ_n100$θ4))))
colnames(datos100) <- c("θ1","θ2","θ3","θ4")
datos1000<-t(data.frame("Promedio"=c(mean(θ_n1000$θ1),mean(θ_n1000$θ2),mean(θ_n1000$θ3),mean(θ_n1000$θ4)),"Desviación"=c(sd(θ_n1000$θ1),sd(θ_n1000$θ2),sd(θ_n1000$θ3),sd(θ_n1000$θ4))))
colnames(datos1000) <- c("θ1","θ2","θ3","θ4")

#Renderizando Tablas
kbl(datos20, caption = "<center><b><big>Tabla 2. Comparación entre los parámetros para n=20</big></b></center>")%>%
  kable_classic(full_width = FALSE, position = "float_left")
kbl(datos50, caption = "<center><b><big>Tabla 3. Comparación entre los parámetros para n=50</big></b></center>")%>%
  kable_classic(full_width = FALSE, position = "center")
kbl(datos100, caption = "<center><b><big>Tabla 4. Comparación entre los parámetros para n=100</big></b></center>")%>%
  kable_classic(full_width = FALSE, position = "float_left")
kbl(datos1000, caption = "<center><b><big>Tabla 5. Comparación entre los parámetros para n=1000</big></b></center>")%>%
  kable_classic(full_width = FALSE, position = "center")
Con la información ahora visible podemos sacar las siguientes conclusiones:
  • El parámetro θ2 es el único que está alejado de la media, pues en los cuatro tamaños de muestra su rango medio estuvo entre 0.95 y 1.01, en tanto que otros parámetros se mantuvieron cercanos a 0.5.
  • El parámetro con mejor comportamiento en insesgadez es θ1, seguido de θ3.
  • En cuanto al parámetro más eficiente, la menor desviación le correspondió a θ3 seguido de θ1.
  • Los cuatro parámetros fueron bastante consistentes con su comportamiento, hubo más bien poca variación en la media y la desviación estándar, desde la gráfica se puede apreciar cierta mejoría en el rango entre cuartiles en los parámetros θ1 y θ4.
  • El parámetro θ3 fue el más eficiente y tuvo un buen comportamiento respecto de la insesgadez, por lo que se le podría considerar como el mejor parámetro evaluado.

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) Realice una simulación

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%.

lote = c(rep("enferma",500),rep("sana",500))
parametro = mean(lote=="enferma")
parametro
## [1] 0.5

b) Genere una función

Genere una función que permita: obtener una muestra aleatoria de la población y calcule el estimador de la proporción muestral pˆpara un tamaño de muestra dado n

# Definir la población y el parámetro global
lote <- c(rep("enferma", 500), rep("sana", 500))
parametro <- mean(lote == "enferma")

# Definir la función para obtener una muestra aleatoria y calcular el estimador de la proporción muestral
obtener_muestra_y_estimar_proporcion <- function(n) {
  # Obtener una muestra aleatoria de la población
  muestra <- sample(lote, size = n)
  
  # Calcular el estimador de la proporción muestral
  proporcion_muestral <- mean(muestra == "enferma")
  
  return(proporcion_muestral)
}

# Ejemplo de uso de la función
tamaño_muestra <- 300
estimacion_proporcion <- obtener_muestra_y_estimar_proporcion(tamaño_muestra)
estimacion_proporcion
## [1] 0.48

c) Repita el escenario anterior(b) 500 veces

Repita el escenario anterior (b) n=500 veces y analice los resultados en cuanto al comportamiento de los 500 resultados del estimador 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

# Definir la función modificada
obtener_muestras_y_estimar_proporcion <- function(n, repeticiones) {
  # Inicializar un vector para almacenar los resultados de las estimaciones
  estimaciones <- numeric(repeticiones)
  
  # Repetir el proceso de obtener una muestra y estimar la proporción 'repeticiones' veces
  for (i in 1:repeticiones) {
    # Obtener una muestra aleatoria de la población
    muestra <- sample(lote, size = n, replace = FALSE)
  
    # Calcular el estimador de la proporción muestral
    proporcion_muestral <- mean(muestra == "enferma")
    
    # Almacenar el resultado de la estimación
    estimaciones[i] <- proporcion_muestral
  }
  
  return(estimaciones)
}

# Repetir el escenario 500 veces para n = 500
repeticiones <- 500
resultados <- obtener_muestras_y_estimar_proporcion(n = 500, repeticiones)

# Analizar los resultados
# Calculamos medidas de tendencia central y dispersión
media <- mean(resultados, na.rm = TRUE)
mediana <- median(resultados, na.rm = TRUE)
desviacion_estandar <- sd(resultados, na.rm = TRUE)


# Graficamos un histograma de los resultados
hist(resultados, breaks = 20, main = "Distribución de estimaciones de proporción muestral",
     xlab = "Proporción muestral", ylab = "Frecuencia")

# Mostramos medidas de tendencia central y dispersión en el gráfico
abline(v = media, col = "red", lwd = 2)
abline(v = mediana, col = "blue", lwd = 2)
legend("topright", legend = c("Media", "Mediana"),
       col = c("red", "blue"), lwd = 2)

# Imprimir resultados
cat("Media de las estimaciones:", media, "\n")
## Media de las estimaciones: 0.500976
cat("Mediana de las estimaciones:", mediana, "\n")
## Mediana de las estimaciones: 0.502
cat("Desviación estándar de las estimaciones:", desviacion_estandar, "\n")
## Desviación estándar de las estimaciones: 0.01671713

Análisis: El análisis de los resultados muestra que la media y mediana de los estimadores de proporción muestral se acercan al valor teórico de 0.50, lo que sugiere que el estimador es insesgado y eficiente. Además, se observa que tanto la media como la mediana se acercan al mismo valor, lo que garantiza poca variación en los resultados y de igual forma el valor de la desviación estándar se acerca a cero, lo cual indica una baja dispersión de los datos.

d) Repita los puntos b y c

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 :shspiro.test()) y métodos gráficos (gráfico de normalidad: qqnorm()). Comente en su informe los resultados obtenidos

# Definir una secuencia de tamaños de muestra
tamanos_muestra <- c(5, 10, 15, 20, 30, 50, 60, 100, 200, 500)

# Función para realizar las pruebas de bondad de ajuste y gráficos QQ
analizar_tamanos_muestra <- function(tamanos_muestra, poblacion) {
  for (n in tamanos_muestra) {
    # Obtener las estimaciones de proporción muestral para el tamaño de muestra actual
    estimaciones <- obtener_muestras_y_estimar_proporcion(n, repeticiones = 500)
    
    # Prueba de Shapiro-Wilk
    shapiro_result <- shapiro.test(estimaciones)
    
    # Gráfico QQ
    qqnorm(estimaciones, main = paste("QQ plot para n =", n))
    qqline(estimaciones, col = 2)
    
    # Imprimir resultados
    cat("Tamaño de muestra:", n, "\n")
    print(shapiro.test(estimaciones))
   
  }
}

# Ejecutar la función para analizar los diferentes tamaños de muestra
analizar_tamanos_muestra(tamanos_muestra, lote)

## Tamaño de muestra: 5 
## 
##  Shapiro-Wilk normality test
## 
## data:  estimaciones
## W = 0.91982, p-value = 1.186e-15

## Tamaño de muestra: 10 
## 
##  Shapiro-Wilk normality test
## 
## data:  estimaciones
## W = 0.96614, p-value = 2.486e-09

## Tamaño de muestra: 15 
## 
##  Shapiro-Wilk normality test
## 
## data:  estimaciones
## W = 0.97049, p-value = 1.721e-08

## Tamaño de muestra: 20 
## 
##  Shapiro-Wilk normality test
## 
## data:  estimaciones
## W = 0.98158, p-value = 5.809e-06

## Tamaño de muestra: 30 
## 
##  Shapiro-Wilk normality test
## 
## data:  estimaciones
## W = 0.98265, p-value = 1.117e-05

## Tamaño de muestra: 50 
## 
##  Shapiro-Wilk normality test
## 
## data:  estimaciones
## W = 0.99014, p-value = 0.001969

## Tamaño de muestra: 60 
## 
##  Shapiro-Wilk normality test
## 
## data:  estimaciones
## W = 0.98713, p-value = 0.0002133

## Tamaño de muestra: 100 
## 
##  Shapiro-Wilk normality test
## 
## data:  estimaciones
## W = 0.99345, p-value = 0.02886

## Tamaño de muestra: 200 
## 
##  Shapiro-Wilk normality test
## 
## data:  estimaciones
## W = 0.99653, p-value = 0.3556

## Tamaño de muestra: 500 
## 
##  Shapiro-Wilk normality test
## 
## data:  estimaciones
## W = 0.99532, p-value = 0.138

Análisis: Al calcular los estimadores ^P cambiando el tamaño de la muestra se puede observar como ^P tiende a estar más cerca del valor del parámetro P = 0.5.Cuando nos apoyamos en las gráficas qq-norm, observamos que en la mayoría de las muestras de tamaños grandes, los puntos (observaciones) están muy cercanos a la línea nominal para una distribución normal, hecho que no ocurre con las muestras de tamaños pequeños.

e) Repita toda la simulación (puntos a - d)

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.

# Función para generar la población con un porcentaje de plantas enfermas dado
generar_poblacion <- function(tamano_lote, porcentaje_enfermas) {
  cantidad_enfermas <- round(tamano_lote * porcentaje_enfermas)
  poblacion <- c(rep(1, cantidad_enfermas), rep(0, tamano_lote - cantidad_enfermas))
  return(poblacion)
}

# Función para obtener una muestra aleatoria y calcular el estimador de la proporción muestral
obtener_muestras_y_estimar_proporcion <- function(poblacion, n, repeticiones) {
  estimaciones <- numeric(repeticiones)
  
  for (i in 1:repeticiones) {
    muestra <- sample(poblacion, size = n, replace = FALSE)
    proporcion_muestral <- mean(muestra)
    estimaciones[i] <- proporcion_muestral
  }
  
  return(estimaciones)
}

# Función para realizar la simulación completa
simulacion_completa <- function(porcentaje_enfermas) {
  # Generar la población
  poblacion <- generar_poblacion(tamano_lote = 1000, porcentaje_enfermas)
  
  # Definir tamaños de muestra
  tamanos_muestra <- c(5, 10, 15, 20, 30, 50, 60, 100, 200, 500)
  
  # Repetir el proceso de obtener estimaciones de proporción muestral para diferentes tamaños de muestra
  for (n in tamanos_muestra) {
    # Obtener estimaciones de proporción muestral
    estimaciones <- obtener_muestras_y_estimar_proporcion(poblacion, n, repeticiones = 500)
    
    # Prueba Shapiro-Wilk
    shapiro_result <- shapiro.test(estimaciones)
    
    # Gráfico QQ
    qqnorm(estimaciones, main = paste("QQ plot para n =", n))
    qqline(estimaciones, col = 2)
    
    # Imprimir resultados
    cat("Porcentaje de plantas enfermas:", porcentaje_enfermas, "\n")
    cat("Tamaño de muestra:", n, "\n")
    print(shapiro.test(estimaciones))
    ##cat("Prueba Shapiro-Wilk p-valor:", shapiro_result$p.value, "\n\n")
  }
}

# Realizar la simulación para lotes con 10% de plantas enfermas
simulacion_completa(porcentaje_enfermas = 0.1)

## Porcentaje de plantas enfermas: 0.1 
## Tamaño de muestra: 5 
## 
##  Shapiro-Wilk normality test
## 
## data:  estimaciones
## W = 0.72146, p-value < 2.2e-16

## Porcentaje de plantas enfermas: 0.1 
## Tamaño de muestra: 10 
## 
##  Shapiro-Wilk normality test
## 
## data:  estimaciones
## W = 0.8474, p-value < 2.2e-16

## Porcentaje de plantas enfermas: 0.1 
## Tamaño de muestra: 15 
## 
##  Shapiro-Wilk normality test
## 
## data:  estimaciones
## W = 0.88974, p-value < 2.2e-16

## Porcentaje de plantas enfermas: 0.1 
## Tamaño de muestra: 20 
## 
##  Shapiro-Wilk normality test
## 
## data:  estimaciones
## W = 0.92924, p-value = 1.272e-14

## Porcentaje de plantas enfermas: 0.1 
## Tamaño de muestra: 30 
## 
##  Shapiro-Wilk normality test
## 
## data:  estimaciones
## W = 0.95607, p-value = 4.774e-11

## Porcentaje de plantas enfermas: 0.1 
## Tamaño de muestra: 50 
## 
##  Shapiro-Wilk normality test
## 
## data:  estimaciones
## W = 0.97478, p-value = 1.374e-07

## Porcentaje de plantas enfermas: 0.1 
## Tamaño de muestra: 60 
## 
##  Shapiro-Wilk normality test
## 
## data:  estimaciones
## W = 0.97838, p-value = 9.237e-07

## Porcentaje de plantas enfermas: 0.1 
## Tamaño de muestra: 100 
## 
##  Shapiro-Wilk normality test
## 
## data:  estimaciones
## W = 0.98758, p-value = 0.0002931

## Porcentaje de plantas enfermas: 0.1 
## Tamaño de muestra: 200 
## 
##  Shapiro-Wilk normality test
## 
## data:  estimaciones
## W = 0.99248, p-value = 0.01289

## Porcentaje de plantas enfermas: 0.1 
## Tamaño de muestra: 500 
## 
##  Shapiro-Wilk normality test
## 
## data:  estimaciones
## W = 0.99049, p-value = 0.002591
# Realizar la simulación para lotes con 90% de plantas enfermas
simulacion_completa(porcentaje_enfermas = 0.9)

## Porcentaje de plantas enfermas: 0.9 
## Tamaño de muestra: 5 
## 
##  Shapiro-Wilk normality test
## 
## data:  estimaciones
## W = 0.71475, p-value < 2.2e-16

## Porcentaje de plantas enfermas: 0.9 
## Tamaño de muestra: 10 
## 
##  Shapiro-Wilk normality test
## 
## data:  estimaciones
## W = 0.85937, p-value < 2.2e-16

## Porcentaje de plantas enfermas: 0.9 
## Tamaño de muestra: 15 
## 
##  Shapiro-Wilk normality test
## 
## data:  estimaciones
## W = 0.89861, p-value < 2.2e-16

## Porcentaje de plantas enfermas: 0.9 
## Tamaño de muestra: 20 
## 
##  Shapiro-Wilk normality test
## 
## data:  estimaciones
## W = 0.92118, p-value = 1.65e-15

## Porcentaje de plantas enfermas: 0.9 
## Tamaño de muestra: 30 
## 
##  Shapiro-Wilk normality test
## 
## data:  estimaciones
## W = 0.95791, p-value = 9.393e-11

## Porcentaje de plantas enfermas: 0.9 
## Tamaño de muestra: 50 
## 
##  Shapiro-Wilk normality test
## 
## data:  estimaciones
## W = 0.97438, p-value = 1.126e-07

## Porcentaje de plantas enfermas: 0.9 
## Tamaño de muestra: 60 
## 
##  Shapiro-Wilk normality test
## 
## data:  estimaciones
## W = 0.97132, p-value = 2.532e-08

## Porcentaje de plantas enfermas: 0.9 
## Tamaño de muestra: 100 
## 
##  Shapiro-Wilk normality test
## 
## data:  estimaciones
## W = 0.98485, p-value = 4.523e-05

## Porcentaje de plantas enfermas: 0.9 
## Tamaño de muestra: 200 
## 
##  Shapiro-Wilk normality test
## 
## data:  estimaciones
## W = 0.99126, p-value = 0.004747

## Porcentaje de plantas enfermas: 0.9 
## Tamaño de muestra: 500 
## 
##  Shapiro-Wilk normality test
## 
## data:  estimaciones
## W = 0.99234, p-value = 0.0114

Análisis: Cuando se modificó la población para generar el Lote A y el Lote B con proporciones de 10% y 90% de individuos enfermos respectivamente, podemos observar que los estimadores ^P para ambos casos están muy cercanos al valor del parámetro P = 0.1 y P = 0.9 respectivamente.

En cuanto a la aplicación de la prueba de normalidad Shapiro-Wilk para muestras de tamaño n > 200, se observa tanto la aceptación como el rechazo de la hipótesis nula H0. Este hecho se respalda con el gráfico qq-plot de normalidad, donde las observaciones se sobreponen cada vez más con muestras de tamaño n > 200 en la línea compuesta por los valores nominales, acercándose más a una distribución gaussiana.

Problema 4

Estimacció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∗2,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:

Intervalo de confianza

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

x=c(7.69, 4.97, 4.56, 6.49, 4.34, 6.24, 4.45) # datos de combustible m/gl de 7 camiones
boot=sample(x,7000,replace=TRUE)   # se extraen n x m muestras
b=matrix(boot,nrow=1000,ncol=7, byrow = TRUE)    # se construye matriz de n x m con k=1000
mx=apply(b,1,mean)                 # se calculan medias por fila
mean(mx)                           # se calculan media población
## [1] 5.505633

Cálculo con primer método

ic1=quantile(mx, probs=c(0.025, 0.975)) # Cálculo IC método 1
ic1
##     2.5%    97.5% 
## 4.732821 6.424929

Cálculo con segundo método

ic2=c(2*mean(mx)-ic1[2], 2*mean(mx)-ic1[1]) # Cálculo IC método 2
ic2
##    97.5%     2.5% 
## 4.586337 6.278444

Comente los resultados. Confiaría en estas estimaciones?

El primer método de cálculo de intervalo de confianza al 95%, representado por las líneas rojas, muestra un intervalo muy cercano al parámetro poblacional (media), con cuantiles bastante precisos. Por otro lado, el segundo método, ilustrado por las líneas verdes, presenta un intervalo ligeramente más amplio en el extremo inferior y más estrecho en el superior. Sin embargo, ambos métodos ofrecen estimaciones confiables, ya que se aproximan al valor real del parámetro poblacional sin diferencias significativas entre sí. Por lo tanto, si confiaría en estas estimaciones ya que ambos métodos son cercanos al valor real del parámetro y no difieren entre si de manera significativa.

# Crea el lienzo del gráfico
hist(mx, las=1, main=" ", ylab = " ", xlab = " ", col="lightblue")

# Agrega las líneas al gráfico
abline(v=ic1, col="red",lwd=2)
abline(v=ic2, col="green",lwd=2)