引き続き, データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)
MCMC
MCMC (Markov chain Monte Carlo methods)については, 詳しい情報が書籍 [1]や Webに溢れていますので, そちらを参照ください。
MCMCは事後分布からのランダムサンプリングを得るための道具ですが, パラメータ推定にも使うことができます。
ベイズ統計では, 普通の統計のようにパラメータをただひとつの真の値が存在するとは考えずに確率変数とします。
単純にベイズ統計によってパラメータを推定するには, 多変量の事後分布の期待値や周辺分布の計算, 多次元積分が必要になってしまいますが, MCMCはこれを現実的に行うための計算手法に関する1990年代に生まれたイノベーションです。
Stan
Stanは, C++で実装された確率的プログラミング言語です。DSLとも呼べるかもしれません。注目はMCMCサンプラーの中でも, Hamiltonian Monte Carloを No-U-Turnで実装している点 [2]です。Stanを使うことで, 容易に確率的プログラミングを行えます。
下記のような特徴があります。
- コンパイルによる高速に動作
- NUTSにより収束が速いことが期待できる
- Docsが丁寧である
- 各言語へのインターフェイスが用意されている
多数のパラメータを扱いたい時に向いてそうです。
RStanはRからStanへのインターフェイスです。他にもPyStanなど各言語へのBindingがあります。
RStanのインストール
環境はOSX 10.10.2です。
最新のインストール方法は RStan Getting Started を覗くのが確実です。RStan 2.8.0はCRANに登録されているので, 下記でインストールできます。
Sys.setenv(MAKEFLAGS = "-j4")
install.packages("rstan", dependencies = TRUE)
下記手順はRStan 2.6.0のインストール方法です。
R 3.0.2未満でRStan 2.6.0をインストールしようとすると, エラーが発生するので先にRをupdateしておきます。
$ ruby -e "$(curl -fsSL https://raw.githubusercontent.com/Homebrew/install/master/install)"
$ brew tap homebrew/science
$ brew install r
Rのconsoleからrstanをインストールします。
> install.packages("inline")
> install.packages("Rcpp")
> Sys.setenv(MAKEFLAGS = "-j4")
> source('https://mc-stan.org/rstan/install.R', echo = TRUE, max.deparse.length = 2000)
> install_rstan()
g++がないと怒られるので, Makeconfでコンパイラ設定でclangを使うように書き換えます。
/bin/sh: llvm-g++-4.2: command not found
$ cd /Library/Frameworks/R.framework/Resources/etc
$ vim Makeconf
CC=clang
CXX=clang++
CXXFLAGS= -O3 -pedantic
# CC = llvm-gcc-4.2 -arch x86_64 -std=gnu99
# CFLAGS = -mtune=core2 -g -O2 $(LTO)
# CPICFLAGS = -fPIC
# CPPFLAGS = -I/usr/local/include
# CXX = llvm-g++-4.2 -arch x86_64
# CXXCPP = $(CXX) -E
# CXXFLAGS = -mtune=core2 -g -O2 $(LTO)
# CXXPICFLAGS = -fPIC
...
再度インストールを試みると今度は成功しました。
> require(rstan)
要求されたパッケージ rstan をロード中です
要求されたパッケージ Rcpp をロード中です
要求されたパッケージ inline をロード中です
次のパッケージを付け加えます: ‘inline’
The following object is masked from ‘package:Rcpp’:
registerPlugin
rstan (Version 2.6.0, packaged: 2015-02-06 21:02:34 UTC, GitRev: 198082f07a60)
Windows環境では, Install Rtools for Windowsから RToolsをインストールし, CRANから rstanをインストールします。私の環境だと Sys.getenv(‘PATH’)で確認したところ, RtoolsにPATHが通ってなかったので追加しました。
rtools <- "C:\\Rtools\\bin"
gcc <- "C:\\Rtools\\gcc-4.6.3\\bin"
path <- strsplit(Sys.getenv("PATH"), ";")[[1]]
new_path <- c(rtools, gcc, path)
new_path <- new_path[!duplicated(tolower(new_path))]
Sys.setenv(PATH = paste(new_path, collapse = ";"))
Example 0 : 正規分布の平均と分散の推定
まずは 二番目に簡単な rstan コード を動かしてみます。
Stanではブロック単位で処理を記述していきます。
- data : RやPythonから入力するデータの定義
- parameters : 推定するパラメータ
- transformed parameters : dataで宣言したデータの変換 / ハイパーパラメータ定義
- model : モデルの定義
require(rstan)
N <- 1000
x <- rnorm(N, mean=50, sd=10)
y <- 10 + 0.8 * x + rnorm(N, mean=0, sd=7)
stancode <- '
data {
int<lower=0> N;
real x[N];
real y[N];
}
parameters { // 推定するパラメータ
real alpha; // 平均の第一項 (10)
real beta; // 平均の第二項の定数 (0.8)
real<lower=0> s; // 標準偏差 (7)
}
model{
for(i in 1:N)
y[i] ~ normal(alpha + beta * x[i], s); // yを正規分布で推定
alpha ~ normal(0, 100); // alphaを正規分布で推定
beta ~ normal(0, 100); // betaを正規分布で推定
s ~ inv_gamma(0.001, 0.001); // sを逆ガンマ分布で推定
}
'
datastan <- list(N=N, x=x, y=y) # x,y によるデータセット
fit <- stan(model_code=stancode, data=datastan, iter=1000, chain=4) # stanを呼ぶ
traceplot(fit, ask=T)
print(fit, digit=2)
</lower=0></lower=0>結果は以下です。alphaだけmeanの推定値のずれが大きい。
Inference for Stan model: stancode.
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
alpha 11.46 0.05 1.10 9.29 10.70 11.50 12.21 13.52 540 1
beta 0.77 0.00 0.02 0.73 0.76 0.77 0.79 0.81 537 1
s 6.79 0.01 0.15 6.51 6.68 6.79 6.89 7.10 572 1
lp__ -2416.25 0.06 1.25 -2419.55 -2416.76 -2415.90 -2415.35 -2414.87 470 1
Rhat[4]は尤度計算が収束したかの指標で1.1未満であることが基準になります。また, DIC 逸脱度モデル選択基準を使うのひとつの手かもしれません。収束しない場合, モデルやサンプリングに問題がありそうです。
95%ベイズ信頼区間 (確信区間, 信用区間)に0を含まないこともパラメータの妥当性を測る上で基準のひとつになります。
また, lp_(log_prob)は事後確率の対数スケールで通常は自動修正されるようですが, increment_log_prob()で修正できるようです。
Example 1 : Eight Schools
続いて, RStanのGetting-Started Example 1を試してみます。
Eight Schools[5]は米国のSAT (大学進学適性試験)スコアのコーチング効果の分析データです。
schools_code が stan codeになります。model_codeに指定していますが, 別ファイルに書いて読み込むこともできます。
require(rstan)
schools_code <- '
data {
int<lower=0> J; // サンプルサイズ
real y[J]; // コーチング効果
real<lower=0> sigma[J]; // 標準偏差
}
parameters {
real mu; // 定数
real<lower=0> tau; // eta[j] の係数
real eta[J]; // 個体差
}
transformed parameters {
real theta[J];
for (j in 1:J)
theta[j] <- mu + tau * eta[j]; // 個体ごとの平均
}
model {
eta ~ normal(0, 1); // eta は標準正規分布に従う
y ~ normal(theta, sigma); // y は N(theta, sigma)に従う
}
'
schools_dat <- list(
J=8,
y=c(28,8,-3,7,-1,1,18,12),
sigma=c(15,10,16,11,9,11,10,18)
)
fit <- stan(model_code = schools_code, data = schools_dat, iter = 1000, chains = 4)
fit2 <- stan(fit = fit, data = schools_dat, iter = 10000, chains = 4)
print(fit2)
plot(fit2)
la <- extract(fit2, permuted = TRUE) # return a list of arrays
mu <- la$mu
### return an array of three dimensions: iterations, chains, parameters
a <- extract(fit2, permuted = FALSE)
### use S3 functions as.array (or as.matrix) on stanfit objects
a2 <- as.array(fit2)
m <- as.matrix(fit2)
print(fit, digits = 1)
</lower=0></lower=0></lower=0>結果は以下です。収束は問題ないようです。
Inference for Stan model: schools_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
mu 7.8 0.3 5.4 -2.3 4.4 7.6 11.0 19.3 253 1
tau 6.8 0.3 6.0 0.2 2.4 5.3 9.3 22.4 380 1
eta[1] 0.4 0.0 0.9 -1.4 -0.2 0.4 1.0 2.2 1148 1
eta[2] 0.0 0.0 0.9 -1.8 -0.6 0.0 0.6 1.7 854 1
eta[3] -0.2 0.0 0.9 -2.0 -0.9 -0.2 0.5 1.7 1197 1
eta[4] 0.0 0.0 0.9 -1.6 -0.6 0.0 0.6 1.7 1161 1
eta[5] -0.3 0.0 0.9 -2.0 -0.9 -0.3 0.2 1.4 948 1
eta[6] -0.2 0.0 0.9 -1.9 -0.7 -0.2 0.4 1.5 856 1
eta[7] 0.4 0.0 0.9 -1.5 -0.2 0.4 0.9 2.0 1456 1
eta[8] 0.1 0.0 0.9 -1.9 -0.5 0.1 0.7 1.9 1248 1
theta[1] 11.4 0.3 8.2 -2.4 6.2 10.2 15.6 31.3 903 1
theta[2] 7.7 0.2 6.3 -4.6 3.7 7.7 11.7 20.6 1525 1
theta[3] 6.0 0.2 7.8 -12.1 1.9 6.7 10.7 20.5 1296 1
theta[4] 7.6 0.2 6.6 -7.3 4.0 7.8 11.4 21.1 1295 1
theta[5] 4.9 0.2 6.6 -10.4 1.1 5.5 9.3 16.3 819 1
theta[6] 5.9 0.2 6.6 -7.9 2.0 6.4 10.4 17.8 1451 1
theta[7] 10.7 0.2 6.9 -1.6 6.0 10.0 14.5 25.9 898 1
theta[8] 8.6 0.4 8.0 -5.9 3.7 8.1 12.9 26.6 476 1
lp__ -4.8 0.1 2.6 -10.5 -6.4 -4.6 -3.0 -0.5 574 1
戻り値のオブジェクト stanfit は可視化を行うplot(), traceplot()などを持っています。
おわりに
階層ベイズモデルでは, 似たようなパラメーターに制約を持たせることで, 一般化線形モデルでは難しかったデータの個数をパラメーターの個数を上回った場合においても上手く表現できるモデルを構築することができます。
複数のランダム効果が含まれる場合のパラメータ推定には, 特に Stan などのMCMCサンプラーを用いたパラメータ推定が良さそうです。
[1] 岩波データサイエンス Vol.1
[2] The No-U-Turn Sampler: Adaptively Setting Path Lengths in Hamiltonian Monte Carlo
[3] 数学セミナー 2007.11 554号
[4] R^ : 3本以上のサンプル列を比較して評価する。1.1以上はサンプル列間のばらつきが大きいので、定常分布は推定できないと判断。
[5] Dataset as seen in Rubin (1981) which was an analysis of coaching effects on SAT scores from eight schools.
[6] Stanで統計モデリングを学ぶ(4): とりあえず階層ベイズモデルを試してみる（基本編） – 銀座で働くデータサイエンティストのブログ
[7] StanTutorial
[8] Hamiltonian Monte Carloによるベイズ推定ソフトウェアStanの紹介
[9] Stanでincrement_log_prob()
