前回の続き。Stanコードに対数尤度を出す部分を追加。。。
(どう見てもカッコ悪い。スマートに書き直したいものだ。ていうか，そもそもこれで合ってるのかな？)
[code]
data{
int<lower =2> K; #number of clusters
int<lower =1> N; #number of observations
real X[N]; #observed data
}
parameters{
vector[K] mu;
vector<lower =0,upper=10>[K] sig2;
simplex[K] theta;
}
transformed parameters{
vector[K] ps[N];
for(n in 1:N){
for(k in 1:K){
ps[n,k] = log(theta[k])+normal_lpdf(X[n]|mu[k],sig2[k]);
}
}
}
model{
sig2 ~ cauchy(0,2.5);
mu ~ normal(0,100);
for(n in 1:N){
target += log_sum_exp(ps[n]);
}
}
generated quantities{
real log_lik[N];
simplex[K] u[N]; #class membership probability
for(n in 1:N){
u[n] = softmax(ps[n]);
log_lik[n] = log_sum_exp(ps[n]);
}
}
[/code]
で，鎖が複数になるとラベルスイッチング問題が生じるんだけど，これを解決するlabel.switching関数を使ってみる。
[code language=”R”]
## どこに分類されたか
allocK < – rstan::extract(fit_sp,pars="u")$u
# 各チェインごとに割り当てられたクラス
zChain <- matrix(ncol=N,nrow=itr*C)
for(i in 1:(itr*C)){
zChain[i,] <- apply(allocK[i,,],1,which.max)
}
## MCMCで得られたパラメタ
mcmc.params <- rstan::extract(fit_sp,pars=c("mu","sig2","theta"))
# mcmc * num.of.clusters * num.of.params
mcmc.pars <- array(data=NA,dim=c(itr*C,K,3))
mcmc.pars[,,1]<- mcmc.params$mu
mcmc.pars[,,2]<- mcmc.params$sig2
mcmc.pars[,,3]<- mcmc.params$theta
# この関数はSJW法を使う時に必要になってくる。completeオプションに関数を渡さないといけないので！
complete.normal.loglikelihood<-function(x,z,pars){
# x: data (size = n)
# z: allocation vector (size = n)
# pars: K\times J vector of normal mixture parameters:
# pars[k,1] = mean of the k-normal component
# pars[k,2] = variance of the k-normal component
# pars[k,3] = weight of the k-normal component
# k = 1,…,K
g <- dim(pars)[1] #K (number of mixture components)
n <- length(x) #this denotes the sample size
logl<- rep(0, n)
logpi <- log(pars[,3])
mean <- pars[,1]
sigma <- sqrt(pars[,2])
logl<-logpi[z] + dnorm(x,mean = mean[z],sd = sigma[z],log = TRUE)
return(sum(logl))
}
# パッケージのお出まし
library(label.switching)
# 全方法試す
set <- c("STEPHENS", "PRA", "ECR", "ECR-ITERATIVE-1", "ECR-ITERATIVE-2", "AIC", "SJW","DATA-BASED")
# 対数尤度
log_liks <- rstan::extract(fit_sp,pars="log_lik")$log_lik
library(loo)
loo(log_liks)
waic(log_liks)
#ピボットは対数尤度が一番大きいやつにしてみた（なんでもいいと思う）
# MCMCサンプルごとに対数尤度が一番大きいやつを算出(apply部)
# 表の形で集計，大きい順に並べかえる(sort(table)部分)
# 表のラベルを取り出して(names部分)，数字にする(as.numeric)
pivot <- as.numeric(names(sort(table(apply(log_liks,1,which.max)),decreasing=T)[1]))
# さあ実行
ls <- label.switching(method = set,
zpivot = zChain[pivot,],
prapivot = mcmc.pars[pivot, , ],
z = zChain,
K = 3,
p = allocK,
mcmc=mcmc.pars,
complete=complete.normal.loglikelihood,
data = X)
# 各推定法で結果が一致したかどうか
ls$similarity
# 推定された所属クラスタ番号
ls$clusters
[/code]
練習ここまで。