非参数统计第二次作业——局部多项式回归与样条回归习题一：一、本题是研究加拿大工人收入情况，即年龄（age）和收入(income)的关系。
此次共调查了205个加拿大工人的年龄和收入，所有工人都是高中毕业。
且本题设定因变量为log.income,协变量为age，运用统计方法来拟合log.income 与age之间的函数关系。
二、模型的建立1.估计方法的选取拟合两个变量之间的函数关系，即因变量和协变量之间的关系，用回归估计的方法，回归估计包括参数回归估计和非参数回归估计。
参数估计是先假定某种数学模型或已知总体的分布，例如总体服从正态分布，其中某些参数未知，如总体均值、方差等，然后利用样本去估计这些未知参数，常用的方法有极大似然估计，Ｂａｙｅｓ估计等，线性模型可以用最小二乘法估计。
非参数估计是不假定具有某种特定的数学模型，或总体分布未知，直接利用样本去估计总体的数学模型，常用的方法有局部多项式回归方法和样条函数回归方法。
本题调查了205个加拿大工人的年龄和收入，但是加拿大工人年龄和收入的具体分布未知，即这两个变量所能建立的数学模型未知，而且由协变量和因变量所形成的散点图可以看出它不符合某种特定的已知模型，需要进一步研究，然后拟合它们之间的函数关系。
因此本题选用非参数回归估计的方法，来拟合因变量和协变量之间的关系。
针对此问题分别采用非参数估计中的局部多项式回归和样条函数回归方法对log.income 与age之间的函数关系进行估计。
２.局部多项式回归方法局部多项式的思想是在某个点x附近，用一个多项式函数来逼近未知的光滑函数g(x)。
选定局部邻域的大小h,对于任意给定某个点x 0，在其小邻域内展开泰勒公式，用一个p阶多项式来局部逼近g(x)，然后再用极大似然估计。
（1）加拿大工人的收入（log.income）与年龄（age）之间的散点图如下所示：注：以下所做的图中用X表示协变量年龄age，用Y表示因变量收入log.income（2）用将X与Y排序的方法拟合的加拿大工人的收入（log.income）与年龄（age）之间函数关系如下图所示：（3）用局部多项式回归方法拟合的加拿大工人的收入（log.income）与年龄（age）之间函数关系如下图所示：（4）用cross-validation的方法选择最佳的smoothing parameter，图形如下：由上图可以大概看出smoothing parameter的取值，使得函数CV.vec达到最小的h.vec取值是7，即最佳的smoothing parameter取值h=7。
（5）结果分析对于最终用局部多项式回归方法拟合的收入（log.income）与年龄（age）之间函数关系图中，黑色线条表示的是将X与Y排序拟合的函数关系；红色线条Local linear estimate1表示的是用Epanechnikov核函数确定的smoothing parameter进行局部多项式回归得到的函数关系；蓝色线条表示用cross-validation方法确定的最佳smoothing parameter进行局部多项式回归得到的函数关系，显然蓝色线条对X与Y拟合的函数关系比较准确。
3. 样条函数回归方法样条函数的思想是在区间[a,b]内等距离选取K个点作为节点，每两个相邻的节点区域内都是一个基函数，且每一个基函数都是分段函数，每一组基函数构成一个线性空间。
在众多基函数选取中，B-样条基函数更稳定，应用更广泛。
对于拟合的函数的光滑程度的控制，P-Spline函数方法更好。
P-Spline函数方法用一些预先定义的节点来定义一组基函数，同时增加一个惩罚函数，来控制拟合函数的光滑程度。
然后用一组B-样条基函数的线性组合来逼近f(x)，最后解最优函数。
（1）加拿大工人的收入（log.income）与年龄（age）之间的散点图如下所示：（2）用penalized-splines方法拟合的加拿大工人的收入（log.income）与年龄（age）之间函数关系如下图所示：（3）用generalized cross-validation的方法选择最佳的smoothing parameter，图形如下：由上图可以大概看出smoothing parameter的取值，最佳的smoothing parameter取值h=0.035。
（4）结果分析上图中红色线条表示的是用generalized cross-validation方法选择的最佳smoothing parameter 进行penalized-splines回归得到的X与Y的函数关系，显然此回归结果与局部多项式回归中蓝色线条所代表的拟合函数相似，而且都充分凸显了散点图中xobs与yobs函数关系的双峰效果，拟合程度较好。
习题二一、本题是对ethanol数据集进行研究，因变量为NOx，协变量为E，运用统计方法来拟合E与NOx之间的函数关系。
二、模型的建立1.估计方法的选取拟合两个变量之间的函数关系，即因变量和协变量之间的关系，用回归估计的方法，回归估计包括参数回归估计和非参数回归估计。
参数估计是先假定某种数学模型或已知总体的分布，例如总体服从正态分布，其中某些参数未知，如总体均值、方差等，然后利用样本去估计这些未知参数，常用的方法有极大似然估计，Ｂａｙｅｓ估计等，线性模型可以用最小二乘法估计。
非参数估计是不假定具有某种特定的数学模型，或总体分布未知，直接利用样本去估计数学模型，常用的方法有局部多项式回归方法，和样条函数回归方法。
本题是针对ethanol数据集进行研究，但是ethanol数据集的具体分布未知，而且由协变量和因变量所形成的散点图可以看出它不符合某种特定的已知模型，需要进一步研究，然后拟合它们之间的函数关系。
因此本题选用非参数回归估计的方法，来拟合因变量和协变量之间的关系。
针对此问题分别采用非参数估计中的局部多项式回归和样条函数回归方法对NOx与E之间的函数关系进行估计。
1.局部多项式回归方法注：以下所绘的图中用X表示协变E，用Y表示因变量NOx。
（1）ethanol数据集中NOx与E之间的函数关系散点图如下所示：（2）用将X与Y排序的方法拟合协变量E与因变量NOx之间函数关系如下图所示：（3）用局部多项式回归方法拟合的协变量E与因变量NOx之间函数关系，如下图所示：（4）用cross-validation的方法选择最佳的smoothing parameter，图形如下：由上图可以大概看出smoothing parameter的取值，使得函数CV.vec达到最小的h.vec取值是0.035，即最佳的smoothing parameter取值h=0.035。
（5）结果分析对于最终用局部多项式回归方法拟合的协变量E与因变量NOx之间函数关系图中，黑色线条表示的是将X与Y排序拟合的函数关系；红色线条Local linear estimate1表示的是用Epanechnikov核函数确定的smoothing parameter进行局部多项式回归得到的函数关系；蓝色线条表示用cross-validation方法确定最佳的smoothing parameter进行局部多项式回归得到的函数关系，显然蓝色线条对X与Y拟合的函数关系比较准确。
2.样条函数回归方法注：以下所绘的图中用xobs表示协变E，用yobs表示因变量NOx。
（1）ethanol数据集中NOx与E之间的函数关系散点图如下所示：（2）用penalized-splines方法拟合的ethanol数据集中NOx与E之间的函数关系如下图所示：（3）用generalized cross-validation的方法选择最佳的smoothing parameter，图形如下：由上图可以大概看出smoothing parameter的取值，使得函数GCV达到最小的横坐标取值是-6，即最佳的smoothing parameter取值h=-6。
（4）结果分析上图中红色线条表示的是用generalized cross-validation方法选择的最佳smoothing parameter 进行penalized-splines回归得到的xobs与yobs的函数关系。
代码:习题一：局部多项式回归library(SemiPar)data(age.income);X<-age.income$age;Y<-age.income$log.income;X2=X^2; X3=X^3; X4=X^4;fit1 <- lm(Y~X+X2+X3+X4);coefE=c(fit1$coeff);resids=fit1$residuals;sigmaE=sqrt(var(resids));CK=1.719temp=cbind(2,3*2*X,4*3*X^2)%*%as.vector(coefE[-(1:2)]);den=sum(temp^2);h.ROT=CK*(sigmaE^2/den)^(1/(2*1+3));h.vec=seq(5,15,by=0.05);CV.vec=0*h.vec;for(k in 1:length(h.vec)){print(k);CV.vec[k] <- CV1.fun(X,Y,h=h.vec[k]);}plot(h.vec,CV.vec,type="l");h.CV=h.vec[which.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)){ypred1[i] <- LLS.fun(xfine[i],X,Y,h=h.ROT);ypred2[i] <- LLS.fun(xfine[i],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);legend(40,12,c("True","Local linear estimate1","Local linear estimate2"),lty=c(1,2,4),col=c(1,2,4))样条回归：library(SemiPar)data(age.income);xobs = age.income$age;yobs = age.income$log.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(min(xobs),max(xobs),h)nquadpts <- length(quadpts)quadwts <- c(1,rep(c(4,2),(nquadpts-1)/2))quadwts[nquadpts] <- 1quadwts <- quadwts*h/3Q2basismat = eval.basis(quadpts, basis,2);Rmat = t(Q2basismat)%*%(Q2basismat*(quadwts%*%t(rep(1,nbasis)))) basismat2 = t(basismat)%*%basismat;lambdaVec = 10^seq(-5,5,1)nlambda = length(lambdaVec)df = rep(0,nlambda)GCV = dffor (s in 1:nlambda){lambda = lambdaVec[s]Bmat = basismat2 + lambda*Rmat;chat = solve(Bmat)%*%t(basismat)%*%yobs;yhat = basismat%*%chat;SSE = t(yhat-yobs)%*%(yhat-yobs)Smat = basismat%*%solve(Bmat)%*%t(basismat)df[s] = sum(diag(Smat))GCV[s] = SSE/(nobs-df[s])^2}plot(seq(-5,5,1),GCV,type = "l")lambda.opt = lambdaVec[which.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$EY<-ethanol$NOx;X2=X^2; X3=X^3; X4=X^4;fit1 <-lm(Y~X+X2+X3+X4);coefE=c(fit1$coeff);resids=fit1$residuals;sigmaE=sqrt(var(resids));CK=1.719temp=cbind(2,3*2*X,4*3*X^2)%*%as.vector(coefE[-(1:2)]); den=sum(temp^2);h.ROT=CK*(sigmaE^2/den)^(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.vec[k] <- CV1.fun(X,Y,h=h.vec[k]);}plot(h.vec,CV.vec,type="l");h.CV=h.vec[which.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)){ypred1[i] <- LLS.fun(xfine[i],X,Y,h=h.ROT);ypred2[i] <- LLS.fun(xfine[i],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);legend(0.8,1,c("True","Local linear estimate1","Local linear 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),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.1 quadpts <- seq(min(xobs),max(xobs),h) nquadpts <- length(quadpts) quadwts <- c(1,rep(c(4,2),(nquadpts-1)/2)) quadwts[nquadpts] <- 1 quadwts <- quadwts*h/3 Q2basismat = eval.basis(quadpts, basis,2); Rmat = t(Q2basismat)%*%(Q2basismat*(quadwts%*%t(rep(1,nbasis)))) basismat2 = t(basismat)%*%basismat; lambdaVec = 10^seq(-10,-1,1) nlambda = length(lambdaVec) df = rep(0,nlambda) GCV = df for (s in 1:nlambda) { lambda = lambdaVec[s] Bmat = basismat2 + lambda*Rmat; chat = solve(Bmat)%*%t(basismat)%*%yobs;yhat = basismat%*%chat; SSE = t(yhat-yobs)%*%(yhat-yobs) Smat = basismat%*%solve(Bmat)%*%t(basismat) df[s] = sum(diag(Smat)) GCV[s] = SSE/(nobs-df[s])^2 } plot(seq(-10,-1,1),GCV,type = "l") lambda.opt = lambdaVec[which.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")。