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の数に合わせる。

続きはこちら

SPONSOR LINK

Spread the love

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です