深夜の盛り場近く, 前後不覚に酔っぱらったおじさん(ということにしておこう)が,あっちこっちと行方を定めることができずに ふらついています。 さて, この人が1000歩あるくと一体どこに行き着いているのでしょうか? ランダムウォーク(Random Walk)というのはこのようなモデルで, そのまま酔歩とも呼ばれています。
最初に線路の上をふらふらと歩く酔っ払いの動きを眺めてみることにしよう。 もちろんこんなことは鉄道運送法違反だかでつかまってしまうから,実際にやってはいけない(って,法律がどうという問題じゃないかもしれないが)。
1次元の完全なランダムウォークでは, 1歩の幅を 1 とすると,$+1$ と $-1$ が等確率でランダムに現れる。 全部で $n$ 歩だけ動いたときに正の方向に進んだ回数を $x$ とすると, $x$ の分布はよく知られた二項分布の形になる。 なお,$n$ が十分に大きければ(30とか50程度以上), 二項分布は正規分布として近似できることにも注意しておこう。 このあたりの議論は統計学の教科書を参考にできる。
$ f(x) = _nC_x (p)^{x}(1-p)^{n-x} = _nC_x (1/2)^n $
前進する確率は後退する確率と等しいから $p = 1/2$ である。
確率論でよく知られているように,$x$ の
期待値 $E[x]$ と分散 $V[x]$ は次のようになる。
$ E[x] = np = n/2 $
$ V[x] = n p (1-p) = n/4$
それではランダムウォークによる前進の距離について考えてみよう。 $n$ 歩のうち,前進の歩数が $x$ であったとすると, 後退の歩数は $n-x$ つまり$-(n-x)$ だけ正の方向に進むわけだから, 差し引きの前進距離 $z$ は次のようになる。
$ z = x - (n-x) = 2x - n $したがって,$z$ の期待値と分散は次のようになる。
$ E[z] = E[x - (n-x)] = E[2x -n] = E[2x] - n = 2 \times n/2 - n = 0$
$ V[z] = V[2x-n] = 2^2 V[x] = 4 \times n/4 = n $
上の式が意味していることを考えてみよう。
まず,酔っ払いの位置の期待値は 0 であるから,平均すると出発点の ところにいることになる。ただし,これは何百人もの酔っ払いが 線路上をさまよったとして, 同じ歩数だけ歩いた後の位置を平均するともとの位置になるということを意味しているのであって, 個々の酔っ払いの位置は広がりをもつことになる。
では,位置の広がりの程度はどうなるかというと,標準偏差は分散の 正の平方根だから,次のようになる。
$ \sqrt{V[z]} = \sqrt{n] $
つまり,酔っ払いが $n$ だけ歩いた後の位置の広がりは $\sqrt{n}$ の標準偏差をもつというわけだ。
上の結論をシミュレーションで確かめてみよう。 繰り返し計算を行って,結果を GNUPLOT で表示するだけでよいので, Ruby でソースを書くことにする。
このプログラムは次のような動作を行うものとする。
nstep = 100 position = 0 nstep.times{ if rand(2) == 1 then position += 1 else position -= 1 end }
位置をどうやって覚えるかを考えよう。最終的には, ある位置に止まった回数はどれだけなのかという統計を取りたいのである。 すると, たとえば 12 のところで止まったとすると,「12という場所」にあるカウンタを 1つだけ増加させるといったことが必要だ。 そのような仕掛けとして有力なのは何だろうか?もちろん配列だ。
たとえば 0という位置に 3 回,1という位置に 5 回,2という位置に 6 回 止まったということは,中身が [3,5,6] という配列で表される(選挙の開票で,候補者ごとの 得票数が並んだボードを想像するとよい)。
そこで位置を記録するために freq という配列を用意することにしよう。これは frequency(頻度,度数分布の「度数」の意味)からとっている。上の ソースで示された1回の試行の 前でこの配列を宣言しておけばよいだろう。試行の後で,その位置の番地に 1を加えればよい。
#!/usr/local/bin/ruby freq = [] nstep = 100 position = 0 nstep.times{ if rand(2) == 1 then position += 1 else position -= 1 end } freq[position] += 1
さて,上のソースを走らせてみよう。すると,次のような エラーメッセージが出力されてしまう。
./test.rb:12: undefined method '+' for nil:NilClass (NoMethodError)
(nil に対して '+' というメソッドは定義されていない)
ソースを見ると,このエラーは freq の position 番地に1を加えようとしたところで 生じている。つまり,この番地には何も代入されていないわけで,それに 1を足すのはできないと言っているのだ。
そもそも, freq はソースの2行目の freq = [] で宣言されているのだが, このときの配列のサイズは 0 である。Java などでは 配列の宣言のときにサイズ(配列要素の数)も同時に指定することになっている わけなので,Ruby でも同じように配列のサイズを同時に宣言し,かつ 初期値として 0 を与えることにしたい。 それには次のようにする。
ref = Array.new(size,0)ここで size は配列のサイズである。
もうひとつ考えておかなければならない問題がある。それは酔っ払いの停止した位置は 正と負の値をとりうるということだ。たとえば10歩歩いたとすると,止まった場所は -10 から 10 までの 21 通りになる。つまり nstep の歩数に対しては 2 * nstep + 1 だけの要素数を確保する必要がある。
以上を考え合わせて, 次の行を freq = [] の代わりに置いてみよう。
freq = Array.new(2*nstep+1,0)
ソースの末尾には,結果を見るために次の2行を追加する。
puts position p freq
さて,ここまでのソースを走らせてみると, position で表された位置と p frep で表される配列の要素の位置が ずれていることに気がつくはずだ。たとえば position が 0 のときには, freq の真ん中の要素が 1 にならなければならないのに, 左端の要素が 1 になってしまう。そこで,真ん中の要素のところに ゼロが来るようにしないといけない。具体的には,
freq[position] += 1
を
freq[position+nstep] += 1
としてやればよい。
ここまでで,酔っ払いが出発点から歩き出して,ある歩数だけ歩いたときに 行き着く先が記録される仕組みができた。 そのソースは下からダウンロードできる。
randomwalk1.rb(酔歩1パスのシミュレーション)
以上で1回の試行によって,酔っ払いがどの位置に止まるかを記録する仕掛けは できた。後は何度も試行しては,その記録を freq に書き込むように,ループを形成すればよい。 具体的には ntrial を試行回数として,次のようなループになる。その結果,ここでは times による二重ループが作られることになる。
nstep = 10 # この値は適宜変更 ntrial = 1000 # 積算回数(最初は少なめに。本番ではなるべく大きく) freq = Array.new(2*nstep+1,0) ntrial.times{ 酔歩 1 パスの処理 ... ... ... } 結果をうまく出力できるように工夫。
最後は結果の出力だ。上のプログラムの最後に
p freqとして,出力を見ると,たとえば次のようになっているはずだ。
[1, 0, 12, 0, 57, 0, 114, 0, 203, 0, 236, 0, 195, 0, 141, 0, 31, 0, 8, 0, 2]
この21個の値は,その位置にたどり着いた回数を示している。たとえば $-10$ の位置にたどりついた回数は 1 である。同様にして,$-8$ の位置には12回たどり着いている。 つまり,この結果を表にすると次のようになる。
-10 | -9 | -8 | -7 | -6 | -5 | -4 | -3 | -2 | -1 | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 |
1 | 0 | 12 | 0 | 57 | 0 | 114 | 0 | 203 | 0 | 236 | 0 | 195 | 0 | 141 | 0 | 31 | 0 | 8 | 0 | 2 |
この結果を GNUPLOT で表示させるためには,次のようなデータが出力されればよい
-10 1 -9 0 -8 12 -7 0 -6 57 ...
実際には, freq の添字は 0 からスタートするわけだから,上のように $-10$ からスタートさせるには, 添字の値から nstep を引いてやる必要がある。
for i in 0 ... freq.size puts "#{i-nstep} #{freq[i]}" end
ここまでの計算で得られているのは,ある位置にたどり着いた回数(頻度)だった。 しかしそれだと,試行回数 ntrial を増やすと値が大きくなってしまう。したがって,回数を試行回数で割った値を 出力したほうがよい。このとき,整数同士で割ってしまうとすべてゼロになってしまうから, 割り算のどちらかの値を .to_f で浮動小数点に変換しておかないといけない。
ここまでのすべてを考慮に入れたプログラムを示しておこう。このプログラムでは, コマンドラインから歩数を入力することになっている。 randWalkStat.rb
完成したプログラムを使って,ランダムウォークの計算機実験をやってみよう。次のように 歩数を変えながら実行して,その結果をファイルにリダイレクトし, GNUPLOT で表示させてみよう。かわいそうな酔っ払いの歩みについていろいろな ことがわかるはずだ。
./randWalkStat.rb 50 > rw0050.dat ./randWalkStat.rb 100 > rw0100.dat ./randWalkStat.rb 200 > rw0200.dat ./randWalkStat.rb 500 > rw0500.dat ./randWalkStat.rb 1000 > rw1000.dat