分位点回帰

分位点回帰#

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

Qτ(yi|Xi)=Fy−1(τ|Xi)

と表す。ここでFy−1(τ|Xi)はyにおいてXiに条件づけられたyiの分布関数である（Fy−1(τ|Xi)=inf{y:Fy(y|Xi)≥τ}）。

例えばτ=0.1のとき、Qτ(yi|Xi)はyiの下位10分位である。

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

E(yi|Xi)=arg minm(Xi)⁡E[(yi−m(Xi))2]

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

Qτ(yi|Xi)=arg minq(Xi)⁡E[ρτ(yi−q(Xi))]

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

ρτ(x)=(τ−1(x≤0))x

あるいは

ρτ(x)={(τ−1)x if x≤0τx if x>0

と書かれる

Hide code cell 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()
../_images/54794bf6727a7695472b5e0372e6c7d57d1fbc5662994c2bf4d08f44988701fd.png

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

ρ0.5(x)={−0.5x if x≤00.5x if x>0=12|x|

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

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

median(yi|Xi)=Q0.5(yi|Xi)=arg minq(Xi)⁡E[ρ0.5(yi−q(Xi))]
Hide code cell 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()
../_images/0963fe9b4ff046432ea5c0e5617a427b32cd329c48c64c4b718f46777d05f849.png
絶対誤差の最適解

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

R(y^)=EY[ℓ(y,y^)]=∫−∞∞|y−y^|f(y)dy

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

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(∵積分の線形性)

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

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^)

よって

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

となる。dR(y^)dy^=0とおいて整理すれば

∫y^∞f(y)dy=∫−∞y^f(y)dy

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

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

∫−∞y^f(y)dy=12

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

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

∫−∞m dF(x)≥12 and ∫m∞dF(x)≥12

というものもある。

参考：定積分の微分

定積分の定義

∫axf(t)dt=F(x)−F(a)

より、この導関数は

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

Pinball Lossを少し表現を変えて

ρτ(x)={τ|x| if x>0(1−τ)|x| if x≤0

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

R(y^)=τ∫y^∞(y−y^)f(y)dy+(1−τ)∫−∞y^(y^−y)f(y)dy

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

dR(y^)dy^=−τ∫y^∞f(y)dy+(1−τ)∫−∞y^f(y)dy

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

F(y^):=∫−∞y^f(y)dy1−F(y^):=∫y^∞f(y)dy

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

dR(y^)dy^=−τ(1−F(y^))+(1−τ)F(y^)=0

これを整理すると

τ(1−F(y^))=(1−τ)F(y^)τ−τF(y^))=(1−τ)F(y^)τ=τF(y^))+(1−τ)F(y^)τ=(τ+(1−τ))F(y^)τ=F(y^)

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

y^=F−1(τ)

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

モデルの評価#

D2 pinball score#

D2はR2の一般化

D2(y,y^)=1−dev(y,y^)dev(y,ynull)

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

このD2に

dev(y,y^)=pinball(y,y^)

を代入したものがD2 pinball score

interval score#

分位点回帰モデルの実践#

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

Quantile regression - statsmodels 0.15.0 (+213)

Hide code cell 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()
../_images/5010b23877a3ca8ed8efeb3e6d37afdf6a18303bca21a065b2acacc2107a04cd.png
Contents