Pythonによるデータ処理17 ～ 1Dランダムウォーク

今回からPythonによるランダムウォークのシミュレーションについて取り上げる。今回は、最も基本的な1次元(1D)のランダムウォークの問題で原点に戻る確率について考える。（理論の方は以下参照）

zakii.la.coocan.jp

1. 1Dランダムウォーク
2. 中心にいる確率の計算値
3. 中心にいる確率の理論値
4. まとめ

1. 1Dランダムウォーク

1Dのランダムウォークをシミュレーションするためには、乱数行列を用いるのが手っ取り早い。以下のような+1と-1の乱数行列を用い、列(横)方向がステップ数、行方向は試行回数としてとらえる。

$\begin{pmatrix} 1 & -1 & \dots & 1 \\ -1 & -1 & \dots & 1 \\ \vdots & \vdots & \ddots & \vdots \\ 1 & 1 & \dots & -1 \end{pmatrix}$

ここで列方向に和を取ったときの値が中心0からのずれになる。

Pythonの場合1と-1をランダムに発生させる関数はないので、0から1まで一様乱数を発生させるnp.ranodm.randomを用いて100ステップまでの100万回試行の1Dランダムウォークを模擬する1000000x100行列を生成する。

import numpy as np
import matplotlib.pyplot as plt

num=1000000
step=100
a=np.random.random((num,step))
a[a>=0.5]=1
a[a<0.5]=-1

ここで生成した行列を列(横)方向に和をとると100ステップ経過したときの位置の分布が得られる。

sa=np.sum(a,axis=1)

100ステップ経過した分布をsa100, 上記でstep=10として再計算した10ステップ経過した分布をsa10としてヒストグラムを表示する。

plt.hist(sa100,bins=40,color='b')
plt.hist(sa10,bins=10,color='r',alpha=0.8)
plt.legend(['step=100','step=10'])

若干崩れているが、0を中心に正規分布状に頻度が分布し、ステップが10から100になると広がっていることがわかる。

2. 中心にいる確率の計算値

　上記100万列の行列のうち列方向の和が0になる比率がランダムウォークであるステップ経過したときに中心0にある確率となる。ランダムウォークで中心に戻ってくる確率は、各ステップ（２ステップごと）の確率を積算する。ただし、まず簡易的に重複を許す（一度中心に戻った列も利用する）ことにする。重複を許すため、各ステップで積算すると確率が1を超えてしまう。

num=1000000
step=100
prop_sim_1d=[]
st_s=[]
nums=[]
a=np.random.random((num,step))
a[a>=0.5]=1
a[a<0.5]=-1

for i in np.arange(2,step+2,2):
 sa=np.sum(a[:,:i],axis=1)
 ind=np.where(sa==0)
 c=np.size(ind)
 p=c/num
 prop_sim_1d.append(p)
 st_s.append(i)
cm_prop_sim_1d=np.cumsum(prop_sim_1d)
plt.plot(st_s,cm_prop_sim_1d)
plt.title('Numerical Simulation')
plt.ylabel('Cumulative probability')
plt.xlabel('Steps')

3. 中心にいる確率の理論値

ところで、あるステップ数経過したときに中心0にいる確率は簡単に計算することができる。まず0にいるためにはステップ数が必ず偶数である必要があるためステップ数を2nとして、nは1以上の自然数とする。次に中心0になるためには、2nのうち+1と-1が半数づつになるため乱数+1,-1がちょうど半分づつになる確率は

$\displaystyle u_{2n}=\frac{1}{{2}^{2n}} \binom{2n}{n}$

上をステップごとに積算していけばよい。具体的な計算コードはいろいろあるがSympyを用いて計算すると以下のようになる。

from sympy import *
init_printing(False)
var('n  k')

step=100
prop_theo_d=[]
st_d=[]
for i in np.arange(1,int(step/2)+1,1):
 p=Sum(1/(2**(2*n))*factorial(2*n)/(factorial(n)*factorial(n)),(n,1,i)).evalf()
 prop_theo_d.append(p)
 st_d.append(int(2*i))

plt.plot(st_d,prop_theo_d)
plt.title('Theoretical calculation')
plt.ylabel('Cumulative probability')
plt.xlabel('Steps')

理論値と乱数行列を使った計算値を比較するとぴったり一致する。

plt.plot(st_d,prop_theo_d)
plt.plot(st_s,cm_prop_sim_1d,'r--')
plt.title('Theoretical vs Simulation')
plt.ylabel('Cumulative probability')
plt.xlabel('Steps')
plt.legend(['Theory','Simulation'])