classification.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# MSAIL Tutorial Series"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "## Classification"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Kaggle Intro and Outline\n",
    "\n",
    "We will be working through the following [Kaggle competition](https://www.kaggle.com/c/titanic) to learn about various classification algorithms."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "First we import the necessary libraries:\n",
    "* pandas for parsing the csv file\n",
    "* numpy for math and arrays\n",
    "* scipy for more scientific computing\n",
    "* matplotlib for generating visualizations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import scipy as sp\n",
    "from scipy import special\n",
    "\n",
    "%matplotlib inline\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Read the training data into a pandas dataframe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "PassengerId      int64\n",
      "Survived         int64\n",
      "Pclass           int64\n",
      "Name            object\n",
      "Sex             object\n",
      "Age            float64\n",
      "SibSp            int64\n",
      "Parch            int64\n",
      "Ticket          object\n",
      "Fare           float64\n",
      "Cabin           object\n",
      "Embarked        object\n",
      "dtype: object\n"
     ]
    }
   ],
   "source": [
    "df = pd.read_csv('train.csv', header=0)\n",
    "print df.dtypes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Assign numerical values for gender"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "df['Gender'] = df['Sex'].map( {'female': 0, 'male': 1} ).astype(int)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Clean the age data by replacing NaN with median given gender/class"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "median_ages = np.zeros((2, 3))\n",
    "for i in range(0, 2):\n",
    "    for j in range (0, 3):\n",
    "        median_ages[i, j] = df[(df['Gender'] == i) & (df['Pclass'] == j+1)]['Age'].dropna().median()\n",
    "\n",
    "df['AgeFill'] = df['Age']\n",
    "        \n",
    "for i in range(0, 2):\n",
    "    for j in range(0, 3):\n",
    "        df.loc[(df.Age.isnull()) & (df.Gender == i) & (df.Pclass == j+1), 'AgeFill'] = median_ages[i, j]\n",
    "\n",
    "df['AgeIsNull'] = pd.isnull(df.Age).astype(int)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Construct the data matrix, x, using class, gender, and age, and ticket price. And split the data into train and test sets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "x = np.array(df[['Pclass', 'Gender', 'AgeFill', 'Fare']])\n",
    "y = np.array(df['Survived'])\n",
    "\n",
    "num_train = 850\n",
    "\n",
    "xtrain = x[0:num_train]\n",
    "ytrain = y[0:num_train]\n",
    "\n",
    "xtest = x[num_train:]\n",
    "ytest = y[num_train:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Define a function that:\n",
    "* takes a model\n",
    "* trains it\n",
    "* uses it to perform classification on the test set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def train_and_test(model, xtrain, ytrain, xtest, ytest):\n",
    "    model = model.fit(xtrain, ytrain)\n",
    "    output = model.predict(xtest)\n",
    "    print (1.0 * sum([i==j for (i, j) in zip(output, ytest)])) / len(ytest)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Define some classification models and get some results!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.780487804878\n"
     ]
    }
   ],
   "source": [
    "import sklearn.linear_model\n",
    "model = sklearn.linear_model.LogisticRegression()\n",
    "train_and_test(model, xtrain, ytrain, xtest, ytest)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Logistic Regression\n",
    "\n",
    "Logistic regression is used for binary classification tasks. It processes a dataset $D = \\{(x_1,t_1),\\dots,(x_n,t_n)\\}$, where $t_i\\in\\{0,1\\}$ and the feature vector of the $i$th example is $\\phi(x_i)\\in\\mathbb{R}^M$.\n",
    "\n",
    "Logistic regression forms a probabilistic model. It estimates probability distributions of the two classes $p(t=1|x;w)$ and $p(t=0|x;w)$. It fits its paramaters $w\\in\\mathbb{R}^M$ to the training data by MLE, ie finds the $w$ that maximizes the probability of the training data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "__Logistic Model__\n",
    "\n",
    "Given features for an example, $\\phi(x)$, logistic regression models the probability of this example belonging to the class 1 as:\n",
    "\n",
    "$$p(t=1|x;w) = \\sigma(w^T\\phi(x))$$\n",
    "\n",
    "and defines the probability of the example belonging to the class 0 as:\n",
    "\n",
    "$$\n",
    "\\begin{eqnarray}\n",
    "p(t=0|x;w) &=& 1 - p(t=1|x;w) \\nonumber\\\\\n",
    "&=& 1 - \\sigma(w^T\\phi(x)) \\nonumber\n",
    "\\end{eqnarray}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "$\\sigma(a)$ is the sigmoid function, which is defined as:\n",
    "\n",
    "$$\\sigma(a) = \\frac{1}{1+e^{-a}}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x10825b950>]"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXUAAAEACAYAAABMEua6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGMZJREFUeJzt3XuYlHXdx/H31wWJkMTCkKOUIqihj5qIZToe0g1MxEri\n6eQhoauwg+Vj5HXV9OiTj6XkY54whZQMPASyBAmojRLmAURBZAVUhAU8g6DJ+fv8cQ86jrs7M7uz\n85u55/O6rvvavWfuHb+u44cf39/v/o25OyIiEg97hC5ARESKR6EuIhIjCnURkRhRqIuIxIhCXUQk\nRhTqIiIxkjPUzWyCmb1iZkuaueZaM1thZk+b2RHFLVFERPKVz0h9IlDb1JNmNgQ40N37AaOAG4tU\nm4iIFChnqLv7PGBDM5ecAdyWvvYxoIuZdStOeSIiUohi9NR7AmsyzhuAXkV4XRERKVCxJkot61x7\nD4iIBNCuCK+xFuidcd4r/dgHmJmCXkSkBdw9e+DcpGKEeh0wBphiZoOBje7+ShOFFeEfJwDJZJJk\nMhm6jFjQ77J57rBmDdTXw4oVsHw5PP88vPQSrF4Nu3ZBz56w337QvTu8+GKSL385yb77QteusM8+\n7x977w2dO0NNTeh/q8phlneeA3mEuplNBk4AuprZGuBXQHsAdx/v7rPMbIiZrQTeAc4tuGoRKQu7\ndkWh/dhj8MQT8PTTsGQJdOwIBx8MBx0E/frBSSfB/vtHR5cukJk7ySRcemmwf4WqlzPU3X1kHteM\nKU45IlJK7vDMM/Dgg/DAAzBvXhTSxxwTHV/5CgwcGI24pTIUo/0iASQSidAlxEa1/S63bYOHHoJ7\n74Xp06FDBzj5ZPjGN+Dmm6M2SmtU2++z3Fip+txm5uqpi4ThDk8+CRMnwpQpcOCBMHw4DBsGAwaE\nrk6aY2YlnygVkTK1ZQtMmgTXXQebNsG558LChVEvXOJJoS4SQxs2wA03RGF+5JEwbhyceCLsoS38\nYk//iUViZMsWuOqqaJXK8uUwdy7MnBn1zBXo1UEjdZEYcIc774SxY+Gww+Dhh6MliFJ9FOoiFa6h\nAUaNgnXr4E9/ghNOCF2RhKS/kIlUKHe45RY44gg49tjoZiEFumikLlKBNm+OVrK8+GJ049DAgaEr\nknKhkbpIhamvh0GD4BOfgEceUaDLBynURSpIXR0cfzz87Gcwfnx0N6hIJrVfRCrExInRRlkzZ8LR\nR4euRsqVQl2kAlxzDfz+9/CPf0D//qGrkXKmUBcpY+7w61/D5MnRDop9+oSuSMqdQl2kjF11Fdx9\nd3QzUTd9nLvkQaEuUqZuvz3au2X+fAW65E+hLlKGZs2Ciy+GVAp69QpdjVQShbpImVmwAL7znWj5\novZvkUJpnbpIGXnjDfjqV+Gmm6Jb/0UKpU8+EikTu3bB0KFw6KHRBKkIFP7JRxqpi5SJyy+Ht9+G\nK64IXYlUMvXURcrAnDlRy2XBAmjfPnQ1UsnUfhEJ7M03o025Jk2Ck04KXY2Um0LbLwp1kcC+/W3o\n0gWuvTZ0JVKOCg11tV9EAvrb3+Cf/4QlS0JXInGhUBcJZONG+N73orZLp06hq5G4UPtFJJDzzoOO\nHeH660NXIuVM7ReRCjBvHtx/PyxdGroSiRutUxcpsV274Mc/hiuvhM6dQ1cjcaNQFymxP/0JPvIR\n+PrXQ1cicaSeukgJbdoEAwbA9On6SDrJj7YJECljv/kNnHqqAl3ajkbqIiXywgswaBAsXgw9eoSu\nRiqFRuoiZSqZhB/+UIEubUtLGkVKYPly+PvfYeXK0JVI3GmkLlICl10GP/oR7L136Eok7nKGupnV\nmlm9ma0ws0saeb6rmd1nZk+Z2TNmdk6bVCpSoZ57DmbPjlovIm2t2YlSM6sBngNOAdYCTwAj3X1Z\nxjVJoIO7jzWzrunru7n7jqzX0kSpVKVvfhMOOQR+8YvQlUglKvZE6SBgpbuvcvftwBRgWNY164GP\npb//GPBGdqCLVKv6+ugDMC68MHQlUi1yTZT2BNZknDcAx2Rd80fgQTNbB3QGzi5eeSKV7fLL4Sc/\n0XYAUjq5Qj2ffskvgKfcPWFmBwBzzexwd9+cfWEymXzv+0QiQSKRKKBUkcqyZk204uWGG0JXIpUk\nlUqRSqVa/PO5euqDgaS716bPxwK73P3KjGtmAf/j7vPT5w8Al7j7gqzXUk9dqsrFF0ebd119dehK\npJIVe+vdBUA/M+sLrANGACOzrqknmkidb2bdgP7AC/kWIBJHmzfDhAnw5JOhK5Fq02you/sOMxsD\nzAZqgFvdfZmZjU4/Px74DTDRzJ4mmnj9L3d/s43rFilrEybAKafA/vuHrkSqjfZ+ESmynTvhwANh\nyhQ4JntZgUiBtPeLSGD33hvt76JAlxAU6iJFNm4cXHRR6CqkWinURYpo4UJYuxbOPDN0JVKtFOoi\nRXTzzXDBBVBTE7oSqVaaKBUpks2boU8fWLpUe6ZL8WiiVCSQyZPhxBMV6BKWQl2kSMaPh1GjQlch\n1U6hLlIECxfCG29EHyotEpJCXaQIxo+PJkj30P9REpgmSkVaafcE6bPPQvfuoauRuNFEqUiJTZkC\niYQCXcqDQl2klW67Dc49N3QVIhG1X0Ra4fnn4dhjoaEB9twzdDUSR2q/iJTQn/8MI0Yo0KV85PqQ\nDBFpgjtMmhTddCRSLjRSF2mhf/0L2rWDz342dCUi71Ooi7TQpEnwrW+B5d3tFGl7migVaYGtW6M9\nXp58Uh9ZJ21LE6UiJTBzJgwcqECX8qNQF2mB3a0XkXKj9otIgd56C3r3htWroUuX0NVI3Kn9ItLG\npk+PtgVQoEs5UqiLFOiuu6IbjkTKkdovIgXYsAH69o22BejcOXQ1Ug3UfhFpQ9OmwcknK9ClfCnU\nRQqg1ouUO7VfRPL0+utwwAGwbh106hS6GqkWar+ItJFp06C2VoEu5U2hLpKnO++Es88OXYVI89R+\nEcnDq6/CQQfB+vXQsWPoaqSaqP0i0gbq6qLWiwJdyp1CXSQPU6fC8OGhqxDJTe0XkRx27/Wydq3W\np0vpqf0iUmQzZ8IJJyjQpTIo1EVymDZNrRepHDlD3cxqzazezFaY2SVNXJMws0Vm9oyZpYpepUgg\n774Lc+bAGWeErkQkP+2ae9LMaoDrgFOAtcATZlbn7ssyrukCXA+c5u4NZta1LQsWKaU5c+Coo6Cr\n3tVSIXKN1AcBK919lbtvB6YAw7Ku+U/gr+7eAODurxe/TJEwpk2Ds84KXYVI/nKFek9gTcZ5Q/qx\nTP2Aj5vZP8xsgZnpQ74kFrZvhxkz4MwzQ1cikr9m2y9APmsQ2wNHAicDHwX+ZWaPuvuK1hYnEtJD\nD8GBB0KvXqErEclfrlBfC/TOOO9NNFrPtAZ43d3fBd41s4eBw4EPhXoymXzv+0QiQSKRKLxikRK5\n916tepHSS6VSpFKpFv98szcfmVk74DmiUfg64HFgZNZE6QCiydTTgA7AY8AId38267V085FUDHfY\nf3+YPRsOPjh0NVLNCr35qNmRurvvMLMxwGygBrjV3ZeZ2ej08+Pdvd7M7gMWA7uAP2YHukilWbQI\nOnSAAQNCVyJSGG0TINKIZBLefhuuuip0JVLttE2ASBFMnw7DshfvilQAhbpIlpdegoYGOPbY0JWI\nFE6hLpJlxgwYOhTa5VobJlKGFOoiWaZP114vUrk0USqSYeNG6NMH1q2DvfYKXY2IJkpFWuW+++D4\n4xXoUrkU6iIZ1HqRSqf2i0ja9u3wyU/C0qXQo0foakQiar+ItNC8edEGXgp0qWQKdZG0GTPUepHK\np1AXIdrAa8YM+PKXQ1ci0joKdRGgvh62bYPDDw9diUjrKNRFiEbpp58Olvd0lEh5UqiLoNaLxIeW\nNErVe+MN+PSn4ZVX4CMfCV2NyAdpSaNIgWbNgpNOUqBLPCjUpeqp9SJxovaLVLVt26Bbt2j1S7du\noasR+TC1X0QK8PDD0L+/Al3iQ6EuVU2tF4kbhbpULd1FKnGkUJeq9eyzsHMnDBwYuhKR4lGoS9Wq\nq4tG6bqLVOJEoS5VS60XiSMtaZSq9Oqr0K9f9LVDh9DViDRNSxpF8jBrFpxyigJd4kehLlVJrReJ\nK7VfpOps3Rp9FumKFdFXkXKm9otIDqkUHHqoAl3iSaEuVUetF4mzdqELECkl92h9+n33ha5EpG1o\npC5VZdGiaMXLwQeHrkSkbSjUparU1cEZZ+guUokvhbpUlenTYdiw0FWItB0taZSqsXo1HHkkvPwy\ntNNsklQILWkUaUJdHQwdqkCXeMsZ6mZWa2b1ZrbCzC5p5rqjzWyHmZ1V3BJFikOtF6kGzbZfzKwG\neA44BVgLPAGMdPdljVw3F/g3MNHd/9rIa6n9IsG89Rb06gXr18Nee4WuRiR/xW6/DAJWuvsqd98O\nTAEaG+tcCNwDvJZ3pSIl9Pe/w/HHK9Al/nKFek9gTcZ5Q/qx95hZT6KgvzH9kIbjUnZ2L2UUibtc\noZ5PQF8D/DzdW7H0IVI2tm2LRuraGkCqQa51AGuB3hnnvYlG65mOAqZYdDdHV+BLZrbd3euyXyyZ\nTL73fSKRIJFIFF6xSIEefDC6g7RHj9CViOSWSqVIpVIt/vlcE6XtiCZKTwbWAY/TyERpxvUTgRnu\nPrWR5zRRKkGMGgUHHQQ/+1noSkQKV+hEabMjdXffYWZjgNlADXCruy8zs9Hp58e3qlqRNrZzZ7SU\n8ZFHQlciUhq6o1Ribd48uPBCeOqp0JWItIzuKBXJMHUqDB8eugqR0tEN0xJb7jBtWvShGCLVQiN1\nia1Fi6B9e/jMZ0JXIlI6CnWJralT4ayztHe6VBeFusTWtGlRqItUE4W6xFJ9PWzcCEcfHboSkdJS\nqEss3XUXfPWrsIfe4VJl9JaXWLrzThgxInQVIqWnUJfYWboUNm2CwYNDVyJSegp1iZ0774Szz1br\nRaqT3vYSK+5qvUh1U6hLrCxeHO2frlUvUq0U6hIru1svuuFIqpX2fpHY2N16ufvu0JWIhKORusTG\nwoXRCP2II0JXIhKOQl1iY/cEqVovUs30IRkSCzt3Qp8+MHcuHHJI6GpEikcfkiFV6cEHoXt3BbqI\nQl1i4fbb4VvfCl2FSHhqv0jFe/tt6NULli+HT34ydDUixaX2i1SdqVPhuOMU6CKgUJcYmDQJvv3t\n0FWIlAe1X6SirV0LAwdGXzt2DF2NSPGp/SJV5Y474CtfUaCL7KZQl4rlrlUvItkU6lKxHn0Utm6F\nL3whdCUi5UOhLhXr5pth1ChtCyCSSROlUpE2boS+fWHFCth339DViLQdTZRKVfjzn6G2VoEukk2h\nLhXHHcaPh9GjQ1ciUn4U6lJxdk+QJhKhKxEpPwp1qTjjx2uCVKQpmiiVirJhA3zqU5ogleqhiVKJ\ntQkTYOhQBbpIUzRSl4qxYwcccAD89a/w2c+GrkakNDRSl9i6555obboCXaRpeYW6mdWaWb2ZrTCz\nSxp5/htm9rSZLTaz+WZ2WPFLlWrmDldfDRddFLoSkfKWM9TNrAa4DqgFDgFGmtnBWZe9ABzv7ocB\nlwE3F7tQqW7z50d3kZ5+euhKRMpbPiP1QcBKd1/l7tuBKcCwzAvc/V/u/lb69DGgV3HLlGo3bhz8\n+MdQUxO6EpHylk+o9wTWZJw3pB9ryvnArNYUJZLp+efh4YfhnHNCVyJS/trlcU3eS1bM7ETgPODz\njT2fTCbf+z6RSJDQLYGSh3Hj4LvfhU6dQlci0vZSqRSpVKrFP59zSaOZDQaS7l6bPh8L7HL3K7Ou\nOwyYCtS6+8pGXkdLGqVgDQ1w2GGwbBl06xa6GpHSa4sljQuAfmbW18z2BEYAdVn/0D5Egf7NxgJd\npKWuuALOP1+BLpKvnO0Xd99hZmOA2UANcKu7LzOz0ennxwO/BPYBbrRoQ47t7j6o7cqWarBmDUye\nDPX1oSsRqRy6o1TK1ve/D3vtBb/9behKRMIptP2iUJeytHo1HHFENErXPi9SzbRNgMTCFVfABRco\n0EUKlc+SRpGSWr4c7r5bvXSRltBIXcrOT38Kl1wCXbuGrkSk8mikLmVlzpxohH7PPaErEalMGqlL\n2dixA37yE7jqKujQIXQ1IpVJoS5l46abYL/94IwzQlciUrm0pFHKwptvwoAB8MADMHBg6GpEyofW\nqUtFOucc6NwZ/vCH0JWIlJdCQ10TpRLcrFnw0EOwZEnoSkQqn0JdgnrrLRg9Gm67LdoSQERaR+0X\nCeq734X27eHGG0NXIlKe1H6RijF7Ntx/v9ouIsWkUJcg1q+H886D22+PJkhFpDi0Tl1Kbvt2GDEi\n6qWffHLoakTiRT11KbmLL4ZnnoGZM2EPDStEmqWeupS1qVOjHRgXLlSgi7QFhbqUzNNPw/e+F43Q\nP/GJ0NWIxJPGSlISL74IQ4fC9dfD0UeHrkYkvhTq0uZeew1OOw3GjoWvfS10NSLxplCXNrV5MwwZ\nEq12+cEPQlcjEn9a/SJt5s03o5bL4YdHd4xa3vP3IrKbPnhaysLLL0MiAZ//vAJdpJQU6lJ0q1bB\nF74QtVx+9zsFukgpKdSlqB56CD73OfjRj+DSSxXoIqWmdepSFO5wzTVw5ZUwaRJ88YuhKxKpTgp1\nabWNG6ObilasgEcfhb59Q1ckUr3UfpFW+dvf4DOfgY9/HP75TwW6SGgaqUuLvP46XHRRFOSTJsGJ\nJ4auSERAI3Up0LvvRn3zAQOi0fnixQp0kXKikbrkZds2+Mtf4Fe/giOPhPnzoX//0FWJSDaFujRr\n0yb44x+jlS39+8Mdd8Bxx4WuSkSaolCXD3GHJ56ACROivc9PPRWmT49G6CJS3hTqAkRBvmwZTJsW\ntVm2boVzz4WnnoLevUNXJyL50oZeVeztt2HePHjgAairgy1b4Mwzo+1xjztOd4OKlINCN/TKGepm\nVgtcA9QAt7j7lY1ccy3wJeDfwDnuvqiRaxTqAblDQwM89hg8/jg88kg0Cj/qKDjpJDj99Ki9oiAX\nKS9FDXUzqwGeA04B1gJPACPdfVnGNUOAMe4+xMyOAf7P3Qc38loK9SJKpVIkEokPPe4e7ZC4fHl0\nh+czz8CSJdEBcMwxMGgQDB4c7aD40Y+Wtu5y1NTvUlpGv8/iKvYHTw8CVrr7qvSLTwGGAcsyrjkD\nuA3A3R8zsy5m1s3dXymocslp69bopp/XXoNbbknxwgsJ1q+PRuCrV8NLL0VHx47Qr190HHoo1NbC\nwIHQo4dG4o1RCBWXfp9h5Qr1nsCajPMG4Jg8rukFVFWou8OOHbB9e3Rs2/b+sXXr+8eWLdENPO++\nC//+d3S8807U337nneiTgjZtio6NG6Njw4boAye2bIGuXaPjnXegpga6d4+Ce8gQ6NMH9t8funQJ\n/dsQkVByhXq+/ZLs8V+jPzdkSPrJjGeb+r6p5zK/NvVY9nPZx65dTX/NPHbu/PDX3ceOHR88du6M\nQrZ9+/ePDh2io337aPS8+7xjx6jt0bFjdOy1F3TqFB377Qcf+xh07hyF8z77vH/svff7I+1kMjpE\nRDLl6qkPBpLuXps+HwvsypwsNbObgJS7T0mf1wMnZLdfzEwNdRGRFihmT30B0M/M+gLrgBHAyKxr\n6oAxwJT0HwIbG+unF1KUiIi0TLOh7u47zGwMMJtoSeOt7r7MzEannx/v7rPMbIiZrQTeAc5t86pF\nRKRRJbv5SERE2l6bbr1rZl8zs6VmttPMjsx6bqyZrTCzejM7tS3riCMzS5pZg5ktSh+1oWuqRGZW\nm34PrjCzS0LXU+nMbJWZLU6/Jx8PXU8lMbMJZvaKmS3JeOzjZjbXzJab2Rwzy7m2ra33U18CDAce\nznzQzA4h6s8fAtQCN5iZ9nYvjAPj3P2I9HFf6IIqTfrmuuuI3oOHACPN7OCwVVU8BxLp9+Sg0MVU\nmIlE78VMPwfmuvtBwAPp82a1aZC6e727L2/kqWHAZHffnr6xaSXRjU5SGE0+t857N9e5+3Zg9811\n0jp6X7aAu88DNmQ9/N7NnemvZ+Z6nVCj4x5ENynt1kB0E5MU5kIze9rMbs3nr2XyIY3dOKf3Yes4\ncL+ZLTCzC0IXEwOZd+e/AnTL9QOt3nrXzOYC+zXy1C/cfUYBL6UZ2yzN/G4vBW4E/jt9fhlwNXB+\niUqLC73niu/z7r7ezPYF5ppZfXoEKq3k7p7P/T6tDnV3/2ILfmwtkLlLd6/0Y5Ih39+tmd0CFPIH\nqESy34e9+eDfIKVA7r4+/fU1M5tG1OJSqLfcK2a2n7u/bGbdgVdz/UAp2y+ZfbY64OtmtqeZfQro\nB2imvADp/8C7DSealJbCvHdznZntSTR5Xxe4poplZh81s87p7zsBp6L3ZWvVAd9Jf/8d4N5cP9Cm\nn3xkZsOBa4GuwEwzW+TuX3L3Z83sLuBZYAfwfe3LW7Arzew/iFoILwKjA9dTcZq6uS5wWZWsGzDN\nog2K2gF3uPucsCVVDjObDJwAdDWzNcAvgf8F7jKz84FVwNk5X0dZKiISH1obLiISIwp1EZEYUaiL\niMSIQl1EJEYU6iIiMaJQFxGJEYW6iEiMKNRFRGLk/wEgL6oU3c1GqAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10815ae10>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "xaxis = np.arange(-10, 10, 0.2)\n",
    "yaxis = sp.special.expit(xaxis)\n",
    "plt.plot(xaxis, yaxis)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "__Maximum Likelihood Estimation__\n",
    "\n",
    "The likelihood function $L(w)$ is defined as the probability that the current $w$ assigns to the training set:\n",
    "\n",
    "$$L(w) = \\prod_{i=1}^Np(t_i|x_i;w)$$\n",
    "\n",
    "We have two separate terms for the probability distributions of the two classes, so we combine these two terms into one like:\n",
    "\n",
    "$$p(t_i|x_i;w) = p(t=1|x_i;w)^{t_i}p(t=0|x_i;w)^{1-t_i}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Logistic regression tries to find the $w$ that maximizes the likelihood $L(w)$, which is the same $w$ that maximizes the log-likelihood $l(w) = \\log{L(w)}$.\n",
    "\n",
    "$$\n",
    "\\begin{eqnarray}\n",
    "\\arg{\\max{L(w)}} &=& \\arg{\\max{\\log{L(w)}}} \\nonumber\\\\\n",
    "&=& \\arg{\\max{l(w)}} \\nonumber\\\\\n",
    "\\nabla_wl(w) &=& \\sum_{i=1}^{N}(t^i - \\sigma(w^Tx_i))x^i \\nonumber\n",
    "\\end{eqnarray}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### SGD (Stochastic Gradient Descent)\n",
    "\n",
    "We are currently calling scikit-learn's `.fit` method to train our model. Let us examine how this actually works.\n",
    "\n",
    "<img src=\"https://upload.wikimedia.org/wikipedia/commons/f/ff/Gradient_descent.svg\" height=\"200\" width=\"200\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "__Gradient descent/ascent__\n",
    "* _descent_: find local minimum of function by taking steps proportional to the _negative_ of the gradient of the function at the current point\n",
    "* _ascent_: find local maximum of function by taking steps proportional to the _positive_ of the gradient of the function at the current point"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "__Pseudocode__\n",
    "```\n",
    "n = learning rate\n",
    "x = random\n",
    "repeat\n",
    "    x = x ± n * f'(x)\n",
    "until convergence\n",
    "return x\n",
    "```\n",
    "\n",
    "We repeatedly run through the training set and for each example, we update the parameters according to the gradient of the function with respect to that single example only."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### SVM (Support Vector Machine)\n",
    "\n",
    "SVM is another discriminative model for classification. Although we are able to develop SVMs that do $K$ class classifications (similar to how logistic regression can be extended to multiclass classification ie softmax), we will restrict ourselves to binary classification in this formulation, where the class label is either +1 (positive) or -1 (negative). SVM is not a probabilistic algorithm; it does not optimize a probability measure as a likelihood, as we did in logistic regression. SVM tries to find the ***hyperplane that maximally separates the positive class from the negative class***."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<img src=\"https://upload.wikimedia.org/wikipedia/commons/b/b5/Svm_separating_hyperplanes_%28SVG%29.svg \" height=\"200\" width=\"200\">\n",
    "\n",
    "Let us now formulate the optimization problem for SVM. The decision boundry is the hyperplane $w^T\\phi(x) + b = 0$. Now the linear discriminant function $y(x) = w^T\\phi(x) + b$ with the maximum margin is a good classifier. Margin is the width that the boundry can be increased before hitting a data point. The maximum margin is good because this will be most robust to outliers and thus has strong generalization ability."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Remember the distance from a point $(x_1^*,x_2^*)$ to a line $a^Tx + b = 0$ in $\\mathbb{R}^2$ is:\n",
    "$$\n",
    "\\frac{a_1x_1^* + a_2x_2^* + b}{\\sqrt{a_1^2 + a_2^2}}\n",
    "$$\n",
    "\n",
    "We can extend this for our maximum margin classifier so the distance from a $\\phi(x)$ to the hyperplane $w^T\\phi(x) + b = 0$ is:\n",
    "$$\n",
    "\\frac{t(w^T\\phi(x) + b)}{||w||}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "So now the margin is:\n",
    "$$\n",
    "\\min_n\\frac{t_n(w^T\\phi(x_n) + b)}{||w||}\n",
    "$$\n",
    "where $n$ is the number of training examples.\n",
    "\n",
    "Our optimization problem then becomes:\n",
    "$$\n",
    "\\arg\\max_{w,b}\\frac{1}{||w||}\\min_n[t_n(w^T\\phi(x_n) + b)]\n",
    "$$\n",
    "which is saying that we would like to maximize the distance of the hyperplane from the closest \"support vectors\"."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Now we can rescale $w$ and $b$ since the scaling factor will cancel when dividing by $||w||$. So, the support vectors will be along $w^T\\phi(x) + b = \\pm1$.\n",
    "\n",
    "<img src=\"https://upload.wikimedia.org/wikipedia/commons/2/2a/Svm_max_sep_hyperplane_with_margin.png \" height=\"200\" width=\"200\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Our optimization problem then becomes:\n",
    "$$\n",
    "\\begin{equation*}\n",
    "\\begin{aligned}\n",
    "& \\underset{w, b}{\\text{minimize}}\n",
    "& & \\arg\\min_{w,b}\\frac{1}{2}||w||^2 \\\\\n",
    "& \\text{subject to}\n",
    "& & t_n(w^T\\phi(x_n) + b) \\geq 1, \\; n = 1, \\ldots, N.\n",
    "\\end{aligned}\n",
    "\\end{equation*}\n",
    "$$\n",
    "since maximizing $\\frac{1}{||w||}$ is the same as minimizing $||w||^2$, and the $\\frac{1}{2}$ is for convenience."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Always cool to look up research papers!\n",
    "http://research.microsoft.com/pubs/68541/dagsvm.pdf\n",
    "This one is especially cool since mixing DAGs and SVMs is a great way for companies to assess a candidate's competency with both core algorithms as well as ML algos!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.731707317073\n"
     ]
    }
   ],
   "source": [
    "import sklearn.svm\n",
    "model = sklearn.svm.SVC(kernel=\"rbf\")\n",
    "train_and_test(model, xtrain, ytrain, xtest, ytest)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.80487804878\n"
     ]
    }
   ],
   "source": [
    "import sklearn.neighbors\n",
    "model = sklearn.neighbors.KNeighborsClassifier(3, weights='distance')\n",
    "train_and_test(model, xtrain, ytrain, xtest, ytest)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### kNN (k Nearest Neighbors)\n",
    "\n",
    "We have a dataset $D = \\{(x_1,t_1),\\dots,(x_n,t_n)\\}$. Given an unseen test example $x_{test}$ to predict the example's class $y_{test}$, in kNN we simply:\n",
    "1. Find the k closest examples closest examples to $x_{test}$ from the dataset, i.e. $kNN(x_{test})=\\{(x_1',t_1'),\\dots,(x_k',t_k')\\}$, such that $x_1',\\dots,x_k'$ are the k best points among all training data at minimizing $d(x', x_{test})$, where $d(x',x)$ is a distance measure between $x'$ and $x$.\n",
    "2. Predict the class label $y_{test}$ as the mode class of the k Nearest Neighbors."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Naive Bayes\n",
    "\n",
    "So far we have been working with discriminative models, namely ones that try to model $P(t_i|x_i)$. Now we will introduce a generative model, namely one that models $P(x_i|t_i)$ (class model) and $P(t_i)$ (prior). Let us write $P(t_i|x_i)$ in terms of the class model and prior using Bayes' rule.\n",
    "\n",
    "$$\n",
    "\\begin{eqnarray}\n",
    "P(t,x) &=& P(x,t) \\nonumber\\\\\n",
    "P(t|x)P(x) &=& P(x|t)P(t) \\nonumber\\\\\n",
    "P(t|x) &=& \\frac{P(x|t)P(t)}{P(x)} \\nonumber\\\\\n",
    "&=& \\frac{P(x|t)P(t)}{\\sum_tP(x,t)} \\nonumber\\\\\n",
    "&=& \\frac{P(x|t)P(t)}{\\sum_tP(x|t)P(t)} \\nonumber\\\\\n",
    "&=& \\frac{P(t)\\prod_mP(x_j|t)}{\\sum_tP(t)\\prod_mP(x_j|t)} \\nonumber\\\\\n",
    "\\end{eqnarray}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can model $P(x_j|t=k)$ using a Gaussian distribution, namely\n",
    "$$\n",
    "\\begin{eqnarray}\n",
    "P(x_j|t=k) &=& \\frac{1}{\\sigma_{jk}\\sqrt{2\\pi}}e^{\\frac{-(x-\\mu_{jk})^2}{2\\sigma_{jk}^2}} \\nonumber\\\\\n",
    "\\end{eqnarray}\n",
    "$$\n",
    "where\n",
    "$$\n",
    "\\begin{eqnarray}\n",
    "\\mu_{jk} &=& \\frac{1}{\\sum_nI(t_i=k)}\\sum_nx_{ij}I(t_i=k) \\nonumber\\\\\n",
    "\\end{eqnarray}\n",
    "$$\n",
    "and\n",
    "$$\n",
    "\\begin{eqnarray}\n",
    "\\sigma_{jk}^2 &=& \\frac{1}{\\sum_nI(t_i=k)}\\sum_n(x_{ij}-\\mu_{jk})^2I(t_i=k) \\nonumber\\\\\n",
    "\\end{eqnarray}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.780487804878\n"
     ]
    }
   ],
   "source": [
    "import sklearn.naive_bayes\n",
    "model = sklearn.naive_bayes.GaussianNB()\n",
    "train_and_test(model, xtrain, ytrain, xtest, ytest)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Random Forests"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.926829268293\n"
     ]
    }
   ],
   "source": [
    "import sklearn.ensemble\n",
    "model = sklearn.ensemble.RandomForestClassifier(n_estimators=100)\n",
    "train_and_test(model, xtrain, ytrain, xtest, ytest)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import csv\n",
    "df = pd.read_csv('test.csv', header=0)\n",
    "\n",
    "df['Gender'] = df['Sex'].map( {'female': 0, 'male': 1} ).astype(int)\n",
    "\n",
    "median_ages = np.zeros((2, 3))\n",
    "for i in range(0, 2):\n",
    "    for j in range (0, 3):\n",
    "        median_ages[i, j] = df[(df['Gender'] == i) & (df['Pclass'] == j+1)]['Age'].dropna().median()\n",
    "\n",
    "df['AgeFill'] = df['Age']\n",
    "        \n",
    "for i in range(0, 2):\n",
    "    for j in range(0, 3):\n",
    "        df.loc[(df.Age.isnull()) & (df.Gender == i) & (df.Pclass == j+1), 'AgeFill'] = median_ages[i, j]\n",
    "\n",
    "df['AgeIsNull'] = pd.isnull(df.Age).astype(int)\n",
    "\n",
    "data = np.array(df[['PassengerId', 'Pclass', 'Gender', 'AgeFill', 'Fare']])\n",
    "\n",
    "with open(\"randomforestbasedmodel.csv\", \"wb\") as pf:\n",
    "    predictions = csv.writer(pf, delimiter=',')\n",
    "    predictions.writerow([\"PassengerId\", \"Survived\"])\n",
    "    for example in data:\n",
    "        example = np.nan_to_num(example)\n",
    "        predictions.writerow([str(example[0].astype(int)), str(model.predict(example[1:])[0])])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}