廿TT

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

rstan で混合二項分布のパラメータ推定

ordered 型は、「小さい順」という制約です。

StanとRでベイズ統計モデリングで解説されている「ラベルスイッチング」を回避するためにこれを使ってます。

StanとRでベイズ統計モデリング (Wonderful R)

StanとRでベイズ統計モデリング (Wonderful R)

inv_logit() で-∞~+∞ の値を 0~1 に変換しています。

なぜか binomial_logit_lpmf を使うとうまくいきませんでした。

// binomial2.stan
data{
  int<lower=0> n;
  int<lower=0> m;
  int<lower=0> y[m];
}
parameters{
  ordered[2] b;
  real<lower=0,upper=1> phi; 
}
transformed parameters{
  matrix[m,2] lp;
  for(i in 1:m){
    lp[i,1] =log(phi)+binomial_lpmf(y[i]|n,inv_logit(b[1]));
    lp[i,2] =log1m(phi)+binomial_lpmf(y[i]|n,inv_logit(sum(b)));
  }
}
model{
  for(i in 1:m){
    target += log_sum_exp(lp[i,]);
  }
}
generated quantities{
  matrix[m,2] label;
  for(i in 1:m){
    label[i,] = softmax(lp[i]')';
  }
}

サンプルが2つの二項分布の山のどちらから得られたかについての正答率は 0.83 でした。

2つの二項分布の成功確率パラメータが近すぎると収束しません。

library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = 4)
prob1 <- 0.4
prob2 <- 0.5
set.seed(1)
clabel <- sample.int(2,100,replace = TRUE,prob = c(0.4,0.6))
y <-rbinom(100,100,ifelse(clabel==1,prob1,prob2))
dat4stan <-list(n=100,m=100,y=y)
fit1 <-stan("binomial2.stan",data=dat4stan)
traceplot(fit1,pars=c("phi","b","lp__"))
b <-get_posterior_mean(fit1,pars=c("b"))[,"mean-all chains"]
phi <-get_posterior_mean(fit1,pars=c("phi"))[,"mean-all chains"]
s <-matrix(get_posterior_mean(fit1,pars="label")[,"mean-all chains"],100,2,byrow = TRUE)
clabel_hat <-apply(s, 1, which.max)
cat(plogis(b))
cat("正答率",mean(clabel==clabel_hat))