沖本計量時系列本の7章GARCHモデルの勉強会用ノートなんですがブラウザから見ると数式の一部及び式番号が見えないようです？
github.com
異常検知と変化検知： 10章メモ（疎構造学習による変化検知）
以下の赤い本の10章を読みます。キャラクターは適当です。誤りがありましたらご指摘いただけますと幸いです。
 | 異常検知と変化検知 (機械学習プロフェッショナルシリーズ) 井手 剛 杉山 将 講談社 2015-08-08 売り上げランキング : 56236 Amazonで詳しく見る by G-Tools |
他の章のメモ： 9章（部分空間法による変化検知）
特異スペクトル変換法は式 (9.8) の過去側の分布／現在側の分布どちらの にもフォンミーゼス・フィッシャー分布（vMF分布）を仮定したものだけど、vMF分布は長さ1のベクトルの分布であるのに対して、部分時系列はそれ自体は長さ1のベクトルじゃない。部分時系列を並べた行列を特異値分解して左特異ベクトルをとる、という操作をすることで長さ1のベクトルにするんだったね。時系列の変化度とは抽象的には「少し前の状況（分布）と今の状況（分布）のKL情報量」と定義されたけど、特異スペクトル変換法は「少し前の状況」及び「今の状況」を「1本（ないしは数本）の単位ベクトル」で表現するという点で、仮定する分布というよりこの変数変換自体が独特であるように感じるね。でもこの「元データ → 元データを特異値分解して出てきた特異値の大きい数本の左特異ベクトル」という変換がノイズ除去に相当するみたいだね。
なんでそんな吹奏楽部がぎすぎすしそうな例を考えたの…。
知らないよ…。
ただ変数間の依存関係に着目するには「直接相関」と「間接相関」の区別が重要ということで…130～131ページにある「都市ごとの教会の数と殺人件数の相関（どちらも都市の人口に比例しているだけ）」のような例って岩波DS3でも見ましたね。
小学生の握力と算数ドリルの点数の疑似相関だったっけ。学年が交絡因子で。
ああそういうのでしたね。でも、算数ってノートにねばり強く書き続ける力だと思うんですよね。つまり、握力が高まるほど長時間ノートに力強く数式を書き続けられるようになって、算数力も向上するのではないでしょうか。
話ややこしくしなくていいよ…。
132ページからもしばらく時系列からは離れて（？）、確率ベクトルの対マルコフグラフをデータからどう描けばいいのか、という話が続きますね。それでもうこの本では多変量正規分布を仮定したガウス型グラフィカルモデルのみ取り扱うようですね。132ページの「多変量正規分布というめがねをかけて」という表現がわかるようでわかりませんが置いておきますね。…あれ、データに多変量正規分布を仮定したときの共分散行列の推定値ってこの (2.4) 式でいいんでしたっけ。最近別件でアンサンブルカルマンフィルタを実装したんですけど、そちらではフィルタを実行するときアンサンブル粒子から以下の下側の式のように分母が N-1 の共分散行列を計算して、これを用いてカルマンゲインを計算したと記憶しているんですが…。
前者が最尤推定量で後者が不偏推定量だね。これらの違いは、雑に言うと以下のような感じかな…。
- 手元にデータがある。色々な多変量正規分布の中から、手元のデータを生成する確率が最も高いものを選ぶ。これは対数尤度を最大化すればよい（計算略）。
- 手元にデータがある。このデータがある多変量正規分布から生成されたということを踏まえる＝つまり、手元のデータはたまたま手元にあるように生成されただけで、別の世界線では同じ多変量正規分布から生成された別のデータになっていたかもしれない。これを踏まえた上で、色々な多変量正規分布の中から、手元のデータを生成する確率が最も高いものを選ぶ。手元のデータを母集団と見做すならば共分散行列は式 (2.4) の標本共分散行列でいいのだが、「標本平均がたまたま観測された量である」という仮定の下では、この「標本平均（データに対して最も都合がよい場所に仮定した平均）の周りのばらつき」でばらつきを見積もるのはやや過小評価になる（証明略）。
- 異常検知におけるガウス型グラフィカルモデルの推定では、手元のデータを生成した最も尤もらしい共分散行列を知りたいので最尤推定量を求める。それを通して手元のデータの変数間の相関構造を知りたいだけなので、サンプリングが違ったらどうだっただろう、という考慮はそぐわない。
- アンサンブルカルマンフィルタでは状態分布を代替的にアンサンブル粒子で表現していて、フィルタのステップではアンサンブル粒子を共分散行列に翻訳したい。このとき、現在のアンサンブル粒子はたまたま今回の計算のために生成しただけなので、現在のアンサンブル粒子を生成した最も尤もらしい共分散行列を知りたいのではない。色々なアンサンブル粒子の選び方があることを踏まえて、不偏推定量を求める。
お疲れさまです。ところで、133ページを読むと、結局 と
の間に直接相関がないことと、精度行列の
成分がゼロであることが同値なんですね？ ということは、「教会の数」「殺人件数」「人口」のデータの場合、「教会の数」と「殺人件数」の組合せに対応する精度行列の要素がゼロになるということですが、本当にちゃんとそうなるんでしょうか…。疑似相関とはいえ相関はあるんですよね？ あとこれって一意に定まらなさそうな気もするんですよね。だって「教会の数」「殺人件数」「人口」のどれがニワトリでタマゴかなんてデータは語らないですよね？
今回のケースでは、132ページのグラフみたいに、「人口」軸で輪切りにしたときに「教会の数」と「殺人件数」は無相関だ、というのはわかるよね。「教会の数」軸での輪切りや「殺人件数」軸での輪切りではそうならないんじゃないかな。雑な例だけど実験してみるね。まず、0番目（ Python だから便宜上 0 から始めるね）の変数から1 番目の変数と2番目の変数それぞれにグラフの辺があって、1番目の変数と2番目の変数の間にはグラフの辺がない例。
# -*- coding: utf-8 -*- import numpy as np # 標本共分散行列 def calc_cov_mat(z_en): z_hat = z_en - np.average(z_en, axis=0) v = np.outer(np.zeros(z_en.shape[1]), np.zeros(z_en.shape[1])) for i_en in range(z_en.shape[0]): v = v + np.outer(z_hat[i_en], z_hat[i_en]) v = v / z_en.shape[0] return v n_en = 100 x0_en = np.random.randn(n_en, 1) x1_en = x0_en + 0.01 * np.random.randn(n_en, 1) x2_en = -1.0 * x0_en + 0.01 * np.random.randn(n_en, 1) x_en = np.concatenate([x0_en, x1_en, x2_en], axis=1) Sigma = calc_cov_mat(x_en) Lambda = np.linalg.inv(Sigma) print(Sigma / Sigma[0, 0]) print(Lambda / Lambda[0, 0])
[[ 1. 0.99916607 -1.00078985] [ 0.99916607 0.99843717 -0.9999501 ] [-1.00078985 -0.9999501 1.00167383]] [[ 1. -0.4735883 0.52634417] [-0.4735883 0.44918519 -0.02475815] [ 0.52634417 -0.02475815 0.50121079]]
精度行列 Lambda の (1, 2) 要素＝(2, 1) 要素は小さいね。試しに変数の役割を入れ替えて、x2―x0、x2―x1 にグラフの辺がある場合もやってみるね。
x2_en = np.random.randn(n_en, 1) x0_en = x2_en + 0.01 * np.random.randn(n_en, 1) x1_en = -1.0 * x2_en + 0.01 * np.random.randn(n_en, 1) x_en = np.concatenate([x0_en, x1_en, x2_en], axis=1) Sigma = calc_cov_mat(x_en) Lambda = np.linalg.inv(Sigma) print(Sigma / Sigma[0, 0]) print(Lambda / Lambda[0, 0])
[[ 1. -0.99804501 0.99911337] [-0.99804501 0.99624298 -0.99724056] [ 0.99911337 -0.99724056 0.99831405]] [[ 1. 0.08005913 -0.92082764] [ 0.08005913 1.15942445 1.07805449] [-0.92082764 1.07805449 1.99854651]]
しかし137ページに の場合は「偏微分＝0」が解析的に解けないとありますね。
だったら解けるんでしょうか…そもそも今回は精度行列を疎にしたいという目的より必ず
なので、
の場合に興味はありませんが…。それで、解析的に解けないのでどうするかというと…
次元目以外固定してしまう？
次元目を一番後ろに寄せて、精度行列、精度行列の逆行列、標本共分散行列を137ページのようにブロック表記して、そうすると各行列が4つの部分に分解されて、うち2つは転置の関係なので、元の「偏微分＝0」は (10.15) (10.16) (10.17) の3式になるわけですね…それで138～139ページの手続きにしたがえば元の精度行列の
次元目の行・列以外を固定したもとで目的関数を極大値にできますが…それを全ての次元について、ぐるぐると繰り返せば目的の精度行列が求まると…このように、元データを多変量正規分布のめがねでみてかつ疎な精度行列を得ようとする手法を「グラフィカルラッソ」というんですね。うーん、手順はわかるんですが、これ最適解に収束するんですか？ ある次元以外固定してしまおう、というのがなんか心配なんですが…。
その辺の議論はこの本にはなさそうだね。参考文献を参照かな。
…まあいいです。何らかの参照データから疎な精度行列が求まり、変数間の依存関係がわかったとしましょう。これを利用する異常検知手法が 10.5 節に2例紹介されています。1例目は、ある多変量の観測値が得られたときに、その 次元目の要素の異常度は、
次元目以外の要素を given とした条件付き確率密度の対数のマイナス1倍とする（外れ値解析）。
つまり、時系列の 次元目だけに着目してよく観測される量だったとしても、
次元目以外がこの状況のときには観測されない量だ、という場合は異常と判定されるようなイメージだね。
2つ目の例は1点の観測値だけでなくまとまった観測値に対する異常判定ですね。このときは、参照データと観測データのそれぞれの分布を求めて、両者間のKL情報量を 次元目以外について積分したものが異常度ということです。どちらの例も特にグラフィカルラッソや、何なら多変量正規分布でなくてもよさそうですね。ただ相関構造をシンプルにしておかなければ有効な解析ができないということなのですかね。多変量正規分布の場合の (10.28) → (10.29)、(10.30) → (10.31) の式変形は後で追っておきましょう。…143ページに「
と
の立場を入れ替えた定義も可能です」とありますが、確かにKL情報量は非対称ですが、実際過去側の分布と現在側の分布のどちらをどちらにすべきなんでしょう？
9章を読んだときにも似た話をしたね。KL情報量は対数密度比の期待値だから、過去側の分布で期待値をとるのか現在側の分布で期待値をとるのかって話だと思うけど…過去側の分布が基準と考えれば (10.30) 式の向きが自然なんだろうけど。異常時に分布がどう変化するかにもよるだろうね。
雑記
機械学習で Leakage ってよく聞きますけど、結局 Leakage の語義というか具体的に言い切ると何なのかよくわからないんですよね。漏れてはいけない情報が漏れていること、のようではあるんですが、何が漏れてはいけないんでしょう？ 我々は総力を挙げてテストデータの正答を予測するモデルの構築に取り組まなければならないのではないんでしょうか？
ここ https://www.kaggle.com/docs/competitions に What is Leakage? とあるから読んでみたらどうかな。「予期せぬ混入した情報」であって以下のようなタイプがあるって書かれているね。
- テストデータに対する正答が漏れてしまっている。
- テストデータであるはずのデータが訓練データに混ざってしまっている。
- 上の2つは言うまでもないね。自分でデータ分析をするときはバグでこれらの漏洩をしないように注意だね。コンペティションではじゅうぶん気を付けられていると思うけど、回答者が調べれば正答が突き止めることができてしまうような出題はできないね。
- モデルが利用できないはずの未来の情報が混ざってしまっている。
- 以下のブログで紹介されている、「売上を予測するコンペティションで、店舗の住所を特定して（予測時点では知りえないはずの）未来の天候情報を紐づける」というのがこれだよね（もっともこれはデータに混ざっていたわけじゃなくて、第三者データから持ち込んだわけだけど）。https://adtech.cyberagent.io/techblog/archives/259
- 意図的に加えたノイズや匿名化が除去されてしまう。
- これはよくわからないけど、研究の目的上とか実用上の制約で匿名化データを使用しなければならないときに、個人名をどこかから探してきてくっつけたら駄目みたいな話だね。モデルが自分の力でノイズを除去したり個人の性質を同定したりするならよさそうな気はするけどね。
- 何らかの事情で利用禁止することにしたために除去した情報が残ってしまっている。
- モデルが運用上入手できない情報が利用されてしまう。
- これらは「未来の情報」「個人情報」以外で目的上「この情報からの学習は意図しないという情報」って感じがするね。コンペティションの出題者も、意図しない情報がどこかから発掘されないように出題しないといけないから大変だよね。
- データがもつ情報が利用目的から外れて歪められてしまう。
- これは以下の記事で挙げられている例のことだと思うんだよね。テニスのルールは詳しくないんだけど、両選手のミスショットの回数とかから勝敗を予測するのが学習の意図なんだけど、プレイヤー1の名前とプレイヤー2の名前を説明変数に入れてしまうともうモデルはこの対戦カードの対戦結果はこれって覚えちゃうよね。でもそんなことをモデルにしてほしかったわけじゃない。
https://tjo.hatenablog.com/entry/2016/01/27/235620 - これらの情報が第三者データから持ち込まれてしまう。
なるほど。その未来の情報という例で思ったんですが、こういうのはどうです？ 未来の天候情報そのものを参照するとズルなので、未来の天候情報を予測するモデルを謎の目的関数で学習して、そのモデルの予測を利用して未来の売上を予測するというのは。
いや私運営じゃないから知らないんだけど…確かにモデルがあからさまにデータを内蔵していなくて自分で天候を予測しているなら問題は少なそうだけど、それで達成できるような精度は直接未来の売上を予測しても達成できるんじゃないの…？ もちろんモデルによるけど…。だいたい中間目標を置くのは多くの場面で得策じゃないよね。天候情報の予測がよほど楽勝だったら別だけど、じゃあやっぱり大してプラスにならないんじゃないの。
ぐぬぬ…。
ちょっとうっかりしすぎじゃないですかね。因果関係ってご存知ですかって感じなんですけど。
いやでも、こんなことも実務で起こらないとはいえないから気を付けないとねってことなんじゃないかな。これは極端な例で、実際の Leakage はもっとわかりづらいって書いてあるよ。あと続く例は、実際に初期の Kaggle のコンペティションで、ソーシャルネットワークの「つながり」の予測問題があったみたいなんだけど、サンプルスクリプトのミスで、ある条件を満たすユーザどうしは無条件で「つながっている」というラベルがついてしまっていたみたいで、これを突いたチームが2位になったみたいだね。さらに1位のチームは、スクレイピングによってデータの匿名性を剥がしちゃったみたいだけど。
…それって出題ミスですよね。2位のチームは、出題者の意図とは違うモデルを構築してしまったとはいえ、データのパターンを見出したのだから Leakage ですらないような。1位のチームは…大会のレギュレーションによるでしょうね。…最後の Leakage in Competitions というパラグラフにはなんて書いてあるんです？
雑記： 沖本本読み勉強会用のメモ
以下の本を読む勉強会用のメモです。もし何か誤りがありましたらご指摘いただけますと幸いです。
 | 経済・ファイナンスデータの計量時系列分析 (統計ライブラリー) 沖本 竜義 朝倉書店 2010-02-01 売り上げランキング : 19813 Amazonで詳しく見る by G-Tools |
今度勉強会で上の本の3章を発表することになりました。
へー。1章と2章は読んだの？
読んでないです。
読もうよ！？
仕方ないですね…1章で思い出しておきたい事項は以下でしょうか。
- 時系列分析では時系列データの背後にデータ生成過程（または単に過程）を仮定する。
- 期待値が時刻tに依存せず一定で、自己共分散も時刻tに依存せずラグkにのみ依存するような過程を弱定常であるという。
- 特に、ある長さkの部分時系列の同時分布が時刻tに依存せず同一の分布にしたがうような過程を強定常であるという。
- 正規過程が弱定常であれば、強定常である。
- 期待値が時刻tに依存せずゼロ、0次の自己共分散（つまりただの分散）が時刻tに依存せず一定、1次以上の自己共分散が時刻tに依存せずゼロであるような過程をホワイトノイズという。
- データの自己相関の有無を判断するには、以下のような検定をする。
とは限らないかな。例を示そうか。
n <- 20000 breaks <- seq(-10,10, 0.5) layout(matrix(c(1,1,2,3,4,4,5,6,7,7,8,9), 4, 3, byrow=FALSE)) MyPlot <- function(x, y) { plot(c(-10,10), c(-10,10), xlab="x", ylab="y", type="n") abline(h=0, v=0, lty=2) points(x, y) hist(x, breaks=breaks) hist(y, breaks=breaks) } ### case 1 ### x <- rnorm(n=n, mean=0, sd=2.0) y <- rnorm(n=n, mean=0, sd=2.0) MyPlot(x, y) ### case 2 ### library(mvtnorm) mu <- c(0, 0) Sigma <- matrix(c(4, 3.6, 3.6, 4), 2, 2) xy <- rmvnorm(n, mu, Sigma) x <- xy[,1] y <- xy[,2] MyPlot(x, y) ### case 2.1 ### x0 <- x[x < -3 | 3 < x] y0 <- y[x < -3 | 3 < x] x1 <- x[!(x < -3 | 3 < x)] y1 <- y[!(x < -3 | 3 < x)] x1a <- x1[x1 < 0] x1b <- x1[!(x1 < 0)] MyPlot(c(x0, x1a, x1b), c(y0, y1))
上の3つの例はどれもXの周辺分布が N(0, 4)、Yの周辺分布も N(0, 4) だけど、一番右のプロットの形状は2変量正規分布にみえる？
一番右のプロットの形状ですか？ 気持ち悪いです。
ごめん…。
次は2章ですね。ここでMA過程とAR過程が降って湧いたように基本的なモデルとして登場しますが、さもそれらが自己相関をもつ確率過程を表現しようとした自然な帰結なのだというような雰囲気ですが、自己相関をもつ確率過程を表現するのに「y_{t} と y_{t-1} が共通の確率変数に依存する」と「y_{t} が y_{t-1} を参照する」以外のモデルは本当にないんでしょうか？ y_{t-1} が y_{t} を参照するとか。
未来の値が現在の値を決定するモデルはやめようよ？
データ生成過程さんにしてみれば、エントロピーの増大方向など知ったことではなくないですか？
仮にデータ生成過程さんが気にしなくても、そのモデルを利用する我々が気にするんだよ！
なるほど。…それで本に戻ると、MA過程の方は次数・パラメータにかかわらず定常であると。これは直感的に理解できますね。有限の過去までのノイズの線形結合で爆発しようがありませんから。それでAR過程の方は (2.15) 式の全ての解の絶対値が1より大きいときに定常となると。なぜそうなるかのやわらかい解釈は 2018-03-07 の記事を参照してください。また、MA過程の問題点として、ある過程を記述する係数とホワイトノイズの分散のセットの候補が複数存在してしまうので、AR(∞)過程に書き直せる（＝反転可能である）ものを選択しようという基準があるということですね。
「同一の期待値と自己相関構造をもつMA(q)過程が2^q通り存在し、その内1つのみが反転可能である」というのをきちんと追えてないけど、「MA過程はAR(∞)過程に書き直せる」「だから少なくともAR(∞)過程に書き直せるようなホワイトノイズの取り方が1つ存在する」というのは理解できるかな。
2章の後半はモデルの推定と解釈の話ですが…43ページに、ARMAモデルやGARCHモデル、マルコフ転換モデルなどは最小2乗推定できないってありますけど、なんででしょう？ 何かモデルを仮定すれば残差平方和を計算することは常にできますよね？
まずARMAモデル（というよりMAモデル）については、ある時刻での残差が次の時刻の説明変数になっちゃうよね。だから「理想的には残差がゼロだろう」っていうスタンスのOLSで推定したら駄目だよね。GARCHとマルコフ転換モデルはもっと後の方で出てくるからすぐ追えないけど、GARCHはぶれの大きさをモデリングするし、マルコフ転換モデルは観測できない真の状態がぶれるモデルだし、やっぱり「理想的にはぶれないだろう」っていうOLSの考え方はそぐわないようにみえる。
あー、確率的な変動が全くないとすると意味をなさないようなモデルだと駄目っぽいですね。それでやっと3章に辿り着きましたが、この章のタイトルは「予測」です。しかし予測とは何なのかと。
大切な問いだね。
AR(1)過程を例に挙げてこれを問うていますが…「AR(1)過程と現時点の観測値を所与とする。次の時点の予測値とは何か」と。「予測A」はAR(1)過程の期待値そのものを予測値としていて、入手している現時点の観測値を完全に無視していくスタイルですね。「予測B」はノイズの最頻値、「予測C」はノイズの平均値が次の時刻に実現するだろうというような予測です。別にどの予測もそういう予測ではあるわけですが、ここでは（一期前までに観測される値が所与の下での）予測誤差の2乗の期待値が小さいほどよい予測とし、予測誤差の2乗の期待値を最小にする予測を最適予測とよぶということです。それで、先の例では「予測C」が最適予測だと…ん、A～Cの3つの予測の中では「予測C」が最もよいというのはわかるんですが、「予測C」よりもよい予測はないことは示されてないですよね？
「予測C」が「一期先の値の期待値を予測値とする」として得られたことは、それが予測誤差の2乗の期待値を最小に予測であることの説明になってないってことかな。 を
について微分してみたらどうだろう？
思い出しました。2次損失は期待値で最小に、線形損失は中央値で最小に、0-1損失は最頻値で最小になるとかDLMの本でやりましたね。…その後は、AR過程の何期も先のMSEを陽に書くのは難しい、区間予測はシミュレーションでやろう、というような話ですね。…なんか3章はそんな感じですね。
