第3章 概率密度函数的估计 3.1引言 在上一章最后的讨论中已经提到,贝叶斯决策的基础是概率密度函数的估计,即根据一定的训练样本来估计统计决策中用到的先验概率P(ωi)和类条件概率密度p(x|ωi)。其中,先验概率的估计比较简单,通常只需根据大量样本计算出各类样本在其中所占的比例,或者根据对所研究问题的领域知识事先确定。因此,本章重点介绍类条件概率密度的估计问题。 这种首先通过训练样本估计概率密度函数,再用统计决策进行类别判定的方法称作基于样本的两步贝叶斯决策。 这样得到的分类器性能与第2章理论上的贝叶斯分类器有所不同。我们希望当样本数目N→∞时,基于样本的分类器能收敛于理论上的结果。为做到这一点,实际只要说明N→∞时,估计的p^(x|ωi)和P^(ωi)收敛于p(x|ωi)和P(ωi)。这在统计学中可通过对估计量性质的讨论来解决。 在监督学习中,训练样本的类别是已知的,而且假定各类样本只包含本类的信息,这在多数情况下是正确的。因此,我们要做的是利用同一类的样本来估计本类的类条件概率密度。为了讨论方便,在本章后面的论述中,除了特别说明外,一概假定所有样本都是来自同一类,不再标出类别标号。由于概率密度估计是统计推断中的一个重要内容,有很多专门的教科书介绍它,所以在本书中只扼要地介绍其中的基本思想和主要方法,并不深入地展开论述。 概率密度函数的估计方法分为两大类: 参数估计参数估计(parametric estimation)与非参数估计非参数估计(nonparametric estimation)。 参数估计中,已知概率密度函数的形式,但其中部分或者全部参数未知,概率密度函数的估计问题就是用样本来估计这些参数。主要方法又有两类: 最大似然估计和贝叶斯估计贝叶斯估计,两者在很多实际情况下结果接近,但从概念上它们的处理方法是不同的。 3.1引言 非参数估计,就是概率密度函数的形式也未知,或者概率密度函数不符合目前研究的任何分布模型,因此不能仅仅估计几个参数,而是用样本把概率密度函数数值化地估计出来。 参数估计是统计推断的基本问题之一,在讨论具体问题之前先介绍几个参数估计中的基本概念。 (1) 统计量统计量。样本中包含着总体的信息,希望通过样本集把有关信息抽取出来,就是说针对不同要求构造出样本的某种函数,这种函数在统计学中称为统计量。 (2) 参数空间参数空间。如上所述,在参数估计中,总是假设总体概率密度函数的形式已知,而未知的仅是分布中的几个参数,将未知参数记为θ,在统计学中,将总体分布未知参数θ的全部可容许值组成的集合称为参数空间,记为Θ。 (3) 点估计点估计、估计量和估计值。点估计问题就是要构造一个统计量d(x1,…,xN)作为参数θ的估计θ^,在统计学中称θ^为θ的估计量。如果x(i)1,…,x(i)N是属于类别ωi的几个样本观察值,代入统计量d就得到对于第i类的θ^的具体数值,这个数值在统计学中称为θ的估计值。 (4) 区间估计区间估计。除点估计外,还有另一类估计,它要求用区间(d1,d2)作为θ可能取值范围的一种估计。这个区间称为置信区间置信区间,这类估计问题称为区间估计。 本章要求估计总体分布的具体参数,显然这是点估计问题。我们将介绍两种主要的点估计方法——最大似然估计最大似然估计和贝叶斯估计,它们都能得到相应的估计值。当然评价一个估计的“好坏”,不能按一次抽样结果得到的估计值与参数真值θ的偏差大小来确定,而必须从平均的和方差的角度出发进行分析。为了表示这种偏差,统计学中做了很多关于估计量性质的定义。 在数理统计中,用来判断估计好坏的常用标准是无偏性无偏性、有效性有效性和一致性一致性,这是本章介绍的估计方法的基本出发点。如果参数θ的估计量θ^(x1,x2,…,xn)的数学期望等于θ,则称估计是无偏的; 如果当样本数趋于无穷时估计才具有无偏性,则称为渐近无偏。如果一种估计的方差比另一种估计的方差小,则称方差小的估计更有效。而如果对于任意给定的正数ε,总有 limn→∞P(|θ^n-θ|>ε)=0 则称θ^是θ的一致估计。显然,无偏性、有效性都只是说明对于多次估计来说,估计量能以较小的方差平均地表示其真实值,并不能保证具体的一次估计的性能; 而一致性则保证当样本数无穷多时,每一次的估计量都将在概率意义上任意地接近其真实值。 在统计学中常见的一些典型分布形式并不总是能够拟合所有实际数据的分布,此外,在很多实际问题中还会遇到多峰分布的情况,在另外一些情况下我们可能无法事先判断数据的分布情况。在这些情况下,都需要直接利用样本去非参数地估计概率密度。本章将介绍其中最基本的直方图法、kN近邻法和Parzen窗法。 3.2最大似然估计最大似然估计 3.2.1最大似然估计的基本原理 3.2最大似然估计 在最大似然估计(maximum likelihood estimation)中,我们做以下基本假设: (1) 我们把要估计的参数记作θ,它是确定但未知的量(多个参数时是向量),这与把它看作随机量的方法是不同的。 (2) 每类的样本集记作Xi,i=1,2,…,c,其中的样本都是从密度为p(x|ωi)的总体中独立抽取出来的,即所谓满足独立同分布条件。 (3) 类条件概率密度p(x|ωi)具有某种确定的函数形式,只是其中的参数θ未知。例如在x是一维正态分布N(μ,σ2)时,未知的参数就可能是θ=[μ,σ2]T,对不同类别的参数可以记作θi,为了强调概率密度中待估计的参数,也可以把p(x|ωi)写作p(x|ωi,θi)或p(x|θi)。 (4) 各类样本只包含本类的分布信息,也就是说,不同类别的参数是独立的,这样就可以分别对每一类单独处理。 在这些假设的前提下,我们就可以分别处理c个独立的问题。即,在一类中独立地按照概率密度p(x|θ)抽取样本集X,用X来估计出未知参数θ。 设样本集包含N个样本,即 X={x1,x2,…,xN}(31) 由于样本是独立地从p(x|θ)中抽取的,所以在概率密度为p(x|θ)时获得样本集X的概率即出现X中的各个样本的联合概率是 l(θ)=p(X|θ)=p(x1,x2,…,xN|θ)=∏Ni=1p(xi|θ)(32) 这个概率反映了在概率密度函数的参数是θ时,得到式(31)中这组样本的概率。现在因为已经得到了式(31)的样本集,而θ是不知道的,式(32)就成为θ的函数,它反映的是在不同参数取值下取得当前样本集的可能性,因此称作参数θ相对于样本集X的似然函数似然函数(likelihood function)。式(32)右边乘积中的每一项p(xi|θ)就是θ相对于每一个样本的似然函数。 似然函数l(θ)给出了从总体中抽出x1,x2,…,xN这样N个样本的概率。为了便于解释,暂且假定θ是已知的,用θ0表示这个已知值。最可能出现的N个样本是使得l(θ)值为最大的样本x′1,x′2,…,x′N。例如,为了简便,假定N=1,x为一维且具有以均值为6,方差为1的正态分布,那么“最可能出现的”样本就是x′1=6,这时似然函数l(θ)=p(x′1|6,1)=maxx∈Ed p(x|6,1)。再回过来看,假定θ为未知,而我们从一次抽样中得到N个样本x1,x2,…,xN,想知道这组样本“最可能”来自哪个密度函数。换句话说,所抽取出的这组样本,来自哪个密度函数(θ取什么值)的可能性最大?即我们要在参数空间Θ中找到一个θ值(用θ^表示),它能使似然函数l(θ)极大化。一般来说,使似然函数的值最大的θ^是样本x1,x2,…,xN的函数,记为θ^=d(x1,x2,…,xN),就把θ^=d(x1,x2,…,xN)叫做θ的最大似然估计量。图31示意了最大似然估计的基本原理。 图31最大似然估计的基本原理 现在可以给出最大似然估计量的定义。 最大似然估计量: 令l(θ)为样本集X的似然函数,X={x1,x2,…,xN},如果θ^=d(X)=d(x1,x2,…,xN)是参数空间Θ中能使似然函数l(θ)极大化的θ值,那么θ^就是θ的最大似然估计量,或者记作 θ^=arg max l(θ)(33) 其中,arg max是一种常用的表示方法,表示使后面的函数取得最大值的变量的取值。有时,为了便于分析,还可以定义对数似然函数 H(θ)=lnl(θ)=ln∏Ni=1p(xi|θ)=∑Ni=1lnp(xi|θ)(34) 容易证明,使对数似然函数最大的θ值也使似然函数最大。 3.2.2最大似然估计的求解 现在再来看式(32)似然函数公式,虽然公式复杂,但其中的函数形式p(·)是已知的,其中的xi也都是已知的,未知量只有θ。在似然函数满足连续、可微的条件下,如果θ是一维变量,即只有一个待估计参数,其最大似然估计量就是如下微分方程的解 dl(θ)dθ=0(35) 或 dH(θ)dθ=0(36) 更一般地,当θ=[θ1,…,θs]T是由多个未知参数组成的向量时,求解似然函数的最大值就需要对θ的每一维分别求偏导,即用下面的梯度算子 θ=θ1,…,θsT(37) 来对似然函数或者对数似然函数求梯度并令其等于零 θl(θ)=0(38) 或 θH(θ)=∑Ni=1θlnp(xi|θ)=0(39) 得到s个方程,方程组的解就是对数似然函数的极值点。需要注意,在某些情况下,似然函数可能有多个极值,此时上述方程组可能有多个解,其中使得似然函数最大的那个解才是最大似然估计量。 并不是所有的概率密度形式都可以用上面的方法求得最大似然估计。例如,已知一维随机变量x服从均匀分布 p(x|θ)=1θ2-θ1θ1<x<θ2 0其他(310) 其中分布的参数θ1、θ2未知。从总体分布中独立抽取了N个样本x1,x2,…,xN,则似然函数为 l(θ)=p(X|θ)=p(x1,x2,…,xN|θ1,θ2)=1(θ2-θ1)N 0(311) 对数似然函数为 H(θ)=-Nln(θ2-θ1)(312) 通过式(39)求 Hθ1=N·1θ2-θ1(313) Hθ2=-N·1θ2-θ1 (314) 从式(313)、式(314)方程组中解出的参数θ1和θ2至少有一个为无穷大,这是无意义的结果。造成这种困难的原因是似然函数在最大值的地方没有零斜率,所以必须用其他方法来找最大值。从式(310)看出,当θ2-θ1越小时,则似然函数越大。而在给定一个有N个观察值x1,x2,…,xN的样本集中,如果用x′表示观察值中最小的一个,用x″表示观察值中最大的一个,显然θ1不能大于x′,θ2不能小于x″,因此θ2-θ1的最小可能值是x″-x′,这时θ的最大似然估计量显然是 θ^1=x′(315) θ^2=x″(316) 3.2.3正态分布正态分布下的最大似然估计 这里仅以单变量正态分布情况下估计其均值和方差为例来说明最大似然估计的用法。我们知道,单变量正态分布的形式为 p(x|θ)=12πσexp-12x-μσ2(317) 其中均值μ和方差σ2为未知参数,即我们要估计的参数为θ=[θ1,θ2]T=[μ,σ2]T,用于估计的样本仍然是X={x1,x2,…,xN}。根据式(311),最大似然估计应该是下面方程组的解 θH(θ)=∑Nk=1θlnp(xk|θ)=0(318) 从正态分布式(317)可以得到 lnp(xk|θ)=-12ln2πθ2-12θ2(xk-θ1)2(319) 分别对两个未知参数求偏导,得到 θlnp(xk|θ)=1θ2(xk-θ1) -12θ2+12θ22(xk-θ1)2(320) 因此,最大似然估计应该是以下方程组的解 ∑Nk=11θ^2(xk-θ^1)=0 -∑Nk=11θ^2+∑Nk=1(xk-θ^1)2θ^22=0(321) 容易解得 μ^=θ^1=1N∑Nk=1xk(322) δ^2=θ^2=1N∑Nk=1(xk-μ^)2(323) 这正是人们经常使用的对均值和方差的估计,它们是对正态分布样本的均值和方差的最大似然估计。 对于多元正态分布,分析原理和上面相同,只是公式略微复杂一些,结论也和单变量情况很相似,即多元正态分布的均值和方差的最大似然估计是 μ^=1N∑Ni=1xi(324) Σ^=1N∑Ni=1(xi-μ^)(xi-μ^)T(325) 从以上结果可以得出结论: 均值向量μ的最大似然估计是样本均值。协方差矩阵Σ的最大似然估计是N个矩阵(xk-μ^)(xk-μ^)T的算术平均。由于真正的协方差矩阵是随机矩阵(x-μ)(x-μ)T的期望值,所以这个结果是非常令人满意的。最大似然估计量是平方误差一致和简单一致估计量,但不一定都是无偏估计量。上例中μ^是无偏的,而Σ^就不是无偏的,Σ^的无偏估计为1N-1∑Nk=1(xk-μ^)(xk-μ^)T,其证明作为习题留给读者。 3.3贝叶斯估计贝叶斯估计与贝叶斯学习贝叶斯学习 3.3贝叶斯估计与贝叶斯学习 贝叶斯估计(Bayesian Estimation)是概率密度估计中的另一类主要的参数估计方法,其结果在很多情况下与最大似然法相同或几乎相同,但是两种方法对问题的处理视角是不同的,在应用上也各自有各自的特点。一个根本的区别是,最大似然估计是把待估计的参数当作未知但固定的量,要做的是根据观测数据估计这个量的取值; 而贝叶斯估计则是把待估计的参数本身也看作是随机变量,要做的是根据观测数据对参数的分布进行估计,除了观测数据外,还可以考虑参数的先验分布。贝叶斯学习(Bayesian learning)则是把贝叶斯估计的原理用于直接从数据对概率密度函数进行迭代估计。 3.3.1贝叶斯估计 可以把概率密度函数的参数估计问题看作一个贝叶斯决策问题,但是这里要决策的不是离散的类别,而是参数的取值,是在连续的空间里做决策。 把待估计参数θ看作具有先验分布密度p(θ)的随机变量,其取值与样本集X有关,我们要做的是根据样本集X=x1,x2,…,xN估计最优的θ(记作θ*)。 在用于分类的贝叶斯决策中,最优的条件可以是最小错误率或者最小风险。在这里,对连续变量θ,我们假定把它估计为θ^所带来的损失为λ(θ^,θ),也称作损失函数损失函数。 设样本的取值空间是Ed,参数的取值空间是Θ,那么,当用θ^来作为估计时总期望风险就是 R=∫Ed∫Θλ(θ^,θ)p(x,θ)dθdx =∫Ed∫Θλ(θ^,θ)p(θ|x)p(x)dθdx(326) 我们定义在样本x下的条件风险为 R(θ^|x)=∫Θλ(θ^,θ)p(θ|x)dθ(327) 那么,式(326)就可以写成 R=∫EdR(θ^|x)p(x)dx(328) 现在的目标是对期望风险求最小。与贝叶斯分类决策时相似,这里的期望风险也是在所有可能的x情况下的条件风险的积分,而条件风险又都是非负的,所以求期望风险最小就等价于对所有可能的x求条件风险最小。在有限样本集合X={x1,x2,…,xN}的情况下,我们所能做的就是对所有的样本求条件风险最小,即 θ*=arg minθ^ R(θ^|X)=∫Θλ(θ^,θ)p(θ|X)dθ(329) 在决策分类时,需要事先定义决策表即损失表,而连续情况下,需要定义损失函数。最常用的损失函数是平方误差损失函数,即 λ(θ^,θ)=(θ-θ^)2(330) 可以证明,如果采用平方误差损失函数,则在样本x条件下θ的贝叶斯估计量θ*是在给定x下θ的条件期望,即 θ*=E[θ|x]=∫Θθp(θ|x)dθ (331) 同样,在给定样本集X下,θ的贝叶斯估计量是 θ*=E[θ|X]=∫Θθp(θ|X)dθ(332) 这样,在最小平方误差损失函数下,贝叶斯估计的步骤是: (1) 根据对问题的认识或者猜测确定θ的先验分布密度p(θ)。 (2) 由于样本是独立同分布的,而且已知样本密度函数的形式p(x|θ),可以形式上求出样本集的联合分布为 p(X|θ)=∏Ni=1p(xi|θ)(333) 其中θ是变量。 (3) 利用贝叶斯公式求θ的后验概率分布 p(θ|X)=p(X|θ)p(θ)∫Θp(X|θ)p(θ)dθ(334) (4) 根据式(332),θ的贝叶斯估计量是 θ*=∫Θθp(θ|X)dθ(335) 在贝叶斯估计中,样本的概率密度函数p(x|θ)的形式是已知的,参数的先验分布密度p(θ)只有在某些特殊形式下才能使式(334)的后验概率形式上方便计算。特别地,对于给定的概率密度函数p(x|θ)模型,如果先验密度p(θ)能够使参数的后验分布p(θ|X)具有与p(x|θ)相同的形式,则这样的先验密度函数形式称作与概率模型p(x|θ)共轭(conjugate)。在Bernardo和Smith的教材(Bernardo,J.M.and Smith,A.F.M.Bayesian Theory.Chichester: Wiley,1994)中给出了一些常用的共轭先验概率密度模型的例子,实际中最常用的是在p(x|θ)为正态分布时p(θ)也为正态分布。 应注意到,我们本来的目的并不是估计概率密度参数,而是估计样本的概率密度函数p(x|X)本身,因为假定概率密度函数的形式已知,才转化为估计密度函数中的参数的问题。实际上,在上面介绍的贝叶斯估计贝叶斯估计框架下,从式(334)得到了参数的后验概率后就可以不必求对参数的估计,而是直接得到样本的概率密度函数 p(x|X)=∫Θp(x|θ)p(θ|X)dθ(336) 可以这样直观地理解式(336): 参数θ是随机变量,它有一定的分布,而要估计的概率密度p(x|X)就是所有可能的参数取值下的样本概率密度的加权平均,而这个加权就是在观测样本下估计出的参数θ的后验概率。在式(334)给出的参数分布估计中,决定分布形状的是p(X|θ)p(θ),即 p(θ|X)~p(X|θ)p(θ)(337) 分母只是对估计出的分布的归一化因子,保证概率密度函数下的积分为1。可以看到,p(θ|X)是由两项决定的,一项就是上一节定义的似然函数p(X|θ),它反映了在不同参数取值下得到观测样本的可能性; 另一项是参数取值的先验概率p(θ),它反映了对参数分布的先验知识或者主观猜测。 极端情况下,如果完全没有先验知识,即认为p(θ)是均匀分布,则p(θ|X)就完全取决于p(X|θ)。与最大似然估计不同的是,这里并没有直接把似然函数最大或者是后验概率最大的值拿来当作对样本概率密度参数的估计,而是根据把所有可能的参数值都考虑进来,用它们的似然函数作为加权来平均出一个对参数的估计(即式(335))或者对样本概率密度函数的估计(式(336))。 另一个极端情况是,如果先验知识非常强,p(θ)就是在某一特定取值θ0上的一个脉冲函数(即p(θ0)=δ(θ0)),则由式(337)可知,除非在θ0的似然函数为0,否则最后的估计就是θ0,样本不再起作用。 通常情况下,p(θ)不是均匀分布也不是脉冲函数,则参数的后验概率就受似然函数和先验概率的共同作用。一种常见的情况是,似然函数在其最大值θ=θ^附近会有一个尖峰,那么如果先验概率在最大似然估计处不为零且变化比较平缓,则参数的后验概率p(θ|X)就会集中在θ=θ^附近,此时式(335)得到的贝叶斯估计就与最大似然估计接近,式(336)估计出的样本密度基本上也就是在最大似然估计下的样本密度。 3.3.2贝叶斯学习 现在来考虑更为一般的情况,即根据观测样本用式(335)来估计样本概率密度函数的参数。为了反映样本的数目,把样本集重新记作XN={x1,x2,…,xN},式(335)重写如下 θ*=∫Θθp(θ|XN)dθ(338) 其中 p(θ|XN)=p(XN|θ)p(θ)∫Θp(XN|θ)p(θ)dθ(339) 当N>1时,有 p(XN|θ)=p(xN|θ)p(XN-1|θ)(340) 把它代入式(339),可以得到如下的递推公式 p(θ|XN)=p(xN|θ)p(θ|XN-1)∫p(xN|θ)p(θ|XN-1)dθ(341) 为了形式统一,把先验概率记作p(θ|X0)=p(θ),表示在没有样本情况下的概率密度估计。根据式(341),随着样本数的增加,可以得到一系列对概率密度函数参数的估计 p(θ),p(θ|x1),p(θ|x1,x2),…,p(θ|x1,x2,…,xN),…(342) 称作递推的贝叶斯估计。如果随着样本数的增加,式(342)的后验概率序列逐渐尖锐,逐步趋向于以θ的真实值为中心的一个尖峰,当样本无穷多时收敛于在参数真实值上的脉冲函数,则这一过程称作贝叶斯学习。 此时,用式(336)估计的样本概率密度函数也逼近真实的密度函数 p(x|XN→∞)=p(x)(343) 3.3.3正态分布时的贝叶斯估计 下面以最简单的一维正态分布模型为例来说明贝叶斯估计的应用。假设模型的均值μ是待估计的参数,方差为σ2为已知,我们可以把分布密度写为 p(x|μ)=12πσexp-12σ2(x-μ)2(344) 假定均值μ的先验分布也是正态分布,其均值为μ0、方差为σ20,即 p(μ)=12πσ0exp-12σ20(μ-μ0)2(345) 用式(334)来对均值μ进行估计 p(μ|X)=p(X|μ)p(μ)∫Θp(X|μ)p(μ)dμ(346) 已经知道,这里的分母只是用来对估计出的后验概率进行归一化的常数项,可以暂时不考虑。现在来计算式(346)右边的分子部分。 p(X|μ)p(μ)=p(μ)∏Ni=1p(xi|μ) =12πσexp-12μ-μ0σ02∏Ni=112πσexp-12xi-μσ2 把所有不依赖于μ的量都写入一个常数中,上式可以整理为 p(X|μ)p(μ)=αexp-12μ-μNσN2(347) 可见p(μ|X)也是一个正态分布,可以得到 p(μ|X)=12πσNexp-12μ-μNσN2(348) 其中的参数满足 1σ2N=1σ20+Nσ2(349) μN=σ2Nμ0σ20+∑Ni=1xiσ2(350) 进一步整理后得 μN=Nσ20Nσ20+σ2mN+σ2Nσ20+σ2μ0(351) σ2N=σ20σ2Nσ20+σ2(352) 其中,mN=1N∑Ni=1xi是所有观测样本的算术平均。 所以,贝叶斯估计告诉我们,待估计的样本密度函数的均值参数服从均值为μN、方差为σ2N的正态分布。显然,可以用式(335)得到参数的贝叶斯估计值,即 μ^=∫μp(μ|X)dμ=∫μ2πσNexp-12μ-μNσN2dμ=μN(353) 在式(351)中,正态分布下贝叶斯估计的结果是由两项组成的,一项是样本的算术平均,另一项是对均值的先验认识。当样本数目趋于无穷大时,第一项的系数趋于1而第二项的系数趋于0,即估计的均值就是样本的算术平均。这与最大似然估计是一致的。当样本数目有限时,如果先验知识非常确定,那么先验分布的方差σ20就很小,此时第一项的系数就很小,而第二项的系数接近1,估计主要由先验知识来决定。一般情况下,均值的贝叶斯估计是在样本算术平均与先验分布均值之间进行加权平均。 从这里可以看到,贝叶斯估计的优势不但在于使用样本中提供的信息进行估计,而且能够很好地把关于待估计参数的先验知识融合进来,并且能够根据数据量大小和先验知识的确定程度来调和两部分信息的相对贡献,这在很多实际问题中是非常有用的。 在得到式(348)的后验分布后,我们也可以直接用式(336)求出样本的密度函数 p(x|X)=12πσ2+σ2Nexp-12x-μNσ2+σ2N2~NμN,σ2+σ2N(354) 其中,μN、σ2N仍然如式(352)、式(353)。可以看到,虽然我们的条件是已知方差σ2,但是由于均值是估计值μN,贝叶斯估计得到的分布密度函数方差增加了,变成σ2+σ2N,而所增加的项σ2N在当样本趋于无穷大时趋于零。 3.3.4其他分布的情况 需要说明的是,在一般情况下,在求出参数的后验概率分布p(θ|X)后,计算式(335)的数学期望和式(336)的积分并不是非常容易的,对于某些概率模型甚至会非常困难。在这种情况下,比较简单的做法是直接根据p(θ|X)选取一个参数值作为估计,例如选择后验概率最大的参数值,但在很多分布下最大值与数学期望的差距可能会很大。一种更完善的做法是,用所谓吉布斯采样(Gibbs Sampling)等方法来对参数的后验概率分布p(θ|X)进行随机采样,用采样得到的参数的算术平均来估算式(335)的数学期望,而这种采样并不需要计算式(334)中的分母部分(在很多分布下这一归一化因子的计算并不容易),而是可以根据与p(θ|X)成正比的p(X|θ)p(θ)进行。有关这方面的内容请参阅Andrew Webb的教材Statistical Pattern Recognition中有关内容或者关于马尔可夫链蒙特卡罗MCMC(Markov Chain Monte Carlo)的有关教材或专著。 3.4概率密度估计的非参数方法 3.4概率密度估计的非参数方法 最大似然方法和贝叶斯方法都属于参数化的估计方法,要求待估计的概率密度函数形式已知,只是利用样本来估计函数中的某些参数。但是,在很多情况下,我们对样本的分布并没有充分的了解,无法事先给出密度函数的形式,而且有些样本分布的情况也很难用简单的函数来描述。在这种情况下,就需要非参数估计非参数估计,即不对概率密度函数的形式作任何假设,而是直接用样本估计出整个函数。当然,这种估计只能是用数值方法取得,无法得到完美的封闭函数形式。从另外的角度来看,概率密度函数的参数估计实际是在指定的一类函数中选择一个函数作为对未知函数的估计,而非参数估计则可以看作从所有可能的函数中进行的一种选择。 3.4.1非参数估计的基本原理与直方图方法 直方图直方图(histogram)方法是最简单直观的非参数估计方法,也是日常人们最常用的对数据进行统计分析的方法。图32给出了一个直方图的例子。 图32直方图举例 进行直方图估计的做法是: (1) 把样本x的每个分量在其取值范围内分成k个等间隔的小窗。如果x是d维向量,则这种分割就会得到kd个小体积或者称作小舱,每个小舱的体积记作V。 (2) 统计落入每个小舱内的样本数目qi。 (3) 把每个小舱内的概率密度看作是常数,并用qi/(NV)作为其估计值,其中N为样本总数。 下面来分析非参数估计的基本原理。我们的问题是: 已知样本集X={x1,…,xN}中的样本是从服从密度函数p(x)的总体中独立抽取出来的,求p(x)的估计p^(x)。与参数估计时相同,这里不考虑类别,即假设样本都是来自同一个类别,对不同类别只需要分别进行估计即可。 考虑在样本所在空间的某个小区域R,某个随机向量落入这个小区域的概率是 PR=∫Rp(x)dx(355) 根据二项分布,在样本集X中恰好有k个落入小区域R的概率是 Pk=CkNPkR(1-PR)N-k(356) 其中CkN表示在N个样本中取k个的组合数。k的期望值是 E[k]=NPR(357) 而且k的众数(概率最大的取值)是 m=[(N+1)PR](358) 其中[]表示取整数。因此,当小区域中实际落入了k个样本时,PR的一个很好的估计是 P^R=kN(359) 当p(x)连续且小区域R的体积V足够小时,可以假定在该小区域范围内p(x)是常数,则式(355)可近似为 PR=∫Rp(x)dx=p(x)V(360) 用式(359)的估计代入到式(360)中,可得: 在小区域R的范围内 p^(x)=kNV(361) 这就是在上面的直方图中使用的对小舱内概率密度的估计。 在上面的直方图估计中,我们采用的是把特征空间在样本取值范围内等分的做法。可以设想,小舱的选择是与估计的效果密切相连的: 如果小舱选择过大,则假设p(x)在小舱内为常数的做法就显得粗糙,导致最终估计出的密度函数也非常粗糙,如图33(a)的例子所示; 而另一方面,如果小舱过小,则有些小舱内可能就会没有样本或者很少样本,导致估计出的概率密度函数很不连续,如图33(b)所示。 图33小窗宽度对直方图估计的影响示意 所以,小舱的选择应该与样本总数相适应。理论上讲,假定样本总数是n,小舱的体积为Vn,在x附近位置上落入小舱的样本个数是kn,那么当样本趋于无穷多时p^(x)收敛于p(x)的条件是 (1) limn→∞Vn=0, (2) limn→∞kn=∞, (3) limn→∞knn=0 (362) 直观的解释: 随着样本数的增加,小舱体积应该尽可能小,同时又必须保证小舱内有充分多的样本,但每个小舱内的样本数又必须是总样本数中很小的一部分。 我们自然可以想到,小舱内有多少样本不但与小舱体积有关,还与样本的分布有关。在有限数目的样本下,如果所有小舱的体积相同,那么就有可能在样本密度大的地方一个小舱里有很多样本,而在密度小的地方则可能一个小舱里只有很少甚至没有样本,这样就可能导致密度的估计在样本密度不同的地方表现不一致。因此,固定小舱宽度的直方图方法只是最简单的非参数估计方法,要想得到更好的估计,需要采用能够根据样本分布情况调整小舱体积的方法。 3.4.2kN近邻估计方法 kN近邻估计就是一种采用可变大小的小舱的密度估计方法,基本做法是: 根据总样本确定一个参数kN,即在总样本数为N时我们要求每个小舱内拥有的样本个数。在求x处的密度估计p^(x)时,我们调整包含x的小舱的体积,直到小舱内恰好落入kN个样本,并用式(360)来估算p^(x),即 p^(x)=kN/NV(363) 这样,在样本密度比较高的区域小舱的体积就会比较小,而在密度低的区域则小舱体积自动增大,这样就能够比较好地兼顾在高密度区域估计的分辨率和在低密度区域估计的连续性。 为了取得好的估计效果,仍然需要根据式(362)来选择kN与N的关系,例如可以选取为kN~kN,k为某个常数。 kN近邻估计与简单的直方图方法相比还有一个不同,就是kN近邻估计并不是把x的取值范围划分为若干个区域,而是在x的取值范围内以每一点为小舱中心用式(363)进行估计,如图34所示。图35给出了两个一维情况下在不同样本数目时kN近邻估计效果的例子。 图34kN法的窗口宽度与样本密度的关系示意 图35不同样本数和不同参数下 kN法估计的效果举例 3.4.3Parzen窗法Parzen窗法 再来看式(361)的基本公式p^(x)=kNV。在采用固定小舱体积的情况下,我们也可以像kN近邻估计那样用滑动的小舱来估计每个点上的概率密度,而不是像直方图中那样仅仅在每个小舱内估计平均密度。 假设x∈Rd是d维特征向量,并假设每个小舱是一个超立方体,它在每一维的棱长都为h,则小舱的体积是 V=hd(364) 要计算每个小舱内落入的样本数目,可以定义如下的d维单位方窗函数 φ([u1,u2,…,ud]T)=1|uj|≤12,j=1,2,…,d 0其他(365) 该函数在以原点为中心的d维单位超正方体内取值为1,而其他地方取值都为0。对于每一个x,要考查某个样本xi是否在这个x为中心、h为棱长的立方小舱内,就可以通过计算φx-xih来进行。现在共有N个观测样本{x1,x2,…,xN},那么落入以x为中心的超立方体内的样本数就可以写成 kN=∑Ni=1φx-xih(366) 把它代入式(361)中,可以得到对于任意一点x的密度估计的表达式 p^(x)=1NV∑Ni=1φx-xih(367) 或者写成 p^(x)=1N∑Ni=11V φx-xih(368) 还可以从另外的角度来理解式(368): 定义核函数核函数(也称窗函数) K(x,xi)=1V φx-xih(369) 它反映了一个观测样本xi对在x处的概率密度估计的贡献,与样本xi与x的距离有关,也可记作K(x-xi)。概率密度估计就是在每一点上把所有观测样本的贡献进行平均,即 p^(x)=1N∑Ni=1K(x,xi)(370) 一个基本的要求是,这样估计出的函数至少应该满足概率密度函数的基本条件,即函数值应该非负且积分为1。显然,这只需核函数本身满足密度函数的要求即可,即 K(x,xi)≥0且∫K(x,xi)dx=1(371) 容易验证,由式(369)和式(365)定义的立方体核函数满足这一条件。 这种用窗函数(核函数)估计概率密度的方法称作Parzen窗方法(Parzen window method)估计或核密度估计核密度估计(kernel density estimation)。Parzen窗估计也可以看作用核函数对样本在取值空间中进行插值。 由式(369)定义的核函数称作方窗函数,还有多种其他核函数。下面列举几种常见的核函数。 (1) 方窗 k(x,xi)=1hd| xj-xji|≤h2,j=1,2,…,d 0其他(372) 其中,h为超立方体的棱长。这种写法与式(369)定义是相同的。 (2) 高斯窗(正态窗) k(x,xi)=1(2π)dρ2d|Q|exp-12(x-xi)TQ-1(x-xi)ρ2(373) 即以样本xi为均值、协方差矩阵为Σ=ρ2Q的正态分布函数。一维情况下则为 k(x,xi)=12πσexp-(x-xi)22σ2(374) (3) 超球窗 k(x,xi)=V-1‖x-xi‖≤ρ 0其他(375) 其中V是超球体的体积,ρ是超球体半径。 在这些窗函数中,都有一个表示窗口宽度的参数,也称作平滑参数,它反映了一个样本对多大范围内的密度估计产生影响。 图36给出用高斯窗估计一个概率密度函数的例子。 图36Parzen窗估计概率密度函数示例 当被估计的密度函数连续时,在核函数及其参数满足一定的条件下,Parzen窗估计是渐近无偏和平方误差一致的。这些条件主要是: 对称且满足式(371)的密度函数条件、有界(不能是无穷大)、核函数取值随着距离的减小而迅速减小、对应小舱的体积应随着样本数的增加而趋于零,但是不能缩减速度太快,即慢于1/N趋于零的速度。 图37给出了两种不同分布情况下用不同的参数和高斯窗进行估计的结果,其中采用σ=h1/N,h1为例子中可调节的参量。 图37不同样本数和不同参数下Parzen窗估计的效果示例 作为非参数方法的共同问题是对样本数目需求较大,只要样本数目足够大,总可以保证收敛于任何复杂的未知密度,但是计算量和存储量都比较大。当样本数很少时,如果能够对密度函数有先验认识,则参数估计方法能取得更好的估计效果。