多項式曲線フィッティング (2)： 過学習を防ぐ「正則化」項とscikit-learn利用の暫定まとめ¶

このページのオリジナルのipynbファイル
Return to home page

(先週話したことの補足）¶

数学的な詳細の理解には時間がかかるが、

2乗誤差最小化問題という多変数関数の極値（最小値）問題が線形代数の（行列）計算に帰着できる

ということに注目してほしい。これは評価関数が２次関数の線形結合になっている場合に限られるが、「計画行列」の活用は後のカーネル法などにも引き継がれる。

(注) 関数の極値問題は勾配法(一変数の場合はニュートン法と呼ばれる)による繰り返し計算を行うことが多い(たとえば、https://fussy.web.fc2.com/algo/math10_extreme.htm)。

多数の局所的な極小がたくさんある場合は困難であり、様々なシミュレーション的方法が行われている。

今週の話¶

"IPython Interactive Computing and Visualization Cookbook" (O'Reilly, 2018)のサンプルプログラム8.1を例に （現在、原文はhttps://ipython-books.github.io/ にて閲覧できる。）

模式図はPRMLより。

2乗誤差最小化 (復習)¶

説明変数$x$に対する目的変数$y$について、$M$次関数でのフィッティングを考える： $$ y(x, {\bf w}) = w_0 + w_1 x + w_2 x^2 + \cdots + w_M x^M = \sum_{j=0}^{M} w_j x^j $$ データ$\{(x_n, t_n)\}$ ($n=1,\cdots N$)があったときに、 2乗誤差 $$ E({\bf w}) = \frac{1}{2} \sum_{n=1}^N \{y(x_n,{\bf w}) - t_n\}^2 $$ を最小にするよう${\bf w}$を決める。

説明用に、PRMLに載っている図を載せておく。

最小２乗法の模式図：

$\sin$カーブにノイズを載せたテストデータをべき関数フィッティングしてみた例：次数$M$をが大きいほうが全データにあうフィッティングができるが、不自然であることがわかるだろう。($M$が次数)

一般に高次項の係数がとても大きくなる。(前回の関数フィッティングのプログラムで係数を表示してみてほしい。function_fitting.htmlの最後のセルにあるプリント文のコメントを外せばよい。）

それを抑える工夫が次に述べるリッジ回帰、ラッソ回帰と呼ばれるものである。

リッジ回帰¶

過学習を防ぐために、次のような係数に対する2次関数を付加し、その最小化を図る。この項は大きな係数の方が不利になるので、高次の項の「暴発」が抑えられる。

$$ \tilde{E}({\bf w}) = \frac{1}{2} \sum_{n=1}^N \{y(x_n,{\bf w}) - t_n\}^2 + \frac{\lambda}{2}||{\bf w}||^2 $$

これはリッジ回帰(ridge regression)を呼ばれる。ニューラルネットでは、荷重減衰(weight decay)として知られている。

また、このことを正則化(regularization)と呼ぶ。

リッジ回帰における重み係数¶

上記のPRMLの事例にある式で係数$\lambda$を固定して最小化を行った例($M=9$)：$\lambda$が大きいと加重がきつくフィッティングされにくくなる。リッジ回帰では$\lambda$を手入力で与える必要がある。（下で試すscikit-learnではalphaという引数で指定）

pythonスクリプトでの実例¶

サンプルデータの作成¶

まずは前回同様、サンプルデータを作っておく。 (sinカーブでなくてもよい）

ここでは、無名関数fを定義し、それにランダムデータを加える形でデータを作成してみる。

In [1]:
# 今まで試してきた多項式近似の復習
import numpy as np
import random
import sklearn.linear_model as lm
import matplotlib.pyplot as plt

f_cos = lambda x: np.cos(2.0*x) # サンプルデータを作るためのベース関数: このような関数定義の記述方法もある

# データ数
n_tr = 20
x_max = 5.0 # xの範囲 [0, x_max]
x = np.linspace(0., x_max, n_tr)
y = f_cos(x) + np.random.randn(len(x))*0.3


# 単純な多項式近似
lrp = lm.LinearRegression()

for deg in [2,5,10]: # 最大べき
    power_matrix_x = np.vander(x, deg+1) # 計画行列の作成
    lrp.fit(power_matrix_x, y)
    # モデルの係数表示
    print('フィッティング関数の次数=' + str(deg))
    print('切片= ' + str(lrp.intercept_))
    print('係数= ' +  ' '. join(['%.3f' % c for c in lrp.coef_]))
    # 予測
    x_lrp = np.linspace(0., x_max*1.2, 100)
    y_lrp = lrp.predict(np.vander(x_lrp, deg+1))
    # 近似曲線の描画
    plt.plot(x_lrp, y_lrp)

# データの描画
plt.plot(x, y, "ob", ms=8)
plt.ylim(-2.0,2.0)

フィッティング関数の次数=2
切片= 0.13568625562071746
係数= -0.020 -0.024 0.000
フィッティング関数の次数=5
切片= 1.2210868946388487
係数= 0.086 -0.993 3.681 -4.329 -0.410 0.000
フィッティング関数の次数=10
切片= 1.1220577110543857
係数= -0.001 0.028 -0.304 1.835 -6.863 16.512 -26.095 27.380 -17.010 2.611 0.000
Out[1]:
(-2.0, 2.0)

scikit-learnのリッジ回帰¶

マニュアルページ： https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.Ridge.html

In [2]:
# リッジ回帰 alpha(正則化係数)の既定値は1.0  = 上の解説文におけるλ
ridge = lm.Ridge(alpha=1.0)

for deg in [2,5,10]: # 最大べき
    ridge.fit(np.vander(x, deg +1), y) # リッジ回帰
    # 予測
    x_lrp = np.linspace(0., x_max*1.2, 100)
    y_lrp = ridge.predict(np.vander(x_lrp, deg+1))
    plt.plot(x_lrp, y_lrp,
            label='degree ' + str(deg))
    plt.plot(x, y , 'ok', ms=10)
    plt.title('Ridge Regression')
    # モデルの係数表示
    # print('フィッティング関数の次数=' + str(deg))
    # print('切片= ' + str(ridge.intercept_))
    # print('  '. join(['%.5f' % c for c in ridge.coef_]))
    # 「真の値」cos(x)との比較による決定係数R2の算出
    # 実際には、真の値はわからない。その場合、トレーニングデータとは別にテストデータを用意し、それで検証する
    #print("alpha=" + str(alpha) + "R2=" + str(ridge.score(np.vander(x, deg+1), f_cos(x))))

# サンプルデータのプロット
plt.ylim(-2.0,2.0)
plt.plot(x, y, "ob", ms=8)
# sin curve
plt.plot(x_lrp, f_cos(x_lrp), label="sin(x)")
plt.legend(loc=2)
Out[2]:
<matplotlib.legend.Legend at 0x7f75a1100fa0>
In [5]:
# リッジ回帰 alpha(正則化係数)の既定値は1.0 

for a in [0.05, 0.1, 0.3, 0.5, 1.0, 2.0]:
    ridge = lm.Ridge(alpha=a)
    
    for deg in [2, 5, 10]: # 最大べき
        ridge.fit(np.vander(x, deg +1), y) # リッジ回帰
        # 予測
        x_lrp = np.linspace(0., x_max*1.2, 100)
        y_lrp = ridge.predict(np.vander(x_lrp, deg+1))
        plt.plot(x_lrp, y_lrp,
                label='degree ' + str(deg))
        plt.legend(loc=2)
        plt.plot(x, y , 'ok', ms=10)
        plt.title('Ridge Regression')
        # モデルの係数表示
        print('フィッティング関数の次数=' + str(deg))
        print('切片= ' + str(ridge.intercept_))
        print('  '. join(['%.5f' % c for c in ridge.coef_]))
        # 「真の値」cos(x)との比較による決定係数R2の算出
        # 実際には、真の値はわからない。その場合、トレーニングデータとは別にテストデータを用意し、それで検証する
        print("alpha=" + str(a) + ", deg=" + str(deg) + ", R2=" + str(ridge.score(np.vander(x, deg+1), f_cos(x))))

# サンプルデータのプロット
plt.ylim(-2.0,2.0)
plt.plot(x, y, "ob", ms=8)

フィッティング関数の次数=2
切片= 0.4164497190283239
0.01799  -0.22651  0.00000
alpha=0.05, deg=2, R2=0.07319032622357546
フィッティング関数の次数=5
切片= 1.4212147654535712
0.06024  -0.68241  2.41715  -2.35013  -1.39076  0.00000
alpha=0.05, deg=5, R2=0.9246990582148639
フィッティング関数の次数=10
切片= 1.1774753605356092
0.00004  -0.00187  0.02135  -0.09331  0.13959  0.01035  0.04015  -0.05689  -0.53787  -0.93996  0.00000
alpha=0.05, deg=10, R2=0.9279894262666633
フィッティング関数の次数=2
切片= 0.41324330712538765
0.01733  -0.22297  0.00000
alpha=0.1, deg=2, R2=0.07365235557690586
フィッティング関数の次数=5
切片= 1.3956163105706856
0.05194  -0.58207  2.00201  -1.68759  -1.68644  0.00000
alpha=0.1, deg=5, R2=0.9099839189650605
フィッティング関数の次数=10
切片= 1.1369338135777762
0.00005  -0.00203  0.02172  -0.08865  0.10783  0.07840  0.02311  -0.15134  -0.52052  -0.82604  0.00000
alpha=0.1, deg=10, R2=0.9273969899665108
フィッティング関数の次数=2
切片= 0.4013772751146817
0.01488  -0.20985  0.00000
alpha=0.3, deg=2, R2=0.075297455097181
フィッティング関数の次数=5
切片= 1.1803628262120314
0.03785  -0.42233  1.41528  -0.99638  -1.61719  0.00000
alpha=0.3, deg=5, R2=0.8613623948014627
フィッティング関数の次数=10
切片= 1.045376268038547
0.00010  -0.00297  0.02865  -0.11117  0.12832  0.12332  -0.03555  -0.24795  -0.47652  -0.62312  0.00000
alpha=0.3, deg=10, R2=0.9219499395256724
フィッティング関数の次数=2
切片= 0.39085105183565494
0.01271  -0.19822  0.00000
alpha=0.5, deg=2, R2=0.07667142623387724
フィッティング関数の次数=5
切片= 1.0120046107744969
0.03038  -0.34092  1.14105  -0.76578  -1.40090  0.00000
alpha=0.5, deg=5, R2=0.8090359941215963
フィッティング関数の次数=10
切片= 0.9911163642329129
0.00010  -0.00308  0.03004  -0.11841  0.14266  0.12584  -0.06313  -0.27105  -0.44588  -0.53172  0.00000
alpha=0.5, deg=10, R2=0.9158860358061822
フィッティング関数の次数=2
切片= 0.369088997747132
0.00823  -0.17418  0.00000
alpha=1.0, deg=2, R2=0.07925710154300392
フィッティング関数の次数=5
切片= 0.7548628842693
0.02008  -0.23028  0.77993  -0.50957  -1.02103  0.00000
alpha=1.0, deg=5, R2=0.6971672203411846
フィッティング関数の次数=10
切片= 0.9027236043260879
0.00006  -0.00236  0.02623  -0.11204  0.14863  0.11277  -0.08989  -0.27591  -0.39200  -0.41894  0.00000
alpha=1.0, deg=10, R2=0.9008057664473685
フィッティング関数の次数=2
切片= 0.3384880212715846
0.00195  -0.14044  0.00000
alpha=2.0, deg=2, R2=0.08230979030066121
フィッティング関数の次数=5
切片= 0.5210910480469141
0.01125  -0.13596  0.47710  -0.31679  -0.65662  0.00000
alpha=2.0, deg=5, R2=0.5614928171903719
フィッティング関数の次数=10
切片= 0.7943102700898788
-0.00005  -0.00053  0.01501  -0.08398  0.13406  0.08279  -0.10212  -0.25133  -0.32304  -0.31727  0.00000
alpha=2.0, deg=10, R2=0.8728212655743505
Out[5]:
[<matplotlib.lines.Line2D at 0x7f4c2dfb1b50>]

Lasso回帰¶

ウェイトの項を単に絶対値$\frac{\lambda}{2}||{\bf w}||$にしたものをlasso回帰と呼ぶ。

In [ ]:
model = lm.Lasso(alpha=0.01)
for deg in [2,5,10]: # 最大べき
    model.fit(np.vander(x, deg +1), y) # リッジ回帰
    # 予測
    x_lrp = np.linspace(0., x_max*1.2, 100)
    y_lrp = model.predict(np.vander(x_lrp, deg+1))
    plt.plot(x_lrp, y_lrp,
            label='degree ' + str(deg))
    plt.legend(loc=2)
    # モデルの係数表示
    # モデルの係数表示
    print('フィッティング関数の次数=' + str(deg))
    print('切片= ' + str(lrp.intercept_))
    print('係数' + ' '. join(['%.5f' % c for c in model.coef_]))
    plt.plot(x, y , 'ok', ms=10)
    plt.title('Lasso Regression')

# サンプルデータのプロット
plt.ylim(-2.0,2.0)
plt.plot(x, y, "ob", ms=8)

Lassoの場合、この程度のばらつきの場合だと、既定の加重係数(alpha=1.0)は大きすぎるようだ。また、データ数を多くしないと収束がわるという警告がでる。

全体として、次数を上げた場合、加重を含む最適化を行っているため、過学習が抑えられている。

係数alpha (上の式では$\lambda$)は外からあたえるパラメータであり、ハイパーパラメータと呼ばれる。

Exercise 3¶

  • データ作成のベース関数やデータ数、ハイパーパラメータ(正則化係数alphaやフィッティング関数の次数)を変えて振る舞いを確かめてみよう。
  • ランダム項を含まないcos(x)を再現(推定)できているかを決定係数R2を用いて評価してみよう。

(注) 決定係数

テストデータの実測値を$(t_1, t_2, \cdots, t_N)$、予測(推定)値を$(y_1, y_2, \cdots, y_N)$としたとき

$$ R^2 = 1 - \frac{\sum_{i=1}^N (y_i - t_i)^2}{ \sum_{i=1}^N (t_i - \bar{t})^2} $$

scikit-learnの利用まとめ Summary on the scikit-learn fitting libraries¶

The regression problem we solve (predict) is to make a numerical machine that predict the value $y$ from the explanation variables $\{\boldsymbol{x}\}$ from a training dataset $\{\boldsymbol{x}_n, y\}$.

$$ y = f(\{\boldsymbol{x}_n\}, \boldsymbol{x}) $$

where the values of the dataset is regarded as parameters of $f(\boldsymbol{x})$.

  • Select a model xxxxx (one of linear, Ridge, Lasso, support vector machine, neural net, and etc)
    import sklearn.linear_model as lm
model = lm.xxxxxx(parameter1, parameter2, ....)

  • Prepare a training dataset: array "x" of explanation variables and list of target "y"

  • Train the machine by the dataset

    model.fit(x,y, ....)

  • Predict $y$ from a set of explanation variables $\{x_1, x_2, \cdots \}$.

    model.predict(x)

補足 : 外れ値について¶

極端な外れ値を含む場合は、別の対応が必要である。

参考：scikit-learnの事例より

https://scikit-learn.org/stable/auto_examples/linear_model/plot_huber_vs_ridge.html#sphx-glr-auto-examples-linear-model-plot-huber-vs-ridge-py

In [14]:
# Authors: Manoj Kumar mks542@nyu.edu
# License: BSD 3 clause

import numpy as np
import matplotlib.pyplot as plt

from sklearn.datasets import make_regression
from sklearn.linear_model import HuberRegressor, Ridge

def huber_vs_ridge(show_ridge=True, show_huber=True):
    # Generate toy data.
    rng = np.random.RandomState(0)
    X, y = make_regression(
        n_samples=20, n_features=1, random_state=0, noise=4.0, bias=100.0
    )
    
    # Add four strong outliers to the dataset.
    X_outliers = rng.normal(0, 0.5, size=(4, 1))
    y_outliers = rng.normal(0, 2.0, size=4)
    X_outliers[:2, :] += X.max() + X.mean() / 4.0
    X_outliers[2:, :] += X.min() - X.mean() / 4.0
    y_outliers[:2] += y.min() - y.mean() / 4.0
    y_outliers[2:] += y.max() + y.mean() / 4.0
    X = np.vstack((X, X_outliers))
    y = np.concatenate((y, y_outliers))
    plt.plot(X, y, "b.")
    x = np.linspace(X.min(), X.max(), 7)
    
    # Fit the huber regressor over a series of epsilon values.
    if show_huber:
        colors = ["r-", "b-", "y-", "m-"]
        epsilon_values = [1, 1.5, 1.75, 1.9]
        for k, epsilon in enumerate(epsilon_values):
            huber = HuberRegressor(alpha=0.0, epsilon=epsilon)
            huber.fit(X, y)
            coef_ = huber.coef_ * x + huber.intercept_
            plt.plot(x, coef_, colors[k], label="huber loss, %s" % epsilon)
    
    # Fit a ridge regressor to compare it to huber regressor.
    if show_ridge:
        ridge = Ridge(alpha=0.0, random_state=0)
        ridge.fit(X, y)
        coef_ridge = ridge.coef_
        coef_ = ridge.coef_ * x + ridge.intercept_
        plt.plot(x, coef_, "g-", label="ridge regression")
    
    plt.title("Comparison of HuberRegressor vs Ridge")
    plt.xlabel("X")
    plt.ylabel("y")
    if show_ridge or show_huber:
        plt.legend(loc=0)
    plt.show()
In [19]:
huber_vs_ridge()
In [18]:
# Using "interact" to switch on or off showing results of two kind of regressions 

from ipywidgets import interact
interact(huber_vs_ridge, show_ridge=False, show_huber=False)

interactive(children=(Checkbox(value=False, description='show_ridge'), Checkbox(value=False, description='show…
Out[18]:
<function __main__.huber_vs_ridge(show_ridge=True, show_huber=True)>
In [20]:
%%html
<link rel="stylesheet" type="text/css" href="custom.css">
In [ ]: