Actividad 2 Ejercicios estadísticos

Problema 1: Estimación del valor de π

La siguiente figura sugiere como estimar el valor de π con una simulación. En la figura, un circulo 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 π.

Solución 1

set.seed(1)  # Establece una semilla para la reproducibilidad
#distribución uniforme con valor mínimo de 0 y valor máximo de 1

x <- runif(1000, min = 0, max = 1)
y <- runif(1000, min = 0, max = 1)

x1 <- runif(10000, min = 0, max = 1)
y1 <- runif(10000, min = 0, max = 1)

# Generar coordenadas para (X a Xn) y (Y a Yn) y realizar la operacion (Xi−0.5)^2+(Yi−0.5)^2 

distancia<- (x - 0.5)^2 + (y - 0.5)^2
puntos_dentro_del_circulo <- distancia < 0.25
num_puntos_dentro_circulo <- sum(puntos_dentro_del_circulo)

distancia1<- (x1 - 0.5)^2 + (y1 - 0.5)^2
puntos_dentro_del_circulo1 <- distancia1 < 0.25
num_puntos_dentro_circulo1 <- sum(puntos_dentro_del_circulo1)

#Como es pi/4 se hace la operacion contraria dividido el numero de puntos totales

estimacionpi <- 4 * (num_puntos_dentro_circulo / 1000)
estimacionpi1 <- 4 * (num_puntos_dentro_circulo1 / 10000)
cat("Número de puntos dentro del círculo:", num_puntos_dentro_circulo, "\n")

## Número de puntos dentro del círculo: 770

## Número de puntos dentro del círculo: 770
cat("Estimación de π:", estimacionpi, "\n")

## Estimación de π: 3.08

cat("Número de puntos dentro del círculo:", num_puntos_dentro_circulo1, "\n")

## Número de puntos dentro del círculo: 7854

## Número de puntos dentro del círculo: 770
cat("Estimación de π:", estimacionpi1, "\n")

## Estimación de π: 3.1416

Grafica para 1000

library(ggplot2)
ggplot(data.frame(x = x, y = y, dentro = puntos_dentro_del_circulo), aes(x = x, y = y, color = dentro)) +
  geom_point() +
  scale_color_manual(values = c("red", "green")) +
  labs(title = "Puntos", x = "X", y = "Y") +
  theme_minimal()

Grafica para 10000

library(ggplot2)
ggplot(data.frame(x = x1, y = y1, dentro = puntos_dentro_del_circulo1), aes(x = x1, y = y1, color = dentro)) +
  geom_point() +
  scale_color_manual(values = c("red", "green")) +
  labs(title = "Puntos", x = "X", y = "Y") +
  theme_minimal()

Conclusión 1

Se evidencia existen diferentes métodos para aproximar el valor de Pi, una de las constantes matemáticas más importantes. Donde uno de ellos es el método de Monte Carlo, que utiliza técnicas aleatorias para estimar el valor de Pi donde tomamos un numero aleatorio para poder estimar y a medida que generamos más puntos aleatorios, nuestra aproximación de Pi se vuelve más precisa, Este método es especialmente útil en problemas donde las soluciones analíticas son difíciles de obtener y en problemas de simulación, optimización y toma de decisiones, entre otros.

Problema 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 X1, X2, X3 y X4, 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:

Image Description

Solución 2

set.seed(5)
#Semilla: fijar un valor inicial para el generador de números aleatorios
theta = 2                                    
#Suponemos un tetha = 2, ya que este es desconocido
n=4                                          
# n: Tamaño de la muestra
m=2000*n                                     
# m: # de replicas del experimento
replicas=rexp(m, rate = theta)                                            
# de replicas exponenciales
datos = data.frame(matrix(replicas, nrow = m, ncol = n, byrow = TRUE))    
# matriz de datos m x n

muestrasn <- function(data, n){
  set.seed(5)
  muestra <- datos[sample(nrow(datos), size=n), ]    
  #Calculo de la muestra n=20, 50, 100, 1000
  muestra
  #Calculo estimadores
  ϑ1 <- ((muestra$X1 + muestra$X2)/6) + ((muestra$X3 + muestra$X4)/3)
  ϑ2 <- ((muestra$X1 + (2*muestra$X2) + (3*muestra$X3) + (4*muestra$X4))/5)
  ϑ3 <- (((muestra$X1) + (muestra$X2) + (muestra$X3) + (muestra$X4))/4)
    Max = apply(muestra,1,max)                      
    #cálculo del valor máximo para las m muestras
    Min = apply(muestra,1,min)                      
    #cálculo del valor minimo para las m muestras
  ϑ4 <- (Min + Max)/2

  T1234=data.frame(ϑ1,ϑ2,ϑ3,ϑ4)
  T1234

  boxplot(T1234, col = c("#FF7256","#C1FFC1","#CD6090","#79CDCD"), ylim = c(0,4.3),
          main = "Estimadores θ" , ylab="Estimación")
  abline(h=0.5,col='red')                                          # línea indicando 1/2
  grid()

  #Calculo medias y variaciones
  cat("Los valores de las medias son: \n")
  print(prom<-apply(T1234,2, mean))        
  # valores de las medias

  cat("Los valores de las desviaciones son: \n")
  print(desv<-apply(T1234,2, sd))          
  # valores de la desviaciones estándar
}

muestrasn(datos,20)

## Los valores de las medias son: 
##        ϑ1        ϑ2        ϑ3        ϑ4 
## 0.4405283 0.8999809 0.4448185 0.5398929 
## Los valores de las desviaciones son: 
##        ϑ1        ϑ2        ϑ3        ϑ4 
## 0.1806392 0.3893014 0.1909691 0.2738864

muestrasn(datos,50)

## Los valores de las medias son: 
##        ϑ1        ϑ2        ϑ3        ϑ4 
## 0.4832036 0.9650741 0.4866900 0.6011739 
## Los valores de las desviaciones son: 
##        ϑ1        ϑ2        ϑ3        ϑ4 
## 0.2407162 0.5085414 0.2173306 0.2920160

muestrasn(datos,100)

## Los valores de las medias son: 
##        ϑ1        ϑ2        ϑ3        ϑ4 
## 0.4877034 0.9705049 0.4898138 0.6039843 
## Los valores de las desviaciones son: 
##        ϑ1        ϑ2        ϑ3        ϑ4 
## 0.2422095 0.4959850 0.2190291 0.2925472

muestrasn(datos,1000)

## Los valores de las medias son: 
##        ϑ1        ϑ2        ϑ3        ϑ4 
## 0.4999912 1.0039324 0.4999569 0.5897709 
## Los valores de las desviaciones son: 
##        ϑ1        ϑ2        ϑ3        ϑ4 
## 0.2631583 0.5558884 0.2429078 0.3230200

Conclusión 2

Analizando los resultados obtenidos en cada una de las muestras encontramos que los mejores estimadores son el 1 y el 3
- T1 y T3 Son los que más se acercan al valor esperado por lo cual son insesgados
- T1 Y T3 estaban inicialmente lejos del valor esperado, pero a medida que incrementó la muestra, empezaron a acercarse al valor esperado por lo cual son consistentes
- T1 Y T3 tienen la menor varianza por lo cual son los más eficiente
- T2 Es el que menos se acerca al valor esperado por lo cual es el más sesgado y el que mayor varianza tiene lo que lo convierte en el menos eficiente
- T4 sobrepasa el valor esperado y la menor varianza por lo cual es sesgado y no muy eficiente

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:

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%.
Genere una función que permita: Obtener una muestra aleatoria de la población Calcule el estimador de la proporciónmuestral pˆpara un tamaño de muestra dado n.
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.
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
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.

Solución 3

#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(123)
poblacion <- rbinom(1000, 1, 0.5)

#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
muestraAleatoria = function(poblacion, n) {
  if (n >= length(poblacion)) {
    return(-1)
  } 
  else {
    muestra = sample(poblacion, n, replace = FALSE)
    mean_m = mean(muestra)
    return(mean_m)
  }
}

escenarios <- replicate(500, muestraAleatoria(poblacion, 500))

mediaEscenarios <- mean(escenarios)
desviacionEscenarios <- sd(escenarios)
min <- min(escenarios)
max <- max(escenarios)
q1 <- quantile(escenarios, probs=0.25)
q3 <- quantile(escenarios, probs=0.75)
var <- var(escenarios)


contenido <- round(as.numeric(c(mediaEscenarios, desviacionEscenarios, min, max, q1, q3, var)),4)

titulos <- c("Media", "Desviacion", "Minimo", "Maximo", "Q1", "Q3", "Varianza")

tabla <- as.data.frame(rbind(titulos, contenido))

tabla

##               V1         V2     V3     V4    V5    V6       V7
## titulos    Media Desviacion Minimo Maximo    Q1    Q3 Varianza
## contenido 0.4939     0.0164  0.448   0.54 0.482 0.504    3e-04

hist(escenarios, main = "Histograma 500 iteraciones", xlab = "Probabilidad de las muestras", breaks=30, ylab = "Frecuencia", las=1,
     font.axis=4)

tamanio_muestras <- c(5, 10, 15, 20, 30, 50, 60, 100, 200, 500)
resultados <- list()

for (i in tamanio_muestras) {
  resultados[[as.character(i)]] <- replicate(500, muestraAleatoria(poblacion, i))
}

for (j in tamanio_muestras){
  cat("Muestra n = ", j, "\n")
  print(shapiro.test(resultados[[as.character(j)]]))
  qqnorm(resultados[[as.character(j)]])
  qqline(resultados[[as.character(j)]])
}

## Muestra n =  5 
## 
##  Shapiro-Wilk normality test
## 
## data:  resultados[[as.character(j)]]
## W = 0.9225, p-value = 2.281e-15

## Muestra n =  10 
## 
##  Shapiro-Wilk normality test
## 
## data:  resultados[[as.character(j)]]
## W = 0.96565, p-value = 2.017e-09

## Muestra n =  15 
## 
##  Shapiro-Wilk normality test
## 
## data:  resultados[[as.character(j)]]
## W = 0.97214, p-value = 3.749e-08

## Muestra n =  20 
## 
##  Shapiro-Wilk normality test
## 
## data:  resultados[[as.character(j)]]
## W = 0.97838, p-value = 9.212e-07

## Muestra n =  30 
## 
##  Shapiro-Wilk normality test
## 
## data:  resultados[[as.character(j)]]
## W = 0.98452, p-value = 3.642e-05

## Muestra n =  50 
## 
##  Shapiro-Wilk normality test
## 
## data:  resultados[[as.character(j)]]
## W = 0.98949, p-value = 0.001197

## Muestra n =  60 
## 
##  Shapiro-Wilk normality test
## 
## data:  resultados[[as.character(j)]]
## W = 0.99039, p-value = 0.002399

## Muestra n =  100 
## 
##  Shapiro-Wilk normality test
## 
## data:  resultados[[as.character(j)]]
## W = 0.99368, p-value = 0.03499

## Muestra n =  200 
## 
##  Shapiro-Wilk normality test
## 
## data:  resultados[[as.character(j)]]
## W = 0.99268, p-value = 0.01519

## Muestra n =  500 
## 
##  Shapiro-Wilk normality test
## 
## data:  resultados[[as.character(j)]]
## W = 0.99617, p-value = 0.2711

Conclusión 3

Los resultados obtenidos son consistentes con el valor real de la población, que en este caso es 0.5. Además, la desviación de 0.0152 indica que la mayoría de los datos se sitúan a 0.0152 unidades por encima o por debajo de la media. Aunque, al observar el histograma, se nota que no es completamente simétrico, evidenciando un coeficiente de asimetría distinto de cero. Estos hallazgos proporcionan una base sólida para demostrar el Teorema del Límite Central. A pesar de la variabilidad en el muestreo, la media de los estimadores se acerca significativamente al valor real de la población, y la desviación estándar ofrece una medida de la variabilidad presente en los resultados.
Las pruebas propuestas nos permiten contrastar si los valores anticipados concuerdan con los valores observados en el orden de la muestra. En la aplicación de la Prueba de Shapiro-Wilk, se destaca que en todos los casos el valor es cercano a 1, lo cual sugiere que los datos exhiben una distribución normal. El valor p asociado a esta prueba se utiliza para evaluar la magnitud de la diferencia obtenida, bajo la suposición de que la población sigue una distribución normal. Cuando este valor es bajo, se rechaza la hipótesis de que los datos no se ajustan a una distribución normal. En las primeras simulaciones, se observa un valor p muy pequeño, indicando que los datos no siguen una distribución normal. A medida que las simulaciones avanzan, el valor de W se aproxima a 1, pero el valor p permanece reducido, indicando persistente falta de normalidad en todas las simulaciones.

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∗1a 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:

Solución 4

set.seed(5)                                                       
#Semilla: fijar un valor inicial para el generador de números aleatorios
muestra_original = c(7.69, 4.97, 4.56, 6.49, 4.34, 6.24, 4.45)    
#Datos muestra de la eficiencia del combustible de 7 camiones
muestras = sample(muestra_original, 7000, replace=TRUE)             
#Se extraen n x m muestras
datos = matrix(muestras,nrow=1000,ncol=7)                         
#Construir matriz de n x m
#head(datos)
medias_boostrap = apply(datos,1,mean)                             
#Calcular media por fila
#medias_boostrap

#Media de la muestra original
mean(muestra_original)

## [1] 5.534286

#Calcular los percentiles 0.025 y 0.975
percentil_inf <- quantile(medias_boostrap, 0.025)
percentil_sup <- quantile(medias_boostrap, 0.975)

#METODO 1
IC1 <- c(percentil_inf, percentil_sup)
IC1

##     2.5%    97.5% 
## 4.725607 6.490536

En el Método 1: Con un 95% de confianza el promedio del rendimiento de combustible(millas/galon) se encuentra entre el intervalo [4.73, 6.49]

#METODO 2
IC2 <- c((2*mean(medias_boostrap) - percentil_sup), (2*mean(medias_boostrap) - percentil_inf))
IC2

##    97.5%     2.5% 
## 4.562576 6.327504

En el Método 2: Con un 95% de confianza el promedio del rendimiento de combustible(millas/galon) se encuentra entre el intervalo [4.56, 6.33]

# Gráfico del histograma de los datos originales
hist(muestra_original, las=1, main = "Histograma Datos Originales",
     xlab = "Eficiencia de Combustible", ylab = "Frecuencia", col="#034A94", border="white") 
# Añadido borde blanco
abline(v=mean(muestra_original), col="red", lty = 2) 
# Línea vertical en la media de la muestra
legend("topright", legend = "Media de la Muestra", col = "red", lty = 2, cex = 1, xjust = 1, yjust = 1)

# Gráfico de las medias de las muestras y los intervalos de confianza
hist(medias_boostrap, las=1, main="Histograma Medias Bootstrap", ylab = "Frecuencia ", xlab = "Promedios eficiencias ", col="#034A94", border="white") 
# Añadido borde blanco
abline(v = mean(medias_boostrap), col = "salmon", lty = 2) 
# Línea vertical en la media de la muestra
legend("topright", legend = "Media de la Muestra", col = "salmon", lty = 2, cex = 1, xjust = 1, yjust = 1)
# Intervalo de confianza 1
abline(v=IC1, col="#FF7F00", lwd=2) 

# Intervalo de confianza 2
abline(v=IC2, col="#0EB0C6", lwd=2)

Conclusión 4

Los intervalos de confianza arrojados por el método 1 y 2 contienen la verdadera media de la eficiencia del combustible de la población, la cual es 5.53, lo cual indica ambos métodos se desempeñan de manera adecuada y consistente el método 1 es un poco más amplio [4.73, 6.49] comparado con el método 2 [4.56, 6.33]

Por lo anterior si confiaría en las estimaciones de estos métodos ya que los métodos están funcionando correctamente y ambas abarcan una gran parte de la población como lo evidencia el histograma.

Métodos y simulación estadística

Juliana Arias Vergara,Jorge Barrera Salgado

2024-03-11