少し考えてみましたが、コスト関数である対数尤度は活性化関数の出力が教師データに近づくにつれ、0に漸近していくようです。
計算機上ではにするとが－∞に発散してしまうので値を求めることができないですが、解析解としては0に漸近します。
したがって、の極限で、コスト関数は最小値0になります。
以下が、ロジスティック回帰を実装したスクリプトになります。
| # -*- coding: utf-8 -*- | |
| import sys | |
| import numpy as np | |
| import pandas as pd | |
| import matplotlib.pyplot as plt | |
| from matplotlib.colors import ListedColormap | |
| def main(): | |
| # ---アヤメデータの取得 | |
| df = pd.read_csv("https://archive.ics.uci.edu/ml/machine-learning-databases/iris/iris.data", header=None) | |
| y = df.iloc[0:100, 4].values | |
| y = np.where(y == "Iris-setosa", 0, 1) | |
| X = df.iloc[0:100, [0, 2]].values | |
| # ---トレーニングサンプルの標準化 | |
| Xstd = np.copy(X) | |
| Xstd[:, 0] = (X[:, 0] - X[:, 0].mean()) / X[:, 0].std() | |
| Xstd[:, 1] = (X[:, 1] - X[:, 1].mean()) / X[:, 1].std() | |
| # ---学習実施(学習率:0.01 最大エポック数:20 (特徴量を標準化,シャッフルオン,バッチ勾配降下法)) | |
| clf = LogisticRegression(eta=0.01, numEpoch=20, shuffle=True) | |
| clf.fit(Xstd, y) | |
| plotResult(clf,Xstd,y,"Standardization","On","BGD") | |
| # ---学習実施(学習率:0.2 最大エポック数:20 (特徴量を標準化,シャッフルオン,バッチ勾配降下法)) | |
| clf = LogisticRegression(eta=0.2, numEpoch=20, shuffle=True) | |
| clf.fit(Xstd, y) | |
| plotResult(clf,Xstd,y,"Standardization","On","BGD") | |
| # ---学習実施(学習率:0.01 最大エポック数:20 (特徴量を標準化,シャッフルオン,確率的勾配降下法)) | |
| clf = LogisticRegression(eta=0.01, numEpoch=20, shuffle=True) | |
| clf.fitSGD(Xstd, y) | |
| plotResult(clf, Xstd, y, "Standardization", "On", "SGD") | |
| # ---学習実施(学習率:0.2 最大エポック数:20 (特徴量を標準化,シャッフルオン,確率的勾配降下法)) | |
| clf = LogisticRegression(eta=0.2, numEpoch=20, shuffle=True) | |
| clf.fitSGD(Xstd, y) | |
| plotResult(clf, Xstd, y, "Standardization", "On", "SGD") | |
| # ---学習実施(学習率:0.01 最大エポック数:20 (特徴量を標準化,シャッフルオン,一発解法)) | |
| clf = LogisticRegression(eta=0.01, numEpoch=20, shuffle=True) | |
| clf.fitLS(Xstd, y) | |
| plotResult(clf, Xstd, y, "Standardization", "On", "LS") | |
| def plotResult(clf,X,y,trainSampleState,shuffeled,fitAlgorism): | |
| #---プロットのプロパティ | |
| markers = ('s', 'o', 'x', '^', 'v') | |
| colors = ('green', 'yellow','red', 'blue', 'lightgreen', 'gray', 'cyan') | |
| cmap = ListedColormap(colors[:len(np.unique(y))]) | |
| labels = ('setosa', 'versicolor') | |
| fig = plt.figure(figsize=(18, 10)) | |
| plt.clf() | |
| ax1 = fig.add_subplot(231) | |
| ax2 = fig.add_subplot(232) | |
| ax3 = fig.add_subplot(233) | |
| ax4 = fig.add_subplot(234) | |
| ax5 = fig.add_subplot(235) | |
| # ---エポック数と重みのプロット | |
| epochTimes = np.array(range(0, clf.numEpoch + 1)) | |
| clf.W_ = np.array(clf.W_) | |
| ax1.plot(epochTimes, clf.W_[:, 0].T[0], color="blue", label="w0") | |
| ax1.plot(epochTimes, clf.W_[:, 1].T[0], color="red", label="w2") | |
| ax1.plot(epochTimes, clf.W_[:, 2].T[0], color="green", label="w1") | |
| ax1.legend(loc="upper left") | |
| ax1.set_xlabel('Epochs') | |
| ax1.set_ylabel('Weight') | |
| ax1.set_title("Epochs and weight") | |
| ax1.set_xlim(0, clf.numEpoch) | |
| # ---ax1.set_xticks(range(0,numEpoch+1,1)) | |
| ax1.grid() | |
| # ---エポック数と正答率のプロット | |
| epochTimes = range(1, clf.numEpoch + 1) | |
| ax2.plot(epochTimes, np.array(clf.correctAnswerRateEachEpoch_) * 100, color="blue", marker="o", | |
| label="Correct answer rate") | |
| # ---ax2.legend(loc="upper left") | |
| ax2.set_title("Epochs and correct answer rate") | |
| ax2.set_xlabel('Epochs') | |
| ax2.set_ylabel('Correct answer rate[%]') | |
| ax2.set_xlim(0, clf.numEpoch) | |
| ax2.grid() | |
| # ---教師データのクラスラベルが1(バージカラー)であるである確率 | |
| numTrainSamples = len(y) | |
| trainSampeleIndex=range(numTrainSamples) | |
| ax3.plot(trainSampeleIndex, clf.probability(X)*100, color="blue", marker="o",label="Probability") | |
| ax3.set_title("TrainSamples and probability(label=1)") | |
| ax3.set_xlabel('TrainSamples') | |
| ax3.set_ylabel('Probability[%]') | |
| ax3.set_xlim(0, numTrainSamples) | |
| ax3.grid() | |
| # ---エポック数とコスト関数のプロット | |
| epochTimes = range(1, clf.numEpoch + 1) | |
| ax4.plot(epochTimes, np.array(clf.J_)[:, 0][:, 0], color="blue", marker="o", | |
| label="Cost function values") | |
| # ---ax3.legend(loc="upper left") | |
| ax4.set_title("Epochs and cost function values") | |
| ax4.set_xlabel('Epochs') | |
| ax4.set_ylabel('Cost function values') | |
| ax4.set_xlim(0, clf.numEpoch) | |
| ax4.grid() | |
| # ---ax3.set_xticks(range(0,numEpoch+1,1)) | |
| # ---決定領域のプロット | |
| x1_min, x1_max = X[:, 0].min() - 1, X[:, 0].max() + 1 | |
| x2_min, x2_max = X[:, 1].min() - 1, X[:, 1].max() + 1 | |
| dx = 0.02 | |
| X1 = np.arange(x1_min, x1_max, dx) | |
| X2 = np.arange(x2_min, x2_max, dx) | |
| X1, X2 = np.meshgrid(X1, X2) | |
| Z = clf.predict(np.array([X1.ravel(), X2.ravel()]).T) | |
| Z = Z.reshape(X1.shape) | |
| # ---決定領域 | |
| ax5.contourf(X1, X2, Z, alpha=0.5, cmap=cmap) | |
| ax5.set_xlim(X1.min(), X1.max()) | |
| ax5.set_ylim(X2.min(), X2.max()) | |
| # ---アヤメデータ | |
| for idx, cl in enumerate(np.unique(y)): | |
| ax5.scatter(x=X[y == cl, 0], y=X[y == cl, 1], | |
| alpha=1.0, c=cmap(idx), | |
| marker=markers[idx], label=labels[idx]) | |
| ax5.set_title("Decision regions") | |
| ax5.set_xlabel("Sepal length [cm]") # がく片の長さ | |
| ax5.set_ylabel("Petal length [cm]") # 花びらの長さ | |
| ax5.legend(loc="upper left") | |
| ax5.grid() | |
| pngFileName=u"ロジスティック回帰学習結果_トレーニングサンプル" +trainSampleState +u"_シャッフル"+shuffeled+u"_Fitアルゴリズム"+fitAlgorism+ u"_学習率"+str(clf.eta)+u"_最大エポック数"+str(clf.numEpoch)+ ".png" | |
| plt.savefig(pngFileName, dpi=300) | |
| class LogisticRegression(object): | |
| def __init__(self, eta=0.01, numEpoch=10, shuffle=True, random_state=None): | |
| self.eta = eta #学習率 | |
| self.numEpoch = numEpoch #最大エポック数 | |
| self.shuffle = shuffle #トレーニングサンプルをシャッフルする？ | |
| if random_state: | |
| np.random.seed(random_state) | |
| self.W_=[] #各エポックごとに重みを保存するリスト | |
| self.J_=[] #各エポックごとにコスト関数を保存するリスト | |
| self.correctAnswerRateEachEpoch_=[] #各エポックごとに正答率を保存するリスト | |
| def __actFunc(self, z): | |
| """活性化関数(__でプライベート関数)""" | |
| return 1.0 / (1.0 + np.exp(-z)) | |
| def __quantizer(self, phi): | |
| """量子化器(__でプライベート関数)""" | |
| return np.where(np.array(phi) >= 0.5, 1, 0) | |
| def __shuffle(self, X, y): | |
| """トレーニングサンプルのシャッフル(__でプライベート関数)""" | |
| r = np.random.permutation(len(y)) | |
| return X[r], y[r] | |
| def predict(self, X): | |
| """予測関数""" | |
| X = np.matrix(X) | |
| m, n = X.shape | |
| X = np.c_[np.matrix(np.ones((m, 1))), X] | |
| z=X*self.w_ | |
| phi=self.__actFunc(z) | |
| return self.__quantizer(phi) | |
| def probability(self, X): | |
| """教師データのクラスラベルが1である確率を求める関数""" | |
| X = np.matrix(X) | |
| m, n = X.shape | |
| X = np.c_[np.matrix(np.ones((m, 1))), X] | |
| z=X*self.w_ | |
| phi=self.__actFunc(z) | |
| return phi | |
| def fit(self, X, y): | |
| """バッチ勾配降下法による学習の実施""" | |
| # ---トレーニングサンプルのシャッフル | |
| if self.shuffle: | |
| X, y = self.__shuffle(X, y) | |
| # ---特徴行列 | |
| X = np.matrix(X) | |
| m, n = X.shape | |
| # ---しきい値用に一番左側の列に1を追加 | |
| X = np.c_[np.matrix(np.ones((m, 1))), X] | |
| # ---教師データベクトル | |
| y = np.matrix(y).T | |
| # ---重みベクトル(初期化) | |
| self.w_ =np.matrix(np.zeros((n+1,1))) | |
| # 重みを学習回数ごとに保存 | |
| self.W_.append(self.w_.tolist()) | |
| for indexEpoch in range(1,self.numEpoch+1): #エポックループ | |
| correctAnswer=[] #各トレーニングセットに対する教師データと予測値の正誤表をイニシャライズ | |
| # ---特徴行列と重みベクトルの掛け合わせ | |
| z = X * self.w_ | |
| # ---活性化関数出力ベクトル | |
| phi =self.__actFunc(z) | |
| # ---コスト関数 | |
| J = -y.T * np.log(phi) - (np.ones((m, 1)) - y).T * np.log(np.ones((m, 1)) - phi) | |
| # ---コスト関数の勾配 | |
| gradJ = -X.T * (y - phi) | |
| # ---重み更新 | |
| dw = - self.eta * gradJ | |
| self.w_ = self.w_ + dw | |
| # 重みをエポックごとに保存 | |
| self.W_.append(self.w_.tolist()) | |
| # コスト関数をエポックごとに保存 | |
| self.J_.append(J.tolist()) | |
| #--正答表 | |
| correctAnswer = np.where(np.array(y == self.__quantizer(phi)) == True, 1, 0) | |
| #各エポックごとに正誤表から正答率算出し保存 | |
| self.correctAnswerRateEachEpoch_.append(float(sum(correctAnswer))/len(correctAnswer)) | |
| def fitSGD(self, X, y): | |
| """確率的勾配降下法による学習の実施""" | |
| # ---トレーニングサンプルのシャッフル | |
| if self.shuffle: | |
| X, y = self.__shuffle(X, y) | |
| # ---特徴行列 | |
| X = np.matrix(X) | |
| m, n = X.shape | |
| # ---しきい値用に一番左側の列に1を追加 | |
| X = np.c_[np.matrix(np.ones((m, 1))), X] | |
| # ---教師データベクトル | |
| y = np.matrix(y).T | |
| # ---重みベクトル(初期化) | |
| self.w_ = np.matrix(np.zeros((n + 1, 1))) | |
| # 重みを学習回数ごとに保存 | |
| self.W_.append(self.w_.tolist()) | |
| for indexEpoch in range(1, self.numEpoch + 1): # エポックループ | |
| correctAnswer = [] # 各トレーニングセットに対する教師データと予測値の正誤表をイニシャライズ | |
| for xi, yi in zip(X, y): # トレーニングセットループ | |
| # ---特徴行列と重みベクトルの掛け合わせ | |
| zi = xi * self.w_ | |
| # ---活性化関数出力ベクトル | |
| phii = self.__actFunc(zi) | |
| # ---コスト関数の勾配 | |
| gradJi = -xi.T * (yi - phii) | |
| # ---重み更新 | |
| dwi = - self.eta * gradJi | |
| self.w_ = self.w_ + dwi | |
| # ---特徴行列と重みベクトルの掛け合わせ | |
| z = X * self.w_ | |
| # ---活性化関数出力ベクトル | |
| phi =self.__actFunc(z) | |
| # ---コスト関数 | |
| J = -y.T * np.log(phi) - (np.ones((m, 1)) - y).T * np.log(np.ones((m, 1)) - phi) | |
| # 重みをエポックごとに保存 | |
| self.W_.append(self.w_.tolist()) | |
| # コスト関数をエポックごとに保存 | |
| self.J_.append(J.tolist()) | |
| # --正答表 | |
| correctAnswer = np.where(np.array(y == self.__quantizer(phi)) == True, 1, 0) | |
| # 各エポックごとに正誤表から正答率算出し保存 | |
| self.correctAnswerRateEachEpoch_.append(float(sum(correctAnswer)) / len(correctAnswer)) | |
| def fitLS(self, X, y): | |
| """一発学習の実施""" | |
| # ---トレーニングサンプルのシャッフル | |
| if self.shuffle: | |
| X, y = self.__shuffle(X, y) | |
| # ---特徴行列 | |
| X = np.matrix(X) | |
| m, n = X.shape | |
| # ---しきい値用に一番左側の列に1を追加 | |
| X = np.c_[np.matrix(np.ones((m, 1))), X] | |
| # ---教師データベクトル | |
| y = np.matrix(y).T | |
| # ---重みベクトル(初期化) | |
| self.w_ = np.matrix(np.zeros((n + 1, 1))) | |
| # 重みを学習回数ごとに保存 | |
| self.W_.append(self.w_.tolist()) | |
| for indexEpoch in range(1, self.numEpoch + 1): # エポックループ | |
| correctAnswer = [] # 各トレーニングセットに対する教師データと予測値の正誤表をイニシャライズ | |
| # ---コスト関数を最小にする重みを一発で算出 | |
| dy = 1e-10 | |
| yHat = np.where(y == 1, 1 - dy, dy) #---log(0)が計算できないのでyをdyだけずらしてやる | |
| self.w_ = np.linalg.inv(X.T*X)*X.T*(np.log(yHat)-np.log(np.ones((m, 1)) - yHat)) | |
| # ---特徴行列と重みベクトルの掛け合わせ | |
| z = X * self.w_ | |
| # ---活性化関数出力ベクトル | |
| phi = self.__actFunc(z) | |
| # ---コスト関数 | |
| J = -y.T * np.log(phi) - (np.ones((m, 1)) - y).T * np.log(np.ones((m, 1)) - phi) | |
| # 重みをエポックごとに保存 | |
| self.W_.append(self.w_.tolist()) | |
| # コスト関数をエポックごとに保存 | |
| self.J_.append(J.tolist()) | |
| # --正答表 | |
| correctAnswer = np.where(np.array(y == self.__quantizer(phi)) == True, 1, 0) | |
| # 各エポックごとに正誤表から正答率算出し保存 | |
| self.correctAnswerRateEachEpoch_.append(float(sum(correctAnswer)) / len(correctAnswer)) | |
| if __name__ == "__main__": | |
| main() | |
ベルヌーイ分布を使うため、教師データのクラスラベルは0と1にしてあります。(セトサが0でバージカラーが1です。)
ロジスティック回帰では、トレーニングサンプルを与えたときの出力として得られるクラスラベルの確率を求めることができるので、probabilityメソットを追加しました。
このスクリプトを実行すると、学習率と勾配効果法を変更した4つの学習結果グラフを出力します。(最大エポック数:20、特徴量の標準化、トレーニングデータのシャッフルは全結果で共通です。)
また、最後に一発解法で求めた結果を出力してみました。
学習率:0.01 バッチ勾配降下法の結果
右上にトレーニングサンプルを与えたときの出力として得られるクラスラベルが1(バージカラー)である確率グラフを追加しました。
最初の50個がセトサで後の50個がバージカラーなので、それなりの確率が出力されています。
しかし、20回のエポックループでも学習が不十分ですね。2こ外してます。
学習率:0.2 バッチ勾配降下法の結果
学習率を0.2にしてやると、正答率が100%になりました。
コスト関数が単調減少なので、学習率を大きくしてやると発散せずに学習スピードが上がるんですね。
しかしこれ以上学習率を上げてやると、コスト関数のの計算で浮動小数点演算の限界を超えるようで、nan出力されてしまいます。
決定境界がADALINEの時と異なっていますね。
学習率:0.01 確率的勾配降下法の結果
確率的勾配降下法を使うと、1エポック目の正答率が96%になり、バッチ勾配降下法の時より良くなっています。
しかし、20回のエポック後では、結果がほとんど変わらないですね。
学習率:0.2 確率的勾配降下法の結果
学習率を0.2では、2回目のエポックで正答率が100%になりました。
上の4つの中では、一番いい結果ですね。
一発解法の結果
ADALINEの時はコスト関数が誤差二乗和だったので最小二乗法と呼べましたが。。今回はなんて言うんでしょうね？
一発解法で求めた重みに、をそのまま代入すると計算できないので、少しだけずらしたを入れてみた結果です。
やはり一発解法は一発でしっかり学習してくれていますね。
これ以上ずらし量を小さくすると浮動小数点演算の限界に達してしまうので、ここがプログラムの限界です。
結局ロジスティック回帰の利点は、トレーニングサンプルを入力したときのクラスラベルを確率で表すことができる事だけなんでしょうかね。。
学習効率はいまいちな感がありますが、テキストには結構実用されているような事が書いてあったので、何か特化した使い道があるんでしょうか？
まあ取りあえずは、ロジスティック回帰の概念が理解できてよかったです。
