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

5. 酔歩問題(時間が余った場合)

5.1. 設計

時間が余った場合は、酔歩問題を考えてみましょう。課題は「(実験書の)3.3の酔歩問題のシミュレーションプログラムを作成し、ワークステーション上で実行せよ。N = 100として、時刻n = 10, 20, ..., 100の時の酔っぱらいの位置xへの到達確率をシミュレーションから求めよ。位置xへの到達頻度分布をヒストグラムとして整理すること。また、各ステップでの変位の平均、分散も合わせて出力し、時間の関数としてプロットせよ。またNについていくつかの場合について試みよ。」です。

この酔歩問題のプログラムをwalkコマンドと呼ぶことにしましょう。実験書の手順を参考にすれば、次のような疑似コードが書けます。(今までと異なり、乱数の発生については処理速度向上のため、randコマンドの出力を読み込むのではなく、関数rand()を直接利用します。)

    for (N回繰り返す)
    {
        x に 0 を代入する;

        for (n回繰り返す)
        {
            0以上1未満の乱数を発生させて r に代入する;

            if (r < 0.5)
            {
                x から1を引く;
            }
            else
            {
                x に1を加える;
            }
        }
        
        /* x はこのとき `n歩後の男の位置' である。*/
        x を出力する; 
    }

人為的なパラメータは、試行回数Nと、一回の試行で動く歩数nです。これらは、walkコマンドへの引数として与えることにします。残るは、位置xに到達した回数Nxを数える作業です。この作業は、プログラムを作成しなくても、後で説明するようにsortコマンドとuniqコマンドで処理することができます。

また、変位の平均と分散を時間の関数としてプロットする作業は、先程作成したaverageコマンドとvarianceコマンドを使用します。

5.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 random-walk Enter

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

foo[bar]% cd random-walk Enter

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

foo[bar]% pwd Enter
/home/appi/foo/random-walk
foo[bar]%

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

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

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

foo[bar]% cp ~yoshiko/skel/random-walk/* . Enter
foo[bar]%

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

5.4. ソースの編集

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

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

5.5. 実行ファイルの構築

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

foo[bar]% make Enter

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

foo[bar]% make
gcc -Wall -O2 -pipe -D__USE_FIXED_PROTOTYPES__ -c walk.c
gcc -o walk walk.o -lm
ln -s ../rand-distribution/average .
ln -s ../rand-distribution/variance .
foo[bar]%

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

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

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

5.6.1. walkコマンド

実行する前に、パラメータN, nを決定しなければなりません。とりあえずそれぞれN = 10n = 4に選んでテストすることにします。ktermで次のように入力します。

foo[bar]% ./walk 10 4 Enter

次のように、酔歩4歩のシミュレーションが10回試行されます(値は異なるかもしれません)。

foo[bar]% ./walk 10 4
0
0
-4
0
4
0
0
4
4
-2
foo[bar]%

4歩の酔歩の場合、最終的な位置は-4, -2, 0, 2, 4のいずれかです。正常に動作しなかったり、値がめちゃくちゃだったりした場合は、Muleで修正して再びmakeします。

5.6.2. sortコマンド

今度は、最終的な位置-4, -2, 0, 2, 4が、それぞれ何回ずつ現れるのか数えなければなりません。最初に数えやすくするために、walkコマンドの出力をsortコマンドで並び換えます。次のように入力してみます。

foo[bar]% ./walk 10 4 | sort -n Enter

ここで、sortコマンドへの引数-nは、入力文字列を数とみなして比較するように指示するものです。結果は、次のように表示されます。

foo[bar]% ./walk 10 4 | sort -n
-4
-2
0
0
0
0
0
4
4
4
foo[bar]%

5.6.2. uniqコマンド

さらに、それぞれの値の出現回数を数えるには、uniqコマンドを使用します。uniqコマンドの本来の用途は、「連続して重複する行をひとつの行にまとめる」ことです。したがって、次のように入力すると、

foo[bar]% ./walk 10 4 | sort -n | uniq Enter

次のように重複した行がひとつになります。

foo[bar]% ./walk 10 4 | sort -n | uniq
-4
-2
0
4
foo[bar]%

しかし、これでは役に立ちません。uniqコマンドには、引数-cを与えることによって、「行が入力中で続けて出現した回数」を各行の先頭に(空白をひとつあけて)表示する機能があります。そこで、次のように入力します。

foo[bar]% ./walk 10 4 | sort -n | uniq -c Enter

次のように出現回数が行の先頭に付加されます。

foo[bar]% ./walk 10 4 | sort -n | uniq -c
1 -4
1 -2
5 0
3 4
foo[bar]%

最終的な位置-4, -2, 0, 2, 4が、それぞれ1, 1, 5, 0, 3回出現したことがわかります。これをグラフに表示すれば、ヒストグラムが作成できます。

5.7. 到達頻度分布のグラフの作成

結局、例えばN = 100として、次のように実行することで、時刻n = 10, 20, ..., 100の時の位置xへの到達頻度分布を得ることができます。

foo[bar]% ./walk 100 nの値 | sort -n | uniq -c Enter

しかし、これについてもシェルスクリプトhistogram.cshを用意してありますので、次のように入力することで、ヒストグラムを表示できます。

foo[bar]% ./histogram.csh Nの値 nの値 [ nの値 ... ] Enter

例として、N = 100n = 10, 50, 100の時のヒストグラムを作成するには、次のように入力します。

foo[bar]% ./histogram.csh 100 10 50 100 Enter

実際には自分でパラメータを適当に選び、ヒストグラムを作成してください。ヒストグラムはxgraphのウィンドウに現れますので、印刷して終了します。

5.8. 変位の平均と分散のグラフの作成

もう一度話を戻しますが、「酔歩4歩のシミュレーション」を10回試行して、得られた10個の(4歩後の)変位の平均を求めてみましょう。次のように入力してみます。

foo[bar]% ./walk 10 4 | ./average Enter

すると次のように表示されます(値は異なるかもしれません)。

foo[bar]% ./walk 10 4 | ./average
1 0.000000
2 0.000000
3 -1.333333
4 -1.000000
5 0.000000
6 0.000000
7 0.000000
8 0.500000
9 0.888889
10 0.600000
foo[bar]%

必要な値は「10個の変位の平均」なので、次のように入力して最終行だけを取り出します。

foo[bar]% ./walk 10 4 | ./average | tail -1 Enter

すると次のように表示されます。

foo[bar]% ./walk 10 4 | ./average | tail -1
10 0.600000
foo[bar]%

あとは、あるNについて、いくつかのnで上記の操作を繰り返し、得られた平均をnについてプロットすれば、必要なグラフが手に入ります。もちろん、分散に関しても同様です(averageコマンドをvarianceコマンドに入れ換えることで、分散も求めることができます)。

しかし、手作業で同じようなことを繰り返すのも面倒です。今回もシェルスクリプトgraph.cshを用意してありますので、それを用いてグラフを作成します。次のように入力することで、平均と分散のグラフを出力できます。

foo[bar]% ./graph.csh Nの値 nの値 [ nの値 ... ] Enter

例として、N = 100n = 10, 20, 30, 40, 50の時の平均と分散のグラフを作成するには、次のように入力します。

foo[bar]% ./graph.csh 100 10 20 30 40 50 Enter

実際には自分でパラメータを適当に選び、グラフを作成してください。グラフはxgraphのウィンドウに現れますので、印刷して終了します。