多変量正規分布をギブスサンプリングで

引き続き「パターン認識と機械学習」(PRML) 11章予習中。

Gibbs サンプリング、これはもう試してみるしか。
syou6162 さんが試してはるの( http://d.hatena.ne.jp/syou6162/20090115/1231965900 )をなぞるだけでもいいんだけど、せっかくだから多次元一般化しよう。

r_mul_norm1 <- function(x, mu, Sig) {
    idx <- 1:length(mu);
    for(a in idx) {
        b <- idx[idx!=a];                    # b = [1,D] - a
        s <- Sig[b,a] %*% solve(Sig[b,b]);   # Σ_ab Σ_bb ^ -1

        # (PRML 2.81) μ_a|b = μ_a + Σ_ab Σ_bb ^ -1 (x_b - μ_b) 
        mu_a_b <- mu[a] + s %*% (x[b] - mu[b]);

        # (PRML 2.82) Σ_a|b = Σ_aa - Σ_ab Σ_bb ^ -1 Σ_ba
        sig_a_b <- Sig[a,a] - s %*% Sig[a,b];

        # [Gibbs] x_a 〜 p(x_a|x_{-a}) = N(μ_a|b, Σ_a|b)
        x[a] <- rnorm(1, mu_a_b, sqrt(sig_a_b));
    };
    x;
}

実装のミソは、条件付き分布を記述するのに必要になる Σ_bb など(詳細は PRML 2章の 2.81, 2.82 式を参照)を求めるところ。
idx <- 1:D (D は空間の次元) に対して b <- idx[idx!=a] とすることで a 以外のインデックスになるので、後は Sig[b,b] で Σ_bb が得られる。
R はこういうのが元の数式に沿って簡単に書けるのがいいよね。

【追記】

ちなみにindexに負数を指定するとその位置の要素だけを取り除いたベクトルが得られます．e.g. c(1, 2, 3)[-2] #=> [1] 1 3 #tsukubar

ってことは、Sig[-a, -a] でいいのか〜。し、知らなかった……orz

実行。
マルコフ連鎖モンテカルロ(MCMC)なので、前回の結果を第１引数に渡して呼び出す。

> mu <- c(1,2);
> Sig <- matrix(c(
+     1.0, 0.7,
+     0.7, 1.0
+ ), 2);
> x <- mu;  # 初期値。(0, 0) とかでもいいんだけど、一応
> (x <- r_mul_norm1(x, mu, Sig));
[1] 1.699274 2.146990
> (x <- r_mul_norm1(x, mu, Sig));
[1] -0.3963841 -0.5353421
> (x <- r_mul_norm1(x, mu, Sig));
[1] -0.7136735  1.2407015
> (x <- r_mul_norm1(x, mu, Sig));
[1] -0.1796579  0.6226741

たくさんサンプリングして、グラフに描いてみる。

N <- 1000;
data <- matrix(numeric(N*length(mu)),N);
x <- mu;
for(i in 1:N) {
    x <- r_mul_norm1(x, mu, Sig);
    data[i,] <- x;
};
plot(data);

それっぽい。
ちゃんとするなら系列の始まり(MC が収束したと判断できるところまで)を捨てないといけない。x の初期値に mu を与えとけば多分大丈夫！　と手抜き。

問題は各パラメータの相関性と、独立な系列、つまり十分 i.i.d とみなしてよいかどうかというあたり。
自己相関をみることで、それを確認することができる。

> acf(data, 6, plot=F)

Autocorrelations of series ‘data’, by lag

, , Series 1

 Series 1    Series 2   
  1.000 ( 0)  0.706 ( 0)
  0.498 ( 1)  0.349 (-1)
  0.265 ( 2)  0.175 (-2)
  0.116 ( 3)  0.038 (-3)
  0.009 ( 4) -0.009 (-4)
  0.002 ( 5)  0.014 (-5)
  0.005 ( 6) -0.008 (-6)

, , Series 2

 Series 1    Series 2   
  0.706 ( 0)  1.000 ( 0)
  0.712 ( 1)  0.502 ( 1)
  0.364 ( 2)  0.237 ( 2)
  0.151 ( 3)  0.102 ( 3)
  0.056 ( 4)  0.009 ( 4)
 -0.005 ( 5) -0.001 ( 5)
 -0.010 ( 6) -0.009 ( 6)

Series 1 と 2 の相関性は 0.706 ということで、共分散行列に与えたとおりになっている。
また、この値が小さければ独立性が高いわけだが、どれくらいの範囲なら「十分独立」と判断していいのかよくわからない。
plot=T にすると描かれるグラフには、0.06 のあたりに閾値が設定されており、acf(runif(1000), plot=T) とかすると、確かにだいたいその線を下回っているので、-0.06〜0.06 ならよさそう。
そうして上の結果を見ると、３回まわしたくらいではまだ相関性が残ってしまっているが、４回まわせば i.i.d としてよい、と考えられそうだ。

もちろんこの回数は条件によって変わる。一番極端な場合、対角行列を分散とすれば、１回で独立になる(当たり前だけど)。
PRML にはこの回数のオーダーは書かれているが、見積もり方は書かれていない。問題にあわせて適宜選ぶしかないのかな。

せっかく多次元に対応できるように作ったので、４次元とか食わせてみよう。
が、毎回逆行列を計算しているので、次元を上げてサンプル数を増やすとめちゃめちゃ遅くなる。
前回の結果を渡さないといけないのとかも気に入らなかったので、クロージャを使って書き直してみた。

r_mul_norm_closure <- function(mu, Sig) {
    D <- length(mu);
    idx <- 1:D;
    s <- list();
    sd_a_b <- numeric(D);
    for(a in idx) {
        b <- idx[idx!=a];
        s[[a]] <- Sig[b,a] %*% solve(Sig[b,b]);
        sd_a_b[a] <- sqrt(Sig[a,a] - s[[a]] %*% Sig[a,b]);
    };
    x <- mu;
    env <- environment();
    function(I=1) {
        for(i in 1:I) for(a in idx) {
            b <- idx[idx!=a];
            mu_a_b <- mu[a] + s[[a]] %*% (env$x[b] - mu[b]);
            env$x[a] <- rnorm(1, mu_a_b, sd_a_b[a]);
        };
        result <- env$x; # necessary
    };
}
r_mul_norm <- function(N, r_mul_norm1, I=1) {
    t(sapply(1:N, function(i) r_mul_norm1(I)));
}

mu <- c(1,2,3,4);
Sig <- matrix(c(
     0.4468, 0.5442, 0.0644,-0.1021,
     0.5442, 1.1502, 0.3061, 0.0189,
     0.0644, 0.3061, 0.2547, 0.1691,
    -0.1021, 0.0189, 0.1691, 0.2112
), 4);
data <- r_mul_norm(1000, r_mul_norm_closure(mu, Sig));