🔖 キーワード索引（追補）

🎨 直感で掴む 📐 ピンボール損失 🧮 実値で計算 🐍 statsmodels 🐍 GBR (sklearn) 🐍 損失検証 🐍 予測区間 ⚠️ 落とし穴 🌐 関連手法 📚 グループ教材

💡 30 秒で分かる結論（追補）

• OLS は条件付き「平均」を、 分位点回帰は条件付き「分位点（中央値・25%/75% など）」を推定する。 名前は「分位点」だが、 実際は「分布のどの部分を捉えるか」を選べる回帰法。
• τ=0.5（中央値回帰）は外れ値に強く、 OLS の代替として有用。 L1 損失と等価で、 LAD（Least Absolute Deviations）と同じ。
• τ=0.05, 0.95 を使えば 予測区間 (90% PI)が直接得られる。 OLS では誤差分布を仮定する必要があるが、 分位点回帰は仮定不要。
• ピンボール損失 $\rho_\tau(u)$ は「過小予測と過大予測を非対称に重み付ける」関数で、 これを最小化することで τ 分位点が得られる。
• SSDSE-B-2026 で「人口 → 出生数」を τ=0.25, 0.5, 0.75 で回帰すると、 傾きが微増する（高人口県では出生数が傾向以上に増える）。 OLS だけでは見えない「分布の裾」が分かる。

📍 文脈ボックス：あなたが今見ているもの

🎨 直感で掴む：「平均回帰」と「分位点回帰」の違い

身長と体重の散布図を想像してほしい。 OLS は「身長 x cm の人の体重の平均」を予測する直線を引く。 分位点回帰 (τ=0.9) は「身長 x cm の人の体重の上位 10% ライン」を引く。 同じデータでも見える線が違う。 SSDSE-B-2026 の「人口 → 出生数」でも同様で、 平均で見ると傾き 0.0061 だが、 上位 25%（τ=0.75）の県では 0.0063 と若干急、 下位 25%（τ=0.25）では 0.0059 と緩い。 大県ほど「出生数のばらつきが平均からポジティブ側に広い」傾向が読み取れる。

📐 ピンボール損失（分位点回帰のコア）

分位点回帰の損失関数 $\rho_\tau(u)$ は「ピンボール（pinball）損失」と呼ばれる非対称な絶対値損失:

$$ \rho_\tau(u) = u \cdot (\tau - \mathbb{1}_{u < 0}) = \begin{cases} \tau \cdot u & \text{if } u \ge 0 \\ (\tau - 1) \cdot u & \text{if } u < 0 \end{cases} $$

$$ \hat{\beta}(\tau) = \arg\min_{\beta} \sum_{i=1}^{n} \rho_\tau(y_i - x_i^\top \beta) $$

🔬 数式を言葉で読み解く

• $u = y - \hat{y}$（残差）。 過小予測 ($u > 0$) には重み $\tau$、 過大予測 ($u < 0$) には重み $1 - \tau$ を掛ける。
• $\tau = 0.5$ なら両側を 等しく重み付け → MAE 最小化 → 中央値が解。
• $\tau = 0.9$ なら過小予測に重み 0.9、 過大予測に重み 0.1 → 過小予測を嫌う → 9 割が予測線以下に来る。
• $\tau = 0.1$ なら過小予測に重み 0.1、 過大予測に重み 0.9 → 過大予測を嫌う → 9 割が予測線以上に来る。
• 係数 $\hat{\beta}(\tau)$ は $\tau$ の関数。 τ を 0.01 から 0.99 まで動かすと、 連続的に変化する係数曲線（quantile process）が得られる。

🧮 実値で計算してみる：SSDSE-B-2026 で 25/50/75 分位回帰

SSDSE-B-2026 (2023 年, 47 県) で、 説明変数 = 人口 (A1101)、 被説明変数 = 出生数 (A4101) として、 OLS と 3 種類の分位回帰を比較する。

🐍 Python 実装

① statsmodels で 3 種類の分位回帰

🎯 このコードでやること: SSDSE-B-2026 の (人口, 出生数) で τ ∈ {0.25, 0.5, 0.75} の分位回帰を statsmodels.formula.api.quantreg で実行し、 OLS と並べて比較。

📥 入力データ (SSDSE-B-2026 2023 年 抜粋):

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15

import pandas as pd
import statsmodels.formula.api as smf

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
d = df[df['SSDSE-B-2026']==2023][['A1101','A4101']]
d.columns = ['pop','birth']

# OLS と分位回帰を 3 つ
ols = smf.ols('birth ~ pop', d).fit()
print(f'OLS: slope={ols.params["pop"]:.6f}, intercept={ols.params["Intercept"]:.2f}')

for tau in [0.25, 0.5, 0.75]:
    qr = smf.quantreg('birth ~ pop', d).fit(q=tau)
    print(f'τ={tau}: slope={qr.params["pop"]:.6f}, '
          f'intercept={qr.params["Intercept"]:.2f}')

📤 実行すると次の出力が得られる:

💬 結果の読み方: τ が増えるほど傾きがおおむね急になる（0.0059 → 0.0063）。 これは「人口の大きい県では出生数の 上限が傾向以上に伸びる」ことを意味する（東京・神奈川の集中効果）。 OLS の傾き 0.0061 はちょうど中央値と平均値の間に来ている。

② sklearn の GradientBoostingRegressor で非線形分位回帰

🎯 このコードでやること: GBR の loss='quantile' オプションで τ=0.1, 0.5, 0.9 の非線形分位回帰を実行。 線形より柔軟に「分布の端」を捉えられる。

📥 入力データ: SSDSE-B-2026 (47 県 × 2023), X=人口、 y=出生数。

 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

import pandas as pd
import numpy as np
from sklearn.ensemble import GradientBoostingRegressor

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
d = df[df['SSDSE-B-2026']==2023]
X = d[['A1101']].values
y = d['A4101'].values

# 3 つの τ で別々にモデルを学習
models = {}
for tau in [0.1, 0.5, 0.9]:
    gbr = GradientBoostingRegressor(
        loss='quantile', alpha=tau,
        n_estimators=100, max_depth=3, random_state=42
    )
    gbr.fit(X, y)
    models[tau] = gbr

# 東京 (人口 14086000) の予測区間
x_tokyo = np.array([[14086000]])
for tau, m in models.items():
    print(f'東京 (人口 14086000) τ={tau}: 予測出生数 = {m.predict(x_tokyo)[0]:,.0f}')

# 平均的県 (人口 200 万) の予測区間
x_avg = np.array([[2000000]])
for tau, m in models.items():
    print(f'仮想県 (人口 2000000) τ={tau}: 予測出生数 = {m.predict(x_avg)[0]:,.0f}')

📤 実行すると次の出力が得られる:

💬 結果の読み方: 東京は人口が突出しており、 訓練データ内で「東京サイズの県は 1 件のみ」のため、 τ=0.5 と τ=0.9 がほぼ同じ値（86,239 と 86,347）になる（過学習）。 仮想県 (人口 200 万) では「下側 10K 〜 上側 12.5K」の予測区間 (約 2300 幅) が得られ、 OLS の単一点予測より「不確実性」を表現できている。 なお τ=0.1 は木モデルの外挿制約で東京・仮想県とも同じ低値 (10,212) を返す点に注意。

③ ピンボール損失を自前で計算 → τ ごとの最小化を確認

🎯 このコードでやること: 「定数予測 c」を様々な値で試し、 ピンボール損失を計算する。 τ=0.5 では中央値、 τ=0.25 では下位 25% 分位点、 τ=0.75 では上位 25% 分位点で損失が最小になることを実値で確認。

📥 入力データ: SSDSE-B-2026 47 県の出生数 A4101 (平均 15,474、 中央値 9,524、 25%=5,472、 75%=14,390)。

 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

import pandas as pd
import numpy as np

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
y = df[df['SSDSE-B-2026']==2023]['A4101'].values  # 出生数

def pinball(y, c, tau):
    u = y - c
    return np.where(u >= 0, tau * u, (tau - 1) * u).mean()

# 候補値: 25%, 中央値, 75%, 平均
candidates = {
    'q25 (5472)': 5472,
    'median (9524)': 9524,
    'q75 (14390)': 14390,
    'mean (15474)': 15474,
}

print('--- τ=0.25 (q25 が最小であるはず) ---')
for name, c in candidates.items():
    print(f'  c={name}: ρ_0.25 = {pinball(y, c, 0.25):.2f}')

print('--- τ=0.50 (中央値が最小であるはず) ---')
for name, c in candidates.items():
    print(f'  c={name}: ρ_0.50 = {pinball(y, c, 0.50):.2f}')

print('--- τ=0.75 (q75 が最小であるはず) ---')
for name, c in candidates.items():
    print(f'  c={name}: ρ_0.75 = {pinball(y, c, 0.75):.2f}')

📤 実行すると次の出力が得られる:

💬 結果の読み方: τ=0.25 のピンボール損失は q25=5472 で最小 (2773)、 τ=0.5 では median=9524 で最小 (4857)、 τ=0.75 では q75=14390 で最小 (5913)。 これが「ピンボール損失最小化 = 分位点推定」の証明的確認。 SSDSE は右裾が長いので、 平均 15474 は中央値 9524 から大きく乖離している。

④ 90% 予測区間の構築（τ=0.05, 0.95）

🎯 このコードでやること: τ=0.05 と τ=0.95 の分位回帰を訓練し、 「90% 予測区間 (PI)」を構築。 OLS の正規分布仮定に依存しない、 ノンパラメトリックな PI が得られる。

📥 入力データ: SSDSE-B-2026 47 県、 X=人口 A1101、 y=出生数 A4101。

 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

import pandas as pd
import numpy as np
from sklearn.ensemble import GradientBoostingRegressor

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
d = df[df['SSDSE-B-2026']==2023]
X = d[['A1101']].values
y = d['A4101'].values

def fit_qr(tau):
    m = GradientBoostingRegressor(
        loss='quantile', alpha=tau,
        n_estimators=200, max_depth=2, random_state=0
    )
    m.fit(X, y)
    return m

m_lo = fit_qr(0.05)
m_md = fit_qr(0.50)
m_hi = fit_qr(0.95)

# 仮想県 (人口 100 万, 300 万, 1000 万) で予測区間
for pop in [1_000_000, 3_000_000, 10_000_000]:
    xq = np.array([[pop]])
    lo, md, hi = m_lo.predict(xq)[0], m_md.predict(xq)[0], m_hi.predict(xq)[0]
    print(f'人口 {pop:>10,}: 90% PI = [{lo:,.0f}, {hi:,.0f}], 中央値={md:,.0f}')

# カバレッジ確認 (訓練データの何割が 90% PI に入るか)
lo_pred = m_lo.predict(X)
hi_pred = m_hi.predict(X)
covered = ((y >= lo_pred) & (y <= hi_pred)).mean()
print(f'\n訓練データのカバレッジ: {covered*100:.1f}% (目標 90%)')

📤 実行すると次の出力が得られる:

💬 結果の読み方: 訓練データの 83.0% が 90% 予測区間に入る → 目標をやや下回るが N=47 の小サンプルでは妥当な範囲。 人口が増えるほど PI の絶対幅は広がる（約 1,800 → 6,800 → 62,000）。 なお τ=0.05 の下限は木モデルの外挿制約で人口 300 万以上では 9,950 に張り付く点に注意。 OLS で同じことをやろうとすると「残差の正規性」を仮定する必要があり、 出生数のような右裾分布では PI が不正確になりがち。 分位点回帰は仮定なしで PI を直接学ぶ。

⚠️ 落とし穴（追補 6 件）

1. 分位点のクロス問題: τ=0.25 の予測が τ=0.5 の予測より大きくなることがある（quantile crossing）。 単調性を保つには cqr ライブラリや QuantileForest を使う。
2. 信頼区間の難しさ: 分位点回帰の係数の標準誤差はブートストラップが基本。 解析解は存在するが、 誤差分布の独立性仮定に依存。 サンプルが小さい時 (n < 100) は注意。
3. サンプルサイズと裾の精度: τ=0.01 や τ=0.99 の極端な分位点は、 サンプル数の少ない裾を推定するため誤差が大きい。 SSDSE 47 件で τ=0.05 の予測は信頼性が低い。
4. 離散変数の問題: y が離散値（出生数の整数値）の場合、 分位点が「同じ値」になることがある（plateau 現象）。 線形補間 (interpolation) で対処。
5. 解釈の罠: 「τ=0.75 の傾きが大きい → 高人口県では出生数の上限が伸びる」と言いたくなるが、 これは「上位 25% の 条件付き分位点」であり「上位 25% の県」ではない。 慎重に言葉を選ぶ必要がある。
6. 「外れ値除去と同じ」という誤解: 分位点回帰は外れ値を 除去しない。 全データを使ったまま、 ピンボール損失が残差の「符号」だけを見る（大きさを二乗しない）ため外れ値 1 点の 影響が有界になる、 というロバスト統計的な仕組み。 外れ値を先に削除（トリミング）すると、 分布の裾情報が失われ τ=0.9 や τ=0.1 の推定そのものが壊れる。 「裾を 捨てる」のではなく「裾を 推定対象にする」のが分位点回帰である。

🌐 関連手法・派生

🔗 関連用語（前提・並列・発展）

前提: 分位点 / OLS / 線形計画 / 損失関数
並列: ロバスト回帰 / GAM / 中央値回帰
発展: 分位点 RF / LightGBM / NGBoost / Conformal Prediction

📚 関連グループ教材

• 回帰分析（グループ教材） — OLS、 ロバスト、 分位点の並列比較
• 機械学習の基礎 — 損失関数の設計と非対称損失
• 不確実性定量化 — 予測区間と Conformal Prediction

📐 補足: 線形計画問題としての分位点回帰

分位点回帰の最適化は線形計画（LP）として定式化できる。 補助変数 $u_i^+, u_i^- \ge 0$ を導入:

$$ \min_{\beta, u^+, u^-} \sum_{i=1}^{n} \left[\tau u_i^+ + (1-\tau) u_i^-\right] $$

$$ \text{s.t.}\quad y_i - x_i^\top \beta = u_i^+ - u_i^-,\ \ u_i^+, u_i^- \ge 0 $$

シンプレックス法や内点法で解ける。 statsmodels の quantreg は内部で内点法（IRLS の派生）を呼んでいる。 n=10,000、 p=50 程度なら 1 秒以内、 n=100,000 でも数秒で解ける。

🐍 Python 実装 (補足): scipy.optimize で線形計画解

🎯 このコードでやること: 上記の LP 定式化を scipy.optimize.linprog で自前で解き、 statsmodels の結果と一致するか確認。

📥 入力データ: SSDSE-B-2026 47 県、 X=人口、 y=出生数 (上記と同じ)。

 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

import pandas as pd
import numpy as np
from scipy.optimize import linprog

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
d = df[df['SSDSE-B-2026']==2023]
y = d['A4101'].values.astype(float)
X = np.column_stack([np.ones(len(d)), d['A1101'].values.astype(float)/1e4])  # [1, pop万人]（well-conditioned）

n, p = X.shape
tau = 0.5  # 中央値回帰

# 変数: [beta (p個), u+ (n個), u- (n個)]
# 目的: minimize sum(tau * u+ + (1-tau) * u-)
c = np.concatenate([np.zeros(p), tau * np.ones(n), (1-tau) * np.ones(n)])

# 制約: X @ beta + u+ - u- = y → [X, I, -I] @ [beta; u+; u-] = y
A_eq = np.hstack([X, np.eye(n), -np.eye(n)])
b_eq = y

# beta は無制約、 u+, u- >= 0
bounds = [(None, None)] * p + [(0, None)] * (2 * n)

res = linprog(c, A_eq=A_eq, b_eq=b_eq, bounds=bounds, method='highs')
beta = res.x[:p]
print(f'LP解 (τ=0.5): intercept={beta[0]:.2f}, slope={beta[1]:.6f}')
print(f'statsmodels比: intercept=-287.53, slope=57.484981 (IRLS 未収束のため LP 解とずれる)')

📤 実行すると次の出力が得られる:

💬 結果の読み方: 自前 LP（HiGHS）は厳密解を返す。 一方 statsmodels の QuantReg は IRLS 反復が上限（max_iter=1000）に達して未収束のまま止まるため、 この n=47 データでは両者が僅かにずれる（ピンボール損失合計: LP 26,872.5 < statsmodels 28,551.0 で LP 解の方が厳密に良い）。 分位点回帰の「本体」は線形計画 1 個に集約されることが確認できる。 OLS が正規方程式 $(X^\top X)^{-1} X^\top y$ の閉形式解を持つのに対し、 分位点回帰は LP の反復解法が必要で計算量は増えるが、 現代の HIGHS / Mosek なら大規模問題でも実用速度。

📐 補足: 「数式を言葉で読み解く」追加版

分位点回帰の条件付き分位点関数 $Q_Y(\tau \mid X = x) = x^\top \beta(\tau)$ の意味を、 もう一段噛み砕く:

• $Q_Y(0.5 \mid X=x) = x^\top \beta(0.5)$ → 「X が x の条件下で Y の 真ん中の値」
• $Q_Y(0.9 \mid X=x) = x^\top \beta(0.9)$ → 「X が x の条件下で Y の 上から 10% 目の値」
• $Q_Y(0.5 \mid X=x) - Q_Y(0.25 \mid X=x)$ → 条件付き IQR の半分（下半分のばらつき幅）
• $Q_Y(0.95 \mid X=x) - Q_Y(0.05 \mid X=x)$ → 条件付き 90% 範囲（予測区間幅）
• $\beta(\tau)$ を τ の関数として描いた図は「coefficient process plot」と呼ばれ、 「説明変数の影響が分布のどこに効くか」を可視化する。

🌐 補足: 業界別ユースケース

📚 確認チェックリスト

• □ OLS と分位点回帰の違い (「平均」vs「分位点」) を 1 文で説明できる
• □ ピンボール損失 $\rho_\tau(u) = u(\tau - \mathbb{1}_{u<0})$ の意味を言葉で読み解ける
• □ τ=0.5 が「中央値 = LAD = L1 損失最小化」と等価であることを理解した
• □ SSDSE-B-2026 で τ=0.25, 0.5, 0.75 の Python コードを動かし、 傾きの推移を確認した
• □ 90% 予測区間 (τ=0.05, 0.95) を構築し、 OLS の正規仮定 PI と比較した
• □ 線形計画 (LP) 定式化と scipy.linprog での解法を確認した
• □ Quantile Crossing 問題と Quantile Random Forest の役割を理解した
• □ 自分の業務で τ をどう選ぶか（VaR=0.01、 SLA=0.95 など）を言える

📐 補足: 分位点回帰 vs Conformal Prediction

両者とも「予測区間」を作る手法だが、 アプローチが全く異なる:

実務では「分位点回帰で適応的な区間を作り、 Conformal で有限サンプル保証を付与する」CQR が主流になりつつある (Romano et al., 2019)。

📚 さらに読むには

• Roger Koenker & Gilbert Bassett, "Regression Quantiles" (Econometrica 1978) — 原典
• Roger Koenker, Quantile Regression (Cambridge University Press, 2005) — 教科書
• Nicolai Meinshausen, "Quantile Regression Forests" (JMLR 2006)
• Yarin Romano et al., "Conformalized Quantile Regression" (NeurIPS 2019)
• Tony Duan et al., "NGBoost: Natural Gradient Boosting for Probabilistic Prediction" (ICML 2020)

🧮 実値で計算 (補足): SSDSE 全 12 年で τ-process プロット相当

SSDSE-B-2026 を 2023〜2012 年（12 年分）に拡張し、 各年で 47 県の (人口, 出生数) に τ=0.1, 0.25, 0.5, 0.75, 0.9 の分位回帰を当てる。 傾き $\beta_1(\tau)$ の年次推移を表で示す。

観察 1: 年が進むにつれて傾き全体が低下 (0.0083 → 0.0061) → 少子化進行で「人口あたり出生数」が全国的に下がっている。
観察 2: 各年内で τ=0.9 / τ=0.1 比 = 約 1.13 倍 → 高分位（高人口県）の方が傾きが急という構造は不変。
観察 3: OLS の傾きはおおむね τ=0.5〜0.75 の間 → 出生数分布の右裾の影響を平均が引っ張られている。

📐 補足: 多変量分位点回帰

複数の説明変数を含めることも可能。 SSDSE-B-2026 で X = {人口 A1101, 年平均気温 B4101, 婚姻件数 A9101} を使い y = 出生数を推定する例:

$$ Q_Y(\tau \mid X) = \beta_0(\tau) + \beta_1(\tau) \cdot \text{pop} + \beta_2(\tau) \cdot \text{temp} + \beta_3(\tau) \cdot \text{marr} $$

係数 $\beta_j(\tau)$ は τ ごとに変わるので、 「婚姻件数が多い県では中央値 (τ=0.5) の出生数は増えるが、 上位 (τ=0.9) では影響が変わる」といった分布依存の関係を抽出できる。 OLS では平均 1 本しか得られないので、 こうした構造は見落とされる。

⚠️ 落とし穴 (補足 5 件)

1. τ の選択ミス: 「とりあえず 0.25, 0.5, 0.75」では業務目的に合わないことが多い。 VaR は τ=0.01, SLA は τ=0.99 など、 「何を保証したいか」から逆算する。
2. カバレッジの誤読: 「90% PI を作ったから 90% 当たる」は周辺的にしか保証されない。 X の特定領域（外挿）ではカバレッジが落ちる。 必ず検証データで 条件付きカバレッジを確認する。
3. τ=0.5 を OLS の代替に使う罠: 中央値回帰は確かに外れ値に強いが、 「平均」を直接知りたい時には誤った結論を導く。 例: 「平均収入」を中央値で報告すると過小評価。
4. 非線形関係の見落とし: 線形分位点回帰は直線しか引けない。 SSDSE のように「東京が突出した非線形性」がある場合、 線形 QR は東京で予測が外れる。 GBR-quantile か Quantile RF を使う。
5. 解釈不可能になりがち: τ ごとに係数が変わるので「人口の影響」を一文で説明できない。 QR は「分布のどこに効くか」を分解する手法であり、 単一の effect size を返さない。 これは欠点ではなく特徴。

🐍 Python 実装 (補足 2): Quantile Crossing の検出と修正

🎯 このコードでやること: SSDSE-B-2026 で τ=0.1, 0.5, 0.9 を別々に推定し、 「クロスしている点」を検出。 単調性を保つ後処理（sorting）を適用。

📥 入力データ: SSDSE-B-2026 47 県、 X=人口、 y=出生数。

 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

import pandas as pd
import numpy as np
import statsmodels.formula.api as smf

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
d = df[df['SSDSE-B-2026']==2023][['A1101','A4101']].copy()
d.columns = ['pop','birth']

# 3 つの τ で線形分位回帰
preds = {}
for tau in [0.1, 0.5, 0.9]:
    qr = smf.quantreg('birth ~ pop', d).fit(q=tau)
    preds[tau] = qr.predict(d)

pred_df = pd.DataFrame(preds)
print('--- 推定結果 (最初の 5 行) ---')
print(pred_df.head())

# クロス検出: τ=0.1 > τ=0.5 や τ=0.5 > τ=0.9 となる行
cross = ((pred_df[0.1] > pred_df[0.5]) | (pred_df[0.5] > pred_df[0.9])).sum()
print(f'\nクロスしている件数: {cross} / {len(d)}')

# 修正: 各行で値をソート
fixed = pred_df.apply(lambda row: pd.Series(sorted(row.values),
                                            index=[0.1, 0.5, 0.9]), axis=1)
print('\n--- ソート後 (クロス解消) ---')
print(fixed.head())

📤 実行すると次の出力が得られる:

💬 結果の読み方: SSDSE のような直線関係が強いデータではクロスはほぼ起きない（0 件）。 一方、 非線形なデータ（医療コスト等）ではクロスが頻発し、 sorting 修正が必須。 より洗練された手法は Chernozhukov et al. (2010) の「rearrangement」で、 連続関数として単調化する。

📐 補足: 推定量の漸近分布

Koenker (2005) によれば、 分位点回帰推定量 $\hat{\beta}(\tau)$ は次の漸近分布を持つ:

$$ \sqrt{n}(\hat{\beta}(\tau) - \beta(\tau)) \xrightarrow{d} N\left(0,\ \tau(1-\tau) \cdot D^{-1} \cdot J \cdot D^{-1}\right) $$

ここで $x_i$ は説明変数ベクトル（計画行列 $X$ の第 $i$ 行）で、 $D = E\left[f_{Y \mid X}\bigl(Q_\tau(Y \mid X) \mid X\bigr) \, x_i x_i^\top\right]$、 $J = E[x_i x_i^\top]$。 $D$ の推定には条件付き密度 $f_{Y \mid X}$ が必要で、 これがブートストラップが好まれる理由でもある。

🐍 Python 実装 (補足 3): ブートストラップで信頼区間

🎯 このコードでやること: τ=0.5 の傾き係数の 95% 信頼区間を、 1000 回のブートストラップで推定。 statsmodels が返すデフォルトの「IID 漸近 SE」と比較。

📥 入力データ: SSDSE-B-2026 47 県、 (人口, 出生数)。

 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

import pandas as pd
import numpy as np
import statsmodels.formula.api as smf

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
d = df[df['SSDSE-B-2026']==2023][['A1101','A4101']].copy()
d.columns = ['pop','birth']

# ブートストラップ 1000 回
B = 1000
slopes = []
for b in range(B):
    sample = d.sample(n=len(d), replace=True, random_state=b)
    qr = smf.quantreg('birth ~ pop', sample).fit(q=0.5)
    slopes.append(qr.params['pop'])
slopes = np.array(slopes)

# 95% 信頼区間
ci_lo, ci_hi = np.percentile(slopes, [2.5, 97.5])
print(f'ブートストラップ平均: {slopes.mean():.6f}')
print(f'ブートストラップ SE: {slopes.std():.6f}')
print(f'95% CI (bootstrap): [{ci_lo:.6f}, {ci_hi:.6f}]')

# statsmodels の漸近 SE と比較
qr_full = smf.quantreg('birth ~ pop', d).fit(q=0.5)
print(f'\nstatsmodels: slope={qr_full.params["pop"]:.6f}, '
      f'SE={qr_full.bse["pop"]:.6f}')
print(f'statsmodels 95% CI: [{qr_full.conf_int().loc["pop", 0]:.6f}, '
      f'{qr_full.conf_int().loc["pop", 1]:.6f}]')

📤 実行すると次の出力が得られる:

💬 結果の読み方: ブートストラップ CI [0.00547, 0.00628] は statsmodels の漸近 CI [0.00558, 0.00592] より明らかに広い。 これは「右裾の長い分布」「N=47 の小サンプル」で漸近近似が楽観的になっているため。 一般に分位点回帰では ブートストラップ CI を信用すべき。

🗺 概念マップ (補足): 損失関数の系譜

🌐 補足: 分位点回帰の歴史

• 1760 年: Boscovich が「絶対値残差最小化」(L1) を初めて提案 → 中央値回帰の原型
• 1789 年: Laplace が L1 推定量を体系化
• 1978 年: Koenker & Bassett, Econometrica で「Regression Quantiles」を発表 → 現代的定式化
• 1990 年代: 経済学（賃金の分位点）、 医療（成長曲線）で標準化
• 2006 年: Meinshausen, Quantile Random Forest で非線形対応
• 2017 年: LightGBM が objective='quantile' を実装、 GBDT 主流に
• 2019 年: CQR (Conformalized Quantile Regression) で有限サンプル保証付き予測区間
• 2020 年〜: 深層学習版（QuantileNeuralNet, DeepAR）で時系列・大規模適用

🐍 Python 実装 (補足 4): LightGBM の quantile regression

🎯 このコードでやること: LightGBM の objective='quantile' で SSDSE-B-2026 の多変量分位回帰。 説明変数 3 つ（人口・気温・婚姻件数）で出生数の τ=0.1, 0.5, 0.9 を推定。

📥 入力データ: SSDSE-B-2026 47 県、 X=[A1101, B4101, A9101]、 y=A4101。

 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

import pandas as pd
import numpy as np
import lightgbm as lgb

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
d = df[df['SSDSE-B-2026']==2023].dropna(subset=['A1101','B4101','A9101','A4101'])
X = d[['A1101', 'B4101', 'A9101']].values
y = d['A4101'].values

# 3 つの τ で別々に学習
preds = {}
for tau in [0.1, 0.5, 0.9]:
    model = lgb.LGBMRegressor(
        objective='quantile', alpha=tau,
        n_estimators=100, max_depth=3,
        learning_rate=0.05, random_state=42, verbose=-1
    )
    model.fit(X, y)
    preds[tau] = model.predict(X)

# 訓練データ上のカバレッジ (90% PI = τ=0.1〜0.9)
covered = ((y >= preds[0.1]) & (y <= preds[0.9])).mean()
print(f'90% PI カバレッジ (訓練): {covered*100:.1f}%')

# 仮想入力（人口 500 万、 気温 15度、 婚姻件数 2 万件）の予測区間
x_new = np.array([[5_000_000, 15.0, 20000]])
for tau in [0.1, 0.5, 0.9]:
    model = lgb.LGBMRegressor(
        objective='quantile', alpha=tau,
        n_estimators=100, max_depth=3,
        learning_rate=0.05, random_state=42, verbose=-1
    )
    model.fit(X, y)
    print(f'τ={tau}: 予測 = {model.predict(x_new)[0]:,.0f}')

📤 実行すると次の出力が得られる:

💬 結果の読み方: LightGBM-quantile は非線形・多変量に対応し、 訓練データではカバレッジ 80.9% と名目 80%（τ=0.1〜0.9）にほぼ一致。 実務では検証データでチューニング (alpha, n_estimators) が必要。 仮想入力の PI 幅 [10K, 52K] は OLS では得られない不確実性表現で、 「人口だけ」より精度が上がる。

📐 補足: Expectile 回帰との対比

Expectile（期待分位点）は 1987 年に Newey & Powell が提案した別の非対称損失。 ピンボール損失の代わりに「非対称二乗損失」を使う:

$$ \rho^E_\tau(u) = u^2 \cdot |\tau - \mathbb{1}_{u<0}| = \begin{cases} \tau u^2 & \text{if } u \ge 0 \\ (1-\tau) u^2 & \text{if } u < 0 \end{cases} $$

Basel III では 2017 年から VaR から ES への移行が進んでおり、 Expectile 回帰の重要性が増している。

🌐 補足: モデル評価指標

分位点回帰の予測精度は「平均二乗誤差 (MSE)」では測れない（OLS 用）。 適切な指標は:

• Pinball Loss (τ ごと): $\frac{1}{n}\sum \rho_\tau(y_i - \hat{q}_i^\tau)$ → 小さいほど良い、 τ ごとに別評価
• Coverage: $\frac{1}{n}\sum \mathbb{1}\{y_i \le \hat{q}_i^\tau\}$ → 理想は τ に近い値（例: τ=0.9 なら 90%）
• Interval Score: 予測区間の幅 + ペナルティ (区間外のサンプルに重い罰)
• CRPS: 全ての τ を統合した「分布ごと評価」、 確率予測の標準指標
• Sharpness: 予測区間の幅。 狭い方が良いが、 カバレッジを犠牲にしてはダメ

📐 補足: 分位点回帰のハイパーパラメータ早見表

🔬 数式を言葉で読み解く (補足): なぜピンボール損失で分位点が出るか

ピンボール損失の期待値 $E[\rho_\tau(Y - c)]$ を $c$ で微分してゼロと置くと:

$$ \frac{\partial}{\partial c} E[\rho_\tau(Y - c)] = -\tau \cdot P(Y \ge c) + (1-\tau) \cdot P(Y < c) = 0 $$

$$ \implies P(Y < c) = \tau \implies c = Q_\tau(Y) $$

つまり「ピンボール損失の期待値を最小化する $c$ は、 分布 $Y$ の τ 分位点に他ならない」が証明される。 これが分位点回帰の数学的核心。 同じロジックで:

• $L_2$ 損失 → 平均（$\frac{\partial}{\partial c} E[(Y-c)^2] = -2E[Y-c] = 0 \Rightarrow c = E[Y]$）
• $L_1$ 損失 → 中央値（ピンボールの τ=0.5 ケース）
• 非対称 L2 → expectile

「損失関数を決めれば、 分布のどの点を予測するかが決まる」。 これは予測モデリングの最も深い原理の 1 つ。

🧮 実値 (補足 2): 損失関数比較表

SSDSE-B-2026 47 県の出生数を「定数 c で予測する」場合、 c の値を変えながら異なる損失を比較:

★ は各列で最小値。 L2 最小 = 平均、 L1 最小 = 中央値、 ρ_0.25 最小 = q25、 ρ_0.75 最小 = q75、 ρ_0.9 最小 = q90 付近 (c=40,000。 実測 q90=38,238)。 「損失と推定対象の対応」が一目で分かる。 SSDSE は右裾が長いので、 平均 (15K) と中央値 (9.5K) が大きく乖離している点に注意。

🐍 Python 実装 (補足 5): カバレッジ・キャリブレーション図

🎯 このコードでやること: τ ∈ {0.1, 0.2, ..., 0.9} の 9 つの分位回帰を学習し、 検証データ上で「予測されたカバレッジ」vs「目標カバレッジ」をテーブル化。 ideally 一致するはず。

📥 入力データ: SSDSE-B-2026 全 12 年 × 47 県 = 564 件、 X=人口、 y=出生数。 ホールドアウト 80/20。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20

import pandas as pd
import numpy as np
import statsmodels.formula.api as smf
from sklearn.model_selection import train_test_split

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
d = df[['A1101','A4101']].dropna()  # 全 12 年分をプールしテスト標本数を確保
d.columns = ['pop','birth']
train, test = train_test_split(d, test_size=0.2, random_state=42)

results = []
for tau in np.arange(0.1, 1.0, 0.1):
    qr = smf.quantreg('birth ~ pop', train).fit(q=tau)
    pred = qr.predict(test)
    coverage = (test['birth'] <= pred).mean()
    results.append({'target_τ': round(tau, 1),
                    'actual_coverage': round(coverage, 3),
                    'gap': round(coverage - tau, 3)})

print(pd.DataFrame(results).to_string(index=False))

📤 実行すると次の出力が得られる:

💬 結果の読み方: 目標 τ と実測カバレッジが多くの点で ±0.02 前後（最大でも gap +0.087 @τ=0.4）に収まり、 分位回帰は おおむねキャリブレーション済み。 OLS + 正規分布仮定の PI は右裾の長いデータでカバレッジが落ちることがあるが、 分位回帰は分布仮定なしで近い水準を保つ。 これが「OLS の代替」として分位回帰が推される最大の理由。

📚 もう一段の検討事項

1. τ を 0.01 間隔で 99 個推定する場合、 LP を 99 回解くのは時間がかかる。 → "Single-pass quantile regression" (Bondell et al., 2010) で同時推定可能。
2. パネルデータ（県 × 年）では固定効果分位点回帰 (Koenker, 2004) が必要。 静的 QR では年内変動が捉えられない。
3. 因果推論との接続: Conditional Quantile Treatment Effect (CQTE) で「処置が分布のどこに効くか」を推定可能。 RCT 不平等分析の標準。
4. ベイズ版: ALD (Asymmetric Laplace Distribution) で尤度を構築するベイズ QR (Yu & Moyeed, 2001)。 事前分布で正則化が自然。
5. 関数データ・縦断データへの拡張: Functional Quantile Regression、 縦断 QR で時系列・空間データ対応。

🌐 補足: 「OLS と QR をどう使い分けるか」フローチャート

📐 補足: 分位点回帰の正則化

高次元（説明変数が多い）QR では正則化が必要。 Belloni & Chernozhukov (2011) の L1-penalized QR:

$$ \hat{\beta}(\tau) = \arg\min_{\beta} \sum_{i=1}^{n} \rho_\tau(y_i - x_i^\top \beta) + \lambda \|\beta\|_1 $$

Lasso と同じく L1 罰則で sparse な係数を得る。 SCAD、 adaptive lasso も適用可能。 Python では statsmodels 単体では未サポートだが、 sklearn-quantile や pyqreg で対応。

📚 さらに深く学ぶには (補足)

• Roger Koenker, "Quantile Regression for Longitudinal Data" (J. Multivar. Anal. 2004)
• Victor Chernozhukov et al., "Quantile and Probability Curves Without Crossing" (Econometrica 2010)
• Alexandre Belloni & Victor Chernozhukov, "L1-Penalized Quantile Regression in High-Dimensional Sparse Models" (Ann. Stat. 2011)
• Keming Yu & Rana A. Moyeed, "Bayesian Quantile Regression" (Stat. & Prob. Lett. 2001)
• Howard Bondell et al., "Noncrossing Quantile Regression Curve Estimation" (Biometrika 2010)

🎨 直感 (補足): 「弓道の的」の比喩

OLS は「的の中心 (黒丸) を狙う矢」と同じ。 ピンボール損失 (τ=0.5) は「的の中央線 (黒丸を通る左右対称線) を狙う」。 τ=0.9 は「的の上から 1 割の位置 (高得点ゾーン) を狙う」。 矢の散らばり方が同じでも、 狙う点が違えば「外れ方の罰則」が違う。 これがピンボール損失の非対称性の本質。

SSDSE で言えば、 「出生数の平均」は的の中央 (15K) だが、 中央値 (9.5K) は的の左寄り、 q75 (14K) は的の右寄り、 q90 (38K) は的の外周。 同じ 47 県データから 5 つの違う「狙い」を導ける。

📈 OLS と 5 つの τ-回帰を 47 都道府県で同時推定

分位点回帰の真価は、 OLS が描く 1 本の回帰線では捉えられない「条件付き分布の形状変化」を可視化できる点にあります。 SSDSE-B-2026 の 人口総数 を説明変数、 出生数 を目的変数として、 τ ∈ {0.10, 0.25, 0.50, 0.75, 0.90} の 5 本の回帰直線を引くと、 大都市圏 (人口総数が右側) で 傾きが拡大することが見えます。 つまり「人口あたり出生数」のばらつきは大都市ほど大きいという事実が、 1 本の OLS では完全に潰れて見えなくなります。

🐍 statsmodels.QuantReg で τ を 5 段階で推定

このコードでやること: SSDSE-B-2026 の data/raw/SSDSE-B-2026.csv から人口と出生数を抽出し、 statsmodels.regression.quantile_regression.QuantReg で 5 本の τ-回帰直線を一気に推定する。 同じ条件で OLS も並走させ、 傾き係数の差を比較する。

📥 入力データ (SSDSE-B-2026 の人口と出生数、 47 行):

1
2
3
4
5
6
7
8
9

🐍 quantreg_5taus.py<button onclick="copyTryitCode(this)">📋 コピー</button></div><div class="pygblock">import pandas as pd
import statsmodels.formula.api as smf
df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
df = df[df['SSDSE-B-2026'] == 2023].rename(columns={'A1101': 'pop', 'A4101': 'births'})
ols = smf.ols('births ~ pop', data=df).fit()
print(f'OLS         slope = {ols.params[1]:.5f}')
for q in [0.10, 0.25, 0.50, 0.75, 0.90]:
    qr = smf.quantreg('births ~ pop', data=df).fit(q=q)
    print(f'tau={q:.2f}    slope = {qr.params[1]:.5f}')</div>

📤 実行すると次の出力が得られる:

💬 τ が大きくなるほど傾きがおおむね増加（0.005695 → 0.006473、 約 1.14 倍）しています。 これは「人口の増加に伴い、 上位 10% の県では出生数の伸びが下位 10% の県よりやや急峻」であることを意味します。 OLS は中央付近 (τ=0.50 付近) の値しか返さず、 この 分布の歪みは永久に見えません。 政策決定では「τ=0.10 の県 (下位群)」をターゲットにすべきか「τ=0.90 の県 (上位群)」をターゲットにすべきかで対策が分かれるため、 分位点回帰が必須となります。

🎨 直感で掴む — 「真ん中」 vs 「上限・下限」

例え話：通勤時間の予測

あなたは「家から駅まで何分かかるか」を予測したいとします。 普通の回帰なら「平均 12 分」と返してきます。 でも、 本当に知りたいのは：

• 「電車に遅れずに着くには何分前に出ればよいか？」 → 90% 分位点（最悪寄りの値）。 例：17 分
• 「だいたい何分かかる？」 → 中央値（50% 分位点）。 例：11 分
• 「最短はどれくらい？」 → 10% 分位点。 例：8 分

これらはすべて「同じデータから抽出できる別々の予測値」です。 分位点回帰は 3 本まとめて出してくれる 強力なツール。

不等分散（heteroskedastic）データに特に強い

「所得が高いほど消費のばらつきも大きい」「人口の多い県ほど延べ宿泊者数のばらつきも大きい」 のような 不等分散 データでは、 通常の OLS が引いた 1 本の直線は 「真ん中だけ」 しか見せてくれません。

OLS との対比

🎨 概念図で押さえる分位点回帰の本質（補遺）

分位点回帰の中核を 3 枚の図で再整理する。 平均回帰との違い、 ピンボール損失の非対称性、 τ ごとの回帰直線の束という 3 つの視点を並べる。

🖼 図 1: 最小二乗（平均）vs 分位点回帰（中央値・上下端）

最小二乗回帰は 条件付き平均 を引く線、 分位点回帰は 条件付き τ 分位点 を引く線である。 外れ値や非対称分布があると、 両者は系統的にずれる。

💬 読み方：赤の OLS 直線は上方の外れ値（赤丸）に引っ張られて中央値直線（青破線, τ=0.5）から離れる。 オレンジ（τ=0.9）と緑（τ=0.1）の 2 本は 予測区間の上下端 を直接与える。

🖼 図 2: ピンボール損失の非対称性

分位点回帰は ρ_τ(u) = u·(τ − 𝟙(u<0)) という非対称損失（ピンボール損失）を最小化する。 τ=0.5 のときは対称の絶対値損失、 τ≠0.5 のときは過小予測と過大予測の罰が異なる。

💬 読み方：青の τ=0.5 は左右対称の V 字（絶対値）。 オレンジ τ=0.9 は 右側（過小予測）の傾きが急 で、 つまり「予測 ŷ が真値 y より小さい状態」を強く罰する → 上方に寄った推定線になる。 緑 τ=0.1 はその逆。

🖼 図 3: τ を変えた予測区間バンド

τ を複数（例: 0.1, 0.5, 0.9）で並行推定すると、 条件付き予測区間が x の関数として自然に得られる。 平均回帰＋等分散仮定の信頼区間と違い、 ヘテロスケダスティック（x で分散が変わる）な状況にも対応できる。

💬 読み方：水色の帯は τ=0.1〜0.9 の 80% 予測区間。 x が大きくなるにつれ帯が広がる場合、 条件付き分散の不均一性 を分位点回帰がそのまま捉えていることが分かる。 OLS の信頼区間（等幅）では再現できない情報である。

以上、 分位点回帰の本質（中央値推定・ピンボール非対称損失・予測区間）を 3 枚の図で整理した。 詳細は本文の各セクションを参照。

🎮 触って理解する

下の散布図は 教材用の擬似データ（x が大きいほど y のばらつきが増える 不均一分散 データ）。 スライダーで分位点 τ を動かすと、 その分位点を通る回帰直線が リアルタイム に引き直される。 散布図を左右にドラッグ（スマホは左右スワイプ）しても τ が変わる。 計算は「最適な分位点回帰直線は少なくとも 2 点を通る」という性質を使い、 全点対を厳密に探索して大域最適解を求めている（近似ではない）。

👀 ここに注目

• τ=0.5（中央値回帰）と OLS の違い：両方を表示して比較する。 データが右上に扇形へ広がると、 OLS（赤破線）は上側に散らばる点へ二乗誤差で引っ張られ、 中央値線（紫）と傾きがずれる。 これが「平均の回帰」と「中央値の回帰」の差。
• τ=0.9 と τ=0.1 で傾きが違う：「同時表示」をオン。 上側（オレンジ, τ=0.9）ほど傾きが急で、 下側（緑, τ=0.1）ほど緩やかなら、 x が増えるほど上位が急伸する 不均一分散のサイン。 3 本が平行に近ければ等分散。 「ばらつきの不均一さ」を 0 にすると 3 本はほぼ平行になる。
• ピンボール損失の非対称性：左下の小図。 τ を上げると「予測が小さすぎる側（u>0）」の傾きが急になり、 その罰を避けるため回帰線が上方へ引き上げられる。 τ=0.5 では左右対称の V 字（絶対値＝L1 損失）になる。

📐 数式 — Koenker & Bassett (1978) の定式化

分位点の定義

pinball loss（チェック関数）

分位点回帰の最適化問題

OLS との対比：損失関数の違い

🔬 数式を「言葉」で読み解く

🧮 SSDSE-B 都道府県データで計算してみる

SSDSE-B 2026 の 2023 年データを用い、 「延べ宿泊者数 $y$ ~ 人口総数 $x$」 の分位点回帰を τ = 0.10, 0.50, 0.90 の 3 本で推定します。

STEP 0: 何が知りたいか

普通の回帰なら 1 本の直線で「人口が 1 万人増えると延べ宿泊者数は何万人泊増えるか（平均）」を見ます。 しかし都道府県では 大都市・観光地と地方で 1 人あたり宿泊者数が大きく違う ので、 「平均だけ」では情報を失います。

STEP 1: 数値の準備（標本値の代表例）

説明変数 $x$ = 人口総数（万人）、 被説明変数 $y$ = 延べ宿泊者数（万人泊）。 5 県を抜粋：

STEP 2: 3 本の分位点回帰直線（47 県全部で推定した典型結果）

STEP 3: 解釈

• 傾きが τ で増加（2.13 → 4.08 → 5.57） → 「上位県ほど、 人口 1 万人あたりの宿泊者数の伸びが大きい」。 観光需要の集中が上位県で強い ことを示唆。
• 3 本の 幅（90% − 10%）が右ほど広がる → 不等分散の証拠。 大都市ほど予測の不確実性が大きい。
• τ=0.50 と τ=0.10/0.90 を比較すると、 OLS（平均直線）は傾き 4.11 で τ=0.50（4.08）とほぼ重なる（東京という高レバレッジ点の影響を受けつつも、 傾き自体は中央値付近に落ち着く）。

🧮 ピンボール損失 (Pinball Loss) の幾何学

分位点回帰が最小化する目的関数は ピンボール損失 (チェック損失とも呼ぶ) です。 τ ∈ (0, 1) のとき、 残差 u = y − ŷ に対し:

$$ \rho_\tau(u) = \begin{cases} \tau \cdot u & (u \geq 0) \\ -(1-\tau) \cdot u & (u < 0) \end{cases} $$

これは「予測が下振れ (u>0) したときのペナルティ τ」と「上振れ (u<0) したときのペナルティ 1-τ」を非対称に課す関数で、 τ=0.5 のとき左右対称となり中央値回帰 (LAD 回帰) と一致します。 τ=0.9 なら下振れ予測を 9 倍重く罰し、 結果として「90% 分位」の上を通る線が得られます。

🐍 ピンボール損失を手書きして 5 本の τ-線を最小化

このコードでやること: sklearn.linear_model.QuantileRegressor を用いて、 ピンボール損失の数値最小化を直接行う。 statsmodels の内部実装 (LP based) と異なり scikit-learn は L1 正則化と組み合わせて使うことで スパース分位点回帰も実現できる。

📥 入力データ (47 県の人口 → 延べ宿泊者数):

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

🐍 sklearn_quantreg_pinball.py<button onclick="copyTryitCode(this)">📋 コピー</button></div><div class="pygblock">import pandas as pd
from sklearn.linear_model import QuantileRegressor
from sklearn.metrics import mean_pinball_loss
df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
df = df[df['SSDSE-B-2026'] == 2023]
X = (df[['A1101']] / 1e4).values
y = (df['G7101'] / 1e4).values
for q in [0.10, 0.50, 0.90]:
    m = QuantileRegressor(quantile=q, alpha=0.0, solver='highs').fit(X, y)
    loss = mean_pinball_loss(y, m.predict(X), alpha=q)
    print(f'tau={q:.2f}  pinball_loss={loss:7.2f}  coef={m.coef_.round(4)}')</div>

📤 実行すると次の出力が得られる:

💬 ピンボール損失は τ=0.5 で最大、 両端 (τ=0.1, 0.9) で小さくなります。 これは「中央値推定は全データの絶対値偏差和を最小化するので大きく、 両端は片側 10% だけ責任を持つので小さい」という幾何学的な帰結です。 業務応用では、 例えば 需要予測の安全在庫 (τ=0.95) や VaR 推定 (τ=0.01) など、 τ を業務 KPI に合わせて設計します。

🧮 数式に値を入れて手で計算する: 分位点回帰の損失

合成データで τ=0.5 (中央値回帰) の損失を計算する。

Step 1: 損失関数

Step 2: 5 点で計算

Step 3: 損失合計

🐍 Python で再現

1
2
3
4
5

import numpy as np
e = np.array([2, 1, 0, -1, -2])
def rho(u, tau): return u * (tau - (u < 0).astype(float))
print(f"τ=0.5: {rho(e, 0.5).sum()}")
print(f"τ=0.9: {rho(e, 0.9).sum()}")

📤 実行結果

💬 手計算 (Step 3) と Python 出力が完全一致。

🐍 Python 実装 — statsmodels と LightGBM の 2 通り

1. statsmodels (古典的・線形分位点回帰)

🎯 このコードでやること：statsmodels の quantreg で SSDSE-B 2026 都道府県人口→延べ宿泊者数の関係を τ=0.10/0.50/0.90 の 3 分位点で線形回帰。 中央値 (50%) と下位 10%・上位 90% の回帰係数の違いから「散らばりが人口で広がるか」を読む。

📥 入力：SSDSE-B-2026.csv の 2023 年データ (47 都道府県)。 列 A1101=人口 → 万人単位 / G7101=延べ宿泊者数 → 万人泊単位に換算。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22

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

# SSDSE-B 2026 を読み込み、2023 年に絞る
df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1]).rename(columns={'SSDSE-B-2026': '年度'})
df = df[df['年度'] == 2023].copy()
df = df.rename(columns={'A1101': '人口', 'G7101': '延べ宿泊者数'})

# 単位を万人・万人泊に揃える（数値の桁を整える）
df['人口万人'] = df['人口'] / 1e4
df['宿泊万人泊'] = df['延べ宿泊者数'] / 1e4

# 3 つの分位点で回帰
taus = [0.10, 0.50, 0.90]
models = {}
for tau in taus:
    mod = smf.quantreg('宿泊万人泊 ~ 人口万人', data=df)
    res = mod.fit(q=tau)
    models[tau] = res
    print(f'τ={tau}: intercept={res.params[0]:.2f}, slope={res.params[1]:.3f}')

📤 実行例：

💬 結果の読み方：slope (人口万人 → 延べ宿泊者数 万人泊) が τ=0.10 で 2.13、 τ=0.90 で 5.57 と増加 → 人口が同じでも上位 (90 パーセンタイル) ほど人口あたりの宿泊者数が多い。 OLS なら 1 本の回帰直線だが、 分位点回帰は「散らばりが人口に依存して扇形に広がる」構造を可視化できる。

2. 結果の可視化（扇形に広がる予測区間）

 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

fig, ax = plt.subplots(figsize=(9, 6))

# 散布図
ax.scatter(df['人口万人'], df['宿泊万人泊'],
           s=70, alpha=0.6, edgecolor='white', linewidth=1.2,
           color='#1565C0', label='47 都道府県 (2023)')

# 3 本の分位点直線
x_grid = np.linspace(df['人口万人'].min(), df['人口万人'].max(), 100)
colors = {0.10: '#D32F2F', 0.50: '#388E3C', 0.90: '#D32F2F'}
styles = {0.10: '--', 0.50: '-', 0.90: '--'}
for tau, res in models.items():
    a, b = res.params
    ax.plot(x_grid, a + b * x_grid,
            color=colors[tau], linestyle=styles[tau], linewidth=2.2,
            label=f'τ = {tau:.2f}')

# 比較のため OLS も
import statsmodels.formula.api as smf
ols = smf.ols('宿泊万人泊 ~ 人口万人', data=df).fit()
ax.plot(x_grid, ols.params[0] + ols.params[1] * x_grid,
        color='#1976D2', linewidth=2.2, linestyle=':', label='OLS (平均)')

ax.set_xlabel('人口総数（万人）')
ax.set_ylabel('延べ宿泊者数（万人泊）')
ax.set_title('分位点回帰：延べ宿泊者数 ~ 人口（10/50/90 パーセンタイル）')
ax.legend(loc='upper left')
plt.tight_layout()
plt.show()

3. LightGBM での非線形分位点回帰

大きなデータや非線形関係には、 LightGBM の objective='quantile' が高速で便利です。

 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

from lightgbm import LGBMRegressor

# 同じデータで 3 本の分位点モデルを学習
X = df[['人口万人']].values
y = df['宿泊万人泊'].values

quantile_models = {}
for tau in [0.10, 0.50, 0.90]:
    m = LGBMRegressor(objective='quantile', alpha=tau,
                       n_estimators=300, learning_rate=0.05,
                       num_leaves=15, random_state=42, verbose=-1)
    m.fit(X, y)
    quantile_models[tau] = m

# 各 τ で予測（同じ x_grid 上）
X_grid = x_grid.reshape(-1, 1)
preds = {tau: m.predict(X_grid) for tau, m in quantile_models.items()}

# 予測区間として fill_between
fig, ax = plt.subplots(figsize=(9, 6))
ax.fill_between(x_grid, preds[0.10], preds[0.90],
                alpha=0.18, color='#D32F2F', label='80% 予測区間')
ax.plot(x_grid, preds[0.50], color='#388E3C', linewidth=2.3, label='中央値')
ax.scatter(df['人口万人'], df['宿泊万人泊'], s=60, alpha=0.6, color='#1565C0')
ax.set_xlabel('人口（万人）'); ax.set_ylabel('延べ宿泊者数（万人泊）')
ax.legend(); plt.tight_layout(); plt.show()

4. pinball loss を直接計算して評価

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

# ── この抜粋で使うモデルを用意します ──
import pandas as pd
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.metrics import mean_pinball_loss

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
df = df[df['SSDSE-B-2026'] == 2023]          # 2023 年の 47 都道府県
X = df[['A1101', 'A1303']].astype(float).values
y = df['A4101'].astype(float).values          # 出生数
quantile_models = {q: GradientBoostingRegressor(loss='quantile', alpha=q,
                                                random_state=0).fit(X, y)
                   for q in (0.10, 0.90)}


# τ=0.9 モデルの平均 pinball loss
y_pred_90 = quantile_models[0.90].predict(X)
loss_90 = mean_pinball_loss(y, y_pred_90, alpha=0.90)
print(f'τ=0.90 の pinball loss: {loss_90:.2f}')

# カバレッジ（予測区間に入った観測値の割合）
y_lo = quantile_models[0.10].predict(X)
y_hi = quantile_models[0.90].predict(X)
coverage = ((y >= y_lo) & (y <= y_hi)).mean()
print(f'80%予測区間のカバレッジ: {coverage*100:.1f}%')  # 理論値 80%

5. 多変量への拡張（複数説明変数）

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20

# ── この抜粋で使うデータを用意します ──
import pandas as pd
import statsmodels.formula.api as smf

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
df = df[df['SSDSE-B-2026'] == 2023]          # 2023 年の 47 都道府県
df = df.copy()
df['宿泊万人泊'] = df['G7101'].astype(float) / 1e4   # 延べ宿泊者数
df['人口万人'] = df['A1101'].astype(float) / 1e4

# 人口 + 気温 + 消費支出（万円）で 3 つの分位点を推定
df['気温'] = df['B4101']
df['消費支出万円'] = df['L3221'] / 1e4

for tau in [0.10, 0.50, 0.90]:
    res = smf.quantreg('宿泊万人泊 ~ 人口万人 + 気温 + 消費支出万円',
                       data=df).fit(q=tau)
    print(f'τ={tau}:')
    print(res.params.round(4))
    print()

⚠️ 5 つの「分位点回帰 落とし穴」

🌐 関連手法・派生

🔗 関連用語（前提・並列・発展）

🟦 前提となる用語

🟨 並列・比較対象

🟩 発展・応用

📚 関連グループ教材

• 回帰分析の総合教材 — OLS、 ロジスティック、 分位点の比較
• 線形回帰入門 — 平均を予測する古典的回帰
• 損失関数 — pinball loss を含む損失関数群の整理
• 外れ値処理 — 分位点回帰を含む頑健統計の章

🧪 SSDSE-B 実データで「OLS と分位点回帰の差」を実感する

「分位点回帰は外れ値に強い」「不等分散データで威力を発揮する」と説明されても、 実感が湧きにくいですよね。 ここでは SSDSE-B の 47 都道府県データで 5 つの比較実験 を行い、 OLS との具体的な差を数値で確認します。

実験 1: 東京を含めるか除くかで OLS の係数はどう変わる？

SSDSE-B 2023 年の「延べ宿泊者数 ~ 人口」回帰で：

解釈：OLS の傾きは東京・大阪など人口も宿泊者数も大きい「梃の効いた点」に強く依存し、 それらを除くと 4.11 → 2.61 と大きく低下する。 一方 中央値回帰（τ=0.5）は 4.08 → 3.50 と変動が小さく、 「典型的な都道府県での関係」 をより安定して表す。

実験 2: 説明変数を増やしたとき τ ごとに係数はどう動く？

$y$ = 延べ宿泊者数（万人泊）、 $x_1$ = 人口（万人）、 $x_2$ = 年平均気温、 $x_3$ = 消費支出（万円）で重回帰：

解釈：

• 人口係数は τ で増加傾向（+2.17 → +5.56）→ 規模（＝観光需要の集中）が上位県ほど強い
• 気温・消費支出の係数は τ ごとに符号が変わり不安定 → これらは宿泊者数の分位を強く規定していないことを示唆
• 結果として 人口が上位分位を最も強く説明する主因 だと読み取れる

OLS 1 本だけでは、 これらの「分布のなかでの違い」は 絶対に 見えません。

実験 3: 残差プロットの読み方

OLS で残差プロット（x 軸：予測値、 y 軸：残差）を描いたとき、 残差の広がりが 「漏斗状」 なら不等分散の証拠。 これは分位点回帰の出番です。

1
2
3
4
5
6
7
8
9

# OLS 残差プロット
import matplotlib.pyplot as plt
ols_resid = ols.resid
fig, ax = plt.subplots(figsize=(8, 5))
ax.scatter(ols.fittedvalues, ols_resid, alpha=0.6, s=60)
ax.axhline(0, color='red', linestyle='--')
ax.set_xlabel('予測値'); ax.set_ylabel('残差')
ax.set_title('OLS 残差プロット — 漏斗状なら不等分散')
plt.show()

実験 4: Breusch-Pagan 検定で不等分散を統計的に確認

1
2
3
4
5

from statsmodels.stats.diagnostic import het_breuschpagan

bp = het_breuschpagan(ols.resid, ols.model.exog)
print(f'BP statistic: {bp[0]:.2f}, p-value: {bp[1]:.4f}')
# p < 0.05 なら不等分散の証拠 → 分位点回帰の出番

実験 5: ブートストラップで分位点係数の信頼区間

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12

import numpy as np

n_boot = 1000
boot_slopes = []
for _ in range(n_boot):
    idx = np.random.choice(len(df), len(df), replace=True)
    sample = df.iloc[idx]
    res = smf.quantreg('宿泊万人泊 ~ 人口万人', data=sample).fit(q=0.9)
    boot_slopes.append(res.params[1])

ci_lo, ci_hi = np.percentile(boot_slopes, [2.5, 97.5])
print(f'τ=0.9 傾きの 95% 信頼区間: [{ci_lo:.2f}, {ci_hi:.2f}]')

📊 応用シナリオ集 — 分位点回帰が真価を発揮する場面

シナリオ 1: 公共政策の格差分析

「教育予算が増えたら学力は上がるか？」を考えるとき、 OLS の平均効果だけ見ても 下位層の子どもに効いたのか、 上位層がさらに伸びたのか が分からない。 分位点回帰で τ=0.1（学力下位 10%）と τ=0.9（学力上位 10%）の係数を比較すると、 教育政策の 平等化効果か格差拡大効果か を識別できる。

SSDSE-B では 47 都道府県の「教育費」と「学力テスト得点」の関係を τ で分解し、 政策の 「底上げ効果」 vs 「トップ層伸長効果」を可視化できる。

シナリオ 2: 金融リスク管理 (VaR の推定)

Value-at-Risk (VaR) は「99% の確率でこれ以上は損しない」金額。 これは τ=0.01 の分位点 そのもの。 リターン分布の左裾を 説明変数（マクロ経済指標）の関数として 推定できる。

• $y$ = ポートフォリオの日次リターン
• $x$ = VIX、 金利スプレッド、 為替変動など
• $\tau = 0.01, 0.05$ の分位点回帰で「悪条件下での最悪リターン」を予測

シナリオ 3: 気象・気候予測の不確実性

「明日の最高気温は 23 ℃」と単一値で予測するより、 「80% の確率で 19〜26 ℃ に入る」と区間で予測する方が意思決定に有用。 分位点回帰で τ=0.1, 0.9 を同時推定すれば、 即座に 80% 予測区間が手に入る。

シナリオ 4: 医療における体重曲線

WHO の 成長曲線（growth chart） は本質的に分位点回帰の応用。 「年齢 $x$ の子どもの体重」の τ=3%, 50%, 97% を引いて、 「あなたの子は同年齢の何 % 位にいる」と判断する。

シナリオ 5: 不平等指標との接続

所得分布の τ=0.1 と τ=0.9 の差は 「10/90 比」（不平等指標）。 これを 地域・年齢・教育水準 の関数として推定すれば、 不平等が「誰のあいだで」「いつ」拡大したかを分解できる。 ピケティ的な分析。

📋 比較チートシート — どの回帰を選ぶ？

分位点回帰を「使うべきサイン」 — 5 つの兆候

1. 残差プロットで 「漏斗状」（広がりが x で変化）の不等分散が見える
2. OLS の残差が 非対称な分布（歪んでいる）
3. 「平均」よりも「下位／上位の動き」 が政策・実務的に重要
4. 「予測値だけでなく 予測区間」を出す必要がある
5. 少数の外れ値が OLS 結果を大きく動かしている

🎓 理論的補足 — 推定量の性質と歴史

歴史的経緯

• 1755 年：Roger Boscovich が天文観測データに 絶対誤差最小化（中央値回帰の祖）を提案。 同年代に Laplace も独立に研究。
• 1789 年：Laplace が「絶対誤差和最小化」の解法を整備。 これが現代の τ=0.5 分位点回帰の原型。
• 1805 年：Legendre が 最小二乗法（OLS）を発表。 計算容易性から L1 系を圧倒し、 19〜20 世紀の主流に。
• 1978 年：Roger Koenker と Gilbert Bassett が論文 “Regression Quantiles”（Econometrica）を発表。 任意 τ の分位点を線形計画問題として推定する現代的フレームワークを確立。
• 1990s：計算機の発達で実用化。 経済学・公衆衛生・気象で急速に普及。
• 2005〜2015：分位点ランダムフォレスト（Meinshausen 2006）、 ニューラル分位点回帰、 LightGBM/XGBoost の objective 統合などで 機械学習版 が確立。

分位点回帰推定量の漸近的性質

適切な正則条件下で、 分位点回帰推定量 $\hat{\boldsymbol\beta}(\tau)$ は次の性質を持つ：

• 一致性：$\hat{\boldsymbol\beta}(\tau) \overset{p}{\to} \boldsymbol\beta(\tau)$
• 漸近正規性：$\sqrt{n}\bigl(\hat{\boldsymbol\beta}(\tau) - \boldsymbol\beta(\tau)\bigr) \overset{d}{\to} \mathcal{N}(0, \Sigma(\tau))$
• 分散共分散行列 $\Sigma(\tau)$ は 誤差の条件付き密度 $f_{u|x}(0)$ に依存。 これが推定誤差の計算を難しくする。

OLS との効率比較

誤差が正規分布のとき、 中央値回帰の漸近相対効率（OLS 比）は：

線形計画問題としての定式化

pinball loss の最小化は、 補助変数 $u_i^+, u_i^- \ge 0$ を導入することで 線形計画問題（LP） として書ける：

これによりシンプレックス法・内点法など 確立された LP ソルバー がそのまま使え、 大規模問題でも実用的に解ける。

🔢 statsmodels.QuantReg / LightGBM 主要パラメータ表

statsmodels.regression.quantile_regression.QuantReg.fit() の主要引数

LightGBM の objective='quantile' 関連パラメータ

📚 ミニケーススタディ — SSDSE-B 死亡率の分位点分析

もう 1 つ実データで使い方を確認しましょう。 SSDSE-B 2026 の 「死亡率 ~ 高齢化率」 の関係を、 通常の回帰と分位点回帰で比較します。

背景

「高齢化率が上がると死亡率も上がる」は当然の関係ですが、 「どの程度の傾きで」「下位／上位の県では別の傾向か」を分位点回帰で詳細に見ます。

OLS の結果（参考）

$y$ = 死亡率（人口千あたり）、 $x$ = 高齢化率（%）。 47 県で OLS：

分位点回帰の結果

解釈

• 傾きが τ で 緩やかに減少（0.646 → 0.583、 約 10% 減）。 高齢化に対する死亡率の感応度は 低死亡率群でやや強い が、 全体としては均一（等分散に近い）。
• OLS と τ=0.5 の傾きはほぼ一致（0.610 vs 0.624）。 これは 外れ値・歪みが小さく、 線形関係が極めて頑健 なことを示す。
• これは 分位点回帰を使ってもあまり情報が増えない ケース。 逆にこういう状況なら OLS を安心して使ってよいと判定できる。

逆に分位点回帰が威力を発揮するケース

SSDSE-B の 「延べ宿泊者数 ~ 人口」 では、 τ=0.10 と τ=0.90 で傾きが 2.13 vs 5.57 と約 2.6 倍も違う。 これは 明確な不等分散 の証拠で、 OLS の 1 本の傾きでは分布の広がりの情報を失う。 こういうデータには分位点回帰が必須です。

使い分けのまとめ

❓ よくある質問 (FAQ)

Q1. τ=0.5 と OLS、 結局どっちを使うべき？

A1. データに 明らかな外れ値や歪み があるなら τ=0.5（中央値回帰）の方が頑健。 誤差が正規・等分散で外れ値もない理想的状況なら OLS の方が 1.57 倍効率がよい（漸近相対効率）。 実務では 両方推定して比較 するのが正攻法。 大きく違うなら外れ値や不等分散を疑う合図。

Q2. 何個 τ を推定すべき？

A2. 探索目的なら τ ∈ {0.05, 0.10, 0.25, 0.50, 0.75, 0.90, 0.95} の 7 点が定番。 予測区間が目的なら 90% 区間 → τ=0.05, 0.95、 80% 区間 → τ=0.1, 0.9。 多すぎると分位点交差問題が増える。

Q3. p 値や信頼区間はどう計算？

A3. statsmodels の summary() はカーネル密度を使った漸近近似で SE を返すが、 小標本では信頼性が低い。 ブートストラップ法（1000 リサンプル）で信頼区間を求めるのが推奨。 上記コード例参照。

Q4. 多変量・カテゴリ変数も扱える？

A4. 完全に扱える。 statsmodels の formula API なら y ~ x1 + x2 + C(category) のように OLS と同じ書き方でカテゴリダミー化も自動。 LightGBM 版は数値・カテゴリ・欠損値もそのまま処理可能。

Q5. 時系列データでも使える？

A5. 使える。 ただし誤差の独立性が崩れるので、 信頼区間にブロックブートストラップを使うこと。 ボラティリティ予測やリスク管理で CAViaR（Conditional Autoregressive VaR、 Engle & Manganelli 2004）という拡張モデルが有名。

Q6. ニューラルネットでもできる？

A6. できる。 損失関数を pinball loss に変えるだけ。 PyTorch なら：

1
2
3
4
5
6
7
8
9

import torch

def pinball_loss(y_pred, y_true, tau):
    err = y_true - y_pred
    return torch.max(tau * err, (tau - 1) * err).mean()

# 学習ループの中で：
#   loss = pinball_loss(model(x), y, tau=0.9)
#   loss.backward(); optimizer.step()

Q7. 解釈は OLS と同じでいい？

A7. ほぼ同じ。 「$x_1$ が 1 単位増えると、 $y$ の τ 分位点 が $\beta_1(\tau)$ だけ増える」と読む。 OLS の「平均が増える」を「τ 分位点が増える」に置き換えるだけ。

Q8. 計算コストは？

A8. 線形分位点回帰は LP として解けるので、 n=10,000、 p=50 程度なら 1 秒以内。 OLS の正規方程式解法に比べれば遅いが、 LightGBM 版なら GBDT と同等。

📚 補足: 分位点回帰の深掘り解説と実務展開

✅ 理解度チェック

以下の問いに自分の言葉で答えられるか確認しよう。 即答できない項目は本文の対応セクションへ戻って復習することを推奨する。

1. 定義: τ=0.5 の分位点回帰と OLS の違いを、 損失関数のレベルで一文で説明できるか。
2. 計算: τ=0.25, 0.5, 0.75 のピンボール損失を 3 つのデータ点で手計算できるか。
3. 使い分け: 外れ値が多い SSDSE-B-2026 の延べ宿泊者数の分布に OLS と分位点回帰のどちらを使うべきか、 理由とともに答えられるか。
4. 解釈: τ=0.9 の回帰係数が正、 τ=0.1 が負だった場合、 何が分かるかを 30 秒で要約できるか。
5. 実装: statsmodels.QuantReg と sklearn.linear_model.QuantileRegressor のどちらをいつ使うか、 違いを 3 行で述べられるか。
6. 応用: 予測区間として τ=0.05 と τ=0.95 を使うとき、 信頼水準は何 % か、 また OLS の予測区間とどう違うか。

💡 ヒント: 5 問以上即答できれば実務レベル。 3 問未満なら 30 秒で分かる結論と Python 実装セクションを優先的に復習。

理論背景 — Koenker と Bassett の発見

分位点回帰 (quantile regression) は 1978 年に Koenker と Bassett が発表した論文 Regression Quantiles によって体系化された。 それ以前にも中央値回帰の発想は Edgeworth (1888) や Boscovich (18 世紀) にまで遡れるが、 これらは LAD (Least Absolute Deviation) 回帰として知られていた。 Koenker と Bassett の革新は 「任意の分位点 τ ∈ (0,1) を一つの統一フレームワークで扱える」 という発見にあった。 τ=0.5 で LAD 回帰が、 τ→0 や τ→1 で極値に近い分位点が、 すべて同じ ピンボール損失 (check loss) $\rho_\tau(u) = u(\tau - \mathbb{1}_{u<0})$ の最小化として解ける。 この発見によって、 「条件付き分布の全体像」 を回帰で記述する道が開けた。 OLS は条件付き期待値しか教えてくれないが、 分位点回帰は分布の任意の点で「説明変数が結果にどう効くか」 を答えてくれる。 SSDSE-B-2026 の延べ宿泊者数のように分布が右に長く裾を引くデータでは、 平均と中央値が大きく乖離するので、 分位点回帰は特に威力を発揮する。

ピンボール損失の数理 — なぜ非対称損失なのか

ピンボール損失の正体は 「下振れと上振れに非対称な重みを置く絶対値損失」 である。 数式で書くと、 残差 $u = y - \hat{y}$ について τ=0.5 なら $\rho_{0.5}(u) = 0.5|u|$ となり、 これは LAD と同じ。 τ=0.9 なら $\rho_{0.9}(u) = 0.9 u \text{ if } u > 0$、 $\rho_{0.9}(u) = -0.1 u \text{ if } u < 0$ となり、 下振れ (予測が大きすぎた) のコストが小さく、 上振れ (予測が小さすぎた) のコストが大きい。 結果として、 推定値は「下振れより上振れを避ける方向」 に引き寄せられ、 0.9 分位点 (上位 10%) が得られる。 この非対称性は、 OLS の二乗損失や LAD の対称絶対値損失とは本質的に異なる発想で、 「リスク管理ではどちら側の誤差をより避けたいか」 という実務的問いを直接表現できる点が魅力である。 SSDSE-B-2026 の県別失業率を予測する際、 「失業率を過小推定する方が政策的に危険」 と考えるなら、 上振れペナルティを大きくした分位点回帰が自然な選択となる。

分位点回帰の頑健性 — 外れ値が結論を変えない

分位点回帰の最大の魅力は 外れ値に対する頑健性 である。 OLS は二乗誤差を最小化するため、 大きな残差を持つ 1 点が回帰係数に過剰な影響を与える (影響関数が無限大)。 一方、 分位点回帰のピンボール損失は線形なので、 影響関数は有界である。 これは Huber や Hampel が研究した ロバスト統計 の伝統に連なる発想で、 SSDSE-B-2026 のように東京・大阪のような巨大都市が外れ値として混じるデータでは特に重要である。 たとえば「県別の延べ宿泊者数」 を OLS で人口に回帰すると、 東京の影響で係数が大きく引っ張られる。 同じデータを τ=0.5 の分位点回帰にかけると、 東京の影響は限定的で、 「典型的な県」 のメカニズムが浮かび上がる。 さらに τ=0.1 と τ=0.9 を比較すれば、 「宿泊需要の小さい県と大きい県でメカニズムがどう違うか」 までわかる。 これは OLS では絶対に得られない構造的洞察である。

SSDSE-B-2026 での実例 — 県別延べ宿泊者数の分位点回帰

発展 — 加法分位点回帰、 ベイズ分位点回帰、 量子勾配ブースティング

分位点回帰は線形に留まらず、 様々な拡張が研究されている。 (1) 加法分位点回帰 (additive quantile regression): スプラインを使って非線形な分位点曲線を推定。 GAM (一般化加法モデル) の分位点版で、 qgam パッケージ (R) が定番。 (2) ベイズ分位点回帰 (Bayesian quantile regression): 非対称ラプラス分布を尤度に使うことで、 ピンボール損失と等価な事後分布を構成できる。 階層構造を入れて都道府県別に異なる分位点回帰を推定する応用が可能。 (3) 量子勾配ブースティング (quantile gradient boosting): LightGBM や XGBoost で objective='quantile' を指定すると、 非線形な分位点回帰を自動的に学習できる。 SSDSE-B-2026 の都道府県データのように小サンプルでも、 量子勾配ブースティングは交差検証で過学習を制御しながら頑健な予測区間を提供してくれる。 (4) 分位点ランダムフォレスト (Meinshausen 2006): RF の葉ノードに落ちた訓練点の分位点を返すだけで、 任意の τ の予測ができる。 これらの拡張は分位点回帰の応用範囲を、 線形仮定の枠を超えて大きく広げる。

落とし穴と回避策

分位点回帰には特有の落とし穴がある。 第一に 分位点クロッシング: 複数の τ で別々に回帰すると、 たとえば τ=0.4 の予測値が τ=0.5 を上回るといった矛盾が生じうる。 これは独立に推定するため避けられず、 同時推定 (joint estimation) や事後のソートで修正する。 第二に 計算コスト: ピンボール損失は微分不可能なため線形計画 (LP) で解く。 大規模データでは内点法や勾配ベース近似 (smoothed quantile regression) を使う。 第三に 標準誤差の計算: 漸近分散の解析式は条件付き密度関数を含むため、 ブートストラップで推定するのが安全。 第四に τ=0 や τ=1 に近い極値分位点: 推定が極めて不安定なので、 τ=0.05 程度を下限に留めるのが実用的。 SSDSE-B-2026 の 47 サンプルというサイズでは、 τ=0.1, 0.5, 0.9 程度を主軸にすると安定した結論が得られる。

まとめ — 平均だけでは見えない世界

分位点回帰は OLS の延長線上にあるが、 哲学が根本的に異なる。 OLS は「平均」 を見る道具、 分位点回帰は「分布全体」 を見る道具である。 SSDSE-B-2026 のような社会経済データでは、 平均は政策の参考にはなるが、 格差や弱者の状況は中央値や下位分位点でしか見えない。 「上位 10% と下位 10% で説明変数の効果がどう違うか」 という問いに答えられるのが分位点回帰だけである。 統計データ分析コンペティションのレポートで、 OLS 結果だけでなく分位点回帰の結果を併記すれば、 「平均的にはこうだが、 弱者層ではこう、 富裕層ではこう」 という多層的な物語が描ける。 これは政策提言の説得力を飛躍的に高める。 分位点回帰は単なる手法というよりも、 「分布全体を真摯に見る」 という分析者の姿勢の宣言である。 この姿勢を身につけた学生は、 平均という単一指標に縛られず、 データが語る豊かな多様性を捉えられる分析者へと成長する。

解釈の落とし穴 — 「分位点回帰係数」 が何を表すか

分位点回帰の係数解釈は OLS よりやや注意が必要である。 τ=0.5 の係数 $\beta_{0.5}$ は「説明変数 $X$ が 1 単位増えたとき、 条件付き中央値が $\beta_{0.5}$ だけ変化する」 と読む。 ここで重要なのは 「条件付き分布の中央値」 であって、 「中央値の人」 の所得変化ではない点である。 SSDSE-B-2026 で 「人口を 1 万人増やすと中央値の延べ宿泊者数が 4.08 万人泊増える」 と言うとき、 これは「中央値だった県」 が 4.08 万人泊増えるのではなく、 「もし全県の人口が 1 万人増えれば、 県別 宿泊者数分布の中央値が 4.08 万人泊シフトする」 という意味になる。 この区別を曖昧にすると、 「中央値の県を 4.08 万人泊押し上げる政策」 という誤った提言が生まれる。 分位点回帰の正しい解釈は「集団全体の分布形状の変化」 であり、 個別の県の動きを保証するものではない。 SSDSE-B-2026 のような横断データでは、 政策含意を語る際にこの集団 vs 個人の混同 (ecological fallacy の一種) に注意を払うことが、 質の高いレポート執筆の決め手となる。

因果推論との接続 — 条件付き分位点と政策評価

近年、 分位点回帰は 因果推論 の文脈でも重要性を増している。 平均処置効果 (Average Treatment Effect, ATE) は OLS や IV で測れるが、 「処置が分布のどこに効くか」 を知りたい場合は 分位点処置効果 (Quantile Treatment Effect, QTE) が必要になる。 たとえば「最低賃金引き上げが所得分布のどこに最も効くか」 を問う研究では、 OLS だと「平均所得は上がった」 という結論しか出ない。 QTE を使えば「下位 25 パーセンタイルでは大きく上昇、 中央値ではほぼ変化なし、 上位では微減」 といった分布全体の動きが見える。 SSDSE-B-2026 を用いて 「観光振興政策が県別経済格差をどう変えたか」 を分析する場合、 単に平均の変化を見るだけでなく、 下位分位点が伸びたか上位分位点だけが伸びたかを QTE で確認すれば、 格差是正効果が定量化できる。 Chernozhukov-Hansen の IV-QR や Firpo の Unconditional QTE といった手法は、 分位点回帰の因果推論版として急速に普及している。

金融・リスク管理での応用 — VaR と CVaR の橋渡し

分位点回帰は金融リスク管理の基幹技術でもある。 銀行や投資ファンドが日常的に計算する Value-at-Risk (VaR) は、 損失分布の極値分位点 (τ=0.95 や τ=0.99) そのものである。 ポートフォリオの収益分布を市場指標 (株価指数、 為替、 金利) に回帰するとき、 OLS は平均的なリスクしか教えてくれない。 一方、 分位点回帰なら「最悪のケースで損失がどう動くか」 が直接モデル化できる。 さらに Expected Shortfall (ES, CVaR) は VaR を超える部分の条件付き期待値で、 ピンボール損失とは異なる ES 損失 (Fissler-Ziegel score) を用いて推定する。 SSDSE-B-2026 の金融指標ではなく、 県別経済データに置き換えても同じ発想が活きる。 「景気後退期に最も大きな損失を被る県はどこか」 を回帰モデルで表現すると、 政策立案者にとって極めて有用な情報になる。 分位点回帰は金融工学から地域経済学まで、 「リスクの非対称性」 が問題になるあらゆる領域で使われている。

学習者へのメッセージ

最後に学習者に伝えたいのは、 「分位点回帰は OLS の代替ではなく補完である」 という点だ。 OLS は平均を見るための最も鋭い道具であり続けるが、 平均だけでは見落とすことが多すぎる。 SSDSE-B-2026 の県別 延べ宿泊者数、 出生数、 消費支出、 気温といったさまざまな指標は、 平均値以上に分布形状そのものが政策的意味を持つ。 「平均値だけで議論を済ませている人」 と 「分位点を語れる人」 では、 同じデータから引き出せる洞察の深さが格段に違う。 統計データ分析コンペティションのレポート執筆では、 OLS と分位点回帰を併記して 「平均ではこう、 下位ではこう、 上位ではこう」 と語れる立体的な分析を目指してほしい。 そうすれば、 単純な相関係数や回帰係数では語れない、 「データの背後にある社会の多様性」 を表現できる。 これは分析者として一段上のレベルに到達する道筋であり、 SSDSE-B-2026 という教材はそのトレーニングに最適である。 分位点回帰は今後あなたの分析人生で何度も助けになってくれるはずだ。

補足として、 分位点回帰の学習リソースを挙げておく。 教科書としては Koenker の Quantile Regression (Cambridge, 2005) が決定版で、 理論から応用まで網羅されている。 日本語文献では小西貞則・北川源四郎『情報量規準』 や金本拓『計量経済学』 が分位点回帰を扱う。 Python での実装は statsmodels.QuantReg、 sklearn.linear_model.QuantileRegressor が標準。 R では quantreg (Koenker 自身がメンテナンス) と qgam (加法分位点回帰) が定番である。 SSDSE-B-2026 を題材にした演習としては、 「県別 延べ宿泊者数を人口・気温・消費支出に回帰」、 「県別 出生数を人口・婚姻件数に回帰」 などが教育効果が高い。 ぜひ自分の興味に近いテーマで τ=0.1, 0.5, 0.9 の比較を試し、 OLS との違いを体感してほしい。 そうすれば、 「平均だけ見ていた頃の自分」 と 「分位点を語れる自分」 の間に決定的な飛躍が起きていることに気付くだろう。 この飛躍こそが、 統計教育の真の成果である。

閉じくくりとして、 分位点回帰の発想は機械学習の損失関数設計にも深く影響を与えている。 LightGBM や Neural Network で予測区間を学ぶ際、 ピンボール損失を使うだけで複数の τ を同時に学習できる。 「Conformal Prediction」 と組み合わせれば、 任意のブラックボックスモデルから理論的に正しい予測区間を構成できる。 SSDSE-B-2026 を題材にした研究では、 こうした最新技術を取り入れた予測区間の評価が、 単純な点推定よりも価値の高い貢献となる。 分位点回帰を起点に、 機械学習・統計推論・因果推論の交差点を探求する旅は、 データサイエンス学習における最も魅力的な旅路の一つである。 ぜひこの旅を、 SSDSE-B-2026 の身近なデータから始めてほしい。 自分の出身県や関心のある地域のデータを τ=0.1, 0.5, 0.9 で回帰してみるだけで、 「平均では見えない物語」 が次々と立ち現れる。 その発見の体験こそが、 統計を学ぶ最大の喜びである。 分位点回帰はその喜びへの入口となる強力な道具であり、 一度習得すればあなたの分析の幅を一気に広げてくれるだろう。 統計データ分析コンペティション 2026 に挑む受講生にとって、 OLS と分位点回帰の併用は差別化要因になりうる。 平均だけを語るレポートではなく、 分布全体を語るレポートを目指せば、 審査員に強い印象を残せる。 分布の声に耳を傾ける分析者であれ。 それが、 分位点回帰が君に教える最大のレッスンである。 平均値だけを見てきた今までの自分から、 分布全体に思考を広げる新しい自分へと脱皮するための鍵が、 ここに、 ピンボール損失というシンプルな式の中に確かに眠っている。 その鍵を取り出して使い始めれば、 君のデータの見方は永遠に変わる。 これは大げさではなく、 多くの分析者が体験してきた事実である。 さあ、 自分のデータで早速試してみよう。 きっと、 今まで気付かなかった世界が広がっているはずだ。 これが分位点回帰の真の力であり、 学習する価値の核心である。 ぜひ自分の研究テーマで実際に活用してみてほしい。

🗺 概念マップ — 分位点回帰の位置づけ

🔗 隣接手法への橋渡し

「分位点回帰」は単独で完結する手法ではなく、 隣接領域と連携することで真価を発揮する。

• 上流: 分位点・回帰分析の基礎の上に成立
• 並列: 最小二乗回帰 (条件付き平均) と比較し、 中央値や上下分位を直接モデル化
• 下流: ロバスト統計・箱ひげ図で結果を可視化

SSDSE-B-2026 の延べ宿泊者数を分位点回帰で予測すると、 中央値 (τ=0.5) と上位 (τ=0.9) で人口の係数が異なり、 「人口が増えると上位県ほど宿泊者数が急増」の不均一性が定量化できる。

🌳 手法選択フロー

「分位点回帰」を実際の課題に当てはめるとき、 状況別に何を選ぶかを 3 段階で判定する。

1. 分布の形は対称か? 対称 → OLS 回帰で十分、 非対称 (人口・所得) → 分位点回帰
2. 関心は中央値か上下端か? 中央値 → 0.5 分位、 上位 → 0.9、 下位 → 0.1 を指定
3. 外れ値はあるか? 多 → 分位点回帰の頑健性が活きる、 少 → OLS と結果が近い

SSDSE-B-2026 の延べ宿泊者数を人口で回帰する際、 OLS は東京の影響で傾きが過大評価される。 0.5 分位回帰なら中央値ベースの「典型的な県の人口-宿泊者数 関係」が得られる。