Pi Estimator in R

Question

The code below estimates pi in R, now I am trying to find the minimum number of terms N_Min you would have to include in your estimate of pie to make it accurate to three decimal places.

pi_Est<- function(NTerms){
  NTerms = 5 # start with an estimate of just five terms
  pi_Est = 0 # initialise the value of pi to zero
  Sum_i = NA # initialise the summation variable to null
  for(ii in 1:NTerms)
  {
    Sum_i[ii] = (-1)^(ii+1)/(2*ii - 1)  # this is the series equation for calculating pi
  }
  Sum_i = 4*Sum_i # multiply by four as required in the formula (see lecture notes)
  
  pi_Est = sum(Sum_i)
  cat('\nThe estimate of pi with terms = ', NTerms ,' is ',pi_Est)
  
}

Answer 1

First of all, I would change some things about your function. Instead of getting it to print out a message, get it to return a value. Otherwise it becomes very difficult to do anything with its output, including testing it for convergence to pi.

Also, no matter what the value of NTerms is you feed this function, you are immediately over-writing NTerms inside the function.

You could rewrite the function like this:

pi_Est <- function(NTerms) {
  
  pi_Est <- 0 
  Sum_i  <- numeric()
  
  for(ii in seq(NTerms))
  {
    Sum_i[ii] <- (-1)^(ii+1)/(2*ii - 1) 
  }
  
  return(sum(4 * Sum_i))
}

And to show it converges to pi, let's test it with 50,000 terms:

pi_Est(50000)
#> [1] 3.141573

Now, if we want to find the first value of NTerms that is correct to 3 decimal places, we are going to need to be able to call this function on a vector of NTerms - at the moment it is only working on a single number. So let's define the function f that vectorizes pi_Est:

f <- Vectorize(pi_Est)

Now, let's create the estimate for all values of NTerms between 1 and 2,000 and store them in a vector:

estimates <- f(1:2000)

We can see that the values of estimates seem to oscillate round and converge to pi if we plot the first 100 values:

plot(estimates[1:100], type = 'l')
abline(h = pi)

Our answer is just the first value which, when rounded to three decimal places, is the same as pi rounded to three decimal places:

result <- which(round(estimates, 3) == round(pi, 3))[1]

result
#> [1] 1103

And we can check this is correct by feeding 1103 into our original function:

pi_Est(result)
#> [1] 3.142499

You will see that this gives us 3.142, which is the same as pi rounded to 3 decimal places.

^{Created on 2022-01-31 by the reprex package (v2.0.1)}

Answer 2

1000 terms are required to make the estimate accurate to within 0.001:

pi_Est1 <- function(n) {
  if (n == 0) return(0)
  neg <- 1/seq(3, 2*n + 1, 4)
  if (n%%2) neg[length(neg)] <- 0
  4*sum(1/seq(1, 2*n, 4) - neg)
}

pi_Est2 <- function(tol) {
  for (i in ceiling(1/tol + 0.5):0) {
    est <- pi_Est1(i)
    if (abs(est - pi) > tol) break
    est1 <- est
  }
  
  list(NTerms = i + 1, Estimate = est1)
}

tol <- 1e-3
pi_Est2(tol)
#> $NTerms
#> [1] 1000
#> 
#> $Estimate
#> [1] 3.140593
tol - abs(pi - pi_Est2(tol)$Estimate)
#> [1] 2.500001e-10
tol - abs(pi - pi_Est1(pi_Est2(tol)$NTerms - 1))
#> [1] -1.00075e-06

Created on 2022-01-31 by the reprex package (v2.0.1)

Answer 3

Perhaps we can try the code below

pi_Est <- function(digits = 3) {
  s <- 0
  ii <- 1
  repeat  {
    s <- s + 4 * (-1)^(ii + 1) / (2 * ii - 1)
    if (round(s, digits) == round(pi, digits)) break
    ii <- ii + 1
  }
  list(est = s, iter = ii)
}

and you will see

> pi_Est()
$est
[1] 3.142499

$iter
[1] 1103


> pi_Est(5)
$est
[1] 3.141585

$iter
[1] 130658