Python機械学習プログラミング第３章(3.3)

f:id:shintaro-football7:20191030160441j:plain
今回のテーマはロジスティック回帰です。
対数尤度関数（交差エントロピー型誤差関数）を最小化するように重みをもとめて、その重みで入力ベクトルを線形写像する。線形変換後のスカラーをロジスティック関数に代入して閾値0.5を超えていてればclass１に分類、という流れ。多クラス分類の場合は、クラスごとに線形写像をする重みｗが違くクラスの数だけ重みが存在しソフトマックス関数で事後確率が最大となるクラスに振り分ける。
(事後確率とは $P(C_k|x) = \pi_{k}(x) = \frac{exp(a_k)}{\sum_{j=1}^{K}exp(a_j)}$ )

[第2版]Python 機械学習プログラミング達人データサイエンティストによる理論と実践 (impress top gear)

作者: Sebastian Raschka,Vahid Mirjalili,福島真太朗,株式会社クイープ
出版社/メーカー: インプレス
発売日: 2018/03/16
メディア: 単行本（ソフトカバー）
この商品を含むブログ (2件) を見る

まずはロジスティック回帰を実装する。

import numpy as np
from sklearn import datasets
from sklearn.model_selection import train_test_split
from matplotlib.colors import ListedColormap
import matplotlib.pyplot as plt
from sklearn.linear_model import LogisticRegression
from sklearn .preprocessing import StandardScaler

class LogisticRegressionGD(object):
    def __init__(self,eta = 0.05,n_iter = 100,random_state = 1):
        self.eta = eta
        self.n_iter = n_iter
        self.random_state = random_state
    def fit(self,X,y):
        rgen = np.random.RandomState(self.random_state)
        self.w_ = rgen.normal(loc = 0,scale = 0.01, size=1+X.shape[1])
        self.cost_ = []
        for i in range(self.n_iter):
            net_input = self.net_input(X)
            output = self.activation(net_input)
            errors = (y - output)
            self.w_[1:] += self.eta * X.T.dot(errors)
            self.w_[0] += self.eta * errors.sum()

            cost = -y.dot(np.log(output)) - ((1-y).dot(np.log(1-output)))
            self.cost_.append(cost)
        return self
    def net_input(self,X):
        return np.dot(X,self.w_[1:]) + self.w_[0]
    def activation(self,z):
        return 1./(1 + np.exp(-np.clip(z,-250,250)))
    def predict(self,X):
        return np.where(self.net_input(X) >= 0.0,1,0)

def plot_decision_regions(X,y,classifier,test_idx = None,resolution = 0.02):
    markers = ('s','x','o', '^', 'v')
    colors = ('red', 'blue', 'lightgreen', 'gray', 'cyan')
    cmap = ListedColormap(colors[:len(np.unique(y))])

    x1_min,x1_max = X[:,0].min()-1,X[:,0].max()+1
    x2_min,x2_max = X[:,1].min()-1,X[:,1].max()+1
    xx1,xx2 = np.meshgrid(np.arange(x1_min,x1_max,resolution),
                          np.arange(x2_min,x2_max,resolution))
    Z = classifier.predict(np.array([xx1.ravel(),xx2.ravel()]).T)
    Z = Z.reshape(xx1.shape)
    plt.contourf(xx1,xx2,Z,alpha = 0.3,cmap = cmap)
    plt.xlim(xx1.min(), xx1.max())
    plt.ylim(xx2.min(), xx2.max())
    
    for idx, cl in enumerate(np.unique(y)):
        plt.scatter(x=X[y == cl, 0], 
                    y=X[y == cl, 1],
                    alpha=0.8, 
                    c=colors[idx],
                    marker=markers[idx], 
                    label=cl, 
                    edgecolor='black')
    if test_idx:
        X_test, y_test = X[test_idx, :], y[test_idx]

        plt.scatter(X_test[:, 0],
                    X_test[:, 1],
                    c='',
                    edgecolor='black',
                    alpha=1.0,
                    linewidth=1,
                    marker='o',
                    s=100, 
                    label='test set')

if __name__ == '__main__':
    iris = datasets.load_iris()
    X = iris.data[:,[2,3]] #3,4列目の特徴量を抽出
    y = iris.target #クラスラベルを取得
    X_train,X_test,y_train,y_test = train_test_split(X,y,test_size = 0.3,random_state = 1,stratify = y)

    X_train_01_subset = X_train[(y_train == 0)|(y_train==1) ]
    y_train_01_subset = y_train[(y_train == 0)|(y_train==1) ]
    lrgd = LogisticRegressionGD(eta=0.05, n_iter=1000, random_state=1)
    lrgd.fit(X_train_01_subset,y_train_01_subset)

    plot_decision_regions(X=X_train_01_subset, 
                      y=y_train_01_subset,
                      classifier=lrgd)

    plt.xlabel('petal length [standardized]')
    plt.ylabel('petal width [standardized]')
    plt.legend(loc='upper left')

    plt.tight_layout()
    plt.show()
    sc = StandardScaler()
    sc.fit(X_train)
    X_train_std = sc.transform(X_train)
    X_test_std = sc.transform(X_test)

    X_combined_std = np.vstack((X_train_std, X_test_std))#行方向に結合
    y_combined = np.hstack((y_train, y_test))#列方向に結合

    lr = LogisticRegression(C=100.0, random_state=1)
    lr.fit(X_train_std, y_train)

    plot_decision_regions(X_combined_std, y_combined,
                      classifier=lr, test_idx=range(105, 150))
    plt.xlabel('petal length [standardized]')
    plt.ylabel('petal width [standardized]')
    plt.legend(loc='upper left')
    plt.tight_layout()
    plt.show()
   
    print(lr.predict_proba(X_test_std[:3,:]))
    # [[3.17983737e-08 1.44886616e-01 8.55113353e-01][8.33962295e-01 1.66037705e-01 4.55557009e-12][8.48762934e-01 1.51237066e-01 4.63166788e-13]]
    weights, params = [], []
    for c in np.arange(-5, 5):
        lr = LogisticRegression(C=10.**c, random_state=1)
        lr.fit(X_train_std, y_train)
        weights.append(lr.coef_[1])
        params.append(10.**c)
    
    weights = np.array(weights)
    plt.plot(params, weights[:, 0],
         label='petal length')
    plt.plot(params, weights[:, 1], linestyle='--',
         label='petal width')
    plt.ylabel('weight coefficient')
    plt.xlabel('C')
    plt.legend(loc='upper left')
    plt.xscale('log')
    
    plt.show()