第20回　ロジスティック回帰の実装

今回は実践編、前回導出した確率的勾配降下法でのロジスティック回帰の学習を実装してみましょう。

環境はこれまでと同じくPython/numpy/matplotlibを用います。インストールなどの準備は第6回を参照してください。

パーセプトロンの実装の復習

今回のロジスティック回帰の実装は、連載第17回のパーセプトロンの実装と非常によく似たものになります。

そこでパーセプトロンの実装を再掲しておきます。パーセプトロンそのものの復習はここではしませんので、必要なら連載第15回をご覧ください。

import numpy as np
import matplotlib.pyplot as plt
import random

# データ点の個数
N = 100

# データ点のために乱数列を固定
np.random.seed(0)

# ランダムな N×2 行列を生成 = 2次元空間上のランダムな点 N 個
X = np.random.randn(N, 2)

def h(x, y):
    return 5 * x + 3 * y - 1  #  真の分離平面 5x + 3y = 1

T = np.array([ 1 if h(x, y) > 0 else -1 for x, y in X])

# 特徴関数
def phi(x, y):
    return np.array([x, y, 1])

w = np.zeros(3)  # パラメータを初期化(3次の 0 ベクトル)

np.random.seed() # 乱数を初期化
while True:
    list = range(N)
    random.shuffle(list)

    misses = 0 # 予測を外した回数
    for n in list:
        x_n, y_n = X[n, :]
        t_n = T[n]

        # 予測
        predict = np.sign((w * phi(x_n, y_n)).sum())

        # 予測が不正解なら、パラメータを更新する
        if predict != t_n:
            w += t_n * phi(x_n, y_n)
            misses += 1

    # 予測が外れる点が無くなったら学習終了(ループを抜ける)
    if misses == 0:
        break

# 図を描くための準備
seq = np.arange(-3, 3, 0.02)
xlist, ylist = np.meshgrid(seq, seq)
zlist = np.sign((w * phi(xlist, ylist)).sum())

# 分離平面と散布図を描画
plt.pcolor(xlist, ylist, zlist, alpha=0.2, edgecolors='white')
plt.plot(X[T== 1,0], X[T== 1,1], 'o', color='red')
plt.plot(X[T==-1,0], X[T==-1,1], 'o', color='blue')
plt.show()

ロジスティック回帰の実装

これをロジスティック回帰の実装に差し替えるにあたって、変更する必要があるのはラベル変数の値、予測確率の計算と、パラメータwの更新式、あとはパラメータの初期化とループの終了条件です。

順に見ていきましょう。

人工データの生成

パーセプトロンは線形分離可能な問題しか解けなかったので（連載第15回参照⁠）⁠、このコードでは線形分離可能なデータを生成しています。

一方、ロジスティック回帰は線形非分離な問題でも大丈夫なのですが、ここではパーセプトロンと同じ問題をきちんと（期待通りに）解けるかを確認しておきましょう。

というわけでデータ生成のコードは基本そのまま使いますが、正解ラベルの値だけはロジスティック回帰にあわせて変更する必要があります。

パーセプトロンの正解ラベルは+1と-1であるのに対し、ロジスティック回帰は1と0を使うため、その点のみ反映したコードは次のようになります。

TT = np.array([ 1 if f(x, y) > 0 else 0 for x, y in X])

これに関連して、出力の図を描くところも少し変更します。

パーセプトロンによって得られる予測は「正か負か⁠」⁠、特に+1または-1であり、描画する必要があるのは分離平面だけでした。一方、ロジスティック回帰の予測は確率、つまり0から1の間で変化する関数として得られます。

よって、なだらかに変化する予測分布を描画したいですから、そのための書き換えのついでに、今回はmatplotlibのimshow()という関数を使ってみましょう。

plt.imshow(zlist, extent=[-3,3,-3,3], origin='lower', cmap=plt.cm.PiYG_r)

imshow関数は、引数cmapに描画する色マップを指定します。PiYG_rの場合は、予測分布が1に近い点はピンク（Pi⁠）⁠、0に近い点は緑（G）に、そしてちょうど中間は白っぽい黄色（Y）となるグラデーションで描画されます（_rは逆順を意味する⁠）⁠。

予測確率とパラメータの更新式

パーセプトロンにしてもロジスティック回帰にしても、「⁠データを選んで、予測して、それを使ってパラメータを更新」という大きい流れは全く同じです。

パーセプトロンのコードでは、訓練データセットからランダムに1つ取りだした点（x_n,y_n⁠）⁠、その点での正解ラベルt_nを使って予測とパラメータの更新をし、それをデータセットを一周するまで繰り返すということをしていますが、その枠組みはロジスティック回帰でもそのまま使えるということですね。

一方、予測やパラメータの更新式はモデルの要ですから、こちらは当然書き換える必要があるでしょう。

ロジスティック回帰の予測確率は次の式で表されるのでした。

ただし

これを実装するには、まずシグモイド関数σ(z)が必要ですね。このシグモイド関数やその多変数版にあたるソフトマックス関数はnumpyかscipyで標準サポートしてくれててもいいのになあとか個人的には思うのですが、無いものは仕方ありません。

せっかくですから関数として実装しておきましょう。

# シグモイド関数
def sigmoid(z):
    return 1.0 / (1.0 + np.exp(-z))

上の式では予測確率にyの文字が当てられていますが、その文字はすでにデータ点のy座標を表すのに使ってしまいました。ここではパーセプトロンのコードにあわせてpredict（予測）という名前の変数を使うことにしましょう。

φ_nもこの後またすぐ使いますから、これは変数feature（特徴、素性）に入れておきます。すると次のようなコードになります。

        # 特徴量と予測確率
        feature = phi(x_n, y_n)
        predict = sigmoid(np.inner(w, feature))

特徴関数phiは定数項にあたる1を付け加えただけの、元のパーセプトロンのコードと同じものを使います。

予測確率がわかると、パラメータを更新できます。

ロジスティック回帰の確率的勾配降下法では、次の式でパラメータwを更新します。

ηは「学習率」と呼ばれる適当な正の値で、一回の更新でパラメータをどれくらい動かすかを調整します。ηが大きいほど学習は速くなりますが、予測確率が安定しません。そこで初期値は0.1くらいにして、学習が進むにつれて徐々に小さくしていく、というのが一般的です。

# 学習率の初期値
eta = 0.1

    # イテレーションごとに学習率を小さくする
    eta *= 0.9

このetaを使ってパラメータを更新します。

        w -= eta * (predict - t_n) * feature

パラメータの初期値と学習の終了条件

後はパラメータwの初期化と学習の終了条件だけです。この2つは実は少々関連があります。

パーセプトロンの解法アルゴリズムは、0で初期化したパラメータで始めるとすべてのデータ点を完全に正しく分離する解にたどり着くことが示されています。

この時の学習ループの終了条件は明確で、「⁠すべてのデータ点での予測が正解と一致するまで」となります。「⁠線形分離可能な問題に限る」という強い条件が付いているだけあって、非常に強い性質を持っているわけですね。

一方の確率的勾配降下法で求まるのは、パラメータwの初期値に依存した局所解になります。そこでランダムな初期値で始めるのが一般的ですが、パラメータが大きくなりすぎるのを嫌って、あえて0から始めることもあります。

w = np.random.randn(3)  # パラメータをランダムに初期化

そうして求めた局所解は、一般に全データ点で正解を予測できるものにはなりませんので、学習ループの終了条件もパーセプトロンのままというわけにはいきません。

そこで確率的勾配降下法では尤度などの指標値があまり変化しなくなったら終了したり、あるいはあらかじめ学習の回数を適当に決めておきます。

今回はシンプルに学習回数を50周と決めておきましょう。

これらをすべて反映したコードは次のようになります。

import numpy as np
import matplotlib.pyplot as plt

# データ点の個数
N = 100

# データ点のために乱数列を固定
np.random.seed(0)

# ランダムな N×2 行列を生成 = 2次元空間上のランダムな点 N 個
X = np.random.randn(N, 2)

def f(x, y):
    return 5 * x + 3 * y - 1  #  真の分離平面 5x + 3y = 1

T = np.array([ 1 if f(x, y) > 0 else 0 for x, y in X])

# 特徴関数
def phi(x, y):
    return np.array([x, y, 1])

# シグモイド関数
def sigmoid(z):
    return 1.0 / (1 + np.exp(-z))

np.random.seed() # 乱数を初期化
w = np.random.randn(3)  # パラメータをランダムに初期化

# 学習率の初期値
eta = 0.1

for i in xrange(50):
    list = range(N)
    np.random.shuffle(list)

    for n in list:
        x_n, y_n = X[n, :]
        t_n = T[n]

        # 予測確率
        feature = phi(x_n, y_n)
        predict = sigmoid(np.inner(w, feature))
        w -= eta * (predict - t_n) * feature

    # イテレーションごとに学習率を小さくする
    eta *= 0.9

# 図を描くための準備
seq = np.arange(-3, 3, 0.01)
xlist, ylist = np.meshgrid(seq, seq)
zlist = [sigmoid(np.inner(w, phi(x, y))) for x, y in zip(xlist, ylist)]

# 散布図と予測分布を描画
plt.imshow(zlist, extent=[-3,3,-3,3], origin='lower', cmap=plt.cm.PiYG_r)
plt.plot(X[T==1,0], X[T==1,1], 'o', color='red')
plt.plot(X[T==0,0], X[T==0,1], 'o', color='blue')
plt.show()

さて、長々と20回続けてきた本連載ですが、次回で最終回となります。

最終回は連載の流れの中では触れられなかった落ち穂拾いをしつつ、これまでのまとめをしようと思います。