regressionExample.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy as sp\n",
    "from scipy import stats\n",
    "from scipy.optimize import curve_fit\n",
    "from scipy.optimize import minimize\n",
    "from scipy.interpolate import *\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### import data "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "url = 'http://apmonitor.com/che263/uploads/Main/stats_data.txt'\n",
    "data = pd.read_csv(url)\n",
    "del data['index']\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = data['x'].values\n",
    "y = data['y'].values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fitting a line to the data. Here we actually get a bunch of useful information \n",
    "### out of the function, but for now we care about the slope and y-intercept"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plotting the line fit\n",
    "plt.plot(x,y,'.') # plot the original data points\n",
    "plt.plot(x,fit,'-') # plot the new fitted data\n",
    "plt.vlines(x,fit,y,color=(0.7, 0.7, 0.7)) # Plot the errors\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('y')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p1 = np.polyfit(x,y,1)\n",
    "p2 = np.polyfit(x,y,2)\n",
    "p3 = np.polyfit(x,y,3)\n",
    "print(p1)\n",
    "print(p2)\n",
    "print(p3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(x,y,'o')\n",
    "xp = np.linspace(16,22,100)\n",
    "plt.plot(xp,np.polyval(p1,xp),'r-')\n",
    "plt.plot(xp,np.polyval(p2,xp),'b--')\n",
    "plt.plot(xp,np.polyval(p3,xp),'m:')\n",
    "yfit = p1[0] * x + p1[1]\n",
    "yresid= y - yfit\n",
    "SSresid = sum(pow(yresid,2))\n",
    "SStotal = len(y) * np.var(y)\n",
    "rsq = 1 - SSresid/SStotal\n",
    "print(yfit)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(rsq)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "slope,intercept,r_value,p_value,std_err = linregress(x,y)\n",
    "print(pow(r_value,2))\n",
    "print(p_value)\n",
    "#show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### New Example http://apmonitor.com/che263/index.php/Main/PythonDataRegression\n",
    "### load data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xm = np.array([18.3447,79.86538,85.09788,10.5211,44.4556, \\\n",
    "               69.567,8.960,86.197,66.857,16.875, \\\n",
    "               52.2697,93.917,24.35,5.118,25.126, \\\n",
    "               34.037,61.4445,42.704,39.531,29.988])\n",
    "ym = np.array([5.072,7.1588,7.263,4.255,6.282, \\\n",
    "               6.9118,4.044,7.2595,6.898,4.8744, \\\n",
    "               6.5179,7.3434,5.4316,3.38,5.464, \\\n",
    "               5.90,6.80,6.193,6.070,5.737])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Preform regression "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate y\n",
    "def calc_y(x):\n",
    "    a,b,c = x\n",
    "    y = a + b/xm + c*np.log(xm)\n",
    "    return y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define objective\n",
    "def objective(x):\n",
    "    return np.sum(((calc_y(x)-ym)/ym)**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# initial guesses\n",
    "x0 = np.zeros(3)\n",
    "\n",
    "# show initial objective\n",
    "print('Initial SSE Objective: ' + str(objective(x0)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# optimize\n",
    "# bounds on variables\n",
    "bnds100 = (-100.0, 100.0)\n",
    "no_bnds = (-1.0e10, 1.0e10)\n",
    "bnds = (no_bnds, no_bnds, bnds100)\n",
    "solution = minimize(objective,x0,method='SLSQP',bounds=bnds)\n",
    "x = solution.x\n",
    "y = calc_y(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# show final objective\n",
    "print('Final SSE Objective: ' + str(objective(x)))\n",
    "\n",
    "# print solution\n",
    "print('Solution')\n",
    "print('a = ' + str(x[0]))\n",
    "print('b = ' + str(x[1]))\n",
    "print('c = ' + str(x[2]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot solution\n",
    "import matplotlib.pyplot as plt\n",
    "plt.figure(1)\n",
    "plt.plot(xm,ym,'ro')\n",
    "plt.plot(xm,y,'bx');\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('y')\n",
    "plt.legend(['Measured','Predicted'],loc='best')\n",
    "plt.savefig('results.png')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}