2013年03月28日

モンテカルロ法と擬似乱数

コンピュータ将棋つながりで最近いろんなCOMの思考ルーチン見てたんですが、
最近(もっと前から?)囲碁の方でモンテカルロ法が流行ってるそうです。

モンテカルロ法ってのは簡単に言えば
「答えを出す方法はあるが答えを出す式が分からない」系統の問題に、
手当たり次第適当に答え出しまくって確率で判断するって理論です。

よくある例えが円周率を求める例えで、

1.(0,0)-(1,1)の正方形に適当に点を打つ
2.0.5,0.5を中心に半径0.5の範囲に入っていればそれは円の中である
3.打った点の位置から2.であるか否かは分かる
4.打った点が円の中である点の数/打った点の総数は円の面積の近似である

こんな乱暴な理論なんですが、試行回数が無限になれば誤差0になることは証明されてます。

こういう単純明快で乱暴でそれでも正しい理論って大好きなので、
もうちょい深く調べてたんですが、モンテカルロ法で超重要な「適当」について
いろいろと思い出したことと知ったことがあります。

「適当」ってのは言い換えれば「乱数」なんですが、コンピュータの乱数
(正確に言うと処理系が提供してる疑似乱数生成関数)は大抵色々制限があります。

例えば範囲(返る乱数の種類)ですと
VBA:0〜2^24-1の整数(を16777216で割った浮動小数点)
VC :0〜2^16-1の整数(でもint型で返される)
しかありません。VCだと65536種類です。今のCPUだと一瞬で重複しますね。

周期も問題でして、こちらは VBA:2^24 VC:2^32 となっています。
周期ってのはその周期分取得すると、また同じものが1から出てくるってことで、
周期より多く乱数を取得すると、もはや乱数ではなくなってしまいます。

つまりどちらも、兆超える試行回数だと使い物になりません。
(シミュレーション系や、まさにモンテカルロ法だと必要になるんですよ)

さて、これらの疑似乱数生成関数は大体全部、線形合同法と言われる式を使っています。

こんなの。

X(n+1) = ( A * X(n) + B ) MOD M

※1)M>A M>B A>0 B>=0
※2)周期はA,Bを上手く設定すれば最大M

こんな単純な式で乱数っぽいのが作れるんだから凄いもんですが、
単純なだけに問題がありまして。

1.下位ビットがランダムではない
 作られた乱数に Mod 2 とかやると、0と1が交互に現れます。
 10進数で言うと偶数の次は必ず奇数になるってことですね。

2.組にして使うとランダムではない
 画面へ(Mod 256, Mod 256)で点を打つとかやると、
 永遠に打てない点が出たりします。

とまあ周期以前にえーって感じですが、
ここら辺をどう対応するかは処理系で変わります。

VBAはSingleで返すので、普通はRnd()*(2^24) Mod 16 じゃなくInt(Rnd()*16)でしょう。
VCだと上位ビットだけ返すようにしてます。(だから2^16しか種類が無いんです)

さて、まともな乱数作る方法がこの時代になっても無いのかって話ですが、
(上記でも数万回あたりまでは十分実用的なんですが)
まともな乱数の定義としては以下あたりでしょうか。

1.再現可能である
2.規則性が見られない
3.処理が速い
4.周期が長い
 
自分で考えたのが、円周率の小数点をn桁で区切って10^nで割ればよくね?っての。
円周率は無理数ですし、○桁から〜って指定すれば1.2.4.はいけそうだと思うんですが、
・・・superPIがベンチマークに使われてる以上、考えるまでもなく遅すぎるでしょうね。

で、今回始めて知ったのがXorShift。ソースこんなん。


// 自分用に書き換えてます。スレッドセーフじゃないです。
unsigned __int32 xss[ 4 ];
void xsinit( unsigned __int32 seed )
{
for( int i = 0; i < 4 ; i++ )
{
seed = 1812433253U * ( seed ^ ( seed >> 30 ) ) + i;
xss[ i ] = seed;
}
}
unsigned __int32 xorshift()
{
unsigned long t = ( xss[ 0 ] ^ ( xss[ 0 ] << 11 ) );
xss[ 0 ] = xss[ 1 ];
xss[ 1 ] = xss[ 2 ];
xss[ 2 ] = xss[ 3 ];
xss[ 3 ] = ( xss[ 3 ] ^ ( xss[ 3 ] >> 19 ) ) ^ ( t ^ ( t >> 8 ) );
return xss[ 3 ];
}


見れば分かりますが、ビットシフトとXORしかしてません。名前そのままです。
アルゴリズム的には線形合同法を128ビットに拡張して改造した感じ?(正直自信ない)
ぱっと見、相当早そうな感じです。Modとか無いですし。

で、このXorShift、周期は2^128-1です。これがどれぐらい長いかというと、
「4GHz16コアのCPUを持つPCを1億台繋げて1クロックで1乱数生成できて
スレッド・プロセス・ネットワーク間ロスはゼロ」だという
ありえない仮定を挙げたとしても、周回するまで15701億年かかります。
Wikipediaの時間の比較でもネタが無くなるレベルですね。

あとはメルセンヌツイスタ。これはアルゴリズムから違います。たぶん。
詳しくはよく分からなかったんですが、
1.周期が2^19937-1である
2.多次元において十分に不規則である
3.かなり凄い発明だと思うんですが普及していない
あたりは分かりました。

1の周期はもう一般的浮動小数点計算範囲から離れてるので諦めました。
たぶん宇宙が再配置とかそういう単位です。

2.3.はまあたぶん需要と供給の関係なんでしょう。
本気で乱数が重要な処理系ならハード乱数作るでしょうしね。
posted by Nick at 20:54| Comment(0) | TrackBack(0) | アルゴリズム