2.logistic_regression/ex2.logistic_regression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "6HPvYgkOjlHg"
   },
   "source": [
    "# 逻辑回归 Logistic Regression\n",
    "\n",
    "此Notebook是配合Andrew Ng \"Machine Leanring\"中[逻辑回归](https://github.com/loveunk/machine-learning-deep-learning-notes/blob/master/machine-learning/logistic-regression.md)部分学习使用。\n",
    "\n",
    "测试用python版本为3.6\n",
    "* 机器学习路径：https://github.com/loveunk/machine-learning-deep-learning-notes/\n",
    "* 内容正文综合参考网络资源，使用中如果有疑问请联络：www.kaikai.ai"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "2bSRbRGTjlHi"
   },
   "source": [
    "在这一次练习中，我们将要实现逻辑回归并且应用到一个分类任务。我们还将通过将正则化加入训练算法，来提高算法的鲁棒性，并用更复杂的情形来测试它。\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "WXk4cLxBjlHj"
   },
   "source": [
    "## 准备数据"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "nfrEaVuajlHl"
   },
   "source": [
    "在训练的初始阶段，我们将要构建一个逻辑回归模型来预测，某个学生是否被大学录取。设想你是大学相关部分的管理者，想通过申请学生两次测试的评分，来决定他们是否被录取。现在你拥有之前申请学生的可以用于训练逻辑回归的训练样本集。对于每一个训练样本，你有他们两次测试的评分和最后是被录取的结果。为了完成这个预测任务，我们准备构建一个可以基于两次测试评分来评估录取可能性的分类模型。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "2zO2vzfujlHm"
   },
   "source": [
    "让我们从检查数据开始。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "meEb58oJjlHn"
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "plt.style.use('fivethirtyeight')\n",
    "from sklearn.metrics import classification_report # 这个包是评价报告"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "zqEIXn3sjlHr"
   },
   "outputs": [],
   "source": [
    "path = 'ex2data1.txt'\n",
    "data = pd.read_csv(path, header=None, names=['exam1', 'exam2', 'admitted'])\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "8r4bOK-K8J4i"
   },
   "outputs": [],
   "source": [
    "data.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "bvQFnaZojlH1"
   },
   "source": [
    "让我们创建两个分数的散点图，并使用颜色编码来可视化，如果样本是正的（被接纳）或负的（未被接纳）。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "4VZ-TIUmjlH3"
   },
   "outputs": [],
   "source": [
    "positive = data[data['admitted'].isin([1])]\n",
    "negative = data[data['admitted'].isin([0])]\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(12,8))\n",
    "ax.scatter(positive['exam1'], positive['exam2'], s=50, c='b', marker='o', label='Admitted')\n",
    "ax.scatter(negative['exam1'], negative['exam2'], s=50, c='r', marker='x', label='Not Admitted')\n",
    "ax.legend()\n",
    "ax.set_xlabel('Exam 1 Score')\n",
    "ax.set_ylabel('Exam 2 Score')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "2dp3JtvpjBYR"
   },
   "outputs": [],
   "source": [
    "def get_X(df):#读取特征\n",
    "#     \"\"\"\n",
    "#     use concat to add intersect feature to avoid side effect\n",
    "#     not efficient for big dataset though\n",
    "#     \"\"\"\n",
    "    ones = pd.DataFrame({'ones': np.ones(len(df))})#ones是m行1列的dataframe\n",
    "    data = pd.concat([ones, df], axis=1)  # 合并数据，根据列合并\n",
    "    return data.iloc[:, :-1].as_matrix()  # 这个操作返回 ndarray,不是矩阵\n",
    "\n",
    "\n",
    "def get_y(df):#读取标签\n",
    "#     '''assume the last column is the target'''\n",
    "    return np.array(df.iloc[:, -1])#df.iloc[:, -1]是指df的最后一列\n",
    "\n",
    "\n",
    "def normalize_feature(df):\n",
    "#     \"\"\"Applies function along input axis(default 0) of DataFrame.\"\"\"\n",
    "    return df.apply(lambda column: (column - column.mean()) / column.std())#特征缩放"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "mMuhs2-hjBYW"
   },
   "outputs": [],
   "source": [
    "X = get_X(data)\n",
    "print(X.shape)\n",
    "\n",
    "y = get_y(data)\n",
    "print(y.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "8R0uQXSZjlH7"
   },
   "source": [
    "看起来在两类间，有一个清晰的决策边界。现在我们需要实现逻辑回归，那样就可以训练一个模型来预测结果。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "R_GK_8JhjlH8"
   },
   "source": [
    "## Sigmoid 函数\n",
    "$g$ 代表一个常用的逻辑函数（logistic function）为S形函数（Sigmoid function），公式为： \\\\[g\\left( z \\right)=\\frac{1}{1+{{e}^{-z}}}\\\\] \n",
    "合起来，我们得到逻辑回归模型的假设函数： \n",
    "\t\\\\[{{h}_{\\theta }}\\left( x \\right)=\\frac{1}{1+{{e}^{-{{\\theta }^{T}}X}}}\\\\] "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "6yd5Ddl2jlH-"
   },
   "outputs": [],
   "source": [
    "def sigmoid(z):\n",
    "    return 1 / (1 + np.exp(-z))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "4SG4dpr3jlID"
   },
   "source": [
    "让我们做一个快速的检查，来确保它可以工作。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "FN4JMcM2jlIE"
   },
   "outputs": [],
   "source": [
    "nums = np.arange(-10, 10, step=1)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(12,8))\n",
    "ax.plot(nums, sigmoid(nums), 'r')\n",
    "ax.set_xlabel('z', fontsize=18)\n",
    "ax.set_ylabel('g(z)', fontsize=18)\n",
    "ax.set_title('sigmoid function', fontsize=18)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "apEHwiQmvoH_"
   },
   "source": [
    "## Cost Function 代价函数"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "Zw3ZmdqVjlII"
   },
   "source": [
    "棒极了！现在，我们需要编写代价函数来评估结果。\n",
    "代价函数：\n",
    "$J\\left( \\theta  \\right)=\\frac{1}{m}\\sum\\limits_{i=1}^{m}{[-{{y}^{(i)}}\\log \\left( {{h}_{\\theta }}\\left( {{x}^{(i)}} \\right) \\right)-\\left( 1-{{y}^{(i)}} \\right)\\log \\left( 1-{{h}_{\\theta }}\\left( {{x}^{(i)}} \\right) \\right)]}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "R0Y4hVpkjlIJ"
   },
   "outputs": [],
   "source": [
    "def cost(theta, X, y):\n",
    "    return np.mean(-y * np.log(sigmoid(X @ theta)) - (1 - y) * np.log(1 - sigmoid(X @ theta)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "XEBSMFxujlIS"
   },
   "source": [
    "让我们来检查矩阵的维度来确保一切良好。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "6Ax8oM2AjlIT"
   },
   "outputs": [],
   "source": [
    "theta = np.zeros(3)\n",
    "theta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "8OY_-JQ8jlIX"
   },
   "outputs": [],
   "source": [
    "X.shape, theta.shape, y.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "ZZxmEmkwjlIa"
   },
   "source": [
    "让我们计算初始化参数的代价函数(theta为0)。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "fXdMxsdrjlIb"
   },
   "outputs": [],
   "source": [
    "cost(theta, X, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "SSGkYiVXjlIf"
   },
   "source": [
    "看起来不错，接下来，我们需要一个函数来计算我们的训练数据、标签和一些参数thata的梯度。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "tI-AX8knjlIi"
   },
   "source": [
    "## Gradient descent 梯度下降\n",
    "* 这是批量梯度下降（batch gradient descent）  \n",
    "* 转化为向量化计算： $\\frac{1}{m} X^T( Sigmoid(X\\theta) - y )$\n",
    "$$\\frac{\\partial J\\left( \\theta  \\right)}{\\partial {{\\theta }_{j}}}=\\frac{1}{m}\\sum\\limits_{i=1}^{m}{\\left({h_{\\theta }}\\left( {{x}^{(i)}} \\right)-{{y}^{(i)}}\\right)x_{_{j}}^{(i)}}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "D1SA0lHXjlIj"
   },
   "outputs": [],
   "source": [
    "# Gradient的实现，利用的是循环的方法\n",
    "def gradient(theta, X, y):\n",
    "    theta = np.matrix(theta)\n",
    "    X = np.matrix(X)\n",
    "    y = np.matrix(y).T\n",
    "    print(\"X.shape = \", X.shape)\n",
    "    print(\"y.shape = \", y.shape)\n",
    "    print(\"theta.shape = \", theta.shape)\n",
    "    \n",
    "    parameters = int(theta.ravel().shape[1])\n",
    "    grad = np.zeros(parameters)\n",
    "    \n",
    "    error = sigmoid(X @ theta.T) - y\n",
    "    print(\"error.shape = \", error.shape)\n",
    "    \n",
    "    for i in range(parameters):\n",
    "        term = X[:,i].T @ error\n",
    "        grad[i] = term / len(X)\n",
    "    \n",
    "    return grad"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "brjL4gt34jMb"
   },
   "outputs": [],
   "source": [
    "# 另外一种实现的方法\n",
    "def gradient(theta, X, y):\n",
    "    return (1 / len(X)) * X.T @ (sigmoid(X @ theta) - y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "6TEJDxaqjlIr"
   },
   "source": [
    "注意，我们实际上没有在这个函数中执行梯度下降，我们仅仅在计算一个梯度步长。在练习中，一个称为“fminunc”的Octave函数是用来优化函数来计算成本和梯度参数。由于我们使用Python，我们可以用SciPy的“optimize”命名空间来做同样的事情。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "FoLQqZg7jlIt"
   },
   "source": [
    "我们看看用我们的数据和初始参数为0的梯度下降法的结果。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "twQhHqydjlIw"
   },
   "outputs": [],
   "source": [
    "gradient(theta, X, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "M5cP5VjE-YdI"
   },
   "source": [
    "## 拟合参数"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "3ye-zhoBjlI1"
   },
   "source": [
    "现在可以用SciPy's truncated newton（TNC）实现寻找最优参数。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "0UHRpG7ijlI2"
   },
   "outputs": [],
   "source": [
    "import scipy.optimize as opt\n",
    "res = opt.minimize(fun=cost, x0=theta, args=(X, y), method='Newton-CG', jac=gradient)\n",
    "print(res)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "o-z0iOvZjlI5"
   },
   "source": [
    "让我们看看在这个结论下代价函数计算结果是什么个样子~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "d7Sc0NlzjlI6"
   },
   "outputs": [],
   "source": [
    "cost(res.x, X, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "jqbJJue1-lK6"
   },
   "source": [
    "## 用训练集预测和验证"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "F5xbkpoXjlJD"
   },
   "source": [
    "接下来，我们需要编写一个函数，用我们所学的参数theta来为数据集X输出预测。然后，我们可以使用这个函数来给我们的分类器的训练精度打分。\n",
    "逻辑回归模型的假设函数： \n",
    "\t\\\\[{{h}_{\\theta }}\\left( x \\right)=\\frac{1}{1+{{e}^{-{{\\theta }^{T}}X}}}\\\\] \n",
    "当${{h}_{\\theta }}$大于等于0.5时，预测 y=1\n",
    "\n",
    "当${{h}_{\\theta }}$小于0.5时，预测 y=0 。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "L5kMfIWijlJE"
   },
   "outputs": [],
   "source": [
    "def predict(theta, X):\n",
    "    probability = sigmoid(X @ theta)\n",
    "    return [1 if x >= 0.5 else 0 for x in probability]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "5yDGG6aXjlJI"
   },
   "outputs": [],
   "source": [
    "theta_min = res.x\n",
    "predictions = predict(theta_min, X)\n",
    "correct = [1 if ((a == 1 and b == 1) or (a == 0 and b == 0)) else 0 for (a, b) in zip(predictions, y)]\n",
    "accuracy = (sum(map(int, correct)) % len(correct))\n",
    "print ('accuracy = {0}%'.format(accuracy))\n",
    "print(classification_report(y, predictions))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "luDLSnm3jlJL"
   },
   "source": [
    "我们的逻辑回归分类器预测正确，如果一个学生被录取或没有录取，达到89%的精确度。不坏！记住，这是训练集的准确性。我们没有保持住了设置或使用交叉验证得到的真实逼近，所以这个数字有可能高于其真实值（这个话题将在以后说明）。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "EO1hfP0j-0cc"
   },
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "xIvApW8HjBZg"
   },
   "source": [
    "## 寻找决策边界\n",
    "决策便捷是一系列的$x$可以满足：$$\\frac{1}{1+e^{-\\theta^T x}}=0.5$$\n",
    "也就是说$$e^{-\\theta^T x} = 1$$\n",
    "意味着 $$\\theta^T X = 0$$\n",
    "\n",
    "这是一条直线。\n",
    "\n",
    "参考这里：http://stats.stackexchange.com/questions/93569/why-is-logistic-regression-a-linear-classifier\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "oErM-l5djBZh"
   },
   "outputs": [],
   "source": [
    "print(res.x) # this is final theta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "Bmtqg98tjBZk"
   },
   "outputs": [],
   "source": [
    "# x0 = 1\n",
    "# 定义x1的range后，根据 θ0 + x1*θ1 + x2*θ2 = 0，可得x2\n",
    "x1 = np.arange(130, step=0.1)\n",
    "\n",
    "coef = -(res.x / res.x[2])  # find the equation\n",
    "x2 = coef[0] + coef[1]*x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "zW-rRiS-jBZu"
   },
   "outputs": [],
   "source": [
    "import seaborn as sns\n",
    "sns.set(context=\"notebook\", style=\"ticks\", font_scale=1.5)\n",
    "\n",
    "sns.lmplot('exam1', 'exam2', hue='admitted', data=data, \n",
    "           size=6, \n",
    "           fit_reg=False, \n",
    "           scatter_kws={\"s\": 25}\n",
    "          )\n",
    "\n",
    "plt.plot(x1, x2, 'grey')\n",
    "plt.xlim(0, 130)\n",
    "plt.ylim(0, 130)\n",
    "plt.title('Decision Boundary')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "qgr6k1ybjlJM"
   },
   "source": [
    "## 正则化逻辑回归"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "raPYPfoXjlJN"
   },
   "source": [
    "在训练的第二部分，我们将要通过加入正则项提升逻辑回归算法。如果你对正则化有点眼生，或者喜欢这一节的方程的背景，请参考在\"exercises\"文件夹中的\"ex2.pdf\"。简而言之，正则化是成本函数中的一个术语，它使算法更倾向于“更简单”的模型（在这种情况下，模型将更小的系数）。这个理论助于减少过拟合，提高模型的泛化能力。这样，我们开始吧。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "8HRDKZsYjlJO"
   },
   "source": [
    "设想你是工厂的生产主管，你有一些芯片在两次测试中的测试结果。对于这两次测试，你想决定是否芯片要被接受或抛弃。为了帮助你做出艰难的决定，你拥有过去芯片的测试数据集，从其中你可以构建一个逻辑回归模型。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "JuiBlOvjjlJR"
   },
   "source": [
    "和第一部分很像，从数据可视化开始吧！"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "PPewmmIgjlJS"
   },
   "outputs": [],
   "source": [
    "data2 = pd.read_csv('ex2data2.txt', names=['test1', 'test2', 'accepted'])\n",
    "data2.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "PH37JgMmjlJV"
   },
   "outputs": [],
   "source": [
    "positive = data2[data2['accepted'].isin([1])]\n",
    "negative = data2[data2['accepted'].isin([0])]\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(12,8))\n",
    "ax.scatter(positive['test1'], positive['test2'], s=50, c='b', marker='o', label='Accepted')\n",
    "ax.scatter(negative['test1'], negative['test2'], s=50, c='r', marker='x', label='Rejected')\n",
    "ax.legend()\n",
    "ax.set_xlabel('Test 1 Score')\n",
    "ax.set_ylabel('Test 2 Score')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "nKZFQDlUjlJY"
   },
   "source": [
    "哇，这个数据看起来可比前一次的复杂得多。特别地，你会注意到其中没有线性决策界限，来良好的分开两类数据。一个方法是用像逻辑回归这样的线性技术来构造从原始特征的多项式中得到的特征。让我们通过创建一组多项式特征入手吧。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "0WuhZaH9PyJv"
   },
   "outputs": [],
   "source": [
    "def feature_mapping(x, y, power, as_ndarray=False):\n",
    "#     \"\"\"return mapped features as ndarray or dataframe\"\"\"\n",
    "    # data = {}\n",
    "    # # inclusive\n",
    "    # for i in np.arange(power + 1):\n",
    "    #     for p in np.arange(i + 1):\n",
    "    #         data[\"f{}{}\".format(i - p, p)] = np.power(x, i - p) * np.power(y, p)\n",
    "\n",
    "    data = {\"f{}{}\".format(i - p, p): np.power(x, i - p) * np.power(y, p)\n",
    "                for i in np.arange(power + 1)\n",
    "                for p in np.arange(i + 1)\n",
    "            }\n",
    "\n",
    "    if as_ndarray:\n",
    "        return pd.DataFrame(data).as_matrix()\n",
    "    else:\n",
    "        return pd.DataFrame(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "sISiRom-jlJZ"
   },
   "outputs": [],
   "source": [
    "degree = 5\n",
    "x1 = np.array(df.test1)\n",
    "x2 = np.array(df.test2)\n",
    "\n",
    "data = feature_mapping(x1, x2, power=6)\n",
    "print(data.shape)\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "UaNs-uSEjlJe"
   },
   "source": [
    "现在，我们需要修改第1部分的成本和梯度函数，包括正则化项。首先是成本函数："
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "YgE_bAV2jlJe"
   },
   "source": [
    "## Regularized cost（正则化代价函数）\n",
    "$$J\\left( \\theta  \\right)=\\frac{1}{m}\\sum\\limits_{i=1}^{m}{[-{{y}^{(i)}}\\log \\left( {{h}_{\\theta }}\\left( {{x}^{(i)}} \\right) \\right)-\\left( 1-{{y}^{(i)}} \\right)\\log \\left( 1-{{h}_{\\theta }}\\left( {{x}^{(i)}} \\right) \\right)]}+\\frac{\\lambda }{2m}\\sum\\limits_{j=1}^{n}{\\theta _{j}^{2}}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "qBnOVnA9RI1v"
   },
   "outputs": [],
   "source": [
    "theta = np.zeros(data.shape[1])\n",
    "X = feature_mapping(x1, x2, power=6, as_ndarray=True)\n",
    "print(X.shape)\n",
    "\n",
    "y = get_y(df)\n",
    "print(y.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "Kkw1EX4DjlJg"
   },
   "outputs": [],
   "source": [
    "def regularized_cost(theta, X, y, l=1):\n",
    "#     '''you don't penalize theta_0'''\n",
    "    theta_j1_to_n = theta[1:]\n",
    "    regularized_term = (l / (2 * len(X))) * np.power(theta_j1_to_n, 2).sum()\n",
    "\n",
    "    return cost(theta, X, y) + regularized_term\n",
    "#正则化代价函数"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "207-7syrRNE_"
   },
   "outputs": [],
   "source": [
    "regularized_cost(theta, X, y, l=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "5nkQHdeGRtF2"
   },
   "source": [
    "## Regularized gradient(正则化梯度) \n",
    "\n",
    "$$\\frac{\\partial J\\left( \\theta  \\right)}{\\partial {{\\theta }_{j}}}=\\left( \\frac{1}{m}\\sum\\limits_{i=1}^{m}{\\left( {{h}_{\\theta }}\\left( {{x}^{\\left( i \\right)}} \\right)-{{y}^{\\left( i \\right)}} \\right)} \\right)+\\frac{\\lambda }{m}{{\\theta }_{j}}\\text{ }\\text{             for  j}\\ge \\text{1}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "1BnWF1tEjlJj"
   },
   "source": [
    "请注意等式中的\"reg\" 项。还注意到另外的一个“学习率”参数。这是一种超参数，用来控制正则化项。现在我们需要添加正则化梯度函数："
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "ycKcQ0ndjlJk"
   },
   "source": [
    "如果我们要使用梯度下降法令这个代价函数最小化，因为我们未对${{\\theta }_{0}}$ 进行正则化，所以梯度下降算法将分两种情形：\n",
    "\\begin{align}\n",
    "  & Repeat\\text{ }until\\text{ }convergence\\text{ }\\!\\!\\{\\!\\!\\text{ } \\\\ \n",
    " & \\text{     }{{\\theta }_{0}}:={{\\theta }_{0}}-a\\frac{1}{m}\\sum\\limits_{i=1}^{m}{[{{h}_{\\theta }}\\left( {{x}^{(i)}} \\right)-{{y}^{(i)}}]x_{_{0}}^{(i)}} \\\\ \n",
    " & \\text{     }{{\\theta }_{j}}:={{\\theta }_{j}}-a\\frac{1}{m}\\sum\\limits_{i=1}^{m}{[{{h}_{\\theta }}\\left( {{x}^{(i)}} \\right)-{{y}^{(i)}}]x_{j}^{(i)}}+\\frac{\\lambda }{m}{{\\theta }_{j}} \\\\ \n",
    " & \\text{          }\\!\\!\\}\\!\\!\\text{ } \\\\ \n",
    " & Repeat \\\\ \n",
    "\\end{align}\n",
    "\n",
    "对上面的算法中 j=1,2,...,n 时的更新式子进行调整可得： \n",
    "${{\\theta }_{j}}:={{\\theta }_{j}}(1-a\\frac{\\lambda }{m})-a\\frac{1}{m}\\sum\\limits_{i=1}^{m}{({{h}_{\\theta }}\\left( {{x}^{(i)}} \\right)-{{y}^{(i)}})x_{j}^{(i)}}$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "cwxZ64B4jlJm"
   },
   "outputs": [],
   "source": [
    "def regularized_gradient(theta, X, y, l=1):\n",
    "#     '''still, leave theta_0 alone'''\n",
    "    theta_j1_to_n = theta[1:]\n",
    "    regularized_theta = (l / len(X)) * theta_j1_to_n\n",
    "\n",
    "    # by doing this, no offset is on theta_0\n",
    "    regularized_term = np.concatenate([np.array([0]), regularized_theta])\n",
    "\n",
    "    return gradient(theta, X, y) + regularized_term"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "m3mVnObqjlJo"
   },
   "source": [
    "就像在第一部分中做的一样，初始化变量。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "SQJGX8C5jlJq"
   },
   "outputs": [],
   "source": [
    "regularized_gradient(theta, X, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "3HDWjlJrjlJs"
   },
   "source": [
    "让我们初始学习率到一个合理值。，果有必要的话（即如果惩罚太强或不够强）,我们可以之后再折腾这个。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "Ra_Rd5sFjlJs"
   },
   "outputs": [],
   "source": [
    "learningRate = 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "c3oj6HgkjlJv"
   },
   "source": [
    "现在，让我们尝试调用新的默认为0的theta的正则化函数，以确保计算工作正常。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "haAzSuP3jlJw"
   },
   "outputs": [],
   "source": [
    "regularized_cost(theta2, X2, y2, learningRate)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "Co_DCQqPjlJ8"
   },
   "source": [
    "现在我们可以使用和第一部分相同的优化函数来计算优化后的结果。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "twdWMetWSeAd"
   },
   "source": [
    "## 拟合参数"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "rb6a5ASMSg09"
   },
   "outputs": [],
   "source": [
    "import scipy.optimize as opt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "AXbKitbDSiqS"
   },
   "outputs": [],
   "source": [
    "print('init cost = {}'.format(regularized_cost(theta, X, y)))\n",
    "\n",
    "res = opt.minimize(fun=regularized_cost, x0=theta, args=(X, y), method='Newton-CG', jac=regularized_gradient)\n",
    "res"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "UyNK7mUXjlJ_"
   },
   "source": [
    "最后，我们可以使用第1部分中的预测函数来查看我们的方案在训练数据上的准确度。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "8D7uUrYGSrW6"
   },
   "source": [
    "## 预测"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "koAbf3UhStRU"
   },
   "outputs": [],
   "source": [
    "final_theta = res.x\n",
    "y_pred = predict(final_theta, X)\n",
    "\n",
    "print(classification_report(y, y_pred))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "_yLVQLTFjlKI"
   },
   "source": [
    "虽然我们实现了这些算法，值得注意的是，我们还可以使用高级Python库像scikit-learn来解决这个问题。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "cih8ekDVjlKK"
   },
   "outputs": [],
   "source": [
    "from sklearn import linear_model#调用sklearn的线性回归包\n",
    "model = linear_model.LogisticRegression(penalty='l2', C=1.0)\n",
    "model.fit(X2, y2.ravel())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "jaT8naLtjlKM"
   },
   "outputs": [],
   "source": [
    "model.score(X2, y2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "A-SkdKuQjlKP"
   },
   "source": [
    "这个准确度和我们刚刚实现的差了好多，不过请记住这个结果可以使用默认参数下计算的结果。我们可能需要做一些参数的调整来获得和我们之前结果相同的精确度。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "ioOB1lqPUlp7"
   },
   "source": [
    "## 使用不同的 $\\lambda$ \n",
    "### 画出决策边界\n",
    "* 我们找到所有满足 $X\\times \\theta = 0$ 的x\n",
    "* instead of solving polynomial equation, just create a coridate x,y grid that is dense enough, and find all those $X\\times \\theta$ that is close enough to 0, then plot them"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "wgDzdelZjBa6"
   },
   "outputs": [],
   "source": [
    "def draw_boundary(power, l):\n",
    "#     \"\"\"\n",
    "#     power: polynomial power for mapped feature\n",
    "#     l: lambda constant\n",
    "#     \"\"\"\n",
    "    density = 1000\n",
    "    threshhold = 2 * 10**-3\n",
    "\n",
    "    final_theta = feature_mapped_logistic_regression(power, l)\n",
    "    x, y = find_decision_boundary(density, power, final_theta, threshhold)\n",
    "\n",
    "    df = pd.read_csv('ex2data2.txt', names=['test1', 'test2', 'accepted'])\n",
    "    sns.lmplot('test1', 'test2', hue='accepted', data=df, size=6, fit_reg=False, scatter_kws={\"s\": 100})\n",
    "\n",
    "    plt.scatter(x, y, c='R', s=10)\n",
    "    plt.title('Decision boundary')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "SK7PmPjkjBa-"
   },
   "outputs": [],
   "source": [
    "def feature_mapped_logistic_regression(power, l):\n",
    "#     \"\"\"for drawing purpose only.. not a well generealize logistic regression\n",
    "#     power: int\n",
    "#         raise x1, x2 to polynomial power\n",
    "#     l: int\n",
    "#         lambda constant for regularization term\n",
    "#     \"\"\"\n",
    "    df = pd.read_csv('ex2data2.txt', names=['test1', 'test2', 'accepted'])\n",
    "    x1 = np.array(df.test1)\n",
    "    x2 = np.array(df.test2)\n",
    "    y = get_y(df)\n",
    "\n",
    "    X = feature_mapping(x1, x2, power, as_ndarray=True)\n",
    "    theta = np.zeros(X.shape[1])\n",
    "\n",
    "    res = opt.minimize(fun=regularized_cost,\n",
    "                       x0=theta,\n",
    "                       args=(X, y, l),\n",
    "                       method='TNC',\n",
    "                       jac=regularized_gradient)\n",
    "    final_theta = res.x\n",
    "\n",
    "    return final_theta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "GqbGlpA7jBbB"
   },
   "outputs": [],
   "source": [
    "def find_decision_boundary(density, power, theta, threshhold):\n",
    "    t1 = np.linspace(-1, 1.5, density)\n",
    "    t2 = np.linspace(-1, 1.5, density)\n",
    "\n",
    "    cordinates = [(x, y) for x in t1 for y in t2]\n",
    "    x_cord, y_cord = zip(*cordinates)\n",
    "    mapped_cord = feature_mapping(x_cord, y_cord, power)  # this is a dataframe\n",
    "\n",
    "    inner_product = mapped_cord.as_matrix() @ theta\n",
    "\n",
    "    decision = mapped_cord[np.abs(inner_product) < threshhold]\n",
    "\n",
    "    return decision.f10, decision.f01\n",
    "#寻找决策边界函数"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "JZWGQjWljBbJ",
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "draw_boundary(power=6, l=1)#lambda=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "CLBBoc6ZjBbN"
   },
   "outputs": [],
   "source": [
    "draw_boundary(power=6, l=0)  # no regularization, over fitting，#lambda=0,没有正则化，过拟合了"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "zkOZ5XMZjBbS"
   },
   "outputs": [],
   "source": [
    "draw_boundary(power=6, l=100)  # underfitting，#lambda=100,欠拟合"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "name": "ML-Exercise2.ipynb",
   "provenance": [],
   "toc_visible": true,
   "version": "0.3.2"
  },
  "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
}