anomaly_detection.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(307L, 2L) (307L, 2L) (307L, 1L)\n[[ 13.04681517  14.74115241]\n [ 13.40852019  13.7632696 ]\n [ 14.19591481  15.85318113]\n [ 14.91470077  16.17425987]\n [ 13.57669961  14.04284944]]\n[[0]\n [0]\n [0]\n [0]\n [0]]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmUXHWd9/H3p5MOgQYCkRDDEoMDSFCHJR3gUfHQQDyA\nbKLw4JFFZaa1nQdQkLiciSK4QCQwo86JooKIIE8YVgEHURoRR6HT7IsMiFHAQFjkIcia5Pv88bv3\nVnWnq7u601XVXfm8zrmn+t6699avTiX1qd/9LVcRgZmZGUBLowtgZmZjh0PBzMwKDgUzMys4FMzM\nrOBQMDOzgkPBzMwKNQsFSZMl3SHpHkkPSPpKtn2qpJskPZI9bl6rMpiZ2fCoVuMUJAloi4iXJLUC\ntwEnA0cAz0fEWZI+D2weEZ+rSSHMzGxYalZTiOSlbLU1WwI4DLgo234RcHitymBmZsMzsZYnlzQB\n6AW2B/4jIm6XND0ilme7PAVMr3BsJ9AJ0NbWNmennXaqZVHNzJpOb2/vsxExbTjH1OzyUZ8XkTYD\nrgJOBG6LiM3KnvtbRAzartDe3h5Lly6tcSnNzJqLpN6IaB/OMXXpfRQRLwDdwAHA05JmAGSPK+pR\nBjMzG1otex9Ny2oISNoQmAf8AbgWOD7b7XjgmlqVwczMhqeWbQozgIuydoUWYElEXCfpd8ASSScA\nfwaOqmEZzMxsGGoWChFxL7DbANufA/ar1euamdnIeUSzmZkVHApmZlZwKJiZWcGhYGZmBYeCmZkV\nHApmZlZwKJiZWcGhYGZmBYeCmZkVHApmZlZwKJiZWcGhYGZmBYeCmZkVHApmZlZwKJiZWcGhYGZm\nBYeCmZkVHApmZlZwKJiZWcGhYGZmBYeCmZkVHApmZlZwKJiZWcGhYGZmBYeCmZkVHApmZlZwKJiZ\nWaFmoSBpW0ndkh6U9ICkk7Ptp0t6UtLd2XJQrcpgZmbDM7GG514FnBoRd0raBOiVdFP23HkRcU4N\nX9vMzEagZqEQEcuB5dnfKyU9BGxdq9czM7N1V5c2BUmzgN2A27NNJ0q6V9IFkjavRxnMzGxoNQ8F\nSRsDVwCfjogXgcXAW4FdSTWJRRWO65S0VNLSZ555ptbFNDMzahwKklpJgXBJRFwJEBFPR8TqiFgD\nfB/YY6BjI+L8iGiPiPZp06bVsph1s3AhdHf33dbdnbabmY0Ftex9JOCHwEMRcW7Z9hllu30AuL9W\nZRhr5s6Fo44qBUN3d1qfO7ex5TIzy9Wy99G7gWOB+yTdnW37IvBhSbsCASwDPlHDMowpHR2wZEkK\ngq4uWLw4rXd0NLpkZmZJLXsf3QZogKduqNVrjgcdHSkQzjwTFixwIJjZ2OIRzXXW3Z1qCAsWpMf+\nbQxmZo3kUKijvA1hyRI444zSpSQHg5mNFQ6FOurp6duGkLcx9PQ0tlxmZjlFRKPLMKT29vZYunRp\no4thZjauSOqNiPbhHOOagpmZFRwKZmZWcCiYmVnBoWBmZgWHgpmZFRwKZmZWcCiYmVnBoWBmZgWH\ngpmZFRwKZmZWcCiYmVnBoWBmZgWHgpmZFRwKZmZWcCiYmVnBoWBmZgWHgpmZFRwKZmZWcCiYmVlh\n4mBPSpoMHAzsDWwFvALcD1wfEQ/UvnhmZlZPFUNB0ldIgXALcDuwApgM7AiclQXGqRFxbx3KaWZm\ndTBYTeGOiPhyhefOlbQlMLMGZTIzswapGAoRcX3/bZJagI0j4sWIWEGqPZiZWZMYsqFZ0qWSNpXU\nRmpPeFDSabUvmpmZ1Vs1vY92jogXgcOBnwPbAccOdZCkbSV1S3pQ0gOSTs62T5V0k6RHssfN1+kd\nmJnZqKkmFFoltZJC4dqIeAOIKo5bRWqI3hnYC/gXSTsDnwd+FRE7AL/K1s3MbAyoJhS+BywD2oBb\nJb0FeHGogyJieUTcmf29EngI2Bo4DLgo2+0iUtiYmdkYUDEUJP0vSYqIb0XE1hFxUEQE8BegYzgv\nImkWsBupa+v0iFiePfUUML3CMZ2Slkpa+swzzwzn5czMbIQGqykcB/RKukzSRyW9GSCSVdW+gKSN\ngSuAT2dtE4UsZAa8FBUR50dEe0S0T5s2rdqXMzOzdTBYl9QuAEk7AQcCP5I0BegG/gv4bUSsHuzk\nWVvEFcAlEXFltvlpSTMiYrmkGbhbq5nZmDFkm0JE/CEizouIA4B9gduAI0mXgiqSJOCHwEMRcW7Z\nU9cCx2d/Hw9cM5KCm5nZ6Bt07qOcpN2B95Au9fw2Ik6s4rB3k7qu3ifp7mzbF4GzgCWSTgD+DBw1\n7FKbmVlNDBkKkr5Eqhnkl38ulHR5RHx1sOMi4jZAFZ7eb1ilNDOzuqimpvARYJeIeBVA0lnA3cCg\noWBmZuNPNeMU/kqaHTW3AfBkbYpjZmaNNNjU2d8mtSH8P+ABSTdl6/OAO+pTPDMzq6fBLh8tzR57\ngavKtt9Ss9KYmVlDDTZO4aJKz5mZWXOqZursgyXdJel5SS9KWilpyLmPzKx+Fi6E7u6+27q703az\n4aimofnfSIPM3hQRm0bEJhGxaY3LZWbDMHcuHHVUKRi6u9P63LmNLZeNP9WEwuPA/dk8RWZNbbz+\n4u7ogCVLUhB86UvpccmStN1sOKoZpzAfuEHSr4HX8o39pq4wawr5L+78CzX/xb1kSaNLNrSODujq\ngjPPhAULHAg2MtXUFL4GvEwaq7BJ2WLWdMbzL+7ubli8OAXC4sVr13jMqlFNTWGriHhHzUtiNkaM\nx1/c5TWajo60jKdAs7GjmprCDZLeV/OSmDVY3p5Q/ov7W9+CT3xi8P3LNar9oaenbwDkNZ6envqX\nxcY3DdV+LGkl6VacrwFvkCa5i3r2QGpvb4+lS5cOvaPZOujuhsMPBwmuyoZrlq/3/8Xd/9d5/3Wz\nRpPUGxHtwzmmmvspbBIRLRGxobuk2nhU7S/6jg44+miIKH3BX311CoSBfnGP5/YHs0oGu0fzrMEO\nVLLNaBfIbKQqffn/8Y9D9+E/6CA491z43vfg5JNTe8Jb3gLf/Gb6kp8/v3KQ5O0Pu+yydpnGQ3dW\nsz4iYsAFuJx0K83jgLcDWwIzSXdfOxP4b2BepeNHc5kzZ06YDeXmmyO22CI99l/P/16woO8+uUWL\nIqSIrq70/Jw5ERBx6KEDnzsi4uyz03H5eTfdNKKtLaKzs/IxZvUELI1hft8O/iTsTOqSegvwMHAX\ncClwDDB5uC820sWhYNUa7Mt/wYL0L37BgtK2s88u7dPVlZ6fMaMUCFLEvHkRra0pAMqPyfd///vT\ncwcfnNYnTkzHtLWVAqmzMx1nVk+jHgpjZXEo2HAM9OXf2RkxZUopLPbYI32pd3aWwuPggyMmTUrH\nbrddOm6HHdL6nDlpv87OdFxra/ry7+qKmDw57dPamoKkpSWt77BDen6DDdJrd3am8CgPh5tvdlhY\n7TgUbL3Xv6aQfxFPmZIu79x8c1qfODH96+/qStva2tJ6HggQseOOqaaQf8kfemj64s/3a22NOPbY\n0npLSwqA/O98+4QJEXvtFbHnnul8ixaVQiivTdjgymt0OQfq0BwKtl4bqE1hypT0Rb1oUQqItraI\nDTdMj3vtlf4HvPnNpS/wefPSfptumta32CIdmz9fvuTh0dLSNwRaWtKXf76eB1AeQnmZ8oCwoQ3W\nXmSVjSQUhhynMBZ4nIJVY+HC1KOovEtodzdcdhlceSUceCBcfHHavmABPP00/OAHsGZN2tZS1hcv\n3waw9dbwZL8b0La09N1nwoQ0nmHVqsr7SWk/SPsdeyz8+Mcjf7/rm7zXWFdXGlzo7r9Dq8k4BUm/\nqmabWaPNn5/GE5R3S+3pgbe9LXUXvfhi2GCD9MX89a/DhRemL2wp7dvamtbLt8HagSD1DYSNNoLV\nq2HiRNh88777lu8XkcJg1SrYe2/4+c89P9FwlHf/7epyINTKYOMUJkuaCmwhaXNJU7NlFrB1vQpo\nVq2FC9MX8yGHpDEHAMuWwamnpi/fCRPSEpG+xN94A6ZNg0mT0r6vvVYKg0oV6Dlz1n7u5ZfT46uv\nwt/+Vrl8+blbWqC3F/baq+/4CRucJ/yrj4qXjySdDHwa2Ar4a9lTLwLfj4jv1L54iS8f2UDyy0U9\nPaWBaB/4AMyeDb//PeywAzz6aAqKN95IX8r9f+W3tqbnqtX/slFeu1i9urpjJ02CLbeEv/wl1TDO\nPDPVHObPr74M6yNPKTIyo3r5KCL+PSK2Az4bEduVLbvUMxDMKsnvfTBxYnq86670K/6++1KN4JFH\noK0tfXFvtVV6rvwLHYYXCND3eCkdX00g5Me+/e0pEGbPhmOOcSBUyxP+1U81E+IdN9D2iKhbE5lr\nClZJ/ovxwAPhJz9JX7RLlqRLQdtuC48/nkIhv55fb1tsAc8+23dbayt87GNpnqX+jeI9PQ4JGz01\naWgG5pYtewOnA4cOu3RmQxjJVNQ9PaVeRdOnp8fXXkuXjh5/HDbdtO+lo3rbay+YPLnvttZW+Otf\nfU9lG5uqmSX1xLLln4HdgY2HOk7SBZJWSLq/bNvpkp6UdHe2HLRuxbdmMpKbz0+cmGoIc+bAU0+V\ntj/2WLqE9OKL6ZLNqlWVG49r6cYb17689MYbcOut8IUveIZVG3uqqSn093dguyr2+xFwwADbz4uI\nXbPlhhG8vjWp4U5F3d0N3/gGfPKTcOedqVYAqUawenVaNtwQHnyw/rWE/PVWrSq1W8yalR7feAM+\n8pH0nLtY2lhTzTiFn0m6NluuJ02Md9VQx0XErcDzo1BGW48Mpy96VxcceWT6sj3mmFQrgFKNYIst\n4JVXBm5gXleVQqa8S+vkySmoWlrSJaNly9KAtba20pTe7mJpY00192g+p+zvVcCfI+KJdXjNE7PG\n66XAqRExYM9uSZ1AJ8DMmTPX4eVsPOnfF72jo9TlNH/MuyRuv33aZ/ZseOaZdAmpt7d0rmefTZeQ\nqu0dNBwRsNlm8MILa2/v6kp/L16cAiAf6XzssWnA2nHHwQUXwPLlpTu6+Z7KNlZU06bwa1LtYAow\nlRQMI7UYeCuwK7AcWDTI654fEe0R0T5t2rR1eEkbq/o3LOe3wzziCDjjjNKlpLzLaf547rnpcd99\nU7//hx5KX/y9vWmfcrUIhFz/QMjdeSdcfnkKqeeeS2U655wUCAceCN/9Luy/f99bfLqLpY0V1Vw+\n+ifgDuAI4EPA7yV9fCQvFhFPR8TqiFgDfB/YYyTnsebQv2H5ssvSr+qjj07r+RflqlXp8RvfSF+q\nn/0szJyZguOjH4U3vzmNJM4HkrW2rh0O9dLSArffDrvtlsJh993T1Bq77ZZqEBdfnC51vfe9a9cI\n8ju8mTXUUDPmkWoJbypbfxPwcDWz7QGzgPvL1meU/f0Z4LJqzuNZUpvXUHdEK5ffJ2HvvdPjhhuW\nbnRTPvU1RMyeHbHZZgPPblqLpXxW1Pzv3Xcvvcd86u5q3qfZaGEEs6RW0/voOWBl2frKbNugJP0U\n+B3wNklPSDoBWCjpPkn3Ah1ZMNh6rNqG5byt4dhj4bbb0mNE2rbjjvD3v6eeRgBTpqRLSq+8Ur/3\nEVGaATUi1RgefbRUC4pINaDyy2JuWLYxaajUAH5Mug3n6cCXgTtJ3U1PAU4ZbgqNZHFNoXlVU1PI\n98nvh5w/7rRT6V4Fc+akX+hTp6b1/H4I5b/g67VMnJhePy9nZ6dvEGONQY1qCn8ErgbyoT/XAH8C\nNskWsxEpn9RssF/Q+bw3edvCKaekx/e+N7Ud7LBDun6///7w/POwzTape2pbW2MGrH3hC3D11aXy\n/sM/uP3Axg/fZMcaptJNcaqZ/6c8UHp60hiA7343Nez29qZgeGJdOk6PUN7t1F1LbSyo1U12dpR0\nvqRfSLo5X0ZeTLN1Vz5r5ty5qQvoOeekaanf9a4UCJv0q8e2jGT8fpn+PZrK748wYULq9XTttaXp\nK9xmYONRNf9NLie1KfwrcFrZYrZORjLXUT62Yf780i/xyy5LYxtOOQVOOw3uuSeNEVi5Mt1KU0pB\nsWbN4MEweXLqPlrJpEnpPP2PmTQplevGG9Plqocf9pgDG7+qCYVVEbE4Iu6IiN58qXnJrOkNd64j\nGDhIrryyNLahpye1T9x7L8yblwaYffKT6fp+foOcPCj6e/XVym0Qs2enOYvuuCOtb7VVCpD994cb\nbkjn7+hIbQl5G4LbDGxcqtQCTRq9PJXU6+hTwIyybVOH26K9Lot7HzW3fPzBggXV7d/Zmfr9l/dY\nynvz5D2VurpSz6OurrQ+e3Z6jR13LPUSammp3INo5sz0OGFCad/W1vT3okWpHPlrecyBjVWMoPfR\nYOM+e0k9jvLfVOWXjII0XYXZOhlorqOhGmiPPjqNDD7zzHQc9G10zmsb22+fRj/PmpXGLXR1pWkv\nHn20dClp4kR4/fV0Cej119O5Ntgg3Ysh37+nJ93Vbc2a1JB8yilpv/KpKdyobE1juCnSiMU1hebU\n/5d2tb+88xHCG24YsdFGaUxApWPy0c/vfGdaz2sZu+2WagEbbRQxb97aNYUZM0pjIqp9LbOxhlqM\nU5B0xADLfpK2rENmWRMbyX13u7vhkENSG8RnPwsvv5yu5991V987tC1cCJ/6VBr9vPfecP/9cNhh\ncMklaSK6O++EE05IbQu/+U06Ju9dNHlymsF0+vTUPhEB118P112X9j/8cPcssiY2VGoA15Pui3BF\ntjwH/AJ4BDh2uCk0ksU1BcudfXbp13s+l9BGG0VMntz3F3w+J1JXV9/1Qw8t7XPzzekc226btkOp\n3SBve9hzz77nvfnmVNvwaGQbD6jRiOaJwOyI+GBEfBDYmdSmsCfwudGPKbPK5s9PM45GlHoQtbam\nNoFyy5al+yNfemmqVVx+eWojeOqpUo2ipyf1FvrLX1LPoq6uNH/SfvulkdFdXTB1at/2go4O+N73\n3LPImlc1obBtRDxdtr4i2/Y88EZtimVWWf5lftJJqbH5pJPSevllpxtugK9/PTUe55PtHXlkunfz\n3Llrj6Y+7bR0aenoo+GXv0yXsS6/PG03W59UEwq3SLpO0vGSjifNfXSLpDagwm1GzGon/5Ve3mup\nfHu5SZPSr/9Fi1JbQPko6Gru5eABaLbeGer6EqlL6oeA87LlQ2RzJtVrcZuClaum11L5tnwcxEYb\nDbyP73FgzYpatClk5/7PiPhMtvxn9mJmDVFNr6V8HyjVKCZOTDWCXLX3cjBbnww5S6qklZSmzZ4E\ntAJ/j4hNa1y2gmdJtZEon0m1o6PyeldXCg7PbGrNZiSzpA55J9uIKOaalCTgMGCv4RfPrL6GqlGU\nB0RHR3VzL5k1uxHdT0HSXRGxWw3KMyDXFGy0rcu9HMzGi5rUFCQdUbbaArQDrw6zbGZjykBf/NXM\nu2TW7IYMBeCQsr9XActIl5DMzKzJVNOm8LF6FMTMzBqvmgnxtpF0laQV2XKFpG3qUTgzM6uvakY0\nXwhcC2yVLT/LtpmZWZOpJhSmRcSFEbEqW34ETKtxuczMrAGqCYXnJB0jaUK2HEOaPtvMxrmFC9e+\nN0R3d997U9j6pZpQ+DhwFPAUsJw095Ebn82aQP+JAfNR3nPnNrZc1jiD9j6SNAE4IiIOrVN5zKyO\n8lHenu7DcoPWFCJiNfDhkZxY0gVZb6X7y7ZNlXSTpEeyx81Hcm4zGz2eGNDKVXP56LeSviNpb0m7\n50sVx/0IOKDfts8Dv4qIHYBfZetm1kDd3X3vTeH7T6/fqhnRvGv2eEbZtgD2HeygiLhV0qx+mw8D\n9sn+vgi4Bd/S06xh+s8c64kBrZoRzaP5T2N6RCzP/n4KmF5pR0mdQCfAzJkzR7EIZpYbbCZZh8L6\nqZr7KWwAfBCYRVmIRMQZlY4pO3YWcF1EvCNbfyEiNit7/m8RMWS7gmdJNTMbvpHMklpNm8I1pMs+\nq4C/ly0j8bSkGQDZ44oRnsesKu6HbzY81bQpbBMR/RuMR+pa4HjgrOzxmlE6r9mA8n74A919zczW\nVk1N4b8lvXO4J5b0U+B3wNskPSHpBFIYzJP0CLB/tm5WM+X98L/0JTeimg2lYk0hG1+wJtvnY5Ie\nA14DBERE/ONgJ46ISuMb9hthWc1GpLwf/oIFDgSzwQx2+WhrSt1Rzcat/v3wfYc1s8oGC4U/RcSf\n61YSsxpwP3yz4RksFLaUdEqlJyPi3BqUx2xUuR++2fAMFgoTgI1JbQhm49L8+Wtv8+Ujs8oGC4Xl\n1QxQMzOz5jFYl1TXEMzM1jODhYK7jpqZrWcqhkJEPF/PgpiZWeNVM6LZzMzWEw4FMzMrOBTMzKzg\nUDAzs4JDwczMCg4FMzMrOBTMzKzgUDAzs4JDwczMCg4FMzMrOBTMzKzgUDAzs4JDwczMCg4FMzMr\nOBTMzKzgUDAzs4JDwczMCg4FMzMrOBTMzKwwsREvKmkZsBJYDayKiPZGlMPMzPpqSChkOiLi2Qa+\nvpmZ9ePLR2ZmVmhUKATwS0m9kjobVAYzM+unUZeP3hMRT0raErhJ0h8i4tbyHbKw6ASYOXNmI8po\nZrbeaUhNISKezB5XAFcBewywz/kR0R4R7dOmTat3Ec3M1kt1DwVJbZI2yf8G3gfcX+9ymJnZ2hpR\nU5gO3CbpHuAO4PqI+K8GlMPMamzhQuju7rutuzttt7Gp7qEQEY9FxC7Z8vaI+Fq9y2Bm9TF3Lhx1\nVCkYurvT+ty5jS2XVdbIcQpm1uQ6OmDJkhQEXV2weHFa7+hodMmsEo9TMLOa6uhIgXDmmenRgTC2\nORTMrKa6u1MNYcGC9Ni/jcHGFoeCmdVM3oawZAmccUbpUpKDYexyKJhZzfT09G1DyNsYenoaWy6r\nTBHR6DIMqb29PZYuXdroYpiZjSuSeoc7C7VrCmZmVnAomJlZwaFgZmYFh4KZmRUcCmZmVnAomJlZ\nwaFgZmYFh4KZmRUcCmZmVnAomJlZwaFgZmYFh4KZmRUcCmZmVnAomJlZwaFgZmYFh4KZmRUcCmZm\nVnAomJlZwaFgZmYFh4KZmRUcCmZmVnAomJlZoSGhIOkASQ9LelTS5xtRBjMzW1vdQ0HSBOA/gAOB\nnYEPS9q53uUwM7O1NaKmsAfwaEQ8FhGvA5cBhzWgHGZm1s/EBrzm1sDjZetPAHv230lSJ9CZrb4m\n6f46lK1RtgCebXQhaqiZ318zvzfw+xvv3jbcAxoRClWJiPOB8wEkLY2I9gYXqWb8/savZn5v4Pc3\n3klaOtxjGnH56Elg27L1bbJtZmbWYI0IhR5gB0nbSZoEHA1c24BymJlZP3W/fBQRqyT9H+BGYAJw\nQUQ8MMRh59e+ZA3l9zd+NfN7A7+/8W7Y708RUYuCmJnZOOQRzWZmVnAomJlZYUyHQrNPhyFpmaT7\nJN09kq5jY42kCyStKB9TImmqpJskPZI9bt7IMq6LCu/vdElPZp/h3ZIOamQZ14WkbSV1S3pQ0gOS\nTs62j/vPcJD31hSfn6TJku6QdE/2/r6SbR/2Zzdm2xSy6TD+B5hHGuDWA3w4Ih5saMFGkaRlQHtE\nNMXgGUnvBV4CfhwR78i2LQSej4izsmDfPCI+18hyjlSF93c68FJEnNPIso0GSTOAGRFxp6RNgF7g\ncOCjjPPPcJD3dhRN8PlJEtAWES9JagVuA04GjmCYn91Yril4OoxxJiJuBZ7vt/kw4KLs74tI/xHH\npQrvr2lExPKIuDP7eyXwEGkGgnH/GQ7y3ppCJC9lq63ZEozgsxvLoTDQdBhN8yFmAvilpN5sWo9m\nND0ilmd/PwVMb2RhauRESfdml5fG3aWVgUiaBewG3E6TfYb93hs0yecnaYKku4EVwE0RMaLPbiyH\nwvrgPRGxK2nG2H/JLk80rUjXKsfm9cqRWwy8FdgVWA4samxx1p2kjYErgE9HxIvlz433z3CA99Y0\nn19ErM6+T7YB9pD0jn7PV/XZjeVQaPrpMCLiyexxBXAV6ZJZs3k6u56bX9dd0eDyjKqIeDr7z7gG\n+D7j/DPMrkdfAVwSEVdmm5viMxzovTXb5wcQES8A3cABjOCzG8uh0NTTYUhqyxq8kNQGvA9oxplg\nrwWOz/4+HrimgWUZdfl/uMwHGMefYdZY+UPgoYg4t+ypcf8ZVnpvzfL5SZomabPs7w1JHXT+wAg+\nuzHb+wgg6x72b5Smw/hag4s0aiS9lVQ7gDTdyKXj/f1J+imwD2k64qeBLwNXA0uAmcCfgaMiYlw2\n1lZ4f/uQLj0EsAz4RNk13HFF0nuA3wD3AWuyzV8kXXsf15/hIO/twzTB5yfpH0kNyRNIP/aXRMQZ\nkt7EMD+7MR0KZmZWX2P58pGZmdWZQ8HMzAoOBTMzKzgUzMys4FAwM7OCQ8HGPUkvDb1Xse8+kt5V\ny/IM8fqflnTcKJznMkk7jEaZzMo5FGx9sw/QkFCQNBH4OHDpKJxuMTB/FM5j1odDwZqSpEMk3S7p\nLkm/lDQ9mwjtk8Bnsrnz985Ggl4hqSdb3p0df3o2Qdotkh6TdFLZuY/LJlC7R9LFkjaR9KdsGgUk\nbVq+XmZf4M6IWJXtd4uk8yQtlfSQpLmSrszmvv9qtk+bpOuz17pf0v/OzvUbYP8saMxGjf9BWbO6\nDdgrIkLSPwHzI+JUSd+lbP58SZcC50XEbZJmAjcCs7Nz7AR0AJsAD0taDOwI/Cvwroh4VtLUiFgp\n6Rbg/aQR3EcDV0bEG/3K9G7SPP7lXo+IdqWbvlwDzCFNz/1HSeeRajZ/jYj3Z+WdAhARayQ9Cuwy\nwDnNRsyhYM1qG+D/ZnPbTAL+VGG//YGd09Q4AGyazaQJcH1EvAa8JmkFadrhfYHL8xsjlU0Z8APS\n5ZyrgY8B/zzAa80gzeNfLp/P6z7ggXyKBUmPkSaEvA9YJOls4LqI+E3ZsSuArXAo2Cjy5SNrVt8G\nvhMR7wQ+AUyusF8LqUaxa7ZsXXazktfK9lvNID+iIuK3wCxJ+wATImKgidVeGaAc+Wus6fd6a4CJ\nEfE/wO6kcPiqpC+V7TM5O6fZqHEoWLOaQmmq9ePLtq8kXQ7K/QI4MV+RtOsQ570ZODKbaAxJU8ue\n+zGpEfnBT5jDAAAA30lEQVTCCsc+BGw/ZMnLSNoKeDkifgJ8kxQQuR0Zp7N62tjlULBmsJGkJ8qW\nU4DTgcsl9QLl98D+GfCBvKEZOAlozxqOHyQ1RFcUEQ8AXwN+LekeoHyK6UuAzYGfVjj858Bwb6T0\nTuCO7I5aXwbyBujpwCsR8dQwz2c2KM+SajZKJH0IOCwijh1kn6tIjd6PrONrfQZ4MSJ+uC7nMevP\nDc1mo0DSt0m3VT1oiF0/T2pwXqdQAF4ALl7Hc5itxTUFMzMruE3BzMwKDgUzMys4FMzMrOBQMDOz\ngkPBzMwK/x9ZT+sMtO/Z9gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x9136438>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#本文是高斯分布用于异常检测\n",
    "#load data set \n",
    "import scipy.io as sio\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "data = sio.loadmat('ex8data1.mat')\n",
    "X = data['X']\n",
    "Xval = data['Xval']\n",
    "yval = data['yval']\n",
    "print X.shape,Xval.shape,yval.shape\n",
    "print X[:5]\n",
    "print yval[:5]\n",
    "\n",
    "plt.plot(X[:,0],X[:,1],'bx')\n",
    "plt.xlim(0, 30)\n",
    "plt.ylim(0, 30)\n",
    "plt.xlabel('Latency (ms)')\n",
    "plt.ylabel('Throughput (mb/s)')\n",
    "plt.show()\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1L, 2L) (1L, 2L)\n[[ 14.11222578  14.99771051]] [[ 1.83263141  1.70974533]]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.06470829]\n [ 0.05030417]\n [ 0.07245035]\n [ 0.05031575]\n [ 0.06368497]]\n[[ 0.04163207]\n [ 0.08190901]\n [ 0.04071578]\n [ 0.06190003]\n [ 0.07118676]]\n"
     ]
    }
   ],
   "source": [
    "#高斯分布\n",
    "from scipy import stats\n",
    "#求得均值和方差\n",
    "def estimateGaussian(x):\n",
    "    #train set x fit the mu,sigma2\n",
    "    m,n = x.shape\n",
    "    mu = np.mean(x,axis=0).reshape(1,-1)\n",
    "    #sigma2 = np.sum(np.square(x - mu),axis=0).reshape(1,-1) / m\n",
    "    sigma2 = np.var(x,axis=0).reshape(1,-1)\n",
    "    \n",
    "    return mu,sigma2\n",
    "#计算高斯概率  \n",
    "def p(x,mu,sigma2):\n",
    "    #x is a new example:[m*n]\n",
    "    m,n = x.shape\n",
    "    p_list=[]\n",
    "    \n",
    "    for j in xrange(m):\n",
    "        p = 1\n",
    "        for i in xrange(n):\n",
    "            #p *= np.exp(-np.power((x[j, i] - mu[0, i]), 2) \\\n",
    "             #       / (2 * sigma2[0, i]**2)) \\\n",
    "             #/(np.sqrt(2 * np.pi * sigma2[0, i]))\n",
    "            p *= stats.norm.pdf(x[j,i], mu[0,i], np.sqrt(sigma2[0,i]))\n",
    "            # stats.norm.pdf(x,mean,sigma)\n",
    "        p_list.append(p)\n",
    "    p_array = np.array(p_list).reshape(-1,1)\n",
    "    return p_array\n",
    "\n",
    "mu,sigma2 = estimateGaussian(X)\n",
    "print mu.shape,sigma2.shape\n",
    "print mu,sigma2\n",
    "\n",
    "p_train = p(X,mu,sigma2)\n",
    "print p_train[:5]\n",
    "\n",
    "p_val = p(Xval,mu,sigma2)\n",
    "print p_val[:5]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "#cross validation for select threshold\n",
    "#交叉验证用于阈值选定,这里用的是F1 score这个评估指标\n",
    "def selectThreshold(y,pval):\n",
    "    bestEpsilon = 0 \n",
    "    bestF1 = 0\n",
    "    \n",
    "    stepSize = (np.max(pval) - np.min(pval)) / 1000\n",
    "    \n",
    "    for epsilon in np.arange(np.min(pval), np.max(pval), stepSize):\n",
    "        predictions = (pval < epsilon)\n",
    "        \n",
    "        fp = np.sum((predictions==1) & (y==0))\n",
    "        fn = np.sum((predictions==0) & (y==1))\n",
    "        tp = np.sum((predictions==1) & (y==1))\n",
    "\n",
    "        if tp + fp == 0:\n",
    "            precision = 0\n",
    "        else:\n",
    "            precision = float(tp) / (tp + fp)  #note!!!!float!!! \n",
    "        if tp + fn == 0:\n",
    "            recall = 0\n",
    "        else:\n",
    "            recall = float(tp) / (tp + fn)\n",
    "        \n",
    "        if precision + recall == 0:\n",
    "            F1 = 0\n",
    "        else:\n",
    "            F1 = 2.0 * precision * recall / (precision + recall) \n",
    "        if F1 > bestF1:\n",
    "            bestF1 = F1\n",
    "            bestEpsilon = epsilon\n",
    "    \n",
    "    return bestEpsilon,bestF1\n",
    "        \n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Best epsilon found using cross-validation: 8.990853e-05\nBest F1 on Cross Validation Set:  0.875000\n(you should see a value epsilon of about 8.99e-05)\nyou should see a Best F1 value of  0.875000)\n"
     ]
    }
   ],
   "source": [
    "#test the cs\n",
    "epsilon,F1 = selectThreshold(yval,p_val)\n",
    "print \"Best epsilon found using cross-validation: %e\"%(epsilon)\n",
    "print \"Best F1 on Cross Validation Set:  %f\"%(F1)\n",
    "print \"(you should see a value epsilon of about 8.99e-05)\"\n",
    "print \"you should see a Best F1 value of  0.875000)\"\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Outliers found: 6\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XucVXW9//HXZy7IQILSwIToiE2aWB0vjDZk9nBUekgl\nXg7NzwywEqkpUrEy6XEw0zL1hJ5T+KAMLNHIBk0jlcJirAwvMN4vdWQMUEIRLS+pE8jn98d3rTV7\nhtkze4bZe8/e834+Huux9lr7u9b6bjazPvt7XebuiIiIAJTkOwMiIjJwKCiIiEhCQUFERBIKCiIi\nklBQEBGRhIKCiIgkshYUzGyomT1gZo+Y2RNm9q1o/ygzu8vMno7We2crDyIi0juWrXEKZmbAcHd/\n3czKgXuAc4HTgJfd/XIzuxDY292/npVMiIhIr2StpODB69FmebQ4cDJwfbT/euCUbOVBRER6pyyb\nJzezUqAFeA9wjbvfb2ZV7r4lSvI8UJXm2NnAbIDhw4dPPPjgg7OZVRGRotPS0rLN3Uf35pisVR91\nuIjZXsCtwJeBe9x9r5T3/uHu3bYr1NbW+rp167KcSxGR4mJmLe5e25tjctL7yN3/CTQDJwIvmNlY\ngGi9NRd5EBGRnmWz99HoqISAmVUAk4G/ACuAM6NkZwK/ylYeRESkd7LZpjAWuD5qVygBmtz9djO7\nF2gys7OAjUBDFvMgIiK9kLWg4O6PAod3sf8l4PhsXVdERPpOI5pFRCShoCAiIgkFBRERSSgoiIhI\nQkFBREQSCgoiIpJQUBARkYSCgoiIJBQUREQkoaAgIiIJBQUREUkoKIiISEJBQUREEgoKIiKSUFAQ\nEZGEgoKIiCQUFEREJKGgICIiCQUFERFJKCiIiEhCQUFERBIKCiIiklBQEBGRhIKCiIgkFBRERCSh\noCAiIgkFBRERSWQtKJjZfmbWbGZPmtkTZnZutP9iM9tsZg9Hy8eylQcREemdsiyeewfwFXd/0Mz2\nBFrM7K7ovavd/XtZvLaIiPRB1oKCu28BtkSvXzOzp4Bx2bqeiIjsvpy0KZjZeOBw4P5o15fN7FEz\nu87M9s5FHgaE1laYOxeqqqC0NKznzg37RUQGgKwHBTN7B3ALcJ67vwosAt4NHEYoSSxIc9xsM1tn\nZutefPHFbGcz+1auhLo6qKiANWugrS2sKyrC/pUr851DERHM3bN3crNy4Hbgt+5+VRfvjwdud/f3\nd3ee2tpaX7duXVbymBOtrVBXx40NKxg3bRL19e1vNTfD5pvvZXrTVLjvPqipyV8+RaSomFmLu9f2\n5phs9j4yYAnwVGpAMLOxKclOBR7PVh4GjIUL4eyzGTdtEg0NIRBAWDc0wLhpk2DWLLjmmvzmU0QG\nvayVFMzsw8CfgMeAndHubwCfIlQdObAB+HzUKJ1WwZcUqqpCVVFNTRIIGhth0SJoaiKUHFpb4eij\n4fnn851bESkSfSkpZLP30T2AdfHWndm65oC1bRvsvz8QAkBjI1x6KcyfT3tVUnV1SCcikkca0ZwL\nlZWwcSMQqowWLQoBYdGi9qokNm0K6URE8khBIRfOOAOWLEmqjpqa4JJLwjppY1i8OKQTEcmjbI5o\nlticOVBXx+ZXTqKpqb33UX19CAybb74XmhaH3kciInmkoJALNTWwdCnTZ06FEbOgelZoQ9i0ifpV\ni0NAWLpU3VFFJO9UfZQrU6aEkkBbW+hlVFER1m1tYf+UKfnOoYhIdgev9ZeC75IqIpIHA2rwmoiI\nFB4FBRERSSgoiIhIQkFBREQSCgoiIpJQUBARkYSCgoiIJBQUREQkoaAgIiIJBQUREUkoKIiISEJB\nQUREEgoKIiKSUFAQEZGEgoKIiCQUFEREJKGgICIiiW6f0WxmQ4FPAMcA+wBvAo8Dd7j7E9nPnoiI\n5FLaoGBm3yIEhLuB+4GtwFDgIODyKGB8xd0fzUE+RUQkB7orKTzg7t9M895VZjYGqM5CnkREJE/S\nBgV3v6PzPjMrAd7h7q+6+1ZC6UFERIpEjw3NZrbMzEaY2XBCe8KTZva17GdNRERyLZPeR4e4+6vA\nKcBK4ABgRk8Hmdl+ZtZsZk+a2RNmdm60f5SZ3WVmT0frvXfrE4iISL/JJCiUm1k5ISiscPftgGdw\n3A5CQ/QhQB3wJTM7BLgQ+L27Hwj8PtoWEZEBIJOg8CNgAzAc+KOZ7Q+82tNB7r7F3R+MXr8GPAWM\nA04Gro+SXU8INiIiMgCkDQpmNsnMzN2/7+7j3P1j7u7AJqC+Nxcxs/HA4YSurVXuviV663mgKs0x\ns81snZmte/HFF3tzORER6aPuSgozgRYzu8nMPmNm7wLwYEemFzCzdwC3AOdFbROJKMh0WRXl7te6\ne627144ePTrTy4mIyG7orktqI4CZHQxMAX5qZiOBZuA3wJ/d/e3uTh61RdwC/MzdfxntfsHMxrr7\nFjMbi7q1iogMGD22Kbj7X9z9anc/ETgOuAf4JKEqKC0zM2AJ8JS7X5Xy1grgzOj1mcCv+pJxERHp\nf93OfRQzsyOADxOqev7s7l/O4LCjCV1XHzOzh6N93wAuB5rM7CxgI9DQ61yLiEhW9BgUzOwiQskg\nrv75iZktd/dvd3ecu98DWJq3j+9VLkVEJCcyKSl8GjjU3d8CMLPLgYeBboOCiIgUnkzGKfydMDtq\nbA9gc3ayIyIi+dTd1Nk/ILQhvAI8YWZ3RduTgQdykz0REcml7qqP1kXrFuDWlP13Zy03ItI3ra2w\ncCEsWwbbtkFlJZxxBsyZAzU1+c6dFJDuxilcn+49ERlAVq6EmTPh7LNhzRrYf3/YuBGWLIG6Oli6\nFKZMyXcupUBkMnX2J8zsITN72cxeNbPXzKzHuY9EJAdaW2HmTG5sWEHz5MtCqaCsDGpqaJ58GTc2\nrAgBo7U13zmVApFJQ/P/EAaZvdPdR7j7nu4+Isv5EpFMLFwIZ5/NuGmTaGiA5uawu7kZGhpg3LRJ\nMGsWXHNNfvMpBSOToPAs8Hg0T5FIUbvyyvYba6y5OewfkJYtg7POor4emppCILjoorBuaoL6ekJQ\nWLYs3zmVApHJOIULgDvN7A9AW7yz09QVIkXhyCM73lDjX9xNTfnOWRrbtoU2BEJ+Gxvh0kth/vwo\nIABUV4d0IhnIpKTwHeANwliFPVMWkaLT7S/ugaiyMjQqEwLYokUhICxalFLi2bQppBPJQCYlhX3c\n/f1Zz4nIAJH2F/dAdMYZsGQJzZMv6xDA6utTAtqqxSGdSAYyCQp3mtlH3X1V1nMjkkdXXhmqj6D9\nF/f3vw8vvAA/+lH69KlBo7kZ1q6FCy7ITZ6ZMwfq6tj8ykk0NU1K8hKXeDbffC80LYb77stRhqTQ\nZVJ91Aj8xszeVJdUKWZHHgmnnAKnntr+i9sdfvGLXRuf4/Rd9fiJA0tO1NTA0qVMb5pK/ap5oevp\n9u3Q2kr9qnlMb5oaxiloAJtkKJPnKezp7iXuXqEuqVKIMu1RVF8Pp58eAkF8g7/tNrj11vDrv7MB\n0/4wZUooCbS1wdFHQ0VFWLe1hf0auCa94e5dLsD4dO9F7xuwb3dp+muZOHGii/TkiivcV6/uuG/1\navfZs90rK9vfW72647a7+5Qp7gsWhNfz57uD+8SJYX/qua64YtfrxumPP77r63d1jEguAOu8l/fb\n7m76ywmP0pwJvA8YA1QTnr52KbAGmNzbC/ZlUVCQTHS+2adux6/nz981ILiHgGDm3tgY3p84Mfx1\nTJ3a9bndw81+wYL2844Y4T58eAhC6Y4RyaV+DQrhfBxC6JJ6N/BX4CFgGTAdGNrbi/V1UVCQTHV3\n849/0c+f374vtXTR2BjeHzu2PSCYuU+e7F5e3l6SiI+J03/84+G9T3wibJeVhWOGD28PSLNnq8Sw\nW9avdz/vPPcxY9xLSsL6vPPCfkmr34PCQFkUFKQ3urr5z57tPnJke7A46qhwU0+tWvrEJ9yHDAnH\nHnBAOO7AAz2pSqqsDOkbG0OQKCsLr4cODWnKy0MgKSkJ2wceGN7fY49w7dmzQ/BIDQ6qXsrAnXeG\nf/x580IQ2L49rOfNC/vvvDPfORywFBRk0OtcUohvxCNHhuqd1avDdllZ+N/f2Bj2DR8etuOAAO4H\nHRRKCvFNfurUcOOP05WXu8+Y0b5dUhICQPw63l9a6l5X5/7BD4bzLVjQHoTi0oSksX69e2Wl3/DF\nNV2219zwxTXhi1aJoUsKCjKoddWmMHJkuFEvWBACxPDh7hUVYV1XF/4C3vWu9hv45Mkh3YgRYbuy\nMhwbv5+6xMGjpKRjECgpCTf/eDsOQHEQivMUBwjpxnnnuc+b1217kV94ofvcuXnN5kCloCCDWk+9\nj1J/1c+fH/Z3vpl3vsGD+7hxuwaEzmlKSzve/LtKZxbSxOlmzMjPv1NBGTMmKQWkbS9av969qip/\neRzA+hIUMnmewu8z2SeSbxdcEMYTpI5JWLsW3vteOPRQuOEG2GMPKC2Fyy6Dn/wEdu4Es5C2vDxs\np+4D2NzpieRmIU1s2DB4++3wGIO99+6YNjWdO+zYEZZjjgnPxulqUJykSDPhX2OjJvzLlrRBwcyG\nmtkooNLM9jazUdEyHhiXqwyKZOrKK8ON+aST4KpoDt8NG+ArXwk339LSsLiHm/j27TB6NAwZEtK2\ntbUHA08zUfzEibu+98YbYf3WW/CPf6TPX3zukhJoaQkPRUsdES1d0IR/OWee5n+/mZ0LnAfsA/w9\n5a1XgR+7+8LsZy+ora31devW9ZxQBpV47qG1a9unljj1VJgwIQzkPfBAWL8+BIrt28NNufOv/PLy\n8F6mSkp2PX7nzhBkMjl2yBAYMybcx4YNC796d+zI4VxJhWbuXKio2GXCv9QpzetXzQsR/SrN5t+Z\nmbW4e21vjklbUnD3/3X3A4CvuvsBKcuhuQwIIunEcw+VlYX1Qw+FX/GPPRZKBE8/DcOHhxv3PvuE\n91Jv6NC7gAAdjzcLx2cSEOJj3/e+EBAmTIDp0xUQejRnDvz4x2y++d4OU4h0mPBv8WL40pfym88i\nksksqa+Y2czOO919aRbyI5Kx1LmHpkyBr3413GibmsKNer/94NlnQ1DYurX/r5+uiilVZWXH6u6W\nlpCfY44J8yzldYbVQhBP+DdzKoyYBdWzQhvCpk1hSvCmxZrwr59lMkvqkSnLMcDFwNQs5kkGqb48\nCnPt2hAQbrgBqqrCuq0tVB09+yyMGNGx6ijX6upg6NCO+8rL4e9/HwAzrBYKTfiXW73trgTsBfwm\ng3TXAVsJz3eO910MbAYejpaPZXJNdUkdHLrti55GPGdRPFdRahfR0tLw+n3v6zhuIJdLeXnHAW/x\nvhEjOs6bpDmSJBvIRpfULvwLOCCDdD8FTuxi/9Xufli03NmH60uR6u1U1M3N8N3vwhe+AA8+GEoF\nEEoEb78dlooKePLJ3JcS4uvt2NHebjF+fFhv3w6f/nR4r8suliJ5lMk4hV+b2YpouYMwMd6tPR3n\n7n8EXu6HPMogkrYvehcaG+GTnww32+nT4dXo0U9xXX9lJbz5ZtcNzLsrXZBJ7dI6dGgIVCUlocpo\nwwaYMSM0fjc3h+fhdNnFUiSPMmlo/l7K6x3ARnd/bjeu+eWo4Xod8BV377Jnt5nNBmYDVFdX78bl\npJB07oteX9/e5TRex10S3/OekGbCBHjxxTCGoKWl/VzbtoVeSJn2DuoNd9hrL/jnP3fd39gYXi9a\nFAKAWSgVzJgRBqzNnAnXXQdbtoQH+OzyTGWVGCSPMnny2h8IpYORwChCYOirRcC7gcOALcCCbq57\nrbvXunvt6NGjd+OSMlB1blhubg6PwzztNLjkkvaqpLjLaby+6qqwPu640O//qafCjb+lJaRJlY2A\nEOscEGIPPgjLl4cg9dJLIU/f+14ICFOmwA9/CCec0B4QoL3qrKsnvInkUibVR7OAB4DTgGnAfWb2\nub5czN1fcPe33X0n8GPgqL6cR4pD52cc33RT+FV9+ulhO75R7tgR1t/9bnvX0+rqEDg+8xl417vC\nSOJ4IFl5+a7BIVdKSuD+++Hww0NwOOKIMLXG4YeHEsQNN4Sqro98ZNcSQX29uqPKANBTSzShlPDO\nlO13An/NpBUbGE/H3kdjU17PBW7K5DzqfVS8enoiWqr4OQnHHBPWFRXtD7pJnfoa3CdMcN9rr9z1\nMkrt3RS/PuKI9s8YT92tnkaSS2Sp99FLwGsp269F+7plZj8H7gXea2bPmdlZwJVm9piZPQrUR4FB\nBrFMG5bjtoYZM+Cee8LaPew76CD4179CTyOAkSNDldKbb+buc7iH9ov4dUlJmGIjLgW5hxJQarWY\nGpZlQOopagBLCY/hvBj4JvAgobvp+cD5vY1CfVlUUihemZQU4jRxv/54ffDB7dNQT5wYfqGPGhW2\n4+ch5GN8QllZx3EIs2d3PaW3nrgm2UaWSgqtwG1A1NGPXwF/A/aMFpE+SZ3UrLtf0GvXdmxbOP/8\nsP7IR0LbwYEHhvr7E06Al1+GffcN3VOHD89sKor+Nm8e3HZbe35ratR+IIUj7SypA4lmSS1O8Syn\nfZn/JzWgrF0bxgD88IehYbelJQSG53an43Qfxd1O1bVUBoJ+nSU15aQHmdm1ZrbKzFbHS9+zKbL7\n4tJDfX0ILMuXh26fY8bAhz4UAsKencqxJX0Zv5+ic4+m1OcjlJaGXk8rVoSSgtoMpFBl8meynNCm\n8F/A11IWkd3SuUtqJpPCxWMbLrig/Zf4TTeFsQ3nnw9f+xo88kgYI/DaazBuXLh5f+hDobtqd4Fh\n6NDQfTSdIUPCeTofM2RIyNdvfxuqq/76V405kMKVSVDY4e6L3P0Bd2+Jl6znTIpeb+c6gq4DyS9/\n2T62Ye3a0D7x6KMweXIYYPaFL4T6/fgBOXGg6Oytt9K3QUyYEOYseuCBsL3PPiGAnHAC3HlnOH99\nfWhLiNsQ1GYgBSldCzRh9PIoQq+jLwJjU/aN6m2L9u4s6n1U3OLxB/PnZ5Z+9uzQ7z+1x1Lcmyfu\nqdTYGHoeNTaG7QkTwjUOOqi9l1BJSfoeRNXVYR3PtFpS0j7b6YIFIR+ZzOIqkk/0ofdRd+M+Wwg9\njuLfVKlVRk6YrkJkt3Q111FPDbSnnx5GBl96aTgOOjY6x6WN97wnjH4ePz6MW2hsDNNerF/fXpVU\nVgb//neoAvr3v8O59tgjPIshTr92bXiq286doSH5/PNDutSpKdSoLEWjt1EkH4tKCsWpL89PiNON\nHBlGNA8bFsYEpDsmHv38gQ+E7biUcfjhoRQwbJj75Mm7lhTGjm0fE5HptUQGGrIxTsHMTutiOd7M\nxuQgZkkRS/1VD5lNCtfcDCedFNogvvpVeOONUJ//0EMdn9B25ZXwxS+G0c/HHAOPPw4nnww/+1mY\niO7BB+Gss0Lbwp/+FI6JexcNHRpmMK2qCu0T7nDHHXD77SH9KaeoZ5EUsZ6iBnAH4bkIt0TLS8Aq\n4GlgRm+jUF8WlRQkdsUV7b/e47mEhg1zHzq04y/4eE6kxsaO21OntqdZvTqcY7/9wv74qWipbQ8f\n/GDH865eHUobGo0shYAsjWguAya4+3+6+38ChxDaFD4IfL3/w5RIehdcEGYcdW/vQVReHtoEUm3Y\nEJ6PvGxZKFUsXx7aCJ5/vr1EsXZt6C20aVPoWdTYGOZPOv74MDK6sRFGjerYXlBfDz/6kXoWSfHK\nJCjs5+4vpGxvjfa9DGzPTrZE0otv5uecExqbzzknbKdWO915J1x2WWg8jifb++Qn4ZlnQrfWzqOp\nv/a1ULV0+unwu9+Faqzly8N+kcEkk6Bwt5ndbmZnmtmZhLmP7jaz4UCax4yIZE/8Kz2111Lq/lRD\nhoRf/wsWhLaA1FHQmTzLQQPQZNDpqX6J0CV1GnB1tEwjmjMpV4vaFCRVJr2WUvfF4yCGDes6jZ5x\nIMWKbLQpROe+2d3nRsvN0cVE8iKTXktxGmgvUZSVhRJBLNNnOYgMJj3Okmpmr9E+bfYQoBz4l7uP\nyHLeEpolVfoidSbV+vr0242NIXBoZlMpNn2ZJbXHJ9m6ezLXpJkZcDJQ1/vsieRWTyWK1ABRX5/Z\n3Esixa5Pz1Mws4fc/fAs5KdLKilIf9udZzmIFIqslBTM7LSUzRKgFnirl3kTGVC6uvFnMu+SSLHr\nMSgAJ6W83gFsIFQhiYhIkcmkTeGzuciIiIjkXyYT4u1rZrea2dZoucXM9s1F5kREJLcyGdH8E2AF\nsE+0/DraJyKFrrUV5s4NU8KWlob13LlhvwxKmQSF0e7+E3ffES0/BUZnOV8ikm0rV4ZZAysqYM0a\naGsL64qKsH/lynznUPIgk6DwkplNN7PSaJlOmD5bRApVayvMnMmNDStonnxZeLB0WRnU1NA8+TJu\nbFgBM2eqxDAIZRIUPgc0AM8DWwhzH6nxWaSQLVwIZ5/NuGmTOkwMGI/yHjdtEsyaBddck998Ss51\n2/vIzEqB09x9ao7yIyK5sGwZrFlDfU0Yxd3ldB/Vs+Doo+Gqq/KdW8mhbksK7v428Km+nNjMrot6\nKz2esm+Umd1lZk9H6737cm4R2U3btsH++wPdTAxYXR3SyaCSSfXRn81soZkdY2ZHxEsGx/0UOLHT\nvguB37v7gcDvo20RybXKSti4EQhVRqnPpkieP71pU0gng0omI5oPi9aXpOxz4LjuDnL3P5rZ+E67\nTwaOjV5fD9yNHukpkntnnAFLltA8+bL0EwOuWhzSyaCSyYjm/pwNpsrdt0Svnweq0iU0s9nAbIDq\n6up+zIIMKq2toVF12bJQFVJZGW50c+aEHjeD1Zw5UFfH5ldOoqlp0i4zyW6++V5oWgz33ZfffErO\nZTKieQ8zO8PMvmFmF8XL7l44elBP2ila3f1ad69199rRozUsQvpA/fDTq6mBpUuZ3jSV+lXzQvDc\nvh1aW6lfNY/pTVNh6dLBHTgHqUzaFH5FqPbZAfwrZemLF8xsLEC03trH84h0T/3wezZlSigJtLWF\nXkYVFWHd1hb2T5mS7xxKHmTy5LXH3f39fTp5aFO4PT7ezP4beMndLzezC4FR7t7j7PV6noL02ty5\nUFGxS515h6evrZoXboDqcilFqi/PU8ikpLDGzD7Qh8z8HLgXeK+ZPWdmZwGXA5PN7GnghGhbpP8t\nWwZnnZXUkTc0wEUXdXq62qxZIZ2IJNI2NEfjC3ZGaT5rZs8AbYARmgT+o7sTu3u68Q3H9zGvIplL\n0w9//nz1wxfpTne9j8bR3h1VpLDE/fBranbph588YU398EV20V1Q+Ju7b8xZTkT6k/rhi/RJd0Fh\njJmdn+5Nd1frnAxc6ocv0ifdBYVS4B2ENgSRwhL3w585FUbMCpO7VVfDpk2hhNC0WP3wRbrQXVDY\n4u6XdPO+yMAW98O/5prQ/z51RPN99ykgiHShu6CgEoIUvpqaMA5BYxFEMtLdOAV1HRURGWTSBgV3\nfzmXGRERkfzLZESziIgMEgoKIiKSUFAQEZGEgoKIiCQUFEREJKGgICIiCQUFERFJKCiIiEhCQUFE\nRBIKCiIiklBQEBGRhIKCiIgkFBRERCShoCAiIgkFBRERSSgoiIhIQkFBREQSCgoiIpJQUBARkURZ\nPi5qZhuA14C3gR3uXpuPfIiISEd5CQqRenfflsfri4hIJ6o+EhGRRL6CggO/M7MWM5udpzyIiEgn\n+ao++rC7bzazMcBdZvYXd/9jaoIoWMwGqK6uzkceRUQGnbyUFNx9c7TeCtwKHNVFmmvdvdbda0eP\nHp3rLIpIf2hthblzoaoKSkvDeu7csF8GpJwHBTMbbmZ7xq+BjwKP5zofIpJlK1dCXR1UVMCaNdDW\nFtYVFWH/ypX5zqF0IR/VR1XArWYWX3+Zu/8mD/kQkWxpbYWZM7mxYQXjJk+ivibaX1ND8+TL2PzK\nSUyfORXuuw9qaro9leRWzksK7v6Mux8aLe9z9+/kOg8ikmULF8LZZzNu2iQaGqC5OexuboaGBhg3\nbRLMmgXXXJPffMouzN3znYce1dbW+rp16/KdDRHJVFVVqCqqqUkCQWMjLFoETU1QX08oTRx9NDz/\nfL5zW7TMrKW3g4M1TkFE+t+2bbD//kAIAI2NcOmlYV1fH6Wprg7pZEBRUBCR/ldZCRs3AqHKaNEi\nmD8/rOOqJDZtCulkQFFQEJH+d8YZsGRJUnXU1ASXXBLWSRvD4sUhnQwo+Zz7SESK1Zw5UFfH5ldO\noqlpUlJlVF8fAsPmm++FpsWh95EMKAoKItL/ampg6dLQ7XTELKieFdoQNm2iftXiEBCWLlV31AFI\n1Ucikh1TpoSSQFtb6GVUURHWbW1h/5Qp+c6hdEFdUkVEipS6pIqIyG5RUBARkYSCgoiIJBQUREQk\noaAgIiIJBQUREUkoKIiISEJBQUREEgoKIiKSUFAQEZGEgoKIiCQUFEREJKGgICIiCQUFERFJKCiI\niEhCQUFERBIKCiIiklBQEBGRhIKCiIgkFBRERCSRl6BgZiea2V/NbL2ZXZiPPIiIyK5yHhTMrBS4\nBpgCHAJ8yswOyXU+RERkV/koKRwFrHf3Z9z938BNwMl5yIeIiHRSlodrjgOeTdl+Dvhg50RmNhuY\nHW22mdnjOchbvlQC2/KdiSwq5s9XzJ8N9PkK3Xt7e0A+gkJG3P1a4FoAM1vn7rV5zlLW6PMVrmL+\nbKDPV+jMbF1vj8lH9dFmYL+U7X2jfSIikmf5CAprgQPN7AAzGwKcDqzIQz5ERKSTnFcfufsOM5sD\n/BYoBa5z9yd6OOza7Ocsr/T5ClcxfzbQ5yt0vf585u7ZyIiIiBQgjWgWEZGEgoKIiCQGdFAo9ukw\nzGyDmT1mZg/3pevYQGNm15nZ1tQxJWY2yszuMrOno/Xe+czj7kjz+S42s83Rd/iwmX0sn3ncHWa2\nn5k1m9mTZvaEmZ0b7S/477Cbz1YU35+ZDTWzB8zskejzfSva3+vvbsC2KUTTYfwfMJkwwG0t8Cl3\nfzKvGetHZrYBqHX3ohg8Y2YfAV4Hlrr7+6N9VwIvu/vlUWDf292/ns989lWaz3cx8Lq7fy+feesP\nZjYWGOvaSASEAAAFNUlEQVTuD5rZnkALcArwGQr8O+zmszVQBN+fmRkw3N1fN7Ny4B7gXOA0evnd\nDeSSgqbDKDDu/kfg5U67Twauj15fT/hDLEhpPl/RcPct7v5g9Po14CnCDAQF/x1289mKggevR5vl\n0eL04bsbyEGhq+kwiuZLjDjwOzNriab1KEZV7r4lev08UJXPzGTJl83s0ah6qeCqVrpiZuOBw4H7\nKbLvsNNngyL5/sys1MweBrYCd7l7n767gRwUBoMPu/thhBljvxRVTxQtD3WVA7O+su8WAe8GDgO2\nAAvym53dZ2bvAG4BznP3V1PfK/TvsIvPVjTfn7u/Hd1P9gWOMrP3d3o/o+9uIAeFop8Ow903R+ut\nwK2EKrNi80JUnxvX627Nc376lbu/EP0x7gR+TIF/h1F99C3Az9z9l9HuovgOu/psxfb9Abj7P4Fm\n4ET68N0N5KBQ1NNhmNnwqMELMxsOfBQoxplgVwBnRq/PBH6Vx7z0u/gPLnIqBfwdRo2VS4Cn3P2q\nlLcK/jtM99mK5fszs9Fmtlf0uoLQQecv9OG7G7C9jwCi7mH/Q/t0GN/Jc5b6jZm9m1A6gDDdyLJC\n/3xm9nPgWMJ0xC8A3wRuA5qAamAj0ODuBdlYm+bzHUuoenBgA/D5lDrcgmJmHwb+BDwG7Ix2f4NQ\n917Q32E3n+1TFMH3Z2b/QWhILiX82G9y90vM7J308rsb0EFBRERyayBXH4mISI4pKIiISEJBQURE\nEgoKIiKSUFAQEZGEgoIUPDN7vedUSdpjzexD2cxPD9c/z8xm9sN5bjKzA/sjTyKpFBRksDkWyEtQ\nMLMy4HPAsn443SLggn44j0gHCgpSlMzsJDO738weMrPfmVlVNBHaF4C50dz5x0QjQW8xs7XRcnR0\n/MXRBGl3m9kzZnZOyrlnRhOoPWJmN5jZnmb2t2gaBcxsROp2iuOAB919R5TubjO72szWmdlTZnak\nmf0ymvv+21Ga4WZ2R3Stx83s/0Xn+hNwQhRoRPqN/kNJsboHqHN3N7NZwAXu/hUz+yEp8+eb2TLg\nane/x8yqgd8CE6JzHAzUA3sCfzWzRcBBwH8BH3L3bWY2yt1fM7O7gY8TRnCfDvzS3bd3ytPRhHn8\nU/3b3WstPPTlV8BEwvTcrWZ2NaFk83d3/3iU35EA7r7TzNYDh3ZxTpE+U1CQYrUv8ItobpshwN/S\npDsBOCRMjQPAiGgmTYA73L0NaDOzrYRph48DlscPRkqZMmAxoTrnNuCzwNldXGssYR7/VPF8Xo8B\nT8RTLJjZM4QJIR8DFpjZFcDt7v6nlGO3AvugoCD9SNVHUqx+ACx09w8AnweGpklXQihRHBYt41Ie\nVtKWku5tuvkR5e5/Bsab2bFAqbt3NbHam13kI77Gzk7X2wmUufv/AUcQgsO3zeyilDRDo3OK9BsF\nBSlWI2mfav3MlP2vEaqDYquAL8cbZnZYD+ddDXwymmgMMxuV8t5SQiPyT9Ic+xTwnh5znsLM9gHe\ncPcbgf8mBIjYQRTorJ4ycCkoSDEYZmbPpSznAxcDy82sBUh9BvavgVPjhmbgHKA2ajh+ktAQnZa7\nPwF8B/iDmT0CpE4x/TNgb+DnaQ5fCfT2QUofAB6Inqj1TSBugK4C3nT353t5PpFuaZZUkX5iZtOA\nk919RjdpbiU0ej+9m9eaC7zq7kt25zwinamhWaQfmNkPCI9V/VgPSS8kNDjvVlAA/gncsJvnENmF\nSgoiIpJQm4KIiCQUFEREJKGgICIiCQUFERFJKCiIiEji/wNkhHEtZpld3AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0xa3746a0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "#可视化一下检测是异常值\n",
    "print \"Outliers found: %d\"%(np.sum(p_train < epsilon))\n",
    "\n",
    "#visualization:Draw a red circle around those outliers\n",
    "outliters = np.where(p_train.ravel() < epsilon)\n",
    "\n",
    "plt.plot(X[:,0],X[:,1],'bx')\n",
    "plt.plot(X[outliters,0],X[outliters,1],'ro',lw=2, markersize=10,\\\n",
    "         fillstyle='none', markeredgewidth=1)\n",
    "plt.xlim(0, 30)\n",
    "plt.ylim(0, 30)\n",
    "plt.xlabel('Latency (ms)')\n",
    "plt.ylabel('Throughput (mb/s)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1000L, 11L) (100L, 11L) (100L, 1L)\n[[  6.65311366  13.68970163  17.55284534 -11.93349739  -5.75640883\n   15.43592131 -12.13675197   8.22412257  -4.88402408   5.34933534\n   17.35137235]]\n"
     ]
    }
   ],
   "source": [
    "# load the data2\n",
    "#接下来用多元高斯分布模型\n",
    "import scipy.io as sio\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "data1 = sio.loadmat('ex8data2.mat')\n",
    "X1 = data1['X']\n",
    "Xval1 = data1['Xval']\n",
    "yval1 = data1['yval']\n",
    "print X1.shape,Xval1.shape,yval1.shape\n",
    "print X1[:1]\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(11L,) (11L, 11L) 60.9748937269 53.2057218642\n4.80014259158e+19\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[  7.43392946e-18]\n [  3.17719963e-18]]\n[[  5.14552429e-22]\n [  1.73600598e-17]]\n"
     ]
    }
   ],
   "source": [
    "# multivariate_normal 多元高斯分布\n",
    "from scipy import stats\n",
    "def get_mean_cov(x):\n",
    "    m,n = x.shape\n",
    "    mu = np.mean(x,axis=0)#平均值\n",
    "    Sigma = np.var(x, axis=0) #方差\n",
    "    #Sigma = np.diagflat(Sigma) #对角矩阵,[[1,0,0],[0,1,0],[0,0,1]]\n",
    "    Sigma = np.diag(Sigma)\n",
    "    detSigma = np.linalg.det(Sigma) #行列式\n",
    "    \n",
    "    return mu,Sigma,detSigma\n",
    "\n",
    "    \n",
    "def mul_P(x,mu,sigma,detSigma):\n",
    "    #x is vector =[n*1]\n",
    "    n = x.shape[0]\n",
    "    xx = x - mu #[n*1]\n",
    "    #dot点乘,note -n/2--->-n/2.0\n",
    "    p = np.exp(-0.5 * np.dot(xx.T.dot(np.linalg.pinv(sigma)), xx)) * \\\n",
    "            ((2*np.pi)**(-n/2.0) * (detSigma**(-0.5)))\n",
    "    return p\n",
    "        \n",
    "mu1,Sigma,detSigma=get_mean_cov(X1)\n",
    "\n",
    "print mu1.shape,Sigma.shape,Sigma[0,0],Sigma[1,1]  #compute correct\n",
    "print detSigma #correct \n",
    "p_train1 = np.array([mul_P(x.reshape(-1,1),mu1.reshape(-1,1),Sigma,detSigma) for x in X1]).reshape(-1,1)\n",
    "p_val1 = np.array([mul_P(x.reshape(-1,1),mu1.reshape(-1,1),Sigma,detSigma) for x in Xval1]).reshape(-1,1)\n",
    "\n",
    "\n",
    "print p_train1[:2]\n",
    "print p_val1[:2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Best epsilon found using cross-validation: 1.377229e-18\nBest F1 on Cross Validation Set:  0.615385\n(you should see a value epsilon of about 1.38e-18)\nyou should see a Best F1 value of  0.615385)\nOutliers found: 117\n"
     ]
    }
   ],
   "source": [
    "#test the cs\n",
    "epsilon1,F1_1 = selectThreshold(yval1,p_val1)\n",
    "print \"Best epsilon found using cross-validation: %e\"%(epsilon1)\n",
    "print \"Best F1 on Cross Validation Set:  %f\"%(F1_1)\n",
    "print \"(you should see a value epsilon of about 1.38e-18)\"\n",
    "print \"you should see a Best F1 value of  0.615385)\"\n",
    "\n",
    "print \"Outliers found: %d\"%(np.sum(p_train1 < epsilon1))\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}