読者です 読者をやめる 読者になる 読者になる

mcmc procedure

せっかくなので触ってみました。

data logistic;
  call streaminit(88745623);
  a = 0; b = 1;
  do i = 1 to 200;
    do x = 0 to 1;
    logit = a +b * x;
    p = exp(logit) / (1 + exp(logit));
    y = rand("Bernoulli", p);
    output;
  end; end;
  keep x y;
run;

ods listing close;
ods pdf file = "c:/mcmc.pdf" style = Journal;
ods graphics on;

/* logistic regression */
proc logistic data = logistic;
  model y(event="1") = x;
run;

/* bayesian logistic regression */
proc mcmc data = logistic
  ntu=500 nbi=1000 nmc=20000 nthin=1 seed=86574
  propcov=quanew diag=all outpost=dpost;
  parms (a b) 0;
  prior a ~ normal(0, var=10000);
  prior b ~ normal(0, var=10000);
  p = logistic(a + b * x);
  model y ~ binomial(1, p);
run;

ods graphics off;
ods pdf close;
ods listing;

サンプリングは Metropolis-Hasting method で実装されている、ntu=500 nbi=1000 nmc=20000 nthin=1 あたりで条件を指定できる。
prior の指定は、nlmixed procedure の様に、general(xxxx) 形式が使えるため、事前分布の確率密度関数を自分で書く事もできる。
hyper statement や hyperprior statement で hyperprior まで指定できてしまう。
文法が nlmixed ライクなので、if 文でモデルを超複雑に分岐したりできてしまう。
これは恐ろしい procedure だ、書き方を覚えたら何でもできちゃうじゃないか。


あと、パラメータが連番の場合、

proc mcmc data = logistic
  ntu=500 nbi=1000 nmc=20000 nthin=1 seed=86574
  propcov=quanew diag=all outpost=dpost;
  array beta[2];
  parms beta: 0;
  prior beta: ~ normal(0, var=10000);
  p = logistic(beta1 + beta2 * x);
  model y ~ binomial(1, p);
run;

こんな感じに略記することもできる。
指示変数を入れて群別にパラメータを出すときには便利ですね。


絵もきれいに出る。


出力例 (Dropbox)