🧠

【Python】ロジスティック回帰を使ってクラスの確率を予測するモデルの構築Part3

2024/04/07に公開

scikit-learn

この記事は株式会社インプレスの「Python機械学習プログラミング Pytorch&scilit-learn編」を読んで、私が学習したことをまとめています。リンク一覧はこちらから。今回は3章3節の「ロジスティック回帰を使ってクラスの確率を予測するモデルの構築」を読んで学んだことをまとめていきます。
※まとめている中で、思った以上にボリュームがあったので、細かく記事を分けています。今回は「3.3.3　ADALLINE実装をロジスティック回帰のアルゴリズムに変換する」を読んで学んだことをまとめています。
用語などの定義については【Python】ADALINE(フルバッチ勾配降下法)と学習の収束にまとめていますので、こちらをご確認ください。

Part1でロジスティック回帰の活性化関数であるシグモイド関数 $\sigma_3$ について学び、それに伴い閾値関数を再定義しました。そしてPart2ではパラメータを学習するための損失関数 $L_2$ として、対数尤度関数について学びました。
今回はいよいよロジスティック回帰の実装を行っていきたいと思います。

ライブラリのインポート

import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn import datasets
from sklearn.metrics import accuracy_score            # 正解率
from sklearn.linear_model import Perceptron           # パーセプトロンモデル
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split  # データセットの分割
from matplotlib.colors import ListedColormap

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

class LogisticRegressionGD:
    """
    パラメータ
    ----------------------------------------------------
    eta          : float.学習率
    n_iter       : int.訓練データの訓練回数
    random_state : int.重みを初期化するための乱数シード
    ----------------------------------------------------

    データ属性
    ----------------------------------------------------
    w_      : 1次元配列.適合後の重み
    b_      : スカラー.適合後のバイアス
    losses_ : リスト.各エポックでの損失関数の値
    ----------------------------------------------------
    """

    def __init__(self, eta = 0.01, n_iter = 50, random_state = 1):
        self.eta          = eta             # 学習率の初期化
        self.n_iter       = n_iter          # 訓練回数の初期化
        self.random_state = random_state    # 乱数シードを設定
    
    def fit(self, X, y):
        """
        パラメータ
        -----------------------------------------------------------------------------------------------
        X : {配列のようなデータ構造}. shape = [n_examples(訓練データの個数), n_features(特徴量の個数)]
        y : {配列のようなデータ構造}. shape = [n_examples(訓練データの個数)]
        -----------------------------------------------------------------------------------------------
        """
        rgen         = np.random.RandomState(self.random_state)
        self.w_      = rgen.normal(loc   = 0.0,
                                   scale = 0.01,
                                   size  = X.shape[1]
                                  )
        self.b_      = np.float_(0.)
        self.losses_ = []

        # 訓練回数分まで訓練データを反復処理
        for i in range(self.n_iter):
            net_input = self.net_input(X)
            output    = self.activation(net_input)
            errors    = y - output                    # 誤差の計算

            """ADALINEの重み、バイアス更新と同じ"""
            # 重み(w_)の更新
            self.w_ += self.eta * (2.0 * X.T.dot(errors) / X.shape[0])
            # バイアス(b_)の更新
            self.b_ += self.eta * (2.0 * errors.mean())
            """---------------------------------"""

            # 損失関数L_2の計算
            loss = (-y.dot(np.log(output)) - ((1 - y).dot(np.log(1 - output))) / X.shape[0])

            # 損失値を格納
            self.losses_.append(loss)
        
        return self
    
    # 総入力を計算
    def net_input(self, X):
        return np.dot(X, self.w_) + self.b_
    
    # ロジスティックシグモイド活性化関数を計算
    def activation(self, z):
        return 1. / (1. + np.exp(-np.clip(z, -250, 250)))
    
    # 1ステップ後のクラスラベルを返す
    def predict(self, X):
        return np.where(self.activation(self.net_input(X)) >= 0.5, 1, 0)

このロジスティック回帰のコードについてですが、ADALINEの重みとバイアスの更新定義と、ロジスティック回帰の重みとバイアスの更新定義は同じであることが知られているようです。まだ自分で確認が終わっていないので、「らしい」にとどめておきます。証明が完了し次第追記します。

それでは次に、データセットを用いてモデルの訓練を行っていきます。

データセットの準備

データセットの準備を行います。

# Irisデータセットをロード
iris = datasets.load_iris()

# データセットの確認
print(iris)

Irisデータセットの準備

# 特徴量の取得
# 150個のデータの2, 3行目(花びらの長さと花びらの幅)を抽出
X = iris.data[:, [2, 3]]

# 目的変数の取得
y = iris.target
# 目的変数のunique値（以後正解ラベルという）を確認する
print('Class labels:', np.unique(y))

実行結果

Class labels: [0 1 2]

訓練データセットとテストデータセットに分割する

# データセットの分割
X_train, X_test, y_train, y_test =\
    train_test_split(X, y,
                     test_size    = 0.3,    # 全体の30%をテストデータとする
                     random_state = 1,      # 乱数を固定
                     stratify     = y       # 層化抽出を実行
                    )

# 層化抽出を指定したため、正解ラベルのバランスを維持してデータセットを分割できる
print('Labels counts in y:', np.bincount(y))
print('Labels counts in y_train:', np.bincount(y_train))
print('Labels counts in y_test:', np.bincount(y_test))

実行結果

Labels counts in y: [50 50 50]
Labels counts in y_train: [35 35 35]
Labels counts in y_test: [15 15 15]

訓練データセットを標準化する

# インスタンス化
sc = StandardScaler()

# 学習し、特徴量毎の平均値と標準偏差を推定する
sc.fit(X_train)

# 推定した平均値と標準偏差を用いて標準化を行う
X_train_std = sc.transform(X_train)
X_test_std  = sc.transform(X_test)

ロジスティック回帰モデルによる学習を行う

品種としてIris-Setosa(クラス0)とIris-Versicolor(クラス1)のみを抽出します。

X_train_01_subset = X_train_std[(y_train == 0) | (y_train == 1)]
y_train_01_subset = y_train[(y_train == 0) | (y_train == 1)]

# 学習率0.3でロジスティック回帰のインスタンスを生成
lrgd = LogisticRegressionGD(eta = 0.3, n_iter = 1000, random_state = 1)

# 学習を実施
lrgd.fit(X_train_01_subset, y_train_01_subset)

決定境界を可視化するための関数の実装

def plot_decision_regions(X, y, classifier, test_idx = None, resolution = 0.02):

    """マーカーとカラーマップの準備"""
    markers = ('o', 's', '^', 'v', '<')
    colors  = ('#3F9877',                  # ジェードグリーン
               '#003F8E',                  # インクブルー
               '#EA5506',                  # 赤橙
               'gray',
               'cyan'
              )
    cmap    = ListedColormap(colors[: len(np.unique(y))])

    """グリッドポイント（格子点）の生成"""    
    # x軸方向の最小値、最大値を定義
    x_min = X[:, 0].min() - 1
    x_max = X[:, 0].max() + 1
    # y軸方向の最小値、最大値を定義
    y_min = X[:, 1].min() - 1
    y_max = X[:, 1].max() + 1

    # 格子点の生成
    xx, yy = np.meshgrid(np.arange(x_min, x_max, resolution),
                         np.arange(y_min, y_max, resolution)
                        )
    # 確認用
    #print(xx1)
    #print(xx2)

    """各特徴量を1次元配列に変換(ravel())して予測を実行"""
    # つまり2つの特徴量から0と予測された格子点と1と予測された値が格子点ごとにlabに格納される
    lab = classifier.predict(np.array([xx.ravel(), yy.ravel()]).T)

    # 確認用
    #print(np.array([xx1.ravel(), xx2.ravel()]).T)
    #print(lab) # 格子点の並びで0,1が格納されているが、この時点では1次元配列なので変換が必要

    """予測結果の元のグリッドポイント(格子点)のデータサイズに変換"""
    lab = lab.reshape(xx.shape)

    """グリッドポイントの等高線のプロット"""
    plt.contourf(xx, yy, lab,
                 alpha = 0.3,                    # 透過度を指定
                 cmap  = cmap
                 )
    plt.xlim(xx.min(), xx.max())
    plt.ylim(yy.min(), yy.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     = f'Class {cl}',
                    edgecolor = 'black'
                    )
        
    """テストデータ点を目立たせる（点を〇で表示）"""
    if test_idx:    # ここはbool(test_idx)と同義。つまりTrueを返す
        # すべてのデータ点を描画
        X_test = X[test_idx, :]
        y_test = y[test_idx]

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

決定領域を可視化

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

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

# グラフを保存
plt.savefig('fig3-5.png')

plt.show()