code/hypothesis-testing_1.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import scipy.stats\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Hypothesis Testing Introduction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- the task of determining the probability that a hypothesis is true\n",
    "- typically, we define a null hypothesis and an alternative hypothesis; then, we determine the probability that the null hypothesis is true\n",
    "- a p-value is a test statistic that can be regarded as evidence against the null hypothesis (for instance, we reject a null hypothesis if the p-value is smaller than a specified threshold or significance level or alpha value) -- the p-value is the probability of a model assuming that the null hypothesis is true\n",
    "- 1-sample test: compare a sample against a number\n",
    "- 2-sample test: compare two samples"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Case 1 -- Simulation under the null hypothesis using permutation sampling (2-sample test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- null hypothesis: data is identically distributed between the two variables\n",
    "- how likely is it to observe this data if the null hypothesis is true?\n",
    "- we can scramble (permute the data) and divide it into 2 new samples and then check if there's any difference (in other words, we generate new data under the null hypothesis)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "df = pd.read_csv('../datasets/iris/iris.csv')\n",
    "x1 = df.loc[df['species'] == 'setosa', 'sepal_width'].values\n",
    "x2 = df.loc[df['species'] == 'versicolor', 'sepal_width'].values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mean: 3.428\n",
      "Std.Dev.: 0.379064369096\n",
      "Mean: 2.77\n",
      "Std.Dev.: 0.313798323378\n"
     ]
    }
   ],
   "source": [
    "for d in (x1, x2):\n",
    "    print('Mean:', np.mean(d))\n",
    "    print('Std.Dev.:', np.std(d, ddof=1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- the difference in means (our test statistic) of permuted samples should be zero if the null hypothesis is true"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def generate_permutation_samples(sample1, sample2, seed=None):\n",
    "    rng = np.random.RandomState(seed)\n",
    "    sample = np.concatenate((sample1, sample2))\n",
    "    sample_perm = rng.permutation(sample)\n",
    "    sample1_perm = sample_perm[:sample1.shape[0]]\n",
    "    sample2_perm = sample_perm[sample1.shape[0]:]\n",
    "    return sample1_perm, sample2_perm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "num_permutations = 1000\n",
    "\n",
    "perm_mean_diffs = np.zeros(num_permutations)\n",
    "for i in range(num_permutations):\n",
    "    x1_perm, x2_perm = generate_permutation_samples(x1, x2, seed=i)\n",
    "    perm_mean_diffs[i] = np.mean(x1_perm) - np.mean(x2_perm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.65799999999999947"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "empirical_mean_diff = np.mean(x1) - np.mean(x2)\n",
    "empirical_mean_diff"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.lines.Line2D at 0x115265518>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAD/tJREFUeJzt3XuMHWd9xvHvUxtSAm3j1Itr4tB1K5cqIBDpNk0BVWkD\nJSQIpxKKTLm41JKFmnKpkKgDEvkrkmkrSis1VFaS4qpRUiukjdVwC4Y0qiCBTQjkYkJccnPqxMul\nQEMV6vDrHzuqtmbtPT5zzh779fcjrc7MO+/s/M5o/fjdd2fmpKqQJLXrpyZdgCRpvAx6SWqcQS9J\njTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuNWTroAgNWrV9f09PSky5Ck5fX9B+Zff/ZFQ+1+\n5513fquqppbqd1wE/fT0NLOzs5MuQ5KW12fPm3999a1D7Z7kkUH6OXUjSY0z6CWpcQa9JDXOoJek\nxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNOy7ujNWJa3rbzQP1e3j7RWOuRNKROKKXpMYZ9JLUOINe\nkhpn0EtS45YM+iTXJDmY5N5Ftr03SSVZvaDtsiT7kjyQ5LWjLliSdGwGGdF/DLjg8MYkZwK/Czy6\noO0sYBPw4m6fK5OsGEmlkqShLBn0VXUb8J1FNv0l8D6gFrRtBK6vqqer6iFgH3DOKAqVJA1nqDn6\nJBuBx6vqq4dtOgN4bMH6/q5tse+xNclsktm5ublhypAkDeCYgz7JqcD7gQ/2OXBV7aiqmaqamZpa\n8iMPJUlDGubO2F8G1gNfTQKwDrgryTnA48CZC/qu69okSRNyzCP6qrqnqp5fVdNVNc389MzZVfUE\nsBvYlOSUJOuBDcCXRlqxJOmYDHJ55XXAF4EXJdmfZMuR+lbVfcAu4H7gU8ClVfXMqIqVJB27Jadu\nqupNS2yfPmz9CuCKfmVJkkbFO2MlqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPo\nJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxg3y4eDXJDmY\n5N4FbX+e5OtJvpbkn5KctmDbZUn2JXkgyWvHVbgkaTCDjOg/BlxwWNstwEuq6qXAN4DLAJKcBWwC\nXtztc2WSFSOrVpJ0zJYM+qq6DfjOYW2fqapD3ertwLpueSNwfVU9XVUPAfuAc0ZYryTpGI1ijv4P\ngU92y2cAjy3Ytr9rkyRNSK+gT/IB4BBw7RD7bk0ym2R2bm6uTxmSpKMYOuiT/AHweuDNVVVd8+PA\nmQu6revafkJV7aiqmaqamZqaGrYMSdIShgr6JBcA7wPeUFU/XLBpN7ApySlJ1gMbgC/1L1OSNKyV\nS3VIch1wHrA6yX7gcuavsjkFuCUJwO1V9Y6qui/JLuB+5qd0Lq2qZ8ZVvCRpaUsGfVW9aZHmq4/S\n/wrgij5FSZJGxztjJalxBr0kNc6gl6TGLTlHr5PX9LabJ12CpBFwRC9JjTPoJalxBr0kNc6gl6TG\nGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY3zztjGDHI368PbL1qGSiQdLxzRS1LjDHpJapxBL0mNM+gl\nqXEGvSQ1bpAPB78GeD1wsKpe0rWdDvwjMA08DFxSVd/ttl0GbAGeAd5VVZ8eS+VqklcNSaM3yIj+\nY8AFh7VtA/ZU1QZgT7dOkrOATcCLu32uTLJiZNVKko7ZkkFfVbcB3zmseSOws1veCVy8oP36qnq6\nqh4C9gHnjKhWSdIQhr1hak1VHeiWnwDWdMtnALcv6Le/a9NJzo8llCan9x9jq6qAOtb9kmxNMptk\ndm5urm8ZkqQjGDbon0yyFqB7Pdi1Pw6cuaDfuq7tJ1TVjqqaqaqZqampIcuQJC1l2KDfDWzuljcD\nNy1o35TklCTrgQ3Al/qVKEnqY5DLK68DzgNWJ9kPXA5sB3Yl2QI8AlwCUFX3JdkF3A8cAi6tqmfG\nVLskaQBLBn1VvekIm84/Qv8rgCv6FCVJGh3vjJWkxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+gl\nqXEGvSQ1zqCXpMYN+5hincB8ZLB0cnFEL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6\nSWqcQS9JjTPoJalxvYI+yZ8kuS/JvUmuS/LTSU5PckuSB7vXVaMqVpJ07IYO+iRnAO8CZqrqJcAK\nYBOwDdhTVRuAPd26JGlC+k7drASek2QlcCrwH8BGYGe3fSdwcc9jSJJ6GDroq+px4C+AR4EDwPeq\n6jPAmqo60HV7Aliz2P5JtiaZTTI7Nzc3bBmSpCX0mbpZxfzofT3wAuC5Sd6ysE9VFVCL7V9VO6pq\npqpmpqamhi1DkrSEPlM3rwYeqqq5qvof4EbgFcCTSdYCdK8H+5cpSRpWnw8eeRQ4N8mpwH8D5wOz\nwFPAZmB793pT3yLlh4VIGt7QQV9VdyS5AbgLOAR8BdgBPA/YlWQL8AhwySgKlSQNp9dHCVbV5cDl\nhzU/zfzoXpJ0HPDOWElqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mN\nM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxvYI+yWlJbkjy9SR7k/xm\nktOT3JLkwe511aiKlSQdu74j+r8CPlVVvwq8DNgLbAP2VNUGYE+3LkmakKGDPsnPAb8FXA1QVT+q\nqv8ENgI7u247gYv7FilJGl6fEf16YA74uyRfSXJVkucCa6rqQNfnCWBN3yIlScPrE/QrgbOBj1bV\ny4GnOGyapqoKqMV2TrI1yWyS2bm5uR5lSJKOpk/Q7wf2V9Ud3foNzAf/k0nWAnSvBxfbuap2VNVM\nVc1MTU31KEOSdDRDB31VPQE8luRFXdP5wP3AbmBz17YZuKlXhZKkXlb23P+dwLVJng18E3g78/95\n7EqyBXgEuKTnMZo2ve3mSZcgqXG9gr6q7gZmFtl0fp/vK0kaHe+MlaTGGfSS1DiDXpIaZ9BLUuMM\neklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1Lj+j69Ulp2gz7x8+HtF425EunE4Ihe\nkhpn0EtS45y6UbMGmeJxekcnA0f0ktQ4g16SGmfQS1Ljegd9khVJvpLkX7r105PckuTB7nVV/zIl\nScMaxYj+3cDeBevbgD1VtQHY061LkiakV9AnWQdcBFy1oHkjsLNb3glc3OcYkqR++o7oPwK8D/jx\ngrY1VXWgW34CWNPzGJKkHoYO+iSvBw5W1Z1H6lNVBdQR9t+aZDbJ7Nzc3LBlSJKW0GdE/0rgDUke\nBq4HfifJPwBPJlkL0L0eXGznqtpRVTNVNTM1NdWjDEnS0Qwd9FV1WVWtq6ppYBPwuap6C7Ab2Nx1\n2wzc1LtKSdLQxnEd/XbgNUkeBF7drUuSJmQkz7qpqluBW7vlbwPnj+L7SpL6885YSWqcQS9JjTPo\nJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDVuJA810+Kmt9086RIk\nyRG9JLXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapyXVw7ByyYlnUiGHtEnOTPJ55Pcn+S+JO/u2k9P\nckuSB7vXVaMrV5J0rPpM3RwC3ltVZwHnApcmOQvYBuypqg3Anm5dkjQhQwd9VR2oqru65R8Ae4Ez\ngI3Azq7bTuDivkVKkoY3kj/GJpkGXg7cAaypqgPdpieANUfYZ2uS2SSzc3NzoyhDkrSI3kGf5HnA\nx4H3VNX3F26rqgJqsf2qakdVzVTVzNTUVN8yJElH0CvokzyL+ZC/tqpu7JqfTLK2274WONivRElS\nH32uuglwNbC3qj68YNNuYHO3vBm4afjyJEl99bmO/pXAW4F7ktzdtb0f2A7sSrIFeAS4pF+JkqQ+\nhg76qvo3IEfYfP6w31eSNFo+AkGSGmfQS1LjDHpJapwPNdNJbdAH1D28/aIxVyKNjyN6SWqcQS9J\njXPqRhrAIFM8Tu/oeGXQH8YPFZHUGqduJKlxBr0kNc6gl6TGnTRz9M69SzpZOaKXpMYZ9JLUOINe\nkhpn0EtS406aP8ZKxwsfpKblZtBLI+KVXTpeOXUjSY0bW9AnuSDJA0n2Jdk2ruNIko5uLFM3SVYA\nfwO8BtgPfDnJ7qq6fxzH81dmtcgnZmpUxjWiPwfYV1XfrKofAdcDG8d0LEnSUYwr6M8AHluwvr9r\nkyQts4lddZNkK7C1W/2vJA9MqpZjsBr41qSLmDDPwXF0DvKhiRz2uHn/EzTic5Bhd/zFQTqNK+gf\nB85csL6ua/s/VbUD2DGm449Fktmqmpl0HZPkOfAcnOzvH068czCuqZsvAxuSrE/ybGATsHtMx5Ik\nHcVYRvRVdSjJHwOfBlYA11TVfeM4liTp6MY2R19VnwA+Ma7vPyEn1FTTmHgOPAcn+/uHE+wcpKom\nXYMkaYx8BIIkNc6gP4okpye5JcmD3euqRfqcmeTzSe5Pcl+Sd0+i1lFa6vEVmffX3favJTl7EnWO\n0wDn4M3de78nyReSvGwSdY7ToI8xSfLrSQ4leeNy1rccBjkHSc5Lcnf37/9fl7vGgVSVX0f4Av4M\n2NYtbwM+tEiftcDZ3fLPAN8Azpp07T3e8wrg34FfAp4NfPXw9wNcCHyS+Yt/zwXumHTdEzgHrwBW\ndcuvOxnPwYJ+n2P+73FvnHTdE/g5OA24H3hht/78Sde92Jcj+qPbCOzslncCFx/eoaoOVNVd3fIP\ngL2c2HcBD/L4io3A39e824HTkqxd7kLHaMlzUFVfqKrvdqu3M3+vSEsGfYzJO4GPAweXs7hlMsg5\n+H3gxqp6FKCqjsvzYNAf3ZqqOtAtPwGsOVrnJNPAy4E7xlvWWA3y+IrWH3FxrO9vC/O/4bRkyXOQ\n5Azg94CPLmNdy2mQn4NfAVYluTXJnUnetmzVHYOT/oNHknwW+IVFNn1g4UpVVZIjXqKU5HnMj2ze\nU1XfH22VOl4l+W3mg/5Vk65lAj4C/GlV/TgZ+hb+E91K4NeA84HnAF9McntVfWOyZf1/J33QV9Wr\nj7QtyZNJ1lbVgW5qYtFfy5I8i/mQv7aqbhxTqctlycdXDNjnRDbQ+0vyUuAq4HVV9e1lqm25DHIO\nZoDru5BfDVyY5FBV/fPylDh2g5yD/cC3q+op4KkktwEvY/5vdccNp26ObjewuVveDNx0eIfM/5Rf\nDeytqg8vY23jMsjjK3YDb+uuvjkX+N6CKa4WLHkOkrwQuBF46/E2ehuRJc9BVa2vqumqmgZuAP6o\noZCHwf4t3AS8KsnKJKcCv8H83+mOKyf9iH4J24FdSbYAjwCXACR5AXBVVV0IvBJ4K3BPkru7/d5f\n83cGn3DqCI+vSPKObvvfMn+FxYXAPuCHwNsnVe84DHgOPgj8PHBlN6I9VCfQQ66WMuA5aNog56Cq\n9ib5FPA14MfM58K9k6t6cd4ZK0mNc+pGkhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS\n1Lj/BSfsNnZ1y1KaAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1152655c0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(perm_mean_diffs, bins=20)\n",
    "plt.axvline(empirical_mean_diff, color='orange')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "P-value 0.0\n"
     ]
    }
   ],
   "source": [
    "p_value = np.sum(perm_mean_diffs >= empirical_mean_diff) / len(perm_mean_diffs)\n",
    "print('P-value', p_value)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- the area under the histogram left (or right) to the vertical bar (difference observed in the original samples) is the p-value\n",
    "- here, the area right to the vertical bar (i.e., values that are at least as extreme as the observed difference of the mean between the two samples) is 0.0, so we can conclude that the probability of observing this data (the difference between the means of the samples) or something more extreme is 0% likely if the null hypothesis was true (remember that the p-value is not the probability that the null hypothesis is true, though)\n",
    "- thus, we can reject the null hypothesis (data is \"statistically significantly different\" from what we would expect to observe if the null hypothesis was true)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Case 2 -- Simulation under the null hypothesis using bootstrapping (1-sample test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.4279999999999999"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df = pd.read_csv('../datasets/iris/iris.csv')\n",
    "x1 = df.loc[df['species'] == 'setosa', 'sepal_width'].values\n",
    "x1_mean = np.mean(x1)\n",
    "x1_mean"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- our null hypothesis is that two variables, x1 and x2, have the same mean (not necessarily same distribution)\n",
    "- assume we don't have a sample from x2, but we have its mean\n",
    "- we set the significance level $\\alpha=0.05$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "x2_mean = 3.55"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "p-value 0.006\n"
     ]
    }
   ],
   "source": [
    "# 1) normalize x1 so that it has the same mean as x2\n",
    "normed_x1 = x1 - np.mean(x1) + x2_mean\n",
    "\n",
    "# 2) draw bootstrap replicates\n",
    "def bootstrap_means(x, n_bootstrap_samples, seed=None):\n",
    "\n",
    "    rng = np.random.RandomState(seed)\n",
    "    sample_means = np.zeros(shape=n_bootstrap_samples)\n",
    "    \n",
    "    for i in range(n_bootstrap_samples):\n",
    "        boot_sample = rng.choice(x, size=x.shape[0], replace=True)\n",
    "        # replicate is a general term for a statistic computed\n",
    "        # from a bootstrap sample\n",
    "        bootstrap_replicate = np.mean(boot_sample)\n",
    "        sample_means[i] = bootstrap_replicate\n",
    "    return sample_means\n",
    "\n",
    "num_replicates = 1000\n",
    "bs_replicates = bootstrap_means(x1, n_bootstrap_samples=num_replicates, seed=123)\n",
    "\n",
    "# 3) compute p-value\n",
    "p_value = np.sum(bs_replicates >= x2_mean) / num_replicates\n",
    "\n",
    "print('p-value', p_value)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- at $\\alpha=0.05$, we can reject the null hypothesis that the two variables are from the same distribtion"
   ]
  }
 ],
 "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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}