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

4. f(x)の求積問題

4.1. 設計

最後も課題を確認しておきます。課題は「求積問題で適当な関数f(x)を選び、適当な区間で定積分し面積を計算するプログラムを作成し、ワークステーション上で実行せよ。試行回数の増加とともにその値はどのように変化するか? 入門的モンテカルロとあたりはずれモンテカルロの2つについて計算し、両者の結果を比較せよ。」です。

乱数の発生については、Buffonの針と同様に、先程作成したrandコマンドがそのまま使えます。あとは乱数を読み込んで定積分するプログラムを作ることになります。今回は入門的モンテカルロ法で定積分するプログラムをcrudeコマンド、あたりはずれモンテカルロ法で定積分するプログラムをhitormissコマンドと呼ぶことにします。テキストにある通り、入門的モンテカルロ法(crudeコマンド)では一つ乱数を発生させることが1ステップに相当します。また、あたりはずれモンテカルロ法(hitormissコマンド)では二つの乱数で二次元座標上の点を表現するので、二つ乱数を発生させることが1ステップに相当します。

4.1.1. crudeコマンド

crudeコマンドで乱数以外に必要になる人為的なパラメータは、被積分関数f(x)とその積分区間です。速度的な問題で被積分関数は直接プログラム中に書くことにします。積分区間に関しては、axbで積分する場合は、次のような動作イメージをもってプログラムを作成します。

foo[bar]% ./rand 10 | ./crude aの値 bの値 Enter

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

foo[bar]% ./rand 10 | ./crude 0 1 Enter
.
.
.
foo[bar]%

上の例では、randコマンドが10個乱数を発生させていますので、crudeコマンドは10ステップ実行されることになります。

4.1.2. hitormissコマンド

hitormissコマンドで乱数以外に必要になる人為的なパラメータは、被積分関数f(x)と矩形領域axb, 0 < ycを表す3つのパラメータa, b, cです。crudeコマンドと同様に速度的な問題で被積分関数は直接プログラム中に書くことにします。したがって、次のような動作イメージをもってプログラムを作成します。

foo[bar]% ./rand 10 | ./hitormiss aの値 bの値 cの値 Enter

例えば、それぞれa = 0b = 1c = 1と選んだ場合は次のようにするわけです。

foo[bar]% ./rand 10 | ./hitormiss 0 1 1 Enter
.
.
.
foo[bar]%

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

4.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]%

求積問題のプログラムを作成するための作業用ディレクトリを次のように作ります。

foo[bar]% mkdir integral Enter

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

foo[bar]% cd integral Enter

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

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

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

求積問題でも雛型を用意してあります。

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

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

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

4.4. ソースの編集

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

ファイルcrude.chitormiss.cの二つを編集して、プログラムを完成させます。

また被積分関数については、ファイルFunction.cを編集して定義します。

4.5. 実行ファイルの構築

求積問題でもMakefileを用意してあります。ファイルcrude.chitormiss.cをコンパイルするには、ktermで次のように入力します。

foo[bar]% make Enter

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

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

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

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

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

4.6.1. crudeコマンド

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

foo[bar]% ./rand 10 | ./crude 0 1 Enter

次のように10ステップ実行されます(値は異なるかもしれません)。

foo[bar]% ./rand 10 | ./crude 0 1
1 0.857868
2 0.921152
3 0.931160
4 0.909657
5 0.791601
6 0.823858
7 0.807871
8 0.828640
9 0.833127
10 0.849034
foo[bar]%

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

4.6.2. hitormissコマンド

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

foo[bar]% ./rand 10 | ./hitormiss 0 1 1 Enter

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

foo[bar]% ./rand 10 | ./hitormiss 0 1 1
1 1.000000
2 1.000000
3 1.000000
4 1.000000
5 1.000000
foo[bar]%

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

4.7. グラフの作成

試行回数nが増加するにしたがって、積分値が収束する様子をグラフにプロットしてみます。これについてもシェルスクリプトgraph.shが用意してあります。

4.7.1. crudeコマンド

試行回数が10000回で、それぞれa = 0, b = 1に選んだならば、ktermで次のように入力します。

foo[bar]% ./rand 10000 | ./crude 0 1 | sh graph.sh -0 crude Enter

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

foo[bar]% ./rand 10000 | ./crude 0 1 | tail Enter
9991 0.784151
9992 0.784167
9993 0.784184
9994 0.784204
9995 0.784208
9996 0.784210
9997 0.784216
9998 0.784232
9999 0.784196
10000 0.784164
foo[bar]%

入力するとxgraphのウィンドウが現れますので、収束している様子を確認したら、印刷して終了します。

4.7.2. hitormissコマンド

試行回数が10000回で、それぞれa = 0, b = 1, c = 1に選んだならば、ktermで次のように入力します。

foo[bar]% ./rand 20000 | ./hitormiss 0 1 1 | sh graph.sh -0 hit-or-miss Enter

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

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

foo[bar]% ./rand 20000 | ./hitormiss 0 1 1 | tail Enter
9991 0.783906
9992 0.783927
9993 0.783949
9994 0.783970
9995 0.783992
9996 0.784014
9997 0.784035
9998 0.784057
9999 0.783978
10000 0.783900
foo[bar]%