【ベンチマーク】分散不均一な回帰

最終更新 2026-03-06 読了時間 6 分

まとめ

このベンチマークでは、\(y = 3x + \epsilon\)、\(\epsilon \sim N(0, (0.5x)^2)\) と \(x\) が大きいほどノイズが広がる 500 点のデータで、OLS・WLS・Ridge・HuberRegressor・QuantileRegressor の 5 モデルを比較します。前処理として「なし」「log(y)」「√y」の 3 通りを掛け合わせ、5 分割交差検証の MAE と RMSE をヒートマップで並べます。OLS の傾きの推定は不偏のまま 3 に近づきますが、\(x\) の大きい区間で残差が扇形に広がり、予測区間の幅を一律に見積もると実態と合わなくなります。WLS に正しい重み \(1/x^2\) を与えるか log 変換で分散を揃えると、残差の広がりが均一になります。

売上が大きい店舗ほど売上のブレも大きい、所得が高いほど消費のばらつきも大きい——実務データでは分散が一定でないケースが多い。等分散を仮定する OLS をそのまま適用すると、係数推定は正しくても信頼区間や予測区間が歪む。ここでは分散不均一（heteroscedasticity）な合成データで最適なアプローチを検証する。

1. データの特徴 #

説明変数は 1 つで、\(x\) は \([1, 20]\) の一様分布から 500 点を引きます。真のモデルは \(y = 3x + \epsilon\) で、ノイズ \(\epsilon \sim N(0, (0.5x)^2)\) の標準偏差が \(0.5x\) と \(x\) に比例して増えます。\(x = 1\) では標準偏差 0.5、\(x = 20\) では 10 なので、散布図は右に行くほど扇形に広がります。

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

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

見るポイント：RMSE は二乗するぶん 1 つの大外れに引きずられる。MAE は線形、中央絶対誤差はほぼ動かない。

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
30
31
32
33
34
35
36
37
38
39
40
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

rng = np.random.default_rng(42)
n = 500

x = rng.uniform(1, 20, size=n)
noise = rng.normal(0, 0.5 * x)
y = 3 * x + noise

fig, axes = plt.subplots(1, 2, figsize=(13, 4))

axes[0].scatter(x, y, alpha=0.4, s=15, color="#2563eb")
axes[0].plot([1, 20], [3, 60], color="#ef4444", linewidth=2, label="真のモデル y=3x")
axes[0].set_title("分散不均一データ（扇形に広がる）")
axes[0].set_xlabel("x")
axes[0].set_ylabel("y")
axes[0].legend()
axes[0].grid(alpha=0.2)

# 残差の分散をxの範囲ごとに可視化
bins = np.linspace(1, 20, 6)
bin_labels = [f"{bins[i]:.0f}-{bins[i+1]:.0f}" for i in range(len(bins)-1)]
residuals = y - 3 * x
bin_idx = np.digitize(x, bins) - 1
bin_idx = np.clip(bin_idx, 0, len(bin_labels) - 1)
bp_data = [residuals[bin_idx == i] for i in range(len(bin_labels))]
bp = axes[1].boxplot(bp_data, labels=bin_labels, patch_artist=True)
for patch in bp["boxes"]:
    patch.set_facecolor("#2563eb")
    patch.set_alpha(0.5)
axes[1].set_title("xの範囲ごとの残差の広がり")
axes[1].set_xlabel("xの範囲")
axes[1].set_ylabel("残差 (y - 3x)")
axes[1].grid(alpha=0.2)

fig.suptitle("不等分散データの構造", fontsize=13)
fig.tight_layout()
plt.show()

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

8-10 行目ノイズの標準偏差を 0.5x にする
rng.normal(0, 0.5 * x) は scale に配列を渡しているので、点ごとに違う幅のノイズが出ます。\(x = 1\) の点は標準偏差 0.5、\(x = 20\) の点は 10 で、これが不等分散そのものです。0.5 * x を定数 5.0 に替えると等分散のデータになり、以降の 5 モデルの差はほぼ消えます。
14-15 行目真の直線 y=3x を基準線として描く
[1, 20] と [3, 60] を結ぶ赤い線が真のモデルで、点がその周りにどう散るかを見るための基準です。左の図で右側ほど点が赤い線から離れて広がる「扇形」が、この後のすべての比較の出発点になります。
23-29 行目x を 5 区間に切って残差の箱ひげ図にする
np.linspace(1, 20, 6) で \(x\) を 5 つの等幅区間に分け、np.digitize で各点の区間番号を求め、真の直線からの残差 y - 3*x を区間ごとに箱ひげ図にしています。右の図で箱の高さが区間ごとに 1〜4 倍と伸びていくのが、分散が \(x\) に比例して増えている数値的な証拠です。

不等分散データの構造

3. 比較するパイプライン #

前処理	説明	モデル
なし（生データ）	そのまま回帰	5モデル
log(y) 変換	対数変換で分散安定化	5モデル
√y 変換	平方根変換で分散安定化	5モデル

モデル	特徴
OLS	最小二乗法（等分散仮定）
WLS	重み付き最小二乗法（weights=1/x²）
Ridge	L2正則化
Huber	ロバスト回帰（外れ値に強い）
Quantile(0.5)	中央値回帰（分布の歪みに強い）

4. 実験と結果 #

 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
95
96
97
98
99
from sklearn.model_selection import cross_val_score, KFold
from sklearn.linear_model import LinearRegression, Ridge, HuberRegressor, QuantileRegressor
from sklearn.preprocessing import FunctionTransformer
from sklearn.pipeline import Pipeline
from sklearn.base import clone
import seaborn as sns

cv = KFold(n_splits=5, shuffle=True, random_state=42)
X = x.reshape(-1, 1)

# 変換関数
def log_transform(y):
    return np.log1p(y)

def sqrt_transform(y):
    return np.sqrt(np.maximum(y, 0))

def inv_log(y):
    return np.expm1(y)

def inv_sqrt(y):
    return y ** 2

models = {
    "OLS": LinearRegression(),
    "WLS": LinearRegression(),  # 手動で重み適用
    "Ridge": Ridge(alpha=1.0),
    "Huber": HuberRegressor(epsilon=1.35),
    "Quantile(0.5)": QuantileRegressor(quantile=0.5, alpha=0.0, solver="highs"),
}

transforms = {
    "なし": (None, None),
    "log(y)": (log_transform, inv_log),
    "√y": (sqrt_transform, inv_sqrt),
}

results = []
for trans_name, (fwd, inv) in transforms.items():
    for mdl_name, mdl in models.items():
        maes, rmses = [], []
        for train_idx, test_idx in cv.split(X):
            X_tr, X_te = X[train_idx], X[test_idx]
            y_tr, y_te = y[train_idx], y[test_idx]

            y_fit = fwd(y_tr) if fwd else y_tr

            m = clone(mdl)
            if mdl_name == "WLS" and fwd is None:
                # 重みを 1/x² に設定
                w = 1.0 / (X_tr.ravel() ** 2 + 1e-6)
                m.fit(X_tr, y_fit, sample_weight=w)
            else:
                try:
                    m.fit(X_tr, y_fit)
                except Exception:
                    maes.append(np.nan)
                    rmses.append(np.nan)
                    continue

            pred = m.predict(X_te)
            if inv:
                pred = inv(pred)

            maes.append(np.mean(np.abs(y_te - pred)))
            rmses.append(np.sqrt(np.mean((y_te - pred) ** 2)))

        results.append({
            "前処理": trans_name, "モデル": mdl_name,
            "MAE": np.nanmean(maes), "RMSE": np.nanmean(rmses),
        })

df_res = pd.DataFrame(results)
pivot = df_res.pivot_table(index="前処理", columns="モデル", values="MAE")
trans_order = ["なし", "log(y)", "√y"]
mdl_order = ["OLS", "WLS", "Ridge", "Huber", "Quantile(0.5)"]
pivot = pivot.reindex(index=[t for t in trans_order if t in pivot.index],
                      columns=[m for m in mdl_order if m in pivot.columns])

fig, axes = plt.subplots(1, 2, figsize=(14, 4))

sns.heatmap(pivot, annot=True, fmt=".2f", cmap="RdYlGn_r",
            linewidths=0.5, ax=axes[0], cbar_kws={"label": "MAE（低いほど良い）"})
axes[0].set_title("前処理 × モデル: MAE ヒートマップ")
axes[0].set_xlabel("")
axes[0].set_ylabel("")

pivot_rmse = df_res.pivot_table(index="前処理", columns="モデル", values="RMSE")
pivot_rmse = pivot_rmse.reindex(index=[t for t in trans_order if t in pivot_rmse.index],
                                columns=[m for m in mdl_order if m in pivot_rmse.columns])
sns.heatmap(pivot_rmse, annot=True, fmt=".2f", cmap="RdYlGn_r",
            linewidths=0.5, ax=axes[1], cbar_kws={"label": "RMSE（低いほど良い）"})
axes[1].set_title("前処理 × モデル: RMSE ヒートマップ")
axes[1].set_xlabel("")
axes[1].set_ylabel("")

fig.suptitle("分散不均一データ: 前処理 × モデル 精度比較", fontsize=13)
fig.tight_layout()
plt.show()

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

12-22 行目目的変数の変換と逆変換を組で用意する
log1p と expm1、sqrt と ** 2 を対にして、学習は変換後の \(y\) で行い、評価は元の尺度に戻してから MAE・RMSE を測ります。逆変換を忘れると log 尺度の誤差と生の誤差を比べることになり、ヒートマップの行どうしが比較できません。np.maximum(y, 0) は √y がノイズで負になった \(y\) で失敗しないための保護です。
24-30 行目5 モデルの設定を 1 か所にまとめる
WLS は LinearRegression そのもので、重みは 49〜52 行目で sample_weight として渡します。HuberRegressor(epsilon=1.35) は既定のしきい値、QuantileRegressor(quantile=0.5, alpha=0.0) は正則化なしの中央値回帰です。Ridge(alpha=1.0) は説明変数が 1 本なので効果は小さく、「正則化は不等分散の対策にならない」ことを示す比較用の 1 行です。
49-52 行目WLS の重みは 1/x²、生データのときだけ
ノイズの分散が \((0.5x)^2\) なので、分散の逆数に比例する 1.0 / X_tr**2 が理論上の最適な重みです。+ 1e-6 はゼロ除算の保護です。fwd is None の条件で、log や √ で変換した \(y\) に対しては重みを付けていません。変換後は分散がすでに均一に近く、重みを掛けると過剰補正になるためです。
54-59 行目学習に失敗した組は NaN で埋める
try / except で fit の失敗を捕まえ、その分割の MAE と RMSE を np.nan にしています。15 組 × 5 分割のループが 1 つの例外で止まらないための保険で、70 行目の np.nanmean が NaN を除いて平均します。ヒートマップに空白のセルがあれば、この経路を通ったと分かります。
73-78 行目前処理 × モデルの表に並べ替える
pivot_table で縦に前処理 3 行、横にモデル 5 列の表を作り、reindex で表示順を「なし → log → √」「OLS → WLS → Ridge → Huber → Quantile」に固定しています。順番を固定しないと pandas が文字コード順に並べ、2 枚のヒートマップで行列の順序が食い違うことがあります。
82-92 行目MAE と RMSE を同じ色尺度で読む
cmap="RdYlGn_r" は小さい値を緑、大きい値を赤にする色で、annot=True で各セルに数値を書き込みます。左が MAE、右が RMSE で、RMSE は残差を 2 乗するため \(x\) の大きい区間の大外れに敏感です。2 枚で順位が入れ替わるセルがあれば、そのモデルは「大外れは少ないが小さな誤差が多い」か、その逆だと読めます。

精度比較ヒートマップ

読み方のポイント #

OLS の傾きの推定は 3 に近く正しいのですが、\(x\) が大きい区間の残差がそのまま MAE と RMSE に乗るため、生データの行では他のモデルと大差がつきません。WLS に正しい重み \(1/x^2\) を与えると、ノイズの小さい左側の点で直線が決まり、傾きの推定の分散が小さくなります。ただし MAE・RMSE は元の尺度で測るので、右側の大きなノイズ自体は消えず、改善幅は限定的です。log(y) の行では変換後の分散がほぼ均一になるため、5 モデルすべてで RMSE が揃う方向に動きます。

5. 誤差パターン分析 #

 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
from sklearn.model_selection import train_test_split
from sklearn.base import clone

X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.3, random_state=42)

configs = [
    ("OLS（生データ）", LinearRegression(), None, None, None),
    ("WLS（1/x²重み）", LinearRegression(), None, None, "wls"),
    ("OLS + log(y)", LinearRegression(), log_transform, inv_log, None),
]

fig, axes = plt.subplots(1, 3, figsize=(15, 4))
for i, (name, mdl, fwd, inv, special) in enumerate(configs):
    m = clone(mdl)
    y_fit = fwd(y_tr) if fwd else y_tr

    if special == "wls":
        w = 1.0 / (X_tr.ravel() ** 2 + 1e-6)
        m.fit(X_tr, y_fit, sample_weight=w)
    else:
        m.fit(X_tr, y_fit)

    pred = m.predict(X_te)
    if inv:
        pred = inv(pred)

    residuals = y_te - pred
    axes[i].scatter(X_te, residuals, alpha=0.4, s=15, color="#2563eb")
    axes[i].axhline(0, color="#ef4444", linewidth=1, linestyle="--")
    axes[i].set_title(f"{name}", fontsize=10)
    axes[i].set_xlabel("x")
    axes[i].set_ylabel("残差")
    axes[i].grid(alpha=0.2)

    # 残差の標準偏差をx範囲ごとに表示
    for xmin, xmax in [(1, 7), (7, 14), (14, 20)]:
        mask = (X_te.ravel() >= xmin) & (X_te.ravel() < xmax)
        if mask.sum() > 0:
            std = residuals[mask].std()
            axes[i].text((xmin + xmax) / 2, axes[i].get_ylim()[1] * 0.85,
                        f"σ={std:.1f}", ha="center", fontsize=8,
                        bbox=dict(boxstyle="round", facecolor="#f0f0f0", alpha=0.8))

fig.suptitle("残差パターンの比較（x範囲ごとのσを表示）", fontsize=13)
fig.tight_layout()
plt.show()

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

6-10 行目3 つの設定をタプルで並べる
各タプルは「表示名・モデル・順変換・逆変換・特別処理」の 5 要素で、special="wls" のときだけ 17〜19 行目で sample_weight を渡します。OLS 生データ、WLS、OLS + log の 3 つに絞っているのは、残差の形の違いがこの 3 つで最もはっきり出るためです。
27-29 行目残差を x に対して散布図にする
y_te - pred をテストデータの \(x\) に対して描き、赤い破線で残差 0 を示します。等分散なら点は横一様な帯になり、不等分散なら右に開く扇形になります。この図は「残差プロット」と呼ばれ、OLS の前提を確認する最初の道具です。
36-42 行目区間ごとの σ を図に書き込む
\(x\) を 1〜7、7〜14、14〜20 の 3 区間に切り、各区間の残差の標準偏差を σ=… として図の上部に表示しています。OLS では左から右へ σ が 3 倍前後に増え、log 変換では 3 つの値が近づきます。数値で見ると「扇形」の印象を定量的に確かめられます。

残差パターン比較

OLS の残差は \(x\) が大きくなるにつれて扇形に広がり、3 区間の σ が左から右へ増えていきます（不等分散の典型）。WLS でも元の尺度の残差は右側で大きいままですが、直線の位置はノイズの小さい左側で決まっています。log 変換では 3 区間の σ がほぼ揃い、分散が均質化されたことが数値で確認できます。

6. よくある失敗パターン #

1 つ目は OLS のまま予測区間を計算することです。等分散を仮定した予測区間は残差全体の標準偏差 1 本で幅を決めるので、\(x = 20\) 付近では実際のばらつき（σ ≈ 10）より狭く、\(x = 1\) 付近では実際（σ ≈ 0.5）より広くなります。2 つ目は対数変換の逆変換でバイアスを無視することで、log(y) で学習した予測を expm1 で戻すと、期待値ではなく中央値の推定になります。期待値が必要なら Duan のスミアリング推定などでバイアス補正を行います。

線形回帰 — OLS の基本と仮定
Ridge・Lasso 回帰 — 正則化による係数安定化
ロバスト回帰 — Huber 損失と M 推定
分位点回帰 — 条件付き分位点の推定
MAE・RMSE — 回帰の基本指標