StatModeling Memorandum

StatModeling Memorandum

StanとRとPythonでベイズ統計モデリングします. たまに書評.

WinBUGSのTrapを止めるには

はじめに

  • モデルがよくない/まちがっていることが多いです。徹底的に可視化してからモデル式を考えましょう。
  • 観測データにモデリングを破綻させるような凄い「外れ値」が含まれているとダメ。
  • 全個体の説明変数が同じ値だとダメ。
  • Linux上でwineを使ってWinBUGSを動かしている場合だと動作が不安定な時があります。

WinBUGSコードの規則によるもの

  • ある変数が左辺(定義)に何度も出てきているとダメ(特殊な例として0-trickや1-trickと呼ばれているものがあります。こういう風に制限をかける場合を参照。)。
  • dnorm(a + b*x[i], tau) といった、stochastic node(~) の右辺の引数に計算式は入れられません。行を分けましょう。
  • inprod(x[i,], exp(z[i,])) といった inprod(,) の中に exp() があるような書き方もダメらしい。行を分けましょう。
  • Truncated sampling distributionとは、ある値以下(or 以上)のものは観測できずデータポイントとして記録にすら残らない確率分布です。これはWinBUGSの I(.,.) では扱えません。
  • Dirichlet分布(ddirch)、Wishart分布(dwish)は共役事前分布としてだけ使ってOK。パラメータは既知である必要があるとのことです。
  • 多項分布は事前分布として使ってはダメで自由度Nは既知である必要があるとのことです。
  • 多変量正規分布や多変量t分布はどこでも(i.e. 事前分布でも尤度でも)使えます。パラメータに制限はないとのことです。

Rのエラー

Error in is.finite(x) : default method not implemented for type 'list'

2次元データをmatrixではなくてdata.frameで渡してしまっている。→as.matrixしてから渡しましょう。

compile時に出るエラー

multiple definitions of node hogehoge

ある変数が何度も左辺に登場する場合や、deterministic nodeにデータを与えている場合に出ます。

multivariate distribution must have more than one component
for (j in 1:N.w) {
   w[j] ~ dmnorm(Zero[], tau[,])
}

はダメ。以下のように書きます。

w[1:N.w] ~ dmnorm(Zero[], tau[,])
vector valued relation w must involve consecutive elements of variable
for (k in 1:K) {
   w[1:N.dim, k] ~ dmnorm(Zero[], tau[,])
}
tau[1:N.dim, 1:N.dim] ~ dwish(R[,], Deg.freedom)

はダメ。以下のように書きます。

for (k in 1:K) {
   w[k, 1:N.dim] ~ dmnorm(Zero[], tau[,])
}
tau[1:N.dim, 1:N.dim] ~ dwish(R[,], Deg.freedom)
expected multivariate node

centerxiはstochastic nodeとすると、

mn[1, 1:N.dim] <- center[] + xi[]
mn[2, 1:N.dim] <- center[] - xi[]

はダメ。以下のように書きます。

for (d in 1:N.dim){
   mn[1, d] <- center[d] + xi[d]
   mn[2, d] <- center[d] - xi[d]
}

data loaded データを読み込む時に出るエラー

undefined variable

WinBUGSに(計算で使わない)余計なデータをRから渡しています。

variable Hogehoge is not defined in model or in data set

WinBUGSで使用する変数がRから渡されていないです。

inits(...) 初期値を与える時に出るメッセージ(エラーではありません)

this chain contains uninitialized variables

一部の初期値が設定してないから、事前分布から乱数を使って発生させますよの意味。モデルによってはたまに変な外れ値的な初期値が渡されてtrap行きになる模様。

inits(...) 初期値を与える時に出るエラー

SplusData1000

Rから渡した初期値の不具合。例として渡すデータの次元とmodel内の変数の次元と一致していない。Rによって作られたWinBUGS用データ(inits1.txt)をよく見てみましょう。

実行時のエラー

undefined real result

計算中にオーバーフローが発生したことを示唆します。対策は以下の通り。

  • 数値的に不可能な値の発生(例:負の値のlog・lowerとupperが逆になっている切断正規分布など)がないかcheckします。
  • もっとinformativeな事前分布を使うようにします(例1:一様分布の範囲を狭くする([0,100]→[0,10]) 例2: Wishart分布の自由度を上げる)。
  • 初期値が全て同一の値の時に発生することがあるらしいので一部を変えます。
  • 分散の大きな事前分布から発生させた初期値が極端に大きいと発生します。適切な初期値を渡すようにします。
  • 計算途中で極端に大きい(or 小さい)値になる場合にも発生します。プロビットモデルでは特にこの問題が起こりやすいため、分布を適当な範囲に制限したほうがよいとのことです。下記の例参照。
probit(p[i]) <- delta[i]
delta[i] ~ dnorm(mu[i], tau) I(-5, 5)

probit()の代わりにphi()を使用する方が問題が起こりにくい場合もあります。

p[i] <- phi(delta[i])
Trap 20 (precondition violated), Trap 21 (precondition violated)

不可能な初期値を与えた時に出てくることがあります。例: Beta分布なのに (0,1) を超える範囲の初期値とか。

Trap 66 (postcondition violated)

事前分布の分散が広すぎの可能性が高いです。残りの対策は上記の「undefined real result」のところと同様です。

incompatible copy

ループの中などで同じ変数を何度も定義してしまっているような場合に出るらしいです。

for (i in 1:N) {
   log.mu <- alpha + b[i]
}
stack overflow

再帰的なdeterministic node(<-)がある場合 →stochastic node(~) に置き換えられるか検討しましょう

NIL deference (read)

配列をスカラ値に変換するような間違ったことやっている時に出てくることがあるらしいです。

収束しない・とても遅い

  • 説明変数の中央化(平均0)と標準化(標準偏差1)、ただしfactorや離散値の場合は逆効果の場合もあります。
  • dlnorm(対数正規分布)は収束しづらいです。→logとってからdnormdpoisを使うことで代替できないか検討しましょう。
  • パラメータが非線形になっているモデルは避けましょう。p1/p2, p1*p2など。
  • パラメータのどれか(どちらか)が大きければOKというようなモデルは避けましょう。p1が大きいところとp2が大きいところでさまよって収束しないことがあります。→p1をゼロに固定する、p1の平均をゼロに固定するなど。  →例: 下記のモデルで tau.xtau を両方とも無情報事前分布にするとダメ。
for (i in 1:N){
   Y[i] ~ dbern(p[i])
   logit(p[i]) <- alpha + beta * x[i]
   x[i] ~ dnorm(0, tau.x)
   alpha ~ dnorm(0, tau)
   beta ~ dnorm(0, tau)
}
  • 収束が悪いパラメータが関わる部分はシンプルにしましょう(グループによって平均が変わる・分散が変わるとしないで、同じ平均・分散を使うとか)。
  • 変数変換を行うことで、各パラメータの独立性をアップさせることができないか考えましょう。
  • 混交のよくないMCMCにおいて、「まれに遠くの値への遷移が発生する」ようにしておくと混交が改善します。上記でprobitの範囲を制限しているのも似た発想です。logitの場合にも使えます。  →例: 下記のように標準偏差を冗長に定義することにより、値の大きなところへの遷移が可能になります(cf. The BUGS Book 10.5のコードを改変)。
for (i in 1:N){
   Y[i] ~ dbern(p[i])
   logit(p[i]) <- alpha + xi * r[i]
   r[i] ~ dnorm(0, tau.eta)
}
alpha ~ dnorm(0, 1E-4)
tau.eta <- 1/(s.eta*s.eta)
s.eta <- s/abs(xi)
xi ~ dnorm(0, 1)
s ~ dunif(0, 100)
  • logit()よりもphi()を使ったほうが収束がよくなることがあります。逆もしかり。
  • Treatが0,1のデータである時に、
model {
   for (i in 1:N) {
      Y[i] <- alpha + beta[Treat[i]+1]
   }
   beta[1] <- 0
   beta[2] <- beta.1
   beta.1 ~ dnorm(0, 1.0E-4)
}

というnest indexを使った書き方は収束がよくない時があります。以下のように書くと大丈夫です。

model {
   for (i in 1:N) {
      Y[i] <- alpha + beta * Treat[i]
   }
   beta ~ dnorm(0, 1.0E-4)
}
  • dcat()の中身のパラメータが 変化する/推定したい 場合、自動で合計1にしてくれるJAGSの方が安定の様子。WinBUGSはすぐTrapになるし、頑張って回避しても収束しないことがある。
  • 最終手段としてWinBUGSにはover relaxオプションというものがあります。{R2WinBUGS}パッケージのbugs関数の引数で over.relax=TRUEを指定すれば使うことができます。この方法では繰り返しのたびに複数のサンプルを生成し、このうち現在の値と負の相関を示す値を選択します。これにより1回の繰り返しに要する時間は少し増加しますが、chain内の相関が減少するため必要な繰り返し数は少なくなります。この方法はいつでも有効というわけではありません。やや自明ですがモデルのパラメータが多い場合にはうまくいかないことが多いです。繰り返し数(n.iter)を十分大きくするなどした後、万策尽きてから試してみてください。chainの混合が改善したかどうかを確認するには自己相関関数が利用できます。