以下のように、入力データ\(x\),\(y\)の値が与えられた場合を考える。
このときの、各点との距離が最小になるような\(y=ax+b\)という直線(=回帰直線)の傾き\(a\), 切片\(b\)を求めていく方法として、最小2乗法がある。
最小2乗法とは、回帰直線と各点との距離の2乗和が最小になるような関係式を求める方法で、以下の計算結果を最小にする\(a\),\(b\)を求めればよい。
\[
\begin{eqnarray}
S &=& (ax_{1}+b-y_{1})^2 + (ax_{2}+b-y_{2})^2 + ･･･ + (ax_{n}+b-y_{n})^2 \\
&=& \sum_{i=1}^{n} (ax_{i}+b-y_{i})^2
\end{eqnarray}
\]
上記\(S\)の値を最小にする\(a\),\(b\)を求めるには、最急降下法で、以下の更新式を繰り返すことで求めることができる。ただし、今回は\(a\),\(b\)それぞれについて偏微分した更新式が必要になる。
\(S\)を\(a\),\(b\)それぞれについて偏微分した式は以下のようになるため、この式を利用することになる。
\[
\begin{eqnarray}
\displaystyle \frac{\partial S}{\partial a} &=& 2\sum_{i=1}^{n} (ax_{i}+b-y_{i})x_{i} \\
\displaystyle \frac{\partial S}{\partial b} &=& 2\sum_{i=1}^{n} (ax_{i}+b-y_{i})
\end{eqnarray}
\]
なお、上記式の導出過程は、以下の記事を参照のこと。
最小2乗法の計算式を偏微分してみた このブログの以下の記事で、最小2乗法と最急降下法を用いて回帰直線を求めている。 https://www.purin-it.co...
今回は、最小2乗法と最急降下法を用いて、各点との距離が最小になるような\(y=ax+b\)という直線(=回帰直線)の傾き\(a\), 切片\(b\)を求めてみたので、そのサンプルプログラムを共有する。
前提条件
下記記事のAnacondaをインストールしJupyter Notebookを利用できること
Python開発用のAnacondaをインストールしJupyter Notebookを利用してみた今回は、Pythonを勉強してみたいと思い、Python開発環境を構築してみたので、その手順を共有する。 Python開発用として...
入力データx,yの値(全20個)を読み込み、散布図として表示すると、以下のようになる。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 | %matplotlib inline import numpy as np import matplotlib.pyplot as plt # 入力データの読み込み input_data = np.array([[33,352], [33,324], [34,338], [34,317], [35,341], [35,360], [34,339], [32,329], [28,283], [35,372], [33,342], [28,262], [32,328], [33,326], [35,354], [30,294], [29,275], [32,336], [34,354], [35,368]]) # x座標、y座標の抜き出し input_data_x = input_data[:, 0] input_data_y = input_data[:, 1] # 入力データを散布図で表示 plt.scatter(input_data_x, input_data_y) plt.title("input_data") plt.xlabel("x", size=14) plt.ylabel("y", size=14) plt.xlim(26, 36) plt.ylim(0, 450) plt.grid() plt.show() |
先ほど読み込んだ入力データの内容を出力した結果は、以下の通り。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 | import numpy as np # 入力データの読み込み input_data = np.array([[33,352], [33,324], [34,338], [34,317], [35,341], [35,360], [34,339], [32,329], [28,283], [35,372], [33,342], [28,262], [32,328], [33,326], [35,354], [30,294], [29,275], [32,336], [34,354], [35,368]]) # 入力データを出力 print("*** 入力データ ***") print(input_data) print("*** 入力データの形状 ***") print(input_data.shape) # x座標、y座標の抜き出し input_data_x = input_data[:, 0] input_data_y = input_data[:, 1] # x座標、y座標を出力 print("*** x座標 ***") print(input_data_x) print("*** y座標 ***") print(input_data_y) |
入力データを読み込み、\(a\),\(b\)に特定の値を設定した後で、\(\displaystyle \frac{\partial S}{\partial a} = 2\sum_{i=1}^{n} (ax_{i}+b-y_{i})x_{i}\)を(2*入力データ数)で除算した結果は、以下の通り。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 | %matplotlib inline import matplotlib.pyplot as plt def da_f(a, b, input_data): ret = 0 input_data_cnt = input_data.shape[0] print("*** da_f(a, b, input_data)の呼び出し ***") print("input_data_cnt=" + str(input_data_cnt)) for tmp in range(input_data_cnt): tmp_x = input_data[tmp, 0] tmp_y = input_data[tmp, 1] print("tmp=" + str(tmp) + ", tmp_x=" + str(tmp_x) + ", tmp_y=" + str(tmp_y)) ret = ret + (( a * tmp_x + b - tmp_y ) * tmp_x) / input_data_cnt return ret # 入力データの読み込み input_data = np.array([[33,352], [33,324], [34,338], [34,317], [35,341], [35,360], [34,339], [32,329], [28,283], [35,372], [33,342], [28,262], [32,328], [33,326], [35,354], [30,294], [29,275], [32,336], [34,354], [35,368]]) # a,bの初期値 a = 10 b = 10 # 関数(da_f(a, b, input_data))を呼び出した場合の結果を確認 da_f(a, b, input_data) |
なお、\(S\)を\(a\)で偏微分した値を(2*入力データ数)で除算しているのは、計算結果をなるべく小さくするためである。また、\(S\)を\(b\)で偏微分した値も似たような式で計算できる。
最小2乗法と最急降下法を用いて、各点との距離が最小になるような\(y=ax+b\)という直線(=回帰直線)の傾き\(a\), 切片\(b\)を求めた後で、その結果をグラフで描画した結果は、以下の通り。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 | %matplotlib inline import numpy as np import matplotlib.pyplot as plt def da_f(a, b, input_data): ret = 0 input_data_cnt = input_data.shape[0] for tmp in range(input_data_cnt): tmp_x = input_data[tmp, 0] tmp_y = input_data[tmp, 1] ret = ret + (( a * tmp_x + b - tmp_y ) * tmp_x) / input_data_cnt return ret def db_f(a, b, input_data): ret = 0 input_data_cnt = input_data.shape[0] for tmp in range(input_data_cnt): tmp_x = input_data[tmp, 0] tmp_y = input_data[tmp, 1] ret = ret + ( a * tmp_x + b - tmp_y ) / input_data_cnt return ret # 入力データの読み込み input_data = np.array([[33,352], [33,324], [34,338], [34,317], [35,341], [35,360], [34,339], [32,329], [28,283], [35,372], [33,342], [28,262], [32,328], [33,326], [35,354], [30,294], [29,275], [32,336], [34,354], [35,368]]) # x座標、y座標の抜き出し input_data_x = input_data[:, 0] input_data_y = input_data[:, 1] # a,bの初期値 a = 10 b = 10 # 学習率η eta = 0.0001 # 最急降下法を1,000回分繰り返した場合を確認 for num in range(1000): a = a - eta * da_f(a, b, input_data) b = b - eta * db_f(a, b, input_data) # 最急降下法で計算後の値を確認 print("*** 最急降下法で計算後の各値 ***") print("*** a,bの値 ***") print("a = " + str(a) + ", b = " + str(b)) print("*** da_f(a, b, input_data)の値 ***") print(da_f(a, b, input_data)) print("*** db_f(a, b, input_data)の値 ***") print(db_f(a, b, input_data)) # 入力データの値を散布図で表示 plt.scatter(input_data_x, input_data_y) plt.title("input_data") plt.xlabel("x", size=14) plt.ylabel("y", size=14) plt.xlim(26, 36) plt.ylim(0, 450) plt.grid() # 算出した直線(y=ax+b)を追加で表示 x = np.linspace(26, 36, 1000) y = a * x + b plt.plot(x, y, label='y=ax+b', color='darkviolet') plt.legend() plt.show() |
また、最急降下法の、更新式を1,000回繰り返した際の、(Sをaで偏微分した式を書く)(Sをbで偏微分した式を書く)の変化は以下の通りで、約50回の繰り返しで最小値に近づいているのが確認できる。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 | %matplotlib inline import numpy as np import matplotlib.pyplot as plt def da_f(a, b, input_data): ret = 0 input_data_cnt = input_data.shape[0] for tmp in range(input_data_cnt): tmp_x = input_data[tmp, 0] tmp_y = input_data[tmp, 1] ret = ret + (( a * tmp_x + b - tmp_y ) * tmp_x) / input_data_cnt return ret def db_f(a, b, input_data): ret = 0 input_data_cnt = input_data.shape[0] for tmp in range(input_data_cnt): tmp_x = input_data[tmp, 0] tmp_y = input_data[tmp, 1] ret = ret + ( a * tmp_x + b - tmp_y ) / input_data_cnt return ret # 入力データの読み込み input_data = np.array([[33,352], [33,324], [34,338], [34,317], [35,341], [35,360], [34,339], [32,329], [28,283], [35,372], [33,342], [28,262], [32,328], [33,326], [35,354], [30,294], [29,275], [32,336], [34,354], [35,368]]) # x座標、y座標の抜き出し input_data_x = input_data[:, 0] input_data_y = input_data[:, 1] # a,bの初期値 a = 10 b = 10 # 学習率η eta = 0.0001 # da_f, db_fの値の変化を格納 repeat_num = [] da_f_value = [] db_f_value = [] # 最急降下法を1,000回分繰り返した場合を確認 for num in range(1000): da_f_num = da_f(a, b, input_data) db_f_num = db_f(a, b, input_data) a = a - eta * da_f_num b = b - eta * db_f_num repeat_num.append(num) da_f_value.append(da_f_num) db_f_value.append(db_f_num) # 最急降下法で計算後の値を確認 print("*** 最急降下法で計算後の各値 ***") print("*** a,bの値 ***") print("a = " + str(a) + ", b = " + str(b)) print("*** da_f(a, b, input_data)の値 ***") print(da_f(a, b, input_data)) print("*** db_f(a, b, input_data)の値 ***") print(db_f(a, b, input_data)) # da_fの値の変化(repeat_num,da_f_value)の値の変化をグラフ化 plt.plot(repeat_num, da_f_value) plt.title("da_f value change") plt.xlabel("repeat_num", size=14) plt.ylabel("da_f_value", size=14) plt.grid() plt.show() # db_fの値の変化(repeat_num,da_f_value)の値の変化をグラフ化 plt.plot(repeat_num, db_f_value) plt.title("db_f value change") plt.xlabel("repeat_num", size=14) plt.ylabel("db_f_value", size=14) plt.grid() plt.show() |
要点まとめ
- 各点との距離が最小になるようなy=ax+bという直線(=回帰直線)の傾きa, 切片bを求めていく方法として、最小2乗法がある。
- 最小2乗法とは、回帰直線と各点との距離の2乗和が最小になるような関係式を求める方法をいい、この計算式を最小にするa,bの値は、最急降下法で求められる。
