02 概率论

2.1 前言

在前面章节中,我们已经了解了概率论在机器学习中所扮演的重要角色。本章将讨论更多关于概率论的细节。我们并没有足够的篇幅展开相关领域的深层次讨论 —— 读者可以自行参考更多的相关书籍。但在后面的章节中,将简明扼要的介绍许多可能会用到的关键思想。

在探讨更多技术细节之前,请思考一个问题:什么是概率?我们对诸如「一个硬币面朝上的概率为 0.5」的表述已经非常熟悉。但这句话到底意味着什么?

关于概率至少有两种不同的解释。一种是 频率 (frequentist) 学派 的解释。在这种观点中,概率代表事件在长时间实验情况下出现的频率。比如,在前面例子中,是指如果投掷一枚硬币很多次,那么我们相信有一半的次数硬币正面朝上。

另一种为 贝叶斯 ( Bayesian ) 学派 的解释。在这种观点中,概率是用来量化人们对某些事件的 不确定性 ( uncertainty ) ,所以本质上与信息而非重复的实验相关。从贝叶斯学派观点看待上述例子,意味着我们相信在下一次投掷硬币时,硬币正面朝上的可能性为 0.5。

贝叶斯解释的一个重要优势在于,它可以用来衡量那些无法进行重复试验的事件的不确定度。比如说估计到 2020 年冰川融化的概率。该事件本身可能只发生一次甚至不会发生,也就是说它是不能被重复的试验。然而,我们可以量化针对该事件发生的不确定度 ( 比如说基于采取的一些抑制全球变暖的行为,可以认为该事件发生的可能性会变小 ) 。再比如第 1 章中提及的垃圾邮件分类任务,我们可能已经收到一个特定邮件信息,希望计算它是垃圾邮件的可能性。或者在雷达屏幕上观察到一个移动光点,我们希望计算这个飞行物身份的概率分布 ( 它是一只鸟还是一个飞机呢? ) 。在上述所有案例中,尽管没有一个事件是可以重复试验的,但贝叶斯观点却是有效且具备可解释性的。所以在本书中我们将采用贝叶斯观点对概率的解释。 不过幸运的是,无论采取哪种观点看待概率,概率论的基本规则都一样。

2.2 关于概率论的简单综述

本节介绍概率论基础知识,仅针对那些对相关知识已经生疏的读者。关于更多相关细节,可以参考其他书籍。已对这块知识比较熟悉的读者可直接跳过本节。

2.2.1 离散型随机变量

表达式 \(p ( A )\) 表示事件 \(A\) 发生的概率。比如 \(A\) 可能表示 “明天会下雨”。其中 \(p ( A )\) 满足 \(0 \le p ( A ) \le1\) ,如果 \(p ( A ) =0\) ,表示事件 \(A\) 不可能发生, \(p ( A ) =1\) 意味着事件 \(A\) 肯定发生。我们使用 \(p ( \bar A )\) 表示事件非 \(A\) 发生的可能性,满足 \(p ( \bar A ) =1-p ( A )\) 。通常将 “\(A\) 发生” 这个事件记为 \(A=1\) ,”\(A\) 不发生” 记为 \(A=0\)

通过定义 离散型随机变量 ( discrete random variable ) \(X\) ,可以扩展出二元事件 ( 即事件只存在两种状态 ) 的符号表达,该离散型变量取值于一个有限集或者可数无限集 \(\chi\)( 译者注:关于可数无限集的例子:比如做一个抛掷硬币的试验,直到第一次出现正面时抛掷硬币的次数 \(\chi\) 的取值所构成的就是一个可数无限集 ) 。

我们将事件 \(X=x\) 发生的概率表示为 \(p ( X=x )\) ,或者直接写成 \(p ( x )\) 。其中符号 \(p ( )\) 称为 概率质量函数 ( probability mass function ) , 满足性质 \(0 \le p ( x ) \le 1,\sum_{x\in\chi} p ( x ) =1\) 2.1 展示了定义在一个有限状态空间 ( state space ) \(\chi = \{1,2,3,4\}\) 上的两种概率质量函数。其中左图属于均匀分布, \(p ( x ) =1/4\) ,右图为一个退化分布 \(p ( x ) =\mathbb {I} ( x=1 )\) ,其中 \(\mathbb {I} ( )\) 为二元指示函数 ( indicator function ) ,该分布意味着 \(X\) 永远等于 1,换句话说,它是一个常数。

图 2.1 (a) 是 \(\{1,2,3,4\}\) 上的均匀分布,其中 \(p(x=k)=1/4\) 。(b) 退化分布, 当 \(x=1\) 时,\(p(x)=1\),当 \(x \in \{2,3,4\}\) 时,\(p(x)=0\)

2.2.2 基本定理

2.2.2.1 两个事件并集发生的概率

给定两个事件 \(A\)\(B\) ,定义事件 \(A\)\(B\) 发生的概率为:

\[\begin{split} \begin {align*} p(A \vee B)&=p(A)+p(B)-p(A \wedge B) ~~\tag {2.1}\label {eqn:2.1}\\\\ &=p(A)+p(B) \text { if } A \text { and } B \text { are mutually exclusive } ~~\tag {2.2}\label {eqn:2.2} \end {align*} \end{split}\]

2.2.2.2 联合概率

定义事件 \(A\)\(B\) 同时发生的概率为:

\[ p ( A,B ) =p ( A \wedge B ) =p ( A|B ) p ( B ) ~~\tag {2.3}\label {eqn:2.3} \]

上式通常又被称为 乘法法则 ( product rule )

给定两个事件的 联合概率分布 ( joint distribution ) \(p ( A,B )\) ,定义 边缘分布 ( marginal distribution ) 如下:

\[ p ( A ) =\sum_b p ( A,B ) =\sum_b p ( A|B=b ) p ( B=b ) ~~\tag {2.4}\label {eqn:2.4} \]

上式针对事件 \(B\) 所有可能的状态进行求和(积分)。类似地,也可以定义 \(p ( B )\) 。 该式被称为 求和法则 ( sum rule ) 或者叫 全概率法则 ( rule of total probability )

可以多次使用乘法法则,进而引出概率论中的 链式法则 ( chain rule )

\[ p ( X_{1: D} ) =p ( X_1 ) p ( X_2|X_1 ) p ( X_3|X_2,X_1 ) ...p ( X_D|X_{1: D-1} ) ~~\tag {2.5}\label {eqn:2.5} \]

式中模仿了 Matlab 中的一种符号写法 \(1: D\) 表示集合 \(\{1,2,…,D\}\)

2.2.2.3 条件概率

我们定义在事件 \(B\) 发生的前提下,事件 \(A\) 发生的概率为 条件概率 ( conditional probability )

\[ p ( A|B ) =\frac {p ( A,B )}{p ( B )} \ {\rm {if}} \ p ( B ) \gt 0 ~~\tag {2.6}\label {eqn:2.6} \]

2.2.3 贝叶斯法则

根据求和法则和求积法则,结合条件概率定义,可以得到 贝叶斯法则( Bayes rule ) 或者 贝叶斯定理 ( Bayes Theorem )

\[ p ( X=x|Y=y ) =\frac {p ( X=x,Y=y )}{p ( Y=y )}=\frac {p ( X=x ) p ( Y=y|X=x )}{\sum_{x^\prime} p ( X=x^\prime ) p ( Y=y|X=x^\prime )} ~~\tag {2.7}\label {eqn:2.7} \]

Sir Harold Jeffreys 认为:「贝叶斯定理之于概率论等价于勾股定理之于几何学」 。下面介绍两个贝叶斯理论的应用,不过在全书后面内容中,我们将会碰到更多例子。

2.2.3.1 案例:医疗诊断

考虑一个关于使用贝叶斯法则的案例:医疗诊断问题。假设你是一个 40 多岁的女性,你决定去做一个名叫 manmogram 的检测,以判断自己是否患有乳腺癌。如果检测结果显示为阳性,那么你患病的概率有多大?

这个问题的答案显然与检测方法的可靠性有关。假设你被告知检测方法的 敏感性 ( sensitivity ) 为 80%,也就是说,如果你患有癌症,那么测试结果显示为阳性的可能性为 0.8。换句话说:

\[ p ( x=1|y=1 ) =0.8 ~~\tag {2.8}\label {eqn:2.8} \]

其中 \(x=1\) 表示检测结果为阳性, \(y=1\) 表示你确实患有乳腺癌。许多人因此就认为他们患有癌症的可能性为 80%。但这是错误的!因为他们忽略了患有乳腺癌的先验概率,幸运的是,这个概率相当低:

\[ p ( y=1 ) =0.004 ~~\tag {2.9}\label {eqn:2.9} \]

忽略先验概率的情况被称为 基率谬论( base rate fallacy ) 。我们同样需要考虑测试结果可能是 假正例 ( false positive ) 或者 虚警 ( false alarm ) 的情况。不幸的是,类似这样的假正例发生的可能性还很高 ( 在现有的筛选技术下 ) :

\[ p ( x=1|y=0 ) =0.1 ~~\tag {2.10}\label {eqn:2.10} \]

将上述三项通过贝叶斯法则进行合并,我们可以计算出正确的答案:

\[\begin{split} \begin {align*} p ( y=1|x=1 ) & = \frac {p ( x=1|y=1 ) p ( y=1 )}{p ( x=1|y=1 ) p ( y=1 ) +p ( x=1|y=0 ) p ( y=0 )} ~~\tag {2.11}\label {eqn:2.11} \\ & = \frac {0.8 \times 0.004}{0.8 \times 0.004 + 0.1 \times 0.996} ~~\tag {2.12}\label {eqn:2.12} \end {align*} \end{split}\]

其中 \(p ( y=0 ) =1-p ( y=1 ) =0.996\) 。换句话说,如果你的检测结果为阳性,你也只有大概 3% 的可能性患有乳腺癌。 ( 该案例中的数据来自于相关文献,基于该项分析,美国政府决定不再推荐女性在 40 岁时进行每年的 mammogram 检测:因为假正例会导致不必要的担心和压力,并最终导致不必要的,昂贵的,并且可能存在潜在伤害的后续检测。5.7 节将介绍当我们在不确定的情况下权衡利弊的最优方法 ) 。

2.2.3.2 案例:生成式分类器

我们可以将医疗诊断的例子一般化,从而实现对任意形式的特征向量 \(\mathbf {x}\) 进行分类:

\[ p ( y=c|\mathbf {x},\mathbf {\theta} ) =\frac {p ( y=c \mid \mathbf {\theta}) p ( \mathbf {x}|y=c,\mathbf {\theta})}{\sum_{c^\prime} p ( y=c^\prime|\mathbf {\theta} ) p ( \mathbf {x}|y=c^\prime ,\mathbf {\theta})} ~~\tag {2.13}\label {eqn:2.13} \]

上式被称为 生成式分类器( generative classifier ) ,因为它通过使用类条件概率密度(似然) \(p ( \mathbf {x}|y=c )\) 和类先验分布 \(p ( y=c )\) 来指定如何生成样本。我们会在第 3 章和第 4 章详细讨论这一类模型。与该模型不同的是,判别式分类器直接对类后验概率分布 \(p ( y=c|\mathbf {x} )\) 进行训练。我们会在 8.6 节讨论两种方法的优缺点。

2.2.4 独立性和条件独立性

如果两个随机变量 \(X\)\(Y\) 的联合概率分布可以表示为边缘分布的乘积,则称 \(X\)\(Y\)无条件独立( unconditionally independent ) 或者 边缘独立 ( marginally independent ) 的,

\[ X \perp Y {\Leftrightarrow} p ( X,Y ) =p ( X ) p ( Y ) ~~\tag {2.14}\label {eqn:2.14} \]

一般情况下,如果一系列变量的联合概率分布等于各边缘分布的乘积,我们称这些变量相互独立。不幸的是,无条件独立很少会发生,因为大部分变量会相互影响,不过这种影响有时可通过其他变量间接产生。因此称已知变量 \(Z\) 的情况下, \(X\)\(Y\) 之间条件独立 ( conditionally independent,CI ) 的充要条件是:条件联合分布等于条件边缘分布的乘积。如下所示:

\[ X \perp Y |Z {\Leftrightarrow} p ( X,Y|Z ) =p ( X|Z ) p ( Y|Z ) ~~\tag {2.15}\label {eqn:2.15} \]

第 10 章会讨论概率图模型,届时可以发现这种关系可以表示为一种 \(X-Z-Y\) 形式的图结构,该结构反映出 \(X\)\(Y\) 之间的所有关联性都需要通过 \(Z\) 实现。举例来说,如果已经知道今天是否会下雨 ( \(Z\) ) ,那么明天将会下雨的概率 ( 事件 \(X\) ) 与今天的地是否潮湿 ( 事件 \(Y\) ) 之间是独立的。直觉上,因为 \(Z\) 同时导致了 \(X\)\(Y\) ,所以如果已经知道了 \(Z\) ,那么为了解 \(X\) ,并不需要关于 \(Y\) 的信息,反之亦然。

条件独立的另一个特性为:

定理 2.2.1 \(X \perp Y | Z\) 当且仅当存在函数 \(g\)\(h\) 满足:

\[ p ( x,y|z ) = g ( x,z ) h ( y,z ) ~~\tag {2.16}\label {eqn:2.16} \]

对于所有的 \(x,y,z\) 成立,且 \(p ( z )>0\) .

条件独立性假设允许通过一些局部信息去构建大规模概率模型。本书将介绍大量相关案例,特别是 3.5 节,将会讨论朴素贝叶斯分类器,在 17.2 节会讨论马尔可夫模型,在第 10 章讨论概率图模型;所有这些模型都充分应用了条件独立性的性质。

图 2.2 计算 \(p(x,y)=p(X)p(Y)\) ,其中 \(X⊥Y\) 。此处 \(X\)\(Y4 是离散型随机变量;\)X\( 有6种可能的状态,\)Y\( 有5种可能状态。正常情况下,两个变量上的联合分布需要 \)(6×5)−1=29\( 个参数来定义( 概率总和等于 1 的约束条件所致)。但通过无条件独立性假设,只需要 \)(6−1)+(5−1)=9\( 个参数来定义 \)p(x,y)$ 。 这种通过假设来简化模型的方法,在本书中会大量出现。

2.2.5 连续型随机变量

截至目前,我们只是讨论了关于离散型变量的情况。本节将介绍有关连续型变量的相关内容。

假设 \(X\) 是某个未知的连续型随机变量。变量 \(X\) 满足 \(a \le X \le b\) 的概率可以通过如下的方式进行计算。 接下来,定义三个事件 \(A= ( X \le a ) , B= ( X \le b ) , W= ( a \lt X \lt b )\)。 显然有 \(B = A \vee W\) , 因为事件 \(A\) 和事件 \(W\) 互斥,根据概率论的求和法则,有:

\[ p ( B ) = p ( A ) + p ( W ) ~~\tag {2.17}\label {eqn:2.17} \]

所以

\[ p ( W ) = p ( B ) - p ( A ) ~~\tag {2.18}\label {eqn:2.18} \]

接下来定义函数 \(F ( q ) \triangleq =p ( X\leq q)\),该函数被称为 累积分布函数 ( cumulative distribution function, cdf ) ,属于单调递增函数。 图 2.3 ( a ) 给出了示意图。基于该定义有:

\[ p ( a \lt X \le b ) =F ( b ) -F ( a ) ~~\tag {2.19}\label {eqn:2.19} \]

定义 \(f ( x ) =\frac {d}{dx} F ( x )\)( 假设导数存在 ) ,该式被称为 概率密度函数 ( probability density function, pdf ) 。 图 2.3 ( b ) 给出了示意图。 在已知概率密度函数的情况下,我们可以计算一个连续型变量属于某个有限区间的概率:

\[ p ( a \lt X \le b ) =\int_a^b f ( x ) dx ~~\tag {2.20}\label {eqn:2.20} \]

如果该区间足够小,可以有:

\[ p ( x \le X \le x+dx ) \approx p ( x ) dx ~~\tag {2.21}\label {eqn:2.21} \]

其中 \(p ( x ) \ge 0\) , 但对于任意给定的 \(X\) , \(p ( x ) \gt 1\) 是存在可能的,只要积分等于 1 即可。举例来说,考虑一个 均匀分布( uniform distribution ) \(Unif ( a,b )\) :

\[ Unif ( x|a,b ) =\frac {1}{b-a}\mathbb {I} ( a \le x\le b ) ~~\tag {2.22}\label {eqn:2.22} \]

如果令 \(a=0,b=\frac {1}{2}\) ,有 \(\forall x\in [0,\frac {1}{2}], p ( x ) =2\)

图 2.3 (a) 标准正态分布 \(N(0,1)\) 的 cdf 曲线图。(b) 相应的 pdf。每一个阴影区域都包含 \(α/2\) 的概率质量。因此,非阴影区域包含 \(1−α\) 的概率质量。如果分布为高斯 \(N(0,1)\),则最左边的截止点为 \(\Phi^{−1}(α/2)\) ,其中 \(\Phi\) 为高斯分布的 cdf。根据对称性,最右边的截止点是 \(\Phi^{−1}(1−α/2)=−\Phi^{−1}(α/2)\)。当 \(α=0.05\) 时,中心区间为 \(95%\) ,左界为 \(-1.96\),右界为 \(1.96\)

2.2.6 分位数

由于累积分布函数 \(F\) 是一个单调递增函数,其逆函数记作 \(F^-\)。如果 \(F\)\(X\) 的累积分布函数,那么 \(F^{-1}(\alpha)\) 就是满足概率 \(P (X\le x_\alpha )=\alpha \) 的值;这也叫做 \(F\)\(\alpha\) 分位数(quantile)。\(F^{-1}(0.5)\) 就是整个分布的中位数(median),左右两侧的概率各一半。而 \(F^{-1}(0.25)\)\(F^{-1}(0.75)\) 则是另外两个分位数。

利用累积分布函数(cdf)的逆函数还可以计算尾部概率。例如,如果有高斯分布 \(N (0,1)\)\(\Phi\) 是这个高斯分布的累积分布函数(cdf),这样在 \(\Phi^{-1}(\alpha/2)\) 左边的点就包含了 \(\alpha/2\) 的概率质量,如图 2.3(b)所示。与之对称,在 \(\Phi^{-1}(1-\alpha/2)\) 右边的点也包含了 \(\alpha /2\) 的概率质量。所以,在区间 \((\Phi^{-1}(\alpha/2),\Phi^{-1}(1-\alpha/2))\) 就包含了 \(1-\alpha\) 的概率质量。如果设置 \(\alpha =0.05\),那么中间 95% 的区间被以下范围覆盖:

\[ (\phi^{-1}(0.025),\phi^{-1}(0.975))=(-1.96,1.96) ~~\tag {2.23}\label {eqn:2.23} \]

如果分布为 \(N (\mu,\sigma^2)\),那么其 95% 区间就位于 \((\mu-1.96\sigma,\mu+1.96\sigma)\)。有时候简写成 \(\mu \pm 2\sigma\)

2.2.7 均值和方差

对正态分布,最常用的性质就是均值(mean),或期望值(expected value),记作 \(\mu\)。对于离散型随机变量,可定义成 \(\mathrm {E}[X] \overset\triangle {=} \sum_{x\in X} x p (x)\);对于连续型随机变量,可以定义为 \(\mathrm {E}[X] \overset\triangle {=} \int_{X} xp (x) dx\)。如果该积分是无穷的,则均值不能定义,更多例子后文会有。

方差(variance)表征的是分布的 “散开程度(spread)”,记作 \(\sigma^2\)。定义如下:

\[\begin{split} \begin {align*} var [X] & \triangleq E [(X-\mu)^2]=\int (x-\mu)^2p (x) dx ~~\tag {2.24}\label {eqn:2.24}\\ & = \int x^2p (x) dx +\mu^2 \int p (x) dx-2\mu\int xp (x) dx=E [X^2]-\mu^2 ~~\tag {2.25}\label {eqn:2.25}\\ \end {align*} \end{split}\]

从上面的式子就可以推导出:

\[ E [X^2]= \mu^2+\sigma^2 ~~\tag {2.26}\label {eqn:2.26} \]

然后就可以定义标准差(standard deviation)了:

\[ std [X]\overset\triangle {=} \sqrt {var [X]} ~~\tag {2.27}\label {eqn:2.27} \]

标准差与 \(X\) 单位一样。

2.3 离散型随机变量的常见分布

本节介绍一些定义在离散状态空间的常用概率分布,都是有限或者无限可列的。

2.3.1 二项分布和伯努利分布

如果抛 \(n\) 次硬币,设 \(X \in \{0,...,n\}\) 是人头朝上的次数。若头朝上的概率是 \(\theta\),就可以说 \(X\) 服从二项分布(binomial distribution),记作 \(X \sim Bin (n,\theta)\)。其 pmf(概率质量函数)可以写作:

\[ Bin (k|n,\theta)\overset\triangle {=} \binom {n}{k} \theta ^k (1- \theta)^{n-k} ~~\tag {2.28}\label {eqn:2.28} \]

上式中的

\[ \binom {n}{k} \overset\triangle {=} \frac {n!}{(n-k)!k!} ~~\tag {2.29}\label {eqn:2.29} \]

是组合数,相当于国内早期教材里面的 \(C_n^k\),从 \(n\) 中取 \(k\) 个样的组合数,也是二项系数(binomial coefficient)。如图 2.4 所示就是一些二项分布。

图 2.4 \(n=10\)\(θ∈\{0.25,0.9\}\) 的二项分布图。

该分布的均值和方差如下所示:

\[ mean=\theta, var =n\theta (1-\theta) ~~\tag {2.30}\label {eqn:2.30} \]

换个思路,如果只抛硬币一次,那么 \(X\in \{0,1 \}\) 就是一个二值化的随机变量,人头朝上的概率就是 \(\theta\)。这时候就说 \(X\) 服从伯努利分布(Bernoulli distribution),记作 \(X \sim Ber (\theta)\),其中概率质量函数 pmf 定义为(\(\mathbb{I}\) 为指示函数):

\[ Ber (x|\theta)=\theta^{\mathbb{I} (x=1)}(1-\theta)^{\mathbb{I} (x=0)} ~~\tag {2.31}\label {eqn:2.31} \]

也可以写成:

\[\begin{split} \begin {align*} Ber (x|\theta)=\begin {cases} \theta &\text { if x =1} ~~\tag {2.32}\label {eqn:2.32}\\ 1-\theta &\text { if x =0} \end {cases} \end {align*} \end{split}\]

很明显,伯努利分布只是二项分布中 \(n=1\) 的特例。

2.3.2 多项分布和多重伯努利分布

(1) 多项分布

二项分布可以用于抛硬币这种情况的建模。要对有 \(K\) 个可能结果的事件进行建模,就要用到多项分布(multinomial distribution)。这个定义如下:设 \(\mathbf{x} =(x_1,...,x_K)\) 是一个随机向量,其中的 \(x_j\) 是第 \(j\) 面出现的次数。这样 \(\mathbf{x}\) 的概率质量函数 pmf 就如下所示:

\[ Mu ( \mathbf{x} |n,\mathbf{\theta}) \triangleq \binom {n}{x_1,..,x_K}\prod^K_{j=1}\theta^{x_j}_j ~~~~\tag {2.33}\label {eqn:2.33} \]

其中的 \(\theta_j\) 是第 \(j\) 面出现的概率,组合的计算如下所示:

\[ \binom {n}{x_1,...,x_K} \overset\triangle {=} \frac {n!}{x_1!x_2!...x_K!} ~~\tag {2.34}\label {eqn:2.34} \]

这样得到的也就是多项式系数(multinomial coefficient),将一个规模为 \(n=\sum^K_{k=1}\) 的集合划分成规模从 \(x_1\)\(x_K\) 个子集的方案数。

(2) 多重伯努利分布(或类别分布)

接下来设 \(n=1\)。这好比将一个 \(K\) 面的骰子投掷一次,所以 \(\mathbf{x}\) 就是由 0 和 1 组成的向量。其中只有一个元素会是 1. 具体来说就是如果 \(k\) 面朝上,就说第 \(k\) 位为 1。这样就可以把 \(\mathbf{x}\) 看做一个用标量分类的有 \(K\) 个状态的离散型随机变量,\(\mathbf{x}\) 就是它的虚拟编码(dummy encoding),即:\(\mathbf{x} =[\mathbb{I} (x=1),...,\mathbb{I}(x=K)]\)\(\mathbb{I}\) 为指示函数。例如,如果 \(K=3\),那么状态 1、2、3 对应的虚拟编码分别是 \((1,0,0),(0,1,0),(0,0,1)\)。该编码也称作独热编码(one-hot encoding),因为只有一个位置是 1. 其对应的概率质量函数 pmf 就是:

\[ Mu (\mathbf{x}|1,\mathbf{\theta})=\prod^K_{j=1 }\theta_j ^{ \mathbb{I} (x_j=1)} ~~\tag {2.35}\label {eqn:2.35} \]

可以参考图 2.1 的(b-c)作为一个例子。这是类别分布(categorical distribution)或者离散分布(discrete distribution)的典型案例,Gustavo Lacerda 建议大家称之为多重伯努利分布(multinoulli distribution),这样与二项分布 / 伯努利分布的区别关系相仿。本书就采用了这个术语,使用下列记号表示这种情况:

\[ Cat (x|\mathbf{\theta})\overset\triangle {=} Mu (\mathbf{x}|1,\mathbf{\theta}) ~~\tag {2.36}\label {eqn:2.36} \]

换句话说,如果 \(x\sim Cat (\mathbb{\theta})\),则 \(p (x=j|\mathbb{\theta})=\theta_j\)。参考表 2.1。

表2.1 多项分布及相关分布汇总。

(3) 应用:DNA 序列模体

图2.5 (a)一些经过对齐的DNA序列。(b)相应的序列标识。

生物序列分析(biosequence analysis)是一个典型应用案例,设有一系列对齐的 DNA 序列,如图 2.5(a)所示,其中有 10 行(序列),15 列(沿着基因组的位置)。图中几个位置是进化的保留位,是基因编码区域的位置,所以对应的列都是 “纯的”,例如第 7 列就都是 g。

如图 2.5(b)所示的可视化方法是序列标识图(sequence logo)。具体方法是把字母 A、C、G 和 T 投到对应位置上,字体大小与其经验概率(empirical probability)成正比,最大概率的字母放在最顶部。

位置 \(t,\hat \theta_t\) 处的经验概率分布通过计数向量的归一化获得,可以参考本书公式 3.48:

\[ \mathbf{N}_t=(\sum^N_{i=1} \mathbb{I}(X_{it}=1),\sum^N_{i=1}\mathbb{I} (X_{it}=2),\sum^N_{i=1}\mathbb{I} (X_{it}=3),\sum^N_{i=1}\mathbb{I} (X_{it}=4)) ~~\tag {2.37}\label {eqn:2.37} \]
\[ \hat\theta_t =\mathbf{N}_t/N ~~\tag {2.38}\label {eqn:2.38} \]

该分布被称作一个模体(motif)。可以计算每个位置上最可能出现的字母,得到共有序列(consensus sequence)

2.3.3 泊松分布

如果一个离散型随机变量 \(X\in \{0,1,2,...\}\) 服从泊松分布,即 \(X\sim Poi (\lambda)\),其参数 \(\lambda >0\),其概率质量函数 pmf 为:

\[ \text{Poi} (x|\lambda )=e^{-\lambda}\frac {\lambda ^x}{x!} ~~\tag {2.39}\label {eqn:2.39} \]

第一项是归一化常数(normalization constant),使用来保证概率密度函数的总和是 1.

泊松分布经常用于对罕见事件的统计建模,比如放射性衰变和交通事故等等,可被视为 \(\theta\) 很小 但 \(n\) 很大的二项分布。图 2.6 是一些泊松分布的示意图。

图 2.6. \(\lambda \in \{ 1,10 \}\) 对应的泊松分布图。为清楚起见,将 \(x\) 轴截断为 \(25\)

2.3.4 经验分布

给定数据集 \(D =\{x_1,...,x_N \}\),可以定义经验分布(empirical distribution),也称经验测度(empirical measure),形式如下:

\[ p_{emp}(A)\overset\triangle {=}\frac 1 N \sum^N_{i=1}\delta _{x_i}(A) ~~\tag {2.40}\label {eqn:2.40} \]

其中 \(\delta_x (A)\)狄拉克测度(Dirac measure),定义为:

\[\begin{split} \begin {align*} \delta_x (A)= \begin {cases} 0 \qquad \text {if } x \notin A \\ 1 \qquad \text {if } x \in A ~~\tag {2.41}\label {eqn:2.41} \end {cases} \end {align*} \end{split}\]

一般来说可以给每个样本关联一个权重(weight)

\[ p (x)=\sum^N_{i=1} w_i\delta_{x_i}(x) ~~\tag {2.42}\label {eqn:2.42} \]

其中要满足 \(0\le w_i \le 1\) 以及 \(\sum^N_{i=1} w_i=1\)。可以将其想象成在数据点 \(x_i\) 处有峰(spike)的一个直方图,而 \(w_i\) 决定了 \(i\) 峰的高低。此分布中所有不在数据集中的点概率均为 0。

2.4 连续性随机变量的常见分布

接下来介绍一些常用的单变量一维连续概率分布。

2.4.1 高斯(正态)分布

不管统计学还是机器学习,最广泛使用的都是高斯分布(Gaussian distribution),也称正态分布(Normal distribution)。其概率密度函数 pdf 为:

\[ \mathcal{N} (x|\mu,\sigma^2) \overset\triangle {=} \frac {1}{\sqrt {2\pi \sigma^2}} e^ {-\frac {1}{2 \sigma^2}(x-\mu)^2} ~~\tag {2.43}\label {eqn:2.43} \]

上式中的 \(\mu=\mathbb{E} [X]\) 是均值(mean),\(\sigma^2=var [X]\) 是方差(variance)。\(\sqrt {2\pi \sigma^2}\) 是归一化常数,用于确保整个密度函数的积分为 1,具体可以参考练习 2.11。

Note

注:高斯分布的均值与众数(mode)相等。

可以用 \(X \sim \mathcal{N}(\mu,\sigma^2)\) 来表示 \(p (X=x)=\mathcal{N}(x|\mu,\sigma^2)\)

\(X \sim \mathcal{N} (0,1)\)标准正态分布(standard normal distribution)。图 2.3(b)是标准正态分布的概率密度函数图,也被称作钟形曲线。

有时用高斯分布的精度(precision) \(\lambda =1/\sigma^2\) 来表达钟形的分散程度,即方差的倒数。精度高,则方差低,图形的大部分对称分布在以均值为中心的狭窄区域。

Note

注意连续型随机变量的概率密度函数(pdf)与离散型随机变量的概率质量函数(pmf)之间的区别:

(1)pmf 定义域不连续且函数值不会大于 1,且所有可能取值的 pmf 总和为 1。

(2) pdf 定义域连续,函数值可能大于 1,但是在整个定义域上 pdf 的积分为 1。 比如在中心位置,\(x=\mu\), 有 \(\mathcal{N} (\mu \mid \mu,\sigma^2)=(\sigma\sqrt {2\pi})^{-1} e^0\), 当标准差 \(\sigma< 1/\sqrt {2\pi}\) 时, 就会出现 \(p (x)>1\) 的情形。

(3)pmf 求和得出事件的概率,而 pdf 求积分得出事件的概率。

高斯分布的累积分布函数 (cdf) 为:

\[ \Phi (x;\mu , \sigma^2)\overset\triangle {=} \int^x_{-\infty} \mathcal{N} (z|\mu,\sigma^2) dz ~~\tag {2.44}\label {eqn:2.44} \]

图 2.3 (a) 所示为 \(\mu=0,\sigma^2=1\) 时的 cdf 函数曲线。该积分没有闭合形式的表达式,但在多数软件包中都内置了。特别地,可以用误差函数 (error function, 缩写为 erf) 来计算 cdf:

\[ \Phi (x;\mu , \sigma^)\overset\triangle {=} \frac 1 2 [1+ \text{erf} (z/\sqrt2)] ~~\tag {2.45}\label {eqn:2.45} \]

其中 \(z=(x-\mu)/\sigma\), 误差函数为:

\[ \text{erf} (x)\overset\triangle {=} \frac {2}{\sqrt\pi}\int^x_0e^{-t^2} dt ~~\tag {2.46}\label {eqn:2.46} \]

高斯分布是统计学里使用最广的分布,有几个原因:

(1)这两个参数很好解释,分别对应着分布中的两个基础特征,均值和方差。

(2)中心极限定理 ( central limit theorem, 参考本书 2.6.3) 表明独立随机变量的和就近似为高斯分布,所以高斯分布很适合用来对残差或者噪音建模。

(3)高斯分布有最少假设 (least number of assumptions)最大熵 (maximum entropy), 适合对特定场景建立有特定均值和方差的约束和假设,如本书 9.2.6 所述,这使其成为很多情况下的默认选择。

(4)高斯分布的数学形式很简单,容易实现,效率也很高。高斯分布更广泛应用参考(Jaynes, 2003)的第七章。

2.4.2 狄拉克函数(退化的高斯分布)

如果让高斯分布的方差趋近于零,即 \(\sigma^2 \rightarrow 0\),则高斯分布就变成高度无穷大而峰值宽度无穷窄的形状了,当然中心还是在 \(\mu\) 位置:

\[ \lim_{\sigma^2\rightarrow 0} N (x| \mu,\sigma^2) =\delta (x-\mu) ~~\tag {2.47}\label {eqn:2.47} \]

此函数 \(\delta\) 被称为狄拉克函数 (Dirac delta function)。 其定义为:

\[\begin{split} \begin{align*} \delta (x)=\begin {cases}\infty \qquad \text{if} x=0\\ 0 \qquad \text{if} x\ne 0 ~~\tag {2.48}\label {eqn:2.48} \end {cases} \end{align*} \end{split}\]

狄拉克函数的积分满足:

\[ \int_{-\infty} ^\infty \delta (x) dx=1 ~~\tag {2.49}\label {eqn:2.49} \]

狄拉克函数的特点就是筛选性 (sifting property), 可以从一系列求和或积分中筛选出某个单一元素项:

\[ \int_{-\infty} ^\infty f (x) \delta (x-\mu ) dx=f (\mu ) ~~\tag {2.50}\label {eqn:2.50} \]

上式中只有当 \(x-\mu=0\) 的时候积分才非零。

2.4.3 学生 \(t\) 分布

高斯分布有个问题就是对异常值很敏感,因为从分布中心往外的对数概率衰减速度和距离成平方关系。有一个更健壮的分布,就是所谓的学生 \(t\) 分布 (student distribution), 其概率密度函数如下所示:

\[ \mathcal{T} (x|\mu,\sigma^2,v)\propto \left [1+\frac 1v (\frac {x-\mu}{\sigma})^2 \right ]^ {-(\frac {v+1}{2})} ~~\tag {2.51}\label {eqn:2.51} \]

上式中的 \(\mu\) 是均值,\(\sigma^2>0\) 是缩放参数 (scale parameter),\(\nu > 0\) 称为自由度 ( degrees of freedom)。 图 2.7 为该函数的曲线。为后文方便,这里介绍几个属性:

\[ mean=\mu, \qquad mode=\mu, \qquad var=\frac {\nu \sigma^2}{\nu -2} ~~\tag {2.52}\label {eqn:2.52} \]

该模型中,当自由度大于 \(\nu > 2\) 时方差才有意义,自由度 \(\nu > 1\) 均值才有意义。

图 2.7 (a) \(\mathcal{N}(0,1)\)\(\mathcal{T}(0,1,1)\)\(\text{Lap}(0,1/\sqrt{2})\) 的 pdf。高斯和拉普拉斯的均值均为 \(0\) ,方差均为 \(1\) 。当 \(\nu =1\) 时,学生分布的均值和方差都未定义。(b) 上述 pdf 的对数。注意,学生分布对任何参数值都不是对数凹的,而拉普拉斯分布总是对数凹的(和对数凸的)。然而,两者都是单峰的。

图2.8 异常值对高斯分布、学生分布和拉普拉斯分布的影响。(a)无异常值的高斯曲线和学生分布曲线重叠。(b)有异常值时,与学生分布和拉普拉斯分布相比,高斯分布更易受异常值影响。根据(Bishop,2006a)的图 2 .16。

\(\mathcal{T}\) 分布的稳定性如图 2.8 所示,左侧是没有异常值的高斯分布和 \(\mathcal{T}\) 分布,右侧是加入了异常值的。很明显异常值对于高斯分布干扰很大,而 \(\mathcal{T}\) 分布则几乎看不出来有影响。因为 \(\mathcal{T}\) 分布比高斯分布更重尾, 至少对于小自由度 \(\nu\) 的时候是这样,如图 2.7 所示。

如果自由度 \(\nu =1\) , 则 \(\mathcal{T}\) 分布就成了柯西分布 (Cauchy distribution) 或者洛伦兹分布 (Lorentz distribution)。 要注意这时候重尾会导致定义均值 (mean) 的积分不收敛。

要确保有限范围的方差, 就需要自由度 \(\nu > 2\)。 一般常用自由度是 \(\nu =4\), 在一系列问题中性能表现不错 (Lange等, 1989)。 如果自由度远超过 5, 即 \(\nu >> 5\)\(\mathcal{T}\) 分布很快近似到高斯分布,也就失去了对异常值的稳健性。

2.4.4 拉普拉斯分布

另外一个常用的重尾分布就是拉普拉斯分布(Laplace distribution),也被称为双面指数分布 (double sided exponential distribution), 其概率密度函数如下所示:

\[ \text{Lap} (x|\mu,b)\overset\triangle {=}\frac1 {2b}\exp (-\frac {|x-\mu|}{b}) ~~\tag {2.53}\label {eqn:2.53} \]

上式中的 \(\mu\) 是位置参数, \(b > 0\) 是缩放参数, 如图 2.7 所示就是其曲线。该分布的各属性如下所示:

\[ mean=\mu,\qquad mode=\mu, \qquad var=2b^2 ~~\tag {2.54}\label {eqn:2.54} \]

该分布的健壮性如图 2.8 中所示,从图中也可以发现拉普拉斯分布比高斯分布在 0 点有更高概率。此特性在 13.3 节要用到,很适合在模型中增强稀疏性。

2.4.5 伽马分布族

此分布很灵活,适用于正实数值的随机变量 \(x>0\) 。 该分布用两个参数定义,分别是形状 (shape) \(a > 0\)速率 (rate) \(b > 0\)

\[ \text{Ga} (T|shape=a, rate=b)\overset\triangle {=} \frac {b^a}{\Gamma (a) } T^{a-1} e^{-Tb} ~~\tag {2.55}\label {eqn:2.55} \]

上式中的 \(\Gamma (a)\) 指伽马函数:

\[ \Gamma (x) \overset\triangle {=} \int_0^{\infty} u^{x-1} e^{-u} du ~~\tag {2.56}\label {eqn:2.56} \]

图 2.9 是一些伽马分布的示例。该分布的各属性如下:

\[ mean=\frac {a}{b},\qquad mode=\frac {a-1}{b},\qquad var=\frac {a}{b^2} ~~\tag {2.57}\label {eqn:2.57} \]

图 2.9 (a)一些 \(Ga(a,b=1)\) 分布。如果为 \(a≤1\) ,则众数为 \(0\) ,否则为 \(a>0\)。我们增加速率 \(b\),减小水平尺度,就把所有东西都往左往上挤了。(b)一些降雨量数据的经验 pdf,叠加了拟合的Gamma分布。

有一些分布实际上是伽马分布的特例,比如下面这几个:

  • 指数分布 (Exponential distribution) 是形状参数 \(a\) 为 1 的伽马分布,定义为 \(\text{Expon}(x|\lambda) \overset\triangle {=} \text{Ga}(x|1,\lambda)\), 其中的 \(\lambda\) 是速率参数 (rate)。 此分布用于建模泊松过程 (Poisson process) 中各事件之间的时间间隔。例如,一个过程可能由一系列事件按照某个固定的平均速率 \(\lambda\) 连续独立发生。

  • 厄兰分布 (Erlang Distribution) 是形状参数 \(a\) 为整数的伽马分布。经常设置 \(a=2\) , 产生一个单参数的厄兰分布,即 \( \text{Erlang} (x|\lambda) = \text{Ga} (x|2, \lambda)\)\(\lambda\) 为速率参数(rate)。

  • 卡方分布 (Chi-squared distribution) 的定义为 \(\chi^2 (x|\nu) \overset\triangle {=} \text{Ga} (x|\frac{\nu}{2},\frac{1}{2})\). 此分布是高斯分布随机变量平方和的分布。更确切地说,如果有 \(Z_i \sim \mathcal{N} (0, 1)\), 而且 \(S=\sum_{i=1}^\nu Z_i^2\), 则 \(S\) 服从卡方分布: \(S \sim \chi_\nu^2\).

另一个有用的结果是:如果一个随机变量服从伽马分布:\(X \sim \text{Ga} (a,b)\) 那么这个随机变量的倒数就服从一个逆伽马分布, 即 \(\frac{1}{X} \sim \text{IG} (a,b)\), 这个在练习 2.10 里面有详细描述。逆伽马分布 (inverse gamma) 定义如下:

\[ \text{IG} (x|shape =a,scale =b)\overset\triangle {=} \frac {b^a}{\Gamma (a)} x^{-(a+1)} e^{-\frac b x} ~~\tag {2.58}\label {eqn:2.58} \]

逆伽马分布的各属性如下:

\[ mean= \frac {b}{a-1}, \qquad mode=\frac {b}{a+1}, \qquad var=\frac {b^2}{(a-1)^2 (a-2)} ~~\tag {2.59}\label {eqn:2.59} \]

均值只在 \(a>1\) 的情况下才存在(收敛),而方差仅在 \(a>2\) 的时候存在(收敛)。

2.4.6 贝塔分布

贝塔分布的定义域为区间 [0,1], 定义如下:

\[ \text{Beta} (x|a,b)=\frac {1}{B (a,b)} x^{a-1 }(1-x)^{b-1} ~~\tag {2.60}\label {eqn:2.60} \]

上式中 \(B (a,b)\) 是贝塔函数,定义如下:

\[ B (a,b)\overset\triangle {=}\frac {\Gamma (a)\Gamma (b)}{\Gamma (a+b)} ~~\tag {2.61}\label {eqn:2.61} \]

该分布见图 2.10 。 需要 \(a\)\(b\) 都大于 0 来确保整个分布可积分,因为需要保证 \(B (a,b)\) 存在。如果 \(a=b=1\) ,得到的是均匀分布 (uniform distirbution), 如图 2.10 中红色虚线所示。如果 \(a\)\(b\) 都小于 1 , 那么得到的就是一个双峰分布 (bimodal distribution), 两个峰值在 0 和 1 位置上,如图 2.10 中的蓝色实线所示。如果 \(a\)\(b\) 都大于 1 了,得到的就是单峰分布 (unimodal distribution) ,如图 2.10 中的另外两条虚线所示。

图 2.10 一些贝塔分布图。

该分布的各属性如下,在练习 2.16 里会用到:

\[ mean= \frac {a}{a+b},\qquad mode=\frac {a-1}{a+b-2},\qquad var=\frac {ab}{(a+b)^2 (a+b+1)} ~~\tag {2.62}\label {eqn:2.62} \]

2.4.7 帕累托分布

帕累托分布 (Pareto distribution)是用来对具有长尾(或重尾)特点的变量进行建模的,特别是符合“80/20”法则特征的随机变量。例如,英语中最常出现的词汇是冠词 the,其出现概率为第二位 of 的 2 倍多,而 of 的出现概率是第四位单词的 2 倍,诸如此类。如果将每个词汇词频和排名作图,得到的是一个幂律 (power law), 也称为齐夫定律 (Zipf's law)。人类财富的分配也有类似特点,尤其是在资本主义国家。

此分布的概率密度函数 (pdf) 如下:

\[ \text{Pareto} (x|k,m)=km^kx^{-(k+1)}\mathbb{I}(x\geq m) ~~\tag {2.63}\label {eqn:2.63} \]

通过定义可知,\(x\) 必须比某一个常数 \(m\) 大,但又不用大特别多,而其中 \(k\) 则控制这个的度,避免 \(x\) 太大。随着 \(k \rightarrow \infty\), 该分布接近于狄拉克分布 \(\delta (x-m)\) 。图 2.11 (a) 就是一些帕累托分布图,如果用对数坐标来制图,会形成一条直线,如图 2.11 (b) 所示(也叫幂律)。这个直线的方程形式为 \(\log p (x)=a\log x+c\), 其中的 \(a\)\(c\) 是某个常数。

帕累托分布的属性如下:

\[ mean=\frac {km} {k-1} \text {if } k>1, \qquad mode=m, \qquad var =\frac {m^2k}{(k-1)^2 (k-2)} \text{if } k > 2 ~~\tag {2.64}\label {eqn:2.64} \]

图 2.11 (a) \(m=1\) 时的帕累托分布 \(\text{Pareto}(x|m,k)\) 。 (b) 对数-对数尺度上的概率密度函数(pdf)。

2.5 联合概率分布

前面都是单变量的概率分布,接下来看更复杂的多变量概率分布,即联合概率分布 ( joint probability distributions)。 其中涉及到多个相关联的随机变量,也是本书核心内容。

对于规模 \(D>1\) 个的随机变量集合,其联合概率分布表示为 \(p (x_1,...,x_D)\), 被用于对这些随机变量间的(随机)关系进行建模。

如果所有变量都是离散型变量,则可以把联合概率分布表示成一个多维数组,每个维度对应一个变量。如果设每个变量的状态数目总共是 \(K\) , 则按照这种表示法,需要定义 \(O(K^D)\) 个参数。当然,通常可以通过条件独立假设,使用更少的参数来定义高维联合分布(第 10 章会做进一步解释)。

对于连续型变量,为了减小复杂度,比较常用的做法是将概率密度函数(pdf )限制为某些特定的形式(如:假设各变量相互独立且呈高斯分布)。

2.5.1 协方差和相关系数

随机变量 \(X\)\(Y\) 之间的协方差 (covariance) 用来衡量两者之间的(线性)相关程度,定义如下:

\[ cov [X,Y]\overset\triangle {=} \mathbb{E} [(X-\mathbb{E}[X])(Y-\mathbb{E}[Y])] =\mathbb{E}[XY]-\mathbb{E}[X] \mathbb{E}[Y] ~~\tag {2.65}\label {eqn:2.65} \]

如果 \(\mathbf{x}\) 是一个 \(d\) 维度的随机向量(每个维度对应一个随机变量),则其协方差矩阵 (covariance matrix) 定义如下,是一个对称正定矩阵 (symmetric positive definite matrix):

\[\begin{split} \begin {align*} cov [x] &\triangleq E [(x-E [x])(x-E [x])^T] ~~\tag {2.66}\label {eqn:2.66}\\ &= \begin {pmatrix} var [X_1] & cov [X_1,X_2] &...& cov [X_1,X_d] \\ cov [X_2,X_1] & var [X_2] &...&cov [X_2,X_d] \\ ...&...&...&...\\ cov [X_d,X_1] & cov [X_d,X_2] &...&var [X_d] \\ \end {pmatrix} ~~\tag {2.67}\label {eqn:2.67}\\ \end {align*} \end{split}\]

协方差可以从 0 到 \(\infty \) 之间取值。但为使用方便,会将其标准化为有上限。

其中比较常用的是皮尔逊相关系数 (correlation coefficient),变量 \(X\)\(Y\) 之间的皮尔逊相关系数为:

\[ \text{corr} [X,Y]\overset\triangle {=} \frac {cov [X,Y]}{\sqrt {var [X] var [Y]}} ~~\tag {2.68}\label {eqn:2.68} \]

相应的,所有维度一起构成相关矩阵 (correlation matrix) :

\[\begin{align*} R= \begin {pmatrix} \text{corr} [X_1,X_1] & \text{corr} [X_1,X_2] &...& \text{corr} [X_1,X_d] \\ ...&...&...&...\\ \text{corr} [X_d,X_1] & \text{corr} [X_d,X_2] &...&\text{corr} [X_d] \\ \end {pmatrix} ~~\tag {2.69}\label {eqn:2.69} \end{align*}\]

从练习 4.3 可知相关系数位于 [-1,1] 区间内。在相关矩阵中,每一个对角线元素值都是 1, 其他的值都在 [-1,1] 区间内。

另外,当且仅当参数 \(a\)\(b\) 满足 \(Y = aX + b\) 时,才有 \(corr [X, Y] = 1\), 也就是说 \(X\)\(Y\) 之间存在线性关系,参考练习 4.3.

直觉可能让人觉得相关系数和线性回归的斜率有关,比如说像 \(Y = aX + b\) 这个表达式当中的系数 \(a\) 一样。然而并非如此,如公式 7.99 所示,回归系数公式中的 \(a = cov [X, Y] /var[X]\),是回归线斜率的量化, 而相关系数则完全不同,可以看做是两个变量之间线性程度的衡量,参考图 2.12.

回想本书 2.2.4 节, 如果 \(X\)\(Y\) 相互独立,则有 \(p(X, Y) = p(X) p(Y)\), 二者的协方差 \(cov [X,Y] = 0\), 相关系数 \(\text{corr}[X,Y] = 0\), 很好理解,相互独立就是不相关了。但反过来不成立,不相关并不能意味着相互独立。例如设 \(X \sim U (-1,1), Y=X^2\) , 很明显,\(X\)\(Y\) 相关,甚至 \(Y\) 就是 \(X\) 唯一决定的,然而如练习 4.1 所示,这两个变量的相关系数等于零,即 \(corr [X,Y]=0\) 。 图 2.12 有更多直观的例子,都是两个变量 \(X\)\(Y\) 具有明显的相关性,但计算出来的相关系数却都是 0 。 实际上更通用的用来衡量两组随机变量之间是否独立的工具是互信息量 (mutual information),这部分在 2.8.3 节中有涉及。如果两个变量真正不相关,互信息量才会等于 0.

图 2.12 几组 \((x,y)\) 点,每组的相关系数为 \(x\)\(y\)。请注意,相关性反映了线性关系的噪声和方向(顶行),但不是该关系的斜率(中间),也不是非线性关系的许多方面(底部)。注:中间的数字斜率为 0 ,但在这种情况下,相关系数未定义,因为 \(Y\) 的方差为 0 。

此处查看原书图 2.13

2.5.2 多元高斯分布

多元高斯分布 (multivariate Gaussian) 或者 多元正态分布 (multivariate normal, 缩写为 MVN), 是连续型随机变量的联合概率分布中应用最广的分布。在第 4 章会对其进行详细说明,这里给出简单定义。

\(D\) 维上的多元正态分布的定义如下所示:

\[ \mathcal{N} (\mathbf{x}|\mu,\Sigma)\overset\triangle {=} \frac {1}{(2\pi )^{\frac D2} |\Sigma|^{\frac12}}\exp [-\frac12 (\mathbf{x}-\mu)^T\Sigma^{-1}(\mathbf{x}-\mu) ] ~~\tag {2.70}\label {eqn:2.70} \]

上式中 \(\mu = \mathbb{E}[x] \in R^D\) 是均值向量,而 \(\Sigma= cov [\mathbf{x}]\) 是一个 \( D\times D\) 的协方差矩阵。

有时候会用到一个叫做精度矩阵 (precision/concentration matrix)的东西, 其实就是协方差矩阵的逆矩阵,也就是 \(\Lambda =\Sigma^{-1 }\). 前面分母为 \((2\pi )^{\frac D2}|\Sigma|^{\frac12}\) 的项为归一化常数,以保证概率密度函数的积分等于 1。

图 2.13 展示了一些多元正态分布的密度图,其中有三个是不同协方差矩阵的下的二维投影,另外一个是立体的曲面图像。一个完整的协方差矩阵有 \(D (D + 1)/2\) 个参数,除以 2 是因为矩阵 \(\Sigma\) 是对称的。对角协方差矩阵的方向有 \(D\) 个参数,非对角线位置的元素的值都是 0。 球面 (spherical) 或者各向同性 (isotropic) 协方差矩阵 \(\Sigma = \delta^2 I_D\) 有一个自由参数。

图2.13 二维高斯的等值线图。(a) 完全协方差矩阵具有椭圆等值线。(b) 对角协方差矩阵是轴对齐的椭圆。(c) 球形协方差矩阵具有圆形形状。(d) c 中球形高斯的曲面图。

2.5.3 多元学生 \(t\) 分布

相比多元正态分布, 多元学生 \(t\) 分布更加健壮,其概率密度函数为:

\[\begin{split} \begin {align*} \mathcal{T} (\mathbf{x}|\mu,\Sigma,\nu)&=\frac {\Gamma (\nu/2+D/2)}{\Gamma (\nu/2)} \frac {|\Sigma|^{-1/2}}{\nu^{D/2}\pi^{D/2}}\times [1+\frac1\nu (\mathbf{x}-\mu )^T\Sigma^{-1}(\mathbf{x}-\mu)]^{-(\frac {\nu+D}{2})} ~~\tag {2.71}\label {eqn:2.71}\\ &=\frac {\Gamma (\nu/2+D/2)}{\Gamma (\nu/2)} |\pi V|^{-1/2}\times [1+(\mathbf{x}-\mu)^TV^{-1}(\mathbf{x}-\mu)]^{-(\frac {\nu+D}{2})} ~~\tag {2.72}\label {eqn:2.72}\\ \end {align*} \end{split}\]

其中的 \(\Sigma\) 叫做范围矩阵 (scale matrix), 而并不是真正的协方差矩阵,\(V=\nu\Sigma\). 这个分布比高斯分布有更重的尾部 (fatter tails)。 参数 \(\nu\) 越小,越重尾;而当 \(\nu \rightarrow \infty\) ,这个分布趋向于高斯分布。

该分布的属性如下:

\[ mean=\mu, \qquad mode=\mu, \qquad Cov=\frac {\nu}{\nu-2}\Sigma ~~\tag {2.73}\label {eqn:2.73} \]

2.5.4 狄利克雷分布

贝塔布扩展到多元就成了狄利克雷分布 (Dirichlet distribution), 支持概率单纯形 (probability simplex), 定义如下:

\[ S_K=\{x:0 \le x_k \le 1, \sum ^K_{k=1} x_k=1\} ~~\tag {2.74}\label {eqn:2.74} \]

其概率密度函数 pdf 如下所示:

\[ \text{Dir} (x|\alpha)\overset\triangle {=} \frac {1}{B (\alpha)} \prod^K_{k=1} x_k^{\alpha_k -1} \mathbb{I} (x\in S_K) ~~\tag {2.75}\label {eqn:2.75} \]

式中 \(B (\alpha_1,...,\alpha_K)\) 是贝塔函数在 \(K\) 个变量上的自然推广, 定义如下:

\[ B (\alpha)\overset\triangle {=} \frac {\prod^K_{k=1}\Gamma (\alpha_k)}{\Gamma (\alpha_0)} ~~\tag {2.76}\label {eqn:2.76} \]

其中的 \(\alpha_0\overset\triangle {=} \sum^K_{k=1}\alpha_k\)

图2.14 (a) 当 \(K = 3\) 时的狄利克雷分布定义了单形上的分布,可以用三角曲面表示。该曲面上的点满足 \(0 ≤ θ_k≤ 1\)\(\sum_{k=1}^3θ_k= 1\) 。(b) 当 \(α = (2,2,2)\) 时的狄利克雷密度图。 (c) \(α = (20,2,2)\) 。 (d) \(α = ( 0.1,0.1,0.1)\)

图 2.15 不同参数值的 5 维对称狄利克雷分布样本。 (a) \(α = ( 0.1,...,0.1)\),这导致分布非常稀疏,有许多 0 。(b) \(α = ( 1,...,1)\),这导致更均匀和密集的分布。

图 2.14 展示的是当 \(K=3\) 时的一些狄利克雷函数图,图 2.15 是一些概率向量样本。很明显其中 \(\alpha_0\overset\triangle {=} \sum^K_{k=1}\alpha_k\) 控制了分布强度,也就是峰值位置。例如 \(\text{Dir} (1, 1, 1)\) 是一个均匀分布,\(\text{Dir} (2, 2, 2)\) 是以 \((1/3, 1/3, 1/3)\) 为中心的宽分布 (broad distribution), 而 \(\text{Dir} (20, 20, 20)\) 是以 \((1/3, 1/3, 1/3)\) 为中心的窄分布 (narrow distribution)。如果对于所有的 \(k\) 都有 \(\alpha_k <1\), 那么峰值位于单纯形的角落。

狄利克雷分布的属性如下:

\[ \mathbb{E} [x_k]=\frac {\alpha_k}{\alpha_0},\qquad mode [x_k]=\frac {\alpha_k-1}{\alpha_0-K},\qquad var [x_k]=\frac {\alpha_k (\alpha_0-\alpha_k)}{\alpha_0^2 (\alpha_0+1)} ~~\tag {2.77}\label {eqn:2.77} \]

其中 \(\alpha_0 = \sum_k \alpha_k\). 通常我们使用形式为 \(\alpha_k=\alpha/K\) 的对称狄利克雷先验,则有均值为 \(1/K\), 方差为 \(var [x_k]=\frac {K-1}{K^2 (\alpha+1)}\)。 此时增大 \(\alpha\) 就能降低方差,提高了模型精度。

2.6 随机变量变换

如果有一个随机变量 \(\mathbf{x} \sim p ()\), 还有个函数 \(\mathbf{y}=f (\mathbf{x})\), 那么 \(\mathbf{y}\) 的分布是什么呢?这就是本节要讨论的内容。

2.6.1 线性变换

\(f()\) 是一个线性函数:

\[ \mathbf{y}=f (\mathbf{x})=A\mathbf{x}+b ~~\tag {2.78}\label {eqn:2.78} \]

这样就可以推导 \(\mathbf{y}\) 的均值和协方差。

均值如下:

\[ \mathbb{E} [\mathbf{y}]=\mathbb{E} [A\mathbf{x}+b]=A\mu+b ~~\tag {2.79}\label {eqn:2.79} \]

上式中的 \(\mu=\mathbb{E}[\mathbf{x}]\). 这就叫线性期望 (linearity of expectation)。 如果 \(f ()\) 是一个标量值的函数 \(f (\mathbf{x})=a^T\mathbf{x}+b \), 那么对应结果就是:

\[ \mathbb{E} [a^T\mathbf{x}+b]=a^T\mu+b ~~\tag {2.80}\label {eqn:2.80} \]

对于协方差,得到的就是(证明略):

\[ cov [\mathbf{y}]=cov [A\mathbf{x}+b]=A\Sigma A^T ~~\tag {2.81}\label {eqn:2.81} \]

其中的 \(\Sigma =cov [\mathbf{x}]\)。如果 \(f ()\) 的值为标量,则结果为:

\[ var [y]=var [a^T\mathbf{x}+b]=a\Sigma a^T ~~\tag {2.82}\label {eqn:2.82} \]

这些结果后面的章节都会多次用到。不过这里要注意,只有 \(\mathbf{x}\) 服从高斯分布时,才能单凭借着均值和协方差来定义 \(\mathbf{y}\) 的分布。通常必须使用一些技巧来对 \(\mathbf{y}\) 的完整分布进行推导,而不能只靠两个属性就确定。

2.6.2 通用变换

如果 \(X\) 是一个离散型随机变量,而 \(Y\)\(X\) 的变换 \(f(x) = y\),可以简单通过将所有 \(X\) 的概率质量相加得出 \(y\) 的概率质量函数(pmf):

\[ p_y (y)=\sum_{x:f (x)=y} p_x (x) ~~\tag {2.83}\label {eqn:2.83} \]

例如,若 \(X\) 是偶数则 \(f(X)=1\), 奇数则 \(f(X)=0\)\(p_x(X)\) 是在 \(\{1, . . . , 10\}\) 上的均匀分布,\(p_y(1) = \sum_{x \in \{2,4,6,8,10\}} p_x(x) = 0.5\) ,并且 \(p_y(0) = 0.5\). 注意这个例子中函数 \(f\) 是多对一的函数。

如果 \(X\) 是连续型随机变量,就可以利用公式 2.83。因为 \(p_x(x)\) 是一个概率密度,而不是概率质量,无法用密度累加方式计算,只能采用累积分布函数 \(P\) 来计算:

\[ P_y (y)\overset\triangle {=} P (Y\le y)=P (f (X)\le y)=P (X\in\{x|f (x)\le y\}) ~~\tag {2.84}\label {eqn:2.84} \]

可以通过对累积分布函数的微分求导得到 \(y\) 的概率密度函数。

如果 \(f\) 是单调函数和可逆函数,则可以改写为:

\[ P_y (y)\overset\triangle {=} P (Y\le y)=P (X\le f^{-1}(y))=P_x (f^{-1}(y)) ~~\tag {2.85}\label {eqn:2.85} \]

求导可得:

\[ p_y (y)\overset\triangle {=} \frac {d}{dy} P_y (y)=\frac {d}{dy} P_x (f^{-1}(y)=\frac {dx}{dy}\frac {d}{dx} P_x(x)=\frac {dx}{dy} p_x (x) ~~\tag {2.86}\label {eqn:2.86} \]

其中 \(x=f^{-1}(y)\) , 可以把 \(dx\) 看作是 \(x\) -空间中的一种体积测度;类似的把 \(dy\) 当作对 \(y\) -空间的体积测度。则 \(\frac {dx}{dy}\) 测量了体积的变化。由于符号无关紧要,所以取绝对值来得到通用表达:

\[ p_y (y)=p_x (x)|\frac {dx}{dy}| ~~\tag {2.87}\label {eqn:2.87} \]

这也叫变量转换公式 (change of variables formula). 按照下面的思路来理解可能更容易。落在区间 \((x, x+\delta x)\) 的观测被变换到区间 \((y, y+\delta y)\) , 其中 \(p_x (x)\delta x \approx p_y (y)\delta y\)。 因此 \(p_y (y) \approx p_x(x)|\frac {\delta x}{\delta y}|\). 例如,假如有个随机变量 \(X \sim U(-1,1)\) , 且 \(Y=X^2\) , 那么则有 \(p_y (y)=\frac{1}{2}y^{-\frac{1}{2}\). 具体看练习 2.10.

2.6.2.1 变量的多重变化

前面的结果可推广到多元分布上。设 \(f\) 是一个函数,从 \(\mathbb{R}^n\) 映射到 \(\mathbb{R}^n\), 设 \( \mathbf{y}=f (\mathbf{x})\) 。则该函数的雅可比矩阵(Jacobian matrix) \(J\) 为:

\[\begin{split} \begin{align*} J_{\mathbf{x}\rightarrow \mathbf{y} } \triangleq \frac {\partial (y_1,...,y_n)}{\partial (x_1,...,x_n)}\overset\triangle {=} \begin {pmatrix} \frac {\partial y_1}{\partial x_1} & ...& \frac {\partial y_1}{\partial x_n} \\ ...&...&...\\ \frac {\partial y_n}{\partial x_1} &...&\frac {\partial y_n}{\partial x_n} \\ \end {pmatrix} ~~\tag {2.88}\label {eqn:2.88} \end{align*} \end{split}\]

矩阵的行列式 \(|\text{det}J|\) 表示在运行函数 \(f\) 时一个单位的超立方体上的体积变化。

如果 \(f\) 是一个可逆映射, 就可以用逆映射 \(\mathbf{y}\rightarrow \mathbf{x}\) 的雅可比矩阵来定义变换后随机变量的概率密度函数 (pdf):

\[ p_y (\mathbf{y})=p_x (\mathbf{x})|det (\frac {\partial \mathbf{x}}{\partial \mathbf{y}})|=p_x (\mathbf{x})|detJ_{\mathbf{y}\rightarrow \mathbf{x}}| ~~\tag {2.89}\label {eqn:2.89} \]

在练习 4.5, 你就要用到上面这个公式来推导一个多元正态分布的归一化常数。

举个简单例子,假如要把一个概率密度函数从笛卡尔坐标系的 \(\mathbf{x}=(x_1,x_2)\) 转换到一个极坐标系 \(y=(r,\theta )\), 有对应关系:\(x_1=r \cos \theta,x_2=r \sin \theta\). 这样则有雅可比矩阵如下:

\[\begin{split} \begin{align*} J_{\mathbf{y}\rightarrow \mathbf{x} }= \begin {pmatrix} \frac {\partial x_1}{\partial r} &\frac {\partial x_1}{\partial \theta} \\ \frac {\partial x_2}{\partial r} &\frac {\partial x_2}{\partial \theta} \\ \end {pmatrix} = \begin {pmatrix} \cos \theta & -r \sin \theta \\ \sin \theta & r\cos \theta\\ \end {pmatrix} ~~\tag {2.90}\label {eqn:2.90} \end{align*} \end{split}\]

矩阵 \(J\) 的行列式为:

\[ |det J|=|r\cos^2\theta+r\sin^2\theta|=|r| ~~\tag {2.91}\label {eqn:2.91} \]

因此:

\[ p_y (y)=p_x (x)|det J| ~~\tag {2.92}\label {eqn:2.92} \]
\[ p_{r,\theta}(r,\theta)=p_{x_1,x_2}(x_1,x_2) r=p_{x_1,x_2}(r\cos\theta,r\sin\theta) r ~~\tag {2.93}\label {eqn:2.93} \]

以几何角度来看,可以参考图 2.16, 其中的阴影部分面积可以用如下公式计算:

\[ P (r \le R \le r + dr, \theta \le \Theta \le \theta + d\theta ) = p_{r,\theta} (r, \theta ) drd\theta ~~\tag {2.94}\label {eqn:2.94} \]

在这个限制范围内,这也就等于阴影中心部分的密度 \(p (r,\theta)\) 乘以阴影部分的面积 \(rdrd\theta\)。 因此则有:

\[ p_{r,\theta} (r, \theta ) drd\theta= p_{x_1,x_2}(r\cos\theta,r\sin\theta) rdrd\theta ~~\tag {2.95}\label {eqn:2.95} \]

图 2.16 变量从极坐标到笛卡尔坐标的变化。阴影块的面积是 \(rdrdθ\)

图 2.17 中心极限定理。 我们绘制了一个直方图 \(\frac{1}{N}\sum^N_{i=1}x_{ij}\) , 其中 \(x_{ij} \sim \text{Beta}(1,5)\) ,对于 \(j = 1 : 10000\)。当 \(N \to \infty\),分布趋于高斯分布。(a) \(N = 1\)。 (b) \(N = 5\)

2.6.3 中心极限定理

现在设想有 \(N\) 个概率密度函数为 \(p(x_i)\) (不一定是正态分布)的随机变量,所有随机变量的均值和方差都是 \(\mu\)\(\sigma^2\),再假设所有随机变量都是独立同分布的,令 \(S_N =\sum^N_{i=1} X_i\) 是所有随机变量的和。随着 \(N\) 的增大,这个和的分布会接近正态分布:

\[ p (S_N=s)=\frac {1}{\sqrt {2\pi N\sigma^2}}\exp (-\frac {(s-N\mu)^2}{2N\sigma^2}) ~~\tag {2.96}\label {eqn:2.96} \]

也就是量 \(Z_N\) 的分布:

\[ Z_N \overset\triangle {=} \frac {S_N-N_{\mu}}{\sigma\sqrt N} = \frac {\bar X-\mu}{\sigma/\sqrt N} ~~\tag {2.97}\label {eqn:2.97} \]

会收敛到标准正态分布了。其中 \(\bar X=\frac{1}{N} \sum^N_{i=1} x_i\) 为样本均值。

上述即为中心极限定理,更多内容参考 (Jaynes 2003, p222) 或者 (Rice 1995, p169).

图 2.17 即是一例,其中计算 \(\beta\) 分布变量均值,右图可见很快收敛到正态分布了。

2.7 蒙特卡罗近似方法

要计算一个随机变量的函数的分布,靠公式变换 通常还挺难的。有另外一个办法,简单又好用。首先就是从分布中生成 S 个样本,就把它们标为 \(x_1,...,x_S\). 生成样本有很多方法,对于高维度分布来说最流行的方法就是马尔科夫链蒙特卡罗方法 (Markov chain Monte Carlo, 缩写为 MCMC), 这部分内容在本书 24 章再行讲解。 还说这个例子,对分布函数 \(f (X)\) 使用经验分布 (empirical distribution)\(\{f (x_s )\}^S_{s=1}\) 来进行近似。这就叫蒙特卡洛近似 (Monte Carlo approximation), 之所以用这个名字是因为欧洲的知名赌城。这种方法首先是在统计物理里面应用发展起来的,确切来说还是在原子弹研究过程中,不过现在已经广泛应用在统计和机器学习领域里面了。 此处查看原书图 2.18

应用蒙特卡罗方法,可以对任意的随机变量的函数进行近似估计。先简单取一些样本,然后计算这些样本的函数的算术平均值 (arithmetic mean). 这个过程如下所示:

\[ E [f (X)]=\int f (xp (x) dx\approx \frac1S\sum^S_{s=1} f (x_s) ~~\tag {2.98}\label {eqn:2.98} \]

上式中 \(x_s \sim p (X)\). 这就叫做蒙特卡罗积分 (Monte Carlo integration), 相比数值积分 (numerical integration) 的一个优势就是在蒙特卡罗积分中只在具有不可忽略概率的地方进行评估计算,而数值积分会对固定网格范围内的所有点的函数进行评估计算。 通过调整函数 \(f ()\), 就能对很多有用的变量进行估计,比如:

  • \(\bar x =\frac 1S \sum^S_{s=1} x_s\rightarrow E [X]\)

  • \(\frac 1 S\sum^S_{s=1}(x_s-\bar x)^2\rightarrow var [X]\)

  • \(\frac 1 S \# \{x_s \le c\}\rightarrow P (X\le c)\)

  • 中位数 (median)\(\{x_1,...,x_S\}\rightarrow median (X)\)

下面是一些例子,后面一些章节有更详细介绍。

2.7.1 样例:更改变量,使用 MC (蒙特卡罗)方法

在 2.6.2,我们讨论了如何分析计算随机变量函数的分布 \(y = f (x)\). 更简单的方法是使用蒙特卡罗方法估计。例如,若 \(x \sim Unif (−1, 1) , y = x^2\). 就可以这样估计 \(p (y)\): 从 \(p (x)\) 中去多次取样,取平方,计算得到的经验分布。如图 2.18 所示。后文中还要广泛应用这个方法。参考图 5.2.

此处查看原书图 2.19

2.7.2 样例:估计圆周率 \(\pi\), 使用蒙特卡罗积分

蒙特卡罗方法还可以有很多种用法,不仅仅是统计学领域。例如我们可以用这个方法来估计圆周率 \(\pi\). 我们知道圆的面积公式可以利用圆周率和圆的半径 r 来计算,就是 \(\pi r^2\), 另外这个面积也等于下面这个定积分 (definite integral):

\[ I=\int _{-r}^r\int _{-r}^r\prod (x^2+y^2\le r^2) dxdy ~~\tag {2.99}\label {eqn:2.99} \]

因此有 \(\pi =I/(r^2)\). 然后就可以用蒙特卡罗积分来对此进行近似了。设 \(f (x,y) =\prod (x^2+y^2\le r^2)\) 是一个指示器函数 (indicator function), 只要点在圆内,则函数值为 1, 反之为 0, 然后设 \(p (x),p (y)\) 都是在闭区间 [-r,r] 上的均匀分布 (uniform distribution), 所以有 \(p (x)=p (y)=\frac {1}{2r}\) 这样则有:

\[\begin{split} \begin {align*} I &= (2r)(2r)\int\int f (x,y) p (x) p (y) dxdy ~~\tag {2.100}\label {eqn:2.100}\\ &= 4r^2 \int\int f (x,y) p (x) p (y) dxdy ~~\tag {2.101}\label {eqn:2.101}\\ &=4r^2\frac1S\sum^S_{s=1} f (x_s,y_s) ~~\tag {2.102}\label {eqn:2.102}\\ \end {align*} \end{split}\]

当标准差为 0.09 的时候,计算得到的圆周率为 \(\hat \pi =3.1416\), 参考本书 2.7.3 就知道什么是标准差了。接受 / 拒绝的点如图 2.19 中所示。 此处查看原书图 2.20

2.7.3 蒙特卡罗方法的精确度

随着取样规模的增加,蒙特卡罗方法的精度就会提高,如图 2.20 所示,在图上部是从一个高斯分布中取样的直方图,底下的两个图使用了核密度估计 (kernel density estimate, 参考本书 14.7.2) 得到的光滑曲线。这种光滑分布函数在密集网格点上进行评估和投图。这里要注意一点,光滑操作只是为了投图看而已,蒙特卡罗方法估计的过程根本用不着光滑。 如果我们知道了均值的确切形式,即 \(\mu =\mathbb{E}[f (X)]\), 然后蒙特卡罗方法近似得到的是 \(\hat\mu\), 那么对于独立取样则有:

\[ (\hat\mu -\mu )\rightarrow N (0,\frac {\sigma^2 }{S}) ~~\tag {2.103}\label {eqn:2.103} \]

其中:

\[ \sigma^2=var [f (X)]=E [f (X)^2]-E [f (X)]^2 ~~\tag {2.104}\label {eqn:2.104} \]

这是由中心极限定理决定的。当然了,上式中的 \(\sigma^2\) 是位置的,但也可以用蒙特卡罗方法来估计出来:

\[ \hat\sigma^2= \frac1S \sum^S_{s=1}(f (x_s)-\hat\mu)^2 ~~\tag {2.105}\label {eqn:2.105} \]

然后则有:

\[ P\{\mu-1.96\frac {\hat \sigma}{\sqrt S}\le \hat\mu \le \mu +1.96\frac {\hat \sigma}{\sqrt S}\}\approx 0.95 ~~\tag {2.106}\label {eqn:2.106} \]

上式中的 \(\frac {\hat \sigma}{\sqrt S} \) 就叫做数值标准差或者经验标准差 (numerical or empirical standard error), 这个量是对我们对 \(\mu\) 估计精度的估计。具体信息查看本书 6.2 有更多讲解。 如果我们希望得到的答案 \(\pm \epsilon\) 范围内的概率至少为 95%, 那就要保证取样数目 S 满足条件 \(1.96\sqrt {\hat\sigma^2/S}\le \epsilon\), 这里的 1.96 可以粗略用 2 替代,这样就得到了 \(S\geq \frac {4 \hat\sigma^2}{\epsilon^2}\).

2.8 信息理论

信息理论 (information theory) 关注的是以紧凑形式进行数据呈现(这种紧凑形式也被称为数据压缩 (data compression) 或者源编码 (source coding)), 以及以能健壮应对错误的方式进行传输和存储(这个过程也叫做纠错 (error correction) 或者信道编码 (channel coding)). 第一眼看上去好像这和概率论以及机器学习没什么关系,不过实际是有联系的。首先,对数据进行紧凑表达需要给高概率的字符串赋予短编码字,而将长编码字留给低概率字符串。就好比自然语言中,特别常用的词汇都往往比少见的词汇短很多,比如冠词 a/the 明显比闪锌矿 sphalerite 短很多。另外,在噪音频道上进行信息解码也需要对人发送的不同信息建立一个良好的概率模型。这就都需要一个能够预测数据可能性的模型,这也是机器学习里面的一个核心问题,关于信息理论和机器学习之间关系的更多内容请参考 (MacKay 2003).

显然这里不可能说太多太深关于信息理论的内容,有兴趣的话去看 (Cover and Thomas 2006). 这里也就是介绍本书中要用到的一些基础概念了。

2.8.1 信息熵

随机变量 X 服从分布 p, 这个随机变量的熵 (entropy) 则表示为 \(H (X)\) 或者 \(H (p)\), 这是对随机变量不确定性的一个衡量。对于一个有 K 个状态的离散型随机变量来说,其信息熵定义如下:

\[ H (X)\overset\triangle {=}-\sum^K_{k=1} p (X=k)\log_2p (X=k) ~~\tag {2.107}\label {eqn:2.107} \]

通常都用 2 作为对数底数,这样单位就是 bit (这个是 binary digits 的缩写). 如果用自然底数 e, 单位就叫做 nats 了。 举个例子,\(X\in \{1,...,5\}\), 柱状分布 (histogram distribution), 概率 \(p=[0.25,0.25,0.2,0.15,0.15]\), 利用上面的公式计算得到 \(H =2.2855\). 熵最大的离散分布就是均匀分布,可以参考本书 9.2.6. 因为对于一个 K 元 (K-ary) 随机变量,如果 \(p (x = k) = 1/K\), 则信息熵最大,这时候的熵为 \(H (X)=\log_2K\). 熵最小的分布就是所有概率质量都在单一状态的 \(\delta\) 分布,这时候熵为 0, 因为只有一个状态有概率,没有任何不确定性。 在图 2.5 当中对 DNA 序列进行了投图,每一列的高度定义为 \(2-H\), 其中的 H 就是这个分布的熵,2 是最大可能熵 (maximum possible entropy). 因此高度为 0 的就表示均匀分布,而高度为 2 就对应着确定性分布 (deterministic distribution).

此处查看原书图 2.21

对于二值化随机变量的特例,\(X\in\{0,1\}\), 则有 \(p (X=1)=\theta, p (X=0)=1-\theta\), 这样熵为:

\[\begin{split} \begin {align*} H (X)&= -[p (X=1)\log_2p (X=1)+p (X=0)\log_2p (X=0)] ~~\tag {2.108}\label {eqn:2.108}\\ &=-[\theta\log_2\theta+(1-\theta)\log_2 (1-\theta)] ~~\tag {2.109}\label {eqn:2.109} \end {align*} \end{split}\]

这也叫做二值熵函数 (binary entropy function), 也写作 \(H (\theta)\), 如图 2.21 所示,可见当 \(\theta=0.5\) 的时候熵值最大为 1, 这时候是均匀分布了。

2.8.2 KL 散度

KL 散度 (Kullback-Leibler divergence), 也称相对熵 (relative entropy), 可以用来衡量 p 和 q 两个概率分布的差异性 (dissimilarity). 定义如下:

\[ KL (p||q)\overset\triangle {=}\sum^K_{k=1} p_k\log\frac {p_k}{q_k} ~~\tag {2.110}\label {eqn:2.110} \]

上式中的求和也可以用概率密度函数的积分来替代。就可以写成:

\[ KL (p||q)=\sum_kp_k\log p_k - \sum_kp_k\log q_k =-H (p)+H (p,q) ~~\tag {2.111}\label {eqn:2.111} \]

上式中的 \(H (p,q)\) 就叫做交叉熵 (cross entropy):

\[ H (p,q)\overset\triangle {=}-\sum_kp_k\log q_k ~~\tag {2.112}\label {eqn:2.112} \]

参考 (Cover and Thomas 2006) 可以证明,当使用模型 q 来定义编码本 (codebook) 的时候,来自分布 p 的待编码数据的平均比特数 (average number of bits) 就是交叉熵。正规熵 (regular entropy)\(H (p)=H (p,p)\), 参考本书 2.8.1 的定义,也就是使用真实模型时候的比特数期望值,因此 KL 散度也就是不同概率分布之间的不同的量度。换个说法,KL 散度就是要对数据编码所需要的额外比特 (extra bits) 的平均数,因为这时候用分布 q 来对数据进行编码,而不是使用分布 p.

既然是额外的比特数,这种表述就很明显说明这个 KL 散度应该是非负的,即 \(KL (p||q)\ge 0\), 等于 0 则意味着两个分布相等,即 \(p=q\). 接下来对此进行一下证明。

定理 2.8.1 信息不等式 (Information inequality)

\(KL (p||q)\ge 0\) 当且仅当 \(p=q\) 的时候,KL 散度为 0.

证明

要证明这个定理,需要用到詹森不等式 (Jensen’s inequality). 这个不等式是说,对于任意的凸函数 (convex function) f, 有以下关系:

\[ f (\sum^n_{i=1}\lambda_i (x_i)) \le \sum^n_{i=1}\lambda_i f (x_i) ~~\tag {2.113}\label {eqn:2.113} \]

其中 \(\lambda_i \ge 0,\sum^n_{i=1}\lambda_i=1\). 由于凸函数的定义,对于 n=2 的时候很显然,对于 n>2 的情况也可以归纳证明 (proved by induction).

对定理的证明参考了 (Cover and Thomas 2006, p28). 设 \(A=\{x:p (x)>0\}\)\( p (x)\) 的支撑集合 (support, 译者注:纯白或许就当做定义域理解好了). 则有:

\[\begin{split} \begin {align*} -KL (p||q)& = -\sum_{x\in A} p (x)\log \frac {p (x)}{q (x)} = \sum_{x\in A} p (x)\log \frac {q (x)}{p (x)} ~~\tag {2.114}\label {eqn:2.114}\\ & \le \sum_{x\in A} p (x)\frac {q (x)}{p (x)} = \log \sum_{x\in A} q (x) ~~\tag {2.115}\label {eqn:2.115}\\ & \le \log \sum_{x\in \chi} q (x) =\log1 =0 ~~\tag {2.116}\label {eqn:2.116}\\ \end {align*} \end{split}\]

当 上面第一个不等式是应用了詹森不等式。因为 \(\log (x)\) 是个严格凸函数,所以在等式 2.115 里面,当且仅当对于某些 c 使 \(p (x)=cq (x)\) 成立的时候,等量关系成立。等式 2.116 中的等量关系当且仅当 \(\sum_{x \in A } q (x)=\sum_{x\in \chi } q (x)=1\) 的时候成立,这时候 c=1. 所以对于所有的 x 来说,当且仅当 \(p (x)=q (x),KL (p||q)=0\).

证明完毕。 这个结果的一个重要推论就是 * 具有最大熵的离散分布就是均匀分布 *. 更确切地说,\(H (X)\le \log|\chi |\),\(|\chi |\) 是 X 的状态数,当且仅当 \(p (x)\) 是均匀分布的时候等号成立。设 \(u (x)=1/|\chi |\), 则有:

\[\begin{split} \begin {align*} 0&\le KL (p||u)= \sum_xp (x)\log\frac {p (x)}{u (x)} ~~\tag {2.117}\label {eqn:2.117}\\ &=\sum_xp (x)\log p (x)-\sum_xp (x)\log u (x)=-H (X)+\log|\chi| ~~\tag {2.118}\label {eqn:2.118}\\ \end {align*} \end{split}\]

上面这个就是公式化的拉普拉斯不充分理由原则 (Laplace’s principle of insufficient reason), 说的是在没理由优先选择某个分布的时候,优先选择均匀分布 (uniform distribution). 关于如何建立满足特定约束条件 (certain constraints) 的分布可以阅读本书 9.2.6, 其他方面尽可能最小化 (as least-commital as possible).(正态分布满足一阶和二阶矩约束条件,但其他方面有最大熵。)

2.8.3 信息量 (Mutual information)

设有两个随机变量 X 和 Y. 如果我们想知道一个变量告诉我们关于另一个变量的多少信息。就可以计算相关系数 (correlation coefficient) 了,可是相关系数只适用于实数值的随机变量。另外相关系数对不相关程度的衡量作用也很有限,如图 2.12 所示。所以更常用的方法是对比联合分布 (joint distribution)\(p (X, Y)\) 和因式分布 (factored distribution)\(p (X) p (Y)\) 的相关性。这就叫互信息量 (mutual information) 或者简写做 MI, 定义如下:

\[ I (X;Y)\overset\triangle {=} KL (p (X,Y)||p (X) p (Y))=\sum_x\sum_yp (x,y)\log\frac {p (x,y)}{p (x) p (y)} ~~\tag {2.119}\label {eqn:2.119} \]

\(I (X;Y)\ge0\) 的等号当且仅当 \(p (X,Y)=p (X) p (Y)\) 的时候成立。也就是如果两个变量相互独立,则互信息量 MI 为 0. 为了深入理解 MI 这个量的含义,咱们用联合和条件熵的方式来重新表述一下。参考练习 2.12 可以得到上面的表达式等价于下列形式:

\[ I (X;Y)=H (X)-H (X|Y)=H (Y)-H (Y|X) ~~\tag {2.120}\label {eqn:2.120} \]

上式中的 \(H (Y|X)\) 就是条件熵 (conditional entropy) 定义为 \(H (Y|X)=\sum_xp (x) H (Y|X=x)\). 这样就可以把 X 和 Y 之间的互信息量 MI 理解成在观测了 Y 之后对 X 的不确定性的降低,或者反过来就是观测了 X 后对 Y 不确定性的降低。本书后面一些内容中还会用到这个概念。参考 2.13 和 2.14 来阅读互信息量 MI 和相关系数之间的联系。 另外一个和互信息量 MI 有很密切联系的量是点互信息量 (pointwise mutual information, 缩写为 PMI), 对于两个事件(而不是随机变量) x 和 y, 其点互信息量 PMI 定义如下:

\[ PMI (x,y)\overset\triangle {=} \log\frac {p (x,y)}{p (x) p (y)}= \log\frac {p (x|y)}{p (x)}= \log\frac {p (y|x)}{p (y)} ~~\tag {2.121}\label {eqn:2.121} \]

这个量衡量的是与偶发事件相比,这些事件之间的差异。很明显 X 和 Y 的互信息量 MI 就是点互信息量 PMI 的期望值。所以就可以把点互信息量 PMI 的表达式写为:

\[ PMI (x,y)= \log\frac {p (x|y)}{p (x)}= \log\frac {p (y|x)}{p (y)} ~~\tag {2.122}\label {eqn:2.122} \]

这个量是通过将先验 (prior) 的 \(p (x)\) 更新到后验 (posterior) 的 \(p (x|y)\) 得到的,也可以是将先验的 \(p (y)\) 更新到后验的 \(p (y|x)\) 得到。

2.8.3.1 连续随机变量的互信息量

上一节中的互信息量 MI 定义是针对离散型随机变量的。对于连续随机变量,可以先对其进行离散化 (discretize) 或者量子化 (quantize), 具体方法可以使将每个随机变量归类到一个区间里面,将变量的变化范围划分出来,然后计算每一段的小区间中的分布数量 (Scott 1979). 然后就可以利用上文的方法公式来计算互信息量 MI 了(代码参考 PMTK3 的 mutualInfoAllPairsMixed, 样例可以参考 miMixedDemo).

此处查看原书图 2.22

然而很不幸,分成多少个小区间,以及小区间边界的位置,都可能对计算结果有很大影响。一种解决方法就是直接对互信息量 MI 进行估计,而不去先进行密度估计 (Learned-Miller 2004)). 另一种办法是尝试很多不同的小区间规模和位置,然后计算得到的最大互信息量 MI. 经过适当的标准化之后,这个统计量就被称为最大信息系数 (maximal information coefficient, 缩写为 MIC)(Reshed et al. 2011). 更确切来说定义如下所示:

\[ m (x,y)=\frac {\max_{G\in G (x,y)} $i$ (X (G);Y (G))}{\log\min (x,y)} ~~\tag {2.123}\label {eqn:2.123} \]

上式中的 \(G (x, y)\) 是一个规模为 \( x\times y\) 的二维网状集合,而 \(X (G),Y (G)\) 表示的是将变量在这个网格上进行离散化得到的结果。在区间位置 (bin locations) 上最大化的过程可以通过使用动态编程 (dynamic programming) 来有效进行 (Reshed et al. 2011). 这样定义了连续型变量互信息量 MIC 如下:

\[ MIC\overset\triangle {=} \max_{x,y:xy<B} m (x,y) ~~\tag {2.124}\label {eqn:2.124} \]

上式中的 B 是一个与取样规模相关的约束条件,用于约束能使用且能可靠估计分布的区间个数。((Reshed et al. 2011) 建议的是 \(B = N^{0.6}\). 显然 MIC 处于区间 [0, 1] 中,其中 - 表示两个变量没关系,而 1 表示二者有无噪音的相关性 (noise-free relationship), 这种相关性可以是任意形式的,不仅仅是线性相关。 图 2.22 给出的是一个实例。其中的数据集包含了 357 个变量,衡量一系列的社会 / 经济 / 健康 / 政治指标,由世界健康组织 WHO 手机。左边的途中看到了对于 65,566 个变量对的相关系数 (CC) 与互信息量 (MIC) 的关系图。有图则投了一些特定变量对的散点图,其中包括了:

  • C 图中的是 CC 和 MIC 都低,相应的散点图很明显表明了这两组变量之间没有关系:因伤致死比例和人群中牙医密度。

  • D 图和 H 图中是 CC 和 MIC 都高,呈现近乎线性的相关性。

  • E/F/G 这三个图都是低 CC 高 MIC. 这是因为这些变量之间的关系是非线性的,例如在 E 图和 F 图中,都是非函数对应关系,比如可能是一对多的对应关系。 总的来说,MIC 这个统计量是基于互信息量的,可以用于发现变量之间的有意义的关系,而这些关系可能是那些简单的统计量,比如相关系数之类无法反应的。由于这个优势,MIC 也被称作是 21 世纪的相关性衡量变量 “a correlation for the 21st century” (Speed 2011).

练习 2