数概率

从抛硬币实验得出沃利斯公式

字数统计: 3.1k字   |   阅读时长: 15分
, , ,

摘要: 从伯努利试验序列中找到沃利斯公式

【对数据分析、人工智能、金融科技、风控服务感兴趣的同学，欢迎关注我哈，阅读更多原创文章】
我的网站：潮汐朝夕的生活实验室
我的公众号：潮汐朝夕
我的知乎：潮汐朝夕
我的github：FennelDumplings
我的leetcode：FennelDumplings

沃利斯公式是我们大学数学分析中的一个著名公式，由沃利斯 1655 年发表在《无穷算术》上，发表时的写法如下：

该公式是第一个把 $\pi$ 表示成容易计算的有理数列的极限的公式，在理论上很有意义。

《无穷算术》第一次把代数学扩展到分析学，发展了卡瓦列里的不可分量原理，开辟了用级数表示函数的道路，把有限算术变成无限算术，成为牛顿、莱布尼兹发明微积分的前奏。

现代数学分析教材中沃利斯公式的表达式如下：

教材中一般是先给出以上公式，然后再用现代分析学的理论体系去证明。这当然没错，但我们更感兴趣的是这个公式是怎么发现的，当年沃利斯是怎么想到去算一系列有理数的乘积，然后发现跟 $\pi$ 的关联的。

从一些数学史的资料上我们可以知道，沃利斯是在研究单位圆第一象限的面积相关的插值问题时得到的沃利斯公式，具体的步骤也很精彩。本文我们不讨论沃利斯当年的推导过程，而从一个似乎八竿子打不着的投硬币的伯努利试验出发，研究 $n$ 次试验中正反面次数之差 $T_{n}$ 的性质，最终推导出了与沃利斯公式吻合的表达式，推导过程中与组合相关的部分比较考验耐心。之后我们给出了模拟的代码，通过模拟投硬币试验的 $T_{n}$ 的结果测出了 $\pi$。

问题背景

考虑投掷硬币的伯努利试验，试验次数为 n，硬币正面的次数 $X \sim B(n, 0.5)$。

在试验中我们可以发现，随着 $n$ 增加，正面的频率逐渐稳定于 0.5。同时我们还可以发现，随着 $n$ 增加，硬币正反面次数的差（记为 $T_{n}$）也是增加的。在试验中我们发现 $T_{n} \rightarrow \infty$，但试验中我们同时发现 $\frac{T_{n}}{n}\rightarrow 0$，也就是说 $T_{n}$ 是比 $n$ 低阶的无穷大，所以我们可以【猜想 $T_{n}^{2}$ 会不会是与 $n$ 同阶的无穷大】，也就是 $\lim\limits_{n\rightarrow\infty}\frac{T_{n}^{2}}{n}$ 收敛于一个常数。

理论推导

$T_{n}$ 的定义

由于问题背景比较简单，研究对象清晰，我们可以进行一些理论推导。

投掷硬币 $n$ 次，正面记为 1，反面记为 0，得到的长为 $n$ 的 01 序列共有 $2^{n}$ 种。其中出现 $m$ 次正面的事件数量为 $\binom{n}{m}$，$0 \leq m \leq n$，其正反面的差为 $|n-2m|$。

记 $Q_{n} = \sum\limits_{m=0}\limits^{n}\binom{n}{m}|n-2m|$，平均正反面的差即为 $T_{n} = \frac{Q_{n}}{2^{n}}$。

根据以上公式，我们进行一些计算，看看 n 比较小的时候的一些值。

1
2
3
4
5
6
7
8
import math

for n in range(5, 13):
    Qn = 0
    for m in range(n + 1):
        Qn += math.comb(n, m) * abs(n - m * 2)
    Tn = Qn / math.pow(2, n)
    print("{}|{}|{}|{}".format(n, Qn, Tn, Tn * Tn / n))

结果如下：

$n$ $Q_{n}$ $T_{n}$ $\frac{T_{n}^{2}}{n}$
5 60 1.875 0.703125
6 120 1.875 0.5859375
7 280 2.1875 0.68359375
8 560 2.1875 0.59814453125
9 1260 2.4609375 0.67291259765625
10 2520 2.4609375 0.605621337890625
11 5544 2.70703125 0.6661834716796875
12 11088 2.70703125 0.6106681823730469

$T_{2n} = T_{2n-1}$

从结果表格中我们发现，平均正反面的差 $T_{n}$ 成对出现，于是可以猜想 $T_{2n} = T_{2n+1}$。证明如下：

首先看 $Q_{2n-1}$：

当 $m$ 取 $0$ 到 $n-1$ 时：

因此：

然后看 $Q_{2n}$：

类似地，当 $m$ 取 $0$ 到 $n-1$ 时：

另外 $m = n$ 时，$|2n-2m|=0$，因此：

将 $Q_{2n-1}$ 和 $Q_{2n}$ 展开：

$$ \begin{align} Q_{2n-1} &= 2\sum\limits_{m=0}\limits^{n-1}(\binom{2n-1}{m}(2n-1-2m)) \\ &= 2(\binom{2n-1}{0}(2n-1) + \binom{2n-1}{1}(2n-3) + \binom{2n-1}{2}(2n-5) + \cdots + \binom{2n-1}{n-2}\cdot 3 + \binom{2n-1}{n-1}\cdot 1) \\ \end{align} $$ $$ \begin{align} Q_{2n} &= 2\sum\limits_{m=0}\limits^{n-1}(\binom{2n}{m}(2n-2m)) \\ &= 2(\binom{2n}{0}(2n) + \binom{2n}{1}(2n-2) + \binom{2n}{2}(2n-4) + \cdots + \binom{2n}{n-2}\cdot 4 + \binom{2n}{n-1}\cdot 2) \\ \end{align} $$

由帕斯卡恒等式 $\binom{n}{k} = \binom{n-1}{k} + \binom{n-1}{k-1}$，可以进一步推导 $Q_{2n}$：

$$ \begin{align} Q_{2n} &= 2(\binom{2n}{0}(2n) + \binom{2n}{1}(2n-2) + \binom{2n}{2}(2n-4) + \cdots + \binom{2n}{n-2}\cdot 4 + \binom{2n}{n-1}\cdot 2) \\ &= 2(\binom{2n-1}{0}(2n) + (\binom{2n-1}{0} + \binom{2n-1}{1})(2n-2) + (\binom{2n-1}{1} + \binom{2n-1}{2})(2n-4) + \cdots \\ & \quad\quad\quad\quad\quad\quad+ (\binom{2n-1}{n-3} + \binom{2n-1}{n-2})\cdot 4) + (\binom{2n-1}{n-2} + \binom{2n-1}{n-1})\cdot 2) \\ &= 2(\binom{2n-1}{0}(4n-2) + \binom{2n-1}{1}(4n-6) + \cdots + \binom{2n-1}{n-2}\cdot 6 + \binom{2n-1}{n-1}\cdot 2) \\ &= 2\times 2(\binom{2n-1}{0}(2n-1) + \binom{2n-1}{1}(2n-3) + \cdots + \binom{2n-1}{n-2}\cdot 3 + \binom{2n-1}{n-1}\cdot 1) \\ \end{align} $$

因此 $Q_{2n} = 2Q_{n-1}$，于是 $T_{2n} = T_{2n-1}$。

$\frac{T_{2n}}{T_{2n+1}} = \frac{2n+1}{2n}$

$T_{2n}$ 与 $T_{2n-1}$ 的关系是从表格计算结果中猜想后证明的。而 $T_{2n}$ 与 $T_{2n+1}$ 的关系就不是很好从表格中猜想的，这里我们进行正面推导。

与推导 $Q_{n-1}$ 的过程一样，首先根据组合数和绝对值的对称性，提出因子 2 并拆掉绝对值，然后用帕斯卡恒等式将 $\binom{2n+1}{i}$ 形式的组合数变为 $\binom{2n}{i}$ 的形式，然后合并同类项，与 $Q_{2n}$ 对比得到结果：

$$ \begin{align} Q_{2n+1} &= \sum\limits_{m=0}\limits^{2n+1}(\binom{2n+1}{m}|2n+1-2m|) \\ &= 2\sum\limits_{m=0}\limits^{n}(\binom{2n+1}{m}(2n+1-2m)) \\ &= 2(\binom{2n+1}{0}(2n+1) + \binom{2n+1}{1}(2n-1) + \binom{2n+1}{2}(2n-3) + \cdots + \binom{2n+1}{n-1}\cdot 3 + \binom{2n+1}{n}\cdot 1) \\ &= 2(\binom{2n}{0}(2n+1) + (\binom{2n}{0} + \binom{2n}{1})(2n-1) + (\binom{2n}{1} + \binom{2n}{2})(2n-3) + \cdots \\ & \quad\quad\quad\quad\quad\quad + (\binom{2n}{n-2} + \binom{2n}{n-1})\cdot 3 + (\binom{2n}{n-1} + \binom{2n}{n})\cdot 1) \\ &= 2(\binom{2n}{0}(4n) + \binom{2n}{1}(4n-4) + \cdots + \binom{2n}{n-1}\cdot 4 + \binom{2n}{n}) \\ &= 2Q_{2n} + 2\binom{2n}{n} \end{align} $$

由 $Q_{2n} = 2Q_{n-1}$、$Q_{2n+1} = 2Q_{2n} + 2\binom{2n}{n}$，以及 $Q_{1}=2$，可以解出 $Q_{2n} = 2n\binom{2n}{n}$、$Q_{2n+1} = (4n+2)\binom{2n}{n}$。

于是：

可以推出 $\frac{T_{2n+1}}{T_{2n}} = \frac{2n+1}{2n}$。

$T_{n}$ 的通项公式

在前面的推导中，我们得到了 $T_{n}$ 的通项公式。用通项公式把 $n = 5,6,\cdots,12$ 时的数据算一下，结果完全一样，代码如下：

1
2
3
4
5
6
7
8
9
10
import math
import numpy as np

for k in range(2, 7):
    n = 2 * k
    Tn = math.comb(2 * k, k) * k * 2 / math.pow(2, 2 * k)
    print("{}|{}|{}".format(n, Tn, Tn * Tn / n))
    n = 2 * k + 1
    Tn = math.comb(2 * k, k) * (k * 2 + 1) / math.pow(2, 2 * k)
    print("{}|{}|{}".format(n, Tn, Tn * Tn / n))

数列 $\{\frac{T_{n}^{2}}{n}\}$ 的收敛性

令 $\{a_{n}\}$ 为 $\{\frac{T_{n}^{2}}{n}\}$ 的奇数项子数列，$a_{n} = \frac{T_{2n-1}^{2}}{2n-1}$。

令 $\{b_{n}\}$ 为 $\{\frac{T_{n}^{2}}{n}\}$ 的偶数项子数列，$b_{n} = \frac{T_{2n}^{2}}{2n}$。

对于 $\{a_{n}\}$，可以推导出 $a_{1} = 1$，$\frac{a_{n+1}}{a_{n}} = \frac{4n^{2}-1}{4n^{2}} < 1$，因此 $\{a_{n}\}$ 单调递减。

对于 $\{b_{n}\}$，可以推导出 $b_{1} = \frac{1}{2}$，$\frac{b_{n+1}}{b_{n}} = \frac{4n^{2}+4n+1}{4n^{2}+4n} > 1$，因此 $\{b_{n}\}$ 单调递增。

因为 $T_{2n} = T_{2n-1}$，所以 $a_{n} > b_{n}$。于是 $n > 1$ 时，$\frac{1}{2} = b_{1} < b_{n} < a_{n} < a_{1} < 1$，于是 $\{a_{n}\}$ 有下界，$b_{n}$ 有上界。

由单调有界定理，$\{a_{n}\}$ 和 $\{b_{n}\}$ 均收敛，记 $A = \lim\limits_{n\rightarrow\infty}a_{n}$，$B = \lim\limits_{n\rightarrow\infty}b_{n}$。

因为 $T_{2n} = T_{2n-1}$，所以 $A = B$。所以数列 $\{\frac{T_{n}^{2}}{n}\}$ 收敛到 $A = B$。

下面求 $B = \lim\limits_{n\rightarrow\infty}\frac{T_{2n}^{2}}{2n}$。

求 $\lim\limits_{n\rightarrow\infty}\frac{T_{n}^{2}}{n}$

利用斯特林公式：$\lim\limits_{n \rightarrow \infty} \frac{n!}{(\frac{n}{e})^{n}\sqrt{2\pi n}} = 1$

因此 $\lim\limits_{n\rightarrow\infty}\frac{T_{n}^{2}}{n} = \frac{2}{\pi}$，由此得到 $\pi$ 的另一个公式：

导出沃利斯公式

由 $T_{2n - 1} = T_{2n}$，$\frac{T_{2n+1}}{T_{2n}} = \frac{2n+1}{2n}$，以及 $b_{n} = \frac{T_{2n}^{2}}{2n}$。有：

于是 $b_{n+1} = \frac{2n+1}{2n}\cdot\frac{2n+1}{2n+2}\cdot b_{n}$：

这正是沃利斯公式：

模拟结果

直接根据 $T_{n}$ 的定义，模拟投硬币试验并记录 $T_{n}$ 和 $\frac{T_{n}^{2}}{n}$ 的结果。代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
import numpy as np
import math

L = int(4e5) # 试验次数
N = int(1e3) # 一次试验的投掷次数
tests = np.random.random_integers(0, 1, size=(L, N))
n = np.arange(N) + 1

def get_Tn(nums):
    sums = np.cumsum(nums)
    Tn = np.abs(n - sums * 2)
    return Tn

Tn = np.zeros((L, N))
for i in range(L):
    Tn[i] = get_Tn(tests[i])

Tn_avg = np.zeros(N)
for i in range(N):
    Tn_avg[i] = np.mean(Tn.T[i])

result = Tn_avg * Tn_avg / n

for i in (5, 6, 7, 8, 9, 10, 11, 12, 100, 200, 400, 800, 1000):
    print("{}|{}|{}|{}".format(i, Tn_avg[i - 1], result[i - 1], 2 / result[i - 1]))

模拟结果如下：

$n$ $T_{n}$ $\frac{T_{n}^{2}}{n}$ $\frac{2n}{T_{n}^{2}}$
5 1.87719 0.70476845922 2.8378114454972816
6 1.879865 0.5889820697041667 3.395689109864682
7 2.191715 0.6862306630321429 2.9144719228413734
8 2.194105 0.601762093878125 3.3235725884805576
9 2.465125 0.6752045850694445 2.962065193609876
10 2.46335 0.6068093222500001 3.295928270488926
11 2.7106 0.6679411236363636 2.9942758863411854
12 2.709445 0.6117576840020834 3.2692682941325977
100 7.96788 0.634871116944 3.1502456902231604
200 11.28695 0.6369762015124999 3.139834730482866
400 15.956695 0.6365402883075625 3.1419849406824087
800 22.594895 0.6381616000762813 3.1340024215824553
1000 25.26087 0.6381115531569 3.134248220558759

总结

在网上可以找到众多通过实验得到圆周率的方法，但是基本上都是预设了圆形然后再投针的方法。本文从伯努利试验中找到了 $\pi$，并没有预设图形，而且吻合上了著名的沃利斯公式。一切的起点是【猜想 $T_{n}^{2}$ 会不会是与 $n$ 同阶的无穷大】，当有了这个想法，后面的推导都是现代分析学和组合数学基础知识的应用。这也启示我们，通过现象提出靠谱的猜想是一个非常牛的能力。

Share