brunch

You can make anything
by writing

C.S.Lewis

by 라인하트 Feb 13. 2021

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

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


Programming Exercise 8: 

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


1. Anomaly detection(이상 탐지)


   In this exercise, you will implement an anomaly detection algorithm to detect anomalous behavior in server computers. The features measure the through- put (mb/s) and latency (ms) of response of each server. While your servers were operating, you collected m = 307 examples of how they were behaving, and thus have an unlabeled dataset {x(1),...,x(m)}. You suspect that the vast majority of these examples are “normal” (non-anomalous) examples of the servers operating normally, but there might also be some examples of servers acting anomalously within this dataset.

   You will use a Gaussian model to detect anomalous examples in your dataset. You will first start on a 2D dataset that will allow you to visualize what the algorithm is doing. On that dataset you will fit a Gaussian dis- tribution and then find values that have very low probability and hence can be considered anomalies. After that, you will apply the anomaly detection algorithm to a larger dataset with many dimensions. You will be using ex8.m for this part of the exercise.

   The first part of ex8.m will visualize the dataset as shown in Figure 1.


   이번 실습은 서버의 이상 동작을 탐지하는 알고리즘을 구현합니다. 각 서버의 응답 속도 (mb/s)와 지연 시간 (ms)를 측정합니다. 학습 예제 m = 307개를 수집하고 레이블이 지정되지 않은 데이터셋 x^(1), x^(2),..., x^(m)이 있습니다. 대부분의 학습 예제는 정상적으로 동작하는 정상 예제이지만 데이터셋 안에 이상 예제가 몇 개 있을 수도 있습니다. 

   가우시안 모델을 적용하여 데이터셋의 비정상적인 예제를 탐지합니다. 먼저 알고리즘의 작업을 시각화하기 위해 2차원 데이터 셋에서 시작합니다. 해당 데이터 셋에서 적합한 가우시안 분포를 계산하고 확률이 매우 낮아 이상 예제로 간주할 수 있는 확률을 찾습니다. 이상 탐지 알고리즘을 차원이 더 높은 큰 데이터셋에 적용합니다. 이번 실습은 ex8.m 파일을 사용합니다.

   ex8.m 파일은 그림 1과 같이 데이터 셋을 시각화합니다. 



<Part 1: 데이터셋 로드>


(1) 데이터 로드 


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

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

clc                   % 터미널을 깨끗이 정리 


load ('ex8data1.mat');

   

   ex8data1.mat 파일은 3개의 데이터셋을 옥타브 프로그램에 로드합니다. 데이터 행렬 X는 피처 x1은 서버의 지연시간 와 피처 x2는 서버의 응답속도에 대한 데이터를 307개를 보유합니다. 또한, 평가를 위한 교차 검증 셋 Xval과 교차 검증 셋의 레이블 yval이 있습니다. 


>> load ('ex8data1.mat');

>> whos

Variables in the current scope:


   Attr Name        Size                     Bytes  Class

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

        X                  307x2                       4912  double

        Xval             307x2                       4912  double

        yval             307x1                       2456  double


Total is 1535 elements using 12280 bytes


>> X

X =

   13.0468   14.7412

   13.4085   13.7633

   14.1959   15.8532

   14.9147   16.1743

   13.5767   14.0428

   ...


>> Xval

Xval =

   15.79026   14.92102

   13.63962   15.32996

   14.86590   16.47387

   13.58468   13.98931

   13.46404   15.63533

   ...


>> yval

yval =

   0

   0

   0

   0

   0

   ...


(2) 데이터 시각화


  데이터를 'X'표로 도식화합니다. 


plot (X(:,1), X(:,2),'x');


axis ([0 30 0 30]);


xlabel('Latency (ms)');

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





1.1 Gaussian distribution (가우시안 분포)


   To perform anomaly detection, you will first need to fit a model to the data’s distribution.

   Given a training set {x(1),...,x(m)} (where x(i) ∈ Rn), you want to estimate the Gaussian distribution for each of the features xi. For each feature i = 1...n, you need to find parameters μi and σi^2 that fit the data in the i-th dimension {x(1), ..., x(m)} (the i-th dimension of each example). 

   The Gaussian distribution is given by


   where μ is the mean and σ2 controls the variance.


   이상 탐지를 수행하기 위해 먼저 모델을 데이터 분포에 적합하게 조정합니다. 학습 셋 x^(10, x^(2),..., x^(m)은 R^(n) 차원이고 피처 xi에 대한 가우스 분포를 추정합니다. 각 피처 i= 1, 2,.., n에 대한 파라미터 μi와 σi^2를 계산합니다. 

    가우스 분포는 다음과 같고, μi는 i번째 피처에 대한 평균이고, σi^2는 i번째 피처에 대한 표준편차의 제곱이자 분산입니다. 



1.2 Estimating parameters for a Gaussian

      (가우시안 분포를 위한 파라미터 추정)


   You can estimate the parameters, (μi, σi^2), of the i-th feature by using the following equations. To estimate the mean and variance, you will use :


   Your task is to complete the code in estimateGaussian.m. This function takes as input the data matrix X and should output an n-dimension vector mu that holds the mean of all the n features and another n-dimension vector sigma2 that holds the variances of all the features. You can implement this using a for-loop over every feature and every training example (though a vectorized implementation might be more efficient; feel free to use a vectorized implementation if you prefer). Note that in Octave/MATLAB, the var function will (by default) use 1/(m-1) , instead of 1/m , when computing σ2. 

   Once you have completed the code in estimateGaussian.m, the next part of ex8.m will visualize the contours of the fitted Gaussian distribution. You should get a plot similar to Figure 2. From your plot, you can see that most of the examples are in the region with the highest probability, while the anomalous examples are in the regions with lower probabilities.

   You should now submit your solutions.


   다음 방정식은 i 번째 피처의 파라미터 평균 μi와 표준 편차의 제곱 σi^2를 추정합니다. 

   먼저 EstimatedGaussian.m 파일의 코드를 완성합니다. 이 파일은 데이터 행렬 X를 입력받아 모든 n개의 피처의 평균을 가지는 n 차원 벡터 mu와  n개의 피처의 분산(표준 편차의 제곱)을 가지는 n 차원 벡터 sigma2를 반환합니다. 모든 피처와 모든 학습 예제에 대해 for 루프를 사용할 수 있지만, 벡터화된 구현이 더 효율적입니다. 옥타브/매트랩 프로그램에서 var() 함수는 σ2를 계산할 때 1/m 이 아닌 1/(m-1)을 사용합니다. 

   EstimatesGaussian.m 파일을 완성한 후 ex8.m 파일은 가우시안 분포를 그림 2와 같이 시각화합니다.  대부분의 예제는 확률이 가장 높은 지역에 위치하고 비정상적인 예제는 확률이 낮은 지역에 있음을 확인합니다. 

   완료한 후에 submit 명령을 실행합니다. 



<Part 2: 데이터셋 통계 추정(상) : 평균과 분산을 계산하기>


(1) 데이터 로드 


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

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

clc                   % 터미널을 깨끗이 정리 


load ('ex8data1.mat');



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


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


[mu sigma2] = estimateGaussian(X);


   estimateGaussian.m 파일은 피처 별 평균 벡터 mu와 피처 별 분산 벡터 sigma2를 반환합니다.


(3) estimateGaussian.m 파일 분석


function [mu sigma2] = estimateGaussian(X)

%ESTIMATEGAUSSIAN  데이터 행렬 X의 피처 별 가우시안 분포의 파라미터를 계산

%   [mu sigma2] = estimateGaussian(X), 

%             입력 X :  데이터셋 X

%              mu : 피처 별 평균을 나타내는 n차원 벡터

%              sigma2 : 피처 별 분산을 나타내는 n차원 벡터 

% 변수 정의 

[m, n] = size(X);


% 반환할 변수를 초기화 

mu = zeros(n, 1);

sigma2 = zeros(n, 1);


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

% Instructions: 

%               데이터의 평균과 분산을 계산 In particular, mu(i) should contain the mean of

%               i번째 피처에 대한 평균 mu(i)와 sigma2(i)

%





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

end


(4) 가우시안 파라미터 mu를 피처 별로 계산


   데이터 행렬 X에 대한 평균을 계산합니다. 데이터 행렬 X는 307 X 2차원입니다. mu에 평균을 입력합니다.  


[m, n] = size(X);              

mu = zeros(n, 1);           % mu를 초기화

mu = mean(X)';              % 열 별로 평균을 계산


>> [m, n] = size(X);

>> mu = zeros(n, 1);

>> mu = mean(X)          

mu =

   14.112   14.998


>> mu = mean(X)'         % 행 벡터를 열 벡터로 전환

mu =

   14.112

   14.998


(5) 가우시안 파라미터 sigma2를 피처 별로 계산


   데이터 행렬 X에 대한 분산을 계산합니다. 데이터 행렬 X는 307 X 2차원입니다. sigma2에 분산을 입력합니다. 분산을 계산하는 함수는 var()가 있지만, 여기서는 직접 아래 공식에 맞추어 계산합니다. var() 함수는 1/(m-1)로 나누기 때문에 약간 값의 차이가 있습니다.



>> mu

mu =

   14.112

   14.998


>> mu'

ans =

   14.112   14.998


>> X-mu'

ans =

  -1.0654e+00  -2.5656e-01

  -7.0371e-01  -1.2344e+00

   8.3689e-02   8.5547e-01

   8.0247e-01   1.1765e+00

  -5.3553e-01  -9.5486e-01

  ...


>> sum((X-mu').^2)

ans =

   562.62   524.89


>> m

m =  307


>> (1/m)* sum((X-mu').^2)

ans =

   1.8326   1.7097


>> (1/m)* sum((X-mu').^2)'

ans =

   1.8326

   1.7097


  따라서, 분산을 계산하는 코드는 다음과 같습니다.

 

sigma2 = (1/m)*sum((X-mu').^2)';



(6) 정답 


function [mu sigma2] = estimateGaussian(X)

%ESTIMATEGAUSSIAN  데이터 행렬 X의 피처 별 가우시안 분포의 파라미터를 계산

%   [mu sigma2] = estimateGaussian(X), 

%             입력 X :  데이터셋 X

%              mu : 피처 별 평균을 나타내는 n차원 벡터

%              sigma2 : 피처 별 분산을 나타내는 n차원 벡터 

% 변수 정의 

[m, n] = size(X);


% 반환할 변수를 초기화 

mu = zeros(n, 1);

sigma2 = zeros(n, 1);


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

% Instructions: 

%               데이터의 평균과 분산을 계산 In particular, mu(i) should contain the mean of

%               i번째 피처에 대한 평균 mu(i)와 sigma2(i)

%


mu = mean(X)';  

sigma2 = (1/m)*sum((X-mu').^2)';


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

end


(7) 결과 확인


>> [mu sigma2] = estimateGaussian(X)

mu =

   14.112

   14.998


sigma2 =

   1.8326

   1.7097



<Part 2: 데이터셋 통계 추정(중) : 가우시안 분포 추정>


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


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


p = multivariateGaussian(X, mu, sigma2);


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


(3)  multivariateGaussian.m 파일 분석


function p = multivariateGaussian(X, mu, Sigma2)

%MULTIVARIATEGAUSSIAN 다변량 가우시안 분포의 확률 밀도 함수를 계산 

%    p = MULTIVARIATEGAUSSIAN(X, mu, Sigma2) 

%    mu와 Sigma2로 다변량 가우시안 분포로 데이터 행렬 X의 확률 밀도 함수를 계산

%    Sigma2가 벡터라면 공분산 행렬의 각 차원의 분산을 sigma^2으로 취급                    

%


k = length(mu);


if (size(Sigma2, 2) == 1) || (size(Sigma2, 1) == 1)

    Sigma2 = diag(Sigma2);

end


X = bsxfun(@minus, X, mu(:)');


p = (2 * pi) ^ (- k / 2) * det(Sigma2) ^ (-0.5) * ...

    exp(-0.5 * sum(bsxfun(@times, X * pinv(Sigma2), X), 2));

end



   파일의 상세 내용을 분석합니다. 


%%% mu의 열의 길이를 계산

k = length(mu);


   피처 별 평균을 가진 벡터 mu의 길이를 계산합니다.

 

>> length (mu)

ans =  2


%%%  If 문의 이해

(size(Sigma2, 2) == 1) || (size(Sigma2, 1) == 1)


옥타브 프로그램에서 계산합니다. 두 가지 조건 중에서 하나라도 참이면 참입니다. 


Sigma2는 2 X 1 벡터이므로  첫 번째 조건은 참입니다. 


>> Sigma2 = sigma2

Sigma2 =

   1.8386

   1.7153


>> size(Sigma2,2)

ans =  1


>> size(Sigma2,2) == 1

ans =  1


Sigma2는 2 X 1 벡터이므로  두 번째 조건은 거짓입니다. 


>> Sigma2

Sigma2 =

   1.8386

   1.7153


>> size(Sigma2,1)

ans =  2


>> size(Sigma2,1) == 1

ans = 0



%%%% diag() 함수의 이해

    Sigma2 = diag(Sigma2);


   diag() 함수는 대각 행렬을 생성하거나 행렬의 대각선 요소를 반환합니다. 


>> v = [1 2 3 4 5]

v =

   1   2   3   4   5


>> diag(v)             % 벡터 v의 성분으로 대각 행렬 생성

ans =

Diagonal Matrix


   1   0   0   0   0

   0   2   0   0   0

   0   0   3   0   0

   0   0   0   4   0

   0   0   0   0   5


>> A = magic(3)

A =

   8   1   6

   3   5   7

   4   9   2


>> diag(A)            % 행렬 A의 대각선 성분을 반환

ans =

   8

   5

   2


>> sigma2

sigma2 =

   1.8386

   1.7153


>> diag(sigma2)          % 벡터 sigma2로 2 X2 정방 행렬로 대각 행렬 생성    

ans =

Diagonal Matrix

   1.8386        0

        0   1.7153



%%% 데이터 행렬 X에서 평균 mu를 계산

X = bsxfun(@minus, X, mu(:)'); 


  데이터 행렬 X에서 평균 mu를 뺀 값으로 데이터 행렬 X를 업데이트합니다. 데이터 행렬 X는 307 X 2차원 행렬이고, mu는 2 X 1차원 행렬입니다. 


>> mu

mu =

   14.112

   14.998


>> mu(:)

ans =

   14.112

   14.998


>> mu(:)'

ans =

   14.112   14.998


   mu(:)' 코드는 mu가 행 벡터 또는 열 벡터에 상관없이 무조건 행 벡터로 전환합니다.


%%% 다변량 가우시안 분포로 확률 밀도를 추정

p = (2 * pi) ^ (- k / 2) * det(Sigma2) ^ (-0.5) * ...

    exp(-0.5 * sum(bsxfun(@times, X * pinv(Sigma2), X), 2));


   다변량 가우시안 분포를 기준으로 확률 밀도를 추정하는 공식은 다음과 같습니다.

  이공식을 옥타브 프로그램 코드로 표현하였습니다. 

 



(4) 다변량 가우시안 분포를 활용한 확률 밀도를 추정


   multivariateGaussian.m 파일로 다변량 가우시안 분포를 활용한 확률 밀도를 추정합니다. 


>> p = multivariateGaussian(X, mu, sigma2)

p =


   6.4708e-02

   5.0304e-02

   7.2450e-02

   5.0316e-02

   6.3685e-02

   ...


<Part 2: 데이터셋 통계 추정(하) : 확률 밀도 추정을 도식화 >


(1) 데이터셋과 추정 분포 시각화


   데이터 분포를 2차원 평면에 그리는 visualizeFit.m 파일을 작성합니다. 이 파일을 호출할 때 데이터 행렬 X, 피처 별 평균 mu, 피처 별 분산 sigma2를 전달합니다. 


visualizeFit(X,  mu, sigma2);


   visualizeFit.m 파일은 다변량 가우시안 분포에 따른 추정 확률을 그립니다.  


(2) visualizeFit.m 파일 분석 

function visualizeFit(X, mu, sigma2)

%VISUALIZEFIT 데이터셋과 추정 분포를 시각화 

%   VISUALIZEFIT(X, p, mu, sigma2) 

%          가우시안 분포의 확률 밀도 함수를 시각화

%           각 예제는 x1과 x2 피처의 데이터를 보유

%


% 0부터 35까지 0.5씩 증가시킨 2차원 격자 구조를 반환

[X1,X2] = meshgrid(0:.5:35);   


% 다변량 가우시안 분포를 그리기

Z = multivariateGaussian([X1(:) X2(:)],mu,sigma2);

Z = reshape(Z,size(X1));


% 데이터 행렬 X를 2차원 평면에 그리기

plot(X(:, 1), X(:, 2),'bx');

hold on;


% Do not plot if there are infinities

if (sum(isinf(Z)) == 0)

    contour(X1, X2, Z, 10.^(-20:3:0)');

end

hold off;


end



%%% meshgrid() 함수의 이해

[X1,X2] = meshgrid(0:.5:35);


   meshgrid() 함수는 벡터를 입력받아 2차원 또는 3차원 격자 구조를 반환합니다.


>> meshgrid(1:5)         % 1에서 5까지 5 X 5 행렬을 반환

ans =

   1   2   3   4   5

   1   2   3   4   5

   1   2   3   4   5

   1   2   3   4   5

   1   2   3   4   5


>> v = 1:5                       

v =

   1   2   3   4   5


>> meshgrid(v)            % 벡터 A를 입력받아 5 X 5 행렬을 반환

ans =

   1   2   3   4   5

   1   2   3   4   5

   1   2   3   4   5

   1   2   3   4   5

   1   2   3   4   5


>> v = 1:3

v =

   1   2   3


>> [X1 X2] = meshgrid(v)   % 벡터 A로 2차원 공간에 맞는 행렬 X1과 X2를 반환

X1 =

   1   2   3

   1   2   3

   1   2   3


X2 =

   1   1   1

   2   2   2

   3   3   3


   [X Y Z ] = meshgrid(v)는 3차원 공간의 좌표를 벡터 v의 값으로 표현합니다. 각 축의 길이는 length(v)입니다.  


>>  meshgrid(1:0.1:1.2)   % 1부터 1.2까지 0.1씩 증가시킵니다. 

ans =

   1.0000   1.1000   1.2000

   1.0000   1.1000   1.2000

   1.0000   1.1000   1.2000


   따라서, meshgrid(0:. 5:35)는 0부터 35까지 0.5씩 숫자를 증가시킨 격자 구조를 반환합니다. 여기서 X1과 X2는 71 X 71차원 행렬입니다.


>> [X1,X2] = meshgrid(0:.5:35);   

>> size(X1)

ans =

   71   71


>> size(X2)

ans =

   71   71



% 격자구조의 다변량 가우시안 분포를 그리기

Z = multivariateGaussian([X1(:) X2(:)],mu,sigma2);

Z = reshape(Z,size(X1));


   격자구조의 데이터를 데이터 행렬로 변환합니다.. 


>> size(X1)         % X1은 71 X 71차원 행렬

ans =

   71   71


>> size(X1(:))     % 행렬 X1을 열 벡터로 변환

ans =

   5041      1


>> size(X2)        % X2은 71 X 71차원 행렬

ans =

   71   71


>> size(X2(:))    % 행렬 X2을 열 벡터로 변환

ans =

   5041      1


>> size([X1(:) X2(:)])    % 하나의 행렬로 합성

ans =

   5041      2


   multivariateGaussian.m 파일은 데이터 행렬 [X1(:) X2(:)]과 이미 계산된 mu와 sigma2를 입력받아 추정 확률을 Z로 반환합니다. 


>> Z = multivariateGaussian([X1(:) X2(:)],mu,sigma2);

>> size(Z)

ans =

   5041      1


   Z는 5041 X1 차원 행렬이므로 원래 값인 71 X 71의 값으로 재배열합니다. 


>> Z = reshape(Z,size(X1));

>> size(Z)

ans =

   71   71



%%% 원래 데이터 행렬 X를 그리기

plot(X(:, 1), X(:, 2),'bx');



%%% 등고선 그리기

  contour(X1, X2, Z, 10.^(-20:3:0)');



(3) 결과 확인


visualizeFit(X,  mu, sigma2);

xlabel('Latency (ms)');

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



브런치는 최신 브라우저에 최적화 되어있습니다. IE chrome safari