{ "cells": [ { "cell_type": "markdown", "id": "b404adc3-5c0b-4b24-a8ec-0b7267820b18", "metadata": {}, "source": [ "# Classification: Logistic Regression" ] }, { "cell_type": "code", "execution_count": null, "id": "8484932e-b20c-43da-a289-fdda602320b5", "metadata": {}, "outputs": [], "source": [ "# initialize packages\n", "%matplotlib ipympl\n", "import matplotlib.pyplot as plt\n", "from mpl_toolkits.mplot3d import Axes3D\n", "import numpy as np\n", "import ipywidgets as widgets\n", "\n", "# load data\n", "from sklearn.datasets import load_iris\n", "iris = load_iris(as_frame=True)\n", "\n", "# initialize data\n", "X_iris = iris.data[[\"petal length (cm)\", \"petal width (cm)\"]].values\n", "y_iris_logical = iris.target_names[iris.target] == 'virginica'\n", "y_iris = y_iris_logical.astype(int)" ] }, { "cell_type": "markdown", "id": "de79f214-413a-4b5f-9b6e-eedd76ee200b", "metadata": {}, "source": [ "# Activity 1: classification boundary" ] }, { "cell_type": "code", "execution_count": null, "id": "68110e18-a64e-4788-9fb7-c2da4787e6f4", "metadata": {}, "outputs": [], "source": [ "# initialize beta ranges\n", "beta0_min, beta0_max = -80, 0\n", "beta1_min, beta1_max = 2, 10\n", "beta2_min, beta2_max = 5, 15\n", "\n", "# initialize x ranges\n", "x1_min, x1_max = 0, 8\n", "x2_min, x2_max = 0, 3\n", "\n", "# set up sliders\n", "slider10 = widgets.FloatSlider(min=beta0_min, max=beta0_max, value=(beta0_min+beta0_max)/2, description=r'$\\beta_0:$')\n", "slider11 = widgets.FloatSlider(min=beta1_min, max=beta1_max, value=(beta1_min+beta1_max)/2, description=r'$\\beta_1:$')\n", "slider12 = widgets.FloatSlider(min=beta2_min, max=beta2_max, value=(beta2_min+beta2_max)/2, description=r'$\\beta_2:$')\n", "\n", "# set up plot figure\n", "fig1, ax1 = plt.subplots(figsize=(6,6))\n", "fig1.canvas.header_visible = False\n", "fig1.canvas.footer_visible = False\n", "fig1.canvas.toolbar_visible = False\n", "\n", "ax1.set_xlabel(\"petal length (cm)\"); ax1.set_ylabel(\"petal width (cm)\")\n", "ax1.set_xlim((x1_min,x1_max)); ax1.set_ylim((x2_min,x2_max))\n", "\n", "# Initialize line object\n", "line1, = ax1.plot([], [], color='red', linestyle='dashed')\n", "\n", "def plot_data_1(X, y_logical):\n", " ax1.scatter(X[y_logical, 0], X[y_logical, 1], color='orange', label='virginica')\n", " ax1.scatter(X[~y_logical, 0], X[~y_logical, 1], color='blue', label='not virginica')\n", " ax1.legend(loc='upper left')\n", "\n", "def plot_boundary_1(beta0, beta1, beta2):\n", " x1line = np.array([x1_min,x1_max])\n", " x2line = -(beta0+beta1*x1line)/beta2\n", " line1.set_data(x1line, x2line)\n", "\n", "plot_data_1(X_iris, y_iris_logical)\n", "\n", "# create initial boundary and monitor sliders\n", "plot_boundary_1(slider10.value, slider11.value, slider12.value)\n", "widgets.interactive(plot_boundary_1, beta0=slider10, beta1=slider11, beta2=slider12)" ] }, { "cell_type": "markdown", "id": "3f87d3bb-0979-4b35-b626-80a8c7a40a0e", "metadata": {}, "source": [ "# Activity 2: probability function" ] }, { "cell_type": "code", "execution_count": null, "id": "ddc5a66a-cfbc-43f4-bf4f-88388dec86f0", "metadata": {}, "outputs": [], "source": [ "# initialize beta ranges\n", "beta0_min, beta0_max = -90, 0\n", "beta1_min, beta1_max = 2, 10\n", "beta2_min, beta2_max = 5, 15\n", "\n", "# initialize x ranges\n", "x1_min, x1_max = -1, 8\n", "x2_min, x2_max = -1, 4\n", "\n", "# set up sliders\n", "slider20 = widgets.FloatSlider(min=beta0_min, max=beta0_max, value=(beta0_min+beta0_max)/2, description=r'$\\beta_0:$')\n", "slider21 = widgets.FloatSlider(min=beta1_min, max=beta1_max, value=(beta1_min+beta1_max)/2, description=r'$\\beta_1:$')\n", "slider22 = widgets.FloatSlider(min=beta2_min, max=beta2_max, value=(beta2_min+beta2_max)/2, description=r'$\\beta_2:$')\n", "\n", "# set up plot figure\n", "fig2 = plt.figure(figsize=(6,6))\n", "fig2.canvas.header_visible = False\n", "fig2.canvas.footer_visible = False\n", "fig2.canvas.toolbar_visible = False\n", "ax2 = plt.axes(projection='3d')\n", "\n", "# create grid of x values\n", "x1s = np.linspace(x1_min, x1_max, 40)\n", "x2s = np.linspace(x2_min, x2_max, 40)\n", "x1sv, x2sv = np.meshgrid(x1s, x2s)\n", "\n", "def plot_pdf_2(beta0, beta1, beta2):\n", " # create grid of heights\n", " def sigmoid(x):\n", " return 1 / (1 + np.exp(-x))\n", " zsv = sigmoid(beta0+beta1*x1sv+beta2*x2sv)\n", " ax2.plot_surface(x1sv, x2sv, zsv, alpha=0.5, color='blue')\n", " ax2.set_xlabel('petal length (cm)')\n", " ax2.set_ylabel('petal width (cm)')\n", " ax2.set_zlabel('probability')\n", " ax2.set_xlim3d((x1_min,x1_max)); ax2.set_ylim3d((x2_min,x2_max)); ax2.set_zlim3d((0,1))\n", "\n", "def plot_data_2(X, y_logical):\n", " ax2.scatter(X[y_logical].T[0], X[y_logical].T[1], color='orange')\n", " ax2.scatter(X[~y_logical].T[0], X[~y_logical].T[1], color='blue')\n", " \n", "def plot_boundary_2(beta0, beta1, beta2): \n", " numpoints=40\n", " xline = np.linspace(x1_min,x1_max,numpoints)\n", " yline = -(beta0+beta1*xline)/beta2\n", " def clip_x(x):\n", " return(x) if x1_min <= x <= x1_max else np.nan\n", " def clip_y(y):\n", " return(y) if x2_min <= y <= x2_max else np.nan\n", " xline = np.array(list(map(clip_x, xline)))\n", " yline = np.array(list(map(clip_y, yline)))\n", " zlinezero = np.full((numpoints,),0)\n", " zlinehalf = np.full((numpoints,),0.5)\n", " ax2.plot(xline,yline, zlinezero, color='red', linestyle='dashed')\n", " ax2.plot(xline,yline, zlinehalf, color='red') \n", "\n", "def update_pdf_plot_2(beta0, beta1, beta2):\n", " ax2.clear()\n", " plot_pdf_2(beta0, beta1, beta2) \n", " plot_data_2(X_iris, y_iris_logical)\n", " plot_boundary_2(beta0, beta1, beta2)\n", " fig2.canvas.draw()\n", "\n", "# create initial plot and monitor sliders\n", "update_pdf_plot_2(slider20.value, slider21.value, slider22.value)\n", "widgets.interactive(update_pdf_plot_2, beta0=slider20, beta1=slider21, beta2=slider22)" ] }, { "cell_type": "markdown", "id": "e2291b57-d0f4-468c-b269-52e0415dd96f", "metadata": {}, "source": [ "# Likelihood function" ] }, { "cell_type": "code", "execution_count": null, "id": "abc61111-8a85-440e-b4ec-2f537a159427", "metadata": {}, "outputs": [], "source": [ "def likelihood(beta0, beta1, beta2, X, y):\n", " def sigmoid(x):\n", " return 1/(1+np.exp(-x))\n", " p = sigmoid(beta0+beta1*X.T[0]+beta2*X.T[1])\n", " return np.prod((p**y)*((1-p)**(1-y)))\n", "\n", "def vectorized_likelihood(beta0, beta1, beta2, X, y):\n", " def sigmoid(x):\n", " return 1/(1+np.exp(-x))\n", " p = sigmoid(beta0[:,:,np.newaxis]+beta1[:,:,np.newaxis]*X.T[0]+beta2[:,:,np.newaxis]*X.T[1])\n", " return np.prod((p**y)*((1-p)**(1-y)), axis=2)\n", "\n", "# initialize beta ranges\n", "beta0_min, beta0_max = -90, 0\n", "beta1_min, beta1_max = 2, 10\n", "beta2_min, beta2_max = 5, 15\n", "\n", "# set up sliders\n", "slider30 = widgets.FloatSlider(min=beta0_min, max=beta0_max, value=(beta0_min+beta0_max)/2, description=r'$\\beta_0:$')\n", "slider31 = widgets.FloatSlider(min=beta1_min, max=beta1_max, value=(beta1_min+beta1_max)/2, description=r'$\\beta_1:$')\n", "slider32 = widgets.FloatSlider(min=beta2_min, max=beta2_max, value=(beta2_min+beta2_max)/2, description=r'$\\beta_2:$')\n", "\n", "# set up plot figure\n", "fig3 = plt.figure(figsize=(6,6))\n", "fig3.canvas.header_visible = False\n", "fig3.canvas.footer_visible = False\n", "fig3.canvas.toolbar_visible = False\n", "ax3 = plt.axes(projection='3d')\n", "\n", "# create grid of beta values\n", "beta1s = np.linspace(beta1_min, beta1_max, 100)\n", "beta2s = np.linspace(beta2_min, beta2_max, 100)\n", "beta1sv, beta2sv = np.meshgrid(beta1s, beta2s)\n", "\n", "def plot_likelihood(beta0):\n", " # create grid of likelihood values\n", " lsv = vectorized_likelihood(np.full(beta1sv.shape, beta0), beta1sv, beta2sv, X_iris, y_iris)\n", " ax3.plot_surface(beta1sv, beta2sv, lsv, alpha=0.4, color='blue')\n", " ax3.set_xlabel(r'$\\beta_1$')\n", " ax3.set_ylabel(r'$\\beta_2$')\n", " ax3.set_zlabel('likelihood')\n", " ax3.set_xlim3d((beta1_min,beta1_max)); ax3.set_ylim3d((beta2_min, beta2_max)); ax3.set_zlim3d((0,4e-05))\n", "\n", "def plot_likelihood_point(beta0, beta1, beta2):\n", " value = likelihood(beta0, beta1, beta2, X_iris, y_iris)\n", " ax3.scatter(beta1, beta2, value, color='red')\n", " ax3.set_title(\"likelihood value = \"+str(value))\n", "\n", "def update_likelihood(beta0, beta1, beta2):\n", " ax3.clear()\n", " plot_likelihood(beta0) \n", " plot_likelihood_point(beta0, beta1, beta2)\n", " fig3.canvas.draw()\n", "\n", "# create initial plot and monitor sliders\n", "update_likelihood(slider30.value, slider31.value, slider32.value)\n", "widgets.interactive(update_likelihood, beta0=slider30, beta1=slider31, beta2=slider32)" ] }, { "cell_type": "markdown", "id": "ee21353e-3d54-47ab-bb3e-edc736c25a54", "metadata": {}, "source": [ "# Activity 3: log-likelihood function" ] }, { "cell_type": "code", "execution_count": null, "id": "74ae29a0-b5e2-4fd0-8ef7-51f925e9f1ff", "metadata": {}, "outputs": [], "source": [ "def log_likelihood(beta0, beta1, beta2, X, y):\n", " def sigmoid(x):\n", " return 1/(1+np.exp(-x))\n", " p = sigmoid(beta0+beta1*X.T[0]+beta2*X.T[1])\n", " # shift boundary values\n", " p[p==0] = 1e-10\n", " p[p==1] = 1-1e-10\n", " return np.sum(y*np.log(p)+(1-y)*np.log(1-p))\n", "\n", "def vectorized_log_likelihood(beta0, beta1, beta2, X, y):\n", " def sigmoid(x):\n", " return 1/(1+np.exp(-x))\n", " p = sigmoid(beta0[:,:,np.newaxis]+beta1[:,:,np.newaxis]*X.T[0]+beta2[:,:,np.newaxis]*X.T[1])\n", " # shift boundary values\n", " p[p==0] = 1e-10\n", " p[p==1] = 1-1e-10\n", " return np.sum(y*np.log(p)+(1-y)*np.log(1-p), axis=2)\n", "\n", "# initialize beta ranges\n", "beta0_min, beta0_max = -90, -20\n", "beta1_min, beta1_max = 2, 10\n", "beta2_min, beta2_max = 5, 15\n", "\n", "# set up sliders\n", "slider40 = widgets.FloatSlider(min=beta0_min, max=beta0_max, value=(beta0_min+beta0_max)/2, description=r'$\\beta_0:$')\n", "slider41 = widgets.FloatSlider(min=beta1_min, max=beta1_max, value=(beta1_min+beta1_max)/2, description=r'$\\beta_1:$')\n", "slider42 = widgets.FloatSlider(min=beta2_min, max=beta2_max, value=(beta2_min+beta2_max)/2, description=r'$\\beta_2:$')\n", "\n", "# set up plot figure\n", "fig4 = plt.figure(figsize=(6,6))\n", "fig4.canvas.header_visible = False\n", "fig4.canvas.footer_visible = False\n", "fig4.canvas.toolbar_visible = False\n", "ax4 = plt.axes(projection='3d')\n", "\n", "# create grid of beta values\n", "beta1s = np.linspace(beta1_min, beta1_max, 100)\n", "beta2s = np.linspace(beta2_min, beta2_max, 100)\n", "beta1sv, beta2sv = np.meshgrid(beta1s, beta2s)\n", "\n", "def plot_log_likelihood(beta0):\n", " # create grid of log-likelihood values\n", " llsv = vectorized_log_likelihood(np.full(beta1sv.shape, beta0), beta1sv, beta2sv, X_iris, y_iris)\n", " ax4.plot_surface(beta1sv, beta2sv, llsv, alpha=0.4, color='blue')\n", " ax4.set_xlabel(r'$\\beta_1$')\n", " ax4.set_ylabel(r'$\\beta_2$')\n", " ax4.set_zlabel('log-likelihood')\n", " ax4.set_xlim3d((beta1_min,beta1_max)); ax4.set_ylim3d((beta2_min, beta2_max)); ax4.set_zlim3d((-2000,0))\n", "\n", "def plot_log_likelihood_point(beta0, beta1, beta2):\n", " value = log_likelihood(beta0, beta1, beta2, X_iris, y_iris)\n", " ax4.scatter(beta1, beta2, value, color='red')\n", " ax4.set_title(\"log-likelihood value = \"+str(value))\n", "\n", "def update_log_likelihood(beta0, beta1, beta2):\n", " ax4.clear()\n", " plot_log_likelihood(beta0) \n", " plot_log_likelihood_point(beta0, beta1, beta2)\n", " fig4.canvas.draw()\n", "\n", "# create initial plot and monitor sliders\n", "update_log_likelihood(slider40.value, slider41.value, slider42.value)\n", "widgets.interactive(update_log_likelihood, beta0=slider40, beta1=slider41, beta2=slider42)" ] }, { "cell_type": "markdown", "id": "d5ecbf61-f2b2-4d1f-a53c-3b9b035ba025", "metadata": {}, "source": [ "# Minimizing the log-likelihood/maximizing the loss function" ] }, { "cell_type": "code", "execution_count": null, "id": "0c2d791c-c800-413c-82c9-137b9d94b80f", "metadata": {}, "outputs": [], "source": [ "from scipy.optimize import minimize\n", "\n", "def minus_log_likelihood(beta):\n", " return -log_likelihood(beta[0], beta[1], beta[2], X_iris, y_iris)\n", "minimize(minus_log_likelihood, [-10,1,1])" ] }, { "cell_type": "markdown", "id": "cbf37bd6-312c-4af3-82ea-55046fbdd684", "metadata": {}, "source": [ "# Activity 4: adjusting the probability threshold" ] }, { "cell_type": "code", "execution_count": null, "id": "c7162d34-56c9-47b6-8bf2-80c989c626e0", "metadata": {}, "outputs": [], "source": [ "# initialize beta ranges\n", "beta0_min, beta0_max = -90, 0\n", "beta1_min, beta1_max = 2, 10\n", "beta2_min, beta2_max = 5, 15\n", "\n", "# maximum likelihood values\n", "beta0_ml, beta1_ml, beta2_ml = -45.27, 5.755, 10.45\n", "\n", "# initialize x ranges\n", "x1_min, x1_max = -1, 8\n", "x2_min, x2_max = -1, 4\n", "\n", "# set up sliders\n", "slider5p = widgets.FloatSlider(min=0.01, max=0.99, step=0.01, value = 0.5, description = r'$p_\\text{threshold}:$')\n", "\n", "# set up plot figure\n", "fig5 = plt.figure(figsize=(6,6))\n", "fig5.canvas.header_visible = False\n", "fig5.canvas.footer_visible = False\n", "fig5.canvas.toolbar_visible = False\n", "ax5 = plt.axes(projection='3d')\n", "\n", "# create grid of x values\n", "x1s = np.linspace(x1_min, x1_max, 40)\n", "x2s = np.linspace(x2_min, x2_max, 40)\n", "x1sv, x2sv = np.meshgrid(x1s, x2s)\n", "\n", "def plot_pdf_5():\n", " # create grid of heights\n", " def sigmoid(x):\n", " return 1 / (1 + np.exp(-x))\n", " zsv = sigmoid(beta0_ml+beta1_ml*x1sv+beta2_ml*x2sv)\n", " ax5.plot_surface(x1sv, x2sv, zsv, alpha=0.5, color='blue')\n", " ax5.set_xlabel('petal length (cm)')\n", " ax5.set_ylabel('petal width (cm)')\n", " ax5.set_zlabel('probability')\n", " ax5.set_xlim3d((x1_min,x1_max)); ax2.set_ylim3d((x2_min,x2_max)); ax2.set_zlim3d((0,1))\n", "\n", "def plot_data_5(X, y_logical):\n", " ax5.scatter(X[y_logical].T[0], X[y_logical].T[1], color='orange')\n", " ax5.scatter(X[~y_logical].T[0], X[~y_logical].T[1], color='blue')\n", " \n", "def plot_boundary_5(p): \n", " numpoints=40\n", " xline = np.linspace(x1_min,x1_max,numpoints)\n", " yline = -(beta0_ml+beta1_ml*xline+np.log(1/p-1))/beta2_ml\n", " def clip_x(x):\n", " return(x) if x1_min <= x <= x1_max else np.nan\n", " def clip_y(y):\n", " return(y) if x2_min <= y <= x2_max else np.nan\n", " xline = np.array(list(map(clip_x, xline)))\n", " yline = np.array(list(map(clip_y, yline)))\n", " zlinezero = np.full((numpoints,),0)\n", " zlinep = np.full((numpoints,),p)\n", " ax5.plot(xline,yline, zlinezero, color='red', linestyle='dashed')\n", " ax5.plot(xline,yline, zlinep, color='red') \n", "\n", "def update_pdf_plot_5(p):\n", " ax5.clear()\n", " plot_pdf_5() \n", " plot_data_5(X_iris, y_iris_logical)\n", " plot_boundary_5(p)\n", " fig5.canvas.draw()\n", "\n", "# create initial plot and monitor sliders\n", "update_pdf_plot_5(slider5p.value)\n", "widgets.interactive(update_pdf_plot_5, p=slider5p)" ] }, { "cell_type": "markdown", "id": "9ed7e5fe-14d8-4640-8b48-8410ffc2affd", "metadata": {}, "source": [ "# Activity 5: confusion matrix" ] }, { "cell_type": "code", "execution_count": null, "id": "6b3f696e-d467-42a4-b2a0-354b9997bde2", "metadata": {}, "outputs": [], "source": [ "# initialize beta ranges\n", "beta0_min, beta0_max = -80, 0\n", "beta1_min, beta1_max = 2, 10\n", "beta2_min, beta2_max = 5, 15\n", "\n", "# maximum likelihood values\n", "beta0_ml, beta1_ml, beta2_ml = -45.27, 5.755, 10.45\n", "\n", "# calculate probabilities of all data points\n", "p_iris = 1/(1+np.exp(-(beta0_ml+beta1_ml*X_iris.T[0]+beta2_ml*X_iris.T[1])))\n", "\n", "# initialize x ranges\n", "x1_min, x1_max = 0, 8\n", "x2_min, x2_max = 0, 3\n", "\n", "# set up sliders\n", "slider6p = widgets.FloatSlider(min=0.01, max=0.99, step=0.01, value=0.5, description=r'$p_\\text{threshold}:$')\n", "\n", "# set up plot figure\n", "fig6, ax6 = plt.subplots(figsize=(6,6))\n", "fig6.canvas.header_visible = False\n", "fig6.canvas.footer_visible = False\n", "fig6.canvas.toolbar_visible = False\n", "\n", "ax6.set_xlabel(\"petal length (cm)\"); ax1.set_ylabel(\"petal width (cm)\")\n", "ax6.set_xlim((x1_min,x1_max)); ax1.set_ylim((x2_min,x2_max))\n", "\n", "# Initialize line object\n", "line6, = ax6.plot([], [], color='red', linestyle='dashed')\n", "\n", "def plot_data_6(X, y_logical):\n", " ax6.scatter(X[y_logical, 0], X[y_logical, 1], color='orange', label='virginica')\n", " ax6.scatter(X[~y_logical, 0], X[~y_logical, 1], color='blue', label='not virginica')\n", " ax6.legend(loc='upper left')\n", "\n", "def plot_boundary_6(p):\n", " xline = np.array([x1_min,x1_max])\n", " yline = -(beta0_ml+beta1_ml*xline+np.log(1/p-1))/beta2_ml\n", " line6.set_data(xline, yline)\n", " TP = np.sum(np.logical_and(p_iris>=p, y_iris_logical==True))\n", " FP = np.sum(np.logical_and(p_iris>=p, y_iris_logical==False))\n", " FN = np.sum(np.logical_and(p_iris