Morten Hjorth-Jensen
2290 lines
72 KiB
Plaintext
2290 lines
72 KiB
Plaintext
{
|
|
"cells": [
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"<!-- dom:TITLE: Data Analysis and Machine Learning: Logistic Regression -->\n",
|
|
"# Data Analysis and Machine Learning: Logistic Regression\n",
|
|
"<!-- dom:AUTHOR: Morten Hjorth-Jensen at Department of Physics, University of Oslo & Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University -->\n",
|
|
"<!-- Author: --> \n",
|
|
"**Morten Hjorth-Jensen**, Department of Physics, University of Oslo and Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University\n",
|
|
"\n",
|
|
"Date: **Sep 21, 2021**\n",
|
|
"\n",
|
|
"Copyright 1999-2021, Morten Hjorth-Jensen. Released under CC Attribution-NonCommercial 4.0 license\n",
|
|
"\n",
|
|
"\n",
|
|
"\n",
|
|
"\n",
|
|
"## Plans for week 38\n",
|
|
"\n",
|
|
"* Thursday: Summary of regression methods and discussion of project 1. Start Logistic Regression\n",
|
|
"\n",
|
|
"* Friday: Logistic Regression and Optimization methods\n",
|
|
"\n",
|
|
"## Ridge and LASSO Regression, reminder\n",
|
|
"\n",
|
|
"The expression for the standard Mean Squared Error (MSE) which we used to define our cost function and the equations for the ordinary least squares (OLS) method, that is \n",
|
|
"our optimization problem is"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"{\\displaystyle \\min_{\\boldsymbol{\\beta}\\in {\\mathbb{R}}^{p}}}\\frac{1}{n}\\left\\{\\left(\\boldsymbol{y}-\\boldsymbol{X}\\boldsymbol{\\beta}\\right)^T\\left(\\boldsymbol{y}-\\boldsymbol{X}\\boldsymbol{\\beta}\\right)\\right\\}.\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"or we can state it as"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"{\\displaystyle \\min_{\\boldsymbol{\\beta}\\in\n",
|
|
"{\\mathbb{R}}^{p}}}\\frac{1}{n}\\sum_{i=0}^{n-1}\\left(y_i-\\tilde{y}_i\\right)^2=\\frac{1}{n}\\vert\\vert \\boldsymbol{y}-\\boldsymbol{X}\\boldsymbol{\\beta}\\vert\\vert_2^2,\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"where we have used the definition of a norm-2 vector, that is"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"\\vert\\vert \\boldsymbol{x}\\vert\\vert_2 = \\sqrt{\\sum_i x_i^2}.\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"By minimizing the above equation with respect to the parameters\n",
|
|
"$\\boldsymbol{\\beta}$ we could then obtain an analytical expression for the\n",
|
|
"parameters $\\boldsymbol{\\beta}$. We can add a regularization parameter $\\lambda$ by\n",
|
|
"defining a new cost function to be optimized, that is"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"{\\displaystyle \\min_{\\boldsymbol{\\beta}\\in\n",
|
|
"{\\mathbb{R}}^{p}}}\\frac{1}{n}\\vert\\vert \\boldsymbol{y}-\\boldsymbol{X}\\boldsymbol{\\beta}\\vert\\vert_2^2+\\lambda\\vert\\vert \\boldsymbol{\\beta}\\vert\\vert_2^2\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"which leads to the Ridge regression minimization problem where we\n",
|
|
"require that $\\vert\\vert \\boldsymbol{\\beta}\\vert\\vert_2^2\\le t$, where $t$ is\n",
|
|
"a finite number larger than zero. By defining"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"C(\\boldsymbol{X},\\boldsymbol{\\beta})=\\frac{1}{n}\\vert\\vert \\boldsymbol{y}-\\boldsymbol{X}\\boldsymbol{\\beta}\\vert\\vert_2^2+\\lambda\\vert\\vert \\boldsymbol{\\beta}\\vert\\vert_1,\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"we have a new optimization equation"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"{\\displaystyle \\min_{\\boldsymbol{\\beta}\\in\n",
|
|
"{\\mathbb{R}}^{p}}}\\frac{1}{n}\\vert\\vert \\boldsymbol{y}-\\boldsymbol{X}\\boldsymbol{\\beta}\\vert\\vert_2^2+\\lambda\\vert\\vert \\boldsymbol{\\beta}\\vert\\vert_1\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"which leads to Lasso regression. Lasso stands for least absolute shrinkage and selection operator. \n",
|
|
"\n",
|
|
"Here we have defined the norm-1 as"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"\\vert\\vert \\boldsymbol{x}\\vert\\vert_1 = \\sum_i \\vert x_i\\vert.\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"<!-- !split -->\n",
|
|
"## Various steps in cross-validation\n",
|
|
"\n",
|
|
"When the repetitive splitting of the data set is done randomly,\n",
|
|
"samples may accidently end up in a fast majority of the splits in\n",
|
|
"either training or test set. Such samples may have an unbalanced\n",
|
|
"influence on either model building or prediction evaluation. To avoid\n",
|
|
"this $k$-fold cross-validation structures the data splitting. The\n",
|
|
"samples are divided into $k$ more or less equally sized exhaustive and\n",
|
|
"mutually exclusive subsets. In turn (at each split) one of these\n",
|
|
"subsets plays the role of the test set while the union of the\n",
|
|
"remaining subsets constitutes the training set. Such a splitting\n",
|
|
"warrants a balanced representation of each sample in both training and\n",
|
|
"test set over the splits. Still the division into the $k$ subsets\n",
|
|
"involves a degree of randomness. This may be fully excluded when\n",
|
|
"choosing $k=n$. This particular case is referred to as leave-one-out\n",
|
|
"cross-validation (LOOCV). \n",
|
|
"\n",
|
|
"<!-- !split -->\n",
|
|
"## How to set up the cross-validation for Ridge and/or Lasso\n",
|
|
"\n",
|
|
"* Define a range of interest for the penalty parameter.\n",
|
|
"\n",
|
|
"* Divide the data set into training and test set comprising samples $\\{1, \\ldots, n\\} \\setminus i$ and $\\{ i \\}$, respectively.\n",
|
|
"\n",
|
|
"* Fit the linear regression model by means of ridge estimation for each $\\lambda$ in the grid using the training set, and the corresponding estimate of the error variance $\\boldsymbol{\\sigma}_{-i}^2(\\lambda)$, as"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"\\begin{align*}\n",
|
|
"\\boldsymbol{\\beta}_{-i}(\\lambda) & = ( \\boldsymbol{X}_{-i, \\ast}^{T}\n",
|
|
"\\boldsymbol{X}_{-i, \\ast} + \\lambda \\boldsymbol{I}_{pp})^{-1}\n",
|
|
"\\boldsymbol{X}_{-i, \\ast}^{T} \\boldsymbol{y}_{-i}\n",
|
|
"\\end{align*}\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"* Evaluate the prediction performance of these models on the test set by $\\log\\{L[y_i, \\boldsymbol{X}_{i, \\ast}; \\boldsymbol{\\beta}_{-i}(\\lambda), \\boldsymbol{\\sigma}_{-i}^2(\\lambda)]\\}$. Or, by the prediction error $|y_i - \\boldsymbol{X}_{i, \\ast} \\boldsymbol{\\beta}_{-i}(\\lambda)|$, the relative error, the error squared or the R2 score function.\n",
|
|
"\n",
|
|
"* Repeat the first three steps such that each sample plays the role of the test set once.\n",
|
|
"\n",
|
|
"* Average the prediction performances of the test sets at each grid point of the penalty bias/parameter. It is an estimate of the prediction performance of the model corresponding to this value of the penalty parameter on novel data. It is defined as"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"\\begin{align*}\n",
|
|
"\\frac{1}{n} \\sum_{i = 1}^n \\log\\{L[y_i, \\mathbf{X}_{i, \\ast}; \\boldsymbol{\\beta}_{-i}(\\lambda), \\boldsymbol{\\sigma}_{-i}^2(\\lambda)]\\}.\n",
|
|
"\\end{align*}\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"## Cross-validation in brief\n",
|
|
"\n",
|
|
"For the various values of $k$\n",
|
|
"\n",
|
|
"1. shuffle the dataset randomly.\n",
|
|
"\n",
|
|
"2. Split the dataset into $k$ groups.\n",
|
|
"\n",
|
|
"3. For each unique group:\n",
|
|
"\n",
|
|
"a. Decide which group to use as set for test data\n",
|
|
"\n",
|
|
"b. Take the remaining groups as a training data set\n",
|
|
"\n",
|
|
"c. Fit a model on the training set and evaluate it on the test set\n",
|
|
"\n",
|
|
"d. Retain the evaluation score and discard the model\n",
|
|
"\n",
|
|
"\n",
|
|
"5. Summarize the model using the sample of model evaluation scores\n",
|
|
"\n",
|
|
"## Code Example for Cross-validation and $k$-fold Cross-validation\n",
|
|
"\n",
|
|
"The code here uses Ridge regression with cross-validation (CV) resampling and $k$-fold CV in order to fit a specific polynomial."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"collapsed": false,
|
|
"editable": true
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"%matplotlib inline\n",
|
|
"\n",
|
|
"import numpy as np\n",
|
|
"import matplotlib.pyplot as plt\n",
|
|
"from sklearn.model_selection import KFold\n",
|
|
"from sklearn.linear_model import Ridge\n",
|
|
"from sklearn.model_selection import cross_val_score\n",
|
|
"from sklearn.preprocessing import PolynomialFeatures\n",
|
|
"\n",
|
|
"# A seed just to ensure that the random numbers are the same for every run.\n",
|
|
"# Useful for eventual debugging.\n",
|
|
"np.random.seed(3155)\n",
|
|
"\n",
|
|
"# Generate the data.\n",
|
|
"nsamples = 100\n",
|
|
"x = np.random.randn(nsamples)\n",
|
|
"y = 3*x**2 + np.random.randn(nsamples)\n",
|
|
"\n",
|
|
"## Cross-validation on Ridge regression using KFold only\n",
|
|
"\n",
|
|
"# Decide degree on polynomial to fit\n",
|
|
"poly = PolynomialFeatures(degree = 6)\n",
|
|
"\n",
|
|
"# Decide which values of lambda to use\n",
|
|
"nlambdas = 500\n",
|
|
"lambdas = np.logspace(-3, 5, nlambdas)\n",
|
|
"\n",
|
|
"# Initialize a KFold instance\n",
|
|
"k = 5\n",
|
|
"kfold = KFold(n_splits = k)\n",
|
|
"\n",
|
|
"# Perform the cross-validation to estimate MSE\n",
|
|
"scores_KFold = np.zeros((nlambdas, k))\n",
|
|
"\n",
|
|
"i = 0\n",
|
|
"for lmb in lambdas:\n",
|
|
" ridge = Ridge(alpha = lmb)\n",
|
|
" j = 0\n",
|
|
" for train_inds, test_inds in kfold.split(x):\n",
|
|
" xtrain = x[train_inds]\n",
|
|
" ytrain = y[train_inds]\n",
|
|
"\n",
|
|
" xtest = x[test_inds]\n",
|
|
" ytest = y[test_inds]\n",
|
|
"\n",
|
|
" Xtrain = poly.fit_transform(xtrain[:, np.newaxis])\n",
|
|
" ridge.fit(Xtrain, ytrain[:, np.newaxis])\n",
|
|
"\n",
|
|
" Xtest = poly.fit_transform(xtest[:, np.newaxis])\n",
|
|
" ypred = ridge.predict(Xtest)\n",
|
|
"\n",
|
|
" scores_KFold[i,j] = np.sum((ypred - ytest[:, np.newaxis])**2)/np.size(ypred)\n",
|
|
"\n",
|
|
" j += 1\n",
|
|
" i += 1\n",
|
|
"\n",
|
|
"\n",
|
|
"estimated_mse_KFold = np.mean(scores_KFold, axis = 1)\n",
|
|
"\n",
|
|
"## Cross-validation using cross_val_score from sklearn along with KFold\n",
|
|
"\n",
|
|
"# kfold is an instance initialized above as:\n",
|
|
"# kfold = KFold(n_splits = k)\n",
|
|
"\n",
|
|
"estimated_mse_sklearn = np.zeros(nlambdas)\n",
|
|
"i = 0\n",
|
|
"for lmb in lambdas:\n",
|
|
" ridge = Ridge(alpha = lmb)\n",
|
|
"\n",
|
|
" X = poly.fit_transform(x[:, np.newaxis])\n",
|
|
" estimated_mse_folds = cross_val_score(ridge, X, y[:, np.newaxis], scoring='neg_mean_squared_error', cv=kfold)\n",
|
|
"\n",
|
|
" # cross_val_score return an array containing the estimated negative mse for every fold.\n",
|
|
" # we have to the the mean of every array in order to get an estimate of the mse of the model\n",
|
|
" estimated_mse_sklearn[i] = np.mean(-estimated_mse_folds)\n",
|
|
"\n",
|
|
" i += 1\n",
|
|
"\n",
|
|
"## Plot and compare the slightly different ways to perform cross-validation\n",
|
|
"\n",
|
|
"plt.figure()\n",
|
|
"\n",
|
|
"plt.plot(np.log10(lambdas), estimated_mse_sklearn, label = 'cross_val_score')\n",
|
|
"plt.plot(np.log10(lambdas), estimated_mse_KFold, 'r--', label = 'KFold')\n",
|
|
"\n",
|
|
"plt.xlabel('log10(lambda)')\n",
|
|
"plt.ylabel('mse')\n",
|
|
"\n",
|
|
"plt.legend()\n",
|
|
"\n",
|
|
"plt.show()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"## To think about\n",
|
|
"\n",
|
|
"When you are comparing your own code with for example **Scikit-Learn**'s library, there are some minor things to keep in mind.\n",
|
|
"The example here shows how one can keep the intercept in order to compare own code."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"collapsed": false,
|
|
"editable": true
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"import numpy as np\n",
|
|
"import pandas as pd\n",
|
|
"import matplotlib.pyplot as plt\n",
|
|
"from sklearn.model_selection import train_test_split\n",
|
|
"from sklearn import linear_model\n",
|
|
"\n",
|
|
"def R2(y_data, y_model):\n",
|
|
" return 1 - np.sum((y_data - y_model) ** 2) / np.sum((y_data - np.mean(y_data)) ** 2)\n",
|
|
"def MSE(y_data,y_model):\n",
|
|
" n = np.size(y_model)\n",
|
|
" return np.sum((y_data-y_model)**2)/n\n",
|
|
"\n",
|
|
"\n",
|
|
"# A seed just to ensure that the random numbers are the same for every run.\n",
|
|
"# Useful for eventual debugging.\n",
|
|
"np.random.seed(3155)\n",
|
|
"\n",
|
|
"n = 100\n",
|
|
"x = np.random.rand(n)\n",
|
|
"y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)\n",
|
|
"\n",
|
|
"Maxpolydegree = 20\n",
|
|
"X = np.zeros((n,Maxpolydegree))\n",
|
|
"X[:,0] = 1.0\n",
|
|
"\n",
|
|
"for polydegree in range(1, Maxpolydegree):\n",
|
|
" for degree in range(polydegree):\n",
|
|
" X[:,degree] = x**degree\n",
|
|
"\n",
|
|
"\n",
|
|
"# We split the data in test and training data\n",
|
|
"X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)\n",
|
|
"\n",
|
|
"# matrix inversion to find beta\n",
|
|
"OLSbeta = np.linalg.pinv(X_train.T @ X_train) @ X_train.T @ y_train\n",
|
|
"print(OLSbeta)\n",
|
|
"# and then make the prediction\n",
|
|
"ytildeOLS = X_train @ OLSbeta\n",
|
|
"print(\"Training MSE for OLS\")\n",
|
|
"print(MSE(y_train,ytildeOLS))\n",
|
|
"ypredictOLS = X_test @ OLSbeta\n",
|
|
"print(\"Test MSE OLS\")\n",
|
|
"print(MSE(y_test,ypredictOLS))\n",
|
|
"\n",
|
|
"p = len(OLSbeta)\n",
|
|
"I = np.eye(p,p)\n",
|
|
"# Decide which values of lambda to use\n",
|
|
"nlambdas = 4\n",
|
|
"MSEOwnRidgePredict = np.zeros(nlambdas)\n",
|
|
"MSEOwnRidgeTrain = np.zeros(nlambdas)\n",
|
|
"MSERidgePredict = np.zeros(nlambdas)\n",
|
|
"MSERidgeTrain = np.zeros(nlambdas)\n",
|
|
"\n",
|
|
"lambdas = np.logspace(-4, 4, nlambdas)\n",
|
|
"for i in range(nlambdas):\n",
|
|
" lmb = lambdas[i]\n",
|
|
" OwnRidgeBeta = np.linalg.pinv(X_train.T @ X_train+lmb*I) @ X_train.T @ y_train\n",
|
|
" # include lasso using Scikit-Learn\n",
|
|
" # Note: we include the intercept\n",
|
|
" RegRidge = linear_model.Ridge(lmb,fit_intercept=False)\n",
|
|
" RegRidge.fit(X_train,y_train)\n",
|
|
" # and then make the prediction\n",
|
|
" ytildeOwnRidge = X_train @ OwnRidgeBeta\n",
|
|
" ypredictOwnRidge = X_test @ OwnRidgeBeta\n",
|
|
" ytildeRidge = RegRidge.predict(X_train)\n",
|
|
" ypredictRidge = RegRidge.predict(X_test)\n",
|
|
" MSEOwnRidgePredict[i] = MSE(y_test,ypredictOwnRidge)\n",
|
|
" MSEOwnRidgeTrain[i] = MSE(y_train,ytildeOwnRidge)\n",
|
|
" MSERidgePredict[i] = MSE(y_test,ypredictRidge)\n",
|
|
" MSERidgeTrain[i] = MSE(y_train,ytildeRidge)\n",
|
|
" print(\"Beta values for own Ridge implementation\")\n",
|
|
" print(OwnRidgeBeta)\n",
|
|
" print(\"Beta values for Scikit-Learn Ridge implementation\")\n",
|
|
" print(RegRidge.coef_)\n",
|
|
"# Now plot the results\n",
|
|
"plt.figure()\n",
|
|
"plt.plot(np.log10(lambdas), MSEOwnRidgeTrain, 'b', label = 'MSE Ridge train')\n",
|
|
"plt.plot(np.log10(lambdas), MSEOwnRidgePredict, 'r', label = 'MSE Ridge Test')\n",
|
|
"plt.plot(np.log10(lambdas), MSERidgeTrain, 'y', label = 'MSE Ridge train')\n",
|
|
"plt.plot(np.log10(lambdas), MSERidgePredict, 'g', label = 'MSE Ridge Test')\n",
|
|
"\n",
|
|
"plt.xlabel('log10(lambda)')\n",
|
|
"plt.ylabel('MSE')\n",
|
|
"plt.legend()\n",
|
|
"plt.show()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"## More complicated Example: The Ising model\n",
|
|
"\n",
|
|
"The one-dimensional Ising model with nearest neighbor interaction, no\n",
|
|
"external field and a constant coupling constant $J$ is given by"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"<!-- Equation labels as ordinary links -->\n",
|
|
"<div id=\"_auto1\"></div>\n",
|
|
"\n",
|
|
"$$\n",
|
|
"\\begin{equation}\n",
|
|
" H = -J \\sum_{k}^L s_k s_{k + 1},\n",
|
|
"\\label{_auto1} \\tag{1}\n",
|
|
"\\end{equation}\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"where $s_i \\in \\{-1, 1\\}$ and $s_{N + 1} = s_1$. The number of spins\n",
|
|
"in the system is determined by $L$. For the one-dimensional system\n",
|
|
"there is no phase transition.\n",
|
|
"\n",
|
|
"We will look at a system of $L = 40$ spins with a coupling constant of\n",
|
|
"$J = 1$. To get enough training data we will generate 10000 states\n",
|
|
"with their respective energies."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"collapsed": false,
|
|
"editable": true
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"import numpy as np\n",
|
|
"import matplotlib.pyplot as plt\n",
|
|
"from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
|
|
"import seaborn as sns\n",
|
|
"import scipy.linalg as scl\n",
|
|
"from sklearn.model_selection import train_test_split\n",
|
|
"import tqdm\n",
|
|
"sns.set(color_codes=True)\n",
|
|
"cmap_args=dict(vmin=-1., vmax=1., cmap='seismic')\n",
|
|
"\n",
|
|
"L = 40\n",
|
|
"n = int(1e4)\n",
|
|
"\n",
|
|
"spins = np.random.choice([-1, 1], size=(n, L))\n",
|
|
"J = 1.0\n",
|
|
"\n",
|
|
"energies = np.zeros(n)\n",
|
|
"\n",
|
|
"for i in range(n):\n",
|
|
" energies[i] = - J * np.dot(spins[i], np.roll(spins[i], 1))"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"Here we use ordinary least squares\n",
|
|
"regression to predict the energy for the nearest neighbor\n",
|
|
"one-dimensional Ising model on a ring, i.e., the endpoints wrap\n",
|
|
"around. We will use linear regression to fit a value for\n",
|
|
"the coupling constant to achieve this.\n",
|
|
"\n",
|
|
"## Reformulating the problem to suit regression\n",
|
|
"\n",
|
|
"A more general form for the one-dimensional Ising model is"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"<!-- Equation labels as ordinary links -->\n",
|
|
"<div id=\"_auto2\"></div>\n",
|
|
"\n",
|
|
"$$\n",
|
|
"\\begin{equation}\n",
|
|
" H = - \\sum_j^L \\sum_k^L s_j s_k J_{jk}.\n",
|
|
"\\label{_auto2} \\tag{2}\n",
|
|
"\\end{equation}\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"Here we allow for interactions beyond the nearest neighbors and a state dependent\n",
|
|
"coupling constant. This latter expression can be formulated as\n",
|
|
"a matrix-product"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"<!-- Equation labels as ordinary links -->\n",
|
|
"<div id=\"_auto3\"></div>\n",
|
|
"\n",
|
|
"$$\n",
|
|
"\\begin{equation}\n",
|
|
" \\boldsymbol{H} = \\boldsymbol{X} J,\n",
|
|
"\\label{_auto3} \\tag{3}\n",
|
|
"\\end{equation}\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"where $X_{jk} = s_j s_k$ and $J$ is a matrix which consists of the\n",
|
|
"elements $-J_{jk}$. This form of writing the energy fits perfectly\n",
|
|
"with the form utilized in linear regression, that is"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"<!-- Equation labels as ordinary links -->\n",
|
|
"<div id=\"_auto4\"></div>\n",
|
|
"\n",
|
|
"$$\n",
|
|
"\\begin{equation}\n",
|
|
" \\boldsymbol{y} = \\boldsymbol{X}\\boldsymbol{\\beta} + \\boldsymbol{\\epsilon},\n",
|
|
"\\label{_auto4} \\tag{4}\n",
|
|
"\\end{equation}\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"We split the data in training and test data as discussed in the previous example"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"collapsed": false,
|
|
"editable": true
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"X = np.zeros((n, L ** 2))\n",
|
|
"for i in range(n):\n",
|
|
" X[i] = np.outer(spins[i], spins[i]).ravel()\n",
|
|
"y = energies\n",
|
|
"X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"## Linear regression\n",
|
|
"\n",
|
|
"In the ordinary least squares method we choose the cost function"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"<!-- Equation labels as ordinary links -->\n",
|
|
"<div id=\"_auto5\"></div>\n",
|
|
"\n",
|
|
"$$\n",
|
|
"\\begin{equation}\n",
|
|
" C(\\boldsymbol{X}, \\boldsymbol{\\beta})= \\frac{1}{n}\\left\\{(\\boldsymbol{X}\\boldsymbol{\\beta} - \\boldsymbol{y})^T(\\boldsymbol{X}\\boldsymbol{\\beta} - \\boldsymbol{y})\\right\\}.\n",
|
|
"\\label{_auto5} \\tag{5}\n",
|
|
"\\end{equation}\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"We then find the extremal point of $C$ by taking the derivative with respect to $\\boldsymbol{\\beta}$ as discussed above.\n",
|
|
"This yields the expression for $\\boldsymbol{\\beta}$ to be"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"\\boldsymbol{\\beta} = \\frac{\\boldsymbol{X}^T \\boldsymbol{y}}{\\boldsymbol{X}^T \\boldsymbol{X}},\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"which immediately imposes some requirements on $\\boldsymbol{X}$ as there must exist\n",
|
|
"an inverse of $\\boldsymbol{X}^T \\boldsymbol{X}$. If the expression we are modeling contains an\n",
|
|
"intercept, i.e., a constant term, we must make sure that the\n",
|
|
"first column of $\\boldsymbol{X}$ consists of $1$. We do this here"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"collapsed": false,
|
|
"editable": true
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"X_train_own = np.concatenate(\n",
|
|
" (np.ones(len(X_train))[:, np.newaxis], X_train),\n",
|
|
" axis=1\n",
|
|
")\n",
|
|
"X_test_own = np.concatenate(\n",
|
|
" (np.ones(len(X_test))[:, np.newaxis], X_test),\n",
|
|
" axis=1\n",
|
|
")"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"collapsed": false,
|
|
"editable": true
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"def ols_inv(x: np.ndarray, y: np.ndarray) -> np.ndarray:\n",
|
|
" return scl.inv(x.T @ x) @ (x.T @ y)\n",
|
|
"beta = ols_inv(X_train_own, y_train)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"## Singular Value decomposition\n",
|
|
"\n",
|
|
"Doing the inversion directly turns out to be a bad idea since the matrix\n",
|
|
"$\\boldsymbol{X}^T\\boldsymbol{X}$ is singular. An alternative approach is to use the **singular\n",
|
|
"value decomposition**. Using the definition of the Moore-Penrose\n",
|
|
"pseudoinverse we can write the equation for $\\boldsymbol{\\beta}$ as"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"\\boldsymbol{\\beta} = \\boldsymbol{X}^{+}\\boldsymbol{y},\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"where the pseudoinverse of $\\boldsymbol{X}$ is given by"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"\\boldsymbol{X}^{+} = \\frac{\\boldsymbol{X}^T}{\\boldsymbol{X}^T\\boldsymbol{X}}.\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"Using singular value decomposition we can decompose the matrix $\\boldsymbol{X} = \\boldsymbol{U}\\boldsymbol{\\Sigma} \\boldsymbol{V}^T$,\n",
|
|
"where $\\boldsymbol{U}$ and $\\boldsymbol{V}$ are orthogonal(unitary) matrices and $\\boldsymbol{\\Sigma}$ contains the singular values (more details below).\n",
|
|
"where $X^{+} = V\\Sigma^{+} U^T$. This reduces the equation for\n",
|
|
"$\\omega$ to"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"<!-- Equation labels as ordinary links -->\n",
|
|
"<div id=\"_auto6\"></div>\n",
|
|
"\n",
|
|
"$$\n",
|
|
"\\begin{equation}\n",
|
|
" \\boldsymbol{\\beta} = \\boldsymbol{V}\\boldsymbol{\\Sigma}^{+} \\boldsymbol{U}^T \\boldsymbol{y}.\n",
|
|
"\\label{_auto6} \\tag{6}\n",
|
|
"\\end{equation}\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"Note that solving this equation by actually doing the pseudoinverse\n",
|
|
"(which is what we will do) is not a good idea as this operation scales\n",
|
|
"as $\\mathcal{O}(n^3)$, where $n$ is the number of elements in a\n",
|
|
"general matrix. Instead, doing $QR$-factorization and solving the\n",
|
|
"linear system as an equation would reduce this down to\n",
|
|
"$\\mathcal{O}(n^2)$ operations."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"collapsed": false,
|
|
"editable": true
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"def ols_svd(x: np.ndarray, y: np.ndarray) -> np.ndarray:\n",
|
|
" u, s, v = scl.svd(x)\n",
|
|
" return v.T @ scl.pinv(scl.diagsvd(s, u.shape[0], v.shape[0])) @ u.T @ y"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"collapsed": false,
|
|
"editable": true
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"beta = ols_svd(X_train_own,y_train)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"When extracting the $J$-matrix we need to make sure that we remove the intercept, as is done here"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"collapsed": false,
|
|
"editable": true
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"J = beta[1:].reshape(L, L)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"A way of looking at the coefficients in $J$ is to plot the matrices as images."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"collapsed": false,
|
|
"editable": true
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"fig = plt.figure(figsize=(20, 14))\n",
|
|
"im = plt.imshow(J, **cmap_args)\n",
|
|
"plt.title(\"OLS\", fontsize=18)\n",
|
|
"plt.xticks(fontsize=18)\n",
|
|
"plt.yticks(fontsize=18)\n",
|
|
"cb = fig.colorbar(im)\n",
|
|
"cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=18)\n",
|
|
"plt.show()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"It is interesting to note that OLS\n",
|
|
"considers both $J_{j, j + 1} = -0.5$ and $J_{j, j - 1} = -0.5$ as\n",
|
|
"valid matrix elements for $J$.\n",
|
|
"In our discussion below on hyperparameters and Ridge and Lasso regression we will see that\n",
|
|
"this problem can be removed, partly and only with Lasso regression. \n",
|
|
"\n",
|
|
"In this case our matrix inversion was actually possible. The obvious question now is what is the mathematics behind the SVD?\n",
|
|
"\n",
|
|
"\n",
|
|
"\n",
|
|
"\n",
|
|
"\n",
|
|
"## The one-dimensional Ising model\n",
|
|
"\n",
|
|
"Let us bring back the Ising model again, but now with an additional\n",
|
|
"focus on Ridge and Lasso regression as well. We repeat some of the\n",
|
|
"basic parts of the Ising model and the setup of the training and test\n",
|
|
"data. The one-dimensional Ising model with nearest neighbor\n",
|
|
"interaction, no external field and a constant coupling constant $J$ is\n",
|
|
"given by"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"<!-- Equation labels as ordinary links -->\n",
|
|
"<div id=\"_auto7\"></div>\n",
|
|
"\n",
|
|
"$$\n",
|
|
"\\begin{equation}\n",
|
|
" H = -J \\sum_{k}^L s_k s_{k + 1},\n",
|
|
"\\label{_auto7} \\tag{7}\n",
|
|
"\\end{equation}\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"where $s_i \\in \\{-1, 1\\}$ and $s_{N + 1} = s_1$. The number of spins in the system is determined by $L$. For the one-dimensional system there is no phase transition.\n",
|
|
"\n",
|
|
"We will look at a system of $L = 40$ spins with a coupling constant of $J = 1$. To get enough training data we will generate 10000 states with their respective energies."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"collapsed": false,
|
|
"editable": true
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"import numpy as np\n",
|
|
"import matplotlib.pyplot as plt\n",
|
|
"from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
|
|
"import seaborn as sns\n",
|
|
"import scipy.linalg as scl\n",
|
|
"from sklearn.model_selection import train_test_split\n",
|
|
"import sklearn.linear_model as skl\n",
|
|
"import tqdm\n",
|
|
"sns.set(color_codes=True)\n",
|
|
"cmap_args=dict(vmin=-1., vmax=1., cmap='seismic')\n",
|
|
"\n",
|
|
"L = 40\n",
|
|
"n = int(1e4)\n",
|
|
"\n",
|
|
"spins = np.random.choice([-1, 1], size=(n, L))\n",
|
|
"J = 1.0\n",
|
|
"\n",
|
|
"energies = np.zeros(n)\n",
|
|
"\n",
|
|
"for i in range(n):\n",
|
|
" energies[i] = - J * np.dot(spins[i], np.roll(spins[i], 1))"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"A more general form for the one-dimensional Ising model is"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"<!-- Equation labels as ordinary links -->\n",
|
|
"<div id=\"_auto8\"></div>\n",
|
|
"\n",
|
|
"$$\n",
|
|
"\\begin{equation}\n",
|
|
" H = - \\sum_j^L \\sum_k^L s_j s_k J_{jk}.\n",
|
|
"\\label{_auto8} \\tag{8}\n",
|
|
"\\end{equation}\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"Here we allow for interactions beyond the nearest neighbors and a more\n",
|
|
"adaptive coupling matrix. This latter expression can be formulated as\n",
|
|
"a matrix-product on the form"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"<!-- Equation labels as ordinary links -->\n",
|
|
"<div id=\"_auto9\"></div>\n",
|
|
"\n",
|
|
"$$\n",
|
|
"\\begin{equation}\n",
|
|
" H = X J,\n",
|
|
"\\label{_auto9} \\tag{9}\n",
|
|
"\\end{equation}\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"where $X_{jk} = s_j s_k$ and $J$ is the matrix consisting of the\n",
|
|
"elements $-J_{jk}$. This form of writing the energy fits perfectly\n",
|
|
"with the form utilized in linear regression, viz."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"<!-- Equation labels as ordinary links -->\n",
|
|
"<div id=\"_auto10\"></div>\n",
|
|
"\n",
|
|
"$$\n",
|
|
"\\begin{equation}\n",
|
|
" \\boldsymbol{y} = \\boldsymbol{X}\\boldsymbol{\\beta} + \\boldsymbol{\\epsilon}.\n",
|
|
"\\label{_auto10} \\tag{10}\n",
|
|
"\\end{equation}\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"We organize the data as we did above"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"collapsed": false,
|
|
"editable": true
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"X = np.zeros((n, L ** 2))\n",
|
|
"for i in range(n):\n",
|
|
" X[i] = np.outer(spins[i], spins[i]).ravel()\n",
|
|
"y = energies\n",
|
|
"X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.96)\n",
|
|
"\n",
|
|
"X_train_own = np.concatenate(\n",
|
|
" (np.ones(len(X_train))[:, np.newaxis], X_train),\n",
|
|
" axis=1\n",
|
|
")\n",
|
|
"\n",
|
|
"X_test_own = np.concatenate(\n",
|
|
" (np.ones(len(X_test))[:, np.newaxis], X_test),\n",
|
|
" axis=1\n",
|
|
")"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"We will do all fitting with **Scikit-Learn**,"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"collapsed": false,
|
|
"editable": true
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"clf = skl.LinearRegression().fit(X_train, y_train)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"When extracting the $J$-matrix we make sure to remove the intercept"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"collapsed": false,
|
|
"editable": true
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"J_sk = clf.coef_.reshape(L, L)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"And then we plot the results"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"collapsed": false,
|
|
"editable": true
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"fig = plt.figure(figsize=(20, 14))\n",
|
|
"im = plt.imshow(J_sk, **cmap_args)\n",
|
|
"plt.title(\"LinearRegression from Scikit-learn\", fontsize=18)\n",
|
|
"plt.xticks(fontsize=18)\n",
|
|
"plt.yticks(fontsize=18)\n",
|
|
"cb = fig.colorbar(im)\n",
|
|
"cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=18)\n",
|
|
"plt.show()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"The results perfectly with our previous discussion where we used our own code.\n",
|
|
"\n",
|
|
"## Ridge regression\n",
|
|
"\n",
|
|
"Having explored the ordinary least squares we move on to ridge\n",
|
|
"regression. In ridge regression we include a **regularizer**. This\n",
|
|
"involves a new cost function which leads to a new estimate for the\n",
|
|
"weights $\\boldsymbol{\\beta}$. This results in a penalized regression problem. The\n",
|
|
"cost function is given by"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"2\n",
|
|
"2\n",
|
|
" \n",
|
|
"<\n",
|
|
"<\n",
|
|
"<\n",
|
|
"!\n",
|
|
"!\n",
|
|
"M\n",
|
|
"A\n",
|
|
"T\n",
|
|
"H\n",
|
|
"_\n",
|
|
"B\n",
|
|
"L\n",
|
|
"O\n",
|
|
"C\n",
|
|
"K"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"collapsed": false,
|
|
"editable": true
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"_lambda = 0.1\n",
|
|
"clf_ridge = skl.Ridge(alpha=_lambda).fit(X_train, y_train)\n",
|
|
"J_ridge_sk = clf_ridge.coef_.reshape(L, L)\n",
|
|
"fig = plt.figure(figsize=(20, 14))\n",
|
|
"im = plt.imshow(J_ridge_sk, **cmap_args)\n",
|
|
"plt.title(\"Ridge from Scikit-learn\", fontsize=18)\n",
|
|
"plt.xticks(fontsize=18)\n",
|
|
"plt.yticks(fontsize=18)\n",
|
|
"cb = fig.colorbar(im)\n",
|
|
"cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=18)\n",
|
|
"\n",
|
|
"plt.show()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"## LASSO regression\n",
|
|
"\n",
|
|
"In the **Least Absolute Shrinkage and Selection Operator** (LASSO)-method we get a third cost function."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"<!-- Equation labels as ordinary links -->\n",
|
|
"<div id=\"_auto12\"></div>\n",
|
|
"\n",
|
|
"$$\n",
|
|
"\\begin{equation}\n",
|
|
" C(\\boldsymbol{X}, \\boldsymbol{\\beta}; \\lambda) = (\\boldsymbol{X}\\boldsymbol{\\beta} - \\boldsymbol{y})^T(\\boldsymbol{X}\\boldsymbol{\\beta} - \\boldsymbol{y}) + \\lambda \\sqrt{\\boldsymbol{\\beta}^T\\boldsymbol{\\beta}}.\n",
|
|
"\\label{_auto12} \\tag{12}\n",
|
|
"\\end{equation}\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"Finding the extremal point of this cost function is not so straight-forward as in least squares and ridge. We will therefore rely solely on the function ``Lasso`` from **Scikit-Learn**."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"collapsed": false,
|
|
"editable": true
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"clf_lasso = skl.Lasso(alpha=_lambda).fit(X_train, y_train)\n",
|
|
"J_lasso_sk = clf_lasso.coef_.reshape(L, L)\n",
|
|
"fig = plt.figure(figsize=(20, 14))\n",
|
|
"im = plt.imshow(J_lasso_sk, **cmap_args)\n",
|
|
"plt.title(\"Lasso from Scikit-learn\", fontsize=18)\n",
|
|
"plt.xticks(fontsize=18)\n",
|
|
"plt.yticks(fontsize=18)\n",
|
|
"cb = fig.colorbar(im)\n",
|
|
"cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=18)\n",
|
|
"\n",
|
|
"plt.show()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"It is quite striking how LASSO breaks the symmetry of the coupling\n",
|
|
"constant as opposed to ridge and OLS. We get a sparse solution with\n",
|
|
"$J_{j, j + 1} = -1$.\n",
|
|
"\n",
|
|
"\n",
|
|
"\n",
|
|
"## Performance as function of the regularization parameter\n",
|
|
"\n",
|
|
"We see how the different models perform for a different set of values for $\\lambda$."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"collapsed": false,
|
|
"editable": true
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"lambdas = np.logspace(-4, 5, 10)\n",
|
|
"\n",
|
|
"train_errors = {\n",
|
|
" \"ols_sk\": np.zeros(lambdas.size),\n",
|
|
" \"ridge_sk\": np.zeros(lambdas.size),\n",
|
|
" \"lasso_sk\": np.zeros(lambdas.size)\n",
|
|
"}\n",
|
|
"\n",
|
|
"test_errors = {\n",
|
|
" \"ols_sk\": np.zeros(lambdas.size),\n",
|
|
" \"ridge_sk\": np.zeros(lambdas.size),\n",
|
|
" \"lasso_sk\": np.zeros(lambdas.size)\n",
|
|
"}\n",
|
|
"\n",
|
|
"plot_counter = 1\n",
|
|
"\n",
|
|
"fig = plt.figure(figsize=(32, 54))\n",
|
|
"\n",
|
|
"for i, _lambda in enumerate(tqdm.tqdm(lambdas)):\n",
|
|
" for key, method in zip(\n",
|
|
" [\"ols_sk\", \"ridge_sk\", \"lasso_sk\"],\n",
|
|
" [skl.LinearRegression(), skl.Ridge(alpha=_lambda), skl.Lasso(alpha=_lambda)]\n",
|
|
" ):\n",
|
|
" method = method.fit(X_train, y_train)\n",
|
|
"\n",
|
|
" train_errors[key][i] = method.score(X_train, y_train)\n",
|
|
" test_errors[key][i] = method.score(X_test, y_test)\n",
|
|
"\n",
|
|
" omega = method.coef_.reshape(L, L)\n",
|
|
"\n",
|
|
" plt.subplot(10, 5, plot_counter)\n",
|
|
" plt.imshow(omega, **cmap_args)\n",
|
|
" plt.title(r\"%s, $\\lambda = %.4f$\" % (key, _lambda))\n",
|
|
" plot_counter += 1\n",
|
|
"\n",
|
|
"plt.show()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"We see that LASSO reaches a good solution for low\n",
|
|
"values of $\\lambda$, but will \"wither\" when we increase $\\lambda$ too\n",
|
|
"much. Ridge is more stable over a larger range of values for\n",
|
|
"$\\lambda$, but eventually also fades away.\n",
|
|
"\n",
|
|
"## Finding the optimal value of $\\lambda$\n",
|
|
"\n",
|
|
"To determine which value of $\\lambda$ is best we plot the accuracy of\n",
|
|
"the models when predicting the training and the testing set. We expect\n",
|
|
"the accuracy of the training set to be quite good, but if the accuracy\n",
|
|
"of the testing set is much lower this tells us that we might be\n",
|
|
"subject to an overfit model. The ideal scenario is an accuracy on the\n",
|
|
"testing set that is close to the accuracy of the training set."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"collapsed": false,
|
|
"editable": true
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"fig = plt.figure(figsize=(20, 14))\n",
|
|
"\n",
|
|
"colors = {\n",
|
|
" \"ols_sk\": \"r\",\n",
|
|
" \"ridge_sk\": \"y\",\n",
|
|
" \"lasso_sk\": \"c\"\n",
|
|
"}\n",
|
|
"\n",
|
|
"for key in train_errors:\n",
|
|
" plt.semilogx(\n",
|
|
" lambdas,\n",
|
|
" train_errors[key],\n",
|
|
" colors[key],\n",
|
|
" label=\"Train {0}\".format(key),\n",
|
|
" linewidth=4.0\n",
|
|
" )\n",
|
|
"\n",
|
|
"for key in test_errors:\n",
|
|
" plt.semilogx(\n",
|
|
" lambdas,\n",
|
|
" test_errors[key],\n",
|
|
" colors[key] + \"--\",\n",
|
|
" label=\"Test {0}\".format(key),\n",
|
|
" linewidth=4.0\n",
|
|
" )\n",
|
|
"plt.legend(loc=\"best\", fontsize=18)\n",
|
|
"plt.xlabel(r\"$\\lambda$\", fontsize=18)\n",
|
|
"plt.ylabel(r\"$R^2$\", fontsize=18)\n",
|
|
"plt.tick_params(labelsize=18)\n",
|
|
"plt.show()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"From the above figure we can see that LASSO with $\\lambda = 10^{-2}$\n",
|
|
"achieves a very good accuracy on the test set. This by far surpasses the\n",
|
|
"other models for all values of $\\lambda$.\n",
|
|
"\n",
|
|
"\n",
|
|
"\n",
|
|
"\n",
|
|
"\n",
|
|
"\n",
|
|
"\n",
|
|
"<!-- !split -->\n",
|
|
"## Logistic Regression\n",
|
|
"\n",
|
|
"In linear regression our main interest was centered on learning the\n",
|
|
"coefficients of a functional fit (say a polynomial) in order to be\n",
|
|
"able to predict the response of a continuous variable on some unseen\n",
|
|
"data. The fit to the continuous variable $y_i$ is based on some\n",
|
|
"independent variables $\\hat{x}_i$. Linear regression resulted in\n",
|
|
"analytical expressions for standard ordinary Least Squares or Ridge\n",
|
|
"regression (in terms of matrices to invert) for several quantities,\n",
|
|
"ranging from the variance and thereby the confidence intervals of the\n",
|
|
"parameters $\\hat{\\beta}$ to the mean squared error. If we can invert\n",
|
|
"the product of the design matrices, linear regression gives then a\n",
|
|
"simple recipe for fitting our data.\n",
|
|
"\n",
|
|
"<!-- !split -->\n",
|
|
"## Classification problems\n",
|
|
"\n",
|
|
"\n",
|
|
"Classification problems, however, are concerned with outcomes taking\n",
|
|
"the form of discrete variables (i.e. categories). We may for example,\n",
|
|
"on the basis of DNA sequencing for a number of patients, like to find\n",
|
|
"out which mutations are important for a certain disease; or based on\n",
|
|
"scans of various patients' brains, figure out if there is a tumor or\n",
|
|
"not; or given a specific physical system, we'd like to identify its\n",
|
|
"state, say whether it is an ordered or disordered system (typical\n",
|
|
"situation in solid state physics); or classify the status of a\n",
|
|
"patient, whether she/he has a stroke or not and many other similar\n",
|
|
"situations.\n",
|
|
"\n",
|
|
"The most common situation we encounter when we apply logistic\n",
|
|
"regression is that of two possible outcomes, normally denoted as a\n",
|
|
"binary outcome, true or false, positive or negative, success or\n",
|
|
"failure etc.\n",
|
|
"\n",
|
|
"## Optimization and Deep learning\n",
|
|
"\n",
|
|
"Logistic regression will also serve as our stepping stone towards\n",
|
|
"neural network algorithms and supervised deep learning. For logistic\n",
|
|
"learning, the minimization of the cost function leads to a non-linear\n",
|
|
"equation in the parameters $\\hat{\\beta}$. The optimization of the\n",
|
|
"problem calls therefore for minimization algorithms. This forms the\n",
|
|
"bottle neck of all machine learning algorithms, namely how to find\n",
|
|
"reliable minima of a multi-variable function. This leads us to the\n",
|
|
"family of gradient descent methods. The latter are the working horses\n",
|
|
"of basically all modern machine learning algorithms.\n",
|
|
"\n",
|
|
"We note also that many of the topics discussed here on logistic \n",
|
|
"regression are also commonly used in modern supervised Deep Learning\n",
|
|
"models, as we will see later.\n",
|
|
"\n",
|
|
"\n",
|
|
"<!-- !split -->\n",
|
|
"## Basics\n",
|
|
"\n",
|
|
"We consider the case where the dependent variables, also called the\n",
|
|
"responses or the outcomes, $y_i$ are discrete and only take values\n",
|
|
"from $k=0,\\dots,K-1$ (i.e. $K$ classes).\n",
|
|
"\n",
|
|
"The goal is to predict the\n",
|
|
"output classes from the design matrix $\\hat{X}\\in\\mathbb{R}^{n\\times p}$\n",
|
|
"made of $n$ samples, each of which carries $p$ features or predictors. The\n",
|
|
"primary goal is to identify the classes to which new unseen samples\n",
|
|
"belong.\n",
|
|
"\n",
|
|
"Let us specialize to the case of two classes only, with outputs\n",
|
|
"$y_i=0$ and $y_i=1$. Our outcomes could represent the status of a\n",
|
|
"credit card user that could default or not on her/his credit card\n",
|
|
"debt. That is"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"y_i = \\begin{bmatrix} 0 & \\mathrm{no}\\\\ 1 & \\mathrm{yes} \\end{bmatrix}.\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"## Linear classifier\n",
|
|
"\n",
|
|
"Before moving to the logistic model, let us try to use our linear\n",
|
|
"regression model to classify these two outcomes. We could for example\n",
|
|
"fit a linear model to the default case if $y_i > 0.5$ and the no\n",
|
|
"default case $y_i \\leq 0.5$.\n",
|
|
"\n",
|
|
"We would then have our \n",
|
|
"weighted linear combination, namely"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"<!-- Equation labels as ordinary links -->\n",
|
|
"<div id=\"_auto13\"></div>\n",
|
|
"\n",
|
|
"$$\n",
|
|
"\\begin{equation}\n",
|
|
"\\hat{y} = \\hat{X}^T\\hat{\\beta} + \\hat{\\epsilon},\n",
|
|
"\\label{_auto13} \\tag{13}\n",
|
|
"\\end{equation}\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"where $\\hat{y}$ is a vector representing the possible outcomes, $\\hat{X}$ is our\n",
|
|
"$n\\times p$ design matrix and $\\hat{\\beta}$ represents our estimators/predictors.\n",
|
|
"\n",
|
|
"## Some selected properties\n",
|
|
"\n",
|
|
"The main problem with our function is that it takes values on the\n",
|
|
"entire real axis. In the case of logistic regression, however, the\n",
|
|
"labels $y_i$ are discrete variables. A typical example is the credit\n",
|
|
"card data discussed below here, where we can set the state of\n",
|
|
"defaulting the debt to $y_i=1$ and not to $y_i=0$ for one the persons\n",
|
|
"in the data set (see the full example below).\n",
|
|
"\n",
|
|
"One simple way to get a discrete output is to have sign\n",
|
|
"functions that map the output of a linear regressor to values $\\{0,1\\}$,\n",
|
|
"$f(s_i)=sign(s_i)=1$ if $s_i\\ge 0$ and 0 if otherwise. \n",
|
|
"We will encounter this model in our first demonstration of neural networks. Historically it is called the ``perceptron\" model in the machine learning\n",
|
|
"literature. This model is extremely simple. However, in many cases it is more\n",
|
|
"favorable to use a ``soft\" classifier that outputs\n",
|
|
"the probability of a given category. This leads us to the logistic function.\n",
|
|
"\n",
|
|
"## Simple example\n",
|
|
"\n",
|
|
"The following example on data for coronary heart disease (CHD) as function of age may serve as an illustration. In the code here we read and plot whether a person has had CHD (output = 1) or not (output = 0). This ouput is plotted the person's against age. Clearly, the figure shows that attempting to make a standard linear regression fit may not be very meaningful."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"collapsed": false,
|
|
"editable": true
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"# Common imports\n",
|
|
"import os\n",
|
|
"import numpy as np\n",
|
|
"import pandas as pd\n",
|
|
"import matplotlib.pyplot as plt\n",
|
|
"from sklearn.linear_model import LinearRegression, Ridge, Lasso\n",
|
|
"from sklearn.model_selection import train_test_split\n",
|
|
"from sklearn.utils import resample\n",
|
|
"from sklearn.metrics import mean_squared_error\n",
|
|
"from IPython.display import display\n",
|
|
"from pylab import plt, mpl\n",
|
|
"plt.style.use('seaborn')\n",
|
|
"mpl.rcParams['font.family'] = 'serif'\n",
|
|
"\n",
|
|
"# Where to save the figures and data files\n",
|
|
"PROJECT_ROOT_DIR = \"Results\"\n",
|
|
"FIGURE_ID = \"Results/FigureFiles\"\n",
|
|
"DATA_ID = \"DataFiles/\"\n",
|
|
"\n",
|
|
"if not os.path.exists(PROJECT_ROOT_DIR):\n",
|
|
" os.mkdir(PROJECT_ROOT_DIR)\n",
|
|
"\n",
|
|
"if not os.path.exists(FIGURE_ID):\n",
|
|
" os.makedirs(FIGURE_ID)\n",
|
|
"\n",
|
|
"if not os.path.exists(DATA_ID):\n",
|
|
" os.makedirs(DATA_ID)\n",
|
|
"\n",
|
|
"def image_path(fig_id):\n",
|
|
" return os.path.join(FIGURE_ID, fig_id)\n",
|
|
"\n",
|
|
"def data_path(dat_id):\n",
|
|
" return os.path.join(DATA_ID, dat_id)\n",
|
|
"\n",
|
|
"def save_fig(fig_id):\n",
|
|
" plt.savefig(image_path(fig_id) + \".png\", format='png')\n",
|
|
"\n",
|
|
"infile = open(data_path(\"chddata.csv\"),'r')\n",
|
|
"\n",
|
|
"# Read the chd data as csv file and organize the data into arrays with age group, age, and chd\n",
|
|
"chd = pd.read_csv(infile, names=('ID', 'Age', 'Agegroup', 'CHD'))\n",
|
|
"chd.columns = ['ID', 'Age', 'Agegroup', 'CHD']\n",
|
|
"output = chd['CHD']\n",
|
|
"age = chd['Age']\n",
|
|
"agegroup = chd['Agegroup']\n",
|
|
"numberID = chd['ID'] \n",
|
|
"display(chd)\n",
|
|
"\n",
|
|
"plt.scatter(age, output, marker='o')\n",
|
|
"plt.axis([18,70.0,-0.1, 1.2])\n",
|
|
"plt.xlabel(r'Age')\n",
|
|
"plt.ylabel(r'CHD')\n",
|
|
"plt.title(r'Age distribution and Coronary heart disease')\n",
|
|
"plt.show()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"## Plotting the mean value for each group\n",
|
|
"\n",
|
|
"What we could attempt however is to plot the mean value for each group."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"collapsed": false,
|
|
"editable": true
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"agegroupmean = np.array([0.1, 0.133, 0.250, 0.333, 0.462, 0.625, 0.765, 0.800])\n",
|
|
"group = np.array([1, 2, 3, 4, 5, 6, 7, 8])\n",
|
|
"plt.plot(group, agegroupmean, \"r-\")\n",
|
|
"plt.axis([0,9,0, 1.0])\n",
|
|
"plt.xlabel(r'Age group')\n",
|
|
"plt.ylabel(r'CHD mean values')\n",
|
|
"plt.title(r'Mean values for each age group')\n",
|
|
"plt.show()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"We are now trying to find a function $f(y\\vert x)$, that is a function which gives us an expected value for the output $y$ with a given input $x$.\n",
|
|
"In standard linear regression with a linear dependence on $x$, we would write this in terms of our model"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"f(y_i\\vert x_i)=\\beta_0+\\beta_1 x_i.\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"This expression implies however that $f(y_i\\vert x_i)$ could take any\n",
|
|
"value from minus infinity to plus infinity. If we however let\n",
|
|
"$f(y\\vert y)$ be represented by the mean value, the above example\n",
|
|
"shows us that we can constrain the function to take values between\n",
|
|
"zero and one, that is we have $0 \\le f(y_i\\vert x_i) \\le 1$. Looking\n",
|
|
"at our last curve we see also that it has an S-shaped form. This leads\n",
|
|
"us to a very popular model for the function $f$, namely the so-called\n",
|
|
"Sigmoid function or logistic model. We will consider this function as\n",
|
|
"representing the probability for finding a value of $y_i$ with a given\n",
|
|
"$x_i$.\n",
|
|
"\n",
|
|
"## The logistic function\n",
|
|
"\n",
|
|
"Another widely studied model, is the so-called \n",
|
|
"perceptron model, which is an example of a \"hard classification\" model. We\n",
|
|
"will encounter this model when we discuss neural networks as\n",
|
|
"well. Each datapoint is deterministically assigned to a category (i.e\n",
|
|
"$y_i=0$ or $y_i=1$). In many cases, and the coronary heart disease data forms one of many such examples, it is favorable to have a \"soft\"\n",
|
|
"classifier that outputs the probability of a given category rather\n",
|
|
"than a single value. For example, given $x_i$, the classifier\n",
|
|
"outputs the probability of being in a category $k$. Logistic regression\n",
|
|
"is the most common example of a so-called soft classifier. In logistic\n",
|
|
"regression, the probability that a data point $x_i$\n",
|
|
"belongs to a category $y_i=\\{0,1\\}$ is given by the so-called logit function (or Sigmoid) which is meant to represent the likelihood for a given event,"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"p(t) = \\frac{1}{1+\\mathrm \\exp{-t}}=\\frac{\\exp{t}}{1+\\mathrm \\exp{t}}.\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"Note that $1-p(t)= p(-t)$.\n",
|
|
"\n",
|
|
"## Examples of likelihood functions used in logistic regression and nueral networks\n",
|
|
"\n",
|
|
"\n",
|
|
"The following code plots the logistic function, the step function and other functions we will encounter from here and on."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"collapsed": false,
|
|
"editable": true
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"\"\"\"The sigmoid function (or the logistic curve) is a\n",
|
|
"function that takes any real number, z, and outputs a number (0,1).\n",
|
|
"It is useful in neural networks for assigning weights on a relative scale.\n",
|
|
"The value z is the weighted sum of parameters involved in the learning algorithm.\"\"\"\n",
|
|
"\n",
|
|
"import numpy\n",
|
|
"import matplotlib.pyplot as plt\n",
|
|
"import math as mt\n",
|
|
"\n",
|
|
"z = numpy.arange(-5, 5, .1)\n",
|
|
"sigma_fn = numpy.vectorize(lambda z: 1/(1+numpy.exp(-z)))\n",
|
|
"sigma = sigma_fn(z)\n",
|
|
"\n",
|
|
"fig = plt.figure()\n",
|
|
"ax = fig.add_subplot(111)\n",
|
|
"ax.plot(z, sigma)\n",
|
|
"ax.set_ylim([-0.1, 1.1])\n",
|
|
"ax.set_xlim([-5,5])\n",
|
|
"ax.grid(True)\n",
|
|
"ax.set_xlabel('z')\n",
|
|
"ax.set_title('sigmoid function')\n",
|
|
"\n",
|
|
"plt.show()\n",
|
|
"\n",
|
|
"\"\"\"Step Function\"\"\"\n",
|
|
"z = numpy.arange(-5, 5, .02)\n",
|
|
"step_fn = numpy.vectorize(lambda z: 1.0 if z >= 0.0 else 0.0)\n",
|
|
"step = step_fn(z)\n",
|
|
"\n",
|
|
"fig = plt.figure()\n",
|
|
"ax = fig.add_subplot(111)\n",
|
|
"ax.plot(z, step)\n",
|
|
"ax.set_ylim([-0.5, 1.5])\n",
|
|
"ax.set_xlim([-5,5])\n",
|
|
"ax.grid(True)\n",
|
|
"ax.set_xlabel('z')\n",
|
|
"ax.set_title('step function')\n",
|
|
"\n",
|
|
"plt.show()\n",
|
|
"\n",
|
|
"\"\"\"tanh Function\"\"\"\n",
|
|
"z = numpy.arange(-2*mt.pi, 2*mt.pi, 0.1)\n",
|
|
"t = numpy.tanh(z)\n",
|
|
"\n",
|
|
"fig = plt.figure()\n",
|
|
"ax = fig.add_subplot(111)\n",
|
|
"ax.plot(z, t)\n",
|
|
"ax.set_ylim([-1.0, 1.0])\n",
|
|
"ax.set_xlim([-2*mt.pi,2*mt.pi])\n",
|
|
"ax.grid(True)\n",
|
|
"ax.set_xlabel('z')\n",
|
|
"ax.set_title('tanh function')\n",
|
|
"\n",
|
|
"plt.show()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"## Two parameters\n",
|
|
"\n",
|
|
"We assume now that we have two classes with $y_i$ either $0$ or $1$. Furthermore we assume also that we have only two parameters $\\beta$ in our fitting of the Sigmoid function, that is we define probabilities"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"\\begin{align*}\n",
|
|
"p(y_i=1|x_i,\\hat{\\beta}) &= \\frac{\\exp{(\\beta_0+\\beta_1x_i)}}{1+\\exp{(\\beta_0+\\beta_1x_i)}},\\nonumber\\\\\n",
|
|
"p(y_i=0|x_i,\\hat{\\beta}) &= 1 - p(y_i=1|x_i,\\hat{\\beta}),\n",
|
|
"\\end{align*}\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"where $\\hat{\\beta}$ are the weights we wish to extract from data, in our case $\\beta_0$ and $\\beta_1$. \n",
|
|
"\n",
|
|
"Note that we used"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"p(y_i=0\\vert x_i, \\hat{\\beta}) = 1-p(y_i=1\\vert x_i, \\hat{\\beta}).\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"<!-- !split -->\n",
|
|
"## Maximum likelihood\n",
|
|
"\n",
|
|
"In order to define the total likelihood for all possible outcomes from a \n",
|
|
"dataset $\\mathcal{D}=\\{(y_i,x_i)\\}$, with the binary labels\n",
|
|
"$y_i\\in\\{0,1\\}$ and where the data points are drawn independently, we use the so-called [Maximum Likelihood Estimation](https://en.wikipedia.org/wiki/Maximum_likelihood_estimation) (MLE) principle. \n",
|
|
"We aim thus at maximizing \n",
|
|
"the probability of seeing the observed data. We can then approximate the \n",
|
|
"likelihood in terms of the product of the individual probabilities of a specific outcome $y_i$, that is"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"\\begin{align*}\n",
|
|
"P(\\mathcal{D}|\\hat{\\beta})& = \\prod_{i=1}^n \\left[p(y_i=1|x_i,\\hat{\\beta})\\right]^{y_i}\\left[1-p(y_i=1|x_i,\\hat{\\beta}))\\right]^{1-y_i}\\nonumber \\\\\n",
|
|
"\\end{align*}\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"from which we obtain the log-likelihood and our **cost/loss** function"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"\\mathcal{C}(\\hat{\\beta}) = \\sum_{i=1}^n \\left( y_i\\log{p(y_i=1|x_i,\\hat{\\beta})} + (1-y_i)\\log\\left[1-p(y_i=1|x_i,\\hat{\\beta}))\\right]\\right).\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"## The cost function rewritten\n",
|
|
"\n",
|
|
"Reordering the logarithms, we can rewrite the **cost/loss** function as"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"\\mathcal{C}(\\hat{\\beta}) = \\sum_{i=1}^n \\left(y_i(\\beta_0+\\beta_1x_i) -\\log{(1+\\exp{(\\beta_0+\\beta_1x_i)})}\\right).\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"The maximum likelihood estimator is defined as the set of parameters that maximize the log-likelihood where we maximize with respect to $\\beta$.\n",
|
|
"Since the cost (error) function is just the negative log-likelihood, for logistic regression we have that"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"\\mathcal{C}(\\hat{\\beta})=-\\sum_{i=1}^n \\left(y_i(\\beta_0+\\beta_1x_i) -\\log{(1+\\exp{(\\beta_0+\\beta_1x_i)})}\\right).\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"This equation is known in statistics as the **cross entropy**. Finally, we note that just as in linear regression, \n",
|
|
"in practice we often supplement the cross-entropy with additional regularization terms, usually $L_1$ and $L_2$ regularization as we did for Ridge and Lasso regression.\n",
|
|
"\n",
|
|
"## Minimizing the cross entropy\n",
|
|
"\n",
|
|
"The cross entropy is a convex function of the weights $\\hat{\\beta}$ and,\n",
|
|
"therefore, any local minimizer is a global minimizer. \n",
|
|
"\n",
|
|
"\n",
|
|
"Minimizing this\n",
|
|
"cost function with respect to the two parameters $\\beta_0$ and $\\beta_1$ we obtain"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"\\frac{\\partial \\mathcal{C}(\\hat{\\beta})}{\\partial \\beta_0} = -\\sum_{i=1}^n \\left(y_i -\\frac{\\exp{(\\beta_0+\\beta_1x_i)}}{1+\\exp{(\\beta_0+\\beta_1x_i)}}\\right),\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"and"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"\\frac{\\partial \\mathcal{C}(\\hat{\\beta})}{\\partial \\beta_1} = -\\sum_{i=1}^n \\left(y_ix_i -x_i\\frac{\\exp{(\\beta_0+\\beta_1x_i)}}{1+\\exp{(\\beta_0+\\beta_1x_i)}}\\right).\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"## A more compact expression\n",
|
|
"\n",
|
|
"Let us now define a vector $\\hat{y}$ with $n$ elements $y_i$, an\n",
|
|
"$n\\times p$ matrix $\\hat{X}$ which contains the $x_i$ values and a\n",
|
|
"vector $\\hat{p}$ of fitted probabilities $p(y_i\\vert x_i,\\hat{\\beta})$. We can rewrite in a more compact form the first\n",
|
|
"derivative of cost function as"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"\\frac{\\partial \\mathcal{C}(\\hat{\\beta})}{\\partial \\hat{\\beta}} = -\\hat{X}^T\\left(\\hat{y}-\\hat{p}\\right).\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"If we in addition define a diagonal matrix $\\hat{W}$ with elements \n",
|
|
"$p(y_i\\vert x_i,\\hat{\\beta})(1-p(y_i\\vert x_i,\\hat{\\beta})$, we can obtain a compact expression of the second derivative as"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"\\frac{\\partial^2 \\mathcal{C}(\\hat{\\beta})}{\\partial \\hat{\\beta}\\partial \\hat{\\beta}^T} = \\hat{X}^T\\hat{W}\\hat{X}.\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"## Extending to more predictors\n",
|
|
"\n",
|
|
"Within a binary classification problem, we can easily expand our model to include multiple predictors. Our ratio between likelihoods is then with $p$ predictors"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"\\log{ \\frac{p(\\hat{\\beta}\\hat{x})}{1-p(\\hat{\\beta}\\hat{x})}} = \\beta_0+\\beta_1x_1+\\beta_2x_2+\\dots+\\beta_px_p.\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"Here we defined $\\hat{x}=[1,x_1,x_2,\\dots,x_p]$ and $\\hat{\\beta}=[\\beta_0, \\beta_1, \\dots, \\beta_p]$ leading to"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"p(\\hat{\\beta}\\hat{x})=\\frac{ \\exp{(\\beta_0+\\beta_1x_1+\\beta_2x_2+\\dots+\\beta_px_p)}}{1+\\exp{(\\beta_0+\\beta_1x_1+\\beta_2x_2+\\dots+\\beta_px_p)}}.\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"## Including more classes\n",
|
|
"\n",
|
|
"Till now we have mainly focused on two classes, the so-called binary\n",
|
|
"system. Suppose we wish to extend to $K$ classes. Let us for the sake\n",
|
|
"of simplicity assume we have only two predictors. We have then following model"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"\\log{\\frac{p(C=1\\vert x)}{p(K\\vert x)}} = \\beta_{10}+\\beta_{11}x_1,\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"and"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"\\log{\\frac{p(C=2\\vert x)}{p(K\\vert x)}} = \\beta_{20}+\\beta_{21}x_1,\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"and so on till the class $C=K-1$ class"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"\\log{\\frac{p(C=K-1\\vert x)}{p(K\\vert x)}} = \\beta_{(K-1)0}+\\beta_{(K-1)1}x_1,\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"and the model is specified in term of $K-1$ so-called log-odds or\n",
|
|
"**logit** transformations.\n",
|
|
"\n",
|
|
"\n",
|
|
"## More classes\n",
|
|
"\n",
|
|
"In our discussion of neural networks we will encounter the above again\n",
|
|
"in terms of a slightly modified function, the so-called **Softmax** function.\n",
|
|
"\n",
|
|
"The softmax function is used in various multiclass classification\n",
|
|
"methods, such as multinomial logistic regression (also known as\n",
|
|
"softmax regression), multiclass linear discriminant analysis, naive\n",
|
|
"Bayes classifiers, and artificial neural networks. Specifically, in\n",
|
|
"multinomial logistic regression and linear discriminant analysis, the\n",
|
|
"input to the function is the result of $K$ distinct linear functions,\n",
|
|
"and the predicted probability for the $k$-th class given a sample\n",
|
|
"vector $\\hat{x}$ and a weighting vector $\\hat{\\beta}$ is (with two\n",
|
|
"predictors):"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"p(C=k\\vert \\mathbf {x} )=\\frac{\\exp{(\\beta_{k0}+\\beta_{k1}x_1)}}{1+\\sum_{l=1}^{K-1}\\exp{(\\beta_{l0}+\\beta_{l1}x_1)}}.\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"It is easy to extend to more predictors. The final class is"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"p(C=K\\vert \\mathbf {x} )=\\frac{1}{1+\\sum_{l=1}^{K-1}\\exp{(\\beta_{l0}+\\beta_{l1}x_1)}},\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"and they sum to one. Our earlier discussions were all specialized to\n",
|
|
"the case with two classes only. It is easy to see from the above that\n",
|
|
"what we derived earlier is compatible with these equations.\n",
|
|
"\n",
|
|
"To find the optimal parameters we would typically use a gradient\n",
|
|
"descent method. Newton's method and gradient descent methods are\n",
|
|
"discussed in the material on [optimization\n",
|
|
"methods](https://compphysics.github.io/MachineLearning/doc/pub/Splines/html/Splines-bs.html).\n",
|
|
"\n",
|
|
"This will be discussed next week. Before we develop our own codes for logistic regression, we end this lecture by studying the functionality that **Scikit-learn** offers. \n",
|
|
"\n",
|
|
"\n",
|
|
"\n",
|
|
"\n",
|
|
"\n",
|
|
"## Wisconsin Cancer Data\n",
|
|
"\n",
|
|
"We show here how we can use a simple regression case on the breast\n",
|
|
"cancer data using Logistic regression as our algorithm for\n",
|
|
"classification."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"collapsed": false,
|
|
"editable": true
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"import matplotlib.pyplot as plt\n",
|
|
"import numpy as np\n",
|
|
"from sklearn.model_selection import train_test_split \n",
|
|
"from sklearn.datasets import load_breast_cancer\n",
|
|
"from sklearn.linear_model import LogisticRegression\n",
|
|
"\n",
|
|
"# Load the data\n",
|
|
"cancer = load_breast_cancer()\n",
|
|
"\n",
|
|
"X_train, X_test, y_train, y_test = train_test_split(cancer.data,cancer.target,random_state=0)\n",
|
|
"print(X_train.shape)\n",
|
|
"print(X_test.shape)\n",
|
|
"# Logistic Regression\n",
|
|
"logreg = LogisticRegression(solver='lbfgs')\n",
|
|
"logreg.fit(X_train, y_train)\n",
|
|
"print(\"Test set accuracy with Logistic Regression: {:.2f}\".format(logreg.score(X_test,y_test)))\n",
|
|
"#now scale the data\n",
|
|
"from sklearn.preprocessing import StandardScaler\n",
|
|
"scaler = StandardScaler()\n",
|
|
"scaler.fit(X_train)\n",
|
|
"X_train_scaled = scaler.transform(X_train)\n",
|
|
"X_test_scaled = scaler.transform(X_test)\n",
|
|
"# Logistic Regression\n",
|
|
"logreg.fit(X_train_scaled, y_train)\n",
|
|
"print(\"Test set accuracy Logistic Regression with scaled data: {:.2f}\".format(logreg.score(X_test_scaled,y_test)))"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"## Using the correlation matrix\n",
|
|
"\n",
|
|
"In addition to the above scores, we could also study the covariance (and the correlation matrix).\n",
|
|
"We use **Pandas** to compute the correlation matrix."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"collapsed": false,
|
|
"editable": true
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"import matplotlib.pyplot as plt\n",
|
|
"import numpy as np\n",
|
|
"from sklearn.model_selection import train_test_split \n",
|
|
"from sklearn.datasets import load_breast_cancer\n",
|
|
"from sklearn.linear_model import LogisticRegression\n",
|
|
"cancer = load_breast_cancer()\n",
|
|
"import pandas as pd\n",
|
|
"# Making a data frame\n",
|
|
"cancerpd = pd.DataFrame(cancer.data, columns=cancer.feature_names)\n",
|
|
"\n",
|
|
"fig, axes = plt.subplots(15,2,figsize=(10,20))\n",
|
|
"malignant = cancer.data[cancer.target == 0]\n",
|
|
"benign = cancer.data[cancer.target == 1]\n",
|
|
"ax = axes.ravel()\n",
|
|
"\n",
|
|
"for i in range(30):\n",
|
|
" _, bins = np.histogram(cancer.data[:,i], bins =50)\n",
|
|
" ax[i].hist(malignant[:,i], bins = bins, alpha = 0.5)\n",
|
|
" ax[i].hist(benign[:,i], bins = bins, alpha = 0.5)\n",
|
|
" ax[i].set_title(cancer.feature_names[i])\n",
|
|
" ax[i].set_yticks(())\n",
|
|
"ax[0].set_xlabel(\"Feature magnitude\")\n",
|
|
"ax[0].set_ylabel(\"Frequency\")\n",
|
|
"ax[0].legend([\"Malignant\", \"Benign\"], loc =\"best\")\n",
|
|
"fig.tight_layout()\n",
|
|
"plt.show()\n",
|
|
"\n",
|
|
"import seaborn as sns\n",
|
|
"correlation_matrix = cancerpd.corr().round(1)\n",
|
|
"# use the heatmap function from seaborn to plot the correlation matrix\n",
|
|
"# annot = True to print the values inside the square\n",
|
|
"plt.figure(figsize=(15,8))\n",
|
|
"sns.heatmap(data=correlation_matrix, annot=True)\n",
|
|
"plt.show()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"## Discussing the correlation data\n",
|
|
"\n",
|
|
"In the above example we note two things. In the first plot we display\n",
|
|
"the overlap of benign and malignant tumors as functions of the various\n",
|
|
"features in the Wisconsing breast cancer data set. We see that for\n",
|
|
"some of the features we can distinguish clearly the benign and\n",
|
|
"malignant cases while for other features we cannot. This can point to\n",
|
|
"us which features may be of greater interest when we wish to classify\n",
|
|
"a benign or not benign tumour.\n",
|
|
"\n",
|
|
"In the second figure we have computed the so-called correlation\n",
|
|
"matrix, which in our case with thirty features becomes a $30\\times 30$\n",
|
|
"matrix.\n",
|
|
"\n",
|
|
"We constructed this matrix using **pandas** via the statements"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"collapsed": false,
|
|
"editable": true
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"cancerpd = pd.DataFrame(cancer.data, columns=cancer.feature_names)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"and then"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"collapsed": false,
|
|
"editable": true
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"correlation_matrix = cancerpd.corr().round(1)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"Diagonalizing this matrix we can in turn say something about which\n",
|
|
"features are of relevance and which are not. This leads us to\n",
|
|
"the classical Principal Component Analysis (PCA) theorem with\n",
|
|
"applications. This will be discussed later this semester ([week 43](https://compphysics.github.io/MachineLearning/doc/pub/week43/html/week43-bs.html)).\n",
|
|
"\n",
|
|
"\n",
|
|
"\n",
|
|
"## Other measures in classification studies: Cancer Data again"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"collapsed": false,
|
|
"editable": true
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"import matplotlib.pyplot as plt\n",
|
|
"import numpy as np\n",
|
|
"from sklearn.model_selection import train_test_split \n",
|
|
"from sklearn.datasets import load_breast_cancer\n",
|
|
"from sklearn.linear_model import LogisticRegression\n",
|
|
"\n",
|
|
"# Load the data\n",
|
|
"cancer = load_breast_cancer()\n",
|
|
"\n",
|
|
"X_train, X_test, y_train, y_test = train_test_split(cancer.data,cancer.target,random_state=0)\n",
|
|
"print(X_train.shape)\n",
|
|
"print(X_test.shape)\n",
|
|
"# Logistic Regression\n",
|
|
"logreg = LogisticRegression(solver='lbfgs')\n",
|
|
"logreg.fit(X_train, y_train)\n",
|
|
"print(\"Test set accuracy with Logistic Regression: {:.2f}\".format(logreg.score(X_test,y_test)))\n",
|
|
"#now scale the data\n",
|
|
"from sklearn.preprocessing import StandardScaler\n",
|
|
"scaler = StandardScaler()\n",
|
|
"scaler.fit(X_train)\n",
|
|
"X_train_scaled = scaler.transform(X_train)\n",
|
|
"X_test_scaled = scaler.transform(X_test)\n",
|
|
"# Logistic Regression\n",
|
|
"logreg.fit(X_train_scaled, y_train)\n",
|
|
"print(\"Test set accuracy Logistic Regression with scaled data: {:.2f}\".format(logreg.score(X_test_scaled,y_test)))\n",
|
|
"\n",
|
|
"\n",
|
|
"from sklearn.preprocessing import LabelEncoder\n",
|
|
"from sklearn.model_selection import cross_validate\n",
|
|
"#Cross validation\n",
|
|
"accuracy = cross_validate(logreg,X_test_scaled,y_test,cv=10)['test_score']\n",
|
|
"print(accuracy)\n",
|
|
"print(\"Test set accuracy with Logistic Regression and scaled data: {:.2f}\".format(logreg.score(X_test_scaled,y_test)))\n",
|
|
"\n",
|
|
"\n",
|
|
"import scikitplot as skplt\n",
|
|
"y_pred = logreg.predict(X_test_scaled)\n",
|
|
"skplt.metrics.plot_confusion_matrix(y_test, y_pred, normalize=True)\n",
|
|
"plt.show()\n",
|
|
"y_probas = logreg.predict_proba(X_test_scaled)\n",
|
|
"skplt.metrics.plot_roc(y_test, y_probas)\n",
|
|
"plt.show()\n",
|
|
"skplt.metrics.plot_cumulative_gain(y_test, y_probas)\n",
|
|
"plt.show()"
|
|
]
|
|
}
|
|
],
|
|
"metadata": {},
|
|
"nbformat": 4,
|
|
"nbformat_minor": 4
|
|
}
|