1. はじめに

本稿では，サンプルごとに観測バイアスが存在するシステムの共分散行列を解析する．
これにより， 共分散行列をと平衡点の固有値の関係を調べる．

2.問題設定

平衡点近傍で線形化された$n$次元の離散システムが，次の用に記述できるとする．

\begin{align}
x(k+1) &= A x(k) + w(k), \\
y(k) &= x(k) + \xi \tag{1}
\end{align}

ここで，$x$は状態，$y$は観測値，$w$はシステムに付加されるノイズ，$\xi_i$はサンプル$i$ごとに定まる観測バイアスとする．

システムのノイズ$w$は，ある半正定行列$C_w$が存在し，

$$ E[w(k)]= 0,\quad E[w(k) w(k)^T] = C_w$$

を満たす．
また，観測バイアス $\xi$ は，ある半正定行列$C_w$が存在し，

$$ E[\xi] = 0,\quad E[\xi \xi^T] = C_\xi $$

を満たす．

システムのノイズ $w$ の期待値は0であり，観測バイアス $\xi$ の期待値も0であるため，状態$x(k)$ と観測値 $y(k)$ の期待値も0となる．

本文献では，観測値 $y$ の共分散行列と，行列 $A$ の固有値の関係を示す．
これにより，観測した共分散行列を元に，元のシステムが分岐する前兆を捉えることができる．

3節では状態$x$の共分散行列を解析し，4節において観測値$y$の共分散行列を解析する．

3.状態共分散行列の解析

状態$x$の共分散行列を

$$ C_x := E[x(k)x(k)^T] $$

と定義する．
このとき，$C_x$は時間に依存しないことに注意する．

状態$x$の共分散行列は時間に依存しないので，以下の関係式が成り立つ．

\begin{align}
C_x &= E\Big[x(k+1) x(k+1)^T\Big]\\
&= E\Big[(Ax(k) + w(k)) (Ax(k) + w(k))^T \Big]\\
&= AE[x(k)x(k)^T] A^T +  AE[x(k)w(k)^T] + E[w(k)x(k)^T] A^T + E[w(k)w(k)^T]\\
&= AC_xA^T + C_w
\end{align}

3行目から4行目への変換で，$x(k)$ と$w(k)$ が独立であることを用いて, $E[x(k)w(k)^T]=0$とした．
上の式は，

$$ AC_xA^T - C_x + C_w =0 $$

と書き直すことがでる．
この行列方程式は，離散システムのリアプノフ方程式と呼ばれており，行列$A$のスペクトル半径が1より小さいならば，

$$ C_x = \sum_{k=0}^\infty A^k C_w A^{kT} \tag{2}$$

と記述される．
よって，$C_w$が正定行列であれば，行列$A$のスペクトル半径が1に近づくほど，$C_x$のスメクトル半径は増大する傾向がある．

次節では，観測値の共分散行列を解析する．

4.観測値の共分散行列の解析

観測値$y$の共分散行列を，

$$ C_y := E[y(k)y(k)^T] $$

を定義すると，次の関係式が成り立つ．

\begin{align}
C_y  &= E[y(k)y(k)^T] \\
&= E[(x(k) + \xi)(x(k)+\xi)^T]\\
&= E[x(k)x(k)^T] + E[x(k)\xi^T] + E[\xi x(k)^T] + E[\xi \xi^T]\\
&= C_x + C_\xi
\end{align}

3行目から4行目への変換で，$x(k)$ と $\xi$ が独立であることを用いて，$E[x(k)\xi^T]=0$とした．
状態の共分散行列の解(2)を用いると，

$$ C_y = \sum_{k=0}^\infty A^k C_w A^{kT} + C_\xi $$

が成り立つ．
状態の共分散行列$C_x$と同様に，行列$A$のスペクトル半径が1に近づくほど，観測値の共分散行列$C_y$のスペクトル半径は増加する傾向がある．

5.共分散行列のスペクトル半径の解析

離散システム(1)の行列$A$が変化した場合おける，観測値の共分散行列のスペクトル半径の変化を考える．
変化前の行列を$A_1$，変化後の行列を$A_2$とする．
また，それぞれの状態の共分散行列を$C_{x1},C_{x2}$，観測値の共分散行列を $C_{y1}, C_{y2}$ と記述する．
この時，以下が成り立つ．

$$
| (\rho(C_{y1}) - \rho(C_{y2})) - (\rho(C_{x1}) - \rho(C_{x2})) | \leq \lambda_\max(C_\xi) - \lambda_\min(C_\xi).
$$

ここで，

$$
C_\xi = \sigma_{\xi} I_n
$$

ならば，

$$
\rho(C_{y1}) - \rho(C_{y2}) = \rho(C_{x1}) - \rho(C_{x2})
$$

が成り立つ．
そのため，
状態の共分散行列$C_x$との差は観測バイアスの共分散行列$C_\xi$のみであるため，現行のDNBの研究と同様に，観測値の共分散行列の増加するとき，システムが分岐する傾向がある．

6.例題

$n=2$ のシステムに対して，その共分散行列のスペクトル半径の解析を行なった．
$A_1$ はスペクトル半径が $0.9$ となるランダムな行列とし，$A_2$ はスペクトル半径が $0.5$ となるランダムな行列とした．
その時の，シミュレーション結果を下記に示す．

上記の図では，$\sigma_\xi$ が大きくなるほど，共分散行列のスペクトル半径が増加しているが，$\rho(C_{y1}) - \rho(C_{y2})$ には変化がないことが確認できた．

7.まとめと今後の課題

本稿では，各サンプルに観測バイアスがあるモデルを提案た．
さらに，観測バイアスがあったとしても，DNB手法には影響がないことを確認した．

今後は，実システムにより近いモデルを提案し，DNBへの影響を示す．

8.謝辞と参考

本稿のモデル(1)のアイディアは，Sysmexの岸様が提案しました．
また，共分散行列の解析手法のアイディアは，2018年4月のNOLTA特集号に掲載予定の奥先生の文献を用いています．