１．乱数の利用−−−モンテカルロ法
計算機において乱数といえば，ゲームや暗号をすぐに思い浮かべるかもしれないが，科学技術計算においても乱数は有用である。乱数を用いる計算を一般にモンテカルロ法と呼んでいる。要するに，サイコロ振って．．．というだけのこと。
(1) モンテカルロ・シミュレーション
熱運動を始めとする確率的現象や不確定要素を除けない現象を，問題にふさわしい確率分布をもつ乱数を利用して行うシミュレーションである。自然現象のみならず，交通流の問題であるとか，商品流通と在庫管理の最適化，ビジネスゲームなど，工学や社会科学においてもモンテカルロシミュレーションは盛んに利用されている。
例１ ランダムウォーク（酔歩）
ｄ次元立方格子上を，いずれかの方向に±１だけ，等確率 1/2ｄ で飛び移る「酔っぱらい」の運動である。少し長い目で見ればブラウン運動 （ △t 間の遷移距離 x がガウス確率 exp[-x2/2△t] をもつような確率的運動） に一致する。（中心極限定理）
無限に拡がった空間におけるランダムウォークの解 （ｔ 秒後の位置の確率分布） の関数形は簡単に与えられる。例えば１次元の直線上を，酔っぱらいが t=0 に x=0 を出発した場合，t 秒後に x にいる確率は
P(t,x)＝tC(t+x)/2 / 2t＝2-t t!/[(t+x)/2]! [(t-x)/2]! （２項分布）
で与えられる。［ t 回のステップのうち，(t + x)/2 回が＋１，(t - x)/2 回が−１のとき，t 秒後の位置は x となる。ｄ次元についても同様な形を与えることができる。］ t が十分大きいとき，これは
P(t,x)≒(2πt)-1/2 exp[-x2/2t]
と，ガウス分布で近似される。
しかしながら，吸収壁とか反射壁とかがあるような場合には，簡単に解の形を与えることは一般にはできない。
実際にシミュレーションを行うときには，例えば３次元であれば，(0,1) 実数区間を 1/6 ずつの部分区間に分けておき，(0,1) 間の一様乱数を発生させて，その値がどの部分区間に属するかで，飛ぶ方向を決める。
0 0.5 1
|---------|--------|--------|--------|--------|--------|
+x方向 -x方向 +y方向 -y方向 +z方向 -z方向
例２ イジングモデル
強磁性体や合金，ときには液体-気体の間の相転移をも表現するのによく用いられる統計力学モデルである。ｄ次元格子点上に±１の値をとる変数（イジングスピン） si が配置されており，全系の配向の組 ｛si｝ により系のエネルギー
E = - H i si- J (ii') si si' (ii') は隣接する格子点の対
が決まる。各スピン si は熱運動により，符号を反転したとした場合に増えるエネルギー
ΔEi = 2si ( H + J i' si' )
の符号に応じて，以下の確率で符号を反転する：
exp(-ΔEi/kT) if ΔEi ≧ 0
P(si→ -si) =
1 if ΔEi ＜0
これを繰り返すことにより，十分に時間がたてば温度 T の熱平衡状態が実現される。全スピン変数が平均して１回更新 （update） される時間を単位にして １MCS (Monte Carlo Step) という。
まずスピンの位置を乱数でランダムに選び，反転した場合のエネルギー差ΔEiを計算する。次に，上の規則に基づいてそのスピンをupdateする。「確率 p である操作を行う」 には，(0,1) 間の一様乱数 r を発生させ，「0＜r≦p」に入るかどうかで判断する。
0 p 1
|---------------------------|-----------------|
DO UNDO
(2) モンテカルロ数値計算法
例１ サイコロで定積分を求める
ｄ次元の多重積分
1 1 1
I = ∫∫・・・∫ f(x) dx
0 0 0
を，ｄ次元一様乱数（ベクトル）を用いて
I' = (1/N) i f(xi)
で近似することができる。図 （円の面積） を見る方が分かりやすいかも知れない。乱数が真に乱数の性質（一様である，相関がない，etc）を持っているならば，N 回の試行による誤差は N-1/2 のオーダーである（中心極限定理）。１億回であれば１万分の１（図の例では 3.1415±0.0003 程度）である。
この方法は多重積分になると俄然，有利になる。素朴な単純区分求積法では，１次元方向の区分数を L （したがってセル数 N は d 次元積分なら Ld） として，誤差は L(d-1)×L-d = L-1 = N-1/d 程度は覚悟しなければならない。これに対してモンテカルロ法では次元によらず N-1/2 程度である。N 個の各点が面積を求めたい領域内に含まれるかどうかを調べるのに，乱数発生のために各点で数倍の計算ステップを要するとしても， N が大きい多重積分では原理的に「サイコロ」の方が有利なのである。
なお，多重積分という目的に限れば，誤差が N-1/2 ではなくて N-1 あるいは N-2 のオーダーになる「準乱数」が開発されており，３重積分以上ではどのような区分求積法やガウス点法と呼ばれる方法よりも効率がよいとされている。
例２ サイコロで静電ポテンシャルを求める
静電ポテンシャルはラプラス方程式を満たす：
▽2u(x,y,z) = 0 ただし ▽2 = ∂2/∂x2 + ∂2/∂y2 + ∂2/∂z2
これは静電気の問題だけでなく，熱工学，流体力学などいろんな場面で出てくる方程式であり，この解を調和関数という。境界条件として
閉曲面Ｂ上でのポテンシャルの値 u = b(xB, yB, zB)
が与えられた問題をディリクレ問題という。一般には，平面や球のような幾何学的対称性のよい，ごく限られた問題しか解の具体的な形はわかっていない。任意の閉曲面境界の問題を，ランダムウォークで解くことができる。
空間を立方格子状のメッシュに分割し，境界面は格子点を通るように近似する。解を求めたい点Ｐから酔っぱらいを出発させ，境界点に達したらその点での境界値 b を賞金として与えて酔歩を止める。これを繰り返したときに酔っぱらいが得る賞金の期待値 u が，点Ｐにおけるポテンシャルの解である。
（理由） ランダムウォークの性質から
点 P から出発したときの期待値 ＝ ｛次の瞬間に P の隣接点 N に飛ぶ確率｝
×｛点 N から出発したときの期待値｝
すなわち，Ｐ点を (i,j,k) として
u(i,j,k) = (1/6) [u(i+1,j,k) + u(i-1,j,k) + u(i,j+1,k) + u(i,j-1,k)
+ u(i,j,k+1) + u(i,j,k-1)]
である （マルコフ性：点 N から出発するときの運動は，それ以前にどこにいたかには依存しない）。これは調和関数の性質
点 P での値 ＝ P を中心とする球面上での値の平均値
の離散化版に等価である。すなわち，上式を書き直せば
[u(i+1,j,k) - 2u(i,j,k) + u(i-1,j,k)] + [u(i,j+1,k) - 2u(i,j,k) + u(i,j-1,k)]
+ [u(i,j,k+1) - 2u(i,j,k) + u(i,j,k-1)] = 0
となり，第２回の講義メモにあるように，第一項は x についての２階微分，第二項は ｙ についての，．．．であって，関数 u はラプラス方程式を満たすことになる。また，境界上の点から出発した場合は，そこでの賞金 b を得てただちに止まるから期待値 u は明らかに b に等しい。（終り）
別の方法も可能である。ディリクレ問題の解は，接地された境界を持つ点電荷のポテンシャル
▽2 GP(x,y,z) = - δ(x-xP)δ(y-yP)δ(z-zP) δはデルタ関数
境界条件： 閉曲面Ｂ上で GP = 0
を用いて
u(xP,yP,zP) = - b(x,y,z) n･▽GP(x,y,z) dS （閉曲面Ｂ上の面積分）
と表すこともできる。（グリーン関数の方法）
このグリーン関数 GP(x,y,z) は，点Ｐから出発して境界Ｂに到達すれば消滅する 「吸収壁を持つランダムウォーク」 が，任意の点 (x,y,z) を訪れる回数の期待値によって与えることができる。
（理由） こちらはブラウン運動で説明する方が簡単である。t=0 に点Ｐを出発した酔っぱらいが，時刻 t に点 (x,y,z) にいる確率を p(x,y,z,t) とすれば，関数 p(x,y,z,t) は次の拡散方程式を満たす（※）。
∂p(x,y,z,t)/∂t = ▽2 p(x,y,z,t)
p(x,y,z,0) = δ(x-xP)δ(y-yP)δ(z-zP)
両辺を t=0 から t=∞ まで積分し，訪問回数の期待値が以下で与えられること
∞
GP(x,y,z) = ∫ p(x,y,z,t) dt
0
また吸収壁のため，酔っぱらいは時間がたてばいずれ必ず消滅してどこにもいなくなる （ p(x,y,z,∞) = 0 ）ことを用いればよい。境界（吸収壁）では酔っぱらいは直ちに消滅するから，境界上にいる確率は常に p = 0 ，したがって GP = 0 である。（終り）
（※）境界がなければガウス分布がこれを満たすことに注目せよ。あるいはランダムウォークの関係式（＝６つの隣接点から1/6ずつの確率でやってくる）
p(i,j,k; t+1) = (1/6)[p(i+1,j,k; t) + p(i-1,j,k; t) + p(i,j+1,k; t)
+ p(i,j-1,k; t) + p(i,j,k+1; t) + p(i,j,k-1; t)]
から，先ほど同じようにして導くこともできる。
(3) 乱数の発生
そもそも規則ずくめの計算機で乱数を発生するなんて，もとより無理難題なのである。いかにして実用的には乱数らしきもの−擬似乱数−を実現するか，である。モンテカルロ法にとって大切なことは，乱数としての性質が良いことはもちろんであるが，計算速度が速く使用メモリもできるだけ少ないことである。算術規則に基づいて発生させるため再現性があることは，シミュレーションにとってはありがたい。
ここではよく使われている発生法を２，３紹介するにとどめ，あとは各種の乱数発生法・検定などについて書かれた比較的読みやすい参考書をあげておく：
津田孝夫「モンテカルロ法とシミュレーション」（培風館）
伏見正則「乱数」（東京大学出版会）
(3-1)線形合同法
数列の過去の要素をいくつか使った一次式で新たな要素を生成し，整数列を発生する方法であるが，比較的簡単でメモリも最小限ですむのは，一つ前の整数一つから次の整数を発生させる方法である：
X(n) = A X(n-1) + C mod M （ mod M は M で割った余りを表す。）
初期値（シード）はなるべく大きな数値を与える。
X を M で規格化した実数数列を (0,1) 間の一様乱数として用いる。
積の上位桁の数字を捨てる（余りを採用する）ことにより規則性が除かれると思えばよい。（⇒『運動の決定論性と予知不能性』） この場合，原理的に周期は高々 M であるため
M をできるだけ大きくすること
できるだけ可能な最大周期を保証すること
について様々な工夫が行われた結果，定数 A や C にある条件が課せられるが，あとは経験的に検定結果が良いと言われている定数の組み合わせを信じて使うしかない。
M としては 2^32 （４バイト整数の限界）や 2^31 （同 正整数の限界）がよく使われる。これは周期をシステムで可能な最大にすることができる候補であると同時に，以下のように mod 計算を簡素化できるという理由からである。
M=2^32 の場合：４バイト表現で，上位へのオーバフローをエラーと
しないシステムでは，何もしなくても mod 計算になる。
M=2^31 の場合：整数 2147483647=[011111111 11111111 11111111 11111111]2
との「ビットごとの AND 演算」をとることにより余りを取り
出せる。（この種の計算をマスク演算という。）
よく使われている例の一つ
M = 2^31 A = 5^11 C = 0 （ C=0 のとき乗算合同法という）
＜Fortranプログラムの例＞
INTEGER, PARAMETER :: a = 48828125, mask = 2147483647
REAL, PARAMETER :: rn = 2.147484E+9
INTEGER :: x
REAL :: r
.....
（x の初期値を与える： 手で入力する，クロックを利用する，など）
.....
(以下，繰り返し）
.....
x = IAND( a*x, mask) ! IAND はビットごとの AND計算をする組込関数
r = x/rn ! (0,1)間の実数に規格化する
.....
（注）定数 mask を「mask = 2**31 -1 」で与えるとエラーとなる。
IAND を使わず単に x = a*x とすれば(-1,1)の一様乱数になる。mask は不要。
＜乗算合同法の使用上の注意＞
(a) このまま整数乱数として使用してはならない。(0,1,2,...,m-1) の一様整数乱数が必要なときには，区間を m 等分し，どの区間に落ちるかで決める。
(b) 高次元乱数として使うときは要注意！ ［例えば，( X(n), X(n+1), X(n+2) )を３次元プロットすると顕著な格子構造が現れ，一様性が失われる。］
(3-2) Ｍ系列法（最大周期列法）
乗算合同法では，ある整数の次に来る整数は一意的に決まってしまい，周期はたかだか M = 2^31 とか 2^32 であり，膨大な乱数列を必要とするシミュレーションでは困ることがある。この周期に関する難点を改良した方法である：
X(n) = X(n-p) .XOR. X(n-q) （.XOR. はビットごとのExclusive-OR 演算。
Fortranには IEOR(i,j) という関数がある。）
初期配列（p個のシード）の与え方 （p＞q とする）
1) 先ず p 個の全てが 0 ではないビット列 {b(i), i=0,1,2,...,p-1} を任意に与える。
2) b(i) = b(i-p) .XOR. b(i-q) (= [b(i-p) + b(i-q)] mod 2) により，ビット列を
あわせて 32p 個まで生成する。
3) 先頭から32ビットずつ用いて p 個の4バイト整数 {X(i)} を生成し，初期値とする。
p=127 の場合：
X(0) = [b(31) b(30) .......... b(1) b(0) ]2
X(1) = [b(63) b(62) ...........b(33) b(32)]2
X(2) = [b(95) b(94) .......... b(65) b(64)]2
X(3) = [b(127)b(126) .........b(97) b(96)]2
X(4) = [b(159)b(158) ......... b(129)b(128)]2
.................................
.................................
X(126) = [b(4063) ................ b(4032)]2
(p,q) の例： (89,38), (127,7), (127,63), (521,32), (521,158)
（参考）ビット演算
(AND) (OR) (XOR)
| 0 1 | 0 1 | 0 1
---+-------- ---+-------- ---+--------
0 | 0 0 0 | 0 1 0 | 0 1
| | |
1 | 0 1 1 | 1 1 1 | 1 0
４バイト整数の各位のビットが並行して生成されるＭ系列ビット列になっており，初期配列を上で述べたように与えることによって最大周期 2p-1 が保証される。周期の点では実用的に問題はないが，使い方によっては有意な歪みが現れることが知られている。（大偏差異常→課題１_３）
(3-3) ＭＴ法（Mersenne Twister 法）
松本眞氏（総合人間学部数理基礎論講座）の開発によるもので，目下，世界中で最も良い性質を持っているとされている。ただし，アルゴリズムが少々複雑で，私自身まだ消化できていないので説明は省略する。興味のある人は以下のページを参照のこと。ソースプログラムも提供されている。 松本先生のページ
第一回課題 −−−少なくとも１題はレポートを提出すること
ソースプログラムと見やすくデザインされた出力結果を一つのファイルにまとめてemailで提出するか、紙に印刷して提出すること。
１_１． 半径１の円の面積をモンテカルロ法によって求めることにより，円周率を推定せよ。
・ 試行回数 N=1000, 10000, 100000, ..., 100000000（1億） について，それぞれ10回ずつ繰り返し，10回の平均値と偏差を求め，偏差が N とともにどのように減っていくか調べよ。［偏差＝平均値からの ずれ の２乗の平均値の平方根］
・ (-1,1) の一様乱数は，M=2^32, A=5^11 の乗算合同法により得られるものを使えばよい。ただし，ｘ座標，ｙ座標は，それぞれ別々の初期値（シード）から出発する２つの独立な乱数列を用いること。
１_２． 以下のいずれかの一次元ポテンシャルの解をモンテカルロ法により求め，解の形を適当な縮尺を用いてヒストグラムで表せ。(2)の出力例 ［注） 一次元の調和関数 （ d2u/dx2=0 の解） は x の一次式であるから，解の形は予め分かっている。］
(1) x=±100 におけるポテンシャルの値（例 u(100)=20, u(-100)=10 ）が与えられたときのポテンシャル u(x)。 ただし，(-100, 100) の直線上に電荷はないとする。
(2) x=0 に点電荷１が置かれ，x=±100 が接地 （u=0）されているときの解 u(x)
１_３． Ｍ系列乱数を用いてコイントスを行い，「 p 回のトスを静観し，p 回のうち，もし 60％ 以上の表（オモテ）が出るのを確認できたときだけ，次のトスで表（オモテ）に賭ける」ことによって，どれくらい儲かるか，(p=127, q=63) を用いて100万回実験し期待値を求めてみよ。［松本眞「コイン投げで一儲けする」情報処理学会誌１９９８年１１月 より］
・ （表，裏）を (0,1) で表す１ビット列でよいから，最初の p 個だけ {0,1} のランダムな列を，例えば乗算合同法で与えればよい。a，b が１ビットのときの XOR演算は算術演算「 (a-b)**2 」で与えることができる。
・ p=127 回のトスで 60％、すなわち76回（10/29訂正）以上の表が出る確率は，２項分布で計算できて，100万回のうち16394回くらいと期待される。また，16394回賭けたときの勝ち負けが 8197±128 程度なら正常な揺らぎである。
・ 配列はｐ×100万個用意する必要はない。｛b(0), b(1), ... ,b(p-1)｝ のｐ個で判定を終え、次の新たなｐ個を用意するには、
「 b(0) と b(p-q) から次のビットを作れば，b(0) はもう用がなくなるから，これを新たに b(0)」 とし，「 b(1) と b(p-q+1) から新たに b(1) を」．．．
とすればよい。ただし「 b(q-1) と b(p-1) から新たに b(q-1) を」から先は数列がなくなるから，先頭へ戻るための工夫がいる。わかりにくければ、先ず b(0)〜b(126) を用いて b(127)〜b(253) を生成してから、全体を 0〜126へ平行移動したと思ってもよい。
（特別付録）コイン投げで、ある回数以上のオモテが出る確率を求めるプログラム（Windowsのexeファイル)