非参数统计模型.docx
《非参数统计模型.docx》由会员分享，可在线阅读，更多相关《非参数统计模型.docx（18页珍藏版）》请在冰豆网上搜索。
非参数统计模型
非参数统计第二次作业
局部多项式回归与样条回归
习题一：
一、本题是研究加拿大工人收入情况，即年龄（age）和收入（income）的关系。
此次共调查了205个加拿大工人的年龄和收入，所有工人都是高中毕业。
且本题设定因变量为log.income,协变量为age，运用统计方法来拟合
log.income与age之间的函数关系。
、模型的建立
1.估计方法的选取拟合两个变量之间的函数关系，即因变量和协变量之间的关系，用回归估计的方法，回归估计包括参数回归估计和非参数回归估计。
参数估计是先假定某种数学模型或已知总体的分布，例如总体服从正态分布，其中某些参数未知，如总体均值、方差等，然后利用样本去估计这些未知参数，常用的方法有极大似然估计，Bayes估计等，线性模型可以用最小二乘法估计。
非参数估计是不假定具有某种特定的数学模型，或总体分布未知，直接利用样本去估计总体的数学模型，常用的方法有局部多项式回归方法和样条函数回归方法。
本题调查了205个加拿大工人的年龄和收入，但是加拿大工人年龄和收入的具体分布未知，即这两个变量所能建立的数学模型未知，而且由协变量和因变量所形成的散点图可以看出它不符合某种特定的已知模型，需要进一步研究，然后拟合它们之间的函数关系。
因此本题选用非参数回归估计的方法，来拟合因变量和协变量之间的关系。
针对此问题分别采用非参数估计中的局部多项式回归和样条函数回归方法对log.income与age之间的函数关系进行估计。
2.局部多项式回归方法
局部多项式的思想是在某个点x附近，用一个多项式函数来逼近未知的光滑函数g（x）。
选定局部邻域的大小h,对于任意给定某个点x0，在其小邻域内展开泰勒公式，用一个p阶多项式来局部逼近g（x），然后再用极大似然估计。
（1）加拿大工人的收入（log.income）与年龄（age）之间的散点图如下所示：
InL3P2euucluu-CTO-
2030405060
age
注：
以下所做的图中用X表示协变量年龄age,用丫表示因变量收入
log.income
（2）用将X与丫排序的方法拟合的加拿大工人的收入（log.income）与年龄（age）之间函数关系如下图所示：
2030405060
X
（3）用局部多项式回归方法拟合的加拿大工人的收入（log.income）与
年龄（age）之间函数关系如下图所示：
2030405060
（4）用cross-validation的方法选择最佳的smoothingparameter，图形如下:
3畧70
nr
丁
6
8
10
12
14
hvec
由上图可以大概看出smoothingparameter的取值，使得函数CV.vec达到最小的h.vec取值是7,即最佳的smoothingparameter取值h=7。
（5）结果分析
对于最终用局部多项式回归方法拟合的收入（log.income）与年龄（age）之间函数关系图中，黑色线条表示的是将X与丫排序拟合的函数关系；红色线条
Locallinearestimatel表示的是用Epanechnikov核函数确疋的smoothing
parameter进行局部多项式回归得到的函数关系；蓝色线条表示用cross-validation方法确定的最佳smoothingparameter进行局部多项式回归得到的函数关系，显然蓝色线条对X与Y拟合的函数关系比较准确。
3.样条函数回归方法
样条函数的思想是在区间［a,b］内等距离选取K个点作为节点，每两个相邻的节点区域内都是一个基函数，且每一个基函数都是分段函数，每一组基函数构成一个线性空间。
在众多基函数选取中，B-样条基函数更稳定，应用更广泛。
对于拟合的函数的光滑程度的控制，P-Spline函数方法更好。
P-Spline函数方法用一些预先定义的节点来定义一组基函数，同时增加一个惩罚函数，来控制拟合函
数的光滑程度。
然后用一组B-样条基函数的线性组合来逼近f（x），最后解最优函数。
（1）加拿大工人的收入（log.income）与年龄（age）之间的散点图如下所
示：
20304050
xot）s
（2）用penalized-splines方法拟合的加拿大工人的收入（log.income）
与年龄（age）之间函数关系如下图所示：
seq（-5T5,1）
由上图可以大概看出smoothingparameter的取值，最佳的smoothing
parameter取值h=0.035。
（4）结果分析
上图中红色线条表示的是用generalizedcross-validation方法选择的最
佳smoothingparameter进行penalized-splines回归得到的X与Y的函数关
而且
系，显然此回归结果与局部多项式回归中蓝色线条所代表的拟合函数相似，
都充分凸显了散点图中xobs与yobs函数关系的双峰效果，拟合程度较好
习题二
NOx协变量为E,运用统计
一、本题是对ethanol数据集进行研究，因变量为方法来拟合E与NOx之间的函数关系。
、模型的建立
1.估计方法的选取
拟合两个变量之间的函数关系，即因变量和协变量之间的关系，用回归估计的方法，回归估计包括参数回归估计和非参数回归估计。
参数估计是先假定某种数学模型或已知总体的分布，例如总体服从正态分布，其中某些参数未知，如总体均值、方差等，然后利用样本去估计这些未知参数，常用的方法有极大似然估计，Bayes估计等，线性模型可以用最小二乘法估计。
非参数估计是不假定具有某种特定的数学模型，或总体分布未知，直接利用样本去估计数学模型，常用的方法有局部多项式回归方法，和样条函数回归方法。
本题是针对ethanol数据集进行研究，但是ethanol数据集的具体分布未知，而且由协变量和因变量所形成的散点图可以看出它不符合某种特定的已知模型，需要进一步研究，然后拟合它们之间的函数关系。
因此本题选用非参数回归估计的方法，来拟合因变量和协变量之间的关系。
针对此问题分别采用非参数估计中的局部多项式回归和样条函数回归方法对NOx与E之间的函数关系进行估计。
1.局部多项式回归方法
注：
以下所绘的图中用X表示协变E,用丫表示因变量NOx
（1）ethanol数据集中NOx与E之间的函数关系散点图如下所示:
寸
（2）用将X与丫排序的方法拟合协变量E与因变量NOx之间函数关系如下图所示：
寸
A
OJ
0.60.70.80.91.01.11.2
X
（3）用局部多项式回归方法拟合的协变量E与因变量NOx之间函数关系,如下图所示：
（4）
GLLO
0.020.030.040.050.06
h.vec
由上图可以大概看出smoothingparameter的取值，使得函数CV.vec达至U最小的h.vec取值是0.035，即最佳的smoothingparameter取值h=0.035。
（5）结果分析
对于最终用局部多项式回归方法拟合的协变量E与因变量N0）之间函数关系
图中，黑色线条表示的是将X与丫排序拟合的函数关系；红色线条Locallinearestimatel表示的是用Epanechnikov核函数确定的smoothingparameter进行局部多项式回归得到的函数关系；蓝色线条表示用cross-validation方法确定
最佳的smoothingparameter进行局部多项式回归得到的函数关系，显然蓝色线
条对X与丫拟合的函数关系比较准确。
2.样条函数回归方法
注：
以下所绘的图中用xobs表示协变E，用yobs表示因变量NOx
（1）ethanol数据集中NOx与E之间的函数关系散点图如下所示：
寸
rt—
_g
ACM—
（2）用penalized-splines方法拟合的ethanol数据集中NOx与E之间
的函数关系如下图所示：
xobs
seq（-103-1T1）
由上图可以大概看出smoothingparameter的取值，使得函数GCVi到最小的横坐标取值是-6，即最佳的smoothingparameter取值h=-6。
（4）结果分析
方法选择的最
xobs与yobs的
上图中红色线条表示的是用generalizedcross-validation
佳smoothingparameter进行penalized-splines回归得到的函数关系。
代码:
习题一：
局部多项式回归library（SemiPar）data（age.income）;X<-age.income$age;Y<-age.income$log.income;
X2=XA2;X3=XA3;X4=XA4;
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*XA2）%*%as.vector（coefE[-（1:
2）]）;den=sum（tempA2）;
h.ROT=CK*（sigmaEA2/den）A（1/（2*1+3））;
h.vec=seq（5,15,by=0.05）;CV.vec=0*h.vec;
for（kin1:
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（iin1:
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）;
linear
legend（40,12,c（"True","Locallinearestimate1","Localestimate2"）,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;basiscreate.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=10Aseq（-5,5,1）
nlambda=length（lambdaVec）
df=rep（0,nlambda）
GCV=df
for（sin1:
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]F2
}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=XA2;X3=XA3;X4=XA4;
fit1<-lm（Y~X+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（0.02,0.06,by=0.0005）;CV.vec=0*h.vec;
for（kin1:
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（iin1:
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）;
linear
legend（0.8,1,c（"True","Locallinearestimate1","Localestimate2"）,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;basiscreate.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=10Aseq（-10,-1,1）nlambda=length（lambdaVec）df=rep（0,nlambda）
GCV=df
for（sin1:
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]F2
}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"）
升级会员 