Skip to content

Estimating Pi

October 16, 2018

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.

From → R

8 Comments
  1. Bruce J. Cochrane permalink

    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. Francesco Giorgetti permalink

    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.

  3. Gareth Owen permalink

    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

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

%d bloggers like this: