一般化線形モデル(GLM)の理論と実装②【ロジスティック回帰_後編】｜スクラッチ実装で理解する機械学習アルゴリズム #3

f:id:lib-arts:20190519144320p:plain

連載の経緯の詳細は#1でまとめましたが、本シリーズではあえてスクラッチ実装を元に機械学習のアルゴリズムを実装していくことで、アルゴリズムの概要を掴んだり理論の流れを掴んだりできるようにできればと思います。

実装のほとんどが車輪の再発明に近くなりますが、ベーシックなアルゴリズムをあえて一から自分で追ってみるというのは引き出しを増やすという意味で非常に有意義です。
#2では一般化線形モデル(GLM; Generalized Linear Model)の具体例の一つとしてロジスティック回帰について実装を通して理解を深める手助けとしました。

#2のスコア法以外でもロジスティック回帰の学習を取り扱えればということで、#3では勾配法を用いたロジスティック回帰のパラメータの学習について取り扱っていければと思います。
以下目次になります。

1. ベルヌーイ分布における最尤法と交差エントロピー(cross entropy)
2. #2の例題再訪＆ロジスティック回帰における勾配法
3. 実装で理解するロジスティック回帰（勾配法）
4. まとめ

1. ベルヌーイ分布における最尤法と交差エントロピー(cross entropy)
1節では問題背景としてベルヌーイ分布における最尤法と交差エントロピーについてまとめていきます。まず最尤法ですが、パラメトリックな確率分布にしたがって得られたデータが生成されたと考えた上で、同時確率を算出しそれを尤度とみなすことで尤度を最大化するようにパラメータを求める手法です。独立の試行においては同時確率は積の形で表すことができるので、手元のサンプルに独立性を仮定し、積の形で表します。
上記の説明だと抽象的なので、ロジスティック回帰において背景に仮定するベルヌーイ分布をまず考えます。ベルヌーイ分布は下記の数式で表すことができます。
$P(y|θ)=θ^y(1-θ)^{1-y}$ $y\in{0,1}$
上記において $θ$ はある事象（ex.コイン投げ）が観測される確率で、二値が前提のため観測されない確率は $1-θ$ となります。上記に沿って2値のそれぞれの確率が計算されます。ここでサンプルが独立的に生成されたと仮定すると、同時確率は積の形で表せるため、尤度 $L(θ|Y)$ は下記のように立式できます。（Yを観測されたサンプルの集合と考えます）
$L(θ|Y)=P(y_{1}|θ)×P(y_{2}|θ)×...×P(y_{n}|θ)=\Pi{P(y_{k}|θ)}$
上記の最大値を求めるのですが、ここで上記の微分を行った際に掛け算の形式になっているのが扱いづらそうです。そのため、単調増加関数である $log$ を用いて変換を行います。
$log(L(θ|Y))=log(\Pi{P(y_{k}|θ)})=\sum{log(P(y_{k}|θ))}=\sum{log(θ^{y_{k}}(1-θ)^{1-y_{k}})}$
この時の $log(L(θ|Y))$ を対数尤度と呼びます。この数式において $θ$ のみが $w$ を用いて定義されているため、 $w$ を変数とみなした際に、
$log(L(θ|Y))=\sum{log(θ^{y_{k}}(1-θ)^{1-y_{k}})}=\sum(y_{k}log(θ)+(1-y_{k})log(1-θ))$
のように書き直すことができ、これは交差エントロピーの式に一致します。
最小二乗法が目的変数の分布に正規分布を仮定して最尤法を用いて計算するのと同様に、交差エントロピーを用いた最適化はベルヌーイ分布を仮定した際に最尤法を用いて計算することに帰着します。
したがって、DeepLearningの計算の際によく用いる交差エントロピーの関数も実は最尤法がベースになっているということが上記のことから言えます。

2. #2の例題再訪＆ロジスティック回帰における勾配法
1節では背景知識としてベルヌーイ分布における最尤法と交差エントロピーについてまとめました。2節では具体的な問題の紹介と、交差エントロピーの勾配の計算を行います。
まず問題のご紹介ですが、#2の例題を引き続き用います。

import numpy as np
import matplotlib.pyplot as plt

x = np.arange(6,8,0.2)
n = np.array([70,74,68,69,74,86,43,56,78,48])
y = np.array([ 7, 14, 19, 28, 38, 50, 30, 45, 70, 48])

n_y = n-y
pi_ = y/n
print(pi_)

plt.scatter(x,pi_)
plt.show()

上記の実行結果は下記のようになります。

f:id:lib-arts:20190417183854p:plain
これに対し、3節ではロジスティック回帰の学習の実装を行っていきます。
問題については一通り把握できたので、次に交差エントロピーの勾配を求めます。まず、ロジスティック回帰は
$θ_{k}=\frac{1}{1+exp(-\beta^Tx_{k})}$
のように線形式部分との連結を行います。これを踏まえた上で交差エントロピー関数の微分を行うと（簡易化のためにlをおきます）、
$∇E(w)=\frac{δl}{δ \beta}=\sum_{k=1}^n(θ_{k}-y_{k})φ_{k}$
上記のようになります。
ここで計算した勾配を元に3節ではパラメータの学習を行なっていきます。

3. 実装で理解するロジスティック回帰（勾配法）
2節で定義した例題に対し、勾配法で学習を実装したのが下記になります。

def logistic(a):
return 1/(1+np.exp(-a))

beta = np.array([0.,0.])
for i in range(100000):
a = beta[0]+beta[1]*x
y_t = (logistic(a)-1)*y+logistic(a)*n_y
grad = np.array([np.sum(y_t),np.sum(y_t*x)])
beta[0] = beta[0] - 0.0001*grad[0]
beta[1] = beta[1] - 0.0001*grad[1]
if (i+1)%20000==0:
print("STEP_"+str(i+1))
print(beta)
print(grad)
print("======")
print(beta)

同じxの値でカウントデータの形式で計測した前提の結果を用いたため少々トリッキーな実装になってしまいましたが、意味的には同じと捉えられるかと思います。
実行結果は下記のようになります。

f:id:lib-arts:20190519154803p:plain
学習率などのパラメータ設定や収束にあたってのタイムステップの設定にやや時間がかかりましたが、シンプルな実装で実現することができました。
また下記を実行することで#2の時と類似のグラフを描くことができます。

import matplotlib.pyplot as plt
x_ = np.arange(5.5,8.5,0.05)
y_ = 1/(1+np.exp(-(beta[0]+beta[1]*x_)))
plt.scatter(x_,y_,color="red")
plt.scatter(x,pi)
plt.show()

実行結果は下記になります。

f:id:lib-arts:20190519154840p:plain

4. まとめ
#3では#2に引き続きロジスティック回帰を題材に学習のアルゴリズムを実装しました。
#2で用いたスコア法よりも、今回用いた勾配法はシンプルなのでより直感的に理解しやすい内容だったのではないかと思われます。
#4ではGLMの他の例も見れればということで、ポアソン回帰について取り扱います。