mahalanobis_distance/anomaly-detection-training-mahalanobis-distance.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "#python -m pip install seaborn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from os import listdir\n",
    "from os.path import join\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import scipy as sp\n",
    "import matplotlib.pyplot as plt\n",
    "import random\n",
    "from scipy import stats\n",
    "from sklearn.metrics import confusion_matrix\n",
    "import seaborn as sn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Python 3.7.6\n",
      "Numpy 1.18.1\n",
      "Pandas 1.0.3\n",
      "SciPy 1.4.1\n"
     ]
    }
   ],
   "source": [
    "# Print versions\n",
    "!python --version\n",
    "print('Numpy ' + np.__version__)\n",
    "print('Pandas ' + pd.__version__)\n",
    "print('SciPy ' + sp.__version__)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Max measurements per file: 128\n",
      "Downsample factor: 2\n"
     ]
    }
   ],
   "source": [
    "# Settings\n",
    "dataset_path = '..\\\\datasets\\\\ceiling-fan'  # Directory where raw accelerometer data is stored\n",
    "normal_op_list = ['fan_0_low_0_weight-moving']\n",
    "anomaly_op_list = ['fan_0_med_0_weight', 'fan_0_high_0_weight',\n",
    "                  'fan_0_low_1_weight', 'fan_0_med_1_weight', 'fan_0_high_1_weight']\n",
    "val_ratio = 0.2             # Percentage of samples that should be held for validation set\n",
    "test_ratio = 0.2            # Percentage of samples that should be held for test set\n",
    "raw_scale = 1               # Multiply raw values to fit into integers\n",
    "sensor_sample_rate = 200    # Hz\n",
    "desired_sample_rate = 50    # Hz\n",
    "sample_time = 0.64           # Time (sec) length of each sample\n",
    "samples_per_file = 128      # Expected number of measurements in each file (truncate to this)\n",
    "max_anomaly_samples = 500   # Max number of (shuffled) anomaly samples to use\n",
    "max_measurements = int(sample_time * sensor_sample_rate)\n",
    "downsample_factor = int(samples_per_file / desired_sample_rate)\n",
    "\n",
    "md_file_name = 'models\\\\md_model-moving'   # Mahalanobis Distance model arrays (.npz will be added)\n",
    "\n",
    "print('Max measurements per file:', max_measurements)\n",
    "print('Downsample factor:', downsample_factor)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create list of filenames\n",
    "def createFilenameList(op_list):\n",
    "    \n",
    "    # Extract paths and filenames in each directory\n",
    "    op_filenames = []\n",
    "    num_samples = 0\n",
    "    for index, target in enumerate(op_list):\n",
    "        samples_in_dir = listdir(join(dataset_path, target))\n",
    "        samples_in_dir = [join(dataset_path, target, sample) for sample in samples_in_dir]\n",
    "        op_filenames.append(samples_in_dir)\n",
    "    \n",
    "    # Flatten list\n",
    "    return [item for sublist in op_filenames for item in sublist]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of normal samples: 299\n",
      "Number of anomaly samples: 2926\n"
     ]
    }
   ],
   "source": [
    "# Create normal and anomaly filename lists\n",
    "normal_op_filenames = createFilenameList(normal_op_list)\n",
    "anomaly_op_filenames = createFilenameList(anomaly_op_list)\n",
    "print('Number of normal samples:', len(normal_op_filenames))\n",
    "print('Number of anomaly samples:', len(anomaly_op_filenames))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Shuffle lists\n",
    "random.shuffle(normal_op_filenames)\n",
    "random.shuffle(anomaly_op_filenames)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate validation and test set sizes\n",
    "val_set_size = int(len(normal_op_filenames) * val_ratio)\n",
    "test_set_size = int(len(normal_op_filenames) * test_ratio)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of training samples: 181\n",
      "Number of validation samples: 59\n",
      "Number of test samples: 59\n"
     ]
    }
   ],
   "source": [
    "# Break dataset apart into train, validation, and test sets\n",
    "num_samples = len(normal_op_filenames)\n",
    "filenames_val = normal_op_filenames[:val_set_size]\n",
    "filenames_test = normal_op_filenames[val_set_size:(val_set_size + test_set_size)]\n",
    "filenames_train = normal_op_filenames[(val_set_size + test_set_size):]\n",
    "\n",
    "# Print out number of samples in each set\n",
    "print('Number of training samples:', len(filenames_train))\n",
    "print('Number of validation samples:', len(filenames_val))\n",
    "print('Number of test samples:', len(filenames_test))\n",
    "\n",
    "# Check that our splits add up correctly\n",
    "assert(len(filenames_train) + len(filenames_val) + len(filenames_test)) == num_samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function: extract specified features (variances, MAD) from sample\n",
    "def extract_features(sample, max_measurements=0, scale=1):\n",
    "    \n",
    "    features = []\n",
    "    \n",
    "    # Truncate sample\n",
    "    if max_measurements == 0:\n",
    "        max_measurements = sample.shape[0]\n",
    "    sample = sample[0:max_measurements]\n",
    "    \n",
    "    # Scale sample\n",
    "    sample = scale * sample\n",
    "    \n",
    "#     # Variance\n",
    "#     features.append(np.var(sample, axis=0))\n",
    "    \n",
    "#     # Kurtosis\n",
    "#     features.append(stats.kurtosis(sample))\n",
    "    \n",
    "#     # Skew\n",
    "#     features.append(stats.skew(sample))\n",
    "    \n",
    "    # Median absolute deviation (MAD)\n",
    "    features.append(stats.median_absolute_deviation(sample))\n",
    "    \n",
    "#     # Correlation\n",
    "#     cov = np.corrcoef(sample.T)\n",
    "#     features.append(np.array([cov[0,1], cov[0,2], cov[1,2]]))\n",
    "    \n",
    "    # Compute a windowed FFT of each axis in the sample (leave off DC)\n",
    "#     sample = sample[::downsample_factor, :]  # Downsample\n",
    "#     sample = np.floor(sample)                # Round down to int\n",
    "#     hann_window = np.hanning(sample.shape[0])\n",
    "#     for i, axis in enumerate(sample.T):\n",
    "#         fft = abs(np.fft.rfft(axis * hann_window))\n",
    "#         features.append(fft[1:])  # Leave off DC\n",
    "    \n",
    "    return np.array(features).flatten()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(3,)\n",
      "[0.02714641 0.01484083 0.01158355]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1a707c7d988>]"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD4CAYAAADlwTGnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8GearUAAAgAElEQVR4nO3de3RV9Z3+8ffnnFyAcE8CcgkEQryAomAgUcSqtRW0lt7GQRRUbmPVTu30ZjvT2llTu7QXOzp1pCBUUdD2Z7VltVrrqC1im0BAUREkAUSCIAlCuJPb5/fH2WIaAzmBJDvJeV5rZZGz93efPCdrc56c797nbHN3REQk8UTCDiAiIuFQAYiIJCgVgIhIglIBiIgkKBWAiEiCSgo7QHNkZGR4dnZ22DFERDqU1atXV7h7ZsPlHaoAsrOzKS4uDjuGiEiHYmZbG1uuKSARkQSlAhARSVAqABGRBKUCEBFJUCoAEZEEpQIQEUlQKgARkQSVEAXw2ra9zPvrprBjiIi0Kx3qjWAn66k1ZSz++1bc4cuX5IQdR0SkXUiIArjz6lFUHq7mnj9tIDlqzJ44POxIIiKhS4gCiEaMn/3TuVTX1vHDP64nJSnCjAuyw44lIhKqhCgAgKRohPumjqGqZg3f//06kqMRrh0/JOxYIiKhSYiDwB9KjkZ44LoxXHJGJt99+g2eXF0WdiQRkdAkVAEApCZFmXf9+UzIyeBbT67l969tDzuSiEgoEq4AALokR1kwI49x2X35t9+s5dk3doQdSUSkzSVkAQB0TYmy6MZxnJfVm688/irPv/V+2JFERNpUwhYAQFpqEr+6aRyjBvbk1iVr+Mvbu8KOJCLSZhK6AAB6dklm8cx8cvt3Z+6jq1lRUhF2JBGRNpHwBQDQq1syj83KZ3hGGrMXr6Jw8+6wI4mItDoVQKBPWgqPzc5ncJ9uzHx4Fau3fhB2JBGRVhVXAZjZJDN728xKzeyORtabmd0frH/dzMYGy7PM7CUzW29m68zsqw22+0pwv+vM7Mct85BOXkb3VJbOzqd/zy7cuGgVa7ftDTuSiEirabIAzCwKPABMBkYC15rZyAbDJgO5wddc4MFgeQ3wdXc/CygAbv1wWzO7FJgCjHb3UcBPT/3hnLp+PbuwdE4+fdJSmL6wiDe3V4YdSUSkVcTzCmA8UOrum929CniC2BN3fVOAxR5TCPQ2swHuvsPd1wC4+35gPTAo2ObLwN3ufjRY325OwRnQqytL5+TTo0sy1y8sYsPOfWFHEhFpcfEUwCBgW73bZXz0JB73GDPLBsYARcGi04GJZlZkZn81s3GN/XAzm2tmxWZWXF5eHkfcljG4TzeWzsmnS1KU6xYUUbprf5v9bBGRthBPAVgjy7w5Y8ysO/Bb4HZ3//DP6SSgD7GpoW8CvzGzj92Pu8939zx3z8vMzIwjbssZmp7Gkjn5mBnTFhSxpeJgm/58EZHWFE8BlAFZ9W4PBt6Ld4yZJRN78l/i7k812OapYNpoJVAHZDQvfuvLyezO43Pyqa1zpi0o5N3dh8KOJCLSIuIpgFVArpkNM7MUYCqwrMGYZcCM4GygAqDS3XcEf9EvBNa7+70NtvkdcBmAmZ0OpADt8l1Yuf178NjsfA5X13LtgkLK9qgERKTja7IA3L0GuA14jthB3N+4+zozu9nMbg6GPQNsBkqBBcAtwfIJwHTgMjN7Lfi6Mli3CBhuZm8SO7B8g7s3nFpqN84a0JPHZuWz70g10xYUsbPySNiRREROibXj59yPycvL8+Li4lAzvPruHqYvXEm/Hqk88S8F9OvRJdQ8IiJNMbPV7p7XcLneCdxMY4b04Vc3jWPnviNct6CI3QeOhh1JROSkqABOwrjsviy8YRzb9hziuoeK2HOwKuxIIiLNpgI4SRfkpLNgRh6bKw4yfVERlYerw44kItIsKoBTMDE3k19efz5v79zPDYtWsv+ISkBEOg4VwCm69Mx+PDBtLG9ur2Tmw6s4eLQm7EgiInFRAbSAT486jfumjmH11j3MemQVh6tqw44kItIkFUALuWr0AH7+z+dRtOUD5j5azJFqlYCItG8qgBY05bxB3PPF0bxcUsEtS9ZQVVMXdiQRkeNSAbSwa/Ky+NHnz+HFDbu4bekaqmtVAiLSPqkAWsG0/CH852dH8ee33uf2J16jRiUgIu1QUtgBOqsbLsymqqaOu55ZT3LU+Nk15xGNNPap2SIi4VABtKI5Fw+nqraOnzz3NsnRCPd8cTQRlYCItBMqgFZ266UjqKqp474XSkhOinDX586mkeveiIi0ORVAG7j98lyqaut48C+bSIlGuPPqkSoBEQmdCqANmBnfuuIMqmrqWLhiCylJEb4z+UyVgIiESgXQRsyM/7jqLKpr65i/fDMp0QjfuOKMsGOJSAKL6zRQM5tkZm+bWamZ3dHIejOz+4P1r5vZ2GB5lpm9ZGbrzWydmX21kW2/YWZuZu3uesAtzcz4wdWjuHZ8Fr94qZT7XygJO5KIJLAmXwGYWRR4APgUsQu5rzKzZe7+Vr1hk4Hc4CsfeDD4twb4uruvMbMewGoze/7Dbc0sK7jfd1vwMbVrkYhx1+fOoarGuff5jSRHI3z5kpywY4lIAopnCmg8UOrumwHM7AlgClC/AKYAi4Nr+haaWW8zG+DuO4AdAO6+38zWA4Pqbftz4FvA71vk0XQQkYjx4y+Nprq2jnv+tIHkqDF74vCwY4lIgomnAAYB2+rdLiP2131TYwYRPPkDmFk2MAYoCm5/Ftju7mtPdDDUzOYCcwGGDBkSR9yOIRox7r3mXKpr6/jhH9eTmhRh+gXZYccSkQQSzzGAxp6dG15J/oRjzKw78FvgdnffZ2bdgH8Hvt/UD3f3+e6e5+55mZmZccTtOJKiEe6bOobLz+rP936/jidWJsxMmIi0A/EUQBmQVe/2YOC9eMeYWTKxJ/8l7v5UsD4HGAasNbN3gvFrzOy05j6Aji4lKcID143hkjMy+c7Tb/Dk6rKwI4lIgoinAFYBuWY2zMxSgKnAsgZjlgEzgrOBCoBKd99hsbmdhcB6d7/3w8Hu/oa793P3bHfPJlYgY919Z0s8qI4mNSnKvOvPZ0JOBt96ci3L1jbsVxGRltdkAbh7DXAb8BywHviNu68zs5vN7OZg2DPAZqAUWADcEiyfAEwHLjOz14KvK1v6QXQGXZKjLJiRR152X77269d49o0dTW8kInIKLHbiTseQl5fnxcXFYcdoVQeO1nDDopWs3baXedefz+Uj+4cdSUQ6ODNb7e55DZfregDtTPfUJH510zhGDezJLUvW8Je3d4UdSUQ6KRVAO9SzSzKLZ+aT2787cx9dzSulFWFHEpFOSAXQTvXqlsyjs/IZnpHGrEdWUbR5d9iRRKSTUQG0Y33TUnhsdj6D+3TjpodXsXrrB2FHEpFORAXQzmV0T2Xp7Hz69+zCjYtWsXbb3rAjiUgnoQLoAPr17MLSOfn0Tktm+sIi3txeGXYkEekEVAAdxIBeXVk6u4AeXWIlsGHnvrAjiUgHpwLoQLL6dmPpnHxSk6Jc/1ARpbv2hx1JRDowFUAHMzQ9jSVz8gFj2oIitlQcDDuSiHRQKoAOKCezO0vn5FNT50xbUMi2Dw6FHUlEOiAVQAd1ev8ePDYrn8PVtUydX8j2vYfDjiQiHYwKoAMbObAnj83KZ9+RaqYtKGRn5ZGwI4lIB6IC6ODOHtSLxTPHs/tAFdMWFLJrv0pAROKjAugExgzpw69uGsfOfUe4bkERuw8cDTuSiHQAKoBOYlx2XxbeMI5tew5x3UNF7DlYFXYkEWnnVACdyAU56SyYkcfmioNMX1RE5eHqsCOJSDsWVwGY2SQze9vMSs3sjkbWm5ndH6x/3czGBsuzzOwlM1tvZuvM7Kv1tvmJmW0Ixj9tZr1b7mElrom5mfzy+vN5e+d+bli0kv1HVAIi0rgmC8DMosADwGRgJHCtmY1sMGwykBt8zQUeDJbXAF9397OAAuDWets+D5zt7qOBjcB3TvGxSODSM/vxi2ljeXN7JTMfXsXBozVhRxKRdiieVwDjgVJ33+zuVcATwJQGY6YAiz2mEOhtZgPcfYe7rwFw9/3Erik8KLj95+B6wwCFwOAWeDwSuGLUadw3dQyrt+5h9iPFHK6qDTuSiLQz8RTAIGBbvdtlwbJmjTGzbGAMUNTIz5gJPBtHFmmGq0YP4Of/fB6FW3Yz99FijlSrBETkI/EUgDWyrOGV5E84xsy6A78Fbnf3f/gYSzP7d2JTRUsa/eFmc82s2MyKy8vL44gr9U05bxD3fHE0L5dUcMuSNVTV1IUdSUTaiXgKoAzIqnd7MPBevGPMLJnYk/8Sd3+q/kZmdgPwGeA6d29YKgC4+3x3z3P3vMzMzDjiSkPX5GVx1+fP5sUNu7ht6Rqqa1UCIhJfAawCcs1smJmlAFOBZQ3GLANmBGcDFQCV7r7DzAxYCKx393vrb2Bmk4BvA591d32aWSu7Ln8oP7h6JH9+631u//Vr1KgERBJeUlMD3L3GzG4DngOiwCJ3X2dmNwfr5wHPAFcCpcAh4KZg8wnAdOANM3stWPZdd38G+AWQCjwf6wkK3f3mFntk8jE3ThhGda1z1zPrSY4YP7vmPKKRxmbvRCQRNFkAAMET9jMNls2r970Dtzay3QoaPz6Au49oVlJpEXMuHk5VbR0/ee5tkqMR7vniaCIqAZGEFFcBSOdy66UjOFpTx/0vlJCSFOGHnzub4FWYiCQQFUCC+trluVTX1vHgXzaRHI1w59UjVQIiCUYFkKDMjG9dcQZVNXUsXLGFlKQI35l8pkpAJIGoABKYmfEfV51FdW0d85dvJiUa4RtXnBF2LBFpIyqABGdm/ODqUVTV1PGLl0pJSYrwr5/MDTuWiLQBFYAQiRg/+vw5VNc69z6/keRohC9fkhN2LBFpZSoAAWIl8OMvjaa6to57/rSB5Kgxe+LwsGOJSCtSAcgx0Yhx7zXnUl1bxw//uJ7UpAjTL8gOO5aItBJdEUz+QVI0wn1Tx3D5Wf343u/X8cTKd8OOJCKtRAUgH5OSFOGB68byidMz+c7Tb/Db1WVhRxKRVqACkEalJkX55fTzmZCTwTefXMuytQ0/AFZEOjoVgBxXl+QoC2bkkZfdl6/9+jWefWNH2JFEpAWpAOSEuqZEWXTjOM4d3IuvPP4q//fW+2FHEpEWogKQJnVPTeLhmeMZNbAntyxZw1/e3hV2JBFpASoAiUvPLsksnplPbv/uzH10Na+UVoQdSUROkQpA4tarWzKPzspneEYasx5ZRdHm3WFHEpFToAKQZumblsJjs/MZ1LsrMx9exeqte8KOJCInKa4CMLNJZva2mZWa2R2NrDczuz9Y/7qZjQ2WZ5nZS2a23szWmdlX623T18yeN7OS4N8+LfewpDVldE/l8TkF9OvZhRsXrWTttr1hRxKRk9BkAZhZFHgAmAyMBK41s5ENhk0GcoOvucCDwfIa4OvufhZQANxab9s7gBfcPRd4IbgtHUS/nl1YOief3mnJTF9YxJvbK8OOJCLNFM8rgPFAqbtvdvcq4AlgSoMxU4DFHlMI9DazAe6+w93XALj7fmA9MKjeNo8E3z8CfO4UH4u0sQG9urJ0dgE9usRKYMPOfWFHEpFmiKcABgHb6t0u46Mn8bjHmFk2MAYoChb1d/cdAMG//Rr74WY218yKzay4vLw8jrjSlrL6dmPpnHxSkiJc/1ARpbsOhB1JROIUTwE0do1Ab84YM+sO/Ba43d2b9Weiu8939zx3z8vMzGzOptJGhqansXROAWBMW1DIloqDYUcSkTjEUwBlQFa924OBhh8Mc9wxZpZM7Ml/ibs/VW/M+2Y2IBgzANC7izqwnMzuLJ2TT02dM21BIds+OBR2JBFpQjwFsArINbNhZpYCTAWWNRizDJgRnA1UAFS6+w6LXWF8IbDe3e9tZJsbgu9vAH5/0o9C2oXT+/fgsVn5HKqqZer8QrbvPRx2JBE5gSYLwN1rgNuA54gdxP2Nu68zs5vN7OZg2DPAZqAUWADcEiyfAEwHLjOz14KvK4N1dwOfMrMS4FPBbengRg7syWOz8tl3pJppCwrZWXkk7Egichzm3nA6v/3Ky8vz4uLisGNIHF59dw/TF66kX89UnphbQL8eXcKOJJKwzGy1u+c1XK53AkurGDOkD7+6aRw7K49w3YIidh84GnYkEWlABSCtZlx2Xx66IY93PzjE9QtXsvdQVdiRRKQeFYC0qgtzMlgwI49Nuw4wfeFKKg9Xhx1JRAIqAGl1F5+eybzpY9mwcx83/mol+4+oBETaAxWAtInLzuzPL6aN5Y2ySmY+vIqDR2vCjiSS8FQA0mauGHUa900dw+qte5j9SDGHq2rDjiSS0FQA0qauGj2Ae685j8Itu5n7aDFHqlUCImFRAUib+9yYQdzzxdG8XFLBLUvWUFVTF3YkkYSkApBQXJOXxV2fP5sXN+zitqVrqK5VCYi0NRWAhOa6/KH84OqR/Pmt97n9169RoxIQaVNJYQeQxHbjhGFU1dbxo2c2kBKN8NN/OpdopLFPFxeRlqYCkNDNvTiH6lrnJ8+9TXLUuPsLo4moBERanQpA2oVbLx3B0Zo67n+hhORohB9+7mxinyYuIq1FBSDtxtcuz6Wqpo55f91EcjTCnVePVAmItCIVgLQbZsa3J51BdW0dC1dsITUpwh2Tz1QJiLQSFYC0K2bGf1x1FtW1dfxy+WZSkiJ8/dNnhB1LpFOK6zRQM5tkZm+bWamZ3dHIejOz+4P1r5vZ2HrrFpnZLjN7s8E255lZYXCVsGIzG3/qD0c6AzPjB1ePYuq4LP7nxVLuf6Ek7EginVKTBWBmUeABYDIwErjWzEY2GDYZyA2+5gIP1lv3MDCpkbv+MfCf7n4e8P3gtggAkYjxo8+fwxfGDuLe5zcy76+bwo4k0unEMwU0Hih1980AZvYEMAV4q96YKcBij11fstDMepvZAHff4e7LzSy7kft1oGfwfS/gvZN8DNJJRSLGT750LjW1zt3PbiA5GmHWRcPCjiXSacRTAIOAbfVulwH5cYwZBOw4wf3eDjxnZj8l9krkwsYGmdlcYq8qGDJkSBxxpTOJRox7rzmX6to6/usPb5ESNaZfkB12LJFOIZ5jAI2dgtHwSvLxjGnoy8DX3D0L+BqwsLFB7j7f3fPcPS8zM7PJsNL5JEUj3Dd1DJef1Y/v/X4dT6x8N+xIIp1CPAVQBmTVuz2Yj0/XxDOmoRuAp4Lv/x+xqSaRRqUkRXjgurF84vRMvvP0G/x2dVnYkUQ6vHgKYBWQa2bDzCwFmAosazBmGTAjOBuoAKh09xNN/0CsID4RfH8ZoFM95IRSk6L8cvr5XJiTzjefXMuytTpsJHIqmiwAd68BbgOeA9YDv3H3dWZ2s5ndHAx7BtgMlAILgFs+3N7MHgf+DpxhZmVmNitYNQf4mZmtBX5EMM8vciJdkqM8NGMcedl9+dqvX+PZN5r6O0NEjsdiJ+50DHl5eV5cXBx2DGkHDhytYcbCIl4vq2Te9edz+cj+YUcSabfMbLW75zVcrusBSIfUPTWJh2eOZ9TAntyyZA1/eXtX2JFEOhwVgHRYPbsks3hmPiP6dedfHl3NK6UVYUcS6VBUANKh9eqWzGOz8xmWkcasR1ZRtHl32JFEOgwVgHR4fdNSeGx2PoN6d2Xmw6tYvXVP2JFEOgQVgHQKGd1TWTqngMweqdy4aCVrt+0NO5JIu6cCkE6jf88uLJ1TQO+0ZKYvLOLN7ZVhRxJp11QA0qkM7N2VpbML6NElVgIbdu4LO5JIu6UCkE4nq283ls7JJyUpwvUPFVG660DYkUTaJRWAdEpD09NYOqcAMKYtKGRLxcGwI4m0OyoA6bRyMruzdE4+NXXOtAWFbPvgUNiRRNoVFYB0aqf378Fjs/I5VFXL1PmFbN97OOxIIu2GCkA6vZEDe/LYrHz2Halm2oJCdlYeCTuSSLugApCEcM7gXjwyczwV+48y7aFCdu1XCYioACRhjB3Sh4dnjmfH3iNc/1ARuw8cDTuSSKhUAJJQxmX3ZeGNeWzdfYjrF65k76GqsCOJhEYFIAnnwpwMFszIY9OuA0xfuJLKw9VhRxIJhQpAEtLFp2cyb/pYNuzcx42/WsmBozVhRxJpc3EVgJlNMrO3zazUzO5oZL2Z2f3B+tfNbGy9dYvMbJeZvdnIdl8J7nedmf341B6KSPNcdmZ/fjFtLG+UVXLTr1ZyqEolIImlyQIwsyjwADAZGAlca2YjGwybDOQGX3OBB+utexiY1Mj9XgpMAUa7+yjgpyeRX+SUXDHqNO6bOobVW/cw6+FiDlfVhh1JpM3E8wpgPFDq7pvdvQp4gtgTd31TgMUeUwj0NrMBAO6+HPigkfv9MnC3ux8NxumafhKKq0YP4N5rzqNwy27mPlrMkWqVgCSGeApgELCt3u2yYFlzxzR0OjDRzIrM7K9mNq6xQWY218yKzay4vLw8jrgizfe5MYO45wujebmkgluWrKGqpi7sSCKtLp4CsEaW+UmMaSgJ6AMUAN8EfmNmH7sfd5/v7nnunpeZmRlHXJGTc824LO76/Nm8uGEXX3l8DdW1KgHp3OIpgDIgq97twcB7JzGmsft9Kpg2WgnUARlx5BFpNdflD+UHV4/kuXXvc/uvX6NGJSCdWDwFsArINbNhZpYCTAWWNRizDJgRnA1UAFS6+44m7vd3wGUAZnY6kAJUNCu9SCu4ccIwvnvlmfzx9R1888nXqa1r6sWsSMeU1NQAd68xs9uA54AosMjd15nZzcH6ecAzwJVAKXAIuOnD7c3sceASIMPMyoA73X0hsAhYFJweWgXc4O76nybtwtyLc6iqqeOnf95IctS4+wujiUQam+kU6biaLAAAd3+G2JN8/WXz6n3vwK3H2fba4yyvAq6PO6lIG7vtslyqap37XyghORrhh587m0YOU4l0WHEVgEii+trluVTV1DHvr5tIjka48+qRKgHpNFQAIidgZnx70hlU1dSx6JUtpCZFuGPymSoB6RRUACJNMDO+95mzqK6t45fLN5OSFOHrnz4j7Fgip0wFIBIHM+M/PzuK6to6/ufFUpKjEf71k7lhxxI5JSoAkThFIsaPPn8OVbV13Pv8RlKSItz8iZywY4mcNBWASDNEIsZPvnQu1bXO3c9uIDkaYdZFw8KOJXJSVAAizRSNGPdecy7VNXX81x/eIiUpwvSCoWHHEmk2XRBG5CQkRyPcf+0YLj+rH9/73ZssKdpKnd4xLB2MdaQ33+bl5XlxcXHYMUSOOVpTy9zFq/nrxnL6pqVw0YgMLsrNYGJuBgN6dQ07nggAZrba3fMaLtcUkMgpSE2KMn/G+fzx9R28XFLByyUVLFsb+xzEEf26MzE3g4tzM8kf3pduKfrvJu2LXgGItCB3Z8PO/awoqWB5STkrt3zA0Zo6kqPG+UP7MDE3k4m5GYwa2IuoPltI2sjxXgGoAERa0ZHqWorf2cPLJeW8XFLBWzv2AdC7WzITRmRwcW4GF+VmMqi3pouk9agARNqB8v1H+dumCpZvrODlknJ27T8KwPDMNC7OzeSiERkU5KTTPVXTRdJyVAAi7Yy7U7LrAMs3lrOitILCzbs5Ul1HUsQYO6QPE3NjB5RHD+6t6SI5JSoAkXbuaE0tq9/Zw8ulsVcHb26PTRf16prMhBHpXDQidvwgq2+3kJNKR6MCEOlgdh84yiubdvNy8AphR+URALLTux07mHxBTjo9uiSHnFTau1MqADObBNxH7IpgD7n73Q3WW7D+SmJXBLvR3dcE6xYBnwF2ufvZjdz3N4CfAJnufsJLQqoAJFG5O5vKDxw71bRw824OVdUSjRhjsnoH7z3I5NzBvUiK6v2d8o9OugDMLApsBD5F7ELuq4Br3f2temOuBL5CrADygfvcPT9YdzFwAFjcsADMLAt4CDgTOF8FIBKfqpo61rwbO7toRUkFr2+vxB16dEniwpz0Y68QhqanhR1V2oFTeSPYeKDU3TcHd/QEMAV4q96YKcSe4B0oNLPeZjbA3Xe4+3Izyz7Off8c+Bbw+/gfioikJEUoGJ5OwfB0vnkF7DlYxSubKlgRvEJ4bt37AAzp242JwTuTL8jJoFdXTRfJR+IpgEHAtnq3y4j9ld/UmEHAjuPdqZl9Ftju7mtPdHUlM5sLzAUYMmRIHHFFEk+ftBQ+M3ognxk9EHdnS8XBY9NFv3t1O0uK3iVicG5WbybmZnJxbgbnZvUmWdNFCS2eAmjs2bnhvFE8Yz4abNYN+Hfg0039cHefD8yH2BRQU+NFEp2ZMTyzO8Mzu3PDhdlU19bx6rt7WVFSzvKSCn7xYgn3v1BC99QkLshJD14hZJKd3k2Xukww8RRAGZBV7/Zg4L2TGFNfDjAM+PCv/8HAGjMb7+4748gkInFKjkYYP6wv44f15d8+fQaVh6pjb0YriZ1u+vxbsemiwX26HiuDC3PS6d0tJeTk0triKYBVQK6ZDQO2A1OBaQ3GLANuC44P5AOV7n7c6R93fwPo9+FtM3sHyGvqILCInLpe3ZKZfM4AJp8zAHdn6+5DsfcebCznD2t38PjKbUQMzhncO/ZRFSMyGDOkDylJmi7qbJosAHevMbPbgOeInQa6yN3XmdnNwfp5wDPEzgAqJXYa6E0fbm9mjwOXABlmVgbc6e4LW/qBiEjzmRnZGWlkZ6QxvWAoNbV1rC3by/KNFaworeB//7KJ/3mxlLSUKAXD04N3J2eSk5mm6aJOQG8EE5Hj2nekmr9v2n3sw+y27j4EwMBeXZiYm8lFuRlMGJFB3zRNF7VneiewiJyyd3cf4uXScl7eWMHfNlWw70gNZnDOoF5cNCJ2/GDs0N6kJkXDjir1qABEpEXV1Nbx+vbK4L0H5bz67l5q6pyuyVEKhvc99ma0Ef26a7ooZCoAEWlV+49UU7j5g2PvTt5ccRCA03p2OXaZzItGZJDePTXkpIlHBSAibWrbB4dYURp7d/KK0goqD1cDMGpgz2OvDs4f2ocuyZouam0qABEJTW2d8+b2ymMHk1dv3UNNndMlOUL+sI/ejHZ6f00XtQYVgEDnSDoAAAmRSURBVIi0GweO1lC0eXfwcRXlbCqPTRf165HKRbkZXJybyYQRGWT20HRRSziVD4MTEWlR3VOT+ORZ/fnkWf0BeG/vYVaUVLC8pJyXNuziqTXbAThrQM9jH2Y3LruvpotamF4BiEi7UlfnrHtv37HTTYu3fkB1rZOaFPtIi9jB5EzOGtBD00Vx0hSQiHRIh6pqKNryAS9vrGBFaTkb3z8AQEb31GNnFk3MzaBfzy4hJ22/NAUkIh1St5QkLj2jH5eeEfv4sJ2VR1gRXDd5+cZynn41Nl10Rv8ewUdVZJA/LJ2uKZouaopeAYhIh1VX56zfue/YweRV7+yhqqaOlGiEccP6cNGI2OmmIwf0JBJJ3OkiTQGJSKd3uKqWle98wIrgdNMNO/cDkJ6WwoRgqmhibian9Uqs6SJNAYlIp9c1JconTs/kE6dnArBr34fTRbGvZWtjlynJ7df92Omm+cP70i0lMZ8K9QpARBKCu7Nh5/5jp5uu3PIBR2vqSI4a5w/tc+zdyWcP7NXppos0BSQiUs+R6lqK39lz7N3Jb+3YB0CfbsnHposuys1kUO+uISc9dSoAEZETKN9/NHapzI2xA8q79h8FYHhmGhfnZnLRiAwKctLpntrxpotOqQDMbBJwH7Ergj3k7nc3WG/B+iuJXRHsRndfE6xbBHwG2OXuZ9fb5ifA1UAVsAm4yd33niiHCkBE2oK7U7LrAMs3lrOitILCzbs5Ul1HUsQYO7QPE0dkMPH0TM4Z1ItoB5guOukCMLMosBH4FLGLv68CrnX3t+qNuRL4CrECyAfuc/f8YN3FwAFgcYMC+DTwYnDJyXsA3P3bJ8qiAhCRMBytqWX1O3ti104uKefN7bHpol5dk5kwIj12dbQRGWT17RZy0sadyllA44FSd98c3NETwBTgrXpjphB7gneg0Mx6m9kAd9/h7svNLLvhnbr7n+vdLAS+FPejERFpQ6lJUS4ckcGFIzL49qQz2X3gKK9s2s3LwSuEZ97YCcCwjLRj706+ICedHl2SQ05+YvEUwCBgW73bZcT+ym9qzCBgR5w5ZgK/jnOsiEio0run8tlzB/LZcwfi7mwqP3DsVNMnV5ex+O9biUaMMVm9Y2cXnZ7B6EG9SIpGwo7+D+IpgMYmuBrOG8UzpvE7N/t3oAZYcpz1c4G5AEOGDInnLkVE2oyZMaJfD0b068FNE4ZRVVPHmnf3HLsy2n+/sJGf/99GenRJYkJOxrH3HwxJD3+6KJ4CKAOy6t0eDLx3EmM+xsxuIHaA+JN+nIMR7j4fmA+xYwBx5BURCU1KUoSC4ekUDE/nm1fAnoNVvLKpIrh2cgV/WhebLhqa3i34ILtMLshJp1fXtp8uiqcAVgG5ZjYM2A5MBaY1GLMMuC04PpAPVLr7Cad/gjOLvg18wt0PNTu5iEgH0Ccthc+MHshnRsemi7ZUHDw2XfS7V7ezpOhdIgbnfThdlJvBuVm9SW6D6aJ4TwO9EvhvYqeBLnL3u8zsZgB3nxecBvoLYBKx00BvcvfiYNvHgUuADOB94E53X2hmpUAqsDv4MYXufvOJcugsIBHpTKpr63j13b2sKClneUkFr5ftpc6hR2oSBTnpXBy8GS07vdspXftAbwQTEWnnKg9Vx96MFny6admewwAM7tOVH39pNBfmZJzU/erD4ERE2rle3ZKZfM4AJp8zAHdn6+5DsfcebCxnQK+W/0gKFYCISDtkZmRnpJGdkcb0gqGt8jPa10mpIiLSZlQAIiIJSgUgIpKgVAAiIglKBSAikqBUACIiCUoFICKSoFQAIiIJqkN9FISZlQNbT3LzDKCiBeO0FOVqHuVqHuVqnvaaC04t21B3z2y4sEMVwKkws+LGPgsjbMrVPMrVPMrVPO01F7RONk0BiYgkKBWAiEiCSqQCmB92gONQruZRruZRruZpr7mgFbIlzDEAERH5R4n0CkBEROpRAYiIJKhOUQBmNsnM3jazUjO7o5H1Zmb3B+tfN7Ox8W7byrmuC/K8bmZ/M7Nz6617x8zeMLPXzKxFr4MZR65LzKwy+Nmvmdn34922lXN9s16mN82s1sz6Buta5fdlZovMbJeZvXmc9WHtW03lCmvfaipXWPtWU7nafN8K7jvLzF4ys/Vmts7MvtrImNbbx9y9Q38Ru1D9JmA4kAKsBUY2GHMl8CxgQAFQFO+2rZzrQqBP8P3kD3MFt98BMkL6fV0C/OFktm3NXA3GXw282Aa/r4uBscCbx1nf5vtWnLnafN+KM1eb71vx5Apj3wruewAwNvi+B7CxLZ+/OsMrgPFAqbtvdvcq4AlgSoMxU4DFHlMI9DazAXFu22q53P1v7r4nuFkIDG6hn31KuVpp25a+72uBx1voZx+Xuy8HPjjBkDD2rSZzhbRvxfP7Op5Qf18NtMm+BeDuO9x9TfD9fmA9MKjBsFbbxzpDAQwCttW7XcbHf4HHGxPPtq2Zq75ZxFr+Qw782cxWm9ncFsrUnFwXmNlaM3vWzEY1c9vWzIWZdQMmAb+tt7i1fl9NCWPfaq622rfi1db7VtzC3LfMLBsYAxQ1WNVq+1hnuCi8NbKs4bmtxxsTz7YnK+77NrNLif0nvaje4gnu/p6Z9QOeN7MNwV8xbZFrDbHPDjlgZlcCvwNy49y2NXN96GrgFXev/xdda/2+mhLGvhW3Nt634hHGvtUcoexbZtadWOnc7u77Gq5uZJMW2cc6wyuAMiCr3u3BwHtxjoln29bMhZmNBh4Cprj77g+Xu/t7wb+7gKeJvdxrk1zuvs/dDwTfPwMkm1lGPNu2Zq56ptLgJXor/r6aEsa+FZcQ9q0mhbRvNUeb71tmlkzsyX+Juz/VyJDW28da48BGW34RexWzGRjGRwdCRjUYcxX/eBBlZbzbtnKuIUApcGGD5WlAj3rf/w2Y1Ia5TuOjNwmOB94Nfneh/r6Ccb2IzeWmtcXvK7jPbI5/ULPN9604c7X5vhVnrjbft+LJFeK+ZcBi4L9PMKbV9rEOPwXk7jVmdhvwHLGj4ovcfZ2Z3Rysnwc8Q+xIeilwCLjpRNu2Ya7vA+nA/5oZQI3HPu2vP/B0sCwJWOruf2rDXF8CvmxmNcBhYKrH9riwf18Anwf+7O4H623ear8vM3uc2JkrGWZWBtwJJNfL1Ob7Vpy52nzfijNXm+9bceaCNt63AhOA6cAbZvZasOy7xAq81fcxfRSEiEiC6gzHAERE5CSoAEREEpQKQEQkQakAREQSlApARCRBqQBERBKUCkBEJEH9f6GihSUOc4WHAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Test with 1 sample\n",
    "sample = np.genfromtxt(filenames_train[0], delimiter=',')\n",
    "features = extract_features(sample, max_measurements, scale=raw_scale)\n",
    "print(features.shape)\n",
    "print(features)\n",
    "plt.plot(features)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function: loop through filenames, creating feature sets\n",
    "def create_feature_set(filenames):\n",
    "    x_out = []\n",
    "    for file in filenames:\n",
    "        sample = np.genfromtxt(file, delimiter=',')\n",
    "        features = extract_features(sample, max_measurements, raw_scale)\n",
    "        x_out.append(features)\n",
    "        \n",
    "    return np.array(x_out)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Extracted features from training set. Shape: (181, 3)\n",
      "Extracted features from validation set. Shape: (59, 3)\n",
      "Extracted features from test set. Shape: (59, 3)\n"
     ]
    }
   ],
   "source": [
    "# Create training, validation, and test sets\n",
    "x_train = create_feature_set(filenames_train)\n",
    "print('Extracted features from training set. Shape:', x_train.shape)\n",
    "x_val = create_feature_set(filenames_val)\n",
    "print('Extracted features from validation set. Shape:', x_val.shape)\n",
    "x_test = create_feature_set(filenames_test)\n",
    "print('Extracted features from test set. Shape:', x_test.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.00769517 0.00414957 0.00490948]\n",
      "[[1.43478461e-05 7.25720447e-06 1.87261642e-06]\n",
      " [7.25720447e-06 1.41098698e-05 3.71337361e-06]\n",
      " [1.87261642e-06 3.71337361e-06 2.22039820e-06]]\n"
     ]
    }
   ],
   "source": [
    "# Use training data (normal) to calculate mean and covariance matrices\n",
    "model_mu = np.mean(x_train, axis=0)\n",
    "model_cov = np.cov(x_train.T)\n",
    "print(model_mu)\n",
    "print(model_cov)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate mahalanobis distance of x from group described by mu, cov\n",
    "# Based on: https://www.machinelearningplus.com/statistics/mahalanobis-distance/\n",
    "def mahalanobis(x, mu, cov):\n",
    "    x_minus_mu = x - mu\n",
    "    inv_covmat = sp.linalg.inv(cov)\n",
    "    left_term = np.dot(x_minus_mu, inv_covmat)\n",
    "    mahal = np.dot(left_term, x_minus_mu.T)\n",
    "    if mahal.shape == ():\n",
    "        return mahal\n",
    "    else:\n",
    "        return mahal.diagonal()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Average MD for normal validation set: 2.5604970923201917\n",
      "Standard deviation of MDs for normal validation set: 6.7613861828484465\n",
      "Recommended threshold (3x std dev + avg): 22.84465564086553\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(array([46.,  7.,  3.,  0.,  0.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,\n",
       "         0.,  0.,  0.,  1.,  0.,  0.,  1.]),\n",
       " array([ 0.08767279,  2.13520302,  4.18273325,  6.23026348,  8.27779371,\n",
       "        10.32532394, 12.37285417, 14.4203844 , 16.46791463, 18.51544486,\n",
       "        20.56297509, 22.61050532, 24.65803555, 26.70556578, 28.75309601,\n",
       "        30.80062624, 32.84815647, 34.8956867 , 36.94321693, 38.99074716,\n",
       "        41.03827738]),\n",
       " <a list of 20 Patch objects>)"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8GearUAAALD0lEQVR4nO3dT6ild33H8c+3kxQFKSbNTRgS2+kiFEVqhMEG7EKMgVQTk01AwTILYTYWIlgkuikWCq7ETTeDigOKJaA0iZsyjAZbkNiJxmoYZUTaVDpkRkXUjUX9dnEf9WbmTu7J/f91Xi8Yznl+95w5X37cvHl47n0y1d0BYJ4/OOgBANgeAQcYSsABhhJwgKEEHGCoG/bzw2655ZY+duzYfn4kwHjPPPPMD7t77cr1fQ34sWPHcu7cuf38SIDxquq/N1t3CQVgKAEHGErAAYYScIChBBxgKAEHGErAAYYScIChBBxgqH29E3MnHnhg++998sndmwPgsHAGDjCUgAMMJeAAQwk4wFACDjCUgAMMJeAAQwk4wFACDjCUgAMMJeAAQwk4wFACDjCUgAMMJeAAQwk4wFACDjCUgAMMJeAAQwk4wFACDjCUgAMMJeAAQ60c8Ko6UlXfqKovLsc3V9WZqrqwPN60d2MCcKWXcwb+SJLzG44fTXK2u+9McnY5BmCfrBTwqrojyTuSfGLD8oNJTi/PTyd5aHdHA+ClrHoG/vEkH0zy6w1rt3X3xSRZHm/d7I1VdbKqzlXVucuXL+9oWAB+Z8uAV9X9SS519zPb+YDuPtXdx7v7+Nra2nb+CgA2ccMKr3lzkndW1duTvCLJH1XVZ5K8UFVHu/tiVR1NcmkvBwXgxbY8A+/uD3X3Hd19LMm7knypu9+T5IkkJ5aXnUjy+J5NCcBVdvJ74B9Ncm9VXUhy73IMwD5Z5RLKb3X3U0meWp7/KMk9uz8SAKtwJybAUAIOMJSAAwwl4ABDCTjAUAIOMJSAAwwl4ABDCTjAUAIOMJSAAwwl4ABDCTjAUAIOMJSAAwwl4ABDCTjAUAIOMJSAAwwl4ABDCTjAUAIOMJSAAwwl4ABDCTjAUAIOMJSAAwwl4ABDCTjAUAIOMJSAAwwl4ABDCTjAUAIOMJSAAwwl4ABDbRnwqnpFVX2tqr5ZVc9V1UeW9Zur6kxVXVgeb9r7cQH4jVXOwH+R5K3d/YYkdyW5r6ruTvJokrPdfWeSs8sxAPtky4D3up8vhzcufzrJg0lOL+unkzy0JxMCsKmVroFX1ZGqejbJpSRnuvvpJLd198UkWR5vvcZ7T1bVuao6d/ny5d2aG+C6t1LAu/tX3X1XkjuSvKmqXr/qB3T3qe4+3t3H19bWtjsnAFd4Wb+F0t0/SfJUkvuSvFBVR5Nkeby069MBcE2r/BbKWlW9enn+yiRvS/KdJE8kObG87ESSx/dqSACudsMKrzma5HRVHcl68B/r7i9W1VeTPFZV703yfJKH93BOAK6wZcC7+z+TvHGT9R8luWcvhgJga+7EBBhKwAGGEnCAoQQcYCgBBxhKwAGGEnCAoQQcYCgBBxhKwAGGEnCAoQQcYCgBBxhKwAGGEnCAoQQcYCgBBxhKwAGGEnCAoQQcYCgBBxhKwAGGEnCAoQQcYCgBBxhKwAGGEnCAoQQcYCgBBxhKwAGGEnCAoQQcYCgBBxhKwAGGEnCAoQQcYKgtA15Vr6mqL1fV+ap6rqoeWdZvrqozVXVhebxp78cF4DdWOQP/ZZIPdPdrk9yd5H1V9bokjyY52913Jjm7HAOwT7YMeHdf7O6vL89/luR8ktuTPJjk9PKy00ke2qshAbjay7oGXlXHkrwxydNJbuvui8l65JPceo33nKyqc1V17vLlyzubFoDfWjngVfWqJJ9P8v7u/umq7+vuU919vLuPr62tbWdGADaxUsCr6sasx/uz3f2FZfmFqjq6fP1okkt7MyIAm1nlt1AqySeTnO/uj2340hNJTizPTyR5fPfHA+BabljhNW9O8jdJvlVVzy5rH07y0SSPVdV7kzyf5OG9GRGAzWwZ8O7+9yR1jS/fs7vjALAqd2ICDCXgAEMJOMBQAg4wlIADDCXgAEMJOMBQAg4wlIADDCXgAEMJOMBQAg4wlIADDCXgAEMJOMBQAg4wlIADDCXgAEMJOMBQAg4wlIADDCXgAEMJOMBQAg4wlIADDCXgAEMJOMBQAg4wlIADDCXgAEMJOMBQAg4wlIADDCXgAEMJOMBQAg4w1JYBr6pPVdWlqvr2hrWbq+pMVV1YHm/a2zEBuNIqZ+CfTnLfFWuPJjnb3XcmObscA7CPtgx4d38lyY+vWH4wyenl+ekkD+3yXABsYbvXwG/r7otJsjzeeq0XVtXJqjpXVecuX768zY8D4Ep7/kPM7j7V3ce7+/ja2tpefxzAdWO7AX+hqo4myfJ4afdGAmAV2w34E0lOLM9PJHl8d8YBYFWr/Brh55J8NcmfV9UPquq9ST6a5N6qupDk3uUYgH10w1Yv6O53X+NL9+zyLAC8DO7EBBhKwAGGEnCAoQQcYCgBBxhKwAGGEnCAobb8PfDfBw88sLP3P/nk7swBsJucgQMMJeAAQwk4wFACDjCUgAMMJeAAQwk4wFACDjCUgAMMJeAAQwk4wFACDjCUgAMMJeAAQwk4wFACDjCUgAMMdV38izw7tZN/0ce/5gPsFWfgAEMJOMBQAg4wlIADDCXgAEMJOMBQAg4wlIADDOVGnt9Tbj7iMJv4/bmTmZO9mdsZOMBQOwp4Vd1XVd+tqu9V1aO7NRQAW9t2wKvqSJJ/SvLXSV6X5N1V9brdGgyAl7aTM/A3Jfled3+/u/8vyT8neXB3xgJgKzv5IebtSf5nw/EPkvzllS+qqpNJTi6HP6+q727z825J8sNtvvfAVO37R+54nw5g5oMw8vvpABy6fTrE358vuVc7nPtPN1vcScA3G6evWug+leTUDj5n/cOqznX38Z3+Pb/v7NNq7NNq7NPqDmKvdnIJ5QdJXrPh+I4k/7uzcQBY1U4C/h9J7qyqP6uqP0zyriRP7M5YAGxl25dQuvuXVfW3Sf41yZEkn+ru53Ztsqvt+DLMdcI+rcY+rcY+rW7f96q6r7psDcAA7sQEGErAAYY69AF3u/61VdWnqupSVX17w9rNVXWmqi4sjzcd5IyHQVW9pqq+XFXnq+q5qnpkWbdXG1TVK6rqa1X1zWWfPrKs26dNVNWRqvpGVX1xOd73fTrUAXe7/pY+neS+K9YeTXK2u+9McnY5vt79MskHuvu1Se5O8r7l+8hevdgvkry1u9+Q5K4k91XV3bFP1/JIkvMbjvd9nw51wON2/ZfU3V9J8uMrlh9Mcnp5fjrJQ/s61CHU3Re7++vL859l/T+622OvXqTX/Xw5vHH507FPV6mqO5K8I8knNizv+z4d9oBvdrv+7Qc0yxS3dffFZD1cSW494HkOlao6luSNSZ6OvbrKclng2SSXkpzpbvu0uY8n+WCSX29Y2/d9OuwBX+l2fVhFVb0qyeeTvL+7f3rQ8xxG3f2r7r4r63dWv6mqXn/QMx02VXV/kkvd/cxBz3LYA+52/Zfvhao6miTL46UDnudQqKobsx7vz3b3F5Zle3UN3f2TJE9l/Wcs9unF3pzknVX1X1m/rPvWqvpMDmCfDnvA3a7/8j2R5MTy/ESSxw9wlkOhqirJJ5Oc7+6PbfiSvdqgqtaq6tXL81cmeVuS78Q+vUh3f6i77+juY1lv0pe6+z05gH069HdiVtXbs3696Te36//jAY90aFTV55K8Jev/G8sXkvx9kn9J8liSP0nyfJKHu/vKH3ReV6rqr5L8W5Jv5XfXLD+c9evg9mpRVX+R9R++Hcn6yd1j3f0PVfXHsU+bqqq3JPm77r7/IPbp0AccgM0d9ksoAFyDgAMMJeAAQwk4wFACDjCUgAMMJeAAQ/0/6gOjl8fVqCwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Calculate the mahalanobis distance for each validation sample\n",
    "normal_mds = mahalanobis(x_val, model_mu, model_cov)\n",
    "print('Average MD for normal validation set:', np.average(normal_mds))\n",
    "print('Standard deviation of MDs for normal validation set:', np.std(normal_mds))\n",
    "print('Recommended threshold (3x std dev + avg):', (3*np.std(normal_mds)) + np.average(normal_mds))\n",
    "fig, ax = plt.subplots(1,1)\n",
    "ax.hist(normal_mds, bins=20, label='normal', color='blue', alpha=0.7)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Extracted features from anomaly set. Shape: (500, 3)\n"
     ]
    }
   ],
   "source": [
    "# Extract features from anomaly test set (truncate to length of X test set)\n",
    "anomaly_features = create_feature_set(anomaly_op_filenames[:max_anomaly_samples])\n",
    "print('Extracted features from anomaly set. Shape:', anomaly_features.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Average MD for anomaly set: 81.32351235879358\n",
      "Standard deviation of MDs for validation set: 83.82816246046214\n"
     ]
    }
   ],
   "source": [
    "# Calculate the mahalanobis distance for each anomaly sample\n",
    "anomaly_mds = mahalanobis(anomaly_features, model_mu, model_cov)\n",
    "print('Average MD for anomaly set:', np.average(anomaly_mds))\n",
    "print('Standard deviation of MDs for validation set:', np.std(anomaly_mds))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([148.,  93.,  60.,   2.,   0.,   0.,   1.,   2.,   7.,  38.,  52.,\n",
       "         26.,  28.,  10.,  16.,   7.,   4.,   2.,   2.,   2.]),\n",
       " array([6.93611983e-02, 1.55682001e+01, 3.10670390e+01, 4.65658779e+01,\n",
       "        6.20647168e+01, 7.75635557e+01, 9.30623946e+01, 1.08561234e+02,\n",
       "        1.24060072e+02, 1.39558911e+02, 1.55057750e+02, 1.70556589e+02,\n",
       "        1.86055428e+02, 2.01554267e+02, 2.17053106e+02, 2.32551945e+02,\n",
       "        2.48050784e+02, 2.63549623e+02, 2.79048462e+02, 2.94547300e+02,\n",
       "        3.10046139e+02]),\n",
       " <a list of 20 Patch objects>)"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8GearUAAAPRklEQVR4nO3df6zdd13H8efLVoZAkM3eLrUddpgGHSiB3FSUxJAUpCpdF3WxU0jFxoZkIhoJrmLcX4tLMCoYwVQ220SypZno2iUgTSNZTIByN0C2lbGGxa2srBeIqJhMO9/+cU+Tk7tze88533Puj0+fj6Q55/v5fD/n+95np6/7ybff7/emqpAkteX7VrsASdLkGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ3auNoFAGzatKm2b9++2mVI0rry0EMPfauqZgb1rYlw3759O3Nzc6tdhiStK0n+bak+T8tIUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGrQmbmLqbM+e1a5AWjtOnFjtCrQGuHKXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBy4Z7kruTXEjyyIC+9yWpJJv62g4lOZvk8SRvm3TBkqTlDbNyPwLsXtyY5DrgrcBTfW03APuA1/TGfCTJholUKkka2rLhXlUPAt8Z0PXnwPuB6mvbC9xbVc9V1ZPAWWDnJAqVJA1vrHPuSW4EvlFVX17UtRV4um/7XK9t0GccTDKXZG5+fn6cMiRJSxg53JO8BPgA8MeDuge01YA2qupwVc1W1ezMzMBf3i1JGtM4Dw77UeB64MtJALYBDyfZycJK/bq+fbcBz3QtUpI0mpFX7lX1laraXFXbq2o7C4H+hqr6JnAc2JfkqiTXAzuA0xOtWJK0rGEuhbwH+Czw6iTnkhxYat+qehQ4BjwGfAq4taqen1SxkqThLHtapqpuWaZ/+6LtO4A7upUlSerCO1QlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBg3zC7LvTnIhySN9bR9M8tUk/5rkH5K8oq/vUJKzSR5P8rZpFS5JWtowK/cjwO5FbSeB11bVTwJfAw4BJLkB2Ae8pjfmI0k2TKxaSdJQlg33qnoQ+M6itk9X1cXe5ueAbb33e4F7q+q5qnoSOAvsnGC9kqQhTOKc+28Cn+y93wo83dd3rtf2AkkOJplLMjc/Pz+BMiRJl3QK9yQfAC4CH7/UNGC3GjS2qg5X1WxVzc7MzHQpQ5K0yMZxBybZD7wd2FVVlwL8HHBd327bgGfGL0+SNI6xVu5JdgN/ANxYVf/d13Uc2JfkqiTXAzuA093LlCSNYtmVe5J7gDcDm5KcA25n4eqYq4CTSQA+V1XvrqpHkxwDHmPhdM2tVfX8tIqXJA22bLhX1S0Dmu+6zP53AHd0KUqS1I13qEpSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNWjbck9yd5EKSR/rarklyMskTvder+/oOJTmb5PEkb5tW4ZKkpQ2zcj8C7F7Udhtwqqp2AKd62yS5AdgHvKY35iNJNkysWknSUJYN96p6EPjOoua9wNHe+6PATX3t91bVc1X1JHAW2DmhWiVJQxr3nPu1VXUeoPe6ude+FXi6b79zvbYXSHIwyVySufn5+THLkCQNMul/UM2Athq0Y1UdrqrZqpqdmZmZcBmSdGUbN9yfTbIFoPd6odd+Driub79twDPjlydJGse44X4c2N97vx+4v699X5KrklwP7ABOdytRkjSqjcvtkOQe4M3ApiTngNuBO4FjSQ4ATwE3A1TVo0mOAY8BF4Fbq+r5KdUuSVrCsuFeVbcs0bVrif3vAO7oUpQkqRvvUJWkBhnuktQgw12SGmS4S1KDlv0HVUnrzJ490z/GiRPTP4Y6ceUuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhrUKdyT/F6SR5M8kuSeJC9Ock2Sk0me6L1ePaliJUnDGTvck2wFfgeYrarXAhuAfcBtwKmq2gGc6m1LklZQ19MyG4EfSLIReAnwDLAXONrrPwrc1PEYkqQRjR3uVfUN4E+Bp4DzwHer6tPAtVV1vrfPeWDzoPFJDiaZSzI3Pz8/bhmSpAG6nJa5moVV+vXADwMvTfKOYcdX1eGqmq2q2ZmZmXHLkCQN0OW0zFuAJ6tqvqr+F/gE8DPAs0m2APReL3QvU5I0ii7h/hTwxiQvSRJgF3AGOA7s7+2zH7i/W4mSpFFtHHdgVX0+yX3Aw8BF4IvAYeBlwLEkB1j4AXDzJAqVJA1v7HAHqKrbgdsXNT/HwipekrRKvENVkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJalCncE/yiiT3JflqkjNJfjrJNUlOJnmi93r1pIqVJA2n68r9Q8CnqurHgNcBZ4DbgFNVtQM41duWJK2gscM9ycuBnwXuAqiq/6mqfwf2Akd7ux0FbupapCRpNBs7jH0VMA/8bZLXAQ8B7wWurarzAFV1PsnmQYOTHAQOArzyla/sUIakFbdnz3jjTpyYbB1aUpfTMhuBNwAfrarXA99jhFMwVXW4qmaranZmZqZDGZKkxbqE+zngXFV9vrd9Hwth/2ySLQC91wvdSpQkjWrscK+qbwJPJ3l1r2kX8BhwHNjfa9sP3N+pQknSyLqccwd4D/DxJC8Cvg68i4UfGMeSHACeAm7ueAxJ0og6hXtVfQmYHdC1q8vnSpK68Q5VSWqQ4S5JDTLcJalBhrskNchwl6QGdb0UUpLa1f+YhXX26ARX7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGXRE3MZ0+vdoVaJJ27lztCqS1z5W7JDXIcJekBhnuktQgw12SGtQ53JNsSPLFJA/0tq9JcjLJE73Xq7uXKUkaxSRW7u8FzvRt3wacqqodwKnetiRpBXUK9yTbgF8EPtbXvBc42nt/FLipyzEkSaPrunL/C+D9wP/1tV1bVecBeq+bBw1McjDJXJK5+fn5jmVIkvqNHe5J3g5cqKqHxhlfVYeraraqZmdmZsYtQ5I0QJc7VN8E3JjkF4AXAy9P8nfAs0m2VNX5JFuAC5MoVJI0vLHDvaoOAYcAkrwZeF9VvSPJB4H9wJ291/snUKckTV//70xd56ZxnfudwFuTPAG8tbctSVpBE3lwWFV9BvhM7/23gV2T+FxJ0ni8Q1WSGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWrQRJ4tI0lXlP6nR544sXp1XIYrd0lqkOEuSQ3ytIwkDWOd/SIPV+6S1CDDXZIaZLhLUoMMd0lq0NjhnuS6JP+c5EySR5O8t9d+TZKTSZ7ovV49uXIlScPosnK/CPx+Vf048Ebg1iQ3ALcBp6pqB3Cqty1JWkFjh3tVna+qh3vv/xM4A2wF9gJHe7sdBW7qWqQkaTQTOeeeZDvweuDzwLVVdR4WfgAAm5cYczDJXJK5+fn5SZQhSerpHO5JXgb8PfC7VfUfw46rqsNVNVtVszMzM13LkCT16RTuSb6fhWD/eFV9otf8bJItvf4twIVuJUqSRtXlapkAdwFnqurP+rqOA/t77/cD949fniRpHF2eLfMm4J3AV5J8qdf2h8CdwLEkB4CngJu7lShJGtXY4V5V/wJkie5d436uJKk771CVpEnZs2fNPD3ScJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUFd7lCVJK2R69oXc+UuSQ0y3CWpQYa7JDXIc+5ad06fnu7n79w53c/XGjKt8+X9n3vixHSOsQxX7pLUIMNdkhpkuEtSgwx3SZqmVXrGu+EuSQ3yahlpkXGvxvEqG60lUwv3JLuBDwEbgI9V1Z3TOpakdWKU0xOTvoRwtR8TsNTxp3Sp5FROyyTZAPwV8PPADcAtSW6YxrEkSS80rZX7TuBsVX0dIMm9wF7gsSkdT2rKUqeG1tqpn0t1Lq5rqfYVsdor9DViWuG+FXi6b/sc8FP9OyQ5CBzsbf5Xksf7un8Q+O6Azx2lfXHb4u1NwLeWqH9Slqp3kmOX2+9y/cPM26A253JQ3wNDfTeHmd+l5/OBy1Q3mnHnc/C4wXUtNSfL1bDQngyz//r/biZdvps/suSIqpr4H+BmFs6zX9p+J/CXI4w/3LV9cduA7blp/LcPU+8kxy633+X6h5k353L8uRxm7tbbfI4ybtz59O/5ZI4xrUshzwHX9W1vA54ZYfxS/8IwSvvittV4wEOXYw47drn9Ltc/zLwNanMuh+8bZu7W03yOMm7c+fTv+QSOkd5PhYlKshH4GrAL+AbwBeDXqurRiR9sTEnmqmp2tetogXM5Wc7n5FzJczmVc+5VdTHJbwP/xMKlkHevpWDvObzaBTTEuZws53Nyrti5nMrKXZK0unz8gCQ1yHCXpAYZ7pLUIMN9gCSvSnJXkvtWu5b1KMlLkxxN8jdJfn2161nP/C5OVpKbet/L+5P83GrXM03NhXuSu5NcSPLIovbdSR5PcjbJbZf7jKr6elUdmG6l68uI8/pLwH1V9VvAjSte7Bo3ylz6XVzeiPP5j73v5W8Av7oK5a6Y5sIdOALs7m9Y6kFmSX4iyQOL/mxe+ZLXhSMMOa8s3LR26fETz69gjevFEYafSy3vCKPP5x/1+pvV3PPcq+rBJNsXNQ98kFlV/Qnw9pWtcH0aZV5ZuEN5G/Al2lxAdDLiXPqwvWWMMp9JzgB3Ap+sqodXtNAVdqX8xRv0ILOtS+2c5IeS/DXw+iSHpl3cOrbUvH4C+OUkH2V1bgdfjwbOpd/FsS313XwP8BbgV5K8ezUKWynNrdyX8ILHywFL3r1VVd8Gmv4fPyED57Wqvge8a6WLWeeWmku/i+NZaj4/DHx4pYtZDVfKyr3rg8w0mPM6Oc7lZF3x83mlhPsXgB1Jrk/yImAfcHyVa2qB8zo5zuVkXfHz2Vy4J7kH+Czw6iTnkhyoqovApQeZnQGOrcEHma1pzuvkOJeT5XwO5oPDJKlBza3cJUmGuyQ1yXCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalB/w9XxfcxWU4otQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot histograms of normal validation vs. anomaly sets (MSEs)\n",
    "fig, ax = plt.subplots(1,1)\n",
    "plt.xscale(\"log\")\n",
    "ax.hist(normal_mds, bins=20, label='normal', color='blue', alpha=0.7)\n",
    "ax.hist(anomaly_mds, bins=20, label='anomaly', color='red', alpha=0.7)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a classifier (0 = normal, 1 = anomaly)\n",
    "def detect_anomaly(x, model_mu, model_cov, threshold=0):\n",
    "    md = mahalanobis(x, model_mu, model_cov)\n",
    "    if md > threshold:\n",
    "        return 1\n",
    "    else:\n",
    "        return 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Choose a threshold\n",
    "anomaly_threshold = 3.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0]\n"
     ]
    }
   ],
   "source": [
    "# Perform classification on test set\n",
    "pred_test = [detect_anomaly(x, model_mu, model_cov, anomaly_threshold) for x in x_test]\n",
    "print(pred_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1, 1, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 0, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 0, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 0, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 0, 0, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 0, 0, 1, 0, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1]\n"
     ]
    }
   ],
   "source": [
    "# Perform classification on anomaly set\n",
    "pred_anomaly = [detect_anomaly(x, model_mu, model_cov, anomaly_threshold) for x in anomaly_features]\n",
    "print(pred_anomaly)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Combine predictions into one long list and create a label list\n",
    "pred = np.array(pred_test + pred_anomaly)\n",
    "labels = ([0] * len(pred_test)) + ([1] * len(pred_anomaly))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 46  13]\n",
      " [102 398]]\n"
     ]
    }
   ],
   "source": [
    "# Create confusion matrix\n",
    "cm = confusion_matrix(labels, pred)\n",
    "print(cm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(33.0, 0.5, 'Actual')"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAEWCAYAAABG030jAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8GearUAAAgAElEQVR4nO3dd5xU1fnH8c9XIagIKKKILaCxxIo/S+w9Go2JmkRDogZbUGM3RtTYUEnsRo3RgA0bilFiS2zE3hAQUOyxIgiCCogG2d3n98e9uw64DDPL3J258H3zuq+dOXPvOWdgePbMc889VxGBmZnlx2LV7oCZmZXHgdvMLGccuM3McsaB28wsZxy4zcxyxoHbzCxnHLhtgUlaUtJ9kqZJunMB6tlf0sOV7Fs1SPq3pN7V7octvBy4FyGSfi1phKQvJE1MA8w2Faj6F0BXYLmI2LellUTErRGxawX6MwdJO0gKSXfPVb5RWv54ifWcLemW+e0XEbtHxKAWdtdsvhy4FxGSTgT+AvyJJMiuBvwN2KsC1X8XeDMi6ipQV1Y+AbaStFxBWW/gzUo1oIT/T1nm/CFbBEjqBJwDHBURd0fEzIiYHRH3RcQf0n3aSfqLpAnp9hdJ7dLXdpA0XtLvJU1OR+sHp6/1A84EfpmO5A+de2QqqXs6sm2TPj9I0juSZkh6V9L+BeVPFxy3laQX0xTMi5K2KnjtcUnnSnomredhSV2K/DV8DfwT6JUevziwH3DrXH9Xl0v6UNJ0SSMlbZuW/wg4reB9jinoR39JzwBfAqunZYelr18t6R8F9V8gaZgklfwPaDYXB+5Fw5bAEsDQIvv8EdgC6AlsBGwOnF7w+opAJ2Bl4FDgKknLRsRZJKP4OyJi6Yi4rlhHJLUHrgB2j4gOwFbA6Gb26ww8kO67HHAp8MBcI+ZfAwcDKwDfAU4q1jZwE/Cb9PFuwDhgwlz7vEjyd9AZuA24U9ISEfHgXO9zo4JjDgT6AB2A9+eq7/fAhukvpW1J/u56h9easAXgwL1oWA6YMp9Uxv7AORExOSI+AfqRBKRGs9PXZ0fEv4AvgLVb2J8GYH1JS0bExIgY18w+PwbeioibI6IuIgYDrwM/Kdjnhoh4MyK+AoaQBNx5iohngc6S1iYJ4Dc1s88tETE1bfMSoB3zf583RsS49JjZc9X3JXAAyS+eW4BjImL8fOozK8qBe9EwFejSmKqYh5WYc7T4flrWVMdcgf9LYOlyOxIRM4FfAkcAEyU9IGmdEvrT2KeVC55/3IL+3AwcDexIM99A0nTQa2l65nOSbxnFUjAAHxZ7MSKGA+8AIvkFY7ZAHLgXDc8B/wP2LrLPBJKTjI1W49tphFLNBJYqeL5i4YsR8VBE/BDoRjKKHlhCfxr79FEL+9ToZuB3wL/S0XCTNJXRlyT3vWxELANMIwm4APNKbxRNe0g6imTkPgE4ueVdN0s4cC8CImIayQnEqyTtLWkpSW0l7S7pwnS3wcDpkpZPT/KdSfLVviVGA9tJWi09MXpq4wuSukr6aZrrnkWScqlvpo5/AWulUxjbSPolsC5wfwv7BEBEvAtsT5LTn1sHoI5kBkobSWcCHQtenwR0L2fmiKS1gPNI0iUHAidLKprSMZsfB+5FRERcCpxIcsLxE5Kv90eTzLSAJLiMAMYCLwOj0rKWtPUIcEda10jmDLaLkZywmwB8ShJEf9dMHVOBPdN9p5KMVPeMiCkt6dNcdT8dEc19m3gI+DfJFMH3Sb6lFKZBGi8umipp1PzaSVNTtwAXRMSYiHiLZGbKzY0zdsxaQj65bWaWLx5xm5nljAO3mVnOOHCbmeWMA7eZWc4UuyCjqnost5HPmtq3fDhjgSeV2EKo7uuPFnjtl9lT3ik55rTtsnpV15rxiNvMLGdqdsRtZtaqGpq7Dqw2OXCbmQHU1/Jy8nNy4DYzAyIaqt2Fkjlwm5kBNDhwm5nli0fcZmY545OTZmY54xG3mVm+hGeVmJnljE9OmpnljFMlZmY545OTZmY54xG3mVnO+OSkmVnO+OSkmVm+RDjHbWaWL85xm5nljFMlZmY54xG3mVnO1M+udg9K5sBtZgZOlZiZ5Y5TJWZmOZOjEfdi1e6AmVlNaGgofStC0hKShksaI2mcpH5p+dmSPpI0Ot32KDjmVElvS3pD0m7z66pH3GZmQFTu5OQsYKeI+EJSW+BpSf9OX7ssIi4u3FnSukAvYD1gJeBRSWtFkSuCPOI2M4Mkx13qVqyaxBfp07bpFkUO2Qu4PSJmRcS7wNvA5sXacOA2M4OyUiWS+kgaUbD1KaxK0uKSRgOTgUci4oX0paMljZV0vaRl07KVgQ8LDh+fls2TA7eZGZQ14o6IARGxacE2YI6qIuojoiewCrC5pPWBq4E1gJ7AROCSdHc115tiXXXgNjODip2cLBQRnwOPAz+KiElpQG8ABvJNOmQ8sGrBYasAE4rV68BtZgYVy3FLWl7SMunjJYFdgNcldSvYbR/glfTxvUAvSe0k9QDWBIYXa8OzSszMAOoqdiOFbsAgSYuTDI6HRMT9km6W1JMkDfIecDhARIyTNAR4FagDjio2owQcuM3MEhW6cjIixgIbN1N+YJFj+gP9S23DgdvMDHJ15aQDt5kZeK0SM7Pc8YjbzCxnPOI2M8uZys0qyZwDt5kZQBS9WLGmOHCbmYFz3GZmuePAbWaWMz45aWaWM/VFrzKvKQ7cZmbgVImZWe44cJuZ5Yxz3GZm+RINnsdtZpYvTpWYmeWMZ5WYmeWMR9xmZjmTo8DtmwXXoMUWW4z7H7uDa2+7sqms929/xbAX7uGhZ+7mlLOOr2LvrBoGDriECePHMPqlYU1l/c7+A6NGPsKIFx/m3w/cRrduXavYw4VAROlblTlw16CDD9+ft998p+n5Fttsxi6778Du2/6C3bb+GQOvuqmKvbNquOmmIfx4z/3nKLv4kqv5v01+yKab7coD/3qU0/94QpV6t5BoaCh9qzIH7hqz4korsOOu23LHLUObyg44eF+uufx6vv56NgBTp3xare5ZlTz19At8+tnnc5TNmPFF0+P27ZciamAkmGsNUfpWhKQlJA2XNEbSOEn90vLOkh6R9Fb6c9mCY06V9LakNyTtNr+uZpLjltS52OsR4cgzD2f2P5nzz76M9ku3byrrscZ32WyL/+OkPx7DrP/N4k9nXcrYl8ZVsZdWK849py8H7P8Lpk2fzi4/3Lfa3cm3ys0qmQXsFBFfSGoLPC3p38DPgGERcb6kU4BTgL6S1gV6AesBKwGPSlorIubZoaxG3COBEenPubcR8zpIUh9JIySNmPG/qRl1rXbttOt2TJnyKa+MeW2O8sXbtKHTMh3ZZ9cD+PPZl/HX6y6qUg+t1pxx5gX0WGMzBg8eylG/O7ja3cm1aGgoeStaT6Lx61DbdAtgL2BQWj4I2Dt9vBdwe0TMioh3gbeBzYu1kcmIOyJ6tPC4AcAAgB7LbbTIfe/b5Ac92eVHO7DjLtvQrl07lu7Qnsuu+RMfT5jEg/cnJ6XGjHqFhoYGOi+3LJ9O/azKPbZaMfj2odx7z030O+eSanclv8q4clJSH6BPQdGANH41vr44yUD1e8BVEfGCpK4RMREgIiZKWiHdfWXg+YK6xqdl85T5dMA0j7MmsERjWUQ8mXW7eXTRuVdw0blXAPCDrTflt0f15oQjTuPXB+3LVttuzgvPjKDHGt+l7XfaOmgb3/teD95++10AfrLnrrzxxn+r3KOcK2OtksJB5jxerwd6SloGGCpp/SLVqbkqirWfaeCWdBhwHLAKMBrYAngO2CnLdhc2d946lAuvPIcHn76L2V/P5qSjzqh2l6yV3XLzVWy/3ZZ06dKZ994ZQb9zLmb33XdirbXWoKGhgQ8++IjfHXVKtbuZbxmsVRIRn0t6HPgRMElSt3S03Q2YnO42Hli14LBVgAnF6lWWZ6IlvQxsBjwfET0lrQP0i4hfzu/YRTFVYvP34Ywp1e6C1aC6rz9qbtRalpln9io55rQ/5/Z5tidpeWB2GrSXBB4GLgC2B6YWnJzsHBEnS1oPuI0kr70SMAxYs9jJyaxTJf+LiP9JQlK7iHhd0toZt2lmVr7KLevaDRiU5rkXA4ZExP2SngOGSDoU+ADYFyAixkkaArwK1AFHFQvakH3gHp/meP4JPCLpM+bzFcDMrCoqlCqJiLHAxs2UTwV2nscx/YH+pbaRaeCOiH3Sh2dLegzoBDyYZZtmZi0xv2l+taS1ZpWsCsxIt/WBUVm3a2ZWFt9IISHpXOAg4B2g8ddZ4FklZlZrHLib7AesERFfZ9yOmdmC8Y0UmrwCLMM38xXNzGqS7zn5jT8DL0l6hWThFQAi4qcZt2tmVh4H7iaDSCaev8w3OW4zs9rjWSVNpkTEFRm3YWa24DzibjJS0p+Be5kzVeLpgGZWWxy4mzRePbRFQZmnA5pZzYl6p0oa16O9NyIuy6oNM7OKydGIO7N7TqaLpHj2iJnlQjREyVu1ZZ0qeVbSX4E7gJmNhc5xm1nNqYGAXKqsA/dW6c9zCsqc4zaz2pOfFHfmqwPumGX9ZmaVEnX5idyZ5bgBJHWSdGnjndslXSKpU5Ztmpm1SEMZW5VlGriB60mWct0v3aYDN2TcpplZ2Xxy8htrRMTPC573kzQ64zbNzMpXAyPpUmU94v5K0jaNTyRtDXyVcZtmZmXL04g768B9JHCVpPckvQf8FTg84zbNzMpXoRy3pFUlPSbpNUnjJB2Xlp8t6SNJo9Ntj4JjTpX0tqQ3JO02v65mnSp5DbgQWINkXe5pwN7A2IzbNTMrS9RVrKo64PcRMUpSB5I1mx5JX7ssIi4u3FnSukAvYD1gJeBRSWsVu9N71oH7HuBzkntMfpRxW2ZmLRYVynFHxERgYvp4hqTXgJWLHLIXcHtEzALelfQ2sDnw3LwOyDpwrxIRP8q4DTOzBVdG4JbUB+hTUDQgIgY0s193ksX2XgC2Bo6W9BtgBMmo/DOSoP58wWHjKR7oM89xPytpg4zbMDNbYNFQxhYxICI2LdiaC9pLA3cBx0fEdOBqkrRxT5IR+SWNuzbXnWJ9zXrEvQ1wkKR3SdbjFhARsWHG7ZqZlaVSqRIASW1JgvatEXE3QERMKnh9IHB/+nQ8sGrB4asAE4rVn3Xg3j3j+s3MKiLqmxv4lk+SgOuA1yLi0oLybmn+G2AfkpupQ3KjmdskXUpycnJNYHixNrJeq+T9LOs3M6uUCo64twYOBF4uuODwNOBXknqSpEHeI50aHRHjJA0BXiWZkXJUsRklkP2I28wsF6KhMiPuiHia5vPW/ypyTH+gf6ltOHCbmVHZHHfWHLjNzICIyoy4W4MDt5kZHnGbmeVOQ4VmlbQGB24zMyp3crI1OHCbmeHAbWaWO1H9ZbZLNs/ALelKilwvHxHHZtIjM7MqWFhG3CNarRdmZlW2UEwHjIhBrdkRM7Nqql+YZpVIWh7oC6wLLNFYHhE7ZdgvM7NWlacRdynrcd9KcguyHkA/ksVRXsywT2ZmrS4aVPJWbaUE7uUi4jpgdkQ8ERGHAFtk3C8zs1YVUfpWbaVMB5yd/pwo6cckC3yvkl2XzMxaXy2MpEtVSuA+T1In4PfAlUBH4IRMe2Vm1srqG7K+k2PlzDdwR0Tj7XWmATtm2x0zs+qohRRIqUqZVXIDzVyIk+a6zcwWCg05mlVSSqrk/oLHS5DcK63ojSzNzPImT9MBS0mV3FX4XNJg4NHMemRmVgULVaqkGWsCq1W6I3PbcunVs27CcujNN4ZWuwu2kKpUqkTSqsBNwIpAAzAgIi6X1Bm4A+hOcj3MfhHxWXrMqcChQD1wbEQ8VKyNUnLcM5gzx/0xyZWUZmYLjQrOKqkDfh8RoyR1AEZKegQ4CBgWEedLOgU4BegraV2gF7AesBLwqKS1it3pvZRUSYcKvBEzs5pWqUxJREwEJqaPZ0h6DVgZ2AvYId1tEPA4ySB4L+D2iJgFvCvpbWBz4Ll5tTHfXzGShpVSZmaWZw2hkjdJfSSNKNj6NFenpO7AxsALQNc0qDcG9xXS3VYGPiw4bHxaNk/F1uNeAlgK6CJpWaAxAdSRZDhvZrbQKGdWSUQMAAYU20fS0sBdwPERMV2aZ/3NvVD0C0CxVMnhwPEkQXpkQeXTgauKVWpmljeVvMm7pLYkQfvWiLg7LZ4kqVtETJTUDZiclo8HVi04fBXmM+V6nqmSiLg8InoAJ0XE6hHRI902ioi/tvgdmZnVoEAlb8UoGVpfB7wWEZcWvHQv0Dt93Bu4p6C8l6R2knqQzNwbXqyNUqYDNkhaJiI+Tzu1LPCriPhbCceameVCXeUuwNkaOBB4WdLotOw04HxgiKRDgQ+AfQEiYpykIcCrJDNSjio2owRKC9y/jYim1EhEfCbpt4ADt5ktNOY3ki65noinaT5vDbDzPI7pD/QvtY1SAvdikhSRXFckaXHgO6U2YGaWB5XMcWetlMD9EMnw/hqSM51HAP/OtFdmZq2sUiPu1lBK4O4L9AGOJBn+vwR0y7JTZmatbaEacUdEg6TngdWBXwKdSaa5mJktNOoXhhG3pLVIrp//FTCVZHEUIsI3UzCzhU6O7lxWdMT9OvAU8JOIeBtAkm9ZZmYLpYYcjbiLrVXyc5KVAB+TNFDSzsx7iouZWa5FGVu1FbtycmhE/BJYh2QVqxOArpKulrRrK/XPzKxVNJSxVdt8VweMiJkRcWtE7ElyDf1oknVkzcwWGg1SyVu1lbVyeER8GhF/j4idsuqQmVk11JexVVtLbl1mZrbQWVhmlZiZLTLyNKvEgdvMjNqYLVIqB24zM5wqMTPLnVqY5lcqB24zM6DeI24zs3zxiNvMLGccuM3McqZyt5zMXllXTpqZLawquVaJpOslTZb0SkHZ2ZI+kjQ63fYoeO1USW9LekPSbvOr3yNuMzMqfin7jcBfgZvmKr8sIi4uLJC0Lsm9D9YDVgIelbRWsTu9e8RtZkYyj7vUbX4i4kng0xKb3gu4PSJmRcS7wNvA5sUOcOA2M6O8VImkPpJGFGx9SmzmaElj01TKsmnZysCHBfuMT8vmyYHbzIzyAndEDIiITQu2ASU0cTWwBtATmAhckpY3N4YvegW+A7eZGdnfASciJkVEfUQ0AAP5Jh0yHli1YNdVgAnF6nLgNjOjsjnu5kjqVvB0H6Bxxsm9QC9J7ST1ANYEhhery7NKzMyo7KwSSYOBHYAuksYDZwE7SOpJMmh/DzgcICLGSRoCvArUAUcVm1ECDtxmZgA0VHBh14j4VTPF1xXZvz/Qv9T6HbjNzPAl72ZmueMbKZiZ5YxH3GZmOVOn/Iy5HbjNzHCqxMwsd5wqMTPLmUpOB8yaA7eZGU6VmJnljlMlZmY5U5+jMbcDt5kZHnGbmeVOeMRtZpYvHnHbHPpcdDQb77Qp06dOo++uxy1wfdv+fEf2OeYXAAy98h88dddjABx1+fH02OB71NfV898xb3HdqVdTX1fhW6BaxdXX1/PLQ49lheW78LeL+i1QXff86xH+Puh2AA7v3Yu99vghAH3PvoBxr79FmzZtWH/dtTjr5GNp28b//QvlaTqgb6TQCp688z9c0Pucso87/fZz6bLK8nOUte+0ND8/fj/O2KsvZ/z0ZH5+/H6079gegGf++SQn7XQ0fXc9ju+0+w479tqlIv23bN1y5z2s3n21so456OiT+WjipDnKpk2fwdU33MbggX9h8MC/cPUNtzFt+gwAfrzrjtw3eCBDb76aWbO+5q77HqxY/xcWWd8Bp5IcuFvB68Nf5YvPZ8xRtsJqK9J30Bn0v/9izryzPyutUfTeoE023L4nLz81hpnTvmDm9Jm8/NQYNtxhYwBGPzaqab//jnmLzt26VO5NWCY+nvwJTz47nJ//ZLemsg/GT+DwE09nv0OO4TdHnsQ7739YpIZvPPPCSLbcbGM6dexAp44d2HKzjXnmhZEAbLfV5khCEht8f20mTZ6SyfvJszqi5K3aMg3cki6WtF6WbeTVYecfyaCzruWPe57Erf1v5ODzDi/puM4rLsfUid/8p/v046l0XnG5OfZZvM3ibPOz7Rnz+Ki5D7cac8Hlf+fE3x2K9M1/xX4XXsFpJxzJkOuv5KSjD+O8i68qqa5Jn0xhxRW++YbWdfkuTPpkzgA9u66O+x4axjY/2LQyb2AhEmX8qbask1yvAwMktQFuAAZHxLR57Zze4r4PwGade/K9pbtn3L3qaLfUEqy1ydoc+7c/NJW1bdcWgO333YndDt4TgBW7r8jJN55B3dd1fPLhJC47/AKkb9/wLmLOD9LB5x3O6y+8yhsvvpbhu7AF9fgzL9B52WVYb501GT5qLABffvkVo19+jRNP/1PTfl/Png3A0Ace5pYh9wDwwUcTOPKkM2jbpi0rr9SVK/58JtFMPJn783LexVexyUbrs0nP9TN6V/nlk5OpiLgWuFbS2sDBwFhJzwADI+KxZvYfAAwA+PV396n+r7WMLLaYmDl9JqftceK3Xnvizv/wxJ3/AZIc9zUnXcGU8Z80vT514hTW3eKb/3SdV1yOV59/pen5z47bj46dO3LZqVdn+A6sEl4a+yqPP/08Tz33IrO+ns3MmV9y6rkX06FDe+4a9O1R9j4/3pV9frwrkOS4+//x96zcrWvT6yuu0IUXXxrb9HzSJ1PYbOMNm57/7fpb+ezzaZz1p9MzfFf5VQsj6VJlnuOWtDiwTrpNAcYAJ0q6Peu2a9VXX3zFJx9O5gd7bNVUttr3u5d07NgnRrPBdj1p37E97Tu2Z4PtejL2idEA7NBrFzbcfmOuPObSb43CrfaccOTBDPvnLTx81yAu6ncKm2+yEZf/+QxW7rYiD/3nKSD5NvX6W++UVN/WP9iEZ4ePYtr0GUybPoNnh49i6x9sAsA/7n2QZ14YyYX9+rLYYj611ZyGMrZqy3TELelS4KfAMOBPEdF4y/kLJL2RZdu15OgrTuT7W65Hh2U7cuXzA7nrstu56rjLOOS8w9n7mF+weNs2PHfv03zw2nvzrWvmtC8YesWdnHvfRQAMvXwIM6d9AcCh/Y9gykef0G/o+QC8+ODzDL1iSGbvy7JxwVknc+7Ff+XvgwZTV1fH7jtvzzprrj7f4zp17MDhB/2KXoclU06POPjXdOrYAYBzL76Sbl1XYP8+ybe8XbbfiiMP2T+7N5FD9RUc7Ei6HtgTmBwR66dlnYE7gO4kd3nfLyI+S187FTiU5Gbzx0bEQ0Xrz3JkJukQ4PaI+LKZ1zoVy3cvzKkSa7lBIy+pdhesBrXtsvq3T/6UqZyYc9v7Q4u2J2k74AvgpoLAfSHwaUScL+kUYNmI6CtpXWAwsDmwEvAosFZEzPMijExG3JL+L304Glhn7hMkETGqWNA2M2ttlcxxR8STkrrPVbwXsEP6eBDwONA3Lb89ImYB70p6mySIPzev+rNKlRQbFgWwU0btmpm1SDm568IZcKkB6eSKYrpGxESAiJgoaYW0fGXg+YL9xqdl85RJ4I6IHbOo18wsK+Vc8l44A64Cmku7FO1M5osVSFofWBdYoqlHETdl3a6ZWTlaYTrgJEnd0tF2N2ByWj4eWLVgv1WACcUqyvrKybOAK9NtR+BCklkmZmY1pT6i5K2F7gV6p497A/cUlPeS1E5SD2BNYHgzxzfJekLnL4CdgY8j4mBgI6Bdxm2amZWtgSh5mx9Jg0lOLq4tabykQ4HzgR9Kegv4YfqciBgHDAFeBR4Ejio2owSyT5V8FRENkuokdST5ajD/CalmZq2skhfWRMSv5vHSzvPYvz/Qv9T6sw7cIyQtAwwERpLMayz6FcDMrBrydMl71muV/C59eI2kB4GOETG22DFmZtWQpxsptMaskg1JLvFskz7/XkTcnXW7ZmblyNP6PlmvVXI9sCEwjm9SSAE4cJtZTan3iLvJFhGxbsZtmJktsDylSrKeDvhcuoCKmVlNi4iSt2rLesQ9iCR4fwzMIrm0MyJiw+KHmZm1rjyNuLMO3NcDBwIvUxvrj5uZNcvTAb/xQUTcm3EbZmYLrJI3Usha5jcLlnQbcB9JqgQATwc0s1rjVMk3liQJ2LsWlHk6oJnVHAfuVLqwlJlZzauF2SKlynpZ11UkDZU0WdIkSXdJWiXLNs3MWqKSqwNmLet53DeQrDW7EsmteO5Ly8zMakqU8afasg7cy0fEDRFRl243Astn3KaZWdnqo6HkrdqyDtxTJB0gafF0OwCYmnGbZmZly9OVk1kH7kOA/YCPgYkkd8Q5JOM2zczKlqccd9azSj7A95g0sxyohdx1qbJe1nV54LcUrMcNEBEedZtZTWmogRRIqbK+AOce4CngUaDozS/NzKqpkiNuSe8BM0jiXl1EbCqpM3AHyUD2PWC/iPisJfVnHbiXioi+GbdhZrbAMpgtsmNETCl4fgowLCLOl3RK+rxF8THrk5P3S9oj4zbMzBZYQ0TJWwvtRbLUNenPvVtaUdaB+ziS4P2VpOmSZkiannGbZmZlK+cCHEl9JI0o2Pp8qzp4WNLIgte6RsREgPTnCi3ta9azSjqkeZ01gSWybMvMbEGUM5KOiAHAgCK7bB0REyStADwi6fUF7V+hrGeVHEYy6l4FGA1sATwL7Jxlu2Zm5arkycmImJD+nCxpKLA5MElSt4iYKKkbMLml9bdGqmQz4P2I2BHYGJhS/BAzs9ZXH/Ulb8VIai+pQ+NjkmWtXyFZt6l3ultvkll3LZL1rJL/RcT/JCGpXUS8LmntjNs0MytbBS9l7woMlQRJjL0tIh6U9CIwRNKhwAfAvi1tIOvAPV7SMsA/SfI8nwETMm7TzKxslbqUPSLeATZqpnwqFUoTZ31ycp/04dmSHgM6AQ9m2aaZWUvUwuJRpcp6xN0kIp5orbbMzMrlS97NzHLGi0yZmeVMLdwgoVQO3GZmOMdtZpY7znGbmeWMR9xmZjlTC7ckK5UDt5kZHnGbmeWOZ5WYmeWMT06ameWMUyVmZjnjKyfNzHLGI24zs5zJU45befots6iS1Ce9x51ZE38uFl1Z37rMKmPuO0ibgT8XiywHbjOznHHgNjPLGQfufHAe05rjz8Uiyicnzcxyxsy/xcgAAAUASURBVCNuM7OcceA2M8sZB+5FgKT3JHWpdj+s+iQ9LmnTavfDFowDd42T5KtbzWwODtytQFJ3Sa9JGihpnKSHJS0pqaek5yWNlTRU0rLp/o9L+pOkJ4Dj0ueXSXoyrWczSXdLekvSeQXt/FPSyLQNX5xRY5r795H0haT+ksakn4Wuafl3JQ1LPxvDJK2Wlt8o6WpJj0l6R9L2kq5PPxc3FrR1taQRaVv9munLoZIuK3j+W0mXZv6XYJUREd4y3oDuQB3QM30+BDgAGAtsn5adA/wlffw48LeC4x8HLkgfHwdMALoB7YDxwHLpa53Tn0sCrxSUvwd0qfbfw6K+NffvAwTwk7T8QuD09PF9QO/08SHAP9PHNwK3AwL2AqYDG5AMwkYWfMYa21o8/fxsWPBZ2hRoD/wXaJuWPwtsUO2/I2+lbR5xt553I2J0+ngksAawTEQ8kZYNArYr2P+OuY6/N/35MjAuIiZGxCzgHWDV9LVjJY0Bnk/L1qzwe7AF09y/z9fA/enrI0l+yQNsCdyWPr4Z2KagnvsiibYvA5Mi4uWIaADGFRy/n6RRwEvAesC6hR2JiJnAf4A9Ja1DEsBfrtD7tIw5f9p6ZhU8rgeWmc/+M+dxfMNcdTUAbSTtAOwCbBkRX0p6HFiixb21iiry7zM7DcKQfC7m9X+y8IKL+X0WegAnAZtFxGdpCqW5z8K1wGnA68AN5b4nqx6PuKtnGvCZpG3T5wcCTxTZf346AZ+lQWEdYIsF7aBVVLn/Ps8CvdLH+wNPl9FWR5Jf/NPSnPnuze0UES+QjPx/DQwuo36rMo+4q6s3cI2kpUhSHgcvQF0PAkdIGgu8QfJ13GpHuf8+xwLXS/oD8AllfDYiYoykl0hSJ+8AzxTZfQhJXvyzUuu36vMl72aLMEn3A5dFxLBq98VK51SJ2SJI0jKS3gS+ctDOH4+4zcxyxiNuM7OcceA2M8sZB24zs5xx4LaKklQvabSkVyTdmU51bGldN0r6Rfr4WknrFtl3B0lbFTw/QtJvWtq2WS1z4LZK+yoiekbE+iSXcx9R+KKkxVtSaUQcFhGvFtllB6ApcEfENRFxU0vaMqt1DtyWpaeA76Wj4cck3Qa8LGlxSRdJejFd/e5wACX+KulVSQ8AKzRWVLiOtKQfSRqVrqg3TFJ3kl8QJ6Sj/W0lnS3ppHT/YqswXiBpuKQ3C65iNatpvnLSMqFkHfHdSa4YBNgcWD8i3k2XNJ0WEZtJagc8I+lhYGNgbZLV7roCrwLXz1Xv8sBAYLu0rs4R8amka4AvIuLidL+dCw67CTgmIp6QdA5wFnB8+lqbiNhc0h5p+S6V/rswqzQHbqu0JSU1roL4FHAdSQpjeES8m5bvCmzYmL8mWcdjTZLVEQdHRD0wQdJ/mql/C+DJxroi4tNinZHUiW+vwnhnwS53pz8LV+Yzq2kO3FZpX0VEz8ICSTDnaociGQE/NNd+ezDnKnjNUQn7lKNxdb1iK/OZ1RTnuK0aHgKOlNQWQNJaktoDTwK90hx4N2DHZo59Dtg+XboUSZ3T8hlAh7l3johKr8JoVnUeYVg1XEuSlhilZDj+CbA3MBTYieQGAW/STICNiE/SHPndkhYDJgM/JLljzD8k7QUcM9dhlVyF0azqvFaJmVnOOFViZpYzDtxmZjnjwG1mljMO3GZmOePAbWaWMw7cZmY548BtZpYz/w8htGLBkE9J7gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Make confusion matrix pretty\n",
    "df_cm = pd.DataFrame(cm, index=['normal', 'anomaly'], columns=['normal', 'anomaly'])\n",
    "plt.figure()\n",
    "sn.heatmap(df_cm, annot=True)\n",
    "plt.title('Confusion Matrix')\n",
    "plt.xlabel('Prediction')\n",
    "plt.ylabel('Actual')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save Mahalanobis Distance model\n",
    "np.savez(md_file_name + '.npz', model_mu=model_mu, model_cov=model_cov)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}