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

ロジスティック回帰とは

ロジスティック回帰は、主に分類問題に用いられる統計学・機械学習の手法です。入力変数（特徴量）から出力ラベル（クラス）への所属確率を予測するため、結果を確率として解釈できます。特に、二値分類（例：スパムか否か、病気の有無など）でよく使われますが、適切な拡張により複数クラスに対応することも可能です。

基本的な考え方

• 目的：入力データに基づいて、対象が特定のクラスに属する確率を出力する。
• 特徴：線形モデルの枠組みを用いつつ、出力を0〜1の範囲に収めるために非線形のシグモイド関数を適用する。
• 応用例：医療診断、マーケティング分析、金融リスク評価など。

数学的な背景

シグモイド関数

ロジスティック回帰では、シグモイド関数（またはロジスティック関数）を用いて、実数値の線形結合を確率に変換します。

シグモイド関数の定義は次の通りです。

$$\sigma(z) = \frac{1}{1 + e^{-z}}$$

• z は線形結合の結果（後述）で、シグモイド関数は入力を「S字状」の曲線に写像し、常に0から1の間の値を返します。

モデルの表現

1. 線形結合 入力ベクトル \(\mathbf{x} = (x_1, x_2, \dots, x_n)\) に対し、パラメータ \(\beta = (\beta_0, \beta_1, \dots, \beta_n)\) を用いて以下のように線形結合を行います。 $$z = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \dots + \beta_n x_n$$
2. 確率の計算 線形結合 z にシグモイド関数を適用して、クラス1に属する確率を求めます。 $$P(y = 1 \mid \mathbf{x}) = \sigma(z) = \frac{1}{1 + e^{-z}}$$ これにより、出力は確率として解釈可能となります。

コスト関数

ロジスティック回帰では、モデルのパラメータを最適化するために交差エントロピー誤差（対数損失）をコスト関数として採用します。

サンプル数を m とすると、コスト関数 J(\beta) は次のように定義されます。

$$J(\beta) = -\frac{1}{m} \sum_{i=1}^{m} \left[ y^{(i)} \log\left(\sigma(z^{(i)})\right) + (1 – y^{(i)}) \log\left(1 – \sigma(z^{(i)})\right) \right]$$

• \(y^{(i)}\) は正解ラベル（0または1）
• \(\sigma(z^{(i)})\) はサンプル i に対する予測確率

この損失関数の値を最小化することで、最も適切なパラメータ \(\beta\) を求めます。

パラメータ推定

• 最適化手法：一般に、勾配降下法や確率的勾配降下法（SGD）を用いてコスト関数を最小化し、パラメータを学習します。
• 収束条件：損失関数の変化が小さくなった時点、またはあらかじめ決めた反復回数に達した場合に学習を終了します。

利点と注意点

利点

• 解釈性が高い：各特徴量の重み（\(\beta_i\)）から、特徴が分類に与える影響が理解できます。
• 確率出力：出力が0〜1の範囲に収まり、確率として解釈できるため、意思決定に役立ちます。
• 計算効率が良い：パラメータが少ない場合、学習や推論が高速です。

注意点

• 線形決定境界：基本的には線形な境界しか表現できないため、複雑なデータに対しては性能が劣る場合があります。
• 特徴量の前処理：入力データのスケールが大きく異なる場合、学習に悪影響を及ぼすことがあるため、標準化などの前処理が重要です。
• 多重共線性：特徴量間の相関が強いと、パラメータの推定が不安定になる可能性があります。

多項ロジスティック回帰

複数クラスを識別する必要がある場合、多項ロジスティック回帰（ソフトマックス回帰）が用いられます。

• 仕組み：各クラスに対して線形結合を計算し、ソフトマックス関数によって各クラスに属する確率を出力します。
• ソフトマックス関数の定義： $$P(y = k \mid \mathbf{x}) = \frac{e^{z_k}}{\sum_{j} e^{z_j}}$$ ここで、\(z_k\) はクラス k に対する線形結合の結果です。

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

では、ここまで学習した理論をプログラムに落とし込んでみましょう。ここではGoogle colaboratoryの環境で実行していきます。

必要なライブラリをインポートしておきます。

from pathlib import Path

import numpy as np
from sklearn.datasets import fetch_openml
from sklearn.metrics import accuracy_score
from sklearn.model_selection import train_test_split
from sklearn.utils import shuffle

DATA_DIR = "/content/drive/MyDrive" #ここにデータがダウンロードされます

• pathlib.Path: ファイルやディレクトリのパスをオブジェクト指向で操作するために使用します。
• numpy (as np): 高速な数値計算、特に多次元配列（ベクトル、行列など）の操作を行うために使用します。
• sklearn.datasets.fetch_openml: OpenML.orgから機械学習用データセットをダウンロードするために使用します。
• sklearn.metrics.accuracy_score: 機械学習モデルの分類精度（正解率）を計算するために使用します。
• sklearn.model_selection.train_test_split: データセットを訓練用とテスト用に分割するために使用します。
• sklearn.utils.shuffle: データセットや配列の要素をランダムに並び替える（シャッフルする）ために使用します。

シグモイド関数の実装

def sigmoid_test(x):
    return 1 / (1 + np.exp(- x))

def sigmoid(x):
    return 0.5 * (1 + np.tanh(x / 2))

本来であればsigmoid_testのように書くのが公式通りですがこれを機会で実行するとOverflowという現象が起こってしまいます（桁が小さくなりすぎて値が消失してしまう現象が起こる）。

• np.tanh: NumPyライブラリに含まれる双曲線正接関数（ハイパボリックタンジェント）を計算する関数です。
• 数式: tanh(x) = (e^x - e^-x) / (e^x + e^-x) で定義されます。
• 特徴:
  • 入力値 x に対して、出力値は必ず -1から1の範囲 に収まります。
  • 形状がシグモイド関数 (1 / (1 + exp(-x))) に似ており、同じく活性化関数などとして利用されます。
• コードでの利用:
  • np.tanh(x * 0.5) * 0.5 + 0.5 という計算を行うことで、出力範囲を 0から1の範囲 に調整し、標準的なシグモイド関数と同じ出力範囲になるようにしています。
  • コメントにあるように、np.exp(-x) を直接計算すると x の値によっては指数関数の計算でオーバーフロー（非常に大きな値になりすぎて計算できなくなること）が発生する可能性があります。np.tanh は内部的に数値的な安定性を考慮した実装になっているため、オーバーフローのリスクを低減できます

データセットの設定と重みづけ

# ORのデータセット
x_train_or = np.array([[0, 1], [1, 0], [0, 0], [1, 1]])
t_train_or = np.array([[1], [1], [0], [1]])
x_valid_or, t_valid_or = x_train_or, t_train_or
x_test_or, t_test_or = x_train_or, t_train_or

# 重み (入力の次元数: 2, 出力の次元数: 1)
W_or = np.random.uniform(low=-0.08, high=0.08, size=(2, 1)).astype("float32")
b_or = np.zeros(shape=(1,)).astype("float32")

• 目的: ORゲートを機械学習モデルで学習させるための準備。
• ライブラリ: 数値計算ライブラリ NumPy (np) を使用。
• データセット:
  • x_train_or: ORゲートの入力データ（全4パターン）。
  • t_train_or: x_train_orに対応する正解出力。
  • x_valid_or, t_valid_or: 検証用データ（訓練データと同じ）。
  • x_test_or, t_test_or: テスト用データ（訓練データと同じ）。
• モデルパラメータ初期化:
  • W_or: 重み行列（入力次元2、出力次元1）。-0.08から0.08の範囲の乱数で初期化。
  • b_or: バイアス項（出力次元1）。ゼロで初期化。
  • データ型: 全てのNumPy配列をfloat32（32ビット浮動小数点数）に指定。

train関数とvalid関数を実装

# Numpyを用いてLogの計算結果が0になるのを防ぐ
def np_log(x):
    return np.log(np.clip(a=x, a_min=1e-10, a_max=1e1000))

def train_or(x, t, eps=1.0):
    """
    :param x: np.ndarray, 入力データ, shape=(batch_size, 入力の次元数)
    :param t: np.ndarray, 教師ラベル, shape=(batch_size, 出力の次元数)
    :param eps: float, 学習率
    """
    global W_or, b_or

    batch_size = x.shape[0]

    # 順伝播
    y = sigmoid(np.matmul(x, W_or) + b_or)  # shape: (batch_size, 出力の次元数)

    # 逆伝播
    cost = (-t * np_log(y) - (1 - t) * np_log(1 - y)).mean()
    delta = y - t  # shape: (batch_size, 出力の次元数)

    # パラメータの更新
    dW = np.matmul(x.T, delta) / batch_size  # shape: (入力の次元数, 出力の次元数)
    db = np.matmul(np.ones(shape=(batch_size,)), delta) / batch_size  # shape: (出力の次元数,)
    W_or -= eps * dW
    b_or -= eps * db

    return cost


def valid_or(x, t):
    y = sigmoid(np.matmul(x, W_or) + b_or)
    cost = (-t * np_log(y) - (1 - t) * np_log(1 - y)).mean()
    return cost, y

• np_log(x) 関数:
  • NumPyの log 関数のカスタム版。
  • 入力 x の値がゼロや非常に小さい値になるのを防ぐため、計算前に値を一定範囲内（1e-10 ～ 1e1000）にクリップ（制限）する。
  • 目的: log(0) によるエラーや計算の不安定性を回避するため。
• train_or(x, t, eps=1.0) 関数:
  • ORゲートモデルの1ステップ分の訓練を実行する関数。
  • 入力 x と正解 t を受け取り、モデルのパラメータ (W_or, b_or) を更新する。
  • 処理フロー:
    1. 順伝播: 入力 x と現在の W_or, b_or から、シグモイド関数 (sigmoid) を使って予測値 y を計算。
    2. 損失計算: 予測値 y と正解 t から交差エントロピー誤差 (cost) を計算（np_log を使用）。
    3. 逆伝播 (勾配計算): 損失を最小化するための W_or と b_or の勾配 (dW, db) を計算。
    4. パラメータ更新: 計算した勾配と学習率 eps を使って W_or と b_or を更新（勾配降下法）。
  • 計算された損失 cost を返す。
  • グローバル変数 W_or, b_or を直接変更する。
• valid_or(x, t) 関数:
  • 現在のモデルの性能を検証する関数（パラメータは更新しない）。
  • 入力 x と正解 t を受け取る。
  • 処理フロー:
    1. 順伝播: 入力 x と現在の W_or, b_or から予測値 y を計算。
    2. 損失計算: 予測値 y と正解 t から交差エントロピー誤差 cost を計算。
  • 計算された損失 cost とモデルの出力 y を返す。

For文で学習を行わせる

for epoch in range(1000):
    # online学習
    x_train_or, t_train_or = shuffle(x_train_or, t_train_or)
    for x, t in zip(x_train_or, t_train_or):
        cost = train_or(x[None, :], t[None, :])
    cost, y_pred = valid_or(x_valid_or, t_valid_or)

print(y_pred)

• 目的: 定義済みの train_or 関数と valid_or 関数を使って、ORゲートモデルを実際に訓練するループ。
• ループ:
  • for epoch in range(1000): で、訓練プロセスを 1000回繰り返す（1000エポック）。
• エポック内の処理:
  1. データシャッフル:
     • x_train_or, t_train_or = shuffle(x_train_or, t_train_or): 各エポックの開始時に、訓練データの順番をランダムに入れ替える（shuffle 関数の実装が必要）。学習の偏りを防ぐため。
  2. オンライン学習 (1サンプルずつの学習):
     • for x, t in zip(x_train_or, t_train_or):: シャッフルされた訓練データから、入力 x と正解 t のペアを1つずつ取り出すループ。
     • cost = train_or(x[None, :], t[None, :]): 取り出した1つのデータサンプルで train_or 関数を呼び出し、モデルの重みとバイアスを更新する。
       • x[None, :] は x を1サンプルだけのバッチ（2次元配列）に変換する。
  3. 検証:
     • cost, y_pred = valid_or(x_valid_or, t_valid_or): 1エポック（全ての訓練データを1回使用）が終了した後、検証データを使って現在のモデルの損失 (cost) と予測値 (y_pred) を計算する。パラメータ更新はしない。
• 最終結果の表示:
  • print(y_pred): 1000エポックの訓練が全て完了した後の、最後の検証ステップで得られた予測値 (y_pred) を表示する。これは、訓練されたモデルが検証データ（ORゲートの全入力パターン）に対してどのような出力をするかを示す。

コード全体のまとめ

from pathlib import Path

import numpy as np
from sklearn.datasets import fetch_openml
from sklearn.metrics import accuracy_score
from sklearn.model_selection import train_test_split
from sklearn.utils import shuffle

DATA_DIR = "/content/drive/MyDrive"

def sigmoid(x):
    return 0.5 * (1 + np.tanh(x / 2))

x_train_or = np.array([[0, 1], [1, 0], [0, 0], [1, 1]])
t_train_or = np.array([[1], [1], [0], [1]])
x_valid_or, t_valid_or = x_train_or, t_train_or
x_test_or, t_test_or = x_train_or, t_train_or

W_or = np.random.uniform(low=-0.08, high=0.08, size=(2, 1)).astype("float32")
b_or = np.zeros(shape=(1,)).astype("float32")

def np_log(x):
    return np.log(np.clip(a=x, a_min=1e-10, a_max=1e1000))

def train_or(x, t, eps=1.0):
    global W_or, b_or

    batch_size = x.shape[0]

    y = sigmoid(np.matmul(x, W_or) + b_or)

    cost = (-t * np_log(y) - (1 - t) * np_log(1 - y)).mean()
    delta = y - t

    dW = np.matmul(x.T, delta) / batch_size
    db = np.matmul(np.ones(shape=(batch_size,)), delta) / batch_size
    W_or -= eps * dW
    b_or -= eps * db

    return cost


def valid_or(x, t):
    y = sigmoid(np.matmul(x, W_or) + b_or)
    cost = (-t * np_log(y) - (1 - t) * np_log(1 - y)).mean()
    return cost, y

for epoch in range(1000):
    x_train_or, t_train_or = shuffle(x_train_or, t_train_or)
    for x, t in zip(x_train_or, t_train_or):
        cost = train_or(x[None, :], t[None, :])
    cost, y_pred = valid_or(x_valid_or, t_valid_or)

print(y_pred)

コードの説明

1. ライブラリのインポートと初期設定:
   • numpyなどの必要なライブラリをインポートします。
   • 乱数シードを設定して結果の再現性を確保します。
   • データディレクトリのパスを設定します（このコードでは直接使用されていませんが、設定されています）。
2. シグモイド関数:
   • sigmoid関数を定義します。これはロジスティック回帰の活性化関数として使われ、入力を0から1の間の確率値に変換します。ここでは数値的安定性のためにtanhを使って実装されています。
3. データセットの準備:
   • ORゲートの入力 (x_train_or) と対応する正しい出力 (t_train_or) を定義します。この例では、訓練データ、検証データ、テストデータはすべて同じものを使用しています。
4. モデルパラメータの初期化:
   • ロジスティック回帰モデルの重み (W_or) とバイアス (b_or) をランダムな値とゼロで初期化します。
5. 補助関数:
   • np_log関数を定義します。これは対数計算時に引数がゼロや非常に小さい値になるのを防ぐためのものです。
6. 学習関数 (train_or):
   • 入力データ (x) と教師ラベル (t) を受け取り、モデルの予測 (y) を計算します（順伝播）。
   • 予測と教師ラベルから損失（交差エントロピー誤差）を計算します。
   • 誤差に基づいて勾配 (dW, db) を計算します（逆伝播）。
   • 計算された勾配と学習率 (eps) を使って、重み (W_or) とバイアス (b_or) を更新します。
   • 計算されたコスト（損失）を返します。
7. 検証関数 (valid_or):
   • 入力データ (x) と教師ラベル (t) を受け取り、現在のモデルパラメータを使って予測 (y) と損失 (cost) を計算します。パラメータの更新は行いません。
   • 計算されたコストと予測値を返します。
8. 学習ループ:
   • 指定されたエポック数（ここでは1000回）だけ学習を繰り返します。
   • 各エポックの開始時に訓練データをシャッフルします。
   • 訓練データを1つずつ（オンライン学習）train_or関数に渡し、モデルのパラメータを更新します。
   • 各エポックの終わりにvalid_or関数を使って検証データに対する損失と予測値を計算します（ただし、このループ内ではcostとy_predは毎エポック上書きされます）。
9. 結果の表示:
   • 学習ループが終了した後、最後に計算された検証データに対する予測値 (y_pred) を表示します。これにより、学習後のモデルがORゲートの各入力に対してどのような確率を出力するかがわかります。

おまけ：多クラス分類用のソフトマックス回帰の実装

先ほどまでのロジスティック回帰分析では2値分類しかできないため複数のクラスに分類したい際にはソフトマックス回帰を用いる必要があります。

$$\mathrm{softmax}({\bf x})_k = \frac{\exp({\bf x}_k)}{\sum^K_{k’=1} \exp({\bf x}_{k’})} \hspace{10mm} \text{for} \, k=1,\ldots, K$$

\({\bf x}\)の要素の中に大きなものがある場合，発散してしまうため，下にあるように普通，

$${\bf x} \leftarrow {\bf x} – (\bf{x}_{max})$$

としてからsoftmax関数にかけることで，オーバーフローを防ぐ。

def softmax(x):
    x -= x.max(axis=1, keepdims=True)  # expのunderflow & overflowを防ぐ
    x_exp = np.exp(x)
    return x_exp / np.sum(x_exp, axis=1, keepdims=True)

• 目的: ソフトマックス関数を計算する。これは、入力された数値ベクトル（スコア）を、合計が1になる確率分布に変換するためによく使われる（特に出力層で）。
• 入力: NumPy配列 x。通常、各行が1つのサンプル、各列がそのサンプルの各クラスに対するスコアを表す（例: shape=(batch_size, num_classes)）。
• 処理フロー:
  1. オーバーフロー/アンダーフロー対策:
     • x -= x.max(axis=1, keepdims=True): 各行（サンプルごと）の最大値を、その行の全ての要素から引く。
     • 理由: 指数関数 np.exp() の計算時に、非常に大きな値（オーバーフロー）や非常に小さな値（アンダーフロー）が発生するのを防ぐため。この操作はソフトマックス関数の結果を変えない。
  2. 指数関数:
     • x_exp = np.exp(x): 値を引いた後の x の各要素に対して、指数関数 exp を適用する。
  3. 正規化:
     • return x_exp / np.sum(x_exp, axis=1, keepdims=True): 指数関数を適用した結果 (x_exp) を、各行（サンプルごと）の合計値で割る。
     • axis=1: 行ごとに合計を計算。
     • keepdims=True: 合計計算後も配列の次元数を維持し、ブロードキャストによる割り算を可能にする。
• 出力: 入力 x と同じ形状のNumPy配列。各行の要素の合計は1になり、各要素は0から1の間の値（確率）となる。