【ベンチマーク】非線形関係の回帰

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

まとめ

このベンチマークでは、\(y = \sin(x_1) + x_2^2 + \log(|x_3| + 1) + \epsilon\) という 3 種類の非線形項でできた 500 件のデータに、Linear・Ridge・SVR(RBF)・KNN・RF・GBM の 6 モデルと「なし・Poly(2)・Scaler・Poly+Scaler」の 4 前処理を掛け合わせ、5 分割交差検証の R² と MAE をヒートマップで比べます。線形モデルは前処理を変えても \(\sin\) や \(\log\) を直線で近似できず R² が低いままで、PolynomialFeatures(2) で \(x_2^2\) の項だけがわずかに改善します。RF・GBM は前処理なしで高い R² を出し、SVR(RBF) は StandardScaler を入れると大きく改善します。

物理現象・生体反応・価格弾力性など、入力と出力の関係が直線ではないケースは実務に多い。線形モデルに多項式特徴量を加えれば部分的に対応できるが、真の関数形が不明な場合は木系モデルやカーネル法が有利になる。ここでは既知の非線形関数で生成した合成データを使い、各アプローチの限界と強みを検証する。

1. データの特徴 #

500 件・5 変数で、各 \(x\) は標準正規分布 \(N(0, 1)\) から独立に引きます。真のモデルは \(y = \sin(x_1) + x_2^2 + \log(|x_3| + 1) + \epsilon\) で、\(\epsilon\) の標準偏差は 0.3 です。\(x_1\) は波、\(x_2\) は放物線、\(x_3\) は対数の形で \(y\) に効き、\(x_4\) と \(x_5\) は \(y\) と無関係なノイズ変数です。

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

コードを読む前に、この手法が図の上で何をしているかを追ってみましょう。「再生」を押すか、スライダーを動かすと 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
30
31
32
33
34
35
36
37
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

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

X = rng.normal(0, 1, size=(n, 5))
y = np.sin(X[:, 0]) + X[:, 1]**2 + np.log(np.abs(X[:, 2]) + 1) + rng.normal(0, 0.3, n)
feature_names = ["x₁", "x₂", "x₃", "x₄(noise)", "x₅(noise)"]

fig, axes = plt.subplots(1, 3, figsize=(15, 4))

# 各特徴量とyの関係
for i, (ax, fname, color) in enumerate(zip(axes, ["x₁ → sin(x₁)", "x₂ → x₂²", "x₃ → log(|x₃|+1)"],
                                            ["#2563eb", "#10b981", "#f97316"])):
    ax.scatter(X[:, i], y, alpha=0.3, s=10, color=color)

    # 真の関数を重ねる
    x_sort = np.sort(X[:, i])
    if i == 0:
        y_true = np.sin(x_sort) + 1.0  # 他の項の期待値を加算
    elif i == 1:
        y_true = x_sort**2 + 0.5
    else:
        y_true = np.log(np.abs(x_sort) + 1) + 1.0

    ax.plot(x_sort, y_true, color="#ef4444", linewidth=2, label="真の関数形")
    ax.set_title(fname)
    ax.set_xlabel(feature_names[i])
    ax.set_ylabel("y")
    ax.legend(fontsize=8)
    ax.grid(alpha=0.2)

fig.suptitle("非線形データの構造", fontsize=13)
fig.tight_layout()
plt.show()

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

8-10 行目3 種類の非線形項と 2 本のノイズ列
X[:, 0] は \(\sin\)、X[:, 1] は 2 乗、X[:, 2] は \(\log(|x| + 1)\) で \(y\) に入り、X[:, 3] と X[:, 4] は式に登場しません。ノイズの標準偏差 0.3 は \(y\) の散らばり（標準偏差 1.5 前後）より十分小さいので、R² の上限はほぼ 0.95 以上と見込めます。feature_names の “(noise)” は、後の図でどの列が無関係かを忘れないための印です。
15-17 行目1 変数ずつ y との散布図にする
3 枚のパネルに \(x_1\)、\(x_2\)、\(x_3\) それぞれと \(y\) の散布図を描きます。他の 2 項が足されているので点は縦に幅を持って散りますが、\(x_1\) では波、\(x_2\) では U 字、\(x_3\) では V 字に近い形が見えます。この「1 変数ずつ眺める」作業が、線形モデルで足りるかを判断する最初の手がかりです。
20-26 行目真の関数を他の項の期待値ぶんだけ持ち上げて重ねる
各パネルの赤い線は、その変数の項だけを描くと点の帯より下に来てしまうので、残り 2 項の期待値（\(E[x^2] = 1\) など）を定数として足して高さを合わせています。np.sort で \(x\) を並べ替えてから描くのは、plot が点を順につなぐので、並べ替えないと線がジグザグになるためです。

非線形データの構造

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

前処理	説明
なし	生の特徴量のみ
PolynomialFeatures(2)	2次の交互作用と二乗項を追加
StandardScaler	標準化のみ
Poly(2) + Scaler	多項式 + 標準化

モデル	特徴
Linear	最小二乗法
Ridge	L2正則化
SVR(RBF)	カーネル法（非線形対応）
KNN	近傍ベース
RF	ランダムフォレスト
GBM	勾配ブースティング

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
from sklearn.model_selection import cross_val_score, KFold
from sklearn.linear_model import LinearRegression, Ridge
from sklearn.svm import SVR
from sklearn.neighbors import KNeighborsRegressor
from sklearn.ensemble import RandomForestRegressor, GradientBoostingRegressor
from sklearn.preprocessing import StandardScaler, PolynomialFeatures
from sklearn.pipeline import Pipeline
from sklearn.base import clone
import seaborn as sns

cv = KFold(n_splits=5, shuffle=True, random_state=42)

def make_pipeline(mdl, poly=False, scale=False):
    steps = []
    if poly:
        steps.append(("poly", PolynomialFeatures(degree=2, include_bias=False)))
    if scale:
        steps.append(("scaler", StandardScaler()))
    steps.append(("reg", mdl))
    return Pipeline(steps) if len(steps) > 1 else mdl

models = {
    "Linear": LinearRegression(),
    "Ridge": Ridge(alpha=1.0),
    "SVR(RBF)": SVR(kernel="rbf", C=10),
    "KNN": KNeighborsRegressor(n_neighbors=10),
    "RF": RandomForestRegressor(n_estimators=100, random_state=42),
    "GBM": GradientBoostingRegressor(n_estimators=100, random_state=42),
}

preprocs = {
    "なし": {"poly": False, "scale": False},
    "Poly(2)": {"poly": True, "scale": False},
    "Scaler": {"poly": False, "scale": True},
    "Poly+Scaler": {"poly": True, "scale": True},
}

results = []
for prep_name, prep_cfg in preprocs.items():
    for mdl_name, mdl in models.items():
        pipe = make_pipeline(clone(mdl), **prep_cfg)
        try:
            mae_scores = -cross_val_score(pipe, X, y, cv=cv, scoring="neg_mean_absolute_error")
            r2_scores = cross_val_score(pipe, X, y, cv=cv, scoring="r2")
            results.append({
                "前処理": prep_name, "モデル": mdl_name,
                "MAE": mae_scores.mean(), "R²": r2_scores.mean(),
            })
        except Exception:
            results.append({
                "前処理": prep_name, "モデル": mdl_name,
                "MAE": np.nan, "R²": np.nan,
            })

df_res = pd.DataFrame(results)

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

pivot_r2 = df_res.pivot_table(index="前処理", columns="モデル", values="R²")
prep_order = ["なし", "Poly(2)", "Scaler", "Poly+Scaler"]
mdl_order = ["Linear", "Ridge", "SVR(RBF)", "KNN", "RF", "GBM"]
pivot_r2 = pivot_r2.reindex(index=[p for p in prep_order if p in pivot_r2.index],
                            columns=[m for m in mdl_order if m in pivot_r2.columns])

sns.heatmap(pivot_r2, annot=True, fmt=".3f", cmap="RdYlGn",
            linewidths=0.5, ax=axes[0], vmin=0, vmax=1,
            cbar_kws={"label": "R²（高いほど良い）"})
axes[0].set_title("前処理 × モデル: R² ヒートマップ")
axes[0].set_xlabel("")
axes[0].set_ylabel("")

pivot_mae = df_res.pivot_table(index="前処理", columns="モデル", values="MAE")
pivot_mae = pivot_mae.reindex(index=[p for p in prep_order if p in pivot_mae.index],
                              columns=[m for m in mdl_order if m in pivot_mae.columns])

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

fig.suptitle("非線形データ: 前処理 × モデル 精度比較", fontsize=13)
fig.tight_layout()
plt.show()

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

13-20 行目前処理の有無で Pipeline を組み替える
poly=True なら PolynomialFeatures(degree=2, include_bias=False)、scale=True なら StandardScaler を手前に差し込み、どちらも無ければモデルをそのまま返します。2 次の多項式特徴量は 5 変数から 20 列（5 本の 1 次 + 5 本の 2 乗 + 10 本の交差項）を作るので、\(x_2^2\) の項だけは線形モデルでも表現できるようになります。
22-29 行目6 モデルのハイパーパラメータは固定
SVR(C=10)、KNeighborsRegressor(n_neighbors=10)、RandomForestRegressor(n_estimators=100)、GradientBoostingRegressor(n_estimators=100) と、すべて手で決めた値で比べています。チューニングすれば差は縮まりますが、ここでは「前処理なし・既定値寄りの設定で、どの手法が非線形を拾えるか」を見るのが目的です。
43-44 行目MAE と R² を同じ分割で測る
cross_val_score を 2 回呼んでいますが、cv が random_state=42 の KFold なので分割は同じです。MAE は neg_mean_absolute_error で返るため先頭にマイナスを付けて正の値に戻しています。計算は 2 回走るので、モデルが重い場合は cross_validate で複数の指標を一度に取る方が速いです。
65-67 行目R² の色尺度を 0〜1 に固定する
vmin=0, vmax=1 で色の範囲を固定しているので、別のデータで同じコードを走らせても色の意味が変わりません。固定しないと、全モデルが低い R² でも相対的に濃い緑が出てしまい、「よく当たっている」と誤読します。MAE 側は尺度が問題ごとに違うので固定していません。

精度比較ヒートマップ

読み方のポイント #

線形モデル（Linear・Ridge）の列は 4 つの前処理のどれでも R² が低いままで、\(\sin\) や \(\log\) は 1 次式では近似できないことが分かります。PolynomialFeatures(2) を足した行では \(x_2^2\) の項を捉えられるため、線形モデルの R² がわずかに上がります。RF・GBM の列は前処理なしでも高い R² を出し、木の分岐が区間ごとの定数で非線形関数を自動的に近似しているためです。SVR(RBF) は Scaler の行で大きく改善し、カーネル法が距離ベースなのでスケーリングが効くことを示します。

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
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 = [
    ("Linear（生データ）", LinearRegression(), False, False),
    ("Ridge + Poly(2)", Ridge(alpha=1.0), True, False),
    ("GBM（生データ）", GradientBoostingRegressor(n_estimators=100, random_state=42), False, False),
]

fig, axes = plt.subplots(1, 3, figsize=(15, 4))
for i, (name, mdl, poly, scale) in enumerate(configs):
    pipe = make_pipeline(clone(mdl), poly=poly, scale=scale)
    pipe.fit(X_tr, y_tr)
    pred = pipe.predict(X_te)
    residuals = y_te - pred

    axes[i].scatter(pred, residuals, alpha=0.4, s=15, color="#2563eb")
    axes[i].axhline(0, color="#ef4444", linewidth=1, linestyle="--")
    r2 = 1 - np.sum(residuals**2) / np.sum((y_te - np.mean(y_te))**2)
    axes[i].set_title(f"{name}\nR²={r2:.3f}", fontsize=10)
    axes[i].set_xlabel("予測値")
    axes[i].set_ylabel("残差")
    axes[i].grid(alpha=0.2)

fig.suptitle("残差パターンの比較（テストデータ）", fontsize=13)
fig.tight_layout()
plt.show()

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

6-10 行目段階の違う 3 設定を選ぶ
生の線形回帰、2 次の多項式特徴量を足した Ridge、前処理なしの GBM の 3 つで、「何もしない → 手で非線形項を足す → モデルに任せる」の 3 段階を並べています。前のブロックのヒートマップで差が大きかった組み合わせを抜き出した形です。
19-22 行目横軸を予測値にした残差プロット
pred を横軸、y_te - pred を縦軸に取っています。モデルが正しければ残差は予測値と無関係に 0 の周りに散るはずで、曲線状の構造が見えたら「捉え損ねた非線形成分」が残っている証拠です。タイトルに書き込む R² は 21 行目で手計算しており、r2_score と同じ値になります。

残差パターン比較

線形モデルの残差には予測値に沿った曲線状の構造が残り、これが「モデルが捉えきれていない非線形成分」の証拠です。Ridge + Poly(2) では \(x_2^2\) の分だけ構造が弱まりますが、\(\sin\) と \(\log\) の成分は残ります。GBM の残差は 0 の周りにランダムに散らばり、系統的なパターンがありません。

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

1 つ目は残差プロットを見ずに R² だけで判断することです。R² = 0.7 でも残差に曲線パターンがあれば非線形性を見落としていて、特徴量を足すか木系モデルに替えれば R² はまだ上がります。2 つ目は PolynomialFeatures の次数を上げすぎることで、5 変数でも degree=3 で 55 列、degree=4 で 125 列と特徴量数が膨らみ、過学習のリスクが急増します。次数は交差検証で選びます。

多項式回帰 — 多項式特徴量の理論と実践
サポートベクター回帰 — カーネルトリックによる非線形回帰
R² — 決定係数の解釈
MAE・RMSE — 回帰の基本指標
学習曲線 — 過学習の診断