MCMC 法
EM
欠測時系列
個体群
最尤推定
1,2
箱山 洋
1
/東京海洋大学2
Hiroshi Hakoyama1,2
Fisheries Research Agency1 /Tokyo University of Marine Science and Technology2
水産総合研究
1
時系列
値 (missing value)
体制・予算
場合
、途中
(図 1)。天候・実行
多
不確実性
、野外調査 継続的
行
難
列
(incomplete time-series data)
当
。欠測値
、
。
従
不完全
最尤推定
、確率過程
場合
欠測
個体群
問題
考
離散時間
、
欠測値
時系
過程
含
不完全
場合
、連続時間
仮定
Population size
個体数
1200
1000
800
600
5
最尤推定
考察
個体群動態
行
。
最尤推定
短
離散時間 t
Xt > 0
陽
仮定
考察
図1
(t = 1, 2, . . . , n)、個体数変動
、X1 , X2 , . . . , Xn
欠測値
25
時系列
、時刻 t − 1
。
、時刻 t
個体数 xt−1
条件付
個体数
分布
、θ = (θ1 , θ2 , . . . , θs )
周辺密度
。
最尤推定量
完全
時系列
ln Lc (θ; x) = ln f (x; θ) =
x = (x1 , x2 , . . . , xn )
n−1
∑
得
、
対数尤度
ln f (xt+1 | xt ; θ) + ln f (x1 ; θ),
(1)
t=1
。
30
s(個
表 。
)
∏n−1
同時密度関数 、f (x1 , x2 , . . . , xn ; θ) =
。
t=1 f (xt+1 | xt ; θ) f (x1 ; θ)
、初期個体数 X1
2.1 完全
欠測値 含 不完全
過程 従
f (xt | Xt−1 = xt−1 ; θ)
。
f (x1 ; θ)
20
過程
確率密度関数
表
15
Time
行 。
2 離散時間
10
、Lc
完全
尤度
表
。式 1
1
右辺第二項
初期個体数
周辺分布
計算
部分
、
報 捨
部分
、式 1
定量
構成
尤度
煩雑
右辺第一項
方法
量n
、θ = arg max ln L (θ; x)
難
。
、個体数初期値
条件付対数尤度関数 (conditional log-likelihood)
。
c
計算
大
定義
、右辺第二項
。以下、解析的
無視
数値的
影響
小
完全
持
情
条件付最尤推
。最尤推定量 θ
対
最尤推定値
θ
計算
場合
2.2 不完全
欠測値
考
最尤推定量
含
不完全
時系列
最尤推定量
data)
場合
。
個体群
使
考
。NA
最尤推定
。例
欠測値
表 。
、
場合、
行
場合、完全
(complete
得
x = (x1 , x2 , NA, x4 , NA, NA, x7 )
不完全
対数尤度
次
：
ln L(θ; x) = ln f (x2 | x1 ; θ) + ln f (x4 | x2 ; θ) + ln f (x7 | x4 ; θ) + ln f (x1 ; θ),
、L
観察
定義
前
不完全
尤度
(incomplete data)
与
簡単
式
表
右辺第二項
3 年前 個体数 条件付 確率密度関数 計算
方程式 (see Goel and Richter-Dyn, 1974)
∫
右辺第一項
。式 2
、欠測値
。
次
(2)
過程
第三項
2年
、
積分式
-
表
：
∞
f (x4 | x2 ; θ) =
f (x4 | x3 ; θ)f (x3 | x2 ; θ)dx3 ,
∫0 ∞ ∫ ∞
f (x7 | x4 ; θ) =
f (x7 | x6 ; θ)f (x6 | x5 ; θ)f (x5 | x4 ; θ)dx5 dx6 .
0
、式 2
条件付
対数尤度
θ
関
定 行
、条件付
間
確率密度関数
(4)
0
t−s年
、t 年
(3)
積
最大化
確率密度関数
積
s − 1 点 欠測値
一重積分
場合
s 重積分 形 表
数値的
、EM
困難
多重積分
。
(Dempster et al., 1977)
式
含
。一般
確率密度関数 f (xt | xt−s ; θ)
条件付
。
行
二重積分
、欠測値
用
項
含
尤度関数
含
工夫
、
用
必要
最尤推
。
2.3 EM
、完全
EM
付
期待値
対
、最大化
対
対数尤度関数
、
推定値
関
確率過程
Xobs 、欠測値
Xmis
、x = {xobs , xmis } = (x1 , . . . , xn )
θ
、欠測値
推定値 θ
(k)
xmis
得
対数尤度
逐次更新
推定値
行
観測
、欠測値
条件
観測
(Dempster et al., 1977)。
最尤推定値 得 方法
個体数 Xt
、
対
。
、Xt
時系列 X1 , . . . , Xn
、実現値
。
確率変数
。完全
xobs , xmis
、EM
、観測値
xobs
条件付期待値 Q
計算
2
観測値
、k 段階目
θ (k)
（E
条件
, expectation step）：
∫
Q(θ|θ
(k)
lc (xobs , xmis ; θ)fXmis (xmis |xobs ; θ (k) )dxmis ,
)=
(5)
Ω
、lc (x; θ) = ln Lc (x; θ)
。Ω
標本空間
xmis
最大化
θ
完全
(k+1)
対
。
計算
対数尤度関数、fXmis
、E
（M
書
条件付
xmis
下
Q(θ|θ
確率密度関数
(k)
)
関
θ
, maximization step）：
θ (k+1) = arg max Q(θ|θ (k) ).
(6)
θ
、EM
定値
更新
E
行
、観察
。θ
(k+1)
尤度 L
尤度関数 凹 単峰
、θ
、L(xobs ; θ
不完全
(k)
完全
(例
簡単
式
尤度 (式 1)
用
、混合正規分布
Q(θ|θ
(k)
推
、
証明
繰返
(Dempster et al., 1977)。
(1)
始
、
十分
行
、煩雑
最大化
難
不完全
対
最尤推定値
、不完全
尤度 (式 2)
得
使
。線形
推定; see McLachlan and Krishnan, 2007)、積分 外
計算
)
返
(Wu, 1983)。
収束
最尤推定
、扱
繰
満
)
空間上 任意 初期値 θ
、EM
場合
) ≥ L(xobs ; θ
(k)
単調増加
、
最尤推定値
手順
M
(k+1)
、比較的簡単
手順 構成
EM
。
積分
2.4
、
考
非線形
評価
、Q(θ|θ
(k)
)
最大化
場合
、式 5
右辺
積分
考
：
Q(θ|θ (k) ) =
(i)
難
煩雑
、一般
。解析的
形
普通
完全
。
書
積分
解析的
積分
置
下
、式 5
η
1∑ c
(i)
l (xobs , xmis ; θ),
η i=1
i 番目 値
(k)
右辺
、式 5
最尤推定値
)
。l
c
書
法
式
、Q(θ|θ
最大化
θ
構成
EM
下
換
(7)
、確率密度関数 fXmis (xmis |xobs , θ (k) )
、xmis
場合、Q(θ|θ
過程
。
(k)
発生
書
)
、Q(θ|θ
(k)
方法 、MCEM(Monte Carlo EM)
下
欠測
式
、多
積分
)
呼
用
(Wei and
Tanner, 1990)。
2.5 乱数 発生
離散時間 非線形
次
問題
従
乱数
乱数
発生
発生
過程
。明示的
、一般
対
最尤推定値 計算 、MCEM
確率密度関数 fXmis (xmis |xobs , θ
簡単
。
3
確率分布
(k)
従
)
用
書
乱数
下
発生
、
、
方法
、分
布関数
逆関数
書
考
下
場合
逆関数法
、確率密度関数
(see Knuth, 1997)、
、fXmis (xmis |xobs , θ
用
1953)
(k)
上
覆
関数
見
場合
棄却法
MCMC(Markov chain Monte Carlo) 法 (Metropolis et al.,
従
)
疑似乱数
発生
考
。
2.5.1 Metropolis-Hastings algorithm
、xmis
問題
一度
考
発生
。三年間 (t = 1, 2, 3)
確率分布 、条件付
確率
f (x2 | x3 , x1 ; θ) =
f (x2 | x3 , x1 ; θ)
従
発生
疑似乱数
分散
過程
乱数
持
方法
x2
疑似乱数
、欠測値 x2
次
発生
条件付
：
満
非負実数）
(
疑似乱数列
′
次 確率
単峰形
、
考
値
提案分布
分散
x2 = x2
(i)
、発生
、新
発生
乱数
固定
目的
平均
提案分布
決
f (x2 | x3 , x1 )
平均 m
分布
提案分布
疑似
疑似乱数
。切断正規分布
場合、提案分布
期待値
疑似乱数
近似
。f (x2 | x3 , x1 ; θ)
現在
）
分布
独立連鎖
)−1
呼
棄却確率
目的
用
受容
簡単
；提案分布
決
分布
d2 f (x2 |x3 ,x1 )
dx2 2
Metropolis-Hastings algorithm
。MCMC 法
、f (x2 | x3 , x1 ; θ)
発生
、x2
方法
(8)
、計算機
正規分布
、提案分布
使
′
得
条件
、f (x2 | x3 , x1 ; θ)
。連続
乱数
合
。
−
単峰形
分布
得
一
。MCMC 法
疑似乱数
詳細釣
採択
、例
考
、MCMC 法
考
従
、分散 σ 2
場合
発生
切断正規分布（定義域
方
欠測値
方程式
-
使
確率分布
、
場合
(x1 , NA, x3 )
確率分布（一様分布
複雑
分布
不完全
定義
発生
、目的
1点
f (x2 , x3 | b1 ; θ)
f (x3 | x2 ; θ)f (x2 | x1 ; θ)
= ∫∞
.
f (x3 | b1 ; θ)
f (x3 | x2 ; θ)f (x2 | x1 ; θ)dx2
0
(Metropolis et al., 1953)
乱数
、欠測
、二峰形（多峰形）
平均
分散
提案分布
発生
疑似
：
(
)
f (x2 ′ | x3 , x1 ; θ)q(x2 )
α(x2 , x2 ) = min
,1 ,
f (x2 | x3 , x1 ; θ)q(x2 ′ )
)
(
−m)2
f (x3 | x2 ′ ; θ)f (x2 ′ | x1 ; θ) exp (x22σ
2
( ′
) , 1 ,
= min
−m)2
f (x3 | x2 ; θ)f (x2 | x1 ; θ) exp (x22σ
2
′
、q(x2 )
x2
(i+1)
= x2
提案分布
′
、棄却
積分項
場合
x2
(i+1)
。
。f
。x2 ′
切断正規分布
受容
、
単峰形
目的分布
f (x2 | x3 , x1 ; θ)
、MCMC 法
分布
乱数 x2 (i+1)
、i + 1 番目
。式 9
= x2
(9)
、採択率
高
従
、効率
提案分布
疑似乱数
疑似乱数
、
分母
発生
発生
。
2.5.2 Metropolis sampling within Gibbs method
(i+1)
疑似乱数
xmis
、i 番目
疑似乱数
一
(i+1)
xmis
一
得
発生
発生
一度 発生
。
・
(i)
欠測値
xmis 与
、i + 1 番目
i+1
i+1
i
i
、(xmis,1 , . . . , xmis,j−1 , . . . , xmis,j+1 , xmis,jmax ) 与
。
、前節
。
方法
考
欠測値
1点
(Geman and Geman, 1984)
欠測値
j 番目
条件付
場合
連続
Metropolis sampling within Gibbs method
4
xi+1
mis,j
適用
呼
。
手順
2.6 MCEM
手順
。
初期値（1 段階目
1. 初期
、
推定
2. 初期欠測値 初期
、k 段階目
仮定
対数尤度
2.
積分中
3. 完全
、最尤推定
手順 、
MCEM
1. 完全
期待値 Q(θ|θ (k) )
条件付
欠測値
疑似乱数
最尤推定量（解析
積分
過程
問題 、MCEM
計算
発生
Metropolis within Gibbs method
数値計算）
、離散時間
推定
推定）
用
、Q(θ|θ
個体数
適用
(k)
)
当
解決
最大化
、θ (k+1)
得
、
。
2.7 分散 推定
最尤推定量
量
少
分散
信頼区間
評価
場合
方法
法
用
、
数
方法
十分
大
仮定
漸近的方法
、
。
2.7.1 漸近的 分散 推定
量
十分大
、最尤推定量
行列
、完全
分散
場合
、
推定
欠測情報行列 除
ˆ x)
観測情報行列 J(θ,
知
。不完全
関係 成立
観測情報行列
、次
表
逆行列
、完全観測情報
(Louis, 1982)：
[ 2
]
ˆ xobs ) =E ˆ − ∂ log fθ (x) | ˆ Xobs = xobs
J(θ,
θ
θ=θ
∂θ∂θ ′
[(
]
)(
)′
∂ log fθ (x)
∂ log fθ (x)
− Eθˆ
|θ=θˆ
|θ=θˆ Xobs = xobs .
∂θ
∂θ ′
不完全
観測情報行列
、
積分 計算
(Wei and Tanner, 1990)。
2.7.2 Monte Carlo Bias Correction
詳
用
述
、Hakoyama and Iwasa (2000)
最尤推定量
修正
漸近的
推定
問題
正確
適用可能
、
数
少
場合
方法 (Monte Carlo Bias Correction, MCBC)
信頼区間
求
方法
開発
。今回
不完全
、
法
、MCBC
同様
個体群
。
2.8 欠測 間隔 情報量
欠測 間隔
量
変化
1 年、2 年 長
、
、欠測値
行
、
異
時間間隔
長
推定
考
。例
n点
時系列
5
情報量
？ 情報
、定常過程
時間間隔
短
平衡個体数
推定
n点
完全
情報量
情報
差
失
。
、観測
、増殖率
2点 間
欠測間隔以上
3 連続時間 確率過程
連続時間
一部
仮定
、個体群動態
確率過程
確率微分方程式
考
観察
。
、観測
、時刻 t
得
難
、数値的
強解
。
個体数 Xt
尤度関数
程式 個体群動態
考慮 値
、
率微分方程式 確率過程 対
、近似尤度 基
Hakoyama and Iwasa (2000)
、強解
時間 非線形
確率過程
予想
長
。
最大化
非線形
仮定
陽
一部
含
場合、
不完全
時系列
、
非線形方程式
強解
。仮定
時間間隔 観察
強解
得
限
式
積分項 含
難
基
構成
解
積分項
含
仮定
確率微分方
(Øksendal, 1998)。確
(Hakoyama and Iwasa, 2000)。
対
、近似尤度
、一般
、解
推定 可能 場合
、最尤推定量
。
場合
、環境変動
得
線形化
、
推定
非線形確率微分方程式
。例
周
欠測値
強解 陽 得
線形方程式
難
、準平衡個体数
採用
、
得
、仮
持
、連続時間
常
最尤推定量 構成
、陽
情報
間隔
最尤推定
時間間隔 観察
確率微分方程式
欠測
、環境変動
推定法
推定値
得
構成
小
仮定
。結論
自体 難
置
、連続
。
参考文献
Dempster, A. P., N. M. Laird, and D. B. Rubin (1977) “Maximum likelihood from incomplete data via
the EM algorithm,” Journal of the Royal Statistical Society Series B-Methodological, Vol. 39, No. 1,
pp. 1–38.
Geman, S and D Geman (1984) “Stochastic relaxation, Gibbs distributions, and the Bayesian restoration
of images,” Ieee Transactions On Pattern Analysis and Machine Intelligence, Vol. 6, No. 6, pp. 721–741.
Goel, N. S. and N. Richter-Dyn (1974) Stochastic models in biology: Academic Press.
Hakoyama, H. and Y. Iwasa (2000) “Extinction risk of a density-dependent population estimated from a
time series of population size,” Journal of theoretical Biology, Vol. 204, No. 3, pp. 337–359.
Knuth, DE (1997) The Art of Computer Programming, Vol. 2: Seminumerical Algorithms: AddisonWesley Professional, 3rd edition.
Louis, T. A. (1982) “Finding the Observed Information Matrix when Using the EM Algorithm,” Journal
of the Royal Statistical Society Series B-Methodological, Vol. 44, No. 2, pp. 226–233.
McLachlan, G. J. and T. Krishnan (2007) The EM algorithm and extensions, Vol. 382: Wiley-Interscience.
Metropolis, N., A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller (1953) “Equation of
state calculations by fast computing machines,” The journal of chemical physics, Vol. 21, p. 1087.
Øksendal, B. (1998) Stochastic differential equations: Springer.
Wei, G. C. G. and M. A. Tanner (1990) “A Monte Carlo implementation of the EM algorithm and the
poor man’s data augmentation algorithms,” Journal of the American Statistical Association, Vol. 85,
No. 411, pp. 699–704, September.
Wu, C. F. J. (1983) “On the convergence properties of the EM algorithm,” Annals of Statistics, Vol. 11,
No. 1, pp. 95–103.
6
© Copyright 2024 ExpyDoc
