多項ロジットの理論と実装

問題設定
訓練データ(x_n, y_n), x_n ∈ R ^D, y_n ∈ {1,2,...C}, n=1,2,...,Nを使って、マルチクラス識別問題を解くモデルを考えます。

多項ロジットのモデル
多項ロジットでは以下のモデルを使います。

分子は訓練データxの特徴量の線形結合をロジスティック関数に入れたものです。
分子は非負の値を取ります。

上の式はxがクラスyに属する確率を表します。
よって、xがクラス1に属する確率、xがクラス2に属する確率、...をすべて足すと1.0になっていて欲しいです。この条件を満たすように分母の値で割って正規化を行います。

パラメータの学習
多項ロジットのパラメータは、

です。ここでCはクラスの個数、Dはデータの次元数です。
つまりクラスの数だけ重みベクトルが存在します。

これらのパラメータを学習するために、以下の対数尤度関数を最大化します。

対数尤度関数の偏微分を計算します。
まず第一項の偏微分は以下のようになります。

[]は中の条件が真のときは1、偽のときは0になるような関数です。

第二項の偏微分は以下のようになります。

よって全体の偏微分は、

となります。

偏微分が求まったので、θを適当な値に初期化して、勾配法で対数尤度関数を最適化することができます。

実装
Pythonで実装しました。

import numpy as np

class MultinomialLogit:

 def __init__(self):
  self.C = 0
  self.N = 0
  self.D = 0
  self.theta = None

 def prob(self, X, y):
  return np.dot(self.theta[y], X) / sum(np.dot(self.theta[c], X) for c in range(self.C))

 def logLiklihood(self, X, y):
  return sum(np.log(self.prob(X[n], y[n])) for n in range(self.N))
 
 def fit(self, X, y, iter=200):
  self.C = len(set(y))
  self.N = len(X)
  self.D = len(X[0])+1   # bias term is added
  self.theta = np.random.rand(self.C, self.D);

  bias = np.ones((len(X), 1))
  X = np.hstack((X, bias))

  def delta(a, b):
   return 1.0 if a == b else 0.0
  
  for i in range(iter):
   for c in range(self.C):
    gradient = sum((delta(c, y[n]) - self.prob(X[n], c)) * X[n] for n in range(self.N)) 
    self.theta[c] += gradient
   print ("i=%d %f" % (i, self.logLiklihood(X, y)))
    
 def predict(self, X):
  bias = np.ones((len(X), 1))
  X = np.hstack((X, bias))
  return [np.argmax([self.prob(x, c) for c in range(self.C)]) for x in X]

 def score(self, X, y):
  Ypred = self.predict(X)
  comp = Ypred == y
  return sum(comp) / len(comp)

Irisで実験
Irisのデータを使って実験してみました。
テスト識別率は80%後半程度です。sklearnのLogisticRegressionだと100%が出るので全然ダメな感じですが、

• iterationごとに対数尤度関数が上昇すること
• 対数尤度関数が上がるとテスト識別率も上がること

を確認できたのでそれっぽい実装は出来てそうです。学習率を設定したり、 ヘッセ行列を計算してニュートン法系の最適化を行ったり、確率的勾配法を使ったりなど改善の余地はたくさんありそうです。

if __name__ == '__main__':
 from sklearn import datasets
 from sklearn import cross_validation
 
 iris = datasets.load_iris()
 X = iris.data
 y = iris.target

 Xtrain, Xtest, Ytrain, Ytest = cross_validation.train_test_split(X, y, test_size=0.3)
 logit = MultinomialLogit()
 logit.fit(Xtrain, Ytrain, iter=1000)
 print (logit.score(Xtest, Ytest))