3.Polynomial Regression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Polynomial Regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A linear graph of averages \n",
    "Recall the graph of averages from the Simple Regression notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%pylab inline\n",
    "import pandas as pd\n",
    "from numpy import arange,array,ones,linalg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "HW=pd.read_csv('data/HW25000.csv')\n",
    "HW=HW.iloc[:,1:]\n",
    "HW.columns=['Height','Weight']\n",
    "HW.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A=np.array(HW['Height'])\n",
    "A=np.array([ones(len(A)),A])\n",
    "y=np.array(HW['Weight'])\n",
    "\n",
    "w1 = linalg.lstsq(A.T,y)[0] # finding the optimal parameters\n",
    "print('w1=',w1)\n",
    "def f(x,w):\n",
    "    return w[0]+w[1]*x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Average_Plot(df,x_name,y_name):\n",
    "    # calculate the mean weight for each 1-inch interval of height.\n",
    "    df['round_'+x_name]=df[x_name].round()\n",
    "\n",
    "    per_height_means=df.groupby('round_'+x_name).mean()[[y_name]]\n",
    "\n",
    "    ax=df.plot(kind='scatter',s=1,x=x_name,y=y_name,figsize=[10,8]);\n",
    "    per_height_means.plot(y=y_name,style='ro',ax=ax,legend=False)\n",
    "\n",
    "    _xlim=xlim()\n",
    "    _ylim=ylim()\n",
    "    for _x in arange(_xlim[0]+0.5,_xlim[1],1):\n",
    "        ax.plot([_x,_x],[_ylim[0],_ylim[1]],'g')\n",
    "    return ax\n",
    "ax=Average_Plot(HW,'Height','Weight');\n",
    "x0,x1=plt.xlim()\n",
    "ax.plot([x0,x1],[f(x0,w1),f(x1,w1)],'k');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## non-linear graph of averages\n",
    "Suppose we have data for which the averages do not fall close to a straight line."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "HW['P2']=HW['Weight']+(HW['Height']-68)**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax=Average_Plot(HW,'Height','P2');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The limits of linear regression\n",
    "We **can** use linear regression for this data, but the fit will be poor."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A=np.array(HW['Height'])\n",
    "A=np.array([ones(len(A)),A])\n",
    "y=np.array(HW['P2'])\n",
    "\n",
    "w1 = linalg.lstsq(A.T,y)[0] # finding the optimal parameters\n",
    "print('w1=',w1)\n",
    "def f(x,w):\n",
    "    return w[0]+w[1]*x\n",
    "ax=Average_Plot(HW,'Height','P2');\n",
    "x0,x1=plt.xlim()\n",
    "ax.plot([x0,x1],[f(x0,w1),f(x1,w1)],'k');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Lets try a second degree polynomial\n",
    "\n",
    "A second degree polynomial can be expressed in the form \n",
    "$$f_2(x)=w_0+w_1 x + w_2 x^2 $$\n",
    "\n",
    "We find the parameters $w_0,w_1,w_2$ that minimize the root-mean-square error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A=np.array(HW['Height'])\n",
    "A=np.array([ones(len(A)),A,A**2])\n",
    "y=np.array(HW['P2'])\n",
    "\n",
    "w2 = linalg.lstsq(A.T,y)[0] # finding the optimal parameters\n",
    "print('w2=',w2)\n",
    "def f2(x,w):\n",
    "    return w[0]+w[1]*x+w[2]*x**2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The second degree polynomial is a pretty good fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax=Average_Plot(HW,'Height','P2')\n",
    "_xmin,_xmax=xlim()\n",
    "_xrange=_xmax-_xmin\n",
    "X=arange(_xmin,_xmax,_xrange/100.)\n",
    "Y=f2(X,w2)\n",
    "ax.plot(X,f2(X,w2),'k');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Over-fitting, Under-fitting and Model selection\n",
    "Often we are given data an we don't know what degree polynomial "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# A function for computing a polynomial of an arbitrary degree\n",
    "# The degree is determined by the length of the parameter vector w\n",
    "def F(X,w):\n",
    "    Accum=w[0]*np.ones(len(X))\n",
    "    for i in range(1,len(w)):\n",
    "        Accum+=w[i]*X**i\n",
    "    return Accum\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(0)\n",
    "X=np.arange(-1,1.6,0.25)\n",
    "Y=X+np.random.rand(len(X))\n",
    "Y\n",
    "\n",
    "Data=pd.DataFrame({'x':X,\n",
    "                   'y':Y})\n",
    "ax=Data.plot(kind='scatter',s=30,c='r',x='x',y='y',figsize=[6,5]);\n",
    "grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "L=len(X)\n",
    "split=[0,1]*L\n",
    "Data['split']=split[:L]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def polyfit(ax,d=8):\n",
    "    Train=Data[Data['split']==1]\n",
    "    Test=Data[Data['split']==0]\n",
    "    \n",
    "    A=np.array([Train['x']])\n",
    "    D=np.ones([1,A.shape[1]])\n",
    "    for i in range(1,d+1):\n",
    "        D=np.concatenate([D,A**i])\n",
    "    w=linalg.lstsq(D.T,Train['y'])[0]\n",
    "    train_RMS=sqrt(mean((Train['y']-F(Train['x'],w))**2))\n",
    "    test_RMS=sqrt(mean((Test['y']-F(Test['x'],w))**2))\n",
    "\n",
    "    Train.plot(kind='scatter',s=30,c='b',x='x',y='y',ax=ax,label='Train');\n",
    "    Test.plot(kind='scatter',s=30,c='r',x='x',y='y',ax=ax,label='Test');\n",
    "    grid()\n",
    "    legend()\n",
    "    _xmin,_xmax=plt.xlim()\n",
    "    _xrange=_xmax-_xmin\n",
    "    X=arange(_xmin,_xmax,_xrange/100.)\n",
    "    ax.plot(X,F(X,w),'k');\n",
    "    plt.title('d=%d, train_RMS=%5.3f, test_RMS=%5.3f'%(d,train_RMS,test_RMS))\n",
    "    return train_RMS,test_RMS\n",
    "fig=plt.figure(figsize=[5,5])\n",
    "ax=plt.subplot(111)\n",
    "polyfit(ax,3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rows=2; cols=3; max_d=6\n",
    "fig=plt.figure(figsize=[14,10])\n",
    "train_RMS=np.zeros(max_d)\n",
    "test_RMS=np.zeros(max_d)\n",
    "for d in range(max_d):\n",
    "    if d==0:\n",
    "        ax=plt.subplot(rows,cols,d+1)\n",
    "        ax0=ax\n",
    "    else:\n",
    "        ax=plt.subplot(rows,cols,d+1,sharex=ax0)\n",
    "    train_RMS[d],test_RMS[d]=polyfit(ax,d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "plot(train_RMS,label='train RMS')\n",
    "plot(test_RMS,label='test_RMS')\n",
    "legend()\n",
    "grid()"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Create Assignment",
  "kernelspec": {
   "display_name": "Python [default]",
   "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.2"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {
    "height": "67px",
    "width": "252px"
   },
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {
    "height": "512px",
    "left": "0px",
    "right": "921px",
    "top": "109px",
    "width": "78px"
   },
   "toc_section_display": "block",
   "toc_window_display": false
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}