《武汉工程大学学报》  2010年07期 99-106   出版日期：2010-07-31   ISSN:1674-2869   CN:42-1779/TQ

---

非参数估计方法


0引言函数估计［1］是一个经典反问题，一般定义为给定输入输出样本对，求未知的系统函数［2］.传统的方法为参数方法，即构建一个参数模型，再定义某个误差项，通过最小化误差项来求解模型的参数［3］.参数方法尽管较为简单，但不够灵活.例如参数模型假设有误，则会导致整个求解流程失败［4］.因此学者们发展出不少新技术，非参数估计就是其中一项较好的方法.该方法无需提前假设参数模型的形式，而是基于数据结构推测回归曲面［5］.本文首先研究了经典的2种参数回归方法：最小二乘法与内插函数法，分析了它们的不足，然后主要讨论8种非参数回归方法：核方法、局部多项式回归、正则化方法(样条估计)、正态均值模型、小波方法、过完全字典、前向神经网络、径向基函数网络，尤其详细介绍了其间的相关性与继承性.最后，研究了高维情况下面临的计算维数诅咒与样本维数诅咒.1回归模型考虑模型yi=r(xi)+εi(1)
式（1）中(xi,yi)为观测样本，假定误差ε具有方差齐性，则r=E(y|x)称为y对x的回归函数，简称回归.一般地，可以假设x取值在［0,1］区间内.定义“规则设计”为xi=i/n(i=1,2,…, n).并定义风险函数为
R=∑ni=1［r(xi)-(xi)］2=∑ni=1［yi-(xi)］2(2)
式（2）中为系统函数r的估计.回归一词源于高尔顿(Galton)，他和学生皮尔逊(Pearson)在研究父母身高和子女身高的关系时，以每对夫妇的平均身高为x，取其一个成年儿子的身高为y，并用直线y=33.73+0.512x来描述y与x的关系.研究发现：如果双亲属于高个，则子女比他们还高的概率较小；反之，若双亲较矮，则子女以较大概率比双亲高.所以，个子偏高或偏矮的夫妇，其子女的身高有“向中心回归”的现象，因此高尔顿称描述子女与双亲身高关系的直线为“回归直线”［6］.然而，并非所有的x-y函数均有回归性，但历史沿用了这个术语.更为精确的表达是“函数估计”.2传统方法理论上描述一个函数需要无穷维数据，因此函数估计本身也可称为“无穷维估计”［7］.传统的估计方法有下列两种极端情形.2.1最小二乘法此时假设(x)=β0+β1x，采用最小二乘法计算权值β=(β0,β1)，得到的解为最小二乘估计［8］，(x)=(XTX)-1XTY(3)
则对给定样本点的估计r=［(x1),(x2),…,(xn)］T可写为r=X=LY(4)
这里Y=(y1,y2,…,yn)T.L=X(XTX)-1XT称为帽子矩阵［9］.以5个样本点的一维规则设计矩阵为例，此时
X=0.2
0.4
0.6
0.8
1.0L=0.0180.0360.0540.0720.090
0.0360.0720.1090.1450.181
0.0540.1090.1630.2180.272
0.0720.1450.2180.2900.363
0.0900.1810.2720.3630.454(5)L满足L=LT,L2=L.另外，L的迹等于输入数据的维数p，即trace(L)=p.这里输入数据是一维的，所以trace(L)=1.2.2内插函数法此时对(x)不加任何限制，得到的是该数据的一个内插函数［10］.同样以5个样本点的一维规则设计矩阵为例，由于样本点的估计r=［(x1),(x2),…,(xn)］T完全等于(y1,y2,…,yn)T，所以帽子矩阵为L=10000
01000
00100
00010
00001(6)2.3两种方法的缺陷图1给出了这两种极端拟合的示意图，数据是被高斯噪声干扰的正弦函数，采用上述两种方法拟合，结果表明：最小二乘法过光滑，未展现数据内部的关系；而内插函数法忽略了噪声影响，显得欠光滑.从帽子矩阵也可看出，式(5)表明最小二乘法对每个数据的估计都利用了所有样本，这显然导致过光滑，且x值越大的数据权重越大，这明显与经验不符；反之，式(6)表明内插函数法仅仅利用了最邻近的样本数据，这显然导致欠光滑.第7期张煜东,等：非参数估计方法
武汉工程大学学报第32卷
图1两种极端拟合
Fig.1Two extreme fitting2.4非参数回归的优势非参数回归(nonparametric regression)作为最近兴起的一种函数估计方法，是一种分布无关(distribution free)的方法，即不依赖于数据的任何先验假设.与此对应的是参数回归(parametric regression)，通常需要预先设置一个模型，然后求取该模型的参数.非参方法的本质在于：模型不是通过先验知识而是通过数据决定.需要注意的是，“非参数”并不表示没有参数，只是表示参数的数目、特征是可变的(flexible).由于非参方法无需数据先验知识，其应用范围较参数方法更广，且性能更稳健.其另一个优点是使用过程较参数方法更为简单.然而，它也存在缺点，一般结构更复杂，需要更多的运算时间.2.5线性光滑器需要说明的是，最小二乘法、内插函数法、核方法、正则化方法、正态均值模型均是线性光滑器.定义为：若对每个x，存在向量l(x)=［l1(x),…,ln(x)］T，使得r(x)的估计可写为(x)=∑ni=1li(x)yi(7)则估计为一个线性光滑器［11］.显然权重li(x)随着x而变化，这与信号处理中的“自适应滤波器”非常相似.3核回归核方法［12］定义为(x)=∑ni=1li(x)Yi(8)权重li由式（9）给出li=Kx-xih∑ni=1Kx-xih(9)这里h是带宽，K是一个核，满足K(x)≥0，以及
∫K(x)dx=1，∫xK(x)dx=0，∫x2K(x)dx＞0，(10)常用的核函数见表1.表1常用的核公式
Table 1Frequentlyused kernel formula
核公式boxcarK(x)=0.5*I(x)GaussianK(x)=12πexp-x22EpanechnikovK(x)=34(1-x2)I(x)TricubeK(x)=7081(1-|x|3)3I(x)以boxcar核为例，帽子矩阵为L=1/21/2000
1/31/31/300
01/31/31/30
001/31/31/3
0001/21/2(11)显然，这可视作最小二乘法与内插函数法的折中.为了估计带宽h，首先必须估计风险函数，一般可采用缺一交叉验证得分CV=(h)=1n∑ni=1［yi--i(xi)］2(12)这里-i(xi)为未用第i个数据所得到的估计，使CV最小的h，即为最佳带宽.为了加速运算，可将式（12）重新写为(h)=1n∑ni=1yi-(xi)1-Lii2(13)这里Lii是光滑矩阵L的第i个对角线元素.另一种方法是采用广义交叉验证法，规定
GCV(h)=(h)=1n∑ni=1yi-(xi)1-v/n2(14)这里v=tr(L).4局部多项式回归采用核回归常会碰到下列2个问题［13］：1）若x不是规则设计的，则风险会增大，称为设计偏倚(design bias)；2）核估计在接近边界处会出现较大偏差，称为边界偏倚(boundary bias).为了解决这2个问题，可采用局部多项式回归.局部多项式回归［14］可视作核估计的一个推广，首先定义权函数ωi(x)=K［(xi-x)/h］，选择a=(x)来使得下面的加权平方和最小∑ni=1ωi(x)(yi-a)2(15)利用高等数学知识，可以看出解为(x)=∑ni=1ωi(x)yi∑ni=1ωi(x)(16)可见式(16)正好是核回归估计.这表明核估计是由局部加权最小二乘得到的局部常数估计.因此，若利用一个p阶的局部多项式而不是一个局部常数，就可能改进估计，使曲线更光滑.定义多项式
Px(u;a)=a0+a1(u-x)+a22!(u-x)2+
…+app!(u-x)p(17)则局部多项式的思想是：选择使下列局部加权平方和∑ni=1ωi(x)［yi-Px(xi;a)］2(18)最小的a，估计=(0,1,p)T依赖于目标值x，最终有(x)=Px(x;)=0(x)(19)当p等于0时，等于核估计；当p=1时，称为局部线性回归(local linear regression)估计［15］，由于其算法简单且性能优越，较为常用.5基于正则化的回归为了描述方便，这里假设数据点为［(x0,y0),(x1,y1),…(xn-1,yn-1)］.在风险函数(2)后增加一项惩罚项，一般设为r(x)的二阶导数J=∑n-1i=0yi-(xi)2+λ∫［r″(x)］2dx(20)λ控制了解的光滑程度：当λ=0时，解为内插函数；当λ→∞时，解为最小二乘直线；当0<λ<∞时，(x)是一个自然三次样条.需要注意下列事项：首先三次样条表示曲线在结点(knot)之间是三次多项式，且在结点处有连续的一阶和二阶导数；其次一个m阶样条为一个逐段m-1阶多项式，所以三次样条是4阶的(m=4)；第三，自然样条表示在边界点处二阶导数为0，即在边界点外是线性的；第四，样条的结点等于数据点.为了加速计算，将数据点重新排序，假设a,b为样本点x的上下界，令a=t1≤t2≤…≤tn-1=b，这里t是x重新排序后的点，称为结点.可用B样条基(Bspline basis)［16］作为该三次样条的基，即
(x)=∑n-m-1i=0Pibi,m(t) t∈［tm-1,tn-m］(21)Pi称为控制点，共n-m个，形成一个凸壳.n-m个B样条基可通过如下计算，首先初始化：bj,0(t)=1if tj＜t＜tj+1
0otherwise(22)然后对i=1，逐步+1，直到i=m-1，重复迭代下式：
bj,i(t)=t-tjtj+m-1-tjbj,i-1(t)+
tj+m-ttj+m-tj+1bj,i+1(t)(23)若结点等距，则称B样条是均匀的(uniform)，否则称为不均匀.如果两个结点相等，计算过程会出现0/0情况，此时默认结果为0.令矩阵B的第(i, j)元素bij=bj(xi)，矩阵Ω的第(i, j)元素Ωij=∫bni(x)bnj(x)dx，则控制点可由式（24）求得P=(BTB+λΩ)-1BTY(24)可见，样条也是一个线性光滑器.表面上看，基于核的估计与基于正则化的估计原理与模型均不一致，但是Silverman证明了如下定理，样条估计(x)可视作如下所示的一种渐近的核估计li(x)≈1f(xi)h(xi)Kxi-xh(xi)(25)式中，f(x)是x的密度函数.h(x)=λnf(x)1/4(26)K(x)=12exp-|x|2sin|x|2+π4(27)显然，若样本x是规则设计，则f(x)=1, h(x)=(λ/n)1/4=h,li(x)∝K［(xi-x)/h］，即此时样条估计可视作形如式(27)的渐近核估计.6正态均值模型令φ1,φ2,…为一个标准正交基，则显然r(x)可以展开为r(x)=∑∞i=1θiφi，定义Zj=1n∑ni=1yiφj(xi)(28)则随机变量Zj是正态分布，且均值与方差满足：E(Zj)=θjV(Zj)=σ2/n(29)可见，若估计出θ，则可近似求得(x)≈∑ni=1θiφi.因此正态均值模型将n个样本的函数估计问题转换为估计n个正态随机变量Zj的均值θ的问题［17］.若直接令=Z，则显然得到一个很差的估计，下面给出风险更小的估计.首先，必须做出一个关于的风险估计，Stein给出下列定理：令Z～N(θ,V),=(Z)为θ的一个估计，并令g(Z1,…,Zn)=-Z.则的风险的一个无偏估计为(z)=tr(V)+2tr(VD)+∑ig2i(z)(30)式中gi=i-zi，且D的第(i, j)个元素为g(z1,…,zn)的第i个元素关于zj的偏导数［18］.假设=bZ=(b1Z1,…,bnZn)，式中b称为调节器，根据b的设置，存在下列3种情况：①b=(b,b,…,b)，称为常数调节器(constant modulator)，此时令式(30)最小的 称为JamesStein估计；②b=(1,…,1,0,…,0)，称为嵌套子集选择调节器(nested subset selection modulator)，此时令式(30)最小的称为REACT方法.需要注意的是，若基选择傅立叶基，则该方法类似于频域低通滤波器方法.③b=(b1,b2,…,bn)满足1≥b1≥b2≥…≥bn≥0，称为单调调节器(monotone modulator)，该方法理论最优，但是需要的运算量太大，几乎不实用.7小波方法小波方法［19］适用于空间非齐次(spatially inhomogeneous)函数，即函数的光滑程度随着x会有本质性的变化.它可视作正态均值模型的推广，但存在两点区别：一是采用小波基代替传统的正交基，因为小波基较一般的正交基具有局部化的优点，能实现多分辨率分析；另一点是采用了一种称为“阈”的收缩方式.不妨假定父小波为φ，母小波为ψ，同时规定下标(j, k)的意义如下：fj,k(x)=2j/2f(2jx-k)(31)为了估计函数r，用n=2J项展开来近似r，r(x)≈∑2J0-1k=0αj0,k(x)φj0,k(x)+
∑Jj=J0∑2j-1k=0βj,kψj,k(x)(32)这里J0是任取常数，满足0≤J0≤J.α称为刻度系数，β称为细节系数.那么如何估计这些系数？首先计算Sk=1n∑iφj0,k(xi)yi(33)Djk=1n∑iψj,k(xi)yi(34)Sk、Djk分别称为经验刻度系数与经验细节系数，可知Sk≈N(αj0,k,σ2/n)，Djk≈N(βj,k,σ2/n)，可估计方差为=nmedian(|Dj-1,k-median(Dj-1,k)
|∶k=0,…,2J-1-1)0.6745(35)然后根据Sk、Djk、可得α与β的估计如下：j0,k=Sk(36)β的估计形式稍许复杂，采用硬阈与软阈的方式分别为jk=0|Djk|＜λ
Djk|Djk|≥λ(37)jk=sign(Djk)|Djk|-λ+(38)之所以采用阈的形式，是因为稀疏性(sparse)的思想［20］：对某些复杂函数，在小波基上展开时系数也是稀疏的.因此，需要采用一种方式来捕获稀疏性.然而，传统的L2范数不能捕捉稀疏性，相反，L1范数与非零基数能够较好地捕捉稀疏性.例如，考虑n维向量a=(1,0,…,0)与b=(1/n1/2,…,1/n1/2)，有‖a‖2=‖b‖2=1，可见，L2范数无法区分稀疏性.反之，‖a‖1=1，‖b‖1= n1/2，因此，L1范数能提取稀疏性；另外，若令非零基数为J(θ)={#(θi≠0)}，则J(a)=1，J(b)=n，因此，非零基数也能提取稀疏性.最后，在正则化估计中若惩罚项分别为L1范数或非零基数，则最优估计恰好对应着软阈估计与硬阈估计.最后，需要解决阈估计中λ的计算问题，这里介绍两种最简单的方式：一是通用阈值(universal threshold)，即对所有水平的分辨率阈值均一致，
λ=2lognn(39)另一种是分层阈值(levelbylevel threshold)，即对不同分辨率采用不同阈值，一般是通过最小化下式求得
S(λj)=∑njk=12n-22nI|jk≤λj+min(2jk,λ2j)(40)λj∈［0,(/nj)2lognj］
式中nj=2j-1为在水平j的参数个数.8超完备字典小波基较标准正交基的改进在于更加局部化，因此能实现对跳跃的捕捉.然而，虽然小波基非常复杂，但面对各种复杂的函数还是不够灵活.这种缺陷的根源在于：小波基是标准正交基，任意两个基函数之间正交，这保证了基函数简单完整的同时，也丧失了灵活性.基追踪(basis pursuit)方法［21］的思想是采用一种超完备(overcomplete)的基，例如对“光滑加跳跃”的函数，传统的傅立叶基能够捕捉光滑部分，但是难以捕捉跳跃部分；采用小波基能轻易捕捉跳跃部分，但是描述光滑部分较为困难.此时若将“傅立叶基”与“小波基”合并成一个新的基，则显然这种基能够轻松地估计“光滑加跳跃”函数.但是，这种新的基不再正交，它以牺牲正交性来获得更好的灵活性［22］，故此时用“字典”来描述更精确，而本文为了简便统一仍采用“基”表述.9前向神经网络以一个双层神经网络为例，记网络的输入神经元个数为m, 隐层神经元个数为n, 输出层神经元个数为q，则网络结构如图2所示.图2前向神经网络
Fig.2Forward neural network与上面几节线性方法不同的是，神经网络属于非线性统计数据建模(nonlinear statistical data modeling)，其隐层暗含了“特征提取”的思想，且可视作输入数据在一种“自适应的非线性非正交的基”上的映射.同样地，此时基牺牲了正交性、线性、不变性，增加了计算负担，但换来了更加强大的灵活性［23］.简而言之，前向神经网络采用了类似基追踪的方法［24］，但基是自适应变化的、非线性的，因此更加灵活.前向神经网络与基追踪相似之处在于，两者的基都不是正交的，都是根据给定数据而自适应选取的最佳基.前向神经网络的优势在于无不需预选字典，字典在算法中自动生成，并可作为特征选择的一种方法.10径向基函数网络首先观察径向基函数（RBF）神经元如图3所示.图3RBF神经元图
Fig.3Neuron of RBF图中输入向量p的维数为R，首先p与输入层权值矩阵IW相减，然后求距离函数dist，再与偏置b1相乘，最后求径向基函数radbas(n)=exp(-n2)，得到神经元的输出为a=radbas(‖IW-p‖b1)(41)整个RBF网络由两层神经元组成，第1层为S1个如图3所示的RBF神经元，第2层为S2个线性神经元，如图4所示.在第2层开始时，第1层的输出a首先经过线性层权值矩阵LW后与偏置b2相加，再通过一个纯线性(purelin)函数purelin(n)=n，得到网络输出y为y=purelin(LW×a+b2)(42)图4RBF神经网络结构图
Fig.4Structure of RNN比较式(41)与式(9)可见，RBF网络与核方法非常类似，不同之处在于RBF网络的LW需要通过求解一个方程组，而核方法的权重是直接通过归一化计算求得，因此RBF网络预测结果更为逼近完全内插函数估计(注意不是未知函数r)，而核方法计算更为简便［25］.11维数灾难将函数估计推广到高维，则会碰到维数诅咒(curse of dimensionality)［26］(图5），它意味着当观测值的维数增加时，估计难度会迅速增大.维数诅咒有两层含义：一是计算的维数诅咒，指的是某些算法的计算量随着维数的增长而成指数增加.解决方法通常采用优化算法，例如遗传算法、粒子群算法、蚁群算法等［27］.二是样本的维数诅咒，指的是数据维数为d时，样本量需要随着d指数增长.在函数估计中，第二层含义更为重要，这里给予详细解释.图5样本的维数诅咒示意图
Fig.5Dimensionality curse of samples假设一个半径r维数为d的超球，被一个边长为2r维数为d的超立方体所包围，假设超立方体内存在一个均匀分布的点，则由于超球的体积为2rdπd/2/［dΓ(d/2)］，超立方体的体积为(2r)d，因此该点同时也落在超球内的概率P为P=πd/2d2d-1Γ(d/2)(43)令维数d由2逐步增长到20，则对应的概率P如图6所示.显然，当d=20时，P仅为2.46×10-8.因此，若在2维空间中1个样本在半径r的意义下能逼近一个正方形，则在20维空间内，则需要1/2.46×10-8=4.06×107个样本才能在半径r的意义下逼近超立方体.图6概率P与维数d的关系
Fig.6The curve of probability P
against dimensionality d因此，在高维问题中，由于数据非常稀少，导致局部邻域中包含极少的数据点［28］，因此估计变得异常困难.目前还没有较好的办法解决.12结语将文中阐述的方法归结并示于图7.图7非参数回归方法
Fig.7Survey of nonparametric regression methods不同类型方法的特点总结如下：a. 核方法、正则化方法、正态均值模型可以视作最基本最原始的方式.另外，正则化方法与正态均值模型可视作一类特殊的核方法.b. 核方法、局部多项式方法、正则化方法、正态均值模型、小波等方法在大多数情况下均非常类似.这些方法都包含了一个偏倚方差平衡，所以都需要选择一个光滑参数.由于这些方法均是线性光滑器，所以均可以采用第4节中基于CV、GCV的方法.c. 小波方法一般面向空间非齐次函数.如果需要一个精确的函数估计，而且噪声水平较低，则小波方法非常有效.但若面对一个标准的非参数回归问题，而且感兴趣于置信集，则小波方法并不比其它方法明显更好.d. 超完备字典缺陷是丧失了基的正交性，因此估计系数变得复杂；优点是更为灵活，能够采用稀疏的系数描述复杂函数.e. 前向神经网络与RBF神经网络是基于不同的模型独立推导出来的，二者不可混淆.另外，神经网络方法的缺点是一般不考虑置信带，并常用训练误差代替风险函数，容易过拟合；优点是面向应用、思想简单且设计灵活.f. 理论上，这些方法没有大的差别，特别在用置信带的宽度来评价时.每种方法都有其拥护者与批评者，没有哪一种方法目前获得应用上的优势.一种解决方案是对每个问题都利用所有可行的方法，如果结果一致，则选择简单者；如果结果不一致，则必须探讨内在的原因.g. 所讨论的方法能够用于高维问题，然而，即使通过智能优化算法解决了计算的维数诅咒，仍然面对样本的维数诅咒.计算一个高维估计相对容易，然而该估计将不如一维情况下那么精确，其置信区间会非常大.但这并不表示方法失效，而是表示问题的固有困难.