Naive_bayes.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Naive Bayes "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 190,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-02-15T13:34:53.921859Z",
     "start_time": "2018-02-15T13:34:53.919761Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "np.set_printoptions(precision=6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-02-15T04:35:15.580155Z",
     "start_time": "2018-02-15T04:35:15.573258Z"
    }
   },
   "source": [
    "## Bayes theorem\n",
    "$$\\underbrace{P(B|A)}_2 = \\underbrace{\\frac{P(A|B)}{P(A)}}_3 \\underbrace{P(B)}_1$$\n",
    "<ol>\n",
    "<li>The prior is what you believed about $B$ before having encountered a new and relevant piece of information (i.e., $A$). In our case, this is the probability that a given document belongs to class $B$(China) and class $\\overline{B}$(not China).</li>\n",
    "\n",
    "<li>The posterior is what you believe (or ought to, if you are rational) about $B$ after having encountered a new and relevant piece of information. With more and more training examples, this should get close to the ground truth. In our case, this means that given information $A$(training data), what is the probability that the it came from class $B$ </li>\n",
    "\n",
    "<li>The quotient of the likelihood divided by the marginal probability of the new piece of information indexes the informativeness of the new information for your beliefs about $B$.  </li> \n",
    "\n",
    "Bayes theorem can be simplified into\n",
    "$${P(B|A)} \\propto {{P(A|B)}}{P(B)}$$\n",
    "This is because the marginal probability $P(A)$ is a normalizing constant, i.e., it's present in calculating posterior of all classes and all it does is give posterior as probabilities (sum to 1). We don’t care what the true posterior value is (might be needed when comparing different models), we just want to know which class has the higher posterior value. Furthermore, computing the marginal probability isn't trivial in most datasets. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-02-15T07:57:33.792473Z",
     "start_time": "2018-02-15T07:57:33.786890Z"
    }
   },
   "source": [
    "## Aim"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The training set below can be interprted as words in a document and the labels represent the class of the document which can be China or not China. <br>\n",
    "<b>Training Set</b> <br>\n",
    "<table class=\"tg\">\n",
    "  <tr>\n",
    "    <th class=\"tg-yw4l\">DocID</th>\n",
    "    <th class=\"tg-yw4l\">Words in Doc<br></th>\n",
    "    <th class=\"tg-yw4l\">Doc class = China?<br></th>\n",
    "  </tr>\n",
    "  <tr>\n",
    "    <td class=\"tg-yw4l\">1</td>\n",
    "    <td class=\"tg-yw4l\">Chinese Beijing Chinese<br></td>\n",
    "    <td class=\"tg-yw4l\">1<br></td>\n",
    "  </tr>\n",
    "  <tr>\n",
    "    <td class=\"tg-yw4l\">2</td>\n",
    "    <td class=\"tg-yw4l\">Chinese Chinese Shanghai<br></td>\n",
    "    <td class=\"tg-yw4l\">1<br></td>\n",
    "  </tr>\n",
    "  <tr>\n",
    "    <td class=\"tg-yw4l\">3</td>\n",
    "    <td class=\"tg-yw4l\">Chinese Macao<br></td>\n",
    "    <td class=\"tg-yw4l\">1<br></td>\n",
    "  </tr>\n",
    "  <tr>\n",
    "    <td class=\"tg-yw4l\">4</td>\n",
    "    <td class=\"tg-yw4l\">Tokyo Japan Chinese<br></td>\n",
    "    <td class=\"tg-yw4l\">0<br></td>\n",
    "  </tr>\n",
    "</table>\n",
    "\n",
    "An array represent a document and its index a word and its value represent the count of appearance. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-02-15T05:54:34.676095Z",
     "start_time": "2018-02-15T05:54:34.672975Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "X = np.array([\n",
    "    [2,1,0,0,0,0],\n",
    "    [2,0,1,0,0,0],\n",
    "    [1,0,0,1,0,0],\n",
    "    [1,0,0,0,1,1]\n",
    "])\n",
    "y = np.array([1,1,1,0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our goal is to find the best class for the document. The best class in NB classification is the most likely or maximum a posteriori ( MAP ) class which is simply MLE(Maximum Likelihood Estimate) multiplied by a prior, where MLE is the relative frequency and corresponds to the most likely value of each parameter given the training data. For instance the MLE of the prior for class $C$ in this example is ${\\frac{C}{C+\\overline{C}}}$.\n",
    "<br><br>\n",
    "$${c}_{MAP} = \\text{argmax}_{c \\in C} \\; \\; [ \\log \\hat P(c) + \\sum_i  \\log \\hat{P}(t_i|c)]$$\n",
    "Where $t_i$s are the unique terms in the document. They are essentially the features of our model. \n",
    "One immediate question is why logs and why are they summed? The answer is that in Eq(2), there maybe many features and many conditional probabilities are multiplied. This can result in a floating point underflow. It is therefore better to perform the computation by adding logarithms of probabilities instead of multiplying probabilities. The class with the highest log probability score is still the most probable; $\\log (x+y) = \\log (x) + \\log (y))$ and the logarithm function is monotonic."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-02-15T09:08:28.856924Z",
     "start_time": "2018-02-15T09:08:28.839425Z"
    }
   },
   "outputs": [],
   "source": [
    "class MultinomailNB(object):\n",
    "    def __init__(self, alpha = 1):\n",
    "        self.alpha = alpha #See [1]\n",
    "    def fit(self, X, y):\n",
    "        #Grouping training data by class\n",
    "        separated = [[x for x,t in zip(X,y) if t == c] for c in np.unique(y)]\n",
    "        print ('Separated Data: ',separated)\n",
    "        count_sample = X.shape[0]\n",
    "        self.class_log_prior = [np.log(len(x)/count_sample) for x in separated]\n",
    "        print ('Priors: ', self.class_log_prior)\n",
    "        count = np.array([np.sum(x, axis=0)+self.alpha for x in separated])\n",
    "        print ('Counts of appearances: ',count)\n",
    "        #log probabilities of each word\n",
    "        self.feature_log_prob = np.log(count/count.sum(axis=1)[np.newaxis].T)\n",
    "        print ('Probabilities of each term: ', self.feature_log_prob)\n",
    "    def predict_log_proba(self, X):\n",
    "        return [(self.feature_log_prob*x).sum(axis=1)+self.class_log_prior for x in X]\n",
    "    def predict(self, X):\n",
    "        return [np.argmax(x) for x in self.predict_log_proba(X)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-02-15T07:45:57.780528Z",
     "start_time": "2018-02-15T07:45:57.764773Z"
    }
   },
   "source": [
    "[1] We use the alpha as a smoothing factor; to fix a problem with the MLE estimate whereby probability is zero for a term-class combination that did not occur in the training data.  If the term Mandarin in the training data only occurred in China documents, then the MLE estimates for not China, will be zero if Mandarin appeared in a test sample with a class not China. The zero probability for Mandarin cannot be \"conditioned away\" no matter how strong the evidence for the class not China from other features. The estimate is 0 because of sparseness: the training data are never large enough to represent the frequency of rare events adequately."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Train"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-02-15T08:59:20.239173Z",
     "start_time": "2018-02-15T08:59:20.227025Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Separated Data:  [[array([1, 0, 0, 0, 1, 1])], [array([2, 1, 0, 0, 0, 0]), array([2, 0, 1, 0, 0, 0]), array([1, 0, 0, 1, 0, 0])]]\n",
      "Priors:  [-1.3862943611198906, -0.2876820724517809]\n",
      "Counts of appearances:  [[2 1 1 1 2 2]\n",
      " [6 2 2 2 1 1]]\n",
      "Probabilities of each term:  [[-1.504077 -2.197225 -2.197225 -2.197225 -1.504077 -1.504077]\n",
      " [-0.847298 -1.94591  -1.94591  -1.94591  -2.639057 -2.639057]]\n"
     ]
    }
   ],
   "source": [
    "model = MultinomailNB()\n",
    "model.fit(X,y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Test"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Data\n",
    "<table class=\"tg\">\n",
    "  <tr>\n",
    "    <th class=\"tg-yw4l\">docID</th>\n",
    "    <th class=\"tg-yw4l\">Words in Document<br></th>\n",
    "  </tr>\n",
    "  <tr>\n",
    "    <td class=\"tg-yw4l\">1</td>\n",
    "    <td class=\"tg-yw4l\">Chinese Chinese Chinese Tokyo Japan</td>\n",
    "  </tr>\n",
    "  <tr>\n",
    "    <td class=\"tg-yw4l\">2</td>\n",
    "    <td class=\"tg-yw4l\">Beijing Shanghai Tokyo Japan</td>\n",
    "  </tr>\n",
    "</table>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-02-15T09:05:02.284692Z",
     "start_time": "2018-02-15T09:05:02.281594Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "X_test = np.array([[3,0,0,0,1,1],[0,1,1,0,1,1]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-02-15T09:09:41.018625Z",
     "start_time": "2018-02-15T09:09:41.014004Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1, 0])"
      ]
     },
     "execution_count": 75,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.predict(X_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Guassian Naive Bayes\n",
    "When the data has continuous variables, then the above model is not suitable. Fitting the model to a gaussian distribution is a good solution for working with such data. Note that not every data follows a Gaussian distribution, but it's a pretty good assumption."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 383,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-02-15T14:21:43.270912Z",
     "start_time": "2018-02-15T14:21:43.268672Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.cross_validation import train_test_split\n",
    "from sklearn import metrics\n",
    "from sklearn.metrics import precision_score"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data\n",
    "HTRU2 pulsar candidate [dataset](http://archive.ics.uci.edu/ml/datasets/HTRU2#)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 500,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-02-15T14:51:45.020003Z",
     "start_time": "2018-02-15T14:51:45.017671Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "path = 'data/HTRU2/HTRU_2.csv'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 501,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-02-15T14:51:45.209581Z",
     "start_time": "2018-02-15T14:51:45.181533Z"
    }
   },
   "outputs": [],
   "source": [
    "full = pd.read_csv(path, header=None)\n",
    "X = full.drop(8, 1)\n",
    "X = X.as_matrix()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 502,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-02-15T14:51:45.373331Z",
     "start_time": "2018-02-15T14:51:45.370284Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "9.157447759526203% of the data is true labels\n"
     ]
    }
   ],
   "source": [
    "y = full.as_matrix(columns=[8])\n",
    "class_names=[0,1]\n",
    "#Number of true labels\n",
    "print ('{}% of the data is true labels'.format(np.count_nonzero(y)/len(full)*100))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 503,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-02-15T14:51:45.560934Z",
     "start_time": "2018-02-15T14:51:45.557862Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#Normalize data\n",
    "scaler = StandardScaler()\n",
    "train = scaler.fit_transform(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 504,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-02-15T14:51:45.749545Z",
     "start_time": "2018-02-15T14:51:45.744008Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=.25, stratify=y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 528,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-02-15T15:02:36.253472Z",
     "start_time": "2018-02-15T15:02:36.227754Z"
    }
   },
   "outputs": [],
   "source": [
    "class GaussianNB(object):\n",
    "    def __init__(self):\n",
    "        pass\n",
    "\n",
    "    def fit(self, X, y):\n",
    "        self.separated = [[x for x, t in zip(X, y) if t == c] for c in np.unique(y)]\n",
    "        self.model = np.array([np.c_[np.mean(x, axis=0), np.std(x, axis=0)] for \n",
    "                              x in self.separated])\n",
    "        count_sample = X.shape[0]\n",
    "        self.class_log_prior = [np.log(len(x)/count_sample) for x in self.separated]\n",
    "        return self\n",
    "    def gaussian(self, x, mean, std):\n",
    "        return np.log((np.exp(-1*(x-mean)**2/(2*std**2)))/(std*np.sqrt(2*np.pi)))\n",
    "    \n",
    "    def predict_log_proba(self, X):\n",
    "        results = []\n",
    "        for x in X:\n",
    "            example = []\n",
    "            for c in self.model:\n",
    "                prob = 0\n",
    "                for e, i in zip(x, c):\n",
    "                    prob += self.gaussian(e, *i)\n",
    "                example.append(prob)\n",
    "            results.append(example)\n",
    "        return np.asarray(results)\n",
    "    def predict(self, X):\n",
    "        return [np.argmax(x) for x in (self.predict_log_proba(X)+self.class_log_prior)] #Watch out for the prior"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 506,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-02-15T14:51:46.490126Z",
     "start_time": "2018-02-15T14:51:46.484834Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def accmetrics(labels, preds, averaging):\n",
    "    prec = precision_score(labels, preds, average=averaging)*100\n",
    "    rand = metrics.adjusted_rand_score(labels.squeeze(), preds.squeeze())*100\n",
    "    return rand, prec "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 507,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-02-15T14:51:47.504859Z",
     "start_time": "2018-02-15T14:51:47.093874Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/anaconda/lib/python3.6/site-packages/ipykernel_launcher.py:13: RuntimeWarning: divide by zero encountered in log\n",
      "  del sys.path[0]\n"
     ]
    }
   ],
   "source": [
    "nb = GaussianNB().fit(X_train, y_train)\n",
    "preds = nb.predict(X_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 508,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-02-15T14:51:49.282725Z",
     "start_time": "2018-02-15T14:51:49.275570Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(64.174310857987322, 62.110726643598611)"
      ]
     },
     "execution_count": 508,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "accmetrics(y_test, np.asarray(preds), 'binary')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The results are not great but the task wasn't easy either. We can confirm that it works on a much simpler classification task."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 521,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-02-15T14:52:49.052528Z",
     "start_time": "2018-02-15T14:52:49.045986Z"
    }
   },
   "outputs": [],
   "source": [
    "iris = datasets.load_iris()\n",
    "X, y = iris.data, iris.target\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=.25)\n",
    "nb = GaussianNB().fit(X_train, y_train)\n",
    "preds = nb.predict(X_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 526,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-02-15T14:52:55.005823Z",
     "start_time": "2018-02-15T14:52:55.001554Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(85.881730049476928, 94.73684210526315)"
      ]
     },
     "execution_count": 526,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "accmetrics(y_test, np.asarray(preds), 'micro')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another thing to note is the importance of meaninful priors. In nebula detection, adding a naive constant prior improved the accuracy by ~4 points whereas it adversely affected the iris dataset. Perhaps a more meaningful choice given the limited dataset is a uniform distribution. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### References\n",
    "<ol><li>[Kenzo Takahashi's blog post](http://kenzotakahashi.github.io/naive-bayes-from-scratch-in-python.html)</li>\n",
    "<li>[Introduction to Information Retrieval](https://nlp.stanford.edu/IR-book/html/htmledition/naive-bayes-text-classification-1.html)\n"
   ]
  }
 ],
 "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.4"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}