打ち切り関数

打ち切りのあるデータ

野外のデータ、あるいは実験データを扱う際、「ある特定の値以上であることはわかるけれども、本当の値はわからない」という状況に出くわすことがあります。このようなデータタイプは「打ち切り」と呼ばれ、統計モデルを組むうえで厄介な問題です。該当する例はとして、以下のようなものが挙げられます。

  • ある一定量のエサを与え、一定時間後に無くなった量を調べる(全部消費された場合、もともと与えたエサ量「以上」を消費したかったことは確か)
  • 標識再捕獲を行い、ある時点で生きていたかどうかを調べる(最後の再捕獲で生きていた場合、その時点より”長く”生きていたことは確か)

他にも例は挙げればたくさんありますが、こういう性質のデータを扱うのは頭を抱えます。

センサリング(Censoring)とは

センサリングとは、そのまま「打ち切り」を意味する言葉なのですが、このテクニックにより上記の問題に対応することが可能です。センサリングを適用する場合、上限値に達した観測値は一時的に欠測、NAとして扱われます。BUGSソフトウェアでは、(通常の線形モデルの場合)応答変数の欠測は自動的に推定されるべきパラメーター(つまり、ほかの偏回帰係数などの推定には用いられない情報)として扱われ、ランダムな値があてがわれます。しかし、打ち切りの場合は部分的な情報(y >= 60など)を持っていますので、上限値に達した値に対しては、こちらで定義する上限値以上の値が(ランダムに?)割り振られます。そうした値が割り振られたうえで、線形モデルの各種パラメーターを推定することになります。

JAGSによる実装

基本的な実装の手順は、普通のBUGSモデルのコードに1行追加するだけで済みます。こちらのモデルファイルを参考にBUGSコードを書いてみます。

model{
#Likelihood
 for(i in 1:N.sample){ # loop for i
   isCensored[i] ~ dinterval(y[i],C[i])
   y[i] ~ dpois(lambda[i]) # poisson model
   log(lambda[i]) <- alpha + beta * X[i] # linear model structure
 }
 
#Prior
 alpha ~ dnorm(0, 0.001)
 beta ~ dnorm(0, 0.001)
}

dintervalが打ち切り関数になっており、isCensored[i](こちらでRから与えるデータ)はデータiが打ち切りかどうかを示す変数(0:打ち切り無し、1:打ち切りあり)、C[i]は打ち切った場合の上限値(こちらデータを与える;任意の変数)になります。なお、打ち切りありの場合、yは”NA”として与えることになります。例えば、観測の上限がy=7だった場合、データは下記のようになります。

> head(D)
   y isCensored C         x
1 NA          1 7 12.056546
2 NA          1 7 11.205303
3  4          0 7  9.865988
4 NA          1 7 12.633539
5 NA          1 7  9.116709
6  6          0 7 12.067588

図にするとこんなふうになります。

打ち切りのあるデータプロット。白丸は真の観察値、グレイは打ち切りされた観察値。

センサリングモデルと普通のポアソンモデルで推定値を比較してみます。

#true values; alpha=0.3, beta=0.15

#with censoring model
      Lower95 Median Upper95 Mean    SD Mode  MCerr MC%ofSD SSeff AC.200 psrf
alpha   0.135   0.59    1.08 0.59 0.241   NA 0.0167     6.9   208  0.100    1
beta    0.083   0.13    0.18 0.13 0.024   NA 0.0017     7.0   202  0.097    1

#with no censoring model
      Lower95 Median Upper95  Mean    SD Mode MCerr MC%ofSD SSeff  AC.200 psrf
alpha   0.587  0.970    1.41 0.972 0.212   NA 0.010       5   407 -0.0027    1
beta    0.026  0.071    0.11 0.071 0.021   NA 0.001       5   403 -0.0037    1

センサリング関数を使わない場合、y>=7の値はy=7とされてしまうので、推定値がかなり過小評価されています。それに対し、センサリング関数を用いたモデルでは、推定値の偏りが解消され、その95%CIも真の値を含んでいます。

注意点:センサリング関数を用いる場合、打ち切りされたデータについては、「打ち切り上限値以上の値」を初期値として与える必要があります。上の例で挙げると、y=NAとなっているデータについては、y=8などの値を初期値として与えます。そうしないとエラーメッセージが帰ってきます。

SPONSOR LINK

Varying slope model: 回帰係数に線形予測子を加える

Varying intercept modelでは切片のみがグループ間で変化することをモデル化していますが、ある変数xに対するyの応答が、グループによって異なることもあります。例えば、魚の種数yは、局所スケールの生息地面積xに従って増えますが、その応答の強さは流域スケールの変数である気温xgで変わるかもしれません。

例:生息地面積xに対する魚の種数yの応答。傾きは流域(色の違い)によって異なっている。

このような「空間スケールをまたがった交互作用」を扱いたい場合は、上の例ので考えると、生息地面積xの効果b(回帰係数)は、流域の平均気温xgによって変わる、という階層構造をもったモデルを構築することになります。

model{
  ninfo <- 0.0001
  
  # Prior
  for(j in 1:2){
    for(k in 1:2){
      B[j,k] ~ dnorm(0,ninfo)
    }
  }
  tau ~ dgamma(ninfo, ninfo); sigma <- sqrt(1/tau)
  TAU[1:2,1:2] ~ dwish(W[,], 3)
  SIGMA[1:2,1:2] <- inverse(TAU[,])
  
  # Likelihood
  for(i in 1:N){
    y[i] ~ dnorm(mu[i], tau)
    mu[i] <- b[group[i],1] + b[group[i],2]*x[i]
  }
  
  for(j in 1:Ng){
    b[j,1:2] ~ dmnorm(b.mu[j,], TAU[,])
    for(k in 1:2){
      b.mu[j,k] <- B[1,k] + B[2,k]*xg[j]
    }
  } 
}

上記のモデルでは、切片と偏回帰係数bがグループ(流域)で変化すると仮定しています。ただし、切片と偏回帰係数において、同じグループ間のばらつきを考えることになるので、それらを独立と考えるのは無理があります。そのため、b[,1]とb[,2]の誤差(グループレベルの説明変数で説明できないばらつき)は、共分散を考慮する必要があります。このスクリプトでは、bの誤差構造には多変量正規分布を用いており(ハイライト箇所21行目)、その分散共分散行列の逆数(TAU[,])の事前分布にウィシャート分布を指定しています(Wはdiag(2)としてRから与えた)。

SPONSOR LINK

Varying intercept model: 切片に線形予測子を加える

ランダム効果入りのベイズ推定モデルでは、グループ間の「ランダムなばらつき」を表現するモデルを作りましたが、そのグループ間のばらつきが何によって変わるのかを知りたい場合もあります。例えば、川にすむ魚の個体数をモデリングするとき、小さい空間スケールの個体数は隠れ場所の量によってよく説明できたとしても、河川ごとの平均的な個体数は河川単位で変化する他の要因、例えば降水量で決まるかもしれません(降水量が多い河川は氾濫が多いため、個体数は全体的に少ないなど)。

色の違いが川の違い;X軸は生息場所の量、Y軸は魚の密度。

このような問題には、ランダム効果に線形予測子を加えることで対処できます。下記スクリプトのハイライト箇所で、グループレベルの線形予測子がモデル化されています。なお、このグループレベルの説明変数だけで説明できないばらつきは、パラメーターsigma.Rで考慮されることになり、これによりランダムなグループ間のばらつきも表現されます。

model{
  ninfo <- 0.0001
  
  # Prior
  for(i in 1:2){ b[i] ~ dnorm(0, ninfo) }
  for(i in 1:2){ B[i] ~ dnorm(0, ninfo) }
  tau ~ dgamma(ninfo, ninfo); sigma <- sqrt(1/tau)
  tau.R ~ dgamma(ninfo, ninfo); sigma.R <- sqrt(1/tau.R)
  
  # Likelihood
  for(i in 1:N){
    y[i] ~ dnorm(mu[i], tau)
    mu[i] <- alpha[group[i]] + b[1]*x[i]
  }
  
  for(j in 1:Ng){
    alpha[j] ~ dnorm(mu.R[j], tau.R)
    mu.R[j] <- B[1] + B[2]*xg[j]
  } 
}

なお、注意点として、グループレベルの説明変数は、繰り返しの単位がグループIDと対応することになります。例えば、ここに10の河川から5地点ずつ集められたデータがあるとします(N=5*10=50)。グループレベルでは「河川」が繰り返しとのなるので、添え字jは1から10(Ng: Number of groups)の繰り返し構文となります。したがって、RからJAGSに渡すデータセットは二つのタイプに分かれます(y、x、groupは長さNのベクター、xgは長さNgのベクター)。

ランダム効果入りのベイズ推定モデル

ここでは一般的なGLMMをJAGSで実装するためのスクリプトを紹介します。GLMM、すなわちグループ(ランダム効果)によって切片の値が異なるとするモデルですが、その中身は以下のようなものになります。

Yi ~ Normal(μi, σ2)
μi = α + β*Xi + Rj(i)

上の式ではRj(i)がランダム効果ですが、これは平均0、分散σR2の正規分布に従う確率変数になります。j(グループID)ごとに正規分布に従う乱数が生成され、その値に応じて切片(グループごとのYの平均値)の値が変化します。平均は0に固定されており、問題となるσR2はデータから推定することになります(グループ間のばらつき)。このようなモデルをJAGSで実装する場合、ランダム効果も数字表記(Numeric)でなければならないので、その点に注意が必要です(lmerなどRの関数では因子型として指定;BUGSでは因子型変数は利用不可)。

データセットを用意します。

> head(D,10)
             x group         y
1   1.89336207     1 3.1029784
2   1.96143183     1 3.0784519
3   1.46765994     1 2.8607315
4  -1.61391038     1 2.1712377
5  -0.57316145     1 2.3753531
6   0.70379683     2 0.6382498
7   2.72372835     2 1.0363066
8   0.97064235     2 0.7026342
9   0.02245384     2 0.6256619
10  1.63627604     2 0.8633311

モデルスクリプトを書きます(model_random.Rとして保存)。ランダム効果は入れ子として表記(R[group[i]])されます。データセットはiをユニットとした繰り返しになるのに対し、ランダム効果はjをユニットとした繰り返しになります。

model{
  ninfo <- 0.0001
  
  # Prior
  for(i in 1:2){ b[i] ~ dnorm(0, ninfo) }
  tau ~ dgamma(ninfo, ninfo); sigma <- sqrt(1/tau)
  tau.R ~ dgamma(ninfo, ninfo); sigma.R <- sqrt(1/tau.R) #ランダム効果の精度の(超)事前分布
  
  # Likelihood
  for(i in 1:N){
    y[i] ~ dnorm(mu[i], tau)
    mu[i] <- b[1] + b[2]*x[i] + R[group[i]]
  }
  
  for(j in 1:Ng){ R[j] ~ dnorm(0, tau.R) } 
}

上記のモデルは正規分布になってますが、該当箇所をdpoisやdnegbinとすることで、ポアソンや負の二項分布のモデルでもランダム効果を加えることができます。とてもフレキシブルです。下にいくつかほかの確率分布の場合の例を挙げてみます。モデルを走らせるRコードも下につけますが、詳細はこちらを参照してください。

model{
  ninfo <- 0.0001
  
  # Prior
  for(i in 1:2){ b[i] ~ dnorm(0, ninfo) }
  tau ~ dgamma(ninfo, ninfo); sigma <- sqrt(1/tau)
  tau.R ~ dgamma(ninfo, ninfo); sigma.R <- sqrt(1/tau.R) #ランダム効果の精度の(超)事前分布
  
  # Likelihood
  for(i in 1:N){
    y[i] ~ dpois(mu[i])
    log(mu[i]) <- b[1] + b[2]*x[i] + R[group[i]]
  }
  
  for(j in 1:Ng){ R[j] ~ dnorm(0, tau.R) } 
}
model{
  ninfo <- 0.0001
  
  # Prior
  for(i in 1:2){ b[i] ~ dnorm(0, ninfo) }
  tau ~ dgamma(ninfo, ninfo); sigma <- sqrt(1/tau)
  tau.R ~ dgamma(ninfo, ninfo); sigma.R <- sqrt(1/tau.R) #ランダム効果の精度の(超)事前分布
  size ~ dnorm(0,ninfo)T(0,1000)

  # Likelihood
  for(i in 1:N){
    y[i] ~ dnegbin(p[i],size)
    p[i] <- size/(mu[i]+size)
    log(mu[i]) <- b[1] + b[2]*x[i] + R[group[i]]
  }
  
  for(j in 1:Ng){ R[j] ~ dnorm(0, tau.R) } 
}
###MCMC setting
n.iter <- 3E+4
n.thin <- max(3, ceiling(n.iter/500))
burn <- ceiling(max(10, n.iter/2))
n.ad <- 100
Sample <- ceiling(n.iter/n.thin)

# data for JAGS
Djags <- list(N = nrow(D), Ng = length(unique(D$group)), y = D$y, x = c(scale(D$x, scale=F)), group=D$group)

# parameters to be monitored
para <- c("b", "sigma", "sigma.R")
  
# set initial values
inits <- replicate(3, list(.RNG.name="base::Wichmann-Hill", .RNG.seed=NA), simplify = FALSE)
for(i in 1:3){ inits[[i]]$.RNG.seed <- 2*i+2 }

# Run JAGS
library(runjags)
m <- read.jagsfile("model_random.R")
post <- run.jags(m$model, monitor = para, data = Djags,
                 n.chains = 3, inits=inits, method = "parallel",
                 burnin = burn, sample = Sample, adapt = n.ad, thin = n.thin,
                 n.sims = 3)

結果の表示

RからBUGSソフトウェアを動かすのスクリプトを入力すると、次のような計算プロセスが表示され、postに計算結果が格納されます。

Compiling rjags model...
Calling the simulation using the rjags method...
Adapting the model for 100 iterations...
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100%
Burning in the model for 500 iterations...
|**************************************************| 100%
Running the model for 1500 iterations...
|**************************************************| 100%
Simulation complete
Calculating summary statistics...
Calculating the Gelman-Rubin statistic for 2 variables....
Finished running the simulation

次のコードを入力し、結果を確認します。

summary(post)
        Lower95    Median   Upper95     psrf
alpha -0.02896274 0.1647402 0.3390248 1.001341
beta   0.63089605 0.7559828 0.9100826 1.001867

Medianが推定結果の中央値で、Lower95とUpper95がいわゆる95%信用区間(CI)と呼ばれるものです(“信頼”区間ではない;意味が異なるので注意)。先の真の値(alpha = 0.1, beta = 0.8)と比べると、いずれ真の値も95%CIに含まれていることが分かります。また、これらの推定結果と合わせて、psrfというものも表示されていますが、これはパラメーターの推定値が十分収束しているかどうかを判定するための値(Gelman-Rubin statistics、もしくはR-hatと呼ばれる)です。この値が1.1を上回る場合は推定値が十分に収束していないので、run.jagas関数のburnin、sample、thinの値を適宜調整し、信頼できる推定値を得る必要があります。

RからBUGSソフトウェアを動かす

BUGSファイルの記述で作った統計モデルを使い、JAGSによるパラメーター推定をしてみます。まずは、真のパラメーターが分かっている確率分布から、適当なテストデータセットを作ります。

###R code###
N.sample <- 100 # sample size
Area <- rpois(N.sample,100) # forest area
X <- scale(Area)  # forest area

alpha <- 0.1 # intercept
beta <- 0.8  # slope

y <- lambda <- NULL
for(i in 1:N.sample){
 lambda[i] <- exp(alpha + beta*X[i])
 y[i] <- rpois(1,lambda[i]) # Poisson error process
}

ここで作られた鳥の個体数のデータセットは、標準化された森林面積Xに対し、切片0.1(α)、傾き0.8(β)という関係性をもっています。このデータをプロットすると、100地点分の鳥個体数と森林面積の関係は右図のようになります。

このデータセットからパラメーターを逆推定し、真の値と比較することで、推定結果の妥当性を判定できます。データセットをJAGSへ渡し、計算を実行します。

###R code###
#①パッケージrunjagsの読み込み(JAGSを動かすための関数がある)
library(runjags)

#②JAGSに渡すためのデータセット
Djags <- list(N.sample = N.sample, y = y, X = X)

#③推定結果を知りたいパラメーター
para <- c("alpha","beta")

#④使う統計モデルの読み込み
m <- read.jagsfile("model_PoissonRegression.R")

#⑤数値計算のための初期値の設定
inits <- replicate(3,
                   list(.RNG.name="base::Wichmann-Hill",
                   .RNG.seed=round(runif(1,1,100))),
                   simplify = FALSE)

#⑥関数run.jagasでJAGSを動かし、計算結果をpostに格納
post <- run.jags(model = m$model, monitor = para, data = Djags,
                 n.chains = 3, inits=inits, method = "parallel",
                 burnin = 500, sample = 500, adapt = 100, thin = 3,
                 n.sims = 3)

①パッケージrunjagsの読み込み

install.packages(runjags)と打つことでインストールできます。このなかに、read.jagsfile, run.jagsという関数が入っています。

②JAGSに渡すためのデータセット

モデルを記述したBUGSファイルで使われているデータセットを入力します。ここで、モデルで使われているデータが入力されていないと、エラーがでて計算がストップします。なお、JAGSの場合(ほかの言語もほぼ同様)、データの型に以下のような制限があります。

  • データセットはlist形式(ただし、BUGSソフトウェアによって異なる)
  • 変数の中身はスカラー、ベクター、マトリクス、あるいはアレイ形式(X = 1, 2, 3, … nのような純粋な数字の羅列;変数名などが入るとエラーが出る)
  • 因子型の変数(“Yes”など)は認識されないため、すべて数字に変換しなければならない

③推定結果を知りたいパラメーター

推定結果を知りたいパラメーターを宣言します。ここに含まれていないパラメーターの推定結果は表示されません。なお、パラメーターはモデルファイル(model_PoissonRegression.R)で定義されているパラメーター名と一致する必要があります。

④使う統計モデルの読み込み

使うモデルファイル(model_PoissonRegression.R)を指定します。

⑤数値計算のための初期値の設定

JAGSのベイズ推定の過程では、MCMC(マルコフ連鎖モンテカルロ)と呼ばれる逐次的な推定手法が用いられています。パラメーターには初期値を与えると、指定した統計モデルのもとであるパラメーターのフィッティングを評価し(尤度)、次のステップではよりいいパラメーターの値を探す、という計算過程を繰り返し行います。ベイズの推定結果は、この(数値)計算過程の結果、パラメーターが収束する「値の範囲」を示したものになります。通常、このMCMCという計算は一度の推定で3つ行うのですが、それぞれ違う初期値からスタートする必要があります(違う初期値から同じ値に落ち着くかどうかを確かめたいため)。ここのコードは、その初期値が違う値となるよう設定する「おまじない」です。

⑥関数run.jagasでJAGSを動かし、計算結果をpostに格納

関数run.jagsが実際に「JAGSを動かせ」という指令を出しており、この関数が入力された時点でJAGSの計算がスタートします。それぞれの引数について説明します。

  • model: 実際に使う統計モデルを指定する。ここでは、「m」というところにモデルファイルが格納されているので、その中のm$modelを渡す。
  • monitor: 推定結果を確認したいパラメータを指定する。③で作ったparaを渡す。
  • data: 推定に使うデータセット。②で作ったDjagsを渡す
  • n.chains: MCMCの数を指定する。通常は3つで、変更する必要はほぼない。
  • inits: MCMCで使う初期値。⑤で作ったinitsを渡す。
  • method: 計算の際の詳細指定。いくつかあるが、parallelとすることでパソコンの複数のコアをフル活用し、計算が高速化される。
  • burnin: MCMCの最初の何試行分を捨てるかを指定する。ここで指定されたMCMCサンプルは推定結果に反映されない。全試行数の30~50%で指定することが多い。
  • sample: 最終的な推定値の計算に用いるMCMC1つ当たりのサンプル数。500とすると、500*3で1500サンプルが推定に用いられる。500~1000くらいが妥当。多いほど推定結果が収束しやすくなるが、その分計算に時間がかかる。
  • adapt: burninよりも先になされる試行計算。推定結果には反映されない。100~1000で設定しておけば問題ない。
  • thin: MCMCサンプルを「何個おき」に取ってくるかを指定する(その間のサンプルは「間引き」され、推定結果に反映されない)。MCMCの計算が自己相関(次のステップの値が前の値と類似する)するときには、間隔を広げたほうがいい。最低3は必要で、複雑なモデルになると1000以上の値を設定する必要がある場合もある。間隔をあけるほど全試行数は増えるので、計算時間は長くなる。
  • n.sims: 計算に使うCPUコア数の指定?のようなもの。chainの数に合わせる。

続きはこちら

BUGSファイルの記述

まずはBUGS言語による統計モデルの記述から説明を進めます。BUGSファイルは「統計モデルの構造」のみを記述するファイルです。ここでは、森林パッチiの鳥の個体数Yiを、森林面積Xiで説明するポアソン線形モデルを例に挙げて考えます。

Eq. 1)
Yi ~ Pois(λi)
log(λi) = α + β*Xi

上記のモデルは、RのGLMでは次のように記述されます。

###R code###
data <- read.csv(“data.csv”)
glm(y ~ X, data, family=poisson)

ここで気を付けなければならないのは、

  • GLMのようにRでパッケージ化されているものだと、eq.1で書かれているほとんど部分が自動化されており、リンク関数(log)などを自分で考える必要がない
  • 特に指定しなくとも、サンプルの繰り返しが“i”として認識されている

ことです。これら2点がBUGSとR言語では大きく異なります。このことを念頭におき、上記のモデルをBUGS言語で書いてみます。

###BUGS code###
model{
#Likelihood
 for(i in 1:N.sample){ # loop for i
   y[i] ~ dpois(lambda[i]) # poisson model
   log(lambda[i]) <- alpha + beta * X[i] # linear model structure
 }

#Prior
 alpha ~ dnorm(0, 0.001)
 beta ~ dnorm(0, 0.001)
}

model{}で囲われている範囲がモデル構造を指定するエリアになります({}の外にテキストを書くとエラーがでるので注意)。「~」は確率論的ノード(Stochastic node)と呼ばれ、左辺が確率変数である場合に使われます。確率論的ノードの右辺には、必ず確率分布の記述が必要です(上記ではdpoisが確率分布、Poisson distributionを指すもの)。「<-」は決定論的ノード(Deterministic node)と呼ばれ、=(イコール)と対応します。こうしてみると、BUGS言語ではEq.1をすべて自分で書き下す必要があると分かります。 Prior以降の記述は「事前分布」と呼ばれるものです。ベイズ統計の枠組みでは、すべてパラメーターに対し、事前分布を設定しなければなりません(事前分布を設定し忘れるとエラーが出る)。ここでは切片αと傾きβに対し、事前の情報は何も得られないので、「無情報事前分布」というものを設定しています(ほとんどケースはこれにあたる)。ここでは、平均0、分散1000の正規分布を事前分布として設定することで、切片αと傾きβは「かなり広い範囲の値をとる可能性がある」としています。「dnorm」が正規分布であることを示しており、一つ目の引数が「平均」、二つ目の引数が「精度(分散の逆数=1/1000)」を表します。ここまで書き終えたら、ファイルを.Rもしくは.txtファイルとして保存します(ここでは“model_PoissonRegression.R”として保存することにします)。続きはこちら

BUGSの使い方

最近ではRパッケージも実装され、容易にRからBUGSソフトウェアを使うことができます。その一連の流れは以下のようになっています。

  • 扱いたい統計モデルをBUGS言語で記述し、.Rもしくは.txtファイルとして保存
  • 解析に使うデータをRで準備(データ整形のしやすさはRのほうが優れているため)
  • RからBUGSソフトウェアにデータを渡し、指定した統計モデルのもとで計算を実行
  • BUGSからRに渡される計算結果を保存

この際、「どんな統計モデルを計算させるのか」という部分だけ、RではなくBUGS言語で記述する必要があります。したがって、統計モデルの実装の際には、RおよびBUGSスクリプトをそれぞれ別個に用意することになります。

ベイズ統計とBUGS

ベイズ統計、と聞くとなにやらすべてが新しいものに聞こえますが、扱うモデルの「構造」自体は大きく変わりません(線形回帰、ANOVAなど)。ただし、回帰係数のような“パラメーター”の扱いと推定手法が大きく異なっており、ベイズ統計の実装にはマルコフ連鎖モンテカルロのような逐次的な推定法が必要です。自前でアルゴリズムを構築することも可能ですが、そんな煩雑なプログラムを自分で書くというのは現実的な選択肢ではないでしょう。

そこで、BUGSソフトウェアの登場です。BUGSとはBayesian inference Using Gibbs Samplerの略で、ギブスサンプリングによるベイズ統計の実装を指します。このギブスサンプリングを実行してくれるのがWinBUGSをはじめとするBUGSソフトウェアで、Rから動かすことができます。そのため、昨今のベイズ統計を用いた研究では、RとBUGSソフトウェアの連携が一般的な手法になっています。なお、WinBUGSのほかにも、OpenBUGSとJAGS(Stanもベイズ実装のためのソフトですが、ギブスサンプリングではない)があります。WinBUGSは開発が終了していますが、OpenBUGSとJAGSは最新バージョンが随時ネット上でダウンロードできるようになっており、計算効率もそれらのほうがより良いものとなりつつあります。ここでは、BUGの中ではもっともクラッシュしにくいJAGS(Just Another Gibbs SamplerなのでBUGSというには語弊がありますが。。。)について解説します。

JAGSのダウンロードは次のページから可能です(JAGS)。ダウンロードはProgram files(デフォルト)もしくはCドライブ直下がお勧めです。