非参数回归模型及半参数回归模型

1、精选优质文档-倾情为你奉上精选优质文档-倾情为你奉上专心-专注-专业专心-专注-专业精选优质文档-倾情为你奉上专心-专注-专业第七章 非参数回归模型与半参数回归模型第一节 非参数回归与权函数法 一、非参数回归概念 前面介绍的回归模型，无论是线性回归还是非线性回归，其回归函数形式都是已知的，只是其中参数待定，所以可称为参数回归。参数回归的最大优点是回归结果可以外延，但其缺点也不可忽视，就是回归形式一旦固定，就比较呆板，往往拟合效果较差。另一类回归，非参数回归，则与参数回归正好相反。它的回归函数形式是不确定的，其结果外延困难，但拟合效果却比较好。设Y是一维观测随机向量，X是m维随机自变量。在第四章

2、我们曾引进过条件期望作回归函数，即称g (X) = E (Y|X) （7.1.1）为Y对X的回归函数。我们证明了这样的回归函数可使误差平方和最小，即（7.1.2） 这里L是关于X的一切函数类。当然，如果限定L是线性函数类，那么g (X)就是线性回归函数了。 细心的读者会在这里立即提出一个问题。既然对拟合函数类L(X)没有任何限制，那么可以使误差平方和等于0。实际上，你只要作一条折线(曲面)通过所有观测点(Yi，Xi)就可以了是的，对拟合函数类不作任何限制是完全没有意义的。正象世界上没有绝对的自由一样，我们实际上从来就没有说放弃对L(X)的一切限制。在下面要研究的具体非参数回归方法，不管是核函数

3、法，最近邻法，样条法，小波法，实际都有参数选择问题(比如窗宽选择，平滑参数选择)。 所以我们知道，参数回归与非参数回归的区分是相对的。用一个多项式去拟合(Yi，Xi)，属于参数回归；用多个低次多项式去分段拟合(Yi，Xi)，叫样条回归，属于非参数回归。 二、权函数方法非参数回归的基本方法有核函数法，最近邻函数法，样条函数法，小波函数法。这些方法尽管起源不一样，数学形式相距甚远，但都可以视为关于Yi的线性组合的某种权函数。也就是说，回归函数g (X)的估计gn(X)总可以表为下述形式：（7.1.3）其中Wi(X)称为权函数。这个表达式表明，gn(X)总是Yi的线性组合，一个Yi对应个Wi。不过W

4、i与Xi倒没有对应关系，Wi如何生成，也许不仅与Xi有关，而且可能与全体的Xi或部分的Xi有关，要视具体函数而定，所以Wi(X)写得更仔细一点应该是Wi(X；X1，,Xn)。这个权函数形式实际也包括了线性回归。如果，则，也是Yi的线性组合。在一般实际问题中，权函数都满足下述条件：（7.1.4）如果考虑在第五章介绍的配方回归与评估模型曾有类似条件，不妨称之为配方条件，并称满足配方条件的权函数为概率权。 下面我们结合具体回归函数看权函数的具体形式。 1核函数法 选定Rm空间上的核函数K，一般取概率密度。如果取正交多项式则可能不满足配方条件。然后令（7.1.5）显然。此时回归函数就是 （7.1.6）

5、2最近邻函数法首先引进一个距离函数，用来衡量Rm空间中两点u = (u1，um) 和v= (v1，vm) 的距离u-v。可以选欧氏距离，也可以选。为了反映各分量的重要程度，可以引进权因子C1，,Cn，使Ci也满足配方条件。然后将距离函数改进为 （7.1.7）（7.1.8） 现在设有了样本(Yi，Xi)，i=1,n，并指定空间中之任一点X，我们来估计回归函数在该点的值g(X)。将X1，Xn按在所选距离意义下与X接近的程度排序： （7.1.9）这表示点与X距离最近，就赋以权函数k1；与X距离次近的就赋予权函数k2。，等等。这里的n个权函数k1,kn也满足配方条件，并且按从大到小排序，即（7.1.1

6、0）就是 （7.1.11）若在Xi-X, i=1,n中有相等的，可将这n个相等的应该赋有的权取平均。比如若前两名相等，X1-X=X2-X, 就令W1 = W2=。这样最近邻回归函数就是（7.1.12）ki尽管是n个常数，事先已选好，但到底排列次序如何与X有关，故可记为ki(X)。 三、权函数估计的矩相合性 首先解释矩相合性的概念。如果对样本 (Yi，Xi)，i=1,n构造了权函数Wi = Wi (X)=WI(X;X1,Xn)，有了回归函数g(X)的权函数估计，当Y的r阶矩存在(E|Y|r0, 当n时， （7.1.15）(3)当n时， （7.1.16）则Wi是矩相合的权函数。 定理条件可以作一些

7、直观解释。条件(1)可以作如下理解，因为权函数是概率权，必有|Wi|1，i=1,n。于是（7.1.17）这里取的是C=1。因此条件(1)可以说不叫做一个条件。条件(2)是说，与X的距离超过一定值的那些Xi，对应算出来的权函数之和很小，也就是说，权函数的值主要取决于那些与X邻近的Xi的值。这个条件合理。条件(3)是说，当n越来越大时，各个权系数将越来越小，这也是合理的要求。 在证明本定理之前，先证两个引理。引理7.1.1 设概率权函数Wi适合定理7.1.1的条件(1)及(2)，又对某个r, E|f(X)|r0，存在0，当u-v时，|f(u)-f(v)|(/2)1/r。于是（7.1.19）其中，此

8、处X表示具体取值。由条件(2)，上式右边第二项依概率收敛于0且不大于1。依控制收敛定理有（7.1.20）故存在n0，使当nn0时，有 （7.1.21）因此当nn0时，有 （7.1.22）于是对这种一致连续的f，引理得证。 证毕对一般的函数f，取一个在Rm上连续，且在一有界域之外为0的函数，使，且，这里是事先指定的。因为（7.1.23）右边括号里第三项等于；第一项根据条件(1)不超过；因为在Rm上有界且一致连续，由前面已证结果知当n时，第二项将趋于0。因此（7.1.24）是任意的，故引理得证。证毕引理7.1.2 设Wi为满足定理7.1.1三个条件的概率权，函数f非负且，则（7.1.25）证明 定

9、义一组新的概率权函数，由于0Wi1, 故01。于是由引理7.1.1，有 （7.1.26）因为01，由条件(3)知 （7.1.27）故由控制收敛定理有 （7.1.28）综合两个极限式可知本引理成立。证毕 下面我们证明定理7.1.1。先设r=2, 则E(Y2)。令（7.1.29）由E(Y2)知E(Z2)，故h (X) = E (Z2|X)（7.1.30）存在。又（7.1.31）还须注意：f (Xi) = E (Yi|Xi) (而非E (Y|Xi)。因此按定义而因为 (X，Y) 与 (Xi，Yi)同分布，有E (Y|X=x) = E (Yi|Xi=x)。故现有（7.1.31）因 E | f (X)

10、|2，依引理7.1.1，有（7.1.32）又若将X固定为x，则有（7.1.33）注意到当X固定为x而X1，, Xn也给定时，Wi (x)成为常数，而Z1，, Zn在给定X1，Xn时，条件相互独立，再注意到E (Zi |Xi)=0，由上式有因此式对一切x都成立，有 （7.1.34）考虑到E (h (X)，h0，由引理6.4及上式，知 （7.1.35）合并考虑 (7.1.31)，(7.1.32) 和上式，得。这证明了定理当r=2的情况。 现在设r1,E |Y| r0，先找K0，使当K K0时，对一切n成立 (7.1.40)。又依 (7.1.41)，找K1，使当KK1时有E | E ( Y | X

11、)- E ( Y(K) | X )| r/3。固定K = max ( K0, K1)。根据上式，存在n0, 使当nn0时（7.1.45）这时由 (7.1.42)推出：当n n0时有 （7.1.46） 这就证明了权函数的矩相合性。证毕关于权函数估计的收敛性质还有更多更深入的讨论，如逐点矩相合性，强相合性等，有兴趣的读者可参看有关专着。这里引述Stone的成果，一是因为它是基本的，可以作为入门的引子；二是因为它是一般的，概括了核估计、最近邻估计、样条估计、小波估计等具体形式。算例7.1.3 一元非参数回归本算例利用核估计给出一元非参数回归。计算过程如下。-一般非参数回归模型计算程序, 例 7.1.

12、4 模型及数据结构说明:本项程序计算一般非参数回归模型: Y(i)=g(t(i) +(i) i=1,2,.,n, 0= t =1 其中函数 g 未知待估.资料准备要点: 因变量 Y 在数据第一列, 自变量 t 是 1 维, 例713.D 数据文件中, n=50要打印原始资料吗? 0= 不打印, 1= 打印 (1)打印 Y 的原始资料 1. 1. 1. 2. 0. 1. 1. 0. 1. 1. 0. 1. 1. 1. 1. 0. 0. 0. 1. 1. 1. 1. 1. 1. 1. 1. 0. 0. 0. -0. 0. 0. 1. 0. 0. 0. 0. 0. 1. 1. 1. 0. 0. 1.

13、 1. 1. 0. 1. 1. 0.打印 X 的原始资料 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.请决定非参数回归的方法: （0） 0= 固定自变量窗宽的核函数法. 这需要事先将自变量变换为 0=t=1. 1= 固定自变量资料点数的平滑法. 这需要自变量资料等距并顺序排列.请键入核函数的窗宽选择h(1/N=h0，由条件 (7.2.7)，存在

14、充分大的T0，使 （7.2.10）这里，并且（7.2.11）于是（7.2.12）由的任意性，可知。这就说明fn(x)是f(x)的渐近无偏估计。再利用X1，,Xn的独立性，有（7.2.13）类似于渐近无偏性的证法可得 （7.2.14）于是（7.2.15）这就说明对一切x，fn(x)均方收敛于f(x)，因此，这就证明密度核估计的相合性。 二、使用正交多项式核的密度及其偏导数核估计的收敛速度 上一段研究的密度核估计的收敛性，针对的是使用概率密度核函数K，它非负，积分为1，从而可以肯定保证密度核估计函数fn(x)非负且积分为1。只是它的收敛速度不会超过。为了提高收敛速度，统计工作者使用正交多项式作理论

15、上的研究，取得不少成果。这里介绍的是本书作者的研究成果，近期发表在国际数学杂志“Communications in Statistics”上。它是直接研究多元密度，并连带一般偏导数的核估计给出收敛速度。记多元密度f(t)的s阶混合偏导数为（7.2.16）这里。使用多元核函数作出f(s) (t)的估计如下：（7.2.17）其中是构造核函数的正交多项式空间维数，可以任意取定。不仅决定于s, 而且决定于s1, sp,且满足：（7.2.18）其中u0是一正常数，u = (u1,up)。我们以C表示某一合适常数，各个C可不相同。 还满足：（7.2.19）这种多元核函数可以如下构造：（7.2.20）其中是

16、普通一元核函数，满足：（7.2.21）及（7.2.22） 这种核函数具体构造及改进我们放到下一段再统一研究。下面研究的收敛性。我们假定偏导函数局部有界，即存在与对t的各分量求偏导次数r1,rp无关的，0, r1+rp = r, 使当，且tXt且t+Xt时，有 （7.2.23）这里Xt是t的样本空间。同理定义f(t)局部有界。En () 表示对n个样本求数学期望。定理7.2.1 设f (r) (t), f(t)局部有界，则 （7.2.24）（7.2.25）证明 由t(1),t(n) 的i.i.d.，令，注意，有（7.2.26）再由多元Taylor展式、多项展式及核函数正交条件得（7.2.27）这

17、里，由f(r) (t)局部有界，核函数有界，积分域有界，可得 (7.2.24)。又 （7.2.28）（7.2.29）可知 (7.2.25)成立。证毕在s=1时，由 (7.2.16) 我们把都记作了f(1) (t), 把它们排成向量得。相应也代表了p种核估计。由 (7.2.17)知它们的核函数构造不同，满足的正交条件不同，也把它们排成向量得。由定理 (7.2.2)有（7.2.30）进一步有（7.2.31）设122，由 Jensen 和不等式有（7.2.32）于是有 推论1 设，（7.2.33） （7.2.34）这就证明了使用正交多项式核的密度及其偏导数核估计的收敛速度。在本书第十章第四节要引用这

18、些结果。 三、密度核估计的连续性及光滑性 这一段介绍本书作者提出的一种正交多项式，用它构造的一元到多元密度及其偏导数的核估计，在样本抽定时，保持连续性，在样本数趋于无穷时可以保持好的收敛速度。密度核估计是一随机函数，它利用随机抽得的历史样本x(1) ，x (n) 构造fn (x)，去估计母体的密度f(x)。它的收敛性是一种大样本性质。对于一个具体的核函数和一个具体的fn (x)的构造，一旦历史样本抽定转入统计计算，fn (x)就是一个普通的函数。这时我们自然要考虑它的分析性质，例如连续性和光滑性。因此，密度核估计的连续性和光滑性是对任意抽定的历史样本而言，它是一种小样本性质。从统计计算的角度，

19、仅仅研究大样本性质是不够的。如果核估计呈跳跃间断，得到的参数估计将随当前样本x的连续变动而发生剧烈跳跃，使其难以进入实用，许多文献要么忽略了核函数的构造，要么给出的核函数不满足连续性光滑性，Lin(1975)构造密度及其(偏)导数核估计如下。（7.2.35）（7.2.36）他进一步具体给出了正交多项式的构造，在n=3时我们画出所给函数式的图像。K1(u)uu120-4811K0(u)09图7.2.3.1. 当0u1时，（7.2.37）（7.2.38）显然这样的fn (x)与都不连续。我们试图寻找截断后仍然连续的正交多项式，从而使密度及其(偏)导数的核估计连续，同时保持较高的收敛速度。 我们先考

20、虑一元密度核估计的连续性、光滑性、收敛性。以下给出的正交多项式与Lin给出的正交多项式区别在于连续性和光滑性，其正交性是一样的。（7.2.39）（7.2.40）（7.2.41）令H是r阶行列式，其第1至r-1行的元素为，第r行元素全为1。将H的第一行换上(u/ut)j得H0，将H的第二行换上(u/ut)j得H1，ut是一常数。显然H0(0) = H0 (ut) = H1(0) = H1(ut)=0。 再令（7.2.42）（7.2.43）则又K0(0) = K0(ut)=0，可见K0是满足正交性及连续性的核函数。在r = 3, ut=1时我们画出它的图形 (图7.2.3.2)。当0uut=1时（

21、7.2.44）同样容易验证K1(u)的正交性及连续性。K1(0)=K1(ut)=0 在这样的构造里，密度核估计的光滑性通过密度导数核估计的连续性实现。下面我们再说明多元密度核估计的连续性、光滑性及收敛性。u-2K04图7.2.3.2设有p元密度f(x), x = (x1, xp), 对于其各阶混合偏导数，我们使用多元核函数，作出它的估计：（7.2.45）其中 （7.2.46）同时要求fn在全空间连续，即（7.2.47）多元核函数要满足： （7.2.48）（7.2.49）且这种多元核函数构造如下：（7.2.50）其中是普通一元核函数，满足（7.2.51）（7.2.52）且。 这种一元核函数构造如

22、下：作行列式，其中第1至r-1行元素为，最后一行元素全为1，将其第ti+1行换上，得r阶行列式（7.2.53）（7.2.54）再令（7.2.55）这样我们已经构造了多元密度及其偏导数的核估计并验证了它们的连续性、光滑性及收敛性。 四、改进多元密度核估计的交互投影迭代算法本段介绍交互投影算法在改进多元密度核估计非负性方面的应用。多元密度核估计是（7.2.56） 这里K ()是核函数，x(j)，j=1,n，是样本。如果K()是一个概率密度函数，(非负、积分为1)，则fn(x)的均方误差的收敛速度不会超过O ( n-4/5)。如果K ()是一个正交多项式，则fn(x)的均方误差的收敛速度可以任意接近

23、O (n-1)，但这时fn(x)不再能保证非负、积分为1。Gajek (1986)利用凸集间的交互投影迭代算法来改进用正交多项式构造的一元密度核估计。该算法可以将fn (x)改进为非负且积分为1的函数，同时保证它的均方误差收敛速度不变，本段将Gajek的方法用之于多元密度核估计，并为它重新写了证明。 多元密度核估计的详细构造上段已述。下面先叙述属于Gajek的迭代算法。对于fn(x)定义加权的均方误差式（7.2.57）这里h是一个非负的权函数。应该说R (fn, f)是一个合适的评价标准。迭代算法是： (1)令，且置k = 0; (2)令，再检查。若Ck+1=1，则令而完成迭代；(3)令； (

24、4)置k=k+2并转向步骤(2)。 从几何直观上看，步骤(2)就是去掉函数的负值而将其改写为零，此时可能函数积分超过1。于是有步骤(3)，就是将函数整体向下拉一点，以使积分为1。此时可能又会有负值出现，于是重复步骤(2)。如此反复。Gajek证明了迭代过程收敛。定理7.2.2 设（7.2.58） （7.2.59）则 (1)上述迭代过程收敛，即存在，且。(2)在 (7.2.57)的加权均方误差意义下，至少保持fn (x)的收敛速度，即（7.2.60）在正式证明定理7.2.2之前，我们先叙述三个引理，定义内积 （7.2.61）令 （7.2.62）满足(7.2.62)的全体p元函数构成内积空间，由内

25、积(7.2.61)导出的距离记作，在空间定义 （7.2.63） （7.2.64） （7.2.65）显然，所有的F*，F+和F1都是凸集。引理7.2.1 设，且F0是一凸集，。则f0是fn在F0上的投影当且仅当，有（7.2.66）证略。引理7.2.2 令（7.2.67）则f+是fn在F+上的投影。利用引理7.2.1，证明是容易的。引理7.2.3 令 （7.2.68）则f1(x)是fn(x)在F1上的投影。 利用引理7.2.2，证明也是容易的。 现在我们叙述定理7.2.2的证明： (1) 因为F*非空，F+和F1之间的距离为零。由引理7.2.2和引理7.2.3我们知定理7.2.2的迭代算法也就是两

26、个凸集间的交互投影。这个迭代过程一定收敛，设收敛于，。 (2)由引理7.2.5 ，有（7.2.69）在k次迭代后，我们有（7.2.70）令k并取数学期望，由Fubini定理，我们有 （7.2.71）证毕实际计算时可以取控制精度，在步骤2中Ck+1=1可用Ck+1-1替代，因为与x无关，所以有如下形式（7.2.72）这里常数。 算例7.2.4 随机数发生、直方图显示与密度核估计本算例程式有4个功能：发生给定密度函数的随机数；作直方图(二维或三维)；作饼图；作密度函数的核估计。其中发生随机数的程序附有常见分布16种，参数也随使用者指定。如果还要另外的函数，也只需改写一行。作直方图与饼图程序是用C语

27、言写的，彩色显示，10个区间或20个区间色彩各不相同，十分绚丽，调用也十分方便。先发生伪随机数。-16 种指定分布的随机数发生程序 最多发生 5000 个随机数请指定需要发生的随机数的分布函数代码 1: 标准正态分布 N(0,1) 2: 一般正态分布 N(,) 3: 卡方分布 2 4: t 分布 5: F 分布 6: 对数正态分布 7: WEIBULL 分布 8: 指数分布 9: 柯西 (CHUCHY) 分布 10: 贝塔分布 (2,2) 11: 均匀连续分布 U(0,1) 12: 均匀离散分布整数 13: 负二项分布 14: 几何分布 15: 超几何分布 16: 泊松分布 请输入需要发生的随

28、机数个数 n (x1,x2,.xn), n = ? (500)请输入发生随机数的种子(任一奇数) NRAN (11)请输入 t 分布的自由度 (10)要显示发生的随机数吗? 0= 不显示, 1= 要显示 (0)资料存在哪个文件中? 0= 不存盘 1= C21.D, 2= C22.D, 3= C23.D, 4= C24.D, 5= C25.D 6= C11.D, 7= C12.D, 8= C13.D, 9= C14.D, 10= C15.D正将数据文件存盘, 请稍侯 资料已存盘,计算结束。 - 再将刚才发生的伪随机数用直方图显示。图7.2.4.1下面我们再对上述伪随机数出出密度核估计。-密度函数

29、核估计计算程序, 例 7.2.5 请输入资料长度(观测点数, li725.d是500) N: （500）要显示原始数据文件吗? 0=不显示; 1=显示. （0）请选择密度估计的核函数: （） 1: K(x)=1-|x|, |x|1; 2: K(x)=exp(-x*x/2)/sqrt(2*3.1416); 3: K(x)=exp(-|x|)/2; 4: K(x)=1/(3.1416*(1+x*x); 5: K(x)=(1-x*x)*3./4., |x|1; 要想得到光滑的密度核估计图像, 样本数要多一些, 比如N=500; 窗宽适当, 比如 h=1-5; 要计算的核函数个数适当, 比如 M=20

30、-100。通过这些参数的调整, 一定可以得到比直方图要好的密度核估计结果。窗宽大致相当于直方图里 X 轴各个分组条形的宽度, 但核估计分组逐点改变.请输入核函数窗宽 h ( h 须为正数, 最好 h 1): （1）请输入您想计算的密度核估计点数 M ( 10=M0。极小化ISE(bt, bu)得到的窗宽记为 (bt,ISE，bu,ISE)极小化MISE(bt, bu)得到的窗宽记为 (bt,MISE，bu,MISE)。对于二元函数g的偏导数，我们记（7.2.80）对于支撑在-1,1上的W，令 （7.2.81）假定g有各二阶连续偏导数，g(2,0)和g(0,2)不完全为0，核函数K是Lipsch

31、itz连续，bt+bu0, nbtbu，(n)。由标准的渐近理论可得（7.2.82）其中 （7.2.83）是渐近MISE。函数Itt，Iuu，Itu，If由下式给定：（7.2.84）（7.2.85）（7.2.86）（7.2.87）令 (bt,AMISE，bu,AMISE)表示极小化AMISE(bt, bu)得到的窗宽选择： （7.2.88）（7.2.89）插入方法需要残差方差2以及函数Itt，Iuu，Itu的估计。如果设计密度f未知，If可由下式估得： （7.2.90）下面我们介绍插入法窗宽的选择。 我们先再回顾一下窗宽选择的意义。大家可以想象直方图，那些长条条的宽度就是窗宽。如果把那些长条条

32、取得特别宽，比如极言之，整个直方图就一个长条条，那就太平滑了，变成了均匀分布。如果把长条条取得特别窄，比如极言之，一个长条条里至多只含一个样本点，那些不含样本点的长条条的高度就等于0，整个直方图就乱起乱落。所以合适的窗宽选择是有必要的。这里介绍的是 (bt, bu)的插入法选择。我们使用下列偏导数的核估计：（7.2.91）（7.2.92）这里L2，M是核函数，t, u, t, u是窗宽，我们称为试验窗宽。将,代入Itt，Iuu，Itu的定义式里就可以得到它们的估计值，记作。选择试验窗宽的迭代法如下。令（7.2.93）这里与是在bt,AMISE与bu, AMISE的定义式里以代替2，以代替而得到

33、的。则获得（）的迭代算法如下： （1）置初值 （2）对i =1,2, 迭代公式是 （3）在i*次迭代后停止并且令 ( 这个方法有些类似于搜寻固定点来确定窗宽。初值的设置取是根据经验，并非必要。膨胀系数是挑选来使渐近最优窗宽不依赖于c, d (当然要c0, d0)。窗宽与都有的收敛速度。实际演算时迭代次数大约是i*在5到9之间。 下面我们更具体给出算例。取核函数（7.2.94） （7.2.95）积分域A=(t, u)|0t1, 0u1,权函数 （7.2.96）窗宽bt, bu有上界，下界取法则使 （7.2.97）覆盖0.05，0.950.05，0.95。膨胀系数里，取c = 1.5，d = 0.

34、25。插入算法可以用于非矩形的设计密度支撑里，不过用于矩形支撑集当然更好。令密度 （7.2.98）设计点集 (ti, uj), i=1, nt , j=1, nu可由下式给出：（7.2.99） （7.2.100）相应的剖分集为 （7.2.101）i = 2, nt-1, j=2, nu-1。回归函数可取二元正态： （7.2.102）或者两个二元正态的混合。ti, ui的设计密度f可取为线性： （7.2.103）或正态 （7.2.104）nt与nu大约在10到20之间。有了这些参考消息，就差不多可以发生资料然后代入插入法迭代公式里计算了。Herrmann, Ward, Engel及Gasser

35、(1995)称他们计算结果十分令人满意。笔者认为这样搞真是把简单的问题弄复杂了。这样一些资料成功，换一些资料不一定成功。窗宽选择最重要的还是要凭经验与直观观察。如果先凭经验选窗宽，拟合一次后看图像，再调整窗宽，那么把窗宽的自动选择搁置一边，单纯的二元核回归并不复杂，可以为成千上万的普通读者掌握。手中资料是Yi, Xi (ti, ui)。要拟合的模型是 (7.2.73)：Yi = g (ti, ui)+i, i=1, n，则由(7.2.77)有： （7.2.105） （7.2.106） （7.1.107）大家看，只剩下调整窗宽问题，简单极了。第三节 非参数回归模型的样条拟合 一、样条回归的基本概

36、念 非参数回归模型样条拟合可以考虑多元函数，几何形象就是作曲面拟合。不过多元样条拟合的理论分析与计算难度要更大一些，Annals of Statistics在前两年才有过一次较大规模的讨论。我们这里只介绍清楚一元样条拟合就可以了。 一元自变量的取值区间本来可以从-到+，经过一个线性变换，总可以变到0,1。为了理论分析方便，我们限于考虑自变量定义域为0,1的情况。并且按一般有关文献惯例，将自变量记作t。这样我们要考虑的模型就是（7.3.1）其中g(t)是一个光滑曲线，(t)是白噪声过程，E(t)=0，E(s)(t)=0，(st)，E(2(t)=2。Y (t)是t = t1, t2,tn, 0t1

37、t21时的观察。我们的目的是根据资料Y (tj)=Yj, j=1,2,n来重构函数g。 根据混杂有白噪声的观察资料Yj, j=1,n去重构光滑函数g，这在统计学中叫非参数回归，在数字信号处理中则叫信噪分离。使用信噪分离这个词能帮助我们清楚理解非参数回归的实质。我们假定未知函数g的m阶导数g(m)是平方可积的，记作绝对连续，v=0,1,m-1,g(m)L20，1，即函数g的估计记作gn，因为gn是根据Y1，Yn来作出的。如果单纯考虑极小化 （7.3.2）怎么样?这样的函数f (t)可以精确通过每一点Yj，也可以光滑，还确实可以是样条函数(由分段低次多项式光滑连接而成)。可是这样的拟合有什么实际意

38、义呢?噪声没有去掉。如果单纯考虑极小化 （7.3.3）怎么样呢?函数f倒是充分光滑，比如取f=c,则S2=0。噪声是给去掉了，可是完全没有考虑经过任何一个Yi。真正是将婴儿连同洗澡水一起泼掉了。于是就有了折衷的办法，极小化（7.3.4）这里称为光滑因子。当=0时，S = S1, 当=时，S = S2。的作用就是在“尽量通过数据点”与“尽量去掉白噪声”之间取调和作用。满足S的极小化的解，就取作非参数回归模型的解gn (t)。可以证明，这个解是2m-1次的平滑样条多项式。由于这个解依赖于参数的选择，故我们有时也记作gn,(t)，不过为省事，经常记作gn(t)。 我们的任务就转到怎样确定参数。Rei

39、nsch (1967)建议，如果2已知的话，可取使（7.3.5）Wahba (1975) 在等距间隔资料及一定的平滑和周期条件下，得到的最佳选择的理论结果。他所选择的是如下误差平方和的极小化结果： （7.3.6）理论分析结果表明，真正最佳的参数选择，应该使(7.3.5)左端比右端稍微小一点。当然这个结果不太现实，小一点，到底小多少?同时2可能并不知道，函数g也是未知的。如果2是已知的，则一个好的可以按如下方法获得。定义nn矩阵A()，它依赖于tj和，且满足 （7.3.7）因为gn,(t)是Y1，,Yn的线性函数，所以A()是存在的。然后令 （7.3.8）这里Y= (Y1,Yn)，g = (g

40、(t1),g (tn)，范数是欧氏范数。经过简单的计算有 （7.3.9）可以证明，由下式定义的是ER()的无偏估计：（7.3.10）即可以证明 （7.3.11）因此，的极小化能得到的好的选择。Mallows(1973)甚至用这样的参数估计去作岭回归。Wahba等提出，当2未知时，的估计可以取作下式的极小化解（7.3.12）并把这个估计称作广义交叉核实 (Generalized Cross-Validation, GCV) 估计。可以证明，在一定条件下，当n充分大时，在ER()的一个邻域内还可以证明，在合适条件下，存在序列，它是EV()的极小化解，且有性质（7.3.13）下面我们叙述一下广义交叉

41、核实GCV的动因。交叉核实的动因很简单：设是使用了除第k个点以外 (也算是一种回避吧) 的所有数据点的平滑样条，我们现在想用去预测Yk，用预测好坏来衡量的选择好坏。用数学语言写就是，设是下式的的函数解： （7.3.14）令 （7.3.15）原始的交叉核实就是取V0()的极小化解作为的估计。对于这个交叉核实办法，在自变量等距离情形，Wahba与Wold曾研究过，并且用Monte Carlo实验显示V0()的极小化解是R()的极小化解的极好估计。实验对各种各样的g与2都显示好的结果。理论研究则对自变量等距（）和为周期函数的情形获得满意的结果。我们把自变量等距且函数为周期函数称为对称场合。注意在对称

42、场合所有数据点都是被对称对待的。因此tk的预测误差与其它点tj的预测误差一样。一般地，令（7.3.16）这里权Wk()是作为资料点可能非等距或者函数g非周期函数的一种补偿。如果（7.3.17）k=1,2,n，其中akk()是A()的主对角元素，则 (7.3.16)中的V()就等于(7.3.12)中的V()，并且 (7.3.13)成立。于是我们得到序列Wk可使 (7.3.13) 成立。在证明 (7.3.13)成立的过程中我们得到了平滑样条gn,。有趣的是，在自变量等距场合，在分段拼接Bernoulli多项式的时候，Gram矩阵是一个循环矩阵。这里的广义交叉核实方法也被用来在求解Fredholm积

43、分方程时计算规则参数。 二、平滑样条的构造上一段讨论什么是样条回归，给出的是隐式解，是满足（7.3.18）极小化的m次导数平方可积函数gn,，是2m-1次多项式。而平滑参数的最佳选择是取广义交叉核实函数（7.3.19）的极小化解。因此，我们应该寻求gn,与A()的显式表达式。这可以借助于贝努里(Bernoulli)多项式。设Br(t),r=0,1,是t0,1上的Bernoulli多项式。它的定义是（7.3.20）其中积分常数的选择应使（7.3.21）例如（7.3.22）下面我们介绍怎样用这样一些贝努里多项式构造平滑样条函数的显式解。令x表示x的小数部分，例如3.14=0.14，因此这个记号不是

44、通常的取整函数，而正好是相反。定义（7.3.23）当然，如果t(0，1)，t= t。令Lk，k=0,1,2,为一族线性函数（7.3.24）（7.3.25）则 （7.3.26）定义贝努里核kr(s, t)为（7.3.27）Abramowitz和Stegun(1964)证明了，贝努里核与贝努里多项式有如下关系： （7.3.28）并且可以验证，对p=1,2,r-2及s, t0,1有（7.3.29） （7.3.30）（7.3.31） （7.3.32）（7.3.33） 有了这些准备，我们现在可以写出用分段贝努里多项式表达的gn,的显式表达式。定理7.3.2 若函数，函数gn,是极小化（7.3.34）的函

45、数f的解，则对nm, gn,是唯一的，并可表达为 （7.3.35）这里= (0,1,m)和= (1,2,n)由下列式子给出：（7.3.36）（7.3.37）（7.3.38）T是n(m+1)维的矩阵，其第jr的元素为（7.3.39）是(m+1)(m+1)维的矩阵，除了一个元素(m+1)(m+1)=1外，其余全为0。矩阵M由下式给出（7.3.40）其K是nn矩阵，第jk的元素为（7.3.41）I是nn单位阵。矩阵A()由下式显式表达：（7.3.42）证明 我们先证明函数gn,属于函数空间。首先定义内积（7.3.43）则可以证明 （7.3.44）是的赋予上述内积的再生核，对一切t, Q (t，)，对

46、f，t0,1,Q(t,)，f= f (t)。 Wahba (1970)证明了，gn,必定属于函数空间kQ：（7.3.45）这里。然而可以看到是由与生成的，故（7.3.46）因此必有gn,kk 。于是存在某个r与j，使 (7.3.35)成立。将 (7.3.35)代入到 (7.3.34)中，再使用 (7.3.33)的求导公式得：（7.3.47）向量和应被选择使上式极小化。对上式右边求导并令导数为0，就可完成定理的证明。证毕三、广义交叉核实普通的交叉核实函数是（7.3.48）其中是下式的f的极小化解（7.3.49）就是每次作函数极小化时空出一个点Yk不参加回归，然后用回归函数在这一点的拟合值去与实际

47、观察作比较求出差值的平方。当这样的k取遍1到n时，就有误差平方和V0()，它是的函数。改变可以计算出不同的V0()，取使V0()较小或最小的参数为的选择。如果我们在极小化 (7.3.4)过程中，用取代Yk，即在 (7.3.4) 中用Y1，Yk-1，Yk+1，Yn去代入使极小化，就可以得到。用引理形式说就是 引理7.3.3 设nm, 令gn,(tj，k, Zk)为下式的f的极小化函数解（7.3.50）则（7.3.51）证明 令，f为属于但不同于h的函数，我们有 (7.3.52）比较左边与最右边的表达式，我们看到h是以Zk代换Yk后的极小化问题的解。证毕引理7.3.4 在上述记号下，有（7.3.5

48、3）证明 令，由于对每-t, gn,(t)是Yk的线性函数(见定理7.3.2)，再由上一引理，我们有（7.3.54）由此可以立即解出引理所证之式。证毕记A()的元素为ajk，j, k=1,n，则（7.3.55）因此 （7.3.56）于是由引理7.3.4，普通交叉核实函数可以写为（7.3.57）为了引进广义交叉核实函数V ()，考虑平滑问题中周期函数的情况，就是要找到g，g是积分为0的周期函数，且是下式的极小化解： （7.3.58）函数g是积分为0的周期函数意味着（7.3.59）可以证明，引理7.3.3证明中设定的hn,可以表为（7.3.60）其中 （7.3.61）这里KM-1就起着矩阵A的作用

49、。如果，则KM-1对每一都是循环的，并且（7.3.62）于是 (7.3.57)之V0()可表为（7.3.63）这个表达式由Wahba与Wold (1975)获得，广义交叉核实正是采取了后面的形式。因为A是对称的，可对A作分解：（7.3.64）这里D2是对角阵，U是正交阵。令（7.3.65）这里矩阵Wnn的元素 （7.3.66）并且可以证明， （7.3.67）其中W*是W的共轭矩阵，而nn的对角阵D的对角元素为 （7.3.68）这里（7.3.69）（7.3.70）于是可以证明 （7.3.71）是循环矩阵。再令（7.3.72）则（7.3.73）如果将普通交叉核实作用于这里的资料，结果就是广义交叉核

50、实函数V ()。注意 (7.3.69)揭示了周期函数与等距间隔下平滑样条的低通滤波性质。hn,的Fourie系数为（7.3.74）它与资料（7.3.75）的Fourier系数通过方程 （7.3.76）相联系，这里 （7.3.77）当vn时，我们可以在 (7.3.69)中取=0项而得的近似式，于是 （7.3.78）这里函数称作Butterworth滤波器，经常使用于电气工程中，其半功率，B*是B的共轭函数。将广义交叉核实函数（7.3.79）与普通交叉核实函数（7.3.80）比较知道，V()是多了一个加权系数。它的作用在于，通过表达式（7.3.81）的获得而知，作V()极小化时再不必知道误差方差2

51、。回忆 (7.3.6)，回归函数估计的误差平方和是（7.3.82）如果想选择来使R()达极小，遇到的困难是函数g是未知的。如果2已知的话，可以通过 (7.3.10)的极小化来选择，从而极小化R(),可是如果2未知呢?那么就只有求助于广义交叉核实函数V()了。为了显示V()的作用，应该区分两种情况。如果g()m-1，这里m-1是次数等于或低于m-1的多项式，我们可以证明ER()和EV()都是对应于=。这是一种情况。第二种情况是一般情形，我们可以证明如果选择是使EV()极小，那么由下式定义的GCV方法的有效性：（7.3.83）在n时，可以趋于1。也就是说，从均方误差意义上讲，的选择已经很好了。 当

52、然，如果gm-1，则I*=1，因为ER()和EV()的极小化参数相同。如果在一般场合，gm-1但，这就要分别求出EV()的极小化参数与ER()的极小化参数*，并且满足。我们可以分几步证明I*单调递减趋于0。让我们首先考虑g ()m-1的情况。令（7.3.84）（7.3.85）（7.3.86）（7.3.87）则（7.3.88）由于g () = (g(t1), g (t2),g (tn)是矩阵的前m列的线性组合，故对一切有(I-A()g=0，b()=0，于是对ER()的极小化减弱到对trA2()的极小化，这就是对应于=的极小化。类似地，EV()变为（7.3.89）现在情况下I-A ()有m个零特征

53、值，余下的n-m个特征值可被表为n(n+vn)-1，v=1,2,，n-m, 这里1n，2n，n-m, n是n-m个正数，不完全相同，这样上式EV()就成为（7.3.90）此时当且仅当=时，EV()达到极小。 让我们再来考虑一般情况。我们有定理7.3.3 在上述记号下（7.3.91）且 （7.3.92）其中 （7.3.93）这个定理的证明主要依赖于下式（7.3.94）证明细节这里略去。从上面这个定理可以推出下面的定理： 定理7.3.4 设*是ER()的极小化参数，则EV()的极小化参数可使如下定义的（7.3.95）满足 （7.3.96）证明 令（7.3.97）因为 （7.3.98）并且ER，EV

54、以及h是的连续函数，是非空闭集。如果0不是的边界点，则EV()有极小值点在的内部(或在点)，我们称之为。由定理7.3.3，（7.3.99）如果包含0，则可能在的边界上，即=0，但是上述不等式依然成立。我们现在需要证明，当n时，h (*)0，h ()0。这可以通过以下几个引理完成。证毕引理7.3.5 如果，则 （7.3.100）引理7.3.6 若满足 （7.3.101）这里(u)是连续正值权函数，则对于gm-1 (不恒等于0的)及非零，当n时，必有b2()不等于0。引理7.3.7 设满足上一引理之假设，0(t)，则 （7.3.102） （7.3.103）这里 （7.3.104）且 （7.3.10

55、5）相反，如果n1/2m异于0，则1()与2()也都异于零，但从这个引理知道。 从上面三个引理知道，如果，则当0时，n1/2m，并且 （7.3.106）并且如果或异于0，则ER()并不趋于0。因此，为了极小化ER()，我们必须使*0，n(*)1/2m，以使h(*)0。现在可以验证EV()2，进一步有EV ()单调递减趋于2 (因为EV()-2ER (*) (1+h (*)0)。如果，则为了达到EV ()单调递减趋于2，必须有0，n ()1/2m，这样当n时，h ()0。 在上述定理与引理的基础上，我们可以得到关于广义交叉核实的主要结果。定理7.3.5 设，序列满足，其中(u)是严格正值递减数，

56、则存在序列（n），它使EV()极小化，且 （7.3.107） 算例7.3.3 样条回归与散乱资料插值设有主信号f (t)混杂了白噪声(t)，即是非参数回归模型y = f (t)+我们要获得f (t)的估计，也就是从y的资料里滤掉噪音。 下面的计算程序采用三次样条函数插值，五点平滑公式作平滑。计算过程中有提示，结果有图像，演习几次就知道样条回归的意义与存在的问题了。程序有配套的信号发生器。非参数回归的三个算例都用C语言编成，以及时显示图像，因为非参数回归的效果与控制更依赖于图像了。-样条回归与散乱资料插值计算程序, 例 7.3.3 请作出工作选择 : （1） 0 = 数据文件准备; 1 = 样条

57、回归与散乱资料插值. 请选择要读入的数据文件名: （0） 1=C21.D, 2=C22.D, 3=C23.D, 4=C24.D, 5=C25.D, 6=C11.D, 7=C12.D, 8=C13.D, 9=C14.D, 10=C15.D, 0=LI733.D, 11=自选本程序读入文本文件, 如果是二进制文件, 请转换之请输入资料长度(观测点数, li733.d是500-550) N: （500）要屏幕显示读入的资料吗 ? 0= 不显示; 1= 显示 （0）请输入样条回归平滑参数 PH (1PH=50): （5）平滑参数 PH 越大, 平滑越厉害; 但可能出现过头, 尤其在两头已经作完初始样条

58、平滑插值已经作完第 1 次样条平滑插值已经作完第 2 次样条平滑插值已经作完第 3 次样条平滑插值已经作完第 4 次样条平滑插值要屏幕显示样条平滑拟合数据吗 ? 0= 不显示; 1= 显示 （1） x(i) y(i) 0.0000 1.9412 1.0000 4.1522 2.0000 6.4290 497.0000 -1.4977 498.0000 2.1833 499.0000 4.8845 t(i) z(i) 0.5000 6.9375 1.5000 8.9958 2.5000 10.3448 489.5000 -5.5044 490.5000 -4.4215491.5000 -8.20

59、43数据存入文件吗？ 0=不存， 1=存 （0）计算结束。-下面显示原始数据与样条拟合结果的图像。图7.3.3.1第四节 非参数回归模型的小波拟合 小波理论可以为非参数回归提供迄今最新的强有力工具。这方面主要是应用它的逼近理论与信号分析理论。前面讲过，非参数回归是将混杂有信号g (t)与噪音(t)的资料Y(t)进行信噪分离。样条回归的思路是在“尽量通过每一点”与“尽量平滑”之间取得折衷。小波回归的思路则是将数据Y(t)接不同频道分解，然后重构信号时舍去噪音频道，就可得到光滑的信号曲线。所以小波回归既可以看作是信噪分离，也可以看作是一种滤波，如果主信号与噪音信号频率相差较大，小波回归比样条回归效

60、果要好得多，这从本章算例图像比较即知。读者可以先不看数学推导，而先看本节算例里的几幅图像，就明白了小波拟合的原理。 不过作为一本书，我们还是先从小波理论起源，与信噪分离有关的小波分析理论准备讲起，搭起一个简明的数学框架，然后编程计算。 一、与信噪分离有关的小波理论准备 我们先从信号函数f (t)的Fourier变换与反变换存在的问题谈起。Fourier分析的本质在于将一个相当任意的函数f (t)表示为一族标准函数的加权求和：（7.4.1）其中权函数 （7.4.2）这两个式子就是经典的付氏变换与反变换，是一种纯频域分析。它的固有缺点是在时域没有任何分辨，变换在任何有限频段上的信息都不足以确定在任