見出し画像

ベイズ統計学 | 事後分布のサンプリング Direct method

gota_morishita

ベイズ統計学でサンプリングが難しい事後分布を近似的にサンプリングするための基本となる方法論をいくつか紹介します．今回は，事後分布自体を近似し，その近似分布を利用してサンプリングするdirect methodについて紹介していきます．

パラメタが高次元の場合は，いわゆる次元の呪い (curse of dimensionality) が発生し，評価する回数が指数関数的に増加するので，実戦では使われないです．実戦では，上の弱点を克服するためにMCMC，特にHamiltonian Monte Carlo (HMC)が使われます．

ベイズ統計学の復習

ベイズ統計学では，事前分布 (prior distribution) と尤度関数 (likelihood function) をモデル化しデータを用いて事後分布を計算します．

例えば，二項分布モデルを考えてみます．つまり，尤度関数を以下のように定義します．

$$
p(y|\theta) = {n \choose x}\theta^y(1-\theta)^{n - y}. 
$$

ここで，$\theta$を二項分布のパラメタ，成功確率のパラメタ，実際に観測された成功回数を$y$とします．試行回数$n$は定数として扱います．

次に事前分布を決めなければいけません．二項分布の共役事前分布 (conjugate prior)であるベータ分布族を用いると，事後分布が事前分布と同じベータ分布族に属します．今回は$\theta$の事前分布$p(\theta)$をベータ分布$\text{Be}(2, 2)$とします．以下に見るように事後分布は，ベータ分布となっていることが確認できます．

$$
\begin{align*}
p(\theta | y) &\propto p(y | \theta)p(\theta) \\
  &=  {n \choose x}\theta^y(1-\theta)^{n - y}\theta(1-\theta) \\
  &=  {n \choose x}\theta^{y+1}(1-\theta)^{n - y + 1}
\end{align*}
$$

よって，事後分布は$\text{Be}(y + 2, n - y + 2)$です．事後分布からサンプリングすることで中央値，期待値，信用区間 (credible interval) などを推定することができます．実際にRであれば `rbeta`関数を用いれば簡単にベータ関数のサンプルを得ることができます．

しかし，実践では，事後分布がよく知られた分布になることは珍しいです．そのため興味のある統計量を簡単に推定することができないです．

そのためベイズ統計学は近似的なサンプリング方法を頻繁に利用することになります．

今回はその事後分布の近似的なサンプリング方法の基礎となる手法を紹介していきます．

言葉の定義をまずします．

サンプルしたい分布$p(\theta | y)$のことを目標分布 (target distribution) と呼びます．$p(\theta |  y)$の定数倍されたものを$q(\theta | y)$とします．ベイズ統計では正規化（積分して１にならない）されていない場合もよく扱うので，$q(\theta | y)$を用いてサンプリングできる方が都合が良く，基本的に$q(\theta | y)$で用いてサンプリングします．

Direct Method

これは数値積分の手法が基になっています．まずパラメタ$\theta$の範囲をグリッドに分割します．その分割点を$\{\theta_i\}_{i=1}^N$とします．そのグリッドの分割点の各点で，$q(\theta_i | y)$を評価します．1次元の場合は下記の図のようになります．仮に$\theta$の範囲が有界でなかった場合は，他の推定法から$\theta$がとりうる範囲の目星を付けるといった工夫が必要です．

画像1

この数値結果を用いて$p(\theta | y)$を以下のように近似します．

$$
p(\theta | y) =
\begin{cases}
\frac{q(\theta_j | y)}{\sum_{i=1}^n q(\theta_i | y)} \quad \big(\exists j = 1, \dots, N \quad \theta \in [\theta_{j-1}, \theta_{j})\big)\\
0 \quad (\text{otherwise})
\end{cases}
$$

このタイミングで正規化するのでグリッドを評価するタイミングでは目標分布ではなく，$q(\theta | y)$を用いても問題ないです．

これで$p(\theta | y)$の近似$\hat{p}(\theta | y)$が求まりました．これをもとに累積分布関数 (cumulative distribution function) を近似します．単純に  $\widehat{\Pr}(\theta < a) = \sum_{\{i : \theta_i < \theta \}} \hat{p}(\theta_i | y)$ とすれば良いです．

累積分布関数が近似できれば，サンプリングは簡単です．$F \sim p$のとき，$[0, 1]$の一様分布から乱数$X$を発生し，$F^{-1}(X)$を計算すれば$F^{-1}(X) \sim p$となります．

証明は簡単で，一様分布$X$の累積分布関数$F_X(x) = x$であることを用いると，

$$
\begin{align*}
\Pr(F^{-1}(X) < a) &= \Pr(X < F(a)) \\
  &=  F_X(F(a)) \\
 & = F(a).
\end{align*}
$$

となり，$F^{-1}(X) \sim F$が証明されました．

よって一様分布の乱数を発生させて，$F^{-1}(X)$してサンプリングができます．

Rで実践してみよう

上の例で紹介した$\text{Be}(y + 2, n - y + 2)$のサンプリングをしてみましょう．$n = 50, y = 14$とします．すると目標分布は$\text{Be}(16, 36)$となります．まずはgridを作成します．ベータ分布の範囲は$[0, 1]$なのでその範囲で分割します．

grid <- seq(0, 1, length = 1000)

次に，分割点ごとにベータ分布の密度関数を評価します．

density <- dbeta(grid, 16, 36)

今回は正規化された密度関数で評価したのですが必要ないですが正規化をします

density <- density / sum(density)

累積密度関数を近似します．

cdf <- cumsum(density)

一様分布の乱数$X$を発生させ，$F^{-1}(X)$を計算します．

sample.u <- runif(10000)
sample.index <- sapply(sample.u, FUN = function(x) which(abs(cdf - x) == min(abs(cdf - x)) ))
sample.dist <- grid[sample.index]

解析的な式からプロットした図と比べてみます．

hist(sample.dist, breaks = 100, freq = FALSE, xlab = expression(theta), main = 'Histogram of Direct approximation')
curve(dbeta(x, 16, 36), 0, 1, xlab = expression(theta), ylab = 'beta', add =TRUE)

画像2

よく近似されていることがわかります．

Direct methodの欠点

すぐにわかるように$\theta$の次元が上がるとgridの数は指数関数的に増加するので次元の呪いが発生します．

実践的な解析では１００次元を超える場合もあるので，使い物になりません．

次回はMCMCのaccept rejectのルールの元となるrejection samplingという手法を紹介したいと思います．

この記事が気に入ったらサポートをしてみませんか？