La estimación del valor de π es una de las aplicaciones de la simulación de Montecarlo, para esto, se toma un cuadrado entre la coordenadas (1,1) (1, -1) (-1,1) (-1,1) y al interior de este cuadrado se encuentra una circcunferencia de radio 1. Tomando el área del círculo igual a A = πr² y el del cuadrado como A= 2r*2r=4r²; si se establece una relación entre estás dos áreas se obtiene que el área del círculo sobre el área del cuadrado que equivale a π/4.
Por lo anterior, se sugiere estimar el valor de π con una simulación. Partiendo del hecho que se encuentra un circulo con un área de π/4, está inscrito en un cuadrado cuya área es igual a 1. Se elige de forma aleatoria 100 puntos dentro del cuadrado . La probabilidad de que un punto esté dentro del círculo es igual a la pracció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, que es 79 para obtener la estimación de π/4≈0.76. De este último resultado se debe encontrar una aproximación para el valor de π.
En primer lugar, se generan n coordenadas x:X1, . . . , Xn.Usando 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)
library(STAT)
x <- runif(1000, 0, 1)
plot(1:1000, x)
hist(x)
En segundo lugar, se genera 1000 coordenadas y : Y1,…,Yn, utilizando nuevamente la distribución uniforme con valor mínimo de 0 y valor máximo de 1.
library(STAT)
y <- runif(1000, 0, 1)
plot(1:1000, y)
hist(y)
En tercer lugar, Cada punto (Xi,Yi) se encuentra dentro del círculo si su distancia desde el centro (0.5,0.5) es menor a 0.5. Para cada par (Xi,Yi) determine si la distancia desde el centro es menor a 0.5. Esto último se puede realizar al calcular el valor (Xi−0.5)²+(Yi−0.5)², que es el cuadrado de la distancia, y al determinar si es menor que 0.25
#En este vector se almacenan las distancias
distancias <- numeric()
#Se calculan las distancias
for (i in 1:length(x)){
distancias[i] = (x[i]-0.5)^2 + (y[i]-0.5)^2
}
#Los puntos se guardan tanto dentro como fuera del circulo
#1 -> esta dentro
#0 -> esta fuera
puntos = ifelse(distancias < 0.25, 1, 0)
Finalmente, ¿Cuántos de los puntos están dentro del círculo? ¿Cuál es su estimación de π?
t <- seq(0, 2*pi, length.out = 100)
plot(0, 0, asp = 1, type = "n",
xlim = c(0, 1), ylim = c(0, 1),
ann = F)
abline(v = seq(0, 1, 0.1), lty = 2, col = "gray50")
abline(h = seq(0, 1, 0.1), lty = 2, col = "gray50")
radio <- 1/2
a <- 0.5 ## origen circunferencia eje x
b <- 0.5 ## origen circunferencia eje y
xx <- a + cos(t)*radio
yy <- b + sin(t)*radio
points(xx, yy, type = "l", col = "purple")
points(x,y,pch = 20)
Con sólo 1000 puntos, es probable que su estimación sea inferior por 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
Ahora, se estiman el número de puntos dentro del círculo para esto sumamos cuantos valores son igual a 1
sum(puntos)
## [1] 783
Se procede a estimar el valor de π
4 * sum(puntos)/length(x)
## [1] 3.132
Este valor es muy cercano al valor real de π solo se alejo en 0.0016