2018-05-18

論文読みメモ：深層自己符号化器＋混合ガウスモデルによる教師なし異常検知（その1）

論文読み

以下の論文を読みます。

Bo Zong, Qi Song, Martin Renqiang Min, Wei Cheng, Cristian Lumezanu, Daeki Cho, Haifeng Chen. Deep Autoencoding Gaussian Mixture Model for Unsupervised Anomaly Detection. International Conference on Learning Representations, 2018. https://openreview.net/forum?id=BJJLHbb0-

※ キャラクターに元ネタはないです。お気付きの点がありましたらお手数ですがコメント等にてご指摘ください。

次回：まだ

f:id:cookie-box:20180513082851p:plain:w60

多次元データの教師なし異常検知をするときは次元削減と密度推定の2段階のアプローチをするのが常套手段ですが、次元削減は次元削減で、密度推定は密度推定で学習しているために局所最適解に陥りやすい、という問題提起です。次元削減の段階で異常検知に必要な情報が失われてしまう可能性があると。

f:id:cookie-box:20180513082908p:plain:w60

そりゃ次元削減で情報が失われるかもって言われたら失われるかもしれないけど…じゃあどうするのかな？

なので次元削減と密度推定をもっと効果的に組合せたいということですが…この方針は何もこの研究が初出ということではなく、過去に既に考案されているらしいです。ただそれらの先行研究では低次元空間での特徴保存や、密度推定モデルの表現力や、密度推定の学習アルゴリズムに難があったと。

まあコンセプトが同じだったら何が優れているのか明記しなきゃだしね。

それで提案手法の深層自己符号化混合ガウスモデル（Deep Autoencoding Gaussian Mixture Model: DAGMM）ですが、このモデルでは深層自己符号化器によって低次元の特徴を生成し、また再構築時の誤差も得て、それらを混合ガウスモデルに導入すると。深層自己符号化器と混合ガウスモデルは別々に学習するのではなく、通常のEMアルゴリズムで学習するのでもなく、深層自己符号化器と混合ガウスモデルのパラメータを同時に最適化すると。

再構築時の誤差も導入するの？低次元の特徴だけではなくて？なんでだろう。

Introduction でこの DAGMM のウリが3つ紹介されているんですが、たぶんその1つ目の内容がそれを説明していますね。1つ目のウリはDAGMM は従来手法と違って入力データ中の必要な情報をちゃんと保持できる。次元削減してから密度推定するとき、密度が小さいところに異常データは位置しますが、そもそも異常なデータは次元削減しにくい＝次元削減した点から再構築したときの誤差が大きくなるものだとこの論文は主張しています。なので、この再構築エラーをも低次元の特徴に concatinate して、この新しい特徴ベクトルの空間で密度推定すると。Figure 1 で例が示されていますね。x 軸近辺のデータはほとんど再構築エラーがなかった＝適切に次元削減されたデータで、この x 軸近辺でも正常データと異常データは分離されていますが、それに加えて再構築エラーが大きいという y 軸方向にも正常データと異常データは分離されています。従来手法はこの y 軸方向の広がりを見逃していたと。でも DAGMM の第一段階「圧縮ネットワーク」は自己符号化器によって次元削減し、低次元の特徴に再構築エラーをくっつけて密度推定のプロセスに渡すと。

そっか。「入力データ中の必要な情報を保持できる」ってどういうことかと思ったけど、従来手法では失われていた情報をサルベージしたって雰囲気なのかな。必要な情報を拾い切ってやったっていうよりは。

2つ目のウリは、密度推定に混合ガウスモデルを採用する。従来手法は混合モデルでないので表現力に乏しいと。ただ、混合モデルを採用する場合、混合モデルは通常EMアルゴリズムでパラメータ推定しますが、これは次元削減と密度推定を同時に最適化するには具合が悪いです。

うーん、EステップとMステップを反復するのが厄介ってことかな。それで完結しちゃうと次元削減側に示唆がないよね。わかんないけど。じゃあどうやって混合分布のパラメータを最適化するのかな？

そこで「推定ネットワーク」を導入してそのサンプル（の特徴）が何番目の分布に帰属するかの予測を出力します。それをもとに混合分布のパラメータを直接推定します。入力データのもとで尤度が最大のパラメータを。Eステップなんて必要なかったんです。

ええ…それって自然な発想なのかな…。まあ、「特徴分布は2つのガウス分布の混合らしい」とか最初からわかってたら、各サンプルがどちらに属すだろうかっていうクラス分類をニューラルネットワークでやっちゃえっていうのはわからなくもない、のかな。あと、それで全体を学習するときって、あくまで「圧縮ネットワーク」は再構築エラーを目的関数に学習して「推定ネットワーク」はその特徴のもとでの尤度を目的関数に学習するの？それとも、全体として尤度を目的関数にする？後者だと自己符号化器部分が自己符号化器じゃなくなっちゃうけど、前者だとそれって「次元削減と密度推定の同時学習」っていえるのかな。

Introduction だけ読むと前者にはみえますが、おいおい数式が出てくればわかるでしょう。3つ目のウリは、もう散々言及されてしまいましたが、次元削減と密度推定の同時学習が実現でき、ネットワークの事前学習も不要ということですね。

（次回があれば）つづく

2018-05-14

雑記：カルマンフィルタとか何とかカルマンフィルタとか

雑記

参考文献：データ同化入門 (予測と発見の科学) | 樋口知之 |本 | 通販 | Amazon
文字の置き方と表式をだいたい上の本に準拠していますが、違うこともあります。導出の仕方は上の本とは異なります。間違っていたらご指摘ください。今回はカルマンフィルタの話しかないです。キャラクターに元ネタはないです。

次回：まだ

カルマンフィルタとその仲間たちみたいな何とかフィルタって色々ありますよね。

逐次ベイズ推定のアルゴリズムのこと？確かに色々あるね。拡張カルマンフィルタとか、粒子フィルタとか。

その辺ってなんかわかった気になってしまいませんか？線形でベイジアンだったらちゃんと解けばよくて、そうじゃなかったら現在の $x_t$ の周りで1次形近似したり、 $x_t$ の分布を粒子のアンサンブルで近似すればいいんだ、って感じで…。

なんで！？いうほどわかった気にならないよ？

というか、そもそも何をしたかったのかもよくわからなくなってしまって…。

目的は状態の推定だよね！？

よく ARIMA のような伝統的な時系列モデルから状態空間モデルへ、って説明されますけど、話が変わってる気がするんですよね。ARIMA では AR や MA の各次数の係数を求めて、状態空間モデルは状態の分布を推定していくんですよね。なんか話が逆になってませんか？モデル $f(x_t)$ に対して時系列データ $x_t$ を所与として $f$ の係数を求めたいのか、 $f$ の形は所与として $x_t$ を推定したいのか…。

うーん…状態って言葉が曖昧なのかな？観測値の真の値であるかのように説明されることも多いけど、というよりは状態って「トラッキングしたいもののベクトル」なんだよね。所与なのはあくまで観測値 $y_t$ で、ARIMA ではむしろ AR や MA の各次数の係数たちが $x_t$ で（ふつう時不変だから $x$ だけど）。状態空間モデルではモデルの形 $f$ は所与とするけどもしそのパラメータが未知ならそれは $x_t$ 側に含めると思うよ。なんていうか結局「確信しているものを所与として不確かなものを推定したい」だね。トートロジーだけど。

なるほど…。じゃあ、このブログの2018年4月23日の記事にある、「（状態空間モデルで）もし『状態』が観測される変数そのもので、『観測モデル』が恒等写像だったら（伝統的な時系列モデルと）一緒の構図になる」というのは対応がおかしかったんですね。より正しくは、「もし『状態』が ARIMA モデルの係数（及び必要なステップだけの過去から現在までのノイズや観測値）で、『観測モデル』がそれらの線形和（ARIMA モデルの形）だったら一緒の構図になる」ですね。観測モデルは観測といっても日常的な意味のように「観測時に入り混じる誤差を考慮したモデル」ではなく、「不確かな変数を知覚できる数値に変換するモデル」なのに、筆者はそこを踏み外して、思考停止でシステムモデルに ARIMA を代入していたんですね。

う、うん。たぶんそうだね…。あくまでそういう状態空間モデルにすれば ARIMA と「係数を求めたい」という目的が大雑把にそろうというだけで、実際には ARIMA モデルの方では逐次的に係数の分布を更新してはいかないと思うけどね。あと、もしそのように修正する前の元の文章の立場に立つなら、「この時系列データはこの係数のARIMAモデルで時間発展する」ということの方に絶対の確信があって、知覚された時系列データがそれに従っていないならそっちがあやしい、観測ノイズによって真の値とずれているんだろう、真の値が知りたい、という姿勢にはなるね。

不確かなものを推定したいという目的はわかったと思います。では早速各種フィルタについてちゃんと追って行きたいと思うんです。まず一般的に、一期先予測とフィルタリングは以下の式で表せます。
　　 $\displaystyle p(x_t | y_{1:t-1}) = \int p(x_t | x_{t-1}) p(x_{t-1} | y_{1:t-1}) dx_{t-1}$ 　（一期先予測分布）
　　 $\displaystyle p(x_t | y_{1:t}) = \frac{p(y_t | x_t) p(x_t | y_{1:t-1})}{p(y_t | y_{1:t-1})} = \frac{p(y_t | x_t) p(x_t | y_{1:t-1})}{\displaystyle \int p(y_t | x_t) p(x_t | y_{1:t-1}) dx_t}$ 　（フィルタ分布）

また端折ったね…。まず前提として、ここでは時刻 $t-1$ までの観測値 $y_{1:t-1}$ が得られたもとでの状態の分布 $p(x_{t-1} | y_{1:t-1})$ を所与とするんだね。かつ、時刻 $t-1$ の状態を $x_{t-1}$ と固定したもとでの次の時刻の状態の分布 $p(x_t | x_{t-1})$ も計算できるとする。これはシステムモデルを決めておきなさいってことだね。これらがわかれば、一期先予測の式は $x_{t-1}$ の分布にわたって次の時刻の $x_t$ がどうなるかの分布を積分するだけ。それで次にフィルタの式は、現時点までの観測値 $y_{1:t}$ が出そろったもとでの状態の分布 $p(x_t | y_{1:t})$ はどうなっているか、ということだけど…これは観測モデル $p(y_t | x_t)$ とさっきの一期先予測の式 $p(x_t | y_{1:t-1})$ の積からベイズの定理により求まるね。ただ、単にこれらの積だと確率分布が正規化されないから、 $p(y_t | y_{1:t-1})$ で割る必要がある。これは、 $p(y_t | y_{1:t-1})$ は $p(y_t | x_t) p(x_t | y_{1:t-1})$ を $x_t$ の分布にわたって積分すれば求まる。あと、ここまで暗に、下図のグラフィカルモデルで表される従属構造が仮定されていることに注意が必要だね。つまり、 $p(x_t | x_{1:t-1}, y_{1:t-1}) = p(x_t | x_{t-1})$ 及び $p(y_t | x_{1:t}, y_{1:t-1}) = p(y_t | x_t)$ が成り立つ。マルコフ性ともいうね。そうじゃなきゃ上の式にならない。

※ ここに状態空間モデルのグラフィカルモデルを貼る。

あ、はい、色々飛ばしてました…補足ありがとうございます。それでまずカルマンフィルタの導出です。ここからは、システムモデルと観測モデルがそれぞれ以下のようにかけると仮定しています。
　　 $\begin{cases} x_t = F_t x_{t-1} + G_t v_t \; \; \; \; & v_t \sim N_m(0, Q_t) \\ y_t = H_t x_t + w_t & w_t \sim N_l(0, R_t)\end{cases}$
上式に出てくる以下の文字は以下のような次元のベクトルや行列です。 $x_t \in \mathbb{R}^k, \; y_t \in \mathbb{R}^l, \; F_t \in \mathbb{R}^{k \times k}, \; G_t \in \mathbb{R}^{k \times m}, \; H_t \in \mathbb{R}^{l \times k}, \; Q_t \in \mathbb{R}^{m \times m}, \; R_t \in \mathbb{R}^{l \times l}$
このうち $Q_t$ と $R_t$ は正規ノイズの分散共分散行列ですね。それで、 $x_t, \; v_t, \; y_t, \; w_t$ が最初から決まっていない、確率変数です。目的は、 $x_t$ の分布をトラッキングすることです。

うん、観測値より状態の次元の方が大きいなら $H_t$ は横長の、逆なら縦長の行列になるね。

そうですね。それで、先に結論を書きます。
　　 $x_{t-1} \, | \, y_{1:t-1} \sim N_k (\mu_{t-1|t-1} , \, V_{t-1|t-1})$
を仮定すると、
　　 $x_{t} \, | \, y_{1:t-1} \sim N_k (\mu_{t|t-1} , \, V_{t|t-1}) = N_k (F_t \mu_{t-1|t-1} , \, F_t V_{t-1|t-1} {F_t}^{\top} + G_t Q_t {G_t}^{\top})$ 　（一期先予測分布）
　　 $x_{t} \, | \, y_{1:t} \sim N_k (\mu_{t|t} , \, V_{t|t}) = N_k (\mu_{t|t-1} + K_t e_t , \, V_{t|t-1} - K_t H_t V_{t|t-1})$ 　（フィルタ分布）
但し、
　　 $e_t = y_t - H_t \mu_{t|t-1} \; \in \mathbb{R}^l$
　　 $K_t = V_{t|t-1} {H_t}^{\top} (H_t V_{t|t-1} {H_t}^{\top} + R_t)^{-1} \; \in \mathbb{R}^{k \times l}$
が成り立ちます。時刻 $t-1$ までの観測値 $y_{1:t-1}$ が得られたもとで推定されている状態の分布がガウシアンならば、一期先予測もガウシアンであり、さらに時刻 $t$ のフィルタ分布もガウシアンであるということです。かつ、その一期先予測とフィルタ分布の平均ベクトルと分散共分散行列は上式のように1つ手前の平均ベクトルや分散共分散行列の式として陽に書き下せます。

有名な結論だね。

でもどうしてこういうことになるのか。

まず一期先予測には以下をつかおうか。つまり、多変量正規分布にしたがう確率ベクトルは線形変換しても多変量正規分布にしたがうし、多変量正規分布にしたがう確率ベクトルどうしを足しても多変量正規分布にしたがう。証明は、前者は $M_{Ax}(t) \equiv E(\exp(t^{\top} Ax))=M_x(A^{\top} t)$ を、後者は、互いに独立な確率ベクトル $x_1, \, x_2$ について $M_{x_1 + x_2}(t) = M_{x_1}(t) M_{x_2}(t)$ を利用すれば導出できるかな。あ、 $M_x(t)$ はモーメント母関数で、これは確率分布に対して一意に定まる関数で、多変量正規分布のモーメント母関数は $M_x(t) = \exp(t ^{\top} \mu + \frac{1}{2} t ^{\top} \Sigma t)$ だね。

$x \sim N_k (\mu_x, \, \Sigma_x)$ で $y = Ax$ のとき、 $y \sim N_l (A \mu_x, \, A \Sigma_x A^{\top})$

$x_1$ と $x_2$ が独立で、 $x_1 \sim N_k (\mu_{x_1}, \, \Sigma_{x_1}), \; x_2 \sim N_k (\mu_{x_2}, \, \Sigma_{x_2})$ のとき、 $x_1 + x_2 \sim N_k (\mu_{x_1} + \mu_{x_2}, \, \Sigma_{x_1} + \Sigma_{x_1})$

それらをつかう……？えっと、いま $x_{t-1} \, | \, y_{1:t-1}$ が多変量正規分布にしたがうと仮定しているので、その線形変換である $F_t x_{t-1} \, | \, y_{1:t-1}$ も多変量正規分布にしたがう。それとは独立にノイズ $G_t v_t$ も多変量正規分布にしたがうから… $x_{t} \, | \, y_{1:t-1} = F_t x_{t-1} \, | \, y_{1:t-1} + G_t v_t$ は多変量正規分布にしたがいますね。それで、その平均ベクトルと分散共分散行列は、
　　 $\mu_{t|t-1} \equiv {\rm E}(x_{t} \, | \, y_{1:t-1}) = {\rm E}(F_t x_{t-1} \, | \, y_{1:t-1}) + E(G_t v_t) = F_t \mu_{t-1|t-1}$
　　 $V_{t|t-1} \equiv {\rm Var}(x_{t} \, | \, y_{1:t-1}) = {\rm Var}(F_t x_{t-1} \, | \, y_{1:t-1}) + {\rm Var}(G_t v_t) = F_t V_{t-1|t-1} {F_t}^{\top} + G_t Q_t {G_t}^{\top}$

一期先予測分布が導出できたね。フィルタ分布はベイズの定理を使うとどうなる？さしあたり正規化因子は考えなくていいよ。

ベイズの定理をつかうと、
$\begin{split} p(x_{t} \, | \, y_{1:t}) &\propto p(y_t | x_t) p(x_t | y_{1:t-1}) \\ & = \exp \left( -\frac{1}{2} (y_t - H_t x_t)^{\top} {R_t}^{-1} (y_t - H_t x_t) \right) \exp \left( -\frac{1}{2} (x_t - \mu_{t|t-1} )^{\top} {V_{t|t-1}}^{-1} (x_t - \mu_{t|t-1} ) \right) \\ & = \exp \left( -\frac{1}{2} (y_t - H_t x_t)^{\top} {R_t}^{-1} (y_t - H_t x_t) -\frac{1}{2} (x_t - \mu_{t|t-1} )^{\top} {V_{t|t-1}}^{-1} (x_t - \mu_{t|t-1} ) \right) \end{split}$
ここからどうすれば…。

$\begin{split} (y_t - H_t x_t)^{\top} {R_t}^{-1} (y_t - H_t x_t) + (x_t - \mu_{t|t-1} )^{\top} {V_{t|t-1}}^{-1} (x_t - \mu_{t|t-1} ) \end{split}$
この式を $x_t$ があらわれない項は無視して展開して整理すると、

$\begin{split} & - {y_t}^{\top} {R_t}^{-1} H_t x_t - {(H_t x_t)}^{\top} {R_t}^{-1} y_t + {(H_t x_t)}^{\top} {R_t}^{-1} H_t x_t + {x_t}^{\top} {V_{t|t-1}}^{-1} x_t - {x_t}^{\top} {V_{t|t-1}}^{-1} \mu_{t|t-1} - {\mu_{t|t-1}}^{\top} {V_{t|t-1}}^{-1} x_t \\ &= {x_t}^{\top} ({H_t}^{\top} {R_t}^{-1} H_t + {V_{t|t-1}}^{-1}) x_t - {x_t}^{\top} ({H_t}^{\top} {R_t}^{-1} y_t + {V_{t|t-1}}^{-1} \mu_{t|t-1}) - ({y_t}^{\top} {R_t}^{-1} H_t - {\mu_{t|t-1}}^{\top} {V_{t|t-1}}^{-1}) x_t \end{split}$

こうなるけど、 $x_t$ に挟まれている ${H_t}^{\top} {R_t}^{-1} H_t + {V_{t|t-1}}^{-1}$ は、逆行列補題より逆行列 $V_{t|t-1} - ( V_{t|t-1} {H_t}^{\top} (H_t V_{t|t-1} {H_t}^{\top} + R_t)^{-1}) H_t V_{t|t-1}$ をもつ。加えて、 ${x_t}^{\top}$ に右側からかかっている行列と $x_t$ に左側からかかっている行列が対称だから、上の式は結局 $\begin{split} & ({x_t} - \alpha)^{\top} ({H_t}^{\top} {R_t}^{-1} H_t + {V_{t|t-1}}^{-1}) (x_t -\alpha) \end{split}$ の形に整理することができる。これを展開して係数比較すれば $\alpha$ が求まるね。

※ ここに逆行列補題をかく。

そして、 $\begin{split} & ({x_t} - \alpha)^{\top} ({H_t}^{\top} {R_t}^{-1} H_t + {V_{t|t-1}}^{-1}) (x_t -\alpha) \end{split}$ の形に整理できるということは、フィルタ分布 $p(x_{t} \, | \, y_{1:t})$ が多変量正規分布にしたがうということに他ならないよね。平均 $\alpha$ で、分散共分散行列 $({H_t}^{\top} {R_t}^{-1} H_t + {V_{t|t-1}}^{-1})^{-1}$ のね。

…ごめん寝てた。

なんで寝るの！？部長のために式展開したのに！

次のフィルタの話に移りたいです。ざっとモンテカルロ手法まで復習したく…というのも、そろそろ部活動の予算申請の時期なので、我らがベイズ統計部も膨大な計算資源の必要性を申請書にしたためなければ。

部活動予算でどんなマシン買う気なの！？運動部じゃないんだからそんな予算下りないよ？

（次回があれば）つづく

2018-04-24

論文読みメモ： Unsupervised Anomaly Detection with GAN（その1）

論文読み

以下の論文を読みます。

Thomas Schlegl, Philipp Seeböck, Sebastian M. Waldstein, Ursula Schmidt-Erfurth, Georg Langs. Unsupervised Anomaly Detection with Generative Adversarial Networks to Guide Marker Discovery. arXiv: 1703.05921, 2017. https://arxiv.org/abs/1703.05921

※ 以下、キャラクターが会話します。原作とは関係ありません。

次回：まだ

f:id:cookie-box:20180405221013p:plain:w60

この論文は、何が異常かって正解を与えることなしに、画像から異常を検知するって話みたいです！「この画像はこういう異常のデータ」って正解ラベルを付けて学習するのって、大量にラベルを付けるのが大変だし、学習データになかった異常に対応できないから、正解なしで学習したいって。GANってモデルと、画像から隠れ変数空間へのマッピング手法？を組み合わせたらしいです。あ、この論文で想定している画像というのは、OCT（光干渉断層撮影）スキャンといって、病院での検査などに用いる、皮膚や眼球の断面が撮影できる技術みたいです。網膜のOCTスキャンデータに対する検証で、病気状態の画像を異常として正しく特定できたそうです。…でも絵理さん！このモデルは、画像を入力すると、異常な画像かどうかを出力してくれるんですよね？何が異常かっていう正解を与えてないのに、どうやってこの画像は異常だろうってわかるんですか？おかしくないですか！？

f:id:cookie-box:20180405222314p:plain:w60

確かに、画像から隠れ空間へのマッピングによって異常を判定するってことは、隠れ空間内のどの座標にマッピングされたら異常っていうのがあるんだよね。じゃあGANっていうモデルは、正常な画像は正常な画像に対応する領域に、異常な画像は異常な画像に対応する領域にうつるように画像から隠れ変数空間への写像を学習するはずだけど、画像に正常か異常かのラベルが貼ってないのにそんな学習はできないんじゃないかな。

f:id:cookie-box:20180405220800p:plain:w60

…この論文では、網膜のOCTスキャン画像のうち、網膜下液という症状のない箇所の画像のみを学習データにしている。だから、モデルは正常な画像と異常な画像をごちゃ混ぜに正解なしで与えられてるわけじゃない。正常な画像のみを学習データとして、正常な画像の特徴を学んでおいて、未知の画像の異常度を判定する。

うーん…それって、教師なし学習っていうんですか？正常な画像か異常な画像かを判定するのに、正常な画像のみで学んでおきましたって…。犬の画像か猫の画像しか含まれていないデータを分類するのに、犬の画像をバッチリ学習しておきました、みたいな話ですよね。それってある意味正解ラベル付きみたいなものだし、予め学習したデータにあてはまりそうかそうでないかを判断するだけじゃないですか。

そう、予め知っているデータ側なのか、そうでない側なのかという話。でもそれを判定する方法は、自明じゃない。例えば、何かデータ間の距離を定義して、予め知っているデータに近いか判定する？未知の画像が犬か猫かを判定したいとして、予め学習した犬の画像のどれにも近くなくても、犬の画像かもしれない。予め学習したある犬の画像に近くても、猫の画像かもしれない。教師なし学習で代表的なクラスタリングは、データの中にクラスタの分布たちを見出す。異常検知もまた、与えられた正常データを表現する正常な分布を見出す。異質なものとの比較によって自分の特徴を把握できないから、後者はきっと前者よりもっと難しい。教師なし学習か教師あり学習かという区別は、いま本質的じゃない。

じゃあ、GAN ってモデルはどうやって正常データの分布を学習するんだろう。

生成モデルと識別モデルを同時に学習するって書いてありますね。ある1つのコスト関数を最小化するのではなく、生成モデルの表現力（識別器に本物と判定される確率）と識別モデルの真贋の識別力のナッシュ均衡を目指すみたいです。

（その2があれば）つづく

2018-04-23

論文読みメモ： Forecasting Stock Price with EKF

論文読み

以下の論文を読みます。

H. Haleh, B. Akbari Moghaddam, and S. Ebrahimijam. A New Approach to Forecasting Stock Price with EKF Data Fusion. International Journal of Trade, Economics and Finance, 2(2):109-114, 2011. http://www.ijtef.org/papers/87-F00046.pdf

※ 以下、キャラクターが会話します。それぞれの原作とは関係ありません。内容の誤りは本ブログ筆者に帰属します。

この論文の内容は、拡張カルマンフィルタで株価を予測する？

株価が予測できるんですか？夢がある話ですね！

何そのテンプレみたいな反応。

株式の真の価値と市場価格に乖離があるために、従来の分析では十分なパフォーマンスが出ないって問題提起してる。

市場価格とは別に株式の真の価値ってのがあるんですか？市場で出回ってる価格が株式の価値そのものなんじゃないんですか？

株式の本来の価値を測る方法っていうと、企業が将来獲得するキャッシュフローを現在価値に割り引く方法とか、配当方針から逆算する方法とかあるけど、そういうのが絶対的に株式の真の価値ってわけでもないし、「株式の真の価値と市場価格の乖離があるために上手くいかない」っていうくだりは「何か真の価値というものがあると考えた方が上手くいく」と読み換えた方がいいんじゃないかな。というかそれでちゃんと上手くいったんだよね？

テクニカル分析とファンダメンタル分析の情報を拡張カルマンフィルタで融合することで、イランのある工業関連の会社の株価について、回帰やニューラルネットよりもよい短期予測ができたとアブストに書いてある。

テクニカル分析とファンダメンタル分析？何かつかみどころがない言葉です。

この時点では株価を過去時系列に基づいて予測する手法をテクニカル分析といい、株式会社の持分であるところの株式の価値を業績予測を通して見積もる手法をファンダメンタル分析というイメージでいいと思う。大丈夫、何もつかみどころがないのは株価だけじゃない。トレーダーでなくたって、あらゆる人は明日どうなるか読めない人生を生きている。だから心配ない。

心配だよ！？

イントロは、株式投資の重要性？こうくるのはちょっと面白いかも。のっけから株価予測は投資家の関心事である、みたいなのが多いから。株価に影響するたくさんの情報を上手く扱う方法がない、だから上手く扱えるようにしたい、と。そして効率的市場仮説の紹介とそれへの反論。ここでいわれている市場の非効率性は元論文までみていないので詳しくわからない。そして既存の手法の欠点を非線形性への未対応や訓練データ不足に帰している。

訓練データが足りないって理由になるんですか？訓練データ足せばいいじゃないですか。

ここでの訓練データ不足は、訓練データが足りないために未知の状況に遭遇するといった意味みたいだから、確かに違和感があるかも。あと予測が困難な理由として、市場価格が理論モデルにしたがわないとも書いてる。理論モデルとは何だったのか。それでこの論文では拡張カルマンフィルタを用いて、非線形性への対応と、テクニカル指標とファンダメンタル指標の融合を実現したって。

そうそう絵理ちゃん、その拡張カルマンフィルタって？

この論文でも説明されているけど…ちゃんとした本を参照した方がよさそうだからとばす。

とばすの！？

111ページの Figure 2. が提案モデルの模式図だけど、本文を読み進めると、ファンダメンタル指標から配当割引モデルのゴードンモデルによって直接観測できない株価 $x_k$ を計算している。でも、テクニカル指標がどう入り込んでいるかがよくわからない。これだと、観測モデルをかませることそのものがテクニカル指標の融合だと言わんばかりのような。あ、違う。 $h(x_k, v_k)$ というのが (28) 式で、これはこれでこういう風に株価を推測する方法？なんで2次関数？

テクニカル指標に詳しくないのでよくわかりません。

ぐだぐだだけどおわり

2018-04-23

隼時系列本：ノート7

本読み

以下の本を読みます。

時系列分析と状態空間モデルの基礎: RとStanで学ぶ理論と実装
馬場真哉

プレアデス出版 2018-02-14
売り上げランキング : 7742

Amazonで詳しく見る by G-Tools

※ 以下、キャラクターが会話します。原作とは関係ありません。上の本からそれる話も多いです。誤りがあればご指摘ください。
※ アイコンまだです。

前回：ノート6 ／次回：まだ

f:id:cookie-box:20180305231302p:plain:w60

後半からはナツキが本読み手伝ってくれるんだな。よろしく。

f:id:cookie-box:20180423175223p:plain:w60

うん。ジュンは学校も仕事も忙しいから、最後まで付き合うのは難しい。だから、後は俺に任せて…。

いや俺も学校も仕事も忙しいけどな！？ …でも、ナツキって前半読んでないんだよな？いきなり後半からで大丈夫なのか？

大丈夫、この本の前半と後半は比較的別々に読めるはず…ハヤト、前半の話をまとめられる？

この本の前半の話は…まず、時系列分析っていうのは、手元の時系列データの将来値の予測をしたりその時系列データについて何かしらの知識を得ることが目的。その目的のために何をするかっていうと、時系列データには各時点での値の背後に確率分布があると仮定して、尤もらしい確率分布列を与えるモデルを特定する。そのモデルの候補として、ARIMAやVAR、GARCHとかが紹介されてた。ただ、このモデルを同定するときの注意点として、予め対象時系列データを定常にしておくことや、あてはめたモデルで説明できない残差の正規性や自己相関の確認が要るんだったな。そのために用いる統計的検定も色々出てきた。

まとめありがとう。後半の「状態空間モデル」も、目的はだいたい一緒だよ（「状態」を明示的に扱うなりの目的も加わってくるかもしれないけど）。でも、目的のために何をするかがもっと一般的になる。時系列データには各時点での値の背後に、直接観測できない「状態」と、状態をその時点まで運んできた「システムモデル」、状態を目に見える観測値の分布に写す「観測モデル」があると考えるよ。もし「状態」が観測される変数そのもので、「観測モデル」が恒等写像だったら、前半までの話と一緒の構図になるけどね。

直接観測できない「状態」があると考えるって？

本の179ページにも魚の数の例があるけど、「状態」を導入する意義の1つは、時系列データって、「観測するときにだけ混ざってくるノイズ」があるかもしれない。例えば、俺たちが持ち回りで部室の気温を毎日記録するとして、ハヤトが記録する日だけはハヤトが温度計に近づくせいで温度計の指示値が少し高めに記録されるかもしれない。でも、そのハヤトのせいで記録値が少し大きくなった分は、翌日以降の気温に影響する成分とは思えないよね。

なんていうか、それはもう俺を記録係から外した方がいいな…。

この本の前半のモデルにも誤差項はあるけど、それはその場限りの誤差ではなくて、その影響はその後ずっと（MAモデルの場合はどこかで打ち切られるけど）続いていくものになっている。そうじゃなくて、現実に何かを測るときはその場限りの誤差ってあるはず。だったら、その場限りの誤差が付加される前の「状態」の推移を考えればいい。それが状態空間モデルの意義。それに加えて、第4部の冒頭には明示的に書いてないと思うんだけど、状態って何も「観測される変数にその場限りのノイズが付加される直前のもの」とは限らなくて、「観測される変数そのものではないけど、その推移を追うと有用な量」だったり「観測できないけど興味がある量」かもしれない。状態空間モデルならそういう設計もできる。例えば、ハルナが毎日食べるドーナツの個数は直接観測できて、ハルナの毎日の空腹度は直接観測できないけど、したがう方程式がわかっているのは後者という場合とか。

だから、ハルナの空腹度がしたがう方程式が研究で明らかにされてたらおかしいだろ！ …でもナツキ、状態空間モデルでは、その状態がしたがう方程式が予めわかってないといけないってこと？前半ではさ、ARIMA だったら ARIMA の次数と係数は AIC なり何なりを指標に計算機にゴリ押しで決めさせるって感じで、「対象時系列はこんなメカニズムをもつんじゃないか」って考えるプロセスがなかった気がするんだけど…。

ARIMA を使ってみる時点である意味「ARIMA で表現できるんじゃないか」っていう考えがあるとも思うけど…状態空間モデルでは、対象時系列に対して「こんなメカニズムをもつんじゃないか」というのがもっと明示的にある場合が多いんだろうね。もちろん、データに対する強い信念なしに、色々な次数の ARIMA を試すような用法もあると思うよ？システムモデルの係数を未知パラメータにして。182ページの最初で「このあたりの手順に関しては」って別の文献を紹介しているね。でもそういうモデル探索よりは、何かこういうシステムモデルだろうという信念が予めある程度あって、直接観測できない「状態」の推移やその分散を捉えたいって向きが強いんじゃないかな。それが180ページの「過去の知見（中略）自由にモデル化できる」「（前略）モデルの解釈が容易」ってことだと思うから。俺も普段から状態空間モデルを使うわけじゃないから、よくわからないけど…。

普段から状態空間モデルを使う高校生はいないだろ！

まあそれで、この後の中心話題は未知パラメータをどうやって推定するかだね。この本では2つの異なるアプローチのそれぞれの代表として、カルマンフィルタとHMC法を大きく取り上げてる。この辺について、ちょっとプロデューサーさんが昔書いた記事を掘り出しておくね。

カルマンフィルタに関する記事は以下など。
- 状態空間モデルのまとめ絵
- R dlm をつかってみるだけ
1時点前までの観測値を得た下での現在の状態の分布は、一時点前の状態を用いて、

だけど、これに現在の観測値が手に入った下での現在の状態はベイズの定理より

になる。もし状態空間モデルのシステムと観測が線形で、システムノイズと観測ノイズがガウシアンなら（線形・ガウス状態空間モデル）、をガウス分布とすると、もガウス分布、さらにもガウス分布、…となって、1つ手前を含む形の式になって、どんどん逐次的に計算していけるんだよね。これがカルマンフィルタだね。
HMC法に関する記事。後者の記事、収束してなさすぎだけどスクリプト合ってるのかなあ…。

もっと具体的な計算手順や、これら以外の状態空間モデルの解き方とかはまた今度にする？

うーん、183～185ページに今ナツキが言ったような2つのアプローチのパラメタ推定があるってあるけどさ、なんでその2つなのかよくわからないんだよな。要するに、真面目に解くか、乱数シミュレーションで解くかって違い？でも、逐次的かどうかも違って…真面目に同時に解くとか、乱数的に逐次的に解くとかはないの？

ないとは言い切れない、けど…真面目に全期間同時に解けることなんてまずないんじゃない？逐次的に解けるような特殊なケースだから真面目に書き下せたってだけで…。あと185ページに、HMC法ではパラメタの推定を同時にやるから逐次的にはしにくいとも書いてある。

そっか、どちらかというと、特殊なケースでは真面目に解けて、副産物的に逐次解法になっている、って感じなのかな…。あとナツキ、209ページの散漫カルマンフィルタってのは？

あまり他では聞かない用語、だけど…状態の分散の初期値をものすごく大きくしておいて、状態の期待値の初期値の影響をすぐ消滅させるテクニックを指しているみたいだね。

（ノート8があれば）つづく