ランダムウォーク

20000回コインを投げて黒板に書いた得点がどのような動きをするか見てみましょう。

> x <- runif(20000)

> x <- 2*x

> x <- floor(x)

> x <- 2*x-1

> for (i in 2:20000) x[i] <- x[i]+x[i-1]

> plot(1:20000,x,type="l")

> abline(0,0,col="red")

5行目は「cumsum({vector})」を使うともっと簡潔に「cumsum(x)」とかけます。 この関数の使用例をちょっと挙げてみましょう。

> w <- 1:10

[1] 1 2 3 4 5 6 7 8 9 10

> cumsum(w)

[1] 1 3 6 10 15 21 28 36 45 55

元のベクトルがa=(a1,a2,a3,...,an)のとき cumsum(a)=(a1,a1+a2,a1+a2+a3,....,a1+a2+a3+...+an)になります。

最後の二行でこの得点の動きをグラフに表示します。実行する場合場合によって 傾向は異なりますがどちらかが勝っている時間が圧倒的に長いという方が 正しそうです。またこのように0からはじめて毎回の偶然の試行に応じた距離ずつ 動く過程をランダムウォークと呼びます。

ところでランダムウォークはこれからも使うことがありますが、毎回この 手続きを書くのは面倒なので一連の手続きをまとめて名前を付けて保存しましょう。 プログラムでもよいのですが関数を定義する方法があります。関数は「{関数の名前} <- function({引数}) {引数に対する処理}」の順で定義します。例えばf(x)=2x^2+3なら

> f <- function(x) 2*x^2+3

> f(5)

[1] 53

使うときは上の二行目のように「{関数の名前}({引数})」です。 ランダムウォークの関数を作ってしまいましょう。

> rwalk <- function(x) {

+ ran <- runif(x)

+ ran <- 2*ran

+ ran <- floor(ran)

+ ran <- 2*ran-1

+ cumsum(ran)

+ }

ではランダムウォークを一つだけ考えるのではなく、いくつも考えて 傾向を探りましょう。フェラーの本ではコイン投げを20000回していますが 計算時間の都合で2000回として、このランダムウォークを1000個分調べましょう。 2000回のコイン投げを1000個ですから2000000回コインを投げることになります Pentium2で数分はかかるでしょう。1000個の各ランダムウォークにおける逆転回数 とピーターのリード時間を調べましょう。その記述場所として1000次元のベクトルy,zを 用意します。

> y <- 1:1000

> z <- y

> for (j in 1:1000){ x <- rwalk(2000)

+ y[j] <- length((1:2000)[x==0])

+ z[j] <- length((1:2000)[x>0])

+ }

ここでy[j]はj個目のランダムウォークが0の位置にいる時間を調べたものであり、 z[j]はj個目のランダムウォークが正の位置つまりピーターのリード領域にいる 時間を調べたものです。リードが逆転するには0を通らなければならず、そこから 元と同じリードに戻ることもあるのでyのだいたい二分の一がリードが逆転する回数です。 リードが逆転する回数を正確に調べるのは難しくありませんが計算時間が 長くなるのでここではやめます。ではyとzのヒストグラムを作ってみましょう。

hist(y)

hist(y,nclass=20)

yからリードが変わらないことがいちばんありそうなことで、変わる回数が低いほうが 可能性が高いことがわかります。

hist(z)

hist(y,nclass=20)

zからはどちらかがずっとリードしている方が伯仲しているより ありえそうなことがわかります。

ホームヘ