非參數(shù)統(tǒng)計(jì)模型(共15頁(yè))

1、非參數(shù)統(tǒng)計(jì)(tngj)第二次作業(yè) 局部(jb)多項(xiàng)式回歸與樣條回歸習(xí)題(xt)一：本題是研究加拿大工人收入情況，即年齡（age）和收入(income)的關(guān)系。此次共調(diào)查了205個(gè)加拿大工人的年齡和收入，所有工人都是高中畢業(yè)。且本題設(shè)定因變量為log.income,協(xié)變量為age，運(yùn)用統(tǒng)計(jì)方法來(lái)擬合log.income 與age之間的函數(shù)關(guān)系。模型的建立1.估計(jì)方法的選取擬合兩個(gè)變量之間的函數(shù)關(guān)系，即因變量和協(xié)變量之間的關(guān)系，用回歸估計(jì)的方法，回歸估計(jì)包括參數(shù)回歸估計(jì)和非參數(shù)回歸估計(jì)。參數(shù)估計(jì)是先假定某種數(shù)學(xué)模型或已知總體的分布，例如總體服從正態(tài)分布，其中某些參數(shù)未知，如總體均值、方差等，然后利

2、用樣本去估計(jì)這些未知參數(shù)，常用的方法有極大似然估計(jì)，估計(jì)等，線性模型可以用最小二乘法估計(jì)。非參數(shù)估計(jì)是不假定具有某種特定的數(shù)學(xué)模型，或總體分布未知，直接利用樣本去估計(jì)總體的數(shù)學(xué)模型，常用的方法有局部多項(xiàng)式回歸方法和樣條函數(shù)回歸方法。本題調(diào)查了205個(gè)加拿大工人的年齡和收入，但是加拿大工人年齡和收入的具體分布未知，即這兩個(gè)變量所能建立的數(shù)學(xué)模型未知，而且由協(xié)變量和因變量所形成的散點(diǎn)圖可以看出它不符合某種特定的已知模型，需要進(jìn)一步研究，然后擬合它們之間的函數(shù)關(guān)系。因此本題選用非參數(shù)回歸估計(jì)的方法，來(lái)擬合因變量和協(xié)變量之間的關(guān)系。針對(duì)此問(wèn)題分別采用非參數(shù)估計(jì)中的局部多項(xiàng)式回歸和樣條函數(shù)回歸方法對(duì)lo

3、g.income 與age之間的函數(shù)關(guān)系進(jìn)行估計(jì)。.局部多項(xiàng)式回歸方法局部多項(xiàng)式的思想是在某個(gè)點(diǎn)x附近，用一個(gè)多項(xiàng)式函數(shù)來(lái)逼近未知的光滑函數(shù)g(x)。選定局部鄰域的大小h,對(duì)于任意給定某個(gè)點(diǎn)x 0，在其小鄰域內(nèi)展開(kāi)泰勒公式，用一個(gè)p階多項(xiàng)式來(lái)局部逼近g(x)，然后再用極大似然估計(jì)。（1）加拿大工人的收入（log.income）與年齡（age）之間的散點(diǎn)圖如下所示：注：以下所做的圖中用(zhngyng)X表示協(xié)變量年齡age，用Y表示因變量收入log.income（2）用將X與Y排序(pi x)的方法(fngf)擬合的加拿大工人的收入（log.income）與年齡（age）之間函數(shù)關(guān)系如下圖所示

4、： （3）用局部多項(xiàng)式回歸方法擬合的加拿大工人的收入（log.income）與年齡（age）之間函數(shù)關(guān)系如下圖所示： （4）用cross-validation的方法(fngf)選擇(xunz)最佳(zu ji)的smoothing parameter，圖形如下：由上圖可以大概看出smoothing parameter的取值，使得函數(shù)CV.vec達(dá)到最小的h.vec取值是7，即最佳的smoothing parameter取值h=7。（5）結(jié)果(ji gu)分析對(duì)于(duy)最終用局部多項(xiàng)式回歸方法(fngf)擬合的收入（log.income）與年齡（age）之間函數(shù)關(guān)系圖中，黑色線條表示的是將X

5、與Y排序擬合的函數(shù)關(guān)系；紅色線條Local linear estimate1表示的是用Epanechnikov核函數(shù)確定的smoothing parameter進(jìn)行局部多項(xiàng)式回歸得到的函數(shù)關(guān)系；藍(lán)色線條表示用cross-validation方法確定的最佳smoothing parameter進(jìn)行局部多項(xiàng)式回歸得到的函數(shù)關(guān)系，顯然藍(lán)色線條對(duì)X與Y擬合的函數(shù)關(guān)系比較準(zhǔn)確。3. 樣條函數(shù)回歸方法 樣條函數(shù)的思想是在區(qū)間a,b內(nèi)等距離選取K個(gè)點(diǎn)作為節(jié)點(diǎn)，每?jī)蓚€(gè)相鄰的節(jié)點(diǎn)區(qū)域內(nèi)都是一個(gè)基函數(shù)，且每一個(gè)基函數(shù)都是分段函數(shù)，每一組基函數(shù)構(gòu)成一個(gè)線性空間。在眾多基函數(shù)選取中，B-樣條基函數(shù)更穩(wěn)定，應(yīng)用更廣泛。

6、對(duì)于擬合的函數(shù)的光滑程度的控制，P-Spline函數(shù)方法更好。P-Spline函數(shù)方法用一些預(yù)先定義的節(jié)點(diǎn)來(lái)定義一組基函數(shù)，同時(shí)增加一個(gè)懲罰函數(shù)，來(lái)控制擬合函數(shù)的光滑程度。然后用一組B-樣條基函數(shù)的線性組合來(lái)逼近f(x)，最后解最優(yōu)函數(shù)。（1）加拿大工人的收入（log.income）與年齡（age）之間的散點(diǎn)圖如下所示：（2）用penalized-splines方法擬合的加拿大工人的收入（log.income）與年齡（age）之間函數(shù)關(guān)系如下圖所示： （3）用generalized cross-validation的方法(fngf)選擇(xunz)最佳(zu ji)的smoothing par

7、ameter，圖形如下：由上圖可以大概看出smoothing parameter的取值，最佳的smoothing parameter取值h=0.035。（4）結(jié)果分析 上圖中紅色線條表示的是用generalized cross-validation方法選擇的最佳smoothing parameter 進(jìn)行penalized-splines回歸得到的X與Y的函數(shù)關(guān)系，顯然此回歸結(jié)果與局部多項(xiàng)式回歸中藍(lán)色線條所代表的擬合函數(shù)相似，而且都充分凸顯了散點(diǎn)圖中xobs與yobs函數(shù)關(guān)系的雙峰效果，擬合程度較好。習(xí)題(xt)二一、本題(bnt)是對(duì)ethanol數(shù)據(jù)集進(jìn)行研究，因變量為NOx，協(xié)變量為E，

8、運(yùn)用(ynyng)統(tǒng)計(jì)方法來(lái)擬合E與NOx之間的函數(shù)關(guān)系。二、模型的建立1.估計(jì)方法的選取擬合兩個(gè)變量之間的函數(shù)關(guān)系，即因變量和協(xié)變量之間的關(guān)系，用回歸估計(jì)的方法，回歸估計(jì)包括參數(shù)回歸估計(jì)和非參數(shù)回歸估計(jì)。參數(shù)估計(jì)是先假定某種數(shù)學(xué)模型或已知總體的分布，例如總體服從正態(tài)分布，其中某些參數(shù)未知，如總體均值、方差等，然后利用樣本去估計(jì)這些未知參數(shù)，常用的方法有極大似然估計(jì)，估計(jì)等，線性模型可以用最小二乘法估計(jì)。非參數(shù)估計(jì)是不假定具有某種特定的數(shù)學(xué)模型，或總體分布未知，直接利用樣本去估計(jì)數(shù)學(xué)模型，常用的方法有局部多項(xiàng)式回歸方法，和樣條函數(shù)回歸方法。本題是針對(duì)ethanol數(shù)據(jù)集進(jìn)行研究，但是ethan

9、ol數(shù)據(jù)集的具體分布未知，而且由協(xié)變量和因變量所形成的散點(diǎn)圖可以看出它不符合某種特定的已知模型，需要進(jìn)一步研究，然后擬合它們之間的函數(shù)關(guān)系。因此本題選用非參數(shù)回歸估計(jì)的方法，來(lái)擬合因變量和協(xié)變量之間的關(guān)系。針對(duì)此問(wèn)題分別采用非參數(shù)估計(jì)中的局部多項(xiàng)式回歸和樣條函數(shù)回歸方法對(duì)NOx與E之間的函數(shù)關(guān)系進(jìn)行估計(jì)。局部多項(xiàng)式回歸方法注：以下所繪的圖中用X表示協(xié)變E，用Y表示因變量NOx。（1）ethanol數(shù)據(jù)集中NOx與E之間的函數(shù)關(guān)系散點(diǎn)圖如下所示： （2）用將X與Y排序(pi x)的方法擬合協(xié)變量E與因變量NOx之間函數(shù)(hnsh)關(guān)系如下圖所示：（3）用局部多項(xiàng)式回歸方法(fngf)擬合的協(xié)變量

10、E與因變量NOx之間函數(shù)關(guān)系，如下圖所示：（4）用cross-validation的方法(fngf)選擇(xunz)最佳(zu ji)的smoothing parameter，圖形如下：由上圖可以(ky)大概看出smoothing parameter的取值，使得(sh de)函數(shù)CV.vec達(dá)到(d do)最小的h.vec取值是0.035，即最佳的smoothing parameter取值h=0.035。（5）結(jié)果分析對(duì)于最終用局部多項(xiàng)式回歸方法擬合的協(xié)變量E與因變量NOx之間函數(shù)關(guān)系圖中，黑色線條表示的是將X與Y排序擬合的函數(shù)關(guān)系；紅色線條Local linear estimate1表示的是

11、用Epanechnikov核函數(shù)確定的smoothing parameter進(jìn)行局部多項(xiàng)式回歸得到的函數(shù)關(guān)系；藍(lán)色線條表示用cross-validation方法確定最佳的smoothing parameter進(jìn)行局部多項(xiàng)式回歸得到的函數(shù)關(guān)系，顯然藍(lán)色線條對(duì)X與Y擬合的函數(shù)關(guān)系比較準(zhǔn)確。2.樣條函數(shù)回歸方法注：以下所繪的圖中用xobs表示協(xié)變E，用yobs表示因變量NOx。（1）ethanol數(shù)據(jù)集中NOx與E之間的函數(shù)關(guān)系散點(diǎn)圖如下所示：（2）用penalized-splines方法擬合的ethanol數(shù)據(jù)集中NOx與E之間的函數(shù)關(guān)系如下圖所示： （3）用generalized cross-v

12、alidation的方法(fngf)選擇(xunz)最佳(zu ji)的smoothing parameter，圖形如下：由上圖可以大概看出smoothing parameter的取值，使得函數(shù)GCV達(dá)到最小的橫坐標(biāo)取值是-6，即最佳的smoothing parameter取值h=-6。（4）結(jié)果(ji gu)分析上圖中紅色(hngs)線條表示的是用generalized cross-validation方法(fngf)選擇的最佳smoothing parameter 進(jìn)行penalized-splines回歸得到的xobs與yobs的函數(shù)關(guān)系。代碼:習(xí)題一：局部多項(xiàng)式回歸library(Se

13、miPar)data(age.income);X-age.income$age;Y-age.income$log.income;X2=X2; X3=X3; X4=X4;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*X2)%*%as.vector(coefE-(1:2);den=sum(temp2);h.ROT=CK*(sigmaE2/den)(1/(2*1+3);h.vec=seq(5,15,by=0

14、.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); ypred2i -

15、 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);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.in

16、come;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(xob

17、s),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 = 10seq(-5,5,1)nlambda = length

18、(lambdaVec)df = rep(0,nlambda)GCV = dffor (s in 1:nlambda) lambda = lambdaVecs 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) dfs = sum(diag(Smat) GCVs = SSE/(nobs-dfs)2plot

19、(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)習(xí)題二：局部(jb)多項(xiàng)式回歸library(locfit);data(ethanol);X-ethanol$EY-ethanol$NOx;X2=X2; X3=X3; X4=X4;

20、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*X2)%*%as.vector(coefE-(1:2);den=sum(temp2);h.ROT=CK*(sigmaE2/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.veck - CV1.fun(X,Y,h

21、=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);lines(xfine,ypr

22、ed1,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),len

23、gth=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) quadwtsnquadpts - 1 quadwts - quadwts*h/3 Q2basismat = eval.basis(quadpts, basis,2); Rmat = t