Tutorial_072_InClass_Exercises.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# In class exercise...\n",
    "* MI is biased in that small sample sizes lead to inaccurate estimates of PDFs, and that can sometimes lead to negative MI values (which should never happen in theory). \n",
    "* A common, and simple, approach, is to compute MI with shuffled condition labels (like randomization tests that we did many weeks back) and then subtract the shuffled MI from the actual MI. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.neighbors import KernelDensity\n",
    "# also define the default font we'll use for figures. \n",
    "fig_font = {'fontname':'Arial', 'size':'20'}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## First set up two arrays of data...make them correlated to some degree so that there is a reasonably high MI..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1.        , 0.59529809],\n",
       "       [0.59529809, 1.        ]])"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n = 1000\n",
    "x = np.round(np.random.rand(n))\n",
    "\n",
    "y = np.zeros(n)\n",
    "\n",
    "for i in np.arange(n):\n",
    "    if np.random.rand(1) < .5:\n",
    "        y[i] = x[i]\n",
    "            \n",
    "np.corrcoef(x,y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Then compute the MI between the arrays. Can do two discrete arrays for simplicity, and import the entropy and conditional entropy functions from the tutorial."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "def entropy(x):\n",
    "    \"\"\"compute entropy of discrete array x\n",
    "\n",
    "    Args:\n",
    "        x (int): array of discrete values\n",
    "\n",
    "    Returns:\n",
    "        Hx (float): entropy of x\n",
    "\n",
    "    \"\"\"\n",
    "    # figure out unique values of x - can be more than just 0s, 1s\n",
    "    uniquex = np.unique(x)\n",
    "\n",
    "    Hx = 0\n",
    "    for i in np.arange(len(uniquex)):\n",
    "        # probability that x==uniquex[i]\n",
    "        px = np.sum(x==uniquex[i])/len(x)    \n",
    "\n",
    "        # check for px==0 because log2(0) = -inf\n",
    "        if px!=0:\n",
    "            Hx += (-np.sum( px * np.log2(px) ))  \n",
    "        else:\n",
    "            print('px is zero for value ', i)\n",
    "        \n",
    "    return Hx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "def condEntropy(x,y):\n",
    "    \n",
    "    \"\"\"\n",
    "    conditional entropy, or the average entropy of x given each y, or Hxy\n",
    "    1) For all Y {i=1:numel(X)}, compute the entropy of X given each Y\n",
    "    2) Multiply H(X|Y==i) with the probability of each Y (i.e. pxi)\n",
    "    3) Sum over all i\n",
    "\n",
    "    Args:\n",
    "        x (int): array of discrete values\n",
    "        y (int): array of discrete values\n",
    "        \n",
    "    Returns:\n",
    "        Hxy (float): average conditional entropy of x given y\n",
    "\n",
    "    \"\"\"\n",
    "\n",
    "    Hxy=0\n",
    "    uniquex = np.unique(x)\n",
    "    uniquey = np.unique(y)\n",
    "\n",
    "    # loop over unique elements of y\n",
    "    for i in np.arange(len(uniquey)): \n",
    "\n",
    "        # probability that y==y(i) (prob of each y)\n",
    "        py = np.sum(y==uniquey[i]) / N\n",
    "\n",
    "        # then loop over all possible x's to compute entropy of x at each y\n",
    "        tmp=0\n",
    "        for j in np.arange(len(uniquex)):\n",
    "            px_y = np.sum((x==uniquex[j]) & (y==uniquey[i])) / np.sum(y==uniquey[i])    # e.g. prob x==1 when y==0\n",
    "            tmp += (-( px_y * np.log2(px_y) ))                                     # entropy      \n",
    "\n",
    "        # then tally up entropy of x given each specific y multiplied by the probability of that y (py)\n",
    "        Hxy += py*tmp\n",
    "\n",
    "    return Hxy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MI is:  0.0002369146743368411\n"
     ]
    }
   ],
   "source": [
    "N = 1000   # number of data points\n",
    "Hxt = entropy(x=x)\n",
    "Hxyt = condEntropy(x=x,y=y)\n",
    "MIt = Hxt-Hxyt\n",
    "print('MI is: ', Hxt-Hxyt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Now repeat the above operations, but shuffle the data arrays and repeat the analysis many times (~500-1000 times). Plot the distribution of MI values that you get."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAEshJREFUeJzt3X+QXWV9x/H3VwL4gx/hx8JgEgyO6w/sDD/c0nSc6aixSqBj0im00CppJm06ljo6OK1Rp9X++ANaRyxV42REWRwUIuok1UilAWprG3SRGMBIWRHJNkhWgVhl1KLf/nGelctms/fu7t29N8+8XzN37jnPee55vvdO7ueefe45N5GZSJLq9axeFyBJml8GvSRVzqCXpMoZ9JJUOYNekipn0EtS5Qx6SaqcQS9JlTPoJalyi3pdAMDJJ5+cy5cv73UZknRYueuuu76fmQPt+vVF0C9fvpyRkZFelyFJh5WI+G4n/Zy6kaTKGfSSVDmDXpIqZ9BLUuUMekmqnEEvSZUz6CWpcga9JFXOoJekyvXFlbFzsXzjF3o29kNXXtizsSWpUx7RS1Ll2gZ9RLwkIna13H4YEW+LiBMj4taIeKDcn1D6R0RcExGjEbE7Is6d/6chSTqUtkGfmfdn5tmZeTbwCuBJ4HPARmBHZg4CO8o6wCpgsNw2AJvmo3BJUmdmOnWzEvh2Zn4XWA0Ml/ZhYE1ZXg1cn42dwOKIOK0r1UqSZmymQX8J8KmyfGpmPgJQ7k8p7UuAvS2PGSttzxARGyJiJCJGxsfHZ1iGJKlTHQd9RBwFvAH4dLuuU7TlQQ2ZmzNzKDOHBgba/m6+JGmWZnJEvwr4emY+WtYfnZiSKff7S/sYsKzlcUuBfXMtVJI0OzMJ+kt5etoGYBuwtiyvBba2tF9Wzr5ZARyYmOKRJC28ji6YiojnAr8J/ElL85XAlohYDzwMXFzatwMXAKM0Z+is61q1kqQZ6yjoM/NJ4KRJbT+gOQtnct8ELu9KdZKkOfPKWEmqnEEvSZUz6CWpcga9JFXOoJekyhn0klQ5g16SKmfQS1LlDHpJqpxBL0mVM+glqXIGvSRVzqCXpMoZ9JJUOYNekipn0EtS5Qx6SaqcQS9JlTPoJalyHQV9RCyOiJsj4lsRsScifj0iToyIWyPigXJ/QukbEXFNRIxGxO6IOHd+n4IkaTqdHtH/I3BLZr4UOAvYA2wEdmTmILCjrAOsAgbLbQOwqasVS5JmpG3QR8RxwG8A1wJk5s8y8wlgNTBcug0Da8ryauD6bOwEFkfEaV2vXJLUkU6O6F8IjAMfj4i7I+KjEfE84NTMfASg3J9S+i8B9rY8fqy0SZJ6oJOgXwScC2zKzHOAH/P0NM1UYoq2PKhTxIaIGImIkfHx8Y6KlSTNXCdBPwaMZeadZf1mmuB/dGJKptzvb+m/rOXxS4F9k3eamZszcygzhwYGBmZbvySpjbZBn5nfA/ZGxEtK00rgm8A2YG1pWwtsLcvbgMvK2TcrgAMTUzySpIW3qMN+bwFuiIijgAeBdTQfElsiYj3wMHBx6bsduAAYBZ4sfSVJPdJR0GfmLmBoik0rp+ibwOVzrEuS1CVeGStJlTPoJalyBr0kVc6gl6TKGfSSVDmDXpIqZ9BLUuUMekmqnEEvSZUz6CWpcga9JFXOoJekyhn0klQ5g16SKmfQS1LlDHpJqpxBL0mVM+glqXIGvSRVzqCXpMp1FPQR8VBE3BMRuyJipLSdGBG3RsQD5f6E0h4RcU1EjEbE7og4dz6fgCRpejM5on91Zp6dmUNlfSOwIzMHgR1lHWAVMFhuG4BN3SpWkjRzc5m6WQ0Ml+VhYE1L+/XZ2AksjojT5jCOJGkOOg36BL4UEXdFxIbSdmpmPgJQ7k8p7UuAvS2PHSttkqQeWNRhv1dm5r6IOAW4NSK+NU3fmKItD+rUfGBsADj99NM7LEOSNFMdHdFn5r5yvx/4HHAe8OjElEy531+6jwHLWh6+FNg3xT43Z+ZQZg4NDAzM/hlIkqbVNugj4nkRcezEMvA64F5gG7C2dFsLbC3L24DLytk3K4ADE1M8kqSF18nUzanA5yJiov8nM/OWiPgasCUi1gMPAxeX/tuBC4BR4ElgXderliR1rG3QZ+aDwFlTtP8AWDlFewKXd6U6SdKceWWsJFXOoJekyhn0klQ5g16SKmfQS1LlDHpJqpxBL0mVM+glqXIGvSRVzqCXpMoZ9JJUOYNekipn0EtS5Qx6SaqcQS9JlTPoJalyBr0kVc6gl6TKGfSSVDmDXpIq13HQR8QREXF3RHy+rJ8REXdGxAMRcVNEHFXajy7ro2X78vkpXZLUiZkc0b8V2NOyfhVwdWYOAo8D60v7euDxzHwRcHXpJ0nqkY6CPiKWAhcCHy3rAbwGuLl0GQbWlOXVZZ2yfWXpL0nqgU6P6D8A/AXwi7J+EvBEZj5V1seAJWV5CbAXoGw/UPpLknqgbdBHxG8B+zPzrtbmKbpmB9ta97shIkYiYmR8fLyjYiVJM9fJEf0rgTdExEPAjTRTNh8AFkfEotJnKbCvLI8BywDK9uOBxybvNDM3Z+ZQZg4NDAzM6UlIkg6tbdBn5jszc2lmLgcuAW7LzD8AbgcuKt3WAlvL8rayTtl+W2YedEQvSVoYczmP/h3AFRExSjMHf21pvxY4qbRfAWycW4mSpLlY1L7L0zLzDuCOsvwgcN4UfX4CXNyF2iRJXeCVsZJUOYNekipn0EtS5Qx6SaqcQS9JlTPoJalyBr0kVc6gl6TKGfSSVDmDXpIqZ9BLUuUMekmqnEEvSZUz6CWpcga9JFXOoJekyhn0klQ5g16SKmfQS1LlDHpJqlzboI+IZ0fEVyPiGxFxX0T8dWk/IyLujIgHIuKmiDiqtB9d1kfL9uXz+xQkSdPp5Ij+p8BrMvMs4Gzg/IhYAVwFXJ2Zg8DjwPrSfz3weGa+CLi69JMk9UjboM/Gj8rqkeWWwGuAm0v7MLCmLK8u65TtKyMiulaxJGlGOpqjj4gjImIXsB+4Ffg28ERmPlW6jAFLyvISYC9A2X4AOGmKfW6IiJGIGBkfH5/bs5AkHVJHQZ+ZP8/Ms4GlwHnAy6bqVu6nOnrPgxoyN2fmUGYODQwMdFqvJGmGZnTWTWY+AdwBrAAWR8SismkpsK8sjwHLAMr244HHulGsJGnmOjnrZiAiFpfl5wCvBfYAtwMXlW5rga1leVtZp2y/LTMPOqKXJC2MRe27cBowHBFH0HwwbMnMz0fEN4EbI+LvgLuBa0v/a4FPRMQozZH8JfNQtySpQ22DPjN3A+dM0f4gzXz95PafABd3pTpJ0px5ZawkVc6gl6TKGfSSVDmDXpIqZ9BLUuUMekmqnEEvSZUz6CWpcga9JFXOoJekyhn0klQ5g16SKmfQS1LlDHpJqpxBL0mVM+glqXIGvSRVzqCXpMoZ9JJUOYNekirXNugjYllE3B4ReyLivoh4a2k/MSJujYgHyv0JpT0i4pqIGI2I3RFx7nw/CUnSoXVyRP8U8PbMfBmwArg8Is4ENgI7MnMQ2FHWAVYBg+W2AdjU9aolSR1rG/SZ+Uhmfr0s/y+wB1gCrAaGS7dhYE1ZXg1cn42dwOKIOK3rlUuSOjKjOfqIWA6cA9wJnJqZj0DzYQCcUrotAfa2PGystE3e14aIGImIkfHx8ZlXLknqSMdBHxHHAJ8B3paZP5yu6xRteVBD5ubMHMrMoYGBgU7LkCTNUEdBHxFH0oT8DZn52dL86MSUTLnfX9rHgGUtD18K7OtOuZKkmerkrJsArgX2ZOb7WzZtA9aW5bXA1pb2y8rZNyuAAxNTPJKkhbeogz6vBN4E3BMRu0rbu4ArgS0RsR54GLi4bNsOXACMAk8C67pasSRpRtoGfWb+B1PPuwOsnKJ/ApfPsS5JUpd4ZawkVc6gl6TKGfSSVDmDXpIqZ9BLUuUMekmqnEEvSZXr5IIpHcLyjV/oybgPXXlhT8aVdHjyiF6SKmfQS1LlDHpJqpxBL0mVM+glqXIGvSRVzqCXpMoZ9JJUOYNekipn0EtS5Qx6Sapc26CPiI9FxP6IuLel7cSIuDUiHij3J5T2iIhrImI0InZHxLnzWbwkqb1OjuivA86f1LYR2JGZg8COsg6wChgstw3Apu6UKUmarbZBn5lfBh6b1LwaGC7Lw8Calvbrs7ETWBwRp3WrWEnSzM12jv7UzHwEoNyfUtqXAHtb+o2VNklSj3T7y9iYoi2n7BixISJGImJkfHy8y2VIkibMNugfnZiSKff7S/sYsKyl31Jg31Q7yMzNmTmUmUMDAwOzLEOS1M5sg34bsLYsrwW2trRfVs6+WQEcmJjikST1Rtv/SjAiPgW8Cjg5IsaA9wBXAlsiYj3wMHBx6b4duAAYBZ4E1s1DzZKkGWgb9Jl56SE2rZyibwKXz7UoSVL3eGWsJFXOoJekyhn0klQ5g16SKmfQS1LlDHpJqpxBL0mVM+glqXJtL5hS/1m+8Qs9G/uhKy/s2diSZscjekmqnEEvSZUz6CWpcga9JFXOoJekyhn0klQ5T6/UjPTq1E5P65RmzyN6SaqcR/Q6LPiXhDR7HtFLUuUMekmq3LwEfUScHxH3R8RoRGycjzEkSZ3petBHxBHAh4BVwJnApRFxZrfHkSR1Zj6+jD0PGM3MBwEi4kZgNfDNeRhLmle9/KXQXvEL6PrMR9AvAfa2rI8BvzYP40iaB/4Mdn3mI+hjirY8qFPEBmBDWf1RRNw/D7UcysnA9xdwvJno59rA+uaqn+vreW1x1bSbe15fG72o7wWddJqPoB8DlrWsLwX2Te6UmZuBzfMwflsRMZKZQ70Yu51+rg2sb676ub5+rg2sby7m46ybrwGDEXFGRBwFXAJsm4dxJEkd6PoRfWY+FRF/BvwLcATwscy8r9vjSJI6My8/gZCZ24Ht87HvLunJlFGH+rk2sL656uf6+rk2sL5Zi8yDvieVJFXEn0CQpNpl5mFxA84H7gdGgY1TbD8auKlsvxNY3rLtnaX9fuD17fYJ3FDa7wU+Bhw5aaxfBX4OXNRv9QGvAnYB9wH/1i+1AccD/wx8o9S2rkev3bWlht3AzcAxHYzRD/VdQXPR4W5gB/CCfqqvZftFNKdTD/VTbcDvltfvPuCT/fTaAacDtwN3l20XdD0/u73D+bjRfKn7beCFwFHlxTpzUp8/BT5Sli8BbirLZ5b+RwNnlP0cMd0+gQtorgcI4FPAmyfVchvNdxAX9VN9wOLyj/n0sn5KH9X2LuCqsjwAPFYeu9D1Hdey3/dT3ojTjNEv9b0aeG5ZfnO/1VfWjwW+DOwEhvqlNmCQJkRPmHhf9NNrRzO3/+aW/T7U7Qw9XKZufvmzCpn5M2DiZxVarQaGy/LNwMqIiNJ+Y2b+NDO/Q/Mpe950+8zM7VkAX6W5FmDCW4DPAPv7sL7fBz6bmQ+Xfvv7qLYEji37PYYm6J/qQX0/BCiPfw5PX8x3qDH6or7MvD0znyxj7Gx5XfuivuJvgb8HftJntf0x8KHMfLz0m3jv9kt9CRxXlo9niuuO5upwCfqpflZhyaH6ZOZTwAHgpGke23afEXEk8CbglrK+BPht4CP9WB/wYuCEiLgjIu6KiMv6qLYPAi+j+Ud8D/DWzPxFL+qLiI8D3wNeCvzTLMdY6PparQe+OHnsXtYXEecAyzLz81ON28vaaN4XL46Ir0TEzog4v8/qey/wxogYo5kpeAtddrgEfSc/q3CoPjNtb/Vh4MuZ+e9l/QPAOzLz531a3yLgFcCFwOuBvwSe3ye1vZ7mu4PnA2cDH4yI4zrcV1fry8x1pY49wO/NcoyFrq8ZKOKNNNMi/9Bm7AWrLyKeBVwNvL3DcResttK8iGb65lXApcBHI2JxH9V3KXBdZi6lmfr8RHlNu+ZwCfpOflbhl30iYhHNn0CPTfPYafcZEe+hmUu+oqXPEHBjRDxE86XThyNiTR/VNwbckpk/zszv08yXHtsnta2jmVbKzBwFvkNzVLPgrx1A+bC+CfidWY6x0PUREa8F3g28ITN/OnnsHtZ3LPArwB3lvbGC5mr4Y/qgtokxtmbm/5Vplvtpgr8fXjto/kLbUrb9F/Bsmt/N6Z5OJ/N7eaP5RH6Q5kuPiS84Xj6pz+U880uTLWX55TzzS5MHab4wOeQ+gT8C/hN4zjQ1XcfTX8b2RX00UyM7ymOfS3Pmy1l9Utsm4L1l+VTgf2j+MS/Ya0dztPWi8tgA3ge8r80Y/VLfOTRf8g326r0xXX2TxruD5qCoL2qjOQtmuCyfTDO1clIf1fdF4A9b3sP7KNc4dS1Dex3iHRfa/Enz3+Uf+7tL29/QHN1A8yn4aZovRb4KvLDlse8uj7sfWDXdPkv7U6VtV7n91RT1XMczT6/si/qAP6c58+Ze4G39UhvNn6tfopmfvxd440K/djR/wX6lpYYbKGdCtBmjH+r7V+DRltd1Wz/VN+m9cQdPn17Z89pogvX9NO+Le4BL+um1oznT5is0Hwq7gNd1Oz+9MlaSKne4zNFLkmbJoJekyhn0klQ5g16SKmfQS1LlDHpJqpxBL0mVM+glqXL/D5of25AG13wYAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1c0a4a2a588>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "num_shuffle = 1000\n",
    "tmp0 = np.zeros(N)\n",
    "tmp1 = np.zeros(N)\n",
    "Hx = np.zeros(num_shuffle)\n",
    "Hxy = np.zeros(num_shuffle)\n",
    "MI = np.zeros(num_shuffle)\n",
    "\n",
    "for i in np.arange(num_shuffle):\n",
    "    for j in np.arange(N):   \n",
    "        if np.random.rand(1) < .5:\n",
    "            tmp0[j] = x[j]\n",
    "            tmp1[j] = y[j]\n",
    "        else:\n",
    "            tmp0[j] = y[j]\n",
    "            tmp1[j] = x[j]\n",
    "    Hx[i] = entropy(x=tmp0)\n",
    "    Hxy[i] = condEntropy(x=tmp0,y=tmp1)\n",
    "    MI[i] = Hx[i] - Hxy[i]\n",
    "\n",
    "#MI = Hx - Hxy\n",
    "\n",
    "plt.hist(MI)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Now subtract the mean of the shuffled MI values from your 'real' MI value...this will help correct for any bias that is introduced by a limited sample size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-1.0329021582514803e-05"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "MIt-np.mean(MI)"
   ]
  }
 ],
 "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}