<jupyter_start><jupyter_text>ML Lab 8-  Anomaly Detection **Objectives**: implement Gaussian distribution model to detect anomalous behavior in server computers. While the servers were operating, 2D dataset regarding the through-put (mb/s) and latency (ms) of response of each server were collected. Vast majority of the examples were normal (non-anomalous) but also a few examples of servers acting anomalously.On that dataset you will ﬁt a Gaussian distribution and then ﬁnd values that have very low probability and hence can be considered anomalies.<jupyter_code># Loading relevant libraries and the dataset
import warnings
warnings.filterwarnings('ignore',category=FutureWarning)
warnings.filterwarnings('ignore',category=RuntimeWarning)

import numpy as np
import matplotlib.pyplot as plt
from scipy.io import loadmat
import matplotlib
%matplotlib inline<jupyter_output><empty_output><jupyter_text>Load Data<jupyter_code>mat = loadmat("ex8data1.mat")

Xtest = mat["X"] # normal and a few (6) abnormal examples
Xval = mat["Xval"] # normal and a few (9) abnormal examples
yval = mat["yval"]

print(Xtest.shape)

Xtrain = Xtest[0:300,:] # extract only normal examples<jupyter_output>(307, 2)<jupyter_text>Plot dataPlot the scater plots of the data subsets (Xtrain,Xval, Xtest) and get figure similar to Fig. 1.  **Fig. 1** : **Latency versus Throughput**<jupyter_code>plt.figure(figsize=(10,4))
plt.subplot(1, 3, 1)
plt.scatter(Xtrain[:,0], Xtrain[:,1], marker="x")
plt.xticks(np.arange(0, 40, 10))
plt.yticks(np.arange(0, 35, 5))
plt.xlabel("Latency (ms)")

plt.subplot(1, 3, 2)
plt.scatter(Xval[:,0], Xval[:,1], marker="x")
plt.xticks(np.arange(0, 40, 10))
plt.yticks(np.arange(0, 35, 5))

plt.subplot(1, 3, 3)
plt.scatter(Xtest[:,0], Xtest[:,1], marker="x")
plt.xticks(np.arange(0, 40, 10))
plt.yticks(np.arange(0, 35, 5))

plt.show()<jupyter_output><empty_output><jupyter_text>Plot DataPlot the histograms of each feature and get figure simular to Fig.2.  **Fig. 1** : **Data distribution**<jupyter_code>plt.subplot(1, 3, 1)
plt.hist(Xtrain[:,0], 20, density=1)

plt.subplot(1, 3, 2)
plt.hist(Xtrain[:,1], 20, density=1)
plt.show()<jupyter_output><empty_output><jupyter_text>Gaussian DistributionTo perform anomaly detection, you will first need to fit a model only to the normal data.Given a training set {x(1), …, x(m)} (where x(i) ∈ R^n, here n = 2), you want to estimate the Gaussian distribution for each of the features. For each feature (i = 1 . . . n), you need to find parameters mean and variance(mu, sigma²). For doing that let’s write down the function that calculates the mean and variance of matrix Xtrain.The Gaussian distribution is given by:  $p(x;\mu,\sigma^2) = \frac{1}{\sqrt{2\pi\sigma^2}}e^{\frac{(x-\mu)^2}{2\sigma^2}}$mean:  $\mu_i = \frac{1}{m}\sum^m_{j=1}x^{(j)}$variance:  $\sigma^2_i = \frac{1}{m}\sum^m_{j=1}(x^{(j)} - \mu_j)^2$<jupyter_code>def estimateGaussian(X):
    """
     This function estimates the parameters of a Gaussian distribution using the data in X
    """
    m = X.shape[0]
    
    # compute mean 
    mu = np.mean(X, axis=0)
     
    # compute variance =sigma^2
    sigma2 = (1/m) * np.sum(np.power(X - mu, 2), axis=0)
    
    return mu, sigma2
# Estimate parameters (mean and variance) for the Gaussian model of the normal data (Xtrain)
# ANSWER:    mu[0] = 14.06411334,      mu[1] = 15.02441479

mu, sigma2 = estimateGaussian(Xtrain)
print(mu)
print(sigma2)<jupyter_output>[14.06411334 15.02441479]
[0.93434723 0.8506542 ]<jupyter_text>Multivariate Gaussian DistributionInside *multivariateGaussian*, $\sigma^2$ vector is converted into the covariance matrix  $\Sigma$.  Assuming the features are independent, each element of vector $\sigma^2$ will be on the diagonal of matrix $\Sigma$ and the rest of the elements are zero. Then the formula for the multivariate distribution to get the probability vector is applied.$p(x;\mu,\Sigma) = \frac{1}{(2\pi)^{n/2}|\Sigma|^{1/2}} exp(-\frac{1}{2}(x-\mu)^T\Sigma^{-1}(x-\mu))$<jupyter_code>def multivariateGaussian(X, mu, sigma2):
    """
    Computes the probability density function of the multivariate Gaussian distribution.
    """
    n = X.shape[1]

    sigma_diag = np.diag(sigma2) # diagonal covariance matrix

    X = X - mu.T

    # cov = X.T * X  # non-diagonal covariance matrix

    p = 1/((2*np.pi)**(n/2)*(np.linalg.det(sigma_diag)**0.5))* np.exp(-0.5* np.sum(X @ np.linalg.pinv(sigma_diag) * X, axis=1))
    
    return p<jupyter_output><empty_output><jupyter_text>Select the threshold $\epsilon$  that will flag an example as anomalyF1 score is used to determine the best parameters i.e best epsilon and best F1 instead of accuracy as the dataset is highly unbalanced. F1 score takes into consideration precision and recall.Precision = true positive/(true positive + false positive)Recall = true positive /(true positive + false negative)F1=2(Recall * Precision) / (Recall + Precision)Best parameters are the ones in which the F1 score value is maximum.Predict anomaly if pval<epsilon that gives a vector of binary values in the variable predictions.Note: We need a try-except block because there can be cases where we divide by zero to calculate precision and recall.<jupyter_code>def selectThreshold(yval, pval):
    """
    Find the best threshold (epsilon) to use for selecting outliers
    """
    best_epi = 0
    best_F1 = 0
    
    # range of threshold values to decide the best one. 
    stepsize = (max(pval) - min(pval))/1000
    epi_range = np.arange(pval.min(), pval.max(), stepsize)
    for epi in epi_range:
        predictions = (pval < epi)[:,np.newaxis]
        
        # compute true positive (TP)
        tp = np.sum(predictions[yval == 1] == 1)
        
        # compute false positive (FP)
        fp = np.sum(predictions[yval == 0] == 1)
        
        # compute false negative 
        fn = np.sum(predictions[yval == 1] == 0)
        
        try:
            # compute precision, recall and F1
            prec = tp / (tp + fp)
            
            # compute recall
            rec = tp / (tp + fn)

            # compute F1
            F1 = 2*(rec * prec) / (rec + prec)
        
            if F1 > best_F1:
                best_F1 = F1
                best_epi = epi
            
        except ZeroDivisionError:
                print('Warning dividing by zero!!')  
        
    return best_epi, best_F1
# Apply multivariateGaussian to compute the probability of the examples in Xval. 
pval = multivariateGaussian(Xval, mu, sigma2)

#Apply selectThreshold to compute the best epsilon and best F1 with the validation data
epsilon, F1 = selectThreshold(yval, pval)

print("Best epsilon found using cross-validation:",epsilon)  # ANSWER: 0.00017841159895769934
print("Best F1 on Cross Validation Set:",F1)  # ANSWER: 0.823529411764706<jupyter_output>Best epsilon found using cross-validation: 0.00017841159895769934
Best F1 on Cross Validation Set: 0.823529411764706<jupyter_text>We need to return the indices of the outliers to identify the faulty servers. This gives us a vector with binary entries where 1 means anomaly and 0 means normal.<jupyter_code>def findIndices(binVec):
    l = []
    for i in range(len(binVec)):
        if binVec[i] == 1:
            l.append(i)
    return l
#Apply multivariateGaussian to detect abnormal examples (outliers) 
p = multivariateGaussian(Xtest, mu, sigma2)

outl = (p < epsilon)
listOfOutliers = findIndices(outl)
count_outliers = len(listOfOutliers)
print("\n\nNumber of outliers (Xtest):", count_outliers)  # ANSWER: 6 
print("\n", listOfOutliers)<jupyter_output>Number of outliers (Xtest): 6

 [300, 301, 303, 304, 305, 306]<jupyter_text>Visualize the anomalies<jupyter_code>plt.figure(figsize=(10,4))
plt.subplot(1, 2, 1)
plt.scatter(Xtest[:, 0], Xtest[:, 1], marker = "x")
plt.xlabel("Latency(ms)")
plt.ylabel("Throughput(mb/s)")
plt.scatter(Xtest[listOfOutliers,0], Xtest[listOfOutliers, 1], facecolors = 'none', edgecolors = 'r')
plt.xlim(-1, 30)
plt.ylim(-1, 30)
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12367402970173468, 'fraction_numerical': 0.04130787470360664, 'mean_word_length': 3.965303593556382, 'pattern_counts': {'":': 0, '<': 36, '<?xml version=': 0, '>': 34, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '644789', 'license': 'MIT', 'n_tokens_mistral': 2940, 'n_tokens_neox': 2621, 'n_words': 828, 'path': 'practical-labs/lab08/Anomaly_Detection.ipynb', 'repo_name': 'vascoalramos/taa', 'size': '8013'}