廿TT

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

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

今日の川柳

abrahamcow.hatenablog.com

コードだけ貼ります。

log_softmax関数の微分が難しかった。

SGDmixPois<-function(y, L, theta, lern_rate, num_iters, batch_size){
  ll<-function(y,theta){
    phi <- c(0,theta[1:(L-1)])
    lam <- theta[-c(1:(L-1))]
    lp <- numeric(L)
    for(l in 1:L){
      lp[l] <- phi[l]-log(sum(exp(phi)))+dpois(y,exp(lam[l]),log = TRUE)
    }
    log(sum(exp(lp)))
  }
  
  dll<-function(y,theta){
    phi <- c(0,theta[1:(L-1)])
    lam <- theta[-c(1:L-1)]
    lp <- numeric(L)
    for(l in 1:L){
      lp[l] <- phi[l]-log(sum(exp(phi)))+dpois(y,exp(lam[l]),log = TRUE)
    }
    softmax_lp <- exp(lp)/sum(exp(lp))
    softmax_phi <- exp(phi)/sum(exp(phi))
    mat <-matrix(NA_real_,L,L)
    for(i in 1:L){
      for(j in 1:L){
        mat[i,j] <- softmax_lp[i]*((i==j)-softmax_phi[j])
      }
    }
    dphi <- colSums(mat)
    dlam <- softmax_lp*(y-exp(lam))
    c(dphi[-1],dlam)
  }
  ll_hist <- numeric(num_iters)
  pb <- txtProgressBar(1,num_iters,style = 3)
  for(i in 1:num_iters){
    setTxtProgressBar(pb,i)
    theta <- theta + lern_rate*rowMeans(sapply(sample(y,batch_size),dll,theta=theta))
    ll_hist[i]  <- sum(sapply(y, ll, theta=theta))
  }
  list(theta=theta, ll_hist=ll_hist)
}
set.seed(123)
y <- sample(c(rpois(100,1),rpois(100,5),rpois(100,10)))
results <- SGDmixPois(y, L = 3, theta=runif(5), lern_rate = 0.1, num_iters =2000, batch_size=100)
plot(results$ll_hist,type = "l",xlab = "iteration",ylab = "log-lik")
print(exp(results$theta[-c(1:2)]))
#[1]  0.9248605  4.9462342 10.0615042
print(softmax(c(0,results$theta[1:2])))
#[1] 0.2969639 0.3973950 0.3056410

f:id:abrahamcow:20181004222054p:plain