ランダムウォークで乱数を生成する - 驚異のアニヲタ社会復帰の予備

MikuHatsune2018-04-25

機械学習のTA なのに機械学習素人なので機械学習と統計の講義を聞いている。
多次元な確率分布から乱数をいい感じに取ってくる方法に、ランダムウォークを使う。
多次元な確率分布としては、単純に2次元正規分布、とする。二次元の正規分布は、x1 とx2 のそれぞれの変数が平均 $\mu_1$ , $\mu_2$ で、変数間の関係として分散共分散行列 $\Sigma=\begin{pmatrix}\sigma_1^2 & \sigma_1\sigma_2 \\ \sigma_2\sigma_1 & \sigma_2^2\end{pmatrix}$ である。
R ならrvtnorm パッケージが使える。
　
アルゴリズムとしては、二次元正規分布の確率密度関数を $g(x)$ として、ある適当な初期値 $x_0$ から開始する。
ある $x$ のとき、どれくらいランダムウォークするかは、一様分布 $U$ からランダムに決める。ここでは歩く量 $w\sim U(-0.5, 0.5)$ とする。つまり、x方向もy方向も -0.5 から0.5 までの一様分布からサンプリングされる。
さて、ここでランダムサンプリングするかしないかは、以下のルールで決める。
ある $x_{t}$ のとき一様分布からランダムウォークする量が決まって、 $v_{t}$ に移動したとすると、次の位置 $x_{t+1}$ は、 $x_{t+1}\leftarrow v_{t}$ として位置を変えるか、 $x_{t+1}\leftarrow x_{t}$ として動かないか、である。
この動く/動かない、を決めるのは、確率密度比 $\frac{g(v_t)}{g(x_t)}$ より、 $a_t=\textrm{min}\{1, \frac{g(v)}{g(x)}\}$ として確率を決めて、
$P(x_{t+1}=v_t)=a_t$ ：動く
$P(x_{t+1}=x_t)=1-a_t$ ：動かない
とする。
確率密度比 $\frac{g(v_t)}{g(x_t)}$ が1 より大きいということは、 $g(v_t)$ のほうが $g(x_t)$ よりもっともらしい、ということで、動くべき、ということになる。逆のときは動かないほうがいい、ということになるが、それでもちょっとは動くことも可能性としては残したほうがよい、という気分がある。というのが、（確率密度の）比は0 以下にはならないからである。
　
という感じで、適当な初期値、確率密度関数(ここではdrvnorm で得られる)、繰り返し回数、を指定すると、乱数ができる。繰り返しは10000回行った。
　
できた。

　
さて、これは初期値依存なので、初期値がうまく設定されないと、いい感じにならない。ここで、初期値を少し離して実行してみる。
　
できた。

　
ここで試している確率密度は単純なもので、初期値が目的の関数から離れていても、サンプリング回数が多くなればそこそこたどり着く、ことが多い。ここで、ある程度離れた初期値で、サンプリング回数も減らしてみる（1500回）。
　
最初の100サンプリングくらいでもう目的の関数部分にたどり着いているので、それなりにいい感じの乱数ができているが、上の図とはサンプル数が段違いなので、カーネル密度推定したものの楕円形（これは二次元正規分布で設定した $\Sigma$ の傾き具合に相当する）がグニャグニャである。

　
上の図では、ランダムウォークの過程を各 $t$ のときで色付けしておいた。ここで、確率密度比 $\frac{g(v_t)}{g(x_t)}$ （の対数）と、確率密度比の確率を用いて点が移動したのか否かの記録をプロットしている。
確率密度比 $\frac{g(v_t)}{g(x_t)}$ （の対数）が非常に大きいもしくは小さいときは、最初の50点くらいで、それ以降はほとんど1、つまり $\frac{g(v_t)}{g(x_t)}$ がほぼ等しい状態になっている。この時点では、ランダムウォークの点は二次元正規分布に近いところをウロウロしている。ランダムウォークの移動量は、 $U(-0.5, 0.5)$ のけっこう狭いグリッド内なので、一発のランダムウォークで二次元正規分布から逸脱することはない。
$\frac{g(v_t)}{g(x_t)}$ が1 以上のときは、確率1 で移動することになるので、移動したことを示すマゼンタが記録されている。 $\frac{g(v_t)}{g(x_t)}$ が1 未満のときは、確率 $\frac{g(v_t)}{g(x_t)}$ で移動して確率 $1-\frac{g(v_t)}{g(x_t)}$ でその場にとどまるので、 $\frac{g(v_t)}{g(x_t)}=1$ （対数では0）より下の領域では、2色出ている。

　
こういうのがRstan でのいわゆるwarm up や収束になる。
　

library(mvtnorm)
library(MASS)
x0 <- c(0.1, 0.2)
x0 <- c(-1, 2)

iter <- 1500 # 繰り返し回数
x <- matrix(0, iter, 2) # ランダムウォークの座標
dratio <- rep(0, iter) # 確率密度比
out <- rep(0, iter) # 確率密度比にしたがって確率的に動いたか、動かなかったかの記録
x[1,] <- x0

rho <- 0.7
sig <- matrix(c(1, rho, rho, 1), 2)
mu <- c(4, 5)

for(i in 2:iter){
  walk <- runif(2, -0.5, 0.5) # 歩く量
  v <- x[i-1,] + walk # 次の候補点
  dratio[i-1] <- dmvnorm(v, mean=mu, sigma=sig)/dmvnorm(x[i-1,], mean=mu, sigma=sig) # 確率密度比
  out[i-1] <- rbinom(1, 1, min(1, dratio[i-1])) # 動いたか、動かなかったか
  x[i,] <- x[i-1,] + out[i-1]*walk
}

kd <- kde2d(x[,1], x[,2], c(bandwidth.nrd(x[,1]), bandwidth.nrd(x[,2])), n=1000)

cols <- jet.colors(iter)
plot(x, type="p", pch=16, cex=0.6, col=2, xlab="x1", ylab="x2")
#lines(x, col=cols)
for(i in 2:iter){
  segments(x[i-1,1], x[i-1,2], x[i,1], x[i,2], col=cols[i])
}
points(x[1,1], x[1,2], pch="★", col=6)
points(mu[1], mu[2], pch="★", col=5)
contour(kd, add=TRUE, col=6)

pcols <- c("black", "magenta")
plot(dratio, log="y", type="n", xlab="sampling number", ylab="log likelihood")
for(i in 2:iter){
  segments(i-1, dratio[i-1], i, dratio[i], col=cols[i])
}
abline(h=1, lty=3)
points(seq(iter), dratio, cex=0.4, col=pcols[out+1], pch=16)
legend("topright", legend=c("動かなかったとき", "動いたとき"), pch=16, col=pcols, cex=2)