分位数回归.

2、不同分位点拟合曲线的比较#散点图attach(engel)#打开engel数据集，直接运行其中的列名，就可以调用相应列plot(income,foodexp,cex=0.25,type=n,#画图，说明①xlab=HouseholdIncome,ylab=FoodExpenditure)points(income,foodexp,cex=0.5,col=blue)#添加点，点的大小为0.5abline(rq(foodexp~income,tau=0.5),col=blue)#画中位数回归的拟合直线，颜色蓝abline(lm(foodexp~income),lty=2,col=red)#画普通最小二乘法拟合直线，颜色红taus=c(0.05,0.1,0.25,0.75,0.9,0.95)for(iin1:length(taus)){#绘制不同分位点下的拟合直线，颜色为灰色abline(rq(foodexp~income,tau=taus[i]),col=gray)}detach(engel)3、穷人和富人的消费分布比较#比较穷人（收入在10%分位点的那个人）和富人（收入在90%分位点的那个人）的估计结果#rq函数中，tau不在[0,1]时，表示按最细的分位点划分方式得到分位点序列z=rq(foodexp~income,tau=-1)z$sol#这里包含了每个分位点下的系数估计结果x.poor=quantile(income,0.1)#10%分位点的收入x.rich=quantile(income,0.9)#90%分位点的收入ps=z$sol[1,]#每个分位点的tau值qs.poor=c(c(1,x.poor)%*%z$sol[4:5,])#10%分位点的收入的消费估计值qs.rich=c(c(1,x.rich)%*%z$sol[4:5,])#90%分位点的收入的消费估计值windows(,10,5)par(mfrow=c(1,2))#把绘图区域划分为一行两列plot(c(ps,ps),c(qs.poor,qs.rich),type=n,#type=”n”表示初始化图形区域，但不画图xlab=expression(tau),ylab=quantile)plot(stepfun(ps,c(qs.poor[1],qs.poor)),do.points=F,add=T)plot(stepfun(ps,c(qs.poor[1],qs.rich)),do.points=F,add=T,col.hor=gray,col.vert=gray)ps.wts=(c(0,diff(ps))+c(diff(ps),0))/2ap=akj(qs.poor,z=qs.poor,p=ps.wts)ar=akj(qs.rich,z=qs.rich,p=ps.wts)plot(c(qs.poor,qs.rich),c(ap$dens,ar$dens),type=n,xlab=FoodExpenditure,ylab=Density)lines(qs.rich,ar$dens,col=gray)lines(qs.poor,ap$dens,col=black)legend(topright,c(poor,rich),lty=c(1,1),col=c(black,gray))上图表示收入（income）为10%分位点处（poor，穷人）和90%分位点处（rich，富人）的食品支出的比较。从左图可以发现，对于穷人而言，在不同分位点估计的食品消费差别不大。而对于富人而言，在不同分位点对食品消费的差别比较大。右图反应了穷人和富人的食品消费分布曲线。穷人的食品消费集中于400左右，比较陡峭；而富人的消费支出集中于800到1200之间，比较分散。（四）模型比较#比较不同分位点下，收入对食品支出的影响机制是否相同fit1=rq(foodexp~income,tau=0.25)fit2=rq(foodexp~income,tau=0.5)fit3=rq(foodexp~income,tau=0.75)anova(fit1,fit2,fit3)结果：QuantileRegressionAnalysisofDevianceTableModel:foodexp~incomeJointTestofEqualityofSlopes:tauin{0.250.50.75}DfResidDfFvaluePr(F)1270315.5572.449e-07***---Signif.codes:0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1其中P值远小于0.05，故不同分位点下收入对食品支出的影响机制不同。（五）残差形态的检验也可以理解为是比较不同分位点的模型之间的关系。主要有两种模型形式：（1）位置漂移模型：不同分位点的估计结果之间的斜率相同或近似，只是截距不同；表现为不同分位点下的拟合曲线是平行的。（2）位置-尺度漂移模型：不同分位点的估计结果之间的斜率和截距都不同；表现为不同分位点下的拟合曲线不是平行的。#残差形态的检验source(C:/ProgramFiles/R/R-2.15.0/library/quantreg/doc/gasprice.R)x=gaspricen=length(x)p=5X=cbind(x[(p-1):(n-1)],x[(p-2):(n-2)],x[(p-3):(n-3)],x[(p-4):(n-4)])y=x[p:n]#位置漂移模型的检验T1=KhmaladzeTest(y~X,taus=-1,nullH=location)T2=KhmaladzeTest(y~X,taus=10:290/300,nullH=location,se=ker)#位置尺度漂移模型的检验T3=KhmaladzeTest(y~X,taus=-1,nullH=location-scale)T4=KhmaladzeTest(y~X,taus=10:290/300,nullH=location-scale,se=ker)结果：运行T1，可以查看其检验结果。其中nullH表示原假设为“location”，即原假设为位置漂移模型。Tn表示模型整体的检验，统计量为4.8。THn是对每个自变量的检验。比较T1和T3的结果（T3的原假设为“位置尺度漂移模型”），T1的统计量大于T3的统计量，可见相对而言，拒绝“位置漂移模型”的概率更大，故相对而言“位置尺度漂移模型”更加合适一些。T1$nullH[1]location$Tn[1]4.803762$THnX1X2X3X41.00031990.53216930.50208340.8926828attr(,class)[1]KhmaladzeTestT3$nullH[1]location-scale$Tn[1]2.705583$THnX1X2X3X41.21028990.69317850.50451630.8957127attr(,class)[1]KhmaladzeTest（六）非线性分位数回归这里的非线性函数为Frankcopula函数。##DemoofnonlinearquantileregressionmodelbasedonFrankcopulavFrank-function(x,df,delta,u)#某个非线性过程，得到的是[0,1]的值-log(1-(1-exp(-delta))/(1+exp(-delta*pt(x,df))*((1/u)-1)))/delta#非线性模型FrankModel-function(x,delta,mu,sigma,df,tau){z-qt(vFrank(x,df,delta,u=tau),df)mu+sigma*z}n-200#样本量df-8#自由度delta-8#初始参数set.seed(1989)x-sort(rt(n,df))#生成基于T分布的随机数v-vFrank(x,df,delta,u=runif(n))#基于x生成理论上的非参数对应值y-qt(v,df)#v对应的T分布统计量windows(5,5)plot(x,y,pch=o,col=blue,cex=.25)#散点图Dat-data.frame(x=x,y=y)#基本数据集us-c(.25,.5,.75)for(iin1:length(us)){v-vFrank(x,df,delta,u=us[i])lines(x,qt(v,df))#v为概率，计算每个概率对应的T分布统计量}cfMat-matrix(0,3,length(us)+1)#初始矩阵，用于保存结果的系数for(iin1:length(us)){tau-us[i]cat(tau=,format(tau),..)fit-nlrq(y~FrankModel(x,delta,mu,sigma,df=8,tau=tau),#非参数模型data=Dat,tau=tau,#data表明数据集,tau分位数回归的分位点start=list(delta=5,mu=0,sigma=1),#初始值trace=T)#每次运行后是否把结果显示出来lines(x,predict(fit,newdata=x),lty=2,col=red)#绘制预测曲线cfMat[i,1]-tau#保存分位点的值cfMat[i,2:4]-coef(fit)#保存系数到cfMat矩阵的第i行cat(\n)#如果前面把每步的结果显示出来，则每次的结果之间添加换行符}colnames(cfMat)-c(分位点,names(coef(fit)))#给保存系数的矩阵添加列名cfMat结果：拟合结果：（过程略）cfMat分位点deltamusigma[1,]0.2514.87165-0.205300410.9134657[2,]0.5016.253620.032325250.9638209[3,]0.7512.098360.119986140.9423476（七）半参数和非参数分位数回归非参数分位数回归在局部多项式的框架下操作起来更加方便。可以基于以下函数。#2-7-1半参数模型----fit-rq(y~bs(x,df=5)+z,tau=.33)#其中bs()表示按b-spline的非参数拟合#2-7-2非参数方法lprq-function(x,y,h,m=50,tau=0.5){#这是自定义的一个非参数计算函数，在其他数据下同样可以使用xx-seq(min(x),max(x),length=m)#m个监测点fv-xxdv-xxfor(iin1:length(xx)){z-x-xx[i]wx-dnorm(z/h)#核函数为正态分布,dnorm计算标准正态分布的密度值r-rq(y~z,weights=wx,tau=tau,#上面计算得到的密度值为权重ci=FALSE)fv[i]-r$coef[1]dv[i]-r$coef[2]}list(xx=xx,fv=fv,dv=dv)#输出结果}library(MASS)data(mcycle)attach(mcycle)windows(5,5)#非参数的结果一般是通过画图查看的plot(times,accel,xlab=milliseconds,ylab=acceleration)hs-c(1,2,3,4)#选择不同窗宽进行估计for(iinhs){h=hs[i]fit-lprq(times,accel,h=h,tau=0.5)#关键拟合函数lines(fit$xx,fit$fv,lty=i)}legend(45,-70,c(h=1,h=2,h=3,h=4),lty=1:length(hs))结果：#2-7-3非参数回归的另一个方法----#考察最大的跑步速度与体重的关系data(Mammals)attach(Mammals)x-log(weight)#取得自变量的值y-log(speed)#