Random Forest Independence Test Algorithm.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Random Forest Independence Test\n",
    "\n",
    "Mutual information is often used as a test for independence for two variables X and Y. Mutual information can be written as follows:  \n",
    "$ I(X, Y) = H(Y) - H(Y|X) $  \n",
    "\n",
    "This can be explained as so: \n",
    "The mutual information of X and Y is the \"amount\" of disorder or entropy in Y that can be explained away by X. Thus, total independence would result in $I(X, Y) = 0$. That is, no information can be explained away by $X$ and thus the conditional entropy of $Y$ knowing $X$ is still just $H(Y)$.  \n",
    "\n",
    "Now, to calculate mutual information, we just have to measure the  the conditional entropy of Y given X (measuring H(Y) is straightforward). We will use random forest to do so.  \n",
    "\n",
    "The algorithm is as follows:\n",
    "![Pseudocode](images/algorithm_v1.png)\n",
    "\n",
    "Note the main point of this algorithm is that random forest is able to capture important features in our data that may classify Y from X. Therefore, after training a random forest, by feeding a new $X_i$ into the forest, we should be able to see how much disorder is left in each partition $X_i$ falls into. This is then our measurement of $H(Y|X_i)$.\n",
    "\n",
    "## Implementation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/lib/python3.7/site-packages/sklearn/ensemble/weight_boosting.py:29: DeprecationWarning: numpy.core.umath_tests is an internal NumPy module and should not be imported. It will be removed in a future NumPy release.\n",
      "  from numpy.core.umath_tests import inner1d\n"
     ]
    }
   ],
   "source": [
    "from sklearn.model_selection import train_test_split\n",
    "import numpy as np\n",
    "from sklearn.ensemble import RandomForestClassifier\n",
    "from sklearn import tree\n",
    "from scipy.stats import entropy\n",
    "\n",
    "def estimate_conditional_entropy_hard_voting(X, y, n_trees = 10, max_depth = None, bootstrap = True):\n",
    "    X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=42)\n",
    "    model = RandomForestClassifier(bootstrap = bootstrap, n_estimators =n_trees, max_depth = max_depth, random_state = 0)\n",
    "    model.fit(X_train, y_train)\n",
    "    class_counts = np.zeros((X_test.shape[0], model.n_classes_))\n",
    "    for tree_in_forest in model:\n",
    "        # get number of training elements in each partition\n",
    "        node_counts = tree_in_forest.tree_.n_node_samples\n",
    "        # get counts w.r.t. testing data now\n",
    "        partition_counts = np.asarray([node_counts[x] for x in tree_in_forest.apply(X_test)])\n",
    "        # get probability\n",
    "        class_probs = tree_in_forest.predict_proba(X_test)\n",
    "        # why are there decimals?!\n",
    "        # bootstrap approximation in sklearn\n",
    "        elems = np.multiply(class_probs, partition_counts[:, np.newaxis])\n",
    "        class_counts += elems\n",
    "    probs = class_counts/class_counts.sum(axis=1, keepdims=True)\n",
    "    entropies = -np.sum(np.log(probs)*probs, axis = 1)\n",
    "    return np.mean(entropies)\n",
    "\n",
    "def estimate_conditional_entropy_soft_voting(X, y, n_trees = 10, max_depth = None, bootstrap = True):\n",
    "    X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=42)\n",
    "    model = RandomForestClassifier(n_estimators = n_trees, max_depth = max_depth, random_state = 0, bootstrap = bootstrap)\n",
    "    model.fit(X_train, y_train)\n",
    "    probs = model.predict_proba(X_test)\n",
    "    entropies = -np.sum(np.log(probs)*probs, axis = 1)\n",
    "    return np.mean(entropies)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are several points we must make about our algorithm. First, we are splitting into training and testing data as we mentioned we would be doing. However, it might be better (and empirically it is better) to run all the data to measure conditional entropy of a dataset. We perform these benchmarks later in this notebook.  \n",
    "\n",
    "Additionally, we have two implementations at the moment. Soft voting and hard voting. Right now the results are very similar however we may want to do further tests later. In future notebooks, we will focus on hard voting just because soft voting does seem to have some approximation errors.\n",
    "\n",
    "## Benchmark\n",
    "We use the following dataset for our benchmarks:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = [0]*20 + [1]*20 + [2]*20 + [3]*20\n",
    "y = [0, 1, 0, 1, 0]*4 + [ 1, 1, 1, 1, 0]*4 + [1, 0, 1, 0, 1]*4 + [0, 0, 0, 0, 1]*4\n",
    "X = np.array(x).reshape(-1, 1)\n",
    "y = np.array(y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "H(X) = 1.386294  \n",
    "H(Y) = 0.693147  \n",
    "H(X, Y) = 1.97300  \n",
    "H(Y|X) = .5867  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.6247088685414731"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "estimate_conditional_entropy_hard_voting(X, y, 10, bootstrap = False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This isn't a very good estimate of our conditional entropy. This is mostly because we can expect our testing data to have different condititional entropy than our training data and both will be different from the conditional entropy of the entire dataset. Indeed, as we increase the size of our dataset, we get slightly better results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5915438951209523"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = [0]*300 + [1]*300 + [2]*300 + [3]*300\n",
    "y = [0, 1, 0, 1, 0]*60 + [ 1, 1, 1, 1, 0]*60 + [1, 0, 1, 0, 1]*60 + [0, 0, 0, 0, 1]*60\n",
    "X = np.array(x).reshape(-1, 1)\n",
    "y = np.array(y)\n",
    "estimate_conditional_entropy_hard_voting(X, y, 10, bootstrap = False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use original dataset\n",
    "x = [0]*20 + [1]*20 + [2]*20 + [3]*20\n",
    "y = [0, 1, 0, 1, 0]*4 + [ 1, 1, 1, 1, 0]*4 + [1, 0, 1, 0, 1]*4 + [0, 0, 0, 0, 1]*4\n",
    "X = np.array(x).reshape(-1, 1)\n",
    "y = np.array(y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we will try our algorithm by running everything on the entire dataset instead of a training testing split. We have to ask ourselves if this counts as \"overfitting\" in the machine learning sense. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "def estimate_conditional_entropy_hard_voting_all(X, y, n_trees = 10, max_depth = None, bootstrap = True):\n",
    "    model = RandomForestClassifier(bootstrap = bootstrap, n_estimators =n_trees, max_depth = max_depth, random_state = 0)\n",
    "    model.fit(X, y)\n",
    "    class_counts = np.zeros((X.shape[0], model.n_classes_))\n",
    "\n",
    "    for tree_in_forest in model:\n",
    "        # get number of training elements in each partition\n",
    "        node_counts = tree_in_forest.tree_.n_node_samples\n",
    "        # get counts for all x (x.length array)\n",
    "        partition_counts = np.asarray([node_counts[x] for x in tree_in_forest.apply(X)])\n",
    "        # get class probability for all x (x.length, n_classes)\n",
    "        class_probs = tree_in_forest.predict_proba(X)\n",
    "        # get elements by performing row wise multiplication\n",
    "        elems = np.multiply(class_probs, partition_counts[:, np.newaxis])\n",
    "        # update counts for that tree\n",
    "        class_counts += elems\n",
    "    probs = class_counts/class_counts.sum(axis=1, keepdims=True)\n",
    "    entropies = -np.sum(np.log(probs)*probs, axis = 1)\n",
    "    # convert nan to 0\n",
    "    entropies = np.nan_to_num(entropies)\n",
    "    return np.mean(entropies)\n",
    "\n",
    "def estimate_conditional_entropy_soft_voting_all(X, y, n_trees = 10, max_depth = None, bootstrap = True):\n",
    "    model = RandomForestClassifier(n_estimators = n_trees, max_depth = max_depth, random_state = 0, bootstrap = bootstrap)\n",
    "    model.fit(X, y)\n",
    "    probs = model.predict_proba(X)\n",
    "    entropies = -np.sum(np.log(probs)*probs, axis = 1)\n",
    "    return np.mean(entropies)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5747883706093118"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "estimate_conditional_entropy_hard_voting_all( X, y, 10, bootstrap = True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is a lot closer to our conditional entropy estimate than the previous iteration. However it is slightly off. After some research, I found out this is because of a bootstrapping approximation that sklearn's random forest does for speed reasons. https://stats.stackexchange.com/questions/130206/sklearn-tree-export-graphviz-values-do-not-add-up-to-samples. \n",
    "\n",
    "Let's turn bootstrapping off then. Note: we may eventually want bootstrapping still. However, in order to do so, we will need to find a different implementation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5867070452737222"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "estimate_conditional_entropy_hard_voting_all( X, y, 10, bootstrap = False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Iris Dataset Test\n",
    "\n",
    "Our final test for this notebook is to measure conditional entropy using this process for the iris datset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,1,'Iris Dataset Scatter Plot')"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWQAAAD7CAYAAABdXO4CAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAHLdJREFUeJzt3X1wVeWdB/DvLze04SrL+3RY0tyr68JqDEoDFcq2cQgsjmx1O9vIbq+IqBNJqg5jdYsbR4Exne3QLdjWUOOMr7kzVFO3O4o7VWEUtcxQYoGUouKsSYiovCaCAZckz/5xTt5zk/Oc3HPuc8/5fmbuxHvuPfc85wR/HM753t8jSikQEVHm5WR6AEREZGFBJiIyBAsyEZEhWJCJiAzBgkxEZAgWZCIiQ7AgB4iIJETk1UyPg7whIkpELsv0OMg7LMhZRESaRGRJqteVUkml1D+4+Nw3ROS8iJwRkc9FpEFE1onIVzU+w5diMdp2ROQrIvKfItIqImftY7YlDdsdcOxFJG6PJXesnz3MZ57tN/Z1Lj7nVhF5O13jIv+wIAdEGgrDXUqpCQBmAPgRgH8B8IqIyJgH568HAMwD8E0AEwBcC+DdTA5oOKP8viYppS4G8K8AHhKR63waFmUYC3KWss+C3hGRzSJyEsD6/mdGYtksIsfss95GEblytM9VSn2hlHoDwA0AFgJYbn/eN0Vkt4i0icgnIvIrEfmK/doue/X99pndChGZLCIvi8hxETlt/3f+oPH/r31W/pGIJPq9dpuIHLLX+72IxFJtZ5hdmA/gv5RSR5WlSSn1bL/P/rqIvGiP66SI/Mpe/jcistNedkJEkiIyyX7tOQAFAF6yt/tvAHrG0mYvWzjS2O3XlIj8UEQOAzjs4HexG8BBAEN+byIyUUSetfejWUQeFJEcEbkcwK8BLLTH1TbadsggSik+suQBoAnAEvu/bwXQCeBuALkAxtvL3rZfXwagAcAkAALgcgAzUnzuGwDuGGb5LgA/tf+7GMACe1txAIcArO33XgXgsn7PpwL4ZwBRWGeqLwD4nf3aRQA+BzDbfj4DQKH93zcC+NAeby6ABwH8IdV2hhnzgwBaAFQCKAIg/V6LANgPYLM9hjwAf2+/dhmApQC+CmC6ve9bhjv29vO4PZbcfsucjP01AFMAjB9m7L2faf/OFgHoAFA6eN8BPAvgv+1jGwfwAYDb+/3ZeDvTf1750H9kfAB8aPyyhhbklkGv9y/Ii+3/SRcAyBnlc1MV5G0AnkixzlpYZ6I9z0crlFcDOG3/90UA2uyCPX7Q+/6np7DYz3PsohRzuJ0IgB8CeAfAlwCOAlhlv7YQwPH+RXSEz/knAH8a7tjbz4cryE7GvniEbfZ8ZhuA07D+0rtn8DG29/H/AFzR77U7Abwx+M8BH9n14CWL7HYk1QtKqZ0AfgXgMQDHRKRWRP5K8/NnAjgFACIyy77s8KmIfA7gJwCmpVpRRKIi8rj9z+nPYZ1xThKRiFLqCwArAKwB8ImIbBeRv7NXjQF41L400mZvX+yxjEop1aWUekwptQjWvw6qATxp/1P+6wCalVKdw4z3ayKyTUQ+tsdbN9L+peBk7Cl/Z/1MU0pNVkpdrpT6xXCvAxgHoLnfsmY4PEZkLhbk7DZiqz6l1C+UUsUArgAwC8D9Tj9YRL4O6zLFW/airQDeA/C3Sqm/AvDvsIpNKj8CMBvANfb7v9Pz0fbYfq+UWgrrcsV7AJ6wXz8C4E6l1KR+j/FKqT84HXsPpdQ5pdRjsM42r7A/uyDFDbWfwDqeRfZ4bx60f4OP9XDH3snY09Fe8QSAC7D+AuhRAODjNG6DMoAFOaBEZL6IXCMi4wB8AeA8gG4H60VFpATW9ck9AF6xX5oA67rvWftstmLQqp8BuLTf8wkAzsG66TUFwMP9tvE1EblRRC6CdVnhbL+x/RrAAyJSaL93ooiUjbCdweNfKyLXish4EckVkVX2WP5k788nAP5DRC4SkTwRWdRvvGcBtIvITAz9y2vwdo/bY+6/bLSxp4VSqgvA8wCqRWSCfePwXlhn9T1jze+56UpZJNPXTPhw/sDQa8hvD3q9dxmAUgAHYBWZEwCSAC5O8blvwCrYZ+zHnwBUAcjr957vwDqTPQvrrHlj/+3DvvwA6/rnTQD+2v7cs7CuZd+JvhtWMwC8CaDdfv8bGHg9dCWARlh/ARwB8GSq7QyzL+Wwbmb2fPYeAP/Y7/UCAL8DcNI+Lr+wlxfa650FsA/WGX5rv/VuhHWzsA3AffayjbAKcxuABQ7GPtr17zgGXZce9Hr/m3qTYRXg4/Z2HoJ9rwDAVwBsh3XJ5ESm/9zy4fwh9i+QiIgyjJcsiIgMwYJMRGQIFmQiIkOwIBMRGYIFmYjIEFodwqZNm6bi8bhHQyEiCqaGhoYTSqnpo71PqyDH43Hs3bvX/aiIiEJIRJpHfxcvWRARGYMFmYjIECzIRESGSNt8YESUXS5cuIDW1lacP38+00MJjLy8POTn52PcuHGu1mdBJgqp1tZWTJgwAfF4HJJ1UyeaRymFkydPorW1FZdccomrz+AlC6KQOn/+PKZOncpinCYigqlTp47pXxwsyJT1ko1JxLfEkbMhB/EtcSQbk5keUtZgMU6vsR5PFmTKasnGJMpfKkdzezMUFJrbm1H+UjmLcsA8/fTTOHr0aKaH4TkWZMpqVTuq0HGhY8CyjgsdqNpRlaERkRdYkImyQEt7i9Zyci+ZBOJxICfH+pkc4z9CvvjiCyxfvhxXXXUVrrzySvzmN79BQ0MDSkpKUFxcjGXLluGTTz5BfX099u7di0Qigauvvhrnzp3Djh07MHfuXBQVFeG2227Dl19+CQBYt24drrjiCsyZMwf33XcfAOCll17CNddcg7lz52LJkiX47LPPxjZwL+lML1JcXKyITBLbHFNYjyGP2OZYpodmvL/85S+O31tXp1Q0qhTQ94hGreVu1dfXqzvuuKP3eVtbm1q4cKE6duyYUkqpbdu2qdWrVyullCopKVF//OMflVJKnTt3TuXn56v3339fKaXUypUr1ebNm9WJEyfUrFmzVHd3t1JKqdOnTyullDp16lTvsieeeELde++97gftwHDHFcBe5aDG8gyZslp1aTWi46IDlkXHRVFdWp2hEQVTVRXQMfDKEDo6rOVuFRUV4bXXXsOPf/xjvPXWWzhy5Aj+/Oc/Y+nSpbj66qvxyCOPoLW1dch677//Pi655BLMmjULALBq1Srs2rULEydORF5eHm6//Xa8+OKLiEatPxetra1YtmwZioqKsGnTJhw8eND9oD3GgkxZLVGUQO13axGbGINAEJsYQ+13a5EoSmR6aIHSkuIKUKrlTsyaNQvvvvsuioqK8OCDD+K3v/0tCgsLsW/fPuzbtw+NjY149dVXHX9ebm4u9uzZg+9///t4+eWXcd111wEA7r77btx1111obGzE448/bvQXYfjFEMp6iaIEC7DHCgqA5mH6lRUUuP/Mo0ePYsqUKbj55psxadIk1NTU4Pjx49i9ezcWLlyICxcu4IMPPkBhYSEmTJiAM2fOAABmz56NpqYmfPjhh7jsssvw3HPPoaSkBGfPnkVHRweuv/56LFq0CJdeeikAoL29HTNnzgQAPPPMM+4H7AMWZDJKsjGJqh1VaGlvQcHEAlSXVrPYGqC6GigvH3jZIhq1lrvV2NiI+++/Hzk5ORg3bhy2bt2K3Nxc3HPPPWhvb0dnZyfWrl2LwsJC3HrrrVizZg3Gjx+P3bt346mnnkJZWRk6Ozsxf/58rFmzBqdOncKNN96I8+fPQymFn//85wCA9evXo6ysDJMnT8bixYvx0UcfjfFoeEes683OzJs3T7EfMnmlJ1PcP8YWHRflJQiPHDp0CJdffrnj9yeT1jXjlhbrzLi6Gkjw1zLEcMdVRBqUUvNGW5fXkMkYzBSbLZEAmpqA7m7rJ4tx+rEgkzGYKaawY0EmYxRMHP4OUarlREHDgkzGYKaYwo4FmYzBTDGFHWNvIWZixIyZYgozniGHFNtWUra7/vrr0dbWpr3e+vXr8bOf/cyDEY0dC3JIMWJG2Uophe7ubrzyyiuYNGlSxsbR2dmZ9s9kQQ4pRsxIV7pnZlm3bh0ee+yx3uc9Z66bNm3C/PnzMWfOHDz88MMAgKamJsyePRu33HILrrzyShw5cgTxeBwnTpwAADz77LOYM2cOrrrqKqxcubJ3ncWLF2POnDkoLS1FyzCNN/bt24cFCxZgzpw5+N73vofTp08DAK699lr0fAnuxIkTiMfjAKy+zDfccAMWL16M0tLSMe3/cFiQQ4oRM9LhxSWuFStW4Pnnn+99/vzzz2P69Ok4fPgw9uzZg3379qGhoQG7du0CABw+fBiVlZU4ePAgYrFY73oHDx7EI488gp07d2L//v149NFHAVhNhVatWoUDBw4gkUjgnnvuGTKGW265BT/96U9x4MABFBUVYcOGDaOO+91330V9fT3efPNN1/ueCgtySDFiRjq8uMQ1d+5cHDt2DEePHsX+/fsxefLk3g5vc+fOxTe+8Q289957OHz4MAAgFothwYIFQz5n586dKCsrw7Rp0wAAU6ZMAQDs3r0bP/jBDwAAK1euxNtvvz1gvfb2drS1taGkpARAXxvP0SxdurR3G+nGlEVI9SQZTEtZkJm8usRVVlaG+vp6fPrpp1ixYgWam5vxwAMP4M477xzwvqamJlx00UVj2paO3NxcdHd3A8CQdp1ejoNnyCGWKEqgaW0Tuh/uRtPaJk+Kse51x8rtlcjdmAvZIMjdmIvK7ZVpHxPp8+oS14oVK7Bt2zbU19ejrKwMy5Ytw5NPPomzZ88CAD7++GMcO3ZsxM9YvHgxXnjhBZw8eRIAcOrUKQDAt771LWzbtg0AkEwm8e1vf3vAehMnTsTkyZPx1ltvAUBvG08AiMfjaGhoAADU19ePaR918AyZPDO4e1vPdUcAwxb/yu2V2Lp3a+/zLtXV+7xmeY0PI6ZUqkurh+3EN9ZLXIWFhThz5gxmzpyJGTNmYMaMGTh06BAWLlwIALj44otRV1eHSCQy4mdUVVWhpKQEkUgEc+fOxdNPP41f/vKXWL16NTZt2oTp06fjqaeeGrLuM888gzVr1qCjowOXXnpp73vuu+8+3HTTTaitrcXy5cvHtI862H6TPBPfEkdz+9Cu5rGJMTStbRqyPHdjLrpU15DlEYmg86H0R4zCTrv9poFfJDLRWNpv8gyZPKN73XG4YjzScvIXv0XpPV5DJs/oXneMyPD/LE21nChoWJDJM7rRuvLicq3lREHDgkye0e3eVrO8BhXzKnrPiCMSQcW8Ct7Q85DOPSQa3ViPJ68hk2N+3NSpWV7DAuyTvLw8nDx5ElOnToWIZHo4WU8phZMnTyIvL8/1Z7AgkyO6ETa365B/8vPz0draiuPHj2d6KIGRl5eH/Px81+sz9kaO6EbY3K5DFEScdZrSys1XZ9lRjkgPCzI54uars+woR6SHBZkccdMdjh3liPSwIJMjbiYg5aSlRHp4U4+IyGO8qUejSveUPEQ0NswhhxQzwkTm4RlySHHWaSLzsCCHFDPCROZhQQ4pZoSJzMOCHFLMCBOZhwU5pJgRJjIPC3KI6c46bWpMztRxEeli7I0cMTUmZ+q4iNzgGTI5YmpMztRxEbnBgkyOmBqTM3VcRG6wIJMjpsbkTB0XkRssyOSIqTE5U8dF5AYLMjliakzO1HERuRH69pt+zKTsxzbIOf4+yG9O22+GOvbmR2SKsSyz8PdBJgv1GbIfsyJz5mWz8PdBmcAG9Q74EZliLMss/H2QyUJdkP2ITDGWZRb+PshkoS7IfkSmGMsyC38fZLJQF2Q/IlOMZZmFvw8yWeBu6jHS5B0eWyJ3Qhl7Y6TJOzy2RN4L1CULdv7yDo8tkfcCVZAZafIOjy2R9wJVkBlp8g6PLZH3AlWQGWnyDo8tkfcCVZAZafIOjy2R9wKVsgCswsEi4Uzl9krUNtSiS3UhIhGUF5ejZnlNWrfBqByRc4EryORM5fZKbN27tfd5l+rqfT5cUXYTe2NUjkhP4L4YQs7kbsxFl+oasjwiEXQ+1DlkuZsuaeysRmRhtzca0XDFeKTlbmJvjMoR6WFBDqmIRLSWu4m9MSpHpIcFOaTKi8u1lruJvTEqR6SHBTmkapbXoGJeRe8ZcUQiqJhXkTJl4Sb2xqgckR7e1CMi8lhob+olG5OIb4kjZ0MO4lviSDYmMz0kX/ix32E9tkC49538E6gcclhzr5w921th3nfyV6AuWYQ198rZs70V5n2n9AjlJYuw5l45e7a3wrzv5K9AFeSw5l45e7a3wrzv5K9AFeSw5l45e7a3wrzv5K9AFeSw5l45e7a3wrzv5K9A3dQj8/jRfpMtPsl0oZx1mszCOB6RnkBdsiCz+DFTNWfDpiBhQSbPMI5HpIcFmTzDOB6RHhZk8gzjeER6WJDJM4zjEelh7M1QulEuP2aQJnIimQSqqoCWFqCgAKiuBhIh//uRsbcsphvl0p1BmsgrySRQXg502MGX5mbrOcCi7ATPkA2k211MdwZpIq/E41YRHiwWA5qa/B6NOULZ7S0odKNcujNIE3mlJUXaMNVyGogF2UC6US7dGaSJvFKQIm2YajkNxIJsIN0ol+4M0kReqa4GogP/6CIatZbT6FiQDaQb5dKdQZrIK4kEUFtrXTMWsX7W1vKGnlMsyAFRs7wGnQ91Qj2s0PlQp6NiXLm9ErkbcyEbBLkbc1G5vXLE97uZ6JOTg5olmbRuvOXkWD+THvw6EgnrBl53t/UzW4uxH8dqMMbeDORHBzPdqJybMbETm1kYSXMuU8eKsTcD+TGppm5Uzs2YODmoWRhJcy7dx4qxtyzmRwcz3aicmzGxE5tZGElzLlPHigXZQH50MNONyrkZEzuxmYWRNOcydaxYkA3kRwcz3aicmzGxE5tZGElzLmPHSinl+FFcXKzIH3UH6lRsc0zJelGxzTFVd6Au7duoeLlCRTZEFNZDRTZEVMXLFWkfkx/7Qc7V1SkViyklYv2s468jpXQeKwB7lYMay5t6REQe4009g/iRxTV1G7pZZ6IwYw7ZY0GZednNNtgWlEgPL1l4zI8srqnbYFtQIgsvWRgiKDMvu9kG24IS6WFB9lhQZl52sw22BSXSw4LssaDMvOxmG2wLSqSHBdljQZl52c022BaUSJOTsLIK8BdD3HxxgV+ooNHwCxjeyrbjC4dfDAl17M2PKBfbVoYP21x6K8jHN9SxNz+iXGxbGT5sc+mtbDy+jL054EeUi20rw4dtLr0V5OMb6oLsR5SLbSvDh20uvRXk4xvqguxHlIttK8OHbS69Fejj6+TOn2LKYgCmLGg02ZYCyDbZdnzhMGUR+oJM3grKXywVFUpFItb/MZGI9TzT/CpKJu57tnFakEMdeyNvBSW+V1kJbO1LOqKrq+95TYa+4+JX9MvEfQ+yUMfeyFtBie/l5lqFaLBIBOjMUNM6v6JfJu57NmLsjTIuKPG94QrSSMv94Ff0y8R9DzIWZPJMUOJ7kRTN6VIt94Nf0S8T9z3IWJDJM0GJ75WnaE6Xarkf/Ip+mbjvQcaCTJ7xowudH2pqgIqKvrPCSMR6nsmbWokEUFtrXTMWsX7W1qa/l4OJ+x5kTFmQpxJFiawrwMOpqQlvEQrzvvuNBZkoCwW541mY8ZIFURaqquorxj06OqzllL1YkImyUJA7noUZCzJRFgpyx7MwY0EmykKB7ngWYizIRFnIr9gb+YsFOSCSjUnEt8SRsyEH8S1xJBuTmR5SoCSTVv+InBzrZ3KUw6v7fjcSCatvRXe39TObi7EfxysbMPYWAEHpqmYq3YgZI2l6eLz6sNtbAASlq5qpdDurZeMknJkUhuPFbm8hEpSuaqbSjZgxkqaHx6sPC3IABKWrmql0I2aMpOnh8erDghwAQemqZirdiBkjaXp4vPqwIAdAULqqmUo3YsZImh4erz68qUdE5DHe1DMIM8JmcZN51V2nstKaj07E+llZOfZxj3VMALBkiTWmnseSJekfF42Bk6mpex7FxcV+zZodGHUH6lS0OqqwHr2PaHVU1R3waM52GlFdnVLRqDWlfc8jGrWWp2udioqB7+15VFRkdj9KS4cfV2lp+sZFwwOwVzmosbxk4TFmhM3iJvOqu44fMzW72Q+R1J+nUQbIBV6yMAQzwmZxk3nVXcePmZqZ3Q0mFmSPMSNsFjeZV911/JipmdndYGJB9hgzwmZxk3nVXcePmZrd7Edpqd5yygAnF5p7Hryp507dgToV2xxTsl5UbHOMN/QyrK5OqVhMKRHr50g3wtyuU1GhVCRi3TSLRNJ7Q8/tmJQaemOPN/T8AYc39Xw9QzYx/uXHmBJFCTStbUL3w91oWtsUqi9s+NFW0Y9tvPMO0NpqlbHWVut5uunuh5v2m6+/PjBj8frrYx83pZGTqq3ScIZsYvzLxDEFiZtoltfbcDMm3Ribm9ibH8eKMgemxd5MjH+ZOKYg8aOtoh+tMXVjbG5ib2FoQRlmxsXeTIx/mTimIPEjmuVHa0zdGJub2BtjbAT4WJBNjH+ZOKYg8SOa5UdrTN0Ym5vYG2NsBPhYkE2Mf5k4piDxo62iH60xdWNsbmJvbEFJAPyNvZkY/zJxTEHiJprl9TbcjEk3xuYm9ubHsaLMgGk39YIk2ZhE1Y4qtLS3oGBiAapLq0MVZdORTAJVVda10IIC64xvpHjWkiXAjh19z0tLGc2i7GfcTb2g6Jnhubm9GQqqd4ZnEzLVpumZTbi52Qpy9cwmnCpfO7gYA9ZztoiksOAZsiZG5ZzTjXKxGxkFFc+QPcKonHOMchHpYUHWxKicc4xyEelhQdbEqJxzulEudiOjsGNB1sQZnp3TnU349deHFl+mLChMWJBdCHP3Nl26XdJWrx5YwFev9meco/Fj0lKi3EwPgIKrshLYurXveVdX3/OamqHv74nJdXRYz3ticoCz1pJe0d0PIrcYeyPP6HY9M7XjmR+TllKwMfZGGafb9czUmJwfk5YSASzI5CHdrmemxuT8mLSUCGBBJg/pdj0zteOZH5OWEgEsyOShmhqgoqLvTDISsZ6nuhGmG5Pzi+5+ELnFm3pERB7jTT0KDWaEKSiYQ6asxowwBQnPkCmr1dbqLScyGQsyZTVmhClIWJApqzEjTEHCgkxZjRlhChLe1KOs1nPjrrbWukwRiVjFmDf0KBvxDJmMkkxaTYZycqyfqSZE7W/RIiA/34q95edbz00Yl4nbIMMppRw/iouLFZFX6uqUikaVsronW49o1FqeznX8GJeJ26DMAbBXOaix/KYeGcNN+00/WnYGZRuUOU6/qceCTMbIybHODQcTAbq707eOH+MycRuUOfzqNGUdN+03/WjZGZRtkPlYkMkYbtpv+tGyMyjboCzg5EJzz4M39chrdXVKxWJKiVg/ndzUcrOOH+MycRuUGeBNPWeSSaCqypomqKDAOiPJdP9dIgoWp9eQQ/3FEFNnOSaicAr1NeSqqr5i3KOjw1pOROS3UBdkU2c5JqJwCnVBZtSIiEwS6oLMqBERmSTUBdnUWY6JKJxCXZABq/g2NVlfT21qYjEeiZtuZOySRuRcqGNv5JybiKAfsUJGFylIQv/FEHImzJ3YiMaKzYUordxEBP2IFTK6SEHCgkyOhLkTG5FfWJDJkTB3YiPyCwsyOeImIuhHrJDRRQoSpizIsUTCzEJn6riIdLEgk2cYSSPSw0sW5Bl20yPSw4JMnmEkjUgPCzJ5hpE0Ij0syOQZRtKI9LAgk2cYSSPSw5QFeYqRNCLneIZMRGQIFmQiIkOwIBMRGYIFmYjIECzIRESGYEEmIjKE1hROInIcwDAT5hAR0QhiSqnpo71JqyATEZF3eMmCiMgQLMhERIZgQSYiMgQLMhGRIViQiYgMwYJMRGQIFmQiIkOwIBMRGYIFmYjIEP8PI0vx+bxvIlkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from sklearn import datasets\n",
    "\n",
    "iris = datasets.load_iris()\n",
    "# take first 2 features\n",
    "X = iris.data[:100, :2]\n",
    "# take first 100 (only two classes)\n",
    "y = iris.target[:100]\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "plt.scatter(X[:50, 1], X[:50, 0], c = \"b\", label = \"setosa\")\n",
    "plt.scatter(X[50:, 1], X[50:, 0], c = \"g\", label = \"vericolour\")\n",
    "plt.legend()\n",
    "plt.xticks([])\n",
    "plt.yticks([])\n",
    "plt.title('Iris Dataset Scatter Plot')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "H(Y): 0.6931471805599453\n",
      "Conditional Entropy: 0.0\n",
      "Mutual Information: 0.6931471805599453\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/Cellar/jupyter/1.0.0_5/libexec/lib/python3.7/site-packages/ipykernel_launcher.py:17: RuntimeWarning: divide by zero encountered in log\n",
      "/usr/local/Cellar/jupyter/1.0.0_5/libexec/lib/python3.7/site-packages/ipykernel_launcher.py:17: RuntimeWarning: invalid value encountered in multiply\n"
     ]
    }
   ],
   "source": [
    "import scipy\n",
    "\n",
    "def calculate_entropy(X):\n",
    "    probs = np.bincount(X)/len(X)\n",
    "    return scipy.stats.entropy(probs)\n",
    "entropy = calculate_entropy(y)\n",
    "cond_entropy = estimate_conditional_entropy_hard_voting_all(X, y, 10, bootstrap = False)\n",
    "print(\"H(Y):\", entropy)\n",
    "print(\"Conditional Entropy:\", cond_entropy)\n",
    "print(\"Mutual Information:\", entropy - cond_entropy)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is good! This is exactly what we wanted to see since the two classes in iris are clearly separable (we are not plotting another dimension) given the features. So, the conditional entropy after knowing X should b 0. Meaning, the features of the flowers and the class labels of the flowers are dependent!"
   ]
  },
  {
   "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.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}