PRML 概率分布

本文地址：https://www.cnblogs.com/faranten/p/15917369.html
转载请注明作者与出处

1 二元变量

1.1 伯努利分布与二项分布

​ 考虑一个最基本的试验：抛硬币试验。在一次实验中只有两个结果，即正面与反面，用随机变量x=1来表示抛掷硬币得到的是正面，x=0来表示抛掷硬币得到的是反面，且先验地猜测得到正面的概率是μ，那么

p(x=1|μ)=μp(x=0|μ)=1−μ

那么x的概率分布可以写作

Bern(x|μ)=μx(1−μ)1−x

这称为二元变量的伯努利分布（Bernoulli distribution），容易得到

E(x|μ)=μvar(x|μ)=μ(1−μ)

​ 如果将试验次数增多到N次，实验数据集为D，并以随机变量x记得到硬币正面的次数，由于各次试验是相互独立的，因此似然函数为

p(D|μ)=N∏n=1p(xn|μ)=N∏n=1μxn(1−μ)1−xn

对于频率主义来说，此时的实验数据集是确定的，因此确定参数μ的方式就是最大化上述形如μa(1−μ)b的似然函数，使用极大似然方法，对数化上述似然函数得到lnp(D|μ)=∑Nn=1{xnlnμ+(1−xn)ln(1−μ)}，进而解得上式中μML=1N∑Nn=1xn，该值称为样本均值（sample mean），其实该值就是N次试验中得到正面的比例mN。值得一提的是，该对数似然函数仅与∑Nn=1xn有关，因此可称该量为充分统计量（sufficient statistic）。后面将看到，二项分布作为伯努利分布的一般化形式，所以只需要对二项分布使用贝叶斯方法，就涵盖了贝叶斯方法处理伯努利分布的情况。

​ 现在考虑N次试验中得到正面的次数m，若将此作为随机变量，则它的分布为

Bin(m|N,μ)=(Nm)μm(1−μ)N−m,(Nm)≡N!(N−m)!m!

该分布称为二项分布（binomial distribution），二项分布可以视作伯努利分布在试验次数上的推广，容易得到

E(m|N,μ)=Nμvar(m|N,μ)=Nμ(1−μ)

​ 在二项分布中，如果给定数据集，在频率方法下用极大似然方法求出μ的估计值仍然是N次试验中得到正面的比例mN，当数据集规模较小的时候，这种思路很容易导致过拟合，这说明了贝叶斯方法的必要性。

1.2 Beta分布

​ 在二项分布中，如果令N=1，那么二项分布就变成了伯努利分布，因此可以认为二项分布是伯努利分布更加一般的形式，所以接下来只讨论二项分布而忽略伯努利分布的情况。对于贝叶斯主义而言，二项分布中的μ是随机变量，我们应该用训练集来找到μ的尽可能精确的分布，由于

p(μ|D)=p(D|μ)p(μ)p(D)

其中p(D)为常数N，而似然函数为

p(D|μ)=N∏n=1p(xn|μ)=N∏n−1μxn(1−μ)1−xn=μm(1−μ)N−m

其中m为得到正面的次数。如果先验分布和后验分布具有相似的函数形式，则这种性质称为共轭性（conjugacy）（这保证了顺序学习过程将会一直进行下去）。观察到似然函数仅与μ与(1−μ)两个因子的幂指数成正比，所以如果我们选择一个正比于μ与(1−μ)​两个因子的幂指数，自然就能保证共轭性。所以选择如下形式的先验分布

Beta(μ|a,b)=Γ(a+b)Γ(a)Γ(b)μa−1(1−μ)b−1

这称为Beta分布，其中参数a与b决定了分布的形态，因此是超参数。容易得到

E(μ|a,b)=aa+bvar(μ|a,b)=ab(a+b)2(a+b+1)

​ 那么，参数μ的后验分布p(μ|m,l,a,b)满足

p(μ|m,l,a,b)∝p(m,l|μ,a,b)⋅p(μ|a,b)∝μm(1−μ)l⋅μa−1(1−μ)b−1=μm+a−1(1−μ)l+b−1

其中l=N−m，和Beta分布的标准形式相比，很快得到归一化系数，于是就有

p(μ|m,l,a,b)=Γ(m+l+a+b)Γ(m+a)Γ(l+b)μm+a−1(1−μ)l+b−1

直观来看：其中m和l是训练集中正面次数和反面次数，a和b是先验知道的正面次数和反面次数，在试验的时候可以简单地认为a的值变大了m、b的值变大了l，可将超参数a与b分别看成是x=1和x=0的有效观测数（effective number of observation）。但是a和b具有更一般化的含义，不仅仅是整数。随机变量μ​的期望和方差分别为

E(μ|m,l,a,b)=m+a(m+a)+(l+b)var(μ|m,l,a,b)=(m+a)(l+b)((m+a)+(l+b))2((m+a)+(l+b)+1)

​ 上述内容暗示了学习过程中的顺序（sequential）方法是合理的，即每次将先验分布乘上似然函数，再进行归一化（找到合适的归一化参数）便得到了后验分布，这个后验分布将在下一次学习过程中扮演先验分布的角色。并且，随着试验次数N→∞，m和l都将趋于正无穷，此时的E(μ)和var(μ)都将趋于各自的极大似然估计，并且与先验的参数a和b无关。

2 多项式变量

2.1 范畴分布与多项式分布

​ 二元变量只能用来描述某两种取值的试验，现在给出一种新的变量：多项式变量（multinomial variable），可以用来描述具有多种离散情况的变量，比如

x=(0,0,1,0,0,0)T

描述的是一个具有六种离散状态的变量，并且此时该变量为第三种状态。多项式变量满足∑Kk=1xk=1。现用参数μk来描述第k个变量xk=1的概率，则有

p(x|μμ)=K∏k=1μxkk

该分布称为范畴分布（multinoulli distribution）或者分类分布（categotical distribution），该分布可以视作伯努利分布在试验“维数”上面的推广，其中μμ=(μ1,μ2,⋯,μk)T，显然μk≥0并且∑Kk=1μk=1。可以得到

E(x|μμ)=(μ1,μ2,⋯,μK)Tvar(x|μμ)=(μ1(1−μ1),μ2(1−μ2),⋯,μK(1−μK))Tcov(x|μμ)=diag(μ1,⋯,μK)−μμμμT

​ 在之前关于伯努利分布和二项分布的内容中，我们将伯努利分布视为单一试验中二元变量的概率分布情况，而将二项分布视为多次试验中二元变量的概率分布情况。现在，我们刚讨论完单一试验中多项式变量的概率分布情况，自然要考虑多次实验中多项式变量的概率分布情况。对于给定的数据集D={x1,x2,⋯,xN}而言，对应的似然函数的形式为

p(D|μμ)=N∏n=1K∏k=1μxnkk=K∏k=1μ∑Nn=1xnkk=K∏k=1μmkk

其中mk描述了在数据集D中取第k种状态的数据点的数量，即xk=1的次数，且这k个值是该似然函数的充分统计量。对此似然函数进行对数化处理，用极大似然方法，注意到约束条件∑Kk=1μk=1，可以解得参数μμ的极大似然估计为μMLk=mkN。

​ 现在考虑每个状态的观测数量在参数μμ和总观测数量N条件下的分布，若将每个状态的观测数量作为一组随机变量，则它们的联合分布为

Mult(m1,m2,⋯,mK|μμ,N)=(Nm1m2⋯mK)K∏k=1μmkk

该分布称为多项式分布（multinomial distribution），约束条件为∑Kk=1mk=N。多项式分布可以视作范畴分布在试验次数上的推广，其中归一化系数就是在PRML 基础知识 6.1节提到的“乘数”概念，具体意义是将N个物体分成大小为m1,m2,⋯,mK的K​组的方案总数。可以得到

E(m1,m2,⋯,mK|μμ,N)=N(μ1,μ2,⋯,μK)Tvar(m1,m2,⋯,mK|μμ,N)=N(μ1(1−μ1),μ2(1−μ2),⋯,μK(1−μK))Tcov(m1,m2,⋯,mK|μμ,N)=Ndiag(μ1,⋯,μK)−NμμμμT

2.2 Dirichlet分布

​ 现在我们考虑多项式分布的先验分布的形式。考虑到多项式分布正比于一系列参数μk的幂指数、或者统一起来说正比于参数μμ中每个元素各自的幂指数，因此为了保证先验分布和后验分布的共轭性，多项式分布的先验分布的形式为

Dir(μμ|αα)=Γ(α0)Γ(α1)Γ(α2)⋯Γ(αK)K∏k=1μαk−1k

其中∑Kk=1μk=1，αα=(α1,α2,⋯,αK)T且α0=∑Kk=1αk，该分布被称为Dirichlet分布。可以得到

E(μμ|αα)=(α1α0,α2α0,⋯,αKα0)T=(~α1,~α2,⋯,~αK)Tvar(μμ|αα)=(α1(α0−α1)α20(α0+1),α2(α0−α2)α20(α0+1),⋯,αK(α0−αK)α20(α0+1))T=(~α1(1−~α1)α0+1,~α2(1−~α2)α0+1,⋯,~αK(1−~αK)α0+1)Tcov(μμ|αα)=1α20(α0+1)diag(α0,α0,⋯,α0)⋅diag(α1,α2,⋯,αK)−μμμμT=1α0+1diag(~α1,~α2,⋯,~αK)−μμμμT

​ 那么，参数μμ的后验分布p(μμ|D,αα)满足

p(μμ|D,αα)∝p(D|μμ)⋅p(μμ|αα)∝K∏k=1μmkk⋅K∏k=1μαk−1=K∏k=1μmk+αk−1k

和Dirichlet分布的标准形式相比，很快得到归一化系数，于是就有

p(μμ|D,αα)=Γ(α0+N)Γ(α1+m1)Γ(α2+m2)⋯Γ(αK+mK)K∏k=1μmk+αk−1k

此时可以分析各项参数实际表示的含义，其含义与Beta分布中各参数的直观解释类似，此处不再讨论。

3 高斯分布

​ 先给出一维变量x和D维变量x情况下的高斯分布通用形式

N(x|μ,σ2)=1(2πσ2)1/2exp{−12σ2(x−μ)2}N(x|μμ,Σ)=1(2π)D/21|Σ|1/2exp{−12(x−μμ)TΣ−1(x−μμ)}

其中Σ是一个D×D的协方差矩阵。高斯分布是十分重要的，由中心极限定理知道，现实生活中很多情形都会推导出高斯分布。下面将从矩阵角度对熟知的高斯分布的一些性质进行重新推导，这不仅有利于将高斯分布从一元情形推广到多元情形，而且有利于理解后续章节的关键概念。

3.1 高斯分布的矩阵视角

​ 从上面多元形式的高斯分布可以看出，多元高斯分布中对x的依赖是通过二次型

Δ2=(x−μμ)TΣ−1(x−μμ)

实现的，Δ被称为马氏距离（Mahalanobis distance），当Σ为单位矩阵时，马氏距离就退化为欧氏距离。如果在一个关于x的空间中，马氏距离是常数，那么此时的多元高斯分布亦是常数。

​ 对于矩阵Σ的形式，我们不妨设其为对称矩阵，这是因为任何非对称项都会从多元高斯分布的指数项中消失，下面给出证明。记Δ2=(x−μ)TA(x−μ)，其中A=Σ−1，接着令A=12(A+AT)+12(A−AT)以及B=12(A+AT)，C=12(A−AT)，那么矩阵B就是对称矩阵（即有bij=bji），而矩阵C是反对称矩阵（即有cij=−cji）且A=B+C。现在将Δ2重新写作

Δ2=(x−μμ)T(B+C)(x−μμ)=(x−μμ)TB(x−μμ)+(x−μμ)TC(x−μμ)

现在如果能证明(x−μμ)TC(x−μμ)=0，则说明了任何非对称项都会从多元高斯分布的指数项中消失，事实证明确实如此，推导过程如下

(x−μμ)TC(x−μμ)=(x1−μ1,⋯,xD−μD)⋅(c11⋯c1D⋮⋮cD1⋯cDD)⋅(x1−μ1⋮xD−μD)=D∑i=1D∑j=1cij(xi−μi)(xj−μj)=D∑i=1D∑j=i+1(cij+cji)(xi−μi)(xj−μj)=0

这就说明了Δ2=(x−μμ)TB(x−μμ)，也就是任何非对称项都会从多元高斯分布的指数项中消失，在后续讨论中，我们默认矩阵Σ是对称矩阵。

​ 对于协方差矩阵Σ，考虑其特征方程Σui=λui，由于Σ是实对称矩阵，那么其特征值也是实数，且其特征向量可以被选为单位正交的（即uTiuj=δij），则协方差矩阵可以写作展开的形式

Σ=D∑i=1λiuiuTi

下面给出该结论的证明。首先构造矩阵U=(u1,⋯,uD)，即其中的每一列是特征向量，该矩阵满足UUT=I（或者等价条件UTU=I或U−1=UT），因此称为正交矩阵（orthogonal matrix），根据线性代数知识有ΣU=UΛ，那么

ΣU=UΛ=(u1,⋯,uD)(λ1⋱λD)=(λ1u1,⋯,λDuD)

因此

Σ=UΛU−1=UΛUT=(λ1u1,⋯,λDuD)(uT1⋮uTD)=D∑i=1λiuiuTi

因此Σ=∑Di=1λiuiuTi得证。而在ΣU=UΛ两侧同时取逆矩阵得到U−1Σ−1=Λ−1U−1，进而

Σ−1=(U−1)−1ΛU−1=UΛ−1UT=D∑i=11λiuiuTi

将此式代入二次型Δ2=(x−μμ)TΣ−1(x−μμ)中得到

Δ2=D∑i=11λi(x−μμ)TuiuTi(x−μμ)=D∑i=1y2iλi,yi=uTi(x−μμ)

如果二次型的值为常数且所有的特征值λi为正数，那么该二次型可以视为一个椭球面，椭球中心位于μμ处，但是该椭球的各轴可能不沿着ui方向，该式的意义就在于给出了一组新的基底y=U(x−μμ)，使得椭球在此坐标系下的中心位于(0,0,⋯,0)T处，且各轴的方向沿着ui方向，缩放因子为λ1/2i。在后面的“高斯分布的矩”一节中，将使用Dirac符号对此结论再给出一个推导。

​ 如果一个矩阵的特征值全为正数，则称此矩阵正定的（positive definite）；如果一个矩阵的特征值全为非负数，则称此矩阵是半正定的（positive semidefinite）。对于任意形式的高斯分布而言，其协方差矩阵必须是正定的。

​ 对于新的基底y=U(x−μμ)而言，高斯分布的形式又当如何？先定义Jacobian矩阵如下

J=(∂x1∂y1⋯∂x1∂yD⋮⋮∂xD∂y1⋯∂xD∂yD)=UT

且有|J|2=|UT|2=|UT||UT|=|UT||U|=|UTU|=|I|=1，以及|Σ|=∏Di=1λ1/2i，所以

p(y)=p(x)|J|=D∏i=11(2πλi)1/2exp{−y2i2λi}

该式也可以理解为p(y)=∏Di=1p(yi)。注意，在任意的坐标变换中，Jacobian矩阵都是十分重要的，它决定了坐标变换后整体的“放缩程度”。

3.2 题外话——Dirac符号

​ 作为拓展内容，现在简要介绍首先Dirac符号的概念。Dirac符号是构成现代量子力学形式体系的重要组成部分，由Dirac在1939年提出，Dirac将“括号（bracket）”一词一分为二，得到的两个单词分别代表左右矢量。对于复向量空间Cn，其中的任意一个元素（列向量）记为

|u⟩=(x1⋮xn)或|u⟩=(x1,⋯,xn)T

（其中xi∈C）称为右矢（ket vector or ket），可以定义右矢的加法和数乘运算、并且可以引入线性相关与线性无关等概念，这和线性代数是一样的，此处不再单独叙述。而任意的行向量

⟨α|=(α1,⋯,αn)

（其中αi∈C）称为左矢（bra vector or bra）。并且右矢和左矢的乘积⟨α|x⟩=∑ni=1αixi称为内积（inner product）。这样，如果保持左矢不变，通过内积运算便将Cn中的任意元素|x⟩映射到了C中的某个元素⟨α|x⟩，这可以看成是一种向量函数。如果在保持左矢不变的情况下，记此时的内积为一个从Cn到C的向量函数f，且f满足线性条件（linearity condition），即

f(c1|x⟩+c2|y⟩)=c1f(|x⟩)+c2f(|y⟩),∀|x⟩,|y⟩∈Cn,∀c1,c2∈C

那么称此时的f为线性函数（linear function）。容易验证，一个左矢和内积运算便构成了一个这样的f。并且，通过选取合适的左矢，可以表示出任何从Cn到C的线性函数，因为当一个线性函数作用在一个右矢上的时候，从内积运算的定义便知道它给出的结果f(|x⟩)总是形如⟨α|x⟩=∑ni=1αixi，那么⟨α|便可以代表此时的线性函数。在线性代数中我们学过：一组正交基（的线性组合）可以用来表示一个线性空间。对于向量空间来说，选取Cn中的单位正交向量组{|u1⟩,⋯,|un⟩}，则其满足⟨ui|uj⟩=δij，且Cn中任意向量|x⟩可以表示为∑ni=1ci|ui⟩，在等式两侧进行内积运算便得到

⟨uj|x⟩=n∑i=1ci⟨ej|ui⟩=n∑i=1ciδij=cj

即有cj=⟨uj|x⟩​，这给出了各个坐标的计算公式，将其代入|x⟩的表达式便有

|x⟩=n∑i=1ci|ui⟩=n∑i=1⟨ui|x⟩|ui⟩=n∑i=1|ui⟩⟨ui|x⟩=(n∑i=1|ui⟩⟨ui|)|x⟩

于是便得到一个重要的关系式

n∑i=1|ui⟩⟨ui|=I

该关系式称为完备性关系（completeness relation）。

​ 最后介绍一个概念：投影算子（projection operator）被定义为Pk=|uk⟩⟨uk|，其含义是对任意向量|v⟩左乘投影算子，必然得到该向量|v⟩在方向|uk⟩上的分量，这是针对向量空间的讨论，向量空间是线性空间最简单的形式。对于矩阵空间来说，该空间仍然是线性空间（因为满足线性空间的八条原则），因此亦可以选择若干基底矩阵来生成该空间。对于D×D矩阵而言，至多用D2个基底矩阵便可以描述该矩阵空间。特别地，在之前的内容中，我们已经证明一个矩阵可以写成其特征向量的展开的形式，因此一组|uk⟩⟨uk|总能够用来表示某一个特定的矩阵，即

Σ=D∑i=1λi|ui⟩⟨ui|

其中λi|ui⟩⟨ui|可认为是基底|ui⟩⟨ui|上的分量。

​ Dirac符号的一个简单用途就是速记代数形式对应的具体结构：⟨⋯⟩对应的是具体的数；⟨⋯|对应的是行向量；|⋯⟩对应的是列向量；|⋯|对应的是矩阵（方阵）。

3.3 题外话——矩阵微分

​ 现在来介绍矩阵微分的概念。给定向量a和b以及标量x、向量x和矩阵A与B，则有以下重要定义

∂标量∂标量略∂标量∂向量(∂x∂a)i=∂x∂ai∂标量∂矩阵∂x∂A=(∂x∂A11∂x∂A12⋯∂x∂A1n∂x∂A21∂x∂A22⋯∂x∂A2n⋮⋮⋮∂x∂An1∂x∂An2⋯∂x∂Ann)∂向量∂标量(∂a∂x)i=∂ai∂x∂向量∂向量(∂a∂b)ij=∂ai∂bj∂向量∂矩阵无∂矩阵∂标量∂A∂x=(∂A11∂x∂A12∂x⋯∂A1n∂x∂A21∂x∂A22∂x⋯∂A2n∂x⋮⋮⋮∂An1∂x∂An2∂x⋯∂Ann∂x)∂矩阵∂向量无∂矩阵∂矩阵无

容易得到下面的结论

(i)∂∂x(xTa)=∂∂x(aTx)=a(ii)∂∂x(AB)=∂A∂xB+A∂B∂x(iii)∂∂x(A−1)=−A−1∂A∂xA−1(iv)∂∂xln|A|=Tr(A−1∂A∂x)(v)∂∂AijTr(AB)=Bji(vi)∂∂ATr(AB)=BT(vii)∂∂ATr(ATB)=B(viii)∂∂ATr(A)=I(ix)∂∂ATr(ABAT)=A(B+BT)(x)∂∂Aln|A|=(A−1)T

3.4 高斯分布的矩

​ 在多元高斯分布中，我们从矩阵视角来分析矩。对于多元高斯分布

N(x|μμ,Σ)=1(2π)D/21|Σ|1/2exp{−12(x−μμ)TΣ−1(x−μμ)}

来说，它的一阶矩为

E(x)=∫N(x|μμ,Σ)xdx=1(2π)D/21|Σ|1/2∫exp{−12(x−μμ)TΣ−1(x−μμ)}xdx=1(2π)D/21|Σ|1/2∫exp{−12zTΣ−1z}(z+μμ)d(z+μμ)=1(2π)D/21|Σ|1/2∫(exp{−12zTΣ−1z}z+exp{−12zTΣ−1z}μμ)dz

由于指数部分是关于z的偶函数且整个积分区间为(−∞,+∞)，因此exp{−12zTΣ−1z}z部分关于z的积分为零，故E(x)=μμ，因此将其称为高斯分布的均值，这和我们的直观感觉是一致的。现在来看二阶矩。在一元情况下，二阶矩由E(x2)给出。对于多元变量而言，有D2个由E(xi⋅xj)给出的二阶矩，可以聚在一起写成矩阵形式，即多元情况下的二阶矩为

E(xxT)=1(2π)D/21|Σ|1/2∫exp{−12(x−μμ)TΣ−1(x−μμ)}xxTdx=1(2π)D/21|Σ|1/2∫exp{−12zTΣ−1z}(z+μμ)(z+μμ)Tdz

注意到(z+μμ)(z+μμ)T=(z+μμ)(zT+μμT)=zzT+zμμT+μμzT+μμμμT，其中涉及到zμμT和μμzT的积分值都将因为对称性（和一元情况类似）而等于零，涉及到常数μμμμT的积分值将等于μμμμT，下面讨论涉及到zzT的积分值。取单位正交基{|u1⟩,|u2⟩,⋯,|uD⟩}使得Σ−1=∑Di=11λi|ui⟩⟨ui|，且|z⟩=∑Di=1yi|ui⟩（其中yi=⟨ui|z⟩​），那么有对于指数部分有

zTΣ−1z=D∑i=1⟨ui|z⟩⟨ui|⋅D∑j=11λj|uj⟩⟨uj|⋅D∑k=1⟨uk|z⟩|uk⟩=D∑i=1D∑j=1D∑k=11λj⟨ui|z⟩⟨ui|uj⟩⟨uj|uk⟩⟨uk|z⟩=D∑i=1D∑j=1D∑k=11λj⟨ui|z⟩δijδjk⟨uk|z⟩=D∑i=11λi⟨ui|z⟩⟨ui|z⟩=D∑i=1y2iλi(=yTΛ−1y)

其中y=(y1,y2,⋯,yD)T，故z=Uy，其中U=(|u1⟩  |u2⟩  ⋯  |uD⟩)且|U|=1，于是有dz=dy，进而

1(2π)D/21|Σ|1/2∫exp{−12zTΣ−1z}zzTdz=1(2π)D/21|Σ|1/2∫exp{−12D∑k=1y2kλk}⋅(D∑i=1⟨ui|z⟩|ui⟩)⋅(D∑j=1⟨uj|z⟩⟨uj|)dy=1(2π)D/21|Σ|1/2∫exp{−12D∑k=1y2kλk}⋅D∑i=1D∑j=1⟨ui|z⟩⟨uj|z⟩|ui⟩⟨uj|dy=1(2π)D/21|Σ|1/2∫exp{−12D∑k=1y2kλk}⋅D∑i=1D∑j=1yiyj|ui⟩⟨uj|dy=1(2π)D/21|Σ|1/2D∑i=1D∑j=1|ui⟩⟨uj|∫exp{−12D∑k=1y2kλk}⋅yiyjdy=1(2π)D/21|Σ|1/2D∑i=1|ui⟩⟨ui|∫exp{−12D∑k=1y2kλk}y2idy=1(2π)D/21|Σ|1/2D∑i=1|ui⟩⟨ui|∬⋯∫exp{−12D∑k=1y2kλk}y2idy1dy2⋯dyD

其中倒数第三个等号到倒数第二个等号的原因是：当i≠j时，被积函数中必定有一项形如exp{−12y2i}yi，而这是奇函数，且积分区间是(−∞,+∞)，因此会等于零。注意，关于矢量的积分应该化成上述多重积分的形式。上述积分可以拆成

∬⋯∫exp{−y212λ1}⋅exp{−y222λ2}⋯exp{−y2D2λD}y2idy1dy2⋯dyD

并且易证下述两个结论

∫exp{−y2k2λk}dyi={0,k≠i(2πλi)1/2,   k=i∫exp{−y2k2λk}y2idyj={0,k≠i or i≠j(2πλi)1/2λi,k=i=j

将这两个结论代入拆解之后的式子，并进行整理，便可以求出

E(xxT)=涉及到μμμμT的积分值+涉及到zzT的积分值=μμμμT+Σ

​ 当然，对于多项式变量情形的二阶矩，我们也可以求二阶中心矩，即为

var(x)=E((x−E(x))(x−E(x)))=Σ

这即为协方差矩阵。

​ 一个协方差矩阵衡量了各变量之间的制约关系，如果协方差矩阵的自由参数越多、则描述的模型越复杂，如果协方差矩阵的自由参数越少、则描述的模型越简单。通常而言，一个任意的对称协方差矩阵有D(D+1)2各自由参数，加之D个自由的μ1,⋯,μD参数，该模型共有D(D+3)2个独立参数。如果协方差矩阵是对角矩阵，那么该模型共有2D个独立参数。如果进一步限制协方差矩阵为单位矩阵的倍数，即Σ=σ2I，则此时的协方差矩阵被称为是各向同性的（isotropic），此时模型共有D+1个独立参数。

3.5 条件高斯分布

​ 条件分布就是已知部分信息的情况下的概率分布，对于多项式变量x而言，将其分为两部分xa和xb，即(xaxb)，分别对应D个变量中的前M个变量后D−M个变量，那么此时的均值划分为μμ=(μμaμμb)，协方差矩阵划分为Σ=(ΣaaΣabΣbaΣbb)，注意Σba=ΣTab，并且Σaa和Σbb都是对称的。在之前我们已经多次使用Σ−1这个矩阵，现在对其命名如下

Λ≡Σ−1

这被称为精度矩阵（precision matrix），并且此时有Λ=(ΛaaΛabΛbaΛbb)，当然Λba=ΛTab并且Λaa和Λbb都是对称的。

​ 现在来寻找条件概率分布p(xa|xb)的表达式。对于二次型Δ2而言

Δ2=−12(x−μμ)TΣ−1(x−μμ)=−12(xa−μμa)TΛaa(xa−μμa)−12(xa−μμa)TΛab(xb−μμb)−12(xb−μμb)TΛba(xa−μμa)−12(xb−μμb)TΛbb(xb−μμb)

若将此二次型视为xa的函数，则这又是一个二次型，因此对应的分布p(xa|xb)亦是一个高斯分布。由于任意二次型可以重新写为

Δ2=−12xTΣ−1x+12xTΣ−1μμ+12μμTΣ−1x+对于x而言是常数=12xTΣ−1x+xTΣ−1μμ+对于x而言是常数

根据Dirac符号，⟨x|ΣΣ−1|μμ⟩和⟨μμ|ΣΣ−1|x⟩均为实数且两式互为转置关系，而实数的转置是它本身，因此两项可以合并。现在，只要将任意二次型整理成上述形式，就能在一次项中直接看出均值μμ、在二次项中直接看出精度矩阵Λ（即协方差矩阵Σ的逆矩阵），而条件概率分布的二次型可以写为

Δ2=−12xTaΛaaxa+xTaΛaa(μμa−Λ−1aaΛab(xb−μμb))+对于xa而言是常数

因此条件概率分布p(xa|xb)的均值和协方差矩阵分别为

μμa|b=μμa−Λ−1aaΛab(xb−μμb)Σa|b=Λ−1aa

​ 从上面的过程中可以看出，使用精度矩阵是十分方便的，对于分块矩阵的逆与各个分块的关系，有如下等式

(ABCD)−1=(M−MBD−1−D−1CMD−1+D−1CMBD−1)

其中M=(A−BD−1C)−1，并且称M−1为左侧矩阵关于子矩阵D的舒尔补（Schur complement）。用该等式可以求得条件概率分布p(xa|xb)的均值和协方差矩阵的另一种形式

μμa|b=μμa+ΣabΣ−1bb(xb−μμb)Σa|b=Σaa−ΣabΣ−1bbΣba

从上面可以看出，条件概率分布的均值仅与xb线性相关且协方差与xb无关，这是线性高斯（linear-Gaussian）模型的一个例子。

3.6 边缘高斯分布

​ 当给定联合分布p(xa,xb)的时候，依照上面的办法可以求出p(xb|xa)，现在来考虑边缘分布

p(xa)=∫p(xa,xb)dxb

显然边缘分布亦是高斯分布，现在的任务就是求出边缘高斯分布的均值和协方差矩阵。先提取出其中关于xb的项，对二次型Δ2进行整理得到

Δ2=−12xTbΛbbxb+xTbm+对于xb而言是常数=−12(xb−Λ−1bbm)TΛbb(xb−Λ−1bbm)+12mTΛ−1bbm+对于xb而言是常数

其中m=Λbbμμb−Λba(xa−μμa)。若取上式右侧第一项，就将被积函数整理成了关于xb的标准形式，此积分值的结果是归一化系数的倒数，从N(x|μμ,Σ)=1(2π)D/21|Σ|1/2exp{−12(x−μμ)TΣ−1(x−μμ)}中可知均值与协方差和归一化系数无关，因此不影响讨论。现在把上式右侧第二项与对于xb而言是常数的项合并得到

12mTΛ−1bbm−12xTaΛaaxa+xTa(Λaaμμa+Λabμμb)+对于xa和xb而言是常数=−12xTa(Λaa−ΛabΛ−1bbΛba)xa+xTa(Λaa−ΛabΛ−1bbΛba)μμa+对于xa和xb而言是常数

与高斯分布的标准形式相比，可以看到边缘分布p(xa)的均值和协方差分别由下面两式给出

Σa=(Λaa−ΛabΛ−1bbΛba)−1Σa(Λaa−ΛabΛ−1bbΛba)μμa=μμa

使用关于分块矩阵逆的恒等式有(Λaa−ΛabΛ−1bbΛba)−1=Σaa，因此

E(xa)=μμacov(xa)=Σaa

​ 在条件高斯分布中，用分块精度矩阵表示均值和协方差更加简便，但是对于边缘高斯分布而言，采用分块协方差矩阵更加简便。在条件高斯分布中，我们只需要分离出唯一变量xa，剩下的所有部分（包括xb）都可以认为是已知的、可以用在均值和协方差表达式中的；但是在边缘高斯分布中，两个变量xa和xb是平等的、都是未知的，因此变量xa的均值和协方差的表达式中不能出现xb，所以先分离出变量xb，再分离出xa就可以求解，注意，这里的变量分离顺序不能改变。

3.7 高斯变量的贝叶斯定理

​ 现在给出如下形式的边缘高斯分布和条件高斯分布

p(x)=N(x|μμ,Λ−1)p(y|x)=N(y|Ax+b,L−1)

其中，如果x的维度为M、y的维度为D，那么矩阵A的大小为D×M。对于联合自变量z=(xy)而言，联合概率分布的对数为

lnp(z)=lnp(x)+lnp(y|x)=−12(x−μμ)TΛ(x−μμ)−12(y−Ax−b)TL−1(y−Ax−b)+对于x和y而言是常数

这是z的分量的一个二次函数，因此亦是高斯分布，为了找到这个高斯分布的均值和协方差，需要将此二次型整理成二次项+一次项+常数的形式，先看二次项，得到

−12xT(Λ+ATLA)x−12yTLy+12yTLAx+12xTATLy=−12(xy)T(Λ+ATLA−ATL−LAL)(xy)=−12zTRz

因此高斯分布的精度矩阵（协方差矩阵的逆矩阵）为

R=(Λ+ATLA−ATL−LAL)

使用关于分块矩阵逆的恒等式可以得到协方差矩阵为

cov(z)=R−1=(Λ−1Λ−1ATAΛ−1L−1+AΛ−1AT)

​ 类似地，找到二次型中的一次项为

xTΛμμ−xTATLb+yTLb=(xy)T(Λμμ−ATLbLb)=zT⋅R⋅E(z)

于是均值为

E(z)=R−1(Λμμ−ATLbLb)=(μμAx+b)

​ 接下来来看边缘分布p(y)​的均值和协方差，使用已经得到的结论，可以推得

E(y)=Aμμ+bcov(y)=L−1+AΛ−1AT

​ 最后来寻找条件分布p(x|y)的表达式，在“条件高斯分布”一节中已经给出了μμa|b和Σa|b的一般表达式，现在直接带入可以得到

E(x|y)=(Λ+ATLA)−1{ATL(y−b)+Λμμ}cov(x|y)=(Λ+ATLA)−1

这给出了先验分布参数与后验分布参数的精确关系。

3.8 高斯分布的最大似然估计

​ 给定数据集X=(x1,x2,⋯,xN)T，且满足多元高斯分布，那么对数似然函数为

lnp(X|μμ,Σ)=−ND2ln(2π)−N2ln|Σ|−12N∑n=1(xn−μμ)TΣ−1(xn−μμ)

经过简单的重新排列，我们可以发现这个对数似然函数的充分统计量为∑Nn=1xn和∑Nn=1xnxTn。该对数似然函数对参数μμ的导数为

∂∂μμlnp(X|μμ,Σ)=N∑n=1Σ−1(xn−μμ)

令此导数为零，便得到参数μμ的极大似然估计为

μμML=1NN∑n=1xn

这和我们的直觉是相符的，并且有E(μμML)=μμ，即这个估计是无偏的。下面我们不加证明地指出参数ΣΣ的极大似然估计为

ΣML=1NN∑n=1(xn−μμML)(xn−μμML)T

这也和我们的直觉是相符的，但是E(ΣML)=N−1NΣ，即这个估计是有偏的，当给定数据集的时候，参数ΣΣ的一个无偏估计为

˜Σ=1N−1N∑n=1(xn−μμML)(xn−μμML)T

3.9 顺序估计

​ 对于上一节中得到的\pmb\mu_{ML}而言，如果我们想定量分析最后一个数据点的贡献时

\begin{aligned} \pmb\mu_{ML}^{(N)}&=\frac1N\sum_{n=1}^N\mathbf{x}_n\\ &=\frac1N\mathbf{x}_N+\frac1N\sum_{n=1}^{N-1}\mathbf{x}_n\\ &=\frac1N\mathbf{x}_N+\frac{N-1}{N}\pmb\mu_{ML}^{(N-1)}\\ &=\pmb\mu_{ML}^{(N-1)}+\frac1N(\mathbf{x}_N-\pmb\mu_{ML}^{(N-1)}) \end{aligned}

该式给出了每次得到一个新数据点之后修正参数\pmb\mu_{ML}的方法：将已经得到的参数\pmb\mu_{ML}^{(N-1)}沿着方向\mathbf{x}_N-\pmb\mu_{ML}^{(N-1)}移动一小段距离，且这段距离会随着数据集的不断扩大而减小。

​ 在上面的例子中，参数\pmb\mu的极大似然估计\pmb\mu_{ML}可以分离出最后一个数据点的贡献，但是在实际应用中，参数的极大似然估计的形式是十分多样的，不能确保一定能够从中分离出最后一个数据点的贡献。下面我们介绍一个更加普适的顺序学习方法：Robbins-Monro方法。一个参数的极大似然估计就是对应的负对数似然函数的一个驻点，即导数值等于零的解，现在先从纯数学角度进行分析：考虑随机变量\theta和z，它们的联合分布为p(z,\theta)，那么在已知\theta的情况下，z的条件期望定义了一个关于\theta的函数f(\theta)

f(\theta)\equiv E(z|\theta)=\int zp(z|\theta)dz

通过这种方式定义的函数被称为回归函数（regression function），它的含义是：每给定一个具体的\theta_0时，z的期望仅由\theta_0表示。我们的目标是寻找根\theta^*使得f(\theta^*)=0，Robbins-Monro方法给出了在顺序观测的情况下找到\theta^*的方法。假设z的条件方差是有穷的，因此

E((z-f)^2|\theta)<\infty

并且设当\theta<\theta^*时f(\theta)<0，当\theta>\theta^*时f(\theta)>0，Robbins-Monro方法通过定义下述序列给出了根\theta^*的估计为

\theta^{(N)}=\theta^{(N-1)}-\alpha_{N-1}z(\theta^{N-1})

其中z(\theta^{N-1})是当\theta取值为\theta^{N-1}时的观测值，系数\{\alpha_n\}表示满足下述三个条件的正数序列

\lim_{N\rightarrow\infty}\alpha_N=0,\quad \sum_{N=1}^\infty\alpha_N=\infty,\quad \sum_{N=1}^\infty\alpha_N^2<\infty

第一个条件保证了根的修正幅度会逐渐减小（因此能够收敛到一个有限值），第二个条件保证了不会收敛不到根的值（因此能够收敛到根的值），第三个条件保证了累计的噪声具有一个有限的方差（因此不会导致收敛失败）。

​ 对于任意一个负对数似然函数而言，它的参数\theta的极大似然估计满足

\frac{\partial}{\partial\theta}\{\frac1N\sum_{n=1}^N-\ln p(x_n|\theta)\}=0

交换导数与求和并取极限N\rightarrow\infty，得到

-\lim_{N\rightarrow\infty}\frac1N\sum_{n=1}^N\frac{\partial}{\partial\theta}\ln p(x_n|\theta)=E_x(-\frac{\partial}{\partial\theta}\ln p(x_n|\theta))

因此我们看到寻找极大似然估计对应于寻找回归函数的根。于是我们可以应用Robbins-Monro方法，此时它的形式为

\theta^{(N)}=\theta^{(N-1)}-\alpha_{N-1}\frac{\partial}{\partial\theta^{(N-1)}}[-\ln p(x_N|\theta^{(N-1)})]

​ 下面以高斯分布为例来看看该方法在实际情况中的应用。随机变量z为

z=-\frac{\partial}{\partial\mu_{ML}}\ln p(x|\mu_{ML},\sigma^2)=-\frac{1}{\sigma^2}(x-\mu_{ML})

因此z的分布仍是高斯分布，将此结果代入Robbins-Monro方法，得到

\mu_{ML}^{(N)}=\mu_{ML}^{(N-1)}-\frac{\sigma^2}{N-1}(-\frac{1}{\sigma^2}(x_N-\mu_{ML}^{(N-1)}))

其中，令\alpha_N=\frac{\sigma^2}{N}。该结果很容易推广到多元情形。

3.10 高斯分布的贝叶斯推断

​ 极大似然估计给出了均值和方差的点估计，现在引入这些参数的先验分布。首先对于一组随机变量\mathbf{x}而言，假设方差已知，需要推断均值，则有似然函数

p(\mathbf{x}|\mu)=\prod_{n=1}^Np(x_n|\mu)=\frac{1}{(2\pi\sigma^2)^{N/2}}\text{exp}\{-\frac{1}{2\sigma^2}\sum_{n=1}^N(x_n-\mu)^2\}

为了保证先验分布和后验分布的共轭性，现令先验分布为p(\mu)=N(\mu|\mu_0,\sigma_0^2)，从而后验概率满足

p(\mu|\mathbf{x})\varpropto p(\mathbf{x}|\mu)p(\mu)=N(\mu|\mu_N,\sigma_N^2)

其中

\mu_N=\frac{\sigma^2}{N\sigma_0^2+\sigma^2}\mu_0+\frac{N\sigma_0^2}{N\sigma_0^2+\sigma^2}\mu_{ML},\qquad\frac{1}{\sigma_N^2}=\frac{1}{\sigma^2}+\frac{N}{\sigma_0^2}

其中\mu_{ML}是\mu的极大似然估计，即为\mu_{ML}=\frac1N\sum_{n=1}^Nx_n。从上式自然可以看出：随着数据集规模N的增加，\mu_N会越来越接近\mu_{ML}。此外，由此式可以发现：精度（方差的倒数）是可以直接叠加的，随着数据集规模N的增加，精度会越来越大（即方差会越来越小）。该结论容易推广到多项式变量的情况。若从顺序角度来看高斯分布的贝叶斯推断，从后验分布中分离出最后一个数据点x_N得到

p(\mu|\mathbf{x})\varpropto[p(\mu)\prod_{n=1}^{N-1}p(x_N|\mu)]p(x_N|\mu)

该式明显揭示了每增加一个数据点所能带来的贡献。

​ 上面假设方差已知并推断均值，下面假设均值已知推断方差，保持似然函数和先验分布的共轭性，得到

\begin{aligned} p(\lambda|\mathbf{x})&\varpropto p(\mathbf{x}|\lambda)\cdot\text{Gam}(\lambda|a,b)\\ &\varpropto\lambda^{N/2}\text{exp}\{-\frac{\lambda}{2}\sum_{n=1}^N(x_n-\mu)^2\}\cdot\frac{1}{\Gamma(a)}b^a\lambda^{a-1}\text{exp}(-b\lambda)\\ &\equiv\text{Gam}(\lambda|a_N,b_N) \end{aligned}

其中\lambda\equiv\frac{1}{\sigma^2}，Gamma分布的形式为

\text{Gam}(\lambda|a,b)=\frac{1}{\Gamma(a)}b^a\lambda^{a-1}\text{exp}(-b\lambda)

并且可以得到

\begin{aligned} E(\lambda|a,b)&=\frac{a}{b}\\ \text{var}(\lambda|a,b)&=\frac{a}{b^2} \end{aligned}

上述参数a_N和b_N为

a_N=a_0+\frac{N}{2},\qquad b_N=b_0+\frac12\sum_{n=1}^N(x_n-\mu)^2=b_0+\frac{N}{2}\sigma_{ML}^2

其中\sigma_{ML}^2是方差的最大似然估计，即为\sigma_{ML}^2=\frac1N\sum_{n=1}^N(x_n-\mu)^2（此时假设\mu已知）。从a_N的表达式知道，当我们观测N个数据点的时候，使得参数a_N增加了\frac{N}{2}，因此可以将先验参数a_0解释为2a_0个先验的有效观测。在Beta分布中我们也将先验参数解释为有效观测，实际上，对于指数族分布而言，把共轭先验视为有效假想数据点是一个很通用的思路。

​ 如果在分析问题的时候考虑的是\sigma^2而非\lambda，那么引入的是逆Gamma（inverse gamma）分布，此处不介绍。

​ 在具体分析问题的时候，一个重要的技巧就是不需要始终注意一个概率的归一化系数，比如此处Gamma分布的归一化系数可以在任何时候与标准Gamma分布的形式相比较而得到。

​ 之前先后讨论了方差已知推断均值、均值已知推断方差两种情况，现在先来看均值和方差都未知时的情况。为了找到共轭先验，先考虑似然函数

p(\mathbf{x}|\mu,\lambda)=\prod_{n=1}^N(\frac{\lambda}{2\pi})^{1/2}\text{exp}\{-\frac{\lambda}{2}(x_n-\mu)^2\}\varpropto[\lambda^{1/2}\text{exp}(-\frac{\lambda\mu^2}{2})]^N\text{exp}\{\lambda\mu\sum_{n=1}^Nx_n-\frac{\lambda}{2}\sum_{n=1}^Nx_n^2\}

故假设先验分布为

p(\mu,\lambda)\varpropto[\lambda^{1/2}\text{exp}(-\frac{\lambda\mu^2}{2})]^\beta\text{exp}\{c\lambda\mu-d\lambda\}=\text{exp}\{-\frac{\beta\lambda}{2}(\mu-\frac{c}{\beta})\}\lambda^{\beta/2}\cdot\text{exp}\{-(d-\frac{c^2}{2\beta})\lambda\}

其中参数c，d和\beta都为常数，由于p(\mu,\lambda)=p(\mu|\lambda)\cdot p(\lambda)，将上式整理为

p(\mu,\lambda)=\mathcal N(\mu|\mu_0,(\beta\lambda)^{-1})\cdot\text{Gam}(\lambda|a,b)=p(\mu|\lambda)\cdot p(\lambda)

其中

\mu_0=\frac{c}{\beta},\quad a=\frac{1+\beta}{2},\quad b=d-\frac{c^2}{2\beta}

该先验分布称为被称为正态-Gamma（normal-gamma）分布或者高斯-Gamma（Gaussian-gamma）分布。

​ 对于D维多项式变量\mathbf{x}的高斯分布\mathcal N(\mathbf{x}|\pmb\mu,\mathbf\Lambda^{-1})，如果精度\mathbf\Lambda已知，那么均值\pmb\mu的先验分布仍为高斯分布；如果均值\pmb\mu已知，那么精度\mathbf\Lambda的先验分布定义如下

\mathcal{W}(\mathbf\Lambda|\mathbf W,\nu)=B|\mathbf\Lambda|^{(\nu-D-1)/2}\text{exp}(-\frac12\text{Tr}(\mathbf W^{-1}\mathbf\Lambda))

该分布称为Wishart分布，其中\nu被称为分布的自由度数量，\mathbf W是一个D\times D的标量矩阵，\text{Tr}(\cdot)表示矩阵的迹，归一化系数B(\mathbf W,\nu)为

B(\mathbf W,\nu)=|\mathbf W|^{-\nu/2}(2^{(\nu D)/2}\pi^{D(D-1)/4}\prod_{i=1}^D\Gamma(\frac{\nu+1-i}{2}))^{-1}

和之前一样，此处也可以用\pmb\Sigma而非\mathbf\Lambda参与讨论，那么引入的是逆Wishart（inverse Wishart）分布。

​ 如果均值\pmb\mu和\mathbf\Lambda都是未知的，那么共轭先验为

p(\pmb\mu,\mathbf\Lambda|\pmb\mu_0,\beta,\mathbf W,\nu)=\mathcal N(\pmb\mu|\pmb\mu_0,(\beta\mathbf\Lambda)^{-1})\cdot\mathcal{W}(\mathbf\Lambda|\mathbf W,\nu)

这被称为正态-Wishart（normal-Wishart）分布或者高斯-Wishart（Gaussian-Wishart）分布。

3.11 学生t分布

​ 高斯-Gamma用于均值和方差都未知时的参数估计，现在令\mathcal N(x|\mu,\tau^{-1})和\text{Gam}(\tau|a,b)，如果对\tau积分则得到

\begin{aligned} p(x|\mu,a,b)&=\int_0^\infty\mathcal N(x|\mu,\tau^{-1})\cdot\text{Gam}(\tau|a,b)d\tau\\ &=\int_0^\infty\frac{b^ae^{-b\tau}\tau^{a-1}}{\Gamma(a)}(\frac{\tau}{2\pi})^{1/2}\text{exp}\{-\frac\tau2(x-\nu)^2\}d\tau\\ &=\frac{b^a}{\Gamma(a)}(\frac{1}{2\pi})^{1/2}[b+\frac{(x-\mu)^2}{2}]^{-a-1/2}\Gamma(a+\frac12) \end{aligned}

现在定义新的参数\nu=2a和\lambda=\frac{a}{b}，则分布p(x|\mu,a,b)为

\text{St}(x|\mu,\lambda,\nu)=\frac{\Gamma(\frac\nu2+\frac12)}{\Gamma(\frac\nu2)}(\frac{\lambda}{\pi\nu})^{1/2}[1+\frac{\lambda(x-\mu)^2}{\nu}]^{-\nu/2-1/2}

该分布被称为学生t分布（Student's t-distribution）。参数\lambda可称为t分布的精度，即使它通常不等于方差的倒数。参数\nu称为t分布的自由度。当\nu=1时，t分布变为柯西分布（Cauchy distribution）；当\nu\rightarrow\infty时，t分布\text{St}(x|\mu,\lambda,\nu)变为高斯分布p(x|\mu,\lambda^{-1})，此时均值为\mu、精度为\lambda。

​ t分布可以视为无限多个同均值不同精度的高斯分布相加得到的，因此t分布比高斯分布具有更好的鲁棒性（robustness），即t分布更加集中于均值附近，因此更少受到少数离群点（outlier）的影响，下图红线表示使用t分布进行拟合，而绿线表示使用高斯分布进行拟合

​ 下面令\nu=2a，\lambda=\frac{a}{b}以及\eta=\frac{\tau b}{a}，从而给出t分布的另一种写法

\text{St}(x|\mu,\lambda,\nu)=\int_0^\infty\mathcal N(x|\mu,(\eta\lambda)^{-1})\cdot\text{Gam}(\eta|\frac\nu2,\frac\nu2)d\eta

之后便容易将此结果推广到D维多元高斯分布的情况并积分得到

\begin{aligned} \text{St}(\mathbf{x}|\pmb\mu,\mathbf\Lambda,\nu)&=\int_0^\infty\mathcal N(\mathbf{x}|\pmb\mu,(\eta\mathbf\Lambda)^{-1})\cdot\text{Gam}(\eta|\frac\nu2,\frac\nu2)d\eta\\ &=\frac{\Gamma(\frac{D}{2}+\frac\nu2)}{\Gamma(\frac\nu2)}\frac{|\mathbf\Lambda|^{1/2}}{(\pi\nu)^{D/2}}[1+\frac{\Delta^2}{\nu}]^{-D/2-\nu/2} \end{aligned}

其中\Delta^2是马氏距离的平方，\Delta^2=(\mathbf{x}-\pmb\mu)^T\mathbf\Lambda(\mathbf{x}-\pmb\mu)。且可以得到

\begin{aligned} E(\mathbf{x})&=\pmb\mu\\ \text{cov}(\mathbf{x})&=\frac{\nu}{\nu-2}\mathbf\Lambda^{-1}\\ \text{mode}(\mathbf{x})&=\pmb\mu \end{aligned}

其中\text{mode}(\mathbf{x})表示众数。一元情况的结论可以类似得到，此处略去。

3.12 周期变量

​ 周期变量常使用极坐标进行描述，这时用高斯分布不能很好地描述该变量。一个容易想到的处理周期变量的思路是在0\leq\theta\leq2\pi中选定一个方向作为原点，然后用传统的概率分布方法进行拟合，事实证明该思路在实际应用中局限颇大。对于周期变量的观测数据集D=\{\theta_1,\cdots,\theta_N\}，其均值\frac{\theta_1+\cdots+\theta_N}{N}严重依赖于坐标系的选择，为了找到均值的一个不变的度量。现将此观测视为单位圆上的点，即用模值为1的二维向量\mathbf{x}来描述变量\theta，对向量\mathbf{x}求平均得到\bar{\mathbf{x}}=\frac1N\sum_{n=1}^N\mathbf{x}_n，注意此时的||\bar{\mathbf{x}}||\leq1，即\bar{\mathbf{x}}通常位于单位圆的内部，对应的\bar{\theta}与极坐标原点的选择无关，这样将样本均值写为\bar{\mathbf{x}}=(\bar{r}\cos\bar\theta,\bar{r}\cos\bar\theta)，设\bar{\mathbf{x}}=(\bar{x}_1,\bar{x}_2)，则

\bar{x}_1=\bar{r}\cos\bar\theta=\frac1N\sum_{n=1}^N\cos\theta_n,\quad\bar{x}_2=\bar{r}\sin\bar\theta=\frac1N\sum_{n=1}^N\sin\theta_n

求两者的比值并代入反三角函数便可得到

\bar\theta=\arctan(\frac{\sum_n\sin\theta_n}{\sum_n\cos\theta_n})

​ 下面介绍高斯分布对周期变量的一种推广，该推广满足以下三个条件

p(\theta)\geq0,\quad\int_0^{2\pi}p(\theta)d\theta=1,\quad p(\theta+2\pi)=p(\theta)

因此对于变量\mathbf{x}=(x_1,x_2)而言，若均值为\pmb\mu=(\mu_1,\mu_2)，协方差矩阵为\mathbf\Sigma=\sigma^2\mathbf I_{2\times 2}，因此有

p(x_1,x_2)=\frac{1}{2\pi\sigma^2}\text{exp}\{-\frac{(x_1-\mu_1)^2+(x_2-\mu_2)^2}{2\sigma^2}\}

显然，如果此处的p(\mathbf{x})为常数，那么对应的轮廓线是一个圆。下面将此分布从笛卡尔坐标(x_1,x_2)转换到极坐标(r,\theta)得到的，即代入x_1=r\cos\theta和x_2=r\sin\theta，并且设\mu_1=r_0\cos\theta_0和\mu_2=r_0\sin\theta_0，从而得到

p(\theta|\theta_0,m)=\frac{1}{2\pi I_0(m)}\text{exp}\{m\cos(\theta-\theta_0)\}

其中m=\frac{r_0}{\sigma^2}。该分布称为von Mises分布或者环形正态（circular normal）分布，参数\theta_0对应分布的均值，参数m称为concentration参数（类似于高斯分布方差的倒数、即精度），参数I_0(m)称为零阶修正的第一类Bessel函数（zeroth-order Bessel function of the first kind），定义为

I_0(m)=\frac{1}{2\pi}\int_0^{2\pi}\text{exp}\{m\cos\theta\}d\theta

​ 现在考虑环形正态分布参数\theta_0和参数m的极大似然估计，对数似然函数为

\ln p(D|\theta_0,m)=-N\ln(2\pi)-N\ln I_0(m)+m\sum_{n=1}^N\cos(\theta_n-\theta_0)

可以解得参数\theta_0的极大似然估计为

\theta_0^{ML}=\arctan(\frac{\sum_n\sin\theta_n}{\sum_n\cos\theta_n})

现设A(m)=\frac{I_0'(m)}{I_0(m)}=\frac{I_1(m)}{I_0(m)}，则参数m的极大似然估计m_{ML}满足下面的式子

A(m_{ML})=(\frac1N\sum_{n=1}^N\cos\theta_n)\cos\theta_0^{ML}+(\frac1N\sum_{n=1}^N\sin\theta_n)\sin\theta_0^{ML}

​ 至于建立周期概率分布的通用方法，最简单的思路是使用观测的直方图：极坐标被划分成了固定大小的箱子，但是该思路具有较大的局限性。另一种方法类似于环形正态分布：先考察欧式空间的高斯分布，但是这会使得概率分布的形式异常复杂。最后⼀种方法的思想是，在实数轴上的任何合法的分布（例如高斯分布）都可以转化成周期分布，转化的方法是连续地把宽度为2\pi的区间映射为周期变量(0,2\pi)，这相当于把实数轴沿着单位圆进行缠绕，该方法最终求出的概率分布在计算上较为复杂。

​ 环形正态分布的⼀个局限性是这个分布是单峰的，但是通过将环形正态分布混合，我们可以得到一些应用性更广的模型。

3.13 混合高斯模型

​ 由于高斯分布只是单峰的，因此许多复杂的模型都不能仅只用朴素的高斯分布来描述，现在引入混合高斯（mixture of Gaussians）分布的概念，混合高斯分布是指高斯分布的线性叠加，即这样一个分布

p(\mathbf{x})=\sum_{k=1}^K\pi_k\mathcal N(\mathbf{x}|\pmb\mu_k,\mathbf\Sigma_k)

其中每个\mathcal N(\mathbf{x}|\pmb\mu_k\mathbf\Sigma_k)称为混合分布的一个成分（component），并且都有自己的均值\pmb\mu_k和协方差\mathbf\Sigma_k，特别地，\pi_k被称为混合系数（mixing coefficient），满足\sum_{k=1}^K\pi_k=1以及0\leq\pi_k\leq1。

​ 这种加权平均的思路很容易让我们想到全概率公式，也就是p(\mathbf{x})=\sum_{k=1}^Kp(k)p(\mathbf{x}|k)，将其中\pi_k视为选择第k个成分的先验概率p(k)，把密度\mathcal N(\mathbf{x}|\pmb\mu_k,\mathbf\Sigma_k)视为k条件下的条件密度，那么此时的p(\mathbf{x})在实际意义上就反映了后验概率p(k|\mathbf{x})，即正比于（但不等于）后验概率p(k|\mathbf{x})。后续章节将讲述该后验概率的重要应用，它也被称为责任（responsibility），根据贝叶斯定理，后验概率可以被表示为

\begin{aligned} \gamma_k(\mathbf{x})&\equiv p(k|\mathbf{x})\\ &=\frac{p(k)p(\mathbf{x}|k)}{\sum_lp(l)p(\mathbf{x}|l)}\\ &=\frac{\pi_k\mathcal N(\mathbf{x}|\pmb\mu_k\mathbf\Sigma_k)}{\sum_l\pi_l\mathcal N(\mathbf{x}|\pmb\mu_l,\mathbf\Sigma_l)} \end{aligned}

可以近似地认为，这就是将之前介绍的p(\mathbf{x})关于参数k做了归一化处理，使其反映随机变量k的概率分布而非随机变量\mathbf{x}的概率分布。

​ 现在，我们已经发现，混合高斯分布由以下三个参数控制

\begin{aligned} \pmb\pi&\equiv\{\pi_1,\pi_2,\cdots,\pi_K\}\\ \pmb\mu&\equiv\{\pmb\mu_1,\pmb\mu_2,\cdots,\pmb\mu_K\}\\ \mathbf\Sigma&\equiv\{\mathbf\Sigma_1,\mathbf\Sigma_2,\cdots,\mathbf\Sigma_K\} \end{aligned}

那么混合高斯分布的对数似然函数就是

\ln p(\mathbf{X}|\pmb\pi,\pmb\mu,\mathbf\Sigma)=\sum_{n=1}^N\ln\{\sum_{k=1}^K\pi_kN(\mathbf{x}_n|\pmb\mu_n,\mathbf\Sigma_n)\}

其中\mathbf{X}=\{\mathbf{x}_1,\mathbf{x}_2,\cdots,\mathbf{x}_N\}。这样就可以使用极大似然法确定各个参数的极大似然估计了，具体内容将在后面的章节讨论。

4 指数族分布

​ 本章到目前为止我们正式接触到的分布均为指数族分布（exponential family）的具体的例子，参数\pmb\eta关于变量\mathbf{x}的一般情形的指数族分布的形式为

p(\mathbf{x}|\pmb\eta)=h(\mathbf{x})g(\pmb\eta)\text{exp}\{\pmb\eta^Tu(\mathbf{x})\}

其中\mathbf{x}可能是标量、也可能是向量，可能是离散的、也可能是连续的。函数g(\pmb\eta)充当了归一化系数的作用。下面来讨论已经遇到的三种主要分布的指数族分布标准形式。

4.1 伯努利分布的指数族分布形式

​ 有如下推导

\begin{aligned} p(x|\mu)&=\text{Bern}(x|\mu)=\mu^x(1-\mu)^{1-x}\\ &=\text{exp}\{x\ln\mu+(1-x)\ln(1-\mu)\}\\ &=(1-\mu)\text{exp}\{\ln(\frac{\mu}{1-\mu})\cdot x\} \end{aligned}

容易得到\eta=\ln(\frac{\mu}{1-\mu})，从中解出\mu=\sigma(\eta)=\frac{1}{1+\text{exp}(-\eta)}，\sigma(\eta)被称为logistic sigmoid函数，那么就有

p(x|\mu)=\sigma(-\eta)\text{exp}(\eta x)\\ u(x)=x,\quad h(x)=1,\quad\eta=\ln(\frac{\mu}{1-\mu}),\quad g(\eta)=\sigma(-\eta)

4.2 多项式分布的指数族形式

​ 有如下推导

\begin{aligned} p(\mathbf{x}|\pmb\mu)&=\text{Mult}(\mathbf{x}|\pmb\mu)=\prod_{k=1}^M\mu_k^{x_k}=\text{exp}\{\sum_{k=1}^Mx_k\ln\mu_k\}\\ &=\text{exp}(\pmb\eta^T\mathbf{x}) \end{aligned}

容易得到

u(\mathbf{x})=\mathbf{x},\quad h(\mathbf{x})=1,\quad\pmb\eta=(\ln\mu_1,\cdots,\ln\mu_M)^T,\quad g(\pmb\eta)=1

​ 由于参数\mu_k受到\sum_{k=1}^M\mu_k=1的限制，因此我们可以用前M-1个参数去表示出最后一个参数\mu_M，那么就有

\begin{aligned} p(\mathbf{x}|\pmb\mu)&=\text{Mult}(\mathbf{x}|\pmb\mu)=\prod_{k=1}^M\mu_k^{x_k}=\text{exp}\{\sum_{k=1}^Mx_k\ln\mu_k\}\\ &=\text{exp}\{\sum_{k=1}^{M-1}x_k\ln\mu_k+(1-\sum_{k-1}^{M-1}x_k)\ln{(1-\sum_{k-1}^{M-1}x_k)}\}\\ &=\text{exp}\{\sum_{k=1}^{M-1}x_k\ln(\frac{\mu_k}{1-\sum_{j=1}^{M-1}\mu_j})+\ln(1-\sum_{k-1}^{M-1}x_k)\} \end{aligned}

可以令\eta_k=\ln(\frac{\mu_k}{1-\sum_{j=1}^{M-1}\mu_j})，在此式两侧对k求和，反解得\mu_k=\frac{\text{exp}(\eta_k)}{1+\sum_{j}\text{exp}(\eta_j)}，这被称为这被称为softmax函数或者归一化指数（normalized exponential），那么就有

p(\mathbf{x}|\pmb\eta)=(1+\sum_{k=1}^{M-1}\text{exp}(\eta_k))^{-1}\text{exp}(\pmb\eta^T\mathbf{x})\\ u(\mathbf{x})=\mathbf{x},\quad h(\mathbf{x})=1\\ \pmb\eta=(\ln(\frac{\mu_1}{1-\sum_{j=1}^{M-1}\mu_j}),\cdots,\ln(\frac{\mu_{M-1}}{1-\sum_{j=1}^{M-1}\mu_j}),0)^T\\ g(\pmb\eta)=(1+\sum_{k=1}^{M-1}\text{exp}(\eta_k))^{-1}

4.3 高斯分布的指数族分布形式

​ 有如下推导

\begin{aligned} p(x|\mu,\sigma^2)&=\mathcal N(x|\mu,\sigma^2)=\frac{1}{(2\pi\sigma^2)^{1/2}}\text{exp}\{-\frac{1}{2\sigma^2}(x-\mu)^2\}\\ &=\frac{1}{(2\pi\sigma^2)^{1/2}}\text{exp}\{-\frac{1}{2\sigma^2}x^2+\frac{\mu}{\sigma^2}x-\frac{1}{2\sigma^2}\mu^2\} \end{aligned}

经过推导可以得到

u(x)=\left(\begin{array}{c}x\\x^2\end{array}\right),\quad h(x)=(2\pi)^{-1/2},\quad\pmb\eta=\left(\begin{array}{c}{\mu}/{\sigma^2}\\{-1}/{2\sigma^2}\end{array}\right),\quad g(\pmb\eta)=(-2\eta_2)^{1/2}\text{exp}(\frac{\eta_1^2}{4\eta_2})

4.4 最大似然估计

​ 下面来看参数\pmb\eta的最大似然估计的一般形式，对指数族分布的一般形式p(\mathbf{x}|\pmb\eta)=h(\mathbf{x})g(\pmb\eta)\text{exp}\{\pmb\eta^Tu(\mathbf{x})\}而言，其似然函数为

p(\mathbf{X}|\pmb\eta)=\prod_{n=1}^Np(\mathbf{x}_n|\pmb\eta)=\prod_{i=1}^Nh(\mathbf{x}_i)\cdot(g(\pmb\eta))^N\cdot\text{exp}\{\pmb\eta^T(\sum_{j=1}^Nu(\mathbf{x}_j))\}

两侧对\pmb\eta取偏导得到

N(g(\pmb\eta))^{N-1}\cdot\nabla g(\pmb\eta)\cdot\text{exp}\{\pmb\eta^T(\sum_{j=1}^Nu(\mathbf{x}_j))\}+(g(\pmb\eta))^N\cdot\text{exp}\{\pmb\eta^T(\sum_{j=1}^Nu(\mathbf{x}_j))\}\cdot(\sum_{k=1}^Nu(\mathbf{x}_k))=0

整理得到

-\frac{\nabla g(\pmb\eta)}{g(\pmb\eta)}=\frac{u(\mathbf{x}_1)+\cdots+u(\mathbf{x}_N)}{N}

也就是

-\nabla\ln g(\pmb\eta_{ML})=\frac1N\sum_{n=1}^Nu(\mathbf{x}_n)=E(u(\mathbf{x})),\quad\text{在等概率情况下最后一个等号成立}

由于参数\pmb\eta极大似然估计仅由\sum_{n=1}^Nu(\mathbf{x}_n)产生，因此这个量就是指数族分布的充分统计量。在具体应用的时候，我们不需要存储整个数据集，只需要存储充分统计量的值即可。例如，在伯努利分布中，u(x)=x，那么充分统计量为\sum_{n=1}^Nx_n，所以我们只需要存储数据点\{x_n\}的和即可。类似地，在高斯分布中，u(x)=(x,x^2)^T，因此我们只需要存储\{x_n\}的和以及\{x_n^2\}的和即可。

4.5 共轭先验

​ 在伯努利分布与二项分布中，共轭先验是Beta分布；在范畴分布与多项式分布中，共轭先验是Dirichlet分布；在高斯分布中，对于不同情况而言，共轭先验分别为高斯分布、Gamma分布、高斯-Gamma分布、Wishart分布或高斯-Wishart分布。实际上，对于指数族分布而言，有如下统一形式的共轭先验

p(\pmb\eta|\chi,\nu)=f(\chi,\nu)g(\pmb\eta)^\nu\text{exp}\{\nu\pmb\eta^T\chi\}

其中f(\chi,\nu)是归一化系数，g(\pmb\eta)充当了归一化系数的作用。将此先验分布与似然函数p(\mathbf{X}|\pmb\eta)相乘得到

p(\pmb\eta|\mathbf{X},\chi,\nu)\varpropto g(\pmb\eta)^{\nu+N}\text{exp}\{\pmb\eta^T(\sum_{n=1}^Nu(\mathbf{x}_n)+\nu\chi)\}

容易发现这满足共轭性。参数\nu可以看成是先验分布中假想观测的有效观测数，在给定\chi的情况下，每个假想观测都对充分统计量有贡献\nu\chi。

4.6 无信息先验

​ 到目前为止，我们得到先验分布的思路都是从结构上根据共轭性而来的，并没有说明先验分布中各个参数具体取值的选择方法。从之前的内容可以知道，先验分布中各个参数的取值会对后验分布产生一定的影响，现在考虑是否可以选择合适形式的先验分布从而尽可能减小先验分布中参数值对后验分布的影响，这种形式的先验分布称为被称为无信息先验（noninformative prior）。

​ 现在假设我们有一个由参数\lambda控制的分布p(x|\lambda)，对于先验分布p(\lambda)，一个最朴素的想法就是取为常数K，如果\lambda的取值为有限个N个离散变量，那么p(\lambda)=K是合理的，这相当于NK=1；如果\lambda的取值为某个有限连续区间(a,b)上的连续值，那么p(\lambda)=K仍然是合理的，这相当于(b-a)K=1；但是，如果\lambda的取值是无限个离散值或者无限区间上的连续值，p(\lambda)=K就是不合理的，此时先验分布无法被正确地归一化，因为对\lambda的积分是发散的，这样的先验分布被称作是反常的（improper）。并且，将p(\lambda)取为常数K可能会在变量替换的时候出现问题，比如令\lambda=\eta^2，则p_\eta(\eta)=p_\lambda(\lambda)\cdot|\frac{d\lambda}{d\eta}|=p_\lambda(\eta^2)\cdot2\eta\varpropto\eta\neq K。

​ 下面来看无信息先验的两个简单例子。第一，如果概率密度的形式为p(x|\mu)=f(x-\mu)，参数\mu被称为位置参数（location parameter）。如果这一类概率具有平移不变性（translation invariance），则说明我们选择的先验分布p(\mu)必须对区间(a,b)和区间(a-c,b-c)赋予相同的概率密度，即

\int_a^bp(\mu)d\mu=\int_{a-c}^{b-c}p(\mu)d\mu=\int_a^bp(\mu-c)d\mu

并且这对任意选择的a与b都成立，因此p(\mu)=p(\mu-c)，故p(\mu)的取值为常数，这是平移不变性的重要性质。以高斯分布为例，均值\mu是一种位置参数，由“高斯分布的贝叶斯推断”一节中的表达式

\mu_N=\frac{\sigma^2}{N\sigma_0^2+\sigma^2}\mu_0+\frac{N\sigma_0^2}{N\sigma_0^2+\sigma^2}\mu_{ML}

可知，随着数据集规模N的不断增加，先验参数\mu_0对后验参数\mu_N的影响不断减小。

​ 第二，如果概率密度的形式为p(x|\sigma)=\frac1\sigma f(\frac{x}{\sigma})​，那么参数\sigma​被称为缩放参数（scale parameter）。如果这一类概率具有缩放不变性（scale invariance），则说明我们选择的先验分布p(\sigma)必须对区间(a,b)和区间(\frac{a}{c},\frac{b}{c})​赋予相同的概率密度，即

\int_a^bp(\sigma)d\sigma=\int_{\frac{a}{c}}^{\frac{b}{c}}p(\sigma)d\sigma=\int_a^bp(\frac1c\sigma)\frac1cd\sigma

并且这对任意选择的a与b都成立，因此p(\sigma)=p(\frac1c\sigma)\frac1c，故p(\sigma)\varpropto\frac1\sigma，这是一个反常先验分布，因为对于0\leq\sigma\leq\infty上的积分是发散的，但p_{\ln\sigma}(\ln\sigma)=p_\sigma(\ln\sigma)\cdot|\frac{d\ln\sigma}{d\sigma}|\varpropto\frac1\sigma\cdot\sigma=1，即p_{\ln\sigma}(\ln\sigma)的值为常数，这是缩放不变性的重要性质。以高斯分布为例，方差\sigma是一种缩放参数，由“高斯分布的贝叶斯推断”一节中的表达式

\frac{1}{\sigma_N^2}=\frac{1}{\sigma^2}+\frac{N}{\sigma_0^2}

可知，随着数据集规模N的不断增加，先验参数\sigma_0对后验参数\sigma_N的影响不断减小。

5 非参数化方法

​ 我们已经详细地从贝叶斯主义的角度介绍了一些概率分布的具体形式，它们通常由一些参数控制，且这些参数通常由数据集推得。但是在实际应用中，这种参数化（parametric）方法有时具有较强的局限性。现在介绍如何从频率主义的角度进行建模。

5.1 直方图方法

​ 直方图方法的核心是分割出数个区域，并分别分析每个区域具有的特征。对于离散数据点，该方法很好理解，即将每个数据点放入对应的区域即可。对于连续的情况，我们将连续变量分割称若干个区域\Delta_i，将数据点放入对应的区域，得到每个区域的概率为

p_i=\frac{n_i}{N\Delta_i}

这就完成了建模，特别地，常令每个区域的覆盖范围的大小相同，即令\Delta_1=\Delta_2=\cdots\equiv\Delta。对于维数较低的情形（比如一元变量x），该方法是有效的；但对于维数较高的情形，该方法会导致维数灾难（小区域的数量为M^D，每一维都分成M个区间）。在直方图方法中，\Delta的选择是十分关键的，这将直接影响到建模的优劣程度，比如

上图中绿线表示真实的分布，蓝色方框表示使用直方图方法建立的模型。另外，在原则上，区域边界的选择也会影响模型的有效性，但影响程度通常小于\Delta的选取。

​ 直方图方法的核心在于用离散形式取考虑连续形式的性质，通过考虑某一点及其邻域的数据点从而得到该区域的特征，并且为了保证各区域之间的连贯性（即尽可能降低离散性的影响），还需要考虑局部区域的空间扩展（即这里的区域大小\Delta）。一方面，较大的\Delta会尽可能保证该区域内性质的准确性，另一方面，较小的\Delta会尽可能保证各个区域之间的连贯性（因为当若干数据点被放入同一个区域时，它们之间的差异性被忽略了，这就导致各个区域之间数据点的差异性被放大了，从而导致各个区域之间的连贯性变差）。接下来将介绍两个常用的密度估计的非参数化方法，这两种方法与直方图方法的核心理念是相通的，但对于维数的放大具有更好的适应性。

5.2 核密度方法

​ 先来讨论核密度方法与近邻方法的公共前提。对于D维欧氏空间而言，一个小区域的概率质量为

P=\int_Rp(\mathbf{x})d\mathbf{x}

现在假设数据集中有服从p(\mathbf{x})的N个数据点，那么位于区域R内部的数据点的数量K满足二项分布

\text{Bin}(K|N,P)=\frac{N!}{K!(N-K)!}P^K(1-P)^{N-K}

当区域R的体积V较小（从而p(\mathbf{x})在R内近似为常数）、且数据点的数量N​较大时，我们有区域R上的概率密度的估计式为

p(\mathbf{x})=\frac{K}{NV}

接下来有两个思路进行下面的处理：一是核密度方法，即固定V然后从数据中确定K；二是近邻方法，即固定K然后从数据中确定V的值。可以证明，在N\rightarrow\infty的情况下，如果V随着N而合适地收缩，并且K随着N而增大，那么两种方法得到的概率密度估计值都会收敛到真实的概率密度。

​ 先来看核密度方法。如果我们取区域R以\mathbf{x}​为中心且边长为1的小超立方体，为了统计落在该小超立方体中的数据点的数量，则可以设

k(\mathbf{u})= \begin{cases} 1,\quad|u_i|\leq\frac12,\quad i=1,\cdots,D\\ 0,\quad\text{others} \end{cases}

这是核函数（kernel function）的一个例子，也被称为Parzen窗。对于规模为N的数据集而言，位于小区域内的数据点的总数为

K=\sum_{n=1}^Nk(\frac{\mathbf{x}_n-\mathbf{x}}{h})

由此式给出的估计称为核密度估计（kernel density estimator）或者Parzen估计。对于以\mathbf{x}为中心且边长为1的小超立方体区域、且选择上述核函数而言，其中h=1表示小超立方体的边长，则\mathbf{x}点处的核密度估计为

p(\mathbf{x})=\frac1N\sum_{n=1}^N\frac{1}{h^D}k(\frac{\mathbf{x}_n-\mathbf{x}}{h})

​ 但是核密度方法也具有直方图方法的一个重要缺陷：人为造成的区域间的连贯性变差，如果选择一个较为平滑的核函数（即：距离\mathbf{x}不同远近的数据点被放入区域内的条件不同），那么就能得到一个较好的模型。一个常见的选择就是高斯核函数，从而得到

p(\mathbf{x})=\frac1N\sum_{n=1}^N\frac{1}{(2\pi h^2)^{D/2}}\text{exp}\{-\frac{||\mathbf{x}_n-\mathbf{x}||^2}{2h^2}\}

其中h表示高斯分布的标准差，如果h过小则会模型对噪声过于敏感，如果h过大则会造成模型过于平滑。

​ 事实上，核函数的选取是任意的，只需要满足下面两个条件

k(\mathbf{u})\geq0,\qquad\int k(\mathbf{u})d\mathbf{u}=1

5.3 近邻方法

​ 核密度方法的一个局限是对于所有区域的数据都给定了同样的核参数h，在数据点密集的区域，这可能导致模型过于平滑，而在数据点稀疏的区域，这又可能导致模型对噪声过于敏感，所以一个合理的思路就是根据数据空间的位置确定不同的核参数h，这通过近邻方法来具体处理。

​ 对于概率密度估计式

p(\mathbf{x})=\frac{K}{NV}

而言，现在考虑固定K然后从数据中确定V的值。考虑以\mathbf{x}为中心的小超球体，其半径可以自由变化直到该区域内包含了K个数据点，但是这样得到的模型并不是真实的概率密度模型，因为它在整个空间的积分是发散的。

​ 近邻方法可以推广到分类问题，对于每个具体的分类C_k而言，有

p(\mathbf{x}|C_k)=\frac{K_k}{N_kV}

并且无条件概率密度为

p(\mathbf{x})=\frac{K}{NV}

并且类先验概率密度为

p(C_k)=\frac{N_k}{N}

则根据贝叶斯定理得到后验概率密度为

p(C_k|\mathbf{x})=\frac{p(\mathbf{x}|C_k)p(C_k)}{p(\mathbf{x})}=\frac{K_k}{K}

现在，如果我们想最小化错误分类的概率，则将测试点\mathbf{x}分配给具有最大后验概率密度的类别，这对应于最大的\frac{K_k}{K}。当K=1时，分类规则被称为最近邻规则（nearest-neighbour rule），因为测试点简单地被分类为训练数据集中距离最近的数据点的类别，下图是一个例子。

图（a）中K=3，此时新输入一个测试点（图中黑色菱形）时，找到与之最近的三个数据点，发现红色数据点多，因此测试点为红色数据点。图（b）中K=1，此时新输入一个测试点时，找到与之最近的两个数据点，观察其距离哪个数据点更近则说明它属于哪一个类别，图中所有类间点对的平分线生成上图中的绿线，在实际应用中，只需要观察它在绿线的哪一侧即可确定其是红色数据点还是蓝色数据点。

​ 核密度方法和近邻方法都需要存储整个训练集。如果训练集很大则会造成很大的计算代价。我们可以建立⼀个基于树的搜索结构，使得（近似）近邻可以高效地被找到，而不必遍历整个数据集。尽管这样，这些非参数化方法仍然有较大的局限性。另外，简单的参数化模型非常受限，因为它们只能表示某一种形式的概率分布。因此，我们 需要寻找⼀种概率密度模型，这种模型既要有灵活性、又要保证它的复杂度可以被控制为与训练数据的规模无关。我们在后续章节中将会看到如何找到这种概率密度模型。

6 参考资料

• Christopher M. Bishop, Pattern Recognition and Machine Learning, Springer, 2006
• Markus Svensen, Christopher M. Bishop, Pattern Recognition and Machine Learning - Solutions to the Exercises: Tutors’ Edition, Springer, 2009
• 马春鹏，《模式识别与机器学习》（本文部分名词翻译来自此书），PRML的网传中文版，2014
• Multinoulli分布与多项式分布
• Multinomial Distribution
• Dirichlet Distribution
• Why: "we note that the matrix Σ can be taken to be symmetric, without loss of generality"?
• Mikio Nakahara, Tetsuo Ohmi, Quantum Computing - From Linear Algebra to Physical Realizations, CRC press, 2008
• 向量或矩阵的微分计算
• Matrix Spaces; Rank 1; Small World Graphs
• Expressing a matrix as an expansion of its eigenvalues
• Expected value of xx^{T} for multidimensional Gaussian
• PRML 模式识别和机器学习 从零开始的公式推导 2.3高斯分布
• PRML 模式识别和机器学习 从零开始的公式推导 2.3.1条件高斯分布
• Kernel Density Estimation