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

本章の概要

乱数を使って最適化問題を解くには2つの方法がある．

  • 確率的探索を用いて関数の最大値を探す（本章前半）
  • シミュレーションを使って最適化する関数を近似する（本章後半，5.4以降）

内容：モンテカルロ最適化で関数の最大値を求める

最適化問題には2種類ある

最適化問題には大きくわけて2種類存在する．

  • 定義域\it{\Theta} で関数h(\theta)の極値を求める
  • 定義域\it{\Theta} で陰的方程式g(\theta) = 0 の解を求める

問題設定はそれぞれ違うのだが，今回は最大化問題だけを取り上げる．つまり，以下の式のように定義域\it{\Theta} で一番大きな値を求めたい．

max_{\theta \in \Theta}=h(\theta)

解き方にも2種類ある

  • 数値的手法（5.2 数値的最適化法）
    • 目標関数の特性に強く依存する
  • 確率的手法（5.3 確率的探索）

さて，ある関数のある定義域内における最大値を求めようと思った場合にじゃあどうするかというと，定義域内の値を片っ端から入れていって大きいものを探すといったアプローチと，別の分布をシミュレーションしてそこから最大値が現れるまでシミュレーションし続けるという方法の2通りがあるということになる．ここで重要になるポイントが定義域で，3章のモンテカルロ積分において積分する範囲を決めていたように，5章モンテカルロ最適化においても定義域内で最大値を探すということに繋がる．ただし，モンテカルロ積分よりもモンテカルロ最適化のほうが，全体の平均を求めるよりもピンポイントで値の大きい極値を求めるという点で，一般的に難しい．

例 5.1 コーシー分布から得られたサンプルの尤度を最大化する

左の図

何かしらのコーシー分布から生成された乱数{x_1, \dots ,x_n}が与えられた時に，それを生成したコーシー分布\mathcal{C}(\theta,1) の\theta を最尤法で推定しようという問題．尤度は以下の式で表される

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

実際に計算するときには以下の対数尤度を使う．

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

今回は実験として，上の尤度と対数尤度の2パターンで計算したときに，どういう違いが現れるかを示す．数式上は違いはないが，計算機上では値の取り扱いにおいて違いが出てくる．

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

xm <- rcauchy(500)
loglik <- function(y){-sum(log(1+(x-y)^2))}
lik <- function(y){prod(1/(1+(x-y)^2))}

mi <- numeric(500)
mi2 <- numeric(500)
for(i in 1:500){
  x <- xm[1:i]
  mi[i] <- optimize(loglik, interval=c(-10,10), maximum=T)$max
  mi2[i] <- optimize(lik, interval=c(-10,10),maximum=T)$max
}

plot(mi, type="l", ylim=c(-5, 10), ylab="theta", xlab="sample size")
par(new=T)
plot(mi2, type="l", ylim=c(-5,10),col="red", ylab="theta", xlab="sample size")

尤度で計算した(赤)と対数尤度で計算した計算(黒)を以下のプロットの左側の図に示してある．これをみると，対数尤度で計算した場合では着実に推定値が収束していくのに対し，尤度で計算した場合ではある地点で値が発散してるのがわかる．これは，Rのoptimizeにおいて，尤度を直接観察すると値が発散することがあるからである．このように，確率を扱う場合には，値が非常に小さくなることに注意しなければならない．

右の図

次に，上記の尤度を，摂動が加えられた尤度に置き換える．これは-[\sin(100y)]^2 という0から1の間の値を取るノイズを加えている．

1
2
3
4
5
6
7
8
9
10
11
12
13

perturbedloglik <- function(y){-sin(y*100)^2 - sum(log(1+(x-y)^2))}
perturbedlik <- function(y){-sin(y*100)^2 + prod(1/(1+(x-y)^2))}
mi <- numeric(500)
mi2 <- numeric(500)
for(i in 1:500){
  x <- xm[1:i]
  mi[i] <- optimize(perturbedloglik, interval=c(-10,10), maximum=T)$max
  mi2[i] <- optimize(perturbedlik, interval=c(-10,10),maximum=T)$max
}

plot(mi, type="l", ylim=c(-5, 10), ylab="theta", xlab="sample size")
par(new=T)
plot(mi2, type="l", ylim=c(-5,10),col="red", ylab="theta", xlab="sample size")

これをみると，尤度に関しては左図で示したように発散してしまっている．対数尤度を見てみると，左図と比べて値がばらつきグラフがギザギザになっている．当然だが，ノイズを加えれば値がばらついてしまう．

(個人的な感想)右図は何が言いたいのかちょっと分からない．ここの解釈は本書では出てこないと思うのだが，尤度関数に何らかの関数が加わる状況といえば，最小二乗法などにおける正則化項を加えて過学習を押さえる場合だろうか？　正則化項は摂動というよりかは，lassoやridgeのように1次/2次の関数になるので，今回の場合みたいにノイズのようにはならないと思うのだが…．

例 5.2 Newton-Raphson法を用いて混合正規分布モデルから得られたサンプルの尤度を最大化する

今回の混合正規モデルでは，

\frac{1}{4}\mathcal{N}(\mu_1,1)+\frac{3}{4}\mathcal{N}(\mu_2,1)

を考える．まずは分布の形を調べるために，\mu_1 = 0 ，\mu_2 = 2.5 における対数尤度の分布を示してみる．そして，初期値(1,1)からスタートして，Newton-Raphson法がどのように推定値に収束していくかを図示する．

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

da <- sample(c(rnorm(10^2), 2.5+rnorm(3*10^2)))
like <- function(mu){ -sum(log((0.25*dnorm(da-mu[1]) + 0.75*dnorm(da-mu[2])))) }

mu1 <- seq(-2, 5, le=250)
mu2 <- seq(-2, 5, le=250)
lli=matrix(0,nco=250,nro=250)
for (i in 1:250){
  for (j in 1:250){
    lli[i,j] <- like(c(mu1[i],mu2[j]))
  }
}
image(mu1, mu2, -lli)
contour(mu1, mu2, -lli,nlevels=100, add=T)

sta <- c(1,1)
mmu <- sta
for(i in 1:(nlm(like, sta)$it)){
  mmu <- rbind(mmu, nlm(like, sta, iterlim=i)$est)
}
lines(mmu, pch=19, lwd=2)

アニメーションにしてみる

例が1つだけでは面白くないので，初期値をバラバラに選んで何回もシミュレーションしてみて，初期値からどのように推定値が更新されていくのかを見てみる．今回は，初期値として\mu_1 ，\mu_2 ともに図に表されている-2から5の範囲内でランダムに値を取るようにして，そこからのNewton-Raphson法による推定値の変化を1ステップごとにアニメーションとして作成した．それを25パターン作成して繋げている．

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

library(animation)

mixtureloglik <- function(Niter){
  for(n in 1:Niter){
    sta <- c(runif(1,2-2,5),runif(1,-2,5))
    mmu <- sta
    da <- sample(c(rnorm(10^2), 2.5+rnorm(3*10^2)))
    for(i in 1:(nlm(like, sta)$it)){
      mmu <- rbind(mmu, nlm(like, sta, iterlim=i)$est)
    }
    for(j in 1:nrow(mmu)){
      image(mu1, mu2, -lli, main=paste("Iter:",n))
      contour(mu1, mu2, -lli,nlevels=100, add=T)
      lines(mmu[1:j,1],mmu[1:j,2], pch=19, lwd=2)
    }
  }
}

saveMovie(mixtureloglik(25), interval=0.05, moviename="",
          movietype="gif", outdir=getwd(),
          width=640, height=480)

Rのanimationというパッケージを使用して作成している．このパッケージでは標準でmpgも出力できるらしいのだが，今回上手くいかなかったので，gifアニメとして出力したのちにImagemagickとffmpegを使用してmpgの動画として出力した．gifアニメは思いの外重いので，今回はYoutubeにアップロードした．

1
2

$ convert animation.gif Rplot%d.png
$ ffmpeg -f image2 -i Rplot%d.png animation.mpg