非参数统计模型.docx

1、非参数统计模型非参数统计第二次作业局部多项式回归与样条回归习题一：一、本题是研究加拿大工人收入情况，即年龄(age)和收入(income)的关系。 此次共调查了 205 个加拿大工人的年龄和收入， 所有工人都是高中毕业。 且 本题设定因变量为 log. in come,协变量为 age，运用统计方法来拟合log.income 与age之间的函数关系。、模型的建立1. 估计方法的选取 拟合两个变量之间的函数关系，即因变量和协变量之间的关系，用回归 估计的方法，回归估计包括参数回归估计和非参数回归估计。参数估计是先 假定某种数学模型或已知总体的分布，例如总体服从正态分布，其中某些参 数未知，如总体

2、均值、方差等，然后利用样本去估计这些未知参数，常用的 方法有极大似然估计，Bayes估计等，线性模型可以用最小二乘法估计。 非参数估计是不假定具有某种特定的数学模型，或总体分布未知，直接利用 样本去估计总体的数学模型，常用的方法有局部多项式回归方法和样条函数 回归方法。本题调查了 205 个加拿大工人的年龄和收入，但是加拿大工人年龄和收 入的具体分布未知，即这两个变量所能建立的数学模型未知，而且由协变量 和因变量所形成的散点图可以看出它不符合某种特定的已知模型，需要进一 步研究，然后拟合它们之间的函数关系。因此本题选用非参数回归估计的方 法，来拟合因变量和协变量之间的关系。针对此问题分别采用非

3、参数估计中的局部多项式回归和样条函数回归方 法对log.income与age之间的函数关系进行估计。2 .局部多项式回归方法局部多项式的思想是在某个点 x 附近，用一个多项式函数来逼近未知的 光滑函数 g(x) 。选定局部邻域的大小 h, 对于任意给定某个点 x 0，在其小邻 域内展开泰勒公式，用一个 p 阶多项式来局部逼近 g(x) ，然后再用极大似然 估计。(1)加拿大工人的收入(log.income )与年龄(age)之间的散点图如 下所示：InL3P2 euucluu- CTO-20 30 40 50 60age注：以下所做的图中用 X表示协变量年龄 age,用丫表示因变量收入log.

4、i ncome（2）用将X与丫排序的方法拟合的加拿大工人的收入（log.income ）与年 龄（age）之间函数关系如下图所示：20 30 40 50 60X（3）用局部多项式回归方法拟合的加拿大工人的收入 （log.income）与年龄（age）之间函数关系如下图所示：20 30 40 50 60(4)用 cross-validation 的方法选择最佳的 smoothing parameter，图形如下:3畧70n r丁68101214hvec由上图可以大概看出smoothing parameter的取值，使得函数CV.vec达到最小 的 h.vec 取值是 7, 即最佳的 smooth

5、ing parameter 取值 h=7。(5)结果分析对于最终用局部多项式回归方法拟合的收入(log.income )与年龄(age) 之间函数关系图中，黑色线条表示的是将X与丫排序拟合的函数关系；红色线条Local linear estimatel 表示的是用 Epanechnikov 核函数确疋的 smoothingparameter进行局部多项式回归得到的函数关系；蓝色线条表示用 cross-validation 方法确定的最佳 smoothing parameter 进行局部多项式回归 得到的函数关系，显然蓝色线条对 X与Y拟合的函数关系比较准确。3.样条函数回归方法样条函数的思想是

6、在区间a,b内等距离选取K个点作为节点，每两个相邻的 节点区域内都是一个基函数，且每一个基函数都是分段函数，每一组基函数构成 一个线性空间。在众多基函数选取中，B-样条基函数更稳定，应用更广泛。对于 拟合的函数的光滑程度的控制，P-Spline函数方法更好。P-Spline函数方法用 一些预先定义的节点来定义一组基函数， 同时增加一个惩罚函数，来控制拟合函数的光滑程度。然后用一组 B-样条基函数的线性组合来逼近f(x)，最后解最优 函数。(1)加拿大工人的收入(log.income )与年龄(age)之间的散点图如下所示：20 30 40 50xot)s(2)用penalized-spline

7、s 方法拟合的加拿大工人的收入(log.income )与年龄(age)之间函数关系如下图所示：seq(-5T 5,1)由上图可以大概看出 smoothi ng parameter 的取值，最佳的 smoothi ngparameter 取值 h=0.035。(4)结果分析上图中红色线条表示的是用 gen eralized cross-validati on 方法选择的最佳 smoothing parameter 进行 penalized-splines 回归得到的 X 与 Y 的函数关而且系，显然此回归结果与局部多项式回归中蓝色线条所代表的拟合函数相似，都充分凸显了散点图中xobs与yobs

8、函数关系的双峰效果，拟合程度较好习题二NOx协变量为E,运用统计一、本题是对 ethanol 数据集进行研究，因变量为 方法来拟合E与NOx之间的函数关系。、模型的建立1. 估计方法的选取拟合两个变量之间的函数关系，即因变量和协变量之间的关系，用回归 估计的方法，回归估计包括参数回归估计和非参数回归估计。参数估计是先 假定某种数学模型或已知总体的分布，例如总体服从正态分布，其中某些参 数未知，如总体均值、方差等，然后利用样本去估计这些未知参数，常用的 方法有极大似然估计，Bayes估计等，线性模型可以用最小二乘法估计。 非参数估计是不假定具有某种特定的数学模型，或总体分布未知，直接利用 样本去

9、估计数学模型，常用的方法有局部多项式回归方法，和样条函数回归 方法。本题是针对 ethanol 数据集进行研究，但是 ethanol 数据集的具体分布 未知，而且由协变量和因变量所形成的散点图可以看出它不符合某种特定的 已知模型，需要进一步研究，然后拟合它们之间的函数关系。因此本题选用 非参数回归估计的方法，来拟合因变量和协变量之间的关系。针对此问题分别采用非参数估计中的局部多项式回归和样条函数回归方 法对NOx与E之间的函数关系进行估计。1. 局部多项式回归方法注：以下所绘的图中用X表示协变E,用丫表示因变量NOx(1)ethanol数据集中NOx与E之间的函数关系散点图如下所示:寸(2)用

10、将X与丫排序的方法拟合协变量E与因变量NOx之间函数关系如下 图所示：寸AOJ0.6 0.7 0.8 0.9 1.0 1.1 1.2X(3)用局部多项式回归方法拟合的协变量 E与因变量NOx之间函数关系, 如下图所示：(4)GLLO0.02 0.03 0.04 0.05 0.06h.vec由上图可以大概看出smoothing parameter 的取值，使得函数CV.vec达至U 最小的 h.vec 取值是 0.035，即最佳的 smoothing parameter 取值 h=0.035。(5)结果分析对于最终用局部多项式回归方法拟合的协变量 E与因变量N0)之间函数关系图中，黑色线条表示的

11、是将X与丫排序拟合的函数关系；红色线条Local linear estimatel 表示的是用 Epanechnikov核函数确定的 smoothing parameter 进行 局部多项式回归得到的函数关系；蓝色线条表示用 cross-validati on 方法确定最佳的smoothing parameter进行局部多项式回归得到的函数关系， 显然蓝色线条对X与丫拟合的函数关系比较准确。2.样条函数回归方法注：以下所绘的图中用xobs表示协变E，用yobs表示因变量NOx(1)ethanol数据集中NOx与E之间的函数关系散点图如下所示：寸rt _gA CM (2)用penalized-s

12、plines 方法拟合的ethanol数据集中NOx与E之间的函数关系如下图所示：xobsseq(-103-1T 1)由上图可以大概看出smoothing parameter的取值，使得函数GCVi到最小 的横坐标取值是-6，即最佳的smoothing parameter 取值h=-6。（ 4）结果分析方法选择的最xobs 与 yobs 的上图中红色线条表示的是用 generalized cross-validation佳 smoothing parameter 进行 penalized-splines 回归得到的 函数关系。代码: 习题一：局部多项式回归 library(SemiPar) d

13、ata(age.income); X-age.income$age; Y-age.income$log.income;X2=XA2; X3=XA3; X4=XA4;fit1 - lm(YX+X2+X3+X4); coefE=c(fit1$coeff); resids=fit1$residuals; sigmaE=sqrt(var(resids);CK=1.719 temp=cbind(2,3*2*X,4*3*XA2)%*%as.vector(coefE-(1:2); den=sum(tempA2);h.ROT=CK*(sigmaEA2/den)A(1/(2*1+3);h.vec=seq(5,1

14、5,by=0.05); CV.vec=0*h.vec;for(k in 1:length(h.vec)print(k);CV.veck - CV1.fun(X,Y,h=h.veck); plot(h.vec,CV.vec,type=l);h.CV=h.vecwhich.min(CV.vec);xfine=seq(20,60,length=50);ypred1 - rep(0,length(xfine);ypred2 - rep(0,length(xfine);for(i in 1:length(xfine)ypred1i - LLS.fun(xfinei,X,Y,h=h.ROT);ypred2

15、i - LLS.fun(xfinei,X,Y,h=h.CV);plot(X,Y) lines(sort(X),sort(Y);lines(xfine,ypred1,lty=2,col=2); lines(xfine,ypred2,lty=4,col=4);linearlegend(40,12,c(True,Local linear estimate1,Local estimate2),lty=c(1,2,4),col=c(1,2,4)样条回归： library(SemiPar) data(age.income);xobs = age.income$age;yobs = age.income$l

16、og.income;nobs = length(yobs); plot(xobs,yobs);library(fda); knots=seq(min(xobs),max(xobs),length=15);nknots = length(knots);norder = 4;nbasis = length(knots) + norder - 2; basis create.bspline.basis(c(min(xobs),max(xobs),nbasis,norder,knots); basismat = eval.basis(xobs, basis);h - 0.1quadpts - seq(

17、min(xobs),max(xobs),h)nquadpts - length(quadpts)quadwts - c(1,rep(c(4,2),(nquadpts-1)/2) quadwtsnquadpts - 1quadwts - quadwts*h/3Q2basismat = eval.basis(quadpts, basis,2);Rmat = t(Q2basismat)%*%(Q2basismat*(quadwts%*%t(rep(1,nbasis) basismat2 = t(basismat)%*%basismat;lambdaVec = 10Aseq(-5,5,1)nlambd

18、a = length(lambdaVec)df = rep(0,nlambda)GCV = dffor (s in 1:nlambda)lambda = lambdaVecsBmat = basismat2 + lambda*Rmat;chat = solve(Bmat)%*%t(basismat)%*%yobs;yhat = basismat%*%chat;SSE = t(yhat-yobs)%*%(yhat-yobs)Smat = basismat%*%solve(Bmat)%*%t(basismat) dfs = sum(diag(Smat)GCVs = SSE/( nobs-dfsF2

19、 plot(seq(-5,5,1),GCV,type = l) lambda.opt = lambdaVecwhich.min(GCV);Bmat = basismat2 + lambda.opt*Rmat; chat = solve(Bmat)%*%t(basismat)%*%yobs; yhat = basismat%*%chat;plot(xobs,yobs);lines(xobs,yhat,type = l,col=red) 习题二：局部多项式回归 library(locfit);data(ethanol); X-ethanol$E Y-ethanol$NOx;X2=XA2; X3=X

20、A3; X4=XA4;fit1 -lm(YX+X2+X3+X4); coefE=c(fit1$coeff);resids=fit1$residuals; sigmaE=sqrt(var(resids);CK=1.719temp=cbind(2,3*2*X,4*3*XA2)%*%as.vector(coefE-(1:2); den=sum(tempA2);h.ROT=CK*(sigmaEA2/den)A(1/(2*1+3);h.vec=seq(0.02,0.06,by=0.0005); CV.vec=0*h.vec;for(k in 1:length(h.vec)print(k);CV.veck

21、 - CV1.fun(X,Y,h=h.veck); plot(h.vec,CV.vec,type=l);h.CV=h.vecwhich.min(CV.vec); xfine=seq(0.5,1.2,length=10); ypred1 - rep(0,length(xfine); ypred2 - rep(0,length(xfine);for(i in 1:length(xfine)ypred1i - LLS.fun(xfinei,X,Y,h=h.ROT);ypred2i - LLS.fun(xfinei,X,Y,h=h.CV); plot(X,Y) lines(sort(X),sort(Y

22、);lines(xfine,ypred1,lty=2,col=2); lines(xfine,ypred2,lty=4,col=4);linearlegend(0.8,1,c(True,Local linear estimate1,Local estimate2),lty=c(1,2,4),col=c(1,2,4)样条回归： library(locfit) data(ethanol);xobs = ethanol$E;yobs = ethanol$NOx;nobs = length(yobs); plot(xobs,yobs);library(fda); knots=seq(min(xobs)

23、,max(xobs),length=15);nknots = length(knots);norder = 4;nbasis = length(knots) + norder - 2; basis create.bspline.basis(c(min(xobs),max(xobs),nbasis,norder,knots); basismat = eval.basis(xobs, basis);h - 0.1quadpts - seq(min(xobs),max(xobs),h)nquadpts - length(quadpts)quadwts - c(1,rep(c(4,2),(nquadp

24、ts-1)/2) quadwtsnquadpts - 1quadwts - quadwts*h/3Q2basismat = eval.basis(quadpts, basis,2);Rmat = t(Q2basismat)%*%(Q2basismat*(quadwts%*%t(rep(1,nbasis) basismat2 = t(basismat)%*%basismat;lambdaVec = 10Aseq(-10,-1,1) nlambda = length(lambdaVec) df = rep(0,nlambda)GCV = dffor (s in 1:nlambda)lambda =

25、 lambdaVecsBmat = basismat2 + lambda*Rmat;chat = solve(Bmat)%*%t(basismat)%*%yobs;yhat = basismat%*%chat;SSE = t(yhat-yobs)%*%(yhat-yobs)Smat = basismat%*%solve(Bmat)%*%t(basismat) dfs = sum(diag(Smat)GCVs = SSE/( nobs-dfsF2 plot(seq(-10,-1,1),GCV,type = l) lambda.opt = lambdaVecwhich.min(GCV);Bmat = basismat2 + lambda.opt*Rmat; chat = solve(Bmat)%*%t(basismat)%*%yobs; yhat = basismat%*%chat;plot(xobs,yobs); lines(xobs,yhat,type = l,col=red)