経済のお勉強：カルマンフィルタ

目次 (2024/3/14 ver.0.02)

1番目のタブが選択された！

1.はじめに

もくじ

本稿は、「 第１章　線形時系列解析と非線形時系列解析の違い 」を参照しています。

1.1 時系列解析の研究の歴史

1.時系列解析の研究の歴史
  • 伝統的な時系列の解析は、Yule (1927) の自己回帰モデル (autoregressive model), 以下 AR と略す）に始まり、移動平均モ デル (moving average model, MA)、自己回帰移動平均モデ ル (autoregressive moving average model,ARMA) などのいわゆる 線形時系列解析よりなされてきた（例えば合原、2000)。
  • また、同解析による時系列 信号の基本的な特性は、自己相関関数 (autocorrelation function)、 自己相関係数 (autocorrelation coefficient)、パワースペクトル (power spectrum) などを用いてなされてきた。
  • ここでは、1.1.1 節でまず最初に線形時系列解析の概要を紹介し、1.1.2 節で一般の時系列信号 の定常性 (stationarity) とエルゴード性 (ergodicity) について述べる。
  • 最後に 1.1.3 節で非線形 時系列解析とりわけカオス時系列解析の例を示す。
2.線形時系列解析の概要
  • 時系列解析の代表的な著書としては、Box and Jenkins (1970, 1976)、Hannan and Deistler (1988)、Hamilton, J. D. (1994) などが知られている。
  • ここでは線形時系列解析モデルを幾つか の側面から最も包括的に紹介している Hannan and Deistler による線形時系列解析の概要をまとめる。
  • Hannan らは線形時系列解析の3つの表現を示している。それらは、以下の3つである。
    1. 線形系の入力ー出力表 現 (the input-output representation)、
    2. 線形系のベクトル差分方程式表現 (the vector difference equation (VDE) representation)、
    3. 線形系の状態空間表現 (the state space-representation)
  • とりわけ、（３）による線形系の表現はもともと時系列解析 用ではなく多変量解析の１つの方法として開発された共分散構造分析 (SEM) と形式的には非常に似通っており、実際最近内外の一部の研究者により脳波の解析などで SEM を 応用した線形解析が行われているが、そのような解析の理解のために役立つ。
  • ここで、それらの表現のもとになるデータは２種類の T 時点の観測時系列信号であり、それらをy(t)、およびz(t)、t=1, .., T とする。
  • これらの信号は共に複数の要素からなり、順に s 個、m 個からなるベクトルである。
3.線形系の入力ー出力表現
  • 線形系 (lenear system) とりわけ線形時不変系 (linear time-invariant system) の第１の表現は、系の入力ー出力表現 (the input-output representation) である。
  • ここでは、離散時系 (discrete-time system) のみを扱い、時点は整数 (integers) Z上にまたがるか、場合により自然数(natural numbers) N上にまたがるものとする。

    • $\Large \displaystyle y(t) = \sum_{n=-\infty}^{\infty}L(n)w(t-n)$

  • ここで、「$y (t)$」 は出力 (the outputs) 、「$w (t)$」 は入力 (the inputs)であり、 「$L(j)\in R^{s\times m}$」は次式で与えられる伝達関数 (the transfer function) の係数(行列）である。

    • $\Large \displaystyle l(t) = \sum_{n=-\infty}^{\infty}L(n)z^n, z\in G$

  • ここで、「$L(n) | n\in Z$」 は (1) 式のウエイト関数 (weight function) と呼ばれる。
  • ウエイト関数は、「$y(t)=L(t-τ)e(τ)$」 が時点 τでのインパルス、e (τ)、から生じる 出力の時間履歴であるので、インパルス応答関数 (impulse response function)、とも呼ばれる。
  • なお、ここでいう線形時不変系における時不変性 (time invariance)と は、L (j) が時点 t に依存しないか、一般的な形で表現するならば (y (t) | tZ) が (w (t) | tZ) に対応する出力ならば (y (t+s) | tZ) は、任意の s Z に 対する (w (t) | t Z) に対応する出力 である、ことを意味する。
  • また、系はもしy (t) が w (s), s
  • 上記の系の入出力表現では、このこ とはL (j)=o, j <0 であることに等しい。
4.線形系のベクトル差分方程式表現
  • ２つ目の系の表現は、ベクトル差分方程式表現 (the vector difference equation (VDE) representation) である。

    • $\Large \displaystyle \sum_{n=0}^{p}A(n)y(t-n) = \sum_{n=1}^{r}D(n)z(t-n)+\sum_{n=0}^{q}B(n)e(t-n)$

1.2 線形システム

1.線形システムの定義
  • 任意の物理的システムは、原因「$f(t)$」が入力として加えられた時、 その結果「$g(t)$」を出力または応答として一意的に生ずる変換器と考える事ができる。
  • 任意の物理的システムは、数学的に表現すると与えられた関数「$f(t)$」から関数「$g(t)$」への変換である。

    • $\Large \displaystyle T\{f(t)\} = g(t)$

  • もしシステムが線形性を有すならば、そのシステムの特性の決定は、比較的、容易である。
2.線形性
  • 線形性は、線形演算子「$L$」を使用して、以下のように表される。

    • $\Large \displaystyle L\{\alpha_1 f_1(t)+\alpha_2 f_2(t) \} = \alpha_1 L\{f_1(t)\}+\alpha_2 L\{f_2(t)\}$

  • 線形システムの例として、以下のような微分回路がある。

    • $\Large \displaystyle L\{f(t)\} = \frac{df}{dt}$

  • もう一つの線形システムの例として、以下のような遅延線路がある。

    • $\Large \displaystyle L\{f(t)\} = f(t-t_1)$

  • システムの入力と出力の関係を考える限り、線形システムと線形演算子と言う用語は、同じものである。
  • 単位インパルス「$\delta(t-t_1)$」に対する線形システムの応答を「$h(t;t_1)$」と書くことにする。

    • $\Large \displaystyle L\{\delta(t-t_1)\} = h(t;t_1)$

3.関数「$f(x)$」は、単位インパルス「$\delta(t-t_1)$」の和で表される。
  • 関数「$f(x)$」は、以下のように単位インパルス「$\delta(t-t_1)$」の和で表される。

    • $\Large \displaystyle f(x)=\int_{-\infty}^{\infty}f(t_1)\delta(t-t_1)dt_1$

  • 関数「$f(x)$」を線形は、以下のように単位インパルス「$\delta(t-t_1)$」の和で表される。

    • $\Large \displaystyle L\{f(x)\}=\int_{-\infty}^{\infty}f(t_1)L\{\delta(t-t_1)\}dt_1 =\int_{-\infty}^{\infty}f(t_1)h(t;t_1)dt_1$

  • 従って、関数「$g(t)$」は、以下のようになる。

    • $\Large \displaystyle g(t) = \int_{-\infty}^{\infty}f(t_1)h(t;t_1)dt_1$

4.不変性
  • システムのパラメータが時間に無関係な時、システムは「不変定数パラメータを持つ」と言う。
  • 入力「$f(t)$」に対するシステムの応答が「$g(t)$」ならば、「$f(t-t_1)$」に対する応答は「$g(t-t_1)$」になる。

    • $\Large \displaystyle L\{f(t-t_1)\} = g(t-t_1)$

  • インパルス「$\delta(t)$」に対するシステムの応答を「$h(t)$」と書けば、

    • $\Large \displaystyle L\{\delta(t)\} = h(t)$

  • となるので「$\delta(t-t_1)$」に対するシステムの応答は以下のように「$h(t-t_1)$」となる。

    • $\Large \displaystyle L\{\delta(t-t_1)\} =h(t:t_1)=h(t-t_1)$

  • 従って、入力「$f(t)$」に対する出力「$g(t)$」は、以下のようになる。

    • $\Large \displaystyle g(t) = \int_{-\infty}^{\infty}f(\tau)h(t;\tau)d\tau =\int_{-\infty}^{\infty}f(\tau)h(t-\tau)d\tau =\int_{-\infty}^{\infty}f(t-\tau)h(\tau)d\tau =f(t)*g(t)$

  • つまり、不変定数パラメータを持つ線形システムの応答は、「$h(t)$」のみの知識から 一意に決定されると言う驚くべき結論に達する。
5.固有関数
  • TODO

1.3 矩形パルス波のフーリエ変換

  • 「$t=-T～+T$」で値が「1」の矩形パルス波「$p_T(t)$」のフーリエ変換は、以下のようになる。

    • $\Large \displaystyle F(\omega)=\int_{-\infty}^{\infty}p_T(t)e^{-j\omega t}dt=\int_{-T}^{T}e^{-j\omega t}dt$

    $\Large \displaystyle =\frac{1}{-j\omega}[e^{-j\omega T}-e^{j\omega T}]=\frac{2\sin(\omega T)}{\omega}$

  • 関数「$f(x)$」が、定数「$t_0$」だけ右方向に時間推移した場合のフーリエ変換は、以下のようになる。

    • $\Large \displaystyle \int_{-\infty}^{\infty}f(t-t_0)e^{j\omega t}dt=\int_{-\infty}^{\infty}f(x)e^{-j\omega(t_0 + x)}dx=F(\omega)e^{-j t_0 \omega} $

  • 従って矩形パルス波「$p_T(t)$」が、定数「$t_0$」だけ右方向に時間推移した場合のフーリエ変換は、以下のようになる。

    • $\Large \displaystyle p_T(t-t_0) \leftrightarrow \frac{2\sin(\omega T)}{\omega} e^{-j\omega t_0}$

1.4 サンプリング定理(サンプリング後)

  • 「$\omega_c$」をサンプリング角周波数とすると、サンプリングパルス列「$f^*(t)$」のフーリエ変換「$F^*(\omega)$」は、以下のようになる。

    • $\Large \displaystyle F^*(\omega) = \sum_{n=-\infty}^{\infty}\frac{\pi}{\omega_c}f_n e^{-jn\pi\frac{\omega}{\omega_c}}$

  • サンプリングパルス列「$f^*(t)$」は、「$F^*(\omega)$」を逆フーリエ変換して以下のようになる。

    • $\Large \displaystyle f^*(t) = \sum_{n=-\infty}^{\infty}\frac{\pi}{\omega_c}f_n \delta(t-\frac{n\pi}{\omega_c})$

1.5 サンプリング定理(サンプリング前)

  • 「$F(\omega)$」が「$f(t)$」のフーリエ積分ならば、以下のような対称性がある。

    • $\Large \displaystyle F(t) \leftrightarrow 2\pi f(-\omega)$

  • この関係は、以下のフーリエ逆変換の式を変形し、「$t$」と「$\omega$」を入れ替える事で得られる。

    • $\Large \displaystyle f(t)=\frac{1}{2\pi}\int_{-\infty}^{\infty} F(\omega)e^{j\omega t}dw$

    $\Large \displaystyle 2\pi f(-t)=\int_{-\infty}^{\infty} F(\omega)e^{-j\omega t}dw$

  • 上式の「$t$」と「$\omega$」を入れ替える事で、以下のように「$F(t)$」のフーリエ積分が得られる。

    • $\Large \displaystyle 2\pi f(-\omega)=\int_{-\infty}^{\infty} F(t)e^{-j\omega t}dt$

1.6 離散フーリエ変換とＺ変換

離散フーリエ変換の定義
  • 離散信号 nx )( の離散フーリエ変換は次式により定義されている。

    • $\Large \displaystyle X(e^{j\omega})=\sum_{n=-\infty}^{\infty}x(n)e^{-j\omega}$

Ｚ変換の定義
  • 離散信号「$x(n)$」のZ変換は次式により定義されている。

    • $\Large \displaystyle X(z)=\sum_{n=-\infty}^{\infty}x(n)z^{-n}$

  • ただし、ｚは複素要素をとる複素変数である。また逆変換は

    • $\Large \displaystyle x(n)=\frac{1}{2\pi i}\oint_C X(z)z^{n-1}dz$

  • と定義される。ただし、積分路Cは収束領域内での原点を内部に含む反時間方向の円周路である。

1.7 離散フーリエ変換とz変換との関係

以上

2.最小二乗法

もくじ

2.1 最小二乗法（一般論）

1.最小二乗法とは
  • 線形回帰モデル(liner regression)は、1つ以上の説明変数「$x$」と被説明変数「$y$」との関係性をモデル化しています。
  • 線形回帰モデルの基底関数(basis function)の違いによって以下のように呼ぶことができます。
    • 単回帰
    • 重回帰
    • 多項式回帰
2.単回帰
  • 単回帰は、説明変数が1変数のみであり、説明変数(1次元)から

    • $\Large \displaystyle X^T=(1,x_1)$

  • と考え、

    • $\Large \displaystyle \hat{y}(x,w)=w_0+w_1 x_1$

  • で、説明変数と被説明変数の関係性をモデル化する。、
  • ここで、線形回帰モデルの基底関数を「$\phi_0(x)=1,\phi_1(x)=x_1$」として、 特徴ベクトル(feature vector)を「$\phi(x)=(\phi_0(x),\phi_1(x))^T$」 とおくと、

    • $\Large \displaystyle \hat{y}(x,w)=w_0+w_1 x_1$

    $\Large \displaystyle =(w_0,w_1) \begin{pmatrix} 1 \\ x_1 \end{pmatrix} $

    $\Large \displaystyle =w^T \phi(x)$

  • 予測値を「$\hat{y}=w_0+w_1 x_n$」を書く事にし、実際の「$y_n$」との誤差は、次のようになる。

    • $\Large \displaystyle y_n - \hat{y_n}=y_n - (w_0+w_1 x_n)$

  • 上記の誤差を最小限にするために、以下の評価関数「$E$」を最小にするパラメータ「$w_0,w_1$」を求める。

    • $\Large \displaystyle E=\sum_{n=1}^{N}(y_n - \hat{y}_n)^2=\sum_{n=1}^{N}(y_n - (w_0+w_1 x_n))^2$

  • 上式が最小になるには、上式の「$w_0$」の偏微分と「$w_1$」の偏微分とが0になれば良い。

    • $\Large \displaystyle E=\sum_{n=1}^{N}(y_n - w_0 - w_1 x_n)^2$

    $\Large \displaystyle =\sum_{n=1}^{N}(y_n^2 + w_0^2 + w_1^2 x_n^2 + 2w_0 y_n -2w_0 w_1 x_n - 2w_1 x_n y_n)$

  • だから、「$w_0$」の偏微分と「$w_1$」の偏微分は、以下のようになる。

    • $\Large \displaystyle \frac{\partial E}{\partial w_0}=\sum_{n=1}^{N}(2w_0 - 2y_n + 2w_1 x_n)=0$

    $\Large \displaystyle \frac{\partial E}{\partial w_1}=\sum_{n=1}^{N}(2w_1 x_n^2 + 2w_0 x_n - 2x_n y_n)=0$

  • 従って、以下の式を得る。

    • $\Large \displaystyle w_0 N + w_1 \sum_{n=1}^{N} x_n=\sum_{n=1}^{N} y_n\ \ (1)$

    $\Large \displaystyle w_0 \sum_{n=1}^{N}x_n + w_1 \sum_{n=1}^{N} x_n^2=\sum_{n=1}^{N} x_n y_n\ \ (2)$

  • (1)式から、以下の「$w_0$」を得る。

    • $\Large \displaystyle w_0= \frac{1}{N}(\sum_{n=1}^{N} y_n - w_1 \sum_{n=1}^{N} x_n)=\bar{y}-w_1\bar{x}$

  • 上式の「$w_0$」を(2)式に代入する。

    • $\Large \displaystyle \sum_{n=1}^{N}x_n(\sum_{n=1}^{N} y_n - w_1 \sum_{n=1}^{N} x_n) + N w_1 \sum_{n=1}^{N} x_n^2=N\sum_{n=1}^{N} x_n y_n$

    $\Large \displaystyle w_1 (-(\sum_{n=1}^{N}x_n)^2 + N \sum_{n=1}^{N} x_n^2)=N\sum_{n=1}^{N} x_n y_n - \sum_{n=1}^{N}x_n\sum_{n=1}^{N} y_n$

  • 上式から、以下の「$w_1$」を得る。

    • $\Large \displaystyle w_1=\frac{N\sum_{n=1}^{N} x_n y_n - \sum_{n=1}^{N}x_n\sum_{n=1}^{N} y_n}{N \sum_{n=1}^{N} x_n^2 - (\sum_{n=1}^{N}x_n)^2}$

    $\Large \displaystyle (w_1=\frac{\bar{xy} - \bar{x}\bar{y}}{\bar{x^2} - \bar{x}^2})$

  • 従って、以下の「$w_0$」を得る。

    • $\Large \displaystyle w_0=\frac{\sum_{n=1}^{N} y_n}{N} - \frac{\sum_{n=1}^{N} x_n}{N} \frac{N\sum_{n=1}^{N} x_n y_n - \sum_{n=1}^{N}x_n\sum_{n=1}^{N} y_n}{N \sum_{n=1}^{N} x_n^2 - (\sum_{n=1}^{N}x_n)^2}$

    $\Large \displaystyle =\frac{\sum_{n=1}^{N} y_n}{N}\frac{N \sum_{n=1}^{N} x_n^2 - (\sum_{n=1}^{N}x_n)^2}{N \sum_{n=1}^{N} x_n^2 - (\sum_{n=1}^{N}x_n)^2} - \frac{\sum_{n=1}^{N} x_n}{N} \frac{N\sum_{n=1}^{N} x_n y_n - \sum_{n=1}^{N}x_n\sum_{n=1}^{N} y_n}{N \sum_{n=1}^{N} x_n^2 - (\sum_{n=1}^{N}x_n)^2}$

    $\Large \displaystyle = \frac{\sum_{n=1}^{N} x_n^2\sum_{n=1}^{N} y_n - \sum_{n=1}^{N}x_n\sum_{n=1}^{N} x_n y_n}{N \sum_{n=1}^{N} x_n^2 - (\sum_{n=1}^{N}x_n)^2}$

    $\Large \displaystyle (w_0=\frac{\bar{x^2}\bar{y} - \bar{x}\bar{xy}}{\bar{x^2} - \bar{x}^2})$

  • 例として「$\{(3,2),(2,4),(-1,1)\}$」を計算して見る。

    • $\Large \displaystyle \sum_{n=1}^{N} x_n=3+2-1=4$

    $\Large \displaystyle \sum_{n=1}^{N} y_n=2+4+1=7$

    $\Large \displaystyle \sum_{n=1}^{N} x_n^2=9+4+1=14$

    $\Large \displaystyle \sum_{n=1}^{N} x_n y_n=6+8-1=13$

  • 従って、パラメータ「$w_0,w_1$」は、以下のようになる。

    • $\Large \displaystyle w_0=\frac{14*7-4*13}{3*14-4*4}=\frac{98-52}{42-16}=\frac{46}{26}=1.77$

    $\Large \displaystyle w_1=\frac{3*13-4*7}{3*14-4*4}=\frac{39-28}{42-16}=\frac{11}{26}=0.42$

3.重回帰
  • 重回帰は、説明変数が多変数となり、説明変数(D次元)から

    • $\Large \displaystyle X^T=(1,x_1,x_2,\cdots,x_D)$

  • と考え、

    • $\Large \displaystyle \hat{y}(x.w)=w_0+w_1 x_1+w_2 x_2+\cdots+w_D x_D$

  • で、説明変数と被説明変数の関係性をモデル化します。、
  • ここで、線形回帰モデルの基底関数を「$\phi_0(x)=1,\phi_1(x)=x_1,\cdots,\phi_D(x)=x_D$」として、 特徴ベクトル(feature vector)を「$\phi(x)=(\phi_0(x),\phi_1(x),\cdots,\phi_D(x))^T$」とおくと、

    • $\Large \displaystyle \hat{y}(x,w)=w_0+w_1 x_1+w_2 x_2+\vdots++w_D x_D$

    $\Large \displaystyle =(w_0,w_1,w_2,\cdots,w_D) \begin{pmatrix} 1 \\ x_1 \\ x_2 \\ \vdots \\ x_D \end{pmatrix} $

    $\Large \displaystyle =w^T \phi(x)$

  • となり、単回帰と同様に内積で表現できます。異なっている部分は基底関数の定義のみです。
  • つまり、重回帰は単回帰を多変数に拡張しています。
4.多項式回帰
  • 多項式回帰は、説明変数に対して非線形な関数「$\phi$」(n次多項式や三角関数等)を基底関数として、

    • $\Large \displaystyle \hat{y}(x,w)=w_0\phi(x_0)+w_1\phi(x_1)+\cdots+w_d\phi(x_D)$

  • で、説明変数と被説明変数の関係性をモデル化します。
  • 例として、多項式

    • $\Large \displaystyle \hat{y}(x,w)=w_0+w_1 x_1+w_2 x_1^2+w_3 \sin(x_1)$

  • に回帰する場合を考えます。
  • この場合、基底関数は「$\phi_0(x)=1,\phi_1(x)=x_1,\phi_2(x)=x_1^2,\phi_3(x)=\sin(x_1)$」、 特徴ベクトルは「$\phi(x)=(\phi_0(x),\phi_1(x),\phi_2(x),\phi_3(x))^T$」となるため、以下のように表現できます。

    • $\Large \displaystyle \hat{y}(x,w)=w_0+w_1 x_1+w_2 x_1^2+w_3 \sin(x_1)$

    $\Large \displaystyle =(w_0,w_1,w_2,w_3) \begin{pmatrix} 1 \\ x_1 \\ x_1^2 \\ \sin(x_1) \end{pmatrix} $

    $\Large \displaystyle =w^T \phi(x)$

  • 多項式回帰では、多項式が「$x$」に対して非線形になっていますが、 「$w」に対しては線形となっているため、線形回帰と呼ばれています。
  • 上記から、基底関数を変更することで、単回帰、重回帰、線形回帰が同一のベクトル表現で表現できることが分かります。
5.ここからが本番の最小二乗法の説明
  • 線形回帰モデルのフィッティングでは、標準的に使われる手法として、最小二乗法があります。
  • 最小二乗法は、実際の測定値「$y$」が「$N$」個ある場合に、 測定値「$y$」と線形回帰モデルでの予測値「$\hat{y}$」の二乗誤差の総和

    • $\Large \displaystyle E=\sum_{n=1}^{N}(y_n-\hat{y}_n)^2$

    $\Large \displaystyle =\sum_{n=1}^{N}(y_n-w^T \phi(x_n))^2$

  • を最小にする「$w$」を求めることで、線形回帰モデルをフィッティングさせます。
  • ここからは、この「$w$」を導出していきます。
6.事前準備
  • 二乗誤差の総和「$E$」を最小にするには、「$E$」をベクトル「$w$」で偏微分して0とすれば、 「$E$」を最小とする「$w$」が導出できます。

    • $\Large \displaystyle \frac{\partial E}{\partial w}=\frac{\partial}{\partial w}\sum_{n=1}^{N}(y_n-\hat{y}_n)^2$

    $\Large \displaystyle =\frac{\partial}{\partial w}\sum_{n=1}^{N}(y_n-w^T \phi(x_n))^2=0$

  • 上記を解くために以下の公式等を準備します。
    • 計画行列
    • ベクトル要素の二乗和
    • 内積の微分
    • 二次形式の微分
6.計画行列
  • 線形回帰モデルでは、単回帰、重回帰、多項式回帰のどれも「$\hat{y}(x,w)=w^T\phi(x)$」で表現できました。
  • そこで、実際の測定値「$y$」が「$N$」個ある場合、「$y$」yの予測値「$\hat{y}$」は、 基底関数を用いて以下の行列で表現できます。

    • $\Large \displaystyle \begin{pmatrix} \hat{y}_1 \\ \hat{y}_2 \\ \vdots \\ \hat{y}_N \end{pmatrix} = \begin{pmatrix} w^T\phi(x_1) \\ w^T\phi(x_2) \\ \vdots \\ w^T\phi(x_N) \\ \end{pmatrix} $

    = $\Large \displaystyle \begin{pmatrix} \phi(x_1)^T \\ \phi(x_2)^T \\ \vdots \\ \phi(x_N)^T \\ \end{pmatrix} w $

    = $\Large \displaystyle \begin{pmatrix} \phi_0(x_1) & \phi_1(x_1) & \cdots & \phi_H(x_1) \\ \phi_0(x_2) & \phi_1(x_2) & \cdots & \phi_H(x_2) \\ \vdots \\ \phi_0(x_N) & \phi_1(x_N) & \cdots & \phi_H(x_N) \\ \end{pmatrix} w $

  • なお、「$\phi_0(x)\equiv 1$」と定義しています。
  • ここで、

    • $\Large \displaystyle \hat{y}= \begin{pmatrix} \hat{y}_1 \\ \hat{y}_2 \\ \vdots \\ \hat{y}_N \end{pmatrix} $

    $\Large \displaystyle \Phi= \begin{pmatrix} \phi_0(x_1) & \phi_1(x_1) & \cdots & \phi_H(x_1) \\ \phi_0(x_2) & \phi_1(x_2) & \cdots & \phi_H(x_2) \\ \vdots \\ \phi_0(x_N) & \phi_1(x_N) & \cdots & \phi_H(x_N) \\ \end{pmatrix} $

  • とおくと、「$\hat{y}=\Phi w$」と表現でき、行列「$\Phi$」を計画行列と呼びます。
  • この計画行列「$\Phi$」は、線形回帰モデルの基底関数を行列で表現しています。
7.ベクトル要素の二乗和
  • ベクトル「$x=(x_1,x_2,\cdots,x_n)^T$」の要素の二乗和は

    • $\Large \displaystyle x_1^2+x_2^2+\cdots x_n^2= (x_1,x_2,\cdots,x_n) \begin{pmatrix} x_1 \\ x_2 \\ \vdots \\ x_n \\ \end{pmatrix} =x^T x $

8.内積の微分
  • 「$w=(w_1,w_2,\cdots, w_n)^T,x=(x_1,x_2,\cdots, x_n)^T$」とすると、

    • $\Large \displaystyle w^T x=w_1 x_1+w_2 x_2+\cdots+w_n x_n$

  • スカラーである内積「$w^T x$」をベクトル「$w$」で微分することを考えます。
  • ベクトルの各要素で微分すれば良いので、

    • $\Large \displaystyle \frac{\partial}{\partial w}w^T x= \begin{pmatrix} \frac{\partial}{\partial w_1}w^T x \\ \frac{\partial}{\partial w_2}w^T x \\ \vdots \\ \frac{\partial}{\partial w_n}w^T x \\ \end{pmatrix} $

    $\Large \displaystyle = \begin{pmatrix} x_1 \\ x_2 \\ \vdots \\ x_n \\ \end{pmatrix} $

    $\Large \displaystyle =x$

  • となります。
9.二次形式の微分
10.二乗誤差の総和を最小にする値の導出

2.2 最小二乗法（カルマンフィルタ用）

1.単回帰
  • カルマンフィルタ用最小二乗法では、観測値「$y$」から状態「$x$」を推定するために使用される。
  • 従って、説明変数が観測値「$y$」となり、被説明変数が状態「$x$」となる。
  • すると説明変数(一次元)は、

    • $\Large \displaystyle X^T=(1,y_1)$

  • と考え、

    • $\Large \displaystyle \hat{x}(y,w)=w_0+w_1 y_1$

  • で、説明変数と被説明変数の関係性をモデル化する。
  • ここで、線形回帰モデルの基底関数を「$\phi_0(x)=1,\phi_1(x)=x_1$」として、 特徴ベクトル(feature vector)を「$\phi(x)=(\phi_0(x),\phi_1(x))^T$」 とおくと、

    • $\Large \displaystyle \hat{x}(y,w)=w_0+w_1 y_1$

    $\Large \displaystyle =(w_0,w_1) \begin{pmatrix} 1 \\ y_1 \end{pmatrix} $

    $\Large \displaystyle =w^T \phi(y)$

  • 予測値を「$\hat{x}=w_0+w_1 y_n$」を書く事にし、実際の「$x_n$」との誤差「$e$」は、次のようになる。

    • $\Large \displaystyle e=x_n - \hat{x_n}=x_n - (w_0+w_1 y_n)$

  • 上記の誤差の平均を0にするために、以下の評価関数「$E_1$」を最小にするパラメータ「$w_0,w_1$」を求める。

    • $\Large \displaystyle E_1=E[e]=\sum_{n=1}^{N}(x_n - \hat{x}_n)=\sum_{n=1}^{N}(x_n - (w_0+w_1 y_n))$

  • ここで、観測値「$y$」は、以下の観測方程式で観測されるものとする。ただし「$w$」は、観測ノイズであり「$N(0,\sigma_x^2)$」に従うものとする。

    • $\Large \displaystyle y_n=cx_n+w$

  • すると、評価関数「$E_1$」は、以下のようになる。

    • $\Large \displaystyle E_1=\sum_{n=1}^{N}(x_n - w_0 - w_1 y_n))$

    $\Large \displaystyle =\sum_{n=1}^{N}(x_n - w_0 - w_1(cx_n + w)))$

    $\Large \displaystyle =\sum_{n=1}^{N}((1- w_1 c)x_n - w_0 - w_1 w)))$

    $\Large \displaystyle =(1- w_1 c)\bar{x} - w_0 - w_1\bar{w}=0$

  • 従って、「$w_0$」は、以下のようになる。

    • $\Large \displaystyle w_0=(1- w_1 c)\bar{x} - w_1\bar{w}=\bar{x}-w_1(c\bar{x}+\bar{w})$

  • 上記の誤差の分散を最小限にするために、以下の評価関数「$E_2$」を最小にするパラメータ「$w_0,w_1$」を求める。

    • $\Large \displaystyle E_2=E[\{e - E(e)\}^2]$

    $\Large \displaystyle =E[[(x_n - w_0 - w_1 y_n) - \{ (1 - w_1 c)\bar{x} - w_0 - w_1 \bar{w}\}]^2]$

    $\Large \displaystyle =E[\{ x_n - w_1(cx_n +w) -(1 - w_1 c)\bar{x} + w_1\bar{w}\}^2]$

    $\Large \displaystyle =E[\{ (1 - w_1 c)(x_n - \bar{x}) - w_1(w - \bar{w}) \}^2]$

    $\Large \displaystyle =(1 - w_1 c)^2 E[(x-\bar{x})^2] - 2(1 - w_1 c)w_1 E[(x-\bar{x})(w-\bar{w})] +w_1^2 E[(w-\bar{w})^2]$

    $\Large \displaystyle =(1 - w_1 c)^2 \sigma_x^2 + w_1^2 \sigma_w^2$

    $\Large \displaystyle =\sigma_x^2 - 2c\sigma_x^2 w_1 + c^2 \sigma_x^2 w_1^2 + w_1^2 \sigma_w^2$

    $\Large \displaystyle =\sigma_x^2 - 2c\sigma_x^2 w_1 + (c^2 \sigma_x^2 + \sigma_w^2) w_1^2$

  • 上式が最小になるには、上式の評価関数「$E_2$」の「$w_1$」についての偏微分が0になれば良い。

    • $\Large \displaystyle \frac{\partial E_2}{\partial w_1}=\frac{\partial E_2}{\partial w_1}((1 - 2w_1 c + w_1^2 c^2) \sigma_x^2 + w_1^2 \sigma_w^2)$

    $\Large \displaystyle =\frac{\partial E_2}{\partial w_1}(\sigma_x^2 - 2w_1 c \sigma_x^2 + w_1^2 c^2 \sigma_x^2 + w_1^2 \sigma_w^2)$

    $\Large \displaystyle =- 2c \sigma_x^2 + 2(c^2 \sigma_x^2 + \sigma_w^2)w_1=0$

  • 従って、以下の式を得る。

    • $\Large \displaystyle w_1=\frac{c \sigma_x^2}{c^2 \sigma_x^2 + \sigma_w^2}$

    $\Large \displaystyle =\frac{c \sigma_w^{-2}}{\sigma_x^{-2} + c^2 \sigma_w^{-2}}$

  • ここで、以下のように「$\sigma^2$」を定義する。

    • $\Large \displaystyle \sigma^2=\frac{1}{\sigma_x^{-2} + c^2 \sigma_w^{-2}}$

  • 「$\sigma^2$」を使うと「$w_1$」は、以下のように表す事ができる2$」を定義する。

    • $\Large \displaystyle w_1=c\frac{\sigma^2}{\sigma_w^2}$

  • 「$\sigma^2$」を使うと評価関数「$E_2$」は、以下のようになる。

    • $\Large \displaystyle E_2=(1 - w_1 c)^2 \sigma_x^2 + w_1^2 \sigma_w^2$

    $\Large \displaystyle =\sigma_x^2 - 2c\sigma_x^2 w_1 + (c^2 \sigma_x^2 + \sigma_w^2) w_1^2$

    $\Large \displaystyle =\sigma_x^2\frac{c^2 \sigma_x^2 + \sigma_w^2}{c^2 \sigma_x^2 + \sigma_w^2} - 2c\sigma_x^2 \frac{c\sigma_x^2}{c^2 \sigma_x^2 + \sigma_w^2} + \frac{c^2 \sigma_x^4}{c^2 \sigma_x^2 + \sigma_w^2}$

    $\Large \displaystyle =\frac{\sigma_x^2 \sigma_w^2}{c^2 \sigma_x^2 + \sigma_w^2}=\sigma^2$

2.3 自己回帰モデル（AR モデル）

1.自己回帰モデル（AR モデル）とは
  • 自己回帰モデル（AR モデル）は、ある時刻「$t$」の値を、時刻「$t$」よりも古いデータを使って回帰するモデルです。
  • 自己回帰モデルは、自己相関の高い時系列データのモデリングに有効です。
  • 回帰式は次のように表せます。「$\phi_i$」は自己回帰係数、「$p$」は次数です。
  • また、「$\epsilon_t$」は誤差項であり、平均「$0$」かつ分散「$\sigma^2$」のホワイトノイズです。

    • $\Large \displaystyle y_t=\phi_0 + \phi_1 y_{t-1} + \phi_2 y_{t-2} + \cdots + \phi_p y_{t-p} + \epsilon_t$

2.自己回帰 AP(1) モデル
  • 自己回帰 AP(1) モデルの回帰式は、以下のように表せます。

    • $\Large \displaystyle y_t=\phi_0 + \phi_1 y_{t-1} + \epsilon_t$

  • 自己回帰 AP(1) モデルの平均は、以下のようになります。

    • $\Large \displaystyle E[y_t]=\phi_0 + \phi_1E[y_{t-1}] + E[\epsilon_t]$

  • 自己回帰モデルを適用できるデータは定常性を持ちます。
  • 定常性が仮定されるとき、「$E[y_t]=\mu$」かつ「$E[\epsilon_t] = 0$」が成り立つので、AP(1) の平均値は次のように計算できます。

    • $\Large \displaystyle \mu=\phi_0 + \phi_1\mu$

  • 従って、AP(1) の平均値は次のようになります。

    • $\Large \displaystyle \mu=\frac{\phi_0}{1-\phi_2}$

  • 次に、AP(1) のモデル式の両辺に「$y_t$」をかけて、AP(1) の自己分散を計算します。

    • $\Large \displaystyle E[{y_t}^2]=\phi_0E[y_t] + \phi_1E[y_t y_{t-1}] + E[y_t\epsilon_t]$

  • 自己回帰モデルのとき、データの定常性が仮定されることから、上式は次のように置き換えられる。

2.4 自己回帰モデルとパワースペクトル

  • ー次 元 時 系 列 用x(t) ,サ ン プ リ ン グ 周 期 を ⊿t,時 刻t=n⊿t(nは 整 数)に お け る サ ン プ ル 値 をx(n)で 表 わ す こ とscす る.
  • 自己回帰モデルとは、時系列x(n)の生成過程用P個の係数 αp1,ap2,...ap,を 用 い て は 式 で 表 わすものである,

    • $\Large \displaystyle x(n)=-a_{p1}x(n-1)-a_{p2}x(n-2)-\cdots-a_{pp}x(n-p)+ \epsilon(n)$

  • 本来であれば，推定すべき分布の解析解が得られれば最も良いはずであるが，粒子フィルタが想定する非線形・非ガウス状態空間モデルでは、 そのような解析解が得られることはほとんどない．
  • そのため，モンテカルロ法で状態ベクトル「$x_t^i$」を発生させ，さらにその尤度「$w_t^i$」を適切に評価して経験分布を構成する．
  • たとえば，フィルタ分布「$p(x_t|y_{1:t})$」を経験分布で近似することを考えよう．
  • 粒子フィルタでは以下のような記法が採用されることが多い．

    • $\Large \displaystyle p(x_t|y_{1:t})\simeq \sum_{i=1}^{M}\frac{w_t^i}{\sum_{i=1}^{M}w_t^i}\delta(x_t-x_{t|t-1}^i)$

2.5 線形・ガウス型の状態空間モデル

以上

3.ベイズ統計

もくじ

3.1 一次元正規分布

1.確率密度関数
  • 以下の確率密度関数「$p(x)$」で与えられる確率変数を正規分布と言う。

    • $\Large \displaystyle p(x)=\frac{1}{\sqrt{2\pi\sigma^2}}exp\{-\frac{1}{2}\frac{(x-\mu)^2}{\sigma^2}\}$

  • 「$\mu$」は平均値であり、「$\sigma$」は標準偏差であり、「$\sigma^2$」は分散である。
  • 上記のような正規分布を「$N(\mu,\sigma^2)$」と表記する。「$N(0,1)$」は「標準正規分布」と呼ばれる。
2.正規分布の性質
  • ①1次モーメント（平均）と2次モーメント（分散）を与える事により、それに従う確率法則を完全に確定できる。
  • ②確率変数「$X$」が正規分布「$N(\mu,\sigma^2)$」に従う時、その線形変換「$Y=aX+b$」は「$N(a\mu+b,a^2\sigma^2)$」に従う。
  • ③
  • ④正規分布の確率密度関数をフーリエ変換すると、やはり正規分布の関数になる。

3.2 多次元正規分布

1.多次元正規分布の性質
  • ①期待値ベクトルと共分散行列を指定すれば分布が定まる。
  • ②相関がないだけで独立だと断言できる。
  • ③多次元正規分布を線形変換したら、また多次元正規分布になる。

3.3 共分散と相関係数

1.期待値と分散
  • 確率変数「$X$」に対する期待値「$\mu$」は、以下のように定義される。

    • $\Large \displaystyle \mu=E[X]$

  • 確率変数「$X$」に対する分散「$V[X]$」は、以下のように定義される。

    • $\Large \displaystyle V[X]=E[(X-\mu)^2]$

  • 確率変数「$X$」に対する標準偏差「$\sigma$」は、以下のように定義される。

    • $\Large \displaystyle \sigma=\sqrt{V[X]}$

2.共分散
  • 確率変数「$X,Y$」に対する共分散「$Cov[X,Y]$」を導入する。
  • 確率変数「$X,Y$」の期待値をそれぞれ「$\mu,\nu$」とする時、「$X,Y$」の共分散を以下のように定義する。

    • $\Large \displaystyle Cov[X,Y]=E[(X-\mu)(Y-\nu)]$

3.共分散の性質
  • 共分散の定義から、以下のような共分散の性質があります。

    • $\Large \displaystyle Cov[X,Y]=Cov[Y,X]$

    $\Large \displaystyle Cov[X,X]=V[X]$

    $\Large \displaystyle Cov[X+a,Y+b]=Cov[X,Y]$

    $\Large \displaystyle Cov[aX,bY]=abCov[X,Y]$

4.相関係数
  • 上記の共分散の性質から、確率変数「$X,Y$」に定数「$a,b$」を掛けたら共分散が「$ab$」倍になると言う問題があった。
  • そこで、「尺度」に惑わされない指標として以下の相関係数「$\rho_{XY}$」を導入する。

    • $\Large \displaystyle \rho_{XY}=\frac{Cov[X,Y]}{\sqrt{\mathstrut V[X]}\sqrt{\mathstrut V[Y]}}$

  • つまり、以下のように確率変数「$X,Y$」を分散の平方根（標準偏差）で割って、正規化する。

    • $\Large \displaystyle \hat{X}=\frac{X}{\sqrt{V[X]}}$

    $\Large \displaystyle \hat{Y}=\frac{Y}{\sqrt{V[Y]}}$

  • つまり、以下のように確率変数「$X,Y$」を分散の平方根（標準偏差）で割って、正規化したものの共分散が相関係数になる。

    • $\Large \displaystyle Cov[\hat{X},\hat{Y}]=Cov[\frac{X}{\sqrt{V[X]}},\frac{Y}{\sqrt{V[Y]}}]=\frac{Cov[X,Y]}{\sqrt{\mathstrut V[X]}\sqrt{\mathstrut V[Y]}}=\rho_{XY}$

5.共分散や相関係数では、測れない事

3.4 共分散行列

1.共分散行列＝分散と共分散の一覧表
  • 例えば、以下のように確率変数「$X_1,X_2,X_3$」の全ての組み合わせの共分散の一覧表を作り、行列化する。
    • $X_1$$X_2$$X_3$
    $X_1$$Cov(X_1,X_1)$$Cov(X_1,X_2)$$Cov(X_1,X_3)$
    $X_2$$Cov(X_2,X_1)$$Cov(X_2,X_2)$$Cov(X_2,X_3)$
    $X_3$$Cov(X_3,X_1)$$Cov(X_3,X_2)$$Cov(X_3,X_3)$
  • ここで、確率変数「$X_1,X_2,\cdots,X_n$」を縦に並べた縦ベクトル「$X$」を導入すると、その分散「$V[X]$」は、以下のように行列化された共分散一覧になる。

    • $\Large \displaystyle V[X]= \begin{pmatrix} V[X_1] & Cov[X_1,X_2] & Cov[X_1,X_3] \\ Cov[X_2,X_1] & V[X_2] & Cov[X_2,X_3] \\ Cov[X_3,X_1] & Cov[X_3,X_2] & V[X_3] \\ \end{pmatrix} $

2.ベクトルでまとめて書く
  • 共分散行列を宣言しただけでは、あまりご利益はないので、確率変数「$X_1,X_1,\cdots,Xn$」を縦に並べた、 以下のような縦ベクトル「$X$」を導入する。

    • $\Large \displaystyle X= \begin{pmatrix} X_1 \\ X_2 \\ \vdots \\ X_n \end{pmatrix} $

  • 縦ベクトル「$X$」を使うと、共分散行列を以下のようにまとめて書ける。

    • $\Large \displaystyle V[X]=E[(X-\mu)(X-\mu)^T]$

  • ただし、「$\mu$」は、以下のように定義される。

    • $\Large \displaystyle \mu=E[X]$

  • つまり、縦ベクトル「$X$」の共分散行列（「$X$」の分散「$V[X]$」）は以下のようになる。

    • $\Large \displaystyle V[X]= \begin{pmatrix} V[X_1] & Cov[X_1,X_2] & Cov[X_1,X_3] \\ Cov[X_2,X_1] & V[X_2] & Cov[X_2,X_3] \\ Cov[X_3,X_1] & Cov[X_3,X_2] & V[X_3] \\ \end{pmatrix} $

3.5 自己共分散の定義

1.自己共分散の定義

    • $\Large \displaystyle \gamma_k =\frac{1}{T}\sum_{t=k+1}^{T}(y_t-\bar{y})(y_{t-k}-\bar{y})$

4.多次元正規分布

もくじ

4.1 ランダムウォーク

1.の近似

5.確率過程から

もくじ

5.1 ランダムウォーク

非ガウス分布の近似

5.2 カルマンフィルタ

1.設定
  • 時刻tにおける的の位置を「$X_t$」で表す。「$X_t$」は、直接、測定できないものとする。
  • 時刻tにおける的の位置「$X_t$」は、時刻t-1の的の位置「$X_{t-1}$」より「$W_t$」だけ移動するものとする。
  • 時刻tにおける的の位置の測定値を「$Y_t$」で表す。「$Y_t$」は、誤差「$Z_t$」を含むものとする。

    • $\Large \displaystyle X_t=X_{t-1}+W_t$

    $\Large \displaystyle Y_t=X_{t}+Z_t$

  • ただし、以下が成立し、「$X_0,W_1,W_2,\cdots,Z_1,Z_2,\cdots$」は、独立とし、「$\sigma_0^2,\alpha^2,\beta^2$」は既知の定数とする。

    • $\Large \displaystyle X_0\sim N(0,\sigma_0^2)$

    $\Large \displaystyle W_t\sim N(0,\alpha^2)$

    $\Large \displaystyle Z_t\sim N(0,\beta^2)$

  • 問題設定は「$Y_1,\cdots,Y_t$」の値を見て「$X_t$」の推定値「$\mu_1$」を当てる事である。
2.導出
  • t=1の場合は、以下のようになる。

    • $\Large \displaystyle X_1=X_0+W_1$

    $\Large \displaystyle Y_1=X_1+Z_1=X_0+W_1+Z_1$

  • 行列を使うと、以下のように表現できる。

    • $\Large \displaystyle \begin{pmatrix} Y_1 \\ X_1 \end{pmatrix} = \begin{pmatrix} 1 & 1& 1 \\ 0 & 1 & 1 \end{pmatrix} \begin{pmatrix} Z_1 \\ W_1 \\ X_0 \end{pmatrix} $

  • 右辺のベクトル「$(Z_1,W_1,X_0)^T$」は、前提から3次元正規分布「$N(0,diag(\beta^2,\alpha^2,\sigma_0^2))$」に従う。
  • 従って、その線形変換である左辺のベクトル「$(Y_1,X_1)^T$」も2次元正規分布「$N(0,V_1)$」に従う。

5.3 マルコフ過程

ServiceData

以上

6.線形カルマンフィルタ

もくじ

6.1 カルマンフィルタリング問題

1.線形離散時間状態空間モデル
  • 以下のようなスカラ時系列データ「$y(k)$」が、線形離散時間状態空間モデルで記述されると仮定する。

    • $\Large \displaystyle x(k+1)=Ax(k)+bv(k)$

    $\Large \displaystyle y(k)=c^Tx(k)+w(k)$

  • ここで、「$x(k)$」は、n次元状態ベクトルである。
  • 「$A$」は、n☓n行列である。
  • 「$v(k)$」は、平均値「$0$」、分散「$\sigma_v^2$」の正規性白色雑音である。システム雑音とも呼ばれる。
  • 「$w(k)$」は、平均値「$0$」、分散「$\sigma_w^2$」の正規性白色雑音である。観測雑音とも呼ばれる。
  • 「$v(k)$」と「$w(k)$」は、お互いに独立であると仮定する。
  • 状態空間モデルの係数行列・ベクトル「$(A,b,c)$」は、既知であると仮定する。
  • この仮定は、時系列の状態空間モデルは、何等かの時系列モデリング法によって、既に得られているとする。
  • また、係数行列・ベクトル「$(A,b,c)$」が時間によらず一定である時、上式で表される線形離散時間状態空間モデルは、 線形時不変（LTI: Linear Time-Invarient）システムとなり、その結果、得られる時系列「$y(k)$」は、定常確率過程になる。
  • 以上の準備の元で、時刻kまでの時系列データ「$y(i),i=1,2,\cdots,k$」を用いて、 n次元状態変数ベクトル「$x(k)$」を推定する問題（フィルタリング問題）を考える。
  • 特に、「$x(k)$」を何らかの意味で最もよく表す最適推定値「$\hat{x}(k)$」を最適フィルタを設計するものとする。
2.最適制御と評価関数
  • 最適制御と言う場合は、必ず評価関数（損失関数）を定義する必要がある。
  • 評価関数（損失関数）の値が最小（又は最大）になる場合を最適制御と呼ぶ。
3.カルマンフィルタリング問題
  • 時系列データ「$y(i),i=1,2,\cdots,k$」に基づいて、状態「$x(k)$」のMSW（平均二乗誤差）の 最小値を与える推定値を見つける事を「カルマンフィルタリング問題」と言う。
  • この推定値の事を「最小平均二乗誤差推定値（MMSE:Minimum Mean Square Error）」と呼ぶ。

6.2 逐次処理

1.カルマンフィルタリング問題とは

6.3 時系列に対するカルマンフィルタ

1.事前推定値と事後推定値
  • 時刻「$k$」における状態「$x(k)$」の推定値として、以下の二つの量を定義する。
  • 事前推定値「$\hat{x}^{-}(k)$」
    • 時刻「$k-1$」までに利用可能なデータに基づいた、時刻「$k$」における予測推定値を意味する。
  • 事後推定値「$\hat{x}(k)$」
    • 時刻「$k$」までに利用可能なデータに基づいた、即ち「$y(k)$」も用いた「$x」のフィルタリング時刻推定値。
    • これが時刻「$k$」において求めるべき推定値である。
  • カルマンフィルタでは、上図のように状態推定値は、時間更新されていく。
  • 上図では、状態推定値（1次モーメント）「$\hat{x}(k-1)$」と共分散行列「$P(k-1)$」について示している。
  • 図において、1時刻前の推定値「$\hat{x}(k-1)$」から現時点の事前推定値「$\hat{x}^{-}(k)$」までの処理を「予測ステップ」と呼ぶ。
  • また、現時点の事前推定値「$\hat{x}^{-}(k)$」から現時点の事後推定値「$\hat{x}(k)$」までの処理を「フィルタリングステップ」と呼ぶ。
2.線形予測器の構成
  • 「最小二乗法」と同じように、ベイズ予測による以下のような事前推定値「$\hat{x}^{-}(k)$」と観測値「$y(k)$」に関する線形予測器を仮定する。

    • $\Large \displaystyle \hat{x}(k)=G(k)\hat{x}^{-}(k)+g(k)y(k)$

  • 上式は、まさにベイズ推定であり、以下のような意味を持つ。

    • 事後推定値 = G(k)・事前推定値 + g(k)・観測値

  • 問題は「$G(k)$」と「$g(k)$」をどうやって決定するかである。
  • ここでは、以下の事後状態推定誤差「$\tilde{x}(k)=x(k)-\hat{x}(k)$」に関する観測値「$y(k)$」との直交性を利用する。

    • $\Large \displaystyle \tilde{x}(k)=x(k)-\hat{x}(k)$

    $\Large \displaystyle E[\tilde{x}(k)y(i)]=0,\ \ i=1,2,\cdots,k-1$

3.カルマンゲインの決定法
2.線形予測器の構成
3.共分散行列の更新

6.4 数値シュミレーション例

6.5 システム制御のためのカルマンフィルタ

以上

7.システム同定

もくじ

7.1 LTIシステムの表現

1.離散時間LTIシステム
  • システムへの入力信号を「$u(k)$」、出力信号を「$y(k)$」とすると、

    • $\Large \displaystyle y(k)=\sum_{l=0}^{\infty}g(l)u(k-l)$

  • なる「畳み込み積分」が成り立つ。ただし、「$g(k)$」は、システムのインパルス応答である。
2.連続時間LTIシステム
  • システムへの入力信号を「$u(t)$」、出力信号を「$y(t)$」とすると、

    • $\Large \displaystyle y_t=\int_0^t g(\tau)u(t-\tau)d\tau=g(t)*u(t)$

  • なる「畳み込み積分」が成り立つ。ただし、「$g(t)$」は、システムのインパルス応答である。

7.2 入出力関係

1.離散時間LTIシステム
  • 上記の離散時間LTIシステムをz変換すると、以下の式が得られる。

    • $\Large \displaystyle y(z)=G(z)u(z)$

  • 但し、「$u(z)$」「$y(z)$」は、それぞれ「$u(k)$」「$y(k)$」のz変換であり、以下のように定義される。

    • $\Large \displaystyle u(z)=\sum_{l=0}^{\infty}u(k)z^{-k}$

    $\Large \displaystyle y(z)=\sum_{l=0}^{\infty}y(k)z^{-k}$

  • 但し、「$G(z)$」は、「$g(k)$」のz変換であり、以下のように定義され、離散時間LTIシステムの伝達関数と呼ばれる。

    • $\Large \displaystyle G(z)=\sum_{l=0}^{\infty}g(k)z^{-k}$

  • 一般に、離散時間LTIシステムは、以下のような線形定係数差分方程式で記述できる。

    • $\Large \displaystyle y(k)+a_1 y(k-1)+\cdots+a_n y(k-n)=b_1 u(k-1)+b_2 u(k-2)+\cdots+b_mu(k-m)$

  • 上式を初期値をゼロとして、z変換すると、以下のようになる。

    • $\Large \displaystyle (1+a_1 z^{-1}+\cdots+a_n z^{-n})y(z)=(b_1 z^{-1}+b_2 z^{-2}+\cdots+b_m z^{-m})u(z)$

  • 従って、以下のような伝達関数「$G(z)$」が得られる。

    • $\Large \displaystyle G(z)=\frac{y(z)}{u(z)}=\frac{b_1 z^{-1}+b_2 z^{-2}+\cdots+b_m z^{-m}}{1+a_1 z^{-1}+\cdots+a_n z^{-n}}$

2.連続時間LTIシステム

7.3 ARXモデル

1.シフトオペレータ「$q$」
  • シフトオペレータ「$q$」は、以下のような性質を持つ。

    • $\Large \displaystyle qy(k)=y(k+1)$

    $\Large \displaystyle q^{-1}y(k)=y(k-1)$

2.ARXモデルの式
  • ARXモデルは、ARモデルに外生（因）入力(eXogenous)が加わったモデルである。
  • 従って、ARXモデルの式は、以下のように定義される。

    • $\Large \displaystyle A(q)y(k)=B(q)u(k)+w(k)$

  • ここで、「$y$」は出力、「$u$」は入力、「$w$」は誤差やノイズを意味し、「$k$」は離散時間系におけるサンプル点を意味する。
  • 「$q$」はシフトオペレータと呼ばれるもので、「$z$」変換における演算子「$z$」と意味は同じである。
  • ARXモデルのブロック図は下記のようになる。
  • すなわちARXモデルは、コントローラによって決定される操作量「$u$」とプラントからの出力「$\frac{w(k)}{A(z)}$」を利用したシステム同定であり、ノイズは正確な同定を阻害する因子と言える。
  • 非常に分かり易いですね。

以上

8.ARモデルの実例

もくじ

8.1 カルマンフィルタ

1.カルマンフィルタの問題設定
  • 測定には誤差がつきものである。誤差を減らす策としては、同一条件で独立な測定を何度も行う事である。
  • 「独立同一分布（i.i.d.）」と言う前提の元でn回測定して平均すれば、大数の法則により標準偏差が「$\frac{1}{\sqrt{2}}$」になる。
  • しかし「止まった的を射る」には、この方法でいいのだが、「動いている的を射る」には、この方法では上手く行かない。
  • 時刻kにおける的の位置を、実数値の確率変数「$x(k)$」で表す事にする。
  • 我々は、的の位置である「$x(k)$」を当てたいが、直接知る事は出来ないとする。
  • 従って、代わりに観測誤差「$w(k)$」を含む観測値「$y(k)$」を得る事が出来るとする。
  • 次に、現在の的の位置である「$x(k)$」」から、次の時間の的の位置である「$x(k+1)$」を予想する。
  • 「k」と「k+1」が短い時間差であれば、「$x(k+1)$」は「$x(k)$」」から少しだけの差「$v(k)$」を足した位置と予想できる。
  • 従って、以下のような式の確率モデルを作る事が出来る。

    • $\Large \displaystyle x(k+1)=x(k)+v(k)$

    $\Large \displaystyle y(k)=x(k)+w(k)$

  • ただし、以下のような確率分布（偏差は定数）を仮定し、お互いに独立であるとする。

    • $\Large \displaystyle x(0)\sim N(0,\sigma_0^2)$

    $\Large \displaystyle v(k)\sim N(0,\sigma_v^2)$

    $\Large \displaystyle w(k)\sim N(0,\sigma_w^2)$

  • つまり、カルマンフィルタの問題設定は、観測値「$y(1),y(2),\cdots, y(k)$」からし、状態変数(位置$x(k)$)を推定する事になる。
2.推定誤差
  • 信号「$x$」は、以下の観測式で観測されるものとする。
  • ただし、「$y$」は、観測値であり、「$w$」は、正規分布「$N(0,\sigma^2)$」に従う観測雑音とする。

    • $\Large \displaystyle y=cx +w$

  • 推定問題は、雑音に汚された観測値「$y$」から、信号「$x$」の推定値「$\hat{x}$」を求める事である。
  • さて、線形推定則を仮定すると、信号「$x$」の推定値「$\hat{x}$」は、以下のような観測値「$y$」の1次式で表される。

    • $\Large \displaystyle \hat{x}=f(y)=\alpha y + \beta$

  • 以下のような推定誤差(error estimate)を定義し、この平均二乗誤差（MSE: Mean Square Error）が 最小になるように、線形推定則のパラメータ「$\alpha$」と「$\beta$」を決定する問題を考える。

    • $\Large \displaystyle e=x-\hat{x}$

  • 以下では、「①推定誤差の平均がゼロにする」と「②推定誤差の分散を最小にする」と言う 二つの条件を使って、未知数「$\alpha,\beta$」を決定する。
3.推定誤差の平均値をゼロにする
  • この条件を計算すると、

    • $\Large \displaystyle E[e]=E[x-\hat{x}]=E[x-\alpha y - \beta]$

    $\Large \displaystyle =E[x-\alpha (cx+w) - \beta]$

    $\Large \displaystyle =(1-\alpha c)\bar{x}-\alpha\bar{w}-\beta=0$

  • となるから、

    • $\Large \displaystyle \hat{\beta}=(1-\alpha c)\bar{x}-\alpha\bar{w}=\bar{x}-\alpha(c\bar{x}+\bar{w})$

  • の時、推定値はバイアスを持たない。この時「$\hat{x}$」は不偏推定値と呼ばれる。
4.推定誤差の分散を最小にする
  • 推定値の平均二乗誤差（MSE）は、以下のようになる。

    • $\Large \displaystyle E[\{e-E(e)\}^2]=E[[(x-\alpha y-\beta)-\{(1-\alpha c)\bar{x}-\alpha \bar{w}-\beta \}]^2]$

    $\Large \displaystyle =E[\{ x-\alpha(cx+w)-(1-\alpha c)\bar{x}+\alpha \bar{w}\}^2]$

    $\Large \displaystyle =E[\{ (1-\alpha c)(x-\bar{x})-\alpha(w-\bar{w}) \}^2]$

    $\Large \displaystyle =E[(1-\alpha c)^2(x-\bar{x})^2+\alpha^2(w-\bar{w})^2 -2\alpha(1-\alpha c)(x-\bar{x})(w-\bar{w}) ]$

    $\Large \displaystyle =(1-\alpha c)^2E[(x-\bar{x})^2]+\alpha^2E[(w-\bar{w})^2] -2\alpha(1-\alpha c)E[(x-\bar{x})(w-\bar{w})] $

    $\Large \displaystyle =(1-\alpha c)^2\sigma_x^2+\alpha^2\sigma_w^2$

  • 上式を「$\alpha$」に関する2次式に書き直して平方完成すると以下のようになる。

    • $\Large \displaystyle E[\{e-E(e)\}^2]=(1-2\alpha c+\alpha^2 c^2)\sigma_x^2+\alpha^2\sigma_w^2$

    $\Large \displaystyle =(c^2 \sigma_x^2 +\sigma_w^2)\alpha^2 - 2c\sigma_x^2 \alpha + \sigma_x^2$

    $\Large \displaystyle =(c^2 \sigma_x^2 +\sigma_w^2)[\alpha^2 -\frac{2c\sigma_x^2}{c^2 \sigma_x^2 +\sigma_w^2}\alpha+\frac{\sigma_w^2}{c^2 \sigma_x^2 +\sigma_w^2} ]$

    $\Large \displaystyle =(c^2 \sigma_x^2 +\sigma_w^2)(\alpha -\frac{c\sigma_x^2}{c^2 \sigma_x^2 + \sigma_w^2})^2 - \frac{c^2 \sigma_x^4}{c^2 \sigma_x^2 + \sigma_w^2} + \sigma_x^2$

    $\Large \displaystyle =(c^2 \sigma_x^2 +\sigma_w^2)(\alpha -\frac{c\sigma_x^2}{c^2 \sigma_x^2 + \sigma_w^2})^2 - \frac{\sigma_x^2 \sigma_w^2}{c^2 \sigma_x^2 + \sigma_w^2}$

    $\Large \displaystyle =(c^2 \sigma_x^2 +\sigma_w^2)(\alpha -\frac{c\sigma_w^{-2}}{c^2 \sigma_w^{-2} + \sigma_x^{-2}})^2 - \frac{1}{c^2 \sigma_w^{-2} + \sigma_x^{-2}}$

  • いま、

    • $\Large \displaystyle \sigma^2=\frac{1}{c^2 \sigma_w^{-2} + \sigma_x^{-2}}$

  • とおくと、上式は、以下のように変形できる。

    • $\Large \displaystyle E[\{e-E(e)\}^2]=(c^2 \sigma_x^2+\sigma_w^2)(\alpha - c\sigma_w^{-2}\sigma^2)^2+\sigma^2$

  • これより、

    • $\Large \displaystyle \alpha=c\frac{\sigma^2}{\sigma_w^2}$

  • の時、推定誤差分散は、最小値、

    • $\Large \displaystyle E[\{e-E(e)\}^2]=\sigma^2$

  • をとる事は分かる。
  • この時、「$\hat{x}$」は、「最小分散推定値」と呼ばれる。
5.カルマンフィルタの導出
  • 以下のようなスカラ時系列データ「$y(k)$」が、線形離散時間状態空間モデルで記述されると仮定する。

    • $\Large \displaystyle x(k+1)=Ax(k)+bv(k)$

    $\Large \displaystyle y(k)=c^Tx(k)+w(k)$

  • ここで、「$x(k)$」は、n次元状態ベクトルである。
  • 「$A$」は、n☓n行列である。
  • 「$v(k)$」は、平均値「$0$」、分散「$\sigma_v^2$」の正規性白色雑音である。システム雑音とも呼ばれる。
  • 「$w(k)$」は、平均値「$0$」、分散「$\sigma_w^2$」の正規性白色雑音である。観測雑音とも呼ばれる。
  • 「$v(k)$」と「$w(k)$」は、お互いに独立であると仮定する。
  • 状態空間モデルの係数行列・ベクトル「$(A,b,c)$」は、既知であると仮定する。
  • この仮定は、時系列の状態空間モデルは、何等かの時系列モデリング法によって、既に得られているとする。
6.カルマンフィルタの導出2
  • まず、k=1の場合を考えて見る。「$y(1)$」を見て、「$x(1)$」を推定するために、両者の同時分布を求める必要がある。
  • 設定より、以下の式が成立する。

    • $\Large \displaystyle x(1)=x(0)+v(0)$

    $\Large \displaystyle y(1)=x(1)+w(1)=x(0)+v(0)+w(1)$

  • 行列を使うと、以下のように書ける。

    • $\Large \displaystyle \begin{pmatrix} y(1) \\ x(1) \end{pmatrix} = \begin{pmatrix} 1 & 1 & 1 \\ 0 & 1 & 1 \end{pmatrix} \begin{pmatrix} w(1) \\ v(0) \\ x(0) \end{pmatrix} $

  • 右辺のベクトル「$(w(1),v(0),x(0))^T$」は、前提から三次元正規分布に従う。
  • 従って、左辺のベクトル「$(y(1),x(1))^T$」も、二次元正規分布に従う。
7.その先
  • 線形離散時間定常確率システムは，下図のように表現される。
  • 線形離散時間定常確率システムの状態方程式は，次のようになる。

    • $\Large \displaystyle x(k+1)=Ax(k)+bv(k)$

  • 「$x(t+1)$」は、1期先の状態変数である。「$x(t)$」は、現在の状態変数である。
  • 「$v(k)$」は、「$v(k)\sim N(0,\sigma_v^2)$」のシステム雑音（神の見えざる手）である。
  • 「$A\in R^{n\times n}$」および「$b\in R^{n}$」である。
  • 線形離散時間定常確率システムの観測方程式は，次のようになる。

    • $\Large \displaystyle y(k)=c^Tx(k)+w(k)$

  • 「$x(t)$」は、現在の状態変数である。
  • 「$w(k)$」は、「$w(k)\sim N(0,\sigma_w^2)$」の観測雑音である。
  • 「$C\in R^{n}$」である。
  • カルマンフィルタでは以下の補正の式を使います。
  • 補正後の状態　＝　補正前の状態 ＋ カルマンゲイン ×（本物の観測値－予測された観測値）
  • カルマンゲイン ＝ 状態の予測誤差の分散　÷（状態の予測誤差の分散＋観測方程式のノイズの分散）