暇人日記

アラフォーおっさんのコーセラの機械学習の課題を解こうと頑張っています!

Coursera Machine Learning Week9 課題 苦闘日記①

今回はこのestimateGaussian.mとselectThreshold.mに取り組みたいと思います。

f:id:omoshiroamericanews:20191130094741p:plain

 

 estimateGaussian.m

 

問題文を読む

 

まずこの2次元で簡単版をこなして、次に大きなデータでやってみようとのことです。

2次元Dataは307個なので、Xは307X2のベクトルですね。

グラフにするとこういうデータということです。

f:id:omoshiroamericanews:20191130103818p:plain

 

 

これをestimateGaussian.mを完成させるとこういうグラフになり、

そうすることができれば正解ということです。

 

f:id:omoshiroamericanews:20191130103908p:plain

estimateGaussian.mに含まれそうな式がこれ↓ですね。

pの方はかなり難しそうです。

f:id:omoshiroamericanews:20191130103844p:plain

 プログラム全体版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を出すのが

目的になります。

絵で表すとこういうことだと思います↓

f:id:omoshiroamericanews:20191130104116p:plain

 

 

Week2のFeatureScaleで

平均を計算するのにはmeanという関数があることを学んだので、

それを使って計算しました↓

 

f:id:omoshiroamericanews:20191130104341p:plain

 

答えは2X1なので、これを転置すれば正解ですね。

 

sigma2を計算するのは

A 307X2の614の要素すべてを平均からひいて

B それを2乗して、

C 1から307まで足す、つまり全要素足して

D m(=307)で割るということなので、

 

H=mean(X)

A=X-H
B=A.^2
C=sum(B)
D=(1/m)*C

 

とすればいいと思います。

計算結果が↓となりました。

 

f:id:omoshiroamericanews:20191130104149p:plain

 

答えは2X1なので、これを転置すれば正解ですね。

Submitしたら見事正解となりました!

 

 

ところで図形はどうなっているか?並べて見てみたいと思います。

f:id:omoshiroamericanews:20191130104427p:plain

 

違う図形な気がします。円の数が違いますね。

なぜ正解になったんだろう?あとなぜ円の一番内側が黄色なのでしょうか?

と思いますが正解だったんで一度前に進むことにします。

 

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でやった内容でした↓

 

f:id:omoshiroamericanews:20191130174324p:plain

 

PrecisionとRecallはTrade-offの関係にあって、片方が上がるともう片方が下がるものでじゃあどうやって一番いいバランスを見つければいいのか、という時に役立つのが

F1Scoreでした。

 

pの予測が1であるか0であるかはε次第なので、F1Scoreが一番大きくなるεを

探しなさいということだと理解しました。

しかしどうやってやるのだろう?

良い関数あったかな。

ループで取り出すのは全部の値はできないから適切でない気がしますと

思いましたがループで計算するようです。

 

f:id:omoshiroamericanews:20191130175131p:plain


 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).とすればv=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になるはずです。

f:id:omoshiroamericanews:20191130234310p:plain

なりました!

 

 

 

試行錯誤する

まずはXvalとYvalを見てみます。

f:id:omoshiroamericanews:20191130183802p:plain

 

 で、yvalに含まれる1は

sum(yval==1)=9なので9個です。

 

便利な式を使って

epsilon=max(pval)

cvpredictions = (pval < epsilon)

cvpredictionsの307X1の要素をすべて1にします。

 

そーするとこの式↓の結果が9になるはずです。

sum*2

 

なりました!

f:id:omoshiroamericanews:20191130235425p:plain

 

次にF1Scoreの出し方を考えます。

 

F1=2*(prec * rec)/(prec + rec)

 

prec=tp/(tp+fp)
rec=tp/(tp+fn)

 

tp=sum*3
fp=sum*4
fn=sum*5

 

これでいいと思います!

Submitしたら正解しました!

 

しかしグラフがうまくいきません。

f:id:omoshiroamericanews:20191206170227p:plain

 違うのは2つです

①円の数

②赤丸の数

 

①はさっきもあった問題ですね。。。

まずは②から取り組んでいきます。

 

こういう式↓からグラフを作ってたので、式の怪しいところを見ていきます。

f:id:omoshiroamericanews:20191206171211p:plain

 

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でした。

 

f:id:omoshiroamericanews:20191206171832p:plain

ふっと思いましたが、bestEpsilonは0.0000899だったわけですが、

それはloopの1番目の値だったってことですね。

 

↑の式の

epsilon -->bestEpsilon

と変更すればできるはずです!

 

f:id:omoshiroamericanews:20191206172427p:plain

できました!

ちなみに

outliers = find(p < bestEpsilon)

はpがbestEpsilonより低い列数を返してくれる関数でした。

f:id:omoshiroamericanews:20191206172530p:plain

続いて①の円の数に取り組みたいと思います。

 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個のpがε以下だったので、

11.7%が異常ということでした。 

 

今回はここまでです。また次回頑張ります!

 

 

*1:cvPredictions == 1) &(yval == 0

*2:cvpredictions==1)&(yval==1

*3:cvpredictions==1)&(yval==1

*4:cvpredictions==1)&(yval==0

*5:cvpredictions==0)&(yval==1