Rediscover the Monte Carlo

僕個人はゲームの思考ルーチンを作ることなどには興味があるので、みんな知っていることだと思っていたのですが、意外と「現在世界最強の囲碁の思考ルーチンはモンテカルロ」ってのは知られてないみたいですね。うっかりすると「そんなわけないだろー」とか言われてしまう。その根底には「モンテカルロはとても収束が遅くて使いものにならない」という過去の記憶があるのかなー。ちょうどJavaScriptが使いものにならないおもちゃ言語だと思われていたように。

囲碁の思考ルーチンを著しく進化させた新しいモンテカルロが昔の単純なモンテカルロとどう違うかというと、UCB1という評価関数で「もっと探索するとヨサゲな局面」を判断して、ヨサゲな局面から優先的に探索するという点なんだけど、そういう定性的な話をしてもピンと来ないよね。同じ発想をモンテカルロで円周率を求めるプログラムに適用したら収束の速さが定量的にはっきり見えて面白いんじゃないかなーと思って、空間を分割して調査する価値の高いところから優先的に調査するプログラムを作ってみました。

結果

理屈やソースを先に説明すると最後まで読まれないんじゃないかと危惧したのでまずは結果から(笑)

単純なモンテカルロでは200万回の試行をしてもまだ3.14までしかわかっていないのに対し、改良版モンテカルロでは30万回で3.141まで到達して、20万回の試行の後には3.141593と、かなり近いところまで円周率を求められている。圧倒的じゃないか、我が軍は！

単純モンテカルロ
0: pi = 3.133040
1: pi = 3.141780
2: pi = 3.142373
3: pi = 3.141050
4: pi = 3.142184
5: pi = 3.142260
6: pi = 3.142771
7: pi = 3.143790
8: pi = 3.143636
9: pi = 3.142952
10: pi = 3.143011
11: pi = 3.143880
12: pi = 3.144342
13: pi = 3.144023
14: pi = 3.143701
15: pi = 3.143102
16: pi = 3.143104
17: pi = 3.142973
18: pi = 3.143078
19: pi = 3.143240

改良モンテカルロ
0: pi = 3.142376
1: pi = 3.142209
2: pi = 3.141998
3: pi = 3.141850
4: pi = 3.141732
5: pi = 3.141667
6: pi = 3.141708
7: pi = 3.141727
8: pi = 3.141715
9: pi = 3.141703
10: pi = 3.141722
11: pi = 3.141702
12: pi = 3.141652
13: pi = 3.141622
14: pi = 3.141634
15: pi = 3.141627
16: pi = 3.141620
17: pi = 3.141604
18: pi = 3.141598
19: pi = 3.141593

実行時間はさすがに4倍ほど遅くなった。しかし試行回数が4分の1になったとしてもまだ改良版モンテカルロの方が収束が速い。

単純モンテカルロ
real	0m0.106s
user	0m0.038s
sys	0m0.003s

改良版モンテカルロ
real	0m0.422s
user	0m0.406s
sys	0m0.005s

理屈

理屈の解説。要するにモンテカルロで探索する空間をWIDTH * WIDTHのメッシュに分割して、それぞれの升目をscoreが高い順に試行するってだけ。scoreは「その升目の中で1回試行をした場合にどれくらい結果が変動するかの平均値」にしてある。問題はそのscoreの計算方法なのだが…

まず、各々の試行ごとに「円の中にある/ない」の0/1の値が返ってくるわけなので、ベルヌーイ過程とみなせる。a + b回の試行をして1がa回、0がb回でたとすると、次の試行で1が出る確率はいくらか？ a / (a + b)じゃあないことは明らかだよね。これだと1回試行して1がでたら「次にも確実に1が出る」という誤った結論を出してしまう。

1が出る確率、つまり元のベルヌーイ分布の平均値μ、の平均値がいくらであるのかという話になる。ここでμの共役事前分布はベータ分布で、a + b回の試行をして1がa回、0がb回でた後での事後分布はBeta(a + 1, b + 1)になるとか教科書に書いてあったので(ぇ)、1が出る確率pは(a + 1) / (a + b + 2)と結論できる。

「このマスに円が占める割合」＝「このマスで試行したときに1が出る確率の最尤推定量」＝ a / (a + b)が次の試行でどうupdateされるか考えると、pの確率で(a + 1) / (a + b + 1)、(1 - p)の確率でa / (a + b + 1)に変わるわけだから普通にその2通りでのupdate量に確率を掛けて足しあわせてやればupdate量の平均値が求まる。これが下記コードのscore関数。

(ここ書いていて、今はa / (a + b)の変化量で計算しているけどもしかして(a + 1) / (a + b + 2)の方が適切なのかな？という気がした。マスに円が占める割合の最尤推定量がa / (a + b)なので今の方針であっているとは思うが一瞬不安に思ったので書き留めておく)

(誤解がないように追記しておくと、この評価関数はUCB1ではない。囲碁では「勝てそうな手」を探すのに対しこの円周率推定では「白黒分かれる領域」を探すのでUCB1は使えなかった。)

コード

おかしなところがあれば容赦なく突っ込んでください。大筋では上で述べたように「探索空間をWIDTH * WIDTHに分割して、scoreの高いところから探索する」をやっているだけです。

#include <stdio.h>
#include <stdint.h>
#include <vector>
#include <queue>
using namespace std;

// http://www.jstatsoft.org/v08/i14/xorshift.pdf
uint32_t xor128() {     // 周期は 2^128-1
    static uint32_t x=123456789,y=362436069,z=521288629,w=88675123;
    uint32_t t;
    t=(x^(x<<11));x=y;y=z;z=w; return( w=(w^(w>>19))^(t^(t>>8)) );
}

const size_t ITER_PER_PRINT = 100000;
const size_t NUM_ITER = 20;

int monte(){
  size_t count = 0;
  for(size_t iter = 0; iter < NUM_ITER; iter++){
    for(size_t i = 0; i < ITER_PER_PRINT; i++){
      uint32_t v = xor128();
      uint32_t u = xor128();
      double x = v / 4294967295.0;
      double y = u / 4294967295.0;
      if(x * x + y * y < 1){
	count++;
      }
    }
    printf("%d: pi = %f\n", iter, double(count) / ITER_PER_PRINT / (iter + 1) * 4);
  }
}


inline double score(size_t a_, size_t b_){
  double a = a_ + 1, b = b_ + 1;
  double n = a + b + 2;
  return double(2 * a * b + n) / n / (n + 1) / (n + 2);
}

int monte2(){
  size_t WIDTH = 100;
  size_t N = WIDTH * WIDTH;
  vector<size_t> denom(N);
  vector<size_t> count(N);
  priority_queue<pair<double, size_t> > queue;
  for(size_t i = 0; i < N; i++){
    queue.push(make_pair(score(0, 0), i));
  }
  for(size_t iter = 0; iter < NUM_ITER; iter++){
    for(size_t i = 0; i < ITER_PER_PRINT; i++){
      uint32_t v = xor128();
      uint32_t u = xor128();
      size_t index = queue.top().second;
      queue.pop();
      double x = v / 4294967295.0 / WIDTH + 1.0 / WIDTH * (index % WIDTH);
      double y = u / 4294967295.0 / WIDTH + 1.0 / WIDTH * (index / WIDTH);
      denom[index]++;
      if(x * x + y * y < 1){
	count[index]++;
      }
      queue.push(make_pair(score(count[index], denom[index] - count[index]), index));
    }
    double area = 0.0;
    for(size_t i = 0; i < N; i++){
      area += 1.0 / N * count[i] / denom[i];
    }
    printf("%d: pi = %f\n", iter, area * 4);
  }
}

int main(){
  monte();
  monte2();
}