Actividad 2
Caris Chia Amaya - Weimar Cortes Montiel
Métodos y Simulación estadística
Maestría en Ciencia de Datos
Pontificia Universidad Javeriana de Cali
12 de septiembre de 2024
Problema 1
Estimación del valor de π
La siguiente figura sugiere como estimar el valor de π con una simulación. En la figura, un círcuito 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 éste, 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 π.
plot(0,0, type = "n", xlim = c(0,1), ylim = c(0,1),
xlab = "Eje X", ylab = "Eje Y", asp = 1,
main = "Círculo Inscrito en un Cuadrado")
symbols(0.5,0.5, circles = 0.5, add = TRUE, inches = FALSE, lwd = 2, col = "blue")
Pasos sugeridos:
Genere n coordenadas \(x: X_1, \ldots, X_n\). 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)\).
Genere 1000 coordenadas \(y: Y_1, \ldots, Y_n\), utilizando nuevamente la distribución uniforme con valor mínimo de 0 y valor máximo de 1.
Cada punto \((X_i, Y_i)\) se encuentra dentro del círculo si su distancia desde el centro \((0.5, 0.5)\) es menor a 0.5. Para cada par \((X_i, Y_i)\) determine si la distancia desde el centro es menor que 0.5. Esto último se puede realizar al calcular el valor \((X_i - 0.5)^2 + (Y_i - 0.5)^2\), que es el cuadrado de la distancia, y al determinar si es menor que 0.25.
¿Cuántos de los puntos están dentro del círculo? ¿Cuál es su estimación de \(\pi\)?
Nota
Con sólo 1000 puntos, es probable que la estimación presente un error de 0.05 o más. Una simulación con 10000 y 100000 puntos tiene mayores probabilidades de dar como resultado una estimación muy cercana al valor verdadero.
Funciones recomendadas:
runif(),function(){}.Entregable: enlace en RPubs con informe 1.
Problema tomado de Navidi (2006).
Análisis
En este análisis, se utiliza el método de Monte Carlo para estimar el valor de \(\pi\) generando puntos aleatorios en un cuadrado de lado 1 y contando cuántos de esos puntos caen dentro de un círculo inscrito. Este enfoque permite aproximar el valor de \(\pi\) a medida que aumenta la cantidad de puntos generados.
A continuación, se muestran gráficos que representan las estimaciones de \(\pi\) para diferentes tamaños de muestra: 100, 1000, 10,000, 100,000 y 1,000,000 puntos. En cada gráfico, los puntos dentro del círculo se muestran en azul, mientras que los puntos fuera del círculo aparecen en rojo. Acompañando los gráficos, se presenta una tabla que muestra el número de puntos dentro del círculo, la estimación de \(\pi\), el error absoluto y el error relativo para cada simulación.
El objetivo es observar cómo la estimación de \(\pi\) mejora conforme aumenta el número de puntos generados, lo que permite una mayor precisión en la estimación y una reducción en los errores.
n <- 100
x <- runif(n, min = 0, max = 1)
y <- runif(n, min = 0, max = 1)
distancia_cuadrada <- (x - 0.5)^2 + (y - 0.5)^2
dentro_del_circulo <- distancia_cuadrada < 0.25
num_dentro_del_circulo <- sum(dentro_del_circulo)
pi_estimado <- 4 * num_dentro_del_circulo / n
plot(x, y, col = ifelse(dentro_del_circulo, "#164C51", "#701c00"), asp = 1, main = "Estimación de π usando Monte Carlo (100)")
symbols(0.5, 0.5, circles = 0.5, add = TRUE, inches = FALSE)
pi_real <- pi
error_absoluto <- abs(pi_estimado - pi_real)
error_relativo <- error_absoluto / pi_real
resultados_df <- data.frame(
'Puntos Dentro Circulo' = num_dentro_del_circulo,
'Estimacion de Pi' = round(pi_estimado, 6),
'Error Absoluto' = round(error_absoluto, 6),
'Error Relativo' = round(error_relativo, 6)
)
pander(resultados_df, caption = "Resultados de la Estimación de π con 100 datos", digits = 6, split.tables = Inf)| Puntos.Dentro.Circulo | Estimacion.de.Pi | Error.Absoluto | Error.Relativo |
|---|---|---|---|
| 71 | 2.84 | 0.301593 | 0.096 |
n <- 1000
x <- runif(n, min = 0, max = 1)
y <- runif(n, min = 0, max = 1)
distancia_cuadrada <- (x - 0.5)^2 + (y - 0.5)^2
dentro_del_circulo <- distancia_cuadrada < 0.25
num_dentro_del_circulo <- sum(dentro_del_circulo)
pi_estimado <- 4 * num_dentro_del_circulo / n
plot(x, y, col = ifelse(dentro_del_circulo, "#164C51", "#701c00"), asp = 1, main = "Estimación de π usando Monte Carlo (1000)")
symbols(0.5, 0.5, circles = 0.5, add = TRUE, inches = FALSE)
pi_real <- pi
error_absoluto <- abs(pi_estimado - pi_real)
error_relativo <- error_absoluto / pi_real
resultados_df <- data.frame(
'Puntos Dentro Circulo' = num_dentro_del_circulo,
'Estimacion de Pi' = round(pi_estimado, 6),
'Error Absoluto' = round(error_absoluto, 6),
'Error Relativo' = round(error_relativo, 6)
)
pander(resultados_df, caption = "Resultados de la Estimación de π con 1000 datos", digits = 6, split.tables = Inf)| Puntos.Dentro.Circulo | Estimacion.de.Pi | Error.Absoluto | Error.Relativo |
|---|---|---|---|
| 778 | 3.112 | 0.029593 | 0.00942 |
n <- 10000
x <- runif(n, min = 0, max = 1)
y <- runif(n, min = 0, max = 1)
distancia_cuadrada <- (x - 0.5)^2 + (y - 0.5)^2
dentro_del_circulo <- distancia_cuadrada < 0.25
num_dentro_del_circulo <- sum(dentro_del_circulo)
pi_estimado <- 4 * num_dentro_del_circulo / n
plot(x, y, col = ifelse(dentro_del_circulo, "#164C51", "#701c00"), asp = 1, main = "Estimación de π usando Monte Carlo (10000)")
symbols(0.5, 0.5, circles = 0.5, add = TRUE, inches = FALSE)
pi_real <- pi
error_absoluto <- abs(pi_estimado - pi_real)
error_relativo <- error_absoluto / pi_real
resultados_df <- data.frame(
'Puntos Dentro Circulo' = num_dentro_del_circulo,
'Estimacion de Pi' = round(pi_estimado, 6),
'Error Absoluto' = round(error_absoluto, 6),
'Error Relativo' = round(error_relativo, 6)
)
pander(resultados_df, caption = "Resultados de la Estimación de π con 10.000 datos", digits = 6, split.tables = Inf)| Puntos.Dentro.Circulo | Estimacion.de.Pi | Error.Absoluto | Error.Relativo |
|---|---|---|---|
| 7809 | 3.1236 | 0.017993 | 0.005727 |
n <- 100000
x <- runif(n, min = 0, max = 1)
y <- runif(n, min = 0, max = 1)
distancia_cuadrada <- (x - 0.5)^2 + (y - 0.5)^2
dentro_del_circulo <- distancia_cuadrada < 0.25
num_dentro_del_circulo <- sum(dentro_del_circulo)
pi_estimado <- 4 * num_dentro_del_circulo / n
plot(x, y, col = ifelse(dentro_del_circulo, "#164C51", "#701c00"), asp = 1, main = "Estimación de π usando Monte Carlo (100000)")
symbols(0.5, 0.5, circles = 0.5, add = TRUE, inches = FALSE)
pi_real <- pi
error_absoluto <- abs(pi_estimado - pi_real)
error_relativo <- error_absoluto / pi_real
resultados_df <- data.frame(
'Puntos Dentro Circulo' = num_dentro_del_circulo,
'Estimacion de Pi' = round(pi_estimado, 6),
'Error Absoluto' = round(error_absoluto, 6),
'Error Relativo' = round(error_relativo, 6)
)
pander(resultados_df, caption = "Resultados de la Estimación de π con 100.000 datos", digits = 6, split.tables = Inf)| Puntos.Dentro.Circulo | Estimacion.de.Pi | Error.Absoluto | Error.Relativo |
|---|---|---|---|
| 78402 | 3.13608 | 0.005513 | 0.001755 |
n <- 1000000
x <- runif(n, min = 0, max = 1)
y <- runif(n, min = 0, max = 1)
distancia_cuadrada <- (x - 0.5)^2 + (y - 0.5)^2
dentro_del_circulo <- distancia_cuadrada < 0.25
num_dentro_del_circulo <- sum(dentro_del_circulo)
pi_estimado <- 4 * num_dentro_del_circulo / n
plot(x, y, col = ifelse(dentro_del_circulo, "#164C51", "#701c00"), asp = 1, main = "Estimación de π usando Monte Carlo (1000000)")
symbols(0.5, 0.5, circles = 0.5, add = TRUE, inches = FALSE)
pi_real <- pi
error_absoluto <- abs(pi_estimado - pi_real)
error_relativo <- error_absoluto / pi_real
resultados_df <- data.frame(
'Puntos Dentro Circulo' = num_dentro_del_circulo,
'Estimacion de Pi' = round(pi_estimado, 6),
'Error Absoluto' = round(error_absoluto, 6),
'Error Relativo' = round(error_relativo, 6)
)
pander(resultados_df, caption = "Resultados de la Estimación de π con 1'000.000 datos", digits = 6, split.tables = Inf)| Puntos.Dentro.Circulo | Estimacion.de.Pi | Error.Absoluto | Error.Relativo |
|---|---|---|---|
| 784851 | 3.1394 | 0.002189 | 0.000697 |
Resultados
Se procede a analizar las tendencias de las muestras en una tabla, para observar cómo, a medida que se incrementa el tamaño de las muestras, la estimación se aproxima cada vez más a π.
options(scipen = 999)
estimacion_pi <- function(n) {
x <- runif(n, min = 0, max = 1)
y <- runif(n, min = 0, max = 1)
distancia_cuadrada <- (x - 0.5)^2 + (y - 0.5)^2
dentro_del_circulo <- distancia_cuadrada < 0.25
num_dentro_del_circulo <- sum(dentro_del_circulo)
pi_estimado <- 4 * num_dentro_del_circulo / n
pi_real <- pi
error_absoluto <- abs(pi_estimado - pi_real)
return(c(pi_estimado, error_absoluto))
}
tamanos_muestra <- c(100, 1000, 10000, 100000, 1000000)
resultados <- sapply(tamanos_muestra, function(n) estimacion_pi(n))
resultados_df <- data.frame(
'Tamaño muestra' = tamanos_muestra,
'Estimacion pi' = round(resultados[1, ], 4),
'Error absoluto' = round(resultados[2, ], 4)
)
kable(resultados_df, col.names = c("Tamaño de Muestra", "Estimación de π", "Error Absoluto"),
digits = 6, align = "c", format = "pandoc")| Tamaño de Muestra | Estimación de π | Error Absoluto |
|---|---|---|
| 100 | 3.1200 | 0.0216 |
| 1000 | 3.1280 | 0.0136 |
| 10000 | 3.1268 | 0.0148 |
| 100000 | 3.1358 | 0.0058 |
| 1000000 | 3.1439 | 0.0023 |
En la siguiente gráfica se puede observar una ligera desviación en la estimación en comparación con el tamaño de muestra anterior, está desviación varia entre las muestras, lo cual puede ser atribuido a la fluctuación inherente en la estimación con muestras grandes. A medida que el tamaño de la muestra aumenta, la estimación de π tiende a acercarse más al valor verdadero de π, reduciendo el error absoluto. Aunque la tendencia general muestra una mejora en la estimación con el aumento del tamaño de la muestra, se observa una ligera fluctuación en el error absoluto entre tamaños de muestra grandes (por ejemplo, entre 10,000 y 100,000), lo cual puede deberse a la variabilidad inherente en la simulación y no necesariamente a un problema en el método. Por lo tanto, se procede a revisar los datos para identificar el número de muestras requerido para minimizar el error absoluto a su mínima expresión.
n_sim <- seq(100000, 1000000, by = 10000)
pi_estimados <- numeric(length(n_sim))
errores_absolutos <- numeric(length(n_sim))
for (i in 1:length(n_sim)) {
n <- n_sim[i]
x <- runif(n, min = 0, max = 1)
y <- runif(n, min = 0, max = 1)
distancia_cuadrada <- (x - 0.5)^2 + (y - 0.5)^2
dentro_del_circulo <- distancia_cuadrada < 0.25
num_dentro_del_circulo <- sum(dentro_del_circulo)
pi_estimado <- 4 * num_dentro_del_circulo / n
pi_estimados[i] <- pi_estimado
errores_absolutos[i] <- abs(pi_estimado - pi)
}
puntos_cercanos <- which(abs(errores_absolutos) < 1e-4)
plot(n_sim, pi_estimados, type = "l", col = "#1f77b4", lwd = 2,
xlab = "Tamaño de la Muestra", ylab = "Estimación de π / Error",
main = "Estimación de π y Error Absoluto")
abline(h = pi, col = "#ff7f0e", lty = 2)
lines(n_sim, errores_absolutos, col = "#2ca02c", lwd = 2)
points(n_sim[puntos_cercanos], pi_estimados[puntos_cercanos], col = "#d62728", pch = 19, cex = 1.5)
points(n_sim[puntos_cercanos], errores_absolutos[puntos_cercanos], col = "#d62728", pch = 19, cex = 1.5)
resultados_puntos_cercanos <- data.frame(
Tamano_muestra = n_sim[puntos_cercanos],
Estimacion_pi = pi_estimados[puntos_cercanos],
Error_absoluto = errores_absolutos[puntos_cercanos]
)
kable(resultados_puntos_cercanos, col.names = c("Tamaño de Muestra", "Estimación de π", "Error Absoluto"),
caption = "Puntos donde la Estimación de π Coincide con el Valor Real", digits = 6, align = "c")| Tamaño de Muestra | Estimación de π | Error Absoluto |
|---|---|---|
| 370000 | 3.141611 | 0.000018 |
| 590000 | 3.141634 | 0.000041 |
| 680000 | 3.141665 | 0.000072 |
| 700000 | 3.141526 | 0.000067 |
| 840000 | 3.141676 | 0.000084 |
Conclusión
epsilon <- 0.00001
confidence_level <- 0.99
z_value <- qnorm((1 + confidence_level) / 2)
pi_estimate <- 0.5
n <- (z_value^2 * (pi_estimate * (1 - pi_estimate))) / epsilon^2
n <- ceiling(n)
n_formatted <- formatC(n, format = "f", big.mark = ",", digits = 0)
cat("El número estimado de muestras necesarias es:", n_formatted, "\n")El número estimado de muestras necesarias es: 16,587,241,503
Este número se tiene en cuenta con una nivel de confianza de 99%, y con un error deseado de 0.00001.