Chapter_9.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy import stats\n",
    "from matplotlib import pyplot as plt\n",
    "from statsmodels.sandbox.stats.runs import mcnemar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 9.1 Fisher's Exact Test: The Tea Experiment\n",
    "\n",
    "# At a party, a lady claimed to be able to tell whether the tea or the milk was added\n",
    "# first to a cup. Fisher proposed to give her eight cups, four of each variety, in random\n",
    "# order. One could then ask what the probability was for her getting the number she\n",
    "# got correct, but just by chance.\n",
    "\n",
    "# The experiment provided the Lady with eight randomly ordered cups of tea—\n",
    "# four prepared by first adding milk, four prepared by first adding the tea. She was to\n",
    "# select the four cups prepared by one method. (This offered the Lady the advantage\n",
    "# of judging cups by comparison.)\n",
    "\n",
    "# The null hypothesis was that the Lady had no such ability. (In the real, historical\n",
    "# experiment, the lady got all eight cups correct.)\n",
    "\n",
    "# Calculate if the claim of the Lady is supported if she gets three out of the four\n",
    "# pairs correct."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The hypergeometric distribution models drawing objects from a bin.\n",
    "# M - total number of cups of tea\n",
    "# n - number of cups of tea with milk poured first\n",
    "# N - number of cups of tea that the lady chooses\n",
    "M, n, N = 8, 4, 4\n",
    "hd = stats.hypergeom(M, n, N)\n",
    "x = np.arange(0, N+1)\n",
    "pmf_tea = hd.pmf(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAGwCAYAAABVdURTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA4s0lEQVR4nO3deXxU1f3/8feQlSUJewgQ2QkoypKgBIzKFr/IA1DU0mJZKvgVlEJIK4LoA6Fq1C9LsC1QZPtSFVBBa1ssRGQJIlRiIltElkBiSMomCYoGkpzfH/llvg5ZmMk2mcvr+XjMQ+6Zc+98Didh3t45c6/NGGMEAABgEXXcXQAAAEBVItwAAABLIdwAAABLIdwAAABLIdwAAABLIdwAAABLIdwAAABL8XZ3ATWtsLBQZ86cUUBAgGw2m7vLAQAATjDG6PLly2rZsqXq1Cn/3MxNF27OnDmj0NBQd5cBAAAqICMjQ61bty63z00XbgICAiQV/eUEBga6uRoAAOCM3NxchYaG2t/Hy3PThZvij6ICAwMJNwAAeBhnlpSwoBgAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFjKTXeFYgDWVFAgJSZKWVlSSIgUFSV5ebm7KgDuQLgB4PE2bZKmTZO+/fb/2lq3lhYvlkaOdF9dANyDj6UAeLRNm6RHHnEMNpKUmVnUvmmTe+oC4D6EGwAeq6Cg6IyNMSWfK26LiSnqB+DmQbgB4LESE0uesfk5Y6SMjKJ+AG4erLkB4LGyshy3A3qfVB3ffBVe9dblL9qX2Q+AtRFuAHiskBDH7cDeafIO+En5l/0dws31/QBYGx9LAfBYUVFF34qy2Up/3maTQkOL+gG4eRBuAHgsL6+ir3tLJQNO8XZ8PNe7AW42hBsAHm3kSOn996VWrRzbW7cuauc6N8DNh3ADwOONHCmdOiU1a1q03ayplJZGsAFuVm4PN0uWLFG7du3k7++v8PBwJZbznc0dO3bIZrOVeHz99dc1WDGA2sjLS/LzK/qznx8fRQE3M7eGmw0bNigmJkazZ89WcnKyoqKiNGTIEKWnp5e739GjR5WVlWV/dOrUqYYqBgAAtZ1bw83ChQs1YcIETZw4UV27dlV8fLxCQ0O1dOnScvdr3ry5WrRoYX94lfO/aHl5ecrNzXV4AAAA63JbuLl69aqSkpIUHR3t0B4dHa09e/aUu2/Pnj0VEhKigQMHavv27eX2jYuLU1BQkP0RGhpa6doBAEDt5bZwc/78eRUUFCg4ONihPTg4WNnZ2aXuExISouXLl2vjxo3atGmTwsLCNHDgQO3atavM15k1a5ZycnLsj4yMjCodBwAAqF3cfoVi23UXpzDGlGgrFhYWprCwMPt2ZGSkMjIyNH/+fN1zzz2l7uPn5ye/4lWGAADA8tx25qZp06by8vIqcZbm7NmzJc7mlKdPnz46duxYVZcHAAA8lNvCja+vr8LDw5WQkODQnpCQoL59+zp9nOTkZIVw4xgAAPD/ufVjqdjYWI0ZM0YRERGKjIzU8uXLlZ6erkmTJkkqWi+TmZmptWvXSpLi4+PVtm1b3Xbbbbp69areeustbdy4URs3bnTnMAAAQC3i1nAzatQoXbhwQfPmzVNWVpa6deumzZs3q02bNpKkrKwsh2veXL16Vb///e+VmZmpunXr6rbbbtM///lPPfDAA+4aAgAAqGVsxhjj7iJqUm5uroKCgpSTk6PAwEB3lwOgCvV5ZZuyc39Si0B/7X1uoLvLAVCFXHn/dvvtFwAAAKoS4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFiK28PNkiVL1K5dO/n7+ys8PFyJiYlO7ffZZ5/J29tbPXr0qN4CAQCAR3FruNmwYYNiYmI0e/ZsJScnKyoqSkOGDFF6enq5++Xk5Gjs2LEaOHBgDVUKAAA8hVvDzcKFCzVhwgRNnDhRXbt2VXx8vEJDQ7V06dJy93vyySc1evRoRUZG1lClAADAU7gt3Fy9elVJSUmKjo52aI+OjtaePXvK3G/16tU6ceKE5syZ49Tr5OXlKTc31+EBAACsy23h5vz58yooKFBwcLBDe3BwsLKzs0vd59ixY5o5c6befvtteXt7O/U6cXFxCgoKsj9CQ0MrXTsAAKi93L6g2GazOWwbY0q0SVJBQYFGjx6tuXPnqnPnzk4ff9asWcrJybE/MjIyKl0zAACovZw7/VENmjZtKi8vrxJnac6ePVvibI4kXb58Wfv371dycrKmTJkiSSosLJQxRt7e3tq6dasGDBhQYj8/Pz/5+flVzyAAAECt47YzN76+vgoPD1dCQoJDe0JCgvr27Vuif2BgoA4ePKiUlBT7Y9KkSQoLC1NKSoruuuuumiodAADUYm47cyNJsbGxGjNmjCIiIhQZGanly5crPT1dkyZNklT0kVJmZqbWrl2rOnXqqFu3bg77N2/eXP7+/iXaAQDAzcut4WbUqFG6cOGC5s2bp6ysLHXr1k2bN29WmzZtJElZWVk3vOYNAADAz9mMMcbdRdSk3NxcBQUFKScnR4GBge4uB0AV6vPKNmXn/qQWgf7a+xwX+QSsxJX3b7d/WwoAAKAqEW4AAIClEG4AAIClEG4AAIClEG4AAIClEG4AAIClEG4AAIClEG4AAIClEG4AAIClEG4AAIClEG4AAIClEG4AAIClEG4AAIClEG4AAIClEG4AAIClEG4AAIClEG4AAIClEG4AAIClEG4AAIClEG4AAIClEG4AAIClEG4AAIClEG4AAIClEG4AAIClEG4AAIClEG4AAIClEG4AAIClEG4AAIClEG4AAIClEG4AAIClEG4AAIClEG4AAIClOB1uPv30U+Xn51dnLQAAAJXmdLgZPHiwLl68aN/u06ePMjMzq6UoAACAinI63BhjHLYPHz6svLy8Ki8IAACgMlhzAwAALMXpcGOz2WSz2crcBgAAqA28ne1ojNHAgQPl7V20y5UrVzRs2DD5+vo69Pvyyy+rtkIAAAAXOB1u5syZ47A9YsSIKi8GAACgsiocbgAAAGojFhQDAABLcfrMzYABA5zq9+mnn1a4GAAAgMpyOtzs2LFDbdq00dChQ+Xj41OdNQEAAFSY0+Hm1Vdf1Zo1a/Tee+/pscce0+OPP65u3bpVZ20AAAAuc3rNzYwZM3TkyBF9+OGHunz5svr166c777xTy5YtU25ubnXWCAAA4DSXFxRHRkbqzTffVFZWlp5++mmtWrVKLVu2JOAAAIBaocLflvryyy+1c+dOpaamqlu3bqzDAQAAtYJL4ebMmTN65ZVX1LlzZz3yyCNq3Lix9u3bp71796pu3brVVSMAAIDTnF5Q/MADD2j79u2Kjo7W//zP/2jo0KH2WzEAAADUFjZjjHGmY506dRQSEqLmzZuXe8PM2n5vqdzcXAUFBSknJ0eBgYHuLgdAFerzyjZl5/6kFoH+2vvcQHeXA6AKufL+ze0XAACApRBuAACApbi0aGbfvn366KOPdO3aNQ0aNEjR0dHVVRcAAECFOB1uPvjgAz366KPy9/eXt7e3FixYoAULFigmJqYaywMAAHCN018Ff+WVVzR+/HhdunRJly5d0ty5c/XSSy9VuoAlS5aoXbt28vf3V3h4uBITE8vsu3v3bvXr109NmjRR3bp11aVLFy1atKjSNQAAAOtwOtwcPXpUM2bMsH/9+5lnntGlS5d0/vz5Cr/4hg0bFBMTo9mzZys5OVlRUVEaMmSI0tPTS+1fv359TZkyRbt27VJqaqqef/55Pf/881q+fHmFawAAANbi0lfBs7Oz1bx5c3tbQECAvvrqK7Vv375CL37XXXepV69eWrp0qb2ta9euevDBBxUXF+fUMUaOHKn69evrr3/9q1P9+So4YF18FRywrmr5KrgkbdmyRUFBQfbtwsJCbdu2TYcOHbK3DR8+3KljXb16VUlJSZo5c6ZDe3R0tPbs2ePUMZKTk7Vnz55yPx7Ly8tTXl6efZt7YAEAYG0uhZtx48aVaHvyySftf7bZbCooKHDqWOfPn1dBQYGCg4Md2oODg5WdnV3uvq1bt9a5c+eUn5+vF198URMnTiyzb1xcnObOnetUTQAAwPM5veamsLDwhg9ng83PXX+1Y2NMuVdAlqTExETt379fy5YtU3x8vNatW1dm31mzZiknJ8f+yMjIcLlGAADgOdx2c6imTZvKy8urxFmas2fPljibc7127dpJkm6//Xb95z//0Ysvvqhf/epXpfb18/OTn59f1RQNAABqPZfuCl6VfH19FR4eroSEBIf2hIQE9e3b1+njGGMc1tQAAICbm1tv6x0bG6sxY8YoIiJCkZGRWr58udLT0zVp0iRJRR8pZWZmau3atZKkP//5z7rlllvUpUsXSUXXvZk/f75++9vfum0MAACgdnFruBk1apQuXLigefPmKSsrS926ddPmzZvVpk0bSVJWVpbDNW8KCws1a9YspaWlydvbWx06dNCrr77qsKgZAADc3Jy+zo1VcJ0bwLq4zg1gXdV2nZvrff/99yosLHRoIzAAAAB3cnlBcVpamoYOHar69esrKChIjRo1UqNGjdSwYUM1atSoOmoEAABwmstnbh577DFJ0qpVqxQcHHzDa9IAAADUJJfDzYEDB5SUlKSwsLDqqAcAAKBSXP5Yqnfv3lzlFwAA1Foun7lZsWKFJk2apMzMTHXr1k0+Pj4Oz99xxx1VVhwAAICrXA43586d04kTJ/Sb3/zG3maz2ez3hKrI/aUAAACqisvh5vHHH1fPnj21bt06FhQDAIBax+Vwc/r0aX300Ufq2LFjddQDAABQKS4vKB4wYIC++uqr6qgFAACg0lw+czNs2DBNnz5dBw8e1O23315iQfHw4cOrrDgAAABXuRxuiu/YPW/evBLPsaAYAAC4m8vh5vp7SQEAANQmLq+5AQAAqM0qFG527typYcOGqWPHjurUqZOGDx+uxMTEqq4NAADAZS6Hm7feekuDBg1SvXr1NHXqVE2ZMkV169bVwIED9c4771RHjQAAAE5zec3Nyy+/rNdff13Tp0+3t02bNk0LFy7UH/7wB40ePbpKCwQAAHCFy2duTp48qWHDhpVoHz58uNLS0qqkKAAAgIpyOdyEhoZq27ZtJdq3bdum0NDQKikKAACgolz+WOp3v/udpk6dqpSUFPXt21c2m027d+/WmjVrtHjx4uqoEQAAwGkuh5vJkyerRYsWWrBggd59911JUteuXbVhwwaNGDGiygsEAABwhcvhRpIeeughPfTQQ1VdCwAAQKVV6Do3ly5d0ooVK/Tcc8/p4sWLkqQvv/xSmZmZVVocAACAq1w+c3PgwAENGjRIQUFBOnXqlCZOnKjGjRvrgw8+0OnTp7V27drqqBMAAMApLp+5iY2N1fjx43Xs2DH5+/vb24cMGaJdu3ZVaXEAAACucjncfPHFF3ryySdLtLdq1UrZ2dlVUhQAAEBFuRxu/P39lZubW6L96NGjatasWZUUBQAAUFEuh5sRI0Zo3rx5unbtmiTJZrMpPT1dM2fO1MMPP1zlBQIAALjC5XAzf/58nTt3Ts2bN9ePP/6oe++9Vx07dlRAQIBefvnl6qgRAADAaS5/WyowMFC7d+/Wp59+qi+//FKFhYXq1auXBg0aVB31AQAAuKRCF/GTpAEDBmjAgAFVWQsAAECluRxu3njjjVLbbTab/P391bFjR91zzz3y8vKqdHEAAACucjncLFq0SOfOndOVK1fUqFEjGWN06dIl1atXTw0aNNDZs2fVvn17bd++nbuEAwCAGufyguJXXnlFvXv31rFjx3ThwgVdvHhR33zzje666y4tXrxY6enpatGihaZPn14d9QIAAJTL5TM3zz//vDZu3KgOHTrY2zp27Kj58+fr4Ycf1smTJ/X666/ztXAAAOAWLoebrKws5efnl2jPz8+3X6G4ZcuWunz5cuWrA2pQQYGUmChlZUkhIVJUlMTSMaDm8DuIquLyx1L9+/fXk08+qeTkZHtbcnKyJk+ebP/21MGDB9WuXbuqqxKoZps2SW3bSv37S6NHF/23bduidgDVj99BVCWXw83KlSvVuHFjhYeHy8/PT35+foqIiFDjxo21cuVKSVKDBg20YMGCKi8WqA6bNkmPPCJ9+61je2ZmUTv/uALVi99BVDWbMcY429kYo/T0dDVr1kwZGRk6evSojDHq0qWLwsLCqrPOKpObm6ugoCDl5OQoMDDQ3eXAzQoKiv7v8Pp/VIvZbFLr1lJaGqfHPUGfV7YpO/cntQj0197nBrq7HDiB30E4y5X3b5fW3Bhj1KlTJx0+fFhhYWEeE2iAsiQmlv2PqiQZI2VkFPW7774aKwu4afA7iOrgUripU6eOOnXqpAsXLqhTp07VVRNQY7KyHLcDep9UHd98FV711uUv2pfZD0DV4HcQ1cHlNTevv/66nnnmGR06dKg66gFqVEiI43Zg7zQ1vPuYAnunldsPQNXgdxDVweWvgv/617/WlStX1L17d/n6+qpu3boOz1+8eLHKigOqW1RU0ef5mZlFp7+vV/x5f1RUzdcG3Az4HUR1cDncxMfHV0MZgHt4eUmLFxd9I8Nmc3yueDs+noWMQHXhdxDVweVwM27cuOqoA3CbkSOl99+Xpk2Tfv4/jq1bF/2jOnKkuyoDbg78DqKqubzmRpJOnDih559/Xr/61a909uxZSdK//vUvHT58uEqLA2rKyJHSqVNSs6ZF282aFn31lH9UgZrB7yCqksvhZufOnbr99tu1b98+bdq0Sd9//70k6cCBA5ozZ06VFwjUFC8vyc+v6M9+fpwGB2oav4OoKi6Hm5kzZ+qll15SQkKCfH197e39+/fX559/XqXFAQAAuMrlcHPw4EE99NBDJdqbNWumCxcuVElRAAAAFeVyuGnYsKGySrmaUnJyslq1alUlRQEAAFSUy+Fm9OjRevbZZ5WdnS2bzabCwkJ99tln+v3vf6+xY8dWR40AAABOczncvPzyy7rlllvUqlUrff/997r11lt1zz33qG/fvnr++eero0YAAACnuXydGx8fH7399tuaN2+ekpOTVVhYqJ49e3KvKQAAUCu4HG6KdejQQe3bF93UzHb9ZSUBAADcpEIX8Vu5cqW6desmf39/+fv7q1u3blqxYkVV1wYAAOAyl8PNCy+8oGnTpmnYsGF677339N5772nYsGGaPn16hdbcLFmyRO3atZO/v7/Cw8OVmJhYZt9NmzZp8ODBatasmQIDAxUZGaktW7a4/JoAAMC6XA43S5cu1Ztvvqm4uDgNHz5cw4cPV1xcnJYvX65ly5a5dKwNGzYoJiZGs2fPVnJysqKiojRkyBClp6eX2n/Xrl0aPHiwNm/erKSkJPXv31/Dhg1TcnKyq8MAAAAW5fKam4KCAkVERJRoDw8PV35+vkvHWrhwoSZMmKCJEydKKrrj+JYtW7R06VLFxcWV6H/9HclfeeUV/e1vf9Pf//539ezZs9TXyMvLU15enn07NzfXpRoBAIBncfnMza9//WstXbq0RPvy5cv12GOPOX2cq1evKikpSdHR0Q7t0dHR2rNnj1PHKCws1OXLl9W4ceMy+8TFxSkoKMj+CA0NdbpGAADgeSr0bamVK1dq69at6tOnjyRp7969ysjI0NixYxUbG2vvt3DhwjKPcf78eRUUFCg4ONihPTg4WNnZ2U7VsWDBAv3www/6xS9+UWafWbNmOdSUm5tLwAEAwMJcDjeHDh1Sr169JEknTpyQVHRfqWbNmunQoUP2fs5+Pfz6fsYYp/Zdt26dXnzxRf3tb39T8+bNy+zn5+cnv+LbzAIAAMtzOdxs3769Sl64adOm8vLyKnGW5uzZsyXO5lxvw4YNmjBhgt577z0NGjSoSuoBAADW4PKamzVr1ujHH3+s9Av7+voqPDxcCQkJDu0JCQnq27dvmfutW7dO48eP1zvvvKOhQ4dWug4AAGAtLoebWbNmKTg4WBMmTHB64W9ZYmNjtWLFCq1atUqpqamaPn260tPTNWnSJPtr/fxmnOvWrdPYsWO1YMEC9enTR9nZ2crOzlZOTk6l6gAAANbhcrj59ttv9dZbb+m7775T//791aVLF7322mtOLwL+uVGjRik+Pl7z5s1Tjx49tGvXLm3evFlt2rSRJGVlZTlc8+Yvf/mL8vPz9fTTTyskJMT+mDZtmsuvDQAArMnlNTdeXl72i/edPXtWb731ltasWaMXXnhB//Vf/6UJEyZo2LBhqlPHudz01FNP6amnnir1uTVr1jhs79ixw9VyAQDATaZC95Yq1rx5c/Xr10+RkZGqU6eODh48qPHjx6tDhw4EEQAA4BYVCjf/+c9/NH/+fN1222267777lJubq3/84x9KS0vTmTNnNHLkSI0bN66qawUAALghlz+WGjZsmLZs2aLOnTvriSee0NixYx2uEFy3bl397ne/06JFi6q0UAAAAGe4HG6aN2+unTt3KjIyssw+ISEhSktLq1RhAAAAFeFyuFm5cuUN+9hsNvs3ngAAAGpShe4ttW3bNm3btk1nz55VYWGhw3OrVq2qksIAAAAqwuVwM3fuXM2bN08REREKCQlx+h5SAAAANcHlcLNs2TKtWbNGY8aMqY56AAAAKsXlr4JfvXq13Hs/AQAAuJPL4WbixIl65513qqMWAACASnPqY6nY2Fj7nwsLC7V8+XJ98sknuuOOO+Tj4+PQd+HChVVbIQAAgAucCjfJyckO2z169JAkHTp0yKGdxcUAAMDdnAo327dvr+46AAAAqkSlbpwJAABQ2xBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApbg93CxZskTt2rWTv7+/wsPDlZiYWGbfrKwsjR49WmFhYapTp45iYmJqrlAAAOAR3BpuNmzYoJiYGM2ePVvJycmKiorSkCFDlJ6eXmr/vLw8NWvWTLNnz1b37t1ruFoAAOAJ3BpuFi5cqAkTJmjixInq2rWr4uPjFRoaqqVLl5bav23btlq8eLHGjh2roKCgGq4WAAB4AreFm6tXryopKUnR0dEO7dHR0dqzZ0+VvU5eXp5yc3MdHgAAwLrcFm7Onz+vgoICBQcHO7QHBwcrOzu7yl4nLi5OQUFB9kdoaGiVHRsAANQ+bl9QbLPZHLaNMSXaKmPWrFnKycmxPzIyMqrs2AAAoPbxdtcLN23aVF5eXiXO0pw9e7bE2ZzK8PPzk5+fX5UdDwAA1G5uO3Pj6+ur8PBwJSQkOLQnJCSob9++bqoKAAB4OreduZGk2NhYjRkzRhEREYqMjNTy5cuVnp6uSZMmSSr6SCkzM1Nr166175OSkiJJ+v7773Xu3DmlpKTI19dXt956qzuGAAAAahm3hptRo0bpwoULmjdvnrKystStWzdt3rxZbdq0kVR00b7rr3nTs2dP+5+TkpL0zjvvqE2bNjp16lRNlg4AAGopt4YbSXrqqaf01FNPlfrcmjVrSrQZY6q5IgAA4Mnc/m0pAACAqkS4AQAAlkK4AQAAlkK4AQAAlkK4AQAAlkK4AQAAlkK4AQAAlkK4AQAAlkK4AQAAlkK4AQAAlkK4AQAAlkK4AQAAlkK4AQAAlkK4AQAAlkK4AQAAlkK4AQAAlkK4AQAAlkK4AQAAlkK4AQAAlkK4AQAAlkK4AQAAlkK4AQAAlkK4AQAAlkK4AQAAlkK4AQAAlkK4AQAAlkK4AQAAlkK4AQAAlkK4AQAAlkK4AQAAlkK4AQAAlkK4AQAAlkK4AQAAlkK4AQAAlkK4AQAAlkK4AQAAlkK4AQAAlkK4AQAAlkK4AQAAlkK4AQAAlkK4AQAAlkK4AQAAlkK4AQAAlkK4AQAAluLt7gKsoqBASkyUsrKkkBApKkry8nJ3VQAA1Jza8l5IuKkCmzZJ06ZJ3377f22tW0uLF0sjR7qvLgAAakptei/kY6lK2rRJeuQRx8mUpMzMovZNm9xTFwAANaW2vRcSbiqhoKAopRpT8rnitpiYon4AAFhRbXwvJNxUQmJiyZT6c8ZIGRlF/QAAsKLa+F7ImptKyMpy3A7ofVJ1fPNVeNVbl79oX2Y/AACsoja+FxJuKiEkxHE7sHeavAN+Uv5lf4cJvb4fAABWURvfC/lYqhKioopWgttspT9vs0mhoUX9AACwotr4Xki4qQQvr6KvuEklJ7V4Oz6e690AAKyrNr4XEm4qaeRI6f33pVatHNtbty5q5zo3AACrq23vhW4PN0uWLFG7du3k7++v8PBwJd5gOfXOnTsVHh4uf39/tW/fXsuWLauhSss2cqR06pTUrGnRdrOmUloawQYAcPOoTe+Fbg03GzZsUExMjGbPnq3k5GRFRUVpyJAhSk9PL7V/WlqaHnjgAUVFRSk5OVnPPfecpk6dqo0bN9Zw5SV5eUl+fkV/9vPjoygAwM2ntrwXujXcLFy4UBMmTNDEiRPVtWtXxcfHKzQ0VEuXLi21/7Jly3TLLbcoPj5eXbt21cSJE/X4449r/vz5NVw5AACordwWbq5evaqkpCRFR0c7tEdHR2vPnj2l7vP555+X6H///fdr//79unbtWqn75OXlKTc31+EBAACsy23XuTl//rwKCgoUHBzs0B4cHKzs7OxS98nOzi61f35+vs6fP6+QUr5EHxcXp7lz51Zd4eWYGNVOl3/KV4A/lw/yVMyhZ2P+PB9z6Plqwxy6/afHdt33xowxJdpu1L+09mKzZs1SbGysfTs3N1ehoaEVLbdcE6Pa37gTajXm0LMxf56POfR8tWEO3RZumjZtKi8vrxJnac6ePVvi7EyxFi1alNrf29tbTZo0KXUfPz8/+RWvbgIAAJbntjU3vr6+Cg8PV0JCgkN7QkKC+vbtW+o+kZGRJfpv3bpVERER8vHxqbZaAQCA53Drt6ViY2O1YsUKrVq1SqmpqZo+fbrS09M1adIkSUUfKY0dO9bef9KkSTp9+rRiY2OVmpqqVatWaeXKlfr973/vriEAAIBaxq1rbkaNGqULFy5o3rx5ysrKUrdu3bR582a1adNGkpSVleVwzZt27dpp8+bNmj59uv785z+rZcuWeuONN/Twww+7awgAAKCWsZniFbk3idzcXAUFBSknJ0eBgYHuLgcAADjBlfdvt99+AQAAoCoRbgAAgKUQbgAAgKUQbgAAgKUQbgAAgKUQbgAAgKUQbgAAgKUQbgAAgKUQbgAAgKW49fYL7lB8Qebc3Fw3VwIAAJxV/L7tzI0Vbrpwc/nyZUlSaGiomysBAACuunz5soKCgsrtc9PdW6qwsFBnzpxRQECAbDZblR47NzdXoaGhysjIsOR9q6w+Psn6Y2R8ns/qY2R8nq+6xmiM0eXLl9WyZUvVqVP+qpqb7sxNnTp11Lp162p9jcDAQMv+0ErWH59k/TEyPs9n9TEyPs9XHWO80RmbYiwoBgAAlkK4AQAAlkK4qUJ+fn6aM2eO/Pz83F1KtbD6+CTrj5HxeT6rj5Hxeb7aMMabbkExAACwNs7cAAAASyHcAAAASyHcAAAASyHcAAAASyHcuGjJkiVq166d/P39FR4ersTExHL779y5U+Hh4fL391f79u21bNmyGqq0YlwZ344dO2Sz2Uo8vv766xqs2Hm7du3SsGHD1LJlS9lsNn344Yc33MfT5s/VMXrSHMbFxal3794KCAhQ8+bN9eCDD+ro0aM33M+T5rAiY/SkOVy6dKnuuOMO+8XdIiMj9fHHH5e7jyfNn6vj86S5K01cXJxsNptiYmLK7eeOOSTcuGDDhg2KiYnR7NmzlZycrKioKA0ZMkTp6eml9k9LS9MDDzygqKgoJScn67nnntPUqVO1cePGGq7cOa6Or9jRo0eVlZVlf3Tq1KmGKnbNDz/8oO7du+tPf/qTU/09bf4k18dYzBPmcOfOnXr66ae1d+9eJSQkKD8/X9HR0frhhx/K3MfT5rAiYyzmCXPYunVrvfrqq9q/f7/279+vAQMGaMSIETp8+HCp/T1t/lwdXzFPmLvrffHFF1q+fLnuuOOOcvu5bQ4NnHbnnXeaSZMmObR16dLFzJw5s9T+M2bMMF26dHFoe/LJJ02fPn2qrcbKcHV827dvN5LMd999VwPVVS1J5oMPPii3j6fN3/WcGaMnz+HZs2eNJLNz584y+3j6HDozRk+eQ2OMadSokVmxYkWpz3n6/BlT/vg8de4uX75sOnXqZBISEsy9995rpk2bVmZfd80hZ26cdPXqVSUlJSk6OtqhPTo6Wnv27Cl1n88//7xE//vvv1/79+/XtWvXqq3WiqjI+Ir17NlTISEhGjhwoLZv316dZdYoT5q/yvLEOczJyZEkNW7cuMw+nj6HzoyxmKfNYUFBgdavX68ffvhBkZGRpfbx5PlzZnzFPG3unn76aQ0dOlSDBg26YV93zSHhxknnz59XQUGBgoODHdqDg4OVnZ1d6j7Z2dml9s/Pz9f58+errdaKqMj4QkJCtHz5cm3cuFGbNm1SWFiYBg4cqF27dtVEydXOk+avojx1Do0xio2N1d13361u3bqV2c+T59DZMXraHB48eFANGjSQn5+fJk2apA8++EC33nprqX09cf5cGZ+nzZ0krV+/Xl9++aXi4uKc6u+uObzp7gpeWTabzWHbGFOi7Ub9S2uvLVwZX1hYmMLCwuzbkZGRysjI0Pz583XPPfdUa501xdPmz1WeOodTpkzRgQMHtHv37hv29dQ5dHaMnjaHYWFhSklJ0aVLl7Rx40aNGzdOO3fuLDMAeNr8uTI+T5u7jIwMTZs2TVu3bpW/v7/T+7ljDjlz46SmTZvKy8urxFmMs2fPlkilxVq0aFFqf29vbzVp0qTaaq2IioyvNH369NGxY8equjy38KT5q0q1fQ5/+9vf6qOPPtL27dvVunXrcvt66hy6MsbS1OY59PX1VceOHRUREaG4uDh1795dixcvLrWvJ86fK+MrTW2eu6SkJJ09e1bh4eHy9vaWt7e3du7cqTfeeEPe3t4qKCgosY+75pBw4yRfX1+Fh4crISHBoT0hIUF9+/YtdZ/IyMgS/bdu3aqIiAj5+PhUW60VUZHxlSY5OVkhISFVXZ5beNL8VaXaOofGGE2ZMkWbNm3Sp59+qnbt2t1wH0+bw4qMsTS1dQ5LY4xRXl5eqc952vyVprzxlaY2z93AgQN18OBBpaSk2B8RERF67LHHlJKSIi8vrxL7uG0Oq3W5ssWsX7/e+Pj4mJUrV5ojR46YmJgYU79+fXPq1CljjDEzZ840Y8aMsfc/efKkqVevnpk+fbo5cuSIWblypfHx8THvv/++u4ZQLlfHt2jRIvPBBx+Yb775xhw6dMjMnDnTSDIbN2501xDKdfnyZZOcnGySk5ONJLNw4UKTnJxsTp8+bYzx/PkzxvUxetIcTp482QQFBZkdO3aYrKws++PKlSv2Pp4+hxUZoyfN4axZs8yuXbtMWlqaOXDggHnuuedMnTp1zNatW40xnj9/ro7Pk+auLNd/W6q2zCHhxkV//vOfTZs2bYyvr6/p1auXw1c0x40bZ+69916H/jt27DA9e/Y0vr6+pm3btmbp0qU1XLFrXBnfa6+9Zjp06GD8/f1No0aNzN13323++c9/uqFq5xR/7fL6x7hx44wx1pg/V8foSXNY2rgkmdWrV9v7ePocVmSMnjSHjz/+uP3fl2bNmpmBAwfa3/iN8fz5c3V8njR3Zbk+3NSWObQZ8/9X9gAAAFgAa24AAIClEG4AAIClEG4AAIClEG4AAIClEG4AAIClEG4AAIClEG4AAIClEG4AAIClEG4AN/r666/Vp08f+fv7q0ePHu4up0otX75coaGhqlOnjuLj491dTqW1bdvWYRw2m00ffvhhpY45fvx4Pfjgg5U6Rm13M4wRtQ/hBnDCuXPn5OPjoytXrig/P1/169dXenp6pY87Z84c1a9fX0ePHtW2bduqoNLaITc3V1OmTNGzzz6rzMxM/fd//3eJPqdOnZLNZlNKSkrNF4gqx3yiNvF2dwGAJ/j888/Vo0cP1atXT/v27VPjxo11yy23VPq4J06c0NChQ9WmTZsqqLL2SE9P17Vr1zR06NBae4fjm5ExRgUFBfL2dvyn/+rVq/L19XVTVUDV48wN4IQ9e/aoX79+kqTdu3fb/1yewsJCzZs3T61bt5afn5969Oihf/3rX/bnbTabkpKSNG/ePNlsNr344otlHue1115Tx44d5efnp1tuuUUvv/yyJGnHjh2y2Wy6dOmSvX9KSopsNptOnTolSVqzZo0aNmyoDz/8UJ07d5a/v78GDx6sjIwM+z5fffWV+vfvr4CAAAUGBio8PFz79+8vc2zp6ekaMWKEGjRooMDAQP3iF7/Qf/7zH/vr3X777ZKk9u3bO9Tyc+3atZMk9ezZUzabTffdd5/9udWrV6tr167y9/dXly5dtGTJEod9n332WXXu3Fn16tVT+/bt9cILL+jatWtl1itJ3377rX75y1+qcePGql+/viIiIrRv3z5JRSFzxIgRCg4OVoMGDdS7d2998skn5R7vepmZmRo1apQaNWqkJk2aaMSIEQ7jLigoUGxsrBo2bKgmTZpoxowZcubWfp999pnuvfde1atXT40aNdL999+v7777TpKUl5enqVOnqnnz5vL399fdd9+tL774wr5v8c/Hli1bFBERIT8/PyUmJuq+++7TlClTFBsbq6ZNm2rw4MGSpCNHjuiBBx5QgwYNFBwcrDFjxuj8+fP245X3s1jefBZbu3atmjRpory8PIf2hx9+WGPHjnXuLxpwRrXfmhPwUKdPnzZBQUEmKCjI+Pj4GH9/fxMUFGR8fX2Nn5+fCQoKMpMnTy5z/4ULF5rAwECzbt068/XXX5sZM2YYHx8f88033xhjjMnKyjK33Xab+d3vfmeysrLM5cuXSz3OjBkzTKNGjcyaNWvM8ePHTWJionnzzTeNMf93F/DvvvvO3j85OdlIMmlpacYYY1avXm18fHxMRESE2bNnj9m/f7+58847Td++fe373HbbbebXv/61SU1NNd9884159913TUpKSqn1FBYWmp49e5q7777b7N+/3+zdu9f06tXLfifgK1eumE8++cRIMv/+979NVlaWyc/PL3Gcf//730aS+eSTT0xWVpa5cOGCMcaY5cuXm5CQELNx40Zz8uRJs3HjRtO4cWOzZs0a+75/+MMfzGeffWbS0tLMRx99ZIKDg81rr71W5lxcvnzZtG/f3kRFRZnExERz7Ngxs2HDBrNnzx5jjDEpKSlm2bJl5sCBA+abb74xs2fPNv7+/ub06dP2Y7Rp08YsWrTIvi3JfPDBB8YYY3744QfTqVMn8/jjj5sDBw6YI0eOmNGjR5uwsDCTl5dnjCm6A3RQUJB5//33zZEjR8yECRNMQECAGTFiRJl1JycnGz8/PzN58mSTkpJiDh06ZP74xz+ac+fOGWOMmTp1qmnZsqXZvHmzOXz4sBk3bpxp1KiR/e+y+OfjjjvuMFu3bjXHjx8358+fN/fee69p0KCBeeaZZ8zXX39tUlNTzZkzZ0zTpk3NrFmzTGpqqvnyyy/N4MGDTf/+/e31lPezWNZ8jhs3zj7GK1eumKCgIPPuu+/aj3nu3Dnj6+trPv300zL/HgBXEW6AMly7ds2kpaWZr776yvj4+JiUlBRz/Phx06BBA7Nz506TlpZmf5MpTcuWLc3LL7/s0Na7d2/z1FNP2be7d+9u5syZU+YxcnNzjZ+fn/0N5HrOhhtJZu/evfY+qampRpLZt2+fMcaYgIAAh/BQnq1btxovLy+Tnp5ubzt8+LA9zJRWQ2nS0tKMJJOcnOzQHhoaat555x2Htj/84Q8mMjKyzGO9/vrrJjw8vMzn//KXv5iAgAD7G64zbr31VvPHP/7Rvl1euFm5cqUJCwszhYWF9ufz8vJM3bp1zZYtW4wxxoSEhJhXX33V/vy1a9dM69atyw03v/rVr0y/fv1Kfe777783Pj4+5u2337a3Xb161bRs2dK8/vrrxpj/+/n48MMPHfa99957TY8ePRzaXnjhBRMdHe3QlpGRYSSZo0eP3vBnsaz5/Hm4McaYyZMnmyFDhti34+PjTfv27R3+7oDKYs0NUAZvb2+1bdtW7777rnr37q3u3bvrs88+U3BwsO65555y983NzdWZM2dKfHzVr18/ffXVV07XkJqaqry8PA0cOLBCYyjm7e2tiIgI+3aXLl3UsGFDpaam6s4771RsbKwmTpyov/71rxo0aJAeffRRdejQocyaQkNDFRoaam+79dZb7cfr3bt3hes8d+6cMjIyNGHCBD3xxBP29vz8fAUFBdm333//fcXHx+v48eP6/vvvlZ+fr8DAwDKPm5KSop49e6px48alPv/DDz9o7ty5+sc//qEzZ84oPz9fP/74o9OLxpOSknT8+HEFBAQ4tP/00086ceKEcnJylJWVpcjISPtzxXNiyvloKiUlRY8++mipz504cULXrl1z+Bnz8fHRnXfeqdTUVIe+P5/7stqSkpK0fft2NWjQoNTXunTpUpX8LD7xxBPq3bu3MjMz1apVK61evVrjx4+XzWar1HGBnyPcAGW47bbbdPr0aV27dk2FhYVq0KCB8vPzlZ+frwYNGqhNmzY6fPhwuce4/h9sY4xL/4jXrVu33Ofr1KljP26xstaelPa6xW0vvviiRo8erX/+85/6+OOPNWfOHK1fv14PPfRQiX3KGoOrYytNYWGhJOnNN9/UXXfd5fCcl5eXJGnv3r365S9/qblz5+r+++9XUFCQ1q9frwULFpR53Bv9PT7zzDPasmWL5s+fr44dO6pu3bp65JFHdPXqVafrDg8P19tvv13iuWbNmjl1jNKUV3fxnDvzM1a/fv0S+1/fVlhYqGHDhum1114r0TckJEQnT550uu7y9OzZU927d9fatWt1//336+DBg/r73/9eJccGirGgGCjD5s2blZKSohYtWuitt95SSkqKunXrpvj4eKWkpGjz5s1l7hsYGKiWLVtq9+7dDu179uxR165dna6hU6dOqlu3bplfEy9+48zKyrK3lfZV3Pz8fIcFwkePHtWlS5fUpUsXe1vnzp01ffp0bd26VSNHjtTq1atLfc1bb71V6enpDguSjxw5opycHJfGVvztnIKCAntbcHCwWrVqpZMnT6pjx44Oj+IFq5999pnatGmj2bNnKyIiQp06ddLp06fLfa077rhDKSkpunjxYqnPJyYmavz48XrooYd0++23q0WLFqUugi5Lr169dOzYMTVv3rxE3UFBQQoKClJISIj27t1r3yc/P19JSUk3rLusue/YsaN8fX0dfsauXbum/fv3uzQPPx/D4cOH1bZt2xJjqF+//g1/Fkubz7JMnDhRq1ev1qpVqzRo0CCHs4BAlXDjR2JArZeVlWX8/PzMjz/+aF9DkZmZ6dS+ixYtMoGBgWb9+vXm66+/Ns8++6zDgmJjbrzmxhhjXnzxRdOoUSPzv//7v+b48ePm888/NytWrDDGFK2xCA0NNY8++qg5evSo+cc//mHCwsJKXVB85513mr1795qkpCQTGRlp+vTpY4wpWuT59NNPm+3bt5tTp06Z3bt3mw4dOpgZM2aUWk/xguKoqCiTlJRk9u3bZ8LDw+0Lio1xbs3NtWvXTN26dc1LL71ksrOzzaVLl4wxxrz55pumbt26Jj4+3hw9etQcOHDArFq1yixYsMAYY8yHH35ovL29zbp168zx48fN4sWLTePGjU1QUFCZr5WXl2c6d+5soqKizO7du82JEyfM+++/b19Q/OCDD5oePXqY5ORkk5KSYoYNG2YCAgLMtGnT7MdwZkHxfffdZ3bt2mVOnjxpduzYYaZOnWoyMjKMMca8+uqrplGjRmbTpk0mNTXVPPHEEzdcUHz06FHj6+trJk+ebL766iuTmppqlixZYl/rNW3aNNOyZUvz8ccfOywovnjxojGm9DVZxhStufn52IwxJjMz0zRr1sw88sgjZt++febEiRNmy5Yt5je/+Y19QXh5P4tlzef1a26MMSYnJ8fUq1fP+Pr6mvXr15c5fqCiCDdAOdatW2fuvvtuY4wxu3btMh07dnR634KCAjN37lzTqlUr4+PjY7p3724+/vhjhz7OhJuCggLz0ksvmTZt2hgfHx9zyy23mFdeecX+/O7du83tt99u/P39TVRUlHnvvfdKhJugoCCzceNG0759e+Pr62sGDBhgTp06ZYwpeuP/5S9/aUJDQ42vr69p2bKlmTJlivnxxx/LrOn06dNm+PDhpn79+iYgIMA8+uijJjs72/68M+HGmKIgExoaaurUqeMQjt5++23To0cP4+vraxo1amTuueces2nTJvvzzzzzjGnSpIlp0KCBGTVqlFm0aFG54cYYY06dOmUefvhhExgYaOrVq2ciIiLsC6rT0tJM//79Td26dU1oaKj505/+VCIAlBdujCkKwmPHjjVNmzY1fn5+pn379uaJJ54wOTk5xpiiN/9p06aZwMBA07BhQxMbG2vGjh1bbrgxxpgdO3aYvn37Gj8/P9OwYUNz//3328PKjz/+aH7729/aX7Nfv372Rd3GuBZujDHmm2++MQ899JBp2LChqVu3runSpYuJiYmxL/a90c9iafNZWrgxxpgxY8aYxo0bm59++qnc8QMVYTPGiQstAPBYa9asUUxMjMO1cAB3Gzx4sLp27ao33njD3aXAglhQDACoMRcvXtTWrVv16aef6k9/+pO7y4FFEW4AADWmV69e+u677/Taa68pLCzM3eXAovhYCgAAWApfBQcAAJZCuAEAAJZCuAEAAJZCuAEAAJZCuAEAAJZCuAEAAJZCuAEAAJZCuAEAAJby/wAh+PBYhDzsUwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111)\n",
    "ax.plot(x, pmf_tea, 'bo')\n",
    "ax.vlines(x, 0, pmf_tea, lw=2)\n",
    "ax.set_xlabel('# of cups of tea called correctly')\n",
    "ax.set_ylabel('hypergeom PMF')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Probability of getting 3 or more cups correct is 0.24, which is greater than 0.05, so we retain the null hypothesis that the lady has no special ability.\n"
     ]
    }
   ],
   "source": [
    "# We're interested in the probability that she gets three or greater correctly if she\n",
    "# selects at random, which is:\n",
    "p_correct = hd.pmf(3) + hd.pmf(4)\n",
    "print(f\"Probability of getting 3 or more cups correct is {p_correct:.2f}, which is greater than 0.05, so we retain the null hypothesis that the lady has no special ability.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 9.2 Chi2 Contingency Test (1 DOF)\n",
    "\n",
    "# A test of the effect of a new drug on the heart rate has yielded the following results.\n",
    "\n",
    "# Does the drug affect the heart rate?\n",
    "# What would be the result if the response of one of the not-treated persons would\n",
    "# have been different? Perform this test with and without the Yates-correction."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "|             | Heart rate |               |       |\n",
    "|-------------|------------|---------------|-------|\n",
    "|             | Increased  | Not-increased | Total |\n",
    "| Treated     | 36         | 14            | 40    |\n",
    "| Not treated | 30         | 25            | 55    |\n",
    "| Total       | 66         | 39            | 105   |"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "p = 0.10 which is greater than 0.05, so columns are not contingent on rows, i.e. no effect of drug on heart rate.\n"
     ]
    }
   ],
   "source": [
    "data = np.array([[36, 14], [30, 25]])\n",
    "V, p, dof, expected = stats.chi2_contingency(data)\n",
    "print(f\"p = {p:.2f} which is greater than 0.05, so columns are not contingent on rows, i.e. no effect of drug on heart rate.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "p2 = 0.04, and p2_yates = 0.07\n"
     ]
    }
   ],
   "source": [
    "# The Yates correction is applied to mitigate the tendency of the chi-square test to overestimate \n",
    "# the significance of the association when the sample size is small. \n",
    "# It is used to reduce the likelihood of type I errors (false positives) by slightly downplaying\n",
    "# the test's sensitivity.\n",
    "\n",
    "data2 = np.array([[36, 14], [29, 26]])\n",
    "_, p2, _, _ = stats.chi2_contingency(data2, correction=False)\n",
    "_, p2_yates, _, _ = stats.chi2_contingency(data2, correction=True)\n",
    "print(f\"p2 = {p2:.2f}, and p2_yates = {p2_yates:.2f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 9.3 One Way Chi2-Test\n",
    "\n",
    "# The city of Linz wants to know if people want to build a long beach along the\n",
    "# Danube. They interview local people, and decide to collect 20 responses from each\n",
    "# of the five age groups: (<15, 15–30, 30–45, 45–60, >60)\n",
    "# The questionnaire states: “A beach-side development will benefit Linz.” and the\n",
    "# possible answers are:"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "| 1              | 2     | 3        | 4                 |\n",
    "|----------------|-------|----------|-------------------|\n",
    "| Strongly agree | Agree | Disagree | Strongly disagree |"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The city council wants to find out if the age of people influenced feelings about\n",
    "# the development, particularly of those who felt negatively (i.e., “disagreed” or\n",
    "# “strongly disagreed”) about the planned development."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "| Age group | Frequency of negative responses |\n",
    "|-----------|---------------------------------|\n",
    "| <15       | 4                               |\n",
    "| 15-30     | 6                               |\n",
    "| 30-45     | 14                              |\n",
    "| 45-60     | 10                              |\n",
    "| >60       | 16                              |"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "p= 0.03 which is less than 0.05, so the differences are statistically significant.\n",
      "Resulting analysis has 4 dof.\n"
     ]
    }
   ],
   "source": [
    "# The categories seem to show large differences of opinion between the groups.\n",
    "# • Are these differences significant?\n",
    "# • How many degrees of freedom does the resulting analysis have? \n",
    "\n",
    "freqs = np.array([4, 6, 14, 10, 16])\n",
    "V, p = stats.chisquare(freqs)\n",
    "print(f\"p= {p:.2f} which is less than 0.05, so the differences are statistically significant.\")\n",
    "print(f\"Resulting analysis has {len(freqs) - 1} dof.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.03420269940871678"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 9.4 McNemar’s Test\n",
    "\n",
    "# In a lawsuit regarding a murder the defense uses a questionnaire to show that the\n",
    "# defendant is insane. As a result of the questionnaire, the accused claims “not guilty\n",
    "# by reason of insanity.”\n",
    "\n",
    "# In reply, the state attorney wants to show that the questionnaire does not work.\n",
    "# He hires an experienced neurologist, and presents him with 40 patients, 20 of whom\n",
    "# have completed the questionnaire with an “insane” result, and 20 with a “sane”\n",
    "# result. When examined by the neurologist, the result is mixed: 19 of the “sane”\n",
    "# people are found sane, but 6 of the 20 “insane” people are labeled as sane by the\n",
    "# expert."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "|        | Sane by expert | Insane by expert | Total |\n",
    "|--------|----------------|------------------|-------|\n",
    "| Sane   | 19             | 1                | 20    |\n",
    "| Insane | 6              | 14               | 20    |\n",
    "| Total  | 22             | 18               | 40    |"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "p = 0.12 and p_correct = 0.03\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "c:\\Users\\rich\\miniconda3\\envs\\pydata-book\\lib\\site-packages\\statsmodels\\sandbox\\stats\\runs.py:542: FutureWarning: Deprecated, use stats.TableSymmetry instead\n",
      "  warnings.warn(\"Deprecated, use stats.TableSymmetry instead\", FutureWarning)\n"
     ]
    }
   ],
   "source": [
    "# Is this result significantly different from the questionnaire? \n",
    "# Would the result be significantly different, if the expert had diagnosed all “sane”\n",
    "# people correctly? \n",
    "\n",
    "# The McNemar test examines the difference between the proportions that derive from the marginal sums\n",
    "# of the table, i.e. (a + b) / N and (a + c) / N.\n",
    "\n",
    "obs = [[19, 1], [6, 14]]\n",
    "obs_correct = [[20, 0], [6, 14]]\n",
    "\n",
    "_, p = mcnemar(obs)\n",
    "_, p_correct = mcnemar(obs_correct)\n",
    "print(f\"p = {p:.2f} and p_correct = {p_correct:.2f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "pydata-book",
   "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.10.10"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}