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 π?

Solución:

Ejercicio planteado con 1.000 puntos, se tiene que:

## Con 1.000 puntos
## Estimación del valor de π:  3.136
## Estimación número de puntos dentro del círculo:  784

Ahora se procede a realizar el mismo ejercicio pero con 10.000 puntos:

## Con 10.000 puntos
## Estimación del valor de π:  3.1012
## Estimación número de puntos dentro del círculo:  7753

Posteriormente se procede a realizar el ejercicio con 100.000 puntos

## Con 100.000 puntos
## Estimación del valor de π:  3.1396
## Estimación número de puntos dentro del círculo:  78490

Bueno y solo por probar con 1.000.000 de puntos

## Con 1.000.000 puntos
## Estimación del valor de π:  3.140908
## Estimación número de puntos dentro del círculo:  785227

A continuación dejo anexo el código usado para la resolución del problema 1:

# Este sería el valor inicial planteado de puntos a generar:
n <- 1000 

# Ahora se procede a la generación de coordenadas en X - Y, con una distribución uniforme con un valor mínimo de 0 y valor máximo de 1.

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

# Se inicializa un contador para poder llevar el conteo de los puntos dentro del círculo.

puntos_dentro <- 0

# Ciclo para calcular la distancia al centro y contar puntos dentro del círculo

for (i in 1:n) {
  distancia_al_centro <- (x[i] - 0.5)^2 + (y[i] - 0.5)^2
  if (distancia_al_centro <= 0.25) {
    puntos_dentro <- puntos_dentro + 1
  }
}

# Formula para calcular la estimación de π

estimacion_pi <- 4 * (puntos_dentro / n)

# Mostrando el resultado de estimación de π

cat("Con 1000 puntos", "\n")

cat("Estimación del valor de π: ", estimacion_pi, "\n")

# Mostrando el resultado de los puntos dentro del círculo

cat("Estimación número de puntos dentro del círculo: ", puntos_dentro, "\n")

Conclusiones:

  • A mayor número de puntos tomados para el ejercicio, el valor de π es mucho más preciso, es decir a mayor muestra mayor exactitud.

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:

θ1ˆ=(X1+X2)/6+(X3+X4)/3

θ2ˆ=(X1+2X2+3X3+4X4)/5

θ3ˆ=((X1+X2+X3+X4)/4)$

θ4ˆ=(min{X1,X2,X3,X4}+max{X1,X2,X3,X4})/2

Solución:

1.- Para el estimador θ1ˆ=(X1+X2)/6+(X3+X4)/3

## Tamaño de muestra: 20 Sesgo: 0.449199 Varianza: 0.2326798 
## Tamaño de muestra: 50 Sesgo: 0.7910028 Varianza: 0.343955 
## Tamaño de muestra: 100 Sesgo: 1.129824 Varianza: 0.3296568 
## Tamaño de muestra: 1000 Sesgo: 1.038689 Varianza: 0.2738798

2.- Para el estimador θ2ˆ=(X1+2X2+3X3+4X4)/5

## Tamaño de muestra: 20 Sesgo: 0.6664907 Varianza: 2.131285 
## Tamaño de muestra: 50 Sesgo: 1.296123 Varianza: 1.022809 
## Tamaño de muestra: 100 Sesgo: 0.6624348 Varianza: 1.009784 
## Tamaño de muestra: 1000 Sesgo: 1.266311 Varianza: 1.318818

3.- Para el estimador θ3ˆ=((X1+X2+X3+X4)/4)$

## Tamaño de muestra: 20 Sesgo: -0.471273 Varianza: 0.2009733 
## Tamaño de muestra: 50 Sesgo: 0.3360119 Varianza: 0.1507182 
## Tamaño de muestra: 100 Sesgo: 0.4925072 Varianza: 0.2357519 
## Tamaño de muestra: 1000 Sesgo: 0.8369656 Varianza: 0.2595523

4.- Para el estimador θ4ˆ=(min{X1,X2,X3,X4}+max{X1,X2,X3,X4})/2

## Tamaño de muestra: 20 Sesgo: 0.9464553 Varianza: 0.6214501 
## Tamaño de muestra: 50 Sesgo: 1.480356 Varianza: 0.1960333 
## Tamaño de muestra: 100 Sesgo: 1.343953 Varianza: 0.4213955 
## Tamaño de muestra: 1000 Sesgo: 1.388852 Varianza: 0.4469479

Anexo el código usado para el problema en mención:

## install.packages("moments")
library(moments)

## Definiendo el tamaño de las muestras 

nume_muestra <- c(20, 50, 100, 1000)

## Definiendo la función 1:

funcion1 <- function(x){(x[1] + x[2])/6 + (x[3] + x[4])/3}
## funcion1 <- function(x){(x[1] + 2 * x[2] + 3 * x[3] + 4 * x[4])/5}
## funcion1 <- function(x){(x[1] + x[2] + x[3] +  x[4])/4}
## funcion1 <- function(x)(min(x) + max(x))/2

## Definiendo un arreglo, tipo lista para guardar los resultados de cada tamaño estipulado en n = 20, 50, 100, 1000

resul <- list()

## Ciclo para iterar sobre los diferentes tamaños de n, trayendo o arrojando como resultado los valores de sesgo y varianza para cada tamaño de muestra.

for (n in nume_muestra) {
  muestras <- matrix(rexp(4 * n), ncol = 4)
  estimaciones <- apply(muestras, 1, funcion1) 
  resul[[as.character(n)]] <- estimaciones
  sesgo <- skewness(estimaciones)
  varianza <- var(estimaciones)

## Impresión del tamaño de la muestra:
  
  cat("Tamaño de muestra:", n, "")
  
## Impresión del valor del sesgo para cada tamaño de muestra.
  
  cat("Sesgo:", sesgo, "") 
  
## Impresión del valor de la varianza para cada tamaño de muestra.  
  
  cat("Varianza:", varianza, "\n")
}

## Grafica de cajas para cada valor de muestra:

boxplot(resul, names = nume_muestra, main = "Distribución de Estimaciones θ1,2,3,4", col = rainbow(ncol(trees)))

Conclusiones:

Se encuentra que entre mayor número de datos de la muestra, tiende a la insesgadez.

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:

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

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

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

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

Solución:

  1. Simulación de población n=1000
set.seed(123)
n_lote <- 1000
plant_enfermas <- 0.5
poblacion <- rbinom(n_lote, 1, plant_enfermas)
  1. Generación de función
estimador <- function(x){
  muestra <- sample(poblacion, x, replace = FALSE)
  propuesta <- mean(muestra)
  return(propuesta)
}
  1. Simulación de n = 500
simulacion <- 500
x <- c(5, 10, 15, 20, 30, 50, 60, 100, 200, 500)

resul <- matrix(nrow = simulacion, ncol = length(x))

for (i in 1:simulacion){
  for (j in 1:length(x)) {
    estimador2 <- (estimador(x[j]))
    resul[i,j] <- (estimador2)
    
  }
}

hist(resul, main = "Histograma prueba", xlab = "Media de muestra", ylab = "Frecuencia")

Resultados: Se aprecia que a medidad que el tamañ de la muestra aumenta, la distribución de los estimadores tiende a normalizarse, es decir tiende a una distribución normal.

  1. Normalidad
library(stats)

analizar_normalidad <- function(x, titulo) {
  shapiro_result <- shapiro.test(x)
  qqnorm(x, main = titulo)
  qqline(x) 

  
  cat("Prueba de Shapiro-Wilk:")
  print(shapiro_result)
}

## Realizando el análisis de normalidad con las diferentes muestras 

for (j in 1:length(x)) { 
  analizar_normalidad(resul[, j], paste("n =", x[j]))
  }

Resultados: A medida que la muestra de enfermos aumenta, la distribución de estimadores tiende a desviarse de la normalidad.

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

  • Método 1 (P2.5;P97.5)
  • Método 2 (2X¯−P97.5;2X¯−P2.5)

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

Solución:

## Intervalo de confianza por el método 1): 4.721393 6.4215
## Intervalo de confianza por el método 2): 4.647071 6.347179

Resultados: Al comparar los resultados viendo la similitud de los intervalos de confianza, uno podría confiar en estas estimaciones.

Cabe resaltar que los datos al no provenir de una distribución normal, no serían precisos respecto a sus intervalos de confianza.

Se anexa código utilizado para la resolución del problema 4.

## Definiendo la muestra con los datos aportados 

data <- c(7.69, 4.97, 4.56, 6.49, 4.34, 6.24, 4.45)

n <- length(data)
k <- 1000
m <- 0.05

## usando boostrap

muestra <- replicate(k,sample(data, replace = TRUE))

## Calcula las medias de la muestra (boostrap)

media <- apply(muestra, 2, mean)

## Aplicando el método 1, para determinar los intervalos de confianza 

metodo1 <- quantile(media, c(m/2, 1 - m/2))

media2 <- mean(data)

## Aplicando el método 2, para determinar los intervalos de confianza 

metodo2 <- c(2*media2 - quantile (media, 1 - m/2),2 * media2 - quantile(media, m/2))


cat("Intervalo de confianza por el método 1):", metodo1, "")
cat("Intervalo de confianza por el método 2):", metodo2, "")