分位点回帰と通常の線形回帰の違いは？

通常の線形回帰は条件付き平均を推定するのに対し、分位点回帰は中央値や任意のパーセンタイル（10%、90% など）を推定します。予測区間や非対称な誤差を扱いたいときに有効です。

どの分位点を選ぶべき？

中央値（0.5）はロバストな代替指標、10%・90% などの両端は予測区間の構築に使われます。需要予測や在庫最適化では、欠品コストと過剰在庫コストの比から最適な分位点が決まります。

勾配ブースティングと組み合わせられる？

LightGBM や XGBoost には quantile 損失があり、非線形な分位点回帰が可能です。複数分位点を学習して予測区間を作るのが実務でよく使われる手法です。

分位点回帰 (Quantile Regression)

中級

02.1.7

分位点回帰 (Quantile Regression)

最終更新 2020-04-22 読了時間 4 分

まとめ

分位点回帰は、平均ではなく「観測値の 10% がこの線より下に来る」のような分位点 $\tau$ の直線を推定する回帰です。損失を二乗誤差からピンボール損失に替えるだけで、上にはずれた残差と下にはずれた残差に $\tau$ と $1-\tau$ の異なる重みをつけます。$\tau = 0.5$ なら中央値回帰になり、外れ値に強い直線が得られます。このページでは右に長いガンマ分布ノイズを足した 400 点に QuantileRegressor で $\tau = 0.1, 0.5, 0.9$ の 3 本を引き、最小二乗法の平均直線と重ねて、点の帯の下端・中央・上端をそれぞれの直線が通る様子を確かめます。

先に線形回帰で最小二乗法が「残差の 2 乗」を最小にしていることを押さえておくと、このページで変わるのが損失関数の形だけだと分かります。

直感 #

配送時間を予測するとき、平均 30 分という 1 本の直線だけでは「90% の荷物は何分以内に届くか」に答えられません。分位点回帰は $\tau = 0.9$ を指定すると、観測点の 90% がその下に来る直線を直接求めます。仕組みは残差の数え方にあり、直線より上の点には $\tau$、下の点には $1 - \tau$ の重みをかけて絶対値を足すので、$\tau = 0.9$ では上にはみ出す点 1 個が下の点 9 個分の罰になります。その結果、$\tau = 0.9$ の直線は点の帯の上端近くまで押し上げられます。$\tau$ を 0.1、0.5、0.9 と 3 本引けば、帯の下端・中央・上端が 3 本の直線として見えます。

何が起きているか（ステップで見る） #

コードを読む前に、この手法が図の上で何をしているかを追ってみましょう。「再生」を押すか、スライダーを動かすと 1 手順ずつ進みます。

見るポイント：誤差の大きい点の重みを下げるだけで、直線は多数派に戻る。分位点回帰の $\tau = 0.5$（中央値回帰）も残差を 2 乗せず絶対値で数えるので、同じように外れ値 1 点に引っぱられません。

式で書くと何を最小にしているか #

二乗誤差を、傾きが $\tau$ と $\tau - 1$ の 2 本の直線でできた V 字の損失に置き換えます。残差を $r = y - \hat{y}$、分位点を $\tau \in (0,1)$ とすると、ピンボール損失は

$$ L_\tau(r) = \begin{cases} \tau \, r & (r \ge 0) \\ (\tau - 1) r & (r < 0) \end{cases} $$

で定義されます。$r \ge 0$（点が直線の上）には傾き $\tau$、$r < 0$（点が直線の下）には傾き $1 - \tau$ の罰がかかり、この非対称さが直線を上下にずらします。この損失の和を最小化すると $\tau$ 分位点に対応する線形予測子が得られ、$\tau = 0.5$ なら両側の傾きが 0.5 で等しい絶対値損失になるので中央値回帰です。

Pythonを用いた実験や説明 #

次のコードは、真の関係 $y = 1.5x + 5$ に右へ長い裾を持つガンマ分布のノイズを足した 400 点を作り、QuantileRegressor で $\tau = 0.1, 0.5, 0.9$ の 3 本の直線を推定します。LinearRegression の平均直線を破線で重ね、ノイズが非対称なときに平均と中央値がどれだけ離れるかを図で確認します。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
from __future__ import annotations

import matplotlib.pyplot as plt
import numpy as np
from sklearn.linear_model import LinearRegression, QuantileRegressor
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

def run_quantile_regression_demo(
    taus: tuple[float, ...] = (0.1, 0.5, 0.9),
    n_samples: int = 400,
    xlabel: str = "input x",
    ylabel: str = "output y",
    label_observations: str = "observations",
    label_mean: str = "mean (OLS)",
    label_template: str = "quantile τ={tau}",
    title: str | None = None,
) -> dict[float, tuple[float, float]]:
    """Fit quantile regressors alongside OLS and plot the conditional bands.

    Args:
        taus: Quantile levels to fit (each in (0, 1)).
        n_samples: Number of synthetic observations to generate.
        xlabel: Label for the x-axis.
        ylabel: Label for the y-axis.
        label_observations: Legend label for the scatter plot.
        label_mean: Legend label for the OLS line.
        label_template: Format string for quantile labels.
        title: Optional title for the plot.

    Returns:
        Mapping of quantile level to (min prediction, max prediction).
    """
    rng = np.random.default_rng(123)

    x_values: np.ndarray = np.linspace(0.0, 10.0, n_samples, dtype=float)
    noise: np.ndarray = rng.gamma(shape=2.0, scale=1.0, size=n_samples) - 2.0
    y_values: np.ndarray = 1.5 * x_values + 5.0 + noise
    X: np.ndarray = x_values[:, np.newaxis]

    quantile_models: dict[float, make_pipeline] = {}
    for tau in taus:
        model = make_pipeline(
            StandardScaler(with_mean=True),
            QuantileRegressor(alpha=0.001, quantile=float(tau), solver="highs"),
        )
        model.fit(X, y_values)
        quantile_models[tau] = model

    ols = LinearRegression()
    ols.fit(X, y_values)

    grid: np.ndarray = np.linspace(0.0, 10.0, 200, dtype=float)[:, np.newaxis]
    preds = {tau: model.predict(grid) for tau, model in quantile_models.items()}
    ols_pred: np.ndarray = ols.predict(grid)

    fig, ax = plt.subplots(figsize=(10, 5))
    ax.scatter(X, y_values, s=15, alpha=0.4, label=label_observations)

    color_cycle = plt.rcParams["axes.prop_cycle"].by_key().get("color", ["#1f77b4", "#ff7f0e", "#2ca02c"])
    for idx, tau in enumerate(taus):
        color = color_cycle[idx % len(color_cycle)]
        ax.plot(
            grid,
            preds[tau],
            color=color,
            linewidth=2,
            label=label_template.format(tau=tau),
        )

    ax.plot(grid, ols_pred, color="#9467bd", linestyle="--", label=label_mean)
    ax.set_xlabel(xlabel)
    ax.set_ylabel(ylabel)
    if title:
        ax.set_title(title)
    ax.legend()
    fig.tight_layout()
    plt.show()

    summary: dict[float, tuple[float, float]] = {
        tau: (float(pred.min()), float(pred.max())) for tau, pred in preds.items()
    }
    return summary

summary = run_quantile_regression_demo(
    xlabel="入力 x",
    ylabel="出力 y",
    label_observations="観測データ",
    label_mean="平均 (OLS)",
    label_template="分位点 τ={tau}",
    title="分位点回帰による条件分布",
)
for tau, (ymin, ymax) in summary.items():
    print(f"τ={tau:.1f}: 予測最小値 {ymin:.2f}, 予測最大値 {ymax:.2f}")

コードの解説番号をクリックすると、コード内の該当行がハイライトされます。

36-39 行目右に長い裾のノイズをわざと使う
rng.gamma(shape=2.0, scale=1.0) - 2.0 は平均 0 ですが、左は −2 で止まり右には長く伸びる非対称なノイズです。正規ノイズだと平均と中央値の直線が一致してしまい、分位点回帰のうれしさが見えません。この非対称さが、図で破線（平均）が $\tau = 0.5$ の直線より上に来る理由です。
43-48 行目τ ごとに 1 本ずつ独立に学習する
QuantileRegressor(quantile=tau) を 3 回作り、それぞれが自分の $\tau$ のピンボール損失を最小化します。3 本は互いを知らないので、データが少ないと $\tau = 0.1$ と 0.5 の直線が交差する「分位点交差」が起きることがあります。alpha=0.001 は L1 正則化の強さで、既定の 1.0 では係数が強く縮むため小さくしています。solver="highs" は線形計画法のソルバー指定で、400 点程度なら一瞬で解けます。
50-51 行目比較用の平均直線
LinearRegression は同じデータで条件付き平均を推定します。ノイズの平均は 0 なので真の直線 $1.5x + 5$ に近い位置に来ますが、分位点 0.5 の直線とは裾の分だけずれます。この 2 行を消すと、平均と中央値の差という本ページの要点が図から消えます。
60-69 行目3 本の直線を色分けして描く
axes.prop_cycle から色を取り出し、$\tau$ の順に 1 本ずつ描いています。label_template.format(tau=tau) で凡例に τ の値を入れるので、どの線が 0.1 でどの線が 0.9 かを図だけで判別できます。3 本の間隔が $x$ に沿って一定なら、ノイズの幅が $x$ に依存しない（等分散）と読めます。
80-83 行目各直線の最小・最大を数値で残す
pred.min() と pred.max() は、$x = 0$ と $x = 10$ での各直線の高さです。$\tau = 0.9$ と 0.1 の差が「80% 予測区間の幅」で、この値を印字しておけば図を見なくても区間の広さを比べられます。

QuantileRegressor を使って 0.1・0.5・0.9 分位点を推定し、線形回帰と比較しますの図

実行結果の読み方 #

3 本の実線は点の帯の下端（$\tau = 0.1$）、中央（0.5）、上端（0.9）をそれぞれ通り、どれもほぼ同じ傾き 1.5 で平行に並びます。紫の破線（最小二乗法の平均）は $\tau = 0.5$ の直線より少し上にあり、これはガンマノイズの右に長い裾が平均を引き上げるためです。$\tau = 0.1$ と 0.9 の 2 本で囲まれた帯には観測点の約 80% が入るので、この 2 本をそのまま予測区間として使えます。印字された予測最小値・最大値の差を 3 本で比べると、間隔が等しいかどうかも数値で確認できます。

使いどころと落とし穴 #

需要予測で「欠品を 10% に抑える在庫量」を出す、配送時間の「90% 以内の到着時刻」を示すなど、平均ではなく特定の分位点が意思決定に直結する場面が本領です。ロバスト回帰の代わりに $\tau = 0.5$ を使うのも有効で、外れ値に対して Huber 損失より単純な性質を持ちます。落とし穴は 2 つあります。1 つ目はデータ量で、$\tau = 0.01$ のような極端な分位点は観測点の 1% しか情報を持たないため、400 点では 4 点に支えられた不安定な直線になります。2 つ目は分位点交差で、$\tau$ ごとに独立に学習するため 0.1 の直線が 0.5 を上回る領域が生じることがあり、必要なら予測後に並べ替えます。

分位点と回帰曲線 #

スライダーで分位点 τ を 0.1 から 0.9 まで変えると、回帰直線が点の帯の下端から上端へ移っていく様子を確認できます。

参考文献 #

Koenker, R., & Bassett, G. (1978). Regression Quantiles. Econometrica, 46(1), 33–50.
Koenker, R. (2005). Quantile Regression. Cambridge University Press.

ロバスト回帰 — 外れ値に強い回帰
線形回帰 — 条件付き平均の回帰
ピンボール損失 — 分位点回帰の損失関数