Pythonで線形回帰（単回帰） - ぺーぺーSEのブログ

線形回帰をPythonで解いてみる。
データセットは下記のBostonデータセットを使用する。

blog.pepese.com

データセットのロードまでは以下。

import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
from sklearn import datasets
%pylab inline --no-import-all

boston = datasets.load_boston()
rooms = boston.data[:, 5]
house_prices = boston.target

『%pylab inline --no-import-all』を入れておく事で、Jupyter Notebook上でグラフが表示できるようになる。（別窓にならない）データの様子が見たい場合は以下。

plt.scatter(rooms, house_prices)
plt.show()

部屋の数は、物件の価格と相関があるはず、という単回帰なモデルを考える。

ライブラリで解く

NumPyで解く

z = np.polyfit(rooms, house_prices, 1)
# rooms: x軸データ
# house_prices: y軸データ
# 1: 回帰モデルの多項式の次数
p = np.poly1d(z) # 引数に入れた（複数の）回帰係数をもつ多項式を生成
# p.c: 求めたモデルのパラメータ
p.c
# array([  9.10210898, -34.67062078])

下記で描画する。

plt.hold(True)
# データセット
plt.scatter(rooms, house_prices)
# 求めた単回帰のデータを3～10間で1刻みで生成
ax = np.arange(3, 10, 1)
ay = p(ax)
# 単回帰のデータをプロット
plt.plot(ax, ay)
# 描画・表示
plt.show()

『hold(True)』を指定することで、plot関数は前回にplotしたデータを破棄せずに、新しいグラフを重ねて表示する。

SciPyで解く

from scipy import stats

# 単回帰を求める
slope, intercept, r_value, _, _ = stats.linregress(rooms, house_prices)
# rooms: x軸データ
# house_prices: y軸データ
# slope: 単回帰直線の傾き
# intercept: 単回帰直線の切片（y = ax + b の b）
# r_value: 相関係数
slope, intercept
# (9.102108981180308, -34.670620776438554)

下記で描画する。

plt.hold(True)
# データセット
plt.scatter(rooms, house_prices)
# 求めた単回帰のデータを3～10間で1刻みで生成
ax = np.arange(3, 10, 1)
ay = slope * ax + intercept
# 単回帰のデータをプロット
plt.plot(ax, ay)
# 描画・表示
plt.show()

scikit-learnで解く

from sklearn import linear_model

clf = linear_model.LinearRegression(fit_intercept=True, normalize=False, copy_X=True, n_jobs=1)
# fit_intercept: False に設定すると切片を求める計算を含めない。(デフォルト値: True)
# normalize: True に設定すると、説明変数を事前に正規化。 (デフォルト値: False)
# copy_X: メモリ内でデータを複製してから実行するかどうか。 (デフォルト値: True)
# n_jobs: 計算に使うジョブの数。-1 に設定すると、すべての CPU を使って計算。 (デフォルト値: 1)

X = rooms
X = X.reshape(-1, 1) # Xの正規化、なぜか正規化しないとコケる
Y = house_prices

# 単回帰を求める
clf.fit(X, Y)
# clf.coef_: 回帰係数、単回帰直線の傾き
# clf.intercept_: 切片（誤差）
# clf.score(X, Y): 決定係数
clf.coef_, clf.intercept_
# (array([ 9.10210898]), -34.670620776438568)

下記で描画する。

plt.hold(True)
# データセット
plt.scatter(rooms, house_prices)
# 求めた単回帰のデータを3～10間で1刻みで生成
ax = np.arange(3, 10, 1)
ay = clf.coef_ * ax + clf.intercept_
# 単回帰のデータをプロット
plt.plot(ax, ay)
# 描画・表示
plt.show()

ライブラリを使わず自力で解く

CourseraのMachine Learningに沿った方法で解く。仮説（Hypothesis）は下記。

${ \displaystyle h_\theta(x) = \theta_0 x_0 + \theta_1 x_1 }$

ただし、 $x_0 = 1$ 。
上記に従って、データセットから $x$ 、 $y$ 、 $\theta$ を作成。

# x_0 = 1を設定
x_0 = np.ones(len(rooms)) # roomsはx_1なので、同じ数だけx_0 = 1を作る
x_1 = rooms
X = np.column_stack((x_0, x_1)) # x_0とx_1の行列作成

y = house_prices.reshape(len(house_prices), 1) # 「reshape」で「1×n」→「n×1」変換

# θは「theta_0 = 0, theta_1 = 0」で初期化
theta = np.matrix(np.array([0,0]))

X.shape, theta.shape, y.shape
# ((506, 2), (1, 2), (506, 1))

また、コスト関数の数式および実装は下記。

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

$m$ はデータセット数（ここでは、506）、 $i$ はデータセットの番号（何番目か、1～506）。

def computeCost(X, y, theta):
    inner = np.power(((X * theta.T) - y), 2)
    return np.sum(inner) / (2 * len(X))

勾配降下法の数式および実装は下記。

$\theta_j = \theta_j - \alpha \frac{\partial J}{\partial \theta_j}$

$\frac{\partial J}{\partial \theta_j} = \frac{1}{m} \sum_{i=1}^m (h(x^{(i)}) - y^{(i)}) x_j^{(i)}$

$j$ はパラメータの番号（何番目の $x_j$ や $\theta_j$ か、 $0 \leqq j \leqq 1$ 、パラメータ2個しかないから0番目か1番目）。

def gradientDescent(X, y, theta, alpha, iters):
    temp = np.matrix(np.zeros(theta.shape))
    parameters = int(theta.ravel().shape[1])
    cost = np.zeros(iters)

    for i in range(iters):
        # h(x) - y
        error = (X * theta.T) - y

        for j in range(parameters):
            # (h(x) - y) * x_j
            term = np.multiply(error, X[:,j])
            # - α 1/m Σ (h(x) - y) * x_j
            temp[0,j] = theta[0,j] - ((alpha / len(X)) * np.sum(term))

        theta = temp
        cost[i] = computeCost(X, y, theta)

    return theta, cost

求める。

# αを0.00001とおく
alpha = 0.00001
# 勾配降下法のループを1000とおく
iters = 1000

# 勾配降下法の実行
g, cost = gradientDescent(X, y, theta, alpha, iters)
g
# matrix([[ 0.5564119 ,  3.49684533]])

最初、 $\alpha = 0.01$ でやっていたが発散してしまったので小さくした。
下記で描画する。

plt.hold(True)
# データセット
plt.scatter(rooms, house_prices)
# 求めた単回帰のデータを3～10間で1刻みで生成
ax = np.arange(3, 10, 1)
ay = g[0, 1] * ax + g[0, 0]
# 単回帰のデータをプロット
plt.plot(ax, ay)
# 描画・表示
plt.show()