モンテカルロシミュレーション

3. Buffonの針

3.1. 設計

また課題を確認しておきます。課題は「Buffonの針の問題から円周率を計算するプログラムを作成し、ワークステーション上で実行せよ。試行回数の増加とともにその値はどのように変化するか。」です。

乱数の発生については、先程作成したrandコマンドがそのまま使えます。あとは乱数を読み込んで円周率を出力するプログラムを作ることになります。今回はそれをbuffonコマンドと呼ぶことにします。Buffonの針では、テキストにある通り、二つ乱数を発生させることが一回針を投げる行為に対応しています(このシミュレーションでは一回針を投げることが1ステップに相当します)。

乱数以外に必要になる人為的なパラメータは、テキストにあるように、aLになります(小文字のエルは字体によっては読みにくいので、以下大文字で表記します)。したがって次のような動作イメージをもってプログラムを作成します。

foo[bar]% ./rand 10 | ./buffon aの値 Lの値 Enter

例えば、それぞれa = 1.0L = 0.5と選んだ場合は次のようにするわけです。

foo[bar]% ./rand 10 | ./buffon 1.0 0.5 Enter
.
.
.
foo[bar]%

上の例では、randコマンドが10個乱数を発生させていますが、buffonコマンドは1ステップで2つの乱数をとるので、5ステップだけ実行されるわけです。

方針は以上です。あとは少し細かい点についてテキストを補足しておきます。まず、針の中心の座標yのとり方ですが、「-aya」ととった方が簡単です。そのかわり、針が平行線と交わっているかどうかはy1y2の符号が異符号か同符号かで判断します(テキストでは「y1 - 2a」と「y2 - 2a」の符号を比較するように書いてあります)。また、θについては、「0 ≦ θ < 2π」ととってください。

3.2. 作業用ディレクトリの作成

次のようにpwdコマンドを入力して、現在自分がホームディレクトリにいることを確認してください。

foo[bar]% pwd Enter
/home/appi/foo
foo[bar]%

もし違っていたらcdコマンドでホームディレクトリに移動します。

foo[bar]% cd Enter
foo[bar]% pwd Enter(確認のため)
/home/appi/foo
foo[bar]%

乱数の度数分布のときと同じように、Buffonの針のプログラムを作成するための作業用ディレクトリを次のように作ります。

foo[bar]% mkdir buffon Enter

次にこのディレクトリに移ります(繰り返しになりますが、このときもtcshのファイル名補完機能を用いればタイプ量を減らすことができます)。

foo[bar]% cd buffon Enter

最後に念のためpwdコマンドでカレントディレクトリを確認します。

foo[bar]% pwd Enter
/home/appi/foo/buffon
foo[bar]%

3.3. ソースの雛型のコピー

乱数の度数分布のときと同様に、Buffonの針でも雛型を用意してあります。

以前と同様に、次のようにcpコマンドを用いて、ファイルを自分のホームディレクトリにコピーします。

foo[bar]% cp ~yoshiko/skel/buffon/* . Enter
foo[bar]%

ファイルをコピーしたら、念のためlsコマンドで確認してみてください。

3.4. ソースの編集

乱数の度数分布で起動したMuleが、まだ動いている場合はそのまま使用しても構いません。終了してしまった場合はもう一度Muleを起動してください。

buffon.cを編集して、プログラムを完成させます。

3.5. 実行ファイルの構築

Buffonの針でもMakefileを用意してあるので、ファイルbuffon.cをコンパイルするには、ktermで次のように入力します。

foo[bar]% make Enter

プログラムに間違いがなければ次のように出力され、自動的に実行ファイルbuffonが構築されます。

foo[bar]% make
gcc -Wall -O2 -pipe -D__USE_FIXED_PROTOTYPES__ -c buffon.c
gcc -Wall -O2 -pipe -D__USE_FIXED_PROTOTYPES__ -c Trial.c
gcc -o buffon buffon.o Trial.o -lm
ln -s ../rand-distribution/rand rand
ln -s ../rand-distribution/graph.sh graph.sh
foo[bar]%

エラーや警告が出た場合は、内容を調べて修正し、もう一度makeしてみます。

正常に構築できたようであれば、次に進みます。

3.6. プログラムの実行テスト

実行する前に、パラメータaLを決定しなければなりません。例えば、それぞれa = 1.0L = 0.5に選んだ場合は、ktermで次のように入力します。

foo[bar]% ./rand 10 | ./buffon 1.0 0.5 Enter

randコマンドは10個の乱数を発生させますが、buffonコマンドは1ステップあたりに2つの乱数をとるので、次のように5ステップだけ実行されます。

foo[bar]% ./rand 10 | ./buffon 1.0 0.5
1 1.600000
2 3.200000
3 4.800000
4 3.200000
5 2.666667
foo[bar]%

正常に動作しなかったり、値がめちゃくちゃだったりした場合は、Muleで修正して再びmakeします。

3.7. グラフの作成

試行回数nが増加するにしたがって、どのように円周率が収束するかグラフにプロットしてみます。これについてもシェルスクリプトgraph.shが用意してあります。試行回数が10000回で、それぞれa = 1.0L = 0.5に選んだならば、ktermで次のように入力します。

foo[bar]% ./rand 20000 | ./buffon 1.0 0.5 | sh graph.sh -0 buffon Enter

randコマンドに与える引数は、試行回数の2倍でなければならないことに注意してください。入力するとxgraphのウィンドウが現れますので、収束している様子を確認したら、印刷して終了します。

また、試行回数10000回目での円周率を数値で知りたい場合は、次のように入力することで知ることができます。

foo[bar]% ./rand 20000 | ./buffon 1.0 0.5 | tail Enter
9991 3.124018
9992 3.123720
9993 3.123423
9994 3.123125
9995 3.122828
9996 3.123140
9997 3.123452
9998 3.123765
9999 3.123467
10000 3.123780
foo[bar]%