廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

R: 確率的勾配降下法でポアソン分布のパラメータの最尤推定

今日の川柳

これはたぶん確率的勾配降下法のもっとも簡単な例題です。

モチベーション:でっかいデータで最尤推定したいときがあって、ふつうに準ニュートン法とか使うと遅すぎていやなので、なんか機械学習の人たちがやってる確率的勾配降下法とかいうやつ使えばいいんじゃないかと思った。

確率的勾配降下法 - Wikipedia

ポアソン分布のパラメータの最尤推定をします。

ポアソン分布のパラメータは正の値しかとらないけど、制約をつけて最適化するのが面倒なので指数関数の中に入れてしまいます。

あるサンプル n についての対数尤度は、
\displaystyle L_n = y_n\theta-\exp(\theta)+\mathrm{Const.}

対数尤度の一階微分は、
\displaystyle \frac{d}{d\theta}L_n = y_n-\exp(\theta)

確率的勾配降下法ではデータをミニバッチと呼ばれる m 個の部分集合に分けて、一階微分を評価しパラメータを更新します。

更新式は、
 \displaystyle \theta_{k+1} = \theta_k+\alpha \frac{1}{m}\sum_{n=1}^{m}\left(\frac{d}{d\theta}L_n\right)
ここでαは学習率と呼ばれる適当な定数。

R による実装。

set.seed(1234)
y <- rpois(100,5)
#hist(y)
ll<-function(theta,y){
  dpois(y,exp(theta),log = TRUE)
}
dll<-function(theta,y){
  y-exp(theta)
}
#dll(1,10) #微分があってるかどうか
#numDeriv::grad(ll,1,y=10) #数値微分と比べて確かめる
SGD<-function(ll, dll, y, theta, lern_rate, num_iters, batch_size=10L){
  ll_hist <- numeric(num_iters)
  for(i in 1:num_iters){
    theta <- theta + lern_rate*mean(dll(theta,sample(y,batch_size)))
    ll_hist[i]  <- sum(ll(theta,y))
  }
  list(theta=theta, ll_hist=ll_hist)
}

results <- SGD(ll, dll, y, theta=0, lern_rate = 0.1, num_iters = 100)
print(exp(results$theta))
#[1] 4.658927
plot(results$ll_hist,type = "l",xlab = "iteration",ylab = "log-lik")
print(mean(y)) #普通の最尤推定量
#[1] 4.57
abline(h=sum(dpois(y,mean(y),log = TRUE)),lty=2)

f:id:abrahamcow:20181004052740p:plain

ちゃんと尤度大きくなってるっぽいですね。

実装かんたんだなあ。

明日以降、機械学習関連のライブラリを無駄に使ってなんでもかんでもMAP推定してやろうかと思ってます。

深層学習 (機械学習プロフェッショナルシリーズ)

深層学習 (機械学習プロフェッショナルシリーズ)

買ったけどぜんぜん読んでなかった『深層学習』をはじめて参照した。

abrahamcow.hatenablog.com