分位数回归

分位数回归（quantile regression）描述响应变量在给定协变量条件下的某个分位数。例如，同样年龄人群的中位数身高与较高身高分位数，可能随年龄以不同速度变化；只拟合条件均值，不能给出这两条曲线。分位数回归直接优化与目标分位数相应的不对称绝对损失。

分位数与检验损失

给定 $0<\tau <1$，若实数 $q$ 满足 $$F(q-)\le \tau \le F(q),$$ 则称它为一个 $\tau$ 分位数，其中 $F$ 是分布函数。这一定义允许离散分布和不唯一的分位数；常用的广义逆约定选取 $inf\{y:F(y)\ge \tau \}$。如果在分位数附近分布连续且严格递增，分位数便唯一且满足 $F(q)=\tau$。

定义检验损失 $${\rho }_{\tau }(u)=u(\tau -{\mathrm{𝟏}}_{\{u<0\}})=\{\begin{array}{cc}\tau u,& u\ge 0,\\ (1-\tau )(-u),& u<0.\end{array}$$ 残差取 $u=y-\hat{y}$：正残差表示预测偏低，负残差表示预测偏高。取 $\tau =0.75$ 时，低估一个单位的损失为零点七五，高估一个单位的损失为零点二五。符号和权重若反过来，求得的是另一端的分位数。

残差为真实值减预测值。分位数水平越高，预测偏低的代价越大；在零点五时两侧对称。

损失最小值与分位数的对应证明

设 $𝔼|Y|<\mathrm{\infty }$，考察 $L(q)=𝔼[{\rho }_{\tau }(Y-q)]$。每个样本的损失对 $q$ 都是凸函数：当 $q<Y$，斜率为 $-\tau$；当 $q>Y$，斜率为 $1-\tau$；在相等处有一个斜率区间。取期望得到左右导数 $$L{\text{'}}_{-}(q)=F(q-)-\tau ,L{\text{'}}_{+}(q)=F(q)-\tau .$$ 凸函数在 $q$ 处达到最小值，当且仅当左导数不大于零、右导数不小于零。因此最小值条件正是 $F(q-)\le \tau \le F(q)$。这解释了检验损失的选择，而不只是在最小二乘法中换一个符号。

条件分位数的线性模型

若 $x_i\in {ℝ}^p$ 已含截距列，线性分位数回归估计 $${\hat{\beta }}_{\tau }\in \arg \underset{\beta }{min}\sum _{i=1}^n{\rho }_{\tau }(y_i-x_i^{𝖳}\beta ).$$ 当真实条件分位数确为 $x^{𝖳}{\beta }_{\tau }$ 时，逐个条件分布的最优性可解释这个目标；若关系并不线性，所得只是按检验损失定义的线性近似。不同 $\tau$ 对应不同系数，不能把某个系数解读成所有个体共同的响应变化，更不能仅凭观测回归推出因果效应。

两组观测的完整计算

设 $x=0$ 组的响应为 $1,2,9$，$x=1$ 组为 $3,5,11$。模型 ${\beta }_0+{\beta }_{1}x$ 在两组分别给出 ${\beta }_0$ 和 ${\beta }_0+{\beta }_1$；这两个组位置可以独立变化，所以分别求组内分位数即可。

取 $\tau =1/2$，两组中位数分别为二、五，故 $\hat{\beta }=(2,3)$。残差为 $(-1,0,7,-2,0,6)$，总检验损失为 $$\frac{1}{2}(1+0+7+2+0+6)=8.$$ 取 $\tau =3/4$，每组三个等权观测的最高值是唯一的这一水平分位数，故 $\hat{\beta }=(9,2)$。残差为 $(-8,-7,0,-8,-6,0)$，总损失为 $(8+7+8+6)/4=7.25$。普通最小二乘给出的组均值却是四和 $19/3$，相应系数 $(4,7/3)$。三个目标估计的不是同一个位置特征。

同一组的数据，其均值、中位数和较高分位数都可以不同。连接两组的线来自指定的线性模型，不是在组间补造观测。

变式。 把第一组的九改成九十。中位数仍为二，第二组仍为五，因此中位数回归系数不变；第一组均值升到三十一，均值回归明显变化。但高分位数回归也会受这个尾部观测影响，不能笼统说分位数回归“完全不怕异常值”。

线性规划的可运行实现

引入非负向量 $u,v$ 使 $y-X\beta =u-v$，得到等价问题 $$\underset{\beta ,u,v}{min}\tau {\mathrm{𝟏}}^{𝖳}u+(1-\tau ){\mathrm{𝟏}}^{𝖳}v,X\beta +u-v=y,u,v\ge 0.$$ 系数 $\beta$ 应允许任意实数；把求解器默认的非负界限直接套给系数会改变问题。以下代码依赖NumPy、SciPy，显式检查求解成功并核对约束残差。

Python 3

import math
import numpy as np
from scipy.optimize import linprog

def quantile_fit(design, response, tau):
    x, y = np.asarray(design, dtype=float), np.asarray(response, dtype=float)
    if not (x.ndim == 2 and y.ndim == 1 and x.shape[0] == len(y)
            and min(x.shape) > 0 and np.all(np.isfinite(x)) and np.all(np.isfinite(y))):
        raise ValueError("finite nonempty paired arrays required")
    if not (math.isfinite(tau) and 0 < tau < 1):
        raise ValueError("tau must lie strictly between zero and one")
    n, p = x.shape
    objective = np.r_[np.zeros(p), np.full(n, tau), np.full(n, 1-tau)]
    equality = np.c_[x, np.eye(n), -np.eye(n)]
    result = linprog(objective, A_eq=equality, b_eq=y,
                     bounds=[(None, None)]*p + [(0, None)]*(2*n), method="highs")
    if not result.success:
        raise RuntimeError(result.message)
    if np.max(abs(equality @ result.x - y)) > 1e-7 * (1 + np.max(abs(y))):
        raise RuntimeError("equality residual too large")
    beta = result.x[:p].copy()
    residual = y - x @ beta
    loss = np.sum(np.where(residual >= 0, tau*residual, (tau-1)*residual))
    return beta, float(loss)

x = np.c_[np.ones(6), [0, 0, 0, 1, 1, 1]]
print(quantile_fit(x, [1, 2, 9, 3, 5, 11], 0.5))
print(quantile_fit(x, [1, 2, 9, 3, 5, 11], 0.75))

输出分别为系数 [2. 3.]、损失八，以及系数 [9. 2.]、损失七点二五。该教学实现构造稠密约束矩阵，存储为 $O(n(p+n))$；大型问题应使用稀疏表示。求解时间取决于线性规划规模、数据及求解器，不能用构造矩阵的时间冒充总求解复杂度。输入数组不被修改；共线设计或离散数据可能使最优系数不唯一。

单独拟合不同分位数的直线还可能相交。主例中 $2+3x$ 与 $9+2x$ 在 $x=7$ 相交，之后中位数线反而更高；这发生在观测范围之外，说明不能把局部拟合无限外推。预测区间、参数置信区间和条件分位数曲线也各有含义，不能互相替换。