異常検知と変化検知: 10章メモ(疎構造学習による変化検知)

以下の赤い本の10章を読みます。キャラクターは適当です。誤りがありましたらご指摘いただけますと幸いです。

異常検知と変化検知 (機械学習プロフェッショナルシリーズ)異常検知と変化検知 (機械学習プロフェッショナルシリーズ)
井手 剛 杉山 将

講談社 2015-08-08
売り上げランキング : 56236

Amazonで詳しく見る
by G-Tools
f:id:cookie-box:20180513082851p:plain:w60

先の9章では「『時系列が変化した』とはどのように判定するのか」という決してトリビアルではない問題について、基本的な発想の例として「正常というより異常らしい度合い(対数尤度)のランプ関数の累積和が閾値に達したら変化したと判定する(累積和法)」という解法と、「部分時系列の、予め定めた部分時系列群のうち最近傍標本までの対数距離が大きい場合に異常部位と判定する(近傍法)」という解法を与えました。また、「少し前の状況と今の状況が違うか否か」という形で変化検知問題を抽象的に定式化し、その最も重要な実装例の1つである特異スペクトル変換法を扱いましたね。

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

特異スペクトル変換法は式 (9.8) の過去側の分布/現在側の分布どちらの p にもフォンミーゼス・フィッシャー分布(vMF分布)を仮定したものだけど、vMF分布は長さ1のベクトルの分布であるのに対して、部分時系列はそれ自体は長さ1のベクトルじゃない。部分時系列を並べた行列を特異値分解して左特異ベクトルをとる、という操作をすることで長さ1のベクトルにするんだったね。時系列の変化度とは抽象的には「少し前の状況(分布)と今の状況(分布)のKL情報量」と定義されたけど、特異スペクトル変換法は「少し前の状況」及び「今の状況」を「1本(ないしは数本)の単位ベクトル」で表現するという点で、仮定する分布というよりこの変数変換自体が独特であるように感じるね。でもこの「元データ → 元データを特異値分解して出てきた特異値の大きい数本の左特異ベクトル」という変換がノイズ除去に相当するみたいだね。

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

続く10章は、多変量時系列の異常検知において各変数の寄与度を知るにはどうすればいいのか、という話のようですね…察するに、ブラスバンドのハーモニーが乱れたときに、トランペットが悪いのかクラリネットが悪いのかチューバが悪いのか、各パートの異常への寄与度を知りたいといった感じですね。

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

なんでそんな吹奏楽部がぎすぎすしそうな例を考えたの…。

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

出力値が常に変動するような時系列から異常を検知するには、変数間の依存関係に着目するとよいみたいですね。トランペットとクラリネットはどちらも高音のメロディ楽器なので、両者の振幅の相関は高いだろうとか、相関が乱れたらどちらかが演奏を誤っているだろうとかそんな感じでしょうか。というかトランペットとクラリネットって実際同じようなタイミングで鳴るんでしょうか。

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

知らないよ…。

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

ただ変数間の依存関係に着目するには「直接相関」と「間接相関」の区別が重要ということで…130~131ページにある「都市ごとの教会の数と殺人件数の相関(どちらも都市の人口に比例しているだけ)」のような例って岩波DS3でも見ましたね。

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

小学生の握力と算数ドリルの点数の疑似相関だったっけ。学年が交絡因子で。

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

ああそういうのでしたね。でも、算数ってノートにねばり強く書き続ける力だと思うんですよね。つまり、握力が高まるほど長時間ノートに力強く数式を書き続けられるようになって、算数力も向上するのではないでしょうか。

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

話ややこしくしなくていいよ…。

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

132ページからもしばらく時系列からは離れて(?)、確率ベクトルの対マルコフグラフをデータからどう描けばいいのか、という話が続きますね。それでもうこの本では多変量正規分布を仮定したガウス型グラフィカルモデルのみ取り扱うようですね。132ページの「多変量正規分布というめがねをかけて」という表現がわかるようでわかりませんが置いておきますね。…あれ、データに多変量正規分布を仮定したときの共分散行列の推定値ってこの (2.4) 式でいいんでしたっけ。最近別件でアンサンブルカルマンフィルタを実装したんですけど、そちらではフィルタを実行するときアンサンブル粒子から以下の下側の式のように分母が N-1 の共分散行列を計算して、これを用いてカルマンゲインを計算したと記憶しているんですが…。

 \displaystyle \hat{\Sigma} \equiv \frac{1}{N} \sum_{n=1}^N (x^{(n)} - \hat{\mu})(x^{(n)} - \hat{\mu})^{\top} \tag{2.4}
 \displaystyle \hat{\Sigma} \equiv \frac{1}{N - 1} \sum_{n=1}^N (x^{(n)} - \hat{\mu})(x^{(n)} - \hat{\mu})^{\top}

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

前者が最尤推定量で後者が不偏推定量だね。これらの違いは、雑に言うと以下のような感じかな…。

  • 手元にデータがある。色々な多変量正規分布の中から、手元のデータを生成する確率が最も高いものを選ぶ。これは対数尤度を最大化すればよい(計算略)。
  • 手元にデータがある。このデータがある多変量正規分布から生成されたということを踏まえる=つまり、手元のデータはたまたま手元にあるように生成されただけで、別の世界線では同じ多変量正規分布から生成された別のデータになっていたかもしれない。これを踏まえた上で、色々な多変量正規分布の中から、手元のデータを生成する確率が最も高いものを選ぶ。手元のデータを母集団と見做すならば共分散行列は式 (2.4) の標本共分散行列でいいのだが、「標本平均がたまたま観測された量である」という仮定の下では、この「標本平均(データに対して最も都合がよい場所に仮定した平均)の周りのばらつき」でばらつきを見積もるのはやや過小評価になる(証明略)。
最尤推定量の方は、データが fix された上でそれを多変量正規分布にぼかして考えたらどうなる、って感じで、不偏推定量の方は、データは真の分布からたまたまこのようにサンプリングされただけって感じがするんだよね。
  • 異常検知におけるガウス型グラフィカルモデルの推定では、手元のデータを生成した最も尤もらしい共分散行列を知りたいので最尤推定量を求める。それを通して手元のデータの変数間の相関構造を知りたいだけなので、サンプリングが違ったらどうだっただろう、という考慮はそぐわない。
  • アンサンブルカルマンフィルタでは状態分布を代替的にアンサンブル粒子で表現していて、フィルタのステップではアンサンブル粒子を共分散行列に翻訳したい。このとき、現在のアンサンブル粒子はたまたま今回の計算のために生成しただけなので、現在のアンサンブル粒子を生成した最も尤もらしい共分散行列を知りたいのではない。色々なアンサンブル粒子の選び方があることを踏まえて、不偏推定量を求める。
長々と書いてみたけど N が大きかったら N と N-1 の違いなんて大したことないし、そもそも10章通して N-1 を N に読み換えても話変わらない気がするんだけどね…。

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

お疲れさまです。ところで、133ページを読むと、結局  x_i x_j の間に直接相関がないことと、精度行列の  (i,j) 成分がゼロであることが同値なんですね? ということは、「教会の数」「殺人件数」「人口」のデータの場合、「教会の数」と「殺人件数」の組合せに対応する精度行列の要素がゼロになるということですが、本当にちゃんとそうなるんでしょうか…。疑似相関とはいえ相関はあるんですよね? あとこれって一意に定まらなさそうな気もするんですよね。だって「教会の数」「殺人件数」「人口」のどれがニワトリでタマゴかなんてデータは語らないですよね?

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

今回のケースでは、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]]

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

確かに今度は (0, 1) 要素=(1, 0) 要素が小さいですね。偏相関係数をきちんととれるものなんですね。それで、134ページからは、この非対角要素に偏相関係数(に比例する量)をもつ精度行列をなるべく疎にしたいという話になるんですね。精度行列において閾値未満の非対角要素をゼロに切り捨てるというアイデアは駄目で、そもそも共分散行列が正則でなく精度行列が求まらなかったり、精度行列が求まっても閾値に敏感だったり切り捨てた結果の精度行列が正定値でなくなったりして多変量正規分布の精度行列ではなくなってしまうので駄目だとありますね。じゃあなんでそのアイデアを出してみたんですかね。

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

いや、素朴な発想が上手くいかないと断っておくことも大事だからね? 何にせよ、そういう要望なら目的関数に正則化項を入れて多変量正規分布をフィッティングするんだろうね。

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

なぜネタバレするんです…。おっしゃる通り、尤度関数に精度行列の要素の和のラプラス分布を乗じた上でこれを最大化するんですね。対数をとると単純に対数尤度関数から  \rho || \Lambda||_1 をマイナスすることになります( \rhoラプラス分布の尺度の逆数で、 || \Lambda||_1 は精度行列の各要素の絶対値の和ですね)。この項のことをしばしばL1正則化項と呼称すると。しかし、このL1正則化項の導入によって、精度行列を推定するのに単に標本共分散行列の逆行列を取ればよいというわけにはいかなくなってしまいました。どうやって精度行列を求めればいいんでしょう…。

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

精度行列に関する(L1正則化項付き)対数尤度の偏微分をとればいいんじゃない。

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

しかし137ページに  \rho>0 の場合は「偏微分=0」が解析的に解けないとありますね。 \rho<0 だったら解けるんでしょうか…そもそも今回は精度行列を疎にしたいという目的より必ず  \rho>0 なので、 \rho<0 の場合に興味はありませんが…。それで、解析的に解けないのでどうするかというと…  j 次元目以外固定してしまう?  j 次元目を一番後ろに寄せて、精度行列、精度行列の逆行列、標本共分散行列を137ページのようにブロック表記して、そうすると各行列が4つの部分に分解されて、うち2つは転置の関係なので、元の「偏微分=0」は (10.15) (10.16) (10.17) の3式になるわけですね…それで138~139ページの手続きにしたがえば元の精度行列の  j 次元目の行・列以外を固定したもとで目的関数を極大値にできますが…それを全ての次元について、ぐるぐると繰り返せば目的の精度行列が求まると…このように、元データを多変量正規分布のめがねでみてかつ疎な精度行列を得ようとする手法を「グラフィカルラッソ」というんですね。うーん、手順はわかるんですが、これ最適解に収束するんですか? ある次元以外固定してしまおう、というのがなんか心配なんですが…。

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

その辺の議論はこの本にはなさそうだね。参考文献を参照かな。

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

…まあいいです。何らかの参照データから疎な精度行列が求まり、変数間の依存関係がわかったとしましょう。これを利用する異常検知手法が 10.5 節に2例紹介されています。1例目は、ある多変量の観測値が得られたときに、その  j 次元目の要素の異常度は、  j 次元目以外の要素を given とした条件付き確率密度の対数のマイナス1倍とする(外れ値解析)。

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

つまり、時系列の j 次元目だけに着目してよく観測される量だったとしても、j 次元目以外がこの状況のときには観測されない量だ、という場合は異常と判定されるようなイメージだね。

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

2つ目の例は1点の観測値だけでなくまとまった観測値に対する異常判定ですね。このときは、参照データと観測データのそれぞれの分布を求めて、両者間のKL情報量を j 次元目以外について積分したものが異常度ということです。どちらの例も特にグラフィカルラッソや、何なら多変量正規分布でなくてもよさそうですね。ただ相関構造をシンプルにしておかなければ有効な解析ができないということなのですかね。多変量正規分布の場合の (10.28) → (10.29)、(10.30) → (10.31) の式変形は後で追っておきましょう。…143ページに「 \mathcal{D} \mathcal{D'} の立場を入れ替えた定義も可能です」とありますが、確かにKL情報量は非対称ですが、実際過去側の分布と現在側の分布のどちらをどちらにすべきなんでしょう?

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

9章を読んだときにも似た話をしたね。KL情報量は対数密度比の期待値だから、過去側の分布で期待値をとるのか現在側の分布で期待値をとるのかって話だと思うけど…過去側の分布が基準と考えれば (10.30) 式の向きが自然なんだろうけど。異常時に分布がどう変化するかにもよるだろうね。

(他の章のメモがあれば)つづく

雑記

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

機械学習で Leakage ってよく聞きますけど、結局 Leakage の語義というか具体的に言い切ると何なのかよくわからないんですよね。漏れてはいけない情報が漏れていること、のようではあるんですが、何が漏れてはいけないんでしょう? 我々は総力を挙げてテストデータの正答を予測するモデルの構築に取り組まなければならないのではないんでしょうか?

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

ここ 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 は「漏れてはいけない情報が漏れた下で予測が実行されること」で、というからには漏れてはいけないものとは何か定められていないといけないけど、その典型例は「テストデータの正答」「モデルを利用する場面ではまだ入手できない未来の情報」「その他のモデルが知り得ない情報」「その他、学術的な目的などでこの情報からは学んでほしくないという情報」とかみたいだね。でも何か予測を実施するときの制約というのはケースバイケースだから、具体的に統一的にこれが漏れていてはならないとは言いづらいんじゃないかな。

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

なるほど。その未来の情報という例で思ったんですが、こういうのはどうです? 未来の天候情報そのものを参照するとズルなので、未来の天候情報を予測するモデルを謎の目的関数で学習して、そのモデルの予測を利用して未来の売上を予測するというのは。

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

いや私運営じゃないから知らないんだけど…確かにモデルがあからさまにデータを内蔵していなくて自分で天候を予測しているなら問題は少なそうだけど、それで達成できるような精度は直接未来の売上を予測しても達成できるんじゃないの…? もちろんモデルによるけど…。だいたい中間目標を置くのは多くの場面で得策じゃないよね。天候情報の予測がよほど楽勝だったら別だけど、じゃあやっぱり大してプラスにならないんじゃないの。

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

ぐぬぬ…。

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

上の Kaggle のページの続きに Leakage の具体例があるね。

  • 各患者が前立腺ガンかどうかを予測するためのデータに「前立腺の手術を受けたか」という変数を含めてしまった。

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

ちょっとうっかりしすぎじゃないですかね。因果関係ってご存知ですかって感じなんですけど。

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

いやでも、こんなことも実務で起こらないとはいえないから気を付けないとねってことなんじゃないかな。これは極端な例で、実際の Leakage はもっとわかりづらいって書いてあるよ。あと続く例は、実際に初期の Kaggle のコンペティションで、ソーシャルネットワークの「つながり」の予測問題があったみたいなんだけど、サンプルスクリプトのミスで、ある条件を満たすユーザどうしは無条件で「つながっている」というラベルがついてしまっていたみたいで、これを突いたチームが2位になったみたいだね。さらに1位のチームは、スクレイピングによってデータの匿名性を剥がしちゃったみたいだけど。

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

…それって出題ミスですよね。2位のチームは、出題者の意図とは違うモデルを構築してしまったとはいえ、データのパターンを見出したのだから Leakage ですらないような。1位のチームは…大会のレギュレーションによるでしょうね。…最後の Leakage in Competitions というパラグラフにはなんて書いてあるんです?

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

コンペティションにも Leakage がない保証があるようになんてできない。以下のような対処はできるけどいつもできるわけじゃないって書いてあるね。

つづかない

雑記: 沖本本読み勉強会用のメモ

以下の本を読む勉強会用のメモです。もし何か誤りがありましたらご指摘いただけますと幸いです。

経済・ファイナンスデータの計量時系列分析 (統計ライブラリー)経済・ファイナンスデータの計量時系列分析 (統計ライブラリー)
沖本 竜義

朝倉書店 2010-02-01
売り上げランキング : 19813

Amazonで詳しく見る
by G-Tools
f:id:cookie-box:20180513082851p:plain:w60

今度勉強会で上の本の3章を発表することになりました。

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

へー。1章と2章は読んだの?

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

読んでないです。

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

読もうよ!?

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

仕方ないですね…1章で思い出しておきたい事項は以下でしょうか。

  • 時系列分析では時系列データの背後にデータ生成過程(または単に過程)を仮定する。
  • 期待値が時刻tに依存せず一定で、自己共分散も時刻tに依存せずラグkにのみ依存するような過程を弱定常であるという。
  • 特に、ある長さkの部分時系列の同時分布が時刻tに依存せず同一の分布にしたがうような過程を強定常であるという。
    • 正規過程が弱定常であれば、強定常である。
  • 期待値が時刻tに依存せずゼロ、0次の自己共分散(つまりただの分散)が時刻tに依存せず一定、1次以上の自己共分散が時刻tに依存せずゼロであるような過程をホワイトノイズという。
  • データの自己相関の有無を判断するには、以下のような検定をする。
    • ある次数の標本自己相関係数に対して「=0」の帰無仮説を「≠0」の対立仮説に対して検定する。
    • ある次数までの全ての標本自己相関係数に対して「全て0」の帰無仮説を「どれかは0でない」の対立仮説に対して検定する。
…9ページの脚注がよくわかりません。正規過程とは、ある長さkの部分時系列が多変量正規分布にしたがう(多変量正規分布の期待値と共分散は時刻に依存してもよい)ような過程なんですよね。任意の時刻tに対して y_t の周辺分布が正規分布ならば、そのような過程は正規過程なのではないでしょうか? それぞれ正規分布にしたがう確率変数をいくつか並べたら、その確率ベクトルは多変量正規分布にしたがいますよね?

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

とは限らないかな。例を示そうか。

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))

f:id:cookie-box:20180618221755p:plain
上の3つの例はどれもXの周辺分布が N(0, 4)、Yの周辺分布も N(0, 4) だけど、一番右のプロットの形状は2変量正規分布にみえる?

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

一番右のプロットの形状ですか? 気持ち悪いです。

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

ごめん…。

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

次は2章ですね。ここでMA過程とAR過程が降って湧いたように基本的なモデルとして登場しますが、さもそれらが自己相関をもつ確率過程を表現しようとした自然な帰結なのだというような雰囲気ですが、自己相関をもつ確率過程を表現するのに「y_{t} と y_{t-1} が共通の確率変数に依存する」と「y_{t} が y_{t-1} を参照する」以外のモデルは本当にないんでしょうか? y_{t-1} が y_{t} を参照するとか。

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

未来の値が現在の値を決定するモデルはやめようよ?

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

データ生成過程さんにしてみれば、エントロピーの増大方向など知ったことではなくないですか?

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

仮にデータ生成過程さんが気にしなくても、そのモデルを利用する我々が気にするんだよ!

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

なるほど。…それで本に戻ると、MA過程の方は次数・パラメータにかかわらず定常であると。これは直感的に理解できますね。有限の過去までのノイズの線形結合で爆発しようがありませんから。それでAR過程の方は (2.15) 式の全ての解の絶対値が1より大きいときに定常となると。なぜそうなるかのやわらかい解釈は 2018-03-07 の記事を参照してください。また、MA過程の問題点として、ある過程を記述する係数とホワイトノイズの分散のセットの候補が複数存在してしまうので、AR(∞)過程に書き直せる(=反転可能である)ものを選択しようという基準があるということですね。

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

「同一の期待値と自己相関構造をもつMA(q)過程が2^q通り存在し、その内1つのみが反転可能である」というのをきちんと追えてないけど、「MA過程はAR(∞)過程に書き直せる」「だから少なくともAR(∞)過程に書き直せるようなホワイトノイズの取り方が1つ存在する」というのは理解できるかな。

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

2章の後半はモデルの推定と解釈の話ですが…43ページに、ARMAモデルやGARCHモデル、マルコフ転換モデルなどは最小2乗推定できないってありますけど、なんででしょう? 何かモデルを仮定すれば残差平方和を計算することは常にできますよね?

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

まずARMAモデル(というよりMAモデル)については、ある時刻での残差が次の時刻の説明変数になっちゃうよね。だから「理想的には残差がゼロだろう」っていうスタンスのOLSで推定したら駄目だよね。GARCHとマルコフ転換モデルはもっと後の方で出てくるからすぐ追えないけど、GARCHはぶれの大きさをモデリングするし、マルコフ転換モデルは観測できない真の状態がぶれるモデルだし、やっぱり「理想的にはぶれないだろう」っていうOLSの考え方はそぐわないようにみえる。

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

あー、確率的な変動が全くないとすると意味をなさないようなモデルだと駄目っぽいですね。それでやっと3章に辿り着きましたが、この章のタイトルは「予測」です。しかし予測とは何なのかと。

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

大切な問いだね。

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

AR(1)過程を例に挙げてこれを問うていますが…「AR(1)過程と現時点の観測値を所与とする。次の時点の予測値とは何か」と。「予測A」はAR(1)過程の期待値そのものを予測値としていて、入手している現時点の観測値を完全に無視していくスタイルですね。「予測B」はノイズの最頻値、「予測C」はノイズの平均値が次の時刻に実現するだろうというような予測です。別にどの予測もそういう予測ではあるわけですが、ここでは(一期前までに観測される値が所与の下での)予測誤差の2乗の期待値が小さいほどよい予測とし、予測誤差の2乗の期待値を最小にする予測を最適予測とよぶということです。それで、先の例では「予測C」が最適予測だと…ん、A~Cの3つの予測の中では「予測C」が最もよいというのはわかるんですが、「予測C」よりもよい予測はないことは示されてないですよね?

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

「予測C」が「一期先の値の期待値を予測値とする」として得られたことは、それが予測誤差の2乗の期待値を最小に予測であることの説明になってないってことかな。 E\bigl((X-a)^2 \bigr) a について微分してみたらどうだろう?

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

思い出しました。2次損失は期待値で最小に、線形損失は中央値で最小に、0-1損失は最頻値で最小になるとかDLMの本でやりましたね。…その後は、AR過程の何期も先のMSEを陽に書くのは難しい、区間予測はシミュレーションでやろう、というような話ですね。…なんか3章はそんな感じですね。

(次回があれば)つづく

異常検知と変化検知: 9章メモ(部分空間法による変化検知)

以下の赤い本の9章を読みます。キャラクターは適当です。誤りがありましたらご指摘いただけますと幸いです。

異常検知と変化検知 (機械学習プロフェッショナルシリーズ)異常検知と変化検知 (機械学習プロフェッショナルシリーズ)
井手 剛 杉山 将

講談社 2015-08-08
売り上げランキング : 56236

Amazonで詳しく見る
by G-Tools
他の章のメモ:まだ
f:id:cookie-box:20180513082851p:plain:w60

「変化を検知したい」ってときありますよね。この本の9章を読みたいです。

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

導入が雑…。

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

最初の累積和法は、「正常状態の分布」「異常状態の分布」を所与として、逐次的に観測される変数が異常状態の方に近い場合にだけ異常さのポイントを積み上げる手法のようですね。…「正常状態の分布」「異常状態の分布」がわかってたら苦労しなくないですか?

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

いやまあそこがわかってたらかなりありがたいけど、じゃあどちらの分布に近いのかというのを対数尤度比という形で定量化して、さらにそれを時系列データにおける継続的な異常状態の検知に適用するために累積和という指標にまとめる、というのも大事な仕事なんじゃないかな。

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

そう、その「継続的な異常状態の検知」という点についてもなんですが、この累積和って異常の継続性を考慮してますか? 「たまに異常になる分には目をつぶって、継続的に異常になったら検知したい」ということで累積和をとったということですが、これって異常が連続していようと数回おきであろうと足し上げていくだけですよね。異常が継続していなくても累積和はどんどん降り積もっていきません?

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

そういう場合は定期的にリセットするんでしょ。結局ある一定の期間内の異常さの和に閾値を置いているのと等しいとは思うよ。ある一定の期間内に多くの異常が発生することがここでの継続的な異常なんだと思う。

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

なるほど。…次に紹介されているのは近傍法ですね。この手法で少なくとも所与とするのは「正常な部分時系列のセット」ですね。それで、検査対象の部分時系列を入力すると、最も近い正常な部分時系列までの距離(の対数の定数倍)を返すといった具合ですね。

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

k=1 の k 近傍法で、明示的に異常標本を与えずにとにかく異常度のみ算出する場合はその手続きになるね。

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

115ページから、ある時点(まで)のデータが異常かどうかを判定するというよりは、「検査データセット」が「参照データセット」と異質か否かという話になっていますね。その1つの抽象的な定式化が、まず「参照データの分布 p」と「検査データの分布 p'」を求めた上で、2つの分布の相違度を出すということですね。「参照データの分布」と「検査データの分布」の相違度は例えば (9.8) 式のように定義できると…対数尤度比の期待値ということですが、どうも既視感がある数式です…これは何かの予兆なんでしょうか…。

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

いや、カルバック・ライブラー情報量ですね。

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

はい。

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

KL情報量は非対称だけど、これは「参照データの分布が厚いところで検査データの分布が薄いことを問題視する」って向きだね。参照データにはよく見られるようなデータが検査データに全然見られなかったらおかしいだろうと。狭い参照データ分布をカバーするように広い検査データ分布よりは、広い参照データ分布をカバーし切らない狭い検査データ分布の方が相違度が大きくなるね。例えば、
   D_{\rm KL}(p=N(0, 1^2),\,  p'=N(0, 2^2)) \approx 0.312 (参照データ分布の方が狭い場合)
   D_{\rm KL}(p=N(0, 2^2), \, p'=N(0, 1^2)) \approx 0.807 (参照データ分布の方が広い場合)

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

参照データ分布としてそんなにだらだら広がった分布を仮定しない方がよさそうですね。検査データ分布側の密度の欠落に厳しくなりそうです。

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

116ページからは、(9.8) 式の具体的な計算手法が1つ紹介されていて、つまり、特異スペクトル変換法だね。

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

フォンミーゼス・フィッシャー分布…扱ったことがないですね…7章を参照するとこれは、M次元超球面上に定義される分布ということですが…なぜこんな分布を導入するのでしょう…。

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

参照データ分布と検査データ分布にフォンミーゼス・フィッシャー分布(長いので以下vMF分布)を仮定するということだけど、(少なくともこの本の)特異スペクトル法の手順の中で「集中度」は推定せず参照データ分布も検査データ分布も適当な同一の集中度を取ると仮定しているから、結局推定対象パラメータは「平均方向」のみで、だから「参照データセットと検査データセットをそれぞれ表現するvMF分布って何だろう」というよりは、実質的には「参照データセットと検査データセットをそれぞれ代表する単位ベクトルって何だろう」ってことだね。そういう単位ベクトルを求めたらそれらを平均方向とするvMF分布どうしのKL情報量を相違度とするわけだけど。ともあれ単位ベクトルの求め方が肝心そうだね。いや…単位ベクトル自体を求めるわけじゃなくて、そのような単位ベクトルどうしの内積に相当する値を直接求めにいってるね。

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

その手順は…参照データを横に並べた行列(M×n)と、検査データを横に並べた行列(M×k)をそれぞれ特異値分解して、それぞれからM次元の左特異ベクトルがM本ずつ出てきますが、特異値上位の何本かずつを抜き出してきて(参照データと検査データから違う本数ずつでもよい)横に並べた行列を  U, \, Q として、 U^{\top}Q の行列2ノルムこそが求める内積に相当する値であると…はて…。

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

参照データの第i成分を全て束ねたM本のベクトルをどのようなM個の係数の組み合わせで線形結合するとそのL2ノルムが極値になるかを考える(ただし、M個の係数の2乗和は1という制約の下で)。極値を与える係数のセット(=左特異ベクトル)がMセット求まる。最も大きい極値を与える左特異ベクトルは参照データを特徴付けるベクトルだろう、って感じだね。だから参照データのそれと検査データのそれの内積を取れば相違度がわかるだろうと。ただ最も大きい極値以外無視するというのもあれだっていうなら特異値上位の何本かずつ左特異ベクトルを選ぶわけだけど、例えば2本ずつ選んだら内積が4パターン計算できちゃう。どうやってトータルの相違度とすべきかわからない。そのときは、ここでまたしても各行ベクトル(列ベクトルでもいいけど)の各成分をどのように線形結合するとL2ノルムが極値になるかという問題を解いて、最大の極値(最大特異値)をトータルの相違度とするんだね。この操作は  U^{\top}Q の行列2ノルムを取ることと等しい。

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

うーん…そうやって求めた左特異ベクトル(たち)がその元データセットならではのベクトル(たち)というのはいいんですが、「元データセットをよく代表するベクトル(たち)です」っていうのはよくわかんないですね…。いやでも…元データセットが仮に2次元のデータたちで、どれも成分どうしが (3x, 4x) のような比率をもつデータだったとしたら、特異値最大の左特異ベクトルとしてはおそらく (3/5, 4/5) のようなものが選ばれるはずで、そういう方向との内積を取った和が最も大きくなるデータたちですよ、というのはデータセットをよく形容するベクトルのようには思えます。

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

それが118ページでいう「一番人気のある方向」だね。

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

ちょっと図示してみたいですね。ランダムな2次元データセットを生成して、左特異ベクトルをプロットしてみましょう。こんな感じでいいんでしょうか。灰色のベクトルたちを代表するのが赤いベクトルになると期待しているんですが…。

n <- 7
x <- c()
for (i in 1:n) {
  t <- runif(1, min=0, max=2)
  x <- c(x, 2*cos(pi*t), 2*sin(pi*t))
}
x <- matrix(x, nrow=2)

plot(c(-2,2), c(-2,2), xlab="", ylab="", type="n")
for (i in 1:n) {
  lines(c(0, x[[1,i]]), c(0, x[[2,i]]), lwd=2, col="dimgray")
}
res <- svd(x)
lines(c(0, res$u[[1,1]]), c(0, res$u[[2,1]]), lwd=4, col="red")
lines(c(0, res$u[[1,2]]), c(0, res$u[[2,2]]), lwd=4, col="darkorange")
print(res)

データセットを全て長さ2のベクトルにしているのは見やすさのためだけです。長さ1だと左特異ベクトルと重なったときに見えなくなっちゃうし、長さばらばらなのも見づらかったんで。
f:id:cookie-box:20180613230237p:plain:w240
122ページ以降は特異スペクトル変換法の高速化の話ですね。ランチョス法という手法を用いるそうです。123ページの定義9.2は…計算せよ? ああ、こういう問題を定義するというわけですか…。なるほど、この問題を効率よく計算する手段があれば、効率よく特異スペクトル変換が実現でき…ないですよね? 「参照データの左特異ベクトルとの内積を取ろう」という時点でその相手である検査データの左特異ベクトルが既にあればいいですが、それを求めるのにやっぱり特異値分解が必要ですよ?

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

だから、特異スペクトル変換の定義を modify するみたいだね。もう検査データ側の特異値分解はしないことにして、代表的なデータを選んで(全て選んで後で相違度を何らかの方法で集約してもいいけど)左特異ベクトルの変わりにするみたいだね。

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

「一番人気のある方向」とは何だったのか…。でもランチョス法は近似的な3重対角化なんですね。ただネットで調べるとループ計算中に a_j が数値誤差で直交しなくなっていくことに注意とかありますね。でも数値計算の授業でランチョス法ってやりましたっけ。QR法を習ったのは覚えているんですがその前後にやってたんですかね。もう10年くらい前ですから記憶がありませんね。

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

いや部長は高校1年生だよね!? 幼稚園の年長に数値計算の授業はないよね!?

(他の章のメモがあれば)つづく