-Python- リッジ回帰

線形回帰で最小化する目的関数に,パラメータの大きさの項を足して回帰を行うのがリッジ回帰です.リッジ回帰では以下の関数を最小化するようなwを決定します.

この計算を実装し,可視化するスクリプトは以下のようになります.

 
# Import Module
import numpy as np
from scipy import linalg
 
# Definition of Class
class RidgeRegression: 
    def __init__(self, lambda_ = 1.): 
        self.lambda_ = lambda_
        self.w_ = None
    
    def fit(self, X, t): 
        Xtil = np.c_[np.ones(X.shape[0]), X]
        c = np.eye(Xtil.shape[1])
        A = np.dot(Xtil.T, Xtil) + self.lambda_ * c
        b = np.dot(Xtil.T, t)
        self.w_ = linalg.solve(A, b)
    
    def predict(self, X): 
        Xtil = np.c_[np.ones(X.shape[0]), X]

 

        return np.dot(Xtil, self.w_)
 
 
# Import Module
import ridge
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d
 
x = np.array([1, 2, 4, 6, 7])
y = np.array([1, 3, 3, 5, 4])
model = ridge.RidgeRegression(1.)
model .fit(x, y)
b, a = model.w_
 
# Drawing Graph
plt.scatter(x, y, color = "k")
xmax = x.max()
plt.plot([0, xmax], [b, b + a * xmax], color = "k")
 
plt.show()
 
# Solving problem
n = 100
scale = 10
 
np.random.seed(0)
X = np.random.random*1 * scale
w0 = 1
w1 = 2
w2 = 3
y = w0 + w1 * X[:, 0] + w2 * X[:, 1] + np.random.randn(n)
 
model = ridge.RidgeRegression(1.)
model.fit(X, y)
 
# Drawing Graph
xmesh, ymesh = np.meshgrid(np.linspace(0, scale, 20), 
                            np.linspace(0, scale, 20))
zmesh = (model.w_[0] + model.w_[1] * xmesh.ravel() +
        model.w_[2] * ymesh.ravel()).reshape(xmesh.shape)
fig = plt.figure()
ax = fig.add_subplot(111, projection = '3d')
ax.scatter(X[:, 0], X[:, 1], y, color = "k")
ax.plot_wireframe(xmesh, ymesh, zmesh, color = "r")
plt.show()
 
実行結果は以下のようになります.
続いて普通の線形回帰との違いを確かめるために以下のようなスクリプトを実行してみます.

 

# Import Module

import ridge

import linearreg

import numpy as np

import matplotlib.pyplot as plt

 

# Solving problem

x = np.arange(12)

y = 1 + 2 * x

y[2] = 20

y[4] = 0

 

# Drawing Graph

xmin = 0

xmax = 12

ymin = -1

ymax = 25

fig, axes = plt.subplots(nrows = 2, ncols = 5)

for i in range(5): 

    axes[0, i].set_xlim([xmin, xmax])

    axes[0, i].set_ylim([ymin, ymax])

    axes[1, i].set_xlim([xmin, xmax])

    axes[1, i].set_ylim([ymin, ymax])

    xx = x[:2 + i * 2]

    yy = y[:2 + i * 2]

    axes[0, i].scatter(xx, yy, color = "k")

    axes[1, i].scatter(xx, yy, color = "k")

    model = linearreg.LinearRegression()

    model.fit(xx, yy)

    xs = [xmin, xmax]

    ys = [model.w_[0] + model.w_[1] * xmin, 

            model.w_[0] + model.w_[1] * xmax]

    axes[0, i].plot(xs, ys, color = "k")

    model = ridge.RidgeRegression(10.)

    model.fit(xx, yy)

    xs = [xmin, xmax]

    ys = [model.w_[0] + model.w_[1] * xmin, 

            model.w_[0] + model.w_[1] * xmax]

    axes[1, i].plot(xs, ys, color = "k")

 

plt.show()

 

実行結果は以下のようになります.

f:id:HidehikoMURAO:20181112212527p:plain

上記の実行結果では,ほぼ線形に並んでいる点列を使いますが,2つだけ大きくずれるような値になるよう人工的に配置したデータを用いています.データには12個の点が含まれているのですが,2個からスタートして2個ずつ増やして回帰の様子がどうなるのかを図示しています.

実行結果の上の行が普通の線形回帰で,下の行がリッジ回帰です.1番左の列が点が2個の場合で,それから右に進むにつれて点が2個ずつ増えていきます.3つ目の点が大きく上に外れているので,その値に影響を受けて左からの2列目の図では回帰の直線が立ち上がっているようになっています.

左から2列目では上の方が下の図より立ち上がり方が大きいように感じます.一方5つ目の点が左に大きく外れていますが,その影響で3列目の図では直線の傾きが小さくなっています.

左から3列目を見ると,上の方が下の図より傾きが小さくなっています.全体として上の線形回帰と下のリッジ回帰を比べると,上の図の方が点が増えるにつれて直線の傾きが大きく変化しているように感じます.下の図では比較的外れ値から受ける影響が小さくなっています.

リッジ回帰はサンプル数が少ないときに例外的なデータからの影響を受けにくいという性質があります.これは,正則化による効果です.

*1:n, 2