06. Regression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6. Regression\n",
    "\n",
    "[Data Science Playlist on YouTube](https://www.youtube.com/watch?v=iOJbOfPHXLg&list=PLLBUgWXdTBDg1Qgmwt4jKtVn9BWh5-zgy)\n",
    "[![Python Data Science](https://apmonitor.com/che263/uploads/Begin_Python/DataScience06.png)](https://www.youtube.com/watch?v=iOJbOfPHXLg&list=PLLBUgWXdTBDg1Qgmwt4jKtVn9BWh5-zgy \"Python Data Science\")\n",
    "\n",
    "Regression is the process of adjusting model parameters to fit a prediction `y` to measured values `z`. There are independent variables `x` as inputs to the model to generate the predictions `y`. For machine learning, the objective is to minimize a loss function by adjusting model parameters. A common loss function is the sum of squared errors between the predicted `y` and measured `z` values.\n",
    "\n",
    "    x = Independent Variable, Input, Feature\n",
    "    y = Dependent Variable, Output, Label\n",
    "    z = Output Measurement\n",
    "\n",
    "![temperature](https://apmonitor.com/che263/uploads/Begin_Python/temperature.png)\n",
    "\n",
    "The is objective is to minimize a loss function such as a sum of squared errors between the measured and predicted values:\n",
    "\n",
    "$Loss = \\sum_{i=1}^{n}\\left(y_i-z_i\\right)^2$\n",
    "\n",
    "where `n` is the number of observations. Regression requires labelled data (output values) for training. Classification, on the other hand, can either be supervised (with `z` measurements, labels) or unsupervised (no labels, `z` measurements). Run the following code to load 30 sample data points with input `x` and measured output `z`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "n = 31\n",
    "x = np.linspace(0,3,n)\n",
    "z = np.array([1.89,-0.12,-0.32,2.11,-0.25,1.01,0.17,2.75,2.01,5.5,\\\n",
    "     0.87,3.31,2.29,2.73,3.67,3.92,4.29,5.27,3.85,4.26,\\\n",
    "     5.75,5.82,6.36,9.13,7.61,9.52,9.53,12.49,12.29,13.7,14.12])\n",
    "data = pd.DataFrame(np.vstack((x,z)).T,columns=['x','z'])\n",
    "plt.plot(x,z,'ro',label='Measured')\n",
    "plt.xlabel('x'); plt.ylabel('z'); plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![idea](https://apmonitor.com/che263/uploads/Begin_Python/idea.png)\n",
    "\n",
    "### Linear Regression\n",
    "\n",
    "There are many model forms such as linear, polynomial, and nonlinear. A familiar linear model is a line with slope `a` and intercept `b`.\n",
    "\n",
    "    y = a x + b\n",
    "    \n",
    "A simple method for linear regression is with `numpy` to fit `p=np.polyfit(x,y,1)` and evaluate `np.polyval(p,x)` the model. Run the following code to determine the slope and intercept that minimize the sum of squared errors (least squares) between the predicted `y` and measured `z` output."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p1 = np.polyfit(x,z,1)\n",
    "\n",
    "print('Slope, Intercept:' + str(p1))\n",
    "\n",
    "plt.plot(x,z,'ro',label='Measured (z)')\n",
    "plt.plot(x,np.polyval(p1,x),'b-',label='Predicted (y)')\n",
    "plt.legend(); plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The $R^2$ value can be calculated with the measured (true) and model (predicted) values for any regression model, not just linear.\n",
    "\n",
    "```python\n",
    "from sklearn.metrics import r2_score\n",
    "meas  = [3.0, 2.0, 1.9, 7.1]\n",
    "model = [2.5, 1.8, 2.0, 8.0]\n",
    "r2_score(meas, model)\n",
    "```\n",
    "\n",
    "Another package is `statsmodels` that performs standard Ordinary Least Squares (OLS) analysis with a nice report summary. The input `x` is augmented with a `np.ones(n)` column so that it also predicts the intercept\n",
    "\n",
    "```python\n",
    "xc = np.vstack((x,np.ones(n))).T\n",
    "```\n",
    "\n",
    "and this is also accomplished more conveniently with `xc=sm.add_constant(x)`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import statsmodels.api as sm\n",
    "xc = sm.add_constant(x)\n",
    "model = sm.OLS(z,xc).fit()\n",
    "predictions = model.predict(xc)\n",
    "model.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![expert](https://apmonitor.com/che263/uploads/Begin_Python/expert.png)\n",
    "\n",
    "### Linear Regression Activity\n",
    "\n",
    "Create a linear model with the data:\n",
    "\n",
    "```python\n",
    "xr = [0.0,1.0,2.0,3.5,5.0]\n",
    "yr = [0.7,0.55,0.34,0.3,0.2]\n",
    "```\n",
    "\n",
    "Calculate the $R^2$ value and show the data and linear fit on a plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![debug](https://apmonitor.com/che263/uploads/Begin_Python/debug.png)\n",
    "\n",
    "### Polynomial Regression\n",
    "\n",
    "A polynomial model may also be quadratic:\n",
    "\n",
    "    y = a x^2 + b x + c\n",
    "    \n",
    "A quadratic model is really just a linear model with two inputs `x` and `z=x^2`.\n",
    "\n",
    "    y = a z + b x + c\n",
    "    \n",
    "This is also called multiple linear regression when there is more than one input `y=f(x,z)` where `f` is a function of inputs `x` and `z`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p2 = np.polyfit(x,z,2)\n",
    "print(p2)\n",
    "plt.plot(x,z,'ro',label='Measured (z)')\n",
    "plt.plot(x,np.polyval(p2,x),'b-',label='Predicted (y)')\n",
    "plt.legend(); plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is more information on the regressed coefficients for the quadratic fit if you use `statsmodels`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import statsmodels.api as sm\n",
    "xc = np.vstack((x**2,x,np.ones(n))).T\n",
    "model = sm.OLS(z,xc).fit()\n",
    "predictions = model.predict(xc)\n",
    "model.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![expert](https://apmonitor.com/che263/uploads/Begin_Python/expert.png)\n",
    "\n",
    "### Polynomial Regression Activity\n",
    "\n",
    "Create a polynomial model with the data:\n",
    "\n",
    "```python\n",
    "xr = [0.0,1.0,2.0,3.5,5.0]\n",
    "yr = [1.7,1.45,1.05,0.4,0.2]\n",
    "```\n",
    "\n",
    "Show the polynomial model on a plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![list](https://apmonitor.com/che263/uploads/Begin_Python/list.png)\n",
    "\n",
    "### Nonlinear Regression\n",
    "\n",
    "Nonlinear regression requires a different tool such as `curve_fit` that requires a function `f` that returns a prediction. It also requires the data `x` and `z`. The unknown parameters `a` and `b` are adjusted so that the predicted output matches the measured output `z`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.optimize import curve_fit\n",
    "def f(x,a,b):\n",
    "    return a * np.exp(b*x)\n",
    "p, pcov = curve_fit(f,x,z)\n",
    "print('p = '+str(p))\n",
    "plt.plot(x,z,'ro')\n",
    "plt.plot(x,f(x,*p),'b-')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![expert](https://apmonitor.com/che263/uploads/Begin_Python/expert.png)\n",
    "\n",
    "### Nonlinear Regression Activity\n",
    "\n",
    "Create a nonlinear model with the data:\n",
    "\n",
    "$y = a \\ln\\left( b \\, x \\right)$\n",
    "\n",
    "```python\n",
    "xr = [0.1,1.0,2.0,3.5,5.0]\n",
    "yr = [0.2,0.4,1.05,1.45,1.7]\n",
    "```\n",
    "\n",
    "Show the nonlinear model on a plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![exercise](https://apmonitor.com/che263/uploads/Begin_Python/exercise.png)\n",
    "\n",
    "### Machine Learning\n",
    "\n",
    "Machine learning is computer algorithms and statistical models that rely on patterns and inference. They perform a specific task without explicit instructions. Machine learned regression models can be as simple as linear regression or as complex as deep learning. This tutorial demonstrates several regression methods with `scikit-learn`.\n",
    "\n",
    "#### Generate Data\n",
    "\n",
    "To make the plot interactive, add the command: \n",
    "\n",
    "```python\n",
    "%matplotlib notebook\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "\n",
    "import math\n",
    "def f(x,y):\n",
    "    return 2*math.cos(x)*y + x*math.cos(y) - 3*x*y\n",
    "\n",
    "n = 500\n",
    "x = (np.random.rand(n)-0.5)*2.0\n",
    "y = (np.random.rand(n)-0.5)*2.0\n",
    "z = np.empty_like(x)\n",
    "for i in range(n):\n",
    "    z[i] = f(x[i],y[i])\n",
    "data = pd.DataFrame(np.vstack((x,y,z)).T,columns=['x','y','z'])\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111, projection='3d')\n",
    "ax.scatter(x,y,z,c=z,cmap='plasma')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Scale Data\n",
    "\n",
    "The data can be scaled with a Standard Scalar or just left unscaled because the values of `x`, `y`, and `z` are already close to acceptable ranges. If scaling is desired, it could be done with a few lines of additional code and changes to the plots.\n",
    "\n",
    "```python\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "s = StandardScaler()\n",
    "ds = s.fit_transform(data)\n",
    "ds = pd.DataFrame(ds,columns=data.columns)\n",
    "```\n",
    "\n",
    "The data is split into training and testing sets with `train_test_split`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# no data scaling\n",
    "ds = data\n",
    "\n",
    "# data splitting into train and test sets\n",
    "from sklearn.model_selection import train_test_split\n",
    "train,test = train_test_split(ds, test_size=0.2, shuffle=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Function for Plotting\n",
    "\n",
    "Run this code so that each of the regressor models will train and display on a 3D scatter and surface plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fit(method):\n",
    "    # create points for plotting surface\n",
    "    xp = np.arange(-1, 1, 0.1)\n",
    "    yp = np.arange(-1, 1, 0.1)\n",
    "    XP, YP = np.meshgrid(xp, yp)\n",
    "\n",
    "    model = method.fit(train[['x','y']],train['z'])\n",
    "    zp = method.predict(np.vstack((XP.flatten(),YP.flatten())).T)\n",
    "    ZP = zp.reshape(np.size(XP,0),np.size(XP,1))\n",
    "\n",
    "    r2 = method.score(test[['x','y']],test['z'])\n",
    "    print('R^2: ' + str(r2))\n",
    "\n",
    "    fig = plt.figure()\n",
    "    ax = fig.add_subplot(111, projection='3d')\n",
    "    ax.scatter(ds['x'],ds['y'],ds['z'],c=z,cmap='plasma',label='data')\n",
    "    ax.plot_surface(XP, YP, ZP, cmap='coolwarm',alpha=0.7,\n",
    "                    linewidth=0, antialiased=False)\n",
    "    plt.show()\n",
    "    return"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Linear Regression with `sklearn`\n",
    "\n",
    "The simplest regressor is a linear model. As expected, this model doesn't perform very well with the nonlinear data but it does predict the slope of the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "from sklearn import linear_model\n",
    "lm = linear_model.LinearRegression()\n",
    "fit(lm)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### K-Nearest Neighbors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from sklearn.neighbors import KNeighborsRegressor\n",
    "knn = KNeighborsRegressor(n_neighbors=2)\n",
    "fit(knn)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Support Vector Regressor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from sklearn import svm\n",
    "s = svm.SVR(gamma='scale')\n",
    "fit(s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Multilayer Perceptron (Neural Network)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.neural_network import MLPRegressor\n",
    "nn = MLPRegressor(hidden_layer_sizes=(3), \n",
    "                  activation='tanh', solver='lbfgs')\n",
    "fit(nn)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![expert](https://apmonitor.com/che263/uploads/Begin_Python/expert.png)\n",
    "\n",
    "### Regressor Activity\n",
    "\n",
    "Find another [regressor in scikit-learn](https://scikit-learn.org/stable/supervised_learning.html) such as *Decision Tree Regressor* or *Passive Agressive Regressor*. Train and test the regressor with the `fit()` function in this notebook by passing in the regressor object such as:\n",
    "\n",
    "*Decision Tree Regressor*\n",
    "\n",
    "```python\n",
    "from sklearn import tree\n",
    "dt = tree.DecisionTreeRegressor()\n",
    "fit(dt)\n",
    "```\n",
    "\n",
    "*Passive Aggressive Regressor*\n",
    "\n",
    "```python\n",
    "from sklearn.linear_model import PassiveAggressiveRegressor\n",
    "par = PassiveAggressiveRegressor(max_iter=2000,tol=1e-3)\n",
    "fit(par)\n",
    "```\n",
    "\n",
    "Change the options of the regressor if you can improve the $R^2$ value such as the `PassiveAggressiveRegressor` with `max_iter=2000`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![gekko](https://apmonitor.com/che263/uploads/Begin_Python/gekko.png)\n",
    "\n",
    "### Deep Learning\n",
    "\n",
    "Deep learning is regression with a neural network with multiple layers. Regression with deep learning has specialized packages such as [Tensorflow that are built for large-scale data](https://www.tensorflow.org) or [Gekko that are built for configurable model structures](https://gekko.readthedocs.io/en/latest/). Below is one of the examples from a [deep learning tutorial with Gekko](https://apmonitor.com/do/index.php/Main/DeepLearning).  This same example with Keras (Tensorflow) is also shown in the link."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from gekko import brain\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n", 
    "\n",
    "x = np.linspace(0.0,2*np.pi)\n",
    "y = np.sin(x)\n",
    "\n",
    "b = brain.Brain(remote=False)\n",
    "b.input_layer(1)\n",
    "b.layer(linear=2)\n",
    "b.layer(tanh=2)\n",
    "b.layer(linear=2)\n",
    "b.output_layer(1)\n",
    "b.learn(x,y,disp=False)      \n",
    "\n",
    "xp = np.linspace(-2*np.pi,4*np.pi,100)\n",
    "yp = b.think(xp)  \n",
    "\n",
    "plt.figure()\n",
    "plt.plot(x,y,'bo',label='meas (label)')\n",
    "plt.plot(xp,yp[0],'r-',label='model')\n",
    "plt.legend(); plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![expert](https://apmonitor.com/che263/uploads/Begin_Python/expert.png)\n",
    "\n",
    "### TCLab Activity\n",
    "\n",
    "### Record Temperatures\n",
    "\n",
    "![connections](https://apmonitor.com/che263/uploads/Begin_Python/connections.png)\n",
    "\n",
    "Set heater 1 to 80% with `lab.Q1(80)` and heater 2 to 60% with `lab.Q2(60)`. Record the temperatures (`lab.T1` and `lab.T2`) every 0.5 seconds (use `time.sleep(0.5)`) for 30 seconds. Store the values for time, temperature 1, and temperature 2 in `numpy` arrays. Use __time.time()__ to get the current time in seconds."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![expert](https://apmonitor.com/che263/uploads/Begin_Python/expert.png)\n",
    "\n",
    "### Linear Regression\n",
    "\n",
    "Create a linear model for `T2` with regression. Report the $R^2$ value. Add the regression line as a black line to a plot with the measured `T2` as blue circles. Add appropriate labels to the plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![expert](https://apmonitor.com/che263/uploads/Begin_Python/expert.png)\n",
    "\n",
    "### Nonlinear Regression\n",
    "\n",
    "Create a nonlinear regression for `T1`. Fit the $T_1$ data with:\n",
    "\n",
    "1. a nonlinear model of the form $T_1 = a + b \\exp{(c \\, t)}$ where `a`, `b`, and `c` are adjustable parameters.\n",
    "2. a nonlinear model using a regressor in `scikit-learn`, `keras (tensorflow)`, or `gekko`\n",
    "\n",
    "Report the $R^2$ value for each."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}