Title:WinBUGS

首先假設看到這個問題的人，有在使用WinBUGS而且已經妥善安裝好了。 使用WinBUGS最麻煩的問題是要在他的software下去調整參數的設定，其實是一件非常麻煩的事情，特別是要做一些大量的模擬研究的時候，所以我個人還是習慣在R底下做事，R2WinBUGS就可以解決這個問題。

以下以Bayesian change point的問題來做個介紹，這個問題以年為單位發生coal mining災難的次數，分析並找出平均災難發生次數變化的年份， 用R2WinBUGS不能避免的還是要另外寫一個.BUG檔來描述整個Bayesian model。

#bugs model
model
{
  for(year in 1:N)
  {
  D[year] ~ dpois(mu[year])
  log(mu[year]) ~ b[1] + step(year - changeyear) * b[2]
  }
  for(j in 1:2) {b[j] ~ dnorm(0.0, 1.0E-6)}
  changeyear ~ dunif(1,N)
}

將以上模型存成"coalmining.bug"，除此之外都在R裡面進行，以下為相關的參數設定與程式執行。 注意參數的名稱要與.bug檔一致

library(R2WinBUGS)

## Loading required package: coda
## Loading required package: lattice
## Loading required package: boot
## 
## Attaching package: 'boot'
## 
## The following object is masked from 'package:lattice':
## 
##     melanoma

library(coda)

 N=112
 D=c(4,5,4,1,0,4,3,4,0,6,
 3,3,4,0,2,6,3,3,5,4,5,3,1,4,4,1,5,5,3,4,2,5,2,2,3,4,2,1,3,2,
 1,1,1,1,1,3,0,0,1,0,1,1,0,0,3,1,0,3,2,2,
 0,1,1,1,0,1,0,1,0,0,0,2,1,0,0,0,1,1,0,2,
 2,3,1,1,2,1,1,1,1,2,4,2,0,0,0,1,4,0,0,0,
 1,0,0,0,0,0,1,0,0,1,0,0)
 data=list("N","D")
 parameters <- c("changeyear","b")
 inits = function() {list(b=c(0,0),changeyear=50)}

 coalmining.sim <- bugs (data, inits, parameters,
 "C:/Users/tai-chi/Documents/coalmining.bug", n.chains=1, n.iter=500,codaPkg=TRUE,debug=FALSE)

 coal.coda=read.bugs(coalmining.sim)

## Abstracting b[1] ... 250 valid values
## Abstracting b[2] ... 250 valid values
## Abstracting changeyear ... 250 valid values
## Abstracting deviance ... 250 valid values

 head(coal.coda)

## [[1]]
## Markov Chain Monte Carlo (MCMC) output:
## Start = 251 
## End = 257 
## Thinning interval = 1 
##       b[1]   b[2] changeyear deviance
## 251 0.9941 -1.229      41.82    338.5
## 252 1.1020 -1.384      39.78    336.5
## 253 1.1500 -1.099      40.50    336.3
## 254 1.1990 -1.329      40.03    334.5
## 255 1.1720 -1.516      41.28    338.7
## 256 1.1500 -1.224      35.78    341.3
## 257 1.2260 -1.352      39.33    335.3
## 
## attr(,"class")
## [1] "mcmc.list"