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": 11,
   "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": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(307, 2)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data = loadmat('data/ex8data1.mat')\n",
    "X = data['X']\n",
    "X.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "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": 14,
   "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": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([14.11222578, 14.99771051]), array([1.83263141, 1.70974533]))"
      ]
     },
     "execution_count": 15,
     "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": 16,
   "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": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((307, 2), (307, 1))"
      ]
     },
     "execution_count": 17,
     "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": 18,
   "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": 19,
   "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(sklearn.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": 20,
   "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": 20,
     "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": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 得到异常值的索引\n",
    "p = multi_normal.pdf(X)\n",
    "outliers = np.where(p < epsilon)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "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": 23,
   "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": 24,
   "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": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = loadmat('data/ex8_movies.mat')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$Y$是包含用户对电影评分（1到5）的数组，其数组大小为电影数量 * 用户数量。\n",
    "\n",
    "$R$是包含指示用户是否给电影评分的二进制值的“指示符”数组。\n",
    "\n",
    "$Y$和$R$具有相同的维度。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((1682, 943), (1682, 943))"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Y = data['Y']\n",
    "R = data['R']\n",
    "Y.shape, R.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "m, u = Y.shape\n",
    "# m: 电影数量\n",
    "# u: 用户数量"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "接下来，我们查看第一步电影（Toy Story）的平均打分："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.8783185840707963"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Y[0, np.where(R[0,:]==1)[0]].mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "还可以通过将用户评分矩阵绘制成成图像来“可视化”数据。 \n",
    "\n",
    "我们不能从这个图里获得很多信息，但它确实给我们了解用户和电影的相对密度。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAARIAAAGoCAYAAACKZmowAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJztvX+0XUd15/ndkYyInDi2ARnr2W5bIKAVG194fsZWx700ckDA0LGnO4mwm2A39vLqAIkTJxOZwT2ku5XuaDptAtMxibGNCUayCAHjYZG8EAsFegToSeHKMnJA6gdjSzIWxPxIMK34vez545562rduVZ06v889d3/W0tK750dVnTp1du3atWsXMTMURVGK8GNNF0BRlPFHBYmiKIVRQaIoSmFUkCiKUhgVJIqiFEYFiaIohRkbQUJEryeirxHRESK6venyKIpyChoHPxIiWgbg6wBeC+AogDkA1zHzoUYLpigKgPHRSC4HcISZ55n5HwA8COCahsukKErC8qYLEMkUgCfF76MAXiMvIKJbANwCAMuwbHolzggm+NyLT8dp3/phycWslpe98ll8/dGV3vPymXzPd/o6xg8Pkfe+Mjh50Y9jxTd+NJLncy8+HQBw2rd+mPosMfBPrgT93bPOPH3Y+abdJ/OIvWeckO/e1Tb+Dt/9DjO/KC2dcREk5Dg2NCZj5rsB3A0AZ9DZ/L2nXoBNq3uYPd4HANxx4hIAwFxvGQBgZnYRc71lS+c3re5hYeM0lu/aj5n+IrauOohNq3uDa5Pfa3ffiDXX95fynD3eX7rGd8z8lmkvlSH52z4GAFtXHcQdJy4ZOj/zkUXstcpl7pPPNNNfxM5D00NlXbpmx+Ca+e29pfMzs4tDdQMAx7asx9S2PVjYOI0r79y7VFf287o4tnlw70x/EV+87fJT93+ljztOXIKdh6axd8P9Q2mZ8tjPZM5tXrd/pAyzX+9j7e5fxprr+5jd08cdJ16zdK9JTz4nAMx8ZBFnyTp76FR+5h3JfCWmfHZd2c8ADN7NQzuuwtS2PafyEs9m16XruX2Ye027BIC1u29cKqO5xrQf+x67jkDDf8t2+Z8v/eT/F1OmcRnaHAVwvvh9HoDjsTdvWt3DXG/Z0IsyL2DT6t7SCzWNaK63bOglG6FiGpJpTOblAacEgP2h2WmbfBc2Ti+Vx7w0c94cN9fK6wAsCT+X8DHHN6879UHIa8w98vzWVQex89D00PXmA7jyzr1L6drMHj/1Ybnunestw9MzKwAMGqp5D3bZgFMfwNZVB0c+4MMb7l8qt12/a67vL9WHK73DG+4fOv7QjquGfss2IYWITOfwhvtxeMP9Sx+Yrz6kwJrrLcNjt941dN53nyyHqz5D99xx4pKltinzt+vJdKSmDLPH+0OC0VWWWMEGjI+xdTkGxtarARzDwNh6PTN/1XX9GXQ2v4aurrGE7cNoFC7snrcujm1Zj2uv+0KmBmqTpefuKraWFUNavfnO/yV/fD8zX5aW/lhoJMy8AOCdAGYBPA7gYz4hogwEha1SS7IKEVtDyHvP1LY9+OJtl+dOE8jWSwJwaiux1x/bst5ZzixaQyyufExZ7PxsIXJsy/qRe+2yp9WbrallZSwECQAw82eY+WXM/BJm/p206+0XYzco81uq5/YLMy9IpiWvMWnI86F8YjDpz2/vRTXY+e29kWddvmt/pkZ0bMv6ocYoh2+zx/tLQxFXffgw9gl5/bEt6/HEzYtLarUP+dyusslzMm37uvntPdxx4pLg+5HpzfQXh4YEJy99Nuqj9T1DSAilCQfg1PDEfn+yjRjsdjC1bc+QvUbeK5/B3HPtdV+IeCo/YzG0yYoObepBGmLros6hTaxxeRzIUm/y2tihjQoSRVG8dMpGkoesQ4oQecfzVRA7PrfL7Bu2ZcnPTnNh47R9ufNaO+8i9ekaWuYltg5CZcibft6886ZXhV1H0llBYlTSIqqpr/KbFCw+/wWbNKu+b0bHlZ9JP8tMwYoDfmcz37RjbHkMRYcdaXXge8+xQ4STlz6b61we0p6l6uFgZwVJWm/l6k1tqe7yuQCyfVA2aUbgLIR6ITsfaZXPWgbz/PbH/8gD93rvkcY720dl9nh/SEDZzB7vB8+5qKrHPbZl/UjaIU1M4nQITNIq0oZcuNqCbYyNJU9ddlaQGAccH1feuXfkmE+q+6R5Hs3EbkAxvarP0m96tZCDl0E+W5p3ro1xYrIdqq5+y03ee2Sd+Z7Z9zFtWt1zPhOQ/l7LZmrbnpHnLuKDk0UziG1fts+QEbYhrTBEHu2ls4JEevC5ejHZOFzTuMCpnse+f2Hj9JLbtrkvVrPw5QUMGoQpr5yWlp6usrc2+ZseX5bBnvLz2YzW7r5xpI5sO4Tv43UJY9fz2lOh9vNJ5DPbfhTz23vYeWh6qfeVZd666qDTniOf3b5Hlsmk6dI+7jhxiXMK2p6Gtet2fntvRFOwy27nJ3+7BK25T6Zz8tJnh7QkI6ilZ63x8DXTv+b60DvIgs7aNEweL0WgmalJu6yx079lTtnmra+25VXncwD528vEz9rkwadKpknoIuPzmMbk6i3NAizAXe6sw64Ye0uWhi/LattIsto/JL4y2GmWYR/wDQ186fm0Dxd1CBFZTqlVVmFP6qwgCalutlpt/r7gHrcHoexNXQ3U/lDyIr0TZT724i9ZNmBYbXUhPVWzMtNf9NornrjZr2nI45vX7Xfmbw8J5apmn8rtew5fHsDowj2TlguXq7gtOHz1YdvYXO8tRMyUsG8o5Jpul+U0bdSnHbq+hyzo0EZRxpzQAs2iTPzQRq4Vkb+B0XUsvr/lvT5jpCuvLOWykcZWX5l860XsZ3SlIf+X9RAy1PpmjXwGU7ts5l9ozYudrm99jS8N33Fbi/Fd5yqPKYcxUNptyTacpzkByvqW56T27GqLrnKGjLSuYzKtMpw0bVQjsfCpfq6l9+baqg2fskx2wJpQmQ0yuJN9nas3s9Mzv+3jCxun8cTNi97ZBV8AJoPP4Biqz9Cz2vUkA1stBUNyBPgJ5S2DKqWVO8tzxFwXY6SePd73BtyS9+c17k68RgKk9/oGnz1F9jhGiMhexbyksnwbfFqBDHrjEgYG234jPyLXfbPH+0tR0CQ+m5B9v0+ImDwN9tQsMBBgPruGvWLXEPqwjm1ZP2Srku/E9XfMLNKxLeux5vr+yLU+G0kMLo0kT/sx95ugRlKDM3Uv7T3SqdDWsO2V03norEZiQi2WMfVY5vRlbFrz23tYcWDlSFwRGQLS7mXsUH8xPW6WMrk0Ep9zVkxZyqSK6VRf+MdYDaJK20UMZdS7aiQ45YrtOu4aX/tsD9I5yqXlZJleTItSZdLfvG7/yHqMh3ZcNTTtK+Nzyt92+e20XWXy2TpkXdkzGk/PrPCmGSpLyAVe2glcznshm4zPyc9n93Ld6zq+89B0lAZhe8BObdsTZY8IOT7a+Gx+LmLafll0XiPJiqtna6JncfUmeW0HrnF+GWWKDdnoKnea/aAtsUCKaDox9qcY0u5xlbFIOE2Z3sRrJLLXyDL+G4o4nkj9LCs1Y6W9r0wmT5c/Q+jDCjU01zg/tjzSzdsOUBzbUF3lThMSZQqRIuP/vOtVAPfarTzC3B4u2VqTS9AVWQ+UR3B2ViMpw4+kjrH9OODrlWN7vbrq0VWeWI3CXFe19pkl/SrqLauGpRHSSnRIm9/ew2FrH5ZxpkzDZFUfXt1rURQ3KkjUs1VRciHtUxNvI0kjzfof8mR1HXf5SfjSLoJrZinNqzEtnZDHqm82IytpnpWumZW0tCSx9RryPg0t/JPkWRMDpHvShu71lSt0rWsGzD7v8so1M4NZZndUI7EYZ7uIDgfy46o7X33m9bvxnSsyPKy6vU68RlLmfLkvcJHxhyiiZdjp2oFrQvdI/xcZaAkY9huR613kfSaftLry9U6he6QPhL12JpSfySu0jsZ1r6++fPVrYwSGSyuxtb7YD/ehHVd5NSCZxtS2PdF74NjlD4Ws9GlF6kcSSYxGUobfRcx6krJxrXtx5R3STuz1J+YjNL2iXL9hesuy/Dp823bK8ub1l8lTxhitoQpNz9d2zLt47Na7vM+SZ+bHtUYL8M9WmeMTb2y1HdLGeciSRsyzFXl+34cU26Bd95f9PsqaPaq6ndgdQUx+aYJspr+InYemo4VdWp7y/MQLknGbtVH7htJGJt5Gkhdj8yg6uxIat7rOuQx9RWZKspTfda1vTUoM9spSibHJ2Hmm2TKybL0Rwn4e+dte2+IL11h0BgtIj8cSk0eW+DdZ7slD5zUSV0/vGy/6aMtamzStRY67yyxzUXW/7OGCKz3bzmOw7UFpFFm17LMzuK6J8QquQkt12bxcMVh0rY1gpu+Ol7FpdS9TA0nbqb1or+3CFQc2rVHJ9UW+9UExPjH28TwxaUOzHFm1Pft617s7Z+6kU3Be/ZabouJ92CuqgcE7MFuPhMpjsPN3rdMxbemJm9PbTKwQybKWyOQvBYctZPMIr84KEhNCQE53Zf2IYsnT2/p6SBlGwCb0Adq9pz0VLHEtTw8JizyBfOwIa3KoltZQ7edMm30CBiENXIR2A/SV17DiwEo8PbMi+sOyy+0K/mzqWaZ5bMt651A2tl2evPTZaOFclSG5s4IEOBVRzHy0duM2mJW2Pi9R+3pg2LckiyBKuzakgh/ecL/TxnJsy3pv7AngVISy+e0DTcz0hjItV6T1+e09fPG2y4fCF0oWNk57NStZxq2rDo6sQA75M5jnNPj8KGRdXXvdF4a8OH0akSuujC/9qW178Nitd420C7uuzLP6FjbK68y9dj6HN9zvjQsTajOhjscXh8UltMy7zOtn0mlBMtdbNrL03SCD0Bh1z9cr7zw0PXRuYeNgqs00ZDugTRHkx29zx4lLsHXVwZGGfO11XxgxEMtIaub3BfcMNDTT4Od6y5b+vvh9bx/Jb831fTxx8ylntDtOXDJUj488cC/W7r7R+Rwu24D9jEC5xkD5jlw7KUpk3cQM3aQwsodKRvC5PnhpB3HZ6oBBHbg6O/MMrvLJvFYcWOlsL64ym+vsjmf5rv1LITl97zREZ42t7zywIVfA3roo0wAZ46dhDGo+42MeZ64iYQR8z1/kPTXpKxNTHpfTnfk/1p3edSyLK/9Mf3Gkkwndp8bWBF8Ywiz4FonZrudZyhNqNCZ4jS/Unm88bGskEtNwNq3ujfQ4MVORWVVeee3WVQdH3Nx9z7/iwMqgcTPm2V1DPxfyvfk+RPs6+dsXmtNgP6N8Nvk+gLBW5KorecwXfMkVdHvrqoOpu+7lEcadFyTSDVkS4yMg/3Y1YDN0yjK0sdViO++53jJMbduDrasOOoP03HHikpFGv3nd/pFe6eL3vX2owct8zPXmfJomsvPQNDat7o0Mfx554F6vcJHCyrjby/x8gsnYC6QNSo7/Q4ZfaciUq5alRiHzDNmVTHonL33W+37tenOFxrSfzZTfPifrxybvZIDLxWHt7huHnqcsrbh2QUJE5xPR54jocSL6KhHdmhw/m4g+S0SHk//PSo4TEb2fiI4Q0aNE9OqYfLauOjjUWO0Kk0YvGarf1zO4jGSGLNsJ2AY1n8+CFATmf9s2Yo5vXXVw5AOb2rZnKb1Nqwdrh0zZjTu1NKLaywkMm1afMhBee90Xhp4/9NxSqBkbgm2vkkMtc9yU5fCG+0fSD2kywLCQkXYgFz4BKPPcvG4/Nq/bj7W7bxwSxK424uposrQL17O5Fvi57jt56bNOYWMHrZ493vcG6JJG6DyCqwmNZAHAbzDzPwVwBYB3ENE6ALcDeISZ1wJ4JPkNAG8AsDb5dwuAD8Rm5FPhFjZOD6mW0gAlPwDzAqUE37S6t7Qni32dq3G6Vn+GVFUj/FwL2+wGYGabjF+MrWWZZ5bq9Ozx/ogqLCPTA8Mfyvz24aGQbcBN+1iObVk/co1dn4a53rKl32t33+jsTUONPEvv6hJiwPCzSw1H1plLK1pzvVvDsoenPu3G1LFrtsnlBSzfqU9LkwtPbQc0u6z2UCsrtQsSZn6Kmf86+fvvADwOYArANQA+nFz2YQDXJn9fA+CPecCXAJxJROfG5PXF2y4HMCrtbd+CtI9h7e4bhyp4atseZ6/kegl5bDObVveG0peNQWo0xunMNAqpsWxetz+6jKEyS03MtnVc/L63O4NUA6cav6krewYhrc5dgjoUxDo05PFpH7YQM3nINM1xOTzylcF+prnespHhqe+5TfldadttyGyKZbBnFV3ltAWFa28j132xNGojIaILAbwKwJcBnMPMTwEDYQNgVXLZFIAnxW1Hk2N2WrcQ0T4i2vccTgI4JTDs3kf2wMaKbbCNggCGhgXyGp8xLkTIqCbzdvl1+HaJM8el5rDz0PTQx2WecfZ4f8SNO1T+O05cgs3r9i/ZSGQDfuzWu5xOV7L8JiSDzGPrqoNLhkCXr8P89mFBKuvZp3rvPDTttHuZ5RDyeULYmplryOIz+Lo8eH12NpfQ8RlvfcZ+eS7Nj0SmdceJS4KG83EZ2gAAiOgnAPwpgF9j5h+ELnUcG5mzZua7mfkyZr7sNKzAptW9oBuweelSnQaGpb/0E7GvyTsLFOo5pW+Kb9tK2QBsO4BdJnP88Ib7lwylLnforasOensh6Yvz2K13DeUf8jcw+Zhe3TdMc/Waa67vDwlS+a7sewyb1+0f0SZMmmlGVYl9rW9bURvfQk+Zniyfa4izafVwGqHntd99mhZhBLN5FyFDcR6/qEYECRGdhoEQ+SgzfyI5/LQZsiT/n0iOHwVwvrj9PADH0/Kwp9mK4DNEhhyGfPgMagZpHJaGWHMsy/PIOrANpTK/tDTNmNoWzisOrEzt4V3emi4hIsucx6fF7tFDgjEWnxB0dSKuzsVlN0nLP62Dcr1Dn6Oba4o4xpU+j0NaE7M2BOBeAI8z853i1MMAbkj+vgHAp8TxtyazN1cA+L4ZAoWQPb+r8mJ8C1z3+3rQWNIWoMl9fe1r7EaUtpjMNhzn2agKOGUjstM/Z+5k0B/GLr/EpxFIrcyF79z89p7znE+bCEUfM8z1ljnbiS9N1/vx5WNPpccOJ2Ke0XWNGc762mqaT00aTWgk/wzALwHYSET95N8bAfwugNcS0WEAr01+A8BnAMwDOALggwBGfbkd2NLYrmy5otfVqOX1xkbii6Ehx+O+tRLSMc4e30pMj+waApkpYXOfa1dAYFRIxjjN2fe7xue2D4TPq9Xl0SntBaHnlr/Typ42Vemazk/rkaX9yNSjHLb4bCZSaPjKZBbnAafa3+zxvnetlK+dyLZontFuw7INmPU15m9Zp3K2Sh6b397LtKq4sy7yob1/FzZO44mbF4eGDiGV03b2Mr+LuHOH1He7fPK47+O1F47Z15pn9D2r7/js8VPxW+3nDbmP234jMu1jW9YHY5K6/Fp8ZTPX5dnrVt7vysP3fmOPxyxdyIsse0w7zNJmNR6JwDXOMxJ2+a79Qz2+7TeRxgX3uI2bWQip7488cK/T7Tn0odhl8V3rm671NW7XYj5DaE9kmZ6tFaStXfH5y4Sue3pmhXNJQdqMlF1e+f5NUGQb8/6LENPOTN6uZ5DvJaYdmmuyXJuFzgoSV2XIBmwa97Et64cahrzP1TCBUx+prU7OHu8v3WO//NgNs4yznOtjW9g4PTIEkWnImQ6fLebkpc86n8sVEtFcb4Suy5iYhsxrYeN0rg29fQJLPuPJS591BguyvWlt25PtqGgLcJN3TLlj6kPOpKVh2oDxhzLMbz/VPky57LZhrnPhm2Ey5HlHnRUkBt8Y27z0qW17UlVi36Iou+fctLq3lJbdw/t6YZ8m4XqZ9odiG9DsGZVY5rf3RvxBTNoX3DM6a5GF5btOTcsu37V/Kd0s7uM+pMaQpiX4psiltrN81/4hoSVVfPn+rrxzbzAvXwckz2WZGbHbp3wGn6PcsS3rvYLNTJX71pvlWe3cWRuJjCJf1rg0Jq08U5dFqTOmbNFwDGW+i7LJY2dpM1nela8NTbyNRJK34dpDh7S0ZvqLuaf1bHzjfVd6rgYQMwaX18aGCIixD8i0pGZ1bMt65+5zaYTqMJRW2lS+/Xv5rlEPZpdmGFv+0BAh7ZlCefhWq9tkEfhT2/bkbqvAhGgkSnuoQmMrqiVVjeuZ69Bcs2h/vjqcaI3kuRefPrLOwF7xKAP2mJ5ezrWHwgGa36bHscMIuvK28zDRyc1xaaGX8TMl8pi9utf2OZB+Dy5/C9eqZJO+rflI/wJTTrsuffea+0L166tLWXbbV0fWnbSV2Nf53p28xtZCTBq2NuXSUEOY640zn3wP0kbk0wZd71Iel23LpBHy1Qlh+y3F+B5JJkIjifFDyHIujay+Gmnnysg7RBP79pRJ1vLLJQeh+qpSayg77arKqlt26tBGUQoz0UObkxf9eOo1Lv8NH1nW5bjSDqXlSyPrcR8xBkdg+PmNi7wvUE+WMsh0zbAN8KvzIRd03zHXMCyNmGGKb5gTm4fL58cegsSUM7SkIvSe7HRiy5qHTgqSFd/4EYBw5dhjQtudXDYwX8yNh3ZcNTKWNFtVyNB189t7Sy9WphUaz7qWcssATdJWIu0CErmGw9g3jH+BLPf89t6Qj8rUtj14aMdVS3FbZPmuvHPv0FDAtjmYY+a4THeud2pPHZ8a7ltZa96XvRTApCuHKzauTsIMhWTcFMlMf3HJ38K2B8Wu9pZlMdt0mtg4treuqx7NfXbIAFsYu+py+a79Q+tqTLvz5WOHJs3KxA5t5NoDGQkr7Z4VB1aOjMdj/A+y2i6KzkTksRv4Fo65njkLWZ89iz9H1rUjIew689VhkXeTpS7Ktl2lldt1fqKHNkD6RkumhzM7nEl8fhVrru+P7AM8e7y/JP1DpPmfuPIKnXf1/hKfFhVTPjm7YJ7ZXukcKpsrbTkrlXaPy3PU97x2zBYXIS9f+Z6lV+tMf3Foha6hiBBZ2Dg9olGGyuUTIln8cEy9zB4fjZdia1qyLrMOdTorSOzd1YGBhLfH7LZaaBqK635gVCBc/Zablo67XrA9tnaNVdOEjCs4kJk+9C2Td8U18Y2T7fG7jBm7sHF6KULcwsbpkTUi9joQ+wMx9W0inbsCBZmhn7nWVR+hYEgmXZmGnE6e2rbHO90s60kKnLnesqU68EXQk+n43M1lGaWWJR0XXR+tmTI209yy3UqnQDmENu3bHEuzi8gYM9Jl3/WO0uikIHnuxac7F8/JHkaGFzQfuPFJsOfvbaQviMzD3GuOz2/vLeU50x/EhjW/pW3D3h9W/i9D6Un/CdPLr7m+v2SrkesrjJHONLqZ/uKQ96I5v7BxemkhnxEWwKnG+vTMiqVG//TMiiUBZsq3fNf+oQZrx8M12poJ6WjulXaoC+4ZBEl2aSKyLg3m2U065oNfvms/lu/aj4WN00v1bP43cUrlYkTpy7Pm+sG2DvKDtPO1/XZMnZp3L+03MqaLS5i7NEa7XkzUtbnesqU9hGb6i7jyzr1LZZDBnKe27cGVd+5dirZvtiSR18nymzzs5zP1rn4klo3ENccux6ry77QYFb40fHnJeB4So/lk9TkxH06aep1lLD6/vefd70SW1QimNPuFVP/lOD80LEg757JjmWeMDRdphJx8165NpGT6D+24aiRYdqxNzZdmmv+K7x7Xs+w8NF3IhmXej6v+Z4/3sezcI5NrI3nuxacv/Z3mdGa/JNmj+oYrBp8KLiW53DrA9jC1Z2GkZmIPWaRqHhrrxvYidnq+uLTAqQ2u1lzfH9rKQ46/j21Zv/R8snznzJ0cOu6yDfgiwklNT2ozdvltLciFHBqZ8hjtylcm4NQaFKMp2Fti2Pn68rfbmpklicFuh/JZNq/b7xQitserLy9zr6v+szi4dVIjWbFmin+0Z2VqRcRqA2m/bWzpnkXLiS1r3rKn5QHELXJ09WA+rcI8v0+oy/t85XXNYMRoJK76AE4Nj05e+qz3OWz13/V8dvppmqa5RqYbwpVOSKvJqokabcb37tWztaBnaxa12XdvXfg+4KwzDFWV2zeMNBQZ8uSZQYl9Tjv9mA83tkyxZUhrf7HD77zvduKnf9NIUytNpae9xNC9Wcg63SaRQydDTIO260DODkm1WF4bqrfQM4S0nVA5Q+dCewQZXOX1OdS50pfnXWW3I7DFCrbYvWNC7U8OjeUxV/ny7FWThU5qJJdd+nw+6+C/WPrd9mXmTROrddWtaTVNGQvhfM51dYQRKIOJ1ki+9o0XDv2O6blsQj1vjEOQHZ/Ut5YnT3xMO5+s6dr3xDZoM80skdOVEtfUdlo5bOyQAr57szhoudL2lccX+DqL9uib4Yqtc1deWdzYY9tXHic0SSc1knN/+iy++PGfdZ4rsyfIOh5uskcvM+8qXMfLuL/KPLsWhjGWidZInvnRSq/ULlOdjG30tjGuSnzpZxEiaT1TKJB13kVf5v48ZO1JQ1PpProuRGybWFY6KUh+7Ec/lltgxCwXNypwmmpuk0Wg2Ns4uDDlk+ezBPuNKWca9nOE/FEkedVo17L+rJqWK+xh3jLlHZrG5lV06Bub31xvWaHI/p0UJP/44/+Y+165KM/XQE3vVKSXSvvgZdrLd+13umqb8uUpR9GA0bJh+hq7Wf/iYuuqg9EfU5pLd+z43lfO0DqeNOxFnJJQffo2/bKR78lln4ol7dmKDn07KUgu+cm/Tb3GfsmmVyrDjhBqQHnVdzv+hi+/2AVbvjzSGqcUYsZbVbrA29fZK6ONYNm0urf0HLYWaKvYZs2JxMRMMWUJhUuU98iy2X8bfHv+uq61p3/lNaF3Hdql0EdM28wahMsQG2PFRyeNrT9+7vn8M09f0XQxxppJNS4C6c9eRt2My1T6RBtbT/vWD4MS2Ce17Z7QtXqzLNIMW1nUVte1sjeNSc8uj/lQfD2x655Q+ewhTqhO89gq0uxIscz0F4eeXf5vcG0N6sIXTgCIG9pUbZgvk05qJBr8uR7GPfq8YVycw5pgojUSIF+v1qYeIDR2D12TJb1YfBrcOXMnc6dZpDxZiNGaYmeaDEXaSZZ729Qe0+isIJHjz9hpWmkcc3ml5p3rwLPfAAAgAElEQVSKyzM/7wphYPf+WcbYVYzHm7ahVPGhxdRTEae5LPeGFiu2DR3aKLkpYnTsyrCoDppcK6ZDG6GitlGCp5HV2Q3Ip/k0RdlCpCzHrRh8wx+7nZVVpjQhUvaz50mvs4JEqqi+gL2u3zGLuerA1dOnCcRYg2Hdz1PHR16WYIoRxr7hj8vXpQ7KzidPeo0JEiJaRkRfIaJPJ78vIqIvE9FhItpJRM9Ljq9Ifh9Jzl8Yk37oo7MrSv52VWKa+l7kQ8lidIyJ0xpzLutwpKhh9LFb78pUR3nzy3Kfr334dg8omzo1KKD4avY0mtRIbgXwuPi9DcB7mXktgO8CuCk5fhOA7zLzSwG8N7kulZCXps9fJOvQwNyXx0vREAo8HItpCKHgNWmBbUIN2+f1KOO32sjyu9ZwhJ4v5LmbB9f9PqFse6ZW1UmUqUXElDGkrZZhf2nE2EpE5wH4MIDfAXAbgH8B4NsAXszMC0R0JYDfZuZNRDSb/P1FIloO4FsAXsSBgp9BZ/M7D2wYcsEu8uKaMHa1ybPU54UZW8ZxNayOa7kNvmj+WWi7sfX3AfwWALO67gUAvsfMC8nvowCmkr+nADwJAMn57yfXR2NXYBU9XNk8PbOiUL4xq4frwqWx+QIiNY3UInwfXmx9tsXeBsC5HqpMahckRPQmACeYWXZl5LiUI87JdG8hon1EtO85nAwuQjLbKwCjK0ddC+DMDmx2XE57oZbZMFxiNma21Vzpem5HApvf3lvamc5WtaXbttxNzvyz0wNGN0h34auDmf6idzc/nzYilxbM9BdHNgA/tmU9tq46GOUnMdN3705ozsk0fUNT31IB+Q6MG798VlmnEpeQl/mYduASoObemKGrK7yBXNQo24e9JEJulibLtnnd/qX2aC9/cC08jKX2oQ0R/WcAvwRgAcDzAZwB4JMANqHEoY3xI3G5P8dEBw9db3Cp9mVFQwttv5B2rWn8rh7VF8U9doOsrBHS7TryDYfy1pcpU+wwy/fuTZ2dM3dyKZ2ikeDqHBqVsfG66/hYbEdBRBsA/CYzv4mI/gTAnzLzg0T0hwAeZea7iOgdAC5h5n9LRG8G8C+Z+RdD6U6aQ5rrI6xidamdZpaPN3Y3uDYG6i5zLU5TYTd99ZrWdtpuI3GxBcBtRHQEAxuImRK4F8ALkuO3Abg9JjF7iCHVtdDsTB6nrirGnVKlTlPZdx6aHlFtdx6aHlKBi5TRt5WEa69eF4c33D8SAMinPqcJkdBzuIYhvhk6H/b5Y1vWOzXaGOT+v4bQ7JkvcFFaezX7IIeQm3zJ4ZHc6dGQR8A1KkiYeTczvyn5e56ZL2fmlzLzLzDzyeT4/0x+vzQ5Px+T9uEN9wPwByxyCZZQGIHQtphZo9THCCv54fkWlZl8N6/bj83r9i89s83WVQdHztljY1+wH9lA7caa5aOwz4XuNXYlV76+Z3QJkZn+4tKG2naZfe9a2kh8dg5fXdr1c/LSZ0fanZkKd32sW1cdHHoOn53P3q4UiJ/CPbzhfmxddbBQWEUXnV1r872nXlCqOgpUs/DNzifrpuWx9wH+8qep7r5tKBc2TuPKO/dmDuNXx7DLN4QwzxIzXJnfXmzDcF+aaVt62sRcG7OdaCzjPrQpFVvi2qt5XbMzgFttDc0AGdWyDOROd65y+By15H2yXD4BYvfG9sbnNj6N64mb8zXUtA208wwvbQ3HVy6jTfiEiHyXWWw1rjK72oVJ0x5SuKaFzTtxaW+hGRqTnq/TcVHUe7mzGok0trbJuatLaL2WR1vrcuI1EtkbpIXGqyvIThayaDl17ZWTxZmqjc5mQPq7LuIw5jLU+rDrx+ebEkPZbSWPht1ZjeS82/9dZ9yb24LtF9HWXnTSqLKtTLxGEtpvJBSI2NUj1b1SMw9Zp0VjkVPH4yaYy3pvTWtXaVqU8bxuks4KEpvQMnopzbP0sFUMiUJpuvwSDLa3pvyIYgRAaJMre0mAuTa2ruR95qP0+UxkxRfpHXA/t/2cMULCNRuSpWwusvizhFwXZvqLmNq2JzWer29KP/aeNDo7tKnDs7Vu70SlXqpycS/SbuqIeK/Tv4JY6Roz1eiT5lUIEVd58vTWMQ53Mcihnu0hm8UwGTvMmD3eL3WFbNa6k89X1VAutt24NIcYIeK6L8vwbOw8W6vErgxf5YRejMsPQ67ANZTZ8E15XJtkp+UjPxq5AZPxEwk9q+9Dv/LOvUPewVLND82G2ep3bPCnTat7uQy4vroxywdc5UrD9/EViY2b5YMuY3vXomnF0llBktfQFtI4zDnbZmAavlzDIHGFF4gtg8TeTNxgbA9SeMjedH77qX12fWk/dutdzuMP7bhqRACZZww1TpfglvUQ+zGG3qM85xM+m9ftz/QRSfd3X4S0tbtvdN4bM/2bZTmFz54UEy7BJ1jt5R6SIja/zgoS1xL8GHx729jnALfW4+r1N612uy+78nDlk3bu5KXPDm3mbdI0ech7fGn7tBW5CZbcADwra67v4+q33JR5hiE0vIjZoMu1KE2SJuBd+fvq0G4PrnvlYkqDr2xZ35U9aeAa3m5etx9rd9/oXNRaZKjeWWNrmWtt2k6TPicxeZdtIKzLyG2MrTHl913jqx9zfZXGU1c9xeSnxlZBTGUVpWmPWHuIYfdsMcOpkKpbFptW96KHdjHX1TVTZjSKmA/dt5rWN5QxaZr7XG2paPty1VPMs6QFCnfRWUGSRhmNMU1trgqTn2kUpsHZjSRGSwlFRoshJo/Z4/3oIU3bvHkNae/Y157S2pk577quaBvNa2TOE2Kgs4Ika0CbWFzp1D2EsvOL6UF8U9hZy57Hy/OOE5eUNpXalJdp2SEpfL/T7nEJNF+dSKFshFJMfqqRCGypWpY67JuNaBIzi5DHoAxkK799b8yHnadh+sgaRKptxLol+K5xCTTfHk6uYzH55RGanRUkvkhYeQh9LGUby4qU1w7s4yu3fTzGl8bVKI9tWe8dsthxVcraaU/GhkmLI1MUVxT3rJRdLt978EVcM2RpV3met7OCxHwcro8krVGbj9Bc5wvvZ9LPUvH2EMPeEsBoUvZx+155TMbjNNjOY6G9kO14rxIZZczOf2rbnqEhi4wxavJbcWAlNq3uLV1nh3X05el6TmD4XUgbUWhnRUNandrtwmh6c71lUQLBFa7ygnvid1L0+Y1IfO/BhfR38XUW8p3JPLI6WXZ2+lfutDeuFJnWjZ3m27rqYJRG5YoT2qYwAlVMgecJi9gmZLmzaM7y2omf/h3HF2/3erEfRtY1Gbbrf5raa9K300wLGFUFMaufY4hR9Y09JlYjyUIdRmPp6Vx2sGebzmokVa/+rWMVZhN5ZaFNGkkVyI3Byg4CPS5MvEYS6nF847/Q2NQev25a3ast4FGMV2Xa8byGSde4XxJrOyiyWC7WsS4LMfYno+Wsub6fS4i40gwZcGPeo4us7zMtvTwG/84KEvPxuT4muQ4hNoyAHf0bGBi5XAa8OnDZLCT2al/jjg2kT+HaQXWM8dmOpr98V9yCuLnesiE1O62+7DRjAgtJo6bsKOxn8+1zFDLwpm1A5Xse36JPwD0l7usw0jqSmPqROw74XPlj83OhQ5uSaeswJETVa3XKSj+v0TM2f9ceyuMWXjJEqP58dTTxQ5syjVlZwvPFCpFQr2zvuxNL3mcOfWShNGOnCH0bVmclr40iVojZQsMnRHxD2ixDsCY8dEP1V1TQd1aQlIndoEylF7GRuPwBDLFDBhvp9ZmloS5snPaWxecHAcTFbHXFEs2joaQJrbzvIiaWqY1PwMR4/Uq7S15i/E1iKCtuLtDhoc0khRHISllDjSKzNl2f8ekKEz+08UWxMpRpIK1LTc2yEDE081HHClt7S0mbuoTITH8xk7Yi67jIECZ0fxNUPSHQWY2kC56tbSdWqxhHA3TTtMWbduI1kjJXnPqoSspXoeE0ucmTS4g0tWI6j03EJlbTKGpDa4o85e6sIElzCc6yCMt2ynJFl89K6EPavG6/U0j5VuK6jtvPZwyxWT8euaAv74fnckjL4rIt841d0WywPwr5ztKM0z6jZmxEfNd1eQylruDPMXUin13me2zL+mD7C+1S6aOzQ5sX7fiN1kbbiiGPD4M0otbhAxE7ZBknf4y6fGpihi5NxuI1TPzQJgtpWxjm2eaxKLG9nkTGBinzwy36vHUJkTL3F6oKl4d02rXjwMQKkiyzGD4/krZgBJ3vgy1qy/EFl776LTcVShfIv37Hhc/wm0XA1PVu615WUXXn11lBkhaSz9dg0iq8rgbg28LAhQwa5CJtGXzsM+WdeXFNW5uP29Uz5/mYQ7aTLOEOyjDGhpBrXmKvDb2fLIv/qqQRGwkRnQngHgAXA2AAbwPwNQA7AVwI4JsAfpGZv0tEBOB9AN4I4FkANzLzX4fSV4e06qYP86Q7v72HFQdWVjrEaYM9oW1MQmCj9wH4c2Z+BYBLATwO4HYAjzDzWgCPJL8B4A0A1ib/bgHwgfqLO35UIUTk1p9AtrU2RohUZceYJCESqymlBY32MRarf4noDAAHAKxhkTkRfQ3ABmZ+iojOBbCbmV9ORH+U/L3Dvs6XRxdmbRSlDbRZI1kD4NsAPkREXyGie4jodADnGOGQ/L8quX4KwJPi/qPJsSGI6BYi2kdE+57DqT1h6xgvNr3jXpPOZjFkdVMvk7bXTYgsG643/ZxNCJLlAF4N4APM/CoAP8SpYYwLchwbUaOY+W5mvoyZLzsNK5aO12EnadqVuYjmFdq2IoaYRp7VyzgUKSxrFLGidVPFtUDcuilT9pg2XKb2PS7bURwFcJSZv5z8/jgGguXpZEiD5P8T4vrzxf3nATheU1k7j71tRVaqENShSGF1GtBNRLgYYWlHjzP4hLT9HE13RpKx2GmPmb8F4Ekienly6GoAhwA8DOCG5NgNAD6V/P0wgLfSgCsAfD9kHzGEpn/LUgPbtLrTpozp3rIoc8tOH2mrvfNQVCM4tmW9N0ZJk6vPqxjuNzVr8ysAPkpEjwLoAfhPAH4XwGuJ6DCA1ya/AeAzAOYBHAHwQQBvL5p51h7Y96JiP446o50ZQs+YpfcLbZ5lE9qVrypM2WLfaZ1lc7UPucVFXvJuW2JIE4xjs2UnM/cTe8YrmflaZv4uM/8tM1/NzGuT/59JrmVmfgczv4SZL2HmfTF5+FRNg69BhaJa5SXP/TH3yMV0Lu0odj/YtHL4IsLZdWU+jiyBlWMICbK0evKVMfb60GK+PKS1y5gymS0yDK5375tmr8oo2/lFe22J61AGac8iz5fpoFU0nkied+C7x3c86/M2vWhuXNplm6d/a8HsEWteVpUqbV3qclrDC+3vC+QvZ8yufSG2rjrojN3qQ+4dbOM7HvrgXTaBshfN2Uv7feQJQZGmWZa9r00eOquRTHqEtDp6vCbirmZ9Lt/1sXsjT3IbAlQjyTR1l4Yt8YvM1sRqBVl6Gdcz2mPxvNqIHGvn9auwHdKKjNPlhx3TK/sEQcxQzdRhmRpnnTN9dRqWOytIzNCmDL8DqeLOb+8VmsqM7eFknnaDsD8a1zOa53flm0W4So3D17O7sK+VdWaXLSa9tDyKbLgek09Zw4GmAjxlKb9u2VkQnwSXx40BV5Kl4n1bSIbytO08MVs0rt1941Dacmy8aXV4m1FfaEm7N33kgXujBPXWVQeHPGhD+btCMMZuPBXTA0v/Itd+yDINea18zphQl7HlzovMy6flSK1Ulj+tXGOxaK8O7EV7GsV8GHvsn7d+2rQ3TRX2jCpnbew6r7qN5k1/4m0kUn32bbJt41v/YPfQNnWNe309oSm3vSeL9DNxaTiGTavdfigyH9dzPz2zwnH1aFld63nyqM+hRWw7D7l3C7RnKczfMX5EJnRlmTYyg90mN60edfwrc/+ckEblO5aFzmokctamaM9Sd9Ccmf4ivnjb5YV6+7ICLs/0F/HQjqswtW3PSJpFNJKyNYg870j20mW8YzuNUJplBsSOKXuW+pZli9VIOitIpq/+raVG3iYVvEtovXafiR/aSLU7pILH0HS8kaYxKnabFylWQZ0xPork1Yb30llB4trOwVS4yxYS2gTqoR1XRefrmtExwY5dC+Dmt/ecvhp5GpZsUL71F3lWBRu12dSpSUNqfPK4z47ksjXFPqdr7UjMzIU5vrBxeug922WW2AsB7bRNWnYZ7DTsvELlDA1NjI3Ll6cZhmRZYGlTVBh1dmjj8mzNMoZVilP1xlhF7CyhYVnaiuI6N/xqevg48UMb2cCMtLUbRtlCJG0I5FuRWWZAZFeP1NTGUfJjK9Jb+ihirDUfp6snliueqyJ2XY4UIj6NugyKtpHOChLDTH8RU9v2BJ3I0hzRbFds3zRsWsP29SzyeOyUsyt/n4a1fNf+oCouVXUbMzQz19nnYpjpL+LwhvuXyibL66rL0M6HvnLaZZPPK9OzyxzSLHzDXTPES3v+0HvLI6hMWdMcCg2+vX9t5rf3Cms9nR3avIauXvpdpyqaRd2ucmiV9ZlDy/OB+PgovlXHc71lpU/5lpFemSEOiqYbS1n7BMWUceKnf6UgGXd0FWo6au+qhom3kRhi1rIUpeppQjsYb0wEcjkkAeIif8Vg5x1j7Q8NV8qkqmn6Iun63pWrXcZ4VtvXmN9ZwmFWsS2taiQFqFJTcKWdpdeto4cOPb885ytLaP1HWXFHiuIq+7hoP2WUU4c2HRrapNFUw17YOI0nbl4s1T07z/XKMGXWnw5taqCJXd9d1CVEbJV3+a79UXlnbdQqROLJG0aybFSQFEAu+Cobn8dlaPwaEmyh+9LKb87bDTTkeyDTbHI7SdujNo2YqGtN+eW4aEt4DBUkEYR8LIBqNAJXmmuuHw6KHBMpzVCklzIBfrJoYLL8oc3KshIbkMkwtW1PpmeX15owAjY+n4s2rHnxUbUw77wgiZm1SXNIs+OfykhfofuL4HMck+dsBy9XGe30sjYomfbVb7lp6NwjD9wblYYvVoivLKEZJp9QmN/ec8Zz9c1SpL1zYHRtUVq5Y/w60tboSEKCyY6rEhsnxzVrJGPX5EWNrYK8RqqmIrDZRlafE1pRY+xMfxFbVx0cecbY55493sfa3Tdm2g0vz3vIcl9M2aswYmdJM2+7MvUg1+nIfNPWnMnfOmsjBEno5RWxcDcxu1BXnkYolR3wx7cIramZGpmvXbamZsPSBEgVdTW/vYfDG+4fyVcFyQRN/zZF0ytTlXiyaDfyWp3+LYmsBrSyLPpZDJty3Osb68akl3U6216Z6ru/6B49gDtCexZ7gysP3yyYL5BTTP3ErnIuyzAbWowqWbv7xmA6sswaRT7BaCRtdmzKs5DQfp4sXpcxanoTe/zGUNYitXGkzgWnLnRoU+HQZlxcpKumi0MbfbfD6NCmIK7Nk3xL6vPGhgiRZcoyy9DBp8rnGZb4hIhcMOhT4fM6ieW5Tg5h7Oe0ny9GiBSZ7i869I11kpMxZmKnel3Dx1gmUiOJVeGbHBqVrdLK9LL2uq7pxDzpSNq6aVmbZuKqKItJM/bdqUYSILYBt82+UqQnlELJbkBpvaSphyvv3Dt0/IJ78k+bly1EQtHTshDyHJZ5AeWttfK1szKCNvnSLHv4llkjIaKzAJzPzI+WWpIS0enfeuiijUQZplSNhIh2E9EZRHQ2gAMAPkREdxYtZJPI3ibNXbquRWdZF7r53OCBUdtEFetAYl3kY9zRq6KI+7dPU/PVZWxA51iKvjNXe/IFwSpKlEZCRF9h5lcR0c0YaCPvIaJHmfmVhUtQAZOikTQ9w9B1jeTYlvU4Z+5kp58xjbJtJMuJ6FwAvwjg04VKBoCIfp2IvkpEjxHRDiJ6PhFdRERfJqLDRLSTiJ6XXLsi+X0kOX9h0fyzEur5Xb/LIG1DKBdV9/BZZhxknbl61rSQCFlJ2/wrFlmmqW17nEIkNp+0jdmL0KZQBkC8RvILAP4dgP+XmX+ZiNYA+C/M/K8yZ0g0BeC/A1jHzD8ioo8B+AyANwL4BDM/SER/COAAM3+AiN4O4JXM/G+J6M0A/jdm3hzKoyyNpOkev+3EziyNUz22yYmxynore9am9unfRJB8CcClAH4A4CEA/zeAjwJ4MTMvENGVAH6bmTcR0Wzy9xeJaDmAbwF4EQcKPilDG0WpmrKNrS8jokeI6LHk9yuJ6I48BWPmYwB+D8ATAJ4C8H0A+wF8j5kXksuOAphK/p4C8GRy70Jy/QscZbyFiPYR0b7ncDJP0RQPTUY4a5JJfe48xNpIPgjgXQCeA4Bk6vfNeTJMpo+vAXARgNUATgfwBselRuOgwLlTB5jvZubLmPmy07AiT9GiqToSVhmBZnzkGZv7VOC2jdOzklbHodATMcd9tqEuEitIVjLzXuvYgvPKdH4WwDeY+dvM/ByATwBYD+DMZOgCAOcBOJ78fRTA+QCQnP8pAM/EZlbFR1/1IqqpbXu8Yf6K0lTMjzbi2xUwjVgHsklaZBgrSL5DRC9BogkQ0c9jMCzJwxMAriCilUREAK4GcAjA5wD8fHLNDQA+lfz9cPIbyfldIfuIjQnOM26MQyOMnRZti/EyhqrLOi5G56zECpJ3APgjAK8gomMAfg3AL+fJkJm/DODjAP4awMGkDHcD2ALgNiI6goENxHg73QvgBcnx2wDcnjXPrr68NgvIkHNcU+VoA2XEZmkjUYKEmeeZ+WcBvAjAK5j5Z5j5m3kzZeb3MPMrmPliZv4lZj6Z5HE5M7+UmX+BmU8m1/7P5PdLk/PzefPNi/kQ0vxJ5LVp17nSS0sHGLZL2PFa24Ts2W3tKjb4jyH0bGnPnVXDSFsdnIVQ2bKkG+vlHBtQaaa/WHp7CU7/EtFbmPkBIrrNdZ6ZW+kmr9O/ilIOZU3/np78/5Oef8oYMk4qszIeLA+dZOY/Sv68i5m/XUN5KqHra0Ky0sY4IMoobfKyTSPW2LqHiP6CiG5K/EDGCjuOxjhThTYRGltL+05s3m0xtLqoIjJbVfiESNvqFMjgIk9El2PghHYtBtO1DzLzAxWWLTdqI+ku47RuJwtNB3n2UXqENGbey8y3AbgcA4ewDxcoX610ySbQdC8pacKztYtCBBgPv6EQsWttziCiG4jozwDswcAZ7fJKS1YiXbIJtGnMrHYnxRCrkRwA0APwH5j5Zcy8hZk704pio2or3SfmndelibXZ4dAmVpCsYeZfB/AYEf1ElQWqC9lgfOryzkPDDSYtXL9rW4fZ4/2lf3mQeaTtQmeH1ottiGkfRlHjnmsnwLKGmzHPGMrLrrMYje/pmfRFoTP9RWfZfOV1Hc+6aNCHq01m2e4khtjARhcD+AiAszFYjfttADcw82O5c66QrMZWn6ErjwHMNgZmDf+fJe0QWae8s0w1ptWLPN8mI2Ko/sy5soy5VT530TJmub/UwEZEtAfAu5n5c8nvDQD+EzO3bx4Ko4JknObjmyBN6MiGl7cRF/XlCeVb5VahVRp327S3j+/9lD1rc7oRIgDAzLtxyuu19RRpZFnHqa4dzwxSdbTTbdI3wDQgVxlmj/eHPqbN6/IJA5NHXvU59EGH3q+v/kOYd1NVKAdDW4QIADxxczF7YKxG8kkMVut+JDn0FgCXMfO1hXKvCPUjUZRyKFsjeRsGK38/AeCTyd//Jn/xlDais1RKXoJrbQzM/F0Av1pxWZSGUTvSgK56z1ZJUJAQ0cOh88z8c+UWR1GaR4VIdtKGNldiED/1CxhEfv+v1r/OMk7OQG0m6+bebVrOII3PZRvD6x5GVp1fmiB5MYD/A8DFAN4H4LUAvsPMf8XMf1VpyRrm8Ib7my5CLtomAOd6yzINmeqYycgjFMr2Cal7GFl1fkFBwsyLzPznzHwDgCsAHAGwm4h+pdJStYA2Tc1lQdXydNriINclUmdtkr13/yWABzAIAv1+DGZvFKV02jRzdO11X1j6u22aXttIi9n6YQyGNX+GQfyRVrrE26gfST1o5LnuU5YfyS8BeBmAWzGIkvaD5N/fEdEPyijoOGH3lq7Fcj5ix+V2HtL46Ita78o/Lb/YHta3wAsYDSOQRZvIEm3fRx4twTbmZs3z2Jb1tRmEfc831hHSxgmXRpJnPUaT/gQxvb36O3SbIu83bdFg7PdQ6qK9cUOHNsVWn5YpoCZ5+NOFZy891KJSnDoMdiaPIjMTRYSIrfaPw4dU1VBFPvu4b7ieRucFSR2zALENMeYDLVJeuVI3ZGvJihFOdpquj8M1bT57vF/auD7ProVp98syy/O+vPI8yzgI1CJ0fmhTRqwKOw21TQzogureFtoUAEqiQ5sEnxDJMsyw02ibEPH1kKH4J1koOiRz9exNqfplzWblxafluIRIFWXxhYAsSucFiY8L7hkVMDHDijY6JknHKYMdf7QK4RcTu9Rgf8BFNBnXMC12SJjW68t6quJdZ9GOy3xnpn7messqaQudH9q0jS4Ni9o0tKki3GKX3lVedGhTI1l6rqYaZuy2nFlYvmt/lOGxDoN3FYvSJl2IZGFiBIn9IUn1uKgK24YGl/ZB5415Cvg9eOe39zIZCJsYFrq2YghRtddoHqEae0/MzFxVQl2HNoqieNGhjaDrzkBVUmYPXXZv38Y1J5PKRAiSqihDTWzLx+ATtmX6NpTtJ9FGvwvgVLuoahjRpihyhsoECRHdR0QniOgxcexsIvosER1O/j8rOU5E9H4iOkJEjxLRq8U9NyTXHyaiG/KUpeieHT7KMPC15WNoevaljdPqeTHtoqqoZG0MulWlRnI/gNdbx24H8AgzrwXwSPIbAN4AYG3y7xYAHwAGggfAewC8BsDlAN5jhE8W5EZHbQqco5wi7+59TSC1yDZqB7GUKbwrEyTM/HkAz1iHrwHw4eTvDwO4Vhz/Yx7wJQBnEtG5ADYB+CwzP5NsifFZjAqnVGSv304MAWAAABTlSURBVOYtF8atV27a9lSFFmUPNV3PKNtTldpB1e2hzNnGum0k5zDzUwCQ/L8qOT4F4Elx3dHkmO94J2nDNHIWmh4OVcHJS58d+t3kM45Te2iLsZUcxzhwfDQBoluIaB8R7XsOJ0stnOLuHZvWSGzKGGbIj3cchsFl+pgUoW5B8nQyZEHy/4nk+FEA54vrzgNwPHB8BGa+m5kvY+bLTkP8GpAsFdyWGZaySRMIPlfxmN7arrMyw1JKZvqLpQ8zHtpxVanpVYFrqD57vD/Sru84ccnIdWW257oFycMAzMzLDQA+JY6/NZm9uQLA95OhzyyA1xHRWYmR9XXJsULISo5pfFmCBdkfSlqPUTT+SNUsbJxeEiLz23uZx+3nzJ3SDhc2TmPFgZWY394LPrddz7485eZbIdtXFs1J1qk9zBkXNq3uYdPq3pCgcNWPrOei2ldlnq1EtAPABgAvBPA0BrMvDwH4GIALADwB4BeY+RkiIgD/DQND6rMA/g0z70vSeRsGm3QBwO8w84fS8i7bs7WKBWFVE1pwlncxmn1fmxbtVbHAThftaczWVEGSpZFkCTpjC522BqxR3FQtPEz6MZ1TGwSZChJda1M5bdJIlGrQtTYBxsEarzRDm9rGOPkVTaQgsenqbExb8G2GpfUepulhTRYmUpDYY1Njw/DNuDTRM7jyTOst5YyDb9q1zGd55IF7o65zhYKc6y2rxXaUVVjZQb7Lpk6NRz57KB5PGaiNpKW4DG2zx/utWrDVhDGwzhk03/O1wQialbxtR20kY46roRYRIlX0rmb2IY0ye+E6p+F9z1f1rI5NjPaQ9n6ztJ0872tiBUnbDVkuJyq7zFnU9qoaf8yHPQ4eoga7juv2H8rbgZT5fvM888QKkrarpq5pVbvMdfmn+IRuyGNU9mquctZlK1CDbj1MrCApk7Y01qq0rM3r9jvV65Cx1e7VbMGRt6fPKoCKCNumtVbfVqltRAVJCWRprFU2iiqHLy71OkblNlqLERw+oWsf99VT2UMNOx9ZhzIglqQuAWPKMg7LM3TWRsmNerZ2H521USony5ad40KbouZX4ftTFaqRTBh1r2RWrWW8UY0kknEwZJVJ3ePtcRYiVRrRxzlotIuJFyRzvWWpqmMR/422UaaaXCTU4jh8SCawkf2cZbz/q99yU+E00qiznerQRgm6fI9jUCelPHRoM4Hk7eXzbjC+sHG600PDcTBytgUVJGNGSFjkXYuT94N5emZFa7SVLgu0cUAFyZhRxerfvI5sMrBz07RFoE0qKkg6RN29cqwfybgOEXyereNCncZWFSQlUvaLi/0AzXV198qxSwPsyPPjQluDdse2izrLr4KkRMp+cbFDjravZJaMk19JHZpUHgN53q1EYslTJhUkkYyret4Gxik+q9SYyhbQrjqoK+Jd1Z2N+pEoiuJF/UgmgJCWVKYGlXcvXtffsbRBA5zf3nM+exnaVZkaWhvqSgVJS4lpHHJPXt+5MvDZfkKGU2n49W10HSJr+asw4vrKEGsLc72XIvtI+6hjCJaGDm2UaGaP93HHiUsyzQ6Ns4t9ldHi27YjgA8d2iils2l1b0goxGgB4ypEgGoNlFUKkSaGOipIKqbLrtuxU7lVrvSN/Whi38M4rEpOI48ALNpOVZBEUMQw1rYeuQnBVmXvG/vRxL6HcRhu2JShgRRtpypIPEjh0VYPxxA+4dc2wdZmmpoNyZpvGxwSVZB4iBEebVaDswi/2BkPu4E34e7e5aGioQ2CISsqSAqwdveNTRehFJ64Oe7jtBt4E+7uOw+5hVcbfCkmGZ3+VXKjgZ27j07/Vkgd60bGaZVsl1DNJh8qSHJQh/G1Sz39OAnFcbRPlE0eYVqZICGi+4joBBE9Jo79FyL6GyJ6lIg+SURninPvIqIjRPQ1Itokjr8+OXaEiG6vqrxZqMPg53qZZeRbZtljbSs+oTguK4JtfB+aXbeu5xsHjSePMK3MRkJE/xzA3wP4Y2a+ODn2OgC7mHmBiLYBADNvIaJ1AHYAuBzAagB/CeBlSVJfB/BaAEcBzAG4jpkPhfIeJxtJjAt5W93M22YjKdvtvEoX+XGhcRsJM38ewDPWsb9g5oXk55cAnJf8fQ2AB5n5JDN/A8ARDITK5QCOMPM8M/8DgAeTa1tBlt7dN1UcIyDqFiKy1ww9Y6wQiU2vaG89Ls5kRZ6zKpeDcfZsfRuAP0v+ngLwpDh3NDnmOz4CEd1CRPuIaN9zaE9QYsO4NHJgWLUtQ4jFpte23r+q8hRJt6p2NJaerUT0bgALAD5qDjku48Dx0YPMdzPzZcx82WmodnPrOmKktsV+ECpHyIg6DraAqmnLO6yD2gUJEd0A4E0A/jWfMtAcBXC+uOw8AMcDx2vF/ijq6Dnb4pY/tW2PVyiEhjZV1lFdQqqout+Wd1gHtQoSIno9gC0Afo6ZnxWnHgbwZiJaQUQXAVgLYC8GxtW1RHQRET0PwJuTa2ulzm0Jqp4qzZr+wsbpEaFQRk9b5COtawhUVONUjaQEiGgHgC8CeDkRHSWimwD8NwA/CeCzRNQnoj8EAGb+KoCPATgE4M8BvIOZFxPD7DsBzAJ4HMDHkmtrpc6exfTyMR9a6BpfI846y7J81/4RDcDUR6xm4DIQlj0srPOjjRWCk6SRqIu8kpu2Tf+WQZlT7UWnj49tWV+6MMr6fI1P/3aZMsforqFGGb1r271J27qKN+Yji61blxDJ8txVaDRVTRCoRqIoCbZDmzqkqUaiKJkxQsRonJMuRLKggiQjbVXJm6KLMxMqQLKjgiQjc71lnfx48rCwcXqiZiYUPypIcqAfz4Cuzdi0iXHTfFWQKIqHJt3827jaO4QKkgI0Gfy5DcOrtk8xF0VtJfGoIClAk8Gf04ZXdfSmOrRRDCpICrDm+n4rNAMXbepNZ/qLE7sauG3PXZUWrQ5piuKhqIt6Fxza1CGtJZgeqYqeadws+2XTdk3HCJFxsCUVbUuqkZRME71QFYu7yqQLPXMTtKHeYjUSFSSKImi7UK4bHdqMOVUZcX1DgXFQv+tgHIRIGw38Kkhagm1Nr6pB+1TlsqZy27yxug9fmbPYX+q0V7VR2KkgaQlNRpkvUxup8jmq6ol9Zc5inxg3T9SyUUGijI1jWRt7YmWAChJFUQqjgqRBylbVJ33tT1m02Telrej0r5KbLgZ/VobR6V+lFEK9c9eFSJe0rKpRQaIEadqzsknUuBuPChJFUQqjgqQh0gx6qlY3g9Z7PtTYqiiKFzW2KopSGypIFEUpjAqSFqArb5VxRwVJC+i6P4bSfVSQKIpSGBUkiqIURgWJoiiFUUGiKEphVJAoilIYFSQNoa7YSpeoTJAQ0X1EdIKIHnOc+00iYiJ6YfKbiOj9RHSEiB4loleLa28gosPJvxuqKm/d6MpSpUnKDoJVpUZyP4DX2weJ6HwArwXwhDj8BgBrk3+3APhAcu3ZAN4D4DUALgfwHiI6q8IyKxlQR7rxpewg3ZUJEmb+PIBnHKfeC+C3AMjVgtcA+GMe8CUAZxLRuQA2AfgsMz/DzN8F8Fk4hJPSDOpIpxhqtZEQ0c8BOMbMB6xTUwCeFL+PJsd8x11p30JE+4ho33M4WWKplUlFY7fGU5sgIaKVAN4N4P90nXYc48Dx0YPMdzPzZcx82WlYkb+gNdN2o+uxLesndne+FQdWNl2EsaFOjeQlAC4CcICIvgngPAB/TUQvxkDTOF9cex6A44HjnaHtRtepbXsmNtxi299NGnV2UrUJEmY+yMyrmPlCZr4QAyHxamb+FoCHAbw1mb25AsD3mfkpALMAXkdEZyVG1tclxzpD2zWSEGojaTd1CsIqp393APgigJcT0VEiuilw+WcAzAM4AuCDAN4OAMz8DID/CGAu+fcfkmOdoaqXXYeAih3a1LkvbpmojSSeKmdtrmPmc5n5NGY+j5nvtc5fyMzfSf5mZn4HM7+EmS9h5n3iuvuY+aXJvw9VVd6uUaaAGmetKSvyWZse0o2TANaYrYqieNGYrUrldH3WJg+TOhxSQdJS2tQg21SWttP0cKgpdGijKIoXHdooSkG6pImVvUjPRgWJonjo0jCl7EV6NipIFEUpjAoSRVEKo4KkQ1Q9Dp40umQjqRoVJB2i6nGwTdc9XrtkI6kaFSRKbs6Zqz/uyzi5jU8SKkg6REgVb1JNLzPvud6y0tJKo+saV5moQ5rSKhY2TrcqPMH89l7jQ5yZ/mKtAlSiDmlK5VSx1qZNQgRoh52kKSGSBRUkSm7a9tErzaGCRFGUwqggURQPVRqou+ajooJEURqgDbaXMlFBoigeuvaxV4kKEkVRCqOCRMmNOmwpBhUkSm6mtu1pjct6W8oxqaggUQrNIKQ5S9WltYyD01aXUUGiVGpUHPdtL5U4VJAo0biGDzqkUAAVJEoGXMMHHVIogAoSpQC6QZZiUEGiKEphVJAoilIYFSSKohRGBYmiKIVRQaIoSmFUkCiKMkQe3yAVJIqiDJHHN0gFiaIohalMkBDRfUR0goges47/ChF9jYi+SkT/lzj+LiI6kpzbJI6/Pjl2hIhur6q8iqLkp0qN5H4Ar5cHiOh/AXANgFcy808D+L3k+DoAbwbw08k9dxHRMiJaBuAPALwBwDoA1yXXKmOCrsVpB1XHiF1eVcLM/HkiutA6/MsAfpeZTybXnEiOXwPgweT4N4joCIDLk3NHmHkeAIjoweTaQ1WVW4knZjsKXYvTDqoOG1m3jeRlAK4ioi8T0V8R0UxyfArAk+K6o8kx33FFqZyuRXqvkso0kkB+ZwG4AsAMgI8R0RoA5LiW4RZ0zj1GiegWALcAwPOxspTCKpONBn+Op25BchTAJ3iw4fBeIvpHAC9Mjp8vrjsPwPHkb9/xIZj5bgB3A4O9f0sut6IoAeoe2jwEYCMAENHLADwPwHcAPAzgzUS0goguArAWwF4AcwDWEtFFRPQ8DAyyD9dcZkVRUqhMIyGiHQA2AHghER0F8B4A9wG4L5kS/gcANyTayVeJ6GMYGFEXALyDmReTdN4JYBbAMgD3MfNXqyqzoij5oMF33C3OoLP5NXR108VQlLHnL/nj+5n5srTr1LNVUZTCqCBRFA/qTBePChJFUQqjgkRRPOw8pMGtY1FBoigeLrhH3ftjUUGiKB6euFltJLGoIFEUD+oiH48KEkVRCqOCRFEaoks7FaogUZSGiInnMi6oIFEUpTAqSBRFKYwKEkVRCtPJ1b9E9G0AP8Qg1sm48EKMV3kBLXMdNF3ef8LML0q7qJOCBACIaF/M8ue2MG7lBbTMdTAu5dWhjaIohVFBoihKYbosSO5uugAZGbfyAlrmOhiL8nbWRqIoSn10WSNRFKUmVJAoilKYzgkSIno9EX2NiI4Q0e1Nl8dAROcT0eeI6HEi+ioR3ZocP5uIPktEh5P/z0qOExG9P3mOR4no1Q2VexkRfYWIPp38vijZcvUwEe1M9htCsifRzqS8X3bs+1xXec8koo8T0d8kdX1lm+uYiH49aQ+PEdEOInp+2+vYRacECREtA/AHAN4AYB2A64hoXbOlWmIBwG8w8z/FYMvSdyRlux3AI8y8FsAjyW9g8Axrk3+3APhA/UUGANwK4HHxexuA9ybl/S6Am5LjNwH4LjO/FMB7k+ua4H0A/pyZXwHgUgzK3so6JqIpAL8K4DJmvhiDvZvejPbX8SjM3Jl/AK4EMCt+vwvAu5oul6esnwLwWgBfA3BucuxcAF9L/v4jANeJ65euq7GM52Hw4W0E8GkM9mj+DoDldn1jsInZlcnfy5PrqObyngHgG3a+ba1jAFMAngRwdlJnnwawqc117PvXKY0Ep16M4WhyrFUkKumrAHwZwDnM/BQAJP+vSi5rw7P8PoDfAvCPye8XAPgeMy84yrRU3uT895Pr62QNgG8D+FAyHLuHiE5HS+uYmY8B+D0ATwB4CoM6249217GTrgkSchxr1fw2Ef0EgD8F8GvM/IPQpY5jtT0LEb0JwAlmlkEzQmVqQ90vB/BqAB9g5ldhsN4qZCdruo7PAnANgIsArAZwOgbDLV+Z2lDHTromSI4COF/8Pg/A8YbKMgIRnYaBEPkoM38iOfw0EZ2bnD8XwInkeNPP8s8A/BwRfRPAgxgMb34fwJlEZPaMlmVaKm9y/qcAPFNjeU0ZjjLzl5PfH8dAsLS1jn8WwDeY+dvM/ByATwBYj3bXsZOuCZI5AGsTq/fzMDBcPdxwmQAMZggA3AvgcWa+U5x6GMANyd83YGA7McffmswsXAHg+0Y9rwNmfhczn8fMF2JQj7uY+V8D+ByAn/eU1zzHzyfX19pbMvO3ADxJRC9PDl2Nwcb0raxjDIY0VxDRyqR9mPK2to69NG2kqcCA9UYAXwfwPwC8u+nyiHL9DAZq6KMA+sm/N2Iwxn0EwOHk/7OT6wmDGaj/AeAgBpb9psq+AcCnk7/XANgL4AiAPwGwIjn+/OT3keT8mobK2gOwL6nnhwCc1eY6BvDvAfwNgMcAfATAirbXseufusgrilKYrg1tFEVpABUkiqIURgWJoiiFUUGiKEphVJAoilIYFSRKKRDRhUT0mHXst4noN5sqk1IfKkiU1iK8O5WWo4JEqRwi+lUiOpTE/HgwOXY6Ed1HRHPJArtrkuM3EtGfENH/A+AviOhcIvo8EfWTmB1XNfowihOV+Eod3A7gImY+SURnJsfejYGL99uSY3uJ6C+Tc1cCeCUzP0NEv4HBMvrfSeLNrKy/+EoaKkiUsvC5SJtlAR8loocwcFsHgNdhsCjQ2FCeD+CC5O/PMrNZjDYH4L5kweNDzNwvv+hKUXRoo5TF32KwrkVyNgbBd/5XDNa0TAPYn9g+CMC/YuZe8u8CZjaR2H5oEmDmzwP45wCOAfgIEb214udQcqCCRCkFZv57AE8R0dXAIBYtgNcD+O8Azmfmz2EQJOlMAD+BQbSvX0lWvYKIXuVKl4j+CQZxUT6IwerpRmLXKmF0aKOUyVsB/AER/dfk97/HYKn854jopzDQQt7LzN8jov+IQXyTRxNh8k0Ab3KkuQHA/05EzwH4+yQPpWXo6l9FUQqjQxtFUQqjgkRRlMKoIFEUpTAqSBRFKYwKEkVRCqOCRFGUwqggURSlMP8/QGLDK/9rnD4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(6,6))\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": [
    "### 2.2 协同过滤学习算法\n",
    "\n",
    "接下来，我们将实施协同过滤的代价函数。 直觉上，“代价”是指一组电影评级预测偏离真实预测的程度。它基于文本中称为$X$和$Theta$的两组参数矩阵。 这些“展开”到“参数”输入中，以便稍后可以使用SciPy的优化包。 请注意，已经在注释中包含数组/矩阵形状（对于我们在本练习中使用的数据），以帮助说明矩阵交互如何工作。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.2.1 代价函数\n",
    "\n",
    "$$\n",
    "J\\left(x^{(1)},...,x^{(n_m)},\\theta^{(1)},...,\\theta^{(n_u)}\\right)\n",
    "=\\dfrac{1}{2}\\sum_{(i,j):r(i,j)=1} \\left((\\theta^{(j)})^Tx^{(i)}-y^{(i,j)}\\right)^2\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "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",
    "    # X和Theta本来是以vector形式传入，下面把它们还原为ndarray\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",
    "    # 初始化\n",
    "    J = 0\n",
    "    \n",
    "    # 计算代价\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": [
    "为了测试这一点，我们提供了一组我们可以评估的预训练参数。\n",
    "\n",
    "为了保持评估时间的少点，我们将只看一小段数据。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((1682, 10), (943, 10))"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "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": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "22.224603725685675"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "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",
    "# 拼接X_sub & Theta_sub为一个vector\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": [
    " #### 2.2.2 梯度下降的协同过滤算法\n",
    " \n",
    " 接下来我们需要实现梯度计算。 就像我们在练习4中使用神经网络实现一样，我们将扩展代价函数来计算梯度。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "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",
    "    # X和Theta本来是以vector形式传入，下面把它们还原为ndarray\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",
    "    # 初始化\n",
    "    J = 0\n",
    "    X_grad = np.zeros(X.shape)  # (1682, 10)\n",
    "    Theta_grad = np.zeros(Theta.shape)  # (943, 10)\n",
    "    \n",
    "    # 计算代价\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",
    "    # 计算梯度\n",
    "    X_grad = error * Theta\n",
    "    Theta_grad = error.T * X\n",
    "    \n",
    "    # 将梯度矩阵解析为vector\n",
    "    grad = np.concatenate((np.ravel(X_grad), np.ravel(Theta_grad)))\n",
    "    \n",
    "    return J, grad"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(22.224603725685675,\n",
       " array([ -2.52899165,   7.57570308,  -1.89979026,  -0.56819597,\n",
       "          3.35265031,  -0.52339845,  -0.83240713,   4.91163297,\n",
       "         -0.76677878,  -0.38358278,   2.26333698,  -0.35334048,\n",
       "         -0.80378006,   4.74271842,  -0.74040871, -10.5680202 ,\n",
       "          4.62776019,  -7.16004443,  -3.05099006,   1.16441367,\n",
       "         -3.47410789,   0.        ,   0.        ,   0.        ,\n",
       "          0.        ,   0.        ,   0.        ]))"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "J, grad = cost(params, Y_sub, R_sub, features)\n",
    "J, grad\n",
    "\n",
    "# J 应该为 22.22460"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "运行梯度检验，确保cost函数实现的没问题。\n",
    "\n",
    "观察梯度检验输出的gradient和上面输出的gradient的对比。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ -2.52899165   7.57570308  -1.89979026  -0.56819597   3.35265031\n",
      "  -0.52339845  -0.83240713   4.91163297  -0.76677878  -0.38358278\n",
      "   2.26333698  -0.35334048  -0.80378006   4.74271842  -0.74040871\n",
      " -10.5680202    4.62776019  -7.16004443  -3.05099006   1.16441367\n",
      "  -3.47410789   0.           0.           0.           0.\n",
      "   0.           0.        ]\n"
     ]
    }
   ],
   "source": [
    "def compute_numerical_gradient(X_sub, Theta_sub, Y_sub, R_sub, features):\n",
    "    e = np.exp(-4)\n",
    "    \n",
    "    theta = np.concatenate((np.ravel(X_sub), np.ravel(Theta_sub)))\n",
    "    numgrad = np.zeros(len(theta))\n",
    "    perturb = np.zeros(len(theta))\n",
    "\n",
    "    for i in range(len(theta)):\n",
    "        perturb[i] = e\n",
    "        \n",
    "        loss1, _ = cost(theta - perturb, Y_sub, R_sub, features)\n",
    "        loss2, _ = cost(theta + perturb, Y_sub, R_sub, features)\n",
    "        \n",
    "        numgrad[i] = (loss2 - loss1) / (2*e)\n",
    "        perturb[i] = 0\n",
    "        \n",
    "        # print(loss1, loss2, numgrad[i])\n",
    "    \n",
    "    return numgrad\n",
    "\n",
    "print(compute_numerical_gradient(X_sub, Theta_sub, Y_sub, R_sub, features))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.2.4 正则化代价函数\n",
    "\n",
    "\n",
    "下一步是在代价和梯度计算中添加正则化。 我们将创建一个最终的正则化版本的功能（请注意，此版本包含一个额外的“学习率”参数，在文本中称为“lambda”）。\n",
    "\n",
    "\n",
    "计算公式为：\n",
    "$$\n",
    "J\\left(x^{(1)},...,x^{(n_m)},\\theta^{(1)},...,\\theta^{(n_u)}\\right)\n",
    "=\\dfrac{1}{2}\\sum_{(i,j):r(i,j)=1} \\left((\\theta^{(j)})^Tx^{(i)}-y^{(i,j)}\\right)^2\n",
    "+ \\left(\\dfrac{\\lambda}{2} \\sum_{j=1}^{n_u} \\sum_{k=1}^{n} (\\theta_k^{(j)})^2 \\right)\n",
    "+ \\left(\\dfrac{\\lambda}{2} \\sum_{i=1}^{n_m} \\sum_{k=1}^{n} (x_k^{(i)})^2\\right)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "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",
    "     # X和Theta本来是以vector形式传入，下面把它们还原为ndarray\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",
    "    # 初始化\n",
    "    J = 0\n",
    "    X_grad = np.zeros(X.shape)  # (1682, 10)\n",
    "    Theta_grad = np.zeros(Theta.shape)  # (943, 10)\n",
    "    \n",
    "    # 计算代价\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",
    "    # 计算代价的正则化部分\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",
    "    # 计算正则化的梯度\n",
    "    X_grad = (error * Theta) + (learning_rate * X)\n",
    "    Theta_grad = (error.T * X) + (learning_rate * Theta)\n",
    "    \n",
    "    # 将梯度矩阵解析为vector\n",
    "    grad = np.concatenate((np.ravel(X_grad), np.ravel(Theta_grad)))\n",
    "    \n",
    "    return J, grad"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(31.34405624427422,\n",
       " array([ -0.95596339,   6.97535514,  -0.10861109,   0.60308088,\n",
       "          2.77421145,   0.25839822,   0.12985616,   4.0898522 ,\n",
       "         -0.89247334,   0.29684395,   1.06300933,   0.66738144,\n",
       "          0.60252677,   4.90185327,  -0.19747928, -10.13985478,\n",
       "          2.10136256,  -6.76563628,  -2.29347024,   0.48244098,\n",
       "         -2.99791422,  -0.64787484,  -0.71820673,   1.27006666,\n",
       "          1.09289758,  -0.40784086,   0.49026541]))"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "J, grad = cost(params, Y_sub, R_sub, features, 1.5)\n",
    "J, grad\n",
    "\n",
    "# J 应该为 31.344056"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "这个结果再次与练习代码的预期输出相匹配，所以看起来正则化是正常的。 \n",
    "\n",
    "### 2.3 电影推荐\n",
    "\n",
    "在我们训练模型之前，我们有一个最后步骤， 我们的任务是创建自己的电影评分，以便我们可以使用该模型来生成个性化的推荐。 为我们提供一个连接电影索引到其标题的文件。 接着我们将文件加载到字典中。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "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": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'Toy Story (1995)'"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "movie_idx[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们将使用练习中提供的评分。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Rated Toy Story (1995) with 4 stars.\n",
      "Rated Twelve Monkeys (1995) with 3 stars.\n",
      "Rated Usual Suspects, The (1995) with 5 stars.\n",
      "Rated Outbreak (1995) with 4 stars.\n",
      "Rated Shawshank Redemption, The (1994) with 5 stars.\n",
      "Rated While You Were Sleeping (1995) with 3 stars.\n",
      "Rated Forrest Gump (1994) with 5 stars.\n",
      "Rated Silence of the Lambs, The (1991) with 2 stars.\n",
      "Rated Alien (1979) with 4 stars.\n",
      "Rated Die Hard 2 (1990) with 5 stars.\n",
      "Rated Sphere (1998) with 5 stars.\n"
     ]
    }
   ],
   "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": [
    "我们可以将自己的评级向量添加到现有数据集中以包含在模型中。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((1682, 944), (1682, 944), (1682, 1))"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "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": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((1682, 10), (944, 10), (26260,))"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "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": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5.507036456515984e-19"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "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": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "     fun: 38965.44283137412\n",
       "     jac: array([ 0.02542188,  0.039419  ,  0.04361765, ..., -0.00874763,\n",
       "        0.00236952,  0.00153375])\n",
       " message: 'Maximum number of iterations has been exceeded.'\n",
       "    nfev: 148\n",
       "     nit: 100\n",
       "    njev: 148\n",
       "  status: 1\n",
       " success: False\n",
       "       x: array([ 4.77228432e-01,  8.25795179e-03, -1.80952067e-03, ...,\n",
       "       -1.90029107e-04,  1.08414339e-02, -9.17601268e-02])"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "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": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((1682, 10), (944, 10))"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "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": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1682, 1)"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "predictions = X * Theta.T \n",
    "my_preds = predictions[:, -1] + Ymean\n",
    "my_preds.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[5.00000018],\n",
       "        [5.00000015],\n",
       "        [5.00000013],\n",
       "        [5.00000008],\n",
       "        [5.00000007],\n",
       "        [5.00000007],\n",
       "        [4.99999978],\n",
       "        [4.99999935],\n",
       "        [4.99999908],\n",
       "        [4.99999822]])"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "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": 48,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[1188],\n",
       "        [ 813],\n",
       "        [1652],\n",
       "        ...,\n",
       "        [1307],\n",
       "        [1574],\n",
       "        [ 857]], dtype=int64)"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "idx = np.argsort(my_preds, axis=0)[::-1]\n",
    "idx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Top 10 movie predictions:\n",
      "Predicted rating of 5.000000183388653 for movie Prefontaine (1997).\n",
      "Predicted rating of 5.0000001535211975 for movie Great Day in Harlem, A (1994).\n",
      "Predicted rating of 5.000000134100087 for movie Entertaining Angels: The Dorothy Day Story (1996).\n",
      "Predicted rating of 5.000000082362455 for movie Someone Else's America (1995).\n",
      "Predicted rating of 5.000000074235328 for movie Santa with Muscles (1996).\n",
      "Predicted rating of 5.000000067666932 for movie They Made Me a Criminal (1939).\n",
      "Predicted rating of 4.999999777214146 for movie Saint of Fort Washington, The (1993).\n",
      "Predicted rating of 4.999999345248747 for movie Marlene Dietrich: Shadow and Light (1996) .\n",
      "Predicted rating of 4.999999079447295 for movie Aiqing wansui (1994).\n",
      "Predicted rating of 4.999998221290874 for movie Star Kid (1997).\n"
     ]
    }
   ],
   "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
}