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

今回は、前回の1Dランダムウォークの問題であるステップ(2n)ではじめて原点に戻る確率を計算する際に階乗の部分をStirling近似を用いたので紹介する。

atatat.hatenablog.com

1. はじめて中心に戻る理論確率
2. Stirling近似の導入
3. 中心に戻る確率計算の実行
4. まとめ

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

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

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

n=1から∞までの和

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

を以下数値計算することで、1Dランダムウォークで中心に戻る確率を計算したが、階乗を計算する関係で100ステップ程度で打ち切っていた。(100ステップで確率0.9程度)

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

今回は階乗の部分をStirling近似に置き換えることで、100ステップより増やして1に近づくことを示す。

2. Stirling近似の導入

Stirling近似はnが大きいときに成立する近似式で以下のように表せる。

$\displaystyle n! \sim \sqrt{2\pi n} \left( \frac{n}{e} \right)^{n}$

上の数値計算のコードの階乗factorialの部分をStirling近似に置き換えれば階乗の計算リソースを大幅に削減できることが期待される。

ただし、Stirling近似で計算は簡単になるが、nⁿで想像できるようにnが増えるに従って大きな桁になり、整数のビット数(16bit/32bit)では足りなくなるという現象が起こる。

解消法としては1度対数に変換して計算を実行し、最後にexpで戻すという操作することで回避できる。

以下はStirling近似のlogを計算する関数を定義する。

def stirf_log(x):
  if x==0:
    s=1 #ダミー, 実際には使わない
  else:
    s=np.log(np.sqrt(2*np.pi*x))+x*np.log(x/np.e)
  return s

3. 中心に戻る確率計算の実行

Stirling近似はnが大きい場合に成立（誤差が小さい）するため、100ステップまでは今まで通り計算し100以上はStirling近似を適用する。 Stirling近似を使用するときは、すべて一度logをとって計算し、expで戻している。

step=20000
st_t_st=[]
prop_theo=[]
for i in np.arange(1,int(step/2)+1,1):
  if i<50:
    p=(1/(2*n*2**(2*n-2))*factorial(2*n-2)/(factorial(n-1)*factorial(n-1))).subs(n,i).evalf()
  else:
    lnp=stirf_log(2*i-2)-(np.log(2*i)+(2*i-2)*np.log(2))-2*stirf_log(i-1)
    p=np.exp(lnp)
  prop_theo.append(p)
  st_t_st.append(int(2*i))

cm_prop_theo_st=np.cumsum(prop_theo)
plt.plot(st_t_st,cm_prop_theo_st)
plt.ylim([0,1])
plt.title('Theoretical Calculation')
plt.ylabel('Cumulative probability')
plt.xlabel('Steps')