[17/08/08 09:44 更新]
最低限のことが出来るような、間に合わせのパッケージを作った。このtar ballには、Perl Package、サンブルスクリプト、それに使うデータが含まれている。
[2017年8月8日追記: ここで紹介する間に合わせのパッケージの次のバージョンが「PerlとJAGSで赤池情報量基準(AIC)」で紹介されている]
このパッケージを使うと、Perl ScriptからJAGSを呼び出してベイズ推定をし、R言語のcodaパッケージを使って結果を解析できる。どちらの機能もフルに使ってるわけではない。codaを使った収束判定は、今のところ、Gewekeの指標しかサポートしてない。このパッケージを叩き台にして、自由に変更・拡張して欲しい。
ここでは、「JAGS入門」でやってることを、おおむね再現することを目指す。
まず、データの準備をする。ここでは、参考ページのデータをそのまま流用している。
x,y
0,72
10,66
20,60
30,58
40,56
50,55
60,52
70,48
これを、regdata.csvという名前で保存しておく。データ形式はCSVである。「#」で初まる行は、コメントとして無視する。コメント以外の最初の行に変数名のリストが書かれる。そのあとに対応するデータがくる。欠損値がある場合は該当個所にNAと書く。変数名等の文字列は、二重引用符でくくってもよいし、くくらなくてもよい。カンマのあとに、一つもしくは複数の空白があってもよい。手抜きの実装なので、変数名にカンマを含むと動作しない。
次に、モデルを用意する。これもそのまま流用。
model{
for(i in 1:N){
y[i] ~ dnorm(mu[i], tau)
mu[i] <- alpha[1] + alpha[2]*x[i]
}
for(i in 1:K){
alpha[i] ~ dnorm(0, 1.0E-6)
}
tau ~ dgamma(1.0E-6, 1.0E-6)
sigma2 <- 1/tau
}
これを、reg_model.jagという名で保存。ほぼ無情報事前分布で直線のモデルになっている。
以下に、サンプルスクリプトtest.plをもとに、本パッケージの使い方を説明する。
パッケージを使うには、まず、
use Makeshift::Jags;
が必要である。ここでは、@INCに含まれるディレクトリにMakeshiftディレクトリがあることを仮定しているが、そうでない場合は、
use lib "$ENV{HOME}/lib/perl5";
などとするか、環境変数PERL5LIBを設定するなどの工夫が必要である。
Jagsインスタンスを作るには、
my $jags = Makeshift::Jags->new;
とする。デフォルトでは一時ディレクトリが/tmpに作られるが、別の場所に作りたい場合は、例えば、
my $jags = Makeshift::Jags->new(tmproot => "$ENV{HOME}/tmp");
などのようにする。指定されるディレクトリは、あらかじめ作っておくこと。一時ディレクトリや一時ファイルの自動削除を防ぎたければ、
my $jags = Makeshift::Jags->new(CLEANUP => 0);
とする。一時ファイルは全てテキストファイルなので、何かうまくいかないときは、覗いてみるとよい。
CSVファイルの読み込みは、
my $n = $jags->read_csv("regdata.csv");
である。返り値は、読み込まれたデータ点数となる。
定数の設定はdata->setで、推定されるパラメータの初期値の設定はparam->setで行う。
$jags->data->set(N => $n, K => 2);
$jags->param->set(alpha => [1, 1], tau => 1);
モデルの読み込みは、
$jags->model_file("reg_model.jag");
である。
デフォルトでは、初期値を指定したパラメータを全てモニターするようになっているが、ここでは、tauのかわりにsigma2をモニターする。
$jags->set_monitor("sigma2");
$jags->del_monitor("tau");
あとは、バーンインのための繰り返し数と本番のサンプリング回数を指定すれば、jagsを実行できる。
$jags->burnin(1000);
$jags->update(5000);
$jags->run;
結果の解析にはcodaを使う。インスタンスの生成は、
my $coda = $jags->coda->[1];
といった形。最後の[1]は、チェインの番号を指定している。ここでは、チェイン数は1(デフォルト)で実行しているので、1番目だけが必要。
収束判定をしたければ、
$coda->diag_method("geweke");
とする。デフォルトでは判定しない。今のところ、Gewekeの指標しかサポートしてない。他の判定方法が必要であれば、Makeshift/Jags/coda以下のファイルを参考に、自作してほしい。
結果の概要は、
$coda->summary;
で得られる。モニターされたパラメータについての情報が簡潔に表示される。収束判定をする場合は、それぞれのパラメータに対する指標と、収束してるかどうかの判断も表示される。
# var, mean, sd, naive_se, time_series_se, geweke
alpha[1], 69.0327, 1.58311, 0.0223886, 0.053710, -1.476
alpha[2], -0.3038, 0.03834, 0.0005421, 0.001338, 1.603
sigma2, 6.0491, 7.10036, 0.1004142, 0.162382, -1.123
All parameters have converged
スクリプトの最後の部分でグラフを描ける。PDLとPDL::Graphics::Gnuplotが必要なのでコメントアウトしてるが、#を外せば使えるはずである。実行結果を活用する例にもなっているので、参考にしてほしい。
他に何ができるかは、パッケージのソースで判断してほしい。足りない部分は、自分で作ってほしい。
[「インストール」に戻る] [「PerlとJAGSでマルコフ連鎖モンテカルロ」に戻る] [「あえてのPerl」に戻る]