notebooks/chapter02_book.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Linear regression and logistic regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Small example of simple linear regression (one independent variable)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The work directory is given from the function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def towdir(s):\n",
    "    return (str('./datasets_book/'+s))\n",
    "\n",
    "import deepglmlib.utils as utils\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import importlib\n",
    "importlib.reload(utils)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this example, the explaining variable is univariate such that $x_{i}=x_{i1}$, hence:\n",
    "\n",
    "$$y_i = \\beta_0 + \\beta_1 x_{i} \\,.$$\n",
    "\n",
    "In python, we generate the noise, the explicative variable $x_i$, compute $y_i$ and show their linear relation with the line."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n      = 10\n",
    "b0     = -0.5\n",
    "b1     = 3.5\n",
    "xmin   = 0\n",
    "xmax   = 1\n",
    "x      = np.random.uniform(xmin,xmax,n)  \n",
    "e      = np.random.randn(n)/3           \n",
    "y      = b0 + b1 * x + e    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "beta = [b0, b1]\n",
    "\n",
    "np.savetxt(towdir(\"./beta_1d_reglinear.txt\"),beta)\n",
    "\n",
    "np.savetxt(towdir(\"./xy_1d_reglinear.txt\"),\n",
    "           np.hstack([x.reshape((n,1)),y.reshape((n,1))]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "beta = np.loadtxt(towdir(\"./beta_1d_reglinear.txt\"))\n",
    "xy   = np.loadtxt(towdir(\"./xy_1d_reglinear.txt\"))\n",
    "x    = xy[:,0].reshape((len(xy),1))\n",
    "y    = xy[:,1].reshape((len(xy),1))\n",
    "x.shape,y.shape, beta.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "xmin   = 0\n",
    "xmax   = 1\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.plot(np.array([xmin,xmax]), beta[0]+beta[1]*np.array([xmin,xmax]), \"b-\")\n",
    "ax.plot(x, y ,'bo')\n",
    "ax.axis([xmin,xmax,0.80*min(y), 1.20*max(y)])\n",
    "ax.set_xlabel(r'$x$')\n",
    "ax.set_xlabel(r'$y$')\n",
    "ax.set_title(r'Sample points and the real linear regression line')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = np.hstack([np.ones((len(x),1)),x.reshape(len(x),1)])\n",
    "y = y.reshape(len(x),1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example of linear regression with the small example"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Implementation from Python with the algebra\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "betahat_np = np.linalg.inv(X.T.dot(X)).dot(X.T).dot(y)\n",
    "print(\"betahat_np=\")\n",
    "print(betahat_np.reshape(2,1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Implementation from **Numpy**\n",
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "betahat_np2 = np.linalg.lstsq(X, y, rcond =None)[0]\n",
    "print(\"betahat_np2=\")\n",
    "print(betahat_np2.reshape(2,1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Implementation from **Scikit-Learn**\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#import sklearn as sklearn\n",
    "from sklearn.linear_model import LinearRegression\n",
    "fit_skl = LinearRegression(fit_intercept=False).fit(X, y)\n",
    "#yhat_skl = fit_skl.predict(X)\n",
    "#betahat_skl = [ fit_skl.coef_, fit_skl.intercept_]\n",
    "betahat_skl = fit_skl.coef_.reshape(2,1)\n",
    "\n",
    "print(\"betahat_skl=\")\n",
    "print(betahat_skl)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Quality indicator with numpy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def MSE_score(y,yhat):\n",
    "    return np.sum((y-yhat)**2)/len(yhat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def R2_score(y, yhat):\n",
    "    return 1 - np.sum((y - yhat) ** 2) / np.sum((y - np.mean(y)) ** 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "yhat = X @ betahat_skl\n",
    "print(f\" R^2 = {np.round(R2_score(y,yhat),2):2.2f} \\\n",
    "        \\n MSE = {np.round(MSE_score(y,yhat),2):2.2f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " #### Solution for $\\sigma$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "sigmahat = np.sqrt( np.sum((y-np.matmul(X,betahat_skl))**2) / len(y) )\n",
    "print(\"sigmahat=\",np.round(sigmahat,5))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "def f_loglik_gauss(y,beta,sigma):\n",
    "    n = len(y)\n",
    "    f = +n/2*np.log(2*np.pi*sigma**2) + np.sum((y-np.matmul(X,beta))**2) /(2*sigma**2)\n",
    "    return -f\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Indicators with sklearn"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mse_yyhat,r2_yyhat = utils.f_metrics_regression(y,yhat,True,False,None,ndec=3,samplename=\"sample\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Logistic Regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\newcommand{\\bH}{H}$\n",
    "$\\newcommand{\\bW}{W}$\n",
    "$\\newcommand{\\bp}{p}$\n",
    "$\\newcommand{\\bX}{X}$\n",
    "$\\newcommand{\\by}{y}$\n",
    "$\\newcommand{\\hy}{\\hat{y}}$\n",
    "$\\newcommand{\\bhy}{\\hat{y}}$\n",
    "$\\newcommand{\\bx}{x}$\n",
    "$\\newcommand{\\bbeta}{\\beta}$\n",
    "$\\newcommand{\\hbeta}{\\hat{\\beta}}$\n",
    "$\\newcommand{\\bepsilon}{\\epsilon}$\n",
    "$\\newcommand{\\bhbeta}{\\hat{\\beta}}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Small example\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f_sigmoid(a):\n",
    "    return np.exp(a)/(1+np.exp(a))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "n1 = n2 = 50\n",
    "n  = 2 * n1\n",
    "\n",
    "beta = np.array([-0.5,3.5,2.0]).reshape((3,1))\n",
    "\n",
    "x = np.vstack([ np.random.normal(1,1,n).reshape((n1,2)) ,\n",
    "                np.random.normal(-1,1,n).reshape((n2,2)) ])\n",
    "\n",
    "p            = f_sigmoid( np.hstack([ np.ones((n,1)), x]) @ beta ) \n",
    "y            = np.random.binomial(1,p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let save the dataset with the python module panda."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt(towdir(\"./xy_2d_reglogistic.txt\"),np.hstack([x,y]))\n",
    "np.savetxt(towdir(\"./beta_2d_reglogistic.txt\"),beta)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "xy   = np.loadtxt(towdir(\"./xy_2d_reglogistic.txt\",))\n",
    "beta = np.loadtxt(towdir(\"./beta_2d_reglogistic.txt\"))\n",
    "\n",
    "x    = xy[:,[0,1]]\n",
    "y    = xy[:,2]\n",
    "n    = len(y)\n",
    "\n",
    "print(xy.shape,x.shape,y.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n0 = np.sum(y==0)\n",
    "n1 = np.sum(y==1)\n",
    "n0,n1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = y.ravel()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let plot the dataset with the known labels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.plot(x[y.ravel()==0,0], x[y.ravel()==0,1], 'bx', \n",
    "         x[y.ravel()==1,0], x[y.ravel()==1,1], 'bo')\n",
    "ax.set_xlabel(r'$x1$')\n",
    "ax.set_ylabel(r'$x2$')\n",
    "ax.set_title(r'Sample points from two classes')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "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",
    "logreg = LogisticRegression(solver='lbfgs')\n",
    "logreg = logreg.fit(x, y.ravel())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "betahat_rg=np.append(logreg.intercept_,logreg.coef_)\n",
    "betahat_rg.reshape((x.shape[1]+1,1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "yhat = logreg.predict(x)\n",
    "print(\"acc=\",np.sum(y==yhat)/len(y))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Indicators with sklearn"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "acc_yyhat, prc_yyhat, rcc_yyhat, cm_yyhat = \\\n",
    "    utils.f_metrics_classification(y,yhat,True,False,None,ndec=3,samplename=\"sample\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "myenv",
   "language": "python",
   "name": "myenv"
  },
  "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}