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

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

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

朝倉書店 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年生だよね!? 幼稚園の年長に数値計算の授業はないよね!?

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

雑記

Keras をよくわかっていないのでおかしな点があればご指摘いただきたいのですが、1つ前までの記事で紹介していた DAGMMメモ1メモ2メモ3メモ4)を Keras で実装しようとすると以下のような手順でできるんでしょうか。よくわかりません。

  • 不整脈データ(274次元)を真ん中の特徴量が2次元の深層自己符号化器にかける。
  • 深層自己符号化器が終わったところで、最初の入力と特徴量を回収してきて、損失に再構築誤差を流し込むと同時に、特徴量に再構築エラー(2次元)を concat してこれを出力する(カスタムレイヤー0)。
  • さっきの出力を何層かの全結合層にかけて各分布への割り当て確率にする(不整脈データは混合数2)。
  • ここで4次元の特徴量を回収してきて、尤度と正則化項を損失に流し込む(カスタムレイヤー1)。

ガワだけ実装すると以下のようなイメージですが、中身を実装してコンパイルできるのか知りませんし、訓練して勾配計算が意図通りに動くのかも知りません。

  # DAGMM
  x_org  = Input(shape=[274])
  h_0    = Dense(10, activation='tanh')(x_org)
  z_c    = Dense(2)(h_0)
  h_1    = Dense(10, activation='tanh')(z_c)
  x_dash = Dense(274)(h_1)
  z      = MyLayer0()([x_org, z_c, x_dash])    # 損失に再構築誤差を付加 & 再構築エラーをconcat
  h_2    = Dense(10, activation='tanh')(z)
  h_3    = Dropout(0.5)(h_2)
  gamma  = Dense(2, activation='softmax')(h_3)
  gamma_ = MyLayer1()([z, gamma])              # 損失にエネルギーと正則化を付加
  DAGMM  = Model(x_org, gamma_)
  DAGMM.compile(optimizer=optimizers.Adam(lr=0.0001), loss=None)
  DAGMM.summary()
_______________________________________________________________________________________
Layer (type)               Output Shape      Param #     Connected to
=======================================================================================
input_1 (InputLayer)       (None, 274)       0
_______________________________________________________________________________________
dense_1 (Dense)            (None, 10)        2750        input_1[0][0]
_______________________________________________________________________________________
dense_2 (Dense)            (None, 2)         22          dense_1[0][0]
_______________________________________________________________________________________
dense_3 (Dense)            (None, 10)        30          dense_2[0][0]
_______________________________________________________________________________________
dense_4 (Dense)            (None, 274)       3014        dense_3[0][0]
_______________________________________________________________________________________
my_layer0_1 (MyLayer0)     (None, 4)         0           input_1[0][0]
                                                         dense_2[0][0]
                                                         dense_4[0][0]
_______________________________________________________________________________________
dense_5 (Dense)            (None, 10)        50          my_layer0_1[0][0]
_______________________________________________________________________________________
dropout_1 (Dropout)        (None, 10)        0           dense_5[0][0]
_______________________________________________________________________________________
dense_6 (Dense)            (None, 2)         22          dropout_1[0][0]
_______________________________________________________________________________________
my_layer1_1 (MyLayer1)     (None, 2)         0           my_layer0_1[0][0]
                                                         dense_6[0][0]
=======================================================================================

カスタムレイヤー0は以下です。再構築エラーをつくって低次元の特徴に concat すればいいと思うんですが、まだつくっていないのでなので低次元の特徴を2つ concat しています(不整脈データの例では低次元の特徴と再構築エラーが同じ2次元なので)。雰囲気です。

from keras import backend as K

# 損失に再構築誤差を付加 & 再構築エラーをconcat
class MyLayer0(Layer):
  def __init__(self, **kwargs):
    self.output_dim = 4
    super(MyLayer0, self).__init__(**kwargs)
  
  def reconstruction_loss(self, x_org, x_dash):
    return metrics.mean_squared_error(x_org, x_dash) 
  
  def reconstruction_error(self, x_org, x_dash):
    z_e = None # IMPLIMENT ME
    return z_e
  
  def call(self, inputs):
    x_org  = inputs[0]
    z_c    = inputs[1]
    x_dash = inputs[2]
    loss = self.reconstruction_loss(x_org, x_dash)
    self.add_loss(loss, inputs=[x_org, x_dash]) # 再構築誤差を付加
    z_e = z_c # self.reconstruction_error(x_org, x_dash) # 再構築エラーを計算
    z = K.concatenate([z_c, z_e], axis=1)
    return z
  
  def compute_output_shape(self, input_shape):
    return(input_shape[1][0], self.output_dim)

カスタムレイヤー1は以下です。ここは損失にエネルギー項と正則化項を足すだけですがまだ書いていません。一応各分布への割り当て gamma をそのまま output していますが output する必要はないです。

# 損失にエネルギーと正則化を付加
class MyLayer1(Layer):
  def __init__(self, **kwargs):
    self.output_dim = 2
    super(MyLayer1, self).__init__(**kwargs)
  
  def custom_loss(self, z, gamma):
    return 0 # IMPLIMENT ME
  
  def call(self, inputs):
    z = inputs[0]
    gamma = inputs[1]
    loss = self.custom_loss(z, gamma)
    self.add_loss(loss, inputs=[z, gamma])
    return gamma

  def compute_output_shape(self, input_shape):
    return(input_shape[1][0], self.output_dim)

まだ実装してないのでわかりませんというか実装すればいいんですが気が向いたらやります。いつ気が向くのかはわかりません。


深層自己符号化器までは何も問題ないはずなので取り出して動かしてみます。

  # 不整脈データ
  df = pd.read_csv('arrhythmia.data', header=None)
  df = df.drop(df.columns[[10,11,12,13,14]], axis=1) # 欠損カラムを除去
  # 異常ラベルを貼り変える
  df.iloc[:,274][(df.iloc[:,274]>=1 ) & (df.iloc[:,274]<=2 )] = 20
  df.iloc[:,274][(df.iloc[:,274]>=3 ) & (df.iloc[:,274]<=5 )] = 21
  df.iloc[:,274][(df.iloc[:,274]>=6 ) & (df.iloc[:,274]<=6 )] = 20
  df.iloc[:,274][(df.iloc[:,274]>=7 ) & (df.iloc[:,274]<=9 )] = 21
  df.iloc[:,274][(df.iloc[:,274]>=10) & (df.iloc[:,274]<=13)] = 20
  df.iloc[:,274][(df.iloc[:,274]>=14) & (df.iloc[:,274]<=15)] = 21
  df.iloc[:,274][(df.iloc[:,274]>=16) & (df.iloc[:,274]<=16)] = 20
  df.iloc[:,274] = df.iloc[:,274] - 20
  
  x = df.drop(df.columns[[274]], axis=1).values
  y = df.iloc[:,274].values
  
  # 深層自己符号化器のみ取り出して学習
  DAE = Model(x_org, x_dash)
  DAE.compile(loss='mean_squared_error', optimizer=optimizers.Adam(lr=0.0001))
  DAE.fit(x, x, shuffle=True, epochs=10000, batch_size=128)
  
  # エンコーダを取り出してプロット
  encoder = Model(x_org, z_c)
  z_trained = encoder.predict(x)
  plt.figure(figsize=(12, 8))
  plt.scatter(z_trained[:,0], z_trained[:,1], c=y, cmap=plt.get_cmap("bwr"))
  cb = plt.colorbar()
  for t in cb.ax.get_yticklabels():
    t.set_fontsize(18)
  plt.tick_params(labelsize=18)
  plt.show()

上で学習すると特徴空間(2次元)は以下のようになりました。青が正常で赤が異常です。452点あったはずなんですが、相当重なっていると思います。

f:id:cookie-box:20180527165625p:plain:w630

あと再構築誤差方向への分布はどうなっているのかも気になると思います。さらに何らかの再構築誤差も含めてプロットしようとするともう3次元プロットになってしまうのでその準備をします。

from mpl_toolkits.mplot3d import Axes3D

def my_relative_distance(x, y):
  ret = []
  for i in range(x.shape[0]):
    e2 = 0
    x2 = 0
    for j in range(x.shape[1]):
      e2 = e2 + (x[i][j] - y[i][j]) * (x[i][j] - y[i][j])
      x2 = x2 + x[i][j] * x[i][j]
    ret.append(e2 / x2)
  return ret

その上でこうすると再構築誤差も含めてプロットできると思います。

  # エンコーダを取り出してプロット
  encoder = Model(x_org, z_c)
  
  x0 = df.loc[df.iloc[:,274]==0,:].drop(df.columns[[274]], axis=1).values
  x1 = df.loc[df.iloc[:,274]==1,:].drop(df.columns[[274]], axis=1).values
  z_trained0 = encoder.predict(x0)
  z_trained1 = encoder.predict(x1)
  r_error0 = my_relative_distance(x0, DAE.predict(x0))
  r_error1 = my_relative_distance(x1, DAE.predict(x1))
  
  fig = plt.figure(figsize=(12, 8))
  ax = Axes3D(fig)
  ax.plot(z_trained0[:,0], z_trained0[:,1], r_error0, "o", color="#0000ff")
  ax.plot(z_trained1[:,0], z_trained1[:,1], r_error1, "o", color="#ff0000")
  plt.tick_params(labelsize=18)
  plt.show()

なんか再構築誤差方向にも異常データと正常データがあまり分離されているようにみえないので、深層自己符号化器のみで分離する場合は適するハイパーパラメータが違うんだと思います。

f:id:cookie-box:20180527171437p:plain:w720

論文読みメモ: 深層自己符号化器+混合ガウスモデルによる教師なし異常検知(その4)

以下の論文を読みます。

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-
※ キャラクターに元ネタはないです。お気付きの点がありましたらお手数ですがコメント等にてご指摘ください。
前回:その3 / 次回:まだ
f:id:cookie-box:20180513082851p:plain:w60

前回までのあらすじですが、今回 DAGMM という異常検知向けに多次元データの密度を推定するモデルを開発したので、通信データや医療データの4データセットを用いて異常検知性能を検証します。DAGMM の比較対象となるのは既存の state-of-the-art な4モデル及び DAGMM のフレームワークの各要素の検証となる variant な7モデルの計11モデルです。DAGMM やデータセット、比較対象11モデルについては前回までのメモを参照してください。

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

そうそう、異常検知って、「そのサンプルに対する確率密度が閾値以下なら異常」みたいな判定をすると思うんだけど、複数のモデルを比較検証するときはどんな風に閾値を決めるんだろ。同じ閾値ではフェアな比較ができなさそうだよね。そもそも密度推定によらない異常検知手法もあるし。

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

それも含めて、検証実験のコンフィギュレーションを説明しましょう。今回は2つの実験を行うのですが、1つ目の実験では、各データセットをまず訓練用とテスト用に2分割します。そして、訓練用データから正常データのみを取り出してモデルを訓練します。そのモデルをテストデータに適用するのですが、今回は教師なし学習なので、モデルの出力は各サンプルが「正常か異常か」ではなく「どれくらい異常か」ですね。そこで、今回は次のようにします。「各データセットの異常データ混入割合をa%としたとき、テストデータ異常度の上位a%が『異常』、その他が『正常』と判定されたとする。」 例えば、3つ目の不整脈のデータセットには異常データが15%含まれているので、異常度の上位15%が『異常』判定ということです。こうすれば、通常の分類問題の評価指標がそのまま利用できますね。つまり、精度、感度、F値で各モデルを比較します。

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

なるほど、それならある程度フェアな比較になるのかな。

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

それでその1つ目の実験の結果ですが、論文10ページに数表がありますが、見やすいようにF値のみグラフにしてみました。

f:id:cookie-box:20180524193826p:plain:w700

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

わざわざ!? 甲状腺データと不整脈データの方が通信データよりも異常データを分離する難易度高めなんだね。あと KDDCUP の E2E-AE は何があったし…。

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

さらに2つ目の実験です。1つ目の実験との違いは、訓練時に正常データのみを利用するのではなく、わざと異常データをコンタミします。KDDCUP データセットに対する結果が以下です。こちらもF値のみですがまた勝手にグラフを描いてみました。ratio は異常データのコンタミ率です。0% のF値は上のグラフの該当するF値と同じです。このグラフから提案手法 DAGMM はコンタミにも強いことがわかります。

f:id:cookie-box:20180524201948p:plain:w640

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

あー、OC-SVM は与えられた訓練データを囲む面を求めようとするから、そこに異常データがコンタミしていたら異常検知としての性能がガタ落ちしちゃうって感じなのかな。

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

4.5 節に低次元表現の可視化がありますね。Figure 3. (a) の DAGMM の低次元表現は、他のモデル (b) (c) (d) よりも正常データ(青)と異常データ(赤)をよく分離できていると思いませんか?

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

まあ、(a) が比較的赤と青が混ざっていないようには見えるけど…3次元プロットじゃ他の角度からも見てみないとよくわかんないような。

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

あとこの節で取り上げられているのは、Figure 3. の (c) や (d) は事前学習をせず結合学習する手法であるにもかかわらず、事前学習する (b) に比べて低次元表現にあまり変化がみられないという点ですね。ということからも、非線形な次元削減を最初から regularization して学習するのが肝要だと。しかも、(a) DAGMM の自己符号化器部分の再構築エラーは事前学習した自己符号化器並みに小さいということなので、DAGMM は自己符号化器の学習のフレームワークとしても優秀なのではないかと書いてあるようなないような気がしますね。

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

実際にベンチマークに用いられたデータセットを見てみましょう。どのデータセットも Web 上で見つかると思いますが、一番サイズが小さいところで不整脈のデータは以下に落ちてました。
UCI Machine Learning Repository: Arrhythmia Data Set
これは一番最後のカラムが不整脈の種類ラベルみたいですね。それで、論文には dimension が 274 とあるのにラベルを除いても 279 カラムあってどういうことだろうかと思ったんですが、11~15カラム目に欠損値があったのでおそらくそこを除いていると思います。なので欠損カラムを取り除いてとりあえず自己符号化器にかけてみました。ニューロン数も活性化関数もバッチサイズもエポック数も全て論文に明記されていますからね。あ、論文には tensorflow で実装したとありますが、さしあたり keras で書きます。

# -*- coding: utf-8 -*-
import numpy as np
import pandas as pd
from keras.models import Sequential
from keras.layers import Dense

if __name__ == '__main__':
  df = pd.read_csv('arrhythmia.data', header=None)
  df = df.drop(df.columns[[10,11,12,13,14]], axis=1)
  y = df.iloc[:,274].values
  x = df.drop(df.columns[[274]], axis=1).values
  print(y.shape)
  print(x.shape)
  
  encoder = Sequential()
  encoder.add(Dense(10, activation='tanh', input_shape=(274,)))
  encoder.add(Dense(2))
  
  decoder = Sequential()
  decoder.add(Dense(10, activation='tanh', input_shape=(2,)))
  decoder.add(Dense(274))
  
  compression_net = Sequential()
  compression_net.add(encoder)
  compression_net.add(decoder)
  compression_net.compile(loss='mean_squared_error', optimizer='RMSProp')
  compression_net.fit(x, x, shuffle=True, epochs=10000, batch_size=128)

これで学習できます。まあ本当は全データを訓練にかけるんじゃなくて、訓練用とテスト用に分割した上で訓練用から異常データを除いて学習しなきゃなんですけど。そもそも、自己符号化器だけ学習しちゃ駄目だし。

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

じゃあもうちょっと頑張ろうよ!?

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

今日はもう遅いので…。ところで、R を書いた直後に Python を書くと混乱しませんか?

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

知らないよ…。

(次回があれば)つづく

論文読みメモ: 深層自己符号化器+混合ガウスモデルによる教師なし異常検知(その3)

以下の論文を読みます。

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-
※ キャラクターに元ネタはないです。お気付きの点がありましたらお手数ですがコメント等にてご指摘ください。
前回:その2 / 次回:まだ
f:id:cookie-box:20180513082851p:plain:w60

前回までのあらすじです。

f:id:cookie-box:20180521234046p:plain:w600

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

あらすじにはなってないよね!?

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

ここで一旦2節の Related Work に戻りましょうか。高次元データに対する異常検知の先行手法は3タイプに分類されると書いています(この辺りの話、全体的に抽象的ですが…)。

  • 「異常サンプルは低次元に圧縮しにくいだろう」という考えのもと、次元削減を試行するタイプ。
    • 次元削減したときにデコードが上手くいかないようなデータは異常なはずだ、ということです。このアプローチでは元々 PCA や kernel PCA やその他の PCA の亜種が用いられ、近年では深層自己符号化器を用いた手法が提案されています。ただ、サンプルの異常性を専ら「再構築エラー」のみから判断しているために限界がある、と書かれていますね。普通に正常サンプル並みの再構築エラーになる異常サンプルもあるだろうと。特にデータが全体的に再構築しにくい場合や、逆に次元削減手法の表現力があらゆるサンプルをカバーする場合はそうなるだろうと言っています。
  • (特徴空間に写像した後に)クラスタリングしてから孤立サンプル/少数クラスタを判別するタイプ。
    • そのままですが、このアプローチでもクラスタリングをする以上、高次元データに対しては次元削減必須です。この手法では次元削減とクラスタリングの学習が分離しているために、クラスタリングに必要な情報が次元削減の段階で失われうると指摘されています。これに対処するために、深層自己符号化器を利用して次元削減とクラスタリングを一貫で学習する手法も提案されていますが、クラスタリング時のモデルが単純すぎたり、自己符号化器の事前学習がクラスタリング側からの modify を妨げていると。
      • 3タイプの分類とは別枠で「次元削減と混合分布の学習の結合にも関心が高まっていた」という記述があるのですが(別枠なのはこれは必ずしも異常検知でないからですね)、話としてはこの2タイプ目に該当するのではないかと思うのでここに書いておきますね。といっても、この種の先行手法は次元削減が線形であったり、次元削減が事前学習を前提としていて柔軟性がなかったり、という主張が繰り返されているだけですね。それらに比べて本手法 DAGMM は「次元削減が事前学習不要で」「かつ非線形で」「かつGMMと結合して学習する」というのがよいのだと。それに加えて、特に異常検知のために再構築エラーも活用していると。
  • (特徴空間に写像した後に)一定割合の外れ値を覗いた正常データを囲む境界面を求めようとするタイプ(=One Class SVM)。
    • これもそのままです。特徴空間で正常データと異常データを分離する超平面を探そうということですね。このアプローチは高次元だと次元の呪いのために局所最適解に陥りやすいと書かれています。あ、Qiita で以下のような記事を見つけました。
      異常検知のための One Class SVM - Qiita
今回の提案手法 DAGMM は2番目のアプローチに1番目のアプローチの要素も取り入れたものといった位置付けなのですかね。そんなこんなで、今回の検証で DAGMM の比較対象にする先行手法は以下です。
OC-SVM
(Chen2001)
これは上の3番目の One Class SVM です。カーネル関数は通常の RBF カーネルを採用したとのことです。
DSEBM-e
(Zhai2016)
これらは上の1番目のアプローチで、ナブラを取ると深層自己符号化器の再構築エラーとなるようなエネルギー関数を定義するようですね。それで、DSEBM-e ではエネルギーが大きいサンプルを異常とみなし、DSEBM-r では再構築エラーが大きいサンプルを異常とみなします。
DSEBM-r
(Zhai2016)
DCN
(Yang2017)
異常検知に限らない次元削減と混合分布の結合学習のところで出てきた先行手法ですが、3タイプの分類だと2番目に該当するでしょうね。これは自己符号化器(深層ではない)の学習を後段の k-means が制御します。今回の異常検知の用途では、所属するクラスタクラスタ中心から距離が離れたデータを異常とみなします。

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

色々なアプローチにつて最新の手法と比較したってことなのかな。One Class SVM はちょっと古いけど、もう最近の研究はないのかな。SVM は異常検知に利用するには柔軟性に欠ける?

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

あと、上の比較対象の先行手法とは別に、DAGMM のフレームワークの各要素の有効性の検証のため、以下の DAGMM の亜種でも異常検知を行って DAGMM の結果と比較してみます。

GMM-EN学習時のみ目的関数から「再構築エラー」部分を取り除いたものです。学習後の異常検知は元の DAGMM で行います。
PAE目的関数から「エネルギー」ごと取り除いたものです。なのでもはやただの深層自己符号化器です。再構築エラーによって異常検知します。また、この場合は深層自己符号化器に事前学習(Vincent2010)を施すようです。
E2E-AEモデルの設定は上の PAE と同じそうなのですが、end-to-end で学習するようです。ちょっとよくわからなかったのですが、Zhai2016 と同様のエネルギー関数を利用して尤度を最大化しようとするということなのでしょうか。再構築エラーによって異常検知するということです。
PAE-GMM-EMまず事前学習済の深層自己符号化器を学習してから、EMアルゴリズムによってGMMのパラメータを推定します。学習後の異常検知は元の DAGMM と同じエネルギー関数で行います。
PAE-GMM上の PAE-GMM-EM とほとんど同じですが、EMアルゴリズムではなく推定ネットワークによってGMMのパラメータを推定します。
DAGMM-p上の PAE-GMM と似ているのですが、まず圧縮ネットワークを学習した後、次にモデル全体を結合学習します。
DAGMM-NVIDAGMM のエネルギー関数を、neural variational inference のそれに置き換えたものです。

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

亜種多いな!

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

なんかもう楽しそうですね。それで、再構築エラーには relative Euclidean distance  || x - x' ||_2 / ||x||_2 及びコサイン類似度  x \cdot  x' / (||x||_2 ||x'||_2) を採用したと書いてあるんですが、これは  f(x, x') \in \mathbb{R}^2 ってことなんですよね? あ、これらの距離関数を選定した理由は Appendix を読むようにと。

(次回があれば)つづく