線形回帰のまとめ（その３、Pythonで実践してみる）

前回は、線形回帰の手法として用いられる、最急降下法について説明しました。

今回は、これまで導いてきた式を使って、実際に最急降下法を試してみたいと思います。
私の環境は、Python3で作った、Jupyter Notebook上で実行しています。

テストデータの生成

ライブラリScikitLearnを使えば、相関のあるランダムなデータセットを作ることができます。

# 必要なライブラリをロード
%matplotlib inline
import numpy as np
import time
from matplotlib import cm
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.axes3d import Axes3D
from sklearn import datasets
import os
font = {'family':'IPAexGothic'}
plt.rc('font', **font)

# テストデータを作成

if not os.path.isdir("./data"):
    os.makedirs("./data")

if os.path.isfile("./data/linear1.dat") == True:
    data = np.loadtxt("./data/linear1.dat", delimiter = ",")
    X = data[:,0][:,np.newaxis]
    Y = data[:,1][:,np.newaxis]

else:
    X, Y = datasets.make_regression(n_samples=100, n_features=1, n_informative=1, noise=50)
    np.savetxt("./data/linear1.dat", np.c_[X, Y], delimiter=",")
    Y = Y[:,np.newaxis]

これで特徴量１、それに相関するデータ１、いくらかのノイズを乗せたデータセットがファイルに出力されました。

グラフにプロットしてみましょう。

# (x, y)データをプロットする関数
def plotData(X, Y):
    plt.xlim([-4,4])
    plt.ylim([-250,250])
    plt.scatter(X,Y)
    plt.xlabel(u'X')
    plt.ylabel(u'Y')
    plt.draw()
plotData(X, Y)

数式で解く

まずは数式にあてはめれば解が出ますので、当てはめてみましょう。シグマだらけの式でしたが、sumは関数でサクッと計算してくれます。

# 数式で解を求める
m = Y.size
w0 = (sum(X**2) * sum(Y) - sum(X) * sum(X*Y)) / (m * sum(X**2) - sum(X)**2)
w1 = (sum(X) * sum(Y) - m * sum(X*Y)) / (sum(X)**2 - m * sum(X**2))
w0,w1

(array([3.40246005]), array([71.96840595]))

解は
w0=3.40246005
w1=71.96840595
となりました。

行列を使わない力技で最急降下させる

次はw0, w1の変位量をそれぞれ計算し移動させていく方法です。
とりあえずα=0.1として、1000回ループしています。

# 目的関数
def computeCost(X, Y, w0, w1):
    m = Y.size
    L = 0
    for i in range(m):
        L = L + (w0 + w1 * X[i] - Y[i])**2 / (2 * m)
    return L

# w0, w1を低い方向へ移動させる関数
def doDescent(X, Y, w0, w1, alpha):
    m = Y.size
    w0 = w0 - sum(w0 + w1 * X -Y) * alpha / m
    w1 = w1 - sum((w0 + w1 * X - Y) * X) * alpha / m
    return w0, w1

alpha = 0.1
for i in range(1000):
    w0, w1 = doDescent(X, Y, w0, w1, alpha)
w0,w1

(array([3.40246005]), array([71.96840595]))

同様の解が出力されました。

行列を使った方法で最急降下させる

続いて行列を使った魔法のような式を使ってベクトルWを移動させていく方法です。
同じくα=0.1として、1000回ループしています。

# 行列を使った最急降下法
# Xはデータ列のみなので、1列目に1の列を追加してやる
# (1, x^(i))の形の行列を作る
Xb=np.c_[np.ones([X.size,1]),X]
# Wの初期値(0,0)を作る
W=np.zeros([2,1])

# 行列版の目的関数
def computeCost(Xb, Y, W):
    m = Y.size
    L = ((Xb.dot(W) - Y)**2).sum() / (2 * m)
    return L

# 下りのベクトル方向へWを動かす関数
def doDescent(Xb, Y, W, alpha):
    m = Y.size

    W = W - (alpha / m) * Xb.T.dot(Xb.dot(W)-Y)

    return W

alpha = 0.1
for i in range(1000):
    W = doDescent(Xb, Y, W, alpha)

(array([3.40246005]), array([[71.96840595]]))

また同じ解が出力されましたね。

ScikitLearnのライブラリを使う

ScikitLearnの中には、線形回帰のライブラリがあります。
αもループ回数も指定する必要がありません…
（上で説明したやつは何だったんだ、って思っちゃいますが…）
LinearRegressionをインスタンス化して、fitの一言で終わってしまいます。最強です。

from sklearn.linear_model import LinearRegression
linear = LinearRegression()

linear.fit(X, Y)
linear.intercept_, linear.coef_

(array([3.40246005]), array([[71.96840595]]))

もちろん解は同じですね。

最後に、プロットして確認

y=w0 + w1xをグラフにプロットして、正しく計算されたか見てみましょう。

def plotDataWithLine(X, Y, w0, w1):
    plotData(X, Y)
    plt.hold(True)
    plt.plot([-4, 4], [(-3)*w1+w0, 3*w1+w0]) 
    plt.show()

ばっちりですね。

これで線形回帰はマスターしました。次は、4通りの方法でのパフォーマンスを比較してみましょうか。

CHAR-IN-CO.SE