MCMC法（メトロポリス・ヘイスティングス法）による一般化モデルの最適化と実装｜スクラッチ実装で理解する機械学習アルゴリズム #5

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

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

#5では最適化についてもう少し見てみようということで、MCMC(Markov Chain Monte Carlo)法による最適化について解説と実装の紹介を行います。
以下目次になります。

1. MCMC法とメトロポリス・ヘイスティングス法
2. ロジスティック回帰復習
3. 実装で理解するMCMCによるロジスティック回帰の実装
4. まとめ

1. MCMC法とメトロポリス・ヘイスティングス法 f:id:lib-arts:20190529161757p:plain

メトロポリス・ヘイスティングス法 - Wikipedia
Wikipediaのメトロポリス・ヘイスティングス法(Metropolis-Hastings algorithm)のページでは上記のように解説がされています。MCMC法はざっくりいうと最適化にあたって最適化する対象のパラメータをサンプリングすることで、数値的に求めていく手法です。勾配法やニュートンラプソン法のように厳密な解が求まらない問題に対し、効力を発揮します。また、MCMC法が他の最適化手法と異なるポイントとしては、パラメータの分布が得られる点です。このことにより、パラメータの分布化を行うベイズの文脈では相性が良いです。
次にメトロポリス・ヘイスティングス法(Metropolis-Hastings algorithm)についてですが、考え方としては焼きなまし法(Simulated Annealing)に近いです。

f:id:lib-arts:20190529162326p:plain
導出については上記でまとまっています。ランダムウォークのようにランダムに移動していく形なのですが、パフォーマンスを向上させる向きに移動する際は必ず移動し、そうでないときは確率的に移動するかどうかを決めます。
出力としては勾配法のように最終結果を出力する形もできますが、MCMCの特徴としてはサンプルが得られるためパラメータの分布として得られることです。サンプルに関しては初期値に近いところは初期値の依存性が強いため棄却(burn-in)されます。
説明だけだとわかりにくいので細かいところについては実際に3節で実装を元に説明していきます。

2. ロジスティック回帰復習
2節ではロジスティック回帰の問題設定について簡単に復習します。

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm

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
p_ = y/n
print(p_)
plt.scatter(x,p_)
plt.show()

上記は#3で動かしたコードと同様のコードです（変数名だけ少々変更を行いました）。こちらを実行した結果のp_をxから予測するというのが今回の問題定義になります。ここで#3では勾配を求めていましたが、今回は下記の対数尤度を用います。
$log(L(θ|Y))=\sum{log(θ^y_{k}(1-θ)^{1-y_{k}})}=\sum(y_{k}log(θ)+(1-y_{k})log(1-θ))$
この時、 $θ_{k}=\frac{1}{1+exp(-w^Tx_{k})}$ のようにリンク関数を設定を行います。
立式ができたので、3節ではこちらを元に実装を行います。

3. 実装で理解するMCMCによるロジスティック回帰の実装
下記の実装を動かすことによってメトロポリス・ヘイスティング法を実行できます。（簡易化のため値が向上しない際の採択率を一律で0.05とおきました。）

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

def log_likelihood(theta):
res = np.sum(y*np.log(theta)+n_y*np.log(1-theta))
return res

np.random.seed(1)
beta = np.array([0.,0.])
beta_smpl = np.zeros([500000,2])
for i in range(520000):
a = beta[0]+beta[1]*x
theta = logistic(a)
e_beta = np.array([norm.rvs(0,0.01),norm.rvs(0,0.01)])
theta_after = logistic*1
print(beta)
print("======")
if (i+1)>20000:
beta_smpl[i-20000,:]=beta
print(beta)

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

STEP_100000
[-17.90291551 2.6360675 ]
======
STEP_200000
[-17.71179173 2.60581579]
======
STEP_300000
[-17.88177012 2.63418973]
======
STEP_400000
[-17.79246473 2.61918144]
======
STEP_500000
[-17.86695219 2.63036266]
======
[-17.8262738 2.62400147]

また、ここで下記のコードを実行することでパラメータの分布を得ることができます。

print(np.mean(beta_smpl[:,1]))
plt.hist(beta_smpl[:,1],bins=50)
plt.show()

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

f:id:lib-arts:20190529162234p:plain
曲線のプロットも#3と同様に行うと下記のようになります。

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

4. まとめ
#5ではMCMC法について取り扱いました。
一部大元のアルゴリズムから変更したりしてしまっていますが、なかなか面白い結果が得られたと思います。
#6以降ではガウス過程について取り扱ってみれればと思います。

*1:beta[0]+e_beta[0])+(beta[1]+e_beta[1])*x)
if log_likelihood(theta)<log_likelihood(theta_after):
beta[0] = beta[0]+e_beta[0]
beta[1] = beta[1]+e_beta[1]
else:
if np.random.rand()<0.05:
beta[0] = beta[0]+e_beta[0]
beta[1] = beta[1]+e_beta[1]
if (i+1)%100000==0: #log
print("STEP_"+str(i+1