Joint distributions of numbers of trials and returns to the

吸収壁に到達するまでの試行数と
原点通過回数の同時分布
城西大 理学研究科 数学専攻 平野勝臣・吉田倫大
統計科学とその周辺/2014/12/7/
於 城西大学紀尾井町キャンパス
発表内容.
結果.d を正の整数とする.原点を出発したランダムウォークが初めて −d または d に到達す
るまでの試行数とその時までに原点を通過する回数の同時分布を導いた.
統計的応用.ランダムウォークが独立な試行に基づいて動くとき,初めて −d または d に到達
するまでの試行数の観測値に基づいて,パラメータ p の推定を試みる.
ランダムウォークが P´
olya urn sampling に基づいて動くとき,初めて −d または d に到達す
るまでの試行数の観測値に基づいて,壺の中の玉の構成を決めているパラメータの推定を試みる.
ランダムウォーク.
数直線上の動点は 1 回の試行で確率 p で右に 1 動き,確率 q = 1 − p で左に 1 動く.この試
行を独立に繰り返す.ランダムウォークは確率モデルの例としてしばしば用いられている.ある
状態を出発して,初めてある状態に到達するまでの推移の回数は興味ある問題である.
本報告では,独立試行に基づいて数直線上の点
{−d, −d + 1, . . . , 0, 1, . . . , d − 1, d},
d は正の整数
を動く動点が,原点を出発し,初めて −d または d に到着するまでの試行数とその時までに原点
を通過する回数の同時分布の問題を扱う.さらに,独立試行を exchangeable sequence とした場合
の問題を調べる.
Feller の 14 章で,{0, 1, . . . , z, . . . , a} 上を動く点が z を出発し,0 に到達するまでの試行数,
a に到達するまでの試行数,0 または a に到達するまでの試行数,の生成母関数を求めている.
{−1, 1}− 値の確率変数の系列 ξ1 , ξ2 , . . . は独立で同一の分布に従う (iid) とする.ここで
P (ξi = 1) = p, P (ξi = −1) = q (1 < p < 1, q = 1 − p), i = 1, 2, . . .
とする.
Xn = ξ1 + · · · + ξn , X0 = 0 と置くと,{Xn , n = 0, 1, 2, . . .} は原点を出発する iid 系列に基
づくランダムウォークである.
ξ1 , ξ2 , . . . が exchangeable ならば,{Xn , n = 0, 1, 2, . . .} は exchangeable sequence に基づくラ
ンダムウォークである.
各 n に対して (ξ1 , ξ2 , . . . , ξn ) の分布が置換に対し不変であるとき,ランダム系列 {ξ1 , ξ2 , . . .}
は exchangeable であるという.
はじめに,iid 系列に基づく場合の同時分布を求める.次に exchangeable 試行に基づく場合の
同時分布を求める.
試行数と原点通過の同時分布
独立試行に基づくランダムウォークの場合の同時分布を調べる.
数直線上を動くランダムウォークを {Xn , n = 0, 1, 2, . . .} とし,
P (Xn+1 = i + 1|Xn = i) = p, P (Xn+1 = i − 1|Xn = i) = q = 1 − p, 0 < p < 1
1
とする.また d を正の整数とする.原点 0 を出発したランダムウォークが d または −d に到達す
るまでの待ち時間を W ,それまでに原点を通る回数を N とする.すなわち,
W
N
= min{n; Xn = d or Xn = −d},
= number of {m; Xm = 0 and 1 5 m < W }
とする.W と N の同時確率母関数を求める.
W と N の conditional joint pgf を次のようにおく.
φi = E[tW sN |X0 = i], for i = −d + 1, ..., d − 1.
φ = φ0 .
このとき,条件付き確率母関数は次の漸化式を満たす.

φi
= ptφi+1 + qtφi−1 , (i = −d + 2, ..., −2 or i = 2, ..., d − 1)




φ
= ptsφ + qtφ−2

−1


φ
= ptφ1 + qtφ−1
(∗) :
φ
= ptφ2 + qtsφ

1



φ
= pt + qtφd−2

d−1


φ−(d−1) = ptφ−(d−2) + qt.
この方程式を φ について解く.
2 変数多項式の列 An (s, t) を,A1 (s, t) = 1, A2 (s, t) = 1 − 2pqt2 s と,漸化式
An (s, t) = An−1 (s, t) − pqt2 An−2 (s, t), (n = 3, 4, . . . , d)
によって定めると,
φ(t, s) =
pd td + q d td
Ad (s, t)
が成り立つ.
さらに,Ad (s, t) を陽に表すと,次の定理を得る.
定理 1. W と N の同時確率生成母関数は
φ(t, s) =
(pd + q d )td
,
Ad (s, t)
d=2
で与えられる.ここに,
1 − 2pqt2 s
pqt2
Ad (s, t) = √
{ξ(t)d−1 − η(t)d−1 } − √
{ξ(t)d−2 − η(t)d−2 },
1 − 4pqt2
1 − 4pqt2
ξ(t)
=
η(t)
=
√
1 − 4pqt2
,
√2
1 − 1 − 4pqt2
.
2
1+
φ(t, 1), φ(1, s) はそれぞれ,W および N の周辺確率母関数である.
2
具体例.各 d に対して Ad (s, t) は次で与えられる.
A2 (s, t) =
1 − 2pqst2
A3 (s, t) = (−2pqs − pq)t2 + 1
A4 (s, t) = 2p2 q 2 st4 + (−2pqs − 2pq)t2 + 1
A5 (s, t) = (4p2 q 2 s + p2 q 2 )t4 + (−2pqs − 3pq)t2 + 1
A6 (s, t) =
A7 (s, t) =
−2p3 q 3 st6 + (6p2 q 2 s + 3p2 q 2 )t4 + (−2pqs − 4pq)t2 + 1
(−6p3 q 3 s − p3 q 3 )t6 + (8p2 q 2 s + 6p2 q 2 )t4 + (−2pqs − 5pq)t2 + 1
A8 (s, t) = 2p4 q 4 st8 + (−12p3 q 3 s − 4p3 q 3 )t6 + (10p2 q 2 s + 10p2 q 2 )t4 + (−2pqs − 6pq)t2 + 1
A9 (s, t) = (8p4 q 4 s + p4 q 4 )t8 + (−20p3 q 3 s − 10p3 q 3 )t6
+(12p2 q 2 s + 15p2 q 2 )t4 + (−2pqs − 7pq)t2 + 1
A10 (s, t) =
−2p5 q 5 st10 + (20p4 q 4 s + 5p4 q 4 )t8 + (−30p3 q 3 s − 20p3 q 3 )t6
+(14p2 q 2 s + 21p2 q 2 )t4 + (−2pqs − 8pq)t2 + 1.
系 1.d = 2 のとき,W = 2N + 2 の関数関係がある.
N の周辺分布.
幾何分布 G(p):成功確率 p のベルヌーイ試行で初めて成功が起こるまでの失敗の回数の分布.
定理 1 の ξ(t) と η(t) で,ξ = ξ(1),η = η(1) と置く.p 6= 1/2 ならば
φ(1, s) =
ここに,
pd =
(ξ − η)(ξ d + η d )
ξ d − ηd
pd
.
1 − qd s
and qd =
2ξη(ξ d−1 − η d−1 )
ξd − ηd
である.η < ξ に注意して 0 < pd < 1, 0 < qd < 1, pd + qd = 1 であることから,次の系 2 を
d
得る.p = 1/2 ならば φ(1, s) = d−(d−1)s
であることから次の系 3 を得る.
系 2. p 6= 1/2 ならば N は G(pd ) に従う.
系 3. p = q = 1/2 ならば N は G( d1 ) に従う.
各 d に対して N の分布は次で与えられる.
d=2
G(q 2 + p2 )
d=3
+p
G( q1−pq
)
d=4
q +p
G( 1−2pq
)
d=5
q +p
G( 1−3pq+p
2 q2 )
d=6
q +p
G( 1−4pq+3p
2 q2 )
d=7
q +p
G( (1−5pq+6p
2 q 2 −3p3 q 3 )
d=8
q +p
G( 1−6pq+10p
2 q 2 −4p3 q 3 )
d=9
G( 1−7pq+15p2qq2+p
−10p3 q 3 +p4 q 4 )
d = 10
G( 1−8pq+21pq2 q2+p
−20p3 q 3 +5p4 q 4 ).
3
3
4
4
5
5
6
6
7
7
8
8
9
10
3
9
10
W の周辺分布
W の確率母関数 φ(t, 1) を展開し,W の確率関数を求めることができる.
Exchangeable case:
ランダムウォークがポイヤサンプリングの結果に従って動くときを考察する.
α 個の赤玉と β 個の黒玉の入った壺から無作為に玉を一つ抽出し,同じ色の玉を 1 個加えて壺
に戻す.この操作を繰り返す.i 番目の抽出で赤なら右へ 1, 黒なら左へ 1 進むとする.抽出された
玉の色の系列は exchangeable である.この試行に基づくランダムウォークを考える.p = α/(α+β)
と置き,p の従う分布の密度関数 π(p) をベータ密度
π(p) =
Γ(α + β) α−1
p
(1 − p)β−1
Γ(α)Γ(β)
とする.このとき,ポイヤサンプリングに基づくランダムウォークが,初めて d または −d に到
達するまでの試行数 W の確率関数は,iid の確率関数を P (W = n) とすると,
∫
1
P (W = n)π(p)dp
0
で与えられる.
例えば,初めて 5 または −5 に到達するまでの試行数 W の確率母関数は次で与えられる.
P (W = 2k + 1) = c5 (k)
(α + k − 3) · · · α(β + k − 3) · · · β
g5 (α, β), k = 2, 3, . . .
(α + β + 2k) · · · (α + β)
ここに,
g5 (α, β) = (α + k + 2) · · · (α + k − 2) + (β + k + 2) · · · (β + k − 2).
注意.系 1 の W と N 関係 W = 2N + 2 は試行列が exchangeable 系列,マルコフ系列でも成立
する.
パラメータ推定
iid 系列
数直線上の原点から出発したランダムウォークが,1 回の試行で確率 p で右に 1 進み,確率 1−p
で左に 1 進む.初めて −d または d に到達するまでの試行数の大きさ n の観測を W1 , W2 , . . . , Wn
とし,この観測に基づいて p を推定する.
p = 0.3 のとき,初めて −5 または 5 に到達するまでの試行数 W のシミュレーションデータから
p の最尤推定を試みる.
W = k が j 回観測されたことを [k, j] と表す.50 個のシミュレーションデータは
[5, 12], [7, 7], [9, 7], [11, 5], [13, 7], [15, 2], [17, 2], [19, 1], [21, 4], [23, 1], [25, 1], [31, 1]
であった.尤度関数は
L(p) = {P (Wi1 = 5)}12 {P (Wi2 = 7)}7 {P (Wi3 = 9)}7 · · · P (Wi12 = 31)
で,対数尤度 log L(p), (0 < p < 1) を最大化する p を求めれば,p の最尤推定値 p = 0.2825 を
得る.
4
Exchangeable 系列
壺の中に赤玉 α 個と黒玉 β 個が入っている.ランダムに γ 個の玉を壺から抽出し,それが赤
玉なら赤玉を γ 個加えて壺に戻し,黒玉ならば黒玉を 1 個加えて壺に戻す.この操作に基づいて
原点を出発したランダムウォークが,抽出した玉が赤玉なら右に 1 進み,黒玉なら左に 1 進む.
初めて d または −d に到達するまでの試行数の大きさ n の観測を W1 , W2 , . . . , Wn とし,この観
測に基づいて α と β の最尤推定を試みる.このような問題を解析的に扱うことは困難であるが,
数値的に扱うことは可能である.
ランダムウォークが吸収壁に到達するまでの観測値から α, β の推定する.このような α, β の
推定の例は Aki [1] にある.
d = 2 のとき φ(t) = (p2 + q 2 )t2 /(1 − 2pqt2 ),π(p) としてベータ密度をとり,t2k の係数
P (W = 2k) を求めると
P (W = 2k) =
2k−1 Γ(α + β)
{Γ(α + k + 1)Γ(β + k − 1) + Γ(α + k − 1)Γ(β + k + 1)}
Γ(α)Γ(β)Γ(α + β + 2k)
である.
この結果を用いて,ランダムウォークが Polya サンプリングに従って移動したとき,初めて d
または −d に到達するまでの試行数の観測値から,壺の中の初期個数 α, β を推定してみる.
α = 2, β = 2 としたとき,大きさ 10 のシミュレーションのデータ 2, 2, 6, 14, 2, 2, 4, 4, 2, 2 か
ら尤度関数は
L(α, β) = P (W1 = 2)P (W2 = 2)P (W3 = 2) · · · P (W10 = 2)
で,対数尤度を最大化した結果,α = 6.061 と β = 6.061 を得る.この最大化は Maxima を用い
た.観測数 10 は推定すべきパラメータ 2 個に対しては小さいだろう.
そこで,β = 2 とし推定すべきパラメータ α ひとつにすれば,同様にして,α = 2.475 を得
る.観測数を大きくすれば計算に支障をきたすことも予想される.n = 50 程度で安定するだろう.
参考文献
[1] Aki, S. (2012). Statistical modeling for discrete patterns in a sequence of exchangeable trials,
Annals of the Institute of Statistical Mathematics, 64 633-655.
[2] Feller, W. (1968). An Introduction to probability theory and its applications, Vol. 1, Third
edition, Wily, New York.
[3] Inoue, K., Aki, S. and Hirano, K. (2011). Distributions of simple patterns in some kind of
exchangeable sequences, Journal of Statistical Planning and Inference, 141 2532-2544.
[4] Ross, S. M. (2014). Introduction to probability models, 10 th edition, Academic Press.
5