ブログ|株式会社セシオス

機械学習入門(回帰モデルその2)

作成者: 株式会社セシオス|Apr 8, 2021 3:20:22 AM
  1.  
  2.  
  3. 機械学習入門(線形回帰)では一つの変数から数値を予測しましたが、今度は二つの変数から数値を予測してみましょう。
  4. 一つの変数から数値を予測するものを単回帰分析、二つ以上の変数から数値を予測するものを重回帰分析と呼んだりします。

 

 

重回帰分析

 

1. 前回のやり方


まず前回との違いが分かるように同じ手順でやってみましょう。
データは前回の身長・体重に加え腹囲を追加して、身長・腹囲から体重を予測します。
身長(x1) 腹囲(x2) 体重(y)
172.5 80.3 70.2
171 74.2 62.1
165.5 74.2 58.1
169.9 75.6 60.4
170.2 78.2 64.1

 

 

データ生成

D = np.array(((172.5,80.3,70.2),(171,74.2,62.1),(165.5,74.2,58.1),(169.9,75.6,60.4),(170.2,78.2,64.1)))
d_avg = np.average(D, axis=0)
x = D - d_avg

データの生成は前回とあまり変わりありません。
 
予測関数

yp = x1 * w1 + x2 * w2

"x2 + w2"が追加されただけです。

 

損失関数 

def lose_score(X, w):
    l = X[:,0] * w[0] + X[:,1] * w[1] - X[:,2]
    lose = np.sum(l**2)
    return lose

損失関数ですが前回と比較して胸囲とそのウェイトが追加されて少し複雑になってます。

 

l = X[:,0] * w[0] + X[:,1] * w[1] - X[:,2]
yd = yp - yt = x1 * w1 + x2 * w2 - yt

前回と同様に代入しただけです。

 

 

勾配関数

def grad(X, w):
    x = X[:,0]
    y = X[:,1]
    z = X[:,2]
    g1 = 2 * x * (x * w[0] + y * w[1] -z)
    grad1 = np.sum(g1)
    g2 = 2 * y * (y * w[1] + x * w[0] -z)
    grad2 = np.sum(g2)
    grad = np.array((grad1, grad2))
    return grad

勾配関数は、変数w1,w2と二つあるのでそれぞれの偏微分の配列が傾きになります。
データの次元数が増える毎にコードが増えてしまいます。

 

g1 = 2 * x * (x * w[0] + y * w[1] -z)
d(l**2) / d(w1) = d((x * w1 + y * w2 - z)**2) / d(w1)
d(l**2) / d(w1) = d(x**2 * w1**2 + y**2 * w2**2 + z**2 + 2 * x * y * w1 * w2 - 2 * y * z * w2 - 2 * z * x * w1) / d(w1)
d(l**2) / d(w1) = x**2 * 2 * w1 + 2 * x * y * w2 - 2 * z * x = 2 * x * (x * w1 + y * w2 - z)

g2 = 2 * y * (y * w[1] + x * w[0] -z)
d(l**2) / d(w2) = d((x * w1 + y * w2 - z)**2) / d(w2)
d(l**2) / d(w2) = d(x**2 * w1**2 + y**2 * w2**2 + z**2 + 2 * x * y * w1 * w2 - 2 * y * z * w2 - 2 * z * x * w1) / d(w2)
d(l**2) / d(w2) = y**2 * 2 * w2 + 2 * x * y * w1 - 2 * y * z = 2 * y * (y * w2 + x * w1 - z)

g1, g2は前回同様l**2をそれぞれw1, w2で偏微分したものです。

 

 

学習関数

def train(X, l_rate=0.00003, max_run=50):
    w = np.ones(2)
    for _ in range(max_run):
        g = grad(X, w)
        if np.linalg.norm(g) == 0:
            break
        else:
            w -= l_rate * g
    return w

データの次元数が増えたためウェイトの次元数が増えているくらいでここはあまり変わりません。

 

できなくはありませんが次元数によってコードが変わってしまったり平方展開がかなり複雑になってますね。データが10次元とかになったらもうやりたくないですね。
そこでより汎用的により簡単に数学的に式を変換して工夫していきます。

 
 

2. 数学的工夫

 

予測・誤差式

 

まず説明しやすいように今回は予測、誤差を先に定義します。

 

yp = x1 * w1 + x2 * w2 = x @ w

 
よくみると(x1, x2)のデータ行列xと(w1, w2)のウェイト行列wの内積になっているので予測式はこのように書きます。

 

yd = yp - yt

これはそのまま誤差配列ですね。

 

 

勾配関数

 

今回は話の流れ上勾配関数から説明します。

d(g・f(w)) / dw = g(x)' * d(f(w)) / dw

 

残差平方和を"yd = yp - yt = w * x - yt" と yd**2 の合成関数とみて、この合成関数の微分法則を使います。

 

d(f(w)) / d(w1) = d(yd) / d(w1)
d(f(w)) / d(w1) = d(yp - yt) / d(w1)
d(f(w)) / d(w1) = d(x1 * w1 + x2 * w2 - yt) / d(w1) = x1
g(yd)' = (yd**2)' = 2 * yd

d(l**2) / d(w1) = d(yd**2) = (yd**2)' * d(yd) / d(w1) = 2 * yd * x1


w1での偏微分

 

d(f(w)) / d(w2) = d(yd) / d(w2)
d(f(w)) / d(w2) = d(yp - yt) / d(w2)
d(f(w)) / d(w2) = d(pred(X, w) - yt) / d(w2)
d(f(w)) / d(w2) = d(x1 * w1 + x2 * w2 - yt) / d(w2) = x2
g(w)' = (yd**2)' = 2 * yd

d(l**2) / d(w2) = d(yd**2) = (yd**2)' * d(yd) / d(w2) = 2 * yd * x2

 

 

w2での偏微分

 

grad(m) = 2 * yd(m) * x(m)

 

ここまで簡単にできました。

損失関数や傾きの値はただの指標として利用するだけで絶対量自体には全く意味がないので微分すると出てくる2で割るのが一般的です。


また差分を2乗したものをデータ数分たしているためデータ数の増減が指標の増減に直接相関してしまうのはよろしくないのでデータ数Mで割るのが一般的です。

 

なので損失関数は一般的に予測値と実数値の残差平方和を2 * Mで割ったものとするのが一般的です。
なんで2 * Mで割るのかと質問されても、単に式や値がきれいになるからとしか言えません。

 

 

grad(m) = yd(m) * x(m) / M

 

 

やっと参考書などでよく見る形になってきました。

ただプログラムでこれを計算するのに、さらに一工夫します。

"yd * x"とは実際にどんな行列でどのように計算すればいいのか考えます。

 

yd * x = yd[0] * x[0] + yd[1] * x[1] + … + yd[4] * x[4]
yd * x = (yd[0] * x1[0] + yd[1] * x1[1] + … + yd[4] * x1[4], yd[0] * x2[0] + yd[1] * x2[1] + … + yd[4] * x2[4])

 

よくみるとパラメータ(身長・腹囲)毎に誤差をかけたものを足した行列になっています。

つまりデータ行列の行と列を入れ替えたもの(転置行列)にydをかけたものになっているのです。

 

grad = x.T @ yd / M

 

一行です。やっとよく公開されてるコードと同じ様なものになりました。

 

 

損失関数

 

勾配関数を説明したときに出たように2 * Mで割ります。

 

lose = np.sum(yd**2) / 2 * M = (np.sum(yd**2) / M) / 2 = np.mean(yd**2) / 2

 

こちらも一行です。よく見る形ですがこういう意味だったのですね。

 

            1. データ生成

            1. コードがn次元のデータにも対応できるように拡張されました。
            2. わざわざ平均で引いてからそのデータで学習していましたが拡張されたことでこれもいらなくなります。
            3.  

              yp = w1 * x + w2 = w1 * x + w2 * 1

               

つまり(w1, w2) の行列w と (x, 1)の行列x の内積ととらえることもできます。

 

yp = w @ x

 

固定値1をいれることで全く同じ様に扱えます。このときの1のほうの変数をダミー変数といいます。

 

d = np.array(((172.5,80.3,1,70.2),(171,74.2,1,62.1),(165.5,74.2,1,58.1),(169.9,75.6,1,60.4),(170.2,78.2,1,64.1)))
x = d[:,:-1]
yt = d[:,-1]

 

 

上記のようにダミー変数をいれることで原点を通るようにデータを変換する措置がいらなくなります。

 

 

全体

 

import numpy as np

d = np.array(((172.5,80.3,1,70.2),(171,74.2,1,62.1),(165.5,74.2,1,58.1),(169.9,75.6,1,60.4),(170.2,78.2,1,64.1)))
x = d[:,:-1]
yt = d[:,-1]

# データ数
M = x.shape[0]
# データ次元数
D = x.shape[1]
w = np.ones(D)
l_rate = 0.01
max_run = 50
loss = 0

for _ in range(max_run):
    yp = x @ w
    yd = yp - yt
    w -= l_rate * x.T @ yd / M
    loss = np.mean(yd**2) / 2

print(w)
print(loss)

 

 

          1. 前回のデータにダミー変数を加えて上記のコードを実行しても結果が同じことに注目してください。
          2. きれいですね。
          3.