Ejercicio N° 1: Estimación de Pi (π)

En este primer ejecicio de la unidad 2 de la materia de Métodos y Simulación Estadistica, se plantea la estimación del número π usando métodos de distribución estadistica, concretamente la distribución uniforme.

Este problema logra solucionarse usando el área del circulo dentro de un cuadrado, en un problema común, el área del un circulo con radio 0.5 unidades dentro de un cuadrado cuyas aristas miden 1 unidad, representa en proporción π/4 del área total del cuadrado. Si se usa una distribución aleatoria uniforme dentro del cuadrado, la cantidad de puntos dentro del circulo tendra la misma propoción respecto al cuadrado que el área.

Teniendo en cuenta lo anterior, se procedió a generar una distribución uniforme aleatoria para la variable X y Y de valores entre 0 y 1, con ellas se crearon parejas de coordenadas que forman puntos con el código que se muestra a continuación. Teoricamente, entre mayor sea el numero de puntos generados, mas cercano será el valor predicho al valor real de π, por lo tanto se determinó usar un valor de n = 1,000,000.

#Importación de la información 

set.seed(1)

n = 1000000

X_ValuesRandom = runif(n = n, min = 0, max = 1)
Y_ValuesRandom = runif(n = n, min = 0, max = 1)
df = data.frame(Coord1 = X_ValuesRandom, Coord2 = Y_ValuesRandom)

Una vez realizada la simulación de los datos, se procede a crear una nueva columna que me diga si cada una de las parejas evaluadas se encuentra adentro (1) o fuera del circulo (0), para ello se crea una función que permita ver a que distancia se encuentra el punto del centro, si este es menor o igual a 0.25, que es el radio al cuadrado, se determinará que el punto está dentro del circulo.A continuación se muestra el código:

Distancia1 = function(n_Row){
  
  Distancia2 = (n_Row["Coord1"]-0.5)^2 + (n_Row["Coord2"]-0.5)^2
  
  if (Distancia2 <= 0.25){
    r = 1
  }else {
    r = 0
  } 
  return(r)
}

df$Distancia = apply(df[c("Coord1", "Coord2")], 1, FUN = Distancia1)

Finalmente, se procede a hacer el conteo de números que hay dentro del circulo, debido a que la proporcion debería ser la misma que entre las áreas, el número de valores dentro se debe dividir entre el número total de parejas, que finalmente despejando la ecuación y multiplicando por 4, arroja el valor más cercano a π (3.141592) que podemos encontrar en la simulación que es 3.142512.

Conteo_Dentro_Circulo = sum(df$Distancia)

Valor_Previsto = Conteo_Dentro_Circulo / n * 4

print(Valor_Previsto)
## [1] 3.142512


Ejercicio N° 2: Calculo caracteristicas de estimadores y determinación del mejor estimador

El siguiente problema pretende evaluar diferentes estimadores estadísticos con el fin de obtener cual es el mejor dentro de una serie de datos generados aleatoriamente basado en una distribución exponencial, además de comprobar la hipótesis que enuncia que entre mayor sea el tamaño de lamuestra, el estimador será más consistente.

Inicialmente se decidió hacer una función que permita realizar una distribución Exponencial teniendo en cuenta parametros como lamda, el cual es igual a 1 / (Media), se determinó un valor de media equivalente a 0.5 con el fin de que los estimadores converjan en ese valor y a partir de esta medida, determinar la insesgadez, eficiencia y consistencia.

Teniendo en cuenta lo anterior, dentro de la función Generador_Boxplot cuyo parametro principal es el numero de filas que se quiere generar de manera aleatoria cada una de las 4 variables con las que se calculan los estimadores. En este caso los estimadores evaluados son los siguientes:

Luego de generado el dataframe con las variables, se crean columnas con el valor de cada uno de lo estimadores por fila. A partir de este punto, el objetivo final de la función es crear diagramas de caja y bigotes que permitan observar que tan lejos se encuentra la media del estimador, de la media real (0.5).

set.seed(11)

Generador_Boxplots = function(x){
n_rows = x
n_cols = 4
Lamda = 1/0.5

# Crear un DataFrame con 4 columnas de datos exponenciales

df_1 = data.frame(
  X1 = rexp(n_rows, rate = Lamda),
  X2 = rexp(n_rows, rate = Lamda),
  X3 = rexp(n_rows, rate = Lamda),
  X4 = rexp(n_rows, rate = Lamda)
)

# Se crean columnas que aplican los estimadores y formulas a cada una de las filas

df_2 = df_1 %>%
  mutate(Theta1 = ((X1 + X2)/6) + ((X3 + X4)/3))

df_2 = df_2 %>%
  mutate(Theta2 = (X1 + 2*X2 + 3*X3 + 4*X4)/5 )

df_2 = df_2 %>%
  mutate(Theta3 = (X1 + X2 + X3 + X4)/4)

df_2 = df_2 %>%
  mutate(Theta4 = ((pmax(X1, X2, X3, X4) + pmin(X1, X2, X3, X4))/2))

df_estimadores <- df_2[c("Theta1", "Theta2", "Theta3", "Theta4")]

#Se crean los boxplots de cada uno de los estimadores

boxplot(df_estimadores, main = paste("Boxplot con estimadores para n =" ,x) , ylab = "Valor del Estimador",
        col = c("orange", "green", "blue", "yellow"))

abline(h = 0.5, col = "red", lty = 2)


print(colMeans(df_estimadores))
varianzas_estimadores = apply(df_estimadores, 2, var)
print(varianzas_estimadores)

}

Una vez creada la función, se procede a evaluarla con diferentes tamaños de muestra, inicialmente se evaluó con un tamaño de 20 filas, en este caso se pudo ver que theta 1, theta 2 y theta 3 presentan valores muy cercanos a la media, lo que implica que son mejores estimadores que 2, el cual se encuentra mas lejano. También se observa que todos los estimadores son sesgados ya que conservan una diferencia con la media real (0.5) sin embargo se denota que theta 1 es el estimador con menor sesgadez y además el de menor varianza (0.0727), por lo tanto se considera el más eficiente.

Generador_Boxplots(20)

##    Theta1    Theta2    Theta3    Theta4 
## 0.5868607 1.1399254 0.6046455 0.7015588 
##     Theta1     Theta2     Theta3     Theta4 
## 0.07272365 0.26780555 0.08168032 0.15674729

En segundo caso evaluado, se presenta un tamaño de muestra de 50, un incremento considerable respecto a la inicial. En este caso se muestra una disminución entre las medias de todos los estimadores, esto podria representar consistencia, sin embargo es pertinente hacer pruebas con números de muestra mayores para comprobar la hipótesis. En esta simulación el mejor y el más eficiente pasa a ser theta 3 ya que conserva la menor diferencia respecto a la media y la menor varianza (0.0553).

Generador_Boxplots(50)

##    Theta1    Theta2    Theta3    Theta4 
## 0.4659159 0.9420917 0.4702298 0.5492470 
##     Theta1     Theta2     Theta3     Theta4 
## 0.05988711 0.28982943 0.05530593 0.08721238

Posteriormente se realiza la simulación con un estimador que representa el doble de muestras que el anterior. En este caso se continua denotando una disminución en el sesgo de theta 1 y theta 3 por lo tanto se puede considerar que estos estimadores son consistentes, ya que a medida que se aumenta el tamaño de la muestra estos se acercan más la media real, Caso diferente a lo que ocurre con theta 2 y theta 4, ya que en una primera prueba disminuyeron su nivel de sesgo, pero en una muestra mayor presentaron una diferencia mayor, por lo tanto no se consideran consistentes. En esta simulación el estimador más eficiente es theta 3 ya que presenta menor varianza (0.0724) respecto a theta 1 (0.0807).

Generador_Boxplots(100)

##    Theta1    Theta2    Theta3    Theta4 
## 0.5180562 1.0510828 0.5104381 0.5830582 
##     Theta1     Theta2     Theta3     Theta4 
## 0.08072583 0.33970688 0.07240393 0.11133134

Finalmente, la última simulación contempla un tamaño de muestra de 1,000. Al igual que la anterior se denota consistencia en los estimadores theta 1 y theta 3 y aunque los niveles de sesgo de theta 2 y theta 4 disminuyeron, previamente se comprobó que no son estimadores consistentes. Con esta simulación se puede concluir que theta 3 es el mejor estimador ya que tiene el menor sesgo, la mayor eficiencia y es el estimador mas consistente durante las simulaciones realizadas. también se pudo comprobar la hipótesis inicial en donde se denota que entre mayor sea el tamaño de la muestra, algunos estimadores disminuyen su nivel de sesgo y por lo tanto se hace consistente, propiciando una mejor aproximación al valor real en caso de implementarse.

Generador_Boxplots(1000)

##    Theta1    Theta2    Theta3    Theta4 
## 0.4971283 0.9946601 0.4942797 0.5753517 
##     Theta1     Theta2     Theta3     Theta4 
## 0.06418810 0.27475953 0.05763378 0.09537147


Ejercicio N° 3: Demostración Teorema del Límite Central (TLC)

Este ejercicio pretende demostrar el teorema del límite central el cual enuncia que a medida que el tamaño de la muestra se incrementa, la media muestral se acercará a la media de la población. Por tanto, mediante el TCL podemos definir la distribución de la media muestral de una determinada población con una varianza conocida. De manera que la distribución seguirá una distribución normal si el tamaño de la muestra es lo suficientemente grande. Para demostrar esta hipótesis, se creó una población aleatoria de 1,000 plantas enfermas con una distribución binomial y una probabilidad de exito del 50% de que estuviera enferma.

set.seed(5)

df_Plantas_1000 = data.frame (Enfermas = rbinom(n = 1000,size = 1, prob = 0.5))
Datos_P05 = df_Plantas_1000$Enfermas

Una vez creada la población, se realizó una función que permitiera tomar una muestra de la población inicial, hacerlo un total de 500 veces, crear un data frame con aquella información y finalmente obtener la media de cada una de las filas con las muestras obtenidas, esta media será el objeto de estudio en esta simulación. El parametro de inicio de la función permite determinar la cantidad de muestras que va a tomar la función por cada fila, siendo esta la parte clave del experimento.

Funcion_Muestreo = function(Datos, Numero_Muestras){
  
  matriz_datos = matrix(nrow = 500, ncol = Numero_Muestras)
  
  for (i in 1:500) {
  
    Vector_Muestras = sample(Datos, size = Numero_Muestras , replace = FALSE, prob = NULL)
    matriz_datos[i, ] = Vector_Muestras
  }
  df_Muestra = data.frame (matriz_datos)
  df_Muestra$Promedio_Muestra = apply(df_Muestra, 1, mean)
  
  return(df_Muestra)
}

poteriormente, después de planteada la función principal, se puso a prueba en un único experimento con un tamaño muestral de n = 10, con este se elaboró un histograma que denota la asimetría en la distribución de la media, lo que implica que no es una distribución normal, esto se ratifica con el diagrama qqnorm, el cual permite determinar que tan similares son los valores de los cuartiles de una distribución normal, con los de la distribución evaluada, denotando que aun no se comporta como una.

df_Muestras_1_10 = Funcion_Muestreo(Datos_P05,10)

hist(df_Muestras_1_10$Promedio_Muestra,  main = "n=10", xlab = "Promedio", ylab = "Frecuencia", freq=FALSE)

qqnorm(df_Muestras_1_10$Promedio_Muestra,  main = paste("n = 10")) ; qqline(df_Muestras_1_10$Promedio_Muestra, col="blue")

Luego de realizada la primera prueba, se quiso realizar una función que permitiera comparar que pasaría si realizo el mismo experimento con diferentes tamaños de muestra (5, 10, 15, 20, 30, 50, 60, 100, 200, 500) y comparar los resultados en una misma gráfica, para eso se realizo la función (Funcion_Diferentes_Tamanos_hist) como se muestra a continuación:

#Funcion histogramas

Funcion_Diferentes_Tamanos_hist = function(Datos) {
  
  set.seed(15)
  
  Tamanos_Muestra = c(5, 10, 15, 20, 30, 50, 60, 100, 200, 500)
  
  resultados_df = data.frame(Promedio = numeric(0), Varianza = numeric(0))
  
  for (n in Tamanos_Muestra) {
    
    df_Muestra = Funcion_Muestreo(Datos, n)
    Promedio = mean(df_Muestra$Promedio_Muestra)
    Varianza = var(df_Muestra$Promedio_Muestra)
    resultados_df = rbind(resultados_df, data.frame(Promedio = Promedio, Varianza = Varianza))
    
  }
  
  par(mfrow = c(3, 4)) 
  
  resultados_df$Muestra = Tamanos_Muestra
  print(resultados_df)
  
  for (n in Tamanos_Muestra) {
    
    df_Muestra = Funcion_Muestreo(Datos, n)
    hist(df_Muestra$Promedio_Muestra, main = paste("n =", n), xlab = "Promedio", ylab = "Frecuencia", freq = FALSE)
  }
  
  par(mfrow = c(1, 1))

}

# Función qqnorm

Funcion_Diferentes_Tamanos_qqnorm = function(Datos) {
  
  par(mfrow = c(3, 4))  
  
  Tamanos_Muestra = c(5, 10, 15, 20, 30, 50, 60, 100, 200, 500)
  
  for (n in Tamanos_Muestra) {
    df_Muestra = Funcion_Muestreo(Datos, n)
    qqnorm(df_Muestra$Promedio_Muestra, main = paste("n =", n)) ; qqline(df_Muestra$Promedio_Muestra, col="blue")
  }
  
  par(mfrow = c(1, 1))
}

Con los resultados obtenidos, y con las caracteristicas inicialmente planteadas de la población inicial, se puede apreciar que a partir de el tamaño de muestra n = 60 se hace más evidente la tranformación de la distribución de las medias de la muestra, evidenciando que a medida que aumenta el tamaño de la muestra se parece más a una distribución normal, lo cual es ratificado por las gráficas qqnorm, algo también a destacar, es que el promedio de la media en cada uno de los tamaños evaluados no tiene mucha variación entre ellos, sin embargo, la varianza si se va haciendo mas pequeña cada vez, haciendo que la estimación sea cada vez mas eficiente.

Funcion_Diferentes_Tamanos_hist(Datos_P05)
##     Promedio     Varianza Muestra
## 1  0.4760000 0.0488817635       5
## 2  0.4916000 0.0233361122      10
## 3  0.4920000 0.0170990871      15
## 4  0.4894000 0.0113403206      20
## 5  0.4979333 0.0083925139      30
## 6  0.4889600 0.0047692569      50
## 7  0.4887667 0.0039745391      60
## 8  0.4920000 0.0021631263     100
## 9  0.4924500 0.0010254985     200
## 10 0.4908880 0.0002789854     500

Funcion_Diferentes_Tamanos_qqnorm(Datos_P05)

En una segunda prueba se hizo el mismo proceso anterior, esta vez con una probabilidad de exito de 30%, obteniendo resultados similares sin embargo, de observa que a partir de n = 50 se comienza a denotar una similitud mayor con una gráfica de distribución normal.

Datos_P03 = data.frame (Enfermas = rbinom(n = 1000,size = 1, prob = 0.3))
Datos_P03 = Datos_P03$Enfermas

Funcion_Diferentes_Tamanos_hist(Datos_P03)
##    Promedio     Varianza Muestra
## 1  0.305600 0.0436559519       5
## 2  0.282400 0.0192086573      10
## 3  0.290800 0.0127319350      15
## 4  0.285100 0.0098226353      20
## 5  0.291600 0.0077404765      30
## 6  0.291440 0.0038848962      50
## 7  0.286600 0.0032280071      60
## 8  0.290160 0.0017250244     100
## 9  0.290290 0.0007751161     200
## 10 0.289712 0.0001988107     500

Funcion_Diferentes_Tamanos_qqnorm(Datos_P03)

Finalmente se hizo la prueba con una probabilidad de exito de 90%, pero contrario al ejemplo anterior, se requirió un tamaño de muestra de n = 200 para obtener un resultado visible de distribución normal.

Datos_P09 = data.frame (Enfermas = rbinom(n = 1000,size = 1, prob = 0.9))
Datos_P09 = Datos_P09$Enfermas

Funcion_Diferentes_Tamanos_hist(Datos_P09)
##     Promedio     Varianza Muestra
## 1  0.8904000 1.954693e-02       5
## 2  0.9030000 8.187375e-03      10
## 3  0.8973333 5.559564e-03      15
## 4  0.8994000 4.398437e-03      20
## 5  0.9004667 3.168341e-03      30
## 6  0.8963200 1.679416e-03      50
## 7  0.8966667 1.515253e-03      60
## 8  0.8980800 8.508152e-04     100
## 9  0.8988000 3.608818e-04     200
## 10 0.8991800 9.733828e-05     500

Funcion_Diferentes_Tamanos_qqnorm(Datos_P09)


Ejercicio N° 4: Estimación con Método de Boostrap

Este ejercicio busca mostrar como utilizando el método de Boostrap se puede reconstruir una población a partir de una muestra, en este caso la muestra representa la eficiencia de combustible en millas/galón de una serie de camiones, en este caso se adquirió una muestra de 7 camiones, con base en esta muestra se pretende reconstruir la población inicial, por esto, inicialmente se crea el vector con las muestras y se hace un muestreo con reemplazo (replace = TRUE), esto implica que se podra obtener de manera aleatoria un número de los 7 iniciales permitiendo repitir la muestras y esto se repetirá un total de 7,000 veces.

set.seed(27)

Muestra_Inicial = c( 7.69, 4.97, 4.56, 6.49, 4.34, 6.24, 4.45) 

Metodo_Boot =sample(Muestra_Inicial,7000,replace=TRUE)

Posteriormente se divide la población obtenida por Bootstrap en 1,000 filas, cada una dividida en 7 columnas, donde cada fila representa una muestra a la que se aplicará la media, esto con el fin de crear un vector de medias con longitud de 1,000.

Division_Boot = matrix( Metodo_Boot, nrow=1000, ncol=7) 

Promedio_Matrix_Boot=apply(Division_Boot, 1, mean)

Se utilizan dos métodos para comprobar la efectividad del experimento, el primero mide el punto donde la probabilidad es del 97.5 % y del 2.5%, encontrando así que se tiene una confianza del 97.5% de que el valor de la media se ubique entre 4.7527 y 6.4602.

Metodo_1 = quantile(Promedio_Matrix_Boot, probs=c(0.025, 0.975)) 

print(Metodo_1)
##    2.5%   97.5% 
## 4.75275 6.46025

El método 2 para construir el intervalo de confianza utiliza la formula [2 * Promedio de la muestra - resultado método 1 (97,5%)] y [2 * Promedio de la muestra - resultado método 1(2.5%)], en este caso se encontró que el intervalo de confianza del 97,5% se encuentra entre 4.60713 y 6.31463.

Metodo_2 = c( 2*mean(Promedio_Matrix_Boot)-Metodo_1[2], 2*mean(Promedio_Matrix_Boot)-Metodo_1[1]) 

print(Metodo_2)
##   97.5%    2.5% 
## 4.60713 6.31463

Finalmente se grafican los resultados obtenidos, encontrando que la distribución de las medias del muestreo se comporta como una distribución normal, además, se muestran las lineas con los respectivos intervalos de confianza, el método 1 en color rojo y el método 2 con color verde.

hist(Promedio_Matrix_Boot, las=1, main=" ", ylab = " ", xlab = " ", col="blue")

abline(v=Metodo_1, col="red",lwd=2)
abline(v=Metodo_2, col="green",lwd=2)