Stanの練習がてら、おなじみのirisデータでロジスティック回帰のパラメータを推定してみた。
今回作成したモデルはirisのSpeciesを予測する多クラスロジスティック回帰。回帰係数の事前分布は平均0の正規分布とした（こうするとL2正則化と同等の効果が得られる）。
なお、正規分布の分散の事前分布には半コーシー分布を使う。サンプル数が少ないときは逆ガンマ分布を無情報事前分布として使うと不都合があるらしいので。（でもirisのサンプル数は150なので逆ガンマ分布のままでも大丈夫かも？） 追記: Gelmanの論文をチラ見したら逆ガンマ分布を使って問題になるケースは「グループレベルの効果をモデル化していて、グループ数が少ないあるいはグループレベルの分散がゼロに近い場合」みたいなので、今回のケースは当てはまらないかも（別にコーシー分布を使っても問題は無いみたいだけど）。
参考: http://www.slideshare.net/KojiKosugi/cauchy20150726
以下コード。
require(rstan)
y <- as.integer(iris[,5])
x <- scale(iris[,-5])
data.list <-
list(
N = NROW(x),
M = NCOL(x),
K = 3,
y = y,
x = x)
stan.code <- "
data {
int<lower=0> N;
int<lower=0> M;
int<lower=0> K;
int<lower=1,upper=3> y[N];
matrix[N,M] x;
}
parameters {
real<lower=0> sigma;
vector[K] alpha;
matrix[K,M] beta;
}
model {
sigma ~ cauchy(0, 1);
for (i in 1:K) {
alpha[i] ~ normal(0, sigma);
}
for (d in 1:M) {
for (i in 1:K) {
beta[i,d] ~ normal(0, sigma);
}
}
for (n in 1:N) {
y[n] ~ categorical(softmax(alpha + beta * x[n]'));
}
}
"
stan.fit <-
stan(
model_code = stan.code,
data = data.list,
iter = 1000,
thin = 10,
chains = 3)
print(stan.fit)
実行結果。
> print(stan.fit)
Inference for Stan model: 45dc1db582ed31d42ac02f63763415d9.
3 chains, each with iter=1000; warmup=500; thin=10;
post-warmup draws per chain=50, total post-warmup draws=150.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
sigma 8.97 0.41 4.45 3.87 5.79 7.74 10.61 19.80 120 0.99
alpha[1] -3.01 0.81 9.98 -24.72 -6.89 -1.33 2.29 9.66 150 0.98
alpha[2] 8.49 0.61 7.48 -5.76 4.50 8.40 12.56 22.48 150 1.00
alpha[3] -6.67 0.66 8.13 -27.79 -9.76 -5.04 -1.95 6.53 150 1.00
beta[1,1] -3.71 0.69 8.44 -23.82 -8.17 -3.68 1.07 10.45 150 1.01
beta[1,2] 5.26 0.65 8.01 -7.40 0.91 4.13 8.77 21.17 150 0.99
beta[1,3] -9.71 0.83 10.18 -41.41 -13.96 -8.47 -3.49 5.19 150 0.99
beta[1,4] -9.22 1.18 10.99 -37.98 -13.46 -8.59 -2.95 8.88 86 1.00
beta[2,1] 2.94 0.53 6.36 -9.38 -0.09 2.51 6.31 15.99 143 0.98
beta[2,2] -1.54 0.51 6.20 -14.04 -4.54 -2.06 1.23 12.04 150 1.00
beta[2,3] -1.29 0.45 5.18 -10.75 -4.87 -1.43 1.76 8.60 134 1.00
beta[2,4] -0.66 0.53 6.50 -14.73 -3.72 -0.65 1.84 13.32 150 1.00
beta[3,1] 1.02 0.51 6.10 -12.27 -1.86 0.66 4.72 14.01 145 0.98
beta[3,2] -3.58 0.52 6.34 -17.62 -6.64 -3.81 -0.08 10.62 150 1.00
beta[3,3] 11.49 0.59 6.91 0.37 6.66 10.36 15.54 26.33 138 0.99
beta[3,4] 10.38 0.64 7.83 -3.65 5.64 9.51 14.63 26.99 150 0.99
lp__ -51.00 0.58 6.13 -63.69 -54.32 -50.18 -46.85 -40.07 112 0.99
Samples were drawn using NUTS(diag_e) at Sun Jan 31 10:59:40 2016.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).
上記コードは、回帰係数の事前分布に平均0の正規分布を置いているので、実質的にL2正則化ロジスティック回帰と同等のことをしてくれる。しかも凄いのは、普通L2正則化ロジスティック回帰を使うときは正則化パラメータをcross-validationとかで決めないといけないのだけれど、上記コードでは正規分布の分散（≒正則化パラメータ）の事前分布に無情報事前分布を指定することで、正則化パラメータまで含めて一度にデータから推定させることができるのだ。これが階層ベイズの力だ！すげー！Σ(ﾟдﾟ;)
でもまぁ、この程度だったら普通にglmnetとか使ったほうがいいかも。MCMCは遅いしね…。