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

例 5.4 数値的最適化法と確率的探索を比較する

例5.1の続き．前回は数値的最適化法によって\theta の値を推定したが，今回は確率的探索による推定を行なって，両者の手法の比較を行う．問題としては，以下の対数尤度を最大化する．

\log l(\theta| x_1,\ldots,x_n) = \sum_{i=1}^{n} \log \frac{1}{1+(x_i-\theta)^2}

この式はn \rightarrow \infty において\theta^* =0 になる．すなわちこれから求める推定値の答えは0ということになる．

ソースコード

以下のソースコードにおける変数の対応は以下のようになっている．

  • 数値的最適化法
    • trul：推定値
    • truv：対数尤度
  • 確率的探索
    • loc：推定値
    • maxv：対数尤度

なお，mcsmパッケージに付随するデモスクリプトでは，数値的最適化法における乱数を正規乱数としていたり対数尤度の関数が微妙に違っているが，今回は例5.1の解法と同様のコーシー乱数と対数尤度の式を用いている．

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

ref <- rcauchy(5001)
f <- function(y){-sum(log(1+(x-y)^2))}

maxv <- NULL
loc <- NULL
truv <- NULL
trul <- NULL

for (i in 10:length(ref)){
  # 数値的最適化
  x <- ref[1:i]
  tru <- optimize(f,c(-5,5),maximum=T)
  trul <- c(trul, tru$max)
  truv <- c(truv, tru$ob)

  # 確率的探索
  prop <- runif(10^3,-5,5)
  vale <- apply(as.matrix(prop),1,f)
  loc <- c(loc, prop[order(-vale)[1]])
  maxv <- c(maxv, max(vale))
}

par(mar=c(4,4,1,1), mfrow=c(2,1))
plot(trul, loc, cex=.5, pch=19, xlab=expression(theta^0), ylab=expression(hat(theta)))
abline(a=0, b=1, col="grey")
plot(10:length(ref), (truv-maxv)/abs(truv), type="l", lwd=2, xlab="Sample size", ylab="Relative error", ylim=c(0,0.02))

まず，上のプロットは，数値的最適化法による推定値を横軸に，確率的探索による推定値を縦軸に取った図となっている．真の値が\theta^* =0 なので，当然ながら両者ともに推定値は(0,0)付近に集まっている．また，推定値が外れる場合では，どちらかの値が大きくなればもう片方も大きくなるといった具合に，そのパターンはグレーの線で示した\theta^0 = \hat{\theta} の線分の上に乗るような形となる．

次に，下のプロットは，横軸にサンプルサイズを取り，縦軸に推定値における対数尤度の差を取った図となっている．今回は相対的誤差なので，(数値的最適化法-確率的探索)/(数値的最適化法)といった形で値を求めている．相対的誤差が大きいということは，数値的最適化法と確率的探索それぞれ求めた対数尤度に差がある，すなわち\theta に差があるということで，上のプロットにおけるグレーの線分から離れれば離れるほど相対的誤差も大きくなるということになる．なお，図示しているサンプルサイズはn=10 からn=5001 までの値となっており，これは両者の手法がどちらも一定サイズのサンプルサイズを必要とするからである（つまり1から10までは図示していない）．

例 5.6 グローバルミニマムのまわりにローカルミニマムが配置されている関数を可視化する

以下の関数の最小化を検討する．

h(x,y)=(x\sin(20y)+y\sin(20x))^2\cosh(\sin(10x)x)+ (x\cos(10y)-y\sin(10x))^2\cosh(\cos(20y)y)

今回は可視化するだけだが，この後の例5.7などで度々登場する関数となる．3Dプロットを描写すると，(0,0)に最小値があり，その周りには山あり谷ありとローカルミニマムがたくさん配置されている関数だとわかる．この図を上から見ると，P.157のプロットのようになる．

1
2
3
4
5
6
7
8
9
10

h <- function(x,y){
  (x*sin(20*y) + y*sin(20*x))^2 * cosh(sin(10*x)*x) +
    (x*cos(10*y) - y*sin(10*x))^2 * cosh(cos(20*y)*y)
}
x <- seq(-3, 3, le=435)
y <- seq(-3, 3, le=435)
z <- outer(x, y, h)
par(bg="wheat", mar=c(1,1,1,1))
persp(x, y, z, theta=155, phi=30, col="green4",
      ltheta=-120, shade=0.75, border=NA, box=FALSE)