カルマンフィルタについて勉強しているのですが，アウトプットの一環で備忘録としてQiitaに投稿します．
理解が不十分なところ，導出の流れが間違えているところがあると思いますので，指摘，補足等をコメント等でしてくれると嬉しいです．
使用した本は統計学 One Point 2 カルマンフィルタです．
#第1章 確率分布と時系列に関する準備事項
##確率ベクトルと期待値
ここで$x$は$m$次元確率ベクトル$x=(x_1,x_2,\cdots,x_{m-1}, x_m)^T$，$y$は$n$次元確率ベクトル$y=(y_1,y_2,\cdots,y_{n-1}, y_n)^T$である．この時，期待値と共分散行列を以下のように定義する．
E[f(x)]=\int f(x)p(x)dx
Cov(x,y)=E\{[x-E(x)][y-E(y)]^T\}\\
また，$Var(x)=Cov(x,x)$を分散共分散行列と呼ぶ．
##確率ベクトルの周辺分布と条件付き分布
2つの確率ベクトルの同時密度関数$p(x,y)$が与えられたとき，$x$単独の確率を周辺分布と呼び，その確率密度関数を$p(x)$で表し，周辺密度関数と呼ぶ．また，$y$が得られた下で$X$が従う確率分布を，$y$が与えられた下での$x$の条件付き分布と呼び，その確率密度関数を$p(x\mid y)$で表し，条件付き密度関数と呼ぶ．
$x$と$y$の周辺密度関数$p(x),p(y)$及び，$y$が与えられた下での$x$の条件付き関数$p(x\mid y)$は以下のように得られる．
p(x)=\int p(x,y)dy,\quad p(y)=\int p(x,y)dx
p(x\mid y)=\frac{p(x,y)}{p(y)}
p(x,y)=p(x\mid y)p(y)=p(y\mid x)p(x)
p(y\mid x)=\frac{p(x\mid y)p(y)}{p(x)}=\frac{p(x\mid y)p(y)}{\int p(x\mid y)p(y)dy}\tag{1}
が得られる．$(1)$式がベイズの定理と呼ばれる．
これは$y$から$x$への因果関係$p(x\mid y)$を利用し，逆に結果$x$から原因$y$が推定できるため強力である．
条件付き分布に基づき，$y$が与えられた下での$x$の実数値関数$f(x)$の**条件付き期待値$E[f(x)\mid y]$を$x$の条件付き密度関数$p(x\mid y)$に対する期待値として以下のように定義する．
E[f(x)\mid y]=\int f(x)p(x\mid y)dx
特に，$y$が与えられた下での$x$の条件付き分布における平均ベクトル$E(x\mid y)$を条件付き平均，分散行列$Var(x\mid y)=E\{[x-E(x)][x-E(x)]^T \mid y\} $を条件付き分散と呼び，$y$が与えられた下での2つの確率ベクトル$x_1,x_2$の条件付き共分散，$Cov(x_1,x_2\mid y)$も同様に定義される．
また，条件付き平均，条件付き分散を用いて，$x$の周辺分布における平均ベクトル$E(x)$，分散行列$Var(x)$を以下のように表せる．
E(x)=E[E(x\mid y)],\quad Var(x)=E[Var(x\mid y)]+Var[E(x\mid y)]
##多変量正規分布の定義と基本的性質
確立ベクトル$x=(x_1,x_2,\cdots,x_{m-1}, x_m)^T$が以下の同時密度関数をもつとき，$x$の同時分布を平均ベクトル$\mu=(\mu_1,\mu_2,\cdots,\mu_{m-1},\mu_m)^T$，分散共分散行列$\Sigma$を持つ多変量正規分布と呼び，記号では$N(\mu,\Sigma)$と表す．
p(x)=(2\pi)^{-m/2}(det\Sigma)^{-1/2}\exp]\left\{-\frac{1}{2}(x-\mu)^T\Sigma^{-1}(x-\mu)\right\}
ここで，$det\Sigma$は$\Sigma$の行列式である．
多変量正規分布における周辺分布と独立性，条件付き分布について示す．
$m$変量正規確率ベクトル$x=(x_1,x_2,\cdots,x_{m-1}, x_m)^T$を先頭から$k$個の成分$y_1$と後ろから$m-k$個の成分$y_2$に分け，$x$の平均ベクトル$\mu$と分散共分散行列$\Sigma$も同じように分けると次のように表すことができる．
x=
\begin{pmatrix}
y_1\\
y_2
\end{pmatrix},\quad
\mu=
\begin{pmatrix}
\mu_1\\
\mu_2
\end{pmatrix},\quad
\Sigma=
\begin{pmatrix}
\Sigma_{11}&\Sigma_{12}\\
\Sigma_{21}&\Sigma_{22}
\end{pmatrix}
この時，平均は$\mu_1=E(y_1),\mu_2=E(y_2)$，分散行列は$\Sigma_{11}=Var(y_1),\Sigma_{12}=\Sigma_{21}^T=Cov(y_1,y_2),\Sigma_{22}=Var(y_2)$となる．
$m$次元単位行列$I_m$の先頭$k$行をとった部分行列を$A$と置くと$y_1=Ax$なので，$y_1$の周辺分布は
y_1\sim N(\mu_1,\Sigma_1)
となり，同様に$y_2$の周辺分布は$y_2\sim (\mu_2,\Sigma_2)$となる．
$y_1$と$y_2$が互いに独立の必要十分条件は$Cov(y_1,y_2)=\Sigma_{12}=O$である．
$y_2$が与えられた下での$y_1$の条件付き分布は多変量正規分布に従い，その条件付き平均，条件付き分散は次式で与えられる．
E(y_1\mid y_2)=\mu_1+\Sigma_{12}\Sigma_{22}^{-1}(y_2-\mu_2)
Var(y_1\mid y_2)=\Sigma_{11}-\Sigma_{12}\Sigma_{22}\Sigma_{12}^T
##定常性とコレログラム
簡単のために単変量時系列を扱う．
時系列$y_1,\cdots,y_n$について，平均と分散を以下のようにおく．
\mu_t=E(y_t),\quad \sigma_t^2=Var(y_t)=E[(y_t-\mu_t)^2],\quad t=1,\cdots,n
同一時系列内の共分散を
Cov(y_s,y_t)=E[(y_s-\mu_s)(y_t-\mu_t)],\quad s=1,\cdots,n,\quad t=1,\cdots,n
と表し，これを自己共分散と呼ぶ．
任意の整数$k$に対して
\mu_{t+k}=\mu_t,\quad \sigma_{t+k}^2=\sigma_t^2,\quad Cov(y_s,y_t)=Cov(y_{s+k},y_{t+k})
であるとき，時系列は弱定常である．さらに，時系列の同時密度関数$f$について任意の時点集合$\{t_1,\cdots,t_m\}$と任意の整数$k$に対して
f(y_{t_1+k},\cdots,y_{t_m+k})=f(y_{t_1},\cdots,y_{t_m})
をみたすとき，その時系列は強定常である．定常性をみたさない時系列は非定常である．
弱定常をみたす時系列では，任意の非負整数$k=0,1,\cdots,n$に対してラグ(時間差)$k$の自己共分$Cov(y_{t},y_{t+k})=C_k$は
時点$t$によらずラグ$k$にのみ依存する関数となり，この$C_k$を自己共分散関数と呼ぶ．
弱定常の簡単な例として自己共分散$Cov(\eta_s,\eta_t)$が$s=t$のとき以外ゼロとなる時系列$\eta_1,\cdots,\eta_n$であり，ホワイトノイズと呼ばれる．
非定常時系列な簡単な例は，ホワイトノイズお累積和として得られる時系列
y_t=\eta_1+\cdots+\eta_n,\quad t=1,\cdots,n
であり，ランダムウォークと呼ばれる．
定常時系列の観測値$y_1,\cdots,y_n$として得られたとき，そこから平均$\mu$，分散$\sigma^2$，自己共分散関数$C_k,\quad k=0,1,\cdots,n$の推定値として，標本平均$\hat{\mu}=\sum^n_{t=1}y_t/n$，標本分散$\hat{\sigma^2}=\sum^n_{t=1}(y_t-\bar{y})^2/n$，標本自己共分散
\hat{C_k}=\frac{1}{n}\sum^{n-k}_{t=1}(y_t-\bar{y})(y_{t+k}-\bar{y})
が求められる．
また，標本自己共分散を用いて，自己共分散$R_k,\quad k=0,1,\cdots,n$の推定値として**標本自己相関関数が
\hat{R}_k=\frac{\hat{C_k}}{\hat{C_0}},\quad k=0,1,\cdots,n
のように得られる．
ラグ$k$に対する自己相関関数または標本自己相関数の推移をプロットしたものをコレログラムと呼ぶ．
##自己回帰移動平均(ARMA)モデル
まず，ARモデルとMAモデルを導入する．
$\eta_1,\cdots,\eta_n$を平均ゼロ，分散$\sigma^2$のホワイトノイズとする．$t=1,\cdots,n$について
y_t=\phi_1y_{t-1}+\phi_2y_{t-2}+\cdots+\phi_py_{t-p}+\eta_t
により定義される$y_t$の時系列モデルは$p$次の自己回帰モデル(ARモデル)と呼ばれ，$AR(p)$と表記される．
このモデルは$y_t$の期待値がラグ$p$まで過去の観測値$y_{t-1},\cdots,y_{t-p}$から定まるモデルであり，その係数$\phi_1,\cdots,\phi_p$は自己回帰係数と呼ばれる．
次に，$\eta_1,\cdots,\eta_n$を平均ゼロ，分散$\sigma^2$のホワイトノイズとする．$t=1,\cdots,n$について
y_t=\eta_t+\lambda_1\eta_{t-1}+\cdots+\lambda_q\eta_{t-q}
により定義される$y_t$の時系列モデルは$q$次の移動平均モデル(MAモデル)と呼ばれ，$MA(q)$と表記される．
このモデルは観測値$y_t$が現時点からラグ$q$の過去までのノイズ$\eta_t,\eta_{t-1},\cdots,\eta_{t-q}$により定まるモデルであり，過去のノイズの影響を表す係数$\lambda_1,\cdots,\lambda_q$は移動平均係数と呼ばれる．
$AR(p)$と$MA(q)$を次式と用に合成することで，$(p,q)$次の自己回帰移動平均モデル(ARMA)を得ることができ，$ARMA(p,q)$と表す．
y_t=\phi_1y_{t-1}+\phi_2y_{t-2}+\cdots+\phi_py_{t-p}+\eta_t+\lambda_1\eta_{t-1}+\cdots+\lambda_q\eta_{t-q}
ARMAモデルの自己共分散関数を導出するために，インパル応答関数を以下のように定義する．
\delta_k=\frac{Cov(\eta_t,y_{t+k})}{Var(\eta_t)}=\frac{E(\eta_ty_{t+k})}{\sigma_{\eta}^2},\quad k=0,1,\cdots
インパルス応答関数は，ある時点$t$に加わったノイズのうち$k$期後に影響する割合を示している．
ノイズ$\eta_t$は過去の観測値とは独立であるため，負の整数$k=-1,-2,\cdots$に対するインパルス応答関数は$\delta_k=0$とする．
$ARMA(p,q)$のインパル応答関数は次の漸化式により逐次的に算出できる．(あってる自信がない)
\begin{eqnarray}
\delta_0&=&\frac{E(\eta_ty_{t+0})}{\sigma_{\eta}^2}\\
&=&\frac{E[\eta_t(\phi_1y_{t-1}+\phi_2y_{t-2}+\cdots+\phi_py_{t-p}+\eta_t+\lambda_1\eta_{t-1}+\cdots+\lambda_q\eta_{t-q})]}{\sigma_{\eta}^2}\\
&=&\frac{E(\eta_t*\eta_t)+0+0+\cdots+0}{\sigma_{\eta}^2}\\
&=&\frac{E\{[\eta_t-E(\eta_t)][\eta_t-E(\eta_t)]\}}{\sigma_{\eta}^2}\\
&=&\frac{Var(\eta_t)}{\sigma_{\eta}^2}=1
\end{eqnarray}
\begin{eqnarray}
\delta_k&=&\frac{E(\eta_ty_{t+k})}{\sigma_{\eta}^2}\\
&=&\frac{E[\eta_t(\phi_1y_{t+k-1}+\phi_2y_{t+k-2}+\cdots+\phi_py_{t+k-p}+\eta_{t+k}+\lambda_1\eta_{t+k-1}+\cdots+\lambda_k\eta_{t}+\cdots+\lambda_q\eta_{t+k-q})]}{\sigma_{\eta}^2}\\
&=&\frac{\phi_1E(\eta_ty_{t+k-1})+\phi_2E(\eta_ty_{t+k-2})+\cdots+\phi_{k-1}E(\eta_ty_{t+1})+\phi_{k}E(\eta_ty_{t})+\lambda_kE(\eta_t^2)}{\sigma_{\eta}^2}\\
&=&\phi_{k}\delta_0+\phi_{k-1}\delta_1+\cdots+\phi_1\delta_{k-1}+\lambda_k\\
&=&\sum_{j=1}^{k}\phi_j\delta_{k-j}+\lambda_k,\quad k=1,2,\cdots,n
\end{eqnarray}
ただし，$j>p$のとき$\phi_j=0$，$k>q$のとき$\lambda_k=0$とする．
得られたインパル応答関数を用いて，自己共分散関数を次の連立方程式を解くことで得られる．
まず，$y_t$に$y_{t-k}$をかけ，期待値をとる．
E(y_t,y_{t-k})=\sum_{i=1}^p\phi_pE(y_{t-i}y_{t-k})+E(\eta_ty_{t-k})+\sum_{j=1}^q\lambda_jE(\eta_{t-j}\eta_{t-k})\tag{2}
自己共分散を利用し，
E(y_{t-j}y_{t-k})=Cov(y_{t-j},y_{t-k})=C_{k-j}\tag{3}
インパルス応答関数より，
E(\eta_ty_{t-k})=\delta_{k}\sigma_{\eta}^2\tag{4}
$(3),(4)$式を$(2)$式に代入すると，($C_k$に関しては自信がない)
\begin{eqnarray}
C_k&=&E(y_t,y_{t-k})\\
&=&\sum_{i=1}^p\phi_pC_{k-j}+\delta_{k}\sigma_{\eta}^2+\sum_{j=1}^q\lambda_j\delta_{j-k}\sigma_{\eta}^2\\
\\
C_0&=&\sum_{i=1}^p\phi_pC_{0-j}+\delta_{0}\sigma_{\eta}^2+\sum_{j=1}^q\lambda_j\delta_{j-0}\sigma_{\eta}^2\\
&=&\sum_{i=1}^p\phi_pC_{j}+\sigma_{\eta}^2+\sum_{j=1}^q\lambda_j\delta_{j}\sigma_{\eta}^2\\
&=&\sum_{i=1}^p\phi_pC_{j}+\sigma_{\eta}^2\left\{1+\sum_{j=1}^q\lambda_j\delta_{j}\right\}
\end{eqnarray}
ただし，$j<0$のとき$C_j=C_{-j}$，$\delta_j=0$とする．
##自己回帰和文移動平均(ARIMA)モデル
ARMAモデルは定常な時系列に適用されるモデルである．
非定常な時系列$y_t$に対して，その階差$\Delta y_t=y_t-y_{t-1}$を時系列と見たときに定常性をみたすため，ARMAモデルが適用できる場合がある．ランダムウォークは非定常時系列であるが，その階差は
\Delta y_t=y_t-y_{t-1}=\eta_t
とホワイトノイズ$\eta_t$そのものであり定常性をみたす．
定常な時系列を得るために時系列の観測値の階差をとるプロセスは和分と呼ばれ，1階の階差で定常にならない場合はさらに階差をとることができる．近似的に定常な時系列が得られるまで$d$階の階差$\Delta^dy_t$をとった後の時系列に$ARMA(p,q)$を当てはめたものを自己回帰和分移動平均モデル(ARIMA)と呼ばれ，
$ARIMA(p,d,q)$と表記する．