2018年5月28日
東京大学 生産技術研究所 岡本有司
1. はじめに
推定論の枠組みの中で，DNBがどのような立ち位置を締めているかを解析する．
本稿のポイントは以下の３つ:
- 時間非依存のシステムか時間依存のシステムか?
- 特徴量が未知か既知か?
- 測定した特徴量によって，現象を十分に再現できるかいかな？(十分統計量かいなか?)
2. 推定論
「推定」とは，得られた観測データから，その元となった現象の「特徴」を導くことであり，推定された特徴の値を「推定値」と呼びます．
例えば，平均 $\mu$ 分散 $\sigma^2$ の正規分布 $\mathcal{N}(\mu,\sigma^2)$ によって, $n$個の観測データ ${Y_1,\ldots ,Y_n}$ が得られた状況を考えます．
$$
{Y_1 ,\ldots , Y_n} \sim \mathcal{N}(\mu,\sigma^2) \Rightarrow \hat{\mu} := \frac{1}{n}\sum_{i = 1}^n Y_i,~\hat{\sigma}^2 := \tfrac{1}{n}\sum_{i = 1}^n (Y_i - \hat{\mu})^2
$$
この時， 観測データ ${Y_1,\ldots ,Y_n}$ から正規分布の平均 $\mu$ や分散 $\sigma^2$ を導くことを「推定」と呼びます．
また，観測データから推定する，平均値 $\hat{\mu} := \tfrac{1}{n}\sum_{i = 1}^n Y_i$ や分散値 $\hat{\sigma}^2 := \tfrac{1}{n}\sum_{i = 1}^n (Y_i - \hat{\mu})^2$ のを 推定値と呼びます．
特に，平均値と分散の組 $(\hat{\mu}, \hat{\sigma}^2)$ は元の現象を十分に再現できる特徴なので， 十分統計量 とも呼ばれます．
3. 回帰解析
説明変数 $x$ と目的変数 $y$ の間に，下記の関係が成り立つ時を考える．
$$
y = f(x) + w, \quad w\sim \mathcal{P}
$$
ここで，$\mathcal{P}$ は平均0の確率分布であり，$\mathcal{N}(0,\sigma^2)$ を仮定することが多い．
式(1)を回帰モデルとよび，得られたデータ $((x_1,y_1),\ldots, (x_N, y_N))$ から，写像 $f$ を推定することを回帰解析と呼ぶ.
特に，
$$
y = \sum_{i = 1}^m a_i \phi_i(x) + w, \quad w\sim \mathcal{N}(0,\sigma^2) \tag{1}
$$
が成り立つ時，(1)を線形回帰モデルと呼び，$a:=[a_1, \ldots , a_m]^T$を重み，$\phi:=[\phi_1,\ldots , \phi_m]^T$を特徴量と呼ぶ.
線形回帰モデルでは，得られたデータから重み$a$を推定することが目的となっている．
得られた説明変数と目的変数の組 $((x_1,y_1),\ldots, (x_N, y_N))$ に対して，
\begin{align}
\Phi :=
\begin{bmatrix}
\phi(x_1)^T \\
\vdots\\
\phi(x_N)^T
\end{bmatrix}, \quad
Y :=
\begin{bmatrix}
y_1\\
\vdots\\
y_N
\end{bmatrix}
\end{align}
と定義すると, 誤差関数
$$
L(\hat{a}) :=| Y - \Phi \hat{a} |^2
$$
を最小化する $\hat{a}$ が重み $a$ の推定量となっている．
誤差関数を最小化する推定値 $\hat{a}$ は
$$
\hat{a} = (\Phi^T \Phi)^{-1}\Phi^T Y
$$
となる．
さらに，
$$
\hat{\sigma}^2 = \frac{1}{n} L(\hat{a})
$$
が，$\sigma^2$ の推定量となる．
以上により， $\hat{a}$ と $\hat{\sigma}^2$ は線形回帰モデル(1)の十分統計量となる.
よって，新たな説明変数 $x$ に対して, 目的変数 $y$ を予測することができる．
次に回帰解析を時系列システムに応用する．
4. 時系列予測
時系列予測は, ある等時間の間隔で与えられた観測値が ${y_1, \ldots, y_D}$ が与えられた時，その未来の観測値 $y_{D+1},y_{D+2}, \ldots$ を予測する問題である．
観測値 $y_{k}$ が時刻 $y_{k-1},\ldots, y_{k-m}$ に依存し,
$$
y_{k} = \sum_{i = 0}^M a_i \phi_i (y_{k-1},\ldots, y_{k-m}) + w_k ,\quad w_k\sim \mathcal{N}(0,\sigma^2) \tag{2}
$$
というダイナミクスに従うことを仮定する．
この問題は，$y_{k-1},\ldots, y_{k-m}$ を説明変数，$y_{k+1}$ を目的変数とした回帰解析と一致している．
よって，重み $a$ と分散 $\sigma^2$ を線形回帰モデルと同様の手法を用いて推定することにより，未来の観測値を予想する．
また，上記の線形回帰モデルの推定は，時系列データ ${y_1,\ldots , y_D}$ に対して，
$$
{ y_{m+1},\ldots ,y_{D} }
$$
を目標変数のデータ，
$$
{(y_1,\ldots ,y_{m}),\ldots ,(y_{D-m},\ldots ,y_{D-1})}
$$
を説明変数のデータとすることにより，重み $a$ と分散 $\sigma^2$ を推定できる．
離散システム(2)の特徴量は，モデルが構造が既知ならば，それに適した設計が行われる．
一方，未知な現象を記述されるために，特徴量を状態の一次項で記述することがある．
ここで，
$$
M = m,\quad
\phi_0(y_1,\ldots, y_m ) = {\bf1},\quad \phi_i(y_1,\ldots, y_m ) = y_i
$$
とすると，離散システム(2)は線形離散システム
$$
y_{k} = a_0 {\bf1} + \sum_{i = 1}^{m} a_i y_{k - i} + w_k \tag{3}
$$
となり, ARモデル と呼ばれている．
ここで，
\begin{align}
A:=
\begin{bmatrix}
0 & I & 0 & \cdots & 0\\
0 & 0 & I & \cdots & 0\\
\vdots & \vdots & & \ddots & \vdots\\
0 & 0 & 0 & \cdots & I\\
a_{m} & a_{m - 1} & a_{m - 2} & \cdots & a_1
\end{bmatrix},\quad
B :=
\begin{bmatrix}
0\\
\vdots\\
0\\
1
\end{bmatrix},\quad
C := [0\cdots 0~1 ]
\end{align}
とし， 行列$A$の最大特異値 $\rho(A)$ が
$$
\rho(A)<1
$$
とする．
このとき，時刻 $k$ の状態 $x_k$ を
\begin{align}
x(k) :=
\begin{bmatrix}
y_{k-m + 1}\\
\vdots\\
y_{k}
\end{bmatrix}
+ (A - I)^{-1}
\begin{bmatrix}
0\\
\vdots\\
0\\
a_0
\end{bmatrix}
\end{align}
と定義すると，ARモデルは
$$
x(k + 1) = A x(k) + B w_k\
y_k = C x(k)
$$
と記述することができる．
さて，一般の回帰解析におけるパラメータ推定とことなり，離散システムのパラメータを推定では，システムの安定性 を保証する必要がある．
そのため，重み $a$ に制約が存在する．
ARモデルでは，行列Aの最大特異値が1より小さいならば，システムが安定となる．
特に，観測値 $y_i$ がスカラーであるARモデル(3)では，
$$
\lambda^m - \sum_{i = 1}^m a_i \lambda^{m - i} = 0
$$
を満たす $\lambda$ が
$$
|\lambda| < 1
$$
ならば，システムの安定性が保証される．
4. Dynamic Network Biomarkerの立ち位置
Dynamic Network Biomarker (DNB) は，確率システムがどれだけ分岐点に近いかを評価する指標である．
詳細は，この資料.
DNBは，モデルのパラメータすべてを推定するのではなく，観測値の共分散行列のみを推定する手法である．
この推定により，異なる分岐パラメータの値を取る確率システムに対して，どちらが分岐に近いかを評価できる．
一般の時系列予測に対するDNBの大きな差異は以下の２つ:
- 時系列予測の分野では，等間隔な時系列データを用いて目標変数のデータと説明変数のデータの組を用意する必要があったが，DNBには必要ない.
- 直接遷移行列 $A$ を推定するわけではないので，分岐点に近いかどうかを正確に調べることはできない．
- DNBでは，十分に収束している状態を仮定した上で，共分散行列と遷移行列の最大特異値の関係を示しているが，過渡状態では解析することができない．
実験対象の時系列データを正確に取得できない対象にはDNBは有効だと考えられる．
一方，分岐点周辺では状態が収束する速度は遅くなるため，時系列データが正確に取得できるならば時系列予測のほうが優位になると考えられる．
6. まとめと今後の課題
DNBは離散システムのパラメータを推定する手法であるが，十分統計量を推定しているわけではない．
実験対象の時系列データを正確に取得できない対象にはDNB有効だと考えられるが，時系列データを正確に取れるのであれば有効であるとは言えない．
今後は，DNBと一般のARモデルのパラメータ推定のどちらが精度がいいか解析する．
参考文献
推定論：
竹村 彰通，「現代数理統計学」，創文社，1991
回帰解析：
C. M. ビショップ, 「パターン認識と機械学習」，丸善出版，2002
時系列予測：
池口 徹 (著), 小室 元政 (著), 山田 泰司 (著), 合原 一幸 (編集), 「カオス時系列解析の基礎と応用」，産業図書，2000
(もう少し読みやすい本があると思います．)
DNB:
Makito Oku, Kazuyuki Aihara, On the covariance matrix of the stationary distribution of a noisy dynamical system, Nonlinear Theory and Its Applications, IEICE, Released April 01, 2018, Online ISSN 2185-4106, Print ISSN, https://doi.org/10.1587/nolta.9.166