8.anomaly_detection_and_recommendation/ex8.Anomaly_Detection_and_Recommendation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 异常检测和推荐系统"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "此Notebook是配合Andrew Ng \"Machine Leanring\"中[异常检测](https://github.com/loveunk/machine-learning-deep-learning-notes/blob/master/machine-learning/anomaly-detection.md)和[推荐系统](https://github.com/loveunk/machine-learning-deep-learning-notes/blob/master/machine-learning/recommender-system.md)部分学习使用。\n",
    "\n",
    "测试用python版本为3.6\n",
    "* 机器学习路径：https://github.com/loveunk/machine-learning-deep-learning-notes/\n",
    "* 内容正文综合参考网络资源，使用中如果有疑问请联络：https://www.kaikai.ai"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在本练习的第一部分，你将使用高斯模型实现异常检测算法，并将其应用于检测网络上的故障服务器。\n",
    "在第二部分，你将使用协作过滤构建一个推荐系统，并将其应用于电影推荐数据集。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. 异常检测（Anomaly detection）"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在这个练习中，你将实现一个异常检测算法，该算法用来检测服务器的异常状态。\n",
    "数据集的特征包括每个服务器的吞吐量（Throughput）、延迟（Latency）。\n",
    "当服务器在工作的时候，你收集了 $m=307$个样本，每个样本记录了服务器的状态，因此你有一个未标记的数据集 $\\{x^{(1)}, ..., x^{(m)}\\}$。\n",
    "你猜测数据集中大部分样本是“正常（normal）”的，但是有少量的样本是“异常（Anomal）”的。\n",
    "\n",
    "你将使用高斯模型来检测数据集。\n",
    "先从一个简单的二维数据集开始，可视化将帮助你感受该算法是如何工作的。\n",
    "在这个数据集上，你将拟合一个高斯分布，然后找到那些很低概率的值，这些值可以被认为是异常的。\n",
    "之后，你将把异常检测算法应用到一个有着更多维度的更大的数据集。\n",
    "\n",
    "下面先加载数据集，看一看。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sb\n",
    "from scipy.io import loadmat\n",
    "from scipy import stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(307, 2)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data = loadmat('data/ex8data1.mat')\n",
    "X = data['X']\n",
    "X.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAD8CAYAAABw1c+bAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAG7tJREFUeJzt3X+QHOWd3/H3R8tytWCKFacFo0Uy5IrSHViHRG3JXJRzgDtLQsEgq2IbzuUjxolMylRZKUdlYVKGQKrgopjLXewykYECJxzmYiSZnIUlFbiKmAocKyQhESELCDbaVZDuQAibTdCPb/6YXml21T07P3pnZqc/r6qpmXm6e+bZUeszPU8//TyKCMzMrDimtboCZmbWXA5+M7OCcfCbmRWMg9/MrGAc/GZmBePgNzMrGAe/mVnBOPjNzArGwW9mVjCnTbSCpFnAD4CPAseBtRHxF5LuBP4FcDBZ9ZsRsTFl+yXAXwBdwAMRce9E7zljxoy48MILq/0bzMwKb+vWrX8XEX3VrDth8ANHga9HxEuSzgK2StqSLPvziPgPWRtK6gK+C3wK2Ae8KOnJiPhfld7wwgsvZHBwsJr6m5kZIOmX1a47YVNPROyPiJeSx+8Du4H+Kl9/AfBaRLwRER8CPwSur7ZyZmaWv5ra+CVdCMwHXkiKbpX0sqSHJE1P2aQfeKvs+T6q/9IwM7NJUHXwS/oI8ASwMiIOA98DfgeYB+wHvp22WUpZ6nCgklZIGpQ0ePDgwbRVzMwsB1UFv6RuSqH/aESsA4iItyPiWEQcB75PqVlnvH3ArLLnFwDDae8REWsjYiAiBvr6qjo/YWZmdZgw+CUJeBDYHRH3lZWfX7baZ4BdKZu/CFws6SJJpwM3AE82VmUzM2tENb16FgJfBHZK2p6UfRO4UdI8Sk03bwJfAZA0k1K3zaURcVTSrcAmSt05H4qIV3L+GwDYsG2INZv2MHxohJm9PaxaPIdl8306wcxsPLXjDFwDAwNRS3fODduGuG3dTkaOHDtR1tPdxT3L5zr8zawQJG2NiIFq1u2IK3fXbNozJvQBRo4cY82mPS2qkZlZ++qI4B8+NFJTuZlZkXVE8M/s7amp3MysyDoi+FctnkNPd9eYsp7uLlYtntOiGpmZta9qevW0vdETuO7VY2Y2sY4IfiiFv4PezGxiHdHUY2Zm1XPwm5kVjIPfzKxgHPxmZgXj4DczKxgHv5lZwTj4zcwKxsFvZlYwDn4zs4Jx8JuZFUw1Uy/OkvQzSbslvSLpa0n5GkmvSnpZ0npJvRnbvylpp6TtkqqfXcXMzCZFNUf8R4GvR8TvAVcAX5V0CbAF+HhE/D7wC+C2Cq9xVUTMq3Z2GDMzmzwTBn9E7I+Il5LH7wO7gf6I2BwRR5PVngcumLxqmplZXmpq45d0ITAfeGHcopuBpzI2C2CzpK2SVlR47RWSBiUNHjx4sJZqmZlZDaoOfkkfAZ4AVkbE4bLy2yk1Bz2asenCiLgcuIZSM9En01aKiLURMRARA319fVX/AWZmVpuqgl9SN6XQfzQi1pWV3wRcC3whIiJt24gYTu4PAOuBBY1W2szM6ldNrx4BDwK7I+K+svIlwDeA6yLig4xtz5R01uhjYBGwK4+Km5lZfao54l8IfBG4OumSuV3SUuA7wFnAlqTsfgBJMyVtTLY9D/i5pB3A3wI/iYif5v9nmJlZtSacejEifg4oZdHGlLLRpp2lyeM3gMsaqaCZmeXLV+6amRWMg9/MrGAc/GZmBePgNzMrGAe/mVnBOPjNzArGwW9mVjAOfjOzgnHwm5kVjIPfzKxgHPxmZgXj4DczKxgHv5lZwTj4zcwKxsFvZlYw1czANUvSzyTtlvSKpK8l5edI2iJpb3I/PWP7m5J19iZTNZqZWQtVc8R/FPh6RPwecAWlCdMvAVYDT0fExcDTyfMxJJ0D3AF8gtJcu3dkfUGYmVlzTBj8EbE/Il5KHr8P7Ab6geuBR5LVHgGWpWy+GNgSEe9ExLvAFmBJHhU3M7P61NTGL+lCYD7wAnBeROyH0pcDcG7KJv3AW2XP9yVlZmbWIlUHv6SPAE8AKyPicLWbpZRFxuuvkDQoafDgwYPVVsvMzGpUVfBL6qYU+o9GxLqk+G1J5yfLzwcOpGy6D5hV9vwCYDjtPSJibUQMRMRAX19ftfU3M7MaVdOrR8CDwO6IuK9s0ZPAaC+dm4Afp2y+CVgkaXpyUndRUmZmZi1SzRH/QuCLwNWStie3pcC9wKck7QU+lTxH0oCkBwAi4h3gbuDF5HZXUmZmZi2iiNQm95YaGBiIwcHBVlfDzGzKkLQ1IgaqWddX7pqZFYyD38ysYBz8ZmYF4+A3MysYB7+ZWcE4+M3MCsbBb2ZWMA5+M7OCcfCbmRWMg9/MrGAc/GZmBePgNzMrmNNaXQEzm5o2bBtizaY9DB8aYWZvD6sWz2HZfE+wNxU4+M1q4LAr2bBtiNvW7WTkyDEAhg6NcNu6nQCF/DymGjf1mFVpNOyGDo0QnAy7DduGWl21pluzac+J0B81cuQYazbtaVGNrBYOfrMqOexOGj40UlO5tZcJm3okPQRcCxyIiI8nZY8Dc5JVeoFDETEvZds3gfeBY8DRaicJMGtHDruTZvb2MJTyd8/s7WlBbaxW1RzxPwwsKS+IiM9HxLwk7J8A1qVtmLgqWdehb1NaVqgVMexWLZ5DT3fXmLKe7i5WLZ6TsYW1kwmDPyKeBVLnyU0mYv8c8FjO9TJrOw67k5bN7+ee5XPp7+1BQH9vD/csn+sTu1NEo716/hB4OyL2ZiwPYLOkAP5zRKzNeiFJK4AVALNnz26wWmb5Gu3NM3LkGF0SxyLoL3CvHiiFf1H/9qmu0eC/kcpH+wsjYljSucAWSa8mvyBOkXwprIXSZOsN1sssN+O7Lh6LOHGkX2vwuTuotYO6e/VIOg1YDjyetU5EDCf3B4D1wIJ638+sVfLqzePuoNYuGunO+cfAqxGxL22hpDMlnTX6GFgE7Grg/cwasmHbEAvvfYaLVv+Ehfc+U3Xg5tWbx91BrV1U053zMeBKYIakfcAdEfEgcAPjmnkkzQQeiIilwHnA+tL5X04D/ioifppv9c1OqtSM0siVpmf3dHNo5Ehq+cJ7n6m62abSF4ibgKyZJgz+iLgxo/yfpZQNA0uTx28AlzVYP7OqTBTslY62KwXshm1DHP6/p4b+NOA3Hx498YVQzRdJVt/3s3u6PfyBNZWv3LWOMFEzSr3NNXc++QrHU7oaHAeOHBu7YKJmm6t+ty+1/NDIETcBWVN5kDbrCBMFey1XmpY3u9TavSyt2eaq3+3jZ68eTH3/SoYOjTD/rs3c8elLfeRvuXLwW0eYKNhXLZ4zpjkFTr34asO2IW5fv5PffHjslNepVgArH99+4vnQoRH+6/O/qvv13v3gCKt+tANgzPkKnw+wRjj4rSNMFOyjwbhm0x6GDo3QJZ3SnLLqRztOab5pB0eOBSsf387X/3oHxyIQnPgl4vMBnaHZX+YOfusI5cGe1atnNPRF6SIsOBmcv3XatLYM/XKjdR5fy2pOUlv7asXcBg5+6xhZQwiM/4+VFpzjT65ONbWeP7D2UW+Ps0a4V491vLT/WK228HfOyfX1ukrXy9gU1Irhvh381vHacbz8515PHfC2bqPNQDb1tGK4bzf1WEdJO0mW1eOnk/QXcE6ATlFNj7O8OfitY/ybDTt59PlfjenxsvLx7fR0T2OaSL0QqxMUdU6ATjFRx4TJ4OC3Kaf8qP7snm6kUn/3LCNHjjexds1V9DkBOkWz5zZw8NuUMr6HTtrgaUUx/Yxunlt9daurYVOQT+7alNKOPXRapdKvHLNKfMRvU0o79tBpBx7GwWrh4LemajSgitBDp1atuPLTprYJm3okPSTpgKRdZWV3ShqStD25Lc3YdomkPZJek7Q6z4rb1JPH1IOrFs+hp7tr8irZRnp7uunt6c5cPi25Zssze1mtqmnjfxhYklL+5xExL7ltHL9QUhfwXeAa4BLgRkmXNFJZm9ryCKhl8/u5Z/ncjr9StXuauPO6Sznzt7J/lP/JJ2YDrbny06a2CYM/Ip4F6rnMcAHwWkS8EREfAj8Erq/jdaxD5BVQy+b3c7zDr1Q9cjxODCqX5d8tmwu05spPm9oa6dVzq6SXk6ag6SnL+4G3yp7vS8qsoPIIqA3bhph/1+aaJ0iZCqaf0U33tJO/ZEZHEk1TfqVuWvOXL+qySuo9ufs94G5KAx3eDXwbuHncOmn7bOb/V0krgBUAs2fPrrNa1s5qvTQ9bSarx198q+2HT67X4ZGjp4y5EzBm/H2S5+XTOLbiyk+b2uoK/oh4e/SxpO8Df5Oy2j5gVtnzC4DhCq+5FlgLMDAw0Jn/swuu2oDasG2IO598ZczFWUOHRsYMx9CJsgZaGx/+ATyxdYiBj51z4rNr9pWfNrXVFfySzo+I/cnTzwC7UlZ7EbhY0kXAEHAD8Cd11dKmnKxumxMF1PiuieU6OfQr6ZJO+VLw5CvWiAmDX9JjwJXADEn7gDuAKyXNo/R/8U3gK8m6M4EHImJpRByVdCuwCegCHoqIVyblr7C20ki/cl+ZO1ZPd1fm5+FeO1avCYM/Im5MKX4wY91hYGnZ843AKV09rbM1MqOQw+ykLol7ls/N7N3jXjtWL4/VY7lrpNtmpTDr7J77Y/V0d/Htz13Gsvn97rVjuXPwW+4a6baZdWXu9DO6C9PG39/bwz3L5445cXvP8rn09/aglOVmtfJYPZa7RmYUqtTzZ+G9z3T8OD2C1KGW3WvH8uTgt9w12q88K+TSvlA6jdvtrRkc/DYpJuMIdfwXyrSUbo5TmdvtrVkc/DallH+hVOrzP9V4CkVrJge/TVnjfwFUOvbvniaOtMFs69OA8hmAu6eJNZ+9zIFvTeVePTYpNmwbYuG9z3DR6p+w8N5nahpzvxbL5vfz3Oqr+d/3/pMxA5eV6+/tYc1nL2v5UM7Tz+jmvs/PG9M7x6FvreAjfstd3jNCVTtrV6XeRMvm9/OvHt9e519Um57uaRw9HmMGk+vp7uKOT1/q3jnWFhRteHJsYGAgBgcHW10Nq1NWt8v+3p7UroqVpLXj93R3ZfZjL/+SOLunmyPHjvObD0vbStCs3X0acPYZ3Rz64IhHy7SmkLQ1IgaqWddH/Ja7PGeEqnX4h9Ej6g3bhlj1ox1jjrrzDv3+3h5+8/+OjhlFdNRx4IzTT2Pbtxbl+6ZmOXAbv+Uuzxmhsr4sJrqQa82mPZnj9ufV1v/c6qt5LyX0R3ncIWtXDn7LXZ5jy2R9WQgqnjCuFLrHI/iPn5+Xy6Ttlb7MfDGWtSsHv+Uuz7FlVi2ekzmVW6VJ2icK5LQ69vZ0V12vM0/vOlG/7q5Ta9g9Tb4Yy9qW2/htUuTVe2XZ/H5WZvTGqdTcs2rxnFPa+GFsII+vYy0XhHV3TTvxGgD/9r+/wrsflJp9enu6ufO6S30y19pWNROxPARcCxyIiI8nZWuATwMfAq8DX4qIQynbvgm8DxwDjlZ7xtmsXH9vT2rIjzb3ZJ3kBfjmupf54EjpkikJPr9gVmYgp40xlPXlUt627y6aNtVU09TzMLBkXNkW4OMR8fvAL4DbKmx/VUTMc+hbvept7imtc3LLiNJctZXODZRfEPbc6qszLwpz+71NZRMGf0Q8C7wzrmxzRBxNnj5PaSJ1s0mxbH5/5nAMlU7iVuoKWi1PgmKdKI+TuzcDT2UsC2CzpK2SVuTwXlZQ9Rx5ZzXT1DKmvydBsU7U0MldSbcDR4FHM1ZZGBHDks4Ftkh6NfkFkfZaK4AVALNnz26kWtaB6pncpStj2OasfvxZQ0O4Dd86Td1H/JJuonTS9wuRMe5DMvk6EXEAWA8syHq9iFgbEQMRMdDX11dvtaxD1XPknTVWf1r5aI+eoWSUz9HxhSZrcDmzVqrriF/SEuAbwD+OiA8y1jkTmBYR7yePFwF31V1TK7xaj7yzegOlNRvVOjSE2VQ24RG/pMeA/wnMkbRP0peB7wBnUWq+2S7p/mTdmZI2JpueB/xc0g7gb4GfRMRPJ+WvMEtRy4nZPMcXMmt3Ex7xR8SNKcUPZqw7DCxNHr8BXNZQ7cwaUMvcv1l99t1t0zqRr9y1jlZt81A9J4/NpioHvxm1/Towm+oc/GYJd9u0ovDonGZmBePgNzMrGAe/mVnBOPjNzArGwW9mVjAOfjOzgnHwm5kVjIPfzKxgHPxmZgXjK3fNrK1kTYhj+XHwm1nbGJ0QZ3SwvNEJcQCHf47c1GNmbaPShDiWHwe/mbUNT4jTHFUFv6SHJB2QtKus7BxJWyTtTe6nZ2x7U7LO3mSeXjOzVFkT33hCnHxVe8T/MLBkXNlq4OmIuBh4Onk+hqRzgDuAT1CaaP2OrC8IM7Napsu0+lUV/BHxLPDOuOLrgUeSx48Ay1I2XQxsiYh3IuJdYAunfoGYmQGlE7j3LJ9Lf28PAvp7e7hn+Vyf2M1ZI716zouI/QARsV/SuSnr9ANvlT3fl5SdQtIKYAXA7NmzG6iWmU1lnhBn8k32yV2llEXaihGxNiIGImKgr69vkqtlZlZcjQT/25LOB0juD6Sssw+YVfb8AmC4gfc0M7MGNRL8TwKjvXRuAn6css4mYJGk6clJ3UVJmZmZtUhVbfySHgOuBGZI2kepp869wF9L+jLwK+CzyboDwC0R8c8j4h1JdwMvJi91V0SMP0lsNuk8DIDZSYpIbXJvqYGBgRgcHGx1NaxDjB8GAEpdBN1bxDqJpK0RMVDNur5y1zqehwEwG8vBbx3PwwCYjeXgt47nYQDMxnLwW8fzMABmY3k8fut4oydw3avHrMTBb4XgYQDMTnJTj5lZwTj4zcwKxsFvZlYwDn4zs4Jx8JuZFYyD38ysYBz8ZmYF4+A3MysYB7+ZWcHUHfyS5kjaXnY7LGnluHWulPRe2TrfarzKZmbWiLqHbIiIPcA8AEldwBCwPmXV/xER19b7PmZmlq+8mnr+CHg9In6Z0+uZmdkkySv4bwAey1j2B5J2SHpK0qU5vZ+ZmdWp4eCXdDpwHfDfUha/BHwsIi4D/hOwocLrrJA0KGnw4MGDjVbLzMwy5HHEfw3wUkS8PX5BRByOiF8njzcC3ZJmpL1IRKyNiIGIGOjr68uhWmZmliaP4L+RjGYeSR+VpOTxguT9/j6H9zQzszo1NBGLpDOATwFfKSu7BSAi7gf+KfAvJR0FRoAbIiIaeU8zM2tMQ8EfER8Avz2u7P6yx98BvtPIe5iZWb585a6ZWcE4+M3MCsbBb2ZWMA5+M7OCcfCbmRWMg9/MrGAc/GZmBePgNzMrGAe/mVnBOPjNzArGwW9mVjAOfjOzgnHwm5kVjIPfzKxgHPxmZgXj4DczK5g8Jlt/U9JOSdslDaYsl6S/lPSapJclXd7oe5qZWf0amoGrzFUR8XcZy64BLk5unwC+l9ybmVkLNKOp53rgB1HyPNAr6fwmvK+ZmaXII/gD2Cxpq6QVKcv7gbfKnu9LysaQtELSoKTBgwcP5lAtMzNLk0fwL4yIyyk16XxV0ifHLVfKNnFKQcTaiBiIiIG+vr4cqmVmZmkaDv6IGE7uDwDrgQXjVtkHzCp7fgEw3Oj7mplZfRoKfklnSjpr9DGwCNg1brUngT9NevdcAbwXEfsbeV8zM6tfo716zgPWSxp9rb+KiJ9KugUgIu4HNgJLgdeAD4AvNfieZmbWgIaCPyLeAC5LKb+/7HEAX23kfczMLD959eM3M+s4G7YNsWbTHoYPjTCzt4dVi+ewbP4pnRKnHAe/mVmKDduGuG3dTkaOHANg6NAIt63bCTDlw99j9ZiZpVizac+J0B81cuQYazbtaVGN8uPgNzNLMXxopKbyqcTBb2aWYmZvT03lU4mD38wsxarFc+jp7hpT1tPdxarFc1pUo/z45K6ZWYrRE7ju1WNmViDL5vd3RNCP56YeM7OCcfCbmRWMg9/MrGAc/GZmBePgNzMrGJUGz2wvkg4Cv8zp5WYAWRPBtwvXMR+uY+PavX7gOmb5WERUNX1hWwZ/niQNRsRAq+tRieuYD9exce1eP3Ad8+CmHjOzgnHwm5kVTBGCf22rK1AF1zEfrmPj2r1+4Do2rOPb+M3MbKwiHPGbmVmZjgl+SW9K2ilpu6TBlOWS9JeSXpP0sqTLm1y/OUndRm+HJa0ct86Vkt4rW+dbTajXQ5IOSNpVVnaOpC2S9ib30zO2vSlZZ6+km5pcxzWSXk3+LddL6s3YtuJ+Mcl1vFPSUNm/59KMbZdI2pPsm6ubWL/Hy+r2pqTtGds26zOcJelnknZLekXS15LyttkfK9SxrfbHCUVER9yAN4EZFZYvBZ4CBFwBvNDCunYB/4dSv9vy8iuBv2lyXT4JXA7sKiv798Dq5PFq4M9StjsHeCO5n548nt7EOi4CTkse/1laHavZLya5jncC/7qKfeF14B8ApwM7gEuaUb9xy78NfKvFn+H5wOXJ47OAXwCXtNP+WKGObbU/TnTrmCP+KlwP/CBKngd6JZ3forr8EfB6ROR1kVrdIuJZ4J1xxdcDjySPHwGWpWy6GNgSEe9ExLvAFmBJs+oYEZsj4mjy9Hnggsl472plfI7VWAC8FhFvRMSHwA8pff65qlQ/SQI+BzyW9/vWIiL2R8RLyeP3gd1AP220P2bVsd32x4l0UvAHsFnSVkkrUpb3A2+VPd+XlLXCDWT/J/sDSTskPSXp0mZWqsx5EbEfSjs6cG7KOu30ed5M6ddcmon2i8l2a/Lz/6GMJop2+Bz/EHg7IvZmLG/6ZyjpQmA+8AJtuj+Oq2O5dt4fgc6aiGVhRAxLOhfYIunV5ChnlFK2aXqXJkmnA9cBt6UsfolS88+vk/bgDcDFzaxfDdrl87wdOAo8mrHKRPvFZPoecDelz+VuSs0pN49bpx0+xxupfLTf1M9Q0keAJ4CVEXG49INk4s1Syibtcxxfx7Lydt4fT+iYI/6IGE7uDwDrKf2ELrcPmFX2/AJguDm1G+Ma4KWIeHv8gog4HBG/Th5vBLolzWh2BYG3R5vBkvsDKeu0/PNMTuBdC3whkgbU8arYLyZNRLwdEcci4jjw/Yz3bunnKOk0YDnweNY6zfwMJXVTCtRHI2JdUtxW+2NGHdt+fyzXEcEv6UxJZ40+pnSiZde41Z4E/lQlVwDvjf58bLLMoytJH03aW5G0gNK/z983sW6jngRGe0XcBPw4ZZ1NwCJJ05MmjEVJWVNIWgJ8A7guIj7IWKea/WIy61h+DukzGe/9InCxpIuSX4M3UPr8m+WPgVcjYl/awmZ+hsm+/yCwOyLuK1vUNvtjVh2nwv44RqvPLudxo9QjYkdyewW4PSm/BbgleSzgu5R6UOwEBlpQzzMoBfnZZWXldbw1qf8OSieI/mET6vQYsB84Qumo6cvAbwNPA3uT+3OSdQeAB8q2vRl4Lbl9qcl1fI1Sm+725HZ/su5MYGOl/aKJdfwvyb72MqXwOn98HZPnSyn1Dnl9suqYVr+k/OHR/a9s3VZ9hv+IUvPMy2X/rkvbaX+sUMe22h8nuvnKXTOzgumIph4zM6ueg9/MrGAc/GZmBePgNzMrGAe/mVnBOPjNzArGwW9mVjAOfjOzgvn/BcjUpDem9WsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.scatter(X[:,0], X[:,1])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "从图中，可以看出，这是一个非常紧密的聚类，几个值远离了聚类。\n",
    "在这个简单的例子中，这些可以被认为是异常的。 "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.1 高斯分布\n",
    "\n",
    "想要做异常检测，你首先需要拟合一个针对该数据的模型。\n",
    "\n",
    "对于给定的数据集$\\{x^{(1)}, ..., x^{(m)}\\}$（其中$x^{(1)} \\in \\mathbb{R}^n$），你想为数据中的每个特征$i = 1, ..., n$估计高斯分布。\n",
    "为此，我们将创建一个返回每个要素的均值($\\mu_i$)和方差的函数$\\sigma^2_i$。\n",
    "\n",
    "高斯分布的公式为：\n",
    "$$\n",
    "p(x; \\mu, \\sigma^2) = \\dfrac{1}{\\sqrt{2\\pi\\sigma^2}} e^{-\\dfrac{(x-\\mu)^2}{2\\sigma^2}}\n",
    "$$\n",
    "\n",
    "其中为$\\mu$为均值，$\\sigma^2$为方差。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.2 估算高斯分布的参数\n",
    "\n",
    "均值($\\mu_i$)公式为：\n",
    "$$\\mu_i = \\dfrac{1}{m}\\sum_{j=1}^m x_i^{(j)}$$\n",
    "\n",
    "方差$\\sigma^2_i$公式为：\n",
    "$$\\sigma^2_i = \\dfrac{1}{m}\\sum_{j=1}^m(x_i^{(j)} - \\mu_i)^2$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def estimate_gaussian(X):\n",
    "    mu = X.mean(axis=0)\n",
    "    sigma = X.var(axis=0)\n",
    "\n",
    "    return mu, sigma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([14.11222578, 14.99771051]), array([1.83263141, 1.70974533]))"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mu, sigma = estimate_gaussian(X)\n",
    "mu, sigma"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "接下来绘制拟合的效果。\n",
    "\n",
    "因为之后将使用多变量高斯模型，因此将计算协方差矩阵，不再使用sigma。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3X2QHPdd5/H3t3se9lm7Wmul1a5j2fJDZGJLwTonEO5wjENxueKcQEJh58EXAg51BJJLDkhxfxDgqIKrkBwVKFPOEw5n5wGcYAe4I8E4GCdEjuVYlmzFD7Jle1crrZ5Wu9qHeej+3h/d215Ju6tZWbOzs/q8qqZmpqdn5tseb3/Uv9+vf23ujoiICEDQ6AJERGTlUCiIiEhGoSAiIhmFgoiIZBQKIiKSUSiIiEimbqFgZi1m9oiZ7TKzJ83s99Lll5rZDjN71sy+YmaFetUgIiJLU88jhRJwo7tvBbYBP2NmbwT+GPiUu18BHAfeX8caRERkCeoWCp44mT7NpzcHbgT+Jl1+F/C2etUgIiJLk6vnh5tZCOwELgf+HNgHjLl7NV1lCBhY4L23A7cDtLe1XffaKzbXs1QRkVVn567dR9x93VLeU9dQcPcI2GZm3cDXgS3zrbbAe+8E7gTYvu1a//43769bnSIiq1Gw/tIXl/yeehRyOncfA74NvBHoNrPZMBoEDixHDSIicnb1HH20Lj1CwMxagZuAvcCDwDvS1W4D7qtXDSIisjT1bD7qB+5K+xUC4Kvu/ndm9hTwZTP7n8APgM/VsQYREVmCuoWCuz8BvH6e5c8D19fre0VE5NzpjGYREckoFEREJKNQEBGRjEJBREQyCgUREckoFEREJKNQEBGRjEJBREQyCgUREckoFEREJKNQEBGRjEJBREQyCgUREckoFEREJKNQEBGRjEJBREQyCgUREckoFEREJKNQEBGRjEJBREQyCgUREckoFEREJKNQEBGRjEJBREQyCgUREckoFEREJKNQEBGRTN1CwcwuNrMHzWyvmT1pZh9Kl3/czIbN7PH09tZ61SAiIkuTq+NnV4GPuvtjZtYJ7DSzb6WvfcrdP1HH7xYRkXNQt1Bw9xFgJH08YWZ7gYF6fZ+IiLx6y9KnYGabgNcDO9JFHzSzJ8zs82bWsxw1iIjI2dU9FMysA7gX+LC7jwN3AJuBbSRHEn+ywPtuN7NHzezRw0eP1rtMERGhzqFgZnmSQLjb3b8G4O6H3D1y9xj4DHD9fO919zvdfbu7b1/X21vPMkVEJFXP0UcGfA7Y6+6fnLO8f85qbwf21KuGlWpovMwjQ5MMjZcbXYqIyCnqOfroTcB7gN1m9ni67HeAW8xsG+DAfuADdaxhxRkaL3PPruOEAUQx3Lq1h8GuQqPLEhEB6jv66GHA5nnpH+r1nc3gwHiFMID+zjwjExUOjFcUCiKyYuiM5mW2sStPFMPIRIUoTp6LiKwU9Ww+knkMdhW4dWsPB8YrbOzK6yhBRFYUhUIDDHYVFAYisiKp+UhERDIKBRERySgUREQko1AQEZGMQkFERDIKBRERySgUREQko1AQEZGMQkFERDIKBRERySgUREQko1AQEZGMQkFERDIKBRERySgUREQko1AQEZGMQkFERDIKBRERySgUREQko1AQEZGMQkFERDIKBRERySgUREQko1AQEZGMQkFERDJ1CwUzu9jMHjSzvWb2pJl9KF2+1sy+ZWbPpvc99apBRESWpp5HClXgo+6+BXgj8GtmdjXwMeABd78CeCB9LiIiK0DdQsHdR9z9sfTxBLAXGABuBu5KV7sLeFu9ahARkaVZlj4FM9sEvB7YAax39xFIggPoW+A9t5vZo2b26OGjR5ejTBGRC17dQ8HMOoB7gQ+7+3it73P3O919u7tvX9fbW78CRUQkU9dQMLM8SSDc7e5fSxcfMrP+9PV+YLSeNYiISO3qOfrIgM8Be939k3Neuh+4LX18G3BfvWoQEZGlydXxs98EvAfYbWaPp8t+B/gj4Ktm9n7gJeCddaxBRESWoG6h4O4PA7bAyz9Vr+8VEZFzpzOaRUQko1AQEZGMQkFERDJnDQUze6CWZSIi0vwW7Gg2sxagDbgonbRuttO4C9i4DLWJiMgyW2z00QeAD5MEwGNzlo8Df17PokREpDEWDAV3/1PgT83s193908tYk4ico6HxMgfGK2zsyjPYVWh0OdKEajlP4YSZvff0he7+xTrUI7IiNOPOdWi8zD27jhMGEMVw69aepqldVo5aQuHfzXncQnLi2WOAQkFWpWbduR4YrxAG0N+ZZ2SiwoHxSlPULSvLWUPB3X997nMzWwP8Vd0qEmmwZt25buzKE8UwMlEhipPnIkt1LtNcTAFXnO9CRFaKZt25DnYVuHVrT9M1e8nKctZQMLNvAJ4+DYEtwFfrWZRIIzXzznWwq9BU9crKU8uRwifmPK4CL7r7UJ3qETnvfN+uJa0/PJXcBtpgoAR++NTXbfPWU543Y6e0yEJq6VP4FzPbAFxPcsSwr+5VibxKSw2CWcNTcPcLAbnAqcbGuy6NGWhb+LOHp+CesQ1N1yktspBaprn4ZeAR4OeAdwDfM7NfqndhIkvl+3Zlt/kMT8EjR5L7hQxPQS5w+luT+8XWnV0/PHGI/plDhEHSSS3SzGppPvpN4PXufhTAzHqB7wKfr2dhIrU6PQROaf5pe2XZQkcA/sSO7L1HT67lByfWc6BvDX0tRmjOI0dO/ay5BtqgGhsj0051cpSNXQcZOg4Hei5Xc5I0pVpCYQiYmPN8Ani5PuWILM18gTDfzn94CsLDQ2zIlRmpFhg6OcbG1olT3rtzupPPjm+k6sahQ9O8t3OEBw63EOJEGLd2HmIgX8KufUP2noE2+KkNMXvG4HXdyXiMu18IyL24L/n+G65SMEhTqSUUhoEdZnYfSZ/CzcAjZvYRgNOuvyyyLBZrIppt/hmZdoZ2P8XG1gk2VopErGekWiDCGMiVT31fpcjfn1xLJTY2F2bYV27hu1Nr6MtXuLIwzUi1wHA12bkP73gqe/+u9dvYcTSgu+D87cvGumLMdOT05eGJcefeh/by85fA4DVbz6hVZCWqJRT2cWrn8n3pfef5L0fk7BbrRJ5tzjmwfyjZ+XcmO++BfIlbOw+xa6YDN8/WH64U+ednTvJwbgM5n2Y0v4bKFBwPjYHKCXZ39QEQufHIVBd7Ku1syJUZj3JMeYAdOcYRL3JdcZwfhgMcKRpDU8ZUBGMleP6ksX8Kbjq+i63btjDYVdBoJVnRahl99HvLUYhILc42qmjjczu4tVJkuLXAQK7MQL4EJDv/x0vt7Biu0O0zPEXAm6sHuT9/JY/kL2EiaOE18Ql+tHIQt5irK0e4Lh5l54kKk+aMe5Hv5NZywopsaCkyE4eYQcFiYjd+WGrHmeCSaIKD1s+JslEI4VgZHj0KhvG9f3qavmLM7mo3GzvztOVDjVY6BwrV+qrl5LUrgf8ObJq7vrvfWL+yRJZutsN4IF/KwgCSQLjnOeOgwUtBL2+OXmTU2vj7/GUctlbWMsMkRUasnYJVuan6As8GF7Ez6GN3uJ6N8Tgvh11MWY6IgKFSQEtcJbKAXD5Pm0VcVxxnZ6mTfzi5liploqDITBUqEZRCaA3h8eMwVgqoMs54bwevXVdsmik0VopmnZeqmdTSfPTXwF8AnwWi+pYjsrDFjhLmjiA6ZfmuXQyF/YThIJt9jJdsDbuDdRwO2tlUPc7+mSKlsMBEPqAlqjBcLfD9k+0Mt7RyLN+OWYWROM90W44KIR1eoeowQ0jojlcqvLZ6kJ4u46VqC5MekjenEJUpe0ghDHHgyTEoR7CuBY6UYWz8JGOduaaZQmOlaNZ5qZpJLaFQdfc76l6JyCKWGgjDlSJDT7/MgHUwEE8QhQGjtLE2mqL34BDdQY6LShOMrNnEgZZuesonMZyJXAv/tvYqZoI8kQXEZrRXZyhUYtqiSczBWlvo9WlOBi2sjaepWsATB2cI2522IKbiULSYSQ+JogiCkJnY2NQBR0qwJgf9rfDui44x2HVxPf+zrTrNOi9VM1nscpxr04ffMLP/CnwdyI7J3f1YnWsTOSc7pzv5Py8VyeVeS9UC3lXew5ur+/nCyU0U42M8376BUpBjZ/dlTAd58h7RUZ3m5dZ1jOdaMAPHMHccpxTksbjCWMsaDCeMY7qmT7C+rUrRq7RT4SQFRsoFiGPiXI4fLZxgvJyniuFxRDEI+JE1jjv0FODG/qRT3PftOmPaDFlYM89L1SwWO1LYSTIEdfbazL855zUHLqtXUSK18id2MFwpMlwtEAKHqgXuG85xMOjgRNjCmniGu/M/wo898whd3dM827GBl1ovYsZyVIOQzqkJTuaKlKoVZjDylQrVfPKvT7ckGMoWMpMvEHiMW0gURBwrdNAzPUnXzAmOd7UxEebBY8pBDoudPZVOOoKIsgdMecjIiQpfmi6wrdfZ1G5A3ND/bs1Mk/7V12KX47x0OQsRWchiTUfDlSL3TKxnMg7YU+6gLyhzMIypYlQ9IDRoGT3C8UI7+1sv4uWWXsZzrcRuRGFI1LYGd2c634JbQDXMQRwDDm7pneHuRASAkY+dchBysKWH0cIaJguttHiFGStgxDgB09WYYmhszM1wOCowQ8hkOeLfDofgnp1xPbt9OlqQlaKW0Uc/N8/iE8Budx89/yWJ1C45QnCKOCFO99ghDgdr6fESURBwWfUYEcbz7RsIPeJ4voPILJn1K46pWggGkUEujshXy8TuWBoQHoR4XIUg2ZljBnHETFhIAscjgjiitTJDtRhQIg84ZsZAvsRPth7nuzPdvFBpISIHDrtPJNNnwJwpOdaV9a9fWRHOOiEe8H6SkUfvSm+fAT4CfMfM3rPQm8zs82Y2amZ75iz7uJkNm9nj6e2tr7J+WaWGxss8MjTJ0Hh50fUGcmUijKNxjrEox5i10EaFy+PjfKj0fa5/9lEKcYVj+Q4qliPnEUEUpVcI8VduFlAN81SDPKE7UZAjazm1EIAwrhJEEZUwaSKaCXJMBgWqQchUWMAip0CZPBE5Yjonx3hH1xHet2aEHE6BmJa4wrXdTuSvTMnx8GHj7m8/fdZtFVkOtYw+ioEt7n4IwMzWA3cAbwAeYuFLc/4l8GeceS3nT7n7J85cXSRxxlj07vkno4PknIQbW49z98R6Ls9P88NgLQUipsI8x4I2Sr3reKF9PccKHQTuOEaAE7tTrFYo5YoQBCTtREBUJQpyyVnPlv6byZJwiIIAcyOMna6ZcSa8nSIxcVwlwMlHFdoMOikzTY6bqi8wkPZPXFs8yd5KO20W01vMM9DmZ0zJoeGVshLUEgqbZgMhNQpc6e7HzGzBeYLd/SEz2/Qq65ML0Olj0ee2v88nAgZzJUoe8DxFQo+JzThAB0HeyXsV3CmbkY8rVD0k9Aprpsc53LEWj0mCwcDzhUVOxgmwuIIBU/lWOkpT9MeTHGrtSZqWPKZ3fIxNHRXavcKW+Cg79gQcvWQLVxanubZlkn3lFt7Q25JtTzbDamwaXikrQi2h8K9m9nckJ7EB/DzwkJm1A2Pn8J0fNLP3Ao8CH3X34/OtZGa3A7cDvGZw4zl8jTSr08eiD2zZAof3Lrj+sSjP4zMddAYRMXAsaCdHTAdlOsNpnmnvpxzmwZ1iVKW1PMlMvoUTrV14mPQBJHxOs5JlRwgZg0JcpW/8CDFGz+QxjvQNMm15ZnJ5+mbG6C8d58biYdb7FA/mNhESM1bqAofQHDdjQ2vycQNt8K5L4zT0nIHDe6FLHc7SWLWEwq+RBMGbSBpZvwjc6+4OvHmJ33cH8Ackf3V/APwJMO8Fe9z9TuBOgO3brvX51pHVab6x6LOXxByegsePJfvrrT0wkk53PRUF7K+28KY1Y4wdG6ObEkWPeKm9I9vlG1DK5anSTugR1SCEOEo+zCxZw2bXhFNGZHtM4DFgTBVaKUQVDq7dyHSujciM1rhCkZiZIM8xa+Vw0Mak5bgqPg5Blb6gzM5SJ2vCKg8cDNjQGmfXaFjsKEhkudUyIZ4Df5PeXpW5zVBm9hng717tZ0rzm2+Cs/nGog9PwV88E/DD8WR3ffUaWDvVRTU2YgsoxwGPz3SwxcdZxzRDdPJi2ENECBg+u68PAiyKyUcRpUJ+zoFC/Eo/ApC0J8XgMSFQqJQxd1oqMxTiChMtHVSCkBijEJcZz7Wwr309o4V+rqweYX/YA0A7Rk+uwsVxif5cmYOBn7VJTKRRahmSOsErfzYFIA9MunvXUr/MzPrdfSR9+nZgz2Lry+pX6wRntnkrwzt2MRU5nflk7z5Zdfo2DHJsX5npOKBCgAXGVJDnomiSg0EbHV5mvU0xQifEMYYRWUCUK2AekytXqOZyyQHB6c1FyTdTjKoUy9Pko2T4KR5zqHM95UIRIyYfR/SVxogsZEthgkPeQa+V6KiOssnHeEtn0vX2VLkjuZ5DbAy06eBXVqZajhROuW6Cmb0NuP5s7zOzLwE3ABeZ2RDwu8ANZraNJGT2Ax9YesmymixlgrOBNmgLjZcqyf9A3QXjUMm4pjDJo9Od9IQVWgInXNPFD8Y2UCZkf9hNiRBwCAI8jtMhpo5bDqNEz9QYx9vWAAbBmf0I1TBHGOSYaO1Kpr6wgFw1GQobuNESldgyPsSJfDuVYivTlqfkAe1UecsV7dmMrbd2HmK4WmDw0o06SpAVq5Y+hVO4+9+a2cdqWO+WeRZ/bqnfJ6vbUiY4G7xmK7/KrqxPwR2enjCuvnItPHOM3aV2IofnKq0MdK3j+GSZQlyhZCF5i4kxIgxwCnFE1QI8NCqFljQMzjxSCKoV2krTTBdbkzObgzA56S2Xp+AROY/JufNy2zomurq5LB7j8vgY2+MRtkWjDOQ3Z581O6W3LRAIOqtZVoKlntEcANt5pTlJ5FVZ6gRncztmh6fgyRPJkM7i2l6uHxulN6jycrXAnnIHcbGNXKVKjogyOWIMgqTzuGoBocdcMnWYAy09WDrqyM1fmZYogFxgTLe2ERO80t8QQBA7sQXEHlMOc5S6ujhpBbp9hnU+zTqfZvDazWfUP/f6ziIrUS1HCj8753GVpNnn5rpUIxekpUxwZpu3ZnMhzR3SGfY5DxzsIzw8RNGcy/PTDFVayBWge+oEx4JWwBgPiqyJp5m2HP0nRxksjdFTmWJX96WYx8QWEMZVcnFEHIZY7DgwnSsSxhFREFKolOmMS1w2OcrMmi5OBgW6vMRRa2WMImspMXjVxcyZVFikadTSp/C+5ShEpFanB8PskcOG1pjh9Rs5+syzfHZ8I21BxDQB790wztOViCeOQiUIWedTjFgHYx09HOm8iKJX6ahMsrYyAW5M51oYy7dTsjxB6AxOH+EoXQREBFVnw8wYQXsrUXcnl0UnklkyzFkfT3JjtJ/XX3nRKVd+y+pe5ChBTUeyUtTSfDQIfJrkPAUHHgY+5O5Dda5NLhDncs3ducEwazYgHuEKrjlgFI6MUMK4ojDNawvTjFb7KVUiDpfamLaQdo8pWciMhXTmy2wIneviYf5v7nKcCnmDgle5ojDJtYxziY+xNp7hcEs7LVGVfUEPN1ZfYGt8mOGgk4F4Im0yWlogiKwktTQffQG4B3hn+vzd6bK31KsouXCc6zV3h8bLHChezsbjz50xkmd2lFKuv5/coQMM5MoMVwtcWZjm0nyJb8S95GJjOmilEkWEBjmcF4M1dPgMffEkHsDaeIbNfpxt8UHMoYcZNsSTHAw6mSHHhniSrfFhBq/dzCAAffPWerZA0FGCrCS1hMI6d//CnOd/aWYfrldBcmE5l2vunhokG7j1kp5kiojU3L6GY139PDDTT9+hZ4kwDkd5esMqfbmYY1GeE4QEZkyHrUx5QIXkqms/WX2Zn4heZr1PcX/+SvYGvWCwpXqUmzc5EckMrXNHF52ulqMDBYKsNLWEwhEzezfwpfT5LcDR+pUkF5JzuebubJC05gP2HS2x6+AUg1ee2pw00AYHp+EzzyUX2on8Kn7lipgYaHt6lBDnRJxjc26Ke0/2MRHncOBkrg0wvlO8jHesi4mAqYk+uuIcuDEZrCNilOtbJxasr9amIgWCrES1hMIvkUyB/SmSPoXvssB8RSJLdS7X3N3YlWdsJuLfXpoE4HsvTbF1QxuD6U52Nhz2jEFoMNAOw5NwaAZu2wxbe/qyaSaGdh/icFxgLMrzz1PdlAkomNNizp5SGze1j9FmES/GLQAMhDEDuTOve7CUPgOFgaxki4aCmYXAz7v7f16meuQCtNRr7g52FXjjYDuT5YjNa4vMVE+9FsHsTvd1R3dx/1ASCJHD67qT958yCd01V9P+QsCawPmRI8aeMSiGMBGHrM9VGMiX+ED3AXbNdODmbCtOMnjdtnPeVgWCrHSLhoK7R2Z2M8lRgsiKsbW/le+9PMm+YyXa8+G8zU7br9/Kr3Sf4JEn93N9L1zXe+bnzO1/uKrTefCQUY6hEOTpveQK7CIYhLQj+dwpDKRZ1NJ89B0z+zPgK8Dk7EJ3f6xuVYnUwsCy6a7PNDReZvehGTr6NrA7hmvWndohPWv2yGH2DOlckFz05tVOWqcgkGZUSyj8eHr/+3OWOXDj+S9HpDYHxit0t4RsWdey4KilhTqk5zq9c/qUi96cw6R1CgJpdrWc0bzUC+mI1F0to5YW7JCeEx7D67ac0sk9t6no9JPjZmnHL6tZLWc0F0muvLZp7vru/vsLvUek3moZtXS2DumznTinnb9ciGppProPOAHsRDN8yQpSy6ilrf2tPDk6w0zVzziiOJcT50RWu1pCYdDdf6bulYjUwWJHFOdy4pzIaldLKHzXzK5x9911r0akDhY6ojiXE+dEVrsFQ8HM9pBcbiQHvM/MnidpPjLA3f3a5SlRpH6WeuKcyGq32JHCAHDup26KiEjTWSwUXnD3F5etEhERabjFQqHPzD6y0Ivu/sk61CMiIg20WCiEQAcLTiIgIiKrzWKhMKIT1ERELizBIq/pCEFE5AKzWCj81LJVISIiK8KCzUfufmw5CxGRxhoaL+tEPqnpjGYRWeXONjmgXDgWaz56Vczs82Y2mp4ZPbtsrZl9y8yeTe976vX9IlK7uZMDhkHyXC5MdQsF4C+B0yfS+xjwgLtfATyQPheRBtPkgDKrbs1H7v6QmW06bfHNwA3p47uAbwO/Xa8aRKQ2mhxQZi13n8J6dx8BcPcRM+tbaEUzux24HeA1gxuXqTyRC5cmBxSob/PRq+Lud7r7dnffvq63t9HliIhcEJY7FA6ZWT9Aej+6zN8vIiKLWO5QuB+4LX18G8mlPkXqbmi8zCNDkwyNlxtdisiKVrc+BTP7Ekmn8kVmNgT8LvBHwFfN7P3AS8A76/X9IrM0Bl+kdvUcfXTLAi9p+gxZVnPH4I9MVDgwXlEoiCxgxXY0i5wvGoMvUjtNcyGrnsbgi9ROoSAXBI3BF6mNmo9ERCSjUBARkYxCQUREMgoFERHJKBRERCSjUBARkYxCQUREMgoFERHJKBRERCSjUBARkYxCQUREMgoFERHJKBRERCSjUBARkYxCQUREMgoFERHJKBRERCSjUBARkYxCQUREMgoFERHJKBRERCSjUBARkYxCQUREMgoFERHJ5BrxpWa2H5gAIqDq7tsbUYeIiJyqIaGQerO7H2ng94uIyGnUfCQiIplGhYID3zSznWZ2e4NqEBGR0zSq+ehN7n7AzPqAb5nZD939obkrpGFxO8BrBjc2okYRkQtOQ44U3P1Aej8KfB24fp517nT37e6+fV1v73KXKCJyQVr2UDCzdjPrnH0M/DSwZ7nrEBGRMzWi+Wg98HUzm/3+e9z9/zWgDhEROc2yh4K7Pw9sXe7vFZHGGhovc2C8wsauPINdhUaXIwto5HkKInKBGBovc8+u44QBRDHcurVHwbBC6TwFEam7A+MVwgD6O/OEQfJcViaFgojU3cauPFEMIxMVojh5LiuTmo9EpO4GuwrcurVHfQpNQKEgIstisKugMGgCaj4SEZGMQkFERDIKBRERySgUREQko1AQEZGMQkFERDIKBRERySgUREQko1AQEZGMQkFERDIKBRERySgUREQko1AQEZGMQkFERDIKBRERySgUREQko1AQEZGMQkFERDIKBRERySgUREQko1AQEZGMQkFERDIKBRERyTQkFMzsZ8zsaTN7zsw+1ogaRETkTMseCmYWAn8O/EfgauAWM7t6uesQEZEzNeJI4XrgOXd/3t3LwJeBmxtQh4iInCbXgO8cAF6e83wIeMPpK5nZ7cDt6dNSsP7SPctQW6NcBBxpdBF1tJq3bzVvG2j7mt1VS31DI0LB5lnmZyxwvxO4E8DMHnX37fUurFG0fc1rNW8baPuanZk9utT3NKL5aAi4eM7zQeBAA+oQEZHTNCIUvg9cYWaXmlkB+EXg/gbUISIip1n25iN3r5rZB4F/BELg8+7+5Fnedmf9K2sobV/zWs3bBtq+Zrfk7TP3M5rzRUTkAqUzmkVEJKNQEBGRzIoOhdU+HYaZ7Tez3Wb2+LkMHVtpzOzzZjZqZnvmLFtrZt8ys2fT+55G1vhqLLB9Hzez4fQ3fNzM3trIGl8NM7vYzB40s71m9qSZfShd3vS/4SLbtip+PzNrMbNHzGxXun2/ly6/1Mx2pL/dV9LBPYt/1krtU0inw3gGeAvJMNbvA7e4+1MNLew8MrP9wHZ3XxUnz5jZfwBOAl9099ely/4XcMzd/ygN9h53/+1G1nmuFti+jwMn3f0TjaztfDCzfqDf3R8zs05gJ/A24L/Q5L/hItv2C6yC38/MDGh395NmlgceBj4EfAT4mrt/2cz+Atjl7ncs9lkr+UhB02E0GXd/CDh22uKbgbvSx3eR/CE2pQW2b9Vw9xF3fyx9PAHsJZmBoOl/w0W2bVXwxMn0aT69OXAj8Dfp8pp+u5UcCvNNh7FqfsSUA980s53ptB6r0Xp3H4HkDxPoa3A99fBBM3sibV5quqaV+ZjZJuD1wA5W2W942rbBKvn9zCw0s8eBUeBbwD5gzN2r6So17UNXcijUNB1Gk3uTu/8oyYyxv5Y2T0hzuQPYDGwDRoA/aWw5r54/vTpnAAADb0lEQVSZdQD3Ah929/FG13M+zbNtq+b3c/fI3beRzBJxPbBlvtXO9jkrORRW/XQY7n4gvR8Fvk7yQ642h9L23Nl23dEG13Neufuh9I8xBj5Dk/+GaXv0vcDd7v61dPGq+A3n27bV9vsBuPsY8G3gjUC3mc2epFzTPnQlh8Kqng7DzNrTDi/MrB34aWA1zgR7P3Bb+vg24L4G1nLeze4sU2+niX/DtLPyc8Bed//knJea/jdcaNtWy+9nZuvMrDt93ArcRNJv8iDwjnS1mn67FTv6CCAdHva/eWU6jD9scEnnjZldRnJ0AMl0I/c0+/aZ2ZeAG0imIz4E/C7wt8BXgdcALwHvdPem7KxdYPtuIGl6cGA/8IHZ9vdmY2Y/AfwrsBuI08W/Q9L23tS/4SLbdgur4Pczs2tJOpJDkn/sf9Xdfz/dz3wZWAv8AHi3u5cW/ayVHAoiIrK8VnLzkYiILDOFgoiIZBQKIiKSUSiIiEhGoSAiIhmFglyQzOzk2dfK1r3BzH68nvWIrBQKBZGzuwFQKMgFQaEgkjKzn03nnv+Bmf2Tma1PJ0/7VeC/pfPt//v07NF7zez76e1N6fs/nk6q9m0ze97MfmPOZ783nXRtl5n9lZl1mtkL6dQLmFmXJdfXyDdi20Vm5c6+isgF42Hgje7uZvbLwG+5+0fTeeizOffN7B7gU+7+sJm9BvhHXpl87LXAm4FO4GkzuwO4EvgfJBMgHjGzte4+YWbfBv4TyVnfvwjc6+6V5dtckTMpFEReMQh8JZ0PpwC8sMB6NwFXJ9PpANA1O48V8PfpNAIlMxsF1pPOaT97MaU5U0R8FvgtklB4H/Ar53l7RJZMoSDyik8Dn3T3+83sBuDjC6wXAD/m7tNzF6YhMXdemYjkb8yYZ8pid/+OmW0ys58EQndvysnYZHVRn4LIK9YAw+nj2+YsnyBpDpr1TeCDs0/MbNtZPvcB4BfMrDddf+2c174IfAn4wjnWLHJeKRTkQtVmZkNzbh8hOTL4azP7V2DudbO/Abx9tqMZ+A1ge9px/BRJR/SC3P1J4A+BfzGzXcDcaanvBnpIgkGk4TRLqkgDmdk7gJvd/T2NrkUE1Kcg0jBm9mmSS7G+tdG1iMzSkYKIiGTUpyAiIhmFgoiIZBQKIiKSUSiIiEhGoSAiIpn/D2q7lvBB2fBQAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 创建多变量高斯模型\n",
    "mu = X.mean(axis=0)\n",
    "cov = np.cov(X.T)\n",
    "multi_normal = stats.multivariate_normal(mu, cov)\n",
    "\n",
    "# 创建 grid\n",
    "x, y = np.mgrid[0:30:0.01, 0:30:0.01]\n",
    "pos = np.dstack((x, y))\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "# 绘制概率密度\n",
    "ax.contourf(x, y, multi_normal.pdf(pos), cmap='Reds')\n",
    "\n",
    "# 绘制原始数据点\n",
    "sb.regplot('Latency', 'Throughput',\n",
    "           data=pd.DataFrame(X, columns=['Latency', 'Throughput']), \n",
    "           fit_reg=False,\n",
    "           ax=ax,\n",
    "           scatter_kws={\"s\":10,\n",
    "                        \"alpha\":0.4})\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.3 选择阈值$\\epsilon$\n",
    "\n",
    "现在我们有了高斯分布的参数。对于给定的分布，你可以评估哪些样本有较高的概率，哪些有较低的概率。\n",
    "低概率的样本意味着更可能是数据集中的异常。\n",
    "\n",
    "决定一个哪些样本是异常的一个办法是基于交叉验证集选择一个阈值。\n",
    "当一个样本的概率小于阈值的时候，这个样本应该被认为是一个异常。\n",
    "\n",
    "在这部分练习中，你将实现一个基于交叉验证集$F_1\\text{ Score}$选择阈值$\\epsilon$的算法。\n",
    "\n",
    "为此，我们需要使用一组标记的验证数据（其中真实异常样本已被标记），并在给出不同阈值的情况下，对模型的性能进行鉴定。\n",
    "\n",
    "首先，加载数据集："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((307, 2), (307, 1))"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Xval = data['Xval']\n",
    "yval = data['yval']\n",
    "\n",
    "Xval.shape, yval.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们还需要为验证集（使用相同的模型参数）执行此操作。 \n",
    "\n",
    "我们将使用与真实标签组合的这些概率来确定将数据点分配为异常的最佳概率阈值。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu = X.mean(axis=0)\n",
    "cov = np.cov(X.T)\n",
    "\n",
    "multi_normal = stats.multivariate_normal(mu, cov)\n",
    "\n",
    "pval = multi_normal.pdf(Xval)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "接下来，我们需要一个函数，找到给定概率密度值和真实标签的最佳阈值。 为了做到这一点，我们将为不同的epsilon值计算F1分数。 F1是真阳性，假阳性和假阴性的数量的函数。\n",
    "\n",
    "方程式如下：\n",
    "$$F_1 = \\dfrac{2 \\cdot prec \\cdot rec}{prec + rec}$$\n",
    "\n",
    "查准率 Precision：\n",
    "$$prec = \\dfrac{tp}{tp + fp}$$\n",
    "\n",
    "查全率 Recall：\n",
    "$$rec = \\dfrac{tp}{tp + fn}$$\n",
    "\n",
    "其中：\n",
    "* $tp$ 是真正例（true positive）\n",
    "* $fp$ 是假正例（false positive）\n",
    "* $fn$ 是假反例（false nagative）"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "计算F1 Score的函数。\n",
    "也可以使用sklearn.metrics.f1_score计算F1 Score，下面是一种实现。\n",
    "\"\"\"\n",
    "def f1_score(yval, preds):\n",
    "    tp = np.sum(np.logical_and(preds == 1, yval == 1)).astype(float)\n",
    "    fp = np.sum(np.logical_and(preds == 1, yval == 0)).astype(float)\n",
    "    fn = np.sum(np.logical_and(preds == 0, yval == 1)).astype(float)\n",
    "\n",
    "    # f1 = (2 * precision * recall) / (precision + recall) \n",
    "    # F1 Score也可以化简为如下公式，避免当tp == 0 && fp == 0时计算出现NaN的情况\n",
    "    f1 = 2*tp / (2*tp + fn + fp)\n",
    "    return f1\n",
    "\n",
    "\n",
    "import sklearn.metrics\n",
    "\"\"\"\n",
    "选择最佳epsilon值，使得F1 Score最大值\n",
    "\"\"\"\n",
    "def select_threshold(pval, yval):\n",
    "    epsilons = np.linspace(np.min(pval), np.max(pval), num=1000)\n",
    "    fs = []\n",
    "\n",
    "    for epsilon in epsilons:\n",
    "        preds = (pval < epsilon).astype('int')\n",
    "\n",
    "        fs.append(sk.metrics.f1_score(yval, preds))\n",
    "\n",
    "    # print(fs)\n",
    "    argmax_fs = np.argmax(fs)\n",
    "\n",
    "    return epsilons[argmax_fs], fs[argmax_fs]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "D:\\Anaconda3\\envs\\python36\\lib\\site-packages\\sklearn\\metrics\\classification.py:1143: UndefinedMetricWarning: F-score is ill-defined and being set to 0.0 due to no predicted samples.\n",
      "  'precision', 'predicted', average, warn_for)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(9.045285962118505e-05, 0.8750000000000001)"
      ]
     },
     "execution_count": 63,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "epsilon, f1 = select_threshold(pval, yval)\n",
    "epsilon, f1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "最后，我们可以将阈值应用于数据集，并可视化结果。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 得到异常值的索引\n",
    "p = multi_normal.pdf(X)\n",
    "outliers = np.where(p < epsilon)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAFyJJREFUeJzt3X9sXeV9x/HPN/Z17Rt7JCyBBReWtrboilAJsxgbXgV0wi5dIWVqXVR1aI2UVioSkTq00E4lXSdBx2DTpKlTOlDZyjJ3K6S0pQmIIlVBAs1AIIkCs2kpYPLDLU1x4ni+dr7745yb3jjn/vT9+fj9kq7u9Tnn+ny5OXzu4+c85znm7gIAtL4VjS4AAFAdBDoABIJAB4BAEOgAEAgCHQACQaADQCAIdAAIBIEOAIEg0AEgEO313NmaNWt8/fr19dwlALS855577hfuvrbYdnUN9PXr12tsbKyeuwSAlmdmPy9lO7pcACAQBDoABKJooJvZhWb2lJkdNLMDZnZbvHybmU2a2d74cX3tywUA5FNKH/q8pC+6+/Nm1iPpOTN7Il73D+7+97UrDwBQqqKB7u6HJB2KX0+b2UFJvbUuDABQnrL60M1svaQNkp6NF91qZi+Z2QNmtrrKtQEAylByoJtZt6TvStri7u9I+oak90m6TFEL/t4879tsZmNmNjY1NVWFkgEASUoKdDNLKQrzh9z9YUly9yPuvuDupyR9U9IVSe919+3uPuDuA2vXFh0XDwCoUNE+dDMzSfdLOuju9+UsXxf3r0vSxyXtr02Ji0xPS6Oj0vi41N8vjYxIPT112TUANLNSRrlcJekzkvaZ2d542Zck3Wxml0lySa9J+lxNKsy1Z48yQ8PKzC8oPTermY5OpW7botTuXdLgYM13DwDNrJRRLnskWcKqx6pfTgHT08oMDSs1c0KpeFF6blaaU7T8yGGpu7uuJQFAM2mdK0VHR5WZX0hclZlfiLphAGAZa51AHx+PWuQJ0nOz0sREnQsCgObSOoHe36+Zjs7EVTMdnVJfX50LAoDm0jqBPjKiVHtb4qpUe1s02gUAlrHWCfSeHqV271ImvfJ0S32mo1OZ9MpolAsnRAEsc3W9wcWSDQ4qdeSwUqOj0sSE0n19UcucMAeAFgt0KQrvTZsaXQUANJ3W6XIBABREoANAIAh0AAgEgQ4AgSDQASAQBDoABIJAB4BAEOgAEAgCHQACQaADQCAIdAAIBIEOAIEg0AEgEAQ6AASCQAeAQBDoABAIAh0AAkGgA0AgCHQACASBDgCBINABIBAEOgAEgkAHgEAQ6AAQCAIdAAJBoANAIAh0AAgEgQ4AgSDQASAQBDoABIJAB4BAFA10M7vQzJ4ys4NmdsDMbouXn2tmT5jZePy8uvblAgDyKaWFPi/pi+7+e5KulPQFM/uApK2SnnT3fklPxj8DABqkaKC7+yF3fz5+PS3poKReSTdKejDe7EFJG2tVJACguLL60M1svaQNkp6VdL67H5Ki0Jd0XrWLAwCUruRAN7NuSd+VtMXd3ynjfZvNbMzMxqampiqpEQBQgpIC3cxSisL8IXd/OF58xMzWxevXSTqa9F533+7uA+4+sHbt2mrUDABIUMooF5N0v6SD7n5fzqpHJd0Sv75F0veqXx4AoFTtJWxzlaTPSNpnZnvjZV+SdLek75jZJkmvS/pEbUoEAJSiaKC7+x5Jlmf1h6tbDgCgUlwpCgCBINABIBAEOgAEgkAHgEAQ6AAQCAIdAAJBoANAIAh0AAgEgQ4AgSDQASAQpczlAmC5mJ6WRkel8XGpv18aGZF6ehpdFUpEoAOI7NmjzNCwMvMLSs/NaqajU6nbtii1e5c0ONjo6lACAh2I7XxhUvfsfkVvHTupC1Z16fahi7VxQ2+jy6qP6WllhoaVmjmhVLwoPTcrzSlafuSw1N3d0BJRHH3ogKIwv+PhfZo8dlIuafLYSd3x8D7tfGGy0aXVx+ioMvMLiasy8wtRNwyaHoEOSLpn9ys6mTkz0E5mFnTP7lcaVFGdjY9HLfIE6blZaWKizgWhEgQ6IOmtYyfLWh6c/n7NdHQmrprp6JT6+upcECpBoAOSLljVVdby4IyMKNXelrgq1d4WjXZB0yPQAUm3D12srtSZgdaVatPtQxc3qKI66+lRavcuZdIrT7fUZzo6lUmvjEa5cEK0JTDKBctednTLycyC2sy04K7e5TbKRZIGB5U6clip0VFpYkLpvr6oZU6YtwwCHctadnRL9oTogvvplnk5YR7MkMfubmnTpkZXgQrR5YJlrRqjW5b9kEc0DVroCEKlLeRqjG4p9KXQkq10tCwCHS2hUGAv7jbJtpAlFQ3Uc7pSOnYyc9byFWZ6z9YflvTlUOhLIZiuGLQEulzQ9Ip1aVTabbLzhUm9M3t2mEtRX3qp3Sf5hjae05WiKwZ1RaCj6RUL7Eq7TbY9ekCnvPj+i305XPP+tYnLj53MLO+rT1F3dLmg6RUL7AtWdWkyYZvclvNf79ynHc++oQUvIcHz7Gtx98k171+rp16eStx3IZPHTmrD3zyuOz92Cd0vqCrzCg/wSgwMDPjY2Fjd9ocwXHX3jxNDs3dVl57eeu1ZfehZq7pS2nbDJRr7+dv69jOv16vcsmXHvEuivz0ENZhT3syec/eBYtvRQkfTu33o4rMCe/FVnO9qX3FWoB87mdEdD+/TbCZ5FsFmMXnspLaM7j1rWakndtFEGjynPH3oaHobN/TqrpsuVe+qLpmiFu1dN12qjRt6T7fOk0aqSFGfdf3+Bq0u+ttbTM6c8tmZK9Nzs0rNnFBmaFg6frzmJdBCR0vYuKE3saWadMI0JMtmtscQxHPKpxJWZeYXoikVanwVLi10tLTQA2/ZzPYYgiaYU54WOlpG0kU6+Ua4hGLZzPYYgnhO+aRQn+nojCY7qzFa6GgJ+S4uuub9a5VaYY0uryZWdaU4IdpKmmBOeVroaCr5LpX/6vcPJF6ks+PZN9SZWqHMXFj96F2pNm274ZJGl4FyZOeUXzzKpb2tbnPKMw4dTSNpPLlJLTtKpVJtZrr3kx+kdd6qjh+PxqFPTES37qvCnPKMQ0fLSRqxstzCXJJOuRPmrayBc8rTh46mEfqIlVKd05U08A0ojhY6mkboI1ZKZTnneJl+F+Uo2kI3swfM7KiZ7c9Zts3MJs1sb/y4vrZlolXsfGFSV939Y71n6w911d0/Lmuq2KQbNYeoK1X4f7tjM9FVr9wJCeUqpcvlW5KGE5b/g7tfFj8eq25ZaEVLDaCNG3r1Z7/fqzAHIUZMUmeRL63sxUTVuD0elpeige7uP5H0dh1qQYurRgA99fJU0CdCP33lRadb4ElyJx2rxu3xsLws5aTorWb2Utwls7pqFaFlLTWAdr4wGXQfeldqhR565nWtsOS/QdrMTk86JuW/7J/pAJBPpYH+DUnvk3SZpEOS7s23oZltNrMxMxubmpqqcHdoBUsJoJ0vTOr2/36x2iU1hWx+n8yckkt5b7LR03nmGIWkcwqLpw0GclU0ysXdj2Rfm9k3Jf2gwLbbJW2XoguLKtkfWkMp85ZnR21MHjupNjMtuKt3VZdO/N+8MgthHh75rt1bYTrjFnjZ+dulM2eXZJQLSlVRoJvZOnc/FP/4cUn7C22PsOQbSlcsgBZfCZptqYbczVJI0v1Ms+ccsp9ZvmmDgSRFA93Mdki6WtIaM3tT0p2SrjazyxRdyPeapM/VsEY0kcWhvPjOOoUCKPS5y6uFk56oVNFAd/ebExbfX4Na0AIKjWQp1pIkqH6jK9Wmd7WvSLzTEic9USku/UdZljKShaCKrE6ndNdNl2rbDZdw0hNVxaX/KEu+y/NLCeukk6bLSW+ek5qc9ES1EOgoSykjWfLJPWmaNMol342eQ2CSnt567VnLOemJaiLQUZalDqXLF2BJc6GHhO4m1AOBjrLVolW5+ItiRdx6DwH94qgXAh1NI/eLIpQWe75+c6AWCHRUpNbzdC9usRdqqzfDbepWSDqV83NqhemeT3AbOdQXwxZRtnrN071xQ6+e3nqtfnb3R9Wbpw+6d1VXw8N8dTql+0YuU++qLllcE2GORqCFjrIt5eKixUpt6RcaXZMdNdMIXak23fmxSxitgqZAoKNs1Zqnu9g0ArkWd8F0pVbo5PyCtozulUlqW2FaSJocpUKpFVLmVPK6NjOdcmfcOJoOXS4oW7Xm6c7X0v/q9w8kbp/tgvn0lRdpJnPq9CyGLmnhlGtlx9JvX2cm/ePIZTrvt/L/t5xy18/u/qie3notYY6mQqCjbNWapztfi/5XM5mC/fE7nn0jcfls5lTevvZSuUdfHIX+2mBMOZoVgY6ybdzQq7tuuvSMk4C5d9opVaFgLHTbunzj0xfc837ZlCtfbSYxphxNiz50VKQaJwFvH7pYW0b3Jq4r1EJuy3PRUZtZ3itZSz1xmr27UNJJWFN0T1C6WdCsCHQ0zMYNvdr26IHEOVzO6Urlfd/Nf3Chvv3M62ctv/K9q0//3qTQLeVCpez3BHcLQisi0NFQ2264RLf/14vKLBqhcmJuXjtfmEwM0L/deKl+NnVcT7/69hnLn3/913nfU+rUArl98AxFRKuhDx0NtXFDr7o7z25XZBa8YD/6a788u/skOxa+0L6yFyrd+8kPMhc5gkMLHQ13bCZ52txC/ej5+sNLvcCILhWEiEBHw1Vy04xCJ0YXK3RTawIcIaHLBQ1Xybj2QkMXc9Vr3hmgGRDoaLhKxrUXmqwrV6F5Z4DQ0OWCplBu90ept8Kr1rwzQCughY6WVGqrvlrzzgCtgBY6WlYprfql3NQaaDUEOoLG8EQsJwQ6gsfwRCwX9KEDQCAIdAAIBIEOAIEg0AEgEAQ6AASCQAeAQBDoABAIAh0AAsGFRQBqa3paGh2Vxsel/n5pZETq6Wl0VUEi0AHUzp49ygwNKzO/oPTcrGY6OpW6bYtSu3dJg4ONri44dLkAqI3paWWGhpWaOaH03KwkKT03q9TMCWWGhqXjxxtcYHgIdAC1MTqqzPxC4qrM/ELUDYOqKhroZvaAmR01s/05y841syfMbDx+Xl3bMgG0nPHx0y3zxdJzs9LERJ0LCl8pLfRvSRpetGyrpCfdvV/Sk/HPAPAb/f2a6ehMXDXT0Sn19dW5oPAVDXR3/4mktxctvlHSg/HrByVtrHJdAFrdyIhS7W2Jq1LtbdFoF1RVpX3o57v7IUmKn8+rXkkAgtDTo9TuXcqkV55uqc90dCqTXhmNcunubnCB4an5sEUz2yxpsyRddNFFtd4dgGYyOKjUkcNKjY5KExNK9/VFLXPCvCYqDfQjZrbO3Q+Z2TpJR/Nt6O7bJW2XpIGBAa9wfwBaVXe3tGlTo6tYFioN9Ecl3SLp7vj5e1WrCCgXVyICkkoIdDPbIelqSWvM7E1JdyoK8u+Y2SZJr0v6RC2LBPLiSkTgtKKB7u4351n14SrXApQn50rEVLwoPTcrzSlafuQwfbVYVrhSFK2LKxGBMxDoaF1ciQicgUBH6+JKROAMBDpaF1ciAmcg0NG6uBIROAM3uEBr40pE4DQCHa2PKxEBSXS5AEAwCHQACASBDgCBINABIBAEOgAEgkAHgEAQ6AAQCAIdAAJBoANAIAh0AAgEgQ4AgSDQASAQBDoABIJAB4BAEOgAEAgCHQACQaADQCAIdAAIBIEOAIEg0AEgEAQ6AASCQAeAQBDoABAIAh0AAkGgA0AgCHQACASBDgCBINABIBAEOgAEgkAHgEAQ6AAQCAIdAALRvpQ3m9lrkqYlLUiad/eBahQFACjfkgI9do27/6IKvwcAsAR0uQBAIJYa6C7pcTN7zsw2J21gZpvNbMzMxqamppa4OwBAPksN9Kvc/XJJH5H0BTP70OIN3H27uw+4+8DatWuXuDsAQD5LCnR3fyt+PirpEUlXVKMoAED5Kg50M1tpZj3Z15Kuk7S/WoUBAMqzlFEu50t6xMyyv+c/3H1XVaoCgGqbnpZGR6Xxcam/XxoZkXp6Gl1VVVUc6O7+U0kfrGItAFAbe/YoMzSszPyC0nOzmunoVOq2LUrt3iUNDja6uqph2CKAsE1PKzM0rNTMCaXnZiVJ6blZpWZOKDM0LB0/3uACq4dABxC20VFl5hcSV2XmF6JumEAQ6ADCNj5+umW+WHpuVpqYqHNBtUOgAwhbf79mOjoTV810dEp9fXUuqHYIdABhGxlRqr0tcVWqvS0a7RIIAh1A2Hp6lNq9S5n0ytMt9ZmOTmXSK6NRLt3dDS6weqox2yIANLfBQaWOHFZqdFSamFC6ry9qmQcU5hKBDmC56O6WNm1qdBU1RZcLAASCQAeAQBDoABAIAh0AAkGgA0AgCHQACASBDgCBINABIBDm7vXbmdmUpJ/X6NevkfSLGv3uaqC+paG+paG+pWl0fb/r7muLbVTXQK8lMxtz94FG15EP9S0N9S0N9S1Ns9eXRZcLAASCQAeAQIQU6NsbXUAR1Lc01Lc01Lc0zV6fpID60AFguQuphQ4Ay1pLBbqZvWZm+8xsr5mNJaw3M/snM5sws5fM7PI61nZxXFf28Y6ZbVm0zdVm9uucbb5Sh7oeMLOjZrY/Z9m5ZvaEmY3Hz6vzvPeWeJtxM7uljvXdY2Yvx/+Gj5jZqjzvLXg81LC+bWY2mfPveH2e9w6b2Svx8bi1jvWN5tT2mpntzfPemn5+ZnahmT1lZgfN7ICZ3RYvb4rjr0B9TXP8lc3dW+Yh6TVJawqsv17SjySZpCslPdugOtskHVY0djR3+dWSflDnWj4k6XJJ+3OW/Z2krfHrrZK+nvC+cyX9NH5eHb9eXaf6rpPUHr/+elJ9pRwPNaxvm6S/LOEYeFXSeyV1SHpR0gfqUd+i9fdK+kojPj9J6yRdHr/ukfS/kj7QLMdfgfqa5vgr99FSLfQS3Cjp3zzyjKRVZrauAXV8WNKr7l6ri6hK5u4/kfT2osU3Snowfv2gpI0Jbx2S9IS7v+3uv5L0hKThetTn7o+7+3z84zOS3l3t/ZYqz+dXiiskTbj7T919TtJ/Kvrcq6pQfWZmkj4paUe191sKdz/k7s/Hr6clHZTUqyY5/vLV10zHX7laLdBd0uNm9pyZbU5Y3yvpjZyf34yX1dunlP9/oj80sxfN7Edmdkk9i8pxvrsfkqKDWtJ5Cds0y2f5WUV/dSUpdjzU0q3xn+QP5OkyaIbP748lHXH38Tzr6/b5mdl6SRskPasmPP4W1ZerWY+/RK12T9Gr3P0tMztP0hNm9nLcQsmyhPfUdRiPmXVIukHSHQmrn1fUDXM87nfdKam/nvWVoRk+yy9Lmpf0UJ5Nih0PtfINSV9T9Hl8TVG3xmcXbdPwz0/SzSrcOq/L52dm3ZK+K2mLu78T/eFQ/G0Jy2ry+S2uL2d5sx5/ebVUC93d34qfj0p6RNGftbnelHRhzs/vlvRWfao77SOSnnf3I4tXuPs77n48fv2YpJSZralzfZJ0JNsVFT8fTdimoZ9lfBLsTyV92uMOy8VKOB5qwt2PuPuCu5+S9M08+23059cu6SZJo/m2qcfnZ2YpRWH5kLs/HC9umuMvT31NffwV0jKBbmYrzawn+1rRiYv9izZ7VNKfW+RKSb/O/mlXR3lbRWb2O3G/pszsCkWf/y/rWFvWo5KyowZukfS9hG12S7rOzFbHXQrXxctqzsyGJf2VpBvcfSbPNqUcD7WqL/e8zMfz7Pd/JPWb2Xviv9o+pehzr5c/kfSyu7+ZtLIen198rN8v6aC735ezqimOv3z1NfvxV1Cjz8qW+lA0WuDF+HFA0pfj5Z+X9Pn4tUn6Z0WjC/ZJGqhzjWlFAX1OzrLc+m6Na39R0cmWP6pDTTskHZKUUdTq2STptyU9KWk8fj433nZA0r/mvPezkibix1/Usb4JRf2ne+PHv8TbXiDpsULHQ53q+/f4+HpJUTitW1xf/PP1ikZOvFrP+uLl38oedznb1vXzkzSoqJvkpZx/y+ub5fgrUF/THH/lPrhSFAAC0TJdLgCAwgh0AAgEgQ4AgSDQASAQBDoABIJAB4BAEOgAEAgCHQAC8f8H5FaOnWKZVAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 画出异常点\n",
    "fig, ax = plt.subplots()\n",
    "ax.scatter(X[:,0], X[:,1])\n",
    "ax.scatter(X[outliers[0],0], X[outliers[0],1], s=50, color='r', marker='o')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "红点是被标记为异常值的点。 这些看起来很合理。 有一些分离（但没有被标记）的右上角也可能是一个异常值，但是相当接近。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.4 高维数据集"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = loadmat('data/ex8data2.mat')\n",
    "\n",
    "X = data['X']\n",
    "Xval = data['Xval']\n",
    "yval = data['yval']\n",
    "\n",
    "mu, cov = X.mean(axis=0), np.cov(X.T)\n",
    "multi_normal = stats.multivariate_normal(mu, cov)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "D:\\Anaconda3\\envs\\python36\\lib\\site-packages\\sklearn\\metrics\\classification.py:1143: UndefinedMetricWarning: F-score is ill-defined and being set to 0.0 due to no predicted samples.\n",
      "  'precision', 'predicted', average, warn_for)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.748247887556513e-18 0.5517241379310345\n",
      "find 122 anamolies\n"
     ]
    }
   ],
   "source": [
    "pval = multi_normal.pdf(Xval)\n",
    "\n",
    "epsilon, f1 = select_threshold(pval, yval)\n",
    "epsilon, f1\n",
    "print(epsilon, f1)\n",
    "\n",
    "p = multi_normal.pdf(X)\n",
    "print('find {} anamolies'.format(np.sum((p <= epsilon).astype('int'))))\n",
    "\n",
    "# 结果是122个异常，接近原Assignment里的117，但有一点差异，目前猜测和Python里的计算精度有关？待验证"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. 推荐系统"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "推荐引擎使用基于项目和用户的相似性度量来检查用户的历史偏好，以便为用户可能感兴趣的新“事物”提供建议。\n",
    "\n",
    "在本练习中，我们将实现一种称为协作过滤的特定推荐系统算法，并将其应用于电影评分的数据集。\n",
    "\n",
    "### 2.1 电影评分数据集\n",
    "我们首先加载并检查我们将要使用的电影评分数据集。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = loadmat('data/ex8_movies.mat')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Y是包含从1到5的等级的（数量的电影x数量的用户）数组。\n",
    "\n",
    "R是包含指示用户是否给电影评分的二进制值的“指示符”数组。 两者应该具有相同的维度。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Y = data['Y']\n",
    "R = data['R']\n",
    "Y.shape, R.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们可以通过平均排序Y来评估电影的平均评级。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Y[1,np.where(R[1,:]==1)[0]].mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们还可以通过将矩阵渲染成图像来尝试“可视化”数据。 我们不能从这里收集太多，但它确实给我们了解用户和电影的相对密度。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(12,12))\n",
    "ax.imshow(Y)\n",
    "ax.set_xlabel('Users')\n",
    "ax.set_ylabel('Movies')\n",
    "fig.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "接下来，我们将实施协同过滤的代价函数。 直觉上，“代价”是指一组电影评级预测偏离真实预测的程度。 代价方程在练习文本中给出。 它基于文本中称为X和Theta的两组参数矩阵。 这些“展开”到“参数”输入中，以便稍后可以使用SciPy的优化包。 请注意，我已经在注释中包含数组/矩阵形状（对于我们在本练习中使用的数据），以帮助说明矩阵交互如何工作。\n",
    "# cost\n",
    "<img style=\"float: left;\" src=\"../img/rcmd_cost.png\">"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cost(params, Y, R, num_features):\n",
    "    Y = np.matrix(Y)  # (1682, 943)\n",
    "    R = np.matrix(R)  # (1682, 943)\n",
    "    num_movies = Y.shape[0]\n",
    "    num_users = Y.shape[1]\n",
    "    \n",
    "    # reshape the parameter array into parameter matrices\n",
    "    X = np.matrix(np.reshape(params[:num_movies * num_features], (num_movies, num_features)))  # (1682, 10)\n",
    "    Theta = np.matrix(np.reshape(params[num_movies * num_features:], (num_users, num_features)))  # (943, 10)\n",
    "    \n",
    "    # initializations\n",
    "    J = 0\n",
    "    \n",
    "    # compute the cost\n",
    "    error = np.multiply((X * Theta.T) - Y, R)  # (1682, 943)\n",
    "    squared_error = np.power(error, 2)  # (1682, 943)\n",
    "    J = (1. / 2) * np.sum(squared_error)\n",
    "    \n",
    "    return J"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "为了测试这一点，我们提供了一组我们可以评估的预训练参数。 为了保持评估时间的少点，我们将只看一小段数据。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "params_data = loadmat('data/ex8_movieParams.mat')\n",
    "X = params_data['X']\n",
    "Theta = params_data['Theta']\n",
    "X.shape, Theta.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "users = 4\n",
    "movies = 5\n",
    "features = 3\n",
    "\n",
    "X_sub = X[:movies, :features]\n",
    "Theta_sub = Theta[:users, :features]\n",
    "Y_sub = Y[:movies, :users]\n",
    "R_sub = R[:movies, :users]\n",
    "\n",
    "params = np.concatenate((np.ravel(X_sub), np.ravel(Theta_sub)))\n",
    "\n",
    "cost(params, Y_sub, R_sub, features)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " 接下来我们需要实现梯度计算。 就像我们在练习4中使用神经网络实现一样，我们将扩展代价函数来计算梯度。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cost(params, Y, R, num_features):\n",
    "    Y = np.matrix(Y)  # (1682, 943)\n",
    "    R = np.matrix(R)  # (1682, 943)\n",
    "    num_movies = Y.shape[0]\n",
    "    num_users = Y.shape[1]\n",
    "    \n",
    "    # reshape the parameter array into parameter matrices\n",
    "    X = np.matrix(np.reshape(params[:num_movies * num_features], (num_movies, num_features)))  # (1682, 10)\n",
    "    Theta = np.matrix(np.reshape(params[num_movies * num_features:], (num_users, num_features)))  # (943, 10)\n",
    "    \n",
    "    # initializations\n",
    "    J = 0\n",
    "    X_grad = np.zeros(X.shape)  # (1682, 10)\n",
    "    Theta_grad = np.zeros(Theta.shape)  # (943, 10)\n",
    "    \n",
    "    # compute the cost\n",
    "    error = np.multiply((X * Theta.T) - Y, R)  # (1682, 943)\n",
    "    squared_error = np.power(error, 2)  # (1682, 943)\n",
    "    J = (1. / 2) * np.sum(squared_error)\n",
    "    \n",
    "    # calculate the gradients\n",
    "    X_grad = error * Theta\n",
    "    Theta_grad = error.T * X\n",
    "    \n",
    "    # unravel the gradient matrices into a single array\n",
    "    grad = np.concatenate((np.ravel(X_grad), np.ravel(Theta_grad)))\n",
    "    \n",
    "    return J, grad"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "J, grad = cost(params, Y_sub, R_sub, features)\n",
    "J, grad"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们的下一步是在代价和梯度计算中添加正则化。 我们将创建一个最终的正则化版本的功能（请注意，此版本包含一个额外的“学习率”参数，在文本中称为“lambda”）。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cost(params, Y, R, num_features, learning_rate):\n",
    "    Y = np.matrix(Y)  # (1682, 943)\n",
    "    R = np.matrix(R)  # (1682, 943)\n",
    "    num_movies = Y.shape[0]\n",
    "    num_users = Y.shape[1]\n",
    "    \n",
    "    # reshape the parameter array into parameter matrices\n",
    "    X = np.matrix(np.reshape(params[:num_movies * num_features], (num_movies, num_features)))  # (1682, 10)\n",
    "    Theta = np.matrix(np.reshape(params[num_movies * num_features:], (num_users, num_features)))  # (943, 10)\n",
    "    \n",
    "    # initializations\n",
    "    J = 0\n",
    "    X_grad = np.zeros(X.shape)  # (1682, 10)\n",
    "    Theta_grad = np.zeros(Theta.shape)  # (943, 10)\n",
    "    \n",
    "    # compute the cost\n",
    "    error = np.multiply((X * Theta.T) - Y, R)  # (1682, 943)\n",
    "    squared_error = np.power(error, 2)  # (1682, 943)\n",
    "    J = (1. / 2) * np.sum(squared_error)\n",
    "    \n",
    "    # add the cost regularization\n",
    "    J = J + ((learning_rate / 2) * np.sum(np.power(Theta, 2)))\n",
    "    J = J + ((learning_rate / 2) * np.sum(np.power(X, 2)))\n",
    "    \n",
    "    # calculate the gradients with regularization\n",
    "    X_grad = (error * Theta) + (learning_rate * X)\n",
    "    Theta_grad = (error.T * X) + (learning_rate * Theta)\n",
    "    \n",
    "    # unravel the gradient matrices into a single array\n",
    "    grad = np.concatenate((np.ravel(X_grad), np.ravel(Theta_grad)))\n",
    "    \n",
    "    return J, grad"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "J, grad = cost(params, Y_sub, R_sub, features, 1.5)\n",
    "J, grad"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "这个结果再次与练习代码的预期输出相匹配，所以看起来正则化是正常的。 在我们训练模型之前，我们有一个最后步骤， 我们的任务是创建自己的电影评分，以便我们可以使用该模型来生成个性化的推荐。 为我们提供一个连接电影索引到其标题的文件。 接着我们将文件加载到字典中。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "movie_idx = {}\n",
    "f = open('data/movie_ids.txt',encoding= 'gbk')\n",
    "for line in f:\n",
    "    tokens = line.split(' ')\n",
    "    tokens[-1] = tokens[-1][:-1]\n",
    "    movie_idx[int(tokens[0]) - 1] = ' '.join(tokens[1:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "movie_idx[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们将使用练习中提供的评分。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ratings = np.zeros((1682, 1))\n",
    "\n",
    "ratings[0] = 4\n",
    "ratings[6] = 3\n",
    "ratings[11] = 5\n",
    "ratings[53] = 4\n",
    "ratings[63] = 5\n",
    "ratings[65] = 3\n",
    "ratings[68] = 5\n",
    "ratings[97] = 2\n",
    "ratings[182] = 4\n",
    "ratings[225] = 5\n",
    "ratings[354] = 5\n",
    "\n",
    "print('Rated {0} with {1} stars.'.format(movie_idx[0], str(int(ratings[0]))))\n",
    "print('Rated {0} with {1} stars.'.format(movie_idx[6], str(int(ratings[6]))))\n",
    "print('Rated {0} with {1} stars.'.format(movie_idx[11], str(int(ratings[11]))))\n",
    "print('Rated {0} with {1} stars.'.format(movie_idx[53], str(int(ratings[53]))))\n",
    "print('Rated {0} with {1} stars.'.format(movie_idx[63], str(int(ratings[63]))))\n",
    "print('Rated {0} with {1} stars.'.format(movie_idx[65], str(int(ratings[65]))))\n",
    "print('Rated {0} with {1} stars.'.format(movie_idx[68], str(int(ratings[68]))))\n",
    "print('Rated {0} with {1} stars.'.format(movie_idx[97], str(int(ratings[97]))))\n",
    "print('Rated {0} with {1} stars.'.format(movie_idx[182], str(int(ratings[182]))))\n",
    "print('Rated {0} with {1} stars.'.format(movie_idx[225], str(int(ratings[225]))))\n",
    "print('Rated {0} with {1} stars.'.format(movie_idx[354], str(int(ratings[354]))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can add our own ratings vector to the existing data set to include in the model.\n",
    "我们可以将自己的评级向量添加到现有数据集中以包含在模型中。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "R = data['R']\n",
    "Y = data['Y']\n",
    "\n",
    "Y = np.append(Y, ratings, axis=1)\n",
    "R = np.append(R, ratings != 0, axis=1)\n",
    "\n",
    "Y.shape, R.shape, ratings.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们不只是准备训练协同过滤模型。 我们只需要定义一些变量并对评级进行规一化。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "movies = Y.shape[0]  # 1682\n",
    "users = Y.shape[1]  # 944\n",
    "features = 10\n",
    "learning_rate = 10.\n",
    "\n",
    "X = np.random.random(size=(movies, features))\n",
    "Theta = np.random.random(size=(users, features))\n",
    "params = np.concatenate((np.ravel(X), np.ravel(Theta)))\n",
    "\n",
    "X.shape, Theta.shape, params.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Ymean = np.zeros((movies, 1))\n",
    "Ynorm = np.zeros((movies, users))\n",
    "\n",
    "for i in range(movies):\n",
    "    idx = np.where(R[i,:] == 1)[0]\n",
    "    Ymean[i] = Y[i,idx].mean()\n",
    "    Ynorm[i,idx] = Y[i,idx] - Ymean[i]\n",
    "\n",
    "Ynorm.mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.optimize import minimize\n",
    "\n",
    "fmin = minimize(fun=cost, x0=params, args=(Ynorm, R, features, learning_rate), \n",
    "                method='CG', jac=True, options={'maxiter': 100})\n",
    "fmin"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = np.matrix(np.reshape(fmin.x[:movies * features], (movies, features)))\n",
    "Theta = np.matrix(np.reshape(fmin.x[movies * features:], (users, features)))\n",
    "\n",
    "X.shape, Theta.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们训练好的参数是X和Theta。 我们可以使用这些来为我们添加的用户创建一些建议。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "predictions = X * Theta.T \n",
    "my_preds = predictions[:, -1] + Ymean\n",
    "my_preds.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sorted_preds = np.sort(my_preds, axis=0)[::-1]\n",
    "sorted_preds[:10]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "这给了我们一个排名最高的评级，但我们失去了这些评级的索引。 我们实际上需要使用argsort函数来预测评分对应的电影。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "idx = np.argsort(my_preds, axis=0)[::-1]\n",
    "idx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Top 10 movie predictions:\")\n",
    "for i in range(10):\n",
    "    j = int(idx[i])\n",
    "    print('Predicted rating of {0} for movie {1}.'.format(str(float(my_preds[j])), movie_idx[j]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> 推荐的电影实际上并不符合练习文本中的内容， 还没有找到原因在哪里。"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python (Python 3.6)",
   "language": "python",
   "name": "python36"
  },
  "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.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}