機械学習アルゴリズム実装シリーズ [線形回帰編]

機械学習アルゴリズム実装　線形回帰編

こんにちわ。
ブレインズテクノロジー最年少樋口です。

、、、と言っても、もうすぐ27なのですが...orz

今回から何回かに分け、「機械学習アルゴリズム実装シリーズ」と題して機械学習アルゴリズムをPythonで実装をしていきます。もちろん、scikit-learnなどの機械学習ライブラリを使うのではなく重要な部分は全て自分で実装するという試みです。なるべく難しい数式などの説明は最小限にして、機械学習アルゴリズムを少しずつ解説・実装していきたいと思いますので、「これから機械学習いじってみたいけど、何したらいいか分からない。。。」って人の助けになれたら嬉しいです！(と言いつつ自分の備忘録代わりにしているのは秘密ですw)

因みにここで書いている事はCourseraという有名なオンライン学習サイトのMachine Learning by Stanford Universityで学んだ事を元に書いています。Machine Learning以外にもAlgorithmなど豊富なコンテンツが揃っているのでぜひ一度受講してみてください。(基本的に全部無料です!)
www.coursera.org

機械学習とは

ということで、早速始めます。が、その前に簡単に機械学習についてほんのちょっとだけ説明を。

定義としてはいろんな言い方がされていますが、自分としては
What is Machine Learning?
で引用として述べられている一文がしっくりきました。

The field of machine learning is concerned with the question of how to construct computer programs that automatically improve with experience.

要するに、

どうやってコンピュータープログラムに自動的改善を行わせるのかを主題としているのが機械学習という分野である

ということみたいです。プログラムが自動的に自身を改善してくれるとは夢の様なことですねｗ

機械学習には大きく分けて教師あり学習と教師なし学習の２種類がありますが、更に教師あり学習は回帰と分類、教師なし学習はクラスタリングという種類にカテゴリー分けができます。

[教師あり学習]

回帰(Regression)
分類(Classification)

[教師なし学習]

クラスタリング(Clustering)

厳密に言うと半教師学習というものもあったりしますが、ここではとりあえず「まぁ、大きく分けたら3つあんのね。日本語でおｋ」という程度で大丈夫かと。勉強していくうちに詳細は分かってきます。
以下簡単に、種類ごとにまとめてみました。

	教師あり	教師なし
回帰	単線形回帰, 重線形回帰	-
分類	Logistic Regression, SVM, Neural Network	K-means, Mean-shift

ちなみにMATLABで有名なMathWorksさんのページにもここらへんの説明や分かりやすい図が載っているので参照してみてください。
jp.mathworks.com

それぞれの特徴としては、

回帰：連続する値を予測する
分類：離散値を予測する

という違いがあります。具体的例としてそれぞれ以下に例を示しておきます。

[回帰]
過去から現在までの商品の売れ行きから将来の売れ行きを予測
cpuリソースなどの時間毎の使用量予測

[分類]
画像認識(猫なのか、犬なのか、車なのか)など
ニュースカテゴリ分け

他にもいろんなところで機械学習が使われています。
もっと機械学習について学びたい場合は以下、専門の本をどうぞ。実践的なテクニックと共に深く学べます。
machinelearningmastery.com

データを見る

ここまで、機械学習とはなんぞや的な部分を書いてきましたが、ここからは今回実装する線形回帰とはなんぞや的なところに入っていきます。まず以下の例を見て下さい。
f:id:mhigu:20151013151606p:plain
X軸、Y軸に (x, y) = (0, 0) (1, 1) (2, 2) (3, 3) の点がそれぞれプロットされています。では、このデータから $x=4$ の時はどのような点を取りうる事が予想できるでしょうか。

単純に考えれば $y=4$ と予想できると思います。直感的にも分かりやすいですよね。数式で表すなら、高校の時に習った直線を表す式 $y = ax+b$ の $a,b$ を求めて $x=4$ を代入するだけですね。

では次の例です。車のエンジンサイズを横軸(X軸)、それに対する車の値段を縦軸(Y軸)として表示しています。 $x=290$ の時、 $y$ の値はいくつになるでしょうか。
f:id:mhigu:20151026122355p:plain
(データソースは以下参照)
UCI Machine Learning Repository: Automobile Data Set

先ほどの例とは違って単純には求まりそうにはありませんが、 $y=ax+b$ の式に $a$ と $b$ の値を適切に設定すればおおよその値は求まりそうな気もします。図で示すとだいたい以下のような赤で示された直線を求めることが出来れば車の値段 $y$ を求めることができそうです。
f:id:mhigu:20151027114013p:plain
このように、今までに観測されたデータから線形モデルを使ってある値を予測することを線形回帰といいます。

線形回帰理論編

では実際に線形回帰の中身( $a$ と $b$ の決め方)はどうなってんの？ってところの説明に移りたいと思います。まずは先ほどの車のエンジンサイズと値段データの表記を以下の様に定めます。

$x_i$ : 車のエンジンのサイズ (i番目) ※入力値
$y_i$ : 車の値段(i番目)
$m$ : サンプルデータの数(今回はm=205)

このデータに対して予測を行うモデルを以下のような数式で表します。因みにモデルとは先ほどの $y=ax+b$ で表した数式の事を指します。英語ではhypothesisという呼び方をするみたいです。

$h_θ(x)=θ_0+θ_1x$

今回の問題はエンジンサイズ $x=290$ の時の車の値段を求めることでした。ただし、それに先立ちまずは $θ=(θ_0, θ_1)$ の値を適切に設定しなければなりません。どういうアプローチをすればよいでしょうか。例として、 $θ_0=0, θ_1=30$ の時を考えてみます。直線を図示すると
f:id:mhigu:20151026135530p:plain
のようになり、直感的にあまり適切ではないように見えます。
ではどれぐらい適切でないのでしょう。一つの方法として、青で図示されている各点に対する $h_θ(x)$ で求まる値との差分の2乗を計算することで表せます。

[誤差] $=\frac{1}{m}\{(y_1 - 30x_1) + (y_2 - 30x_2) + .... + (y_m - 30x_m)\}^2$
(具体的な値の計算はここでは省きます。)

これを一般化すると,

$j(θ_0, θ_1)=\frac{1}{m}\sum_{k=1}^{m} \{y_k-(θ_0 + θ_1x_k)\}^2$

となります。この $j(θ)$ を目的関数(Cost function)と呼びます。誤差は原則として0に近ければ近いほど良いので(※Over fittingを考慮しない)、目的はこの目的関数を最小にする $θ_0, θ_1$ の値を求める事になります。

いろいろと定義して来たので一度整理しておくと、

問題

　　エンジンの大きさが $x=290$ の時、車の値段 $y$ の値を求める

モデル

　　 $h_θ(x)=θ_0+θ_1x$

パラメータ

　　 $θ_0, θ_1$

目的関数

　　 $J(θ_0, θ_1)=\frac{1}{2m}\sum_{k=1}^{m} \{y_k-(θ_0 + θ_1x_k)\}^2$
　　※後の計算単純化のため全体を2で割っておく

目的

　　 $J(θ_0, θ_1)$ を最小化する $θ_0, θ_1$ を見つける。
以上がここまでのまとめです。

最急降下法

ここからは目的である $J(θ_0, θ_1)$ の最小化をするためのアルゴリズムの説明をしていきます。この部分が教師データ(これまでの観測データ)から実際に自動的改善を行う部分、つまりは機械学習の肝となる部分です。今回は最急降下法(Gradient Descent)というアルゴリズムを使います。これ以外にも、目的関数を最小化するためのアルゴリズムはありますが、シンプルであり線形回帰以外でも使われることのあるアルゴリズムなので覚えておくのも良いかと思います。

最急降下法は以下の手順で処理を進めます。
1. 任意の $θ_0, θ_1$ を設定する
2. $J(θ_0, θ_1)$ が最小となるように $θ_0, θ_1$ を動かす
3. $θ_0, θ_1$ を更新して2.のプロセスへ戻る

これを収束まで繰り返します。

式で表すと、 $θ_0, θ_1$ は以下計算によりにより更新されていきます。
$θ_j-α\frac{∂}{∂θ_j}J(θ_0, θ_1)$ 　※αは学習率

今回の例ではモデルが $h_θ(x)=θ_0+θ_1x$ なので、パラメータは $θ_0, θ_1$ の２つでした。そのため $j=0, 1$ の場合を考えればよく、 $J(θ_0, θ_1)$ をそれぞれで偏微分した結果は以下の様になります(具体的な計算は省略)。

$j=0$ の時
(次の) $θ_0=θ_0-α\sum_{k=1}^{m}\frac{1}{m}\{y_k-(θ_0 + θ_1x_k)\}$
$j=1$ の時
(次の) $θ_1=θ_1-α\sum_{k=1}^{m}\frac{1}{m}\{y_k-(θ_0 + θ_1x_k)\}x_k$

「はて、偏微分てなんぞや」という人は難しいことは考えずに、とりあえず最終的にはこんな式になるのかと思ってもらえればよいです。
これを繰り返し計算していくことで目的関数 $J(θ_0,θ_1)$ を最小とする、 $θ_0, θ_1$ の値が求まります。あとはこれをコードに落として、データを入れてあげれば線形回帰が動くはずです。。。

実装！！！

さてここまでで一通りの線形回帰の処理の流れを解説してきましたので、いよいよいよいよ実装に移りたいと思います。といっても実装自体はここまでに説明してきた数式をコードに落とすだけですので、説明はコードに書いてしまいたいと思います。最急降下法をメソッドの部分は行列計算をしているため少し分かりにくいかもしれませんが、数式とコードを見比べてどう実装されているのか、データがどう処理されているかに注意して追っていけば分かるかと思います。
※ちなみに今回はデータをcsvの形で読み、データクレンジングは事前にしてあります。csvの形式は以下の構造になっています。

エンジンサイズ	車の値段
130	13495
~	~
~	~

[コード]

# -*- coding:utf-8 -*-
import numpy as np


def gradient_descent(x, y, iter_num, alpha=0.01):
    """
    最急降下法メソッド
    :param x: テストデータarray
    :param y: 目標データarray
    :param iter_num: イテレーションの回数
    :param alpha: 学習率
    :return: 学習後のθ
    """
    theta = np.array([0, 0])  # 1x2行列
    m = len(x)  # データ数

    # ----以下で目的関数Jの偏微分を計算して繰り返しθの値を更新する
    for i in range(iter_num):
        h = np.dot(x, theta)  # 1x205行列   h(θ) = θ0 + θ1X を計算
        J = np.sum(np.square(h-y))/(2*m)  # 目的関数。Jはスカラー値になる。
        print("J=", J)
        tmp1 = theta[0] - alpha*np.sum(h-y)/m   # θ0の計算
        tmp2 = theta[1] - alpha*np.sum((h-y)*x[:, 1])/m  # θ1の計算
        theta[0] = tmp1   # θ0の更新
        theta[1] = tmp2   # θ1の更新

    return theta

if __name__ == "__main__":

    print("Linear Regression for Car Price Evaluation!")

    # -----データの読込、変数初期化
    alpha = 0.0001  #学習率α この値によっては適切に学習できない場合がある
    iter_num = 30  #イテレーションの回数
    data = np.loadtxt("car_2.csv", delimiter=",")
    engine_size = data[:, 0]    # 1x205行列
    x = np.c_[np.ones((len(engine_size),)), engine_size]  # 205x2 最急降下法関数で行列計算を行うため列を足している。
    y = data[:, 1]  # 1x205行列

    # -----最急降下法処理をiter_num分実施
    theta = gradient_descent(x, y, iter_num, alpha)
    print("[theta_0 theta_1] = ", theta)

    # -----後処理
    ans = np.dot(theta, np.array([1, 290]))  # x=290の時のyを計算
    print("Answer(x=290): ", ans)

実行時の結果を以下に示します。

$ python sample.py
Linear Regression for Car Price Evaluation!
J= 118634960.46
J= 84567747.8682
J= 60867761.8781
J= 44969646.0721
J= 33868970.9527
J= 26389948.1617
J= 21570569.0423
J= 18134456.6841
J= 15760941.5348
J= 14182798.9677
J= 13186249.2811
J= 12350034.3706
J= 11881632.6393
J= 11484490.8085
J= 11312642.3557
J= 11158608.8781
J= 11022390.3756
J= 10903986.8483
J= 10803398.296
J= 10803398.296
J= 10803398.296
J= 10803398.296
J= 10803398.296
J= 10803398.296
J= 10803398.296
J= 10803398.296
J= 10803398.296
J= 10803398.296
J= 10803398.296
J= 10803398.296
[theta_0 theta_1] = [  0 105]
Answer(x=290):  30450

結果を見ると最終的に、目的関数の値が収束して動いていないことが確認できます。これにより適切な $J(θ_0, θ_1)$ の値が定まりました。
求めた $J(θ_0, θ_1)$ で結果を図示してみました。

[横軸:エンジンサイズ( $x$ ), 縦軸:車の値段( $y$ )]
f:id:mhigu:20151027194528p:plain
直感的には、 $x$ の大きい時にあまり正確に予測できていないように見えますが、アルゴリズム的にはこれが最適解という事になります。
ただし、この結果は学習率やイテレーションの回数を調整した結果ですので、学習率やイテレーションの回数を変えると結果は変わってきます。特に学習率を変えると収束せずに発散したりします。どうしてそうなるのかは今回は説明しませんが是非考えて見て下さい。