Problema 1

Se ha solicitado que se calcule un valor aproximado de pi (π) utilizando la función function() y una función de distribución uniforme la cual es continua con valores de mínimo y máximo de 0 y 1 respectivamente.

Creación de la función:

La función uniforme en r es runif, se dejarán n valores con mínimo 0 y máximo 1, se crearán los valores aleatorios de xi y yi

La distancia al centro corresponde (xi - 0.5)^2 + (yi - 0.5)^2

Los puntos dentro del círculo son aquellos donde la distacia al centro son menores a 0.25

Al final la función suma todos los puntos dentro del círculo cuya distancia sea menor a 0.25, este valor se multiplica por 4 y se divide entre n para obtener una aproximación a pi

# Función para estimar el valor de π
set.seed(123)
estimar_pi <- function(n) {
  # Generar las coordenadas x, y
  x <- runif(n, min = 0, max = 1)
  y <- runif(n, min = 0, max = 1)
  
  # Determinar si los puntos están dentro del círculo
  distancia_al_centro <- (x - 0.5)^2 + (y - 0.5)^2
  puntos_dentro <- distancia_al_centro < 0.25
  
  # Contar los puntos dentro del círculo y estimar π
  puntos_dentro_del_circulo <- sum(puntos_dentro)
  estimacion_pi <- 4 * puntos_dentro_del_circulo / n
  
  return(estimacion_pi)
}

Estimando pi a partir de n conjunto de puntos

Se estima el valor de π con n = 1.000, 10.000 y 100.000.

## Estimación de π con mil puntos es de: 3.2
## Estimación de π con diez mil puntos es de: 3.1528
## Estimación de π con cien mil puntos es de: 3.144

Elaborar errores

Se puede ver que entre más puntos hayan más precisa será la aproximación a pi, de igual forma todas las aproximaciones tienen un error relativo menor al 2%.

n Estimacion de pi Error Total Error Relativo(%)
1000 3.2000 0.0584073 1.8592
10000 3.1528 0.0112073 0.3567
100000 3.1440 0.0024073 0.0766

Elaborar el gráfico

A continuación, veremos representado gráficamente como se ve la simulación con 100.000 puntos, se puede apreciar que a pesar de que esta aproximación sea la más cercana al valor real el circulo se ve irregular por lo que harían falta muchos más datos para tener una aproximación realmente precisa, sin embargo es un buen acercamiento.

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.

set.seed(123)
# Crear la función exponencial
estim <- function(θ){
  x <- rexp(4, rate = 1)
  θ1 <- (x[1] + x[2] / 6) + (x[3] + x[4] / 3)
  θ2 <- (x[1] + 2 * x[2] + 3 * x[3] + 4 * x[4]) / 5
  θ3 <- (x[1] + x[2] + x[3] + x[4]) / 4
  θ4 <- (min(x[1], x[2], x[3], x[4]) + max(x[1], x[2], x[3], x[4])) / 2
  θ <- c(θ1, θ2, θ3, θ4)
  return(θ)
}

Generar muestras de n=20, 50, 100 y 1000 para cada uno de los estimadores planteados.

set.seed(123)

# Función para calcular el ECM
ecm <- function(estimaciones, valor_real) {
  mean((estimaciones - valor_real)^2)
}

# Valor real de θ para el cálculo del ECM
valor_real_θ <- 1

# Generar muestras de n=20, 50, 100 y 1000 para cada uno de los estimadores planteados.
n_values <- c(20, 50, 100, 1000)

# Inicializar los dataframes con el número correcto de filas y nombres de columnas
df_medias <- data.frame(matrix(ncol = 4, nrow = length(n_values)))
names(df_medias) <- c("θ1", "θ2", "θ3", "θ4")

df_varianzas <- data.frame(matrix(ncol = 4, nrow = length(n_values)))
names(df_varianzas) <- c("θ1", "θ2", "θ3", "θ4")

df_medianas <- data.frame(matrix(ncol = 4, nrow = length(n_values)))
names(df_medianas) <- c("θ1", "θ2", "θ3", "θ4")

df_ecm <- data.frame(matrix(ncol = 4, nrow = length(n_values)))
names(df_ecm) <- c("θ1", "θ2", "θ3", "θ4")

# Encontrar los límites globales para la escala del eje y
all_results <- matrix(unlist(lapply(n_values, function(n) replicate(n, estim(1)))), ncol = 4)
y_lim <- range(all_results)

# Dividir la ventana gráfica en múltiples paneles
par(mfrow = c(2, 2))

for (i in 1:length(n_values)) {
  n <- n_values[i]
  result <- matrix(replicate(n, estim(1)), ncol = 4)  # Aplicar estim a cada elemento de 1 a n y transponer el resultado
  
  # Generar boxplot con colores aleatorios y la misma escala en el eje y
  boxplot(result, las = 1, main = paste("n =", n), col = rainbow(4), ylim = y_lim)
  
  # Calcular y almacenar los resúmenes en dataframes
  df_medias[i, ] <- colMeans(result)
  df_varianzas[i, ] <- apply(result, 2, var)
  df_medianas[i, ] <- apply(result, 2, median)
  df_ecm[i, ] <- apply(result, 2, ecm, valor_real_θ)
}

# Asignar los nombres de las filas a los dataframes
rownames(df_medias) <- paste("n", n_values, sep = "_")
rownames(df_varianzas) <- paste("n", n_values, sep = "_")
rownames(df_medianas) <- paste("n", n_values, sep = "_")
rownames(df_ecm) <- paste("n", n_values, sep = "_")

Evaluar las propiedades de insesgadez, eficiencia y consistencia

Insesgadez

Un estimador se considera insesgado si su valor esperado es igual al parámetro que se está estimando. Comprobaremos si los estimadores son insesgados calculando su sesgo, que se define como la diferencia entre el valor esperado del estimador y el verdadero valor del parámetro.

Valor del sesgo de cada estimador con la simulación n = 20, 50, 100 y 1000.

# Mostrar el dataframe de medias
df_sesgo <- df_medias - valor_real_θ
knitr::kable(df_sesgo)%>%
  kable_styling(bootstrap_options = Form.Basic)%>%
  column_spec(1, bold = T)%>%
  row_spec(0,color = "red")
θ1 θ2 θ3 θ4
n_20 0.7695005 0.3455291 0.2576369 1.4003349
n_50 1.0505006 0.8248018 0.8827098 0.4491871
n_100 0.9030025 0.6976104 0.5751681 0.6418539
n_1000 0.6326381 0.6676144 0.6430772 0.6816256

Como se observa en la tabla, para n=20 y n=100, el estimador con menor sesgo es θ3, para n=50, θ4 refleja menor sesgo. En conclusión para tamaños de muestra pequeños, el valor del sesgo varía entre los 4 estimadores, mientras que para valores grandes como n=1000 el valor del estimador tiende a estabilizarse en 0.6. Así las cosas, θ3 tiene menor insesgadez.

Eficiencia

La eficiencia de un estimador se refiere a la precisión con la que estima el parámetro desconocido. Se presenta a continuación la varianza de cada estimador para evaluar su eficiencia.

Varianza de cada estimador con la simulación n = 20, 50, 100 y 1000.

# Mostrar el dataframe de varianzas
knitr::kable(df_varianzas)%>%
  kable_styling(bootstrap_options = Form.Basic)%>%
  column_spec(1, bold = T)%>%
  row_spec(0,color = "red")
θ1 θ2 θ3 θ4
n_20 1.452332 0.5218423 0.7564330 2.0375534
n_50 2.098443 2.5370587 1.6479542 0.5623773
n_100 1.532540 1.8122101 0.9601632 1.3022048
n_1000 1.316660 1.3776755 1.4272373 1.1635706

Observando la tabla, para n=20, θ2 tiene la menor varianza, mientras que θ4 tiene menor varianza para n=50 y n=1000, por último, para n=100 la menor varianza corresponde a θ3. En ese orden de ideas, θ4 es más eficiente.

Consistencia

Un estimador se considera consistente si converge al verdadero valor del parámetro a medida que el tamaño de la muestra aumenta.

Media de cada estimador con la simulación n = 20, 50, 100 y 1000.

knitr::kable(df_medias)%>%
  kable_styling(bootstrap_options = Form.Basic)%>%
  column_spec(1, bold = T)%>%
  row_spec(0,color = "red")
θ1 θ2 θ3 θ4
n_20 1.769500 1.345529 1.257637 2.400335
n_50 2.050501 1.824802 1.882710 1.449187
n_100 1.903002 1.697610 1.575168 1.641854
n_1000 1.632638 1.667614 1.643077 1.681626

Como se aprecia en la tabla, θ1 y θ2 con consistentes, es decir, a partir de n = 50 cuando aumenta el tamaño de n, el valor esperado de los parámetros se acercan al valor del parámetro (θ=1). En contraste, θ3 y θ4 presentan valores erráticos. θ4 tiene el valor mas cercano al parámetro (θ=1) cuando n tiene el valor más alto.

Así las cosas, no hay un único estimador que sea consistentemente el mejor en todas las simulaciones.

Problema 3

a. Generar la población de tamaño n=1000 con un porcentaje de individuos enfermos del 50%:

pob <- c(rep(1, 500), rep(0, 500))

b. Generar una función para obtener una muestra aleatoria de la población y calcular el estimador de la proporción muestral p-hat para un tamaño de muestra dado n:

set.seed(123)
muestra_proporcion <- function(poblacion, n) {
  muestra <- sample(poblacion, size = n, replace = FALSE)
  proporcion <- sum(muestra) / n
  return(proporcion)
}

c. Repetir el escenario anterior (b) n=500 veces y analizar los resultados en cuanto al comportamiento de los 500 resultados del estimador p-hat:

Tomamos una muestra de tamaño n = 100.

set.seed(123)
resultados <- data.frame(x = numeric())
for (i in 1:500){
  z <- data.frame(x = muestra_proporcion(pob, 100))
  resultados <- rbind(resultados, z)
}
# Comprobar simetría y sesgo de los resultados
media <- mean(resultados$x)
mediana <- median(resultados$x)
# Calcular la variabilidad de los resultados
desv <- sd(resultados$x)

Si realizamos la simulación vemos que el promedio nos da 0.49998 lo cual es muy cercano a la distribución 50-50 original de la población, la mediana nos da exactamente 0.5 lo que es poco mayor que la media por lo que podemos decir que la distribución es bastante simétrica y finalmente la desviación estándar es de 0.0475807 por lo que los datos no se alejaron mucho de la media en general.

d. Repetir los puntos b y c para diferentes tamaños de muestra y evaluar la normalidad de los resultados:

# Establecer la semilla aleatoria
set.seed(123)

# Crear una lista para almacenar los resultados de Shapiro
shapiro <- list()

# Valores de n
ns <- c(5, 10, 15, 20, 30, 50, 60, 100, 200, 500)

# Aumentar el tamaño de la ventana gráfica
options(repr.plot.width = 12, repr.plot.height = 10)

# Dividir el espacio de la figura en una matriz de 5x2
par(mfrow = c(5, 2))

# Iterar sobre los valores de n
for (i in 1:10) {
  n <- ns[i]
  
  # Crear un dataframe temporal para almacenar las muestras
  temp <- data.frame(x = numeric())
  
  # Generar muestras y realizar el test de Shapiro-Wilk
  for (j in 1:500) {
    z <- data.frame(x = muestra_proporcion(pob, n))
    temp <- rbind(temp, z)
  }
  
  # Almacenar el p-valor del test de Shapiro-Wilk
  shapiro[[i]] <- shapiro.test(temp$x)$p
  
  # Ajustar los márgenes
  par(mar = c(2, 2, 1, 1))
  
  # Graficar histograma
  hist(temp$x, main = paste0("Histograma n=", n), xlab = "Valor", ylab = "Frecuencia", col = "skyblue")
  
  # Graficar Q-Q normal
  par(mar = c(2, 2, 1, 1))
  qqnorm(temp$x, main = paste0("Q-Q Normal n=", n), cex = 0.8)
  qqline(temp$x, col = "red")
}

 p_valor<- as.data.frame(shapiro)
  colnames(p_valor)<-c("n=5", "n=10", "n=15", "n=20", "n=30", "n=50", "n=60", "n=100", "n=200", "n=500" )
  p_valor <- round(p_valor, digits = 4)
  
  kable(p_valor)%>%
  kable_styling(bootstrap_options = Form.Basic)%>%
  column_spec(1, bold = T)%>%
  row_spec(0,color = "red")
n=5 n=10 n=15 n=20 n=30 n=50 n=60 n=100 n=200 n=500
0 0 0 0 0 0.0013 0.0043 0.0069 0.1709 0.0368

Si hacemos la prueba de Shapiro-Wilk con una significancia del 0.05 casi ninguna muestra se podría categorizar como una distribución normal, solo la que tiene n=200 que es la única con un p-valor mayor al 0.05 el cual es 0.17. Las muestras n=500 se acercan teniendo un p-valor de 0.036 sin embargo tenemos suficiente información para rechazar que esa distribución es normal, el resto de valores están muy por debajo del 0.05 por lo que ninguno tiene una distribución normal.

e. Repetir toda la simulación para porcentajes de plantas enfermas del 10% y 90%:

10% enfermas

a10. Generar la población de tamaño n=1000 con un porcentaje de individuos enfermos del 10%:
pob <- c(rep(1, 100), rep(0, 900))

c10. Repetir el escenario anterior (b) n=500 veces y analizar los resultados en cuanto al comportamiento de los 500 resultados del estimador p-hat:

Tomamos una muestra de tamaño n = 100.

Si realizamos la simulación vemos que el promedio nos da 0.10146 lo cual es muy cercano a la distribución 10-90 original de la población, la mediana nos da exactamente 0.1 lo que es poco menor que la media por lo que podemos decir que la distribución es bastante simétrica y finalmente la desviación estándar es de 0.0301045 por lo que los datos no se alejaron mucho de la media en general.

d10. Repetir los puntos b y c para diferentes tamaños de muestra y evaluar la normalidad de los resultados:

n=5 n=10 n=15 n=20 n=30 n=50 n=60 n=100 n=200 n=500
0 0 0 0 0 0 0 0 0.0118 0.0883

En este caso solo el n=500 tuvo un valor p superior a 0.05 siendo 0.088, el resto de distribuciones entonces no son normales.

90% enfermas

a90. Generar la población de tamaño n=1000 con un porcentaje de individuos enfermos del 90%:
pob <- c(rep(1, 900), rep(0, 100))
c90. Repetir el escenario anterior (b) n=500 veces y analizar los resultados en cuanto al comportamiento de los 500 resultados del estimador p-hat:

Tomamos una muestra de tamaño n = 100.

Si realizamos la simulación vemos que el promedio nos da 0.90008 lo cual es muy cercano a la distribución 90-10 original de la población, la mediana nos da exactamente 0.9 lo que es poco menor que la media por lo que podemos decir que la distribución es bastante simétrica y finalmente la desviación estándar es de 0.028984 por lo que los datos no se alejaron mucho de la media en general.

d90. Repetir los puntos b y c para diferentes tamaños de muestra y evaluar la normalidad de los resultados:

n=5 n=10 n=15 n=20 n=30 n=50 n=60 n=100 n=200 n=500
0 0 0 0 0 0 0 1e-04 0.002 0.0062

En esta ocasión ninguna distribución es normal ya que todos los valores están por debajo del valor de significancia 0.05, el valor p que más se aproxima es el de n=500 que es de 0.006.

Problema 4: Estimación boostrap

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.

Método 1 percentiles P2.5 y P97.5

Se calculan los percentiles P2.5 y P97.5

#MÉTODO 1

intervalo_metodo1 <- c(P_2.5, P_97.5)
intervalo_metodo1
##    2.5%   97.5% 
## 4.82895 6.27805

Método 2: 2X¯−P97.5;2X¯−P2.5

Uso de la media y los percentiles

#MÉTODO 2

intervalor_metodo2 <- c(2*X_bar - P_97.5, 2*X_bar -P_2.5)
intervalor_metodo2
##    97.5%     2.5% 
## 4.783085 6.232185

Histograma

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

En la gráfica, se muestran los intervalos en color coral para el método 1 y en verde para el método 2. ## Conclusiones:

Ambos métodos proporcionan intervalos de confianza con un ancho similar. La diferencia entre los límites de los intervalos es pequeña (0.006). En este caso, parece que ambos métodos son confiables para estimar el intervalo de confianza de la media de la eficiencia de combustible.

El tamaño de la muestra (n) y el nivel de confianza (nivel_confianza) pueden afectar el ancho del intervalo de confianza. La distribución de los datos también puede afectar la confiabilidad de los intervalos de confianza.