概率统计随机过程之如何推导得到正态分布—正态分布的理解角度

概率统计随机过程之如何推导得到正态分布————正态分布的理解角度

原文:Kim, Kiseon, 和Georgy Shevlyakov. 《Why Gaussianity?》 IEEE Signal Processing Magazine 25, 期 2 (2008年3月): 102–13. https://doi.org/10.1109/MSP.2007.913700.

中文部分靳志辉正态分布章节推导完善。

正态分布表达式形式的出现

正态分布是概率论和数理统计中的一个重要的分布。在最开始的工作中(1733年),法国数学家棣莫佛为了求二项分布在数值较大时的近似计算值,推导出正态分布的表达形式。但是,棣莫佛并没有把它当作是一个概率分布,只是认为这是二项分布在\(N→∞\)时的一个近似表达式。40年后,法国另一个数学家拉普拉斯在推导最初版本的中心极限定理时也得出了正态分布表达式的形式,但当时他也没有意识这是一个概率的分布,也把它当作单纯的一种数学表达式。但是,这是第一次正态密度函数被数学家刻画出来,而且是以二项分布的极限分布的形式被推导出来的。熟悉基础概率统计的人都知道这个结果其实叫棣莫弗-拉普拉斯中心极限定理。当然当时还没有这个称呼,直到1920年才有数学家波利亚提出中心极限定理这个名称。

棣莫弗-拉普拉斯中心极限定理:若\(X\thicksim B(n,p)\)是n次伯努利试验中出现事件\(A\)的次数,每次试验成功的概率为\(p\),且\(q=1-p\),则对任意有限区间\([a,b]\),令\(x_k=\frac{k-np}{\sqrt{npq}}\),当\(n→∞\)时:

  1. \(P(X=k)→\frac{1}{\sqrt{npq}}\cdot\frac{1}{\sqrt{2\pi}}e^{-{1\over2}x^2_k}\)
  2. \(P(a≤\frac{k-np}{\sqrt{npq}}≤b)→\int_a^b\varphi(x)\mathrm{d}x\),其中\(\varphi(x)=\frac{1}{\sqrt{2\pi}}e^{-\frac{x^2}{2}}(-∞<x<∞)\)

棣莫佛-拉普拉斯中心极限定理指出:参数\(n,p\)的二项分布以\(np\)为均值,\(np(1-p)\)为方差的正态分布为极限。当\(n\)很大时,可以用作二项分布的近似值。

总结一句:这个阶段正态分布表达式只是作为一个函数近似的结果,并没有被认为是一种概率分布。

高斯——误差与正态分布

高斯所在的18-19世纪是天文学迅速发展的时期,积累了大量观测数据。但是,这些数据或多或少都是存在误差的。如何处理这些数据得到最准确的结果,或者这些误差有着什么样的数学规律?在当时,这是一个亟待解决的问题。

辛普森和拉普拉斯认为误差应该服从这样一些普遍规律:

  1. 误差分布是对称的,即\(f(x)=f(-x)\)
  2. 小误差应该比大误差更常见,即\(f(x)\)\(x\)的减函数。且误差应该随\(|x|\)增大而逐渐趋于0。
  3. 每次测量产生的误差应该是独立的。
  4. 误差分布函数应是连续函数。

为此,辛普森和拉普拉斯分别设计了三角概率分布函数和重指数分布(拉普拉斯分布)作为误差的分布函数。并用这些分布证明观测值取平均数比恰当选择观测值更能接近真实值。拉普拉斯在重指数分布基础上,希望使用贝叶斯法(当时拉普拉斯称之为不充分推理原则)逐渐修正后验概率得到的误差分布。但是,拉普拉斯的工作到此基本进展就很小了,在18世纪70-80年代沿着这条路径磕磕绊绊走了十几年,依旧进展甚微。

1809年,高斯发表了其数学和天体力学的名著《绕日天体运动的理论》。在此书末尾,他写了一节有关“数据结合”的问题,实际涉及的就是这个误差分布的问题。

高斯采用了拉普拉斯设计的函数\(L(\theta)\),设真值为\(\theta\),n个独立的测量值\(X_1,X_2,\dotsb,X_n\)\(L(\theta)\)为出现这一现象的联合概率密度: \[\begin{aligned} L(\theta)=&L(\theta;X_1,X_2,\dotsb,X_n)\\ \stackrel{\text{独立性}}{=}&f(X_1-\theta)f(X_2-\theta)\dotsb f(X_n-\theta)\tag{1} \end{aligned}\] 其中,\(f(\cdot)\)待定的误差密度函数。当时接下来的处理高斯与拉普拉斯完全不同。

一是高斯没有采取贝叶斯的推理方式,而是径直让\((1)\)达到最大值的\(\hat{\theta}=\hat{\theta}(X_1,\dotsb,X_n)\)作为\(\theta\)的估计。即 \[L(\hat{\theta})=\max_{\theta} L(\theta)\tag{2}\] 当然现在我们把这种方法叫做最大似然估计(Maximum Likelihood Estimation, MLE),\(L(\theta)\)叫做似然函数(正式方法由Fisher推广于1922年)。最大似然估计的思想其实就是很直接:最合理的事情(即真值)最有可能发生。举个简化的例子,如果在一个不透明的袋子里有10个球,黑白两种颜色,但是不知道各有几个;我们袋子中有放回地取100次,其中71次是白球,那么大家认为白球、黑球各有多少个?一般人都会猜7个白球,3个黑球吧。但如果是黑白球各有5个,会不会抽样出上述结果呢?当然有可能啊,只不过概率没有7个白球,3个黑球大。这就是我们在不知不觉中直观地使用了最大似然估计。

二是高斯把辛普森和拉普拉斯的思路倒了过来。高斯先认可算数平均\(\bar{X}\)是应取的估计值,然后去找误差密度函数\(f(\cdot)\)以迎合这一点,即找这样的\(f(\cdot)\),使由\((2)\)式决定的\(\hat{\theta}\)就是\(\bar{X}\)

高斯发现,使用最大似然估计的情况下,只有当 \[f(x)=\frac{1}{\sqrt{2\pi}h}e^{\frac{-x^2}{2h^2}},h>0\] 时,才有\(\hat{\theta}=\bar{X}\)。注意,高斯在原文中,并没有让\(h=\sigma\),即方差。下面,我们来推导这个过程。

我们从似然函数\((1)\)式开始,要求关于\((1)\)式的最小值,而其定义域\(x\in (-∞,∞)\),由费马引理(皮埃尔·德·费马 1601-1665,法国律师和业余数学家,业余数学家之王)可知,\((1)\)式的极值点必然是驻点,即导数为0的点。但是,连乘表达式的求导难以处理,因此高斯采用了化乘为加,且不改变驻点位置(证明见附1)的取\(\log(ln)\)法,即 \[\frac{\mathrm{d}ln L(\theta)}{\mathrm{d}\theta}=0\] 展开\(\ln L(\theta)\)求导可得: \[\sum_{i=1}^n\frac{f'(x_i-\theta)}{f(x_i-\theta)}=0\] 为了方便计算,我们令\(g(x)=\frac{f'(x)}{f(x)}\),则有 \[\sum_{i=1}^n g(x_i-\theta)=0\] 由于高斯假设极大似然估计的解就是算术平均\(\bar{x}\),把解直接带入上式可得 \[\sum_{i=1}^n g(x_i-\bar{x})=0 \tag{3}\] 下面就是高斯精彩的数学表演,由于\(n,x_i\)的任意性,在\((3)\)式中取\(n=2\)有: \[g(x_1-x)+g(x_2-x)=0\] 又因为\(\bar{x}=\frac{x_1+x_2}{2}\),所以\(x_1-x=-(x_2-x)\),即\(g(x_1-x)+g(-(x_1-x))=0\),由于\(x_i\)的任意性可得 \[g(-x)=-g(x),g(x)是奇函数\]\((3)\)式中再取\(n=m+1\),并要求\(x_1=\dotsb=x_m=-x,x_{m+1}=mx\),则有\(\bar{x}=0\),并且 \[\sum_{i=1}^n g(x_i-\bar{x})=mg(-x)+g(mx)=0\] 又因为\(g(x)\)是奇函数,所以有 \[g(mx)=mg(x)\] 而满足上式的唯一的连续奇函数解就是\(g(x)=cx\)(可用两边求导数证明),从而进一步可以解微分方程: \[\frac{f'(x)}{f(x)}=cx\Rightarrow \frac{1}{y}\frac{dy}{dx}=cx\stackrel{\text{分离变量法}}{\Rightarrow}\frac{1}{y}dy=cx dx\] 两边求补丁积分,且y=f(x)作为概率值恒大于等于0,可得: \[\ln y=\frac{1}{2}cx^2+K\Rightarrow y = e^Ke^{\frac{1}{2}cx^2}\] 我们令\(e^K=M\),且\(f(x)\)作为概率密度函数,必然有\(\int_{-∞}^∞f(x)dx=1\),所以 \[\int_{-∞}^∞f(x) dx= \int_{-∞}^∞ Me^{\frac{1}{2}cx^2}dx =\int_{-∞}^∞\frac{M}{\sqrt{0.5|c|}} e^{-(\sqrt{0.5|c|}x)^2} d(\sqrt{0.5|c|}x)= 1\] 这里需要注意,为了保证积分收敛到1,那么常数\(c<0\)是必然的,在\(e\)的指数部分要加上负号,同时开方时需要加绝对值符号。此外,令\(t=\sqrt{0.5|c|}x\),则t的积分区间为\((-∞,∞)\),上式可转换为 \[\frac{M}{\sqrt{0.5|c|}}\int_{-∞}^∞ e^{-t^2}dt=1\] 而这个广义积分正好有\(\int_{-∞}^∞ e^{-t^2}dt=\sqrt{\pi}\)(即高斯积分证明附2),所以可得 \[M=\frac{\sqrt{|c|}}{\sqrt{2\pi}}\]\(M\)代入\(f(x)\)可得 \[ f(x)=\frac{\sqrt{|c|}}{\sqrt{2\pi}}e^\frac{cx^2}{2}\tag{4} \] 这就是高斯推导的正态分布的表达式,它满足以算数平均作为最大似然估计的结果。如果我们用方差定义式计算一下\(\sigma^2\),可以得到\(\sigma^2=-\frac{1}{c}\)(证明见附3),由此可得,\(f(x)\)就是\(N(0,\sigma^2)\)的表达式。正态分布在误差分布计算上有着极大的理论优势,对后世影响极大,以至于正态分布有了高斯分布的名称。

但是,高斯的解释在逻辑上有点循环论证的味道。由于取算术平均是合理的方式,所以用MLE导出正态分布;又因为误差服从正态分布,所以算术平均能够获得更小的误差值。但是,拉普拉斯很快得知了高斯的工作,并将正态分布与中心极限定理联系起来。在1810年,拉普拉斯提出了元误差学说:如误差可以看成是大量的,有各种原因导致的微小的量的叠加,则根据中心极限定理,误差的分布应服从正态分布。拉普拉斯的补充使误差的正态分布理论有了一个更自然、更令人信服的解释。

赫歇尔和麦克斯韦——正态分布的\(\pi\)从哪里来?

每当公式中出现\(\pi\)的时候,我都不禁会问自己?圆在哪里?而正态分布的表达式中正好有一个\(\frac{1}{\sqrt{2\pi}}\)系数,这个怎么解释呢?我们可以从赫歇尔和麦克斯韦的推导过程中一窥究竟。

1850年,天文学家赫歇尔在对星星位置进行测量时,需要考虑二维误差的分布,为了推导这个误差的概率密度分布\(p(x,y)\),赫歇尔设置了两个准则:

  1. x轴,y轴的误差时相互独立的,即随机误差在正交方向上相互独立
  2. 误差的概率分布在空间上具有旋转对称性,即误差的概率分布和角度没有关系

这两个准则对于进行了大量实际测量的赫歇尔都非常合理,而空间的旋转对称性,在二维空间里,那不正是让我们想到了“圆”吗?(三维空间就是球~)

由第一条准则,我们可以把\(x,y\)联合概率分布,分离成两部分: \[p(x,y)=f(x)\cdot f(y)\] 而面对旋转对称性,和角度相关的,我们不禁就和极坐标联系起来: \[p(x,y)=p(r\cos\theta,r\sin\theta)=g(r,\theta)\] 由第二条准则,旋转不变性可知,函数\(g(r,\theta)\)实际上与角度\(\theta\)无关,即\(g(r,\theta)=g(r)\)。实际上,正是这个旋转不变性,使得\(\theta\)的积分域是\([0,2\pi]\),也就是一个圆,即正态分布中\(\pi\)的来历。综上所述,我们可得: \[f(x)f(y)=g(r)=g(\sqrt{x^2+y^2})\tag{5}\] 由于\(x,y\)的任意性,我们取\(y=0\),得到\(g(x)=f(x)f(0)\),所以上式中 \[g(\sqrt{x^2+y^2})=f(\sqrt{x^2+y^2})f(0)\tag{6}\] 联合式\((5)(6)\)并在两边同时除以\(f^2(0)\)有: \[\frac{f(x)f(y)}{f(0)f(0)}=\frac{f(\sqrt{x^2+y^2})f(0)}{f(0)f(0)}\Rightarrow \frac{f(x)}{f(0)}\frac{f(y)}{f(0)}=\frac{f(\sqrt{x^2+y^2})}{f(0)}\] 两边取\(\log\)则有 \[\log\left[\frac{f(x)}{f(0)}\right]+\log\left[\frac{f(y)}{f(0)}\right]=\log\left[\frac{f(\sqrt{x^2+y^2})}{f(0)}\right]\]\(h(x)=\log[\frac{f(x)}{f(0)}]\),则上式可简化为函数方程: \[h(x)+h(y)=h(\sqrt{x^2+y^2})\tag{7}\] 函数方程与代数方程、微分方程不同,并没有通用的解法。所以这个分支也没能发展起来。但是,我们不难猜出\(h(x)=ax^2\)是符合\((7)\)的一个连续函数解(我觉得可以用此函数方程的二阶导数为常数和\(h(0)=0\)证明其解的唯一性)。同时,由于概率积分的收敛性,需要\(a<0\),不妨令\(\alpha=-a\)\[h(x)=ax^2\Rightarrow f(x)=f(0)e^{-\alpha x^2}\] 对于概率密度函数定义域积分为1,则有 \[\int_{-\infty}^\infty f(0)e^{-\alpha x^2}dx = 1\Rightarrow \frac{f(0)}{\sqrt{\alpha}}\int_{-\infty}^\infty e^{-(\sqrt{\alpha} x)^2}d\sqrt{\alpha}x=1\] 不妨令\(\sqrt{\alpha}x = t, t\in (-\infty,\infty)\),根据附2的高斯积分\(\int_{-∞}^∞ e^{-x^2}dx=\sqrt{\pi}\),可推得\(f(0)=\sqrt{\frac{\alpha}{\pi}}\),综上所述有: \[f(x)=\sqrt{\frac{\alpha}{\pi}}e^{-\alpha x^2}\] 这也是正态分布的表达形式。这也是在0均值正态分布,当\(\sigma^2=\frac{1}{2\alpha}\)的结果(方差证明的方式类似附3)。 而满足准则1,2的二维表达式,就是两个独立同分布\(f(x)\)的乘机,即\(p(x,y)\)为二维i.i.d正态分布的密度函数: \[p(x,y)=\frac{\alpha}{\pi}e^{-\alpha(x^2+y^2)}\]

1860年麦克斯韦在考虑气体分子运动速度分布的时候,基于类似的准则推出了著名的气体分子运动速率分布的麦克斯韦-玻尔兹曼气体速率分布定率: \[ F(v)=(\frac{m}{2\pi kT})^{3/2}e^{-\frac{mv^2}{2kT}} \] 这拆开看其实就是三维空间的i.i.d正态分布嘛,其中\(\sigma^2=\frac{kT}{m}\) \[\begin{aligned} F(v)&=(\frac{m}{2\pi kT})^{3/2}e^{-\frac{mv^2}{2kT}}\\ &=(\frac{m}{2\pi kT})^{1/2}e^{-\frac{mv_x^2}{2kT}}\times(\frac{m}{2\pi kT})^{1/2}e^{-\frac{mv_y^2}{2kT}}\times(\frac{m}{2\pi kT})^{1/2}e^{-\frac{mv_z^2}{2kT}} \end{aligned}\] 其中,\(v^2=v_x^2+v_y^2+v_z^2,v_x,v_y,v_z\)是三维空间正交的三个方向向量。赫歇尔和麦克斯韦的工作神奇之处在于,他们没有利用任何概率论只是,只是基于空间不变性,就推导出了正态分布。

为什么正态分布那么普遍?——稳定性,兰登的推导

1941年电子工程师兰登通过分析经验数据他发现噪声电压的分布模式很相似,不同的是分布的层级,而这个层级可以使用方差\(σ^2\)来刻画。基于这些经验,兰登提出以下两个准则:

  1. 随机噪声具有稳定的分布模式。
  2. 累加一个微小的随机噪声,不改变其稳定的分布模式,只改变分布的层级。(用方差度量)
  3. 微小的随机噪声关于0对称分布。(本质上是相位(\([0-2\pi]\))的随机性,也是正态分布中\(\pi\)的来源)

用数学的语言描述: 如果 \[X \sim p(x;\sigma^2),\epsilon\sim q(e),X'=X+\epsilon\] 则有 \[ X'\sim p(x;\sigma^2+\text{var}(\epsilon)) \] 按照两个随机变量和的分布的计算方式,\(X′\)的分布密度函数将是\(X\)的分布密度函数和\(\epsilon\)的分布密度函数的卷积,即有 \[ f(x')=\int_{-\infty}^\infty p(x'-e;\sigma^2)q(e)de \] 通过泰勒级数展开和解二阶偏微分方程(也是著名的扩散方程),我们可以解得 \[p(x;\sigma^2)=\frac{1}{\sqrt{2\pi\sigma^2}}e^{-\frac{-x^2}{2\sigma^2}} \] (这里详情可以扩展写)

杰恩斯(概率论沉思录作者)指出这个推导这基本上就是中心极限定理的增量式版本,相比于中心极限定理是一次性累加所有的因素,兰登的推导是每次在原有的分布上去累加一个微小的扰动。而在这个推导中,我们看到,正态分布具有相当好的稳定性;只要数据中正态的模式已经形成,他就容易继续保持正态分布,无论外部累加的随机噪声 q(e) 是什么分布,正态分布就像一个黑洞一样把这个累加噪声吃掉。这也是稳定性也是正态分布普遍存在的原因。

基于最大熵的推导

如果给定一个分布密度函数\(p(x)\)的均值\(µ\)和方差\(σ^2\)(给定均值和方差这个条件,也可以描述为给定一阶原点矩和二阶原点矩,这两个条件是等价的),则在所有满足这两个限制的概率分布中,熵最大的概率分布\(p(x|µ; σ^2)\)就是正态分布\(N(µ; σ^2)\)

此部分的证明用到了古典泛函分析的变分法。

附录

附1 似然函数取对数

命题:对于定义域\([x_a,x_b]\)上的正函数\(f(x)\),取对数不影响驻点位置。

证明: 驻点要求函数导数为0,即\(f'(x)=0\);而\((\log_a f(x))'=\frac{f'(x)}{f(x)\ln a}\)。由于\(f(x)>0\),且\(a \neq 1\),所以只要当\(f'(x)=0\)时,\((\log_a f(x))'=0\),即二者驻点位置一致。

附2 高斯积分

命题:证明高斯积分\(I=\int_{-∞}^∞ e^{-x^2}dx=\sqrt{\pi}\)

证明:令两个独立的变量\(x,y\) \[ I^2=\int_{-∞}^∞ e^{-x^2}dx\times \int_{-∞}^∞ e^{-y^2}dy=\int_{-∞}^∞\int_{-∞}^∞ e^{-(x^2+y^2)} dx dy \] 对于含有\(x^2+y^2\)的形式,我们一般使用极坐标变换:\(x=r\cos\theta,y=r\sin\theta\);元面积变换:\(dxdy=rdrd\theta\) \[ \begin{aligned} I^2 &= \int_0^{2\pi}\int_0^∞ e^{-r^2}rdrd\theta=\int_0^{2\pi}d\theta\int_0^∞ e^{-r^2}rdr \\ &=2\pi\times (-\frac{1}{2}e^{-r^2}\bigg|_0^∞)=2\pi\times{1\over2}=\pi\\ &\Rightarrow I=\sqrt{\pi} \end{aligned} \] 得证。

补充:高斯积分与Gamma函数(第二类欧拉积分)的关系。在\(I\)中,因为\(e^{-x^2}\)是一个偶函数,因此积分域可以从\(\int_{-∞}^∞dx\)变成\(2\cdot\int_0^∞dx\),此时\(x\in(0,∞)\)。如果我们做一个积分变量替换\(x=\sqrt{t},dx={1\over 2}t^{-{1\over 2}}dt\),则\(t\in(0,∞)\),且\(I\)有如下变换: \[ I= \sqrt{\pi}=\int_{-∞}^∞ e^{-x^2}dx =2\int_0^∞ e^{-t}\cdot{1\over 2}t^{-{1\over 2}}dt=\int_0^∞ e^{-t}t^{-{1\over 2}}dt=\Gamma({1\over2}) \]\[I=\underbrace{\int_{-∞}^∞ e^{-x^2}dx}_{\text{高斯积分}}=\underbrace{\int_0^∞ e^{-x}x^{-{1\over 2}}dx}_{\text{第二类欧拉积分}}=\Gamma({1\over2})=\sqrt{\pi}\]

附3 高斯在推导正态分布时系数\(c\)和方差\(\sigma^2\)的关系

命题:式\((4)\)中的系数\(c=-\frac{1}{\sigma^2}\)

证明:首先,从\((4)\)中,我们不难观察出分布函数是一个偶函数,所以有\(E[x]=\mu=0\),则方差\(E[(x-\mu)^2]=E[x^2]\),则有: \[ \begin{aligned} \sigma^2=&\frac{\sqrt{|c|}}{\sqrt{2\pi}}\int_∞^∞ x^2e^\frac{cx^2}{2}dx\\ \stackrel{\text{分部积分}}{=}&\frac{\sqrt{|c|}}{\sqrt{2\pi}}\left\{x\cdot\frac{1}{c}e^\frac{cx^2}{2}\bigg |_{-\infty}^\infty-\int_∞^∞ 1\cdot \frac{1}{c}e^\frac{cx^2}{2}dx\right\}\\ \stackrel{\text{等价无穷小}}{=}&\frac{\sqrt{|c|}}{\sqrt{2\pi}}\left\{0-\int_∞^∞ 1\cdot \frac{1}{c}e^\frac{cx^2}{2}dx\right\}\\ \stackrel{\text{第一类积分变量替换}}{=}&\frac{\sqrt{|c|}}{\sqrt{2\pi}}\left\{-\int_∞^∞ \frac{1}{c}\cdot\sqrt{\frac{2}{|c|}} e^{-(\sqrt{\frac{|c|}{2}}x)^2}d(\sqrt{\frac{|c|}{2}}x)\right\}\\ \stackrel{t=\sqrt{\frac{|c|}{2}}x}{=}&\frac{-1}{\sqrt{\pi}c}\left\{\int_∞^∞ e^{-t^2}dt\right\}\\ \stackrel{\text{附录2高斯积分}}{=}&\frac{-1}{\sqrt{\pi}c}\cdot \sqrt{\pi}=\frac{-1}{c} \end{aligned} \] 即系数\(c=-\frac{1}{\sigma^2}\),代入\((4)\)可得正态分布\(N(0,\sigma^2)\)