ドラゴンボールの勝敗結果から戦闘力を推定する

以下のような個人データと勝敗データを扱います。

CID	name	power.level
1	ベジータ	18000
2	悟空	8000
3	ナッパ	4000
4	クリリン	1770
5	栽培マン	1200

winner	loser
1	2
1	3
1	5
2	3
3	4
4	5

モデルを記述する部分であるBUGSコードは以下になります。

model{
   for (game in 1:N.game) {
      winner.performance[game] ~ dnorm(skill[Winner[game]], 1)
      loser.performance[game]  ~ dnorm(skill[Loser[game]], 1)
      z[game] <- 1
      z[game] ~ dbern(constraint[game])
      constraint[game] <- step(winner.performance[game] - loser.performance[game])
   }

   for (i in 1:N.member) {
      skill[i] ~ dnorm(0, 1.0E-4)
   }
}

真の実力をskill、試合ごとに実力を発揮した力をperformanceとします。performanceの大小で勝敗が決まると考えるわけです。そしてperformanceは平均がskillの正規分布に従うと考えます。本当は実力を発揮する“ムラ”は個人ごとにだいぶ違うのは周知の事実です。例えば、ある人は安定していて格下にはキッチリ勝つけど格上には勝てないとか、ある人は格下にも負けるけど格上にも勝つときがあるとかありますよね。しかし、今回は勝敗データが少ないのでそこを推定するのは厳しいです。そこで、その“ムラ”の分散は「1」と固定しました。

BUGSコード的に少し難しい箇所は『performanceが高いから勝ち』という制限をうまく入れる必要がありまして、BUGS Bookの9.7節を参考に入れています（5行～7行）。このようにzを<-と~の二つで定義することにより、<-が成り立った時のみサンプリングするようにできます。詳しくはBUGS Bookを参照。

実行するRコードは以下になります。

source("R2WBwrapper.R")  # kubo-san's script: [http://hosho.ees.hokudai.ac.jp/~kubo/ce/RtoWbwrapper.html:title]

d <- read.delim("input/game.txt", as.is=TRUE, sep="\t")

N.member <- 5
N.game <- nrow(d)

clear.data.param()

set.data("N.member", N.member)
set.data("N.game", N.game)
set.data("Winner", d$winner)
set.data("Loser", d$loser)

set.param("winner.performance", rep(1, N.game))
set.param("loser.performance", rep(0, N.game))
set.param("skill", rep(0, N.member))

post.bugs <- call.bugs(
   # debug = FALSE,
   file = "model/model.bugs",
   n.iter = 10000000, n.burnin = 20000, n.thin = 2000
)
post.list <- to.list(post.bugs)
post.mcmc <- to.mcmc(post.bugs)
save.image("output/result.RData")