いつの間にかシリーズ化して、今回はロジスティック回帰をやる。自分は行列計算ができないクラスタ所属なので、入力が3次元以上 / 出力が多クラスになるとちょっときつい。教科書を読んでいるときはなんかわかった感じになるんだが、式とか字面を追ってるだけだからな、、、やっぱり自分で手を動かさないとダメだ。
また、ちょっとした事情により今回は Python でやりたい。Python のわかりやすい実装ないんかな？と探していたら 以下の ipyton notebook を見つけた。
http://nbviewer.ipython.org/gist/mitmul/9283713
こちらのリンク先に2クラス/多クラスのロジスティック回帰 (確率的勾配降下法) のサンプルがある。ありがたいことです。理論的な説明も書いてあるので ロジスティック回帰って何？という方は上を読んでください (放り投げ)。
この記事では、"自分的ロジスティック回帰わからないポイント"である勾配法を中心に書く。
補足 この記事ではデータを pandas で扱う / matplotlib でアニメーションさせる都合上、元のプログラムとは処理/変数名などだいぶ変わっているため見比べる場合は注意。
補足 実際の分析で使いたい場合は scikit-learn の SGDClassifier 。
サンプルデータのロードと準備
みんな大好き iris で、"Species" を目的変数にする。
補足 rpy2 を設定している方は rpy2から、そうでない方は こちら から .csv でダウンロードして読み込み (もしくは read_csv のファイルパスとして直接 URL 指定しても読める)。
import numpy as np import pandas as pd import matplotlib.pyplot as plt import matplotlib.animation as animation # iris の読み込みはどちらかで # rpy2 経由で R から iris をロード # import pandas.rpy.common as com # iris = com.load_data('iris') # csv から読み込み # http://aima.cs.berkeley.edu/data/iris.csv names = ['Sepal.Length', 'Sepal.Width', 'Petal.Length', 'Petal.Width', 'Species'] iris = pd.read_csv('iris.csv', header=None, names=names) np.random.seed(1) # 描画領域のサイズ FIGSIZE = (5, 3.5)
ロジスティック回帰 (2クラス)
まず iris データの一部を切り出して、説明変数 2次元 / 目的変数 2クラスのデータにする。
# 2 クラスにするため、setosa, versicolor のデータのみ抽出 data = iris[:100] # 説明変数は 2つ = 2 次元 columns = ['Petal.Width', 'Petal.Length'] x = data[columns] # データ (説明変数) y = data['Species'] # ラベル (目的変数)
ロジスティック回帰での目的変数は 0, 1 のみからなるベクトルである必要があるので、以下のようにして変換。論理演算で bool の pd.Series を作り、int 型に変換すればよい。
# ラベルを0, 1の列に変換 y = (y == 'setosa').astype(int) y # 1 1 # 2 1 # ... # 99 0 # 100 0 # Name: Species, Length: 100, dtype: int64
クラス数は変数 d として保存しておく。Series.unique で当該データの ユニークな要素数が取れる。d はユニークな要素の長さになるので、
data['Species'].unique() # array(['setosa', 'versicolor'], dtype=object) # クラス数 = 2 d = len(data['Species'].unique())
プロットしてみる。線形分離できるのでサンプルとしては OK だろう。
def plot_x_by_y(x, y, colors, ax=None): if ax is None: # 描画領域を作成 fig = plt.figure(figsize=FIGSIZE) # 描画領域に Axes を追加、マージン調整 ax = fig.add_subplot(1, 1, 1) fig.subplots_adjust(bottom=0.15) x1 = x.columns[0] x2 = x.columns[1] for (species, group), c in zip(x.groupby(y), colors): ax = group.plot(kind='scatter', x=x1, y=x2, color=c, ax=ax, figsize=FIGSIZE) return ax plot_x_by_y(x, y, colors=['red', 'blue']) plt.show()
最適化法
ロジスティック回帰における予測値 は、
で計算される。それぞれの意味は、
: 真のラベル
の予測値ベクトル。次元は (入力データ数, 1)
: シグモイド関数。計算結果を (0, 1) 区間に写像する
: 入力データ。次元は (入力データ数, 説明変数の数)
: 係数ベクトル。次元は(クラス数 = 2, 1)
: バイアス (スカラー)
ロジスティック回帰では、予測値 と真のラベル
の差から計算される損失関数を最小化する回帰直線を求めることになる。この回帰直線を求める方法として勾配法が知られている。勾配法とは、損失関数から 勾配 (Gradient) を求めることによって損失関数をより小さくする
,
を逐次更新していく手法 (詳細は上記リンク)。
ここでは以下 3 種類の勾配法を試す。上のリンクでやっているのは 3つ目の ミニバッチ確率的勾配降下法 / MSGD になる。
- 勾配降下法 (Gradient Descent)
- 確率的勾配降下法 (Stochastic Gradient Descent / SGD)
- ミニバッチ確率的勾配降下法 (Minibatch SGD / MSGD)
※ MSGD は SGD と特に区別されない場合も多い。
まず、共通で使う関数を定義する。
def p_y_given_x(x, w, b): # x, w, b から y の予測値 (yhat) を計算 def sigmoid(a): return 1.0 / (1.0 + np.exp(-a)) return sigmoid(np.dot(x, w) + b) def grad(x, y, w, b): # 現予測値から勾配を計算 error = y - p_y_given_x(x, w, b) w_grad = -np.mean(x.T * error, axis=1) b_grad = -np.mean(error) return w_grad, b_grad
勾配降下法 (Gradient Descent)
勾配法の中ではもっともシンプルな考え方で、全データ (x, y) を使って勾配を求める。実装は以下のようになる。
引数の意味は、
x,y,w,b: 上記数式に対応eta: 学習率 (勾配をどの程度 w, b に反映させるか)num: イテレーション回数
補足 num の既定値はそれぞれの手法がある程度収束する値にしている。ちゃんとやる場合は収束判定すべき。
def gd(x, y, w, b, eta=0.1, num=100): for i in range(1, num): # 入力をまとめて処理 w_grad, b_grad = grad(x, y, w, b) w -= eta * w_grad b -= eta * b_grad e = np.mean(np.abs(y - p_y_given_x(x, w, b))) yield i, w, b, e
実行してみる。w, b の初期値は全て 0 のベクトル / スカラーにした。
# w, b の初期値を作成 w, b = np.zeros(d), 0 gen = gd(x, y, w, b) # gen はジェネレータ gen # <generator object gd at 0x11108e5f0> # 以降、gen.next() を呼ぶたびに 一回 勾配降下法を実行して更新した結果を返す。 # タプルの中身は (イテレーション回数, w, b, 誤差) gen.next() # (1, array([-0.027 , -0.06995]), 0.0, 0.47227246182037463) gen.next() # (2, array([-0.04810306, -0.12007078]), 0.0054926687253766763, 0.45337584157628485) gen.next() # (3, array([-0.06522081, -0.15679859]), 0.014689105435764377, 0.44035208019270661) gen.next() # (4, array([-0.07963026, -0.18443466]), 0.026392079742058178, 0.43101983106700503) # ...
実行のたびに w, b が更新されるとともに、誤差の値が小さくなっていることがわかる。ここでの誤差 は 誤判別率 ではなく、真のクラスと クラス所属確率の差。
確率的勾配降下法 (Stochastic Gradient Descent / SGD)
勾配降下法では勾配の計算を全データを行列とみて行うため、データ量が増えると計算が厳しい。確率的勾配降下法では、入力をランダムに (確率的に) 並べ替えて、データ一行ずつから勾配を求めて係数/バイアスを更新する。一行ずつの処理になるため計算量の問題は解消する。実装は以下のようになる。
def sgd(x, y, w, b, eta=0.1, num=4): for i in range(1, num): for index in range(x.shape[0]): # 一行ずつ処理 _x = x.iloc[[index], ] _y = y.iloc[[index], ] w_grad, b_grad = grad(_x, _y, w, b) w -= eta * w_grad b -= eta * b_grad e = np.mean(np.abs(y - p_y_given_x(x, w, b))) yield i, w, b, e
確率的に処理するため、データをランダムに並べ替える。ランダムに並べ替えられた結果の表示から、x, y の index が一致している = 対応関係が維持されていることを確認。
# index と同じ長さの配列を作成し、ランダムにシャッフル indexer = np.arange(x.shape[0]) np.random.shuffle(indexer) # x, y を シャッフルされた index の順序に並び替え x = x.iloc[indexer, ] y = y.iloc[indexer, ] x.head() # Petal.Width Petal.Length # 81 1.1 3.8 # 85 1.5 4.5 # 34 0.2 1.4 # 82 1.0 3.7 # 94 1.0 3.3 y.head() # 81 0 # 85 0 # 34 1 # 82 0 # 94 0 # Name: Species, dtype: int64
実行してみよう。
# w, b の初期値を作成 w, b = np.zeros(d), 0 gen = sgd(x, y, w, b) # 以降、gen.next() を呼ぶたびに 一行分のデータで更新した結果を返す。 # タプルの中身は (イテレーション回数, w, b, 誤差) gen.next() # (1, array([-0.055, -0.19 ]), -0.050000000000000003, 0.43367318357380596) gen.next() # (1, array([-0.09571092, -0.31213277]), -0.07714061571720722, 0.4071586323043157) gen.next() # (1, array([-0.08310602, -0.22389845]), -0.014116100082250033, 0.42197958122210588)
gen.next によって 一行処理した結果が返ってくる。一行あたりの処理をみると損失は増加していることもある (が徐々に減っていく)。
補足 イテレーション回数は全データに対する回数のため、データが一回りするまで増えない。
ミニバッチ確率的勾配降下法 (Minibatch SGD / MSGD)
確率的勾配降下法について、一行ずつではなく 計算できる程度の固まり (Minibatch) 単位で処理していく方法。実装は以下のようになる。
※ 事前にデータをランダムに並べ替えるのは確率的勾配降下法と同様。
引数はこれまでから一つ増えて、
batch_size: Minibatch として切り出すデータ数 (行数)
def msgd(x, y, w, b, eta=0.1, num=25, batch_size=10): for i in range(1, num): for index in range(0, x.shape[0], batch_size): # index 番目のバッチを切り出して処理 _x = x[index:index + batch_size] _y = y[index:index + batch_size] w_grad, b_grad = grad(_x, _y, w, b) w -= eta * w_grad b -= eta * b_grad e = np.mean(np.abs(y - p_y_given_x(x, w, b))) yield i, w, b, e
実行。
# w, b の初期値を作成 w, b = np.zeros(d), 0 gen = msgd(x, y, w, b) # 読み方/補足は 確率的勾配降下法と同じ。 gen.next() # (1, array([ 0.011 , 0.0725]), 0.050000000000000003, 0.52633544830099133) gen.next() # (1, array([ 0.02251519, 0.13792904]), 0.096115503575706834, 0.54788728897754557)
可視化
以上 3 つの処理を見比べてわかるとおり、差異は勾配を計算するのに利用するデータの切り出し処理のみ。それぞれどんな動きをしているのか可視化してみる。
def plot_logreg(x, y, fitter, title): # 描画領域を作成 fig = plt.figure(figsize=FIGSIZE) # 描画領域に Axes を追加、マージン調整 ax = fig.add_subplot(1, 1, 1) fig.subplots_adjust(bottom=0.15) # 描画オブジェクト保存用 objs = [] # 回帰直線描画用の x 座標 bx = np.arange(x.iloc[:, 0].min(), x.iloc[:, 0].max(), 0.1) # w, b の初期値を作成 w, b = np.zeros(d), 0 # 勾配法の関数からジェネレータを生成 gen = fitter(x, y, w, b) # ジェネレータを実行し、勾配法 1ステップごとの結果を得る for i, w, b, e in gen: # 回帰直線の y 座標を計算 by = -b/w[1] - w[0]/w[1]*bx # 回帰直線を生成 l = ax.plot(bx, by, color='gray', linestyle='dashed') # 描画するテキストを生成 wt = """Iteration = {0} times w = [{1[0]:.2f}, {1[1]:.2f}] b = {2:.2f} error = {3:.3}""".format(i, w, b, e) # axes 上の相対座標 (0.1, 0.9) に text の上部を合わせて描画 t = ax.text(0.1, 0.9, wt, va='top', transform=ax.transAxes) # 描画した line, text をアニメーション用の配列に入れる objs.append(tuple(l) + (t, )) # データ, 表題を描画 ax = plot_x_by_y(x, y, colors=['red', 'blue'], ax=ax) ax.set_title(title) # アニメーション開始 ani = animation.ArtistAnimation(fig, objs, interval=1, repeat=False) plt.show() plot_logreg(x, y, gd, title='Gradient descent') plot_logreg(x, y, sgd, title='Stochastic gradient descent') plot_logreg(x, y, msgd, title='Minibatch SGD')
勾配降下法 (Gradient Descent)
勾配降下法では 全データから勾配をもとめるので、イテレーション時に損失が増えることはない。ただしデータ数が多くなると使えない。
確率的勾配降下法 (Stochastic Gradient Descent / SGD)
確率的勾配降下法では 一行ずつで勾配をもとめるので、処理される行によっては損失が増えることもある。が、全体を通じて徐々に損失は減っていく。
ミニバッチ確率的勾配降下法 (Minibatch SGD / MSGD)
ミニバッチ確率的勾配降下法では ある程度の固まりから勾配をもとめるので、確率的勾配降下法と比べると 更新処理ごとの損失のばらつきは小さい。
多項ロジスティック回帰 (多クラス)
続けて、目的変数のクラス数が 2より大きい場合 = 多項ロジスティック回帰を行う。データは同じく iris を使う。説明変数は 2次元のままとする。
data = iris # クラス数 d = len(data['Species'].unique()) x = data[columns] y = data['Species']
線形分離できないクラスもあるが、まあやってみよう。
多項ロジスティック回帰の場合は の次元が変わる。伴って、予測値を求める式は
となり関連する変数の次元も以下のように変わる。
: ソフトマックス関数。
の計算結果が各クラスに所属する確率を計算する。詳細は上k(略)
上記のとおり y を多項ロジスティック回帰の形式に変換する。
# ラベルを カテゴリに対応した 0, 1 からなる3列に変換 y = pd.DataFrame({'setosa': (y == 'setosa').astype(int), 'versicolor': (y == 'versicolor').astype(int), 'virginica': (y == 'virginica').astype(int)}) y.head() # setosa versicolor virginica # 1 1 0 0 # 2 1 0 0 # 3 1 0 0 # 4 1 0 0 # 5 1 0 0
続けて、多項ロジスティック回帰で使う関数を定義する。
def p_y_given_x(x, w, b): def softmax(x): return (np.exp(x).T / np.sum(np.exp(x), axis=1)).T return softmax(np.dot(x, w.T) + b) def grad(x, y, w, b): error = p_y_given_x(x, w, b) - y w_grad = np.zeros_like(w) b_grad = np.zeros_like(b) for j in range(w.shape[0]): w_grad[j] = (error.iloc[:, j] * x.T).mean(axis=1) b_grad[j] = error.iloc[:, j].mean() return w_grad, b_grad, np.mean(np.abs(error), axis=0)
ミニバッチ確率的勾配降下法 (Minibatch SGD / MSGD)
勾配法については考え方は一緒なので、MSGD のみ。ほかの手法の場合は ループ箇所とデータ切り出し処理を "2クラスのロジスティック回帰" の場合と同様に書き換えればよい。
def msgd2(x, y, w, b, eta=0.1, num=800, batch_size=10): for i in range(num): for index in range(0, x.shape[0], batch_size): # index 番目のバッチを切り出して処理 _x = x[index:index + batch_size] _y = y[index:index + batch_size] w_grad, b_grad, error = grad(_x, _y, w, b) w -= eta * w_grad b -= eta * b_grad e = np.sum(np.mean(np.abs(y - p_y_given_x(x, w, b)))) yield i, w, b, e
可視化
可視化部分も似たようなコードなので、全体は gist で。各クラスを分割する位置に少しずつ回帰直線が寄っていく。
補足 行列 w の出力がずれていて済まぬ、、済まぬ、、。 tex 記法がうまく動かなかったんや、、。
まとめ
ロジスティック回帰、とくに勾配法についてまとめた。目的変数のクラス数によってデータの次元 / 予測値算出のロジックは若干異なる。が、おおまかな流れは一緒。