習作;ラベルスイッチング問題をなんとかする関数

前回の続き。Stanコードに対数尤度を出す部分を追加。。。
(どう見てもカッコ悪い。スマートに書き直したいものだ。ていうか,そもそもこれで合ってるのかな?)


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]);
  }

}

で,鎖が複数になるとラベルスイッチング問題が生じるんだけど,これを解決するlabel.switching関数を使ってみる。

## どこに分類されたか
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

練習ここまで。