R言語でデータ拡大(Data Augmentation)サンプリング法を実装してみた

データ拡大 (Data Augmentation)サンプリング法ってのをお勉強したのでそのまとめ。

概要＆アルゴリズム

説明は２つ目の参考LINKに書いてある方法でずばり良いと思う。
要するに

分布 $f(x)$ から直にサンプリングしたいんだけど、それが難しいという状況の時、適当なダミー変数 $y$ 的なものをかませて
$\int f(x,y) dy = f(x)$ となるような $f(x,y)$ を（適当に・何とかして）構築。
この $f(x,y)$ に従う乱数をギブスサンプラー

$f(x|y)$ から $x$ を取得

$f(y|x)$ から $y$ を取得

を何度も繰り返し、系列データ $\{(x_1,y_1),(x_2,y_2),\dots\}$ を生成する事が出来たならば、
ここで得られた $x$ の系列 $x=\{x_1,x_2,\dots\}$ はまさに望む $f(x)$ からの分布となっている

というもの。
元の $x$ が作る空間に加えて $y$ が作る空間も考えているわけなのでデータ"拡大"サンプリングということですかね。
元の $f(x)$ からサンプルできない状況でも条件付の $f(x|y),f(y|x)$ からサンプリングできればいいってのがアイディアの部分なのだろうか。
１つ前のエントリ
R言語でスライスサンプリング(Slice Sampling)を実装してみた - My Life as a Mock Quant
で書いたスライスサンプリングよりも、条件付分布をサンプリング対象の分布ごとに手計算してやらなければならないという意味で汎用性が落ちるように思える。

ちなみに↓の本では10章をデータ拡大 (Data Augmentation)の話に充てている。

R言語での実装＆動かす

とまぁ、つらつら書いてきたけど、自分の中では「この手法、要するにギブスサンプラーだろ」と思ってしまっているので、あまり新規性を感じない。
というわけで、 $f(x)= \frac{1}{\sqrt{2\pi}} \e^{-\frac{x^2}{2}}$ と一番簡単に生成・結果チェック出来るであろう正規分布としてやってみる*1。

この時、同時分布として
$f(x,y) = \frac{1}{\sqrt{2}\pi} \e^{-(x^2-\sqrt{2}xy+y^2)}$
を考えると、それぞれの条件付分布は
$f(x|y) = \frac{1}{\sqrt{\pi}} \e^{-(x-\frac{y}{\sqrt{2}})^2},Y|X \sim N( \frac{x}{\sqrt{2}},\frac{1}{2})$
$f(y|x) = \frac{1}{\sqrt{\pi}} \e^{-(y-\frac{x}{\sqrt{2}})^2},X|Y \sim N( \frac{y}{\sqrt{2}},\frac{1}{2})$
となり、 $x$ の周辺分布は
$f(x) = \int f(x,y) dy = \frac{1}{\sqrt{2\pi}} \e^{-\frac{x^2}{2}}$
となるので、データ拡大アルゴリズムを適用してサンプリングを行えば $f(x)$ に従うような $x$ が手に入る。

#f(x|y)とf(y|x)。この例だとたまたま同じ関数でいける
rand.conditional <- function(val){1/sqrt(2)*rnorm(1)+val/sqrt(2)}
#Data Augmentationアルゴリズムで生成（要するにギブスサンプラー）
SIZE <- 10^5
points <- numeric(SIZE)
for(i in 1:SIZE)
{
  y <- rand.conditional(points[i])
  points[i+1] <- rand.conditional(y)
}
#yは完全に無視してx(points)のみの分布をチェック
hist(points, SIZE^0.5, freq=FALSE)
x<-seq(-3 ,3 ,0.01)
lines(x, dnorm(x), col=2, lwd=3)