Time Series Analysis： ノート5章 最尤推定（その2）

本読み

以下の本の5章を読みます。キャラクターは架空のものです。私の誤りは私に帰属します。お気付きの点がありましたらコメントでご指摘いただけますと幸いです。

Time Series Analysis

Time Series Analysis

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

前回の復習ですが、ノイズが正規ホワイトノイズであるような定常 AR(1) 過程 Y_t = c + \phi Y_{t-1} + \varepsilon_t の実現値 (y_1, \cdots, y_T) の同時分布は以下になりました。この分布は「最初に出てくる y_1 はどんな分布だろう」「y_1 を観測した上で次の y_2 はどんな分布だろう」 と積み上げていけば導くことができます。

定常ガウシアン AR(1) 過程の実現値の同時分布

\displaystyle f_{Y_T, \cdots, Y_1}(y_T, \cdots, y_1 \, ; c, \phi, \sigma^2) = \frac{\sqrt{1 - \phi^2}}{\left( \sqrt{2 \pi \sigma^2} \right)^T} \exp \left( - \sum_{t=2}^T \frac{ \bigl( y_t - c - \phi y_{t-1} \bigr)^2}{2\sigma^2} - \frac{ \bigl( y_1 - c/ (1 - \phi) \bigr)^2}{2\sigma^2/ (1 - \phi^2)}\right)
実現値 (y_1, \cdots, y_T) からパラメータ (c, \phi, \sigma^2) を推定するには上の同時分布が最大になるような (c, \phi, \sigma^2) を選べばいいわけです（実際にはこれの対数を取って最大化しますが）。これを最尤推定といいます。

119~122ページでは上の同時分布を求める別解として、「y \equiv (y_1, \cdots, y_T)^{\top} はどんな多変量正規分布にしたがうだろうか？」というアプローチで同じ同時分布を導きました。

122~123ページでは「y_1 が所与」という考えやすい場合で実際に最尤推定値を求めますが、前回は流してしまったので改めて追ってみます。y_1 が所与なら上の同時分布は以下のように化けますね。

定常ガウシアン AR(1) 過程の実現値の同時分布（y_1 が所与の場合）

\displaystyle f_{Y_T, \cdots, Y_1}(y_T, \cdots, y_2 \, ; y_1, c, \phi, \sigma^2) = \frac{1}{\left( \sqrt{2 \pi \sigma^2} \right)^{T-1}} \exp \left( - \sum_{t=2}^T \frac{ \bigl( y_t - c - \phi y_{t-1} \bigr)^2}{2\sigma^2} \right)
対数を取るとこうですね。
\displaystyle \mathcal{L}(c, \phi, \sigma^2) \equiv \log f_{Y_T, \cdots, Y_1}(y_T, \cdots, y_2 \, ; y_1, c, \phi, \sigma^2) = - \frac{T-1}{2} \log \left( 2 \pi \sigma^2 \right) - \sum_{t=2}^T \frac{ \bigl( y_t - c - \phi y_{t-1} \bigr)^2}{2\sigma^2}
これを (c, \phi, \sigma^2) について最大化したいわけですから、これらで偏微分してみましょう。
\displaystyle \frac{\partial \mathcal{L}(c, \phi, \sigma^2)}{\partial c} = \sum_{t=2}^T \frac{ 2 \bigl( y_t - c - \phi y_{t-1} \bigr)}{2\sigma^2}
\displaystyle \frac{\partial \mathcal{L}(c, \phi, \sigma^2)}{\partial \phi} = \sum_{t=2}^T \frac{ 2 y_{t-1} \bigl( y_t - c - \phi y_{t-1} \bigr)}{2\sigma^2}
\displaystyle \frac{\partial \mathcal{L}(c, \phi, \sigma^2)}{\partial \sigma^2} = - \frac{T-1}{2\sigma^2} + \sum_{t=2}^T \frac{ \bigl( y_t - c - \phi y_{t-1} \bigr)^2}{2\sigma^4}
これらの式をイコールゼロとしたものを満たす (c, \phi, \sigma^2) が最尤推定値 (\hat{c}, \hat{\phi}, \hat{\sigma}^2) になりますが、まず上の2式より (\hat{c}, \hat{\phi}) が導かれ、次に最後の式より \hat{\sigma}^2 が導かれますね。

y_1 が所与の定常ガウシアン AR(1) 過程の最尤推定値 (\hat{c}, \hat{\phi}, \hat{\sigma}^2)

\displaystyle \left(\begin{array}{cc} T-1 & \sum_{t=2}^T y_{t-1} \\ \sum_{t=2}^T y_{t-1} & \sum_{t=2}^T y_{t-1}^2 \end{array}\right) \left(\begin{array}{c} \hat{c} \\ \hat{\phi} \end{array}\right) = \left(\begin{array}{c} \sum_{t=2}^T y_{t-1} \\ \sum_{t=2}^T y_{t-1}y_t \end{array}\right)
\displaystyle \hat{\sigma}^2 = \frac{1}{T-1} \sum_{t=2}^T \bigl( y_t - \hat{c} - \hat{\phi} y_{t-1} \bigr)^2
y_1 が所与なら最尤推定値を求めることは容易ですね。ただ、「言い換えると、これは最小2乗回帰（OLS）である」とありますが、どういうことなんでしょう？

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

いま、以下の前者の問からその最尤推定値にたどり着いたよね。でもそれって、以下の後者の問から出発したのと同じになってるんだよね。

  • y_1 が所与の定常ガウシアン AR(1) 過程の、実現値 (y_2, \cdots, y_T) に対する最尤推定値は何だろうか？
  • y_1 が所与の AR(1) 過程の、実現値 (y_2, \cdots, y_T) に対する最小2乗推定値は何だろうか？
だって (\hat{c}, \hat{\phi}) を求めるのに \bigl( y_t - c - \phi y_{t-1} \bigr)^2 を最小化しているけど、これは各 y_t を c + \phi y_{t-1} で予測したときの2乗誤差に他ならないし、\hat{\sigma}^2 は残差の分散になっているしね。ちなみに、後者は定常性とガウス性を仮定していない。「誤差2乗和を最小化し、残差の分散を求める」という作業は定常性もガウス性も要求しないからね。

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

あっと、最小2乗推定と等価ということはわかりました。しかし、なら最初から最小2乗推定を紹介していただければよかったのではありませんか？ 定常性とガウス性を仮定して頑張って尤度関数を求めて、「実は単に最小2乗推定するのと同じだったんだけどね」では私が報われなくありませんか？

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

そんなことないよ。確かに最小2乗法はいつでもできるよ。最小2乗推定値自体はいつでも求まると思う。でも、それがまともな推定値になっていることを支持する材料ってないからね。定常性とガウス性をもつ場合には最小2乗推定値は最尤推定値になっていることがいまわかったけど、逆にいうとそうでない場合は変な値になっているかもしれないってことだからね。まあそう言い出すと「なら最尤推定値ならまともなのか」って話になるけど、単なる最小2乗推定値よりはまともだと信じられる場面が多いってことだと思ってる。

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

な、なるほど…？ それで、123ページからは AR(p) 過程の場合ですね。AR(p) 過程の同時分布は…

  • まず「最初の p 点 (y_1, \cdots, y_p)^{\top} はどんな多変量正規分布にしたがうだろうか？」を考えた上で、
  • その後「(y_1, \cdots, y_p) を観測した上で次の y_{p+1} はどんな分布だろう」を積み上げていっています。
AR(1) 過程でみた2つアプローチの合わせ技をやっていますね。素朴なアプローチではできないんでしょうか？

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

…5.3.8 式に AR(2) の場合が具体的に書き下してあるから、AR(2) で求めて検算してみればいいんじゃないかな…と思ったけど、素朴にはできない気がするかな。まず、何も観測していない下での最初の1点の分布 f_{Y_1}(y_1) は AR(1) のときと一緒だよね。期待値と分散は AR(2) のそれにしなければならないけど。あと、直近の p 点 (y_1, y_2) を観測した下での y_3 の分布 f_{Y_3 | Y_2, Y_1}(y_3 \, | y_2, y_1) も AR(1) のときと同じ要領だね。

\displaystyle f_{Y_1}(y_1 \, ; c, \phi_1, \phi_2, \sigma^2) = \frac{1}{\sqrt{2 \pi \sigma^2/ (1 - \phi_1^2 - \phi_2^2)}} \exp \left( - \frac{ \bigl( y_1 - c/ (1 - \phi_1 - \phi_2) \bigr)^2}{2\sigma^2/ (1 - \phi_1^2 - \phi_2^2)}\right)
\displaystyle f_{Y_3 | Y_2, Y_1}(y_3 \, | y_2, y_1 \, ; c, \phi_1, \phi_2, \sigma^2) = \frac{1}{\sqrt{2 \pi \sigma^2}} \exp \left( - \frac{ \bigl( y_3 - c - \phi_1 y_2 - \phi_2 y_1 \bigr)^2}{2\sigma^2}\right)
でもこれらの中間の f_{Y_2 | Y_1}(y_2 \, | y_1) がこんな感じじゃ求められない気がするんだよね。だって、「何も観測していない下での」でもないし、「直近の p 点を観測した下での」でもないから、どっちのアプローチも使えない。結局 AR(2) の共分散構造を満たすように (y_1, y_2) の同時分布を真面目に求めないといけない気がする。その分散共分散行列が125ページの一番上の式だと思う。

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

素朴にはいかなさそうなんですね。…あれ、126ページには「ガウス性をもたない場合の最尤推定」とありますね。5章ではガウス性をもたない場合は扱わないはずでは？

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

…ここにかいてあるのは、ガウス性をもたない過程にガウス性をもつと仮定した最尤推定をすることを「準最尤推定」ということ。生データがガウス性をもたない場合でもシンプルな変換でガウス性をもつようにできることがあること。とかかな。その変換の例が 126ページの下の方にある Box-Cox 変換だね。

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

この変換は \lambda というパラメータを与えて用いるんですね。どのように決定すればいいんでしょうか。

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

無論ガウス分布に近くなるように決めるんだけど、scipy.stats.boxcox は \lambda を与えなければ勝手に決めてくれるって。

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

文明の利器…127ページの一番上には非ガウスの場合の ARMA 過程の推定に関する研究が紹介されていますね。この本自体が17年前の本なのでそれを考慮して参考にしなければなりませんが…。5.4~5.6節は MA 過程、ARMA 過程の話ですが一旦とばします。5.6 節の最後には、ガウシアン ARMA 過程を最尤推定するのに最もシンプルなアプローチはカルマンフィルタであるとありますね。また、他の参考文献も紹介されています。

それで、尤度 \mathcal{L}(\theta) を出す式が特定できた下で、それを最大化する \hat{\theta} を如何に求めるかという話が残っていますが（y_1 が所与の場合の最小2乗法はやりましたが）、5.7 節は数値的に最大化する方法ですね。紹介されているのは以下でしょうか。

  • グリッドサーチ
  • 最急降下法
  • Newton-Raphson 法
  • Davidon-Fletcher-Powell 法（準 Newton 法）
それで、5.8 節が最尤推定を用いた統計的検定で、5.9 節の「不等式制約」というのは、パラメータ探索は「過程が定常である」「分散は正である」のような制約の下で行わなければならないので、そのために使えるテクニックですかね。

つづいたらつづく

論文読みメモ： Trellis Networks for Sequence Modeling（その1）

論文読み

以下の論文を読みます。

Shaojie Bai, J. Z. Kolter, V. Koltun. Trellis Networks for Sequence Modeling. In 7th International Conference on Learning Representations, ICLR 2019, 2019. [1810.06682] Trellis Networks for Sequence Modeling
※ キャラクターは架空のものです。私の誤りは私に帰属します。お気付きの点がありましたらコメント等で指摘ください。
f:id:cookie-box:20190101155733p:plain:w60

これは TCN を提案された方々の論文ですが、この論文では系列データを学習する新しいモデル構造である Trellis Networks（「格子垣」ネットワーク）を提案するようですが、この TrellisNet はアブストラクトによると

  • ちょっと特殊な Temporal Convolutional Networks (TCN) ともいえるし、
  • ちょっと特殊な Recurrent Neural Networks (RNN) ともいえる
とのことですね…しかし、RNN と TCN はどちらも系列データを学習するモデルですがその仕組みは全く異なりますよね？ RNN は現時点の入力と前の時点の出力から現時点の出力をつくりますし、TCN は現時点までの入力たちをどんどん畳込んで現時点の出力をつくりますし、そもそも下図をみても、どちらもお世辞にも「格子」にはみえないと思います。

f:id:cookie-box:20201010183627p:plain:w570
f:id:cookie-box:20190101160814p:plain:w60

じゃあ RNN のようでもあって、TCN のようでもあって、格子にみえるモデルを考えてみたらいいんじゃないかな。そもそもどんなモデルだったら格子にみえるだろう？

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

どんなって、「格子」というからには…あれ、格子ってなんでしょう…あ、そうです！ 複素数平面の格子点です！ つまり、格子点のように、何か横方向のインデックスと縦方向のインデックスをもつものがあれば、格子にみえるでしょう。むしろそれが格子の定義かもしれませんが…まあいいです。

それで、RNN や TCN には縦方向と横方向のインデックスをもつものがあるかという話です。この場合、横方向のインデックスは時間ステップでしょう。縦方向のインデックスは何層目の特徴かといったようなものが考えられると思います。入力を0層目、出力を1層目とすれば RNN も TCN も格子とはいえますね。でも、それではちょっとお粗末だと思います。それでは格子というより梯子ですから。

…あっ、でも、RNN も TCN も構造上2つ以上積み重ねることができますね。そうすれば、時刻 t の i 層目の出力が格子にみえると思います。

…いや、よく考えると、RNN に TCN を重ねることも、TCN に RNN を重ねることも妨げられるものではないですね。ステップ毎に出力を出すタイプの系列モデルであれば中身がどうあっても積み重ねらることができますから。であれば、TCN → RNN → TCN → RNN のような積み重ねだって可能ですし、かつ格子にみえます。TrellisNet の正体はもうこれでは？

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

いい線なんじゃないかな。でも TCN → RNN → TCN → RNN というモデルではないみたいだね。それよりは TCN → TCN → TCN → TCN が近いと思う。3ページ目の Figure 1 をみる限りは。たぶん TrellisNet には純粋な再帰路はない。

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

うおおおい！ 私には考えさせておいて自分はちゃっかりカンニングしてるじゃないですか！ それに、TCN → TCN → TCN → TCN だと RNN 要素なくないですか？ TrellisNet はちょっと特殊な RNN なんでしょう？

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

どう RNN なのか解説してくれるのが 4.2 節なんだと思う。

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

…おとなしくイントロダクションから読んでいきます。って、最初の段落、どばっと文献の引用がありますが、個々の内容には何も言及されていないですね…。近年の系列データの学習の研究では、RNN も改良が続けられているし、TCN も高いパフォーマンスを発揮しているし、セルフアテンションに基づくネットワークも快進撃をみせているということですが…セルフアテンションって RNN や TCN の並びに入るんですか？ RNN や TCN は将来のデータをカンニングしてはならないという縛りがあるはずですが、文章に対する品詞分類や翻訳などのタスクでは、ふつう将来のデータ（単語）をカンニングしてもいいですよね？ むしろカンニングしろって感じですが…。

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

セルフアテンションが RNN や TCN の並びに入るかということについては、セルフアテンション層は未来の単語も過去の単語も参照するということはしていないと思う。だから分の途中まででも出力は同じだよ。Positional Encoding とかあるから文章中の座標は単語に貼り付けてある必要はあるし、ネットワークを学習するときに次の単語を説明しようとしたり前の単語の説明しようとしたりということをしないということではないよ。

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

なるほど？ …イントロダクションの3段落目には、アブストラクトにあったように、TrellisNet はちょっと特殊な TCN であるということがもう少し詳しくかかれているようですね。

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

そこを読む限り、TrellisNet は TCN → TCN → TCN → TCN （多層 TCN）なんだけど、全ての TCN は重みを共有しているし、どの層にも入力データが injection されるんだね（右図）。なるほど格子にみえる。z_t^{(0)} はゼロベクトルね。

でも「これをみせられてもよい性能をもつとは思われないでしょう」と。まあこれだけぱっとみせられても何がいいのか、全ての TCN が重みを共有している意義は何なのかわからないね。でもこれが記憶をトランケートされた RNN の一般化になっていることをこの論文で示すらしい。それが優れている根拠になっているといいたいのかな。

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

「RNN とみなせるから優れているのだ」ですか。RNN に一定の信頼があるのは確かでしょうが。イントロダクションの5段落目にあるのは TrellisNet の性能を検証したタスクですね。おそらく先の先行研究にならって選定されていると思われます。個々のタスクの出自やデータサイズなんかの詳細は Appendix C にあるようです。

  • 文字レベル Penn Treebank（PTB）
  • 単語レベル Penn Treebank（PTB）
  • 単語レベル WikiText-103（WT103） ― 大規模データとして採用されているようですね。
  • sequential MNIST
  • permuted MNIST
  • sequential CIFAR-10
そしてこのタスクのどれもで SOTA を達成したようです（厳密には WT103 については同時期に出た論文の transformer モデルに劣っているようですが）。これまでそれぞれ異なるモデルが SOTA を達成していた6タスクを TrellisNet が一気に塗り替えたのだというのが強調されていますね。

ところで sequential MNIST 以下はもう「系列モデルに対するスタンダードなストレステスト」とされていますが、手書き数字や写真って本来全てのピクセルを同時に認識でき、縦横のつながりも明らかなものですよね。明らかに CNN さんに学習させるべきですよね。それを出鱈目な順序でピクセルごとに流すような、無駄な曲芸のようなタスクを系列モデルにさせるのってパワハラに当たらないんですか？

TrellisNet
f:id:cookie-box:20201011140713p:plain:w200
f:id:cookie-box:20190101160814p:plain:w60

TrellisNet（乙）は私たち（甲）との間に「乙は甲に自然に系列としての構造をもつデータを学習する労務を提供すること」などという雇用契約は結んでいないからパワハラには当たらないね。雇用契約にない仕事をさせるのが問題になるからね。

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

いや法律的な見解ではなく、本来は2次元としての構造をもつデータがベンチマークとなっているのはどうなのかといいたかったんですが…。2節の「背景」は短いですね。かいてある内容としては、

  • 再帰的ニューラルネットって自然言語以外にも色々なドメインに利用されているよね（という感じでまたどばっと引用だけありますね）。
  • 畳込みニューラルネットも長期記憶が必要なタスクで優れているんだよね（ここでこの方たちの TCN を含んでどばっと引用がありますね）。
  • self-attention もあるよね（どばっと略）。
  • 今回は RNN と TCN の架け橋となるモデルを提案するよ。
といった感じでまあ何にせよ RNN と TCN を組み合わせるといった感じですが、RNN と TCN を組み合わせた先行研究も引用されていますね。ただどのようなものだったのか詳しくわかりませんね。TrellisNet の方が両者をより深く結びつけたのだということではあるんでしょうが。

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

「convolutional LSTMs combine convolutional and recurrent units」というのは RNN ユニットと CNN ユニットを組み合わせたというようにみえるけどね。なら RNN の結果を CNN したり、CNN の結果を RNN したり、それこそ TCN → RNN → TCN → RNN のようなものなのかな。知らないけど。

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

3節の SEQUENCE MODELING AND TRELLIS NETWORKS に進みます。冒頭では「系列モデルとは長さ T の系列を受け取って長さ T の系列を出すものであって t 番目の出力は未来の入力に依存してはいけないものとする」という TCN の論文同様の定義が示されています。その後にはもう TrellisNet が以下のように定式化されていますね。

  • 系列モデルの入力を x_{1:T} = (x_1, \cdots, x_T) \in \mathbb{R}^{T \times p} とする。
  • 時刻 t の i 層目の出力を z_{t}^{(i)} \in \mathbb{R}^{q} とする（つまり何層目の出力であっても q 次元である）。 z_{t}^{(0)} はゼロベクトルとする。
  • TrellisNet の z_{t+1}^{(i + 1)} \in \mathbb{R}^{q} のつくり方としては、まず \tilde{z}_{t+1}^{(i + 1)} \in \mathbb{R}^{r} を下図のように線形結合で生成し（Conv1D とも解釈できますね）、 \tilde{z}_{t+1}^{(i + 1)} と z_{t}^{(i)} に適当な非線形関数 f を適用して z_{t+1}^{(i + 1)} とする。
f:id:cookie-box:20201011173741p:plain:w360
W_1 と W_2 はすべての i 層目で共通ですがら、x_t と x_{t+1} を変換した \tilde{x}_{t+1} = W_1^{x} x_t + W_2^{x} x_{t+1} は使い回せるとありますね（W_j^x \in \mathbb{R}^{r \times p}）。

…いっていることはわかりますが、腑に落ちませんね。特に最終出力 z_{t+1}^{(i + 1)} = f \bigl( \tilde{z}_{t+1}^{(i + 1)}, z_{t}^{(i)} \bigr) について、なぜ z_{t}^{(i)} を使うのか、f は具体的にどういった関数なのかという点です。f がわからなければ TrellisNet が何をするモデルなのかイメージすることは難しいと思います。もっといえばすべての層でウェイトを共有する理由だってわかりません。本文中には Vogel & Pock (2017) が画像処理で層方向のウェイトの共有をやっているとありますがだから何やら。

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

3節の最後に LSTM セルに基づく活性化をするとあるよ。4節にその理論的根拠があるとあるけど、先に活性化の式がある 5.1 節をみてみようか。気になるしね。

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

LSTM で活性化？ LSTM は活性化関数ではないと思うんですが…5.1 節をみると、ゲート付き活性化というのがそもそもあるんですね？ 5.1 節を読む限り、その LSTM に基づくゲート付き活性化というのは以下になるはずです。

f:id:cookie-box:20201011230518p:plain:w560
えっと、これは何なんでしょう。何やら LSTM めいてはいますが…。

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

…もし「前回の入力」と「前回の前の層の出力」を使用しないなら、これは時刻 i から時刻 i+1 に進む LSTM だよね。ただし、実際には i とは何層目かであって時刻を表すインデックスではないけど。それに、実際には「前回の入力」と「前回の前の層の出力」が混ぜ込まれるし、いわば「入力とゲートに混ぜ物があり、他人の記憶を植え付けられた LSTM」かな。微塵も LSTM ではないな…。

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

読み飛ばした 4 節に戻りますね。4.1 節は、TrellisNet と TCN の共通点と相違点について、「各層で未来のデータを参照せず畳込みを行うこと」「多層化するほど参照範囲が広がること」の2点は TCN と同様ですが、「各層の重みが共有されていること」「各層に入力データがインジェクションされること」は TCN と異なるといっていますね。それはまあそうですが。「各層の重みが共有されていること」は学習の安定化とモデルサイズの節約に寄与する正則化になっているとのことです。

その2があればつづく

雑記： NeurIPS 2020 pre-proceedings の「シーケンス」を含むタイトル（※）

雑記

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

NeurIPS 2020 pre-proceedings をみるとタイトルに sequence/sequential を含む論文が…23件もありますね。多いので、特に「学習対象が系列データである」といえるものに絞りましょう（※ このセリフの最下部にある12件をスキップします）。それでも11件ありますね…。以下にあやしい理解をメモしておきます。

  • 粒子フィルタによる一般化ベイズフィルタリング ― βダイバージェンスを組み込んでモデルや観測の誤りにもロバストなフィルタリングを実現した。
  • 出力間の関連性を考慮した Sequence to Multi-Sequence モデル ― 出力間の関連性を考慮することで従来より高精度な音声分離を実現した。
  • Funnel-Transformer ― エンコーダ内で入力系列を徐々に短い系列に圧縮する「漏斗」transformer を導入し、入力系列を一つのベクトルに埋め込むタスクで transformer より高性能を発揮した。
  • ペアワイズシーケンスアラインメントのための Rank-One モデルの適応的学習 ― を考案した。
  • COT-GAN ― 「因果最適輸送距離（COT）」を導入し時系列を生成するGANを実現した。
  • 2つのBERTをアダプターでつないだ seq2seq モデル ― で効果的なファインチューニングを実現し、機械翻訳タスクで推論速度を半減させつつ SOTA と同等のスコアを達成した。
  • Temporal Spike Sequence Learning Backpropagation（TSSL-BP） ― スパイキングニューラルネットを精度よく学習するための新しい誤差逆伝播法を考案した（系列データの学習ではない）。
  • スパイク列を検出するための点過程モデル ― を考案した。
  • 長文学習のための Big Bird ― 普遍的な表現力を維持しつつアテンションを疎にすることで同じ計算資源で8倍の長さの系列を学習できるようにした。
  • 解釈可能な Covid-19 予測 ― SEIR のような感染症モデルに解釈可能なままに機械学習を導入した。
  • ユーザごとの特徴を考慮したイベントの種別と発生タイミングの予測 ― を考案した。
タイトルを sequence/sequential で検索すると transformer（BERT）の論文が中途半端にヒットしますが、これらを含む論文自体はさらにたくさんあることに留意が必要です。
Generalised Bayesian Filtering via Sequential Monte Carlo
Ayman Boustati, Omer Deniz Akyildiz, Theodoros Damoulas, Adam Johansen
状態空間モデルで分布がガウシアンでない一般的な場合であって、モデルが必ずしも正しく特定されていない場合（？）に、如何にフィルタリングするかという話なのでしょうか？ 通常の逐次モンテカルロ法（粒子フィルタ）を何か改良しているのではないかと思いますが、ベータダイバージェンスとは何でしょうか。以下をみるとKLダイバージェンスの拡張のようですね。
Sequence to Multi-Sequence Learning via Conditional Chain Mapping for Mixture Signals
Jing Shi, Xuankai Chang, Pengcheng Guo, Shinji Watanabe, Yusuke Fujita, Jiaming Xu, Bo Xu, Lei Xie
seq2seq モデルならぬ、1つの系列を複数系列にマッピングするモデルでしょうか。AさんとBさんとCさんの会話が重なった音声データから、Aさんの声、Bさんの声、Cさんの声を分離するイメージのようですね。出力系列間の関連性を明示的にモデリングするということですが…Aさんの声を解読できたら、もとの音声データとAさんの声を元にBさんの声を出力するということですかね。それはそうした方がいいですよね。効果的な停止則ももっているようですが、この音声データにはもうしゃべっている人はいないよ、というのを何らかの方法で検知するのでしょうね。
後から大量のラベルなしデータで効果的にチューニングできるような、スケーラブルな言語モデルが望まれているということでしょうか。それはそうですよね。ここでは特に文章を単一のベクトルにエンコーディングするタスクを考えていると。そして、スケーラビリティのためにわざと冗長性をもたせるのでしょうか。隠れ状態の系列を徐々に圧縮していくということですか。funnel とは「漏斗」なんですね。エンコーダ層内で1層を経るごとに16単語→8単語→4単語になっていくんでしょうか。それで空いた計算資源を deeper なモデルや wider なモデルの学習に再投資すると。確かにメモリの節約になりますが再投資の仕方もアルゴリズムに含まれているんでしょうか。とにかく系列に対して1つの判定を下すようなモデルで通常の transformer より高性能であったということです。
シーケンスアラインメントとは2系列の類似箇所を検出する作業なのでしょうか。計算が高コストなんですね。そして、実際には多数のリード（調査対象の塩基配列？の1本？1本をリードというのでしょうか）のうちどの2本が似ているかにしか興味がないんですかね。それで提案手法では Rank-One モデルと多腕バンディットアルゴリズムを導入したということですが。アブストラクトにクラウドソーシングと出てきていますが、クラウドソーシングでも複数の人のアノテーションを統合しますが、Rank-One モデルも複数のセグメントの類似性を統合するのでクラウドソーシングの問題設定で用いられているんでしょうか。
COT-GAN: Generating Sequential Data via Causal Optimal Transport
Tianlin Xu, Wenliang Le, Michael Munn, Beatrice Acciaio
系列データを生成するGANということです。GAN（WGAN）は訓練データの分布と生成データの分布の Wasserstein 距離（最適輸送距離）を小さくしようとするんですよね。しかし、生成するデータが系列データの場合にどうやって分布を似せるべきかは自明でないですね。そこで「因果最適輸送距離（Causal Optimal Transport: COT）」を導入したということですが、どのようなものかはアブストラクトから推測することは難しいですね…。学習のアルゴリズムには Sinkhorn アルゴリズムという Wasserstein 距離を近似するアルゴリズムを適用しているんですね。
Incorporating BERT into Parallel Sequence Decoding with Adapters
Junliang Guo, Zhirui Zhang, Linli Xu, Hao-Ran Wei, Boxing Chen, Enhong Chen
BERT を如何に効果的に seq2seq モデルに組み込むかは自明ではないと。そこで本手法では2つのBERTモデルをエンコーダとデコーダとして採用し、その間に軽量な「アダプター」（これをタスクに応じてチューニングする）を挟む構成をとると。"catastrophic forgetting" というのはファインチューニング時に忘れるべきでない特徴抽出まで全て忘れてしまうことをいうんでしょうか。とりあえずそれを回避できるようです。機械翻訳タスクで推論速度を2倍にしながらアブストラクトに記述の BLEU スコアを達成したようですが、これはよいスコアなんでしょうか…あ、最後に SOTA と同等とありますね。…関係ないですが、英独翻訳と英仏翻訳の BLEU スコアに随分差がありますが、英独翻訳は英仏翻訳よりそんなに困難なものなのでしょうか。
そもそも Spiking neural networks（SNNs）をよく知らないですね…ニューラルネットワークだがスパイクを伝達させてさらに神経細胞に寄せたものなんでしょうか。しかし既存の誤差逆伝播法ではスパイクの不連続性を適切に扱えないと。また、まともな性能を達成するのに何ステップも要するのでスケーラブルでないと。なので Temporal Spike Sequence Learning Backpropagation（TSSL-BP）なる新しい誤差逆伝播法を考案したという話でしょうか。…これ系列データの学習の話ではなかったですね。
Point process models for sequence detection in high-dimensional neural spike trains
Alex Williams, Anthony Degleris, Yixin Wang, Scott Linderman
さっきの論文と紛らわしいですがこれは高次元のスパイク列が学習対象ですね（動物の神経データが学習データであるようです）。教師なしでパターンを発見するという話でしょうか。既存手法は時間ステップが離散化されていることが必要であるとか色々好ましくない点があるんですね。なので点過程モデルを提案したということですが、どのようなモデルなんでしょうか。
Big Bird: Transformers for Longer Sequences
Manzil Zaheer, Guru Guruganesh, Kumar Avinava Dubey, Joshua Ainslie, Chris Alberti, Santiago Ontanon, Philip Pham, Anirudh Ravula, Qifan Wang, Li Yang, Amr Ahmed
長い系列（文章）のための BERT、ですよね。BERT はアテンションがあるために系列長の2乗のメモリを消費する、ということでしょうか。そうなんでしたっけ？ まあそうとして、それを回避する疎なアテンションを導入したということですね。それでいて表現力をもつと。それで同じ計算資源で従来の8倍の長さの系列を扱えるようになったということです。長文を扱えるようになるのは可能性が広がりますね。
Interpretable Sequence Learning for Covid-19 Forecasting
Sercan Arik, Chun-Liang Li, Jinsung Yoon, Rajarishi Sinha, Arkady Epshteyn, Long Le, Vikas Menon, Shashank Singh, Leyou Zhang, Martin Nikoltchev, Yash Sonthalia, Hootan Nakhost, Elli Kanal, Tomas Pfister
感染症数理モデル（SEIR）のようなモデルに機械学習を組み込んで性能向上させかつ解釈は可能なようにしたということですかね。説明変数のエンコード、コンパートメント毎の時間変化、というコンポーネントを明示的に用意しているのが肝なのでしょうか。政策をサポートするモデルは特に解釈性が重要ということですね。アメリカの州単位や郡単位で検証しているようですね。
User-Dependent Neural Sequence Models for Continuous-Time Event Data
Alex Boyd, Robert Bamler, Stephan Mandt, Padhraic Smyth
ここでは金融商品の取引データや医学的データが意識されているのでしょうか。様々なタイプの出来事がいつ起きるかを予測するのは困難と。現在の SOTA は強度関数（いまこそこのイベントが起きるぞという強度を出力するのでしょうね）を RNN でモデリングしたものであると。しかし既存のモデルは"come from the same data distribution"？ これは色々な性格の個々人がいることをきちんと考慮していないという意味でしょうか…。提案手法はユーザ行動のモデルを償却変分推論で学習するんですね。
※ 以下の12件をスキップします。sequential だと強化学習の論文が多くヒットしますね…。

つづかない

雑記： NeurIPS 2020 pre-proceedings の「時系列」を含むタイトル

雑記

キャラクターは架空のものです。何かありましたらご指摘いただけますと幸いです。
参考文献

関連記事

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

NeurIPS 2020 pre-proceedings をみるとタイトルに time series を含む論文が…9件もありますね。上から順にメモしていきますか…面倒ですが…。なお、現時点でクリック先で閲覧できるアブストラクトを一通り読んだのみの妄想的な理解であることに留意ください。先に手短にまとめておきます。

  • ODE-LSTMs ― 不定期的サンプルに強い ODE-RNNs を改良し長期依存性を学習できるようにした。
  • Adversarial Sparse Transformer（AST） ― ジェネレータで疎なアテンションを生成しディスクリミネータで予測性能を高めることで長期予測しても誤差が蓄積しないモデルを実現した。
  • Deep Rao-Blackwellised Particle Filters ― 提案分布の factorization と decoder を組み込んだディープなラオ・ブラックウェル化粒子フィルタを提案しスイッチングするシステムを予測できるようにした。
  • 自己相関をもつ時系列どうしの因果関係を高感度に発見する手法 ― を提案した。
  • NeuralCDE ― ODE モデルの途中経過を後から修正できるようにした。
  • temporal saliency rescaling（TSR）― 時系列モデルの変数重要度を上手く出せるようにした。
  • STRIPE ― 非定常な幅広い時系列の分布予測を実現した。
  • Normalizing Kalman Filters ― ディープニューラルネットを用いた Normalizing Flow 付きのカルマンフィルタで複雑に分布する多変量時系列の分布予測を実現した。
  • 自己符号化器を用いた時系列からのストレンジアトラクタの再構築 ― カオティックな時系列から主要な成分を抜き出す自己符号化器を学習した。
さらに目的別に整理してみます。
  • 予測：
    • 観測が不規則： ODE-LSTMs、NeuralCDE
    • 長期予測： AST
    • スイッチングするシステムの分布予測： Deep Rao-Blackwellised Particle Filters
    • 非定常な幅広い時系列の分布予測： STRIPE
    • 複雑に分布する多変量時系列の分布予測： Normalizing Kalman Filters
  • 因果推論： 自己相関をもつ時系列どうしの因果関係を高感度に発見する手法
  • 予測モデルの解釈： TSR
  • 主成分の特定： 自己符号化器を用いた時系列からのストレンジアトラクタの再構築
大多数の目的が予測であることは昨年と変わりません。昨年に引き続き ODE モデルが複数提案されていますが、昨年と違い再帰 and/or 畳込みのみのモデルがありませんね。そして、予測に分類される論文の半数は昨年なかった（はず）の分布予測です。そして予測に分類されない発表は昨年に比してかなり実践的な興味に寄っているようにみえます。無論、昨年も今年もタイトルを "time series" でフィルタしているので ”sequential” や "recurrent" やその他のワードも確認しなければなりませんが疲れたので今回はこれで…。
「不規則にサンプリングされた時系列からの長期依存性の学習」…興味深いタイトルです。えっと、
  • 常微分方程式による連続的な隠れ状態をもつ RNN は不規則にサンプリングされた時系列のモデリングに有利であると（昨年 ODE-RNNs がありましたよね）。
  • しかし長期依存性の学習が難しいと。それは RNN の勾配消失／爆発に起因することを示すと。
  • なので ODE-LSTMs にすると。
素人考えでは ODE-RNNs が提案されたら一瞬で ODE-LSTMs が提案されそうな気もするんですが、LSTM への移行がそんなに単純ではないのか、LSTM にする意義が非自明なのか、ODE-RNNs 自体がまだ精査される段階であったんでしょうか。いっても翌年に提案されるのは一瞬でしょうか。昨年2つあった ODE 系が今年も2つありますね。
Adversarial Sparse Transformer for Time Series Forecasting
Sifan Wu, Xi Xiao, Qianggang Ding, Peilin Zhao, Ying Wei, Junzhou Huang
「時系列予測のための敵対的で疎なトランスフォーマー」でしょうか。「点予測モデルはデータの確率的な性質を捉えない」というのはわかりますが、確率的なモデルも同様だとはどういう意味でしょうか…。また、多くのモデルは長いステップ先の予測をすると誤差が蓄積されるとありますね。AR モデルのようなのを想定されているはずです。それらのような問題を GAN をベースにした Adversarial Sparse Transformer（AST） で解決すると。ええ…話がいきなりですね…いまの話のいったいどこから GAN が…。…ジェネレータが生成する疎なアテンションは「次のステップを予測するのに重要なのはこれまでのステップのこのステップとこのステップだけだ」みたいな雰囲気なのでしょうか、わかりませんが。ディスクリミネータの役割は、モデルが生成した次ステップと真の次ステップを見分けて「まだ全然見分けられるんだけど？」と叱咤する係なんですかね、妄想ですが。普通は単にモデルの予測値と真の値との誤差を小さくすればよさそうですが、特に自己回帰しながら長期予測する時系列モデルでは中間予測が重要なので特別にこういう工夫を入れてもよさそうです。…これは結構新鮮なアイデアな感じがしますね。
Deep Rao-Blackwellised Particle Filters for Time Series Forecasting
Richard Kurle, Syama Sundar Rangapuram, Emmanuel de Bézenac, Stephan Günnemann, Jan Gasthaus
「時系列予測のためのディープなラオ・ブラックウェル化粒子フィルタ」…粒子フィルタのラオ・ブラックウェル化がわからないので詰んだんですが…検索したら以下が出てきましたね。「粒子フィルタだが一部の状態については解析的にできるので解析的にやる」というのに「ラオ・ブラックウェル化」という名が付いているんでしょうか…？ それでアブストラクトに入ると、スイッチングガウス線形動的システム？ これは響き的に以下の図のようなシステムでいいんでしょうか。現実だとどんなデータがこうなるんでしょうか、コマースサイトのログデータはセール期間中と期間外でモードがスイッチする気はしますが。とりあえずこのようなシステムが生成する時系列をいい感じに学習できるんですね。これは元々ラオ・ブラックウェル化粒子フィルタで扱う問題という認識でいいんでしょうか。スイッチングしたらその先は解析的といった感じがしますし。…あ、ディープ化することで非線形にも対応できるんですかね。何やら提案分布を factorization するようですが…？
未知の交絡因子がある場合でも自己相関をもつ時系列どうしの因果関係を高感度に発見するということでいいんでしょうか。既存手法（FCI and variants？）は感度が低いことを示すと。なので新手法を提案するんですが、このアブストラクトは既存手法の知識がないと読むのが困難ですね…。以下にコードが公開されているようです。この論文はモデリングではなく因果関係の検知ですが（広義にはモデリングなんでしょうか）、その検知にもニューラルネットを使うということでもなさそうですね？ アブストラクトにもそのようにありませんしリポジトリの setup.py でもそれらしいライブラリを入れていませんし。
Neural Controlled Differential Equations for Irregular Time Series
Patrick Kidger, James Morrill, James Foster, Terry Lyons
また常微分方程式ですね。不規則にサンプリングされた時系列のモデリングに有利だが途中経過が初期条件に縛られてしまい後から修正できないと。なので controlled differential equations なるものを導入して解決するんでしょうか？ 提案する NeuralCDE は観測をまたいだ adjoint-based backpropagation を利用できるのも既存手法と異なるようです。これ、途中経過を後から修正できるから一気に複数時点のデータから誤差逆伝播できるという話なんでしょうか？ わかりませんが。NeuralCDE が他の ODE モデルを包含していることも示されているようですね。あ、controlled differential equations で検索したときに本論文の GitHub が出てきたので張っておきます。
Benchmarking Deep Learning Interpretability in Time Series Predictions
Aya Abdelsalam Ismail, Mohamed Gunady, Hector Corrada Bravo, Soheil Feizi
Saliency 法とは何でしょうか…とりあえず入力特徴の重要度を出す手法のようですが、時系列データにはあまり適用されていないようですね。RNN、TCN、transformers などの多様なアーキテクチャ間で Saliency 法に基づく解釈手法のパフォーマンスを比較するようです？ 解釈手法のパフォーマンスとは何でしょうか…？ 精度と感度について以下のようにありますが。
  • 精度 ― 検知された特徴が重要な信号を発していたか。
  • 感度 ― （本当に重要な信号を発していたうち？）検知された特徴の数。
…しかし、これであるネットワーク構造が解釈しやすいとか解釈しにくいとか結論付けられるんでしょうか。Saliency 法がわかっていないからかもしれませんが、手法が悪いかどうかと区別できなさそうですし、どの特徴を重視するかもモデル依存なような…って、基本的には重要な特徴の検知に失敗するんですね？ 2ステップの temporal saliency rescaling（TSR） を用いると改善されるようですが…これはモデルに手を加えなくてもいいんでしょうか。よくわかりません。
非定常な様々な時系列の将来のある時点の分布を予測するのでしょうか。STRIPE というのがモデル名なのですかね。何の略なのか…STRIPE モデルは DPP（決定的点過程）カーネルというのを組み込んで用いるのでしょうか。寡聞にしてアブストラクトからどんなモデルの系譜なのかよくわかりません…非定常な時系列をそのまま分布予測しようとする自体が取り組まれてこなかったとかなんでしょうか。
Normalizing Kalman Filters for Multivariate Time Series Analysis
Emmanuel de Bézenac, Syama Sundar Rangapuram, Konstantinos Benidis, Michael Bohlke-Schneider, Richard Kurle, Lorenzo Stella, Hilaf Hasson, Patrick Gallinari, Tim Januschowski
カルマンフィルタを Normalizing Flow で強化したということですが、Normalizing Flow とは…以下の記事を参考にすると、ガウス分布にしたがう確率変数に非線形変換を繰り返して複雑な分布を表現するようですが、ディープニューラルネットを用いているのはこの非線形変換なんでしょうか。確かに「KF なので大規模な時系列に適応できるし Normalizing Flow を用いるので分布が複雑でもいい」となりそうな気もしますが。
「時系列からのストレンジアトラクタの再構築（ディープニューラルネットを用いた）」でしょうか。ストレンジアトラクタとは？ ローレンツアトラクタを描いたことはありますが、あれはストレンジアトラクタの一種なんですね？ 奇妙であればよいと。まあそれで、問題設定としては、普通はたくさんの変数を観測して主要な成分を突き止めるが、観測が少ないときに隠れた支配的な軸を見つけ出すことをしたいんですかね。入力空間を自己符号化器で理想の空間に埋め込むのでしょうか。損失関数が気になりますね。これが novel とありますが。しかし、これって少ない次元の観測から欠けた次元を復元できるといった話なんでしょうか。普通はできない気がしますが、カオスに動き回るのでその動きさえきちんと捉えればできるということなんでしょうか。わかりませんが。アブストラクト内にあるように、色々なデータについて主要な成分を取り出すのに応用できそうですね。

つづかない

Time Series Analysis： ノート5章 最尤推定（その1）

本読み

以下の本の5章を読みます。キャラクターは架空のものです。私の誤りは私に帰属します。お気付きの点がありましたらコメントでご指摘いただけますと幸いです。

Time Series Analysis

Time Series Analysis

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

5章はそのタイトルも「最尤推定」ですね。これまでは ARMA 過程

Y_t = c + \phi_1 Y_{t-1} + \cdots + \phi_p Y_{t-p} + \varepsilon_t + \theta_1 \varepsilon_{t-1} + \cdots + \theta_q \varepsilon_{t-q}
のパラメータ \theta = ( c, \, \phi_1, \, \cdots, \, \phi_p, \, \theta_1, \, \cdots, \, \theta_q, \, \sigma^2 )^\top が所与の下で、Y_t がどんな性質をもつかとか、あるステップまで観測した下で先のステップの値がどのように予測できるかを扱ったのですよね。対してこの章ではパラメータを未知として、Y_t の観測値 (y_1, \cdots, y_T) からそれを推定しにいくということです。
  • これまで： \theta が所与の下で Y_t の期待値や自己共分散、条件付き期待値を調べた。
  • ここでは： \theta が未知の下で Y_t の観測値 (y_1, \cdots, y_T) から \theta を推定する。
そして \theta の推定のやり方としては、(Y_1, \cdots, Y_T) の確率密度関数 f_{Y_T, \, \cdots, \, Y_1}(y_T, \, \cdots, \, y_1 \, ; \, \theta) の観測値における値を最大にするような \hat{\theta} を推定値とするということです。なんと、そのような考え方があったとは…。

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

それが「最尤推定」だからね。その確率密度関数は、パラメータ \theta の関数と考えると \theta の尤度関数というね。尤度関数の最大点 \hat{\theta} は \theta の最尤推定量だね。最尤推定量は観測値の関数 \hat{\theta}(y_1, \cdots, y_T) だね。

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

それで、「最尤推定をするためにはホワイトノイズ \varepsilon_t の分布の形を特定しなければならない」？ なぜですか？

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

例えば AR(1) 過程 Y_t = c + \phi Y_{t-1} + \varepsilon_t のある1ステップだけの観測値 y_1 が手に入ったとき、尤度関数 L(c, \phi, \sigma^2) = f_{Y_1}(y_1 \, ; c, \phi, \sigma^2) ってどうなる？ 過程は定常とするよ。

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

y_1 だけ観測したとき、ですか。パラメータがある (c, \phi, \sigma^2) だったときにその y_1 がどれくらい観測されやすいのかを考えなければなりませんね。うーん…そもそも定常 AR(1) 過程には期待値があったはずです。パラメータが (c, \phi, \sigma^2) であるような定常 AR(1) 過程の期待値 E(Y_t) を計算して、y_1 がその期待値に近ければその (c, \phi, \sigma^2) の尤度は大きいでしょうし、期待値から遠ければその (c, \phi, \sigma^2) の尤度は小さいのではないかと思います。定常 AR(1) 過程の期待値を復習します。AR(1) 過程の式の両辺の期待値をとります。

\begin{split} E(Y_t) &= E(c + \phi Y_{t-1} + \varepsilon_t) \\ \Longleftrightarrow \; E(Y_t) &= c + \phi E(Y_{t-1}) \end{split}
いま過程が定常であることを仮定していますから、期待値は時刻によらず E(Y_t) = E(Y_{t-1}) のはずです。これより、E(Y_t) = c/ (1 - \phi) を得ます。…AR(1) 過程の期待値はわかりましたが、この期待値に照らし合わせて y_1 がどれくらい観測されやすいかは、AR(1) 過程の分散にもよりますね。分散が大きければ期待値から離れた値も観測されやすいですから。分散も先ほどと同様に求められます。
\begin{split} V(Y_t) &= V(c + \phi Y_{t-1} + \varepsilon_t) \\ \Longleftrightarrow \; V(Y_t) &= V(\phi Y_{t-1}) + V(\varepsilon_t)+ 2{\rm Cov}(\phi Y_{t-1}, \varepsilon_t) \\ \Longleftrightarrow \; V(Y_t) &= \phi^2 V(Y_{t-1}) + \sigma^2 \end{split}
分散もまた時刻によらず V(Y_t) = V(Y_{t-1}) ですから、V(Y_t) = \sigma^2/ (1 - \phi^2) を得ます。過程の期待値と分散がわかったので、これに照らし合わせて y_1 がどれくらい観測しやすいかを求めればいいです。

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

期待値と分散だけで分布の形って決まる？

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

！ き、決まりませんでした…。で、では、歪度や尖度も求めれば…。

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

きりがないな…というか、定常（弱定常）って歪度や尖度まで一定とはいっていないからね。そもそもそのやり方で求まらないよ。

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

なんと…詰みました…。

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

だから扱いやすいホワイトノイズを選べばいいんだよ。正規ホワイトノイズをね。それなら任意の t で Y_t の分布は正規分布だよ。正規分布は再生性をもつから、正規分布にしたがうノイズを足していっても正規分布にしたがうからね。なら、期待値と分散だけで分布の形が決まるね。

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

確かに正規分布であれば期待値 \mu と分散 \sigma^2 だけで確率密度関数が特定できます。

\displaystyle f(x) = \frac{1}{\sqrt{2 \pi \sigma^2}} \exp \left( - \frac{(x-\mu)^2}{2\sigma^2}\right)
でしたら Y_t の分布は先ほど求めた期待値と分散を上式に当てはめて、以下になるはずです。
\displaystyle f_{Y_1}(y_1 \, ; c, \phi, \sigma^2) = \frac{1}{\sqrt{2 \pi \sigma^2/ (1 - \phi^2)}} \exp \left( - \frac{ \bigl( y_1 - c/ (1 - \phi) \bigr)^2}{2\sigma^2/ (1 - \phi^2)}\right)
これが尤度関数 L(c, \phi, \sigma^2) に他なりません。観測値 y_1 が手に入ったとき、 L(c, \phi, \sigma^2) を最も大きくする (c, \phi, \sigma^2) が y_1 を生成した AR(1) 過程のパラメータだったのだろうと推定すればいいんですね。正規分布のお陰で詰みを回避できました…しかし、ホワイトノイズが正規分布でなかったらどうなるんです？ 確率密度関数が特定できなくてやはり詰みませんか？

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

この章では扱われていないと思う。5.2～5.6 節の全てのタイトルに明示的に Gaussian って付いてるし。正規分布じゃなかったら尤度関数をかき下すことはできないからね。じゃあ正規分布じゃなかったらどうするのかは、maximum likelihood estimates of non-Gaussian ARMA processes とかで調べたら色々文献が出てくるんじゃないかな。

ところで、いまはある1ステップだけの観測値 y_1 が手に入ったときの尤度関数を計算したけど、普通は連続する観測値 (y_1, \cdots, y_T) からパラメータを推定するよね。じゃあ、連続する2ステップ (y_1, y_2) を観測したときの尤度関数はどうなる？

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

2点観測したのであれば、その確率密度関数は、1点観測したときの確率密度関数を2乗すればいいのでは？

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

Y_1 と Y_2 は独立じゃないからね…Y_2 が Y_1 に依存しているからこそ AR（自己回帰）過程なんだけど…。

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

！ す、すみません…AR(1) 過程は Y_t = c + \phi Y_{t-1} + \varepsilon_t ですから、 Y_2 = c + \phi Y_1 + \varepsilon_2 ですね。Y_2 はこれにしたがって生成されなければなりません。これは、Y_1 = y_1 が観測されている下では、c + \phi y_1 にホワイトノイズ \varepsilon_2 を足したものですから、

\displaystyle f_{Y_2 | Y_1}(y_2 \, | y_1 \, ; c, \phi, \sigma^2) = \frac{1}{\sqrt{2 \pi \sigma^2}} \exp \left( - \frac{ \bigl( y_2 - c - \phi y_1 \bigr)^2}{2\sigma^2}\right)
です。でしたら (y_1, y_2) の同時分布は、 f_{Y_2, Y_1}(y_2, y_1 \, ; c, \phi, \sigma^2) =f_{Y_2 | Y_1}(y_2 \, | y_1 \, ; c, \phi, \sigma^2) f_{Y_1}(y_1 \, ; c, \phi, \sigma^2) より、以下ですね。これが (y_1, y_2) を観測したときの尤度関数です。
\displaystyle f_{Y_2, Y_1}(y_2, y_1 \, ; c, \phi, \sigma^2) = \frac{\sqrt{1 - \phi^2}}{\left( \sqrt{2 \pi \sigma^2} \right)^2} \exp \left( - \frac{ \bigl( y_2 - c - \phi y_1 \bigr)^2}{2\sigma^2} - \frac{ \bigl( y_1 - c/ (1 - \phi) \bigr)^2}{2\sigma^2/ (1 - \phi^2)}\right)

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

うん。さらに連続する3ステップ (y_1, y_2, y_3) を観測したときはこうなるね。

\begin{split} f_{Y_3, Y_2, Y_1}(y_3, y_2, y_1 \, ; c, \phi, \sigma^2) &= f_{Y_3 | Y_2, Y_1}(y_3 \, | y_2, y_1 \, ; c, \phi, \sigma^2) f_{Y_2, Y_1}(y_2, y_1 \, ; c, \phi, \sigma^2) \\ &= \; f_{Y_3 | Y_2, Y_1}(y_3 \, | y_2, y_1 \, ; c, \phi, \sigma^2) \\ & \quad \cdot f_{Y_2 | Y_1}(y_2 \, | y_1 \, ; c, \phi, \sigma^2) \\ & \quad \cdot f_{Y_1}(y_1 \, ; c, \phi, \sigma^2) \\ &= \; f_{Y_3 | Y_2}(y_3 \, | y_2 \, ; c, \phi, \sigma^2) \\ & \quad \cdot f_{Y_2 | Y_1}(y_2 \, | y_1 \, ; c, \phi, \sigma^2) \\ & \quad \cdot f_{Y_1}(y_1 \, ; c, \phi, \sigma^2) \end{split}
最後の式への式変形は、AR(1) 過程が1ステップ前までにしか依存しないことによったよ。具体的に確率密度関数をかき下すとこうだね。
\displaystyle f_{Y_3, Y_2, Y_1}(y_3, y_2, y_1 \, ; c, \phi, \sigma^2) = \frac{\sqrt{1 - \phi^2}}{\left( \sqrt{2 \pi \sigma^2} \right)^3} \exp \left( - \sum_{t=2}^3 \frac{ \bigl( y_t - c - \phi y_{t-1} \bigr)^2}{2\sigma^2} - \frac{ \bigl( y_1 - c/ (1 - \phi) \bigr)^2}{2\sigma^2/ (1 - \phi^2)}\right)
これはこのまま一般化できて、(y_1, \cdots, y_T) を観測したときの確率密度関数は以下になるね。
\displaystyle f_{Y_T, \cdots, Y_1}(y_T, \cdots, y_1 \, ; c, \phi, \sigma^2) = \frac{\sqrt{1 - \phi^2}}{\left( \sqrt{2 \pi \sigma^2} \right)^T} \exp \left( - \sum_{t=2}^T \frac{ \bigl( y_t - c - \phi y_{t-1} \bigr)^2}{2\sigma^2} - \frac{ \bigl( y_1 - c/ (1 - \phi) \bigr)^2}{2\sigma^2/ (1 - \phi^2)}\right)
これが尤度関数 L(c, \phi, \sigma^2) であり、これを最大化する (c, \phi, \sigma^2) を探すことになるけど、通常はこれの対数をとった対数尤度関数 \mathcal{L}(c, \phi, \sigma^2) を最大化するから、テキストではこちらで表記しているね。それで、上の式は（対数を取れば）予測誤差 y_t - c - \phi y_{t-1} の2乗を足していっているようにみえる。これを尤度関数の予測誤差分解（prediction-error decomposition）というみたいだね。

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

おお、AR(1) 過程の尤度関数はこれで求まりましたね。…ん？ 119ページからは AR(1) 過程の尤度関数を求める別の方法が紹介されていますね？ y \equiv (y_1, \cdots, y_T)^{\top} と定義するようです。このようなベクトルを定義してどうするんでしょうか…。

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

y が多変量正規分布にしたがうと考えるんだね。期待値ベクトル \mu と分散共分散行列 \Omega は以下になるね。

\displaystyle \mu = \left(\begin{array}{c} E(Y_1) \\ E(Y_2) \\ \vdots \\ E(Y_T) \end{array}\right)= \left(\begin{array}{c} c/ (1 - \phi) \\ c/ (1 - \phi) \\ \vdots \\ c/ (1 - \phi) \end{array}\right)
\displaystyle \Omega = \left( \begin{array}{c} V(Y_1) & {\rm Cov}(Y_1, Y_2) & \cdots & {\rm Cov}(Y_1, Y_T) \\ {\rm Cov}(Y_1, Y_2) & V(Y_2) & \cdots & {\rm Cov}(Y_2, Y_T) \\ \vdots & \vdots & \ldots & \vdots \\ {\rm Cov}(Y_1, Y_T) & {\rm Cov}(Y_2, Y_T) & \cdots & V(Y_T) \end{array} \right)= \frac{\sigma^2}{1-\phi^2}\left(\begin{array}{c} 1 & \phi & \cdots & \phi^{T-1} \\ \phi & 1 & \cdots & \phi^{T-2} \\ \vdots & \vdots & \ldots & \vdots \\ \phi^{T-1} & \phi^{T-2} & \cdots & 1 \end{array}\right)
ちなみに AR(1) 過程の自己共分散は1時点ずれるごとに分散に \phi が1つかかってくるからね。
{\rm Cov}(Y_t, Y_{t-1}) = {\rm Cov}(c + \phi Y_{t-1} + \varepsilon_t, Y_{t-1}) = \phi {\rm Cov}(Y_{t-1}, Y_{t-1}) = \phi V(Y_{t-1})
そうすると直接 y の確率密度関数が出せる。
\displaystyle f_Y(y \, ; c, \phi, \sigma^2) = \frac{1}{ \sqrt{ (2 \pi)^T | \Omega | } } \exp \left( - \frac{1}{2} (y - \mu)^\top \Omega^{-1} (y - \mu)\right)
これはさっき求めた確率密度関数と一致する。計算には \Omega = \sigma^2 V とおいたときに L^{\top} L = V^{-1} を満たす L を利用しているね。

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

いやいやいや、L が降って湧いていますよね？ L^{\top} L = V^{-1} が成り立つことは検算して確かめよというようにありますが、この L はどうやって考えついたものなんです？

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

うーん、V^{-1} から L を求めるのは LU 分解ともちょっと違うかな？ でも、L にはちゃんと意味があるよ。だって、\Omega^{-1} = \sigma^{-2} L\top L を確率密度関数に代入したら、

\displaystyle f_Y(y \, ; c, \phi, \sigma^2) = \frac{1}{ \sqrt{ (2 \pi)^T | \Omega | } } \exp \left( - \frac{1}{2 \sigma^2} (y - \mu)^\top L^\top L (y - \mu)\right)
となるから、\hat{y} = L(y-\mu) とおくと \hat{y} は各成分が独立で各成分の分散が \sigma^2 になるわけだよね。y の各成分に対して上から順に、「前の成分と独立にした上で、平均を 0 にして、分散を \sigma^2 にする」という操作をやっていけば \hat{y} になるよ。 L はその変換を行列表示したものに他ならない。AR(1) 過程は1ステップ前に依存するだけだからそんな変換は求まるし、そんな変換を求めておけば、分散共分散行列の逆行列を計算する必要もないよね。

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

なるほど？ それで、122ページは実際に対数尤度関数の最大点を求める方法ですね。あれ？ 対数尤度関数そのものではなくて、y_1 が所与の場合の対数尤度関数を最大化しているんですね？

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

その場合は最小2乗法をかき下せるからね。長い系列を観測している場合、最初の観測値の影響は無視できるし。

つづいたらつづく