微不足道——但很神奇：

BBP <- function(n = 13) {
  sum(sapply(seq_len(n), function(k) {
    ((sample.int(8*k+1, 1) <= 4) -
       (sample.int(8*k+4, 1) <= 2) -
       (sample.int(8*k+5, 1) <= 1) -
       (sample.int(8*k+6, 1) <= 1)) / 16^k
  })) + (4 - 2/4 - 1/5 - 1/6)
}

正如您在此代码中看到的，R仅对使用. 默认情况下，仅次绘制（值从不大于）——但结果的预期值是到完全双精度！sample.int$13*4=52$$110$$\pi$

这是次迭代的示例运行（需要一秒钟的时间）：$10,000$

x <- replicate(1e4, BBP())
mu <- mean(x)
se <- sd(x) / sqrt(length(x))
signif(c(Estimate=mu, SE=se, Z=(mu-pi)/se), 4)

它的输出是

 Estimate        SE         Z 
3.1430000 0.0004514 2.0870000

换句话说，的这个（随机）估计值为的小 Z 值的真实值没有显着偏差。$\pi$$3.143\pm 0.00045$$2.08$$\pi.$

---

这是微不足道的，因为我希望代码很明显，像sample.int(b,1) <= a（当整数a不超过时b）这样的计算只是估计有理分数的愚蠢方法a/b。因此，此代码估计Bailey Borwein Plouffe 公式

$$\pi = \sum_{k=0}^\infty \frac{1}{16^k}\left(\frac{4}{8k+1}-\frac{2}{8k+4}-\frac{1}{8k+5}-\frac{1}{8k+6}\right)$$

通过明确表达项并通过 由于公式中的每一项在额外的位，因此在该点终止采样会在二进制点之后给出位，这略高于最大位的总精度在 . 使用的 IEEE 双精度浮点数中可用。$k=0$$k=13.$$4$$\pi,$$4*(13)=52$$52$ R

尽管我们可以通过分析计算出方差，但前面的示例已经为我们提供了一个很好的估计，因为标准误差仅为每次迭代的方差为$0.0045,$$0.002$

var(x)

[1] 0.002037781

因此，如果您想将BBP估计在 sigma 的标准误差内您将需要大约次迭代。例如，以这种方式估计到小数点后六位将需要大约 20 亿次迭代（大约需要三天的计算）。$\pi$$\sigma,$$0.002/\sigma^2$$\pi$

减少方差（极大）的一种方法是一劳永逸地计算 BBP 和中的一些初始项，仅使用蒙特卡罗模拟来估计结果的最低有效位:-)。

生成一对整数并查看它是否在象限内，即让数字表示为。如果，则该点位于象限内。如果数字是连续的，则概率将为。因此，增加会增加估计的精度。下面是一个 Python one 班轮：$[0,M)$$m_1, m_2$$m_1^2+m_2^2<M^2$$\pi/4$$M$

N = int(1e8)
M = int(1e9)
4 * np.mean(np.sum(np.random.randint(0, M, (2, N))**2, axis=0) < M**2)

为了产生概率，可以使用以下算法（Flajolet et al. 2010），该算法基于 Ramanujan 的级数展开：$1/\pi$

1. 将设置为 0。$t$
2. 翻转两个公平的硬币。如果两个都显示头，则将 1 添加到并重复此步骤。否则，转到步骤 3。$t$
3. 翻转两个公平的硬币。如果两个都显示头，则将 1 添加到并重复此步骤。否则，转到步骤 4。$t$
4. 以 5/9 的概率，将 1 添加到。（例如，在 [1, 9] 中生成一个统一的随机整数，如果该整数小于或等于 5，则将 1 添加到。）$t$$t$
5. 掷一枚公平的硬币次，如果正面出现的次数多于反面，则返回 0，反之亦然。再执行此步骤两次。$2t$
6. 返回 1。

然后，运行上面的算法，直到得到 1，然后让为包含最后一次的运行次数。的期望值为。$X$$X$$\pi$

请注意，该算法根本不涉及分数。

另请参阅：https ://math.stackexchange.com/questions/4189867/obtaining-irrational-probabilities

参考：

• Flajolet, P.、Pelletier, M.、Soria, M.，“关于布冯机器和数字”，arXiv:0906.5560 [math.PR]，2010。