線形回帰による直線フィッティング

パターン認識と機械学習（PRML）まとめ（2010/8/29）のつづき。

PRML3章の線形回帰（Linear Regression）を実装してみます。そういえば、3章の実装はまるまる無視していた。何でだろう？

今回は、受講している Courseraの機械学習コースの第一週目の課題を参考にしているのでPRMLと若干表記が違います。Courseraの課題ではOctaveというプログラミング言語で書いたのですが、ここではPythonを使います。

今回は、一変数の線形回帰で直線フィッティングする例を取り上げます。線形回帰は、「線形」という名前がついていますが、これはパラメータと変数の線形結合をモデルとして仮定するという意味で、直線でフィッティングするという意味ではないのだ！ここけっこう勘違いしやすい？というか講義で曲線フィッティングしているのに線形回帰ですとか言ってて「あれ？」となってしまった・・・

今回は、Courseraの講義で使用したデータファイル（ex1data1.txt）をそのまま使います。プログラムとともにgitで公開しているので必要ならダウンロードしてください。

https://github.com/sylvan5/PRML/tree/master/ch3

直線フィッティングの結果はこんな感じになります。

f:id:aidiary:20140401223614p:plain

PRMLでは、

正規方程式による解法
逐次更新による解法

の両方とも説明がありますが、今回は勾配降下法（Gradient Descent）を使った逐次更新によるパラメータ推定を試します（正規方程式の方も後で試す予定）。最適化するのは予測値と真の値の二乗誤差関数です。

$\displaystyle J(\theta) = \frac{1}{2m} \sum_{i=1}^m (h_{\theta} (x^{(i)}) - y^{(i)})^2$

ここで、mは訓練データ数で $x^{(i)}$ と $y^{(i)}$ はi番目のデータの組です。どちらも一次元データ。仮説は次のような線形モデルです。

$h_{\theta} (x) = \theta^{T} x = \theta_0 + \theta_1 x_1$

これは、直線の方程式なので直線でのフィッティング仮定しています。ここに、変数xの高次の項やシグモイド関数などの基底関数を導入することで曲線のフィッティングもできます。線形回帰による曲線フィッティングは次回やってみます。先の二乗誤差関数を最小化するパラメータ $\theta$ を勾配降下法で求めるのが今回の課題です。

勾配降下法

勾配降下法では、次式でパラメータ $\theta$ を収束するまで繰り返し更新します。

$\displaystyle \theta_j = \theta_j - \alpha \frac{1}{m} \sum_{i=1}^m (h_{\theta} (x^{(i)}) - y^{(i)}) x_j^{(i)}$

ここで、 $\alpha$ は学習率です。

Pythonによる実装

というわけで、Pythonによる実装です。線形回帰の結果の他に、更新を繰り返すたびに二乗誤差が最小化される様子をプロットした図、コスト関数 $J(\theta)$ の3次元プロット、コスト関数 $J(\theta)$ の等高線の図も描画してみました。

#coding:utf-8
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

"""3.1 勾配降下法で線形回帰モデル（一変数）のパラメータ推定"""

def plotData(X, y):
    plt.scatter(X, y, c='red', marker='o', label="Training data")
    plt.xlabel("Population of city in 10,000s")
    plt.ylabel("Profit in $10,000s")
    plt.xlim(4, 24)
    plt.ylim(-5, 25)

def computeCost(X, y, theta):
    m = len(y)  # 訓練データ数
    tmp = np.dot(X, theta) - y
    J = 1.0 / (2 * m) * np.dot(tmp.T, tmp)
    return J

def gradientDescent(X, y, theta, alpha, iterations):
    m = len(y)      # 訓練データ数
    J_history = []  # 各更新でのコスト
    for iter in range(iterations):
        theta = theta - alpha * (1.0 / m) * np.dot(X.T, np.dot(X, theta) - y)
        J_history.append(computeCost(X, y, theta))
    return theta, J_history

if __name__ == "__main__":
    # 訓練データをロード
    data = np.genfromtxt("ex1data1.txt", delimiter=",")
    X = data[:, 0]
    y = data[:, 1]
    m = len(y)  # 訓練データ数

    # 訓練データをプロット
    plt.figure(1)
    plotData(X, y)

    # 訓練データの1列目に1を追加
    X = X.reshape((m, 1))
    X = np.hstack((np.ones((m, 1)), X))

    # パラメータを0で初期化
    # yとthetaはわざわざ縦ベクトルにreshapeしなくてもOK
    # np.dot()では自動的に縦ベクトル扱いして計算してくれる
    theta = np.zeros(2)
    iterations = 1500
    alpha = 0.01

    # 初期状態のコストを計算
    initialCost = computeCost(X, y, theta)
    print "initial cost:", initialCost

    # 勾配降下法でパラメータ推定
    theta, J_history = gradientDescent(X, y, theta, alpha, iterations)
    print theta

    # コストの履歴をプロット
    plt.figure(2)
    plt.plot(J_history)
    plt.xlabel("iteration")
    plt.ylabel("J(theta)")

    # 直線をプロット
    plt.figure(1)
    plt.plot(X[:,1], np.dot(X, theta), 'b-', label="Linear regression")
    plt.legend()

    # 35,000と70,000のときの予測値を計算
    predict1 = np.dot(np.array([1, 3.5]), theta)
    predict2 = np.dot(np.array([1, 7.0]), theta)
    print "3.5 => %f" % predict1
    print "7.0 => %f" % predict2

    # J(theta)の三次元プロット
    fig = plt.figure(3)
    ax = fig.gca(projection='3d')
    gridsize = 200
    theta0_vals = np.linspace(-10, 10, gridsize)
    theta1_vals = np.linspace(-1, 4, gridsize)
    theta0_vals, theta1_vals = np.meshgrid(theta0_vals, theta1_vals)
    J_vals = np.zeros((gridsize, gridsize))
    for i in range(gridsize):
        for j in range(gridsize):
            t = np.array([theta0_vals[i, j], theta1_vals[i, j]])
            J_vals[i, j] = computeCost(X, y, t)
    ax.plot_surface(theta0_vals, theta1_vals, J_vals)
    plt.xlabel("theta_0")
    plt.ylabel("theta_1")

    # 等高線を描画
    fig = plt.figure(4)
    plt.contour(theta0_vals, theta1_vals, J_vals, np.logspace(-2, 3, 20))
    plt.plot(theta[0], theta[1], 'rx', markersize=10, linewidth=2)
    plt.xlabel("theta_0")
    plt.ylabel("theta_1")
    plt.show()

実行結果は、

initial cost: 32.0727338775
theta: [-3.63029144  1.16636235]
final cost: 4.48338825659
3.5 => 0.451977
7.0 => 4.534245

となります。最適なパラメータ $\theta$ は、(-3.63, 1.17) であることがわかります。また、そのときの二乗誤差は4.48です。下のような図も表示されます。

f:id:aidiary:20140401223621p:plain

この図からパラメータの更新を繰り返すたびに二乗誤差がどんどん小さくなることがわかります。学習率の値を大きくするとより速く収束する場合もあるけど、逆に発散してしまう場合もあるので要注意！このような図を書いておくと学習が順調に進んでいることが確認できるので、計算量は多少大きくなるけどオススメとのこと。最終的に4.48に収束しています。

f:id:aidiary:20140401223626p:plain

コスト関数の3次元プロットです。matplotlibだと描くのが少し大変。下に凸の関数で、(-3.63, 1.17) の辺りで最小値をとることがわかります。今回は、パラメータ2つだから3次元プロットできたけど、これ以上大きくなると可視化は難しいな。

f:id:aidiary:20140401223630p:plain

最後に等高線です。こちらも (-3.63, 1.17) で最適になることがわかります。

宿題で解いた Octave（Matlab）版のプログラムもあるけど、PythonよりOctaveの方がシンプルでわかりやすいかも。Pythonは少しごちゃごちゃしている感じがする？まあ数値計算専用言語でないから仕方ないけど。次回は、線形回帰で曲線フィッティングしてみます。一度でも自分で試しておけば「線形回帰は直線フィッティングしかできない」という誤解しなくて済みそうだから。