用R语言进行分位数回归Word文档格式.docx

1、98/100） ）fit2 = summary（ rq（foodexp income, tau = c,） ）windows（5,5） # 新建一个图形窗口，可以去掉这句plot（fit1）plot（fit2）# 散点图attach（engel） # 打开engel数据集，直接运行其中的列名，就可以调用相应列plot（income,foodexp,cex=,type=n, # 画图，说明 xlab=Household Income, ylab=Food Expenditure）points（income,foodexp,cex=,col=blue）# 添加点，点的大小为abline（ rq（f

2、oodexp income, tau=, col= ）# 画中位数回归的拟合直线，颜色蓝abline（ lm（foodexp income）, lty = 2, col=red# 画普通最小二乘法拟合直线，颜色红taus = c, , , , , for（i in 1:length（taus） # 绘制不同分位点下的拟合直线，颜色为灰色abline（ rq（foodexp income, tau=tausi）, col=gray ）detach（engel）# 比较穷人（收入在10%分位点的那个人）和富人（收入在90%分位点的那个人）的估计结果# rq函数中，tau不在0,1时，表示按最细的分

3、位点划分方式得到分位点序列z = rq（foodexp income, tau=-1）z$sol # 这里包含了每个分位点下的系数估计结果 = quantile（income, # 10%分位点的收入 = quantile（income, # 90%分位点的收入ps = z$sol1, # 每个分位点的tau值 = c（ c（1, %*% z$sol4:5, ）# 10%分位点的收入的消费估计值 # 90%分位点的收入的消费估计值windows（10,5）par（mfrow=c（1,2） # 把绘图区域划分为一行两列plot（c（ps,ps）,c,type= # type=”n”表示初始化图形

4、区域，但不画图 xlab=expression（tau）, ylab=quantileplot（stepfun（ps,c1,）, =F, add=T）plot（stepfun（ps,c1,）, =F, add=T, =, = = （ c（0,diff（ps） + c（diff（ps）,0） ）/2ap = akj, z=, p=ar = akj, z=, p=plot（c, c（ap$dens, ar$dens）, type=, xlab=Densitylines,ar$dens,col=lines,ap$dens,col=blacklegend（topright, c（poor,rich）,

5、 lty=c（1,1）, col=c（）# 比较不同分位点下，收入对食品支出的影响机制是否相同fit1 = rq（foodexp income, tau = fit2 = rq（foodexp income, tau = fit3 = rq（foodexp income, tau = anova（fit1,fit2,fit3）# 残差形态的检验source（C:/Program Files/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 = xp:n#

6、位置漂移模型的检验T1 = KhmaladzeTest（yX, taus = -1, nullH=location） T2 = KhmaladzeTest（yX, taus = 10:290/300, nullH=, se=ker# 位置尺度漂移模型的检验T3 = KhmaladzeTest（yX, taus = -1, nullH=location-scaleT4 = KhmaladzeTest（yX, taus = 10:# Demo of nonlinear quantile regression model based on Frank copulavFrank - function（

7、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*zn - 200# 样本量df - 8 # 自由度delta - 8 # 初始参数（1989） x - sort（rt（n,df）# 生成基于T分布的随机数v - vFrank（x, df,

8、delta, u = runif（n）# 基于x生成理论上的非参数对应值y - qt（v, df） # v 对应的T分布统计量windows（5,5）plot（x, y, pch=o, col=, cex = .25） # 散点图Dat - （x = x, y = y） # 基本数据集us - c（.25,.5,.75）length（us）- vFrank（x, df, delta, u = usi）lines（x, qt（v,df）# v为概率，计算每个概率对应的T分布统计量cfMat - matrix（0, 3, length（us）+1） # 初始矩阵，用于保存结果的系数length（u

9、s） tau - usicat（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=# 绘制预测曲线cfMati,1 - tau #

10、 保存分位点的值cfMati,2:4 - coef（fit） # 保存系数到cfMat矩阵的第i行n # 如果前面把每步的结果显示出来，则每次的结果之间添加换行符colnames（cfMat） - c（分位点,names（coef（fit）# 给保存系数的矩阵添加列名cfMat# 2-7-1 半参数模型-fit-rq（ybs（x,df=5）+z,tau=.33）# 其中bs（）表示按b-spline的非参数拟合# 2-7-2 非参数方法lprq -function（x,y,h,m=50,tau=# 这是自定义的一个非参数计算函数，在其他数据下同样可以使用xx-seq（min（x）,max（x）

11、,length=m） # m个监测点fv-xxdvlength（xx） z-x-xxi wx-dnorm（z/h） # 核函数为正态分布,dnorm计算标准正态分布的密度值 r-rq（yz,weights=wx,tau=tau, # 上面计算得到的密度值为权重 ci=FALSE） fvi-r$coef1 dvi-r$coef2list（xx=xx,fv=fv,dv=dv）# 输出结果library（MASS）data（mcycle）attach（mcycle） # 非参数的结果一般是通过画图查看的plot（times,accel,xlab=milliseconds,ylab=accelerat

12、ionhs-c（1,2,3,4） # 选择不同窗宽进行估计for（i in hs）h=hsi-lprq（times,accel,h=h,tau=# 关键拟合函数lines（fit$xx,fit$fv,lty=i）legend（45,-70,c（h=1h=2h=3h=4）, lty=1:length（hs）# 2-7-3 非参数回归的另一个方法-# 考察最大的跑步速度与体重的关系data（Mammals）attach（Mammals）x-log（weight）# 取得自变量的值y-log（speed）# 取得因变量的值plot（x,y,xlab=Weightinlog（Kg）Speedinlog

13、（Km/hour）,points（xhoppers,yhoppers,pch=h,col=points（xspecials,yspecials,pch=sothers0 & dat$foodexp0,dat2,lnincome = log（ dat2, ）lnfoodexpdat2 = （dat2）# 去掉log变换以后的缺失值行dim（dat2）# 散点图，观察是否需要log变换attach（dat）# 不变换plot（income,foodexp,cex=,main=散点图detach（dat）# 变换attach（dat2）plot（lnincome,lnfoodexp,cexp=,ma

14、in=散点图（对数变换以后）detach（dat2）# 建立基础的分位数回归模型taus = c,fit1 = rq（ lnfoodexp lnincome, data=dat2, tau=taus, method=fnfit2 = rq（ lnfoodexp lnincome, data=dat2dat2$urtype=1, tau=taus, method=fit3 = rq（ lnfoodexp lnincome, data=dat2dat2$urtype=2,# 所有数据放在一起是参数结果s1 = summary（fit1）# 男性s2 = summary（fit2）# 女性s3 =

15、summary（fit3）# 参数结果的直接比较tabs（s1）# 这里的tabs函数是笔者自己写的，方便显示结果tabs（s2）tabs（s3）# 具体的系数估计值及假设检验tabcoef（s1）# 这里的tabcoef函数是笔者自己写的，方便显示结果tabcoef（s2）tabcoef（s3）Quantiles（Intercept）lnincomeValueStd. Errort valuePr（|t|）# 参数结果的检验T1 = KhmaladzeTest（lnfoodexp lnincome, taus=seq（.1,.9,by = .1）, T2 = KhmaladzeTest（lnfoo