物理を題材にしたC その3
乱数を用いたシミュレーションを行ってみます。乱数は自然界のランダム性に支配される現象を計算機で再現しようとするとき、特によく使われる道具です。ランダム性とは、自然ですら制御できない、ましてや人間の制御の利かない、いわばサイコロを振って次の行き先を決めているような性質で、自然界の根本部分に関わっています。
まず乱数とは、普通一様乱数のことで、その名の示すとおり「乱れた；みだれた」数です。乱れたという意味は、何の規則性もなく乱れたということで、ある数が作られ次の数を予言できないほどでたらめだという意味です。しかし実は計算機では全く規則のない乱数は造れません。つまり順番は同じだが、発生された（作られた）数つまり、乱数としては十分に乱れている数の列を作り出すことができます。これを計算機乱数、あるいは疑似乱数とよびます。疑似とついているのは、順番に関してはある種の規則性があるので、完全な乱数ではないという意味です。計算機乱数では呼び出す度に異なる（疑似でたらめな）値を返す程度に理解して下さい。計算機により作られた乱数の列はいつかは同じ値を取ります。この後は同じ繰り返しとなります。繰り返しまでの回数が充分にお大きければ実用上、「真の乱数」とみなすことができます。
良く知られた線形合同法による疑似乱数の生成プログラムを次に示します。
線形合同法では初期値X0から初めて、Xnを次の漸化式に従って生成（計算）します。
ここで%は左側の（AXn-1+C）を右側のMで割った余りを計算する演算子です。A,C,Mは適当な整数ですが、Mが2^n （２のｎ乗）で、Aが８で割ると５余る数で、Cが奇数の時、この乱数の周期（もう一回同じ数に戻るまでの回数）は最大になることが分かっているのでこの条件を用います。ちなみに周期はMです。つまり、0から(M-1)までの整数が全部でM回発生する乱数のうち１回だけ現れます。 この条件を満たすA=109,C=1021,M=32767を使ってプログラムは書かれています。ちなみにXの初期値X0は13を用いました。
---------------------------------------------
#include <stdio.h>
int rndnum=13; /* initial number for the rand */
int irand();
main()
{
int j;
for (j=0;j<100;j++){
printf(" %8d ",irand());
}
}
irand() /* generate 0-32767 int rand */
{
rndnum=(rndnum*109+1021) % 32768;
return(rndnum);
}
---------------------------------------------
gcc irand.c -o irand でコンパイルをして実行プログラムirandを作ります。
irand として実行した結果は、-------------------
1021 14006 20347 23388 27177 14194 8071 28792 26389 26606 17491 6996 9921 1066 18911 30704 5421 2086 31787 25164 24153 12258 26423 30312 28229 30558 22275 4164 28913 6810 22415 19424 21085 5526 13531 1340 16009 9298 31463 22616 8565 17102 30131 8500 10017 11530 12607 31696 15245 24326 31115 17452 2745 5314 23191 5704 165 19006 8291 20004 18769 15226 22255 1984 20669 25718 19003 7964 17129 306 1607 12344 3029 3502 22291 5908 22401 17898 18591 28592 4589 9702 9963 5644 26393 27042 32247 9768 17157 3358 6595 31748 20913 19546 1615 13216 32541 9046 3995 10492
------------------------------------------------のようになります。1021 の次に14006 さらに20347 その次が23388と言う風です。非常に簡単な計算ではあるが、それらしい「乱数」が作られていることが見えます。もう一回同じプログラムを走らせてみて下さい。結果は同じですね。つまり同じ乱数列を発生します。これは int rndnum=13という初期値を用いているためです。プログラムをeditしてこの値を変更してみてください。できれば素数が望ましいです。今度は別の乱数列を発生しましたね。このように計算機では、シード（種）と呼ばれる初期値に依存した乱数の列を決めて発生させることで動かします。これは計算機で乱数を用いて計算させたとき、結果が違う場合、同じ乱数列で計算結果をチェックするために大変有用です。ただし自然界ではこの様な事は起こりませんから、毎回異なるシードを見つける工夫が必要です。例えば実行時の時刻とか、、、、皆さんも考えて使えるようにプログラムを変更してみて下さい。
しかしこれでは、本当に一様かどうかは良く分かりませんね。そこで、ヒストグラムを描いて見てみましょう。
ヒストグラムは、物理の世界では良く利用されるグラフです。同じようなことが何度も繰り返し起こる場合にその結果の値を集計して、値が適当に同じ大きさ範囲に入る時その繰り返しの数を縦軸に積み上げて行きます。横軸には値の大きさを取ります。ですから、適当な大きさの範囲を大きく取りすぎても、小さく取りすぎても意味のないグラフとなります。下のプログラムではヒストグラムのデータはUNIX-likeにファイルとして読みとります。またヒストグラム化は１０個のビンについて行うように作られています。
このプログラムをコンパイルして gcc hist.c -o hist -lm で出力したとき、上で作ったデータを
irand > irand.d で出力しておくと、 hist < irand.d として入力してやると、ヒストグラムが見れます。こんな原始的なプログラムを走らせているのは、配列の勉強とグラフィカルな環境を仮定しないでテキスト画面に表示できるものを敢えて使うためです。
-----------------------------------hist < irand.dの結果の画面表示
(..............inputs are .... という入力データの表示は省略しました。)
.
inputs are 99 10492.000000
min and max values 300 61.000000 300 32732.000000
＜＜＜入力データ数、その最大値、最小値が表示されます。
bin width 3267.100098 61.000000 32732.000000
＜＜＜この結果１０個のビンに分けるビン幅（値を適当な同じ大きさとみなす幅）
Ndata= 300, average= 16322.43, std = 9812.12
＜＜＜平均値と標準偏差を計算してあります。 ここからしたがヒストグラムです。データの領域、個数、これを棒グラフで長さに換算してあらわしたものという順番で表示されます。
61.00 - 3328.10 39 |||||||||||||||||||||||||||||||||||||||
3328.10 - 6595.20 31 |||||||||||||||||||||||||||||||
6595.20 - 9862.30 25 |||||||||||||||||||||||||
9862.30- 13129.40 22 ||||||||||||||||||||||
13129.40- 16396.50 26 ||||||||||||||||||||||||||
16396.50- 19663.60 35 |||||||||||||||||||||||||||||||||||
19663.60- 22930.70 23 |||||||||||||||||||||||
22930.70- 26197.80 34 ||||||||||||||||||||||||||||||||||
26197.80- 29464.90 36 ||||||||||||||||||||||||||||||||||||
29464.90- 32732.00 28 ||||||||||||||||||||||||||||
＜＜＜＜これがヒストグラムです。
-----------------------------------------------------------------
このヒストグラムをみると、とても一様分布とは言えません。この場合３００個の乱数を発生させました。もう一桁あるいは２桁ふやして実行してみてください。数が多くなると一定らしく見えてくるでしょう。これを大数の法則と言います。
上のヒストグラムをデータ（テキストファイル）から読みとり画面に表示させるプログラムの中心部分です。
-----------------------------------------
/* hist draw */
for(i=0;i<nbin;i++) /* clear histgram buffer */
{ h[i]=0; } ヒストグラムバッファh[100]という配列が用意されています。これの全成分をゼロにします。
for(i=0;i<nbin;i++) /* loop in the bins */ nbinはヒストグラムの棒の数です、各棒のことを
ビンと呼びます。
{
b=bins+i*binw; /* lower edge of the bin */
bn=b+binw; /* higher edge of the bin */
for (j=0;j<nmax;j++) /* loop in the all data */
{
if(d[j]>=b && d[j]<bn){ /* find data in THE bin */ 読み込んだデータはd[]という配列に
保存してあります。特定のビンに該当するとき
h[i]=h[i]+1; /* accumulate in the hist buffer */ ヒストグラムの値を１増やします。
}
} これをすべてのデータについて回します。
/* final histgram printing */ 最後に画面に印刷します。グラフらしく見せるために
printf(" %5.2f- %5.2f %3d ",b,bn,h[i]); /* data values */
for(j=0;j<h[i];j++) { printf("%c",charh);} /* markers */ このループがあります。
printf("\n"); /* to the next line */
----------------------------------------------------------------
汎用性があると思いますので、いろいろな計算結果を簡単にヒストグラムを表示させたい時に、データをファイルに出力し（redirection)、これをhistが拾って表示させると状況がわかります。また変数は自分で理解して変更する、あるいは入力時に変更可能にすると良いでしょう。
さて乱数を使った物理は、統計力学が得意とするところです。そこでまず手始めに、ブラウン運動やランダムウオークあるいは酔歩と呼ばれるシミュレーションを実行してみましょう。これは、酔っぱらいが歩く姿にも見えるため、この名が付けられました。
まず簡単な１次元で考えます。つまり酔っぱらいは右か左へ行くかの選択枝しかありません。一様乱数も上で作った乱数は０〜３２８６８までの整数でしたが、このような場面では0〜1.0の間を一様に分布する乱数が望ましいので、irand()を３２７６８で割ってやります。ただし、整数を整数で割って
ゼロか１しかでない、なんてやらないでくださいね。こういうときは
float randn=(float)irand()/(float)32768
としてください。実数randnが現れます。さて選択肢は同じ確率の２つですから、0から1.0の一様乱数が出たとき、0〜0.5までならたとえば「右へゆく」、反対に0.5から1.0の間なら「左へゆく」と決めれば、よろしい。ｚを酔っぱらいの座標としましょう。
プログラムの重要な部分は、右や左へ１歩を１の大きさにとると、
for(i=0;i<nmax;i++) { nmaxは最大歩数です。すでによそで与えられています。
if(randn>0.5) z=z-1; 左はマイナス方向なので、
else z=z+1; それ以外は必ず右で、右はプラス方向です。
}
で十分です。まずprintf(" %d %d \n",i,z)文で位置座標を刻々と打ち出してやりましょう。さてこれを動かしてみますが、数値を見ても全然楽しくありません。これを実行して、redirection > で別のファイルへ出力してグラフで楽しんでみましょう。
たとえば、この酔歩プログラムをrandw.c という名前で作り、コンパイルしたとしましょう。
gcc randw.c -o randw
としてコンパイルを実行し、実行ファイル名をrandwとして作ります。
実行を randw > randw.d として、ファイル名randw.dへ出力します。
ファイルの中身はたとえば、more randw.dで見てみると
0 0
1 1
2 0
3 1
4 0
5 -1
6 -2
7 -1
8 0
9 1
10 2
11 3
12 4
13 3
14 4
15 3
..........のようになっていたとします。これは第１列が何歩目か表し、第２列が位置を示します。
この例の場合、始めは原点z=0にいました。次に右z=1へ移動します。次の１歩は左へ動くので、また原点の戻ってしまいます。。。
さてこれをgnuplot でグラフ化してやりましょう。
plot 'randw.d' とすると、つぎのようなグラフ（プロット）をえることができます。
次に気になるのは当然２次元平面内での酔っぱらいの足取りですね。その場合、酔っぱらいがとれる選択肢は（２次元平面を簡単のため格子点上だけの世界にしてしまいましょう、そうすると）４つあります。つまり、前、後、右、左です。0から1.0の一様乱数を使って、0〜0.25までならたとえば「右へゆく」、0.25から0.5までなら「前に行く」、0.5〜0.75なら「左に行く」、0.5から1.0の間なら「後ろへゆく」と決めれば、良いでしょう。酔っぱらいの座標は２次元ですから、(x,y)とかかれます。プログラムの部分は、前後、左右への１歩を１の大きさにとると、
if(randn<0.25) x+=1;
else if(0.25<=randn<0.5) y+=1;
else if(0.5<=randn<0.75) x-=1;
else if(0.75<=randn<1.0) y-=1;
と書けば、ゼロから１の間の一様乱数でどこかへ１歩踏み出す、酔っぱらいを作ることができます。
その結果の例を２次元平面で図にしてみると、
なかなかおもしろい絵でしょう。この場合、酔っぱらいはなんと5000歩も歩くという酔いも醒めてしまうあるきっぷりでした。
統計力学では、このようなランダムな動きにもある種の法則性を見つけることができます。すなわち酔っぱらいの歩みは粒子の拡散現象の１種と考えることもできます。またフラクタルとも関係が深い。 粒子が全体としてN歩進んだとしましょう。これによって作られるみちすじの持つ性質を考えみましょう、しかし１粒子を追っかけた実験だけでは不十分であることは明らかですね。このような実験を多数回行いその平均（統計力学の用語ではアンサンブル）をとります。そしてN歩目とゼロ歩目（つまり初期位置）との距離R(N)を計算させてみてください。この値のゆらぎは 次の量に比例します。
<R(N)*R(N)> - <R(N)><R(N)>
これはよく使われる量ですね。分散（標準偏差は分散の平方根です）とも言われます。ちなみに<x>はｘの平均を表し、
<x>=(x1+x2+x3+.....+xN)/N
xという量の１回目の値がx1,２回目がx2,,,というとき、N回目まで測定すると、平均値<x>はこの式でかけるというわけです。
ゆらぎは物理の言葉では２次のモーメントです。（ちなみに１次のモーメントは平均のことですね）
統計力学を勉強してもらうと、R(N)はNの0.5乗に比例することを示すことができます。 これは、酔っぱらいが無限に歩いたなら２次元の全平面をくまなく歩きすべての点に足跡を付けることができることを意味します。これを計算物理的には、計算して比べてみましょう。歩く歩数を決めて、何人かの酔っぱらいさんに登場してもらい、到達距離を計算します。そして酔っぱらい人分の到達距離の平均値を得ることができます。
歩数を変えながら、到達距離の平均値の２乗をプロットしたものが、次の図です。縦軸を２乗にしたのはR(N)はNの0.5乗という表現よりは、R(N)の２乗はNに比例するというほうが、線形でわかりやすいためです。
もう少し乱数を使った話をしてみましょう。たとえばつぎの問題を考えてください。
いまあなたは銀行の窓口で待つ係りの人になってください。１時間に平均N人のお客さんが来るとしましょう（来客人数分布はポアソン分布するんですが）。またお客さんとの対応時間はゼロ（つまり書類を受け取るだけ）と仮定します。さてあなたの休んでいられる時間（あるお客さんとの対応が終わって次のお客さんが来るまでの時間）分布はどうなるでしょうか？
この問題は、バス停でバスを待つ人の待ち時間分布とも一致します。
ゼロから１までの一様乱数をN個生成します。さらにこれを小さい順に並べ替えます。これを英語ではsort （ソート）といいます。そして、この順番の隣同士の値の差が上の「待ち時間」となります。そして、これをヒストグラム化してやればよろしい。UNIX-likeに３つのプロセスに分けて実行してみました。（１）一様乱数をN個生成、（２）小さい順に並べ替え、（３）ヒストグラム。それぞれの実行ファイル名を (1) randgen (2)sort (3) dhistとし、
（１）randgen > rang としてrangなる一様乱数のファイルを作り出します。その結果は次のようになります。ここでは並べ替えによる順番が正しくなったことを示すために１０個で行います。
0.087708
0.591278
0.480469
0.402252
0.876648
0.585785
0.881714
0.137970
0.069885
0.648651
（２）sortは(1)の出力ファイル rangを読み込み、順序をそろえたファイルを作り出します。
それは次のようにして 実行されます sort < rang > rano
出力ファイルranoを見てみると、 more rano
0.069885
0.087708
0.137970
0.402252
0.480469
0.585785
0.591278
0.648651
0.876648
0.881714
このようにちゃんと小さい順に並んでいます。
（３）これを入力して、隣同士の差のヒストグラムをつくるプログラムを実行して、簡単なヒストグラムをつけて見ました。２０個のビンに分けました。このままでは、２０番目のビンの値は表示されていませんが、0.67ですから、１ビンあたり0.07位となります。注目すべきは最初のビンが一番多い？。これは待ち時間分布ですから、最小待ち時間で次のお客があるいはバスが来る！！！？？ってことを意味します。みなさんの経験と一致しますでしょうか？自然界の全くでたらめな現象の間の時間差ならば、その差が小さい確率がもっとも大きい。と言っています。
0 3
1 2
2 1
3 1
4 0
5 0
6 1
7 1
8 0
9 0
10 0
11 0
12 0
13 0
14 0
15 0
16 0
17 0
18 0
19 0
これで、この結論に到達するには１０個では数が不足していますね。こういうとき、統計不足といいます。そこで、１０００個の場合の結果を示します。
0 254
1 221
2 143
3 114
4 74
5 60
6 36
7 37
8 13
9 18
10 10
11 5
12 5
13 4
14 1
15 3
16 0
17 0
18 1
19 1
ただし、軸の倍率は変更しました。なぜなら０−１の間に１０個乱数を落とした、距離分布と、１０００個落とした距離分布は大きく異なるはずだからです。やっぱり最初の最小のビンが最大値を取っています。さてどんな分布をしている、あるいはどんな関数形をしているのでしょうか？これを gnuplotで図にしてみましょう。１０００個の場合の出力ファイル名を
randwait.d として、
gnuplot "randwait.d" with steps で描いてみましょう。
ふむふむですね。右下がりの関数形ですね。この関数形と特定するためには縦軸をlogスケールにすると物事が見えてきます。そのためには set logscale y とgnuplot内でやってください。再びgnuplot "randwait.d" with steps とすると次の絵が見えます。
やはり予想どおりexponetial（指数関数）ですね。
指数関数はみなさんのもっとも身近な関数ですね。自然界のたとえば物理学実験でよく用いられるミュー粒子の寿命を測定するときに現れる関数です。そこには、１個１個のミュー粒子ははランダムなタイミングで崩壊します。これを多数のミュー粒子を測定対象として、測定するとき、短い時間間隔dt内でに崩壊確率は一定値aとしましょう。 [時間tの間生き延び（崩壊することなく）]かつ[t〜t+dtの間に崩壊する確率] dP(t,a)は、二つの確率の積で表されます。すなわち時刻0からt までは崩壊しないので、(1-P(t,a))とt〜t+dtの間に崩壊する確率dt*aであり、
dP(t,a)=(1-P(t,a))*dt*a この微分方程式は変数分離されて、
dP/(1-P)=adt とかけ、これを積分すると、
log(1-P)=-at+c すなわち P=1-exp(-at) となります。これが上の図で見えているわけです。