ここでは，前のページで説明したリッジ回帰を行うプログラムを自作してみよう．次のページで説明するように，scikit-learn を使えば短いコードで簡単に分析を行うことができますが，ブラックボックスとして scikit-learn を使うよりもある程度中身を理解して scikit-learn を利用することをおすすめします．ここの手順でプログラムを記述することでより理解が深まることでしょう．
まずは必要なライブラリのインポートを行います．ここでは目的関数を最小化するために scipy.optimize を利用することにします．
ライブラリの読み込み
import pandas as pd
import numpy as np
import scipy.optimize as optimize
次に，回帰式とリッジ回帰の目的関数を求めるための関数を定義します．
関数の定義
"""
回帰式の定義
"""
def my_func(w, x):
return np.dot(x, w)
"""
リッジ回帰の目的関数
"""
def get_ridge(w, x, y, alpha=1.0):
y_pred = my_func(w, x)
error = (y - y_pred)**2
return np.sum(error) + alpha * np.dot(w, w)
次に GitHub のリポジトリから CSV ファイルを読み込みます．
データファイルの読み込み
url = "https://github.com/rinsaka/sample-data-sets/blob/master/ridge-train.csv?raw=true"
df = pd.read_csv(url)
print(df)
no x1 x2 y 0 0 9.852 19.712 205.060 1 1 8.516 16.936 179.572 2 2 0.416 0.816 20.035 3 3 7.852 15.671 168.983 4 4 5.855 11.800 127.056 5 5 4.463 8.936 99.450 6 6 0.623 1.342 22.804 7 7 8.490 17.023 178.753 8 8 4.143 8.298 93.451 9 9 0.800 1.531 27.555
次に，行列 \(\bf X\) のデータを Pandas データフレームから NumPy 配列に変換します．このとき，推定したいパラメータベクトル \(\bf{w}\) は3行の行ベクトルであることから，行列 \(\bf X\) が3列になるように，先頭に値1の列を追加します．これによって，回帰式の関数 my_func(w, x) の行列演算 \({\bf X}\bf{w}\) が5行目の np.dot(x, w) で計算できるようになります．
NumPy 配列に変換し列を追加
x_data = df.loc[:,['x1', 'x2']].values
x_data = np.insert(x_data, 0, 1, axis=1) # 先頭に列を追加する
print(x_data)
[[ 1. 9.852 19.712] [ 1. 8.516 16.936] [ 1. 0.416 0.816] [ 1. 7.852 15.671] [ 1. 5.855 11.8 ] [ 1. 4.463 8.936] [ 1. 0.623 1.342] [ 1. 8.49 17.023] [ 1. 4.143 8.298] [ 1. 0.8 1.531]]
さらに，\(y\) データも NumPy 配列に変換しておきます．
NumPy 配列に変換
y_data = df.loc[:, 'y'].values
print(y_data)
[205.06 179.572 20.035 168.983 127.056 99.45 22.804 178.753 93.451 27.555]
ある程度の準備ができたので，my_func() 関数がどのような動作をするのか確認しておこう．まず，パラメータベクトル \(\bf{w}\) に [3, 4, 5] という値を与えて， \(\bf X\) の先頭行のデータに対する \(y\) の値を計算してみよう．
y の計算例
w = np.array([3, 4, 5])
my_func(w, x_data[0])
140.96800000000002
上の例では次の計算を行ったことになります．ただし，\(x_0 = 1\) です．
また，my_func() に行列 \(\bf X\) を与えれば，\(\bf X\) のすべての行について上と同じ計算をまとめて実行できます．
y の計算例
my_func(w, x_data)
array([140.968, 121.744, 8.744, 112.763, 85.42 , 65.532, 12.202,
122.075, 61.062, 13.855])
次に，関数 get_ridge() にデータを与えて目的関数値を計算してみよう．このとき，目的関数の定義において，alpha のデフォルト値が 1.0 と定義されているので，引数の alpha を省略することが可能です．
目的関数値の計算例
w = np.array([3, 4, 5])
get_ridge(w, x_data, y_data)
18235.559757999996
引数 alpha に値を指定することも可能です．
alpha を指定した目的関数値の計算例
get_ridge(w, x_data, y_data, 0.1)
18190.559757999996
パラメータベクトル \(\bf{w}\) を [10, 10, 5] に変更して，alpha = 0.1 のままで目的関数値を計算します．
目的関数値の計算例
w = np.array([10, 10, 5])
get_ridge(w, x_data, y_data)
242.58465000000007
この結果，目的関数値が 18189.56 から 242.58 まで小さくなりました．つまり，問題はこの目的関数値が最も小さくなるような重みパラメータの値を求めることになります．
重みパラメータの最適化についてはここで説明したとおりです．ここでも同じ方法で最適化を行います．初期値は [1.0, 1.0, 1.0] としておきます．alpha は 0.1 で固定することにします．最適化のアルゴリズムにはネルダーミード法を利用します．
最適化
w = np.array([1.0, 1.0, 1.0])
alpha = 0.1
results = optimize.minimize(
get_ridge,
w,
args=(x_data, y_data, alpha),
method='Nelder-Mead'
)
print(results)
final_simplex: (array([[11.00937185, 5.73576888, 7.0439693 ],
[11.00938021, 5.73582026, 7.04394338],
[11.00938364, 5.73583333, 7.04393587],
[11.0094044 , 5.73571055, 7.04399628]]), array([33.68924498, 33.68924498, 33.68924498, 33.68924498]))
fun: 33.689244982421954
message: 'Optimization terminated successfully.'
nfev: 372
nit: 202
status: 0
success: True
x: array([11.00937185, 5.73576888, 7.0439693 ])
この結果，目的関数の最小値は 33.69 となり，重みパラメータの最適解は [11.0, 5.7, 7.0] となりました．重回帰分析で重みパラメータの推定値が大きくなってしまった問題が解決されました．
なお，前のページで説明したとおり，重回帰分析の解とリッジ回帰の解はそれぞれ
となるので，これらの式を用いても同じ結果が得られることを確認しよう．
まず，x_data に行列のデータが格納されているので，これを X に代入してから，行列 \(\bf X\) の転置行列 \({\bf X}^{\rm T}\) を transpose() で求める．
行列の転置
X = x_data
X.transpose()
array([[ 1. , 1. , 1. , 1. , 1. , 1. , 1. , 1. ,
1. , 1. ],
[ 9.852, 8.516, 0.416, 7.852, 5.855, 4.463, 0.623, 8.49 ,
4.143, 0.8 ],
[19.712, 16.936, 0.816, 15.671, 11.8 , 8.936, 1.342, 17.023,
8.298, 1.531]])
なお，行列 \(\bf X\) の転置行列 \({\bf X}^{\rm T}\) は X.T でも求めることができます．このほうが表記が短くなるので以降では X.T を使うことにします．
行列の転置
X.T
array([[ 1. , 1. , 1. , 1. , 1. , 1. , 1. , 1. ,
1. , 1. ],
[ 9.852, 8.516, 0.416, 7.852, 5.855, 4.463, 0.623, 8.49 ,
4.143, 0.8 ],
[19.712, 16.936, 0.816, 15.671, 11.8 , 8.936, 1.342, 17.023,
8.298, 1.531]])
行列の積 \({\bf X}^{\rm T} {\bf X}\) は np.dot() で求めることができます．
行列の積
np.dot(X.T, X)
array([[ 10. , 51.01 , 102.065 ],
[ 51.01 , 375.883192, 751.752866],
[ 102.065 , 751.752866, 1503.513491]])
逆行列は np.linalg.inv() で求めることができるので，\(\left( {\bf X}^{\rm T} {\bf X} \right)^{-1} \) は次のようになります．
逆行列
np.linalg.inv(np.dot(X.T, X))
array([[ 0.33172528, 0.84799139, -0.44651246],
[ 0.84799139, 117.12690377, -58.62071499],
[ -0.44651246, -58.62071499, 29.34118254]])
上の結果に転置行列をかけて \(\left( {\bf X}^{\rm T} {\bf X} \right)^{-1} {\bf X}^{\rm T} \) を求めると次のようになります．
np.dot(np.linalg.inv(np.dot(X.T, X)), X.T)
array([[-0.11551704, -0.00891496, 0.32013554, -0.00714299, 0.02786791,
0.12627557, 0.26080421, -0.06980932, 0.13979327, 0.32650783],
[-0.74928651, 5.50027486, 1.73827993, 1.88321522, -5.09842389,
-0.24934621, -4.85094707, -2.64502684, -0.32993926, 4.80119976],
[ 0.39559374, -2.73825377, -0.89032494, -0.93069493, 2.55515528,
0.12204374, 2.40864908, 1.3385677 , 0.16099808, -2.42173398]])
最後にベクトル \(\bf y\) をかけることで重回帰分析の解 \({\bf w} = \left( {\bf X}^{\rm T} {\bf X} \right)^{-1} {\bf X}^{\rm T} {\bf y}\) を求めることができます．
重回帰分析の解
np.dot(np.dot(np.linalg.inv(np.dot(X.T, X)), X.T), y_data)
array([11.54645642, 32.55806264, -6.40309939])
上の結果はここで求めた解と一致しました．
さらに，リッジ回帰の解
も同じように求めてみよう．
まず，単位行列は np.eye() で定義できるので，行数および列数が未知パラメータ数 (3) と一致するように単位行列 \(\bf I \)を生成します．
単位行列
I = np.eye(3)
I
array([[1., 0., 0.],
[0., 1., 0.],
[0., 0., 1.]])
次に，alpha に 0.1 を代入してから \( \left( {\bf X}^{\rm T} {\bf X} + \alpha {\bf I} \right) \) を求めると次のようになります．
行列の和
alpha = 0.1
np.dot(X.T, X) + alpha * I
array([[ 10.1 , 51.01 , 102.065 ],
[ 51.01 , 375.983192, 751.752866],
[ 102.065 , 751.752866, 1503.613491]])
上の式の逆行列 \( \left( {\bf X}^{\rm T} {\bf X} + \alpha {\bf I} \right)^{-1}\) を求める．
逆行列
np.linalg.inv(np.dot(X.T, X) + alpha * I)
array([[ 0.31553789, 0.04447059, -0.04365235],
[ 0.04447059, 7.485877 , -3.74568891],
[-0.04365235, -3.74568891, 1.87633841]])
更に，転置行列 \( {\bf X}^{\rm T} \) をかけると， \( \left( {\bf X}^{\rm T} {\bf X} + \alpha {\bf I} \right)^{-1} {\bf X}^{\rm T}\) は次のようになります．
転置行列をかける
np.dot(np.linalg.inv(np.dot(X.T, X) + alpha * I), X.T)
array([[-0.10681306, -0.04504683, 0.29841734, -0.01935508, 0.06081542,
0.1239327 , 0.28466161, -0.05000082, 0.13755231, 0.28428261],
[-0.03968902, 0.35721172, 0.10211327, 0.12488587, -0.32484873,
-0.01753647, -0.31854256, -0.16329602, -0.02326758, 0.29852247],
[ 0.04020329, -0.16427176, -0.0707668 , -0.05070241, 0.16613235,
0.0062981 , 0.1408296 , 0.0963576 , 0.00781464, -0.16752937]])
最後にベクトル \(\bf y\) をかけることでリッジ回帰の解 \({\bf w} = \left( {\bf X}^{\rm T} {\bf X} + \alpha {\bf I} \right)^{-1} {\bf X}^{\rm T} {\bf y}\) を求めることができます．
リッジ回帰の解
np.dot(np.dot(np.linalg.inv(np.dot(X.T, X) + alpha * I), X.T), y_data)
array([11.00938331, 5.7357478 , 7.04397916])
上の解はやはり最適化によって求めた解と一致しました．
このページではリッジ回帰のコードを記述してその内容を確認しました．Python の scikit-learn を使えば短いコードで簡単にリッジ回帰も行えるので，次のページではその方法を確認します．