本文参考自 R 语言求圆周率：基于蒙特卡洛方法。主要区别是代码用 base R 实现。

我们先画一个示例图，如下图所示。图中包含一个半径为 1 的单位圆，外接一个正方形。

# 图形参数设置
opar = par(no.readonly = TRUE)
on.exit(par(opar), add = TRUE)
par(pty = "s")  # 保持 1:1 比例

# 绘制单位圆
theta = seq(0, 2 * pi, length.out = 100)
plot(
  x = cos(theta), 
  y = sin(theta), 
  col = "red", 
  lwd = 2, 
  type = "l", 
  bty = "n",
  xlab = "", 
  ylab = "",
)

# 绘制外接正方形
rect(
  xleft = -1, 
  ybottom = -1, 
  xright = 1, 
  ytop = 1,
  border = "blue", 
  lwd = 2
)

# 生成并绘制随机点
set.seed(123)
get_points = \(n) {
  df = data.frame(
    x = runif(n, min = -1, max = 1), 
    y = runif(n, min = -1, max = 1)
  )
  return(df)
}

points_data = get_points(200)

points(
  points_data$x, 
  points_data$y, 
  pch = 19, 
  col = "green"
)

图上，圆的面积就是 $\pi$，正方形的面积为 4。如果我们往正方形内随意撒豆子，豆子落在圆中的概率为 $p$，那么对 $p$ 的估计的期望值就是圆在正方形中的面积占比，公式表示如下：

$$ p = E(\hat{p}) = \frac{\pi}{4} $$

根据上面的公式很容易推导出：

$$ \pi = p * 4 = E(\hat{p}) * 4 $$

此时估计 \(\pi\) 就转化为估计 $p$ 的问题:

$$ \hat{\pi} = \hat{p} * 4 $$

根据大数定理和中心极限定理，假设我们要撒 $n$ 个豆子，当 $n$ 趋于无穷大时，豆子落在圆中个数除以豆子落在正方形内的个数以 $p$ 为极限。换句话说，$n$ 越大，估计的 $\pi$ 就会相对越准确。我们可以通过增加 $n$ 逐步逼近真实的 \(\pi\)。

# 编写估计圆周率的函数
set.seed(123)
pi_hat = \(n) {
  df = get_points(n)
  mean(df$x^2 + df$y^2 < 1) * 4
}

# 做 9 次撒豆子实验
n = list(
  n1 = 1, 
  n2 = 10, 
  n3 = 100, 
  n4 = 1000, 
  n5 = 10000,
  n6 = 100000,
  n7 = 1000000,
  n8 = 10000000,
  n9 = 100000000
)

# 生成估计的圆周率
lapply(n, pi_hat)

#> $n1
#> [1] 4
#> 
#> $n2
#> [1] 2.4
#> 
#> $n3
#> [1] 3.48
#> 
#> $n4
#> [1] 3.18
#> 
#> $n5
#> [1] 3.156
#> 
#> $n6
#> [1] 3.14364
#> 
#> $n7
#> [1] 3.142608
#> 
#> $n8
#> [1] 3.141214
#> 
#> $n9
#> [1] 3.141456

从估计的结果来看，随着撒豆子数量的增加，估计值总体上越来越接近真实的 $\pi = $ 3.1415927。