Estimating Pi
I came across this post which gives a method to estimate Pi by using a circle, it’s circumscribed square and (lots of) random points within said square. Booth used Stata to estimate Pi, but here’s some R code to do the same thing…
x <- 0.5 # center x y <- 0.5 # center y n <- 1000 # nr of pts r <- 0.5 # radius pts <- seq(0, 2 * pi, length.out = n) plot(sin(pts), cos(pts), type = 'l', asp = 1) # test require(sp) xy <- cbind(x + r * sin(pts), y + r * cos(pts)) sl <- SpatialPolygons(list(Polygons(list(Polygon(xy)), "polygon"))) plot(sl, add=FALSE, col = 'red', axes=T ) # the square xy <- cbind(c(0, 1, 1, 0), c(0, 0, 1, 1)) sq <- SpatialPolygons(list(Polygons(list(Polygon(xy)), "polygon"))) plot(sq, add = TRUE) N <- 1e6 x <- runif(N, 0, 1) y <- runif(N, 0, 1) sp <- SpatialPoints(cbind(x, y)) plot(sp, add = TRUE, col = "green") require(rgeos) (sim_pi <- (sum(gIntersects(sp, sl, byid = TRUE))/N) *4) sim_pi - pi
Note the use of sp and rgeos packages to calculate the intersections.
Insights of a PhD 
I remember doing this in Fortran nearly 50 years ago. Thanks for posting the R version.
There’s a second method in Francesco’s comment too… should be quite a bit faster, particularly for large Ns, as it doesn’t go down the spatial route
Hi, here a version I wrote some month ago which applies the same concept but is far simpler in term of coding.
Imagine a circle centered in (0,0) and with radius = 1: the distance of each of the random points from the center tells you if the point is inside the circle or not.
# load required libraries
library(dplyr)
library(ggplot2)
library(purrr)
# number of points of the experiment
n <- 1e6
# experiment data frame
experiment %
mutate(X = runif(n, -1, 1), # random number X axis
Y = runif(n, -1, 1), # random number Y axis
D = sqrt(X^2 + Y^2), # distance from center
circle = if_else(D <= 1, TRUE, FALSE) # if distance is %
summarize(pi_approx = sum(circle) / n() * 4,
pi_diff = pi_approx – pi,
pi_err = (pi_approx – pi) / pi * 100)
# ————————
# check how the simulation converge to the exact solution as the size of the experiment increases
# summarize the process in a function
pi_approximator <- function(n = 10) {
experiment %
mutate(X = runif(n, -1, 1),
Y = runif(n, -1, 1),
D = sqrt(X^2 + Y^2),
circle = if_else(D %
summarize(N = n(),
pi_approx = sum(circle) / n() * 4,
pi_diff = pi_approx – pi,
pi_err = (pi_approx – pi) / pi * 100)
}
# map the results for different sizes of the experiment
output %
ggplot(aes(x = repetitions, y = approx_error)) +
geom_line() +
geom_hline(aes(yintercept = 0), color = “red”, linetype = 2) +
scale_x_log10()
Regards,
Francesco
Thanks for that… it is indeed a nice solution. Using the spatial stuff can get a bit slow, particularly when n gets large.
Spatial polygons??? Talk about over engineered.
geom_line()?? Yech.
# Set n
n <- 10000
pi_estimate <- 4*sum(runif(n,-1,+1)^2 + runif(n,-1,+1)^2 < 1)/n
I never said it was the easiest nor only way to do it (there’s always a multitude of methods in R). I’ve been doing some spatial stuff recently which is probably why that came to mind first. And plotting is a good way to get the idea across too.
Your post has given me the idea to reinterpret quadrature of the circle with gganimate and R code. It’s a geometric and visual construction, not random, but I find it’s good to see what is Pi in an more pragmatic way.
it intuits infinity, segments become more and more tiny as long as n becomes larger, then we can see Pi as a perfection.
https://github.com/GuillaumePressiat/under_pi
Interesting…thanks!