La siguiente figura sugiere como estimar el valor de π con una simulación. En la figura, un circuito 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 puede encontrar una aproximación para el valor de π.
Genere n coordenadas de:
\[x : X1, . . . , Xn\]
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).
## [1] 0.293958078 0.702404311 0.001782669 0.293773957 0.602512982 0.805944872
Genere 1000 coordenadas de:
\[y: Y1,...,Yn \]
utilizando la distribución uniforme con valor mínimo de 0 y valor máximo de 1
## [1] 0.6700729 0.5277298 0.9204307 0.1761803 0.5988496 0.9322926
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)^2+(Yi−0.5)^2\]
que es el cuadrado de la distancia, y al determinar si es menor que 0.25.
set.seed(1)
distancia <- numeric()
for (i in 1:length(x))
{
distancia[i] = (x[i]-0.5)^2 + (y[i]-0.5)^2
}
punto = ifelse(distancia< 0.25, 1, 0) # 1 = adentro, 0 = afuera## [1] 0.07137806 0.04173645 0.42498249 0.14738837 0.02028016 0.28047916
## [1] 1 1 0 1 1 0
#graficamos los puntos
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)
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 = "darkred")
points(x,y,pch = 20)(Nota: 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
## [1] 3.184
## [1] -0.04240735
set.seed(1)
c = runif(10000,0,1)
d = runif(10000,0,1)
distancia_2 <- numeric()
for (i in 1:length(c))
{
distancia_2[i] = (c[i]-0.5)^2 + (d[i]-0.5)^2
}
punto_2 = ifelse(distancia_2< 0.25, 1, 0)
estimacion_pi_2 <- ((4*sum(punto_2))/length(c))
(estimacion_pi_2)## [1] 3.1224
## [1] 0.01919265
set.seed(1)
n = runif(100000,0,1)
m = runif(100000,0,1)
distancia_3 <- numeric()
for (i in 1:length(n))
{
distancia_3[i] = (n[i]-0.5)^2 + (m[i]-0.5)^2
}
punto_3 = ifelse(distancia_3< 0.25, 1, 0)
estimacion_pi_3 <- ((4*sum(punto_3))/length(n))
(estimacion_pi_3)## [1] 3.13648
## [1] 0.005112654
## [1] 3.18400 3.12240 3.13648
barplot(vector_pi,
main="Comparativa de estimaciones de π",
xlab = "estimaciones de π ",
ylim = c(0,5.5),
col = c("lightyellow", "lightblue", "darkred"),
ylab = "Valores para π ",
legend.text=c("π_mil","π_10mil","π_100mil"))
abline (h = 3.14, col = 'blue', lwd = 3)La gráfica anterior muestra como a medida que aumentamos los puntos tiene mayores probabilidades de dar como resultado una estimación muy cercana al valor real de π, como lo muestra la línea azul que contiene el valor 3.1416
\[ π = 3.1416 \]