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

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

実装のほとんどが車輪の再発明に近くなりますが、ベーシックなアルゴリズムをあえて一から自分で追ってみるというのは引き出しを増やすという意味で非常に有意義です。
#1では主成分分析について実装を通して理解を深める手助けとしました。

#2ではよく用いられる一般化線形モデル(GLM; Generalized Linear Model)の例としてロジスティック回帰を取り扱います。
以下目次になります。

1. 指数型分布族と一般化線形モデル(GLM)
2. ロジスティック回帰で学ぶGLMの具体例
3. 実装で理解するロジスティック回帰分析
4. まとめ

1. 指数型分布族と一般化線形モデル(GLM)
一般化線形モデル(GLM; Generalized Linear Model)は、通常の線形回帰のモデル(LM; Linear Model)の拡張と考えられます。具体的な拡張については、主に二点があります。
まず一点目は目的変数の分布がLMでは正規分布を仮定しているのに対し、GLMでは指数型分布族(Exponential Family of Distributions)を仮定します。ちなみに正規分布も指数型分布族に含まれるため、この観点においてGLMはLMの拡張となっています。
次に二点目は目的変数の期待値 $\mu$ （推論時は予測値となる）のパラメータと線形式部分との関係性です。数式で表すと $g(\mu_{i})=x^T_{i}\beta$ を考え、このg(x)のことをリンク関数(link function)と呼んでいます。LMではg(x)=xの恒等関数を仮定していたのに対し、GLMでは $g(\pi)=log(\frac{\pi}{1-\pi})$ （ロジスティック回帰）や $g(\mu)=log(\mu)$ （ポワソン回帰）など様々なリンク関数があります。
上記の二点がLMからGLMへの拡張にあたっての具体的な拡張になります。
指数型分布族とGLMについては下記でもまとめていますので、こちらも参考にしていただけたらと思います。

2. ロジスティック回帰で学ぶGLMの具体例
GLMの理解にあたって、概念だけわかっても通常の回帰分析とは具体的にどのように違うのかはなかなかイメージがつきません。そのため実際に具体例を取り扱っていこうということで、今回はロジスティック回帰を取り扱います。最初にロジスティック回帰を取り扱う理由としては、回帰分析の出力を確率的にするにあたってのシンプルな拡張として用いることができることに加え、近年猛威を振るうDeepLearningもロジスティック回帰の拡張と見なせなくもなく、一つ理解するだけで多方面で役に立つことが期待できるからです。
ロジスティック回帰の概要についてまとめたので、具体的な問題について解いてみましょう。

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])
pi = y/n
print(pi)

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

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

f:id:lib-arts:20190417183854p:plain
単なるx,yなどだとイメージがつきにくいので上記の具体的な問題設定として、部屋の家賃設定を考えた際に、6万〜7.8万まで2,000円単位で値上げした際に(x)、顧客数(n)、成約した顧客数(y)、成約率(pi)のデータが得られたとします。この際に、家賃金額から成約率を予測するロジスティック回帰のモデルを作成するとします。
解法としては、まずリンク関数として $g(\pi)=log(\frac{\pi}{1-\pi})$ 、確率分布としてはベルヌーイ分布を仮定します。これらの仮定の下での最尤法により、対数尤度関数は下記のように計算できます。
$l=\sum_{i=1}^{N}(y_{i}(β_{1}+β_{2}x_{i})-n_{i}log(1+exp(β_{1}+β_{2}x_{i}))+\binom{n_{i}}{x_{i}})$
以下ニュートン・ラプソン法をベースにした手法であるスコア法(method of scoring)を用いて解いて行きます。この際にパラメータ $β_{1}$ 、 $β_{2}$ に関するスコア関数は、
$U_{1}=\frac{δl}{δβ_{1}}=\sum_{i=1}^{N}(y_{i}-n_{i}π_{i})$
$U_{2}=\frac{δl}{δβ_{2}}=\sum_{i=1}^{N}x_{i}(y_{i}-n_{i}π_{i})$
となり、情報行列は
$I=\left(\begin{array}{cc} \sum_{i=1}^{N}n_{i}π_{i}(1-π_{i}) \sum_{i=1}^{N}n_{i}x_{i}π_{i}(1-π_{i}) \\ \sum_{i=1}^{N}n_{i}x_{i}π_{i}(1-π_{i}) \sum_{i=1}^{N}n_{i}x_{i}^2π_{i}(1-π_{i}) \end{array}\right)$

と計算できます。ここで最尤推定値は逐次式の
$b_{m}=b_{m-1}+I_{m-1}^{-1}U_{m-1}$
を解くと得られます。

3. 実装で理解するロジスティック回帰
2節ではGLMの具体例としてロジスティック回帰を挙げ、例題の紹介と解法の提示を行いました。3節ではこれに対応した実装を取り扱っていければと思います。

from scipy import linalg

beta = np.array([0.,0.])
for i in range(0,10):
pi_ = np.exp(beta[0]+beta[1]*x)/(1+np.exp(beta[0]+beta[1]*x))
U = np.array([np.sum(y-n*pi_),np.sum(x*(y-n*pi_))])
I = np.array([[np.sum(n*pi_*(1-pi_)),np.sum(n*x*pi_*(1-pi_))],[np.sum(n*x*pi_*(1-pi_)),np.sum(n*(x**2)*pi_*(1-pi_))]])
beta = beta+linalg.inv(I).dot(U)
print(beta)

上記を実行すると、下記のように逐次的なパラメータの推定量を得ることができます。

f:id:lib-arts:20190417184128p:plain
また、piとxの関係性を散布図とともにプロットすると下記のようになります。

import matplotlib.pyplot as plt

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

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

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

上記を確認することで、xの値が大きくなったり小さくなったりしても、piの値域が0~1の範囲におさまるというのが確認できます。

4. まとめ
ロジスティック回帰自体はscikit-learnなどのライブラリに実装されているものですが、あえて自前で実装してみるというのもなかなか面白いのではと思います。
導出の過程から最終的なグラフを実際に描いてみることで、スコア法によるアルゴリズムのイメージがつかめて非常に有意義な印象でした。