第3章
CHAPTER 3


概率密度函数的参数估计








第2章讨论了贝叶斯决策理论，在采用贝叶斯决策理论进行分类决策时，需要计算后验概率P(ωi|X)，或者需要事先知道各类的先验概率P(ωi)和样本的类条件概率密度函数p(X|ωi)，但实际应用中先验概率和类条件概率密度函数往往是未知的。通常，对研究的对象只有一些模糊性的知识，或者通过实验采样而得到的一些样本。这就需要根据已有的样本，利用统计推断中的估计理论对样本的分布做出估计，然后将估计值当成真实值来使用。在模式识别问题中，先验概率的估计相对比较容易，它可以由各类样本在总体样本集中所占的比例进行估计。但类条件概率密度函数的估计却比较困难，从样本出发估计其函数形式和参数，这就是本章要讨论的参数估计问题。

3.1概率密度函数估计概述

所谓的概率密度函数估计是已知某类别ωi的样本Xi(i=1,2,…,N)，采用某种规则估计出样本所属类的概率函数p(X|ωi)。从估计的方法来讲，可分为参数估计和非参数估计。参数估计是先假定样本的类条件概率密度函数p(X|ωi)的类型已知，如服从正态分布、二项分布，再用已知类别的学习样本估计函数里的未知参数θ，这项工作也叫训练或学习。参数估计的方法通常采用的是最大似然估计方法和贝叶斯估计方法。非参数估计则是类条件概率密度函数的形式未知，直接用已知类别的学习样本去估计函数的数学模型。非参数估计的方法通常采用的是Parzen窗法、kN近邻法等。

为了便于理解，首先介绍参数估计中的一些基本概念。

(1) 统计量。假如概率密度函数的形式已知，但表征函数的参数θ未知，则可将θ的估计值构造成样本Xi(i=1,2,…,N)的某种函数，这种函数称为统计量。参数估计的任务，就是利用样本求出参数θ的估计值θ^=θ(X1,X2,…,XN)。

(2) 参数空间。参数θ的取值范围称为参数空间，书中用Θ来表示。

(3) 点估计、估计量和估计值。构造一统计量作为未知参数θ的估计，称为点估计，由样本(X1,X2,…,XN)作为自变量计算出来的θ^值称为估计值，θ^称为估计量。

(4) 区间估计。通过从总体中抽取的样本，根据一定的正确度与精确度的要求，构造出适当的区间，作为未知参数的真值所在范围的估计。

下面分别介绍最大似然估计、贝叶斯估计、贝叶斯学习三种参数估计方法，以及Parzen窗法和kN近邻法两种非参数估计方法。






第5集
微课视频



3.2最大似然估计

对c类问题，设类别ωi的概率密度函数p(X|ωi)的形式已知，但表征该函数的参数未知，记为θi。从ωi中独立抽取N个样本，如果能从这N个样本中推断出θi的估计值θ^i，则完成了概率密度函数p(X|ωi)的估计。为了强调p(X|ωi)与参数θi的关联性，也可把概率密度函数写成p(X|ωi,θi)。例如，如果已知某一类别ωi概率密度函数服从正态分布，则未知参数θi包含了表征该函数的均值 μi和协方差Σi的全部信息，对参数θi的估计，实质上就是对正态函数的均值 μi和协方差Σi的估计。下面首先给出似然函数的定义，然后从似然函数出发，讨论最大似然估计的原理。






1. 似然函数

从ωi类中抽取N个样本X(N)=X1,X2,…,XN,由于这N个样本均来自ωi类，因此可将其概率密度函数p(X|ωi,θi)简化为p(X|θ)，则称这N个样本的联合概率密度函数p(X(N),θ)为相对于样本集X(N)的θ的似然函数。由于θ是概率密度函数的一个确定性参数集，因此概率密度函数p(X(N),θ)实际上就是条件概率p(XN|θ)。如果N个样本为独立抽取，似然函数可表示为

p(X(N)|θ)=p(X1,X2,…,XN|θ)=∏Nk=1p(Xk|θ)(31)

式(31)是在参数θ下观测到的样本集X(N)的概率(联合分布)密度。

2. 最大似然估计

从ωi类中独立抽取N个样本X(N)=X1,X2,…,XN，那么这N个样本最有可能来自哪个概率密度函数，或者说与这N个样本最匹配的未知参数θ是什么。这是最大似然估计要解决的问题，它的主要思想是，给定样本集X(N)=X1,X2,…,XN，通过极大化似然函数p(X(N)|θ)去求与样本匹配的参数θ，θ的最大似然估计量θ^就是使似然函数达到最大的估计量，图31所示为θ为一维时的最大似然估计。由dp(X(N)|θ)dθ=0可求得解。



图31θ为一维时的最大似然估计


由于对数函数具有单调性，为了便于分析，对似然函数取对数

H(θ)=lnp(X(N)|θ)(32)

显然，当估计量θ^使数函数取最大值时，似然函数达到最大值，θ的最大似然估计是下面微分方程的解： 

dH(θ)dθ=0(33)

设ωi类的概率密度函数包含q个未知参数，则θ为q维向量

θ=［θ1,θ2,…,θq］T(34)

此时

H(θ)=lnp(X(N)|θ)=∑Nk=1lnp(Xk|θ)(35)

式(33)可表示为

θ∑Nk=1lnp(Xk|θ)=0(36)

即


∑Nk=1θ1lnp(Xk|θ)=0
∑Nk=1θ2lnp(Xk|θ)=0

∑Nk=1θqlnp(Xk|θ)=0
(37)

求解式(37)微分方程组，可得到θ的最大似然估计值 θ^。

【例3.1】设从ωi中抽取了N个样本，表示为X(N)，这N个样本是从一维正态分布概率密度函数pX|ωi［或p(X(N)|θ)］总体中独立抽取的，用最大似然估计方法，估计正态分布的均值和协方差。

【解】p(X|ωi)可表示为pX|θ~Nμ,σ2，其中，θ=［θ1,θ2］T，θ1=μ，θ2=σ2。因为X(N)是从ωi中独立抽取的N个样本，则θ的似然函数为

p(X(N)|θ)=∏Nk=1p(Xk|θ)

式中，p(Xk|θ)=12πσexp－Xk－μ22σ2。

取似然函数的对数，得

lnp(Xk|θ)=－12ln(2πσ2)－Xk－μ22σ2

函数lnp(Xk|θ)分别对θ1和θ2求导，并令导数为零，即


∑Nk=1θ1lnp(Xk|θ)=∑Nk=1Xk－θ1θ2=0
∑Nk=1θ2lnp(Xk|θ)=∑Nk=1
－12θ2+(Xk－θ1)22θ22=0


由以上方程组解得均值和方差的估计量为


μ^=θ^1=1N∑Nk=1Xk(38)
σ^2=θ^2=1N∑Nk=1(Xk－μ^)2(39)


对于一般的多维正态分布情况，用类似例3.1的方法，可以求得其最大似然估计值为


μ^i=1N∑Nk=1Xk(310)
Σ^i=1N∑Nk=1(Xk－μ^i)(Xk－μ^i)T(311)


式(310)与式(311)表明，在多元正态分布情况下，均值向量的最大似然估计是样本的算术平均值，而协方差矩阵的最大似然估计是N个矩阵的(Xk－μ^i)(Xk－μ^i)T的算术平均值。

3.3贝叶斯估计与贝叶斯学习

1. 贝叶斯估计

贝叶斯估计可描述为给定样本集X(N)=X1,X2,…,XN，对样本的概率密度函数的真实参数θ进行估计，使其估计值θ^带来的贝叶斯风险最小。回顾上一章的最小风险贝叶斯决策，可以看出贝叶斯决策和贝叶斯估计都是以贝叶斯风险最小为基础，只是要解决的问题不同，前者是要判决样本X的类别归属，而后者是估计样本集X(N)所属总体分布的参数，本质上二者是统一的。贝叶斯决策和贝叶斯估计各变量的对应关系如表31所示。


表31贝叶斯决策和贝叶斯估计各变量的对应关系




决 策 问 题估 计 问 题


样本x样本集X(N)
决策αi估计量θ^
真实状态ωj真实参数θ
状态空间A是离散空间参数空间Θ
先验概率P(ωj)参数的先验分布p(θ)


在第2章研究分类问题时，我们用式(211)定义了条件平均风险： 

R(αi|X)=E［L(αi|ωj)］=∑cj=1L(αi|ωj)·P(ωj|X),i=1,2,…,a

参考上式，并对照表31中贝叶斯决策和贝叶斯估计各变量的对应关系，可以定义在观测样本集X(N)=X1,X2,…,XN的条件下，用θ^作为θ的估计的期望损失为

R(θ^|X(N))=∫ΘL(θ^,θ)p(θ|X(N))dθ(312)

式中，L(θ^,θ)为用θ^代替θ所造成的损失，Θ为参数空间。

考虑X(N)的各种取值，应该求R(θ^|X(N))在空间ΩN=Ω×Ω×…×Ω中的期望，即

R=∫ΩNR(θ^|X(N))p(X(N))dX(N)(313)

将式(312)代入上式，得

R=∫ΩN∫ΘL(θ^,θ)p(θ|X(N))p(X(N))dθdX(N)(314)

使R最小的参数θ的估计值θ^即贝叶斯估计。显然，损失函数L(θ^,θ)对θ^的求解有重要影响，当选用不同形式的损失函数时，所得到的贝叶斯估计值θ^也不同。当损失函数为二次函数时,有

L(θ^,θ)=θ－θ^Tθ－θ^(315)

可证明θ^的求解公式如下： 


θ^=∫Θθp(θ|X(N))dθ(316)


上式表明，θ的最小方差贝叶斯估计是观测样本集X(N)条件下的θ的条件期望。

综上所述，观测到一组样本X(N)，通过似然函数pX(N)|θ并利用贝叶斯公式将随机变量θ的先验概率密度p(θ)转换为后验概率密度，然后根据θ的后验概率密度求出估计量θ^，具体步骤如下。

(1) 确定θ的先验概率密度p(θ)。

(2) 由样本集X(N)={X1,X2,…,XN}求出样本的联合概率密度函数也就是θ的似然函数pX(N)|θ。

(3) 利用贝叶斯公式求出θ的后验概率密度：

p(θ|X(N))=p(X(N)|θ)p(θ)∫Θp(X(N)|θ)p(θ)dθ(317)

(4) 根据式(316)求贝叶斯估计量θ^。

在步骤(2)涉及pX(N)|θ的求解，当样本的类概率密度函数的类型已知时，由于样本X1,X2,…,XN为独立抽取，因此有

pX(N)|θ=pX1,X2,…,XNθ=∏Ni=1p(Xiθ)

2. 贝叶斯学习

贝叶斯学习是指利用θ的先验概率密度p(θ)及样本提供的信息递推求出θ的后验概率密度p(θ|X(N))，根据后验概率密度直接求出类概率密度函数p(X|X(N))。因此，贝叶斯学习和贝叶斯估计的前提条件完全相同，区别在于当求出后验概率密度p(θ|X(N))后，贝叶斯学习没有对参数θ进行估计，而是直接进行总体概率密度的推断得到p(X|X(N))。所以，贝叶斯学习的前三步与贝叶斯估计完全一致，最后p(X|X(N))可由迭代计算完成。迭代计算式的推导如下。

p(X|ωi)由未知参数θ确定，可写为p(X|ωi)=pX|θ，假定X(N)={X1,X2,…,XN}是独立抽取的ωi类的一组样本，设θ的后验概率密度函数为p(θ|X(N))，根据贝叶斯公式有

p(θ|X(N))=p(X(N)|θ)p(θ)∫Θp(X(N)|θ)p(θ)dθ(318)

由条件独立可知，当N>1时

p(X(N)|θ)=p(XN|θ)p(X(N－1)|θ)(319)

式中，X(N－1)表示除样本XN以外其余样本的集合。

将式(319)代入式(318)得

p(θ|X(N))=p(XN|θ)p(X(N－1)|θ)p(θ)∫Θp(XN|θ)p(X(N－1)|θ)p(θ)dθ(320)

类似地，由式(318)也可推导出

p(θ|X(N－1))=p(X(N－1)|θ)p(θ)∫Θp(X(N－1)|θ)p(θ)dθ(321)

将式(321)代入式(320)得

p(θ|X(N))=p(XN|θ)p(θ|X(N－1))∫Θp(XN|θ)p(θ|X(N－1))dθ
(322)

式(322)就是利用样本集X(N)估计p(θ|X(N))的迭代计算方法。对于参数估计的递推贝叶斯方法，其迭代过程是贝叶斯学习的过程。下面简述迭代式的使用。

(1) 根据先验知识得到θ的先验概率密度函数的初始估计p(θ)。相当于N=0时(X(N)=X(0))密度函数的一个估计。

(2) 用X1对初始的p(θ)进行修改，根据式(322)，令N=1，得

p(θ|X(1))=p(θ|X1)=p(X1|θ)p(θ)∫Θp(X1|θ)p(θ)dθ(323)

pX1|θ根据式p(X|ωi)=pX|θ计算得到。

(3) 给出X2，对用X1估计的结果进行修改

p(θ|X(2))=p(θ|X1,X2)=p(X2|θ)p(θ|X(1))∫Θp(X2|θ)p(θ|X(1))dθ(324)

(4) 逐次给出X3,X4,…,XN，每次在前一次的基础上进行修改，p(θ|X(N－1)) 可以看成p(θ|X(N))的先验概率。最后，当XN给出后，得

p(θ|X(N))=p(XN|θ)p(θ|X(N－1))∫Θp(XN|θ)p(θ|X(N－1))dθ

(5) p(X|ωi)直接由pθ|X(N)计算得到，此时p(X|ωi)可以写为pX|X(N)，由一般概率公式得

pX|X(N)=∫p(X,θ|X(N)) dθ=∫p(X|θ)p(θ|X(N))dθ(325)

这就是贝叶斯学习。下面通过两个例子，讨论正态分布密度函数的贝叶斯估计和贝叶斯学习问题。

【例3.2】对一个单变量正态分布，已知方差σ2，试用贝叶斯估计方法估计均值μ。

【解】设X(N)={X1,X2,…,XN}是ωi类的N个独立抽取的样本，pX|μ~N(μ,σ2)，μ为未知随机参数。假定μ服从正态分布，且μ的先验概率密度p(μ)也为正态分布，即p(μ)~Nμ0,σ20。利用贝叶斯公式求μ的后验概率密度函数p(μ|X(N))，有

p(μ|X(N))=p(X(N)|μ)p(μ)∫p(X(N)|μ)p(μ)dμ

式中，μ的似然函数p(X(N)|μ)可以表示为p(X(N)|μ)=∏Nk=1p(Xk|μ)，且有p(μ|X(N))=α∏Nk=1p(Xk|μ)p(μ)，α=1/∫p(X(N)|μ)p(μ)dμ，α是与μ无关的比例因子，不影响p(μ|X(N))的形式。

因为

pX|μ~Nμ,σ2pμ~Nμ0,σ20

所以


p(μ|X(N))=α∏Nk=1p(Xk|μ)p(μ)
=α∏Nk=112πσexp－Xk－μ22σ2·12πσ0exp－μ－μ022σ20
=α′exp－12∑Nk=1μ－Xk2σ2+μ－μ02σ20
=α″exp－12Nσ2+1σ20μ2－21σ2∑Nk=1Xk+μ0σ20μ(326)


式中，α′和α″是与μ无关的项。

将p(μ|X(N))写为正态分布密度函数的标准形式N(μN,σ2N)： 

p(μ|X(N))=12πσNexp－12μ－μNσN2(327)

比较式(326)与式(327)，可求得μN和σ2N分别为


μN=Nσ20Nσ20+σ2mN+σ2Nσ20+σ2μ0(328)
σ2N=σ20σ2Nσ20+σ2(329)


式中，mN=1N∑Nk=1Xk。

将所求的μN和σ2N代入式(327)就得到了μ的后验概率密度函数p(μ|X(N))。这时，由式(316)计算μ的贝叶斯估计为

μ^=∫μp(μ|XN)dμ=∫μ12πσNexp－12μ－μNσN2dμ=μN

将式(328)结果代入上式，得

μ^=Nσ20Nσ20+σ2mN+σ2Nσ20+σ2μ0(330)

当Nμ0,σ20=N(0,1)且σ2=1时，有

μ^=NN+1mN=1N+1∑Nk=1Xk(331)

也就是说，此时μ的贝叶斯估计与最大似然估计有类似的形式，只是分母不同。

【例3.3】同例3.2，用贝叶斯学习方法计算。

【解】递推求解后验概率密度p(μ|X(N))： 

p(μ|X(N))=12πσNexp－12μ－μNσN2

式中，μN为观察了N个样本后对μ的最好估计，σ2N为估计的不确定性，分别由式(328)和式(329)求得。

由后验概率密度p(μ|X(N))计算类概率密度函数pX|X(N)： 


pX|X(N)=∫p(X|μ)p(μ|X(N))dμ
=∫12πσexp－X－μ22σ2·12πσNexp－μ－μN22σ2Ndμ
=∫12πσσNexp

－(X－μN)22(σ2+σ2N)·12πσN
∫exp

－σ2+σ2N2σ2σ2N
μ－σ2NX+σ2μN
σ2N+σ2dμ
=12πσ2+σ2Nexp－X－μN22σ2+σ2N(332)



可见pX|X(N)是正态分布，均值为μN，方差为(σ2+σ2N)。均值与贝叶斯估计的结果是相同的，原方差σ2增加到(σ2+σ2N)，这是由于用μ的估计值代替了真实值，引起了不确定性的增加。

对于多维正态分布，可以采用与一维情况类似的方法估计均值向量，但计算比较复杂。

3.4非参数估计

以上内容讨论了最大似然估计、贝叶斯估计和贝叶斯学习这三种参数估计方法，其共同的特点是样本概率密度函数的分布的形式已知，而表征函数的参数未知，所需要做的工作是从样本估计出参数的最优取值。但在实际应用中，上述条件往往并不能得到满足，人们并不知道概率密度函数的分布形式，或者函数分布并不典型，或者不能写出某些参数的函数。为了设计贝叶斯分类器，仍然需要获取概率密度函数的分布知识，所以非常有必要研究如何从样本出发，直接推断其概率密度函数。于是，人们提出一些直接用样本来估计总体分布的方法，称为估计分布的非参数法。

非参数估计方法的任务是从样本集X(N)={X1,X2,…,XN}中估计样本空间Ω中任何一点的概率密度p(X)。如果样本集来自某个确定类别(如ωi类)，则估计的结果为该类的类条件概率密度p(X|ωi)。如果样本集来自多个类别，且不能分清哪个样本来自哪个类别，则估计结果为混合概率密度。

3.4.1非参数估计的基本方法

下面从一个例子说明非参数估计的基本思想。假如样本集X(N)={X1,X2,…,XN}由N个一维样本组成，每个样本Xi在以Xi为中心，宽度为h的范围内，对分布的贡献为a。显然，可以把每个样本在Xi点的“贡献”相加作为这点的概率密度p(Xi)的估计。对所有的样本X都这么做，就可以得到总体分布p(X)的估计值。通常，采用某种函数表示某一样本对某点概率密度的贡献，则某点概率密度p(X)的估计为所有样本所做贡献的线性组合。非参数估计的原理如图32所示。



图32非参数估计的原理


当然，也可以认为每个样本对自己所在位置的分布“贡献”最大，距离越远，贡献越小。下面讨论如何估计概率密度函数。

设一个随机向量X落入特征空间区域R的概率P为

P=∫Rp(X)dX
(333)

式中，p(X)是X的概率密度函数，P是概率密度函数的一种平均形式，对P做估计就是估计出p(X)的这个平均值。

设N个样本X(N)={X1,X2,…,XN}，它们是从概率密度为p(X)的总体分布中独立抽取的，则N个样本中有k个样本落入区域R的概率最大，则可以得到


kNP^(334)
P^k/N(335)


式中，P^希望是X落入区域R中的概率P的一个理想的估计，但我们要估计的是类概率密度函数p(X)的估计p^(X)。为此，设p(X)连续，且区域R足够小，以致p(X)在这样小的区域中没有什么变化，由此可得

P=∫Rp(X)dX≈p(X)V(336)

式中，X是R中的一个点，V是R的“体积”。

由式(335)和式(336)可得

kNP^=∫Rp^(X)dX=p^(X)V(337)

所以

p^(X)=k/NV(338)

式（338）就是X点概率密度的估计，它与样本数N、包含X的区域R的体积V和落入R中的样本数有关。

从理论上来讲，如果使p^(X)趋近p(X)，就必须让体积V趋近于零，同时k、N趋向于无穷大。但事实上，V不可能无穷小，样本也总是有限的，N不可能无穷大，所以p^(X)总是存在误差。如果碰巧有一个或几个样本重合于X出现在R，则会使估计发散，甚至到无穷大。因此要是采用这种估计，我们在使用式(338)时就必须注意V,k,k/N随N变化的趋势，使得当N适当增大时能保持式(338)的合理性。

从理论上考虑，假设有无穷多个样本，可以采取如下措施去提高X处的概率密度p(X)的估计精度。构造一个区域序列R1,R2,…，对R1采用一个样本进行估计，对R2采用两个样本进行估计，以此类推，对RN采用N个样本进行估计。设VN是区域RN的体积，kN是落入RN的样本个数。第N次估计的总体概率密度为

p^N(X)=kN/NVN(339)

为了保证上述估计的合理性，应满足以下三个条件： 

(1) limN→∞VN=0(340)
(2) limN→∞kN=∞(341)
(3) limN→∞kNVN=0(342)


此时，总体概率密度p^N(X)的估计值收敛于实际值p(X)。

在上述条件中，条件(1)保证了空间平均式的收敛性； 条件(2)保证了频数比的收敛性； 条件(3)保证了估计式的收敛性。以上三个条件说明当N增大时，落入RN的样本数也增加； VN不断减少，以使p^N(X)趋于p(X)； 尽管在一个小区域RN中落入了大量的样本，但它的数目与样本总数相比还是可以忽略的。满足上述三个条件的区域序列一般有两种选择方法，从而得到两种非参数估计法。

(1) Parzen窗函数法： 选定一个中心在X处的区域RN，其体积VN以N的某个函数(如VN=1/N)的关系不断缩小，同时需对kN和kN/N加以限制，以使p^N(X)收敛于p(X)，然后计算落入RN的样本数kN，用来估计局部密度p^N(X)的值。


(2) kN近邻法： 令kN为N的某个函数(例如，kN=N)，以X为中心构造一个体积为VN的区域RN，使RN恰好包含kN个样本，用这时的体积来估计p^N(X)的值。


3.4.2Parzen窗法

假设区域RN为d维超立方体，向量X为d维特征空间中的一个点，超立方体RN以原点为中心，侧棱长为hN,则其体积VN为

VN=hdN(343)

为了计算RN中包含的样本数kN,定义d维空间的基本窗函数：

φ(u)=1,uj≤12,j=1,2,…,d0,其他(344)

式中，u=(u1,u2,u3,…,ud)，φ(u)称为Parzen窗函数，它是以原点为中心的超立方体。

利用函数φ(u)可以实现对落在区域RN的样本进行计数，当Xi落在以X为中心、体积为VN的超立方体内时，计数为1，即

φ(u)=φX－XihN=1(345)

否则取值为零。因此，落入该立方体的样本数kN为

kN=∑Ni=1φX－XihN(346)

将式(346)代入式(339)，可得概率密度估计： 

p^N(X)=kN/NVN=1N∑Ni=1φX－XihNVN(347)

式(347)是Parzen窗法估计的基本公式。上式表明，Parzen窗的估计函数实质上由一系列的基函数叠加而成，式(344)定义的窗函数即叠加的基函数，每个样本点处作为叠加节点，使用kN个以样本Xi为中心的窗函数叠加作为X处的概率密度函数的估计，每一个样本对估计所起的作用依赖于它到Xi的距离。显然，概率密度函数与样本的密集程度相关，样本在该区域越密集，叠加函数的值越大，也意味着概率密度函数的估计值越大。


式(344)定义了Parzen窗法估计法的窗函数，从上面的讨论可以看出估计结果和窗函数密切相关。下面讨论窗函数需要满足什么条件，以及如何去选择合适的窗函数。为了使p^N(X)成为一个概率密度函数，则其必须满足概率密度函数的一般要求，即p^N(X)非负且积分为1，相应地要求窗函数φ(u)满足下面两个条件： 


φ(u)≥0(348)
∫φ(u)du=1(349)



上述两式表明，窗函数本身满足密度函数的要求。式(347)窗函数φ(u)的非负性保证了p^N(X)的非负性，进一步有


∫pN(X)dX=
∫1N∑Ni=11VNφ
X－XihNdX
=1N∑Ni=1∫1VNφ
X－XihNdX
=1N∑Ni=1∫φ(u)du=1


从而证明了p^N(X)是一个概率密度函数。

由此可见，一个函数只要满足式(348)和式(349)，它就可以作为窗函数。除了上面选择的超立方体窗函数以外，还可以选择其他的窗函数形式。以一维窗函数为例，常用的方窗函数、正态窗函数和指数窗函数的定义如下。

(1) 方窗函数： 

φ(u)=1,|u|≤12
0,其他(350)

(2) 正态窗函数： 

φ(u)=12πexp-12u2(351)

(3) 指数窗函数： 

φ(u)=12exp－|u|(352)

以上三种窗函数如图33所示。



图33三种窗函数


对密度函数的估计有影响的另一个因素是窗函数的宽度hN，下面分析hN对估计量的影响。定义函数

δN(X)=1VNφXhN(353)

由式(347)可以表示为

p^N(X)=1N∑Ni=1φX－XihNVN=1N∑Ni=1δN(X－Xi)VN(354)

由VN=hdN 可知，hN既影响δN(X)的幅度又影响它的宽度。如果hN较大，则δN(X)的幅度就较小，从而δN(X)宽度较大，只有当Xi距离X较远时才使得δN(X－Xi)与δN(0)相差较大。此时，p^N(X)就变成N个宽度较大的缓慢变化函数的叠加，造成估计值p^N(X)较平滑，跟不上函数p(X)的变化，估计分辨率较低。反过来，如果hN选得较小，则δN(X－Xi)的幅度就较大，δN(X)宽度较小。此时p^N(X)就是N个以Xi为中心的尖脉冲的叠加，p^N(X)波动较大，从而使估计不稳定。

综上所述，hN的选取对p^N(X)的影响很大，hN太大或太小，都对估计的精度不利。如何选择hN，需要一定的经验，当样本数目有限时，可做适当的折中，通过试探的方法选出合理的结果； 当样本数目无限时，可让VN随N的增大而缓慢地趋于零，从而使p^N(X)收敛于p(X)。

【例3.4】设待估计的p(X)是均值为0、方差为1的正态密度函数，用Parzen窗法估计p(X)，即求估计式p^N(X)。

【解】考虑X是一维模式向量的情况。选择正态窗函数： 

φ(u)=12πexp－12u2

并设hN=h1/N，h1为可调节的参数。有

X－XihN=12πexp－12X－XihN2=12πexp－12X－Xih1/N2

这样，估计值p^N(X)是一个以样本为中心的正态密度的平均值： 


p^N(X)=1N∑Ni=11hNX－XihN=1N∑Ni=1Nh1X－XihN
=1h1N12πexp－12X－Xih1/N2




当采集到正态分布样本后，就可求得估计值p^N(X)，结果如图34所示。图34给出了h1分别取0.25、1和4，样本量N取1、16、256和∞时p^N(X)的估计情况。从图中可以看出，样本量越大，估计结果越精确； 同时，当样本量较小时，窗宽的选择对估计结果也有一定的影响。如N=1时，h1=1的估计结果明显好于其他两种情况。


【例3.5】已知一维随机变量x，假设待估计的概率密度函数p(x)为两个均匀分布密度函数的混合，用Parzen窗法估计p(x)，即求估计式p^N(x)。

p(x)=1,－2.5＜x＜－2
0.25,0＜x＜2
0，其他


【解】仍选择正态窗函数，hN的定义同例3.4，估计结果如图35所示。从图中可以得出，当N较小时，估计结果与真实分布相差很大，当N=1时，p^N(x)只是反映出窗函数本身，当h1=0.25，N=16时还能看到单个样本的作用，当h1=1和h1=4时，就显得比较平滑了。当N=16时，还无法确定哪一种估计比较好，但当N=256和h1=1时，估计基本满足要求； 当N增大时，估计结果与真实分布越来越接近。


通过以上的例子可以归纳出非参数估计的优缺点。其优点是该方法具有一般性，对规则或不规则分布、单峰或多峰分布都可以用这个方法得到概率密度估计，而且只要样本足够多，总可以保证收敛于任何复杂的未知概率密度函数。其缺点是要想得到较为满意的估计结果，就需要比参数估计方法所要求的样本数多得多的样本，因此就需要大量的计算时间和存储量。而且随着样本特征维数的增加，用于估计的样本数量也需要相应地增加。



图34正态分布的Parzen窗估计




图35均匀分布的Parzen窗估计




3.4.3kN近邻估计法

Parzen窗估计方法是基于落入体积VN中的样本来做总体估计。由于体积VN中的样本是样本数N的函数，故无论窗函数怎样选取，Parzen窗估计方法都难以是最佳的。此外，体积序列V1,V2,…,VN的选择也是一个难题，由前面的分析，我们知道估计的结果对体积的选择比较敏感。为了解决这一问题提出了kN近邻估计方法(简称kN近邻法)。


kN近邻估计方法的基本思想是，使包含X点的序列体积V1,V2,…,VN受落入RN中的样本数kN控制，而不是作为实验样本数N的函数。具体来说，可以预先确定N的一个函数kN，在X点近邻选择一个体积，并使其不断增大直到捕获kN个样本为止，这kN个样本就是X点的近邻。很显然，如果X点附近的概率密度较大，则包含kN体积较小； 反之，如果X点附近的概率密度较小，则包含kN体积较大。kN近邻估计方法使用的基本估计公式仍为

p^N(X)=kN/NVN

当满足式(340)、式(341)和式(342)三个条件时，p^N(X)收敛于概率密度p(X)。kN近邻估计具有以下特点。

(1) kN大小的选择会影响估计的结果。kN可以选择为样本容量N的某种函数，在kN=k1N,kN≥1条件下，当样本容量N→∞时，可以保证p^N(X)收敛于真实分布p(X)。但是在有限样本容量条件下，k1的选择也会影响估计结果的正确性。

(2) kN近邻估计方法的计算量大。与Parzen窗法一样，为保证估计结果的正确性，所需样本量N一般要很大，尤其当样本的特征维数比较高时更是如此，因此存储量和计算量都很大。经验数据表明，当样本的特征维数为一维时，所需样本容量一般为数百个。但是，当样本特征维数为二维时，样本容量就需要几千个。随着样本特征维数的增加，样本容量会急剧增长，带来超大计算量与超大存储量的问题。目前已有一些针对该问题的解决方法，有兴趣的读者可以参考相关资料。

图36所示为概率密度函数分别为正态分布和双峰均匀分布时用kN近邻法估计p(X)的结果。




图36kN近邻估计p(X)的结果


3.5Python示例

【例3.6】生成100万个服从标准正态分布的随机数，用最大似然估计这些随机数服从的正态分布的参数值。


import numpy as np

import matplotlib.pyplot as plt

from scipy import stats

from scipy.stats import norm



# 生成正态分布样本

testdata = np.random.normal(loc=0, scale=1.0, size=1000000)

print('样本均值的无偏估计：', np.mean(testdata))

print('样本标准差的无偏估计：', np.std(testdata, ddof=1))




运行结果：


样本均值的无偏估计： -3.8854235939791765e-06

样本标准差的无偏估计： 1.0008638816906816




程序生成了100万个服从标准正态分布的随机数，并直接计算了它们的均值和标准差。然后，绘制出这些样本点的直方图（如图37所示），观察其与正态分布函数的逼近程度。



图37样本点的直方图




# 绘制正态分布图

plt.hist(testdata, density=True, bins=30)# 归一化直方图(用出现频率代替次数)，将划分
# 区间变为 20(默认 10)

x = np.linspace(-3, 3, 50)# 在(-3,3)返回均匀间隔的50个数字

plt.plot(x, norm.pdf(x), 'r-')

plt.title('正态分布图')

plt.show()

# 正态分布均值区间估计

def confidence_interval_u(data, sigma=-1, alpha=0.05, side_both=True):

xb = np.mean(data)  # 均值

# s = np.std(data, ddof=1)  # 无偏标准差

# sigma已知，服从标准正态分布

z = stats.norm(loc=0, scale=1.0)

if side_both:  # 求双侧置信区间

tmp = sigma/np.sqrt(len(data))*z.ppf(1-alpha/2)

else:  # 单侧置信下限或单侧置信上限

tmp = sigma/np.sqrt(len(data))*z.ppf(1-alpha)

bottom_limit = xb-tmp

top_limit = xb+tmp

return bottom_limit, top_limit



# 正态分布标准差区间估计

def confidence_interval_sigma(data, mu=-1, alpha=0.05, side_both=True):

sum_tmp = 0.0

for i in data:

sum_tmp = sum_tmp + (i - mu) ** 2

if side_both:  # 双侧置信区间

bottom_limit = sum_tmp / stats.chi2.ppf(1 - alpha / 2, df=len(data))

top_limit = sum_tmp / stats.chi2.ppf(alpha / 2, df=len(data))

else:  # 单侧置信下限或单侧置信上限

bottom_limit = sum_tmp / stats.chi2.ppf(1 - alpha, df=len(data))

top_limit = sum_tmp / stats.chi2.ppf(alpha, df=len(data))

return np.sqrt(bottom_limit), np.sqrt(top_limit)



paramhat = confidence_interval_u(testdata, 1)





paramint = confidence_interval_sigma(testdata, 0)

print('样本均值的区间估计:', np.round(paramhat, 6))

print('样本标准差的区间估计:', np.round(paramint, 6))




运行结果：


样本均值的区间估计: ［-0.001964  0.001956］

样本标准差的区间估计: ［0.999478 1.002252］




程序可以获得两个结果： paramhat和paramint。它们分别对应均值和标准差的区间估计。从结果来看，均值和标准差的区间估计都很窄，而且很接近生成随机数时使用的参数值。

【例3.7】设置两类为2×60的样本，第一类样本是一些随机设定的0.6~1.2的数，第二类样本是第一类样本的第一行加0.1，第二行减0.1而生成的。用贝叶斯方法分别估计这两类样本的条件概率密度的均值参数u。其中，设两类样本的先验估计都为［0.93,0.94］，p(u)的方差为［0.05,0.05］。


import numpy as np

def bayes_estimation(class_id, X, sigma0, mu0):

# 求样本均值

m1, m2 = np.mean(X, axis=1)

print(f'第{class_id}类样本均值：', ［m1, m2］)



# 求解样本方差

s1, s2 = np.var(X, axis=1, ddof=1)

print(f'第{class_id}类样本方差：', ［s1, s2］)



# 利用贝叶斯公式(3-30)

N = 60

mu_1 = N * sigma0［0］ * m1 / (N * sigma0［0］ + s1) + s1 * mu0［0］ / (N * sigma0［0］ + s1)

mu_2 = N * sigma0［1］ * m2 / (N * sigma0［1］ + s2) + s2 * mu0［1］ / (N * sigma0［1］ + s2)

print(f'第{class_id}类样本均值估计值：', ［mu_1, mu_2］)



# 设置方差

c = ［0.05, 0.05］

# 设置先验估计参数

d = ［0.93, 0.94］



# 生成第一类样本

# 生成范围在(0.6,1.2)的2*60的样本，并保留一位小数

X = np.random.uniform(0.6, 1.2, size=(2, 60)).round(1)

print('第一类样本数据：', X)

bayes_estimation(1, X, c, d)


运行结果：

第一类样本数据： ［［0.8 0.7 1.  1.  1.  0.6 1.1 1.2 1.1 0.7 0.8 0.8 0.6 0.6 1.  1.  0.7 0.7 0.8 0.7 1.  0.7 1.1 1.  1.  0.7 1.  0.9 0.8 1.  0.8 1.  0.7 1.  1.1 1. 0.8 0.7 1.1 1.2 0.8 0.9 0.9 0.9 1.1 0.9 0.8 0.8 0.6 1.2 0.9 1.  0.8 0.6 1.1 0.7 1.  0.7 1.1 0.8］

［1.1 1.1 0.7 1.1 0.9 0.9 0.6 0.8 1.1 0.9 1.  0.7 0.9 0.6 0.9 1.  1.1 0.8 1.  1.2 0.8 0.6 1.2 1.  0.6 1.1 1.  0.6 1.  1.  0.9 1.  0.8 1.1 0.8 0.9 1.1 0.8 0.9 1.1 0.7 1.  1.  1.1 0.6 1.  0.6 1.1 1.  1.1 0.8 0.8 0.7 1.2 1.  1.1 0.9 0.7 1.2 0.8］］






第1类样本均值： ［0.885, 0.9183333333333333］

第1类样本方差： ［0.029432203389830512, 0.03237005649717515］

第1类样本均值估计值： ［0.8854371938579977, 0.918564621471337］



# 生成第二类样本

x1 = X［0, :］ + 0.1  # 生成第一行数据

x2 = X［1, :］ - 0.1  # 生成第二行数据

X2 = np.vstack((x1, x2))  # 将两行数据进行垂直拼接

# print('第二类样本数据：', X2)

bayes_estimation(2, X2, c, d)



运行结果：

第2类样本均值： ［0.985, 0.8183333333333335］

第2类样本方差： ［0.029432203389830533, 0.032370056497175136］

第2类样本均值估计值： ［0.9844656519513362, 0.8196321051852005］




【例3.8】产生1、16、256和16384个服从一维标准正态分布的样本。

(1) 用窗宽h1分别为0.25、1、4，hN=h1/N，窗函数为高斯函数的情形估计所给样本的密度函数并画出图形。

(2) 当kN=N时，用kN近邻法估计所给样本的密度函数并画出图形。

【解】Python代码及运行结果如下。

(1) 样本量N分别取1、16、256、16384，h1分别取0.25、1、4时采用Parzen窗法的估计，结果如图38所示。



图38用Parzen窗法估计的结果




import numpy as np

import math

import matplotlib.pyplot as plt



def Paezen_window(u):# 正态分布窗函数

s = np.exp(-(u * u) / 2.0) / math.sqrt(2 * math.pi)

return s



def Parzen(N, x, h1):

n = len(N)

hn = h1 / (math.sqrt(n))

# hn = h1/(math.log(n) + 1)

Px = ［］

for i in x:

p = Paezen_window((N - i) / hn)

Px.append(np.sum(p) / (h1 * math.sqrt(n)))# 概率密度公式

return np.array(Px)



x = np.arange(-5, 5, 0.1).reshape(［-1, 1］)

N_number = ［1, 16, 256, 16384］

h1 = ［0.25, 1, 4］

result = ［］



for i in range(4):

N_i = np.random.normal(0, 1, N_number［i］).reshape(［-1, 1］)

for k in range(0, 3):

result_i = Parzen(N_i, x, h1［k］)

result.append(result_i)



plt.figure(figsize=(8, 6), dpi=300)

for i in range(1, 13):

plt.subplot(4, 3, i)

plt.plot(x, result［i-1］)

if i % 3 == 1:

plt.ylabel('N={0}'.format(N_number［(i-1)//3］))

if i <= 3:

plt.title('h1={}'.format(h1［i-1］))

plt.tight_layout()

plt.show()




(2) 样本量N分别取1、16、256、16384时采用kN近邻法估计，结果如图39所示。




from scipy.spatial.distance import cdist

import numpy as np

import math

import matplotlib.pyplot as plt



def knn_estimate(n, x, d):

# n 随机样本, x 横坐标，d 特征维度

Kn = int(math.sqrt(len(n)))

Px = ［］

q = Kn/len(n)

for i in x:

hn = cdist(np.array(［i］), n, metric='euclidean').reshape(-1)





hn.sort()

if Kn == 1:

Vn = (hn［Kn-1］ * 2) ** d

else:

Vn = (hn［Kn］ * 2) ** d

Px.append(q/Vn)# 计算公式

return np.array(Px)



X = np.arange(-4, 4, 0.005).reshape(［-1, 1］)# reshape(［-1,1］)将矩阵转换为一列



result = ［］

N_number = ［1, 16, 256, 16384］



for i in range(4):

N_i = np.random.normal(0, 1, N_number［i］).reshape(［-1, 1］)

result_i = knn_estimate(N_i, X, d=1)

result.append(result_i)



fig, axes = plt.subplots(2, 2, figsize=(8, 6), dpi=200)

ax = axes.ravel()

for i in range(0, 4):

ax［i］.plot(X, result［i］)

ax［i］.set_title('N={0}＼nKn={1}'.format(

N_number［i］, math.sqrt(N_number［i］)))

plt.tight_layout()

plt.show()






图39用kN近邻法估计的结果


在实际应用中，类条件概率密度通常是未知的。那么，在先验概率和类条件概率密度都未知或者其中之一未知的情况下，该如何来进行类别判断呢？其实，只要能收集到一定数量的样本，根据统计学的知识，我们是可以从样本集来推断总体概率分布的。


监督参数估计就是由已知类别的样本集对总体分布的某些参数进行统计推断； 非监督参数估计是已知总体概率密度函数形式但未知样本所属的类别，要求推断出概率密度函数的某些参数。通常采用最大似然估计方法和贝叶斯估计方法。最大似然估计把参数看成确定(非随机)而未知的，最好的估计值是在获得实际观察样本的概率为最大的条件下得到的； 而贝叶斯估计则是把参数当成具有某种分布的随机变量，样本的观察结果使先验分布转换为后验分布，再根据后验分布修正原先对参数的估计。非参数估计是已知样本所属的类别，但未知总体概率密度函数的形式，要求我们直接推断概率密度函数本身。统计学中常见的一些典型分布形式不总是能够拟合实际中的分布。此外，在许多实际问题中经常遇到多峰分布的情况，这就迫使我们必须用样本来推断总体分布，常见的总体类条件概率密度估计方法有Parzen窗法和kN近邻法两种。


习题及思考题


3.1证明： 按照贝叶斯决策理论进行分类时，其结果满足分类错误率最小。

3.2一元正态分布的最大似然估计： 假设样本X服从正态分布N(μ,σ2),已获得一组样本X1,X2,…,XN，用Python语言设计一程序片段，计算估计参数μ,σ2。

3.3简述参数估计、非参数估计和非参数分类器等概念间的关系。

3.4证明： 对正态总体的期望μ的最大似然估计是无偏的，对方差σ2的最大似然估计是有偏的。

3.5设样本X1,X2,…,XN是来自p(X|θ)的随机样本，其中0≤x≤θ时，p(X|θ)=1θ，否则为0。证明θ的最大似然估计是maxXk。

3.6设p(X)~N(μ,σ2)，窗函数φ(X)~N(0,1)，指出Parzen窗估计p^N(X)=1NhN ∑Ni=1φX－XihN，对于小的hN有如下性质。

(1) Ep^N(X)~Nμ,σ2+h2N。

(2) varp^N(X)=1NhN2πp(X)。

3.7设总体X的概率密度函数为f(X,θ)=(θα)Xα－1e－θxa，求参数θ的最大似然估计。

3.8在掷硬币的游戏实验中，正面出现的概率是q，反面出现的概率是1－q。设Xi，i=1,2,…,N是这个实验的结果，Xi∈(0,1)。

(1) 证明q的最大似然估计是qML=1N∑Ni=1Xi。

提示： 似然函数是p(X,q)=∏Ni=1qXi(1－q)(1－Xi)。

(2) 证明最大似然估计结果是下列方程的解： 

qΣiXi(1－q)(N－ΣiXi)ΣiXiq－N－ΣiXi1－q=0

3.9证明： 对于对数正态分布p(X)=1σX2πexp－lnx－θ22σ2,X>0,
最大似然估计θ^ML=1N∑Nk=1lnXk。