Problema 1: Estimación del valor de π

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 cuadr ado 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 π.

Pasos sugeridos:

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

n=1000
x=runif(n,0,1)
#print(x)

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(n,0,1)
#print(y)

De esta forma ya tenemos las variables X & Y, usando la función runif que genera random deviates

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.

# Distancia de cada punto al centro (0.5, 0.5).
distance = numeric()

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

#vector con valor 1 si distance es menor que 0.25.
dots = ifelse(distance < 0.25, 1, 0)

#print(dots)

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

# Creamos una secuencia de valores desde 0 hasta 2π con 100 puntos con el mismo espaciado.
t = seq(0, 2*pi, length.out = 100)

# Gráfico  con relación de aspecto 1 y límites de x & y, entre 0 & 1
plot(0, 0, asp = 1, type = "n", xlim = c(0, 1), ylim = c(0, 1),  ann = F)

# Dibujamos líneas de cuadrícula con intervalos de 0.1
abline(v = seq(0, 1, 0.1), lty = 2, col = "gray")
abline(h = seq(0, 1, 0.1), lty = 2, col = "gray")

# Definimos el centro de x, centro de y, y el radio del círculo.
xcent = 0.5
ycent = 0.5
radius = 0.5

# Calculamos las coordenadas xx e yy
xx <- xcent + cos(t)*radius
yy <- ycent + sin(t)*radius

points(xx, yy, type = "l", col = "green")
points(x,y,pch = 1)

sum(dots)
## [1] 767

A mayor número de n, mayor será la aproximación a π

4 * sum(dots)/length(x) 
## [1] 3.068

Para 10.000:

[1] 3.1564