{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Lecture 8: Logistic regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "![](https://www.tensorflow.org/images/colab_logo_32px.png)\n",
    "[Run in colab](https://colab.research.google.com/drive/1NSiPlX-azox9No8o7PgcF4o70bwMUH3A )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-01-10T00:19:56.779513Z",
     "iopub.status.busy": "2024-01-10T00:19:56.779124Z",
     "iopub.status.idle": "2024-01-10T00:19:56.787591Z",
     "shell.execute_reply": "2024-01-10T00:19:56.787055Z"
    },
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Last executed: 2024-01-10 00:19:56\n"
     ]
    }
   ],
   "source": [
    "import datetime\n",
    "now = datetime.datetime.now()\n",
    "print(\"Last executed: \" + now.strftime(\"%Y-%m-%d %H:%M:%S\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Estimating probabilities\n",
    "\n",
    "Estimate the probability of an instance belonging to a particular class.\n",
    "\n",
    "Can adapt linear regression algorithm for this purpose to perform *logistic regression*."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Sigmoid function\n",
    "\n",
    "Consider linear weighted sum of inputs $\\theta^{\\rm T} x$ again but then apply sigmoid function $\\sigma$:\n",
    "\n",
    "$$\\hat{p} = h_\\theta(x) = \\sigma(\\theta^{\\rm T} x), $$\n",
    "\n",
    "where\n",
    "\n",
    "$$\\sigma(t) = \\frac{1}{1+\\exp{(-t)}}. $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": [
     "inclass_exercise"
    ]
   },
   "source": [
    "#### Plot the sigmoid function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-01-10T00:19:56.827687Z",
     "iopub.status.busy": "2024-01-10T00:19:56.827194Z",
     "iopub.status.idle": "2024-01-10T00:19:57.238220Z",
     "shell.execute_reply": "2024-01-10T00:19:57.237523Z"
    },
    "tags": [
     "solution",
     "inclass_exercise"
    ]
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "%matplotlib inline\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-01-10T00:19:57.241818Z",
     "iopub.status.busy": "2024-01-10T00:19:57.241170Z",
     "iopub.status.idle": "2024-01-10T00:19:57.489933Z",
     "shell.execute_reply": "2024-01-10T00:19:57.489276Z"
    },
    "tags": [
     "solution",
     "inclass_exercise"
    ]
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 27976 characters]",
      "text/plain": [
       "<Figure size 800x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "t = np.linspace(-10, 10, 100)\n",
    "sig = 1 / (1 + np.exp(-t))\n",
    "plt.figure(figsize=(8, 4))\n",
    "plt.plot([-10, 10], [0, 0], \"k-\")\n",
    "plt.plot([-10, 10], [0.5, 0.5], \"k:\")\n",
    "plt.plot([-10, 10], [1, 1], \"k:\")\n",
    "plt.plot([0, 0], [-1.1, 1.1], \"k-\")\n",
    "plt.plot(t, sig, \"b-\", linewidth=2, label=r\"$\\sigma(t) = \\frac{1}{1 + e^{-t}}$\")\n",
    "plt.xlabel(\"t\")\n",
    "plt.legend(loc=\"upper left\", fontsize=20)\n",
    "plt.axis([-10, 10, -0.1, 1.1]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Predictions\n",
    "\n",
    "Can then make class predictions depending on whether the predicted probability $\\hat{p}$ is greater than 0.5, i.e.\n",
    "\n",
    "$$ \\hat{y} = \\biggl \\{\n",
    "\\begin{split}\n",
    "0,\\ \\text{if}\\ \\hat{p} < 0.5\\\\\n",
    "1,\\ \\text{if}\\ \\hat{p} \\geq  0.5\n",
    "\\end{split}\n",
    ",\n",
    "$$\n",
    "\n",
    "where we recall \n",
    "$\\hat{p} = h_\\theta(x) = \\sigma(\\theta^{\\rm T} x) \\quad \\text{and} \\quad \\sigma(t) = \\frac{1}{1+\\exp{(-t)}}.$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Note that $\\sigma(t) < 0.5$ when $t<0$ and $\\sigma(t) \\geq 0.5$ when $t \\geq 0$.  \n",
    "\n",
    "That is, logistic regression predicts model 1 when $\\theta^{\\rm T} x$ is positive, and model 0 when it is negative.  \n",
    "\n",
    "The decision boundary is defined by $\\theta^{\\rm T} x=0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Cost functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider the cost function:\n",
    "\n",
    "$$ C(\\theta) = \\biggl \\{\n",
    "\\begin{split}\n",
    "-\\log(\\hat{p}),\\ \\text{if}\\ y=1\\\\\n",
    "-\\log(1 - \\hat{p}),\\ \\text{if}\\ y=0\n",
    "\\end{split}\n",
    ".\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Exercise: Plot the cost function for $y=1$ as a function of $\\hat{p}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "inclass_exercise"
    ]
   },
   "source": [
    "Consider the cost function:\n",
    "\n",
    "$$ C(\\theta) = \\biggl \\{\n",
    "\\begin{split}\n",
    "-\\log(\\hat{p}),\\ \\text{if}\\ y=1\\\\\n",
    "-\\log(1 - \\hat{p}),\\ \\text{if}\\ y=0\n",
    "\\end{split}\n",
    ".\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-01-10T00:19:57.493545Z",
     "iopub.status.busy": "2024-01-10T00:19:57.492986Z",
     "iopub.status.idle": "2024-01-10T00:19:57.715480Z",
     "shell.execute_reply": "2024-01-10T00:19:57.714807Z"
    },
    "tags": [
     "solution",
     "inclass_exercise"
    ]
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 28640 characters]",
      "text/plain": [
       "<Figure size 800x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ph = np.linspace(0.01, 0.99, 100)\n",
    "cost_one = -np.log(ph)\n",
    "plt.figure(figsize=(8, 4))\n",
    "plt.plot([0, 0], [-1.1, 1.1], \"k-\")\n",
    "plt.plot(ph, cost_one, \"b-\", linewidth=2, label=r\"$-\\log{(\\hat{p})}$\")\n",
    "plt.xlabel(\"$\\hat{p}$\")\n",
    "plt.legend(loc=\"upper right\", fontsize=20)\n",
    "plt.axis([0, 1, 0, 4]);\n",
    "plt.title('Cost for $y=1$');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": [
     "inclass_exercise"
    ]
   },
   "source": [
    "What can you say intuitively about the cost function at the edges of the domain?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": [
     "solution",
     "inclass_exercise"
    ]
   },
   "source": [
    "- For $\\hat{p}=1$, $C(\\hat{p}) = 0$.\n",
    "- For $\\hat{p}=0$, $C(\\hat{p}) \\rightarrow \\infty$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": [
     "inclass_exercise"
    ]
   },
   "source": [
    "### Exercise: Plot the cost function for $y=0$ as a function of $\\hat{p}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-01-10T00:19:57.718928Z",
     "iopub.status.busy": "2024-01-10T00:19:57.718470Z",
     "iopub.status.idle": "2024-01-10T00:19:57.934523Z",
     "shell.execute_reply": "2024-01-10T00:19:57.933872Z"
    },
    "tags": [
     "inclass_exercise",
     "solution"
    ]
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 28208 characters]",
      "text/plain": [
       "<Figure size 800x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "cost_zero = -np.log(1-ph)\n",
    "plt.figure(figsize=(8, 4))\n",
    "plt.plot([0, 0], [-1.1, 1.1], \"k-\")\n",
    "plt.plot(ph, cost_zero, \"b-\", linewidth=2, label=r\"$-\\log{(1-\\hat{p})}$\")\n",
    "plt.xlabel(\"$\\hat{p}$\")\n",
    "plt.legend(loc=\"upper left\", fontsize=20)\n",
    "plt.axis([0, 1, 0, 4]);\n",
    "plt.title('Cost for $y=0$');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": [
     "inclass_exercise"
    ]
   },
   "source": [
    "What can you say intuitively about the cost function at the edges of the domain?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": [
     "solution",
     "inclass_exercise"
    ]
   },
   "source": [
    "- For $\\hat{p}=0$, $C(\\hat{p}) = 0$.\n",
    "- For $\\hat{p}=1$, $C(\\hat{p}) \\rightarrow \\infty$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Log-loss function for logistic regression\n",
    "\n",
    "Cost function can be written by the single expression\n",
    "\n",
    "$$\n",
    "C(\\theta) = -\\frac{1}{m} \\sum_{i=1}^m \n",
    "\\left [ \n",
    "y^{(i)} \\log(\\hat{p}^{(i)})\n",
    "+\n",
    "(1 - y^{(i)}) \\log(1 - \\hat{p}^{(i)})\n",
    "\\right],\n",
    "$$\n",
    "\n",
    "since $y^{(i)}$ is always 0 or 1 and we thus recover the separate cases considered above."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "#### Aside: statistical interpretation\n",
    "\n",
    "Interpret $p$ as probability of target $y$:\n",
    "\n",
    "$$\n",
    "{\\rm P}(y \\: \\vert \\: p ) = p^{y} (1-p)^{1-y}\n",
    "$$\n",
    "\n",
    "$$\n",
    "{\\rm log}{\\rm P}(y \\: \\vert \\: p, x ) = y {\\rm log}(p) + (1-y) {\\rm log}(1-p)\n",
    "$$\n",
    "\n",
    "See [MacKay](http://www.inference.org.uk/itila/book.html) [Chapter 41] for further details."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Minimising the cost function\n",
    "\n",
    "No closed form solution like linear regression.\n",
    "\n",
    "But since the cost function is convex guaranteed to find global minimum by gradient descent."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Derivative of the cost function\n",
    "\n",
    "\\begin{align*}\n",
    "\\frac{\\partial C}{\\partial \\theta} \n",
    "&=\n",
    "\\frac{1}{m} \\sum_{i=1}^m \n",
    "\\left[ \\sigma\\left(\\theta^{\\rm T} x^{(i)} \\right) - y^{(i)} \\right] \n",
    "x^{(i)} \\\\\n",
    "&=\n",
    "\\frac{1}{m} \n",
    "X^{\\rm T}\n",
    "\\left[ \\sigma\\left(X \\theta \\right) - y \\right] \\\\\n",
    "&=\n",
    "\\frac{1}{m} \n",
    "X^{\\rm T}\n",
    "\\left[ h_\\theta\\left(X\\right) - y \\right]\n",
    "\\end{align*}\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Similarity with linear regression\n",
    "\n",
    "Identical to linear regression (up to factor of 2 depending on conventions adopted) but with a different prediction function:\n",
    "\n",
    "$$h_\\theta(x) = \\sigma(\\theta^{\\rm T} x), $$\n",
    "\n",
    "instead of\n",
    "\n",
    "$$h_\\theta(x) = \\theta^{\\rm T} x. $$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Example of logistic regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider [Iris flower data](https://en.wikipedia.org/wiki/Iris_flower_data_set) again."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-01-10T00:19:57.939374Z",
     "iopub.status.busy": "2024-01-10T00:19:57.937962Z",
     "iopub.status.idle": "2024-01-10T00:19:58.295225Z",
     "shell.execute_reply": "2024-01-10T00:19:58.294555Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['data',\n",
       " 'target',\n",
       " 'frame',\n",
       " 'target_names',\n",
       " 'DESCR',\n",
       " 'feature_names',\n",
       " 'filename',\n",
       " 'data_module']"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sklearn import datasets\n",
    "iris = datasets.load_iris()\n",
    "list(iris.keys())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-01-10T00:19:58.299485Z",
     "iopub.status.busy": "2024-01-10T00:19:58.298085Z",
     "iopub.status.idle": "2024-01-10T00:19:58.303784Z",
     "shell.execute_reply": "2024-01-10T00:19:58.303189Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "%matplotlib inline\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Train model\n",
    "\n",
    "Use petal width to classify whether Virginica or not."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-01-10T00:19:58.307016Z",
     "iopub.status.busy": "2024-01-10T00:19:58.306450Z",
     "iopub.status.idle": "2024-01-10T00:19:58.311305Z",
     "shell.execute_reply": "2024-01-10T00:19:58.310697Z"
    }
   },
   "outputs": [],
   "source": [
    "# Set up training data\n",
    "X_1d = iris[\"data\"][:, 3:]  # petal width\n",
    "y = (iris[\"target\"] == 2).astype(int)  # 1 if Iris-Virginica, else 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-01-10T00:19:58.315162Z",
     "iopub.status.busy": "2024-01-10T00:19:58.313802Z",
     "iopub.status.idle": "2024-01-10T00:19:58.370892Z",
     "shell.execute_reply": "2024-01-10T00:19:58.370148Z"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LogisticRegression\n",
    "log_reg = LogisticRegression(random_state=42)\n",
    "log_reg.fit(X_1d, y);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that Scikit-Learn automatically adds $\\ell_2$ regularizer to cost function."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Prediction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-01-10T00:19:58.375263Z",
     "iopub.status.busy": "2024-01-10T00:19:58.373897Z",
     "iopub.status.idle": "2024-01-10T00:19:58.555888Z",
     "shell.execute_reply": "2024-01-10T00:19:58.555274Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 43048 characters]",
      "text/plain": [
       "<Figure size 800x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "X_1d_new = np.linspace(0, 3, 1000).reshape(-1, 1)\n",
    "y_1d_proba = log_reg.predict_proba(X_1d_new)\n",
    "\n",
    "plt.figure(figsize=(8,4))\n",
    "plt.plot(X_1d_new, y_1d_proba[:, 1], \"g-\", linewidth=2, label=\"Iris-Virginica\")\n",
    "plt.plot(X_1d_new, y_1d_proba[:, 0], \"b--\", linewidth=2, label=\"Not Iris-Virginica\")\n",
    "plt.xlabel(\"Petal width (cm)\", fontsize=14)\n",
    "plt.ylabel(\"Probability\", fontsize=14)\n",
    "plt.legend(loc=\"center left\", fontsize=14);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Decision boundary"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Recall the decision boundary is given by $\\hat{p}=0.5$ or, equivalently, $\\theta^{\\rm T} x=0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-01-10T00:19:58.559105Z",
     "iopub.status.busy": "2024-01-10T00:19:58.558652Z",
     "iopub.status.idle": "2024-01-10T00:19:58.565651Z",
     "shell.execute_reply": "2024-01-10T00:19:58.565067Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1.66066066])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "decision_boundary = X_1d_new[y_1d_proba[:, 1] >= 0.5][0]\n",
    "decision_boundary"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Updating plot with decision boundary and training data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-01-10T00:19:58.568895Z",
     "iopub.status.busy": "2024-01-10T00:19:58.568352Z",
     "iopub.status.idle": "2024-01-10T00:19:58.801548Z",
     "shell.execute_reply": "2024-01-10T00:19:58.800895Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 54460 characters]",
      "text/plain": [
       "<Figure size 800x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(8, 4))\n",
    "plt.plot(X_1d[y==0], y[y==0], \"bs\")\n",
    "plt.plot(X_1d[y==1], y[y==1], \"g^\")\n",
    "plt.plot([decision_boundary, decision_boundary], [-1, 2], \"k:\", linewidth=2)\n",
    "plt.plot(X_1d_new, y_1d_proba[:, 1], \"g-\", linewidth=2, label=\"Iris-Virginica\")\n",
    "plt.plot(X_1d_new, y_1d_proba[:, 0], \"b--\", linewidth=2, label=\"Not Iris-Virginica\")\n",
    "plt.text(decision_boundary+0.02, 0.15, \"Decision  boundary\", fontsize=14, color=\"k\", ha=\"center\")\n",
    "plt.arrow(decision_boundary[0], 0.08, -0.3, 0, head_width=0.05, head_length=0.1, fc='b', ec='b')\n",
    "plt.arrow(decision_boundary[0], 0.92, 0.3, 0, head_width=0.05, head_length=0.1, fc='g', ec='g')\n",
    "plt.xlabel(\"Petal width (cm)\", fontsize=14)\n",
    "plt.ylabel(\"Probability\", fontsize=14)\n",
    "plt.legend(loc=\"center left\", fontsize=14)\n",
    "plt.axis([0, 3, -0.02, 1.02]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Predictions depend on what side of decision boundary fall."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-01-10T00:19:58.804779Z",
     "iopub.status.busy": "2024-01-10T00:19:58.804533Z",
     "iopub.status.idle": "2024-01-10T00:19:58.811352Z",
     "shell.execute_reply": "2024-01-10T00:19:58.810755Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1, 0])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "log_reg.predict([[1.7], [1.5]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercises:** *You can now complete Exercise 1 in the exercises associated with this lecture.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Extending to two features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-01-10T00:19:58.814372Z",
     "iopub.status.busy": "2024-01-10T00:19:58.813922Z",
     "iopub.status.idle": "2024-01-10T00:19:58.827801Z",
     "shell.execute_reply": "2024-01-10T00:19:58.827253Z"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LogisticRegression\n",
    "\n",
    "X = iris[\"data\"][:, (2, 3)]  # petal length, petal width\n",
    "y = (iris[\"target\"] == 2).astype(int)  # 1 if Iris-Virginica, else 0\n",
    " \n",
    "C = 1000 # inverse regularization (smaller values correspond to stronger regularization)\n",
    "log_reg = LogisticRegression(C=C, random_state=42)\n",
    "log_reg.fit(X, y)\n",
    "\n",
    "x0, x1 = np.meshgrid(\n",
    "        np.linspace(2.9, 7, 500).reshape(-1, 1),\n",
    "        np.linspace(0.8, 2.7, 200).reshape(-1, 1),\n",
    "    )\n",
    "X_new = np.c_[x0.ravel(), x1.ravel()]\n",
    "\n",
    "y_proba = log_reg.predict_proba(X_new)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-01-10T00:19:58.830821Z",
     "iopub.status.busy": "2024-01-10T00:19:58.830145Z",
     "iopub.status.idle": "2024-01-10T00:19:59.072590Z",
     "shell.execute_reply": "2024-01-10T00:19:59.071979Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 129460 characters]",
      "text/plain": [
       "<Figure size 1000x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(10, 5))\n",
    "plt.plot(X[y==0, 0], X[y==0, 1], \"bs\")\n",
    "plt.plot(X[y==1, 0], X[y==1, 1], \"g^\")\n",
    "\n",
    "zz = y_proba[:, 1].reshape(x0.shape)\n",
    "contour = plt.contour(x0, x1, zz, cmap=plt.cm.brg)\n",
    "\n",
    "# Solve theta^T x = 0 to determine boundary\n",
    "left_right = np.array([2.9, 7]) \n",
    "boundary = -(log_reg.coef_[0][0] * left_right + log_reg.intercept_[0]) / log_reg.coef_[0][1]\n",
    "\n",
    "plt.clabel(contour, inline=1, fontsize=12)\n",
    "plt.plot(left_right, boundary, \"k--\", linewidth=3)\n",
    "plt.text(3.5, 1.5, \"Not Iris-Virginica\", fontsize=14, color=\"b\", ha=\"center\")\n",
    "plt.text(6.5, 2.3, \"Iris-Virginica\", fontsize=14, color=\"g\", ha=\"center\")\n",
    "plt.xlabel(\"Petal length\", fontsize=14)\n",
    "plt.ylabel(\"Petal width\", fontsize=14)\n",
    "plt.axis([2.9, 7, 0.8, 2.7]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercises:** *You can now complete Exercise 2 in the exercises associated with this lecture.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Softmax regression\n",
    "\n",
    "Can generalise logistic regression to classify multiple classes."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Softmax score\n",
    "\n",
    "Consider the softmax score function for class $k$:\n",
    "\n",
    "$$s_k(x) = \\left(\\theta^{(k)}\\right)^{\\rm T} x .$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "**Important note**: each class $k$ has its own score and set of parameters $\\theta^{(k)}$, for $K$ classes (i.e. $k=1,\\ldots,K$).\n",
    "\n",
    "Define:\n",
    "- Parameter matrix: $\\Theta_{K \\times n} = [ \\theta^{(1)},\\ \\theta^{(2)},\\ ...,\\ \\theta^{(K)}]^{\\rm T}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Softmax function\n",
    "\n",
    "Predictions are then given by the softmax function $\\sigma_k(s(x))$ for each $k$:\n",
    "\n",
    "$$ \n",
    "\\hat{p}_k = \\sigma_k(s(x)) = \\frac{\\exp\\left(s_k(x)\\right)}{\\sum_{k^\\prime=1}^K \\exp\\left(s_{k^\\prime}(x)\\right)}\n",
    ".\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Normalised such that\n",
    "- $\\sum_k \\hat{p}_k = 1$\n",
    "- $0 \\leq \\hat{p}_k \\leq 1$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Predictions\n",
    "\n",
    "Can then make class predictions based on which class has the highest predicted probability, i.e.\n",
    "\n",
    "$$  \n",
    "\\hat{y}\n",
    "= \\arg\\max_k \\hat{p}_k \n",
    "= \\arg\\max_k s_k(x)\n",
    "= \\arg\\max_k \\left(\\theta^{(k)}\\right)^{\\rm T} x\n",
    ",\n",
    "$$\n",
    "\n",
    "where we recall \n",
    "$\\hat{p}_k = \\sigma_k(s(x)) = \\frac{\\exp\\left(s_k(x)\\right)}{\\sum_{k^\\prime=1}^K \\exp\\left(s_{k^\\prime}(x)\\right)}\n",
    "\\quad \\text{and} \\quad \n",
    "s_k(x) = \\left(\\theta^{(k)}\\right)^{\\rm T} x.$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Cost function\n",
    "\n",
    "Generalization of the logistic regression cost function is given by the *cross-entropy* (measure of similarity of probability distributions):\n",
    "\n",
    "$$\n",
    "C(\\Theta) = -\\frac{1}{m} \\sum_{i=1}^m \\sum_{k=1}^K y_k^{(i)} \\log \\left(\\hat{p}_k^{(i)}\\right)\n",
    ".\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "For the case $K=2$, the cost functions reduces to the standard cost function for logistic regression."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Minimising the cost function\n",
    "\n",
    "Can solve by gradient descent.\n",
    "\n",
    "Derivative of cost function given by\n",
    "\n",
    "$$\n",
    "\\frac{\\partial C}{\\partial \\theta^{(k)}}\n",
    "= \\frac{1}{m} \\sum_{i=1}^m \\left(\\hat{p}_k^{(i)} - y_k^{(i)} \\right) x^{(i)} .\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Example of softmax regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-01-10T00:19:59.076992Z",
     "iopub.status.busy": "2024-01-10T00:19:59.076312Z",
     "iopub.status.idle": "2024-01-10T00:19:59.088101Z",
     "shell.execute_reply": "2024-01-10T00:19:59.087548Z"
    }
   },
   "outputs": [],
   "source": [
    "X = iris[\"data\"][:, (2, 3)]  # petal length, petal width\n",
    "y = iris[\"target\"] # consider all three target classes\n",
    "\n",
    "C = 10\n",
    "softmax_reg = LogisticRegression(multi_class=\"multinomial\", solver=\"lbfgs\", C=C, random_state=42)\n",
    "softmax_reg.fit(X, y);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-01-10T00:19:59.091075Z",
     "iopub.status.busy": "2024-01-10T00:19:59.090368Z",
     "iopub.status.idle": "2024-01-10T00:19:59.115277Z",
     "shell.execute_reply": "2024-01-10T00:19:59.114655Z"
    }
   },
   "outputs": [],
   "source": [
    "x0, x1 = np.meshgrid(\n",
    "        np.linspace(0, 8, 500).reshape(-1, 1),\n",
    "        np.linspace(0, 3.5, 200).reshape(-1, 1),\n",
    "    )\n",
    "X_new = np.c_[x0.ravel(), x1.ravel()]\n",
    "\n",
    "y_proba = softmax_reg.predict_proba(X_new)\n",
    "y_predict = softmax_reg.predict(X_new)\n",
    "\n",
    "# Select contours to plot\n",
    "# zz1 = y_proba[:, 0].reshape(x0.shape)\n",
    "zz1 = y_proba[:, 1].reshape(x0.shape)\n",
    "# zz1 = y_proba[:, 2].reshape(x0.shape)\n",
    "zz = y_predict.reshape(x0.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-01-10T00:19:59.118315Z",
     "iopub.status.busy": "2024-01-10T00:19:59.117863Z",
     "iopub.status.idle": "2024-01-10T00:19:59.387990Z",
     "shell.execute_reply": "2024-01-10T00:19:59.387264Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 132636 characters]",
      "text/plain": [
       "<Figure size 1000x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(10, 5))\n",
    "plt.plot(X[y==2, 0], X[y==2, 1], \"g^\", label=\"Iris-Virginica\")\n",
    "plt.plot(X[y==1, 0], X[y==1, 1], \"bs\", label=\"Iris-Versicolor\")\n",
    "plt.plot(X[y==0, 0], X[y==0, 1], \"yo\", label=\"Iris-Setosa\")\n",
    "\n",
    "from matplotlib.colors import ListedColormap\n",
    "custom_cmap = ListedColormap(['#fafab0','#9898ff','#a0faa0'])\n",
    "\n",
    "plt.contourf(x0, x1, zz, cmap=custom_cmap)\n",
    "contour = plt.contour(x0, x1, zz1, cmap=plt.cm.brg)\n",
    "plt.clabel(contour, inline=1, fontsize=12)\n",
    "plt.xlabel(\"Petal length\", fontsize=14)\n",
    "plt.ylabel(\"Petal width\", fontsize=14)\n",
    "plt.legend(loc=\"center left\", fontsize=14)\n",
    "plt.axis([0, 7, 0, 3.5]);"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.13.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
