ガッキーより竹内結子

なのでダンダリン見てる。


あっ, どーも僕です。

リアップは劇薬じゃなくてよかったです。


rstanでちょこちょこ

こちらのブログにrstanの導入&入門編が紹介されていました。

MCMCの計算にStanを使ってみた(超基礎・導入編) - 銀座で働くデータサイエンティストのブログ

普段から読んでいるブログで毎度毎度参考になります。
ですが、はじめてrstanを触る人には誤解を与えるかもな内容でしたので、ちょっと付け足しておきます。

付け足しポイントは次の2つ

  • 回帰係数にはベクトル演算をつかうstanはいっぱい関数があってよくわからないので, まかせます
  • codaを使わなくてもトレースと事後分布は出せる

ベクトル演算は公式リファレンスで推奨されていますし、ひとつづ係数を書くのは大変なのでオススメしません。

書いて説明するのは大変なので、とりあえずRのスクリプトはこちら。リファレンスに似たような例題が載っていますので、それを真似ています。
rstanをロードしてからこれをそのまま貼れば、前述のブログと同じ結果を得られます。
細かいメモはスクリプトに書いてあるので気になったら見てください。

# *************************************************************************
# http://tjo.hatenablog.com/entry/2013/11/06/201735
# ロジスティック回帰の例題
# *************************************************************************
# データ
dat <- read.table ("https://raw.github.com/ozt-ca/tjo.hatenablog.samples/master/r_samples/public_lib/jp/conflict_sample.txt", header = T)
cv  <- c(dat[,8]) - 1
X   <- dat[,-8]; 
N   <- nrow (X)
M   <- ncol (X)

stan_code <- '
    data {
        int<lower=0> N;
        int<lower=0> M;
        matrix[N, M] X;
        int<lower=0, upper=1> y[N];
    }
    parameters {
        real beta0;
        vector[M] beta; // 7をMに修正(2013/11/7)
    }
    model {
        for (i in 1:N)
            // X[i] は row_vector, beta は vector だが, dot_product が吸収してくれる
            y[i] ~ bernoulli(inv_logit (beta0 + dot_product(X[i] , beta)));
    }
'

# *************************************************************************
# MCMC
fitchan <- stan (model_code = stan_code, # file でファイル読み込み, fitでコンパイル済みもの
                 data = list (N = N, 
                              M = M,
                              X = X, 
                              y = cv), 
                 iter  = 1000, 
                 chain = 4)

# *************************************************************************
# グラフ
# トレース
# codaのtraceplotとブッキングするので、codaは読み込ませないように
traceplot (fitchan)

# 事後分布
library (reshape)
library (ggplot2)
post   <- extract (fitchan, permuted = F)
m.post <- melt (post)
graph  <- ggplot (m.post, aes(x = value))
graph  <- graph + geom_density () + facet_grid(. ~ parameters, scales = "free") + theme_bw() #facet_wrapのほうがいいかも
plot (graph)

グラフは全部の変数をだすとおっきくて見にくいので好き変数を抽出してください。

結果の統計量はこちら。

 fitchan
Inference for Stan model: stan_code.
4 chains, each with iter=1000; warmup=500; thin=1; 
post-warmup draws per chain=500, total post-warmup draws=2000.

          mean se_mean  sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
beta0     -1.4     0.0 0.2   -1.9   -1.6   -1.4   -1.2   -0.9  1229    1
beta[1]    1.1     0.0 0.2    0.7    1.0    1.1    1.2    1.4  1615    1
beta[2]   -0.6     0.0 0.2   -0.9   -0.7   -0.6   -0.4   -0.2  1724    1
beta[3]    0.1     0.0 0.2   -0.2    0.0    0.1    0.2    0.4  1659    1
beta[4]   -3.0     0.0 0.2   -3.4   -3.2   -3.0   -2.9   -2.6  1295    1
beta[5]    1.5     0.0 0.2    1.2    1.4    1.5    1.6    1.9  1580    1
beta[6]    5.4     0.0 0.2    5.0    5.3    5.4    5.5    5.8  1255    1
beta[7]    0.1     0.0 0.2   -0.2    0.0    0.1    0.2    0.4  1724    1
lp__    -525.9     0.1 1.9 -530.3 -526.9 -525.6 -524.5 -523.2   744    1

Samples were drawn using NUTS(diag_e) at Wed Nov 06 22:31:29 2013.
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).

わたしも一つ入門編を

わたしもstan勉強し始めて三週間あまりなので、初心者の域をでていないのですが、入門ネタを一つ。

ユーザー定義の関数

わたしはrstanを勉強するにあたり、ネットに落ちているコードや、公式HPの例題(BUGSコードをstanに翻訳したもの)をとりあえず動かせるように試してみました。

そのなかで、たぶんリファレンスを読んでいないと詰まるところを紹介。
コードは石倉究さんの逆ガウス分布から
データ解析環境Rの整備と利用 - 2012年度 データ解析環境Rの整備と利用

細かい式は資料を見てください。
石倉さんの資料ではユーザー定義の対数尤度は次の通りと書いてあります。

model {
lp__ <- lp__ + 0.5 * N * log(lambda / (2.0 * pi())) - 1.5 * sum_log_y;
for(i in 1:N)
    lp__ <- lp__ - lambda * pow((y[i] - mu) / mu, 2) / (2.0 * y[i]);
}


Stan2ではどうやらこの書き方だと動作が安定しないので次のようにincrement_log_probを使います。
かなり汚いですが、安定して動きます。increment_log_probの引数には対数尤度をいれておきます。

model {
    for (i in 1:N)
        increment_log_prob (0.5 *  log (lambda / (2 * pi() * pow(y[i], 3))) - lambda * pow ((y[i] - mu) / mu, 2) / (2 * y[i]));
}

おまけ

ちなみに、上のロジスティック回帰の結果をグラフにするとこんな感じです。


f:id:aaaazzzz036:20131106225328j:plain
f:id:aaaazzzz036:20131106225346j:plain