머신러닝 옥타브 실습(8-2):이상 탐지 알고리즘 구현

   온라인 강의 플랫폼 코세라의 창립자인 앤드류 응 (Andrew Ng) 교수는 인공지능 업계의 거장입니다. 그가 스탠퍼드 대학에서 머신 러닝 입문자에게 한 강의를 그대로 코세라 온라인 강의 (Coursera.org)에서 무료로 배울 수 있습니다. 이 강의는 머신러닝 입문자들의 필수코스입니다. 인공지능과 머신러닝을 혼자 공부하면서 자연스럽게 만나게 되는 강의입니다. 

Programming Exercise 8: 

Anomaly Detection and Recommender Systems (이상 탐지와 추천 시스템)

1. Anomaly detection(이상 탐지)

1.3 Selecting the threshold, ε (임계치 결정하기)

   Now that you have estimated the Gaussian parameters, you can investigate which examples have a very high probability given this distribution and which examples have a very low probability. The low probability examples are more likely to be the anomalies in our dataset. One way to determine which examples are anomalies is to select a threshold based on a cross validation set. In this part of the exercise, you will implement an algorithm to select the threshold ε using the F1 score on a cross validation set.

   You should now complete the code in selectThreshold.m. For this, we will use a cross validation set {(x(1), y(1)), . . . , (x(mcv), y(mcv))}, where the label y = 1 corresponds to an anomalous example, and y = 0 corresponds to a normal example. For each cross validation example, we will compute p(x(i)). The vector of all of these probabilities p(xcv(1)), . . . , p(x(mcv) ) is passed to selectThreshold.m in the vector pval. The corresponding labels ycv(1), . . . , ycv(mcv) is passed to the same function in the vector yval.

    The function selectThreshold.m should return two values; the first is the selected threshold ε. If an example x has a low probability p(x) < ε, then it is considered to be an anomaly. The function should also return the F1 score, which tells you how well you’re doing on finding the ground truth anomalies given a certain threshold. For many different values of ε, you will compute the resulting F1 score by computing how many examples the current threshold classifies correctly and incorrectly.

   지금까지 가우시안 분포의 파라미터를 계산하였으므로 확률이 매우 낮거나 매우 높은 예제를 찾을 수 있습니다. 정상 범위를 벗어난 예제는 이상일 가능성이 높습니다. 어떤 예제가 이상이 있는지 확인하는 방법은 교차 검증 셋에서 임계치를 선택합니다. 실습에서 교차 검증 셋의 F1 스코어를 사용하여 임계값 ε을 선택하는 알고리즘을 구현합니다.

   selectThreshold.m 파일을 작성합니다. 교차 검증 셋은 다음과 같습니다.

   여기서 레이블 y=1은 이상 에제이고, y = 0은 정상 예제입니다. 각 교차 검증 셋의 예제에 대한 추정 확률은 다음과 같습니다.

   교차 검증 셋의 추정 확률 pval과  교차 검증 셋에 대응하는 레이블  yval을 selectThreshold.m 파일에 전달합니다.

   selectTehreshold.m 파일은 두 가지 변수를 반환합니다. 하나는 임계값 ε입니다. 학습 예제 x가 p(x) < ε 라면 이상 예제로 간주합니다.  또 하나는 F1 score입니다. F1 score는 특정 임계값에서 실제로 이상 예제를 얼마나 잘 찾는 지를 알려줍니다. 여러 ε 값에 대해 현재 임계값이 잘못 분류하는 예제 수를 계산합니다. 

  The F1 score is computed using precision (prec) and recall (rec):

   F1 스코어는 정확도(prec)와 재현율(rec)을 사용하여 계산합니다. 

 

•. tp is the number of true positives: the ground truth label says it’s an anomaly and our algorithm correctly classified it as an anomaly.

• fp is the number of false positives: the ground truth label says it’s not an anomaly, but our algorithm incorrectly classified it as an anomaly.

• fn is the number of false negatives: the ground truth label says it’s an anomaly, but our algorithm incorrectly classified it as not being anoma- lous.

   tp (True positives의 수) :  yval에서 이상 예제를 pval에서도 이상 예제로 예측한 경우입니다.  

   fp(False Positive의 수) : yval에서 정상 예제를 pval에서는 이상 예제로 예측한 경우입니다.  

   fn( False Negative의 수) : yval에서 이상 예제를 pval에서는 정상 예제로 예측한 경우입니다.  

   

   In the provided code selectThreshold.m, there is already a loop that will try many different values of ε and select the best ε based on the F1 score. You should now complete the code in selectThreshold.m. You can im- plement the computation of the F1 score using a for-loop over all the cross validation examples (to compute the values tp, fp, fn). You should see a value for epsilon of about 8.99e-05.

   selectThreshold.m 파일에서 이미 다양한 ε값을 시도하고 F1 score에 따라 최상의 ε을 선택하는 루프를 제공합니다. selectThreshold.m 파일의 코드를 완료합니다. 모든 교차 검증 예제에 대한 for 루프를 사용하여 F1 score를 계산을 구현합니다. 각각의 tp, fp, fn을 계산할 수 있습니다. epsilon은 약 8.99e-0.5의  값을 표시합니다. 

   Implementation Note: 

   In order to compute tp, fp and fn, you may be able to use a vectorized implementation rather than loop over all the examples. This can be implemented by Octave/MATLAB’s equality test between a vector and a single number. If you have several binary values in an n-dimensional binary vector v ∈ {0,1}n, you can find out how many values in this vector are 0 by using: sum(v == 0). You can also apply a logical and operator to such binary vectors. For instance, let cvPredictions be a binary vector of the size of your number of cross validation set, where the i-th element is 1 if your algorithm considers xcv^(i) an anomaly, and 0 otherwise. You can then, for example, compute the number of false positives using: fp = sum((cvPredictions == 1) & (yval == 0)).

   구현 노트:

   tp, fp, fn을 계산하기 위해 모든 예제에 대한 루프보다 벡터화 구현을 활용합니다. 이것은 옥타브/ 매트랩의 벡터와 단일 숫자 간의 동등성 테스트로 구현합니다. n차원 이진 벡터 v가 0 또는 1의 값이 여러 개 있는 경우 논리 연산자 sum( v == 0) 코드로 벡터의 값이 0인지 확인합니다.  예를 들어, cvPredictions를 교차 거증 셋의 이진 벡터로 설정합니다. 여기서 알고리즘이 xcv^(i)를 이상 예제로 간주하면 i번째 성분이 1이고 정상 예제로 간주하면 0입니다. 예를 들어,

 fp = sum((cvPredictions == 1) & (yval == 0))

   위의 코드를 사용하면 간단하게 실제 레이블은 y = 0이지만 예측은 1인 fp (False Negative)의 수를 쉽게 계산할 수 있습니다. 

 

   Once you have completed the code in selectThreshold.m, the next step in ex8.m will run your anomaly detection code and circle the anomalies in the plot (Figure 3).

   You should now submit your solutions.

   selectThreshold.m 파일의 코드를 완료합니다. ex8.m 파일은 이상 탐지 코드를 실행하고 그림 3과 같이 이상 에제에 빨간색 동그라미로 표시합니다.

   완료 후 submit 합니다.

<Part 3 : 아웃라이어 찾기(상) : 앱실론 계산히기>

(1) 교차 검증 셋을 다변량 가우시안 분포로 확률 추정  

   교차 검증 셋 Xval을 가우시안 분포로 확률을 추정하는 multivariateGaussian.m 파일을 작성했습니다. 이 파일을 호출할 때 데이터 행렬 Xval, 피처 별 평균 mu, 피처 별 분산 sigma2를 전달합니다. 

pval = multivariateGaussian(Xval, mu, sigma2);

   multivariateGaussian.m 파일은 다변량 가우시안 분포에 따른 추정 확률 pval을 반환합니다. 

   데이터셋 X로 계산한 피처 별 평균 mu와 분산 sigma2의 값을 확인합니다. 

>> mu

mu =

   14.112

   14.998

>> sigma2

sigma2 =

   1.8326

   1.7097

   multivariateGaussian.m 파일을 호출하여 pval을 계산합니다. 

>> pval = multivariateGaussian(Xval, mu, sigma2);

>> size(pval)

ans =

   307     1

>> pval

pval =

   4.1632e-02

   8.1909e-02

   4.0716e-02

   6.1900e-02

   7.1187e-02

   ...

(2) 임계값 epsilon 계산하기

   추정 확률에서 이상 예제 여부를 판단하는 임계값을 계산하는 selectThreshold.m 파일을 작성합니다. 이 파일을 호출할 때 교차 검증 셋의 레이블 yval과 교차 검증 셋의 추정 확률 pval을 전달합니다.  

[epsilon F1] = selectThreshold(yval, pval);

   selectThreshold.m 파일은 임계값 epsilon과 F1 score를 반환합니다. 두 변수의 값은 여러 가지 다른 값들과 비교하여 가장 최적의 값입니다  

(3) selectThreshold.m 파일 분석

 

function [bestEpsilon bestF1] = selectThreshold(yval, pval)

%SELECTTHRESHOLD 아웃라이어를 찾기 위한 최적의 임계값 epsilon을 계산

%       아웃라이어(outliers)는 정상의 범위를 크게 벗어난 값을 의미

%   [bestEpsilon bestF1] = SELECTTHRESHOLD(yval, pval) 

%        교차 검증 셋의 추정 확률 pval과 실제 레이블 yval에 기반한  최고의 임계값 계산

%

% 반환할 변수 초기화

bestEpsilon = 0;

bestF1 = 0;

% 변수 초기화

F1 = 0;

stepsize = (max(pval) - min(pval)) / 1000;  % epsilon의 증가량

% 최고의 epsilon을 찾기 위한 루프 

for epsilon = min(pval):stepsize:max(pval)    

    % ====================== YOUR CODE HERE ======================

    % Instructions: 

    %        F1 값을 지정하고 임계치로 apsilon을 선택하기 위한 F1 스코어를 계산

    %        루프의 마지막 코드는 현재 계산한 F1 스코어와 기존에 계산한 최고의 값과 비교

    %         최적의 F1 스코어와 epsilon을 남김 

    %               

    % Note: 

     %       아웃라이어 예측 0 또는 1의 이진 벡터를 얻기 위해 아래 코드 활용

     %           predictions = (pval < epsilon)

    % =============================================================

  % F1 스코어가 기존의 최고값 bestF1보다 좋으면 값을 대체    

    if F1 > bestF1             

       bestF1 = F1;

       bestEpsilon = epsilon;

    end

end

end

(4) F1 score  계산하기 : tp, fp, fn 계산하기

  필요한 변수를 초기화합니다.

bestEpsilon = 0;

bestF1 = 0;

F1 = 0;

%%% 최고의 epsilon을 찾는 for 루프 만들기

   for 루프는 epsilon 값을 pval의 최소값에서 최대값까지 stepsize 간격으로 증가시킵니다. 

>> min(pval)

ans =    4.5133e-36

>> max(pval)

ans =  0.089909

>> stepsize = (max(pval) - min(pval)) / 1000

stepsize =    8.9909e-05

   for 루프는 다음과 같이 만듭니다.

stepsize = (max(pval) - min(pval)) / 1000;

for epsilon = min(pval):stepsize:max(pval)    

end

%%% predictions의 값 확인

   추정 확률 pval이 임계값 epsilon 보다 작으면 이상 예제입니다. predicion는 추정 확률 pval이 epsilon보다 작으면  참(1)이고 이상 예제입니다. 그리고 추정 확률 pval이 epsilon 보다 크면 거짓(0)이고 정상 예제입니다.  따라서, selectThreshold.m 파일에서 예측치 predictions을 다음의 코드로 제시합니다. 

 predictions = pval < epsilon;

%%% True positives의 수 tp 계산하기 

   관계 연산자' =='은 벡터의 성분과 같은 값이면 1의 값을 반환합니다. 

>> v

v =

   1   2   3   4   5

>> v == 1

ans =

   1   0   0   0   0

   논리 연산자를 활용하여 F1 score를 계산하기 위한 tp, fp, tn을 계산합니다. 

   tp(True positives의 수)는 yval에서 이상 예제를 pval에서도 이상 예제로 예측한 경우입니다. 따라서, preditions의 값도 1이고, yval의 값도 1인 예제를 찾는 코드는 다음과 같습니다.

predictions == 1 & yval ==1

   predictions와 yval은 모두 307 X 1차원 행렬입니다. 모두 참일 경우에만 1이므로 결과값을 모두 합하면 True positives의 수와 같습니다.

tp = sum(predictions == 1 & yval ==1);

    

%%% False positives의 수 fp 계산하기 

   fp(False Positive의 수)는 yval에서 정상 예제를 pval에서는 이상 예제로 예측한 경우입니다.  따라서, predictions의 값은 1이지만 yval의 값은 0인 예제를 찾는 코드는 다음과 같습니다.

predictions == 1 & yval ==0

   predictions와 yval은 모두 307 X 1차원 행렬입니다. 모두 참일 경우에만 1이므로 결과값을 모두 합하면 False positives의 수와 같습니다.

fp = sum(predictions == 1 & yval ==0);

%%% False Negative의 수 fn 계산하기 

   fn(False Negative의 수)는 yval에서 이상 예제를 pval에서는 정상 예제로 예측한 경우입니다.  따라서, predictions의 값은 0이지만 yval의 값은 1인 예제를 찾는 코드는 다음과 같습니다.

predictions == 0 & yval ==1

   predictions와 yval은 모두 307 X 1차원 행렬입니다. 모두 참일 경우에만 1이므로 결과값을 모두 합하면 False Negatiges의 수와 같습니다.

fn = sum(predictions == 0 & yval ==1);

%%%% for 루프에서 제대로 연산하는 지를 확인

stepsize = (max(pval) - min(pval)) / 1000;

for epsilon = min(pval):stepsize:max(pval)    

    predictions = pval < epsilon;

   tp = sum(predictions == 1 & yval ==1);

   fp = sum(predictions == 1 & yval ==0);

   fn = sum(predictions == 0 & yval ==1);

end

결과는 다음과 같습니다.

>> epsilon

epsilon =  0.089909

>> tp

tp =  9

>> fp

fp =  297

>> fn

fn = 0

(5) F1 score  계산하기 : 정확도(prec)와 재현율 (recall) 계산하기

정확도를 계산하는 공식에 따라 코드는 다음과 같습니다.

prec = tp/(tp +fp);

재현율을 계산하는 공식에 따라 코드는 다음과 같습니다.

rec = tp/(tp+fn);

(5) F1 score  계산하기 

   F1 Score를 계산하는 공식은 다음과 같습니다.

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

(6) 정답

function [bestEpsilon bestF1] = selectThreshold(yval, pval)

%SELECTTHRESHOLD 아웃라이어를 찾기 위한 최적의 임계값 epsilon을 계산

%       아웃라이어(outliers)는 정상의 범위를 크게 벗어난 값을 의미

%   [bestEpsilon bestF1] = SELECTTHRESHOLD(yval, pval) 

%        교차 검증 셋의 추정 확률 pval과 실제 레이블 yval에 기반한  최고의 임계값 계산

%

% 반환할 변수 초기화

bestEpsilon = 0;

bestF1 = 0;

% 변수 초기화

F1 = 0;

stepsize = (max(pval) - min(pval)) / 1000;  % epsilon의 증가량

% 최고의 epsilon을 찾기 위한 루프 

for epsilon = min(pval):stepsize:max(pval)    

    % ====================== YOUR CODE HERE ======================

    % Instructions: 

    %        F1 값을 지정하고 임계치로 apsilon을 선택하기 위한 F1 스코어를 계산

    %        루프의 마지막 코드는 현재 계산한 F1 스코어와 기존에 계산한 최고의 값과 비교

    %         최적의 F1 스코어와 epsilon을 남김 

    %               

    % Note: 

     %       아웃라이어 예측 0 또는 1의 이진 벡터를 얻기 위해 아래 코드 활용

     %           predictions = (pval < epsilon)

predictions = pval < epsilon;

tp = sum(predictions == 1 && yval ==1);

fp = sum(predictions == 1 && yval ==0);

fn = sum(predictions == 0 && yval ==1);

prec = tp/(tp +fp);

rec = tp/(tp+fn);

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

    % =============================================================

  % F1 스코어가 기존의 최고값 bestF1보다 좋으면 값을 대체    

    if F1 > bestF1             

       bestF1 = F1;

       bestEpsilon = epsilon;

    end

end

(7) 결과 확인

   옥타브 프로그램에 다음 코드를 입력합니다.

clear ; close all; clc

load('ex8data1.mat');

[mu sigma2] = estimateGaussian(X);

p = multivariateGaussian(X, mu, sigma2);

pval = multivariateGaussian(Xval, mu, sigma2);

[epsilon F1] = selectThreshold(yval, pval);

   selectThreshold.m 파일이 반환한 epsilon과 F1의 값을 확인합니다. 

>> epsilon

epsilon =    8.9909e-05

>> F1

F1 =  0.87500

<Part 3 : 아웃라이어 찾기(하) : 아웃라이어 표시하기>

(1) 아웃라이어 찾기

   selectThreshold.m 파일에서 교차 검증 셋 xval과 교차 검증 셋 레이블 yval에서 임계값 epsilon을 계산했습니다. 아웃라이어(outliers)는 임계값보다 작은 값입니다. 학습 셋 X의 추정 확률 p가 epsilon 보다 작은 값을 찾습니다. 행렬 성분에서 특정 값을 찾을 때 find() 함수를 사용합니다. 

   find () 함수는 조건에 해당하는 행렬 성분의 인덱스를 반환합니다. 

>> outliers = find(p < epsilon)

outliers =

   301

   302

   304

   305

   306

   307

(2) 기존 데이터에서 아웃라이어 표시하기

   데이터셋과 다변량 가우시안 분포에 따른 확률 밀도 함수에 대한 등고선을 표시하는 visualizeFit.m 파일을 작성했습니다. 이 파일을 호출할 때 데이터셋 X, 평균 mu, 분산 sigma2를 전달하면 그래프를 그립니다.  

visualizeFit(X,  mu, sigma2);

xlabel('Latency (ms)');

ylabel('Throughput (mb/s)');

axis([0 30 0 30]);

 

   아웃라이어를 찾는 코드를 입력합니다.

outliers = find(p < epsilon);

   기존 그림에 아웃라이어를 표시합니다.

hold on

plot(X(outliers, 1), X(outliers,2), 'ro', 'LineWidth', 2,'MarkerSize', 10);

hold off

1.4 High dimensional dataset (고차원 데이터 셋)

   The last part of the script ex8.m will run the anomaly detection algorithm you implemented on a more realistic and much harder dataset. In this dataset, each example is described by 11 features, capturing many more properties of your compute servers.

   The script will use your code to estimate the Gaussian parameters (μi and σi2), evaluate the probabilities for both the training data X from which you estimated the Gaussian parameters, and do so for the the cross-validation set Xval. Finally, it will use selectThreshold to find the best threshold ε. You should see a value epsilon of about 1.38e-18, and 117 anomalies found.

   ex8.m 파일의 마지막 실습은 더 현실적이고 더 어려운 데이터 셋에서 이상 탐지 알고리즘을 실행합니다. 데이터셋의 각 학습 예제는 11 가지 피처로 구성되어 서버의 더 많은 특성을 포착합니다. 

   ex9.m은 코드를 사용하여 가우스 파라미터 μi and σi^2를 추정하고, 학습 데이터 X에 대한 확률을 추정합니다. 교차 검증 셋 Xval에서 확률을 추정합니다. 마지막으로  selectThreshold를 사용하여 최적의 임계값 ε을 찾습니다. epsilon은 약 1.28e-18이고 117개의 이상을 발견합니다. 

<Part 4 : 다차원 아웃라이어>

(1) 데이터 로드 

clear ;             % 옥타브 프로그램에 모든 변수를 제거

close all;         % 터미널 이외의 창을 닫음

clc                   % 터미널을 깨끗이 정리 

load ('ex8data2.mat');

  데이터를 옥타브 프로그램에 로드합니다.

>> load ('ex8data2.mat');

>> whos

Variables in the current scope:

   Attr Name        Size                     Bytes  Class

   ==== ====        ====                     =====  =====

        X               1000x11                     88000  double

        Xval            100x11                      8800  double

        yval            100x1                        800  double

Total is 12200 elements using 97600 bytes

(2) 피처 별 평균과 분산 계산   

   데이터 행렬 X에서 피처 별로 평균과 분산을 계산하는 estimateGaussian.m 파일을 작성했습니다. 파일을 호출할 때  데이터 행렬 X를 전달합니다. 

[mu sigma2] = estimateGaussian(X);

   estimateGaussian.m 파일은 피처 별 평균 벡터 mu와 피처 별 분산 벡터 sigma2를 반환합니다. 벡터 mu와 벡터 sigma2는 11 X 1차원 벡터입니다. 

>> [mu sigma2] = estimateGaussian(X)

mu =

    4.9394

   -9.6373

   13.8147

  -10.4645

   -7.9562

   10.1995

   -6.0194

    7.9698

   -6.2532

    2.3245

    8.4737

sigma2 =

   60.975

   53.206

   58.515

   84.204

   65.269

   89.575

   55.633

   87.162

   29.629

   70.785

   50.504

(3) 다변량 가우시안 분포로 확률 추정 

   데이터셋 X을 다변량 가우시안 분포로 확률을 추정하는 multivariateGaussian.m 파일을 작성했습니다. 이 파일을 호출할 때 데이터 행렬 X, 피처 별 평균 mu, 피처 별 분산 sigma2를 전달하면, 다변량 가우시안 분포에 따른 추정 확률 p를 반환합니다. 

p = multivariateGaussian(X, mu, sigma2);

   이 파일을 호출할 때 교차 검증 셋 Xval, 피처 별 평균 mu, 피처 별 분산 sigma2를 전달하면,  다변량 가우시안 분포에 따른 추정 확률 pval을 반환합니다. 

pval = multivariateGaussian(Xval, mu, sigma2);

   

(4) 임계값 epsilon 계산하기

   추정 확률에서 이상 예제 여부를 판단하는 임계값을 계산하는 selectThreshold.m 파일을 작성했습니다. 이 파일을 호출할 때 교차 검증 셋의 레이블 yval과 교차 검증 셋의 추정 확률 pval을 전달합니다.  

[epsilon F1] = selectThreshold(yval, pval);

   selectThreshold.m 파일은 임계값 epsilon과 F1 score를 반환합니다. 두 변수의 값은 여러 가지 다른 값들과 비교하여 가장 최적의 값입니다  

   

>> [epsilon F1] = selectThreshold(yval, pval)

epsilon =    1.3772e-18

F1 =  0.61538

(5) 아웃라이어 찾기

   아웃라이어(outliers)는 임계값보다 작은 값입니다. 학습 셋 X의 추정 확률 p가 epsilon 보다 작은 값을 찾습니다. 행렬 성분에서 특정 값을 찾을 때 find() 함수를 사용합니다. find() 함수는 조건에 맞는 행렬 성분의 인덱스를 반환합니다. 

>> outliers = find(p < epsilon)

outliers =

    10

    21

    22

    31

    ...

>> size(outliers)

ans =

   117     1