前回に引き続きPythonによるランダムウォークのシミュレーションについて取り上げる。今回は、1Dランダムウォークの問題であるステップ(2n)ではじめて原点に戻る確率について考える。
1. はじめて中心に戻る理論確率
今回はランダムウォークでの2nステップ目ではじめて中心0にいる確率を求める。理論式は複雑なので以下を参照する。
2nステップ目で0になる確率f2nは
n=1から∞までの和が1に収束するため、1Dランダムウォークは必ず中心0に戻ることが数学的に得られる。
上の式をSymPyを用いて以下のように数値計算してみる。
from sympy import * init_printing(False) var('n k') import numpy as np import matplotlib.pyplot as plt step=100 st_t=[] prop_theo=[] for i in np.arange(1,int(step/2)+1,1): p=(1/(2*n*2**(2*n-2))*factorial(2*n-2)/(factorial(n-1)*factorial(n-1))).subs(n,i).evalf() prop_theo.append(p) st_t.append(int(2*i)) cm_prop_theo=np.cumsum(prop_theo) plt.plot(st_t,cm_prop_theo) plt.ylim([0,1]) plt.title('Theoretical Calculation') plt.ylabel('Cumulative probability') plt.xlabel('Steps')
100ステップでだいたい0.9程度になっている。階乗の直接計算のリソースの関係でステップ数を増やすと時間がかかる。(Stirlingの近似など使うとよいかもしれない)
2. はじめて中心に戻るシミュレーション確率
前回のコードを継承し、1000000x100の+1,-1を含むランダム行列を用いる。前回は行列の各ステップでの部分和が0になる行数を確率としていたが、今回は部分和が0になった行をnp.delete
を用いて削除する。ステップが進むごとに行数が減っていき、乱数シミュレーションとしてはサンプル数が少なくなるという懸念があるが、理論計算と同じく100ステップまで計算すると同じような曲線が得られた。
num=1000000 step=100 prop_sim=[] st_s=[] 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) a=np.delete(a,ind,0) prop_sim.append(c/num) st_s.append(i) cm_prop_sim=np.cumsum(prop_sim) plt.plot(st_s,cm_prop_sim) plt.ylim([0,1]) plt.title('Numerical Simulation') plt.ylabel('Cumulative probability') plt.xlabel('Steps')
そこで、理論値と比較するとぴったり一致していた。
plt.plot(st_t,cm_prop_theo) plt.plot(st_s,cm_prop_sim,'r--') plt.title('Theoretical vs Simulation') plt.ylabel('Cumulative probability') plt.xlabel('Steps') plt.legend(['Theory','Simulation'])
3. まとめ
今回は2nステップ目ではじめて中心0に戻る確率についての計算について紹介した。