せっかくなので触ってみました。
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;
こんな感じに略記することもできる。
指示変数を入れて群別にパラメータを出すときには便利ですね。
絵もきれいに出る。
