Enunciado

La siguiente figura sugiere como estimar el valor de π con una simulación. En la figura, un circulo 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 π.

Desarrollo

La solución del problema propuesto gira entorno a una tecnica numerica conocida como “Método de Montecarlo”.

Este método utiliza números aleatorios para estimar resultados en situaciones complejas o probabilísticas. Se basa en la generación de un gran número de simulaciones mediante números aleatorios controlados, aplicados a un modelo matemático del problema en cuestión. Luego, se calculan estadísticas sobre los resultados de estas simulaciones para obtener una estimación del resultado deseado. Este método es ampliamente utilizado en áreas como la física, la estadística, la ingeniería y la toma de decisiones cuando las soluciones analíticas son difíciles de obtener o costosas.

Para la solución del problema se siguen los siguientes pasos:

  1. Se genera el par de coordenadas que se utilizaran en la simulación en un rango de 0 - 1. En este caso el valor de n es una variable arbitraria que refleja la cantidad de datos aleatorios generados mediante una distribución uniforme, los cuales son almacenados en un dataframe para una manipulación más eficaz.
n = 1000
x <- runif(n, 0, 1)
y <- runif(n, 0, 1)

datos <- data.frame(x, y)
head(datos, 5)
##           x        y
## 1 0.0236046 0.256794
## 2 0.1893151 0.899770
## 3 0.7043433 0.551076
## 4 0.6312421 0.977738
## 5 0.3636267 0.795067
  1. Teniendo en cuenta que el circulo circunscrito en el cuadrado tiene radio de 0.5, creamos una función auxiliar que nos permite calcular la distancia Euclideana que hay desde el centro del circulo hasta el punto en cuestión. Esta distancia se almacena en el dataframe “datos” en una nueva columna para procesos posteriores.
dist_func <- function(x, y){
  return(sqrt((x-0.5)**2+(y-0.5)**2))
}

datos$distancia <- apply(datos, 1, function(row) dist_func(row["x"], row["y"]))
head(datos)
##           x        y distancia
## 1 0.0236046 0.256794  0.534885
## 2 0.1893151 0.899770  0.506301
## 3 0.7043433 0.551076  0.210630
## 4 0.6312421 0.977738  0.495437
## 5 0.3636267 0.795067  0.325057
## 6 0.1733114 0.651249  0.360003
  1. Si la distancia calculada es mayor a 0.5 significa que la coordenada en cuestión se encuentra por fuera del circulo, caso contrario si esta es menor o igual a 0.5 la coordenada esta contenida en el circulo. Con base en estas premisas, creamos una nueva columna que nos permite validar dicha condición donde 1 indica que el punto esta en el circulo y 0 en otro caso.
datos$condicion <- as.numeric(datos$distancia < 0.5)
head(datos)
##           x        y distancia condicion
## 1 0.0236046 0.256794  0.534885         0
## 2 0.1893151 0.899770  0.506301         0
## 3 0.7043433 0.551076  0.210630         1
## 4 0.6312421 0.977738  0.495437         1
## 5 0.3636267 0.795067  0.325057         1
## 6 0.1733114 0.651249  0.360003         1
  1. Ahora bien, la relación que hay entre el área del circulo, el área del cuadrado y la cantidad de puntos utilizados en la simulación se muestra a continuación:

\[ \frac{A_o}{A_c} = \frac{n_o}{n} \] Donde,

\(A_o\) es el área del circulo.

\(A_c\) es el área del cuadrado.

\(n_o\) es la cantidad de puntos que hay dentro del circulo.

\(n\) es la cantidad de puntos total usados en la simulación.

Luego, el área del circulo y el cuadrado estan definidas como sigue:

\[ A_o = πr^2 = π(0.5)^2 = \frac{π}{4} \\ A_c = l^2 = 1^2 = 1 \\ \frac{A_o}{A_c} = \frac{π}{4} \] Reemplazando este resultado en la ecuación anterior, podemos establecer una relación entre π y el número de puntos con el cual se realizó la simulación. Dicho resultado establece una ecuación que nos permite estimar el valor de π.

\[ \frac{π}{4} = \frac{n_o}{n} \\ π \approx 4\frac{n_o}{n} \]

  1. La sumatoria de la columna de “Condicion” nos indica el número de datos contenidos en el circulo durante la simulación (\(n_o\)). Asimismo, la cantidad total de datos (\(n\)) se definió de 1000`, por lo cual, el valor estimado de π en este caso es el siguiente:
# Cantidad de datos dentro del circulo
n_o <- sum(datos$condicion)
pi_approx <- 4*n_o/n
pi_approx_string <- sprintf("%.6f", pi_approx)

\[ π \approx 4\frac{n_o}{n} \approx 4\frac{798}{1000} \approx 3.192000 \] 6. El error porcentual obtenido para la estimación es el siguiente:

error_porcentual <- abs((pi-pi_approx)/pi)*100

\[ error = 1.604516 \]

Resultados

A partir de la metodología expuesta anteriormente, se consiguio un encontrar una aproximación de π mediante el método de Montecarlo, con un error de error. Con el objetivo de disminuir el error y encontrar un valor de (\(n\)) óptimo, se ejecuta la metodología antes expuesta para diferentes valores de n, donde en cada caso se realizan 50 simulaciones obteniendo los siguientes resultados:

est_pi <- function(n){
  # n corresponde a la cantidad de datos aleatorios
  
  # x, y corresponden a las variables aleatorias
  x <- runif(n, 0, 1)
  y <- runif(n, 0, 1)
  
  datos <- data.frame(x, y)
  
  # Se calcula la distancia euclideana en relación al centro del circulo
  datos$distancia <- apply(datos, 1, function(row) dist_func(row["x"], row["y"]))
  
  # Se valida si los datos estan dentro del circulo
  datos$condicion <- as.numeric(datos$distancia < 0.5)
  
  # Se calcula una aproximación de pi
  n_o <- sum(datos$condicion)
  pi_approx <- 4*n_o/n
  return(pi_approx)
}

data <- c(1000, 5000, 10000, 50000, 100000)

n <- c()
for (value in data){
  n_values <- rep(value, times = 50)
  n <- append(n, c(n_values))
}

results <- data.frame(n)

results$estimacion <- apply(results, 1, function(row) est_pi(row['n']))
results$error <- abs((pi-results$estimacion)/pi)*100

results$n <- factor(results$n, levels = data, labels = c("1000", "5000", "10000", "50000", "100000"))

# Create a box plot of the grouped data
plot <- ggplot(results, aes(x= n, y = estimacion)) +
geom_boxplot() +  # Add box plots
labs(x = "Valores de n", y = "Estimado", title = "Estimaciones de π con 50 simulaciones") +
geom_hline(yintercept = pi, color = "red", linetype = "dashed")

# Display the plot
print(plot)

A partir del gráfico anterior se puede evidenciar que cuando aumenta la cantidad de variables aleatorias, disminuye la dispersión de los datos, y de esta forma la estimación de π es más precisa. Asimismo, se presenta el resumen de los parametros para cada una de las simulaciones, en este se puede evidenciar que la media genera una buena aproximación a partir de los 5.000 datos. No obstante, teneindo en cuenta la dispersión de los datos, desde los 50.000 datos la diferencia entre el 1er y 3er cuartil se reduce considerablemente. Por lo cual, las simulaciones realizadas a partir de los 50.000 datos generan una buena aproimación de π.

summary_data <- aggregate(estimacion ~ n, data = results, FUN = summary)
summary_data
##        n estimacion.Min. estimacion.1st Qu. estimacion.Median estimacion.Mean
## 1   1000         3.04400            3.10400           3.14200         3.13960
## 2   5000         3.08960            3.12300           3.14200         3.14163
## 3  10000         3.11040            3.13330           3.14360         3.14130
## 4  50000         3.12432            3.13604           3.14004         3.14147
## 5 100000         3.13032            3.13731           3.14158         3.14053
##   estimacion.3rd Qu. estimacion.Max.
## 1            3.17500         3.28800
## 2            3.15620         3.20080
## 3            3.15340         3.16440
## 4            3.14420         3.16688
## 5            3.14287         3.15280

Conclusiones

Se realizaron 50 simulaciones en cada uno de los diferentes conjuntos de datos, generados a partir de variables aleatorias que obedecen a una distribución uniforme, utilizando el método de Montecarlo. Dentro de los hallazgos encontrados se tienen los siguientes datos: