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

Para la estimación de \(\pi\) usando el método de Monte Carlo. La idea es generar una gran cantidad de puntos aleatorios dentro de un cuadrado que circunscribe a un círculo y luego contar cuántos de esos puntos caen dentro del círculo. Finalmente, calculando la proporción entre los puntos dentro y fuera del círculo, se puede obtener una aproximación del valor de \(\pi\).

En este caso, vamos a generar un cuadrado de lado 1, usando puntos aleatorios de una distribución uniforme con valor mínimo 0 y valor máximo 1. En el caso del círculo, recordemos que la ecuación canónica de la circunferencia con centro \(c(h, k)\) y radio r está dada por:

\[ (x-h)^2 + (y - k)^2 = r^2 \] En nuestro caso, la circunferencia tiene centro en el punto \(c \left( \frac{1}{2}, \frac{1}{2} \right)\) y radio \(r= 0.5\). Luego, su ecuación canónica será:

\[ \left( x- \frac{1}{2} \right)^2 + \left(y - \frac{1}{2} \right)^2 = \left(\frac{1}{2} \right)^2 \]

Una vez definido el espacio donde vamos a trabajar, procederemos creando una función teniendo en cuenta las siguientes consideraciones:

  • la función tendrá como entrada el número de iteraciones y el valor del parámetro \(\theta\).

  • En esta función, se deben generar dos muestras aleatorias con distribución uniforme tomando un valor mínimo de 0 y valor máximo de 1.

  • Queremos que nos clasifique los puntos de acuerdo a si están o no dentro de la circunferencia. Para esto, se debe tener en cuenta que un punto (x,y) está en la circunferencia si se cumple que

\[ \left( x- \frac{1}{2} \right)^2 + \left(y - \frac{1}{2} \right)^2 \leq \left(\frac{1}{2} \right)^2 \].

  • Finalmente, necesitamos que toda la información se almacene en un data frame, para poder visualizar los resultados obtenidos y entender mejor el proceso de convergencia.
Teniendo en cuenta dichas consideraciones, procederemos a crear la función.

Función que retorna el valor de \(\pi\)

La siguienete función retorna el valor estimado de pi, de acuerdo con el número de iteraciones.

pi_estimacion1 <- function(n){
  
  set.seed(1645) # semilla para que los datos no cambien
  
  x <- runif(n, 
            min = 0,     # Coordenadas en x
            max = 1)

  y <- runif(n, 
            min = 0,    # coordenadas en y
            max = 1)
  

  puntos_circ1 <- which(sqrt((x-0.5)^2+(y- 0.5)^2) <= 0.5) 
  
  pi_estimado1 <- 4*length(puntos_circ1)/n
  
  return(pi_estimado1)
}

Veamos algunas estimaciones usando: 10, 100, 1000, 10000, 1000000 y 10000000 puntos aleatorios.

n Estimación de \(\pi\)
10 3.2
100 3.04
1000 3.148
10000 3.176
100000 3.13284
1000000 3.137124
10000000 3.140863

Para una muestra de tamaño 1000, tenemos que 775 puntos se encuentran dentro de la circunferencia, con lo cual se obtiene una aproximaxion de \(\pi= 3.148\).

Note que la estimación mejora cada vez que incrementamos el número de puntos aleatorios. Por ejemplo, usando 10000000 nos muestra una aproximación con 3 cifras significativas.

Función que retorna el valor de \(\pi\) en cada iteración

La anterior función cumple su proposito a la perfección. Sin embargo, sería interesante ver como mejora o empeora la aproximación de acuerdo con el número de puntos seleccionado. A continuación, se muestra otra alternativa para la estimación de \(\pi\), la cual retorna un data frame con las estimaciones para cada iteración. 
pi_estimacion2 <- function(iteraciones){
  set.seed(1645)
  
  x <- runif(iteraciones, # Coordenadas en x
            min = 0,
            max = 1)

  y <- runif(iteraciones, # coordenadas en y
            min = 0,
            max = 1)
  
  df <- data.frame(x, y) # Data frame con los puntos aleatorios uniformes (x, y)
  
  df$iteraciones <- 1:iteraciones # Añadir columna con el número de iteraciones
  
  df$puntos_circ <- ifelse( sqrt((x-0.5)^2+(y- 0.5)^2) <= 0.5, 1, 0)  #Añadir columnas para identificar puntos en el círculo
  
  df$pi_estimado <- 4*cumsum(df$puntos_circ)/ df$iteraciones # Cálcular el valor de pi
  
  return(df)
}

Esta función, nos permite ver las estimaciones de \(\pi\) de acuerdo al número de iteraciones. A continuación, haremos una visualización de dicho comportamiento.

pi_valores <- pi_estimacion2(5000)

ggplot(pi_valores, aes(x = iteraciones, y = pi_estimado))+
  geom_line(col = "blue")+
  geom_hline(yintercept = pi, col = "red")+
  labs(title = expression(paste("Aproximación de ", pi)),
       x = "Número de puntos",
       y = "Valor estimado")

En este caso podemos observar que las primeras aproximaciones muestran valores cercanos o iguales a 4. Sin embargo, cuando hacemos n lo sufcientemente grande, es decir \(n \to \infty\), la precisión de las estimaciones mejora cada vez más.

En la última gráfica se observa que a medida que aumenta el número de puntos la aproximación mejora. Sin embargo, se siguen observando picos y valles en la gráfica. Esto es debido a que los métodos de Monte Carlo son un proceso estocástico numerico, es decir, una secuencia de estados cuya evolucion viene determinada por los sucesos aleatorios, los cuales pueden ser irregulares en muchos casos.