將分位數回歸公式化為線性規劃問題?
如何將分位數回歸表述為線性規劃問題?在查看中位數分位數問題時,我知道它是
minimize n∑i=1|β0+Xiβ1−Yi| transforms into minimize n∑i=1ei s.t. ei≥β0+Xiβ1−Yi ei≥−(β0+Xiβ1−Yi)
但是如何轉換其他分位數的最小化?
您使用分位數回歸估計器
ˆβ(τ):=argminθ∈RKN∑i=1ρτ(yi−x⊤iθ).
在哪裡 τ∈(0,1) 是根據需要估計的分位數和函數選擇的常數 ρτ(.) 定義為
ρτ(r)=r(τ−I(r<0)).
為了看清目的 ρτ(.) 首先考慮將殘差作為參數,當它們被定義為 ϵi=yi−x⊤iθ . 因此,最小化問題中的總和可以重寫為
N∑i=1ρτ(ϵi)=N∑i=1τ|ϵi|I[ϵi≥0]+(1−τ)|ϵi|I[ϵi<0]
使得與觀察相關的正殘差 yi 高於建議的分位數回歸線 x⊤iθ 被賦予權重 τ 而與觀察相關的負殘差 yi 低於建議的分位數回歸線 x⊤iθ 加權 (1−τ) .
直觀地說:
和 τ=0.5 正負殘差以相同的權重被“懲罰”,並且相同數量的觀察值在最優的“線”上方和下方,因此該線 x⊤iˆβ 是中位數回歸“線”。
什麼時候 τ=0.9 每個正殘差的權重是負殘差的 9 倍 1−τ=0.1 因此對於“線”上方的每個觀察都是最佳的 x⊤iˆβ 大約 9 將放置在該線下方。因此,“線”代表 0.9 分位數。(有關這方面的確切陳述,請參閱 Koenker (2005)“分位數回歸”中的 THM.2.2 和推論 2.1)
這些圖中說明了這兩種情況。左面板 τ=0.5 和右面板 τ=0.9 .
線性程序主要使用標準形式進行分析和求解
(1) minz c⊤z subject to Az=b,z≥0
要達到標準形式的線性程序,第一個問題是在這樣的程序中 (1) 執行最小化的所有變量 z 應該是積極的。為了實現這一點,使用鬆弛變量將殘差分解為正負部分:
ϵi=ui−vi
積極的部分 ui=max(0,ϵi)=|ϵi|I[ϵi≥0] 和 vi=max(0,−ϵi)=|ϵi|I[ϵi<0] 是負面的部分。由檢查函數分配權重的殘差總和被視為
N∑i=1ρτ(ϵi)=N∑i=1τui+(1−τ)vi=τ1⊤Nu+(1−τ)1⊤Nv,
在哪裡 u=(u1,…,uN)⊤ 和 v=(v1,…,vN)⊤ 和 1N 是向量 N×1 所有坐標等於 1 .
殘差必須滿足 N 約束
yi−x⊤iθ=ϵi=ui−vi
這導致公式為線性程序
minθ∈RK,u∈RN+,v∈RN+τ1⊤Nu+(1−τ)1⊤Nv|yi=xiθ+ui−vi,i=1,…,N,
如 Koenker (2005)“分位數回歸”第 10 頁方程 (1.20) 所述。
然而值得注意的是 θ∈R 仍然不限於標準形式(1)上的線性程序中要求的為正。因此再次使用分解成正負部分
θ=θ+−θ−
又在哪裡 θ+=max(0,θ) 是積極的部分並且 θ−=max(0,−θ) 是負面的部分。這 N 約束可以寫成
y:=[y1 ⋮ yN]=[x⊤1 ⋮ x⊤N](θ+−θ−)+INu−INv,
在哪裡 IN=diag1N .
接下來定義 b:=y 和設計矩陣 X 將自變量上的數據存儲為
X:=[x⊤1 ⋮ x⊤N]
重寫約束:
b=X(θ+−θ−)+INu−INv=[X,−X,IN,−IN][θ+ θ− u v]
定義 (N×2K+2N) 矩陣
A:=[X,−X,IN,−IN]
並介紹 θ+ 和 θ− 作為要最小化的變量,因此它們是 z 要得到b=A[θ+ θ− u v]=Az
因為 θ+ 和 θ− 僅通過約束 a 影響最小化問題 0 維度的 2K×1 必須作為係數向量的一部分引入 c 可以適當地定義為
c=[0 τ1N (1−τ)1N],
從而確保 $ c^\top z = \underbrace{\mathbf 0^\top(\theta^+ - \theta^-)}{=0}+\tau \mathbf 1_N^\top u + (1-\tau)\mathbf 1_N^\top v = \sum{i=1}^N \rho_\tau(\epsilon_i). $
因此 c,A 和 b 然後定義和程序如 (1) 完全指定。
這可能最好用一個例子來消化。要在 R 中解決這個問題,請使用 Roger Koenker 的 quantreg 包。這裡還說明瞭如何設置線性程序並使用線性程序求解器求解:
base=read.table("http://freakonometrics.free.fr/rent98_00.txt",header=TRUE) attach(base) library(quantreg) library(lpSolve) tau <- 0.3 # Problem (1) only one covariate X <- cbind(1,base$area) K <- ncol(X) N <- nrow(X) A <- cbind(X,-X,diag(N),-diag(N)) c <- c(rep(0,2*ncol(X)),tau*rep(1,N),(1-tau)*rep(1,N)) b <- base$rent_euro const_type <- rep("=",N) linprog <- lp("min",c,A,const_type,b) beta <- linprog$sol[1:K] - linprog$sol[(1:K+K)] beta rq(rent_euro~area, tau=tau, data=base) # Problem (2) with 2 covariates X <- cbind(1,base$area,base$yearc) K <- ncol(X) N <- nrow(X) A <- cbind(X,-X,diag(N),-diag(N)) c <- c(rep(0,2*ncol(X)),tau*rep(1,N),(1-tau)*rep(1,N)) b <- base$rent_euro const_type <- rep("=",N) linprog <- lp("min",c,A,const_type,b) beta <- linprog$sol[1:K] - linprog$sol[(1:K+K)] beta rq(rent_euro~ area + yearc, tau=tau, data=base)