Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

分位点回帰

分位点τ\tauにおける条件付分位関数を

Qτ(yi∣Xi)=Fy−1(τ∣Xi)Q_\tau(y_i | X_i) = F_y^{-1}(\tau | X_i)
(1)

と表す。ここでFy−1(τ∣Xi)F_y^{-1}(\tau | X_i)はyyにおいてXiX_iに条件づけられたyiy_iの分布関数である（Fy−1(τ∣Xi)=inf⁡{y:Fy(y∣Xi)≥τ}F_y^{-1}(\tau | X_i) = \inf \{ y: F_y(y|X_i) \geq \tau \}）。

例えばτ=0.1\tau = 0.1のとき、Qτ(yi∣Xi)Q_\tau(y_i | X_i)はyiy_iの下位10分位である。

標準的な回帰モデルは二乗誤差(yi−m(Xi))2(y_i - m(X_i))^2の和や期待値を最小化するようにモデルm(Xi)m(X_i)を学習して条件付き期待値E(yi∣Xi)E(y_i|X_i)を予測する

\newcommand{\argmin} attempting to redefine \argmin; use \renewcommand

\newcommand{\argmin}{\mathop{\rm arg~min}\limits}
E(y_i|X_i) = \argmin_{m(X_i)}
E\big[ (y_i - m(X_i))^2 \big]

分位点回帰 （quantile regression）モデルはpinball lossρτ(yi−q(Xi))\rho_{\tau}(y_i - q(X_i))の和や期待値を最小化するようにモデルq(Xi)q(X_i)を学習させ、条件付き分位関数Qτ(yi∣Xi)=Fy−1(τ∣Xi)Q_{\tau}(y_i|X_i) = F^{-1}_y(\tau|X_i)を予測する

Qτ(yi∣Xi)=arg min⁡q(Xi)E[ρτ(yi−q(Xi))]Q_{\tau}(y_i|X_i) = \argmin_{q(X_i)} E\big[ \rho_{\tau}(y_i - q(X_i)) \big]
(3)

pinball lossは τ\tau-tiled absolute value function や 検定関数（check function）とも呼ばれる（グラフを描くとチェックマークに似てるため）

ρτ(x)=(τ−1(x≤0))x\rho_{\tau} (x) = \big(\tau - \mathbb{1}(x \leq 0) \big) x
(4)

あるいは

ρτ(x)={(τ−1)x if x≤0τx if x>0\rho_{\tau} (x) = \begin{cases} (\tau - 1) x & \text{ if } x \leq 0\\ \tau x & \text{ if } x > 0\\ \end{cases}
(5)

と書かれる

Source
import numpy as np
import matplotlib.pyplot as plt

def pinball_loss(x, tau):
    return (tau - 1 * (x <= 0)) * x

x = np.linspace(-1, 1, 100)
fig, axes = plt.subplots(figsize=[10, 2], ncols=3)
for i, tau in enumerate([0.1, 0.5, 0.9]):
    y = pinball_loss(x, tau=tau)
    axes[i].plot(x, y)
    if i == 0:
        axes[i].set(title=f"τ={tau}", xlabel=r"$x$", ylabel=r"$y = (\tau - 1(x <= 0)) x$")
    else:
        axes[i].set(title=f"τ={tau}", xlabel=r"$x$")
fig.show()
<Figure size 1000x200 with 3 Axes>

なお、pinball lossはτ=0.5\tau=0.5のとき

ρ0.5(x)={−0.5x if x≤00.5x if x>0=12∣x∣\begin{align} \rho_{0.5} (x) &= \begin{cases} -0.5 x & \text{ if } x \leq 0\\ 0.5 x & \text{ if } x > 0\\ \end{cases} \\ &= \frac{1}{2} |x| \end{align}
(10)

と、絶対誤差と比例する形になる。

絶対誤差の和を目的関数にとった線形モデルは統計学においてleast absolute deviations (LAD) と呼ばれ、その解は条件付き中央値になる

median(yi∣Xi)=Q0.5(yi∣Xi)=arg min⁡q(Xi)E[ρ0.5(yi−q(Xi))]\text{median}(y_i|X_i) = Q_{0.5}(y_i|X_i) = \argmin_{q(X_i)} E\big[ \rho_{0.5}(y_i - q(X_i)) \big]
(11)
Source
import numpy as np
import matplotlib.pyplot as plt

def pinball_loss(x, tau):
    return (tau - 1 * (x <= 0)) * x

x = np.linspace(-3, 3, 100)
fig, ax = plt.subplots()
ax.plot(x, pinball_loss(x, tau=0.5), label=r"$\rho_{0.5}(x)$")
ax.plot(x, abs(x), label="|x|")
ax.legend()
fig.show()
<Figure size 640x480 with 1 Axes>
絶対誤差の最適解

誤差関数ℓ(y,y^)\ell(y, \hat{y})を絶対誤差∣y−y^∣|y - \hat{y}|とする。予測損失（期待予測誤差）は

R(y^)=EY[ℓ(y,y^)]=∫−∞∞∣y−y^∣f(y)dy\begin{aligned} R(\hat{y}) &= \mathbb{E}_Y[\ell(y, \hat{y})]\\ &= \int_{-\infty}^{\infty} |y - \hat{y}| f(y) dy\\ \end{aligned}
(12)

絶対値の中身の符号で場合分けすると

R(y^)=∫−∞∞∣y−y^∣f(y)dy=∫y^∞(y−y^)f(y)dy+∫−∞y^(y^−y)f(y)dy(絶対値を符号ごとに場合分けした)=∫y^∞yf(y)dy−y^∫y^∞f(y)dy+y^∫−∞y^f(y)dy−∫−∞y^yf(y)dy(∵積分の線形性)\begin{aligned} R(\hat{y}) &= \int_{-\infty}^{\infty} |y - \hat{y}| f(y) dy\\ &= \int_{\hat{y}}^{\infty} (y - \hat{y}) f(y) dy + \int_{-\infty}^{\hat{y}} (\hat{y} - y) f(y) dy \quad (絶対値を符号ごとに場合分けした)\\ &= \int_{\hat{y}}^{\infty} y f(y) dy - \hat{y} \int_{\hat{y}}^{\infty} f(y) dy + \hat{y} \int_{-\infty}^{\hat{y}} f(y) dy - \int_{-\infty}^{\hat{y}} y f(y) dy \quad (\because 積分の線形性) \end{aligned}
(13)

予測損失を微分するとそれぞれの項は

ddy^∫y^∞yf(y)dy=−y^f(y^)ddy^(−y^⋅∫y^∞f(y)dy)=ddy^(−y^)∫y^∞f(y)dy+(−y^)ddy^(∫y^∞f(y)dy)(積の微分)=−1⋅∫y^∞f(y)dy−y^⋅−f(y^)=−∫y^∞f(y)dy+y^f(y^)ddy^(y^⋅∫−∞y^f(y)dy)=ddy^y^⋅∫−∞y^f(y)dy+y^⋅ddy^∫−∞y^f(y)dy(積の微分)=∫−∞y^f(y)dy+y^⋅f(y^)ddy^(−∫−∞y^yf(y)dy)=−y^f(y^)\begin{aligned} \frac{d}{d \hat{y}} \int_{\hat{y}}^{\infty} y f(y) dy &= - \hat{y} f(\hat{y}) \\ \frac{d}{d \hat{y}} \left( - \hat{y} \cdot \int_{\hat{y}}^{\infty} f(y) dy \right) &= \frac{d}{d \hat{y}} \left( - \hat{y} \right) \int_{\hat{y}}^{\infty} f(y) dy + (- \hat{y}) \frac{d}{d \hat{y}} \left( \int_{\hat{y}}^{\infty} f(y) dy \right) \quad (積の微分)\\ &= -1 \cdot \int_{\hat{y}}^{\infty} f(y) dy - \hat{y} \cdot - f(\hat{y})\\ &= -\int_{\hat{y}}^{\infty} f(y) dy + \hat{y} f(\hat{y}) \\ \frac{d}{d \hat{y}} \left( \hat{y} \cdot \int_{-\infty}^{\hat{y}} f(y) dy \right) &= \frac{d}{d \hat{y}} \hat{y} \cdot \int_{-\infty}^{\hat{y}} f(y) dy + \hat{y} \cdot \frac{d}{d \hat{y}} \int_{-\infty}^{\hat{y}} f(y) dy \quad (積の微分)\\ &= \int_{-\infty}^{\hat{y}} f(y) dy + \hat{y} \cdot f(\hat{y}) \\ \frac{d}{d \hat{y}} \left( - \int_{-\infty}^{\hat{y}} y f(y) dy \right) &= - \hat{y} f(\hat{y}) \end{aligned}
(14)

よって

dR(y^)dy^=−y^f(y^)−∫y^∞f(y)dy+y^f(y^)+∫−∞y^f(y)dy+y^f(y^)−y^f(y^)=−∫y^∞f(y)dy+∫−∞y^f(y)dy\begin{aligned} \frac{d R(\hat{y})}{d \hat{y}} &= - \hat{y} f(\hat{y}) -\int_{\hat{y}}^{\infty} f(y) dy + \hat{y} f(\hat{y}) + \int_{-\infty}^{\hat{y}} f(y) dy + \hat{y} f(\hat{y}) - \hat{y} f(\hat{y})\\ &= -\int_{\hat{y}}^{\infty} f(y) dy + \int_{-\infty}^{\hat{y}} f(y) dy \end{aligned}
(15)

となる。dR(y^)dy^=0\frac{d R(\hat{y})}{d \hat{y}}=0とおいて整理すれば

∫y^∞f(y)dy=∫−∞y^f(y)dy\int_{\hat{y}}^{\infty} f(y) dy = \int_{-\infty}^{\hat{y}} f(y) dy
(16)

となる点が予測損失を極小化することがわかる。これはy^\hat{y}が中央値となるときである。

f(y)f(y)は確率密度関数なので、∫−∞∞f(y)dy=1\int_{-\infty}^{\infty} f(y) dy=1になる。−∞-\inftyからy^\hat{y}への積分とy^\hat{y}から∞\inftyへの積分が等しくなるのはその半分、すなわち

∫−∞y^f(y)dy=12\int_{-\infty}^{\hat{y}} f(y) dy = \frac{1}{2}
(17)

である。なお、これは累積分布関数P(y^)\mathrm{P}(\hat{y})に等しい。よってy^\hat{y}は中央値である。

なお、中央値の定義には、以下の式を満たすmm

∫−∞m dF(x)≥12 and ∫m∞dF(x)≥12\int_{-\infty}^m \mathrm{~d} F(x) \geq \frac{1}{2} \text { and } \int_m^{\infty} \mathrm{d} F(x) \geq \frac{1}{2}
(18)

というものもある。

参考：定積分の微分

定積分の定義

∫axf(t)dt=F(x)−F(a)\int_a^x f(t) d t=F(x)-F(a)
(19)

より、この導関数は

F′(x)−0=f(x)F'(x)-0 = f(x)
(20)
Pinball lossの最適解

Pinball Lossを少し表現を変えて

ρτ(x)={τ∣x∣ if x>0(1−τ)∣x∣ if x≤0\rho_{\tau} (x) = \begin{cases} \tau |x| & \text{ if } x > 0\\ (1 - \tau) |x| & \text{ if } x \leq 0\\ \end{cases}
(21)

と表すと、さきほどの絶対誤差の場合分けした項にτ,(1−τ)\tau, (1 - \tau)を掛けた形になる

R(y^)=τ∫y^∞(y−y^)f(y)dy+(1−τ)∫−∞y^(y^−y)f(y)dy\begin{aligned} R(\hat{y}) &= \tau \int_{\hat{y}}^{\infty} (y - \hat{y}) f(y) dy + (1 - \tau) \int_{-\infty}^{\hat{y}} (\hat{y} - y) f(y) dy \\ \end{aligned}
(22)

絶対誤差の場合と同様に導関数は

dR(y^)dy^=−τ∫y^∞f(y)dy+(1−τ)∫−∞y^f(y)dy\begin{aligned} \frac{d R(\hat{y})}{d \hat{y}} &= - \tau \int_{\hat{y}}^{\infty} f(y) dy + (1 - \tau)\int_{-\infty}^{\hat{y}} f(y) dy \end{aligned}
(23)

となる。ここで累積分布関数F(y^)F(\hat{y})と補累積分布関数をそれぞれ

F(y^):=∫−∞y^f(y)dy1−F(y^):=∫y^∞f(y)dyF(\hat{y}) := \int_{-\infty}^{\hat{y}} f(y) dy\\ 1 - F(\hat{y}) := \int_{\hat{y}}^{\infty} f(y) dy\\
(24)

とおき、導関数を0とおく。

dR(y^)dy^=−τ(1−F(y^))+(1−τ)F(y^)=0\begin{aligned} \frac{d R(\hat{y})}{d \hat{y}} &= - \tau (1 - F(\hat{y})) + (1 - \tau) F(\hat{y}) = 0 \end{aligned}
(25)

これを整理すると

τ(1−F(y^))=(1−τ)F(y^)τ−τF(y^))=(1−τ)F(y^)τ=τF(y^))+(1−τ)F(y^)τ=(τ+(1−τ))F(y^)τ=F(y^)\begin{aligned} \tau (1 - F(\hat{y})) &= (1 - \tau) F(\hat{y})\\ \tau - \tau F(\hat{y})) &= (1 - \tau) F(\hat{y})\\ \tau &= \tau F(\hat{y})) + (1 - \tau) F(\hat{y})\\ \tau &= (\tau + (1 - \tau)) F(\hat{y})\\ \tau &= F(\hat{y})\\ \end{aligned}
(26)

となり、累積分布F(y^)F(\hat{y})で表される確率（分布のうちYYが−∞-\infty から y^\hat{y}までの面積）がτ\tauになることがわかる。 逆関数をとると

y^=F−1(τ)\hat{y} = F^{-1}(\tau)\\
(27)

となる。累積分布関数の逆関数は分位点なので、y^\hat{y}は確率τ\tauに対応する分位点である。

モデルの評価

D2 pinball score

D2D^2はR2R^2の一般化

D2(y,y^)=1−dev(y,y^)dev(y,ynull)D^2(y, \hat{y}) = 1 - \frac{\text{dev}(y, \hat{y})}{\text{dev}(y, y_{\text{null}})}
(28)

ここでynully_{\text{null}}は切片のみのモデルの最適解（例：二乗誤差ならyyの平均値、絶対誤差ならyyの中央値、pinball lossならyyの指定されたquantile）

このD2D^2に

dev(y,y^)=pinball(y,y^)\text{dev}(y, \hat{y}) = \text{pinball}(y, \hat{y})
(29)

を代入したものがD2D^2 pinball score

interval score

分位点回帰モデルの実践

statsmodelsでは quantreg() で実行できる

Quantile regression - statsmodels 0.15.0 (+213)

Source
import numpy as np
import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf
import matplotlib.pyplot as plt

data = sm.datasets.engel.load_pandas().data

fig, ax = plt.subplots()
ax.scatter(data["income"], data["foodexp"])
ax.set(xlabel="income", ylabel="foodexp", title="Quantile Linear Regression")

x = np.linspace(data["income"].min(), data["income"].max(), 10)
model = smf.quantreg("foodexp ~ income", data)
for q in [0.1, 0.5, 0.9]:
    res = model.fit(q=q)
    y_hat = res.predict(pd.DataFrame({"income": x}))
    ax.plot(x, y_hat, label=fr"$\tau = {q}$")
ax.legend()
fig.show()
/tmp/ipykernel_1497/2518785238.py:2: DeprecationWarning: 
Pyarrow will become a required dependency of pandas in the next major release of pandas (pandas 3.0),
(to allow more performant data types, such as the Arrow string type, and better interoperability with other libraries)
but was not found to be installed on your system.
If this would cause problems for you,
please provide us feedback at https://github.com/pandas-dev/pandas/issues/54466
        
  import pandas as pd
<Figure size 640x480 with 1 Axes>
統計モデリング
ガウス過程回帰
統計モデリング
生存分析