モンテカルロ法で条件付期待値計算をする際の試行錯誤−２（Rao-Blackwellization法（条件付きモンテカルロ法））

問題設定

以前試行錯誤した内容
モンテカルロ法で条件付期待値計算をする際の試行錯誤−１ - My Life as a Mock Quant
の続き。R言語の書籍「Rによるモンテカルロ法入門」

の4.6「Rao-Blackwellization法と脱条件化」をよく読んでみると私が抱えている問題に対する回答があった（ように思える）ので、メモがてらまとめる。間違っていたらコメント欄等でご指摘ください。

Rao-Blackwellの定理

まず「Rao-Blackwellization法」と名前にある通り数理統計学の１定理である「Rao-Blackwellの定理」を活用した計算技法なので、
Rao-Blackwellの定理について書いておく。

$g(X)$ は $\theta$ に対する不偏推定量であり、 $T(X)$ は十分統計量で
このとき、 $\mathbb{E} \left[\left. g(X) \right| T(X) = t \right]$ は以下を満たす、
(1) $\mathbb{E} \left[\left. g(X) \right| T(X) = t \right]$ は $\theta$ に関数する不偏推定量
(2) $Var(\mathbb{E} \left[\left. g(X) \right| T(X) = t \right]) \leq Var(g(X))$
となる。

というものがあって、これがモンテカルロ計算に応用できるぞというお話。
ここでいう十分統計量と言うのは”その統計量を所与のものとした場合の確率分布はパラメータ $\theta$ に依存しない”という統計量。
数式で書くと十分統計量 $T(X)$ とは
$Pr\{X |T(X) = t, \theta \} = Pr\{X |T(X) = t\}$
ということ。つまり、”パラメーターに関する情報はすべてこの十分統計量に集約されていて、それ以上の情報は必要ない”ということを意味している。

Rao-Blackwellization法

上述のRao-Blackwellの定理は十分統計量で条件付けた上での定理だったが、モンテカルロ計算への応用としては一般の確率変数（の作るσ集合体）で条件付けしたものを使う。実際、Rao-BlackWellの定理に類似した関係式、というか、条件付分散に対して一般に成り立つ関係式があって
$\theta := \mathbb{E}_X \left[X\right]$
$Var(X) = \mathbb{E}_X \left[ \left( X - \theta \right)^2 \right] = \mathbb{E}_Y \left[ \mathbb{E}_X \left[ \left. \left( X - \theta \right)^2 \right| Y \right] \right] \geq \mathbb{E}_Y \left[ \mathbb{E}_X \left[ \left. \left( X - \theta \right) \right| Y \right]^2 \right] = \mathbb{E}_Y \left[ \left( \mathbb{E}_X \left[ \left. X \right| Y \right] - \theta \right)^2 \right] = Var\left( \mathbb{E}_X \left[ \left. X \right| Y \right] \right)$
$\therefore Var\left( \mathbb{E}_X \left[ \left. X \right| Y \right] \right) \leq Var(X)$
となる。
ここで
$\mathbb{E}_X\left[ f(X) \right] = \int_{-\infty}^{\infty} p(x) f(x) dx = \int_{-\infty}^{\infty} \left( \int_{-\infty}^{\infty} p(x,y) dy \right) f(x) dx = \int_{-\infty}^{\infty} \left( \int_{-\infty}^{\infty} p(x|y)p(y) dy \right) f(x) dx = \int_{-\infty}^{\infty} p(y) dy \left( \int_{-\infty}^{\infty} p(x|y)f(x) dx \right)=\mathbb{E}_Y \left[ \mathbb{E}_X \left[ \left. f(X) \right| Y \right] \right]$
$\therefore \mathbb{E}_X\left[ f(X) \right] = \mathbb{E}_Y \left[ \mathbb{E}_X \left[ \left. f(X) \right| Y \right] \right]$
及び
$\mathbb{E}_X \left[ X^2 \right] \geq \mathbb{E}_X \left[ X \right]^2$ （Jensenの不等式より）
という２つの関係式を使用している。また $\mathbb{E}_Z$ と書いた場合、確率変数 $Z$ に対する期待値操作とする。

結局、
$Var\left( \mathbb{E}_X \left[ \left. X \right| Y \right] \right) \leq Var(X)$
というこの不等式が何を言っているのかと言うと、同じモンテカルロパス数で真のパラメーター $\theta$ を
・ $\mathbb{E}_X \left[ X \right]$
・ $\mathbb{E}_Y \left[ \mathbb{E}_X \left[ \left. X \right| Y \right] \right]$
の２通りのやり方でモンテカルロ推定・計算する際には、モンテカルロ誤差(モンテカルロ推定量の誤差)が、各確率変数 $X,\mathbb{E}_X \left[ \left. X \right| Y \right]$ の標準偏差に比例するので、条件付での計算の方がモンテカルロ誤差が小さくなるということを言っており、これが「Rao-Blackwellization法」と呼ばれているモンテカルロ誤差を減少させる計算手法となっている*1。

実際にやってみる

ここでは練習問題4.15の(C)を試しに解いてみる。設定としては
・ $X|y \sim \mathcal{Bin}(n,y)$
・ $Y \sim \mathcal{Be}(\alpha, \beta)$
と、確率変数 $X$ の $y$ での条件付分布が二項分布、 $Y$ の分布がベータ分布に従うものとする。
この時、確率変数 $X$ の（条件なしの）分布は周辺化することで
$p(X) = \int_0^1 P(x|n,y) P(y|\alpha,\beta) dy = \int_0^1 _nC_x y^x (1-y)^{n-x} \frac{y^{\alpha-1}(1-y)^{\beta-1}}{B(\alpha,\beta)} dy = \frac{_nC_x}{B(\alpha,\beta)} \int_0^1 y^{x + \alpha-1}(1-y)^{n-x+\beta-1} dy = \frac{_nC_xB(x+\alpha,n-x+\beta)}{B(\alpha,\beta)}$
と求められる。ただし $B(\alpha,\beta)$ はベータ関数 $B(\alpha,\beta)=\int_0^1 y^{\alpha-1}(1-y)^{\beta-1}dy$
この時の $X$ の分布はベータ二項分布という分布になっている。この分布に従う乱数生成関数はR言語にデフォルトで搭載されてはいないが、二項分布とベータ分布の混合分布となっているので、混合分布からの乱数生成アルゴリズムをそのまま適用できる。つまり

$y$ をベータ分布からサンプルとして引く
サンプリングした $y$ の値を二項分布の確率パラメータ(試行の成功確率)として用い、その二項分布からサンプリングする

とすればベータ二項分布からのサンプルを生成できる。

以下、R言語での実験結果。それぞれ
$\mathbb{E}_Y \left[ \mathbb{E}_X \left[ \left. X \right| Y \right] \right] \sim \frac{1}{N} \sum_{i=1}^N nY_i, \left( Y_i \sim \mathcal{Beta}(\alpha,\beta) \right)$
$\mathbb{E}_X \left[ X \right] \sim \frac{1}{N} \sum_{i=i}^N X_i, \left( X_i \sim \mathcal{BetaBin}(n,\alpha,\beta) \right)$
となる。ここで $X|y \sim \mathcal{Bin}(n,y)$ の場合、 $\mathbb{E}_X \left[ \left. X \right| Y \right] = nY$ と計算可能であることを使っている。

#各種パラメーターの設定
SIZE.SAMPLE <- 10000
SIZE.TRIAL <- 10
ALPHA <- 0.5
BETA  <- 0.5
#ベータ二項分布に従う乱数の生成
rbetabinom <- function(n, size, alpha, beta)
{
  rbinom(n, size = size, prob = rbeta(n, alpha, beta))
}
#ふつーにモンテカルロ計算する法
x <- rbetabinom(SIZE.SAMPLE, SIZE.TRIAL, ALPHA, BETA)
#Rao-Blackwellizationで分散を軽減させた法
y <- rbeta(SIZE.SAMPLE, ALPHA, BETA) * SIZE.TRIAL

結果を見ると

> c(mean(x), sd(x))
[1] 4.975400 3.708931
> c(mean(y), sd(y))
[1] 5.044833 3.534325

ちょっとだけだけど、分散が減少しているのがわかる*2。

参考

Rao–Blackwell theorem - Wikipedia
http://www.imes.boj.or.jp/research/papers/japanese/kk26-b2-3.pdf:title:"CDOプライシングの離散高速アプローチ(pdf)"（条件付きモンテカルロ法についての解説が少しある）

の５章６節「条件付きモンテカルロ法」に”分散減少法”の１つとして紹介があった。

*1:Rao-Blackwellの定理自体を直接使ってるものではないと思ってる。名前・人物にあやかっているだけということかな？？？

*2:正直、「え、こんなもんなの・・・」レベルの印象