廿TT

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

崩壊型ギブスサンプリングによるポアソン混合分布の学習をRで

機械学習スタートアップシリーズ ベイズ推論による機械学習入門 (KS情報科学専門書)

機械学習スタートアップシリーズ ベイズ推論による機械学習入門 (KS情報科学専門書)

ベイズ推論による機械学習入門』の例題です。

R の dnbinom のパラメータの置き方がテキストと微妙に違うのに気づかずつまづいてた。

set.seed(1234)
N <- 100
ind <- sample.int(2,N,replace = TRUE)
lambda <- c(2,10)
y <- rpois(N,lambda[ind])
a <- 1 #事前分布のパラメータ
b <- 1 #事前分布のパラメータ
alpha <-c(2,2) #事前分布のパラメータ
softmax <- function(x){
  maxx <- max(x)
  exp(x-maxx)/sum(exp(x-maxx))
}
maxit <- 1000
shist <- array(NA_integer_,dim = c(N,2,maxit))
lphist <- numeric(maxit-1)
s <- t(rmultinom(N,1,c(0.5,0.5))) #潜在変数の初期化
alphahat <- colSums(s)+alpha
ahat <- colSums(s*y)+a
bhat <- colSums(s)+b
shist[,,1]<-s
for(i in 2:maxit){
  for(n in 1:N){
    ahat <- ahat - y[n]*shist[n,,i-1]
    alphahat <- alphahat - shist[n,,i-1]
    bhat <- bhat - shist[n,,i-1]
    unnorm <-dnbinom(y[n],ahat,1-1/(1+bhat),log = TRUE)+log(alphahat)
    lphist[i-1] <- lphist[i] + log(sum(exp(unnorm)))
    snew <- drop(rmultinom(1,1,softmax(unnorm)))
    shist[n,,i] <- snew
    alphahat <- alphahat + snew 
    ahat <- ahat + snew*y[n]
    bhat <- bhat + snew
  }
}
plot(lphist,type="l")
shat <-apply(shist[,,-c(1:100)], 1:2, mean)
cluster <-factor(apply(shat, 1, which.max))
mean(ind==cluster)
#94%くらい正解
library(ggplot2)
df <- data.frame(y,cluster)
ggplot(df,aes(x=y,fill=cluster,colour=cluster))+
  geom_histogram(position = "identity",alpha=0.5,bins=15)

条件付き予測尤度(っていうのか)のトレースラインです。

f:id:abrahamcow:20181110022813p:plain

クラスタリングの結果です。

f:id:abrahamcow:20181110022729p:plain