11._斯特林公式与概率 - Probability Theory and Mathematical Statistics

斯特林公式

在没有阶乘的情况下研究概率是可行的，但你必须做很多工作来避免它们．记住 $n!$ 是前 $n$ 个正整数的乘积，而且我们约定 $0!=1$ ．因此 $3!=3 \cdot 2 \cdot 1=6$ 且 $4!=4 \cdot 3!=24$ ．有个很好的组合学解释： $n!$ 表示把 $n$ 个人有序排成一列的方法数（第一个人有 $n$ 种选择，第二个人有 $n-1$ 种选择，依此类推）．有了这个观点，我们就可以解释 $0!=1$ 意味着什么都不安排的方法只有一种！

但是，这里的阶乘函数并不十分明显．最隐蔽的一种情况可能发生在标准正态分布的概率密度函数中，即 $\frac{1}{\sqrt{2 \pi}} \exp \left(-x^2 / 2\right)$ ．事实上， $\sqrt{\pi}=(-1 / 2)$ ！．当读到这里时，你心里应该响起了警报声．毕竟，我们已经定义了整数输入的阶乘函数，但现在取阶乘的不仅仅是负数，还是一个负的有理数！这是什么意思呢？我们该如何理解 $-1 / 2$ 的阶乘？对把 $-1 / 2$ 个人排成一列的方法数提问是什么意思？

对这个问题的回答要用到伽马函数，即 $\Gamma(s)$ ．

\Gamma(s)=\int_0^{\infty} e^{-x} x^s \frac{d x}{x}

关于伽马函数可以到《高等数学伽玛函数》里查看证明,但是解决需要记住： $\Gamma(\frac{1}{2})=\sqrt{\pi}$

斯特林公式

斯特林公式：当 $n \rightarrow \infty$ 时，我们有

n!\approx n^n e^{-n} \sqrt{2 \pi n} ;

这意味着

\lim _{n \rightarrow \infty} \frac{n!}{n^n e^{-n} \sqrt{2 \pi n}}=1

更准确地说，有下列级数展开式：

n!=n^n e^{-n} \sqrt{2 \pi n}\left(1+\frac{1}{12 n}+\frac{1}{288 n^2}-\frac{139}{51840 n^3}-\cdots\right) .

斯特林公式与概率

在深入了解斯特林公式的证明细节之前，我们先来看看如何使用它．第一个问题的灵感来源于一种让许多学生感到惊讶的现象。假设有 1 枚均匀的硬币，它掷出正面和反面的概率各占一半。如果抛郑 $2 n$ 次，我们预计会得到 $n$ 个正面．所以，当抛掷 200 万次时，预计会出现 100 万个正面．但事实证明，当 $n \rightarrow \infty$ 时，抛掷一枚均匀的硬币 $2 n$ 次，恰好出现 $n$ 个正面的概率趋近于 0 ！

这个结果看起来很惊人。在 $2 n$ 次抛掷中，出现正面的期望数是 $n$ ，得到 $n$ 个正面的概率却趋向于 0 ？这是怎么回事？如果这是期望值，那么它不应该很有可能发生吗？理解这个问题的关键是要注意，在 $2 n$ 次抛郑中，虽然出现正面的期望值是 $n$ ，但标准差是 $\sqrt{n / 2}$ ．回到抛掷 200 万次的情况，我们预计会出现 100 万个正面，但其波动幅度是 700 ．如果现在把硬币抛郑 2 万亿次，那么预计有 1 万亿个正面，而波动幅度约为 700000 。伴随着 $n$ 的增加，描述可能出现结果的均值＂窗口＂会像标准差那样增长，也就是说，它像 $\sqrt{n}$ 那样增长（直到某些常数）．因此，概率所对应的取值会越来越多，这样就使得单独一个特定结果（比如恰好出现 $n$ 次正面）的概率不断下降。总而言之：随着概率覆盖的集合越来越大， $2 n$ 次抛郑中恰好出现 $n$ 个正面的概率就会不断下降。现在，我们用斯特林公式来量化这个概率的下降速度．

对于取值较大的 $n$ ，我们可以利用斯特林公式来轻松地近似 $n!$ ．记住，我们正在努力回答下面这个问题：抛掷一枚均匀的硬币 $2 n$ 次，恰好得到 $n$ 个正面的概率是多少？

答案可以很轻松地写出来：就是

\operatorname{P}(\text { 在 } 2 n \text { 次抛掷中恰好出现 } n \text { 个正面 })=\binom{2 n}{n}\left(\frac{1}{2}\right)^n\left(\frac{1}{2}\right)^n=\binom{2 n}{n} \frac{1}{2^{2 n}} \text {. }

原因是，在 $2^{2 n}$ 个可能的结果中，每一个结果出现的概率都是相等的，并且恰好出现 $n$ 个正面的结果共有 $\binom{2 n}{n}$ 个．另外，也可以把第二个 $(1 / 2)^n$ 看作 $(1-1 / 2)^{2 n-n}$ ．

那么 $\binom{2 n}{n}$ 是多大？是时候使用斯特林公式了．我们有

\begin{aligned} \binom{2 n}{n} & =\frac{(2 n)!}{n!n!} \\ & \approx \frac{(2 n)^{2 n} e^{-2 n} \sqrt{2 \pi \cdot 2 n}}{n^n e^{-n} \sqrt{2 \pi n} \cdot n^n e^{-n} \sqrt{2 \pi n}} \\ & =\frac{2^{2 n}}{\sqrt{\pi n}} \end{aligned}

因此，恰好得到 $n$ 个正面的概率是

\binom{2 n}{n} \frac{1}{2^{2 n}} \approx \frac{1}{\sqrt{\pi n}}

这意味着，如果 $n=100$ ，那么得到一半正面的概率略小于 $6 \%$ ．但是，如果 $n=$ 1000000 ，这个概率就会小于 $0.06 \%$ ．

实际上, 上面的练习类似于对 2008 年美国总统初选的讨论.2008 年美国总统初选时．希拉里和奥巴马分别在纽约州锡拉丘兹的民主党初选中获得了 6001 张选票．虽然初选共有 12346 张选票，但为了简单起见，我们假设只有 12002 张选票．请问出现平局的概率是多少？如果假设每位候选人都等可能地获得任何 1 张选票，答案就是 $\binom{12002}{6001} / 2^{12002}$ ．精确答案约等于 0.0072829 ．利用上述近似，结果是

\frac{1}{\sqrt{\pi \cdot 6001}} \approx 0.00728305

这与真实值非常接近．

这是一份非常有意思的结论, 12002张选票，希拉里和奥巴马每个人各获得6001张，理论上双方平局的结果应该是100%，但是实际上平局的结果只有0.728%

虽然我们得到的概率略低于 $1 \%$ ，但一些新闻媒体报道说这个概率大约是百万分之一，有些人甚至说这种情况＂几乎不可能＂发生。为什么会有如此不同的答案？这完全取决于你对问题如何建模。如果假设两位候选人会等可能地获得每 1张选票，那么得到的答案就会略少于 $1 \%$ ．但如果考虑其他信息，情况就会发生变化．例如，希拉里当时是纽约州的一名参议员，所以不难想到她在家乡的选票会更多．事实上，她以 $57.37 \%$ 的得票率在这个州获胜，而奥巴马的得票率为 $40.32 \%$ ．同样地，为了便于分析，我们假设希拉里获得了 $57 \%$ 的选票，而奥巴马得到了 $43 \%$ 的选票．如果把这些考虑进来，就不能假设每个选民都等可能地投票给希拉里或奥巴马，此时出现平局的概率只有 $\binom{12}{6001} \cdot 0.57^{6001} \cdot 0.43^{6001}$ ．由斯特林公式可得， $\binom{2 n}{n} \approx 2^{2 n} / \sqrt{\pi n}$ ．因此，在这些假设下，出现平局的概率约等于

\frac{2^{12002}}{\sqrt{\pi \cdot 6001}} \cdot 0.57^{6001} 0.43^{6001} \approx 1.877 \cdot 10^{-54}

请注意，概率之间的差别取决于我们对预期的不同假设！

从斯特林公式到中心极限定理

中心极限定理是概率论中非常重要的结论之一。它指出随着相互独立的随机变量个数不断增加，这些变量的和会趋向于服从正态分布。作为斯特林公式的一个强有力的应用，我们将证明斯特林公式蕴含着中心极限定理的一种特殊情形，即随机变量 $X_1, \cdots, X_{2 N}$ 都服从参数为 $p=1 / 2$ 的伯努利分布的情形。从技术角度考虑，最简单的情形是下列标准化形式

\operatorname{Prob}\left(X_i=n\right)= \begin{cases}1 / 2 & \text { 若 } n=1 \\ 1 / 2 & \text { 若 } n=-1 \\ 0 & \text { 其他. }\end{cases}

把1看作硬币郑出了正面，-1 看作掷出了反面。从这个标准化形式中可以看出， $X_1+\cdots+X_{2 N}$ 的期望值是 0 。由于 $X_i$ 表示第 $i$ 次掷出了正面或反面，所以这个和就是正面比反面多出的个数（期望值为 0 ）．设 $X_1, \cdots, X_{2 N}$ 是相互独立的随机变量，且均服从概率密度函数为式（18．1）的二项分布．那么，每个变量的均值都是 0 ，因为 $1 \cdot(1 / 2)+(-1) \cdot(1 / 2)=0$ ，而方差均为

\sigma^2=(1-0)^2 \cdot \frac{1}{2}+(-1-0)^2 \cdot \frac{1}{2}=1

最后，我们令

S_{2 N}=X_1+\cdots+X_{2 N} .

由

E \left[S_{2 N}\right]= E \left[X_1\right]+\cdots+ E \left[X_{2 N}\right]=0+\cdots+0=0

可知， $S_{2 N}$ 的均值是 0 ．类似地，不难看出 $S_{2 N}$ 的方差是 $2 N$ ．因此，我们希望 $S_{2 N}$ 约等于 0 ，而波动幅度约为 $\sqrt{2 N}$ ．

现在看一下 $S_{2 N}$ 服从的分布．我们首先注意到， $S_{2 N}=2 k+1$ 的概率是 0 ．这是因为 $S_{2 N}$ 等于出现正面的次数减去出现反面的次数，而这个值始终是偶数：如果出现 $k$ 次正面和 $2 N-k$ 次反面，那么 $S_{2 n}$ 就等于 $2 N-2 k$ ． $S_{2 N}$ 等于 $2 k$ 的概率是 $\binom{2 N}{N+k}\left(\frac{1}{2}\right)^{N+k}\left(\frac{1}{2}\right)^{N-k}$ ．理由如下：当 $S_{2 N}$ 等于 $2 k$ 时， 1 （正面）比 -1 （反面）多了 $2 k$ 个，并且 1 和 -1 的总个数是 $2 N$ 。因此，我们有 $N+k$ 个正面（1）和 $N-k$ 个反面（ -1 ）。所有可能的结果一共有 $2^{2 N}$ 种，恰好郑出 $N+k$ 个正面和 $N-k$ 个反面的结果共有 $\binom{2 N}{N+k}$ 种，而且每种结果出现的概率均为 $\left(\frac{1}{2}\right)^{2 N}$ ．为了说明该如何处理更一般的情形（即掷出正面的概率为 $p$ ，掷出反面的概率为 $1-p$ ），我们把上面的概率写成 $\left(\frac{1}{2}\right)^{N+k}\left(\frac{1}{2}\right)^{N-k}$ ．现在用斯特林公式来近似 $\binom{2 N}{N+k}$ 。我们有

\begin{aligned} \binom{2 N}{N+k} & \approx \frac{(2 N)^{2 N} e^{-2 N} \sqrt{2 \pi \cdot 2 N}}{(N+k)^{N+k} e^{-(N+k)} \sqrt{2 \pi(N+k)}(N-k)^{N-k} e^{-(N-k)} \sqrt{2 \pi(N-k)}} \\ & =\frac{(2 N)^{2 N}}{(N+k)^{N+k}(N-k)^{N-k}} \sqrt{\frac{N}{\pi(N+k)(N-k)}} \\ & =\frac{2^{2 N}}{\sqrt{\pi N}} \frac{1}{\left(1+\frac{k}{N}\right)^{N+\frac{1}{2}+k}\left(1-\frac{k}{N}\right)^{N+\frac{1}{2}-k}} . \end{aligned}

剩下的部分就是做一些代数运算来证明这个随机变量收敛于正态分布．不幸的是，在处理因子时，人们经常会掉进一个常见的陷阱。为了避免将来出现这样的问题，我们先描述这个常见错误，然后完成证明。

我们想利用 $e ^x$ 的定义来推导：当 $N \rightarrow \infty$ 时， $\left(1+\frac{w}{N}\right)^N \approx$ $e ^w$ 。然而必须小心一点，因为 $k$ 值会随着 $N$ 增长．例如，我们可能会认为 $\left(1+\frac{k}{N}\right)^N \rightarrow$$e ^k$ 且 $\left.\left(1-\frac{k}{N}\right)\right)^N \rightarrow e ^{-k}$ ，这样两个因子就相互抵消了。因为 $k$ 相对于 $N$ 较小，所以我们可能会忽略 $1 / 2$ ，并给出

\left(1+\frac{k}{N}\right)^k=\left(1+\frac{k}{N}\right)^{N \cdot \frac{k}{N}} \rightarrow e^{k^2 / N}

类似地，有 $\left(1-\frac{k}{N}\right)^{-k} \rightarrow e ^{k^2 / N}$ 。于是，我们可能会断言（稍后在引理18．3．1中会看到，这一说法是错误的！）

\left(1+\frac{k}{N}\right)^{N+\frac{1}{2}+k}\left(1-\frac{k}{N}\right)^{N+\frac{1}{2}-k} \rightarrow e^{2 k^2 / N}

我们来证明 $\left(1+\frac{k}{N}\right)^{N+\frac{1}{2}+k}\left(1-\frac{k}{N}\right)^{N+\frac{1}{2}-k} \rightarrow e ^{k^2 / N}$ 。这一计算的重要性在于它强调了收敛速度的重要性。虽然主要部分 $\left(1 \pm \frac{k}{N}\right)^N$ 等于 $e ^{ \pm k}$ 是正确的，但误差项（在极限中）也相当重要。当 $k$ 是 $N$ 的方幂时，误差会产生一个很大的项。这里的情况是，由这两个因子生成的项会相互强化。另一种说法是，一个因子会趋向于无穷大，而另一个则趋向于 0 。记住 $\infty \cdot 0$ 是无定义的表达式之一，根据各项增长或衰减速度的快慢，它可以是任何值．本节最后会给出更多说明．

简而言之，不能只利用 $\left(1+\frac{w}{N}\right)^N \approx e ^w$ 。我们要更加小心。正确的做法是对两个因子取对数，并用泰勒公式展开，然后做指数运算．这样我们就能更好地处理误差项．在做这些之前，我们需要大致了解 $k$ 的重要取值范围是什么．因为标准差是 $\sqrt{2 N}$ ，所以预计唯一真正重要的 $k$ 是那些与 0 相距几个标准差的值．也就是说， $k$ 的最大值会比 $\sqrt{2 N}$ 略大。我们可以利用切比雪夫不等式（定理 17．3．1）来仔细地量化需要研究多大的 $k$ 。由此可知，只需考察那些 $|k|$ 不大于 $N^{\frac{1}{2}+\epsilon}$ 的 $k$ ，这是因为 $S_{2 N}$ 的标准差是 $\sqrt{2 N}$ ．于是，由 $(2 N)^{1 / 2+\epsilon}=(2 N)^\epsilon \operatorname{StDev}\left(S_{2 N}\right)$ 可得

\operatorname{Prob}\left(\left|S_{2 N}-0\right| \geqslant(2 N)^{1 / 2+\epsilon}\right) \leqslant \frac{1}{(2 N)^{2 \epsilon}}

因此，现在只需要考察当 $|k| \leqslant N^{1 / 2+1 / 9}$ 时， $S_{2 N}=2 k$ 的概率．现在来考察前面所说的引理，它会给出这个乘积的正确极限，而这个证明将告诉我们该如何解决这类问题。

引理 18．3．1 对于任意的 $\epsilon \leqslant 1 / 9$ ，当 $|k| \leqslant(2 N)^{1 / 2+\epsilon}$ 时，如果 $N \rightarrow \infty$ ，那么有

\left(1+\frac{k}{N}\right)^{N+\frac{1}{2}+k}\left(1-\frac{k}{N}\right)^{N+\frac{1}{2}-k} \longrightarrow e^{k^2 / N} e^{O\left(N^{-1 / 6}\right)}

证明：回顾一下，当 $|x|<1$ 时，

\log (1+x)=\sum_{n=1}^{\infty} \frac{(-1)^{n+1} x^n}{n}

因为假设 $k \leqslant(2 N)^{1 / 2+\epsilon}$ ，所以任何小于 $k^2 / N^2, k^3 / N^2$ 或 $k^4 / N^3$ 的项都可以忽略不计。因此，如果定义

P_{k, N}:=\left(1+\frac{k}{N}\right)^{N+\frac{1}{2}+k}\left(1-\frac{k}{N}\right)^{N+\frac{1}{2}-k}

可以得到

\begin{aligned} \log P_{k, N}= & \left(N+\frac{1}{2}+k\right) \log \left(1+\frac{k}{N}\right)^{N+\frac{1}{2}+k} \\ & +\left(N+\frac{1}{2}-k\right) \log \left(1-\frac{k}{N}\right)^{N+\frac{1}{2}-k} \\ = & \left(N+\frac{1}{2}+k\right)\left(\frac{k}{N}-\frac{k^2}{2 N^2}+O\left(\frac{k^3}{N^3}\right)\right) \\ & +\left(N+\frac{1}{2}-k\right)\left(-\frac{k}{N}-\frac{k^2}{2 N^2}+O\left(\frac{k^3}{N^3}\right)\right) \\ = & \frac{2 k^2}{N}-2\left(N+\frac{1}{2}\right) \frac{k^2}{2 N^2}+O\left(\frac{k^3}{N^2}+\frac{k^4}{N^3}\right) \\ = & \frac{k^2}{N}+O\left(\frac{k^2}{N^2}+\frac{k^3}{N^2}+\frac{k^4}{N^3}\right) . \end{aligned}

当 $\epsilon<1 / 9$ 时， $k \leqslant(2 N)^{1 / 2+\epsilon}$ ，所以大 $O$ 项完全由 $N^{-1 / 6}$ 来确定．最后，我们得到了

P_{k, N}=e^{k^2 / N} e^{O\left(N^{-1 / 6}\right)}

结论得证．现在，我们完成了 $S_{2 N}$ 收敛于高斯分布的证明．结合起来可得

\binom{2 N}{N+k} \frac{1}{2^{2 N}} \approx \frac{1}{\sqrt{\pi N}} e^{-k^2 / N}

（通过仔细地分析引理，我们会注意到 $e ^{-k^2 / N}$ 的存在，这是我们快速而松散的计算所遗漏的）。在这个例子中，对中心极限定理的证明完全由一些简单的代数运算来完成。因为我们研究的是 $S_{2 N}=2 k$ ，所以应该把 $k^2$ 写成 $(2 k)^2 / 4$ 。类似地，因为 $S_{2 N}$ 的方差是 $2 N$ ，所以应该用 $(2 N) / 2$ 来代替 $N$ 。虽然这些看起来是不重要的代 18.4 积分判别法与较弱的斯特林公式 481

数技巧，但能轻松地做到这一点会非常有用．通过这些小的调整，我们可以更轻松地将表达式与猜测值进行比较。现在有

\operatorname{P }\left(S_{2 N}=2 k\right)=\binom{2 N}{N+k} \frac{1}{2^{2 N}} \approx \frac{2}{\sqrt{2 \pi \cdot(2 N)}} e^{-(2 k)^2 / 2(2 N)}

数技巧，但能轻松地做到这一点会非常有用．通过这些小的调整，我们可以更轻松地将表达式与猜测值进行比较。现在有

\operatorname{Prob}\left(S_{2 N}=2 k\right)=\binom{2 N}{N+k} \frac{1}{2^{2 N}} \approx \frac{2}{\sqrt{2 \pi \cdot(2 N)}} e^{-(2 k)^2 / 2(2 N)} .

回忆一下， $S_{2 N}$ 不可能是奇数．在上述标准化常数的分子中，因子 2 反映了这样的事实： $S_{2 N}$ 为偶数的概率是我们所期望的 2 倍，原因在于＂ $S_{2 N}$ 为奇数的概率是 0 ＂也必须被考虑到．所以，这看起来像是一个均值为 0 且方差为 $2 N$ 的高斯分布．对于较大的 $N$ ，这样的高斯分布是缓慢变化的，从 $2 k$ 到 $2 k+2$ 的积分基本上是 $2 / \sqrt{2 \pi(2 N)} \cdot \exp -(2 k)^2 / 2(2 N)$ 。

由于我们的证明很长，现在花点时间来回顾一下要点．非常幸运，我们有一个明确的概率公式，其中还包含了二项式系数．利用切比雪夫不等式，我们确定了概率的考察范围．然后利用斯特林公式做展开，并进行了一些代数运算，从而使我们的表达式看起来像一个高斯函数。

这里有个不错的挑战：你可以把上面的论述推广到 $p \neq 1 / 2$ 的情形吗？