今回はこのestimateGaussian.mとselectThreshold.mに取り組みたいと思います。
estimateGaussian.m
問題文を読む
まずこの2次元で簡単版をこなして、次に大きなデータでやってみようとのことです。
2次元Dataは307個なので、Xは307X2のベクトルですね。
グラフにするとこういうデータということです。
これをestimateGaussian.mを完成させるとこういうグラフになり、
そうすることができれば正解ということです。
estimateGaussian.mに含まれそうな式がこれ↓ですね。
ｐの方はかなり難しそうです。
プログラム全体版ex8.mを読む
estimateGaussian.mとmultivariateGaussian.mがカギのようです。
estamateGaussian.mは
[mu sigma2] = estimateGaussian(X)
なので、Xを入れるとmuとsigma2の答えを返す式ですね。
multivariateGaussian.mは
p = multivariateGaussian(X, mu, sigma2)
なので、X,mu,sigma2を入れるとpを返す式ですね。
こちらが難しそうな式ですが、幸いなことに今回の課題にはなっていません。
少し見てみましたが、いきなりkがでてきたりdatという新しい関数が出てきましたので
よくわかりません。。。課題じゃなくてよかった笑
プログラム詳細版estimateGaussian.mを読む
答えのmuとsigma2の答えは次の形式で指定されている点が
ポイントだと思いました。
mu = zeros(n, 1);
sigma2 = zeros(n, 1);
今回はn=2なので、
muは
0
0
で最初規定されていることになります。
試行錯誤する
それでは試行錯誤開始です。
Xは307X2のベクトルなので縦行ごとにμ（平均）とsigma2を出すのが
目的になります。
絵で表すとこういうことだと思います↓
Week2のFeatureScaleで
平均を計算するのにはmeanという関数があることを学んだので、
それを使って計算しました↓
答えは2X1なので、これを転置すれば正解ですね。
sigma2を計算するのは
A 307X2の614の要素すべてを平均からひいて
B それを2乗して、
C 1から307まで足す、つまり全要素足して
D ｍ(=307)で割るということなので、
H=mean(X)
A=X-H
B=A.^2
C=sum(B)
D=(1/m)*C
とすればいいと思います。
計算結果が↓となりました。
答えは2X1なので、これを転置すれば正解ですね。
Submitしたら見事正解となりました！
ところで図形はどうなっているか？並べて見てみたいと思います。
違う図形な気がします。円の数が違いますね。
なぜ正解になったんだろう？あとなぜ円の一番内側が黄色なのでしょうか？
と思いますが正解だったんで一度前に進むことにします。
selectThrehold.m
問題文を読む
Goalはタイトルの通り、Threholdであるεを選ぶことです。
Dataはcross validation setを使う。
cross validation setはラベリングされているDataでy=0はNormal,y=1はanomalyを示す。
どうやってεが大きいとか小さいとか決めるんだろうと思ってましたが、
こういう↓ことのようです。
you will implement an algorithm to select
the threshold " using the F1 score on a cross validation set.
εを決めるのにF1 scoreを用いなさいということです。
ですがF1 scoreってなんだっけな？と思いましたが、
Week6でやった内容でした↓
PrecisionとRecallはTrade-offの関係にあって、片方が上がるともう片方が下がるものでじゃあどうやって一番いいバランスを見つければいいのか、という時に役立つのが
F1Scoreでした。
pの予測が1であるか0であるかはε次第なので、F1Scoreが一番大きくなるεを
探しなさいということだと理解しました。
しかしどうやってやるのだろう？
良い関数あったかな。
ループで取り出すのは全部の値はできないから適切でない気がしますと
思いましたがループで計算するようです。
F1Scoreの計算の仕方が書いてありました。
tp; yval=1 かつpval=1
fp;yval=0かつpval=1
fn;yval=1かつpval=0
ということですね。
個数を数える式はpos=find(yval==1)とpos=find(pval==!)をうまく使う感じかなと
思っていたら、
sum(v == 0).とすればｖ＝0の個数が出る、
fp = sum*1とすればpval=1かつyval=1の
個数が出るという素晴らしいアドバイスコメントがありました！
εは8.99e-0.5(=0.0000899)が正解と問題文に書いてあるので、
これを目指します！
プログラム全体版(ex8.m)を読む
問題文以上のことは書いてありませんでした。
プログラム詳細版(selectThrehold.m)を読む
[bestEpsilon bestF1] = selectThreshold(yval, pval)
yvalとpvalがあれば、EpsilonとF1を返してくれる関数ですね。
yvalとpvalでF1を計算できるので、一番大きいF1を実現させるεを
探し出すということです。
ただそれをloopで探すとなると、
全部の値を検証できないんじゃないかと思っていましたが、
この式↓でかなり値をカバーするからよいという考え方をしていると思います。
stepsize = (max(pval) - min(pval)) / 1000;
for epsilon = min(pval):stepsize:max(pval)
end
stepsizeは0.000089で奇しくもGoalと同じです。
max(pval)は0.089
min(pval)はほぼ0
なので0.089を1000で割った値でloopを刻んでいくということです。
もう一つ便利な式が登場します。
cvpredictions = (pval < epsilon)
これはpvalがepsilonより小さいものは1として表記してくれます。
つまりAnomalyとして表記してくれるということです。
実験してみました。
epsilon=max(pval)としてみます。
そうするとpvalはすべてepsilonより小さくなるので
cvpredictionsはすべて1になるはずです。
なりました！
試行錯誤する
まずはXvalとYvalを見てみます。
で、yvalに含まれる1は
sum(yval==1)=9なので9個です。
便利な式を使って
epsilon=max(pval)
cvpredictions = (pval < epsilon)
cvpredictionsの307X1の要素をすべて1にします。
そーするとこの式↓の結果が9になるはずです。
sum*2
なりました！
次にF1Scoreの出し方を考えます。
F1=2*(prec * rec)/(prec + rec)
prec=tp/(tp+fp)
rec=tp/(tp+fn)
これでいいと思います！
Submitしたら正解しました！
しかしグラフがうまくいきません。
違うのは2つです
①円の数
②赤丸の数
①はさっきもあった問題ですね。。。
まずは②から取り組んでいきます。
こういう式↓からグラフを作ってたので、式の怪しいところを見ていきます。
epsilonが怪しいです。
bestEpsilonに紐づいて異常かどうか見てほしいのですが、
これだとloopの最後のepsilonになってそうです。
loopの最後のepsilonは0.089のはずです。
というのは、
-----------------------------------------------------------------------------------------
stepsizeは0.000089で奇しくもGoalと同じです。
max(pval)は0.089
min(pval)はほぼ0
なので0.089を1000で割った値でloopを刻んでいくということです。
-------------------------------------------------------------------------------------------
となっているので、
0--＞0.00089（10番目）---＞0.0089（100番目）--＞0.089（1000番目）
とloopさせているはずだからです。
計算結果を細かく見てみると、やはりそうです。
espilonは0.089でした。
ふっと思いましたが、bestEpsilonは0.0000899だったわけですが、
それはloopの1番目の値だったってことですね。
↑の式の
epsilon -->bestEpsilon
と変更すればできるはずです！
できました！
ちなみに
outliers = find(p < bestEpsilon)
はpがbestEpsilonより低い列数を返してくれる関数でした。
続いて①の円の数に取り組みたいと思います。
visualizefit.mに関係しているはずですがこれ以上はわからないので、
一旦諦めます。
課題ではありませんが2次元ではなく11次元のデータが用意されていました。
Xは1000X11なので11個の特徴を持つ1000のサンプルです。
Anomaly DetectionなのでYはありません。
ただしvalidation SetはXval；100X11、yval;100X1で
yvalの中に1が10個入っています。
2次元セットでやったようにXvalを使って
100X1のpvalを計算します。
（pval = multivariateGaussian(Xval, mu, sigma2)）
あとはpval(max)とpval(min)の間を1000刻みで
εを動かしていってF1が最大になるεを探します。
こうしてValidation Setで見つけた最適εを
X；1000X11を基にしたp；1000X1に適用する。
1000個中117個のｐがε以下だったので、
11.7％が異常ということでした。
今回はここまでです。また次回頑張ります！