{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<link rel=\"stylesheet\" href=\"berkeley.css\">\n",
    "\n",
    "<h1 class=\"cal cal-h1\">Lecture 09: Logistic Regression – CS 189, Fall 2026</h1>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import plotly.graph_objects as go\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.linear_model import LinearRegression, LogisticRegression\n",
    "from sklearn.metrics import accuracy_score, classification_report, mean_squared_error as mse\n",
    "from scipy import stats\n",
    "import plotly.express as px\n",
    "import sklearn.linear_model as lm\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "\n",
    "import plotly.io as pio\n",
    "pio.renderers.default = \"notebook_connected\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.datasets import load_breast_cancer\n",
    "\n",
    "# Load the dataset\n",
    "data_dict = load_breast_cancer()\n",
    "data = pd.DataFrame(data_dict['data'], columns=data_dict['feature_names'])\n",
    "data['malignant'] = (data_dict['target'] == 0)  # 1 for malignant, 0 for benign\n",
    "\n",
    "# Display the first few rows\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Split the data into training and testing sets\n",
    "X = data[['mean radius']].to_numpy()\n",
    "y = data['malignant'].to_numpy()\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)\n",
    "\n",
    "print(\"Training set size:\", len(X_train))\n",
    "print(\"Testing set size:\", len(X_test))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "toy_df = pd.DataFrame({\n",
    "        \"x\": [-4, -2, -0.5, 1, 3, 5],\n",
    "        \"y\": [0, 0, 1, 0, 1, 1]\n",
    "})\n",
    "toy_df[\"str_y\"] = toy_df[\"y\"].astype(str)\n",
    "toy_df.sort_values(\"x\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = px.scatter(toy_df, x=\"x\", y=\"y\", color=\"str_y\", width=800)\n",
    "fig.update_traces(marker_size=20)\n",
    "# fig.write_image(\"toy_data.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sigmoid(z):\n",
    "    return 1/(1+np.e**-z)\n",
    "\n",
    "def ss_error_on_toy_data(theta):\n",
    "    p_hat = sigmoid(toy_df['x'] * theta)\n",
    "    return 1/2.0 * np.sum((toy_df['y'] - p_hat)**2)\n",
    "\n",
    "w_error = pd.DataFrame({\"w\": np.linspace(-10, 10, 100)})\n",
    "w_error[\"Error\"] = w_error[\"w\"].apply(ss_error_on_toy_data)\n",
    "fig = px.line(w_error, x=\"w\", y=\"Error\", width=800,\n",
    "    title=\"Error on Toy Classification Data\")\n",
    "fig.update_traces(line=dict(width=4))\n",
    "fig.show()\n",
    "# fig.write_image(\"toy_data_error.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set the initial guess as w = 0\n",
    "from scipy.optimize import minimize\n",
    "\n",
    "best_w = minimize(ss_error_on_toy_data, x0 = 0)[\"x\"][0]\n",
    "best_w"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "non_optimal_w = minimize(ss_error_on_toy_data, x0 = -5)[\"x\"][0]\n",
    "non_optimal_w"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = px.scatter(toy_df, x=\"x\", y=\"y\", color=\"str_y\", width=800)\n",
    "xs = np.linspace(-10, 10, 100)\n",
    "fig.add_trace(go.Scatter(\n",
    "    x=xs, y=sigmoid(xs * best_w), \n",
    "    mode=\"lines\", line_color=\"black\", \n",
    "    name=f\"LR Model: w = {best_w:.2f}\"))\n",
    "fig.update_traces(line=dict(width=4))\n",
    "fig.update_traces(marker_size=20)\n",
    "fig.show()\n",
    "# fig.write_image(\"toy_data_with_model_optimal.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = px.scatter(toy_df, x=\"x\", y=\"y\", color=\"str_y\", width=800)\n",
    "xs = np.linspace(-10, 10, 100)\n",
    "fig.add_trace(go.Scatter(\n",
    "    x=xs, y=sigmoid(xs * non_optimal_w), \n",
    "    mode=\"lines\", line_color=\"black\", \n",
    "    name=f\"LR Model: w = {best_w:.2f}\"))\n",
    "fig.update_traces(line=dict(width=4))\n",
    "fig.update_traces(marker_size=20)\n",
    "fig.show()\n",
    "# fig.write_image(\"toy_data_with_model_non_optimal.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p_hat_error = pd.DataFrame({\"p\": np.arange(0.001, 0.999, 0.01)}) \n",
    "p_hat_error[\"Squared Error\"] = 1/2.0* (1 - p_hat_error[\"p\"])**2\n",
    "fig = px.line(p_hat_error, x=\"p\", y=\"Squared Error\", width=800,\n",
    "        title=\"Squared Loss for One Individual when y=1\")\n",
    "fig.update_traces(line=dict(width=4))\n",
    "fig.show()\n",
    "# fig.write_image(\"squared_loss_y1.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p_hat_error[\"Neg Log Error\"] = -np.log(p_hat_error[\"p\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = px.line(p_hat_error.melt(id_vars=\"p\", value_name=\"Error\"), \n",
    "        x=\"p\", y=\"Error\", color=\"variable\", width=800,\n",
    "        title=\"Error Comparison for One Observation when y = 1\")\n",
    "fig.update_traces(line=dict(width=4))\n",
    "fig.show()\n",
    "# fig.write_image(\"loss_comparison_y1.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p_hat_error = pd.DataFrame({\"p\": np.arange(0.001, 0.999, 0.01)}) \n",
    "p_hat_error[\"Squared Error\"] = 1/2.0 * (1 - (1-p_hat_error[\"p\"]))**2\n",
    "p_hat_error[\"Neg Log Error\"] = -np.log(1 - p_hat_error[\"p\"])\n",
    "fig = px.line(p_hat_error.melt(id_vars=\"p\", value_name=\"Error\"), \n",
    "        x=\"p\", y=\"Error\", color=\"variable\", width=800,\n",
    "        title=\"Error Comparison for One Observation when y = 0\")\n",
    "fig.update_traces(line=dict(width=4))\n",
    "fig.show()\n",
    "# fig.write_image(\"loss_comparison_y0.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cross_entropy(y, p):\n",
    "    return - y * np.log(p) - (1 - y) * np.log(1 - p)\n",
    "def mean_cross_entropy_on_toy_data(w):\n",
    "    p = sigmoid(toy_df[\"x\"] * w)\n",
    "    return np.mean(cross_entropy(toy_df[\"y\"], p))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "w_error[\"Cross-Entropy\"] = w_error[\"w\"].apply(mean_cross_entropy_on_toy_data).dropna()\n",
    "fig = px.line(w_error, x=\"w\", y=\"Cross-Entropy\", width=800,\n",
    "           title=\"Cross-Entropy on Toy Classification Data\")\n",
    "fig.update_xaxes(range=[w_error[\"w\"].min(), 4])\n",
    "fig.update_traces(line=dict(width=4))\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "toy_model = lm.LogisticRegression(C=10)\n",
    "\n",
    "# We fit to two data points: (-1, 0) and (1, 1).\n",
    "toy_model.fit([[-1], [1]], [0,1])\n",
    "\n",
    "# Generate estimated probabilities across a range of x-values.\n",
    "xtest = np.linspace(-5, 5, 1000)[:, np.newaxis]\n",
    "p = toy_model.predict_proba(xtest)[:,1]\n",
    "\n",
    "fig = px.scatter(toy_df, x=\"x\", y=\"y\", \n",
    "         color=\"str_y\", symbol=\"str_y\", \n",
    "         symbol_sequence=[\"circle\", \"cross\"],\n",
    "         title=f\"LR Fit (slope = {toy_model.coef_[0][0]}, intercept = {toy_model.intercept_[0]})\",\n",
    "         render_mode=\"svg\")\n",
    "fig.update_traces(marker=dict(size=15))\n",
    "fig.update_layout(\n",
    "  xaxis_title=dict(font=dict(size=22)),\n",
    "  yaxis_title=dict(font=dict(size=22))\n",
    ")\n",
    "fig.add_scatter(x=np.ravel(xtest), y=p, mode=\"lines\", name=\"LR Model with C=10\", \n",
    "                line_color=\"black\", opacity=0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "toy_model = lm.LogisticRegression(C=1000)\n",
    "\n",
    "# We fit to two data points: (-1, 0) and (1, 1).\n",
    "toy_model.fit([[-1], [1]], [0,1])\n",
    "\n",
    "# Generate estimated probabilities across a range of x-values.\n",
    "xtest = np.linspace(-5, 5, 1000)[:, np.newaxis]\n",
    "p = toy_model.predict_proba(xtest)[:,1]\n",
    "\n",
    "fig = px.scatter(toy_df, x=\"x\", y=\"y\", \n",
    "         color=\"str_y\", symbol=\"str_y\", \n",
    "         symbol_sequence=[\"circle\", \"cross\"],\n",
    "         title=f\"LR Fit (slope = {toy_model.coef_[0][0]}, intercept = {toy_model.intercept_[0]})\",\n",
    "         render_mode=\"svg\")\n",
    "fig.update_traces(marker=dict(size=15))\n",
    "fig.update_layout(\n",
    "  xaxis_title=dict(font=dict(size=22)),\n",
    "  yaxis_title=dict(font=dict(size=22))\n",
    ")\n",
    "fig.add_scatter(x=np.ravel(xtest), y=p, mode=\"lines\", name=\"LR Model with C=1000\", \n",
    "                line_color=\"black\", opacity=0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Build a Logistic Regression Model\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = lm.LogisticRegression()\n",
    "model.fit(X_train, y_train)\n",
    "\n",
    "print(\"Slope:\", model.coef_[0][0])\n",
    "print(\"Intercept:\", model.intercept_[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, rather than predict a numeric output, we predict the *probability* of a datapoint belonging to Class 1. We do this using the `.predict_proba` method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Preview the first 10 rows\n",
    "model.predict_proba(X_train)[:10]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By default, `.predict_proba` returns a 2D array. \n",
    "\n",
    "One column contains the predicted probability that the datapoint belongs to Class 0, and the other contains the predicted probability that it belongs to Class 1 (notice that all rows sum to a total probability of 1). \n",
    "\n",
    "To check which is which, we can use the `.classes_` attribute."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model.classes_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This tells us that the first column contains the probabilities of belonging to Class 0 (benign), and the second column contains the probabilities of belonging to Class 1 (malignant). Let's grab just the probabilities of Class 1.\n",
    "\n",
    "We then apply a decision rule: Predict Class 1 if the predicted probability of belonging to Class 1 is 0.5 or higher. Otherwise, predict Class 0.\n",
    "\n",
    "- Remember that 0.5 is a common threshold, but we are not required to always use 0.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Obtain P(Y=1|x) from the output.\n",
    "p = model.predict_proba(X_train)[:, 1]\n",
    "\n",
    "# Apply decision rule: predict Class 1 if P(Y=1|x) >= 0.5.\n",
    "(p >= 0.5).astype(int)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `.predict` method of `LogisticRegression` will apply a 0.5 threshold to classify data, by default"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# .predict will automatically apply a 0.5 threshold for a logistic regression model.\n",
    "classes = model.predict(X_train).astype(int)\n",
    "\n",
    "classes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The point where the sigmoid function outputs 0.5 is the **decision boundary**. \n",
    "\n",
    "- This is the point where the model is indifferent between predicting Class 0 and Class 1.  \n",
    "\n",
    "- This is also the point where $\\theta_0 + \\theta_1 x = 0$. \n",
    "\n",
    "For this one dimensional case we can solve for the $x$ value of the decision boundary:\n",
    "\n",
    "$$\n",
    "x = - \\frac{\\theta_0}{\\theta_1} = - \\frac{\\text{intercept}}{\\text{slope}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's visualize our predictions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Convert X_train to a DataFrame for compatibility\n",
    "X_train_df = pd.DataFrame(X_train, columns=[\"mean radius\"])\n",
    "X_train_df[\"Predicted Class\"] = pd.Categorical(model.predict(X_train))\n",
    "\n",
    "test_points = pd.DataFrame({\"mean radius\": np.linspace(5, 30, 100)})\n",
    "test_points[\"Predicted Prob\"] = model.predict_proba(test_points[[\"mean radius\"]])[:, 1]\n",
    "\n",
    "fig = px.scatter(X_train_df, x=\"mean radius\", y=y_train.astype(int), color=\"Predicted Class\", opacity=0.6)\n",
    "\n",
    "# Add the logistic regression model predictions\n",
    "fig.add_trace(go.Scatter(x=test_points[\"mean radius\"], y=test_points[\"Predicted Prob\"], \n",
    "                         mode=\"lines\", name=\"Logistic Regression Model\", \n",
    "                         line_color=\"black\", line_width=5, line_dash=\"dash\"))\n",
    "fig.add_vline(x = -model.intercept_[0]/model.coef_[0][0], line_dash=\"dash\", \n",
    "              line_color=\"black\",\n",
    "              annotation_text=\"Decision Boundary\", \n",
    "              annotation_position=\"right\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Any time the predicted probability $p$ is less than 0.5, the model predicts Class 0. Otherwise, it predicts Class 1."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A decision boundary describes the **line** that splits the data into classes based on the *features*.\n",
    "\n",
    "For a model with one feature, the decision boundary is a *point* that separates the two classes. The number of dimensions of the decision boundary plot is the number of features.\n",
    "\n",
    "- We visualize this using a 1D plot to plot all data points in terms of *just* the feature.\n",
    "\n",
    "- We cannot define a decision boundary in terms of the predictions, so we remove that axis from our plot. \n",
    "\n",
    "Notice that all data points to the right of our decision boundary are classified as Class 1, while all data points to the left are classified as Class 0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = px.scatter(X_train_df, x=\"mean radius\", y=np.zeros(len(X_train_df)),\n",
    "                 symbol=\"Predicted Class\", symbol_sequence=[ \"circle-open\", \"cross\"], \n",
    "                 color=\"Predicted Class\", height=300, opacity=0.6)\n",
    "# fig.update_traces(marker_symbol='line-ns-open')\n",
    "fig.update_traces(marker_size=8)\n",
    "fig.update_layout(\n",
    "    yaxis=dict(showticklabels=False, showgrid=False, zeroline=False, title=\"\"),\n",
    ")\n",
    "\n",
    "decision_boundary =  -model.intercept_[0]/model.coef_[0][0]\n",
    "fig.add_vline(x = decision_boundary, line_dash=\"dash\", \n",
    "              line_color=\"black\",\n",
    "              annotation_text=\"Decision Boundary\", \n",
    "              annotation_position=\"top right\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Decision Boundaries\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `LogisticRegression` class of `sklearn.linear_model` behaves very similarly to the `LinearRegression` class. As before, we:\n",
    "\n",
    "1. Initialize a model object, and\n",
    "2. Fit it to our data.\n",
    "\n",
    "You find it helpful to recall the model formulation of a fitted logistic regression model with one input:\n",
    "\n",
    "$$\n",
    "\\hat{P}_{\\hat{w}}(Y=1 \\mid X) = \\sigma \\left( \\hat{w}_0 + \\hat{w}_1 X \\right) = \\frac{1}{1 + e^{-(\\hat{w}_0 + \\hat{w}_1 X)}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 🎯 Performance Metrics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's return to our data. We'll compute the **accuracy** of our `model` on this data.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def accuracy(X, Y):\n",
    "    return np.mean(model.predict(X) == Y)\n",
    "\n",
    "print(model.predict(X_train)[:5].astype(int))\n",
    "print(y_train[:5].astype(int))\n",
    "accuracy(X_train, y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model.score(X_train, y_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Important Note: `model.predict` and `model.score` use a threshold of 0.5.\n",
    "To use a different threshold, you must use `model.predict_proba` and work with probabilities directly."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Confusion matrix\n",
    "\n",
    "`scikit-learn` has an built-in `confusion_matrix` method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics import confusion_matrix\n",
    "\n",
    "# Be careful – confusion_matrix takes in y_true as the first parameter and y_pred as the second.\n",
    "# Don't mix these up!\n",
    "cm = confusion_matrix(y_train, model.predict(X_train))\n",
    "cm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = px.imshow(cm, x=[\"0\", \"1\"], y=[\"0\", \"1\"],\n",
    "          labels=dict(x=\"Predicted\", y=\"Actual\"), \n",
    "          text_auto=True, \n",
    "          color_continuous_scale=\"Blues\", \n",
    "          width=400, height=400)\n",
    "fig.update_xaxes(side=\"top\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Precision and Recall\n",
    "\n",
    "We can also compute the number of TP, TN, FP, and TN for our classifier, and then its precision and recall."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Y_hat = model.predict(X_train)\n",
    "tp = np.sum((Y_hat == 1) & (y_train == 1))\n",
    "tn = np.sum((Y_hat == 0) & (y_train == 0))\n",
    "\n",
    "fp = np.sum((Y_hat == 1) & (y_train == 0))\n",
    "fn = np.sum((Y_hat == 0) & (y_train == 1))\n",
    "\n",
    "\n",
    "print(\"True Positives: \", tp)\n",
    "print(\"True Negatives: \", tn)\n",
    "print(\"False Positives:\", fp)\n",
    "print(\"False Negatives:\", fn)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These numbers match what we see in the confusion matrix above. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Precision and Recall\n",
    "\n",
    "**Precision** -- How precise are my positive predictions? In other words, what fraction of the things the model predicted positive are actually positive?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "precision = tp / (tp + fp)\n",
    "precision"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Recall** -- What proportion of actual positives did my model recall in its predictions? In other words, what proportion of actual positive cases that were correctly identified by the model?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "recall = tp / (tp + fn)\n",
    "recall"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### True and False Positive Rates\n",
    "\n",
    "The TP, TN, FP, and TN we just calculated also allow us to compute the true and false positive rates (TPR and FPR). Recall that TPR is the same as recall."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fpr = fp/(fp + tn)\n",
    "fpr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tpr = tp/(tp + fn)\n",
    "tpr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It's important to remember that these values are all for the threshold of $T = 0.5$, which is `scikit-learn`'s default."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 🎛️ Adjusting the Classification Threshold"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before, we used a threshold of 0.5 in our decision rule: If the predicted probability was greater than 0.5 we predicted Class 1, otherwise, we predicted Class 0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_predictions(threshold = 0.5):\n",
    "    # Convert X_train to a DataFrame for compatibility\n",
    "    X_train_df = pd.DataFrame(X_train, columns=[\"mean radius\"])\n",
    "    X_train_df[\"Predicted Class\"] = model.predict_proba(X_train)[:, 1] >= threshold\n",
    "\n",
    "    X_train_df[\"Predicted Class\"] = pd.Categorical(X_train_df[\"Predicted Class\"].astype(int))\n",
    "    fig = px.scatter(X_train_df, \n",
    "            x=\"mean radius\", y=y_train.astype(int), color=\"Predicted Class\", \n",
    "            title=f\"Logistic Regression Predictions (Threshold = {threshold})\")\n",
    "    # Add the logistic regression model predictions\n",
    "    # Make the data points for the LR model curve\n",
    "    test_points = pd.DataFrame({\"mean radius\": np.linspace(5, 30, 100)})\n",
    "    test_points[\"Predicted Prob\"] = model.predict_proba(test_points)[:, 1]\n",
    "    fig.add_trace(go.Scatter(x=test_points[\"mean radius\"], y=test_points[\"Predicted Prob\"], \n",
    "                            mode=\"lines\", name=\"Logistic Regression Model\", \n",
    "                            line_color=\"black\", line_width=5, line_dash=\"dash\"))\n",
    "    decision_boundary = (-np.log(1/threshold - 1) - model.intercept_[0])/model.coef_[0][0]\n",
    "    fig.add_vline(x = decision_boundary, line_dash=\"dash\", line_color=\"black\",\n",
    "                  annotation_text=\"Decision Boundary\", annotation_position=\"right\")\n",
    "    return fig\n",
    "\n",
    "plot_predictions(0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_predictions(0.25)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When we **lower the threshold**, we require a lower predicted probability before we predict Class 1. We can think of this as us telling our model that it needs to be less \"confident\" about a data point being Class 1 before making a positive prediction. The total number of data points predicted to be Class 1 **either stays the same or increases**.\n",
    "\n",
    "The converse happens if we raise the threshold. Consider setting $T=0.75$. Now, we require a higher predicted probability before we predict Class 1. The total number of data points predicted to be Class 1 decreases."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_predictions(0.75)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Thresholds and Performance Metrics\n",
    "\n",
    "How does changing the threshold impact our performance metrics?\n",
    "\n",
    "Let's run an experiment: we'll test out several different possible thresholds. \n",
    "\n",
    "- For each threshold $T$, we'll make a decision rule where we classify any point with a predicted probability equal to or greater than $T$ as being in Class 1. \n",
    "\n",
    "- Otherwise, we'll predict Class 0. \n",
    "\n",
    "- We'll then compute the overall accuracy of the classifier when using that threshold."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define performance metrics dependent on the threshold value.\n",
    "def predict_threshold(model, X, T): \n",
    "    prob_one = model.predict_proba(X)[:, 1]\n",
    "    return (prob_one >= T).astype(int)\n",
    "\n",
    "def accuracy_threshold(X, Y, T):\n",
    "    return np.mean(predict_threshold(model, X, T) == Y)\n",
    "\n",
    "def precision_threshold(X, Y, T):\n",
    "    Y_hat = predict_threshold(model, X, T)\n",
    "    denominator = np.sum(Y_hat == 1)\n",
    "    if denominator == 0:\n",
    "        denominator = np.nan\n",
    "    return np.sum((Y_hat == 1) & (Y == 1)) / denominator\n",
    "    \n",
    "def recall_threshold(X, Y, T):\n",
    "    Y_hat = predict_threshold(model, X, T)\n",
    "    return np.sum((Y_hat == 1) & (Y == 1)) / np.sum(Y == 1)\n",
    "\n",
    "def tpr_threshold(X, Y, T): # Same as recall\n",
    "    Y_hat = predict_threshold(model, X, T)\n",
    "    return np.sum((Y_hat == 1) & (Y == 1)) / np.sum(Y == 1)\n",
    "\n",
    "def fpr_threshold(X, Y, T):\n",
    "    Y_hat = predict_threshold(model, X, T)\n",
    "    return np.sum((Y_hat == 1) & (Y == 0)) / np.sum(Y == 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "metrics = pd.DataFrame()\n",
    "metrics[\"Threshold\"] = np.linspace(0, 1, 1000)\n",
    "metrics[\"Accuracy\"] = [accuracy_threshold(X_train, y_train, t) for t in metrics[\"Threshold\"]]\n",
    "metrics[\"Precision\"] = [precision_threshold(X_train, y_train, t) for t in metrics[\"Threshold\"]]\n",
    "metrics[\"Recall\"] = [recall_threshold(X_train, y_train, t) for t in metrics[\"Threshold\"]]\n",
    "metrics.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = px.line(metrics, \n",
    "        x=\"Threshold\", y=\"Accuracy\",\n",
    "        title=\"Accuracy vs. Threshold\",\n",
    "        render_mode=\"svg\", width=600, height=600)\n",
    "\n",
    "fig.add_scatter(x=[metrics.loc[metrics['Accuracy'].idxmax(), 'Threshold']], y=[metrics.loc[metrics['Accuracy'].idxmax(), 'Accuracy']], \n",
    "                mode='markers', marker=dict(size=10, color='red'),\n",
    "                name=f\"Accuracy Max {metrics.loc[metrics['Accuracy'].idxmax(), 'Accuracy']:.5f}\",)\n",
    "fig.update_layout(\n",
    "  xaxis_title=dict(font=dict(size=22)),\n",
    "  yaxis_title=dict(font=dict(size=22))\n",
    ")\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The threshold that maximizes accuracy.\n",
    "metrics.sort_values(\"Accuracy\", ascending=False).head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It turns out that setting $T=0.5$ does not always result in the best performance! Part of the model design process for classification includes **choosing an appropriate threshold value**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Precision-Recall Curves\n",
    "In the lecture, we noted that there is a tradeoff between precision and recall.\n",
    "\n",
    "Precision $=\\frac{TP}{\\text{Positive Predictions}}=\\frac{TP}{TP+FP}$ increases as the number of false positives decreases, which occurs as the threshold is raised, since raising the threshold tends to reduce the number of positive predictions.\n",
    "\n",
    "Recall $=\\frac{TP}{\\text{Actual Class 1s}}=\\frac{TP}{TP+FN}$ increases as the number of false negatives decreases, which occurs as the threshold is lowered, since lowering the threshold tends to decrease number of negative predictions.\n",
    "\n",
    "We want to keep both precision and recall high. To do so, we'll need to strategically choose a threshold value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = px.line(metrics, \n",
    "        x=\"Threshold\", y=[\"Accuracy\", \"Precision\", \"Recall\"],\n",
    "        title=\"Performance Metrics vs. Threshold\",\n",
    "        render_mode=\"svg\", height=600, width=600)\n",
    "fig.update_layout(\n",
    "  xaxis_title=dict(font=dict(size=22)),\n",
    "  yaxis_title=dict(font=dict(size=22))\n",
    ")\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A **precision-recall** curve tests out many possible thresholds. Each point on the curve represents the precision and recall of the classifier for a *particular choice of threshold*.\n",
    "\n",
    "We choose a threshold value that keeps both precision and recall high (usually in the rightmost \"corner\" of the curve)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = px.line(metrics, x=\"Recall\", y=\"Precision\",\n",
    "        title=\"Precision vs. Recall\",\n",
    "        width=600, height=600,\n",
    "        render_mode=\"svg\")\n",
    "fig.update_layout(\n",
    "  xaxis_title=dict(font=dict(size=22)),\n",
    "  yaxis_title=dict(font=dict(size=22))\n",
    ")\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One way to balance precision and recall is to compute the **F1 score**. The F1 score is the harmonic mean of precision and recall:\n",
    "\n",
    "$$F1 = 2 \\times \\frac{\\text{precision} \\times \\text{recall}}{\\text{precision} + \\text{recall}}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "metrics[\"F1\"] = (2 * metrics[\"Precision\"] * metrics[\"Recall\"] \n",
    "                     / (metrics[\"Precision\"] + metrics[\"Recall\"]))\n",
    "ind = metrics['F1'].idxmax()\n",
    "metrics.loc[ind,:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = px.line(metrics, x=\"Threshold\", y=\"F1\",\n",
    "              title=\"Finding F1 Score Maximum\",\n",
    "              render_mode=\"svg\",\n",
    "              height=600, width=600)\n",
    "\n",
    "fig.add_scatter(x=[metrics.loc[ind, 'Threshold']], y=[metrics.loc[ind, 'F1']], \n",
    "                mode='markers', marker=dict(size=10, color='red'),\n",
    "                name=f\"F1 Max {metrics.loc[ind, 'Threshold']:.5f}\",)\n",
    "\n",
    "fig.update_layout(\n",
    "  xaxis_title=dict(font=dict(size=22)),\n",
    "  yaxis_title=dict(font=dict(size=22))\n",
    ")\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = px.line(metrics, x=\"Recall\", y=\"Precision\",\n",
    "              title=\"Precision vs. Recall\", width=600, height=600,\n",
    "              render_mode=\"svg\")\n",
    "fig.add_scatter(x=[metrics.loc[ind, 'Recall']], y=[metrics.loc[ind, 'Precision']], \n",
    "                mode='markers', marker=dict(size=10, color='red'),\n",
    "                name=f\"F1 Max {metrics.loc[ind, 'Threshold']:.5f}\")\n",
    "fig.update_layout(legend=dict(x=.5, y=.1))\n",
    "fig.update_layout(\n",
    "  xaxis_title=dict(font=dict(size=22)),\n",
    "  yaxis_title=dict(font=dict(size=22))\n",
    ")\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ROC Curves\n",
    "\n",
    "We can repeat a similar experiment for the FPR and TPR. Remember that we want to keep FPR *low* and TPR *high*. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "metrics[\"TPR\"] = [tpr_threshold(X_train, y_train, t) for t in metrics[\"Threshold\"]]\n",
    "metrics[\"FPR\"] = [fpr_threshold(X_train, y_train, t) for t in metrics[\"Threshold\"]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = px.line(metrics, x=\"Threshold\", y=[\"TPR\", \"FPR\", \"Accuracy\"],\n",
    "        render_mode=\"svg\", width=600, height=600)\n",
    "fig.update_layout(\n",
    "  xaxis_title=dict(font=dict(size=22)),\n",
    "  yaxis_title=dict(font=dict(size=22))\n",
    ")\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A **ROC curve** tests many possible decision rule thresholds. For each possible threshold, it plots the corresponding TPR and FPR of the classifier.\n",
    "\n",
    "\"ROC\" stands for \"Receiver Operating Characteristic\". It comes from the field of signal processing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = px.line(metrics, x=\"FPR\", y=\"TPR\", title=\"ROC Curve\", \n",
    "        width=600, height=600,\n",
    "        render_mode=\"svg\")\n",
    "fig.update_layout(\n",
    "  xaxis_title=dict(font=dict(size=22)),\n",
    "  yaxis_title=dict(font=dict(size=22))\n",
    ")\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " Ideally, a perfect classifier would have a FPR of 0 and TPR of 1. The area under the perfect classifier is 1. \n",
    " \n",
    " We often use the area under the ROC curve (abbreviated \"AUC\") as an indicator of model performance. The closer the AUC is to 1, the better."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = px.line(metrics, x=\"FPR\", y=\"TPR\", title=\"ROC Curve\", \n",
    "              width=600, height=600,\n",
    "              render_mode=\"svg\")\n",
    "fig.add_scatter(x=[0,0,1], y=[0,1,1], mode='lines', \n",
    "                line_dash='dash', line_color='black',\n",
    "                name=\"Perfect Classifier\")\n",
    "# move the legend inside the plot\n",
    "fig.update_layout(legend=dict(x=.5, y=.1))\n",
    "fig.update_layout(\n",
    "  xaxis_title=dict(font=dict(size=22)),\n",
    "  yaxis_title=dict(font=dict(size=22))\n",
    ")\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.9.6"
  },
  "toc": {
   "colors": {
    "hover_highlight": "#DAA520",
    "navigate_num": "#000000",
    "navigate_text": "#333333",
    "running_highlight": "#FF0000",
    "selected_highlight": "#FFD700",
    "sidebar_border": "#EEEEEE",
    "wrapper_background": "#FFFFFF"
   },
   "moveMenuLeft": true,
   "nav_menu": {
    "height": "512px",
    "width": "252px"
   },
   "navigate_menu": true,
   "number_sections": true,
   "sideBar": true,
   "threshold": 4,
   "toc_cell": false,
   "toc_section_display": "block",
   "toc_window_display": false,
   "widenNotebook": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
