勾配降下法で重回帰分析してみた
Pythonで勾配降下法では、単回帰(回帰直線)で実験しました。
今回は、単回帰含め、重回帰分析まで広げてみます。
使うデータセットは、skleanのボストン住宅価格です。
モデルの説明
計画行列
$$
\boldsymbol { X } = \left( \begin{array} { c c c c c } { 1 } & { a _ { 11 } } & { a _ { 21 } } & { \dots } & { a _ { n 1 } } \\ { 1 } & { a _ { 12 } } & { a _ { 22 } } & { \dots } & { a _ { n 2 } } \\ { \vdots } & { \vdots } & { \vdots } & { \ddots } & { \vdots } \\ { 1 } & { a _ { 1 m } } & { a _ { 2 m } } & { \dots } & { a _ { n m } } \end{array} \right)
$$
被説明変数
$$
\boldsymbol { y } = \left( \begin{array} { c } { y _ { 1 } } \\ { y _ { 2 } } \\ { \vdots } \\ { y _ { m } } \end{array} \right)
$$
パラメータ
$$
\boldsymbol { \omega } = \left( \begin{array} { c } { \omega _ { 0 } } \\ { \omega _ { 1 } } \\ { \vdots } \\ { \omega _ { m } } \end{array} \right)
$$
回帰モデル
$$
f = X \omega
$$
勾配法による更新
評価関数
$$
E \left( \omega _ { j } \right) = \frac { 1 } { 2 } \sum _ { i = 1 } ^ { n } \left( f _ { i } – y _ { i } \right) ^ { 2 }
$$
パラメータ更新
$$
\omega _ { j } = \omega _ { j } – \alpha \frac { \partial E } { \partial \omega _ { j } }
$$
このとき微分を予め計算しておくと、数値微分をしなくても済む。
合成微分(チェーンルールとかなんちゃら)
$$
\begin{aligned} \frac { \partial E \left( \omega _ { j } \right) } { \partial \omega _ { j } } & = \frac { \partial E \left( \omega _ { j } \right) } { \partial f } \frac { \partial f } { \partial \omega _ { j } } \\ & = \sum \left\{ \left( f _ { i } – y _ { i } \right) x _ { j i } \right\}\end{aligned}
$$
\(x_{0i}\)は、1であることに注意が必要。
iはデータのインデックスである。
練習:回帰直線で検証
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
np.random.seed(1)
#学習率
alpha = 1e-3
#数値微分のh
h = 1e-3
x = np.arange(10).reshape(-1, 1) +1
seg = np.ones(10).reshape(-1,1)
X = np.c_[seg, x]
beta = np.array([[5], [3]])
y = np.dot(X, beta)
print("モデルのパラメータ:\n", beta)
#パラメータの初期値
beta = np.random.normal(0, 1, size = (X.shape[1], 1))
print("パラメータの初期値:\n", beta)
def E(x, y, beta):
return 0.5 * ((y - np.dot(x, beta))**2).sum()
def diff(X, y, beta, n, h=1e-3):
beta_a = beta.copy()
beta_b = beta.copy()
beta_a[n] = beta_a[n] + h
beta_b[n] =beta_b[n] - h
a = E(X, y, beta_a)
b = E(X, y, beta_b)
return (a-b)/(2*h)
grad = 10
while abs(grad) > 0.01:
grad = 0
for i in range(beta.shape[0]):
beta[i] = beta[i] - alpha * diff(X, y, beta, i)
grad = grad + diff(X, y, beta, i)
print("勾配法のパラメータ推定:\n", beta)
''' 計算結果 ''' モデルのパラメータ: [[5] [3]] パラメータの初期値: [[ 1.62434536] [-0.61175641]] 勾配法のパラメータ推定: [[4.99478088] [3.00074812]]
本番:重回帰分析
パラメータの推定がうまくいかなかった。
今回、パラメータの推定を10回行い、その推定値をparamリストに保持するようにしたが、１つに推定値が定まることはなかった。
失敗しているのか…
多重共線性も関係してくるのか…
多重共線性とVIF統計量の求め方
検討しなければ…
sklearnや、最小二乗法の行列計算による重回帰分析は下のリンクで実行しています。
重回帰分析をPythonで実装する
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
from sklearn.datasets import load_boston
# サンプルデータを用意
dataset = load_boston()
# 説明変数データを取得
X = pd.DataFrame(dataset.data, columns=dataset.feature_names)
X = X.values
#被説明変数データを取得
y = pd.DataFrame(dataset.target, columns=['target'])
y = y.values
seg = np.ones((X.shape[0], 1))
X = np.c_[seg, X]
#パラメータ初期値
beta = np.random.normal(0, 1, size = (X.shape[1], 1))
def E(x, y, beta):
return 0.5 * ((y - np.dot(x, beta))**2).sum()
def diff(X, y, beta, n, h=1e-3):
beta_a = beta.copy()
beta_b = beta.copy()
beta_a[n] = beta_a[n] + h
beta_b[n] =beta_b[n] - h
a = E(X, y, beta_a)
b = E(X, y, beta_b)
return (a-b)/(2*h)
alpha = 1e-5
'''
for i in range(10000):
grad = 0
for i in range(beta.shape[0]):
beta[i] = beta[i] - alpha * diff(X, y, beta, i)
grad = grad + diff(X, y, beta, i)
print(grad)
beta
'''
#推定を10回行う
param = []
for i in range(10):
#パラメータ初期値
beta = np.random.normal(0, 1, size = (X.shape[1], 1))
grad = 10
while abs(grad) > 0.001:
grad = 0
for i in range(beta.shape[0]):
beta[i] = beta[i] - alpha * diff(X, y, beta, i)
grad = grad + diff(X, y, beta, i)
param.append(beta)
beta