Actividad 2: Métodos y simulación estadística

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:

1. Estimación del valor de π

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

a. Genere n coordenadas:

Se genera la función create_sample la cual recibe como parámetro n y genera una distribución uniforme de coordenadas usando runif:

# Function to create a sample of coordinates
create_sample <- function(n) {
  x <- runif(n)
  y <- runif(n)
  return(data.frame(x, y))
}

La función recibe el número de puntos a tomar, y devuelve un dataframe con dos columnas (x e y), las cuales contienen los puntos generados con la distribución uniforme.

b. Genere 1.000 coordenadas

Se usa la función create_sample pasando como parámetro la cantidad de coordenadas deseadas, en este caso n = 1.000. Además, se genera un seed para que los resultados pierdan aleatoriedad.

# Generate a sample of coordinates
set.seed(1) # Set seed for reproducibility
n = 1000
sample <- create_sample(n)

c. calcular distancia del centro del círculo

Se usa la función center_distance para realizar el cálculo por pitágoras para hallar la distancia del punto hasta el centro del círculo de diámetro 1. La función se aplica al dataframe generado en el paso anterior, el resultado es un vector de distancias, 1.000 distancias dado que es el número de n asignado.

# Function to calculate the distance from each point to the center
center_distance <- function(sample) {
  return((sample$x - 0.5)^2 + (sample$y - 0.5)^2)
}

d. ¿Cuántos de los puntos están dentro del círculo? ¿Cuál es su estimación de π?

# Calculate distances to the center for each point
distances <- center_distance(sample)

points_in <- sum(distances < 0.25)
cat("cantidad de puntos dentro del círculo: ", points_in)
## cantidad de puntos dentro del círculo:  770
perc_points_in <- points_in/n
cat("Estimación de π para muestra con 1.000 puntos: ", perc_points_in)
## Estimación de π para muestra con 1.000 puntos:  0.77
perc_points_in <- points_in/n
cat("Valor real π: ", pi/4)
## Valor real π:  0.7853982

e. Estimación con 10.000 puntos

# Generate sample
n = 10000
sample <- create_sample(n)

# Calculate distances to the center for each point
distances <- center_distance(sample)

points_in <- sum(distances < 0.25)
cat("cantidad de puntos dentro del círculo: ", points_in)
## cantidad de puntos dentro del círculo:  7854
perc_points_in <- points_in/n
cat("Estimación de π para muestra con 10.000 puntos: ", perc_points_in)
## Estimación de π para muestra con 10.000 puntos:  0.7854
perc_points_in <- points_in/n
cat("Valor real π: ", pi/4)
## Valor real π:  0.7853982

En conclusión, se encuentra que el valor estimado (con muestra de 1.000: 0.77 y con muestra de 10.000: 0.7854) tiende a ser muy cercano al valor real de 0.7853. Mientras mayor es la muestra, más se acerca el valor estimado al valor real.

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.

2.1 Función de cálculo estimadores

El primer paso es generar una función para aplicarla sobre las muestras de la distribución con las características solicitadas.

set.seed(4)

calc_estimators <- function(x1, x2, x3, x4) {
  t1 <- (x1 + x2)/6 + (x3 + x4)/3
  t2 <- (x1 + 2*x2 + 3*x3 + 4*x4)/5
  t3 <- (x1 + x2 + x3 + x4)/4
  t4 <- (min(x1, x2, x3, x4) + max(x1, x2, x3, x4))/2

  return(c(t1, t2, t3, t4))
}

2.2 Aplicación de la función y generación de datos

Lo siguiente es realizar la aplicación de la función en un dataframe de las características solicitadas, en este caso será una con un vector de 4.000 datos de la distribución exponencial distribuidos en cuatro columnas. Luego de la aplicación de la función df contendrá la información del cálculo de los indicadores.

rate <- 1 ## Rate = 1/lambda => lambda = 1
random_rexp <- matrix(rexp(4000, rate = rate), ncol = 4, nrow = 1000)
df_rexp <- data.frame(random_rexp)
series <- apply(df_rexp, 1, function(col) calc_estimators(col[1], col[2], col[3], col[4]))
df <- data.frame(t(series))
colnames(df) <- c("t1", "t2", "t3", "t4")
head(df, 3)
##          t1       t2        t3        t4
## 1 0.5910686 1.359468 0.4732722 0.7988039
## 2 1.2439094 1.906450 1.5126452 2.3130521
## 3 0.7849282 1.607025 0.9413257 1.1906450

2.3 Revisión de los estimadores generados y sus propiedades

Se realiza un proceso iterativo para realizar revisión del comportamiento de los estimadores conforme aumenta la cantidad de estimadores tomados de la muestra.

sample_n = c(20, 50, 100, 1000)

complete_stats = data.frame()

for (n in sample_n) {
  cat(sprintf("Boxplot de estimadores para %s observaciones", n))
  df_sample <- head(df, n)
  boxplot(df_sample, main = sprintf("%s observaciones", n))
  abline(h=1,  col="red")
  stats = sapply(df_sample, function(x) c(mean=mean(x), var=var(x), sd=sd(x)))
  colnames(stats) <- c(sprintf("t1_%s", n), sprintf("t2_%s", n), sprintf("t3_%s", n), sprintf("t4_%s", n)) 
  complete_stats = rbind(complete_stats, t(stats))
}
## Boxplot de estimadores para 20 observaciones

## Boxplot de estimadores para 50 observaciones

## Boxplot de estimadores para 100 observaciones

## Boxplot de estimadores para 1000 observaciones

También se genera una tabla que contiene los datos de la media, la varianza y la desviación estándar de los indicadores, la cual tiene como nombre complete_stats. La tabla se ve de la siguiente manera:

complete_stats[rownames(complete_stats)[order(rownames(complete_stats))], ]
##              mean       var        sd
## t1_100  0.9724227 0.2788806 0.5280915
## t1_1000 1.0082160 0.2781248 0.5273754
## t1_20   0.8084008 0.2585561 0.5084841
## t1_50   0.9119348 0.3122715 0.5588126
## t2_100  1.9710698 1.2589164 1.1220144
## t2_1000 2.0115496 1.1766746 1.0847463
## t2_20   1.6195049 1.0121450 1.0060542
## t2_50   1.8412903 1.3891100 1.1786051
## t3_100  0.9748315 0.2292699 0.4788214
## t3_1000 1.0129165 0.2500660 0.5000660
## t3_20   0.8238081 0.2317945 0.4814505
## t3_50   0.9097480 0.2559280 0.5058932
## t4_100  1.1151086 0.3372509 0.5807331
## t4_1000 1.1880112 0.4126039 0.6423425
## t4_20   0.9471879 0.2869552 0.5356820
## t4_50   1.0829970 0.4187947 0.6471435

2.4. Conclusiones

Realizando la revisión de la tabla resultado y las gráficas generadas se encuentra que:

  • Estimador 1: En promedio t1 se acerca al valor lambda definido como 1, por lo tanto se podría definir que el estimador es insesgado. También presenta una varianza pequeña, lo que permite clasificarlo como eficiente. Por último, a medida que fue subiendo la cantidad de la muestra, su promedio empezó a acercarse a lamda, por lo tanto se puede decir que también es un estimador consistente.
  • Estimador 2: Este estimador en promedio se encuentra en un valor lambda de 2, lo que permite decir que es sesgado. Tiene una varianza grande con respecto al resto de estimadores, por lo tanto es ineficiente y a medida que crece la muestra se aleja más del valor esperado de lambda, por lo tanto es inconsistente
  • Estimador 3: Este estimador en promedio se encuentra orbitando el valor 1, correspondiente a lambda, por lo tanto es es insesgado. Tiene una varianza en promedio más pequeña que la de t1, lo que permite decir que es eficiente. Y mientras la cantidad de la muestra sube, su promedio empezó a acercarse a lamda, por lo tanto se puede decir que también es un estimador consistente
  • Estimador 4: Este estimador en promedio se encuentra orbitando el valor 1, correspondiente a lambda, por lo tanto es es insesgado. Tiene una varianza similara la de t3 en promedio, lo que permite decir que es eficiente. Sin embargo, mientras aumenta la muestra se va a alejando del valor esperado de lambda, por lo tanto es inconsistente

3. Teorema del Límite Central

3.1. Generar una población de n=1000

Población de 1.000 plantas, 500 enfermas representadas por el 1 en la secuencia.

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

3.2. Genere una función para muestreo y cálculo de p

# Función para muestreo y cálculo de p_estimados
get_p_estimados <- function(population, n) {
  p_estimados = c()
  for (i in seq(500)) {
    p_sample <-  sample(population, n)
    p_estimados <- append(p_estimados, (sum(p_sample)/n))  
  }
  return (p_estimados)
}

3.3. Analice los resultados en cuanto al comportamiento de los 500 resultados

p_estimados <- get_p_estimados(population, 500)
boxplot(p_estimados)
abline(h=0.5,  col="blue")

plotn(p_estimados)

qqnorm(p_estimados, main = sprintf("qqnorm para muestra de tamaño %s", 500))

cat("Promedio:", mean(p_estimados), "Varianza:", var(p_estimados))
## Promedio: 0.501136 Varianza: 0.0002696648

  • ¿Qué tan simétricos o sesgados son los resultados obtenidos?_: Los resultados tienden a ser insesgados dado que la media de los promedios encontrados se acerca mucho al valor poblacional definido de 0.5. El valor muestral es de 0.501136

  • ¿Qué se puede observar en cuanto a la variabilidad?_: Se tiene una varianza relativamente pequeña, lo que indica que la distribución va a tender a orbitar a la media definida de 0.5.

3.4. Pruebas de bondad y ajuste

size_sample = c(5, 10, 15, 20, 30, 50, 60, 100, 200, 500)
for (m in size_sample) {
  p_estimados <- get_p_estimados(population, m)
  cat("\nShapiro test para muestra de tamaño", m)
  print(shapiro.test(p_estimados))
}
## 
## Shapiro test para muestra de tamaño 5
##  Shapiro-Wilk normality test
## 
## data:  p_estimados
## W = 0.92986, p-value = 1.498e-14
## 
## 
## Shapiro test para muestra de tamaño 10
##  Shapiro-Wilk normality test
## 
## data:  p_estimados
## W = 0.96177, p-value = 4.139e-10
## 
## 
## Shapiro test para muestra de tamaño 15
##  Shapiro-Wilk normality test
## 
## data:  p_estimados
## W = 0.97766, p-value = 6.23e-07
## 
## 
## Shapiro test para muestra de tamaño 20
##  Shapiro-Wilk normality test
## 
## data:  p_estimados
## W = 0.98045, p-value = 2.986e-06
## 
## 
## Shapiro test para muestra de tamaño 30
##  Shapiro-Wilk normality test
## 
## data:  p_estimados
## W = 0.98277, p-value = 1.2e-05
## 
## 
## Shapiro test para muestra de tamaño 50
##  Shapiro-Wilk normality test
## 
## data:  p_estimados
## W = 0.98929, p-value = 0.001029
## 
## 
## Shapiro test para muestra de tamaño 60
##  Shapiro-Wilk normality test
## 
## data:  p_estimados
## W = 0.98776, p-value = 0.0003343
## 
## 
## Shapiro test para muestra de tamaño 100
##  Shapiro-Wilk normality test
## 
## data:  p_estimados
## W = 0.99454, p-value = 0.07206
## 
## 
## Shapiro test para muestra de tamaño 200
##  Shapiro-Wilk normality test
## 
## data:  p_estimados
## W = 0.99579, p-value = 0.2022
## 
## 
## Shapiro test para muestra de tamaño 500
##  Shapiro-Wilk normality test
## 
## data:  p_estimados
## W = 0.99529, p-value = 0.1347

En cuanto a la normalidad, realizando la revisión de los gráficos generados y los test de shapiro aplicados, se puede evidenciar, en primera medida, un aumento paulatino del valor de p_value al ir aumentando la cantidad de muestra procesada. Esto puede indicar que mientras más datos son tomados, más se va acercando la muestra a tener una distribución normal. El p-value mayor a 0.05 indica una mayor probabilidad de tener una distribución normal.

A la misma conclusión se puede llegar realizando la revisión de los gráficos qqnorm. Mientras más aumenta el valor de n tomado para realizar la muestra, mayor va siendo la correspondencia de los datos con los cuantiles. Esto permite, al final, develar la distribución normal de la población.

3.5. Repetición de la simulación:

3.5.1. Población de 10 % plantas enfermas:

population_10 <- c(rep(1, 100), rep(0, 900))

p_estimados <- get_p_estimados(population_10, 500)
boxplot(p_estimados)
abline(h=0.1,  col="blue")

cat("Promedio:", mean(p_estimados), "Varianza:", var(p_estimados))
## Promedio: 0.099804 Varianza: 8.726411e-05
plotn(p_estimados)

qqnorm(p_estimados, main = sprintf("qqnorm para muestra de tamaño %s", 500))

size_sample = c(5, 10, 15, 20, 30, 50, 60, 100, 200, 500)
for (m in size_sample) {
  p_estimados <- get_p_estimados(population_10, m)
  cat("\nShapiro test para muestra de tamaño", m)
  print(shapiro.test(p_estimados))
}
## 
## Shapiro test para muestra de tamaño 5
##  Shapiro-Wilk normality test
## 
## data:  p_estimados
## W = 0.70902, p-value < 2.2e-16
## 
## 
## Shapiro test para muestra de tamaño 10
##  Shapiro-Wilk normality test
## 
## data:  p_estimados
## W = 0.84755, p-value < 2.2e-16
## 
## 
## Shapiro test para muestra de tamaño 15
##  Shapiro-Wilk normality test
## 
## data:  p_estimados
## W = 0.90031, p-value < 2.2e-16
## 
## 
## Shapiro test para muestra de tamaño 20
##  Shapiro-Wilk normality test
## 
## data:  p_estimados
## W = 0.92119, p-value = 1.652e-15
## 
## 
## Shapiro test para muestra de tamaño 30
##  Shapiro-Wilk normality test
## 
## data:  p_estimados
## W = 0.94929, p-value = 4.594e-12
## 
## 
## Shapiro test para muestra de tamaño 50
##  Shapiro-Wilk normality test
## 
## data:  p_estimados
## W = 0.97061, p-value = 1.818e-08
## 
## 
## Shapiro test para muestra de tamaño 60
##  Shapiro-Wilk normality test
## 
## data:  p_estimados
## W = 0.97806, p-value = 7.73e-07
## 
## 
## Shapiro test para muestra de tamaño 100
##  Shapiro-Wilk normality test
## 
## data:  p_estimados
## W = 0.9823, p-value = 8.966e-06
## 
## 
## Shapiro test para muestra de tamaño 200
##  Shapiro-Wilk normality test
## 
## data:  p_estimados
## W = 0.99045, p-value = 0.00252
## 
## 
## Shapiro test para muestra de tamaño 500
##  Shapiro-Wilk normality test
## 
## data:  p_estimados
## W = 0.99289, p-value = 0.01803

3.5.2. Población de 90 % plantas enfermas:

population_90 <- c(rep(1, 900), rep(0, 100))

p_estimados <- get_p_estimados(population_90, 500)
boxplot(p_estimados)
abline(h=0.9,  col="blue")

cat("Promedio:", mean(p_estimados), "Varianza:", var(p_estimados))
## Promedio: 0.899592 Varianza: 8.624603e-05
plotn(p_estimados)

qqnorm(p_estimados, main = sprintf("qqnorm para muestra de tamaño %s", 500))

size_sample = c(5, 10, 15, 20, 30, 50, 60, 100, 200, 500)
for (m in size_sample) {
  p_estimados <- get_p_estimados(population_90, m)
  cat("\nShapiro test para muestra de tamaño", m)
  print(shapiro.test(p_estimados))
}
## 
## Shapiro test para muestra de tamaño 5
##  Shapiro-Wilk normality test
## 
## data:  p_estimados
## W = 0.69951, p-value < 2.2e-16
## 
## 
## Shapiro test para muestra de tamaño 10
##  Shapiro-Wilk normality test
## 
## data:  p_estimados
## W = 0.85444, p-value < 2.2e-16
## 
## 
## Shapiro test para muestra de tamaño 15
##  Shapiro-Wilk normality test
## 
## data:  p_estimados
## W = 0.90255, p-value < 2.2e-16
## 
## 
## Shapiro test para muestra de tamaño 20
##  Shapiro-Wilk normality test
## 
## data:  p_estimados
## W = 0.92294, p-value = 2.544e-15
## 
## 
## Shapiro test para muestra de tamaño 30
##  Shapiro-Wilk normality test
## 
## data:  p_estimados
## W = 0.94829, p-value = 3.305e-12
## 
## 
## Shapiro test para muestra de tamaño 50
##  Shapiro-Wilk normality test
## 
## data:  p_estimados
## W = 0.96989, p-value = 1.306e-08
## 
## 
## Shapiro test para muestra de tamaño 60
##  Shapiro-Wilk normality test
## 
## data:  p_estimados
## W = 0.97811, p-value = 7.939e-07
## 
## 
## Shapiro test para muestra de tamaño 100
##  Shapiro-Wilk normality test
## 
## data:  p_estimados
## W = 0.97968, p-value = 1.91e-06
## 
## 
## Shapiro test para muestra de tamaño 200
##  Shapiro-Wilk normality test
## 
## data:  p_estimados
## W = 0.98789, p-value = 0.0003663
## 
## 
## Shapiro test para muestra de tamaño 500
##  Shapiro-Wilk normality test
## 
## data:  p_estimados
## W = 0.99172, p-value = 0.006878

3.5. Conclusiones

Realizando la revisión de los shapiro test generados, se encuentra que para este tipo de distribuciones que poseen colas hacia la izquierda o hacia la derecha no pasan métodos de normalidad. Los p-values del test nunca superan el umbral de 0.05, lo que quiere decir que estadísticamente se rechaza la hipótesis nula, lo que permite corroborar que los datos no siguen una distribución normal. Sin embargo, dentro de los qq norm se encuentra una tendencia hacia la normalidad, la cual no se confirma dentro de los test estadísticos para las variaciones de poblaciones no balanceadas.

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.Construya el intervalo de confianza por los dos métodos y compare los resultados obtenidos.

4.1. Generación de muestras y promedios

Se realiza cálculo de los promedios muestrales, suponiendo muestras de 4 para cada uno de las iteraciones. Se realizan 1.000 iteraciones con el vector de promedios que el ejercicio indica. Así mismo, se calcula el promedio de los promedios y se almacena en X_mean, además de los percentiles 2.5 y 97.5

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

means_sample = c()
for (i in 1:1000) {
  means_sample <- append(means_sample, mean(sample(x, 4)))
}

X_mean <- mean(means_sample)
q025 <- quantile(means_sample, 0.025) 
q975 <- quantile(means_sample, 0.975)

4.2 Método 1: (P2.5; P97.5)

El primer método asume que los percentiles calculados sobre los promedios de las muestras son suficientes para realizar la estimación. Por lo tanto, se muestran los percentiles 2.5 y 97.5

cat("Intervalo de confianza desde", q025)
## Intervalo de confianza desde 4.889563
cat("Intervalo de confianza hasta", q975)
## Intervalo de confianza hasta 6.3475

4.3 Método 2: (2X−P97.5; 2X−P2.5)

En el segundo método se realiza el cálculo de dos veces el promedio de los promedios de las muestras, operado contra los percentiles anteriormente calculados.

liminf <- 2 * X_mean - q975
limsup <- 2 * X_mean - q025

cat("Intervalo de confianza desde", liminf)
## Intervalo de confianza desde 4.721665
cat("Intervalo de confianza hasta", limsup)
## Intervalo de confianza hasta 6.179602

4.4 Conclusiones:

El médodo de estimación bootstrap posee ventajas claves en su flexibilidad para estimar estadísticas de interés a partir de datos de muestra. Al generar muestras de la misma distribución que los datos observados, el método bootstrap permite estimar intervalos de confianza con pequeñas muestras.

Además, se identifica que al reutilizar y generar múltiples muestras a partir de los datos originales, se obtienen distribuciones empíricas de las estadísticas de interés, lo que permite calcular intervalos de confianza entendiendo que la distribución de los datos es desconocida.

Sin embargo, si los datos de muestra son pequeños, es posible que las estimaciones bootstrap no reflejen con precisión la verdadera distribución de los datos. Por ejemplo, en el caso de este ejercicio, para los dos métodos, los intervalos de confianza dejarían por fuera tres de siete datos disponibles y que son reales. Es algo importante y a tener en cuenta al realizar estimaciones por este tipo de modelación.