abcのABC (R Advent Calendar 2011 11日目)

近似ベイズ計算 (ABC: Approximate Bayesian Computation) に基づいたパラメータ推定のできるパッケージabcを紹介します．

ABCとは

ABCは，尤度関数の計算が不可能・困難なときに，尤度関数を求めることなく事後分布を推定できる手法です．比較的新しい手法で，集団遺伝学・生態学，疫学などさまざまな分野で応用と開発が進んでいます．
ABCに関しては，kos59125 さんの以下のエントリなどがとても参考になります．

近似ベイズ計算によるベイズ推定 (slideshare)
近似ベイズ計算によるベイズ推定 (捨てられたブログ)

パッケージabc

パッケージabcは，

観察データから計算した要約統計量 (引数 target)
事前分布からサンプリングしたパラメータ (引数 param)
それらをもとにシミュレートした要約統計量 (引数 sumstat)

を与えると，観察データと十分に近いシミュレーションセットを許容して，それらのパラメータを事後分布として返してくれます．要約統計量やパラメータが複数ある場合には，それぞれベクターやマトリックスの型にします (行が各シミュレーション，列が各要約統計量・パラメータ)．
許容範囲の大きさは，引数tolで指定します．
引数methodでは3種のアルゴリズムを選択できます．内訳は以下のとおり．

"rejection" 単純にもっとも近いn%のセットを許容
"loclinear" 局所線形回帰によりパラメータを重み付け
"neuralnet" 非線形回帰によりパラメータを重み付け

引数はほかにもありますが，とりあえず重要なのはこれくらいでしょうか．

作者自身の手による論文 [1] もありますので，ヘルプドキュメントとあわせてご参照ください．

デモ

では，実際にデモをしてみましょう！
5つの値からなる観察データ $\{D_{obs} | 2.2, 3.5, 2.9, 3.6, 2.8\}$ があったとします．簡便のため，正規分布にしたがうことがわかっており，分散 $\sigma ^2$ は既知で~~0.025~~ 0.25 (訂正 111213) とします．近似ベイズ計算のフレームワークで，平均値を算出してみます．真の平均値 $\hat{\theta}$ は3.0としました．

まずこの場合，データを生じるモデルは $f(D|\theta) \sim ND(\theta, 0.25)$ ということになります．
モデルのパラメータ $\theta$ には，事前分布として，平均値0で分散100の正規分布でも与えておきましょう．
データをあらわす要約統計量 $S$ には，5つの値の平均値を用います．
シミュレーションデータと観察データの要約統計量は，それぞれ $S_i$ ， $S_{obs}$ とします．
それらの距離を $\rho(S_i, S_{obs})$ とおきましょう．

試行回数は50万回で，以下のような手順になります．

事前分布から50万個のパラメータ $\theta _i$ をサンプリング
それらの $\theta_i$ をもとに，モデル $f(D|\theta)$ を利用して50万セットのシミュレーションデータ $D_i$ を生成
$D_i$ の要約統計量 $S_i$ を計算
観察データ $D_{obs}$ の要約統計量 $S_{obs}$ に十分近いシミュレーションを採用
$\rho$ が十分小さければ，採用されたシミュレーションのパラメータ群が事後分布に近似できる

abc は，このうち4, 5を計算・評価してくれるパッケージになっています．

まず手順1-3をやってみます．

# Set observed data, iteration number, and others.
data.obs <- c(2.2, 3.5, 2.9, 3.6, 2.8)
iteration <- 500000
sigma.hat <- 0.5
num <- 5
library(abc)

# Set prior distribution for theta.
SetPrior <- function(n){
  rnorm(n, 0, 10)
}

# Simulate data from sampled theta.
SetNormData <- function(theta){
  rnorm(num, theta, sigma.hat)
}

# Model
Model <- function(iteration){
  theta <- SetPrior(iteration)
  data <- sapply(theta, SetNormData)
  ss.sim <- apply(data, 2, mean)
  return(list(theta = theta, ss = ss.sim))
}

# Calculate summary statistics for observed data.
ss.obs <- mean(data.obs)

# Simulate theta and  summary statistics.
sim <- Model(iteration)

この段階で $\theta$ の分布を見ると，当たり前ですが，きちんと事前分布 $ND(0, 10^2)$ にしたがっています．

hist(sim$theta, breaks = 50, xlab = expression(theta),
  main = "All sampled theta")

次に，パッケージabcを利用して手順4-5をやってみます．許容範囲は1%としておきましょう．

# Rejection algorithm (tolerance = 1%).
abc.rej <- abc(target = ss.obs, param = matrix(sim$theta, dimnames = list(NULL, "theta")),
  sumstat = sim$ss, tol = 0.01, method = "rejection")

# Local linear regression algorithm (tolerance = 1%).
abc.lin <- abc(target = ss.obs, matrix(sim$theta, dimnames = list(NULL, "theta")),
  sumstat = sim$ss, tol = 0.01,method = "loclinear", transf = "none")

# Non-linear regression algorithm (tolerance = 1%).
abc.net <- abc(target = ss.obs, matrix(sim$theta, dimnames = list(NULL, "theta")),
  sumstat = sim$ss, 0.01, method = "neuralnet", transf = "none")

abcの適用結果は，print, summary, hist, plotを用いて表すことができます．
hist()でそれぞれ許容されたパラメータの分布を見てみると，以下のような感じになります．
許容範囲が十分小さければ，これらは，事後分布からのサンプルと見ることができます．

range.x <- seq(-10, 10, x <- 0.05)
oldpar <- par(mfrow = c(3, 1))
hist(abc.rej, breaks = range.x)
points(range.x, sapply(range.x, function(y){dnorm(x = y, sd = 10)}) * iteration * 0.01 * x,
  type = "l", lty = "dotted")
hist(abc.lin, breaks = range.x)
points(range.x, sapply(range.x, function(y){dnorm(x = y, sd = 10)}) * iteration * 0.01 * x,
  type = "l", lty = "dotted")
hist(abc.net, breaks = range.x)
points(range.x, sapply(range.x, function(y){dnorm(x = y, sd = 10)}) * iteration * 0.01 * x,
  type = "l", lty = "dotted")

(※ 事前分布の式が間違っていたので訂正しました．誤: dnorm()，正: dnorm(sd = 10)．図は未訂正です．20111211)

上から順に，"rejection"，"loclinear"，"neuralnet"の結果で，点線は事前分布です．

summary()で，以下のように95%信用区間や中央値を見ることができます．

summary(abc.rej)
Call: 
abc(target = ss.obs, param = matrix(sim$theta, dimnames = list(NULL, 
    "theta")), sumstat = sim$ss, tol = 0.01, method = "rejection")
Data:
 abc.out$unadj.values (5000 posterior samples)

              theta
Min.:        2.1308
2.5% Perc.:  2.5357
Median:      2.9951
Mean:        2.9951
Mode:        2.9780
97.5% Perc.: 3.4587
Max.:        3.8066

summary(abc.lin)
Call: 
abc(target = ss.obs, param = matrix(sim$theta, dimnames = list(NULL, 
    "theta")), sumstat = sim$ss, tol = 0.01, method = "loclinear", 
    transf = "none")
Data:
 abc.out$adj.values (5000 posterior samples)
Weights:
 abc.out$weights

                        theta
Min.:                  2.1389
Weighted 2.5 % Perc.:  2.5542
Weighted Median:       2.9940
Weighted Mean:         2.9936
Weighted Mode:         2.9915
Weighted 97.5 % Perc.: 3.4140
Max.:                  3.7821

summary(abc.net)
Call: 
abc(target = ss.obs, param = matrix(sim$theta, dimnames = list(NULL, 
    "theta")), sumstat = sim$ss, tol = 0.01, method = "neuralnet", 
    transf = "none")
Data:
 abc.out$adj.values (5000 posterior samples)
Weights:
 abc.out$weights

                        theta
Min.:                  2.1388
Weighted 2.5 % Perc.:  2.5542
Weighted Median:       2.9941
Weighted Mean:         2.9936
Weighted Mode:         2.9920
Weighted 97.5 % Perc.: 3.4150
Max.:                  3.7763

重みをつけると，結果が多少変わっているのがわかります．
全体的に， $\theta$ の事後分布として，3.0のまわりを最頻値に2.5-3.5の範囲が推定できており， $\hat{\theta} = 3.0$ と整合的な結果です :D

パラメータの重み付け

abcの"loclinear"と"neuralnet"アルゴリズムが何をしているのか，簡単に説明しておきます．

これらは，許容範囲内のパラメータを許容した後， $S_i$ と $S_{obs}$ が完全には一致していないセットに対して，以下のような補正を行ないます．
$\theta_i = m(S_i) + \epsilon _i$
$m$ は回帰式で， $\epsilon _i$ は共通の分散から生成される誤差項です．
さらに，許容されたシミュレーションには， $S_{obs}$ との一致に応じて重みが与えられます．重みは，カーネル $K$ を用いて， $K(\rho(S_i , S_{obs}))$ で与えられます．
”loclinear”は $m$ に線形回帰を，"neuralnet"は非線形回帰を仮定します．"neuralnet"を利用すれば要約統計量セットの次元を減らすことができます．

さらに，引数hcorrにTRUEを指定すると (デフォルト)，以下の式を用いて $S_i$ と $S_{obs}$ の不等分散を補正してくれます．
$\theta_i = m(S_i) + \frac{\sigma(S_{obs})}{\sigma(S_i)} \epsilon _i$

最適な許容範囲の検討

ABCでは，許容範囲の設定が肝になります．大きすぎると事後分布を正確に近似しませんし，小さすぎるとサンプリング効率が低下します．この問題に対し，パッケージabcのcv4abc関数を利用すれば，最適な許容範囲の検討が可能です．

cv4abc関数は，任意の許容範囲について，leave-one-out cross-validationによって，パラメータ推定の正確さとロバストネスを評価します．アルゴリズムは以下のような感じ．

j番目のシミュレーションがランダムに選ばれ，その $S_j$ が $S_{obs}$ のダミーとして使われる
j番目をのぞいたすべてのシミュレーションから，abcでパラメータ推定
これをn回繰り返して，正確に推定できたか評価

評価には，以下の式を用います (summaryで返されるのがこの値)．
$E_{pred} = \frac{\Sigma (\theta_j - \hat{\theta_j})^2}{Var(\hat{\theta_j})}$
$\hat{\theta_j}$ はj番目のシミュレーションの真のパラメータで， $\theta_j$ は推定されたパラメータです．

全シミュレーションの回数分繰り返すのが理想ですが，時間がかかりすぎるので，ランダムにいくつか選んで評価します．回数は，引数nvalで与えます．

rejection samplingで，許容範囲を1%と50%にとって，実際にやってみるとこんな感じになります．

# Rejection sampling (tolerance = 1%).
cv.rej.1 <- cv4abc(param = matrix(sim$theta, dimnames = list(NULL, "theta")),
   sumstat = sim$ss, tols = 0.01, nval = 10, method = "rejection")
summary(cv.rej.1)
Prediction error based on a cross-validation sample of 10

           theta
0.01 0.000571858

# Rejection sampling (tolerance = 50%).
cv.rej.50 <- cv4abc(param = matrix(sim$theta, dimnames = list(NULL, "theta")),
   sumstat = sim$ss, tols = 0.5, nval = 10, method = "rejection")
summary(cv.rej.50)
Prediction error based on a cross-validation sample of 10

        theta
0.5 0.2077764

真の値と推定値の適合をみると以下のようになります (左は許容範囲1%，右は許容範囲50%)．

oldpar <- par(mfrow = c(1, 2))
plot(cv.rej.1)
plot(cv.rej.50)

summaryの結果と図から，許容範囲1%だとよく適合しているが，50%だと推定誤差が大きくなっていることがわかりますね．

考察

ABCの最大の難点のひとつは，計算が非常に非効率なところです．許容範囲を十分に小さくするほど，また要約統計量の次元が増えるほど，より大部分のシミュレーションが棄却されることになり，計算時間が必然的に増大します．
この問題を解決するため，MCMCやSMC (Sequential Monte Carlo) を利用した効率的なサンプリングが考案されていますが [2-5]，パッケージabcではそれらを利用できません．
実際にABCを実装しようとすると，パッケージabcで採用されているrejection samplingベースのアルゴリズムは効率が悪すぎて使いものにならないことが多いように思いますが，ABCのフレームワークで試しにちょっと解析してみたい，というようなときには手軽にできて良いのではないかと思います．

参考

[1] Csillery et al. 2010. abc: an R package for Approximate Bayesian Computation (ABC). arXiv:1106.2793.

[2] Marjoram et al. 2003. Markov chain Monte Carlo without likelihoods. PNAS 100:15324-15328.

[3] Sisson et al. 2007. Sequential Monte Carlo without likelihoods. PNAS 104:1760-1765.

[4] Beaumont et al. 2009. Adaptive approximate Bayesian computation. Biometrika 96:983-990.

[5] Toni et al. 2009. Approximate Bayesian computation scheme for parameter inference and model selection in dynamical systems. J. R. Soc. Interface 6:187-202.