Ejercicio

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

JuveYell

Desarrollo

Coordenadas n=1000

set.seed(123) # Establece una semilla para la reproducibilidad
n=1000
x <- runif(n, min = 0, max = 1)
y <- runif(n, min = 0, max = 1)

DistanciaCentro <- (x-0.5)^2+(y-0.5)^2

# Crear el plot base sin ejes
plot(x, y, col="blue", xlim=c(0,1), ylim=c(0,1), xlab="X", ylab="Y", asp=1, axes=FALSE)

# Añadir manualmente los ejes
axis(1, at=seq(0, 1, by=0.5), pos=0)  # Eje X
axis(2, at=seq(0, 1, by=0.5), pos=0)  # Eje Y


# Dibujar una circunferencia de radio 0.5 alrededor del punto (1,1)
draw.circle(0.5,0.5 , radius=0.5, border="red")

# Valores menores que 0.25
puntosenCirculo <- DistanciaCentro <= 0.25

# Suma los valores TRUE (1) para contar cuántos son menores que 0.25
TotalpuntosenCirculo <- sum(puntosenCirculo)

# Imprime el resultado
cat("Con n=1000, el total de puntos en el círculo es:", TotalpuntosenCirculo)
## Con n=1000, el total de puntos en el círculo es: 800
estimadopi = (4*TotalpuntosenCirculo)/n
cat("Con n=1000, el estimador de π es:", estimadopi)
## Con n=1000, el estimador de π es: 3.2
cat("Con n=1000, Error Absoluto (πEstimado - πReal) es:", estimadopi - pi)
## Con n=1000, Error Absoluto (πEstimado - πReal) es: 0.05840735
cat("Con n=1000, Error Relativo (πEstimado - πReal)/(πReal) es:", ((estimadopi - pi)/pi))
## Con n=1000, Error Relativo (πEstimado - πReal)/(πReal) es: 0.01859164

Coordenadas n=10000

set.seed(123) # Establece una semilla para la reproducibilidad
n=10000
x <- runif(n, min = 0, max = 1)
y <- runif(n, min = 0, max = 1)

DistanciaCentro <- (x-0.5)^2+(y-0.5)^2

# Crear el plot base sin ejes
plot(x, y, col="blue", xlim=c(0,1), ylim=c(0,1), xlab="X", ylab="Y", asp=1, axes=FALSE)

# Añadir manualmente los ejes
axis(1, at=seq(0, 1, by=0.5), pos=0)  # Eje X
axis(2, at=seq(0, 1, by=0.5), pos=0)  # Eje Y


# Dibujar una circunferencia de radio 0.5 alrededor del punto (1,1)
draw.circle(0.5,0.5 , radius=0.5, border="red")

# Valores menores que 0.25
puntosenCirculo <- DistanciaCentro <= 0.25

# Suma los valores TRUE (1) para contar cuántos son menores que 0.25
TotalpuntosenCirculo <- sum(puntosenCirculo)

# Imprime el resultado
cat("Con n=10000, el total de puntos en el círculo es:", TotalpuntosenCirculo)
## Con n=10000, el total de puntos en el círculo es: 7894
estimadopi = (4*TotalpuntosenCirculo)/n
cat("Con n=10000, el estimador de π es:", estimadopi)
## Con n=10000, el estimador de π es: 3.1576
cat("Con n=10000, Error Absoluto (πEstimado - πReal) es:", estimadopi - pi)
## Con n=10000, Error Absoluto (πEstimado - πReal) es: 0.01600735
cat("Con n=10000, Error Relativo (πEstimado - πReal)/(πReal) es:", ((estimadopi - pi)/pi))
## Con n=10000, Error Relativo (πEstimado - πReal)/(πReal) es: 0.005095297

Coordenadas n=100000

set.seed(123) # Establece una semilla para la reproducibilidad
n=100000
x <- runif(n, min = 0, max = 1)
y <- runif(n, min = 0, max = 1)

DistanciaCentro <- (x-0.5)^2+(y-0.5)^2

# Crear el plot base sin ejes
plot(x, y, col="blue", xlim=c(0,1), ylim=c(0,1), xlab="X", ylab="Y", asp=1, axes=FALSE)

# Añadir manualmente los ejes
axis(1, at=seq(0, 1, by=0.5), pos=0)  # Eje X
axis(2, at=seq(0, 1, by=0.5), pos=0)  # Eje Y


# Dibujar una circunferencia de radio 0.5 alrededor del punto (1,1)
draw.circle(0.5,0.5 , radius=0.5, border="red")

# Valores menores que 0.25
puntosenCirculo <- DistanciaCentro <= 0.25

# Suma los valores TRUE (1) para contar cuántos son menores que 0.25
TotalpuntosenCirculo <- sum(puntosenCirculo)

# Imprime el resultado
cat("Con n=100000, el total de puntos en el círculo es:", TotalpuntosenCirculo)
## Con n=100000, el total de puntos en el círculo es: 78658
estimadopi = TotalpuntosenCirculo/n*4
cat("Con n=100000, el estimador de π es:", estimadopi)
## Con n=100000, el estimador de π es: 3.14632
cat("Con n=100000, Error Absoluto (πEstimado - πReal) es:", (estimadopi - pi))
## Con n=100000, Error Absoluto (πEstimado - πReal) es: 0.004727346
cat("Con n=100000, Error Relativo (πEstimado - πReal)/(πReal) es:", ((estimadopi - pi)/pi))
## Con n=100000, Error Relativo (πEstimado - πReal)/(πReal) es: 0.001504761

Conclusión

A lo largo de este experimento, hemos implementado un método de estimación de π utilizando un enfoque probabilístico, específicamente, el método de Montecarlo. Generamos puntos aleatorios dentro de un cuadrado unitario y determinamos la proporción de estos puntos que caen dentro de un círculo inscrito de radio 0.5. A medida que incrementamos el número de puntos generados aleatoriamente, nuestra estimación de π se ha vuelto más precisa, acercándose al valor real de π.