線形回帰でデータに直線を当てはめる - 薬剤師のプログラミング学習日記

これまでに何度か片対数プロットした血中濃度データから直線のパラメータを求めるということをしましたが、これは線形回帰の中でも説明変数が一つの単回帰に当たります。
今まではnumpyやscipyといったライブラリから既成のモジュールを呼び出してデータを与えれば結果（傾きと切片）が返ってきましたが、今回は最急降下法というアルゴリズムを実装してこれをやってみようと思います。

なお、この記事はCourseraのMachine Learning（機械学習）講座『Linear Regression with One Variable』の内容を基に書いています。

線形回帰のモデル表現
目的関数を定義する
最急降下法（Gradient Descent）のアルゴリズム
Pythonで実装する

線形回帰のモデル表現

f:id:enokisaute:20200412161225p:plain
例えば、あるレストランの経営会社が都市の人口 $x$ からその都市のレストランの利益 $y$ を線形回帰で予測する場合、xからyを出力する関数は以下のように表現します。

${y = h_\theta (x) = \theta_0 + \theta_1 x}$

$x$ ：説明変数（目的変数に作用する変数）
$y$ ：目的変数（予測したい変数）
$\theta_0 、 \theta_1$ ：モデルのパラメータ

この $h_\theta (x)$ を機械学習では仮説関数(hypothesis function)と呼び、このモデルは説明変数が一つの時の線形回帰（単回帰）となります。
データに最もうまく当てはまるようなパラメータ $\theta_0 、 \theta_1$ を見つけることが実装するアルゴリズムの目標です。

目的関数を定義する

仮説関数 $h_\theta (x)$ がデータに対してどれくらいうまく当てはまっているかを測るために、次のような関数*1を定義します。

$J(\theta_0, \theta_1) = \dfrac {1}{2m} \displaystyle \sum _{i=1}^m \left (h_\theta (x_{i}) - y_{i} \right)^2$

数式では関数 $h_\theta (x)$ の出力と実際の値との差の二乗の合計を $ｍ$ （データ数）で割って平均を求めています。この $J(\theta_0, \theta_1)$ が返すのは、 $h_\theta (x)$ と実際の値yとのずれの大きさであり、 $h_\theta (x)$ がどれくらいデータに当てはまっているかの指標となります。
予測 $h_\theta (x)$ とデータとのずれが大ききれば $J(\theta_0, \theta_1)$ も大きく、両者のずれが小さければ小さい値を返します。
ちなみに、２で割っていますがこれは後の計算が少し楽になるというだけで、結果には影響しません。

最急降下法（Gradient Descent）のアルゴリズム

上の目的関数が最小となるような $\theta$ を見つける方法である最急降下法とは、以下のようなアルゴリズムです。

$\theta_j$ (ここではj=0, 1)に適当な初期値を設定する。
現在のθにおける目的関数 $J(\theta)$ の勾配とは逆方向に各θの値を更新する。
$\theta$ の更新幅が充分小さくなるまで２を繰り返す。

目的関数の勾配とは、各θの偏微分 $\frac{\partial}{\partial \theta_j}J(\theta)$ を並べたもので、関数の値が一番大きくなる方向を示しています。
したがって、この勾配とは逆方向にθを更新することにより $J(\theta)$ を減少させることができ、それを繰り返すことによって目的関数の最小値を見つける、というのが最急降下法のやり方です。

今回のモデルにおける２のステップを数式で表すと、次のように書くことができます。

${\theta_j := \displaystyle \theta_j - \alpha\frac{\partial}{\partial \theta_j}J(\theta_0, \theta_1) = \displaystyle \theta_j - \alpha \frac{1}{m} \sum_{i=1}^m (h_\theta (x^{(i)}) - y^{(i)}) x_j ^{(i)}}$

ここで、 $α$ は学習率（learning rate）と呼ばれ、θを更新する際の更新幅の大きさを制御します。

また、ポイントとして

$\theta_0 、 \theta_1$ の更新は同時に行うこと。
$α$ は小さすぎると更新の１ステップが少しずつになり、なかなか収束しない。また、大きすぎても最小値を通り越してしまい、収束しないこともあるため、「ちょうど良い値」を設定しなくてはいけない。

があります。

Pythonで実装する

データは講座のプログラミング課題のものを使っていますが、１０個だけランダムに抽出したものを使用*2し、これら１０個のデータに対し最急降下法を使ってパラメータ $\theta_0 、 \theta_1$ を求めています。

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D
　
# 線形回帰でデータに直線を当てはめる
　
# 目的関数
def compute_cost(x, y, theta):
    m = len(y)
    hypothesis = np.dot(x, theta)
    diff = hypothesis - y.reshape(m, 1)
    sqr_errors = np.square(diff)
    j = 1 / (2 * m) * np.sum(sqr_errors)
    return j
　

# 最急降下法でJ(θ)の値を最小化する.
# イテレーション毎にcompute_cost()を呼び出し, コストを出力する.
def gradient_descent(x, y, theta, alpha, num_iters):
    m = len(y)
    j_history = np.zeros((num_iters, 1))
    for i in range(0, num_iters):
        hypothesis = np.dot(x, theta)
        diff = hypothesis - y.reshape(m, 1)
        theta = theta - alpha / m * np.dot(x.T, diff)
        j_history[i] = compute_cost(x, y, theta)
    return theta, j_history
　

if __name__ == '__main__':
    # X: 利益(1万ドル単位), y: 都市の人口(1万人単位)
    X = np.array([8.51, 18.95, 6.35,  5.71, 11.7, 10.23,  13.17, 7.22, 10.95, 12.82])
    y = np.array([4.24, 17.05, 5.49,  3.25,  8.00, 7.77, 14.69, 3.34, 7.04, 13.5])
    m = len(X)      # トレーニングデータの数
　
    intercept_term = np.ones((m, 1))
    X = np.hstack((intercept_term, X.reshape(m, 1)))   # 1列目に切片項を追加
    theta = np.zeros((2, 1))      # thetaを0で初期化
　
    iterations = 1500
    alpha = 0.01        # 学習率
　
    theta, j_history = gradient_descent(X, y, theta, alpha, iterations)
    print('theta: \n', theta)
　
    # 勾配降下が正しく機能していることを確認するため各ステップで値が減少していることを確認する.
    plt.figure()
    plt.plot(j_history)
    plt.xlim(0,)
    plt.xlabel('iterations')
    plt.ylabel('J(theta)')
　
    plt.figure()
    plt.scatter(X[:, 1], y, color='red', marker='x', alpha=0.5)
    plt.ylabel('利益(1万ドル単位)')
    plt.xlabel('都市の人口(1万人単位)')
    plt.plot(X[:, 1], np.dot(X, theta))
    plt.show()
　

    # コスト関数J(θ)はボウル型でグローバルな最小値を持つ(等高線プロットで確認).
    # 勾配降下の各ステップはこのポイントに近づく.
    theta0_vals = np.linspace(-10, 10, 100)
    theta1_vals = np.linspace(-1, 4, 100)
    t0, t1 = np.meshgrid(theta0_vals, theta1_vals)
　
    j_vals = np.zeros((len(theta0_vals), len(theta1_vals)))
　
    t = np.array([t0, t1])
　
    for i in range(len(theta0_vals)):
        for j in range(len(theta1_vals)):
            j_vals[i, j] = compute_cost(X, y, t[:, i, j].reshape(2, 1))
　

    fig = plt.figure(figsize=(9, 4))
    fig.suptitle('Cost function J(θ)')
    plt.subplots_adjust(wspace=0.5, hspace=0.5)
　
    # J(theta)の3次元プロット
    ax1 = fig.add_subplot(121, projection='3d')
    ax1.plot_surface(t0, t1, j_vals, cmap=cm.coolwarm)
    ax1.set_title('Surface')
    ax1.set_xlabel(r'$\rmθ_0$')
    ax1.set_ylabel(r'$\rmθ_1$')
    ax1.set_zlabel(r'$\rmJ(θ_0, θ_1)$')

    # 等高線の描画
    ax2 = fig.add_subplot(122)
    ax2.contour(t0, t1, j_vals, levels=np.logspace(-2, 3, 20))
    ax2.plot(theta[0], theta[1], 'x', color='red')
    ax2.set_title('Contour, showing minimum')
    ax2.set_xlabel('θ0')
    ax2.set_ylabel('θ1')
    plt.show()

compute_cost()が目的関数、gradient_descent()が最急降下法です。上の数式をそのまま実装していますが、こういったコードを書く場合はデータ1つ1つについてループを回すのではなく、ベクトルや行列としてまとめて計算するのがコツだそうです。このようなスタイルで書くと、コードがシンプルになるだけでなく、より効率的にもなるとのこと。

gradient_descent()ではθの更新を繰り返すごとに $J(\theta)$ の値を出力しこれをj_historyに記録して、後でプロットしています。 $J(\theta)$ の値は繰り返し毎に減少していくはずなので、この変化を確認することで、プログラムが正しく動いているかどうか、 $α$ の値が適切かどうかを知ることができます。
f:id:enokisaute:20200415010755p:plain
最終的にパラメータ $\theta_0 、 \theta_1$ は

theta: 
 [[-3.04638561]
 [ 1.09496772]]

となり、この値を使って $: h_\theta (x)$ を描画すると次の図のようになります。
f:id:enokisaute:20200415010800p:plain

下図は目的関数 $J(\theta_0, \theta_1)$ の３次元プロットとそれを等高線で描画したものです。
線形回帰の場合の目的関数は必ず凸関数（弓型）となり、この関数は大域的な最小値以外の局所最小値は持たないのだそう。局所的な最小値ではないか、との心配はいらないということですね。
等高線プロットの方は赤い×点が最終的に求まった $\theta_0 、 \theta_1$ を示しています。

f:id:enokisaute:20200415010757p:plain

講義を受けているときはOctaveというプログラミング言語でコードを書いていましたが、今回Pythonで書き直して良い復習になったと思います。