線形回帰による曲線フィッティング - 人工知能に関する断創録

線形回帰による直線フィッティング（2014/4/1）のつづき。

今回は、線形回帰で曲線フィッティングをしてみます。PRMLによると、線形回帰で曲線のフィッティングをするためには、入力変数を非線形の基底関数（basis function）で変換すればよいそうです。今回は、元の変数を二乗、三乗した項を加えます。こういうのは多項式回帰（polynomial regression）っていいますね。多項式回帰も線形回帰の一種だったのだ。仮説は下のような3次多項式になります。入力変数はx、ただ一つである点に注意！

$h_{\theta} (x) = \theta^{T} x = \theta_0 + \theta_1 x + \theta_2 x^2 + \theta_3 x^3$

実はこれと同じことは前にもやったことがあった。PRMLの1章にある多項式曲線フィッティング（2010/3/27）の仮説がまさしくこの形だ。PRML1章では解析的に解いて直接パラメータを計算していたけど、今回は前回と同様に勾配降下法による逐次更新で求めてみます。学習データは、多項式曲線フィッティング（2010/3/27）と同様にサイン関数に誤差を加えたデータを作成しました。

https://github.com/sylvan5/PRML/blob/master/ch3/linear_regression_curve_fitting.py

#coding:utf-8
import numpy as np
import matplotlib.pyplot as plt

"""
勾配降下法で線形回帰による曲線フィッティング
PRML1章の多項式曲線フィッティングと同じだが、
正規方程式ではなく、勾配降下法で解く点が異なる
"""

def plotData(X, y):
    plt.scatter(X, y, c='red', marker='o', label="Training data")
    plt.xlabel("x")
    plt.ylabel("y")

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__":
    # 訓練データを作成
    # sin(2 * pi * x) の値にガウス分布に従う小さなランダムノイズを加える
    X = np.linspace(0, 1, 100)
    y = np.sin(2 * np.pi * X) + np.random.normal(0, 0.2, X.size)

    # 訓練データ数
    m = len(y)

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

    # 訓練データの1列目に1を追加
    # 元の変数の二乗、三乗の項も加える
    # 多項式回帰になり曲線フィッティングも可能
    X = X.reshape((m, 1))
    X = np.hstack((np.ones((m, 1)), X, X**2, X**3))

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

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

    # 勾配降下法でパラメータ推定
    theta, J_history = gradientDescent(X, y, theta, alpha, iterations)
    print "theta:", theta
    print "final cost:", J_history[-1]

    # コストの履歴をプロット
    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.xlim((0, 1))
    plt.legend()
    plt.show()

前回の線形回帰による直線フィッティング（2014/4/1）からほとんど変わりがありません。訓練データ集合 X に元の入力変数の二乗、三乗の値を追加した点とパラメータ theta を4次元ベクトルに拡張した点が追加の変更です。コスト関数や勾配降下法の関数は変更しなくてOK。出力は下のような感じ。

initial cost: 0.260955056593
theta: [ -0.10204345  10.19472864 -29.85813642  19.77507164]
final cost: 0.0271042753428

二乗誤差関数を最適なパラメータ theta が求まったので、これを先の3次多項式に代入してプロットすると下のような曲線が得られます。データに対して曲線でフィッティングしていることがわかります。

f:id:aidiary:20140402221552p:plain

勾配降下法の各繰り返しでコスト関数が減少して収束することも確認できます。繰り返し回数は10万回とかなり多いです。学習率 alpha によっては収束しない場合もあるので微調整が面倒でした。あとで学習率を指定しなくてもよい最適化アルゴリズムであるConjugate gradient、BFGSなども試してみたいと思います。

f:id:aidiary:20140402221559p:plain

次回は多変数の線形回帰、いわゆる重回帰をやってみます。