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

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.

From → R

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

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

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.

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

4. 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!