生態学のデータ解析 - BUGS コーディング

• BUGS 言語を使ったベイズモデル のかきかた解説
  • 参照: WinBUGS 雑, WinBUGS 多変量正規分布, bugs() 関数の返り値, coda 雑, JAGS 雑
• とりあえず manual をよく読んでください: WinBUGS User Manual
  • 余計なこと (WinBUGS GUI の使いかたとか) も書かれていますが，有用な情報もいろいろ掲載されています

[もくじ]

• 概要
• 構造
• ルール
• こつとか?

概要 

• BUGS はベイズ統計モデルを記述するための言語である
• みかけはプログラミング言語っぽいけれど，ふつーのプログラミング言語ではないと考えたほうがよい
  • 「宣言型」プログラミング言語と分類されることもある
  • 「式」みたいなものが並んでいるだけのように見える
  • これらの「式」のならび順は結果にほとんど影響を与えない

構造 

• 変数 (node)
  • 数値 (浮動小数点数) だけを格納する
  • 文字とか数値ではないものは格納できない
  • 変数名として使える文字は A-Z,a-z,0-9,.(ピリオド)
• 変数のタイプ
  • 定数やデータを格納しておくもの (定数的な使いかた)
  • 決定論的関係 (<-) で定義される変数
  • 確率論的関係 (
) で定義される変数
• 配列 (array)
  • y[1],y[2],... と書く
    • 添字は #{1, 2, 3, ...}
  • 多次元配列も可能 x[1:3,1:4] など

• for ループ

(例)
for (i in 1:3) {
	Y[i] ~ dpois(lambda)
}

(上記は以下のように展開されると考える)
	Y[1] ~ dpois(lambda)
	Y[2] ~ dpois(lambda)
	Y[3] ~ dpois(lambda)

• if ~ then ~ else や for 以外の制御構造はない
• 組みこみ関数 step() (ステップ関数) は擬似 if として使える

ルール 

• 各行の # 以降はコメント
• 式の区切りは ; セミコロンだけど，省略可 (JAGS では省略できない)
• (上と矛盾するようだが) 複数の行にまたがって記述する場合はカッコでくくったほうが安全かも

(例)
lamba[i] <- (
	a
	+ b[1] * x[1, i]
	+ b[2] * x[2, i]
	+ b[3] * x[3, i]
	+ b[4] * x[4, i]
)

• 定数は決定論的関係 <- で定義する

(例)
Tau.noninformative <- 1.0E-4

• 変数を決定論的関係 <- で定義するときには， 同じ内容なのに複数のかきかたでかく場合がある

(例)
lambda[i] <- exp(a + b * x[i])
log(lambda[i]) <- a + b * x[i]

q <- 1 / (1 + exp(-a)) # ロジスティック関数
logit(q) <- a

• 変数を確率論的関係 ~ で定義するときには， 右辺を複雑にできない

(だめな例)
Y[i] ~ dpois(exp(a))

(正しい例: 2 行にわける)
Y[i] ~ dpois(lambda)
lambda <- exp(a)

• 変数を確率論的関係 <- で定義するときには， 右辺が複雑でもよい

(例)
lambda <- exp(a + b * log(x))

• 各変数は決定論的・確率論的関係で少なくとも一度ずつ定義できる

(しかしあまり見かけない例: 応答変数の変換?)
Z[i] <- log(Y[i])
Z[i] ~ dnorm(mu, tau)

• データが格納されておらず， しかも定義されない変数があるとエラーがでる

こつとか? 

• 線形予測子内の説明変数は中央化または標準化しておくと MCMC の収束が良くなる

• 階層ベイズモデルは定義のしかたで計算速度・収束の良さが異なる

(遅くなる書きかた)
for (i in 1:N) {
	Y[i] ~ dpois(lambda[i])
	lambda[i] <- exp(a + b * x[i] + r[i])
	r[i] ~ dnorm(0, tau)
}
tau ~ dgamma(P.gamma, P.gamma)

(速くなる書きかた)
for (i in 1:N) {
	Y[i] ~ dpois(lambda[i])
	lambda[i] <- exp(log.lambda[i])
	log.lambda[i] ~ dnorm(log.mean, tau)
	log.mean <- a + b * x[i]
}
tau ~ dgamma(P.gamma, P.gamma)