PONTIFICIA UNIVERSIDAD JAVERIANA DE CALI

Estimación del valor de π

En la figura, un círcuito 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 éste, 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 π.

  1. Se genera n coordenadas 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).
#generar las coordenadas
n1=1000
n2=10000
n3=100000
n4=10000000
x1=runif(n1,0,1)
x2=runif(n2,0,1)
x3=runif(n3,0,1)
x4=runif(n4,0,1)
  1. 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.
n1=1000
n2=10000
n3=100000
n4=10000000
y1=runif(n1,0,1)
y2=runif(n2,0,1)
y3=runif(n3,0,1)
y4=runif(n4,0,1)
  1. 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.
d=function(a){(a[1]-0.5)^2+(a[2]-0.5)^2}
xy1=data.frame(x1,y1)
xy2=data.frame(x2,y2)
xy3=data.frame(x3,y3)
xy4=data.frame(x4,y4)
dist1=apply(xy1,1,d)
dist2=apply(xy2,1,d)
dist3=apply(xy3,1,d)
dist4=apply(xy4,1,d)
cont1=as.numeric(dist1<=0.25)
cont2=as.numeric(dist2<=0.25)
cont3=as.numeric(dist3<=0.25)
cont4=as.numeric(dist4<=0.25)

pi1=4*(sum(cont1)/n1)
pi2=4*(sum(cont2)/n2)
pi3=4*(sum(cont3)/n3)
pi4=4*(sum(cont4)/n4)

error1=4*(sum(cont1)/n1)-pi
error2=4*(sum(cont2)/n2)-pi
error3=4*(sum(cont3)/n3)-pi
error4=4*(sum(cont4)/n4)-pi

cat(" grafica de puntos n= 1000")
 grafica de puntos n= 1000
#graficar los puntos
t <- seq(0, 2*pi1, 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 = "black")
abline(h = seq(0, 1, 0.1), lty = 2, col = "green3")

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 = "red")
points(x1,y1,pch = 20, col="green4")

  1. ¿Cuántos de los puntos están dentro del círculo? ¿Cuál es su estimación de π?
cat(" Estimación de π a distintos tamaños de muestra")
 Estimación de π a distintos tamaños de muestra
tabla <- data.frame(
   tamaño_muetra=c("1000","10000","100000","10000000"),
  pi = c(pi1, pi2, pi3, pi4),
  Error = c(error1, error2, error3, error4)
)

# Visualizar la tabla
print(tabla)
  tamaño_muetra       pi         Error
1          1000 3.108000 -0.0335926536
2         10000 3.170000  0.0284073464
3        100000 3.144680  0.0030873464
4      10000000 3.142413  0.0008201464
#numero de puntos dentro del circulo
cat(" Número de puntos según el tamaño de muestra")
 Número de puntos según el tamaño de muestra
cat(" n= 1000, puntos:",sum(cont1))
 n= 1000, puntos: 777
cat("n= 10000, puntos: ",sum(cont2))
n= 10000, puntos:  7925
cat("n= 100000, puntos: ",sum(cont3))
n= 100000, puntos:  78617
cat("n= 10000000, puntos: ",sum(cont4))
n= 10000000, puntos:  7856032

Conclusión: Se observa que los puntos van aumentando a medida que el tamaño de muestra es más grande y hace que la estimación se acerque al valor real del parámetro.