將分位數回歸公式化為線性規劃問題？
如何將分位數回歸表述為線性規劃問題？在查看中位數分位數問題時，我知道它是
$$ \begin{align} \text{minimize } & \sum_{i=1}^n |\beta_0 + X_i \beta_1-Y_i|\ \text{transforms into } & \ \text{minimize } & \sum_{i=1}^n e_i\ \text{s.t.} & \ & e_i\geq \beta_0 + X_i\beta_{1}-Y_i\ & e_i\geq -(\beta_0 + X_i\beta_{1}-Y_i) \end{align} $$ 但是如何轉換其他分位數的最小化？
您使用分位數回歸估計器
$$ \hat \beta(\tau) := \arg \min_{\theta \in \mathbb R^K} \sum_{i=1}^N \rho_\tau(y_i - \mathbf x_i^\top \theta). $$
在哪裡 $ \tau \in (0,1) $ 是根據需要估計的分位數和函數選擇的常數 $ \rho_\tau(.) $ 定義為
$$ \rho_\tau(r) = r(\tau - I(r<0)). $$
為了看清目的 $ \rho_\tau(.) $ 首先考慮將殘差作為參數，當它們被定義為 $ \epsilon_i =y_i - \mathbf x_i^\top \theta $ . 因此，最小化問題中的總和可以重寫為
$$ \sum_{i=1}^N \rho_\tau(\epsilon_i) =\sum_{i=1}^N \tau \lvert \epsilon_i \lvert I[\epsilon_i \geq 0] + (1-\tau) \lvert \epsilon_i \lvert I[\epsilon_i < 0] $$
使得與觀察相關的正殘差 $ y_i $ 高於建議的分位數回歸線 $ \mathbf x_i^\top \theta $ 被賦予權重 $ \tau $ 而與觀察相關的負殘差 $ y_i $ 低於建議的分位數回歸線 $ \mathbf x_i^\top \theta $ 加權 $ (1-\tau) $ .
直觀地說：
和 $ \tau=0.5 $ 正負殘差以相同的權重被“懲罰”，並且相同數量的觀察值在最優的“線”上方和下方，因此該線 $ \mathbf x_i^\top \hat \beta $ 是中位數回歸“線”。
什麼時候 $ \tau=0.9 $ 每個正殘差的權重是負殘差的 9 倍 $ 1-\tau= 0.1 $ 因此對於“線”上方的每個觀察都是最佳的 $ \mathbf x_i^\top \hat \beta $ 大約 9 將放置在該線下方。因此，“線”代表 0.9 分位數。（有關這方面的確切陳述，請參閱 Koenker (2005)“分位數回歸”中的 THM.2.2 和推論 2.1）
這些圖中說明了這兩種情況。左面板 $ \tau=0.5 $ 和右面板 $ \tau=0.9 $ .
線性程序主要使用標準形式進行分析和求解
$$ (1) \ \ \min_z \ \ c^\top z \ \ \mbox{subject to } A z = b , z \geq 0 $$
要達到標準形式的線性程序，第一個問題是在這樣的程序中 (1) 執行最小化的所有變量 $ z $ 應該是積極的。為了實現這一點，使用鬆弛變量將殘差分解為正負部分：
$$ \epsilon_i = u_i - v_i $$
積極的部分 $ u_i = \max(0,\epsilon_i) = \lvert \epsilon_i \lvert I[\epsilon_i \geq 0] $ 和 $ v_i = \max(0,-\epsilon_i) =\lvert \epsilon_i \lvert I[\epsilon_i < 0] $ 是負面的部分。由檢查函數分配權重的殘差總和被視為
$$ \sum_{i=1}^N \rho_\tau(\epsilon_i) = \sum_{i=1}^N \tau u_i + (1-\tau) v_i = \tau \mathbf 1_N^\top u + (1-\tau)\mathbf 1_N^\top v, $$
在哪裡 $ u = (u_1,…,u_N)^\top $ 和 $ v=(v_1,…,v_N)^\top $ 和 $ \mathbf 1_N $ 是向量 $ N \times 1 $ 所有坐標等於 $ 1 $ .
殘差必須滿足 $ N $ 約束
$$ y_i - \mathbf x_i^\top\theta = \epsilon_i = u_i - v_i $$
這導致公式為線性程序
$$ \min_{\theta \in \mathbb R^K,u\in \mathbb R_+^N,v\in \mathbb R_+^N}{ \tau \mathbf 1_N^\top u + (1-\tau)\mathbf 1_N^\top v\lvert y_i= \mathbf x_i\theta + u_i - v_i, i=1,…,N}, $$
如 Koenker (2005)“分位數回歸”第 10 頁方程 (1.20) 所述。
然而值得注意的是 $ \theta\in \mathbb R $ 仍然不限於標準形式（1）上的線性程序中要求的為正。因此再次使用分解成正負部分
$$ \theta = \theta^+ - \theta^- $$
又在哪裡 $ \theta^+=max(0,\theta) $ 是積極的部分並且 $ \theta^- = \max(0,-\theta) $ 是負面的部分。這 $ N $ 約束可以寫成
$$ \mathbf y:= \begin{bmatrix} y_1 \ \vdots \ y_N\end{bmatrix} = \begin{bmatrix} \mathbf x_1^\top \ \vdots \ \mathbf x_N^\top \end{bmatrix}(\theta^+ - \theta^-) + \mathbf I_Nu - \mathbf I_Nv , $$
在哪裡 $ \mathbf I_N = diag{\mathbf 1_N} $ .
接下來定義 $ b:=\mathbf y $ 和設計矩陣 $ \mathbf X $ 將自變量上的數據存儲為
$$ \mathbf X := \begin{bmatrix} \mathbf x_1^\top \ \vdots \ \mathbf x_N^\top \end{bmatrix} $$
重寫約束：
$$ b= \mathbf X(\theta^+ - \theta^-) + \mathbf I_N u- \mathbf I_N v= [\mathbf X , -\mathbf X , \mathbf I_N ,\mathbf - \mathbf I_N] \begin{bmatrix} \theta^+ \ \theta^- \ u \ v\end{bmatrix} $$
定義 $ (N \times 2K + 2N ) $ 矩陣
$$ A := [\mathbf X , -\mathbf X , \mathbf I_N ,\mathbf - \mathbf I_N] $$ 並介紹 $ \theta^+ $ 和 $ \theta^- $ 作為要最小化的變量，因此它們是 $ z $ 要得到
$$ b = A \begin{bmatrix} \theta^+ \ \theta^- \ u \ v\end{bmatrix} = Az $$
因為 $ \theta^+ $ 和 $ \theta^- $ 僅通過約束 a 影響最小化問題 $ \mathbf 0 $ 維度的 $ 2K\times 1 $ 必須作為係數向量的一部分引入 $ c $ 可以適當地定義為
$$ c = \begin{bmatrix}\mathbf 0 \ \tau \mathbf 1_N \ (1-\tau) \mathbf 1_N \end{bmatrix}, $$
從而確保 $ 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)
