Actividad 2, Simulación R
Jesus Olier
2023-09-17
Problema 1
Estimación del valor de π La siguiente figura sugiere como estimar el valor de π con una simulación. En la figura (imagen 1), 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 π.
set.seed(1) # Establecer una semilla para reproducibilidad
nRowGen = 100000 # Número total de puntos a generar
X_Values = runif(n = nRowGen, min = 0, max = 1)
Y_Values = runif(n = nRowGen, min = 0, max = 1)
R_Values = runif(n = nRowGen, min = 0, max = 0)
dataFrame = data.frame(X = X_Values, Y = Y_Values, R = R_Values)
dataFramelibrary(ggplot2) # Carga la librería ggplot2
# Crear un gráfico de dispersión
ggplot(data = dataFrame, aes(x = X, y = Y)) +
geom_point(size = 1, alpha = 0.5) + labs(title = "Estimación de π con Puntos Aleatorios", x = "Eje X", y = "Eje Y") +
# Agregar una circunferencia unitaria
geom_path(data = data.frame(x = cos(seq(0, 2 * pi, length.out = 100)),
y = sin(seq(0, 2 * pi, length.out = 100))),
aes(x, y), color = "blue", size = 1)
# Función para establecer R en 1 para puntos dentro de la circunferencia
Puntos_Dentro <- function(df) {
df$R[df$X^2 + df$Y^2 <= 0.25] <- 1
return(df)
}
# Llamar a la función con el dataframe y Modificandolo asignando los puntos dentro y fuera de la circunferencia
dataFrame <- Puntos_Dentro(dataFrame)
dataFrame# Crear un gráfico de dispersión
ggplot(data = df_PuntosDentro, aes(x = X, y = Y)) +
geom_point(size = 1, alpha = 0.5) + labs(title = "Estimación de π con Puntos dentro de la Estimacion", x = "Eje X", y = "Eje Y") +
# Agregar una circunferencia unitaria
geom_path(data = data.frame(x = cos(seq(0, 2 * pi, length.out = 100)),
y = sin(seq(0, 2 * pi, length.out = 100))),
aes(x, y), color = "blue", size = 1)
# Crear un gráfico de dispersión
ggplot(data = df_PuntosFuera, aes(x = X, y = Y)) +
geom_point(size = 1, alpha = 0.5) + labs(title = "Estimación de π con Puntos fuera de la Estmacion", x = "Eje X", y = "Eje Y") +
# Agregar una circunferencia unitaria
geom_path(data = data.frame(x = cos(seq(0, 2 * pi, length.out = 100)),
y = sin(seq(0, 2 * pi, length.out = 100))),
aes(x, y), color = "blue", size = 1)
nRowIn <- nrow(df_PuntosDentro)
nRowOut <- nrow(df_PuntosFuera)
Estimacion = (nRowIn/nRowGen*4)
cat("Número de registros en la circunferencia:", nRowIn, "\n")## Número de registros en la circunferencia: 19569## Número de registros Fuera circunferencia: 80431## Estimación de π: 0.78276Problema 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
Genere una muestras de n=20, 50, 100 y 1000 para cada uno de los estimadores planteados. En cada caso evalue las propiedades de insesgadez, eficiencia y consistencia Suponga un valor para el parámetro θ
## Loading required package: dplyr##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionGenerador_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(Teta1 = ((X1 + X2)/6) + ((X3 + X4)/3))
df_2 = df_2 %>%
mutate(Teta2 = (X1 + 2*X2 + 3*X3 + 4*X4)/5 )
df_2 = df_2 %>%
mutate(Teta3 = (X1 + X2 + X3 + X4)/4)
df_2 = df_2 %>%
mutate(Teta4 = ((pmax(X1, X2, X3, X4) + pmin(X1, X2, X3, X4))/2))
df_estimadores <- df_2[c("Teta1", "Teta2", "Teta3", "Teta4")]
#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("red", "blue", "orange", "purple"))
print(colMeans(df_estimadores))
varianzas_estimadores = apply(df_estimadores, 2, var)
print(varianzas_estimadores)
}
## Teta1 Teta2 Teta3 Teta4
## 0.5351541 1.0472595 0.5313409 0.6298271
## Teta1 Teta2 Teta3 Teta4
## 0.06189673 0.29816264 0.06120700 0.10850040
## Teta1 Teta2 Teta3 Teta4
## 0.5247243 1.0246726 0.5192176 0.6052899
## Teta1 Teta2 Teta3 Teta4
## 0.06920607 0.28154585 0.06329876 0.08109426
## Teta1 Teta2 Teta3 Teta4
## 0.5169422 1.0301765 0.5173151 0.6083772
## Teta1 Teta2 Teta3 Teta4
## 0.06437965 0.26179119 0.05933459 0.08841087
## Teta1 Teta2 Teta3 Teta4
## 0.4964142 0.9914359 0.4930187 0.5767244
## Teta1 Teta2 Teta3 Teta4
## 0.06847236 0.28928552 0.06028784 0.09760957cat("En resumen, Teta1 y Teta3 tienden a ser estimadores insesgados y consistentes para θ, independientemente del tamaño de muestra, y tienen errores estándar más bajos a medida que n aumenta. Teta2 tiende a ser un estimador sesgado con errores estándar más altos. Teta4 también es insesgado y consistente, pero con un error estándar más bajo que Teta2. El tamaño de muestra tiene un impacto en la precisión de las estimaciones y en el grado de sesgo en el caso de θ̂2.")## En resumen, Teta1 y Teta3 tienden a ser estimadores insesgados y consistentes para θ, independientemente del tamaño de muestra, y tienen errores estándar más bajos a medida que n aumenta. Teta2 tiende a ser un estimador sesgado con errores estándar más altos. Teta4 también es insesgado y consistente, pero con un error estándar más bajo que Teta2. El tamaño de muestra tiene un impacto en la precisión de las estimaciones y en el grado de sesgo en el caso de θ̂2.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 y Calcule el estimador de la proporción muestral 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.
#xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
set.seed(1) # Establecer una semilla para reproducibilidad
# Generar una población de 1000 plantas con un 50% de probabilidad de estar enfermas
poblacion <- rbinom(n = 1000, size = 1, prob = 0.5)
# Función para obtener muestra aleatoria y calcular pˆ
calcular_p_muestra <- function(poblacion, n) {
muestra <- sample(poblacion, size = n, replace = FALSE)
p_muestra <- mean(muestra)
return(p_muestra)
}
n_simulaciones <- 500
tamanos_muestra <- c(5, 10, 15, 20, 30, 50, 60, 100, 200, 500)
resultados <- matrix(NA, nrow = n_simulaciones, ncol = length(tamanos_muestra))
for (i in 1:n_simulaciones) {
for (j in 1:length(tamanos_muestra)) {
n <- tamanos_muestra[j]
p_muestra <- calcular_p_muestra(poblacion, n)
resultados[i, j] <- p_muestra
}
}
# Función para realizar pruebas de normalidad y generar gráfico QQ
evaluar_normalidad <- function(vector) {
shapiro_test <- shapiro.test(vector)
qq_plot <- qqnorm(vector)
return(list(ShapiroTest = shapiro_test, QQPlot = qq_plot))
}
# Aplicar las pruebas de normalidad y gráficos QQ a los resultados
resultados_normalidad <- lapply(1:length(tamanos_muestra), function(i) {
evaluar_normalidad(resultados[, i])
})
# Función para realizar toda la simulación
simulacion <- function(prob_enfermas) {
# Generar población con la probabilidad dada
poblacion <- rbinom(n = 1000, size = 1, prob = prob_enfermas)
# Repetir el proceso y evaluar normalidad
resultados <- matrix(NA, nrow = n_simulaciones, ncol = length(tamanos_muestra))
for (i in 1:n_simulaciones) {
for (j in 1:length(tamanos_muestra)) {
n <- tamanos_muestra[j]
p_muestra <- calcular_p_muestra(poblacion, n)
resultados[i, j] <- p_muestra
}
}
resultados_normalidad <- lapply(1:length(tamanos_muestra), function(i) {
evaluar_normalidad(resultados[, i])
})
return(list(Resultados = resultados, Normalidad = resultados_normalidad))
}
# Realizar la simulación para probabilidad de enfermas del 10%
simulacion_10 <- simulacion(0.1)
# Crear un solo gráfico de Normal Q-Q Plot para mostrar los resultados por tamaño de muestra
par(mfrow = c(2, 5), mar = c(4, 4, 2, 1))
for (i in 1:length(tamanos_muestra)) {
n <- tamanos_muestra[i]
qqnorm(simulacion_10$Resultados[, i], main = paste("n =", n),
xlab = "Cuantiles teóricos", ylab = "Cuantiles observados", col = "blue")
qqline(simulacion_10$Resultados[, i], col = "red")
qqnorm(simulacion_90$Resultados[, i], main = paste("n =", n),
xlab = "Cuantiles teóricos", ylab = "Cuantiles observados", col = "blue")
qqline(simulacion_90$Resultados[, i], col = "red")
}
# Restaurar el diseño de la ventana gráfica
par(mfrow = c(1, 1))
# Crear un solo gráfico de cajas para mostrar los resultados por tamaño de muestra
par(mfrow = c(2, 5), mar = c(4, 4, 2, 1))
for (i in 1:length(tamanos_muestra)) {
n <- tamanos_muestra[i]
boxplot(simulacion_10$Resultados[, i], main = paste("n =", n),
ylab = "pˆ", ylim = c(0, 1), col = "lightblue", border = "blue")
boxplot(simulacion_90$Resultados[, i], main = paste("n =", n),
ylab = "pˆ", ylim = c(0, 1), col = "lightpink", border = "red")
}
Problema 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. 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:
datos <- c(7.69, 4.97, 4.56, 6.49, 4.34, 6.24, 4.45)
n <- length(datos)
k <- 1000
alpha <- 0.05
muestras_bootstrap <- replicate(k, sample(datos,replace=TRUE))
medias_bootstrap <- apply(muestras_bootstrap,2,mean)
intervalo_metodo1 <- quantile(medias_bootstrap,c(alpha/2,1-alpha/2))
media_muestral <- mean(datos)
intervalo_metodo2 <- c(2 * media_muestral-quantile(medias_bootstrap,1-alpha/2),2* media_muestral-quantile(medias_bootstrap,alpha/2))
cat("Intervalo de confianza (Método 1):", intervalo_metodo1, "")## Intervalo de confianza (Método 1): 4.78 6.5265## Intervalo de confianza (Método 2): 4.542071 6.288571ambos métodos proporcionan resultados coherentes y confiables para la estimación del intervalo de confianza de la media de los datos. El Método 2 tiene la ventaja de tener en cuenta la distribución t-Student, lo que lo hace más robusto en situaciones en las que la normalidad de los datos no está garantizada. En cualquier caso, ambos métodos indican que podemos tener un alto grado de confianza en que la media poblacional se encuentra dentro de los intervalos calculados.