整数格子上のランダムウォークは2次元以下だとスタート地点に戻ってくるが、3次元だと戻ってこない、という話です。
Pythonでシミュレーションを行い、escape probabilityの収束(する/しきらない)様子を確認します。
Foundation of Data ScienceのCahpter4がマルコフ連鎖~MCMC~ランダムウォークの話になっていて面白かったので。
問題設定
d次元整数格子において、原点から出発してランダムウォーク(=2d個の隣接点からランダムに選んで移動)を無限回行います。このとき、原点に再び戻ってくることがあるのか、いずれ遠くに行って2度と帰ってこなくなるのかが気になります。
これをシミュレーションによって近似的に求めるために、nステップの長さのランダムウォークをN回行い、1回以上原点に戻ってきた経験的確率を求めます。
実験のイメージ
3次元のランダムウォーク(100000ステップ)の様子をこちらを参考に可視化してみました。この例だと原点(左上)のあたりから始まって右下の方に行って戻ってこれなさそうです。
この試行を何度も繰り返して、原点に戻ってくる試行が現れる確率を求めます。
コード
numpy配列演算や日本語表示などの工夫をしていますが、説明は省略します。
% matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
import seaborn as sns
sns.set()
matplotlib.rcParams['font.family'] = 'sans-serif'
matplotlib.rcParams['font.sans-serif'] = ['Hiragino Maru Gothic Pro', 'Yu Gothic',
'Meirio', 'Takao', 'IPAexGothic', 'IPAPGothic', 'VL PGothic', 'Noto Sans CJK JP']
def try_mc(n, N, dim):
"""dim次元格子において、原点から出発してnステップ移動するランダムウォークをN回行う
"""
index = np.random.randint(dim, size=(N, n))
sign = np.random.choice([-1,1], size=(N, n))
ps = np.zeros((N, n,dim), dtype=int)
ps[np.arange(N).reshape(N,1), np.arange(n), index] += sign
ps = ps.cumsum(axis=1)
return ps
n = 100000
N = 1000
dims = [1,2,3]
fig, axes = plt.subplots(1,3, sharey=True, sharex=True, figsize=(15,5))
for i, dim in enumerate(dims):
ps = try_mc(n, N, dim)
ds = np.abs(ps).sum(axis=2)
ax = axes[i]
ax.set_xlabel("ランダムウォークの長さ")
if i==0:
ax.set_ylabel("一回以上原点に戻ってくる確率")
prob = ((ds == 0).cumsum(axis=1)>0).sum(axis=0) / N
ax.plot(prob)
シミュレーション結果
ランダムウォークを続ければいつかは戻ってくる、のであれば確率が1に収束してほしいところですが...
* 1次元ではすぐに1に収束しています。
* 2次元も本当は1に収束するはずですがかなり遅く、100000ステップだとまだまだという感じです。参考文献を見ると1次元のescape probabilityが $O(1/n)$ で上からバウンドできているのに対して、2次元は $O(1/\log n)$ でしか抑えられていないので仕方なさそうです(ここの $n$ は原点からの距離の意味です)。
* 3次元は2次元に比べると傾きが0に近く、収束しているようにも見えます。
戻ってくる確率が0.34くらいです(escape probabilityが0.66)。
1,2次元だと戻ってくる確率が1なので十分長い時間続ければいずれ原点に戻ってきますが、3次元では「戻ってこない確率が0より大きい」、よって十分長い時間続けているといずれ必ず二度と帰ってこないということになります(途中で1度や2度戻ってきても再出発を繰り返せばいずれ戻ってこない試行を引くので)。
(厳密には戻ってこない、のではなく原点に戻ってくる前に無限遠に行く、なので表現がよく分かりませんが。。)
背景の話
数学的には、escape probabilityという値を
$$p_{escape} (i,j) := (iから出発して、jに到達する前にiに戻る確率)$$
と定義します。今回の格子グラフでは $i=原点$, $j=$ {原点から(マンハッタン)距離がnの点}として、$n\to \infty$ のときの $p_{escape}$ を求めていたイメージです。
$dim=1,2$ のときは $p_{escape} = 0$ で、$dim=3$ のときは
$$1/6 <= p_{escape} <= 5/6$$
であることが比較的簡単に証明できます。
この証明にはグラフを各エッジに1Ωの抵抗を挟んだ電気回路において、escape probablityが
$$ p_{escape} = c_{eff} / c_i $$
($c_{eff}$ はi,j間の合成コンダクタンス、$c_i$ は $i$ に接続するエッジのコンダクタンスの和)
と表されることを用いるなどテクニカルに面白いので、興味を持ったらぜひ参考文献を御覧ください。
escape probability以外にも、電圧、電流とランダムウォークの関係、commute timeと合成抵抗の関係などの話題があり、面白いです。
(ちなみに今回の「計算機実験でescape probabilityを求める」というのが演習問題に含まれています。)