読者です 読者をやめる 読者になる 読者になる

廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

R によるシミュレーテッド・アニーリング法の練習

アルゴリズム

関数 h(\theta) を最大にする \theta を探す問題を最適化問題という.

シミュレーテッド・アニーリング法は以下のようなアルゴリズム最適化問題を解く.

反復 t で,

  1. \xi \sim g(\xi) をシミュレートする.
  2. 確率 p_t = \min\{\exp(\Delta h_t/T_t), 1\}\theta _{t+1}=\theta_t + \xi を受理する. さもなくば  \theta_{t+1} = \theta_t とする.

ここで  \Delta h_t = h(\theta_t + \xi) - h(\theta_t). また密度 g は原点 0 を中心に対称であること以外は任意とする.

温度と呼ばれるパラメータ T_t は適当に決めてやる必要がある.

例題

人工的な関数だとなんとなく気分がのらないので, ゼロ切断ポアソン分布のパラメータの最尤推定ゼロ切断ポアソン分布のパラメータの最尤推定 - 廿TT)をやることにする.

密度 g を一様分布とする.

g のスケールは大きくしすぎると,( \Delta h_t が負になるので)推移のほとんどが棄却されてしまい効率が悪い.

一方で小さすぎると, 遷移が小さくなり, 初期値に依存しやすくなる.

受理率が 0.2 から 0.6 の間くらいになるようにするとよいとのこと.

温度の系列は, 今回 1/(1+t)^2 とした. こうするとなんかいいらしい.

温度とスケールはうまいぐあいにチューニングする必要があるが, どうすればうまくいくかはよくわからない.

シミュレーテッド・アニーリング法を実行する R のコードは以下のようになった.

#尤度関数の定義
lambda <- 1
logL0 <- function(ysum,n){
  function(loglambda) {
    lambda = exp(loglambda)
    log(lambda)*ysum - n*lambda - n*log(1-exp(-lambda))
    }
}
y <- rpois(300,lambda) 
y <- y[y!=0]
h <- logL0(sum(y),length(y))

#シミュレーテッド・アニーリング法を実行する関数
mySA <- function(theta0,h,Nsim,Temp,scale=1){ 
  hcur = h(theta0)
  theta = numeric(Nsim)
  theta[1] = theta0
  xis = runif(Nsim,-1,1)*scale
  accept = 0
  for(t in 2:Nsim){
    prop = theta[t-1] + xis[t]
    hprop = h(prop)
    if(Temp[t]*log(runif(1)) < hprop-hcur){
      theta[t] = prop
      hcur = hprop
      accept = accept + 1
    }else{
      theta[t] = theta[t-1]
    }
  }
  list(theta=theta[Nsim],theta_all=theta,accept_rate=accept/Nsim)
}
Nsim <- 10000
Temp1 = 1/((1+1:Nsim)^2)
res1 <-mySA(-20,h,Nsim,Temp1,scale = 0.01)

ポアソン分布のパラメータ λ は負の値を取らないが, そういった制約をつけるのはむずかしいので, 指数関数の中に入れることで, 負の値を取らないようにしている.

Temp[t]*log(runif(1)) < hprop-hcur の部分は, u を一様乱数として, \exp(\Delta h_t/T_t) < u と同じ意味.

初期値はわざとちょっと遠くにおいている.

受理率は 0.4 くらいになった.

> res1$accept_rate
[1] 0.402

\theta_t の系列をプロットしてみる.

plot(res1$theta_all,type="l")

f:id:abrahamcow:20160918033603p:plain

わりと一直線に解に向かっている.(もっとギザギザしてくれたほうが確率的最適化っぽい気がする. )

得られた最適解を尤度関数のグラフに重ねてみる.

curve(h(x),-2,2,n=1000)
abline(v=res1$theta,lty=2)
abline(h=h(res1$theta),lty=2)

f:id:abrahamcow:20160918034004p:plain

optimise 関数(黄金分割法:一次元の最適化、黄金比による探索 - 廿TT)で求めた結果と比較する.

> res1$theta
[1] 0.008970653
> h(res1$theta)
[1] -107.7183
> optimise(f=h,interval = c(-2,2),maximum = TRUE)
$maximum
[1] 0.008967951

$objective
[1] -107.7183

ほぼ同じ値が求まっている.

参考文献

Rによるモンテカルロ法入門

Rによるモンテカルロ法入門