「Rによるモンテカルロ法入門」読書ノート：5章モンテカルロ最適化確率的探索

読書ノートアーカイブ：http://yagays.github.com/blog/2012/10/20/archive-introducing-monte-carlo-methods-with-r/

例 5.3 確率的探索を使って，一様分布からのシミュレーションにより関数の最大値を求める

例 3.3で積分を考えた関数 $h(x) = [\cos(50x) + \sin(20x) ]^2$ について，今度はその最大値について考える．

今回の実験は1000回シミュレーションを行う．それぞれのシミュレーションには，まず一様乱数 $\mathcal{U}_{[0,1]} = u_1, \ldots , u_m$ を1000個作成し，それぞれに関数を掛けた値の中から最大値 $h_{m}^{*} = \max(h(u_1), \ldots ,h(u_m))$ を求める．

確率的探索といっても，今回の場合にやっていることは結局のところ定義域内で乱数を発生させて関数の値が大きいところを探しているだけとなる．今回は一様乱数で定義域内をまんべんなく探すという方法を取っているが，この後の確率的探索の例では「密度」を考えて，関数の値が大きい領域でのシミュレーションの確率を高くし，その他の関数の値の小さい領域での確率を小さくするという方法が取られる．

今回のシミュレーションの仕組み

今回のコードは少しややこしいので，まずは試しにシミュレーションを1回だけ行なってみよう．今回は推定値がどのように更新されていくかを見るために，単純にmax()関数を使うのではなく，cummax()関数を使用している．これにより，ベクトルの前から順に値を見ていって，最大値を更新していく様子を見ることができる．以下の図の黒の線が推定値 $h_{m}^{*}$ ，赤の線がoptimize()関数で求めた定義域内における関数 $h(x)$ の最大値を示している．

h <- function(x){(cos(50*x) + sin(20*x))^2}
plot(cummax(h(runif(10^3))),type="l")
abline(h=optimize(h,int=c(0,1), maximum=T)$ob,col="red")

実際に例5.3のシミュレーションを行う

では実際に1000回のシミュレーションを行なって，推定値がどのように更新されていくのか，そしてその値がどれほどバラけるのかを見ていく．

以下の図が今回の結果となっている．まず，本書に載っている左図ではなく，右の赤い線が重なった図から説明する．この図は1000回シミュレーションした中から10回だけ抜き出してプロットしたものとなっている．これを見ると，最初のうちは0に近いような値から始まっていたりと値が安定しないものの，200回もすると実際の値と殆ど変わらないような推定値になり落ち着いているのがわかる．このように，値がきちんと実際の値に近づくかどうかを見るために，各イテレーションごとの最小値を求めてプロットしたのが左図となっている．灰色の部分は，イテレーションごとの推定値の幅を示しており，これを見ると，800回程度実験すれば一応どのシミュレーションでも実際の値と推定値がほぼ同じになることが分かる．

h <- function(x){(cos(50*x) + sin(20*x))^2}
rangom <- h(matrix(runif(10^6), ncol=10^3))
monitor <- t(apply(rangom, 1, cummax))

par(mfrow=c(1,2))
plot(monitor[1,], type="l", col="white", ylim=c(0, optimize(h,int=c(0,1), maximum=T)$ob),
     xlab="Iterations", ylab="h(theta)") # 描写する領域をplotしているだけ
par(new=T)
polygon(c(1:10^3, 10^3:1), c(apply(monitor,2,max), rev(apply(monitor, 2, min))), col="gray")
abline(h=optimize(h, int=c(0,1), maximum=T)$ob)

plot(monitor[1,],type="l",col="red", ylim=c(0,optimize(h, int=c(0,1), maximum=T)$ob),
     xlab="Iterations", ylab="h(theta)")
for(i in 2:10){
  lines(monitor[i,],col="red")
}

練習問題 5.3 一様分布を用いたシミュレーションのパフォーマンスを向上させるために使用する定義域における制約を可視化する

まず最初に注意．この問題は解けないわけではないが，厳密には間違っているらしい．訂正は解答(pdf)に書かれているが，つまりは定義域 $\it{\Theta}$ の制約は，

$x^2(1+\sin(3y)\cos(8x)) + y^2(2+\cos(5x)\cos(8y)) \leq 1$

ではなく，

$x^2(1+\sin(\frac{3}{y})\cos(8x)) + y^2(2+\cos(5x)\cos(8y)) \leq 1$

としなければいけないようだ．今回はどちらの式も試してみて，本に書かれている式の何がおかしかったのを見てみる．

まず，左図は解答に記載されている修正後の式を使った場合，右図は本書に載っているオリジナルの式を使った場合となっている．これを見ると，左図では黒の円の中に式で表された制約の定義域が収まっているのに対し，右図では黒の円から一部制約の定義域がはみ出ているように見える．

theta=runif(10^5)*2*pi
rho <- runif(10^5)
xunif <- rho*cos(theta)/.77
yunif <- rho*sin(theta)

par(mfrow=c(1,2))
plot(xunif, yunif, pch=19, cex=.4, xlab="x", ylab="y")
const <- (xunif^2*(1+sin(yunif/3)*cos(xunif*8))+
  yunif^2*(2+cos(5*xunif)*cos(8*yunif))<1)
points(xunif[const], yunif[const], col="cornsilk2", pch=19,cex=.4)

plot(xunif, yunif, pch=19, cex=.4, xlab="x", ylab="y")
const2 <- (xunif^2*(1+sin(yunif*3)*cos(xunif*8))+
  yunif^2*(2+cos(5*xunif)*cos(8*yunif))<1)
points(xunif[const2], yunif[const2], col="red", pch=19, cex=.4)

また，この方法のパフォーマンスを測定するために，棄却された点の平均数で評価してみる．生成した乱数の総数が10⁵ 個なので，74%は受理され，26%が棄却されたことがわかる．

summary(const)
# Mode   FALSE    TRUE    NA's 
# logical   26285   73715       0

それにしても，この図に何とも言えない可愛らしさを感じる，宇宙の侵略者を撃ち落とす某ゲームで見たことがあるような…．