3
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 3 years have passed since last update.

基礎から分かる時系列分析 2.カルマンフィルタ (数理的な導出)

Last updated at Posted at 2021-02-21

2.1カルマンフィルタの適用条件

シミュレーションと観測演算子が共に行列で与えられた状態空間モデル

\begin{align}
{\bf x}_t=F_t{\bf x}_{t-1}+{\bf v}_t,\ \ \ \ \ {\bf v}_t\sim \mathcal{N}({\bf 0},Q_t)\\
{\bf y}_t=H_t{\bf x}_t+{\bf w}_t,\ \ \ \ \ \ {\bf w}_t\sim \mathcal{N}({\bf 0},R_t)
\end{align}

において、全ての確率変数が正規分布に従う場合(線形ガウス状態空間モデル)、状態を逐次推定するカルマンフィルタを適用することができる。

2.2一期先予測

昨日までのデータから今日の状態を推定することが一期先予測である。これは

\begin{align}
p({\bf x}_t|{\bf y}_{1:t-1})&=\int p({\bf x}_t,{\bf x}_{t-1}|{\bf y}_{1:t-1})d{\bf x}_{t-1}\\
&=\int p({\bf x}_{t-1}|{\bf y}_{1:t-1})p({\bf x}_t|{\bf x}_{t-1})d{\bf x}_{t-1}\tag{1}
\end{align}

と書ける。(1)に対応するアルゴリズムを導出したい。今、k期までのデータが与えられた時のj期の状態が従う分布を

p({\bf x}_j|{\bf y}_{1:k})\equiv\mathcal{N}({\bf x}_{j|k},V_{j|k})

と表記する。すると(1)のうち

p({\bf x}_{t-1}|{\bf y}_{1:t-1})=\mathcal{N}({\bf x}_{t-1|t-1},V_{t-1|t-1})

は既知(一期前にフィルタ分布によって導出ずみ)であるから、

p({\bf x}_t|{\bf x}_{t-1})

を求めれば良い。これを${\bf v}_t$に関する周辺化積分形に展開すると、

p({\bf x}_t|{\bf x}_{t-1})=\int p({\bf x}_t|{\bf x}_{t-1},{\bf v}_t)p({\bf v}_t)d{\bf v}_t

システムモデルより、x_{t-1}とv_tの実現値が与えられればx_tは唯一の確定値$F_t{\bf x}_{t-1}+{\bf v}_t$を取るので、

p({\bf x}_t|{\bf x}_{t-1},{\bf v}_t)=\delta ({\bf x}_t-F_t{\bf x}_{t-1}-{\bf v}_t)

今、正規分布の分散共分散行列を$\epsilon^2 I$とすると、

p({\bf x}_t|{\bf x}_{t-1},{\bf v}_t)=\mathcal{N}(F_t{\bf x}_{t-1}+{\bf v}_t,\epsilon^2I),\ \ \ \ \ as\ \ \ \epsilon\rightarrow 0

すなわち

p({\bf x}_t-F_t{\bf x}_{t-1}|{\bf v}_t)=\mathcal{N}({\bf v}_t,\epsilon^2I)

これを周辺化することができれば、x_{t-1}を定数扱いして線形変換することでp(x_t|x_{t-1})も求まる。

(2)をよく見ると、条件付き分布の期待値が条件付けている変数${\bf v}_t$に依存している。そこで補題Aを用いれば

p({\bf x}_t-F_t{\bf x}_{t-1})=\mathcal{N}({\bf 0},Q_t+\epsilon^2I)

すなわち

p({\bf x}_t|{\bf x}_{t-1})=\mathcal{N}(F_t{\bf x}_{t-1},Q_t+\epsilon^2I)

さらに補題Aを

\begin{align}
p({\bf x}_{t-1}|{\bf y}_{1:t-1})&=\mathcal{N}({\bf x}_{t-1|t-1},V_{t-1|t-1})\\
p({\bf x}_t|{\bf x}_{t-1})&=\mathcal{N}(F_t{\bf x}_{t-1},Q_t+\epsilon^2I)
\end{align}

に適用すれば、

p({\bf x}_t|{\bf y}_{1:t-1})=\mathcal{N}(F_t{\bf x}_{t-1|t-1},F_tV_{t-1|t-1}\ ^TF_t+Q_t)

となる。

2.2 フィルタ

$t$期までのデータから$t$期の状態を推定するのがフィルタ分布である。これは

\begin{align}
p({\bf x}_t|{\bf y}_{1:t})&=\frac{p({\bf x}_t,{\bf y}_{1:t})}{p({\bf y}_{1:t})}\\
&= \frac{p({\bf x}_t,{\bf y}_{1:t})}{\int p({\bf x}_t,{\bf y}_{1:t})d{\bf x}_t}\\
&=\frac{p({\bf y}_t|{\bf x}_t)p({\bf x}_t|{\bf y}_{1:t-1})}{\int p({\bf y}_t|{\bf x}_t)p({\bf x}_t|{\bf y}_{1:t-1})d{\bf x}_t}\tag{2}
\end{align}

と書ける。(2)に対応するアルゴリズムを導出したい。観測モデルより一時点尤度は

\begin{align}
p({\bf y}_t|{\bf x}_t)=\mathcal{N}(H_t{\bf x}_t,R_t)
\end{align}

である。また

p({\bf x}_t|{\bf y}_{1:t-1})=\mathcal{N}({\bf x}_{t|t-1},V_{t|t-1})

は一期先予測より既知なので、補題Bを適用すれば

p({\bf x}_t|{\bf y}_{1:t})=\mathcal{N}({\bf x}_{t|t},V_{t|t})

が得られ、

\begin{align}
{\bf x}_{t|t}&={\bf x}_{t|t-1}+K_t({\bf y}_t-H_t{\bf x}_{t|t-1})\nonumber\\
V_{t|t}&=V_{t|t-1}K_tH_tV_{t|t-1}\nonumber\\
K_t&=V_{t|t-1}\ ^TH_t(H_tV_{t|t-1}\ ^TH_t+R_t)^{-1}
\end{align}

となる。$K_t$をカルマンゲインと呼ぶ。

補題A,B

補題A:

p({\bf z})=\mathcal{N}({\bf \mu}_{{\bf z}},\Sigma_{{\bf z}}),\ p({\bf x}|{\bf z})=\mathcal{N}(F{\bf z},\Sigma_{{\bf x}|{\bf z}})

ならば、

p({\bf x})=\mathcal{N}(F{\bf \mu}_{{\bf z}},F\Sigma_{{\bf z}}\ ^TF+\Sigma_{{\bf x}|{\bf z}})

補題B:

p({\bf z})=\mathcal{N}({\bf \mu}_{{\bf z}},\Sigma_{\bf {z}}),\ p({\bf x}|{\bf z})=\mathcal{N}(F{\bf z},\Sigma_{{\bf x}|{\bf z}})

ならば、

\begin{align}
p({\bf z}|{\bf x})&=\mathcal{N}({\bf \mu}_{{\bf z}|{\bf x}},\Sigma_{{\bf z}|{\bf x}})\\
{\bf \mu}_{{\bf z}|{\bf x}}&={\bf \mu}_{{\bf z}}+\Sigma_{{\bf z}}\ ^TF(F\Sigma_{{\bf z}}\ ^TF+\Sigma_{{\bf x}|{\bf z}})^{-1}({\bf x}-F{\bf \mu})\\
\Sigma_{{\bf z}|{\bf x}}&=\Sigma_{{\bf z}}-\Sigma_{{\bf z}}\ ^TF(F\Sigma_{{\bf z}}\ ^TF+\Sigma_{{\bf x}|{\bf z}})^{-1}F\Sigma_{{\bf z}}\nonumber
\end{align}

(証明)

補題A、Bを証明する。

\begin{eqnarray}
p({\bf z})=\mathcal{N}({\bf \mu}_{{\bf z}},\Lambda_{{\bf z}}^{-1}),\ p({\bf x}|{\bf z})=\mathcal{N}(F{\bf z},\Lambda_{{\bf x}|{\bf z}}^{-1})\ \ \ \ \ \Lambda_{{\bf z}}=\Sigma_{{\bf z}}^{-1},\ \Lambda_{{\bf x}|{\bf z}}=\Sigma_{{\bf x}|{\bf z}}^{-1}
\end{eqnarray}

なので、まずは

{\bf y}=\begin{pmatrix}{\bf z}\\{\bf x}\end{pmatrix}

と置いて、z,xの同時分布$p({\bf y})$を求める。

\begin{align}
\log p({\bf y})&=\ ^T({\bf z}-{\bf \mu}_{{\bf z}})\Lambda_{{\bf z}}({\bf z}-{\bf \mu}_{{\bf z}})+\ ^T({\bf x}-F{\bf z})\Lambda_{{\bf x}|{\bf z}}({\bf x}-F{\bf z})+const\\
&&\\
&=\ ^T{\bf z}\Lambda_{{\bf z}}{\bf z}+\ ^T{\bf z}\ ^TF\Lambda_{\bf {x}|{\bf z}}F{\bf z}-\ ^T{\bf z}\ ^TF\Lambda_{{\bf x}|{\bf z}}{\bf x}-\ ^T{\bf x}\Lambda_{\bf {x}|{\bf z}}F{\bf z}+\ ^T{\bf x}\Lambda_{{\bf x}|{\bf z}}{\bf x}\\
&\ \ -\ ^T{\bf z}\Lambda_{{\bf z}}{\bf \mu}_{{\bf z}}-\ ^T{\bf \mu}_{{\bf z}}\Lambda_{{\bf z}}{\bf z}+const\\
&&\\
&=\begin{pmatrix}{\bf z} & {\bf x}\end{pmatrix}\begin{pmatrix}
\Lambda_{{\bf z}}+\ ^TF\Lambda_{{\bf x}|{\bf z}}F & -\ ^TF\Lambda_{{\bf x}|{\bf z}}\\
-\Lambda_{{\bf x}|{\bf z}}F &\Lambda_{{\bf x}|{\bf z}}
\end{pmatrix}\begin{pmatrix}{\bf z} \\ {\bf x}\end{pmatrix}\\
&\ \ -\begin{pmatrix}{\bf z}&{\bf x}\end{pmatrix}\begin{pmatrix}\Lambda_{{\bf z}}{\bf \mu}_{{\bf z}}\\ {\bf 0}\end{pmatrix}-\begin{pmatrix}\ ^T{\bf \mu}_{{\bf z}}\Lambda &{\bf 0}\end{pmatrix}\begin{pmatrix}{\bf z}\\{\bf 0}\end{pmatrix}+const\\
&&\\
&=\ ^T{\bf y}\Lambda_{{\bf y}}{\bf y}-\ ^T{\bf y}A-\ ^TA{\bf y}+const\\
&&\\
&=\ ^T({\bf y}-\Lambda_{{\bf y}}^{-1}A)\Lambda_{{\bf y}}({\bf y}-\Lambda_{{\bf y}}A)+const
\end{align}

ただし、

\begin{align}
\Lambda_{{\bf y}}=\begin{pmatrix}
\Lambda_{{\bf z}}+\ ^TF\Lambda_{{\bf x}|{\bf z}}F & -\ ^TF\Lambda_{{\bf x}|{\bf z}}\\
-\Lambda_{{\bf x}|{\bf z}}F &\Lambda_{{\bf x}|{\bf z}}
\end{pmatrix},\ \ A=\begin{pmatrix}\Lambda_{{\bf z}}{\bf \mu}_{{\bf z}}\\ {\bf 0}\end{pmatrix}
\end{align}

と置いた($\Lambda_{{\bf y}}$は同時分布の精度行列)。したがって、同時分布の分散共分散行列と平均ベクトルは、

\begin{align}
\Sigma_{{\bf y}}&=\Lambda_{{\bf y}}^{-1}=\begin{pmatrix}
\Sigma_{{\bf z}} & \Sigma_{{\bf z}}\ ^TF\\
F\Sigma_{{\bf z}} & \Sigma_{{\bf x}|{\bf z}}+F\Sigma_{{\bf z}}\ ^TF
\end{pmatrix}\\
{\bf \mu}_{{\bf y}}&=\Lambda_{{\bf y}}^{-1}A=\begin{pmatrix}{\bf \mu}_{{\bf z}}\\ F{\bf \mu}_{{\bf z}}\end{pmatrix}
\end{align}

補題Cにより、${\bf x}$の周辺分布は

p({\bf x})=\mathcal{N}({\bf \mu}_{{\bf z}},\Sigma_{{\bf x}|{\bf z}}+F\Sigma_{{\bf z}}\ ^TF)

なので、補題Aは示された。

また補題Cにより、xが与えられた時zは正規分布に従い、その分散は

\begin{align}
\Sigma_{{\bf z}|{\bf x}}=(\Lambda_{{\bf z}}+\ ^TF\Lambda_{{\bf x}|{\bf z}}F)^{-1}=\Sigma_{{\bf z}}-\Sigma_{{\bf z}}\ ^TF(F\Sigma_{{\bf z}}\ ^TF+\Sigma_{{\bf x}|{\bf z}})^{-1}F\Sigma_{{\bf z}}
\end{align}

平均は

\begin{align}
{\bf \mu}_{{\bf z}|{\bf x}}={\bf \mu}_{{\bf z}}+\{\Sigma_{{\bf z}}-\Sigma_{{\bf z}}\ ^TF(F\Sigma_{{\bf z}}\ ^TF+\Sigma_{{\bf x}|{\bf z}})^{-1}F\Sigma_{{\bf z}}\}\ ^TF\Sigma_{{\bf x}|{\bf z}}^{-1}({\bf x}-F{\bf \mu}_{{\bf z}})
\end{align}

ここで、

\begin{align}
&\ \ \Sigma_{{\bf z}}-\Sigma_{{\bf z}}\ ^TF(F\Sigma_{{\bf z}}\ ^TF+\Sigma_{{\bf x}|{\bf z}})^{-1}F\Sigma_{{\bf z}}\ ^TF\Sigma_{{\bf x}|{\bf z}}^{-1}\\
&= \Sigma_{{\bf z}}\ ^TF\{\Sigma_{{\bf x}|{\bf z}}^{-1}+(F\Sigma_{{\bf z}}\ ^TF+\Sigma_{{\bf x}|{\bf z}})^{-1}F\Sigma_{{\bf z}}\ ^TF\Sigma_{{\bf x}|{\bf z}}\}\\
&= \Sigma_{{\bf z}}\ ^TF(F\Sigma_{{\bf z}}\ ^TF+\Sigma_{{\bf x}|{\bf z}})^{-1}\{(F\Sigma_{{\bf z}}\ ^TF+\Sigma_{{\bf x}|{\bf z}})\Sigma_{{\bf x}|{\bf z}}^{-1}-F\Sigma_{{\bf z}}\ ^TF\Sigma_{{\bf x}|{\bf z}}^{-1}\}\\
&=\Sigma_{{\bf z}}\ ^TF(F\Sigma_{{\bf z}}\ ^TF+\Sigma_{{\bf x}|{\bf z}})^{-1}
\end{align}

なので、

\begin{align}
{\bf \mu}_{{\bf z}|{\bf x}}={\bf \mu}_{{\bf z}}+\Sigma_{{\bf z}}\ ^TF(F\Sigma_{{\bf z}}\ ^TF+\Sigma_{{\bf x}|{\bf z}})^{-1}({\bf x}-F{\bf \mu})
\end{align}

以上より、補題Bが示された。

補題C

x_a,x_bの同時分布が

\begin{align}
p(\begin{pmatrix}{\bf x}_a\\ {\bf x}_b\end{pmatrix})&=\mathcal{N}({\bf \mu},\Sigma)\\
{\bf \mu}&=\begin{pmatrix}{\bf \mu}_a\\ {\bf \mu}_b\end{pmatrix}\\
\Sigma&=\begin{pmatrix}\Sigma_{aa} & \Sigma_{ab}\\ \Sigma_{ba} &\Sigma_{bb}\end{pmatrix}\\
\Lambda=\Sigma^{-1}&=\begin{pmatrix}\Lambda_{aa} & \Lambda_{ab}\\ \Lambda_{ba} &\Lambda_{bb}\end{pmatrix}
\end{align}

で与えられているとき、x_aの周辺分布は

\begin{align}
p({\bf x}_a)=\mathcal{N}({\bf \mu}_a,\Sigma_{aa})
\end{align}

x_b$が与えられた時のx_aの条件付き分布は

\begin{align}
p({\bf x}_a|{\bf x}_b)&=\mathcal{N}({\bf \mu}_{a|b},\Lambda_{aa}^{-1})\\
{\bf \mu}_{a|b}&={\bf \mu}_a-\Lambda_{aa}^{-1}\Lambda_{ab}({\bf x}_b-{\bf \mu}_b)
\end{align}

(証明)

補題Cを証明する。まずは周辺分布p(x_a)を求める。

\begin{align}
p({\bf x}_a)&=\int p({\bf x}_a,{\bf x}_b)d{\bf x}_b\\
&\propto\int\exp\{-\frac{1}{2}\ ^T({\bf x}-{\bf \mu})\Lambda({\bf x}-{\bf \mu})\}d{\bf x}_b
\end{align}

指数関数の中身をx_bで平方完成する。

\begin{align}
&\ \ ^T({\bf x}-{\bf \mu})\Lambda({\bf x}-{\bf \mu})\\
&&\\
&=\ ^T({\bf x}_a-{\bf \mu}_a)\Lambda_{aa}({\bf x}_a-{\bf \mu}_a)+\ ^T({\bf x}_a-{\bf \mu}_b)\Lambda_{ab}({\bf x}_b-{\bf \mu}_b)\\
&\ \ +\ ^T({\bf x}_b-{\bf \mu}_b)\Lambda_{ba}({\bf x}_a-{\bf \mu}_a)+\ ^T({\bf x}_b-{\bf \mu}_b)\Lambda_{bb}({\bf x}_b-{\bf \mu}_b)\\
&&\\
&=\ ^T{\bf x}_b\Lambda_{bb}{\bf x}_b+\ ^T{\bf x}_b\Lambda_{ba}({\bf x}_a-{\bf \mu}_a)\\
&\ -\ ^T{\bf x}_b\Lambda_{bb}{\bf \mu}_b-\ ^T{\bf \mu}_b\Lambda_{bb}{\bf x}_b+\ ^T({\bf x}_a-{\bf \mu}_a)\Lambda_{ab}{\bf x}_b\\
&\ \ +\ ^T{\bf x}_a\Lambda_{aa}{\bf x}_a-\ ^T{\bf \mu}_a\Lambda_{aa}{\bf x}_a-\ ^T{\bf x}_a\Lambda_{aa}{\bf \mu}_a-\ ^T{\bf x}_a\Lambda_{ab}{\bf \mu}_b-\ ^T{\bf \mu}_b\Lambda_{ba}{\bf x}_a+const\\
&&\\
&=\ ^T{\bf x}_b\Lambda_{bb}{\bf x}_b-\ ^T{\bf x}_b\{\Lambda_{bb}{\bf \mu}_b-\Lambda_{ba}({\bf x}_a-{\bf \mu}_a)\}-\{\ ^T{\bf \mu}_b\Lambda_{bb}-\ ^T({\bf x}_a-{\bf \mu}_a)\Lambda_{ab}\}{\bf x}_b\\
&\ \ +\ ^T{\bf x}_a\Lambda_{aa}{\bf x}_a-\ ^T{\bf x}(\Lambda_{aa}{\bf \mu}_a+\Lambda_{ab}{\bf \mu}_b)-(\ ^T{\bf \mu}_a\Lambda_{aa}+\ ^T{\bf \mu}_b\Lambda_{ba}){\bf x}_a+const\\
&&\\
&=\ \ ^T{\bf x}_bA^{-1}{\bf x}_b-\ ^T{\bf x}_bB-\ ^TB{\bf x}_b\\
&\ \ +\ ^T{\bf x}_a\Lambda_{aa}{\bf x}_a-\ ^T{\bf x}(\Lambda_{aa}{\bf \mu}_a+\Lambda_{ab}{\bf \mu}_b)-(\ ^T{\bf \mu}_a\Lambda_{aa}+\ ^T{\bf \mu}_b\Lambda_{ba}){\bf x}_a+const\\
&&\\
&=\ ^T({\bf x}_b-AB)A^{-1}({\bf x}-AB)-\ ^TB\ ^TAB\\
&\ \ +\ ^T{\bf x}_a\Lambda_{aa}{\bf x}_a-\ ^T{\bf x}(\Lambda_{aa}{\bf \mu}_a+\Lambda_{ab}{\bf \mu}_b)-(\ ^T{\bf \mu}_a\Lambda_{aa}+\ ^T{\bf \mu}_b\Lambda_{ba}){\bf x}_a+const
\end{align}

ただし、

\begin{align}
A=\Lambda_{bb}^{-1},\ \ B=\Lambda_{bb}{\bf \mu}_b-\Lambda_{ba}({\bf x}_a-{\bf \mu}_a)
\end{align}

と置いた。ゆえ、

\begin{align}
p({\bf x}_a)&\propto\exp[\\
&\ \ -\frac{1}{2}\{-\ ^TB\ ^TAB+\ ^T{\bf x}_a\Lambda_{aa}{\bf x}_a-\ ^T{\bf x}(\Lambda_{aa}{\bf \mu}_a+\Lambda_{ab}{\bf \mu}_b)-(\ ^T{\bf \mu}_a\Lambda_{aa}+\ ^T{\bf \mu}_b\Lambda_{ba}){\bf x}_a\}]
\end{align}

右辺の指数部分の中身を変形してまとめると、

\begin{align}
&\ \ ^T{\bf x}_a(-\ ^T\Lambda_{ba}\Lambda_{bb}^{-1}\Lambda_{ba}+\Lambda_{aa}){\bf x}_a\\
&\ \ -\ ^T{\bf x}_a(-\ ^T\Lambda_{ba}{\bf \mu}_b-\ ^T\Lambda_{ba}\Lambda_{bb}^{-1}\Lambda_{ba}{\bf \mu}_a+\Lambda_{aa}{\bf \mu}_a+\Lambda_{ab}{\bf \mu}_b)\\
&\ \ -(\ ^T{\bf \mu}_b\Lambda_{ba}{\bf x}_a-\ ^T{\bf \mu}_a\ ^T\Lambda_{ba}\Lambda_{bb}^{-1}\Lambda_{ba}+\ ^T{\bf \mu}_a\Lambda_{aa}+\ ^T{\bf \mu}_b\Lambda_{ab}){\bf x}_a\\
&=\ ^T{\bf x}_aB^{-1}{\bf x}_a-\ ^T{\bf x}_aC-\ ^C{\bf x}_a\\
&=\ ^T({\bf x}_a-BC)B^{-1}({\bf x}_a-BC)+const
\end{align}

ただし、

\begin{align}
B&=(\Lambda_{aa}-\Lambda_{ab}\Lambda_{bb}^{-1}\Lambda_{ba})^{-1}\\
C&=(\Lambda_{aa}-\Lambda_{ab}\Lambda_{bb}^{-1}\Lambda_{ba}){\bf \mu}_a
\end{align}

したがって、

\begin{align}
{\bf \mu}_a&=BC={\bf \mu}_a\\
\Sigma_a&=B=\Sigma_{aa}
\end{align}

つぎに条件付き分布p(x_a|x_b)を求める。

\begin{align}
^T({\bf x}-{\bf \mu})\Sigma^{-1}({\bf x}-{\bf \mu})&=
^T\begin{pmatrix} {\bf x}_a-{\bf \mu}_a \\ {\bf x}_b-{\bf \mu}_b\end{pmatrix}
\begin{pmatrix}\Lambda_{aa} & \Lambda_{ab}\\ \Lambda_{ba} & \Lambda_{bb}\end{pmatrix}
\begin{pmatrix}{\bf x}_a-{\bf \mu}_a\\ {\bf x}_b-{\bf \mu}_b\end{pmatrix}\\
&&\\
&=\ ^T({\bf x}_a-{\bf \mu}_a)\Lambda_{aa}({\bf x}_a-{\bf \mu}_a)+\ ^T({\bf x}_a-{\bf \mu}_b)\Lambda_{ab}({\bf x}_b-{\bf \mu}_b)\\
&\ \ +\ ^T({\bf x}_b-{\bf \mu}_b)\Lambda_{ba}({\bf x}_a-{\bf \mu}_a)+\ ^T({\bf x}_b-{\bf \mu}_b)\Lambda_{bb}({\bf x}_b-{\bf \mu}_b)\\
&&\\
&=\ ^T{\bf x}_a\Lambda_{aa}{\bf x}_a-\ ^T{\bf x}_a\Lambda_{aa}{\bf \mu}_a+\ ^T{\bf x}_a\Lambda_{ab}({\bf x}_b-{\bf \mu}_b)\\
&\ \ -\ ^T{\bf \mu}_a\Lambda_{aa}{\bf x}_a+\ ^T({\bf x}_b-{\bf \mu}_b)\Lambda_{ba}{\bf x}_a+const\\
&&\\
&=\ ^T{\bf x}_a\Lambda_{aa}{\bf x}_a-\ ^T{\bf x}\{\Lambda_{aa}{\bf \mu}_a-\Lambda_{ab}({\bf x}_b-{\bf \mu}_b)\}\\
&\ \ -\{\ ^T{\bf \mu}_a\Lambda_{aa}-\ ^T({\bf x}_b-{\bf \mu}_b)\Lambda_{ba}\}{\bf x}_a+const\\
&&\\
&=\ ^T{\bf x}_aA^{-1}{\bf x}_a-\ ^T{\bf x}_aB-\ ^TB{\bf x}_a+const\\
&&\\
&=\ ^T({\bf x}_a-AB)A^{-1}({\bf x}_a-AB)+const
\end{align}

ただし、

\begin{align}
A=\Lambda_{aa}^{-1},\ B=\Lambda_{aa}{\bf \mu}_a-\Lambda_{ab}({\bf x}_b-{\bf \mu}_b)
\end{align}

と置いた。したがって、条件付き分布p(x_a|x_b)の平均と分散は、

\begin{align}
{\bf \mu_{a|b}}&=&\Lambda_{aa}^{-1}\{\Lambda_{aa}{\bf \mu}_a-\Lambda_{ab}({\bf x}_b-{\bf \mu}_b)\}\\
&={\bf \mu}_a-\Lambda_{aa}^{-1}\Lambda_{ab}({\bf x}_b-{\bf \mu}_b)\\
\Sigma_{a|b}&=\Lambda_{aa}^{-1}
\end{align}

となる。(証明終わり)

つぎに$\Lambda_{aa}^{-1},\Lambda_{ab}$を求める(このままの方が簡潔なので普通は求めない)。今、

\begin{align}
\begin{pmatrix}\Sigma_{aa} &\Sigma_{ab}\\ \Sigma_{ba} & \Sigma_{bb}\end{pmatrix}^{-1}=\begin{pmatrix}\Lambda_{aa} &\Lambda_{ab}\\ \Lambda_{ba} & \Lambda_{bb}\end{pmatrix}
\end{align}

であるから、

\begin{align}
\Lambda_{aa}&=(\Sigma_{aa}-\Sigma_{ab}\Sigma_{bb}^{-1}\Sigma_{ba})^{-1}\\
\Lambda_{ab}&=-(\Sigma_{aa}-\Sigma_{ab}\Sigma_{bb}^{-1}\Sigma_{ba})^{-1}\Sigma_{ab}\Sigma_{bb}^{-1}
\end{align}

より、

\begin{align}
{\bf \mu}_{a|b}&={\bf \mu}_a+\Sigma_{ab}\Sigma_{bb}^{-1}({\bf x}_b-{\bf \mu}_b)\\
\Sigma_{a|b}&=\Sigma_{aa}-\Sigma_{ab}\Sigma_{bb}^{-1}\Sigma_{ba}
\end{align}

となった。条件付き分布の平均はx_bの線形関数。

2.3 尤度

尤度$p({\bf y}_{1:T})$は

\begin{align}
p({\bf y}_{1:T})=\Pi_{t=1}^Tp({\bf y}_{t}|{\bf y}_{1:t-1})
\end{align}

ここで、

\begin{eqnarray}
p({\bf x}_t|{\bf y}_{1:t-1})&=\mathcal{N}({\bf x}_{t|t-1}|V_{t|t-1})\\
p({\bf y}_t|{\bf x}_t)&=\mathcal{N}(H_t{\bf x}_t,R_t)
\end{eqnarray}

を補題Aに適用して、

\begin{align}
p({\bf y}_t|{\bf y}_{1:t-1})=\mathcal{N}(H_t{\bf x}_{t|t-1},H_tV_{t|t-1}\ ^TH_t+R_t)
\end{align}

となる。したがって、y_tの次元数をmと置くと対数尤度は

\begin{align}
\log p({\bf y}_{1:T})&=\sum_{t=1}^T\log p({\bf y}_t|{\bf y}_{1:t-1})\\
&=-\frac{1}{2}\sum_{t=1}^T\log|D_{t|t-1}|-\frac{Tm}{2}\log(2\pi)\\
&\ \ -\frac{1}{2}\sum_{t=1}^T\ ^T({\bf y}_t-H_t{\bf x}_{t|t-1})D_{t|t-1}^{-1}({\bf y}_t-H_t{\bf x}_{t|t-1})
\end{align}

ただし、

\begin{align}
D_{t|t-1}=H_tV_{t|t-1}\ ^TH_t+R_t
\end{align}

と置いた。パラメータの次元数を$M$と置くと、AICは

\begin{align}
AIC=\log p({\bf y}_{1:T})-M
\end{align}

なので、これをもとにモデル選択を行う。

次回:基礎から分かる時系列分析 3.カルマンフィルタ(ローカルレベルモデルの実装)

3
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
3
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?