楼主: 耕耘使者
3947 3

[问答] R的分位回归中,分位点的选择有限制?? [推广有奖]

贵宾

已卖:5006份资源

学术权威

39%

还不是VIP/贵宾

-

威望
4
论坛币
1811021 个
通用积分
171.7911
学术水平
109 点
热心指数
173 点
信用等级
87 点
经验
93428 点
帖子
4549
精华
0
在线时间
2848 小时
注册时间
2006-4-6
最后登录
2025-7-12

楼主
耕耘使者 发表于 2012-5-5 10:05:33 |AI写论文
200论坛币
y=c(2167,2509,3530,4669,5868,6763,6820,2866,8248,8868,9336,10464,11040,12631,13773,14762,17255,19398,20992,23200)
income=c(2486,3009,4277,5868,7172,8159,8439,8773,10932,11718,12883,13250,14867,16683,18645,20668,23623,26675,28838,21754)
cpi=c(110.5,110.0,120.2,123.9,118.7,109.2,102.8,100.0,101.5,102.5,100.0,100.5,100.1,102.2,101.0,101.2,103.2,105.8,99.6,103.1)
library(quantreg)
fit1=rq(y~income+cpi,tau=c(0.2,0.5,0.8,0.99))
summary(fit1,se="nid")

提示:
错误于summary.rq(xi, ...) : tau - h < 0:  error in summary.rq>


最佳答案

ywh19860616 查看完整内容

使者兄,根据软件提示,应该是带宽h-tau
关键词:分位回归 Summary Library Income Summa cpi library income

沙发
ywh19860616 发表于 2012-5-5 10:05:34
使者兄,根据软件提示,应该是带宽h-tau<0所导致的。summary(fit1)命令其实是通过呼叫summary.rq,然后看下summary.rq内部codes:
else if (se == "nid") {
        h <- bandwidth.rq(tau, n, hs = hs)
        if (tau + h > 1)
            stop("tau + h > 1:  error in summary.rq")
        if (tau - h < 0)
            stop("tau - h < 0:  error in summary.rq")
        bhi <- rq.fit.fnb(x, y, tau = tau + h)$coef
        blo <- rq.fit.fnb(x, y, tau = tau - h)$coef
        dyhat <- x %*% (bhi - blo)
        if (any(dyhat <= 0))
            warning(paste(sum(dyhat <= 0), "non-positive fis"))
        f <- pmax(0, (2 * h)/(dyhat - eps))
        fxxinv <- diag(p)
        fxxinv <- backsolve(qr(sqrt(f) * x)$qr[1:p, 1:p, drop = FALSE],
            fxxinv)
        fxxinv <- fxxinv %*% t(fxxinv)
        cov <- tau * (1 - tau) * fxxinv %*% crossprod(x) %*%
            fxxinv
        scale <- mean(f)
        serr <- sqrt(diag(cov))
    }
    else if (se == "ker") {
        h <- bandwidth.rq(tau, n, hs = hs)
        if (tau + h > 1)
            stop("tau + h > 1:  error in summary.rq")
        if (tau - h < 0)
            stop("tau - h < 0:  error in summary.rq")
        uhat <- c(y - x %*% coef)
        h <- (qnorm(tau + h) - qnorm(tau - h)) * min(sqrt(var(uhat)),
            (quantile(uhat, 0.75) - quantile(uhat, 0.25))/1.34)
        f <- dnorm(uhat/h)/h
        fxxinv <- diag(p)
        fxxinv <- backsolve(qr(sqrt(f) * x)$qr[1:p, 1:p, drop = FALSE],
            fxxinv)
        fxxinv <- fxxinv %*% t(fxxinv)
        cov <- tau * (1 - tau) * fxxinv %*% crossprod(x) %*%
            fxxinv
        scale <- mean(f)
        serr <- sqrt(diag(cov))
    }
    else if (se == "boot") {
        B <- boot.rq(x, y, tau, ...)
        cov <- cov(B)
        serr <- sqrt(diag(cov))
    }


从你取的tau可以看出,你的tau+h很可能会超过1.所以解决方法有2种方法:
第一是换过一种方法产生标准误:
y=c(2167,2509,3530,4669,5868,6763,6820,2866,8248,8868,9336,10464,11040,12631,13773,14762,17255,19398,20992,23200)
income=c(2486,3009,4277,5868,7172,8159,8439,8773,10932,11718,12883,13250,14867,16683,18645,20668,23623,26675,28838,21754)
cpi=c(110.5,110.0,120.2,123.9,118.7,109.2,102.8,100.0,101.5,102.5,100.0,100.5,100.1,102.2,101.0,101.2,103.2,105.8,99.6,103.1)
library(quantreg)
fit1=rq(y~income+cpi,tau=c(0.2,0.5,0.8,0.99))
summary(fit1,se="boot")

第二种是改变tau值,我只尝试tau=0.5可以运行


y=c(2167,2509,3530,4669,5868,6763,6820,2866,8248,8868,9336,10464,11040,12631,13773,14762,17255,19398,20992,23200)
income=c(2486,3009,4277,5868,7172,8159,8439,8773,10932,11718,12883,13250,14867,16683,18645,20668,23623,26675,28838,21754)
cpi=c(110.5,110.0,120.2,123.9,118.7,109.2,102.8,100.0,101.5,102.5,100.0,100.5,100.1,102.2,101.0,101.2,103.2,105.8,99.6,103.1)
library(quantreg)
fit1=rq(y~income+cpi,tau=c(0.5))
summary(fit1,se="nid")




已有 3 人评分学术水平 热心指数 信用等级 收起 理由
耕耘使者 + 1 + 2 热心帮助其他会员
kk22boy + 2 + 2 + 2 观点有启发
epoh + 1 + 1 + 1 热心帮助其他会员

总评分: 学术水平 + 4  热心指数 + 5  信用等级 + 3   查看全部评分

一份耕耘,一份收获。

藤椅
ywh19860616 发表于 2012-5-5 12:24:59
tau=c(0.2,0.5,0.8,0.99)
n=20
h=bandwidth.rq(tau, n, hs=TRUE, alpha=0.05)

#h
[1] 0.21062684 0.35792541 0.21062684 0.02586753
所以产生了tau + h > 1 和tau-h<0错误。
一份耕耘,一份收获。

板凳
耕耘使者 发表于 2012-5-7 19:54:27
谢谢!

您需要登录后才可以回帖 登录 | 我要注册

本版微信群
加好友,备注cda
拉您进交流群
GMT+8, 2026-1-4 02:29