山内の授業のページ/07春/研究室ゼミ?

数値計算のプログラムで遊ぼう2

課題2 もう1度、円の面積を求めて円周率を計算してみよう

第2回は、ぜんぜん違う方法で円の面積を計算し、円周率πの値を求めてみます。

確率の復習

確率論を習ったとき、さいころを振って1の目の出る確率が1/6、2の目の出る確率も1/6なんて話を聞いてますよね。これを使おうと思います。

たとえば、区間0から1までの間に、均等に分布するような確率分布があるとします。さいころなら、目の種類が0から1までの数すべて、という状況です。

区間0〜1の間のある数 x が出る確率は?といいたいですが、これはいけません。実数の稠密性の話があるので、1つの実数xの出る確率は限りなく0に近くなってしまいます。 だから、区間で考えましょう。 区間0〜0.5の数の出る確率は? これなら0.5と答えられます。

ところで、2つの実数の対(ベクトルと呼んでもよいですが) (x, y) があるとしましょう。この数xとyがいずれも0〜1の区間の実数で、ランダムに値を取り、確率分布が一様であるとします。 つまり、0〜1の値の出るさいころを2個振っていると思ってください。この(x, y)が、2次元の座標を表すと読むことにします。 
ということは?

そう、点(x, y)は、0<=x<=1, 0<=y<=1 が占める正方形の上に、一様にばらまかれることになります。 さいころを振るたびに、この正方形の中の1点がランダムに書かれるわけです。図の上で「×」で書いたのが、さいころを振ったときの点です。

1/4円の面積

さて、正方形を区切ります。前と同じように、半径1の1/4円にしましょう。 前と同じように図形を書くことにします。ちょうど(0, 0)〜(1, 1)の正方形の中に入りますよね。
random1x1.gif

さて、上の節で、区間0〜1のさいころでランダムに発生した数 x が区間0〜0.5に落ちる確率は、0.5だと言いました。同じように、区間(0<=x<=1, 0<=y<=1)のさいころでランダムに発生したベクトルv=(x, y)が、この1/4円の内側に落ちる確率は、どうなるでしょうか?一次元の場合は全体に対する長さの比(全体が長さ1に対して、区間0〜0.5の長さは0.5、この比率が0.5/1=0.5)でした。二次元だと面積の比になります。つまり、1/4円の内側に落ちる確率は、全体の面積(=1)に対する1/4円の面積(=π/4)になるはずです。

ここでは逆に、たくさんの回数さいころを振って、ランダムに発生する点が1/4円の内側に落ちるか外側に落ちるかを数えてみましょう。たとえば10000回さいころを振って、もしそのさいころが公平であれば、1/4円の内側に落ちる回数は、(π/4)*10000 のはず(正確には回数を増やすとπ/4に近づくはず)ですね。 だから、(内側に落ちた回数)/10000*4 を計算すれば、πになるはずではないか、という話です。
random1x1_2.gif

計算機上のさいころ

計算機上でさいころを実現するには、ちょっとした仕掛けが要ります。これは別にプログラミング技術の話ではないので、サンプルを見てしまうことにします。 (乱数発生の原理の話は割愛します。キーワード「擬似乱数」や「混合合同法」で検索してみると出てくると思います。)

C言語では、乱数を1つ発生させる関数を呼び出すことができます。この関数は、1回呼び出すと1つ乱数を返します。もう1度呼び出すと、異なる乱数を返します。また呼び出すと、また異なる乱数を返します。繰り返して呼び出すことによって、乱数の列を発生することができます。 乱数といっても、本当に「ランダム」=次の数が前から予測できない=というのはプログラムである限り実現ができないので、見かけ上ランダム〜擬似乱数〜とみなされる列です。このあたりの議論はいろいろできますが、今はちょっと置いておきます。

乱数発生は、関数randを使ってみます。このあたりはシステムによって違っていて、関数randomを勧めているケースも多いです。このあたりの情報は、オンラインマニュアルページ man rand や man random として、よく読んでみてください。RedHat 7.2 Linuxでは、次のように書いてありました。


RAND(3) Linux Programmer's Manual RAND(3)

名前 rand, srand - 乱数を生成する関数

書式 #include <stdlib.h>
 int rand(void);
 void srand(unsigned int seed);

説明
rand() 関数は 0 から RAND_MAX の間の疑似乱数整数を返す。

srand()関数は、 rand()関数で作られる疑似乱数整数系列の新しい種として、その引数の値を使用する。

これらの関数を使用して作られた疑似乱数系列は、同じ値を引数に使用してsrand()を呼ぶことに
よって、再現することが可能である。

種の値が与えられない場合には、rand()関数は自動的に 1 を 種とする。

返り値
rand() 関数は 0 と RAND_MAX の間の数を返す。 srand()関数はどのような値も返さない。

注意
rand()とsrand()の Linux C Library 版は random()とsrandom()のふたつの関数と同じ乱数生成の
アルゴリズムを使用している。
そのため、下位のビットは上位のビットと同じくらいにランダムである。しかし旧版の rand() の
実装は、下位のビットが上位のビットほどランダムになってはいない。

Numerical Recipes in C: The Art of Scientific Computing
(William H. Press, Brian P. Flannery, Saul A. Teukolsky,
William T. Vetterling; New York: Cambridge University
Press, 1992 (2nd ed., p. 277)) では、次のようなコメントがなされている。

「1 と 10 の間の乱数を生成したいのであれば、常に
j=1+(int) (10.0*rand()/(RAND_MAX+1.0));
のように上位ビットを用いること。
j=1+(rand() % 10);
のような、下位のビットを用いるような方法は一切止めること。」
 
乱数の生成は複雑な話題である。前述の Numerical Recipes in C でも実用的な乱数生成を論点と
した優れた議論を第7章(乱数)で展開している。

より 理論的な議論については Donald E. Knuth の The Art of Computer Programming, 
volume 2 (Seminumerical Algorithms), 2nd ed.; Reading, Massachusetts: AddisonWesley
Publishing Company, 1981 の第 3 章(乱数)を見よ。これは、たくさんの実用的な話題についても
深く網羅している。

準拠
SVID 3, BSD 4.3, ISO 9899

関連項目
random(3), srandom(3), initstate(3), setstate(3)

GNU 18 May 1995 RAND(3)
------------------------------------------------------------------------

では、早速プログラムを書いてみましょう。

まず、0〜1の間の一様な乱数を発生させて見ましょう。 乱数が一様であるかどうかを確認しておきましょう。たとえば、N個の乱数を発生させ、0〜0.1、0.1〜0.2、0.2〜0.3...のように、0.1ずつの区間にいくつ属するかを、数えてみましょう。N/10個ぐらいずつになっていればよいわけですが、実際やってみるとかなりばらついていることがわかるでしょう。

いよいよ、面積を求めてみます。区間0〜1に分布する乱数を2回発生し(さいころを2回振ることに相当する)、一方をx、他方をyとしましょう。この2つの値(x, y)が点Pの座標になります。 

次に、この点Pが、1/4円の内側に属するのか、外側に属するのかを判定します。これは、x^2+y^2が1より大きいかどうかで判定すればいいでしょう。

全部でN個の点 (=2N個の乱数) を発生します。 このN個の点のうち、いくつが内側で、いくつが外側だったかを数えておきます。内側だった個数がW個だったとすると、W/Nは、1×1の正方形の中に落ちるN個の点のうちで、1/4円の内側に落ちた点の個数の比率を示しています。 もしN個の点が一様に分布しているとすれば、W/Nは四角に対する円の面積の比率になっているはずです。

この原理で、1/4円の面積を推定します。 このようにランダムなプロセスを使って値を推定する方法は、いわゆる「モンテカルロ法」と呼ばれるシミュレーションの1つです。

値の精度は、点の総数Nを増やしても、あまり上がらないと思います。 おそらく乱数の一様性に問題があるからだと思います。


添付ファイル: filerandom1x1_2.gif 400件 [詳細] filerandom1x1.gif 430件 [詳細]

トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2007-02-25 (日) 19:41:23 (3921d)