つれづれなる備忘録

日々の発見をあるがままに綴る

Pythonによるデータ処理18 ~ 1Dランダムウォーク2

前回に引き続きPythonによるランダムウォークのシミュレーションについて取り上げる。今回は、1Dランダムウォークの問題であるステップ(2n)ではじめて原点に戻る確率について考える。

atatat.hatenablog.com

1. はじめて中心に戻る理論確率

今回はランダムウォークでの2nステップ目ではじめて中心0にいる確率を求める。理論式は複雑なので以下を参照する。

zakii.la.coocan.jp

2nステップ目で0になる確率f2n

 \displaystyle f_{2n}=\frac{1}{2n 2^{2n-2}} \binom{2n-2}{n-1} \,(n\geq 1)

n=1から∞までの和が1に収束するため、1Dランダムウォークは必ず中心0に戻ることが数学的に得られる。

 \displaystyle \sum_{n=1}^{\infty} f_{2n}= 1

上の式を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')

"はじめて0に戻る確率"
はじめて0に戻る確率

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')

"はじめて0に戻るシミュレーション確率"
はじめて0に戻るシミュレーション確率

そこで、理論値と比較するとぴったり一致していた。

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'])

"はじめて0に戻る確率比較"
はじめて0に戻る確率比較

3. まとめ

今回は2nステップ目ではじめて中心0に戻る確率についての計算について紹介した。

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

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

zakii.la.coocan.jp

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になると広がっていることがわかる。

"1Dランダムウォーク分布"
1Dランダムウォーク分布

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

 上記100万列の行列のうち列方向の和が0になる比率がランダムウォークであるステップ経過したときに中心0にある確率となる。 ランダムウォークで中心に戻ってくる確率は、各ステップ(2ステップごと)の確率を積算する。ただし、まず簡易的に重複を許す(一度中心に戻った列も利用する)ことにする。重複を許すため、各ステップで積算すると確率が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')

"1Dランダムウォーク計算値"
1Dランダムウォーク計算値

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')

"1Dランダムウォーク理論値"
1Dランダムウォーク理論値

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

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'])

"1Dランダムウォーク比較"
1Dランダムウォーク比較

4. まとめ

今回は1Dランダムウォークのシミュレーション、原点に戻ってくる確率について紹介した。今回は簡単のため重複を許すケースとしたが、次回は重複しないケースについて紹介する。

今週のお題「わたしは○○ナー」

今週のお題「わたしは○○ナー」

今週のお題「わたしは○○ナー」ということで、例にあるのであまり面白くはないが、ずばりランナー(正確にはジョギングに近い)です。

10月に入ってきて、気温・湿気かなり運動に適してきており、最近は走っていても大量に汗をかくことがなくなってきた。

ところで、ランニングのオンラインイベントというものに初めて参加してみた。

vrwc.runtrip.jp

距離はとりあえず10km, いつも走っている距離なので手持ちのGarminデータのアップロードだけがうまくいくか心配だったが難なく完走証をもらうことができた。

毎月やっているようなので次はハーフマラソンでもよいかもしれない。

あとジムで筋トレも習慣的にやっているので、トレーナーといいたいところだが、筋トレをする人はトレーニーで、教えるコーチがトレーナーになる。