Stan つかってますかー．収束してますかー？収束してませんねー？って事でチェーンが収束しないことありますよね．

そんな方に収束するモデルの書き方を伝授します．それはずばり“収束しないモデルを書かない事”です．つまりですよ，収束しないモデルを自由自在にかけるようになれば，自由自在に収束するモデルが書ける様になるという訳です．前置きはいいとしてやってみましょー

library(rstan) # 相棒を召還
## Loading required package: Rcpp
## Loading required package: inline
## 
## Attaching package: 'inline'
## 
## The following object is masked from 'package:Rcpp':
## 
##     registerPlugin
## 
## rstan (Version 2.4.0, packaged: 2014-07-21 17:16:17 UTC, GitRev: c995bc8b9d67)

初期値の取り方系

理屈はさておきコード

stancode <- "
parameters{
  real<lower=0,upper=200> a;
}
model{
  a ~ uniform(0,1);
}
"
fit <- stan(model_code = stancode)
## 
## TRANSLATING MODEL 'stancode' FROM Stan CODE TO C++ CODE NOW.
## COMPILING THE C++ CODE FOR MODEL 'stancode' NOW.
## 
## SAMPLING FOR MODEL 'stancode' NOW (CHAIN 1).
## Initialization failed after 100 attempts.  Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.
## error occurred during calling the sampler; sampling not done

はい，収束しませんねー．てかまずシミュレーションしてませんね．これで引っかかる人多いはずです．こうやればいける

stancode <- "
parameters{
  real<lower=0,upper=1> a;
}
model{
  a ~ uniform(0,1);
}
"
fit <- stan(model_code = stancode)
traceplot(fit)

plot of chunk unnamed-chunk-3

詳細はStan Reference P155です．ヒントは初期値．

つぎー

無情報事前分布？そりゃー臆病だよ．あんたは臆病だ！

stancode <- "
parameters{
  real a;
  real b;
}
model{
  a ~ normal(0,10000);
  b ~ normal(0,10000);

  0.0 ~ normal(a + b,1.0);
}
"
fit <- stan(model_code = stancode)
traceplot(fit)

plot of chunk unnamed-chunk-4

a + b = 0となるaとbの組み合わせ，つまり最尤値が無限個ある．そして事前情報が緩すぎる．

符号をいれるとうまくいく

stancode <- "
parameters{
  real<lower=0> a;
  real<lower=0> b;
}
model{
  a ~ normal(0,10000);
  b ~ normal(0,10000);

  0.0 ~ normal(a + b,1.0);
}
"
fit <- stan(model_code = stancode)
traceplot(fit)

plot of chunk unnamed-chunk-5

そもそも何もしらないもんモデリングしてんじゃねーよってことで 知ってる事があったら教えてあげてください．

stancode <- "
parameters{
  real a;
  real b;
}
model{
  a ~ uniform(0,100); # 例えば年齢ならこんなんで十分無情報
  b ~ normal(15,10); # 気温ならこんなんでも無情報

  0.0 ~ normal(a + b,1.0);
}
"
fit <- stan(model_code = stancode)
traceplot(fit)

plot of chunk unnamed-chunk-6

もちろんデータがまだあるんなら入れてあげてもうまくいく

stancode <- "
parameters{
  real a;
  real b;
}
model{
  a ~ normal(0,10000);
  b ~ normal(0,10000); 

  0.0 ~ normal(a + b,1.0);
  2.0 ~ normal(2.0*a + b,1.0);
}
"
fit <- stan(model_code = stancode)
traceplot(fit)

plot of chunk unnamed-chunk-7

あとなんか色々原因思いつくけどうまく再現できないです． またなんか事例増えたら記事か来ます．