モンテカルロ法とは，乱数を利用してシミュレーションを行う技法の総称である．第２次大戦中，ロスアラモス研究所のフォン・ノイマンやウラムは，中性子が物質の中を動く様を知るために，乱数を用いたコンピュータ実験を行った．これがモンテカルロ法が実用化された最初の応用例である．

　

　彼らはこの方法の暗号名をモンテカルロ法と呼んだが，この呼称はカジノで有名な地：モンテカルロから由来している．もっとも，ラスベガス法という確率論的技法もあるのだが，迂闊（粗忽？）にも両者の違いについては失念してしまった．

　

　ところで，コンピュータを用いた数値実験とは，モデルとなるシステムを設定してさまざまな現象に関する模擬実験を行なうことであるが，実験が不可能な状況下における現象の情報を提供してくれる唯一の手法であり，今後ますます重要性が増すものと考えられている．たとえば，人体のように臓器構成が複雑で直接実験の対象とすることに危険を伴うもの，飛行機の翼のように大規模で実験をするには経費のかさむものに対して，これと類似の構造をコンピュータ内部に実現できればもっとも望ましい条件設定が得られるものと期待できるからである．

　

　現在，コンピュータシミュレーションは実験を行なわなくともシミュレーションで済む問題に汎用されていて，実験と理論につぐ第３の科学（計算機実験学）と呼ばれる分野を開拓し，実験と理論が互いに密接に関係しあって発展してきたこれまでの自然探究の歴史に変革を迫るものとなった．実際，理論を検証する実験が不可能な分野ではコンピュータを利用して極限状況を模擬実験し，理論と数値シミュレーションによって観測結果を近似することが可能になってきたのである．

　

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

　

　さて，ある論文に掲載されていたモンテカルロ法の適用にふと疑問を感じたことが，このコラムを書くきっかけとなった．今回のコラムでは，モンテカルロ法の誤用例を取り上げる．要するにデタラメな論文だったというわけであるが，読者諸賢もこの論文のウソを見破れるだろうか．

　

　論文で扱っている問題は，母数θが正規分布Ｎ（ｍ，ｓｅ^2）に従うとき，関数ｙ＝ｆ（ｘ，θ）の積分値

　　∫ｆ（ｘ，θ）ｄｘ

の標準誤差（ＳＥ）求めよという形に定式化できる．

　

　この問題に対するモンテカルロ法のプロトコールは，以下のごとくである．

　　(1)平均と標準誤差が所与の正規分布Ｎ（ｍ，ｓｅ^2）に従う正規乱数ｒを得る．

　　(2)ｙ＝ｆ（ｘ，ｒ）をｎ通り作る．

　　(3)これを数値積分して∫ｙｄｘをｎ通り得る．

　　(4)∫ｙｄｘの平均値Ｍと標準偏差ＳＤを求める．

　　(5)ＳＥ＝ＳＤ／√ｎとする．

　

　何の問題もないように見えるが，このようにして得られた積分値の誤差は，シミュレーションの回数ｎによって，いくらでも０に近づけることができる．ところが，誤差伝播の理論によると，

　　ＳＥ＝ｓｅ×∫∂ｙ／∂θｄｘ

であり，積分値の誤差はこれ以上小さくできない限界値として一意に定まる．これは明らかに矛盾である．

　

　データ解析には，テクノロジカルな側面とフィロソフィカルな側面があるわけだが，誤差をいくらでも０に近づけられるならば，その方法自体が意味のないことになるから，この論文に取り上げられた方法のフィロソフィカルな間違いは，モンテカルロの適用にあると考えられた．

　

　この間違いは，N(mean,SD^2)に従って乱数を発生させると，SEはいくらでも小さくできることに起因している．乱数発生のときのinvariantはSEではなくてSDなのである．このことから，求めたい誤差ＳＥは母数θの誤差ｓｅに対応した積分値の誤差ということであって，ステップ(4)のＳＤが目的とするＳＥなのである．すなわち，ステップ(5)は不要ということになる．

　

　ただし，モンテカルロ法で得られた結果は最適なものとはいいがたく，多くの場合，粗雑である．実際に計算すると，

　　ＳＥ＞ｓｅ×∫∂ｙ／∂θｄｘ

となるが，これがモンテカルロ法の限界であるから満足せざるを得ないであろう．

　

　一方，この論文のテクノロジカルな面については，なぜ誤差伝播の法則の積分型（コラム「積分値の誤差と誤差の積分値」参照）を導き出せなかったのか？という疑問が残る．

　

　統計的な手法を用いて，integrated form of propagation of errorを導き出せれば，計算精度，計算速度ともモンテカルロ法よりも優れた結果が得られるのだが，しかし，それは誤差伝播の法則に精通してはじめて可能になることであって，一般的には到底無理なことであろう．

　

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

　

　モンテカルロ法・ラスベガス法を実行する場合，不可欠な道具は乱数である．そこで，乱数発生法をいくつか紹介してみたい．

　

【１】逆関数法による乱数発生法

　

　確率変数xが確率密度f(x)をもつとすれば，y=F(x)は区間(0,1)で一様分布します．

（証明）

ｙの密度関数をg(y)とする．

　　f(x)dx=g(y)dy

また，y=F(x)よりdy={F(x)}'dx=f(x)dx

これより，g(y)=1

　

　このことより，一様乱数列{yi}を発生させて，

　　yi=F(xi)すなわちxi=F^(-1)(yi)

を満たす数列{xi}は与えられた確率密度f(x)をもつ乱数列となります．例えば，変換y=-logxは一様乱数R(0,1)を平均値１の指数分布をもつ乱数に変換するのに用いられます．

　

（例）指数乱数

F(x)=1-exp(-x)　 F^(-1)(x)=-log(1-x)

したがって，x=-log(1-u)

1-uも(0,1)上の一様乱数であるから，x=-loguとする．

これより，平均値１の指数分布（自由度１のガンマ分布）が得られます．

　

　また，尺度母数θの指数分布，すなわち，平均値θの指数分布乱数は，変換

　　y=-1/θlogx

によって得ることができます．たとえば，y=-2logxは一様乱数を平均値２の指数分布に変換します．

　

（例）ガンマ乱数

　　y=-1/θlog(u1*u2*･･･*uk)

とすると平均値θ，自由度ｋのガンマ乱数が得られますが，ｋが大きいときには効率的ではなく，またｋが整数でないときにはこの方法は適用できません．

　

（例）レイリー乱数・ワイブル乱数

　レイリー分布やワイブル分布は指数分布に基づく分布と考えることができますが，

　　z=(-2logx)^(1/2)

と変数変換すると

x=exp(-z^2/2),dx=zdz,p(x)dx=zexp(-z^2/2)dzよりｚの分布はレイリー分布となります．

　

　レイリー分布は形状母数２のワイブル分布であり，同様にして尺度母数β，形状母数αのワイブル乱数は

　　ｘ＝−β{log(u)}^(1/α)

として発生させることができます．

　

（例）ガンベル分布

　　F(x)=exp(-exp(-x))　 x=-log(-logu)

（例）コーシー分布

　　F(x)=1/πtan-1x　 x=tan(πu)

（例）ロジスティック分布

　　F(x)=1/(1+exp(-x))　 x=log(u/(1-u))

（例）最大値分布

　　G(x)={F(x)}^n　 x=F^(-1)(u^(1/n))

（例）最小値分布

　　G(x)=1-{1-F(x)}^n　 x=F^(-1)(1-u^(1/n))

　

　このように，逆関数による乱数発生法は原理的に簡単なので，指数分布，ワイブル分布，二重指数分布，コーシー分布，ロジステック分布などの乱数発生に用いられています．ただし，正規分布のように逆関数が簡単に計算できない場合も多く，また，逆関数法による乱数発生は高速でないという欠点もあります．そのため，任意の統計分布を対象とした高速乱数発生法が多数提案されています．

　

【２】ボックス・ミューラーの正規乱数発生法

　

　正規乱数を発生させる方法に，ボックス・ミューラー（Box-Muller）法があります．

　

　z1,z2が正規分布N(0,1)にしたがい，独立のとき(z1^2+z2^2)^(1/2)はレイリー分布にしたがいます．したがって，レイリー乱数を発生させることができると正規乱数に変換することができます．ｒを一様乱数とすると(-2lnr)は平均値２の指数分布，したがって(-2lnr)^(1/2)はレイリー乱数，また2πrは(0,2π)の一様乱数になります．

　

　r^2=-2logr1,θ=2πr2ですから，これにより，２個の一様乱数ｒ1，ｒ2から互いに独立に標準正規分布に従う２個の正規乱数ｚ1，ｚ2

　　ｚ1＝（−２ｌｎｒ1）^(1/2)ｃｏｓ（２πｒ2）

　　ｚ2＝（−２ｌｎｒ1）^(1/2)ｓｉｎ（２πｒ2）

を作りだすことができます．このように，ボックス・ミューラー法は，標的問題の解であるレイリー分布を応用していると考えることができます．

　

【３】中心極限定理を利用した正規乱数発生法

　

　特性関数は積率母関数同様に和の分布を求めるときなど利用されていますが，ここでは，区間（０，１）の一様分布の特性関数が

　　φ(t)=exp(it/2)sin(t/2)/(t/2)

となることを利用して，一様乱数ｒi（０〜１）をｎ個合計したものの分布が，ｎ→∞の極限で正規分布になることを示してみましょう．

　

　一様乱数をｎ個の合計のしたものの分布の特性関数は

　　[φ(t)]^n=exp(int/2){sin(t/2)/(t/2)}^n

一方，シンク関数

　　ｓｉｎｘ／ｘ＝Σ(-1)^mx^2m/(2m+1)!

　　　　　　　　＝１−１／３！ｘ^2 ＋１／５！ｘ^4 −・・・

の解が±π，±２π，±３π，±４π，・・・となることを利用して，無限積表示すると

　　sinx/x=(1-x2/π2)(1-x2/4π2)(1-x2/9π2)(1-x2/16π2)･･･

　　　　　=Π(1-x2/k2π2)

　

　ｋが大きいとき

　　(1-x2/k2π2)〜exp(-x2/k2π2)

　　Π(1-x2/k2π2)〜exp(-x2/π2Σ1/k2)

ここで，(Σ1/k2=π2/6)より

　　Π(1-x2/k2π2)〜exp(-x2/6)

これより，

　　{sin(t/2)/(t/2)}^n→exp(-nt2/24)　

ですから，

　　[φ(t)]^n→exp(int/2-nt2/24)

　

　正規分布Ｎ（μ，σ^2）の特性関数はexp(iμt-σ2t2/2)ですから，この結果はｎ個の独立した一様乱数の和の分布は平均値ｎ／２，分散ｎ／１２の正規分布に近づくことを示しています．

　

【補】(1+x/n)^n→exp(x)

【補】Σ1/k2=π2/6＝ζ(2)

　　sinx/x=1-1/6x2+120x4-･･･（ベキ級数表示）

また，

　　sinx/x=Π(1-x2/k2π2)　　（無限積表示）

　　　　　=1-1/π2(Σ1/k2)x2+･･･

の両辺を比較することにより，Σ1/k2=π2/6，Σ1/k4=π4/90が計算される．Σ1/k2はゼータ関数ζ(2)に，Σ1/k4はゼータ関数ζ(4)に相当する．

　

　次に，一様乱数をもとに正規乱数を発生させる具体的な方法を述べてみましょう．中心極限定理によって，区間（０，１）の一様乱数ｒi（平均１／２，分散１／１２）をｎ個合計したものの分布は平均値μがｎ／２，分散σ^2がｎ／１２の正規分布に近くなりますから，正規化して

　　Ｚ＝√１２／ｎ（Σｒi−ｎ／２）とおけば，

Ｚの分布は標準正規分布Ｎ（０，１）となります．

　

　そこで，１２個の一様乱数を加えることにすると平方根の計算をしないで済みます．

　　Ｚ＝（Σｒi−６）　　　（ｉ＝１〜１２）

実際，一様分布に従う確率変数１２個ずつの平均をとり１００個のデータから構成したヒストグラムは元の一様分布とは似ても似つかない滑らかな分布となります．

　

　なお，中心極限定理を利用した正規乱数発生法は，１２個の一様乱数から１個の正規乱数が得られる効率のよくない方法で，それに比べ，ボックス・ミューラー法はかなり効率がよくなっています．

　

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝