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

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!