R 中的 Pi 估计器答案

【问题标题】：Pi Estimator in RR 中的 Pi 估计器
【发布时间】：2023-02-08 17:48:38
【问题描述】：

下面的代码估计 R 中的 pi，现在我试图找到最少的项数 N_Min 你必须在你对馅饼的估计中包括在内，以使其精确到小数点后三位。

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

【问题讨论】：

如果你在函数内部设置NTerms = 5，输入参数将被覆盖，你将始终得到NTerms = 5的结果。也许将其定义为默认值：pi_Est <- function(NTerms = 5){...}

标签： r

【解决方案1】：

首先，我会改变一些关于你的功能的事情。与其让它打印出一条消息，不如让它返回一个值。否则，很难对其输出做任何事情，包括测试它是否收敛到 pi。

此外，无论您将 NTerms 的值提供给此函数，您都会立即在函数内部覆盖 NTerms。

您可以像这样重写函数：

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

为了证明它收敛于 pi，让我们用 50,000 项来测试它：

pi_Est(50000)
#> [1] 3.141573

现在，如果我们想找到 NTerms 的第一个精确到小数点后 3 位的值，我们将需要能够在向量NTerms - 目前它只处理一个号码。那么让我们定义向量化pi_Est的函数f：

f <- Vectorize(pi_Est)

现在，让我们为 NTerms 在 1 到 2,000 之间的所有值创建估计并将它们存储在一个向量中：

estimates <- f(1:2000)

如果我们绘制前 100 个值，我们可以看到 estimates 的值似乎在振荡并收敛到 pi：

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

我们的答案只是第一个值，当四舍五入到小数点后三位时，它与四舍五入到小数点后三位的圆周率相同：

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

result
#> [1] 1103

我们可以通过将 1103 输入到我们的原始函数中来检查它是否正确：

pi_Est(result)
#> [1] 3.142499

你会看到这给了我们 3.142，这与四舍五入到小数点后三位的圆周率相同。

^{由 reprex package (v2.0.1) 创建于 2022-01-31}

【讨论】：

我觉得可以使用 cumsum 代替 Vectorize 调用。

【解决方案2】：

需要1000条款才能使估计准确到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)

【讨论】：

【解决方案3】：

也许我们可以试试下面的代码

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

你会看到

> pi_Est()
$est
[1] 3.142499

$iter
[1] 1103


> pi_Est(5)
$est
[1] 3.141585

$iter
[1] 130658

【讨论】：

【解决方案4】：

为什么不用一行代码来计算呢？

   Pi <- tail(cumsum(4*(1/seq(1,4*50000000,2))*rep(c(1,-1), 50000000)),1)

【讨论】：