MLエンジニアへの道 #3 - ロジスティック回帰

ロジスティック回帰を学ぶ前に、前回の記事「MLエンジニアへの道 #2 - ロジスティック回帰の前提知識」で紹介した前提知識を理解しておく必要があります。まだ読んでいない場合は、この記事を進める前にそちらをチェックしてください。

ステップ 1. データ探索

ここでは、線形回帰で使用したのと同じデータセットであるアヤメデータセットを使用します。このデータセットには、アヤメ属の3つの種、つまりアヤメ・セトサ、アヤメ・ヴァージカラー、アヤメ・ヴァージニカの各50行のデータが含まれています。データセットを読み込むには、以下のコードを使用します。

# sklearnからデータセットをインポート
from sklearn.datasets import load_iris
 
# アヤメデータを読み込む
iris = load_iris()
 
# pd DataFramesを作成
# 以下は、種のカテゴリを数値から文字に変換するコードです
iris_df = pd.DataFrame(data= iris.data, columns= iris.feature_names)
target_df = pd.DataFrame(data= iris.target, columns= ['species'])
def converter(specie):
    if specie == 0:
        return 'setosa'
    elif specie == 1:
        return 'versicolor'
    else:
        return 'virginica'
target_df['species'] = target_df['species'].apply(converter)
 
# DataFramesを結合する
iris_df = pd.concat([iris_df, target_df], axis= 1)
 
# データを表示する
iris_df

上記のコードをGoogle Colabで実行すると、以下のようなものが表示されるはずです：

私たちのタスクは、セトサと他の種を分類することだと仮定します。前回と同様に、データを視覚化して情報を得ましょう。これを行うには、次のコードを使用します。

sns.pairplot(iris_df, hue= 'species')

異なる種に対して、すべての特徴量が異なる分布を持っていることが観察できます。この情報を使用して、どの変数をモデルに使用するかを選択できます。

ステップ 2. データ前処理

前のステップで得た情報に基づいて、データ前処理を進めましょう。すべての変数を使用すれば種をより良く分類できることがわかったので、データから列を削除する必要はありません。しかし、まだいくつかの作業があります。

まず、species変数を数値変数に戻す必要があります。その際、ヴァージカラーとヴァージニカを同じに扱う必要があります。また、species列を分離する必要があります。

iris_df.drop('species', axis= 1, inplace= True)
target_df = pd.DataFrame(columns= ['species'], data= iris.target)
 
# Set Setosa = 1, Others = 0
def converter(specie):
    if specie == 0:
        return 1
    else:
      return 0
 
X = np.array(iris_df)
y = np.array(target_df['species'].apply(converter))

説明変数にはX、分類したい結果変数にはyを使用するのが一般的です。さらに、データセットを訓練データとテストデータに分割する必要があります。

from sklearn.model_selection import train_test_split
 
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size= 0.33, random_state= 101)

上記のコードでは、データセットをランダムに分割して、訓練データとテストデータがそれぞれ元のデータセットの67%と33%になるようにしています。

ステップ 3. モデル

問題の理解

モデルを構築する前に、セトサとがく片の幅の関係をプロットして、どのようなものかを見てみましょう。

以下に示すように、データにロジスティック関数を適合させることができることがわかります。

データにロジスティック関数を適合させることができれば、がく片の幅の値を入力し、推定された確率をしきい値（通常は0.5）と比較することで、種がセトサであるかどうかを分類できます。では、どのようにしてロジスティック関数をデータに適合させるのでしょうか？

モデルの定義

ロジスティック関数の方程式を思い出してください：

f(x) = \frac{L}{1+e^{-k(x-x_0)}}

がく片の幅だけを使用する場合、 $x$ はがく片の幅です。種がセトサであるかどうかの確率を予測したいので、 $L$ を1に設定します。この方程式を整理すると、次のように書き換えることができます。

p(x) = \frac{1}{1+e^{-\beta_1x-\beta_0}}

ロジスティック関数についての議論を思い出すと、項 $\beta_1x + \beta_0$ はロジットに対応します。したがって、データにロジスティック関数を適合させるためには、単純に次のように適合させることができます。

log(\frac{p(x)}{1-p(x)}) = \beta_1x+\beta_0

つまり、がく片の幅とセトサであるかないかの確率に対応するロジットの線形回帰を単に行えばロジスティック関数を適合させることができるということです。複数の説明変数を使用する場合は、次のように簡単に適合させることができます。

log(\frac{p(x)}{1-p(x)}) = \beta_1sl+\beta_2sw+\beta_3pl+\beta_4pw+\beta_0

線形回帰を扱ったとき、最適化のためにMSEをコスト関数として使用しました。しかし、ここではMSEを使用することはできません。なぜなら、セトサとその他のラベルは0と1であり、そのロジットは $-\infty$ と $\infty$ だからです。どのセットの $\beta$ sでも、MSEは $\infty$ になってしまい使い物になりません。その代わりに、尤度を使用することができます。前回の記事で学んだように、最大尤度推定（MLE）を行って、データを生成するモデルに最も近いパラメータを得られます。では、この状況の尤度関数を決定してみましょう。

モデルが正しく分類したデータが多ければ多いほど、モデルの尤度は高くなるはずです。つまり、ラベルが1（ $y_i = 1$ ）とされるサンプルでは、ラベルが1であると予測される確率が高くなければならず（ $p(x_i)$ が高くなければならない）、ラベルが0（ $y_i = 0$ ）とされるサンプルでは、ラベルが1であると予測される確率が低くなければなりません（ $(1-p(x_i))$ が高くなる）。したがって、次のように尤度関数を定義できます。

L(\beta) = \prod_{y_i=1}(p(x_i)) * \prod_{y_i=0}(1-p(x_i))

この式は次のように簡略化できます。

L(\beta) = \prod_{i=1}^{n}(p(x_i))^{y_i} * (1-p(x_i))^{1-y_i}

これは、サンプルのラベルが1の場合、 $1-y_i$ は0になり、全項が $p(x_i)$ になるからです。サンプルのラベルが0の場合、 $y_i$ は0になり、 $1-y_i$ は1になり、全項が $1-p(x_i)$ になります。したがって、モデルを評価するために、上記の尤度関数を使用できます。

しかし、最良のパラメータセットを見つけるために、コスト関数を $\beta$ sに関して偏微分する必要があります。上記の尤度関数ではそれが難しいため、通常は尤度の対数(log)を取ってコスト関数をさらに簡略化します。

log(L(\beta)) = \sum_{i=1}^{n}y_ilog(p(x_i))+(1-y_i)log(1-p(x_i))

また、通常はコスト関数は最大化ではなく最小化したいと考えます。対数尤度を最大化することは、負の対数尤度を最小化することと同じであるため、負の対数尤度をコスト関数として使用します。

J(\beta) = \sum_{i=1}^{n}-y_ilog(p(x_i))-(1-y_i)log(1-p(x_i))

上記の式はどこかで見覚えがあるでしょう。それもそのはず、上の式はクロスエントロピー関数と非常に似ています。実際、負の対数尤度の最初の項は、ラベルと予測確率のクロスエントロピーと同じです！

H(y, p) = -\sum_{i=1}^{n}y_ilog(p(x_i))

二つ目の項はもう一方のラベルのクロスエントロピーと捉えることができるため、負の対数尤度とクロスエントロピーが同一であることを確認できます。

学習メカニズム

最大尤度推定を行う方法や、コスト関数を最小化する方法はいくつかありますが、ここでは勾配降下法を使用しましょう。勾配降下法を実装するためには、コスト関数をパラメータ $\beta$ sに関して偏微分する必要があります。そのためにまず、 $\beta$ sを使って負の対数尤度を表しましょう。

J(\beta) = \sum_{i=1}^{n}-y_ilog(\frac{1}{1+e^{-\beta x_i}})-(1-y_i)log(1-\frac{1}{1+e^{-\beta x_i}}) \\ = \sum_{i=1}^{n}-y_ilog(\frac{1}{1+e^{-\beta x_i}})-(1-y_i)log(\frac{e^{-\beta x_i}}{1+e^{-\beta x_i}}) \\ = \sum_{i=1}^{n}y_i[log(\frac{e^{-\beta x_i}}{1+e^{-\beta x_i}})-log(\frac{1}{1+e^{-\beta x_i}})]-log(\frac{e^{-\beta x_i}}{1+e^{-\beta x_i}}) \\ = \sum_{i=1}^{n}y_i[log(e^{-\beta x_i})-log(1+e^{-\beta x_i})-log(1)+log(1+e^{-\beta x_i})]-log(\frac{e^{-\beta x_i}}{1+e^{-\beta x_i}}*\frac{e^{\beta x_i}}{e^{\beta x_i}}) \\ = \sum_{i=1}^{n}y_ilog(e^{-\beta x_i})-log(\frac{1}{1+e^{\beta x_i}}) \\ = \sum_{i=1}^{n}-y_i\beta x_i+log(1+e^{\beta x_i})

ここで、 $\beta x_i$ はロジットの線形方程式の略です。このコスト関数を傾きと切片に関してチェーンルール(連鎖律)を使って微分しましょう。傾きに対しては：

\frac{d}{d\beta} = \sum_{i=1}^{n}-y_ix_i+\frac{1}{1+e^{\beta x_i}}e^{\beta x_i}x_i \\ = \sum_{i=1}^{n}-y_ix_i+\frac{e^{\beta x_i}}{1+e^{\beta x_i}}x_i \\ = \sum_{i=1}^{n}-y_ix_i+\frac{1}{1+e^{-\beta x_i}}x_i \\ = \sum_{i=1}^{n}-y_ix_i+p(x_i)x_i \\ = \sum_{i=1}^{n}(p(x_i)-y_i)x_i

切片に対しては：

\frac{d}{d\beta_0} = \sum_{i=1}^{n}-y_i+\frac{1}{1+e^{\beta x_i}}e^{\beta x_i} \\ = \sum_{i=1}^{n}-y_i+\frac{e^{\beta x_i}}{1+e^{\beta x_i}} \\ = \sum_{i=1}^{n}-y_i+\frac{1}{1+e^{-\beta x_i}} \\ = \sum_{i=1}^{n}p(x_i)-y_i

上記から、切片と傾きの偏微分は、予測確率と実際の確率の差を合計し、予測確率と実際の確率の差に対応する $x$ を掛けたものの合計を取ることで簡単に計算できることがわかります。これにより、コードの実装が非常に簡単になります。

コードの実装

上記の情報を使用して、LogisticRegressionGDモデルを作成しましょう。

class LogisticRegressionGD():
  def __init__(self, lr=0.001):
    self.W = np.zeros(X.shape[1])
    self.b = 0
    self.lr = lr # Learning rate
    self.history = [] # History of loss
 
  def predict(self, X):
    log_odds = np.sum(self.W*X + self.b, axis=1)
    return 1 / (1 + np.exp(-log_odds))
 
  def fit(self, X, y, epochs=200):
    for i in range(epochs):
      pred = self.predict(X)
 
      self.history.append(log_loss(y, pred))
 
      diff = pred - y
      grad_W = np.sum(diff[:, np.newaxis]*X, axis=0)
      grad_b = np.sum(diff)
 
      self.W -= self.lr * grad_W
      self.b -= self.lr * grad_b
    return self.history

コードを詳しく見ると、LinearRegressionGDと非常に似ていることに気づくかもしれません。実際、異なるのはpredictメソッド、損失計算、および勾配だけで、他はすべて同じです。勾配計算の唯一の違いは、LogisticRegressionGDが勾配をデータポイントの数で割らないことのみです。驚きですね！

LogisticRegressionGDを初期化し、以下のように訓練データに適合させることができます。

lr = LogisticRegressionGD()
history = lr.fit(X_train, y_train, epochs=200)

各エポックにおける負の対数尤度損失を追跡することで、勾配降下法がどのように機能しているかをプロットできます。

plt.plot(history)
plt.title("Loss vs Epoch")
plt.ylabel("Loss")
plt.xlabel("Epoch")
plt.show()

上記のプロットからわかるように、損失は各エポックで徐々に減少しています。これは、勾配降下法が機能している良い兆候です。そろそろ記事が長くなってきたので、パイプラインの最終ステップである評価については、今後の記事に残しておきます。

シグモイド関数についてのトリビア

Waite, E.（2020）のビデオで、バイナリ分類にシグモイド関数を使用する理由について興味深い解釈が紹介されていました。まず、実用的な観点から、シグモイド関数が選ばれる理由として、微分が容易であり、勾配降下法に適していることが挙げられます。これは、上記の議論からも十分に納得できます。

次に、彼は、2つのクラスの2つの正規分布が同じ分散を持つ場合に、右側の分布からサンプルを引く確率分布を計算すると、シグモイド関数と同じ曲線が得られることを示しています。シグモイド関数に至るいくつかの数学的導出が示されているので、興味がある方は以下のリンクから動画をチェックしてみてください。この視点は考えたことがなく、個人的には、シグモイド関数を確率的観点から正当化できることが美しいと思います。

リソース

Addagatla, A. 2021. Maximum Likelihood Estimation in Logistic Regression YouTube.
Waite, E. 2020. Why Do We Use the Sigmoid Function for Binary Classification? YouTube.