Estimación del valor de Pi

Problema 1

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 encontrar una aproximación para el valor de π.

Punto a

Genere 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).

z<-1000
x= runif(z,0,1)

Punto b

Genere 1000 coordenadas y : Y1,…,Yn, utilizando nuevamente la distribución uniforme con valor mínimo de 0 y valor máximo de 1

y= runif(1000,0,1)

Punto c

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

dis_p_c <- numeric()

for (i in 1:length(x)){
  dis_p_c[i] = (x[i]-0.5)^2 + (y[i]-0.5)^2
}

puntos = ifelse(dis_p_c < 0.25, 1, 0)

Punto d

¿Cuántos de los puntos están dentro del círculo? ¿Cuál es su estimación de π?

#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)

abline(v = seq(0, 1, 0.1), lty = 2, col = "gray")
abline(h = seq(0, 1, 0.1), lty = 2, col = "gray")

radio <- 1/2

a <- 0.5
b <- 0.5 

xx <- a + cos(t)*radio
yy <- b + sin(t)*radio

points(xx, yy, type = "l", col = "red")
points(x,y,pch = 20)

Número de puntos dentro del círculo

sum(puntos)

## [1] 793

Estimación de pi

4 * sum(puntos)/length(x)

## [1] 3.172

En 1000 puntos el error absoluto es de:

4 * sum(puntos)/length(x) - pi

## [1] 0.03040735

Simulando 10000 puntos

x=runif(10000,0,1)
y= runif(10000,0,1)

dis_p_c <- numeric()

for (i in 1:length(x)){
  dis_p_c[i] = (x[i]-0.5)^2 + (y[i]-0.5)^2
}

puntos = ifelse(dis_p_c < 0.25, 1, 0)

Número de puntos con simulación de 10000

sum(puntos)

## [1] 7839

Estimación de Pi con 10000 puntos

4 * sum(puntos)/length(x)

## [1] 3.1356

Con un error abosluto de:

4 * sum(puntos)/length(x) - pi

## [1] -0.005992654

Simulando 100000 puntos

x=runif(100000,0,1)
y= runif(100000,0,1)

dis_p_c <- numeric()

for (i in 1:length(x)){
  dis_p_c[i] = (x[i]-0.5)^2 + (y[i]-0.5)^2
}

puntos = ifelse(dis_p_c < 0.25, 1, 0)

Número de puntos con simulación de 100000

sum(puntos)

## [1] 78676

Estimación de Pi con 100000 puntos

4 * sum(puntos)/length(x)

## [1] 3.14704