ロジスティック回帰
今回は、ロジスティック回帰です。この方法はPRMLで初めて知りましたが、統計学の方では一般的な方法のようです。回帰という名前がついてますが、実際は分類のためのモデルとのこと。ロジスティック回帰では、クラス1の事後確率が特徴ベクトルの線形関数のロジスティックシグモイド関数として書けることを利用しています。
ここで、σ(a)は式(4.59)のロジスティックシグモイド関数です。
訓練データ集合 {x_n, t_n} (今度は、クラス1のときt_n=0, クラス1のときt_n=1なので注意)からパラメータwを最尤推定で求めます。尤度関数は、
と書けるので、誤差関数(尤度関数の負の対数)は、
となります。誤差関数を最小化するようなwを求めたいってことですね。で、普通だったら今までのようにwで偏微分して0とおいてwを解析的に求めるところですが、yにロジスティックシグモイド関数が入っているせいで解析的に解けないようなのです。そこで4.3.3に出てくるニュートン・ラフソン法(Newton-Raphson method)を応用したwの反復更新手法を使います。このwの更新手法は、反復重み付き最小二乗法(iterative reweighted least squares method: IRLS)として知られているとのこと。ニュートン・ラフソン法によるwの更新式は、(4.92)で与えられます。
ここで、E(w)は式(4.90)です。∇がついているのでwで偏微分した行列になります。HはE(w)をwで二回微分した要素を持つヘッセ行列です。実際に、計算してみると
となります。これらを(4.92)に代入すると、最終的にwの更新式は、
となります。ちなみに決定境界は、
のときです。これは、図4.9のロジスティックシグモイド関数のグラフから分かりますが、ロジスティックシグモイド関数の引数が0のときです。
今回の例題では、2次元データ(バイアスの1が追加なので3次元)でφは恒等関数なので
整理すると
という直線の方程式が得られます。では、これをプログラムして試してみます。
#coding:utf-8 # 4.3.2 ロジスティック回帰 import numpy as np from pylab import * import sys N = 100 # データ数 def sigmoid(a): return 1 / (1 + np.exp(-a)) def f(x1, w): # 決定境界の直線の方程式 # p(C1|X) > 0.5 -> XをC1に分類 # p(C1|X) < 0.5 -> XをC2に分類 # p(C1|X) = 0.5が決定境界 <-> シグモイド関数の引数が0のとき(図4.9) a = - (w[1] / w[2]) b = - (w[0] / w[2]) return a * x1 + b if __name__ == "__main__": # 訓練データを作成 cls1 = [] cls2 = [] mean1 = [-1, 2] mean2 = [1, -1] mean3 = [8, -6] cov = [[1.0,0.8], [0.8,1.0]] # ノイズありデータ作成 cls1.extend(np.random.multivariate_normal(mean1, cov, N/2)) cls2.extend(np.random.multivariate_normal(mean2, cov, N/2-20)) cls2.extend(np.random.multivariate_normal(mean3, cov, 20)) # データ行列Xを作成 temp = vstack((cls1, cls2)) temp2 = ones((N, 1)) X = hstack((temp2, temp)) # ラベルTを作成(1-of-K表現ではないので注意) t = [] for i in range(N/2): t.append(1.0) for i in range(N/2): t.append(0.0) t = array(t) # パラメータwをIRLSで更新 turn = 0 w = array([0.0, 0.0, 0.0]) # 適当な初期値 while True: # Φを計算(恒等式とするのでデータ行列Xをそのまま使う) phi = X # Rとyを計算 R = np.zeros((N, N)) y = [] for n in range(N): a = np.dot(w, phi[n,]) y_n = sigmoid(a) R[n, n] = y_n * (1 - y_n) y.append(y_n) # ヘッセ行列Hを計算 phi_T = phi.T H = np.dot(phi_T, np.dot(R, phi)) # wを更新 w_new = w - np.dot(np.linalg.inv(H), np.dot(phi_T, y-t)) # wの収束判定 diff = np.linalg.norm(w_new - w) / np.linalg.norm(w) print turn, diff if diff < 0.1: break w = w_new turn += 1 # 訓練データを描画 x1, x2 = np.array(cls1).transpose() plot(x1, x2, 'rx') x1, x2 = np.array(cls2).transpose() plot(x1, x2, 'bo') # 識別境界を描画 x1 = np.linspace(-6, 10, 1000) x2 = [f(x, w) for x in x1] plot(x1, x2, 'g-') xlim(-6, 10) ylim(-10, 6) show()
結果は、
分類における最小二乗(2010/4/24)と違ってノイズがあっても決定境界が影響を受けないことが確認できます。PRMLの図4.4の緑色の決定境界がロジスティック回帰のものなので同じ結果が得られました。
次に、4.3.4の多クラスロジスティック回帰を試してみます。これは、IRLSの更新式が書いてなかったので2クラスと同じだろうと試行錯誤してみました。一応、うまく動いているように見えますが、正しいかわからないのでそこらへんは気をつけてください。特にp.209にヘッセ行列のブロックj, kとか書いてあるのがよくわからなかったのですが、クラスごとにヘッセ行列Hが必要なのか?と解釈しました。
#coding:utf-8 # 4.3.4 多クラスロジスティック回帰 import numpy as np from pylab import * import sys N = 300 # データ数 M = 3 # パラメータの次数 K = 3 # クラス数 def sigmoid(a): return 1 / (1 + np.exp(-a)) def f(x1, W_t, c1, c2): # クラスc1とクラスc2の決定境界の直線の方程式 a = - ((W_t[c1,1]-W_t[c2,1]) / (W_t[c1,2]-W_t[c2,2])) b = - ((W_t[c1,0]-W_t[c2,0]) / (W_t[c1,2]-W_t[c2,2])) return a * x1 + b if __name__ == "__main__": # 訓練データの作成 cls1 = [] cls2 = [] cls3 = [] mean1 = [-2, 2] # クラス1の平均 mean2 = [0, 0] # クラス2の平均 mean3 = [2, -2] # クラス3の平均 cov = [[1.0,0.8], [0.8,1.0]] # 共分散行列(全クラス共通) # データを作成 cls1.extend(np.random.multivariate_normal(mean1, cov, N/3)) cls2.extend(np.random.multivariate_normal(mean2, cov, N/3)) cls3.extend(np.random.multivariate_normal(mean3, cov, N/3)) # データ行列Xを作成 temp = vstack((cls1, cls2, cls3)) temp2 = ones((N, 1)) X = hstack((temp2, temp)) # ラベル行列Tを作成 T = [] for i in range(N/3): T.append(array([1, 0, 0])) for i in range(N/3): T.append(array([0, 1, 0])) for i in range(N/3): T.append(array([0, 0, 1])) T = array(T) # パラメータ行列WをIRLSで更新 turn = 0 W = np.zeros((M, K)) while True: # Φを計算(恒等式とするのでデータ行列Xをそのまま使う) phi = X # 予測値の行列Yを計算 Y = np.zeros((N, K)) for n in range(N): denominator = 0.0 for k in range(K): denominator += np.exp(np.dot(W[:,k], X[n,:])) for k in range(K): Y[n, k] = np.exp(np.dot(W[:,k], X[n,:])) / denominator print Y # ヘッセ行列Hを計算 I = np.identity(K) H = np.zeros((K*K, M, M)) for j in range(K): for k in range(K): # (4.110)に従った計算 for n in range(N): temp = Y[n, k] * (I[k, j] - Y[n,j]) H[k+j*K] += temp * matrix(phi)[n].reshape(M,1) * matrix(phi)[n].reshape(1,M) # 縦ベクトルx横ベクトル # もしくは下のように行列計算でもできる # Ijk = 1 if j == k else 0 # R = np.diag(Y[:,k] * (Ijk - Y[:,j])) # H[k+j*K] = np.dot(phi.T, np.dot(R, phi)) # 重み行列を更新 W_new = np.zeros((M, K)) phi_T = phi.T for i in range(K): temp = np.dot(phi_T, Y[:,i] - T[:,i]) W_new[:,i] = W[:,i] - np.dot(np.linalg.inv(H[i+i*K]), temp) # Wの収束判定 diff = np.linalg.norm(W_new - W) / np.linalg.norm(W) print turn, diff if diff < 0.1: break W = W_new turn += 1 W_t = np.transpose(W) print W_t # 訓練データを描画 x1, x2 = np.transpose(np.array(cls1)) plot(x1, x2, 'rx') x1, x2 = np.transpose(np.array(cls2)) plot(x1, x2, 'g+') x1, x2 = np.transpose(np.array(cls3)) plot(x1, x2, 'bo') # クラス1とクラス2の間の識別境界を描画 x1 = np.linspace(-6, 6, 1000) x2 = [f(x, W_t, 0, 1) for x in x1] plot(x1, x2, 'r-') # クラス2とクラス3の間の識別境界を描画 x1 = np.linspace(-6, 6, 1000) x2 = [f(x, W_t, 1, 2) for x in x1] plot(x1, x2, 'b-') xlim(-6, 6) ylim(-6, 6) show()
結果は、下のようになり、図4.5のロジスティック回帰の3クラス分類の結果と同じになりました。