Estimación del valor de \(\pi\)

Se requiere estimar el valor de \(\pi\) con una simulación. Para ello se hará uso del método de Monte Carlo, que consiste en dibujar un cuadrado de área igual a 1 y dibujar inscrito en el cuadrado un círculo con un área igual a \(\pi /4\), Luego se dibujan puntos de manera aleatoria sobre la superficie dibujada. Los puntos que están fuera del círculo y los que están dentro, sirven como estimadores de las áreas internas y externas del círculo. como se muestra en la figura.

El procedimiento es elegir de forma aleatoria n puntos dentro del cuadrado. La probabilidad de que un punto esté dentro del circulo es igual a la fracción del área del cuadrado que abarca a este, la cual es \(\pi/4\). Por tanto, se puede estimar el valor de \(\pi/4\) al contar el número de puntos dentro del círculo, para obtener la estimación de \(\pi /4\). De este último resultado se puede encontrar la estimación para el valor de \(\pi\).

t <- proc.time()
# Función para estimar pi mediante el método de Monte Carlo
est_pi <- function(num_puntos) {
  
  x1 = runif(num_puntos, 0, 1)
  y1 = runif(num_puntos, 0, 1)
  
  d = (x1 - 0.5)^2 +(y1 - 0.5)^2
  y = as.numeric(d <= 0.25)
  
  pi_est = sum(y)/num_puntos*4
  
  return(pi_est)
}

# Valores de puntos para la simulación
num_ptos_simulacion <- c(1000, 10000, 100000, 1000000, 10000000, 100000000, 1000000000)

# Crear una lista para almacenar los resultados
result <- list()

# Realizar la simulación para cada valor de puntos
for (num_ptos in num_ptos_simulacion) {
  # Estimar pi
  pi_est <- est_pi(num_ptos)
  
  # Valor real de pi
  pi_real <- pi
  
  # Diferencia entre el estimado y el valor real de pi
  dif_pi <- abs(pi_est - pi_real)
  
  # Almacenar los resultados en la lista
  result[[as.character(num_ptos)]] <- c(num_ptos, pi_est, dif_pi)
}

# Crear la tabla resumen
res <- do.call(rbind, result)
colnames(res) <- c("Núm. Puntos", "Estimado de Pi", "Diferencia con Pi Real")

# Mostrar la tabla resumen
print(res)
##       Núm. Puntos Estimado de Pi Diferencia con Pi Real
## 1000        1e+03       3.132000           9.592654e-03
## 10000       1e+04       3.128000           1.359265e-02
## 1e+05       1e+05       3.132680           8.912654e-03
## 1e+06       1e+06       3.140224           1.368654e-03
## 1e+07       1e+07       3.141767           1.741464e-04
## 1e+08       1e+08       3.141657           6.458641e-05
## 1e+09       1e+09       3.141636           4.335441e-05
# Calcular tiempo de ejecución
proc.time() - t
##    user  system elapsed 
##  66.082  12.265  78.349

Como se aprecia en el resultado, la estimación de pi con dos cifras decimales se encuentra a partir de los 1e+05 puntos.

Conclusiones

En el caso de la estimación de \(\pi\), se usan valores aleatorios con distribución uniforme para aproximar un valor determinístico.

El principal desafío al utilizar este tipo de simulaciones es su dependencia de los valores de entrada y de la distribución, de ahí que se deba entender muy bien el problema antes de escoger la distribución. En el caso de la estimación de pi, es claro que se requiere la distribución uniforme dado que buscamos que cada variable tenga la misma probabilidad, la probabilidad que el “disparo” caiga dentro del círculo o fuera de el. Otro aspecto a tener en cuenta cuando se utilizan este tipo de simulaciones es la necesidad de potencia computacional robusta, ya que tiende a requerir mucho tiempo en una sola computadora y a mayor número de puntos de simulación mayor precisión en la estimación.

Eficiencia del algoritmo

Se ha dicho que en las simulaciones Monte Carlo, el requerimiento de potencia computacional es un punto en contra de este método, sin embargo, la eficiencia de R también afecta el tiempo de simulación, esto debido al tipo de algoritmo, las estructuras de control y flujo utilizadas en su implementación. Como validación se incluye esta segunda versión del código basado en ciclos for, y se comparará con la versión presentada inicialmente.

t <- proc.time()
# Función para estimar pi mediante el método de Monte Carlo
estimar_pi <- function(num_puntos) {
  puntos_dentro_circulo <- 0
  
  for (i in 1:num_puntos) {
    x <- runif(1, -0.5, 0.5)  # Coordenada x aleatoria dentro del cuadrado
    y <- runif(1, -0.5, 0.5)  # Coordenada y aleatoria dentro del cuadrado
    distancia_al_origen <- x^2 + y^2  # Distancia al origen (centro del cuadrado)
    
    if (distancia_al_origen <= 0.25) {  # Si la distancia está dentro del círculo (radio = 0.5)
      puntos_dentro_circulo <- puntos_dentro_circulo + 1
    }
  }
  
  # Estimación de pi
  pi_estimado <- 4 * puntos_dentro_circulo / num_puntos
  return(pi_estimado)
}

# Valores de puntos para la simulación
num_puntos_simulacion <- c(1000, 10000, 100000, 1000000, 10000000, 100000000)

# Crear una lista para almacenar los resultados
resultados <- list()

# Realizar la simulación para cada valor de puntos
for (num_puntos in num_puntos_simulacion) {
  # Estimar pi
  pi_estimado <- estimar_pi(num_puntos)
  
  # Valor real de pi
  pi_real <- pi
  
  # Diferencia entre el estimado y el valor real de pi
  diferencia_pi <- abs(pi_estimado - pi_real)
  
  # Almacenar los resultados en la lista
  resultados[[as.character(num_puntos)]] <- c(num_puntos, pi_estimado, diferencia_pi)
}

# Crear la tabla resumen
resumen <- do.call(rbind, resultados)
colnames(resumen) <- c("Número de Puntos", "Estimado de Pi", "Diferencia con Pi Real")

# Mostrar la tabla resumen
print(resumen)
##       Número de Puntos Estimado de Pi Diferencia con Pi Real
## 1000             1e+03       3.228000           8.640735e-02
## 10000            1e+04       3.145200           3.607346e-03
## 1e+05            1e+05       3.143160           1.567346e-03
## 1e+06            1e+06       3.140256           1.336654e-03
## 1e+07            1e+07       3.141588           5.053590e-06
## 1e+08            1e+08       3.141305           2.877336e-04
proc.time() - t
##    user  system elapsed 
## 443.226   1.443 444.658

Analizando el desempeño de los dos algoritmos se encuentra que la versión 1 es más intensiva en consumo de RAM llegando hasta los 50GB durante la simulación, mientras que el segundo algoritmo si bien es más intensivo en uso de CPU no consume más allá de los 16GB, sin embargo le toma casi seis veces el tiempo que le toma al algoritmo 1 completar simulaciones hasta para 1e+08 puntos.

Se concluye también que son más eficientes las implementaciones que hacen mayor uso de RAM.

Algoritmo 1 - mayor consumo de RAM (hasta 60GB)

Algoritmo 2 - menor consumo de RAM (hasta 16GB), tarda más de 5 veces lo que tarda el algorimo 1.