正态分布(Normal Distribution)公式为什么长这样？

By Long Luo

相信大家或多或少都听过六西格玛( \(\text{6 Sigma}\) ) ¹ 这个词，六西格玛是指生产的产品中， \(99.99966\%\) 的产品是没有质量问题的，即只有 \(3.4ppm\) 的不良率。

假如一家工厂生产某型号零件，零件的长度要求是 \(100mm\) ，允许的标准差是 \(0.1mm\) 。根据 \(6 \sigma\) 原则，零件规格允许的偏差范围是： \(100 \pm 6 \times 0.1 = 100 \pm 0.6\) 。

这意味着，零件长度超过 \(100.6mm\) 或低于 \(99.4mm\) 的概率是非常低的，约为 \(0.00034\%\) 。如果工厂每天生产 100 万个零件，只允许有 \(3.4\) 个零件会超出 \(6 \sigma\) 的范围，几乎可以忽略不计。因此，生产过程是极其稳定和可靠的，达到了六西格玛水平。

那么 \(6 \sigma\) 中 \(3.4ppm\) 的不良率来自哪里呢？

学过中学数学都知道，在正态分布( \(\text{Normal Distribution}\) ) ² 中， \(68.27\%\) 的数据位于平均值的一个标准差内， \(95.45\%\) 位于两个标准差内， \(99.73\%\) 位于三个标准差内，这也是著名的 68-95-99.7 Rule ³ ，如下图 1 所示：

什么是正态分布？

数据可以用不同的方式“分布”，比如数据可以向左散布的多一些，也可以向右散布的多一些，或者分布的乱七八糟，如下图 2 - 4 所示，

但数据经常会集中在一个中心值的附近，而不向左或右偏斜，像一个钟形，如下图 5 所示。

正态分布，又称高斯分布（ \(\text{Gaussian Distribution}\) ），是一种重要的概率分布，数学王子高斯 ⁴ 在正态分布的研究和应用上做出了巨大贡献。有很多日常现象都符合这种分布，如人的身高、考试成绩等。正因为它几乎无处不在，所以叫 \(\text{Normal Distribution}\) 。德国曾经发行的一款 10 马克的纸币上就印着高斯和正态分布曲线，如下图 6 所示。

这个曲线的数学公式大家在中学里都早已见过，如下所示：

\[ f(x; \mu ,\sigma ) = {\frac {1}{\sigma {\sqrt {2\pi }}}}\;e^{-{\frac {\left(x - \mu \right)^{2}}{2 \sigma ^{2}}}} \tag{1} \label{1} \]

\(\mu\) 是正态分布的数学期望值，可解释为位置参数，决定了分布的位置，表示曲线中心在哪里；方差 \(\sigma ^{2}\) 为尺度参数，决定了分布的幅度，表示曲线的扁平情况。均值 \(\mu\) 和方差 \(\sigma ^{2}\) 不同，曲线形状也有所不同，如下图 7 所示：

正态分布的公式看起来非常复杂，里面有 \(\pi\) 、 \(e\) 、 \(\mu\) 、 \(\sigma\) ，组合起来非常复杂。在学习时，课本介绍正态分布时就直接给出这个公式，却从来不说明这个概率密度函数是怎么推导来的，来龙去脉是什么。最近看了 3Blue1Brown 关于 概率论的系列视频 ，我知道了正态分布曲线公式为什么是这样，我们将在下一章节中推导出这个公式。

正态分布公式怎么来的？

有很多种方法都可以推导出正态分布公式，这里将介绍一种既优雅又直观的推导方式，由天文学家赫歇尔（ \(\text{John Herschel}\) ） ⁵ 在 1850 年给出的。3Blue1Brown 的视频 Why π is in the normal distribution (beyond integral tricks) 中详细介绍了这种方式。不过视频中有一些不够严谨的地方，下面会先介绍视频中的推导方法，然后再介绍严谨的数学分析法。

考虑向一个镖盘投掷飞镖，过镖盘中心作 \(x\) 轴和 \(y\) 轴。每次投掷飞镖都会受到随机因素而偏离目标，故每次飞镖在镖盘的落点 \((x, y)\) 是 \(2\) 维随机变量。

假设满足以下 \(2\) 个条件：

1. 落点的 \(x\) 轴和 \(y\) 轴坐标是相互独立的；
2. 落点的概率密度函数仅与其到原点的距离有关，即分布在空间上具有旋转对称性。

3Blue1Brown Herschel 推导方法

如下图 9 所示，设箭头所示落点区域 \(P\) 的概率密度函数为 \(f_2(x, y)\) ， \(f_2\) 表示有 \(2\) 个输入参数，落点的坐标 \(x\) 和 \(y\) 。

由条件 \(1\) 我们知道每个落点区域 \(P\) 的概率密度，可以表示为 \(x\) 轴方向上概率密度函数与 \(y\) 轴方向上概率密度函数的乘积，每个方向上的概率密度函数只有对应方向上一个参数。

设 \(x\) 轴方向上概率密度函数为 \(g(x)\) ， \(y\) 轴方向上概率密度函数为 \(h(y)\) ，则有以下公式：

\[ f_2(x, y) = g(x) h(y) \tag{2} \label{2} \]

又因为条件 \(2\) 我们知道可以通过旋转对称性，可知 \(x\) 轴和 \(y\) 轴方向上概率密度函数相同，所以 \(g(y) = h(y)\) 。同时落点 \(P\) 距离原点的距离为 \(r = \sqrt{x^2 + y^2}\) ，如下图 10 所示：

则 \(f_2(x, y)\) 概率密度函数表示为半径（距原点的距离）的单变量函数 \(f(r)\) ，则有：

\[ f_2(x, y) = f(r) = f(\sqrt {x^2 + y^2}) \tag{3} \label{3} \]

综合 \(\eqref{2}\) 和 \(\eqref{3}\) 可得：

\[ f(\sqrt {x^2 + y^2}) = g(x) g(y) \tag{4} \label{4} \]

假设 \(x\) 轴上距离原点距离为 \(r\) 的点 \(P\) ，坐标为 \((r, 0)\) ，如下图 11 所示：

点 \(P(r, 0)\) 处的概率密度函数可以写成：

\[ f_2(r, 0) = f(r) = g(r) g(0) \tag{5} \label{5} \]

从上式 \(\eqref{5}\) 可知 的 \(f(r)\) 等于 \(g(r)\) 乘以某个常数 \(g(0)\) ，所以 \(f(r) = C g(r)\) ，其中 \(C\) 为某个常数。

由于最终都需要进行归一化，这里不妨设 \(C = 1\) ，则公式 \(\eqref{4}\) 可以写成：

\[ f(\sqrt {x^2 + y^2}) = f(x) f(y) \tag{6} \label{6} \]

至此我们得到了最重要的关系式，问题转变为如何求解函数 \(f\) 。公式 \(\eqref{6}\) 是一个函数方程 ⁶ ，熟悉函数方程的同学可能一眼就知道满足公式 \(\eqref{6}\) 的函数解是指数函数 ⁷ 。

\[ f({r}) = e^{-{r}^2} \tag{7} \label{7} \]

如果我们不知道这个函数方程的解怎么办呢？下面我们就来求解下。

设函数 \(h(x) = f(\sqrt {x})\) ，则有 \(h(x^2) = f(x)\) ，那么公式 \(\eqref{6}\) 可以写成：

\[ h(x^2 + y^2) = h(x^2) h(y^2) \tag{8} \label{8} \]

从公式 \(\eqref{8}\) 我们知道任取 \(2\) 个正整数，先相加再代入函数 \(h\) ，结果等于先分别代入函数 \(h\) 再相乘。这意味着，函数 \(h\) 把加法变成了乘法。

实际上到了这里，还记得中学时学过的指数函数吗？

对于以某个常数 \(c\) 为基数的指数函数 \(h(x) = c^x\) ，可以将加法变成乘法。

但这个只是我们的猜测，我们还需要严谨的证明我们的猜测。

如果这个性质适用于任意 \(2\) 个正整数，那么很容易扩展至任意 \(n\) 个正整数：

\[ h(x_1 + x_2 + \cdots + x_n) = h(x_1)h(x_2) \cdots h(x_n) \tag{9} \label{9} \]

不妨代入 \(x = 5\) ，则有：

\[ \begin{aligned} h(5) & = h(1 + 1 + 1 + 1 + 1) \\ & = h(1)h(1)h(1)h(1)h(1) \\ & = h(1)^5 \end{aligned} \]

这里 \(x\) 可以为任意整数 \(n\) ，则有：

\[ h(n) = h(1 + \cdots + 1) = h(1) \cdots h(1) = h(1)^n \tag{10} \label{10} \]

由公式 \(\eqref{6}\) ，我们知道：

\[ f(0) = f(0)f(0) \]

因为 \(f(0) \ne 0\) ，因此 \(f(0) = 1\) ，那么有：

\[ h(0) = h(0)h(0) = h(0)^2 \]

对于 \(p \in \mathbb{Z}^+\) ，有：

\[ 1 = h(0) = h(p - p) = h(p)h(-p) = h(1)^p h(-p) \]

易得对于任意负整数 \(-p\) ，有：

\[ h(-p) = h(p)^{-1} = h(1)^{-p} \]

这样我们就将函数 \(h\) 推广至整数域 \(\mathbb{Z}\) 。

同理，对于任意正整数 \(q \in \mathbb{Z}^+\) ，可以得到以下表达式

\[ h(q) = \prod_{k=1}^q h(1) = h(1)^q \]

对公式适当变形，可得：

\[ h(1) = h(\frac {q}{q}) = h(\frac {1}{q})^q \]

可以得到：

\[ h(\frac {1}{q}) = h(1)^{\frac{1}{q}} \]

综合上述可得，对于任意有理数 \(\dfrac{p}{q} \in \mathbb{Q}\) ，都满足：

\[ h(\frac{p}{q}) = h(1)^\frac{p}{q} \]

同时由于概率密度函数连续，那么对于无理数也成立。则根据上述推导我们知道公式 \(\eqref{10}\) 在实数域 \(\mathbb{R}\) 都成立。

不妨令 \(a = h(1)\) ，则我们得到了最终表达式：

\[ h(x) = a^x \tag{11} \label{11} \]

为了方便和统一，任意指数函数都可以写成以 \(e\) 为底的指数函数：

\[ h(x) = e^{\ln a x} = e^{c x} \tag{12} \label{12} \]

至此我们就求出了满足函数方程 \(\eqref{6}\) 函数 \(f\) 为：

\[ f(x) = e^{cx^2} \tag{13} \label{13} \]

那么常数 \(c\) 是多少呢？根据概率论我们知道：

\[ \int f(x) \mathrm{d}x = \int e^{cx^2} \mathrm{d}x = 1 \tag{14} \label{14} \]

上式就是大名鼎鼎的高斯积分（ \(\text{Gaussian Integral}\) ） ⁸ ，如何求解高斯积分可以参考这篇文章： 数学之美：几何视角下的高斯积分(Gaussian Integral) 。

最终我们得出了 \(2\) 维情况下的未归一化的概率密度函数：

\[ f_2(x, y) = f_1(x)f_1(y) = e^{-x^2} e^{-y^2} = e^{-(x^2 + y^2)} \tag{15} \label{15} \]

更严谨的数学分析法

上一节我们使用了不那么严谨的方法得到了正态分布的概率密度函数，下面我们使用另外一种方法求出正态分布的概率密度函数。

由落点分布分布在空间上具有旋转对称性，我们可知 \(x\) 轴和 \(y\) 轴具有相同且连续的概率密度函数。

设落点 \(P\) 的概率密度函数为 \(\rho (x, y)\) ， \(x\) 轴方向上概率密度函数为 \(f(x)\) ，则 \(y\) 轴方向上的概率密度函数为 \(f(y)\) ，那么考虑如下图 12 所示的一个充分小的黄色区域 \(\Box ABCD\) ：

飞镖落在黄色区域 \(\Box ABCD\) 的概率为：

\[ \rho (x, y) \mathrm{d}x \mathrm{d}y = f(x) \mathrm{d}x \cdot f(y) \mathrm{d}y \tag{16} \label{16} \]

将等式左边转换为极坐标形式，

\[ \begin{cases} x = r \cos \theta \\ y = r \sin \theta \end{cases} \]

在极坐标下的概率密度函数设为 \(g(r, \theta )\) , 则有：

\[ \rho (x, y) = \rho (r \cos \theta , r \sin \theta) = g(r, \theta ) \]

由条件 \(2\) , \(g(r, \theta )\) 具有旋转对称性，也就是和 \(\theta\) 无关，所以

\[ \frac{\mathrm{d} g(r, \theta)}{\mathrm{d} \theta} = 0 \]

对公式 \(\eqref{16}\) 两边对 \(\theta\) 求导，可得：

\[ \frac{\mathrm{d} f(x)}{\mathrm{d} \theta } f(y) + f(x) \frac{\mathrm{d} f(y)}{\mathrm{d} \theta } = 0 \]

利用链式法则，有：

\[ -r \sin \theta f'(x)f(y) + f(x) f'(y) r \cos \theta = 0 \]

上式移项可得：

\[ \frac{f'(x)}{f(x) x} = \frac{f'(y)}{f(y) y} = C \]

我们则有：

\[ \frac{f'(x)}{f(x)} = C x \]

对上式进行积分，可得：

\[ \int \frac{f'(x)}{f(x)} \mathrm{d}x = \int C x \mathrm{d}x \]

求解上式可得：

\[ \ln f(x) = \frac{C}{2}x^2 + C' \]

则求得函数 \(f(x)\) 为：

\[ f(x) = A e^{\frac{1}{2}Cx^2} \tag{17} \label{17} \]

同理 \(f(y)\) 为：

\[ f(y) = A e^{\frac{1}{2}Cy^2} \tag{18} \label{18} \]

由概率论我们知道 \(C < 0\) ，同时 \(\int_{0}^{\infty } f(x) \mathrm{d}x = \dfrac{1}{2}\) ，则：

\[ \int_{0}^{\infty } e^{\frac{C}{2} x^2} \mathrm{d}x = \frac{1}{2A} \]

对高斯核函数进行升维，对于 \(2\) 维正态分布，从笛卡尔坐标系转换为极坐标系则有：

\[ \begin{aligned} \int_{0}^{\infty} \int_{0}^{\infty} e^{\frac{C}{2}(x^2 + y^2)} \mathrm{d}x \mathrm{d}y & = \int_{0}^{\infty} \int_{0}^{\infty} e^{\frac{C}{2}r^2} r \mathrm{d}r \mathrm{d}\theta \\ & = \int_{0}^{\frac{\pi }{2}} \int_{0}^{\infty} e^{\frac{C}{2}r^2} r \mathrm{d}r \mathrm{d} \theta \\ & = \frac{1}{4A^2} \end{aligned} \]

令 \(u = r^2\) ，则有：

\[ \int_{0}^{\infty} e^{\frac{C}{2}r^2} r \mathrm{d}r = \int_{0}^{\infty} e^{\frac{C}{2}u} \frac{\mathrm{d}u}{2} = \frac{1}{2} \left [ \frac{2}{C} e^{\frac{C}{2} u} \right ]_{0}^{\infty } = -\frac{1}{C} \]

可得：

\[ -\frac{1}{C} \int_{0}^{\frac{\pi }{2}} \mathrm{d}\theta = -\frac{1}{C} \frac{\pi}{2} = \frac{1}{4A^2} \]

所以可求得：

\[ A = \sqrt{\frac {-C}{2 \pi}} \]

代入公式 \(\eqref{17}\) ，即可得 \(f(x)\) 表达式为：

\[ f(x) = A e^{\frac{1}{2}Cx^2} = \sqrt{\frac {-C}{2 \pi}} e^{\frac{1}{2}Cx^2} \tag{19} \label{19} \]

至此，我们只剩下一个未知参数 \(C\) 就得到所求概率密度函数。

根据期望 \(\mu\) 和方差 \(\sigma^2\) 定义，则有：

\[ \mu = \int_{-\infty }^{\infty } x f(x) \mathrm{d}x \tag{20} \label{20} \]

\[ \sigma ^2 = \int_{-\infty }^{\infty }x^2 f(x) \mathrm{d}x \tag{21} \label{21} \]

由公式 \(\eqref{17}\) 已经得到的 \(f(x)\) 表达式代入公式 \(\eqref{20}\) ，易得期望 \(\mu\) 表达式为：

\[ \mu = \int_{-\infty }^{\infty }x \sqrt {\frac {-C}{2\pi}} e^{\frac{1}{2}Cx^2} \mathrm{d}x \tag{22} \label{22} \]

分析公式 \(\eqref{22}\) ， \(f(x)\) 为偶函数， \(x\) 为奇函数，所以上式 \(\mu = 0\) 。

将 \(f(x)\) 表达式代入方差 \(\sigma^2\) 公式，则有：

\[ \sigma ^2 = \int_{-\infty }^{\infty }x^2 f(x) \mathrm{d}x = 2 \sqrt{\frac {-C}{2 \pi}} \int_{0}^{\infty }x^2 e^{\frac{C}{2} x^2} \mathrm{d}x \tag{23} \label{23} \]

令 \(u = x\) ， \(v = \dfrac{1}{C} e^{\frac{C}{2} x^2}\) ，则有：

\[ \mathrm{d}v = x e^{\frac{C}{2} x^2} \]

根据分部积分公式：

\[ \int u \frac{\mathrm{d}v}{\mathrm{d}x}\,\mathrm{d}x = uv - \int \frac{\mathrm{d}u}{\mathrm{d}x}v\,\mathrm{d}x \]

则对公式 \(\eqref{19}\) 进行分部积分求解可得：

\[ \begin{aligned} \int_{0}^{\infty } u \frac{\mathrm{d}v}{\mathrm{d}x}\,\mathrm{d}x & = \lim_{x \to \infty} \frac{x}{C} e^{\frac{C}{2} x^2} - \frac{0}{C} e^{\frac{C}{2} 0} - \frac{1}{C} \int_{0}^{\infty } e^{\frac{C}{2} x^2} \mathrm{d}x \\ & = 0 - 0 - \frac{1}{C} \frac{\sqrt {2 \pi }}{2 \sqrt {-C}} \\ & = - \frac{1}{C} \frac{\sqrt {2 \pi }}{2 \sqrt {-C}} \end{aligned} \]

所以我们求得 \(\sigma^2\) 为：

\[ \sigma^2 = 2 \sqrt{\frac {-C}{2 \pi}} \frac{-1}{C} \frac{\sqrt {2 \pi }}{2 \sqrt {-C}} = -\frac{1}{C} \]

即：

\[ C = -\frac{1}{\sigma^2} \]

将求得 \(A\) 和 \(C\) 代入公式 \(\eqref{17}\) ，最终我们求得概率密度函数为：

\[ f(x) = {\frac {1}{\sigma {\sqrt {2\pi }}}}\;e^{-{\frac {x^{2}}{2\sigma ^{2}}}} \tag{24} \label{24} \]

公式 \(\eqref{24}\) 是期望 \(\mu = 0\) 的特殊情况。

而当期望 \(\mu \ne 0\) 时，将 \(x\) 平移 \(\mu\) ，则为 \(x - \mu\) ，所以更一般的公式为：

\[ f(x) = {\frac {1}{\sigma {\sqrt {2\pi }}}}\;e^{-{\frac {\left(x - \mu \right)^{2}}{2\sigma ^{2}}}} \tag{25} \label{25} \]

特别地，当 \(\mu = 0\) 且 \(\sigma = 1\) ，这个分布被称为标准正态分布：

\[ f(x)={\frac {1}{\sqrt {2\pi }}}\, e^{-{\frac {x^{2}}{2}}} \tag{26} \label{26} \]

正态分布公式的几何意义

通过 \(\text{Herschel}\) 给出的优雅直观方法，仅仅依靠那 \(2\) 个假设条件，我们居然最终求出了正态分布的公式。有没有感觉到数学的美感？

最初看到 3Blue1Brown 的这个视频，感觉非常美，正态分布那么复杂的公式居然有这么优雅直观的方式自然而然的出来了！

分析正态分布公式，公式中的 \(\pi\) 意味着空间上的对称性，即点分布距离中心是对称的。而 \(e\) 的出现意味着取了时间上的极限，而这和中心极限定理（ \(\text{Central limit theorem}\) ） ⁹ 有关，我会在下一篇文章详细解释，敬请期待！

参考文献