（2021-01-16 追記: 回帰係数の推定量 (\(\hat{\beta}\)) についての記述を、予測値についての記述と勘違いしたまま、 的外れな妄言を書き連ねていた。以下に書いている話は、ほぼ価値が無い）

色々あってスパース回帰分析を勉強しており、講談社データサイエンス入門シリーズの 「スパース回帰分析とパターン認識」 を読んでいる。まったく理解に苦労しており、全然読めていないのだけれど。

ともあれ、その p.16 に次のような記述があった:

リッジ推定量は x に対して傾き \(\frac{1}{1-λ}\) の直線であることから， 大きな x に対しては軟しきい値作用素や硬しきい値作用素と比べて， 最小２乗推定量との乖離が大きくなることが見て取れる．

リッジ回帰では推定量が過少になるということなのだろうか。 取り急ぎ、プログラムを書いて少し様子を見ていた。

X の大きさに応じて大きくなるようなノイズを直線に載せて作ったトイデータに対して statsmodels と scikit-learn の両方で、 Ridge と Lasso のモデルを学習させて回帰直線を描いてみる。 すると、statsmodels では次のような結果になったが:

scikit-learn では次のような結果になった:

x が大きくなるほど最小二乗法 (OLS) の結果から乖離していく様子は λ が大きいほど強く表れているように見える。 また、Ridge の場合は乖離具合が数理的に求められるのかな、 scikit-learn の Ridge では λ や x の大小に関わらず OLS と同等の値が予測されている。

若干話がそれるけれど、こうしてみると scikit-learn の Ridge は statsmodels の実装と動きが違うということだね。 両者は異なる数式でモデルを定義しているので、ちゃんと見ておくべきだな。

使ったソースは以下のようなもの。

import matplotlib.pyplot as plt
import numpy as np
import sklearn.linear_model
import statsmodels.api as sm


# %% データ生成
def sample_data(X):
    X = np.asarray(X)
    if X.ndim == 1:
        (n_features,) = X.shape
        w = np.array([1.3])
        y = X @ w
        y += np.random.rand() * 4.0 * X[0]
        return y
    elif X.ndim == 2:
        n_samples, n_features = X.shape
        return np.array([sample_data(X[i, :]) for i in range(n_samples)])


X = np.random.rand(1, 256).T * 10
assert X.flags.f_contiguous
y = sample_data(X)

fig, ax = plt.subplots()
ax.scatter(X[:, 0], y)
ax.set_xlabel("x")
ax.set_ylabel("y")
fig.tight_layout()

# %% statsmodels の Ridge と Lasso で、λ と推定量の様子を見る
ols = sm.OLS(y, X, hasconst=False).fit()
ridge2 = sm.OLS(y, X, hasconst=False).fit_regularized(alpha=2, L1_wt=0)
ridge8 = sm.OLS(y, X, hasconst=False).fit_regularized(alpha=8, L1_wt=0)
ridge32 = sm.OLS(y, X, hasconst=False).fit_regularized(alpha=32, L1_wt=0)
lasso2 = sm.OLS(y, X, hasconst=False).fit_regularized(alpha=2, L1_wt=1)
lasso8 = sm.OLS(y, X, hasconst=False).fit_regularized(alpha=8, L1_wt=1)
lasso32 = sm.OLS(y, X, hasconst=False).fit_regularized(alpha=32, L1_wt=1)

fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(10, 5))
x_ = np.linspace(0, 10).reshape(-1, 1)
ax1.scatter(X[:, 0], y, marker=".")
ax1.plot(x_, ols.predict(x_), label="OLS")
ax1.plot(x_, ridge2.predict(x_), label="Ridge (λ=2.0)")
ax1.plot(x_, ridge8.predict(x_), label="Ridge (λ=8.0)")
ax1.plot(x_, ridge32.predict(x_), label="Ridge (λ=32.0)")
ax2.scatter(X[:, 0], y, marker=".")
ax2.plot(x_, ols.predict(x_), label="OLS")
ax2.plot(x_, lasso2.predict(x_), label="Lasso (λ=2.0)")
ax2.plot(x_, lasso8.predict(x_), label="Lasso (λ=8.0)")
ax2.plot(x_, lasso32.predict(x_), label="Lasso (λ=32.0)")
ax1.legend()
ax2.legend()
fig.tight_layout()
plt.show()

# %% scikit-learn の Ridge と Lasso で、λ と推定量の様子を見る
ols = sklearn.linear_model.LinearRegression(fit_intercept=False).fit(X, y)
ridge2 = sklearn.linear_model.Ridge(2, fit_intercept=False).fit(X, y)
ridge8 = sklearn.linear_model.Ridge(8, fit_intercept=False).fit(X, y)
ridge32 = sklearn.linear_model.Ridge(32, fit_intercept=False).fit(X, y)
lasso2 = sklearn.linear_model.Lasso(2, fit_intercept=False).fit(X, y)
lasso8 = sklearn.linear_model.Lasso(8, fit_intercept=False).fit(X, y)
lasso32 = sklearn.linear_model.Lasso(32, fit_intercept=False).fit(X, y)

fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(10, 5))
x_ = np.linspace(0, 10).reshape(-1, 1)
ax1.scatter(X[:, 0], y, marker=".")
ax1.plot(x_, ols.predict(x_), label="OLS")
ax1.plot(x_, ridge2.predict(x_), label="Ridge (λ=2.0)")
ax1.plot(x_, ridge8.predict(x_), label="Ridge (λ=8.0)")
ax1.plot(x_, ridge32.predict(x_), label="Ridge (λ=32.0)")
ax2.scatter(X[:, 0], y, marker=".")
ax2.plot(x_, ols.predict(x_), label="OLS")
ax2.plot(x_, lasso2.predict(x_), label="Lasso (λ=2.0)")
ax2.plot(x_, lasso8.predict(x_), label="Lasso (λ=8.0)")
ax2.plot(x_, lasso32.predict(x_), label="Lasso (λ=32.0)")
ax1.legend()
ax2.legend()
fig.tight_layout()
plt.show()

以上。