StatModeling Memorandum

StatModeling Memorandum

StanとRとPythonでベイズ統計モデリングします. たまに書評.

蟻本シリーズ 2 ランダムウォーク

今回は以下のランダムウォークの問題を考えます。

I×Jの大きさのグリッドがあります。(1,1)からスタートして、1ターンに上下左右4マスのうち移動できる方向にそれぞれ確率p1,p2,p3,p4で移動します。いくつかのマスには石が置いてあり、通行不可能になっています。(I,J)にはじめて辿り着くまでにかかるターン数の期待値を求めなさい。ただし、(1,1)から(I,J)に移動するパスが少なくとも1つは存在すると仮定します。

例:I = 3, J = 10

     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,]    1    0    1    1    1    0    1    1    1     0
[2,]    1    0    1    0    1    0    1    0    1     0
[3,]    1    1    1    0    1    1    1    0    1     1

0が石があるマスで、1が移動できるマスです。以降ではこのグリッドを「グリッドA」と呼びます。

蟻本の4-1節(p.257)では、上下左右へ移動する確率が等確率の場合を扱っています。

この問題は期待値についての連立方程式を立てて解くと速いようです。

ここで、{ E (i,j)}{ (i,j)}からスタートした場合にゴールするまでにかかるターン数の期待値とします。ただし、{i}は行方向(上下方向)、{j}は列方向(左右方向)の添え字です。{ (i,j)}から上下左右に移動可能な場合には、次の関係式が成り立ちます。

f:id:StatModeling:20201106190147p:plain

なぜなら、{ p_1}の確率で{ (i-1,j)}へ移動する場合は期待値は{ (i-1,j)}からの期待値{ E (i-1,j)}に今移動に使った1ターンを加えたもの、{ p_2}の確率で{ (i+1,j)}へ移動する場合は期待値は{ (i+1,j)}からの期待値{ E (i+1,j)}に今移動に使った1ターンを加えたもの、…となり、各場合は排他的なので和になるからです。この式を整理すると以下になります。

f:id:StatModeling:20201106190153p:plain

ここで、例えば{ p_3 = 0}(左のマスに石があっていけない)の場合を考えます。{ p_3 = 0}になった影響で{ p_1}{ p_4}が正規化されます。すると、同様の考察で以下の式が成り立ちます。

f:id:StatModeling:20201106190157p:plain

つまり、{ p_k = 0}の場合にも拡張すると、以下の関係式が成り立ちます。

f:id:StatModeling:20201106190201p:plain

これをすべての{(i,j)}の組み合わせについて行列形式で書き下します。{ E (i,j)}から移動する場合だけを抜き出して書くと以下になります。

f:id:StatModeling:20201106190205p:plain

_0_は0が続くことを表しています。この式は期待値に関する連立方程式になっています(以下★式と呼びます)。これは行列演算で簡単に解くことができます。


以上は自分用の備忘録でした。以下では次のような問題を考えます。

グリッドAがあります。あるロボットが(1,1)からスタートして、1ターンに上下左右4マスのうち移動できる方向にそれぞれ未知の確率p1,p2,p3,p4で移動します。いくつかのマスには石が置いてあり、通行不可能になっています。そして、(I,J)にはじめて辿り着くまでにかかるターン数の期待値を10回測定しました。結果は測定誤差を含めてY <- c(68.8, 75.4, 111.2, 81.4, 82.6, 114.3, 89.2, 54.7, 66.3, 71.1)でした。この結果からp1,p2,p3,p4を推定することを考えます。

Stanでの実装例は以下になりました。

  • 2~36行目: 方程式を解いて期待値を返す関数です。
  • 3~5行目: Aは★式の左辺の行列、bは★式の右辺のベクトル、resは★式の左辺の期待値のベクトルです。
  • 7~10行目: はじめにゼロで埋めています。
  • 14~15行目: ゴール地点と石があるマスの期待値はゼロです。そのため、このように{ E (i,j)}に掛かる係数を1にしておけば、右辺はゼロ埋めのままなので、{ E (i,j) = 0}となります。
  • 17~30行目: Abの係数を冒頭の議論をもとに埋めています。
  • 34行目: Stanの行列演算を使って方程式を解いています。内部ではC++のEigenの関数が呼ばれています。
  • 35行目: res[1,1]が(1,1)から出発した場合のターン数の期待値になります。
  • 65行目: 期待値に測定誤差が加わってYになっています。

実行するRコードは以下の通りです。

推定結果は以下になりました。

meanse_meansd2.5%25%50%75%97.5%n_effRhat
p[1]0.250.000.120.040.150.230.350.487101.00
p[2]0.260.000.130.050.150.250.360.529211.00
p[3]0.080.000.050.000.030.070.120.1710321.00
p[4]0.420.010.150.200.290.390.520.777431.00
sigma11.790.103.437.159.4211.1813.4220.0110771.00
mu75.950.063.8168.3473.5875.8578.3683.5440001.00
lp__-33.240.051.54-37.11-34.05-32.90-32.10-31.289321.00

このグリッドAにおいて、上下左右へ移動する確率が等確率の場合には期待値は361.0になります。今回与えた測定データはおよそ平均が80なので、比較するとだいぶ小さい期待値です。そのため、右へ移動してゴールに近づく確率(p[4])がかなり高く、左へ移動する確率(p[3])がかなり低く推定されています。