第十章非参数密度估计密度估计的参数解是首先假设一个参数模型,X1,…,X n~i.i.d. f Xθ,其中θ为低维参数向量。
然后通过一些估计方法得到θ,如极大似然估计,矩估计等等。
然后到处密度函数。
此方法的危险性在于初始假设模型的不正确可能导致严重的推断错误。
一种常见的非参数密度估计是直方图,他是一种分段常数的密度估计。
另一种基本的密度估计可通过考虑密度函数如何将概率分配到各区间上受到启发,如果f 足够光滑,我们假设f将某概率不但赋予给x i点,而且赋予给x i周围的一个区域。
因此,要从X1,…,X n~i.i.d.f估计f,将X i周围区域的概率密度累加起来时合理的。
10.1 绩效度量绩效度量是为了评价密度估计量的性质。
令f为整个支撑区域上f的估计量,引入积分平方误差ISE h= f x−f x 2 dx∞−∞如果我们想讨论估计量的一般性质,那么在所有可能的样本上对ISE h进行平均是比较合理的。
积分平均误差为MISE h=E{ISE h}其中的期望是关于分布f。
因此MISE h可以看成是误差(ISE h)关于抽样密度的整体度量的平均值。
又由期望和积分的可交换性,MISE h=MSE f x dx其中MSE f x=E f x−f x 2=var f x+ bias f x2bias f x=E f x−f(x)MISE和ISE都可用来研究选择h值的准则。
两者的好坏已知都有争论,详见Birgit Grunda; Peter Hallb; J. S. Marronc.Loss and risk in smoothing parameter selectionPeter Hall and J. S. Marron.lower bounds for bandwidth selection in density estimation10.2 核密度估计一元核密度估计允许采取灵活的加权方案,即拟合f x=1nhK(x−X i)ni=1(10.6)其中K为核密度,h为固定值,通常称为窗宽。
一些常见的核为:(10.6)的估计量为固定窗宽核密度估计。
而窗宽的大小对估计量有很大的影响,小的窗宽会将密度分配得太局限于观测数据附近,导致估计密度函数有很多错误的峰值;而大的窗宽会将密度贡献分布得太开,从而会因光滑而遗失掉f的一些特征。
10.2.1 窗宽的选择MISE等于积分均方误差。
这表明窗宽的选择是偏差和方差的折衷例10.1(双峰密度)实际上,我们只需对h试一串值,然后选择一个比较合适的。
当然,我们希望得到一个相对正规的窗宽选择程序:如自动算法。
假设K是连续对称的概率密度函数,均值为0,方差0<σK2<∞.令R(g)表示给定函数g的粗超度的度量,定义为R g=g2(z)dz然后假设R K<∞且f足够光滑。
即有二阶有界连续导数。
MISE h=var f x+ bias f x2dxE f x=1Kx−uf u du=K(t)f(x− t)dt在上式中用Taylor级数展开f x− t=f x− tf′x+ 2t2f′′x2+o( 2)因此bias f x2dx= 4σK4R f′′4+o( 4)同样可以计算得到:var f x=1f x R K+o(1)将其对x积分得var f x=R Kn+o(1n)因此MISE =AMISE +o 1nh+h4,其中AMISE h =R K nh + 4σK 4R f′′ 4称为渐进均方误差,h 最小化上式可得= R KK 4 15很多窗宽的选择方法依赖于优化或者找到关于h 的函数的根,例如最小化AMISE(h)的一个近似量。
1、交叉验证许多窗宽的选择是把fx 作为 f 的估计量而与h 联系起来,用某个量Q(h)量化,如果Q 表示根据对在某种意义上对观测数据的拟合程度,那么观测数据在计算fx 和计算拟合程度时候用了两次,这样会对观测提供一个过于乐观的观点,为纠正这一问题,可以采用交叉验证,计算f x 在第i 个点的质量时,模型采用除去第i 个点之外的所有数据拟合,令f−i X i =1n −1 K X i −X jj ≠i表示X i 点处核密度估计量用除X i 外所有数据估计的密度。
交叉验证中一种常见的Q 的选择是伪似然PL h = f−i X i ni =1尽管此方法简单,但其得到的密度估计常常有太多的摆动且对异常值过于敏感。
且其估计量很多时候是不相合的 另一种方法是将积分平方误差写成ISE h = f2 x dx −2E f x + f 2 x dx =R f −2E f x +R (f ) 组后一项是常数,中间项可以通过2n f−i X i n i =1来估计,因此通过关于h 最小化 UCV h =R f−2f −i X i ni =1得到窗宽,此方法称为无偏交叉验证准则。
如果不可能解析计算R(f ),那么计算上式的最好的方式是寻找一个核来简化解析,对于正态核UCV h =R ∅ +1 [1 8π 1∅12 X i −X j −2∅(X i −X j )]j ≠ini =1 详见Nonparametric Econometrics :Theory and Practice例10.2 (鲸的洄游)2001年春天在阿拉斯加巴罗岛附近的海冰边缘对弓头鲸幼仔做了一个目测调查,为了得到摸清其洄游规律,我们估计其洄游可能会按照一个大致的节奏出现。
我们用正态核对鲸鱼幼仔洄游数据进行核密度估计。
通过PL 算法得到h 的最优值为17.4,我们采用了121个样本中的20个,而书本上的结果为9.75.而UCV 的结果同样可以计算。
2、插入法对于一维核密度估计我们知道,最小化AMISE 得到的窗宽为=R KK 4 15(10.24) 因此h 的估计会依赖于未知函数f ,现可以提出多种方法来估计R f ′′Silverman 提出一种初等的方法:把f 替换成正态分布,该正态分布的均值为0,方差为样本方差。
于是有R f ′′ =38π−12σ5,当k 为标准正态核时。
有h pilot = 4π−11038π−12−1σ n−15= 43n15σ此方法称为Silverman 大拇指法,作为产生近似的窗宽的一种方法,此方法是很有价值的。
(10.24)中的R f ′′ 的经验估计是比Silverman 大拇指更好的方法。
基于核的估计量为f ′′ x =d 22 10 L x −X i 0 ni =1=103 L′′(x −X i)ni =1(10.26) 其中 0为窗宽,L 为用来估计f′′的充分可微的核函数。
R (f′′)的估计直接从(10.26)式可得。
Sheather-Jones 方法为:(1)用简单的大拇指法计算窗宽 0,该窗宽用来估计R (f′′) (2)然后用10.24式计算窗宽h 并产生最后的核密度估计。
对用导频核L =∅的一元核密度估计,Sheather-Jones 窗宽的计算如下:R KnσK 4R αf ′′ 15− =0 其中R α f ′′ =1 5 ∅ 4 (X i −X j )nj =1ni =1 α =6 2h 5Ra f ′′ Rb f ′′′ 17R a f ′′ =1 ∅ 4 (X i −X j )nj =1ni =1R b f ′′′ =1n n −1 b 7 ∅ 6 (X i −X j b)nj =1ni =1a=0.920(IQR)/n17b=0.912(IQR)/n1IQR为数据的四分位间距3、极大光滑原则思想:对所有的f均计算h值,然后选择其中的最大值窗宽的选择: =3R K35n 1 5σ详见:The maximal smoothing principleindensity estimation G.R. Terrell例(鲸鱼洄游,续)10.2 核的选择1、艾氏核假设K为各阶矩有限、方差为1的有界对称密度,Epanechnikov证明了关于K最小化AMISE等价于在这些条件下关于K最小化R(K)。
该问题的解是5∗(z/5)的核,其中K∗为艾氏核K∗z=31−z2若z<1 0, 其它2、典则核由(10.29)式可得h K L =δK其中,δK=R KσK41.要想达到与核为K时的窗宽h同样的光滑度,那么核L的窗宽应该取 δL/δK,进一步,如果我们希望对给定的h不同的核可以达到相同的光滑度,我们可以将核进行改进,使得h=1相当于δK的窗宽。
核密度改写为f X x=1nK δK(x−X i)ni=1,其中K δK z=1hδKK(zhδK),按照此方式定义可以给出每种形式的典则核,这样的好处是:单独的h值对每个典则核交换使用不影响其光滑程度。
10.3 对数样条通过三次样条估计f的对数。
令S为包含节点在t1,…,t M上的三次样条,且在L,t1和[t M,U)上为线性的M-维空间,令S的基表示为函数{1,B1,…,B M−1},现在考虑用如下的参数化定义的密度f X|θ对f建模log f X|θxθ=θ1B1x+⋯+θM−1B M−1x−cθ.其中exp{c(θ)}=exp{θ1B1x+⋯+θM−1B M−1x}dxUL该模型我们通过两个条件来保证。
(1)L>−∞或θ1<0(2)U<∞或θM−1<0对给定的数据值x1,…,x n,该模型的对数似然为lθx1,…,x n=log f X|θ(x i|θ)ni=1在c(θ)的限制下最大化上式可以得到极大似然估计θ。
为估计模型,我们取f x=f X|θ(x|θ)作为f(x)的极大似然对数样条密度估计关于节点的摆放:令x i表示数据的第i个次序统计量,定义一个近似的分位数函数为q i−1n−1=x i,1≤i≤n,对一列数,0<r2<⋯<r M−1<1,M个节点将放在x1,x n及由q r2,…,q(r M−1)标记的次序统计量的位置上。
内部节点的放置由下面决定,对1≤i≤M2,n r i+1−r i=4∗max4−ϵ,1∗max4−2ϵ,1∗…∗max{4−i−1ϵ,1}ϵ的选择满足当M为奇数时r M+1=1/2,或当M为偶数时,r M+r M+1=1/2.其余节点保持分位数对称,于是对M2≤i≤M−1,r M+1−i−r M−i=r i+1−r ir M=1上面假定的M是预先给定的,其实实际上有很多选择M的方法:概括如下,首先把少量节点放在给定的位置上,建议的最小值为min{2.5n 15,n4,n∗,25},其中的n∗为不同数据点的个数,然后其他的节点一个个计入到现存的集合中,每次循环中,在该节点不存在的模型满足Rao检验统计量的最大值的位置上增加一个节点,直到总节点数达到min{4n 15,n4,n∗,30},或是没有节点可以继续添加为止然后各节点依次逐个删除。