ケィオスの時系列解析メモランダム

時系列解析,生体情報学,数学・物理などの解説です.

【Rでシミュレーション】ブラウン運動のランダムウォークモデル

講義用の参考資料です.今回は,ランダムウォークのサンプルパスを描きます.簡単化のためにx軸上の動きのみを考えます.

ランダムウォークのサンプルパス (p = 0.5, q=0.5)

 ランダムウォークは,ブラウン運動の確率モデルです.物理学におけるブラウン運動の重要性を知っていますか.今では奇跡の年と呼ばれる1905年に,26歳のアインシュタインは3本のノーベル賞級の論文を発表しました.その中の1本がブラウン運動(分子の大きさの決め方)についてのものです.その他の2本は,光電効果ノーベル賞受賞)と特殊相対論( E=mc^2の前段)です.特殊相対論の論文を博士論文として大学に提出したところ,理解されず受理されなかったので,ブラウン運動のネタの方を提出したそうです.

 アインシュタインの論文では,ブラウン運動を観測することで,分子の実在を示す理論が導かれています.当時はまだ分子の実在に否定的な意見がありました.なぜなら,「分子が存在し,その運動が古典力学で記述できる」と仮定してみても,マクロな熱現象を説明できないからです.背理法に基づけば,最初の仮定が間違っています.つまり,「分子」です.しかし,本当に間違っていたのは「分子」の部分ではなく,「古典力学」でした.当時はまだ量子力学が確立していなかったのです.

 後に,ペランが実験的にこの理論を検証して,分子の実在を示したことで1926年にノーベル賞をもらっています.興味があれば,米沢 富美子先生の本『ブラウン運動』(共立出版)を読んでみてください.ということで,ランダムウォークは物理学では重要なモデルです.さらに,時系列のフラクタル解析でも重要な基礎となりますので,ランダムウォークのイメージをもつために,各自で数値実験をしてみてください.

ランダムウォークの100本のサンプルパス(時間ステップ100まで)

 以下のRスクリプトで数値実験を実行できます.確率pで+1,確率q=1-pで-1になります.式で書くときは,確率を {\rm Pr}で表すことにして,

  {\rm Pr}( \Delta X_i = +1)=p
  {\rm Pr}( \Delta X_i = -1)=q = 1-p

であり,ランダムウォークの軌道は,

  \displaystyle X_i = \sum_{j=1}^{i} \Delta X_j

です.今回は確率変数を大文字Xで書きました.実現値を表すときは小文字xを使います.

グラフの上段は, \Delta X_iの実現値 \Delta x_iで,下段はその累積(積分 x_iです.pの値を0.5から変えてみたり,Nを大きい値にしたりして,ランダムウォークの振る舞いを観察してみてください.

# +1が出る確率
p <- 0.5
# 総ステップ数
N <- 100
###############
i <- 1:N
q <- 1-p
# -1,+1の乱数の生成
tmp <- runif(N)
x <- (tmp > p)*2 - 1
###############
par(mfrow=c(2,1),mar=c(5,5,2,2))
plot(i,x,"p",pch=16,col=2,xaxs="i",las=1,xlab=expression(italic(i)),ylab=expression(paste(Delta, italic(x)[italic(i)])),
  lwd=2,xlim=c(0,N+1),main=paste("p = ",p,"; q = ",q,sep=""))
arrows(i,numeric(N),i,x,col=2,lwd=2,length=0)
abline(h=0,lty=2,col=8)
y <- cumsum(x)
plot(c(0,i),c(0,y),"o",pch=16,col=2,xaxs="i",las=1,xlab=expression(italic(i)),ylab=expression(italic(x)[italic(i)]),
  lwd=2,xlim=c(0,N+1),main=expression(paste(italic(x)[italic(i)], "=",Sigma, Delta, italic(x)[italic(i)])))
abline(h=0,lty=2,col=8)

このスクリプトを実行すると,下のような図が描かれます.

Rスクリプトの実行例

デモ用の動画も作っておきました.

www.youtube.com