最急降下法による単回帰分析

はじめに

データ分析を学ぶためにcouseraのコースを受講しています。

その中で最小二乗法、最急降下法を用いた単回帰分析を受講し終わりました。

  • 更に理解を深める目的
  • 同じように学習している人の理解の助けになれば嬉しい

上記の理由から実際にpythonによる実装を行い解説を記載しました。

参考になれば幸いです。

対象読者

  • 単回帰分析について仕組みを知りたい人
  • 単回帰分析についておおよそわかったのでプログラムでの実現方法を知りたい人

ざっくり概要

  • 単回帰分析の解説
  • 回帰式の求め方の紹介
  • pythonでの動作確認

単回帰分析とは

「既知の情報から知りたい情報(情報)を予測する」ための手段の1つです。

例えば家の価値を例に考えます。

前提

  • リビングの広さ(sqft)と家の価値のデータ(既知の情報)が手元にある。

この場合価値が不明の家の価値(知りたい情報)を予測するために回帰分析が活用できます。

仮にリビングの広さが1,1000の家の価値を回帰分析で求めると下記になります。

f:id:campusista-takahiro:20200524064635p:plain
図1. リビングの広さが11000の場合の家の価値(赤点)

図1中の青線で表されている1次関数を回帰式と呼びます。回帰分析とはこの回帰式を求めることを言います。

回帰式がどうして予測につながるのでしょうか。

図1の青線の式は下記です。(小数点4桁以降切り捨て)

$$ f(x) = 263.089x -0.101 $$

この式中の「x」と「f(x)」は下記に読み替えることができます

$$家の価値 = 263.089*家の広さ$$

こうするとわかってくるのではないでしょうか。家の広さに11,000を入力し、計算して得られた値は家の価値です。既知の情報には家の広さが11,000のときの家の価値はありませんでした。

つまり回帰式を利用すれば指定した家の広さの家の価値を予測することができているというわけです。

このように回帰分析を活用することで既知の情報から未知の情報を予測できることが可能になります。

この手法が役に立つ場合は「厳密に求めるとコストが掛かってしまうようなことを予測したい」場合に役に立ちます。

例に出した家の価値を求める場合では以下のようなコストが発生することが考えられます。

  • 実際に家を訪れ人が確認しなければならない(人件費)

人ではなくロボットかもしれませんがいずれにしろコストが掛かっていることがわかると思います。

そこで回帰分析を活用して家の価値を予測するとどうでしょうか。

求めた回帰式に調べたい家の広さを入力するだけでおおよそ尤もらしい値が返ってきます。

回帰分析を求めるのとそうでないのとでは労力が雲泥の差です。このような場合に回帰分析は価値を発揮すると言えます。

まとめると

  • 単回帰分析は「既知の情報から知りたい情報(情報)を予測する」手段の1つ
  • 既知の情報から大きく外れないような回帰式を求めることを回帰分析と呼ぶ
  • その情報を求めるときのコストが大きい場合に役に立つ

回帰式の求めを考える

回帰式をどのように求めるか考えてみましょう。仮に人間が目で見て予測する場合を考えてみます。

図1から、リビングの広さ(以降広さ)と価値は相関していることが伺えます。

同様に広さが11000の家の価値はどれほどかグラフの見た目から予測してみましょう。

自分は赤点のあたりだと予想しました。(図2

f:id:campusista-takahiro:20200524064535p:plain
図2. グラフの見た目から広さ11000の場合の値段を予想した図

図2を見ている人の殆どは納得できるのではないかと思います。

図2の赤点を選択する時何を意識していたでしょうか。恐らく大多数の方は下記のようなことを意識したのではないかと思います

  • 既知の情報から大きく外れないようにする

この「大きく外れないようにする」が回帰式を求める上で重要な考えになってきます。

シンプルな例で回帰式を求めることを考えてみましょう。

f:id:campusista-takahiro:20200524064646p:plain
図3. シンプルデータでの回帰式予想

図3のような情報が既知で回帰式を3通り予想してみました。どれが予想として尤もらしい回帰式でしょうか。正解は青線です。その理由は「最も既知のデータから外れていない」からです。

「外れている」を数学的に表現してみましょう。

赤、青、緑はそれぞれ以下の1次関数で表せます。

$$赤線(x)=1.5x - 0.5$$

$$青線(x)=x$$

$$緑線(x)=0.5x + 0.5$$

それぞれの式で求めた家の価値と既知のデータとの誤差を比較します。誤差は以下の式で求めます。

$$誤差 = \Sigma{i=1}^{n}(既知の情報i - 回帰式(x_i))2$$

自分自身数学が不得意です。細かいところが間違っていると思いますがご容赦ください。また、類似記事ではnで割っていたりしていますが、数値は異なるものの誤差という意味自体は同じです。

例えば上記の数式で i=3 の赤線の誤差を求めると下記のようになります

$$赤線の誤差(3) = (3-(1.5*3-0.5)) ^ 2 = 1$$

上記の具合でそれぞれの誤差を求めてみましょう。結果は下記です

  • 赤線の誤差 : 7.5
  • 青線の誤差 : 0
  • 緑線の誤差 : 7.5

結果より、青線が「最も誤差が小さい」ので青線の式が回帰式として最適と言えます。

以上からどのようにして回帰式を求めると良いでしょうか。

赤・青・緑の式を見比べると係数と切片に差があることが見て取れます。シンプルな例では係数1で切片が0だと誤差が最小でした。

つまり、式の「係数」と「切片」をああでもないこうでもないと調整し誤差が最も小さくなる「係数」と「切片」が求まれば回帰式が求められたと言えそうです。

まとめると

  • 誤差が最小になる「係数」と「切片」を求めれば回帰式を求めたと言える

でした。

回帰式を求める

この作業は人力で行うと途方もく時間がかかりそうであることが想像できるかと思います。

ここからはPCの計算能力を利用して係数と切片を求める方法を理解したいと思います。

今までの解説で回帰式を求めるために「誤差が最小になる係数と切片を求める」事が必要であることがわかりました。ここでは関数の結果を最小にするようなパラメータを求める「最急降下法」について紹介していきます。

例えば下記関数について結果が最小となるパラメータを求めるとします。

$$f(x) = 10x2+3$$

上記式でのパラメータは x です。これをグラフで表すと図のようになります。

f:id:campusista-takahiro:20200524064650p:plain
図4. 式のグラフ

グラフを見ると関数の値が最小になるパラメータxが0だとわかります。

これを最急降下法で求めると以下のような式で求めることができます。

$$x = x - a*{df(x)\over dx}$$

$${df(x)\over dx} = 20x$$

これはつまりどういうことかというと、右辺の計算結果でxを更新し続けるということになります。実際に更新し続けてみましょう。aの意味は後ほど説明します。

loop_count = 10000
a = 0.01

def x(x):
  return 20*x

result_x = 11
for i in range(loop_count):
  result_x = result_x - a * x(result_x)

print(result_x)

# 1e-323

上記の式をpythonで表してみました。結果は0.000...1(小数点以下0が323個)と限りなく0に近づいています。上記プログラムから理解できる情報はだいたい下記です。

  • 10000回繰り返しxの更新を行う
  • ステップ幅は 0.01 (ステップ幅については後ほど説明します
  • xの初期値は11

どういう経緯で0に近づいていっているのか、わかりやすい図を示そうと思います。

f:id:campusista-takahiro:20200524064657p:plain
図5. 最急降下法の解説

図5は下記条件で各ステップで求めた傾きの接線を描画したものです

  • 5回繰り返しxの更新を行う
  • ステップ幅は0.02
  • xの初期値は10
  • 赤・青・緑・黄・茶の順番

ループ1回目の計算を詳しく考えてみます。最急降下法の式に当てはめると下記のようになります。

$$ x_1 = 10 - 0.02 * 20 * 10 = 6 $$

これが2回、3回になると図5のようにだんだん0に近づいていくわけです。このように最小値に近づいていくことを「収束する」といいます。

ステップ幅が大きいとどうなるのでしょうか。

f:id:campusista-takahiro:20200524064701p:plain
図6. ステップ幅を大きくした場合の図

図6はステップ幅を0.08にした場合の図です。今度は左に行ったり右に行ったり忙しそうです。

以上から、ステップ幅は「最小値へ近づく際の1計算あたりの距離の大きさ」といえます。今度は初期値を1、ステップ幅を0.13にしてみましょう。すると以下のようになります

f:id:campusista-takahiro:20200524064704p:plain
図7. ステップ幅を大きすぎる値にした場合

今度は逆に最小値からどんどん離れていってしまっています。これはステップ幅が大きすぎるために起こります。このように値がどんどん大きくなってしまうことを「発散する」と言います。

以上より、最急降下法の仕組みと特徴を知っていただけたかと思います。

今までわかっていることを整理すると以下です。

  • 最急降下法を用いて関数の最小値を求めることができる
  • 「係数」と「切片」の最小値が求まれば回帰式が求められる

つまり、2点目の最小値を最急降下法を用いて求めることができるというわけです。以上理解できたことを駆使して図1データに最もフィットする1次関数を求めて行きたいと思います。

まとめると

  • 最急降下法を用いて誤差が最小となる係数と切片を求めることができる
  • ステップ幅は最小値へ近づく際の1計算あたりの大きさ
  • ステップ幅が大きすぎれば「発散」する場合があり、小さすぎれば「収束」するまでの計算回数が膨大になってしまう。

家の価値データセットで単回帰分析を行う

今までに理解した内容から最初に紹介したデータを元にリビングの広さ(以降広さ)から家の価値を求めるために単回帰分析を行ってみます。今回求めたい回帰式は以下のようになります。

$$家の価値 = 係数 * 広さ + 切片$$

最急降下法を用いて誤差が最小となるような「係数」、「切片」を求めていきます。式を整理していきましょう。

$$家の価値 = price   $$

$$家の広さ = size   $$

$$回帰式 = exp   $$

$$係数 = coefficient  $$

$$切片 = intercept  $$

$$誤差 = \sum_{i=1}^{n} (price_i + exp(coefficient,intercept)) ^ 2   $$

$$x = x - a*{df(x)\over dx}$$

今回最小化させたいのは「誤差」なので最急降下法の式中 f(x) は誤差を求める関数が入ります。それを踏まえ導き出せる式は以下の2つです。

$$coefficient = coefficient - a*{d\over dcoefficient} \sum_{i=1}^{n} (price_i - exp(coefficient, intercept)) ^ 2    $$

$$intercept = intercept - a*{d\over dintercept} \sum_{i=1}^{n} (price_i - exp(coefficient, intercept)) ^ 2    $$

求めたい値が2つなので式も2つになり、それぞれの値で偏微分した式となります。そこから微分計算を行った結果が以下の式です。

$$coefficient = coefficient - 2 * a * \sum_{i=1}^{n} (price_i-(coefficient*size_i + intercept)) * size_i$$

$$intercept = intercept - 2 * a * \sum{i=1}^{n} (price_i-(coefficient*size_i + intercept))$$

ここまで来たら後はプログラムに表現して実行するだけです。早速下記プログラムを実行してみましょう。

import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline

csv = pd.read_csv('/path/to/home_data.csv', usecols=['price', 'sqft_living'])

size_of_house = csv['sqft_living']
value_of_house = csv['price']
step_width = 0.000000000001
itera = 1000
errors = np.zeros(itera)
intercept = 0.1
coefficient = 1

for i in range(itera):
    errors[i] = 2 * np.sum(np.square(value_of_house - (intercept + coefficient * size_of_house)))
    temp0 = intercept - 2 * step_width * np.sum((intercept+coefficient * size_of_house) - value_of_house)
    temp1 = coefficient - 2 * step_width * np.sum(((intercept+coefficient * size_of_house) - value_of_house) * size_of_house)
    
    intercept = temp0
    coefficient = temp1

print(intercept)
print(coefficient)
plt.plot(errors)
plt.xlabel("Iteration")
plt.ylabel("Cost")

plt.show()

上記のポイントは temp0, temp1 の部分です。求めた式を直感的にプログラムに実装すると

下記のように考える方もいらっしゃると思います。

intercept = intercept - 2 * step_width * np.sum((intercept+coefficient * size_of_house) - value_of_house)

が、これだとそこからcoefficientを求める際のinterceptは更新後の値となってしまうため都合が悪いのです。そのために別の変数を用意しています。

プログラム実行結果は下記となっています。

f:id:campusista-takahiro:20200524064707p:plain
図8. プログラム実行結果

これは横軸が「繰り返し回数」縦軸が「誤差」を表しており、計算を重ねる毎に誤差が小さくなっていることが伺えます。また、上の数値はそれぞれ切片・係数を表しており、最初のあたりで記載した値が求められていることがわかります。

以上が単回帰分析についての解説でした。おかしな点などありましたら指摘いただけるととても嬉しいです。最後まで読んでいただきありがとうございました。

参考資料