BDA 3.10.4.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.stats import beta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Problem 3.10.4\n",
    "\n",
    "Inference for a 2 × 2 table: an experiment was performed to estimate the effect of betablockers on mortality of cardiac patients. A group of patients were randomly assigned to treatment and control groups: out of 674 patients receiving the control, 39 died, and out of 680 receiving the treatment, 22 died. Assume that the outcomes are independent and binomially distributed, with probabilities of death of p0 and p1 under the control and treatment, respectively. We return to this example in Section 5.6. \n",
    "\n",
    "(a) Set up a noninformative prior distribution on ( p0 , p1 ) and obtain posterior simulations. \n",
    "\n",
    "(b) Summarize the posterior distribution for the odds ratio, ( p1 /(1−p1 )) /( p0 /(1−p0 )). \n",
    "\n",
    "(c) Discuss the sensitivity of your inference to your choice of noninformative prior density.\n",
    "\n",
    "Gelman, Andrew; Carlin, John B.; Stern, Hal S.; Dunson, David B.; Vehtari, Aki; Rubin, Donald B.. Bayesian Data Analysis, Third Edition (Chapman & Hall/CRC Texts in Statistical Science) (Page 80). CRC Press. Kindle Edition. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Here we use the beta(1,1) uniform prior"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEICAYAAABPgw/pAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAGo5JREFUeJzt3X+4XVV95/H3hySABTRAAg0J8aIyFPoHPyalTLGaBmmBKKAFrDotauZJ+0x/YNFCxM6MMOgE20ejHW0nDygoaoioFUmrZlLS1ilGAoQARhtIIwkJJEhSQRQIfOePta4cLufcc+6959x9ztqf1/Oc5579e+299v7utdY+625FBGZmNvj2qzoBZmbWHQ7oZmaFcEA3MyuEA7qZWSEc0M3MCuGAbmZWCAf0Ckm6XtLVVafDrGqS5kvaPsr0vrhWJF0h6dqq09FKLQK6pLdLWi/pSUk7Jf29pNdOcJ19cYKVIufN8Od5ST9tGH5HD7Z3o6QPdnu9HW77NZKK7gAi6Z2S7pX0lKRHJP21pOkVpCMk/SSfRw9L+qikKR0u+5KbTER8OCL+S29SO3HFB3RJlwLLgA8DRwJzgU8B5/V4u1N7uf7SRMTBwx/gIeBNDeM+P3J+H9/+Jem9wDXAnwGvAE4DXgmslrR/BUk6MZ9XrwfeCry7gjRMjogo9kM6mZ4ELmwx/QBSsN+RP8uAA/K0+cB24L3ALmAn8K48bTHwLPBMXv/X8/itwOXARuBpYCpwPLAW2AvcD5zbsP3rgaurPk799snH8Q0jxl0N3AR8EXgCeCepQHIF8CDwGLACODTPvx9wM/BIPvZrgePztP86Iv++msdvB94H3JfHLycVAr4J/Bj4FjC9IU2nA9/J698AvK5h2reBK4F/yen9BnBYnrYDiLyNJ4FfqfqYdzHvXp736aIR4w/O19G78/DL8vm/B/geKfhvb5j/ZOCufOxuynl7dZ42A7g1H/fHgX8G9muRngBe0zC8Evhkw/C7gE15O1uA38/jDwJ+CjzfkE9HAR8EbmxY/tx8Xb/oHMvTLgcezuv+AXBGz49/1SdAj0+us4B9wNQW06/KF+QRwMx88f3PPG1+XvYqYBpwDvAULwSM6xkRjEmBaANwdD5hpwEPkILO/sCCnLnHtVqHP6MG9GeAN5GC9ctIwff/AbOBA4HrgM/l+fcjBf1D8rT/DaxvWN+NwAdHbGN7PgeOAOYAPwLWAyfmdfwj8IE879F5+m/lbZ1Fuqkcnqd/G9gMHAv8AinoDAek1wBR9XHuUd61vOaAG4Av5u9L8zE5LB/L+8gBPV8rPwT+NF9DF5BuwMPH738Bf5OnTQN+HVCL9Pw8oAO/RCqY/WnD9IXAqwGRSvBPAafkafNpuMnkcR8kB3TgPwA/Ac7M6bgsX+/7A8cB24Cj8rxDwKt7ffxLb3I5HHgsIva1mP4O4KqI2BURu0klqt9tmP5snv5sRPwd6S59XJttfiIitkXET0lVzYOBpRHxTET8A6lk8bYJ7FOdfTsivh4Rz+fj+/vAFRHxcET8jHSxXSRpvzzP9RHxRMO0/yjpoDbb+Hg+H7aTgvLtEXFPXsffkkqOAL8H3BIR38zb+gZwDymgDbsuIjZHxFPAl4CTunMY+toMWl9zO/N0gIuAD0XE4xGxDfhEw3ynkQLksnzt3Qzc0TD9WWAW8Mo8/Z8jR80W7pL0E1JJfC2pyRWAiFgVEQ9G8o+kWtivd7ivbwVWRcTqiHgW+EtSQePXgOdILQAnSJoWEVsj4sEO1ztupQf0HwEzRmlvPYpUEhj2wzzu58uPODGfIgXo0Wwbsf5tEfH8iG3MbrMOa27biOG5wNcl7ZW0F7iXVCI7QtIUSR+RtEXSj0klJ3ghoLTyaMP3nzYZHs7/VwJvG9523v5pvPj8eaTheyfnTgkeo/U1NytPh3xtNExrvA6PAh4eEaQbp/8FKT+/lfN3SZs0nUI69m8FfpXUnAKApLMlfUfS4zkPz6H9OdKYzp+nK1/n24DZEfEA8B5SQWKXpBWSjmq6li4qPaDfDvwMOL/F9B2kC3PY3DyuE61KBI3jdwBHS2o8znNJ7Wo2diOP+XbgzIiY3vA5MCIeIZWgzyE1c72C1MwBqWrdbF1jtQ34zIhtHxQRfzGO/SjJ7aTnR29pHJlrRmcDa/KonaSmlmFzG77vBGZLUrPpudb13oh4FakJ7lJJZ4yWqFwCX5nT999zmg4AvkwqWR8ZEdOBv6Pzc+RF8SOn92jy9R0RX4iI1+Z5gvSguKeKDugR8e+kzPukpPMl/YKkafmu/BHSA7Y/lzRT0ow8740drv5R4FVt5llHamO7LG93PukEXDGe/bGX+Bvgw5LmAkg6QtK5edohpMDyI1Ib9odGLNtJ/o3mc8CbJZ2ZawMHSvqNDkthu4CQNJHt96V8zV0J/JWks/J5P0RqctpOOm6QHk6+X9KhkuYAf9ywmttJ7fB/ImmqpLcApw5PlPTG/NNPkR5WP5c/nVgKLJb0i6S27gOA3cA+SWcDv9kw76PA4ZJe0WJdK4GFks6QNI30A4qngX+RdJykBfmm8TNS7a7TNI5b0QEdICI+ClwK/Dkp47YBf0RqD72a9NBrI6m6flce14nrSO1jeyX9bYttP0N6Cn42qar5KeD3IuL7494ha/RR0q9H1kh6gvRA81fytM/wwq+X7s/TGl0LnChpj6Sbx7rhiNgKvBn4b6Tz6iHSBd32moqIJ0gP9tbl82feWLffzyLiI6QfAvwlKeCuI113Z0TE03m2K0nNFf9Garf+XMPyz5BK+O8k/QrmrcBXGjZxLPB/Sc+0bgc+FRFrO0zbvaSH23+W8+FPSIF5D/B24JaGeb9PKvRtyfl01Ih1/QD4z8Bfka7vN5F+bvsM6UaxNI9/hPSg/YpO0jgRGv1ZgpmZDYriS+hmZnXhgG5mVggHdDOzQjigm5kVYlL/wdGMGTNiaGhoMjdpTdx5552PRcTMbq3P+dofnK/l6jRvJzWgDw0NsX79+sncpDUh6Yft5+qc87U/OF/L1WneusnFzKwQDuhmZoVwQDczK4QDuplZIRzQzcwK4YBuZlYIB3Qzs0I4oJuZFcIB3cysEA7oLQwtWVV1EmzA+JyZXD7eL+WAbmZWCAd0M7NCOKCbFUbSdEk3S/q+pE2S/pOkwyStlrQ5/z206nR2w9CSVW56aeCAblaejwPfiIhfAk4ENgFLgDURcSywJg9bYRzQzQoi6eXA64DrACLimYjYC5wH3JBnuwE4v5oUWi85oHfAVTobIK8CdgOfkXS3pGslHQQcGRE7AfLfI5otLGmxpPWS1u/evXvyUj0KN6t0zgHdbIL6LNhMBU4B/joiTgZ+whiaVyJieUTMi4h5M2d27eVHNkkc0GusTg/PamQ7sD0i1uXhm0kB/lFJswDy310Vpc96yAG9QQ2rdn54VpiIeATYJum4POoM4HvALcDFedzFwNcqSJ712KS+U9T6R8PDs3dCengGPCPpPGB+nu0GYC1w+eSn0Cbgj4HPS9of2AK8i1R4WylpEfAQcGGF6bMecUCvr8aHZycCdwKXMOLhmaSWD8+AxQBz586dnBRbRyJiAzCvyaQzJjstNrnc5JLVrKkF/PDMrDgO6PXlh2dmhXFAryk/PDMrj9vQqWVzyzA/PDMriAN6jfnhmZVuuLC2denCilMyOdzkYmZWCAd0M7NCOKCbdUGz5zA1fjZjFXEbupkNFN8oW3MJ3cysEA7oZl3mEqRVxU0uZjYQfKNszyV0M7NCOKCbmRXCAd1snNwE0H9q+JKaF3FANzMrRMcBXdKU/BbxW/PwMZLW5XdP3pT/wZOZmVVkLCX0S0jvnBx2DfCx/O7JPcCibibMzMzGpqOALmkOsBC4Ng8LWEB6KQKkd0+e34sEmplZZzotoS8DLgOez8OHA3sjYl8e3g7M7nLaKlPnhypmNrjaBnRJbwR2RcSdjaObzBotll8sab2k9bt37x5nMs3MrJ1OSuinA+dK2gqsIDW1LAOmSxruaToH2NFsYb9M2MyqVpdad9uu/xHxfuD9AJLmA++LiHdI+hJwASnID+S7J+uSyWZWDxP5HfrlwKWSHiC1qV/XnSSZmdl4jOmfc0XEWmBt/r4FOLX7STIz6746vF/U/23RrAeqbM7Lz7ueAJ4D9kXEPEmHATcBQ8BW4KKI2FNVGq033PXfrEy/EREnRcS8PLwEWJM7Aq7Jw1YYB3SzejiP1AEQ3BGwWG5yqTFXzYsVwLckBfB/ImI5cGRE7ASIiJ2Sjmi2oKTFwGKAuXPnTlZ6J8y/WEtcQjdXzctzekScApwN/KGk13W6oPuNDDYHdBvJVfMBFxE78t9dwFdJv0Z7VNIsgPx3V3UptF5xQK+34ar5nbmqDSOq5kDLqrn/pUP/kXSQpEOGvwO/CdwH3ELqAAgD2hFwLOr6ogu3odfb6RGxI7enrpb0/U4XzO2yywHmzZvX9P/4WCWOBL6a/iEqU4EvRMQ3JN0BrJS0CHgIuLDCNFqPOKDXWGPVXNKLqub5wZmr5gMmd/g7scn4HwFnTH6KbDK5yaWmXDU3K49L6PXlqrlZYRzQa8pVc7PyuMnFzKwQDuhmZoVwQDczK4QDulkXjezMUsfOLVYdPxQdhS9GMxskLqGbmRXCAd3MrBAO6GZmhXBANzMrhAO6mVkhHNA7VNf/r2xmg8MB3WwCOrnJuyBgk8UB3cysEA7oZmaFcEA3MyuEA7qZWSEc0M3MCuGAbmZWCAd0M7NCOKCbmRXCAd3MrBAO6GYFkjRF0t2Sbs3Dx0haJ2mzpJsk7V91Gq37HNDNxmEAuvNfAmxqGL4G+FhEHAvsARZVkirrqbYBXdKBkr4r6R5J90u6Mo8f6Dv+eC/IAbiQreYkzQEWAtfmYQELgJvzLDcA51eTOuulTkroTwMLIuJE4CTgLEmn4Tt+EVw1L9Iy4DLg+Tx8OLA3Ivbl4e3A7GYLSlosab2k9bt37+59Sq2r2gb0SJ7Mg9PyJ/AdvxSumhdE0huBXRFxZ+PoJrNGs+UjYnlEzIuIeTNnzuxJGq13OmpDz6W4DcAuYDXwIL7jDzxXzYt0OnCupK3AClJ+LgOmS5qa55kD7KgmedZLHQX0iHguIk4inQinAsc3m63Fsr7j9y9XzQsTEe+PiDkRMQT8DvAPEfEO4DbggjzbxcDXKkqi9dDU9rO8ICL2SloLnEa+4+eL33f8AdNYNZc0f3h0k1lb3qiB5QDz5s1rOo/1lcuBFZKuBu4Grqs4PS35hwfj1zagS5oJPJuD+cuAN5DaWYfv+CvwHX8QDVfNzwEOBF5OQ9XcN+rBFxFrgbX5+xZS7doK1kmTyyzgNkkbgTuA1RFxK+mOf6mkB0hV9b6949tLuWo+fi5BWr9qW0KPiI3AyU3G+45fpoGpmpvZi42pDd3K5Kq5WRnc9d/MrBAO6GZmhXBANzMrhAO6mVkhHNDNzArhgG5mtTK0ZFWxfQmKDuilZpqZWTNFB3QzszpxQDczK0RtAno3283clGNm/ag2Ad3MrHQO6GZmhXBANzMrhAO6mVkhHNDNzArhgG5mVggHdDOzQjigm5kVonavoHOnILP+1C/X5tCSVWxdurDqZIyLS+hmZoVwQDczK4QDullBJB0o6buS7pF0v6Qr8/hjJK2TtFnSTZL2rzqt1n0O6GZleRpYEBEnAicBZ0k6DbgG+FhEHAvsARZVmEbrEQd0sw71y0O70UTyZB6clj8BLABuzuNvAM6vIHnWYw7oNeWqebkkTZG0AdgFrAYeBPZGxL48y3ZgdlXps95xQK8vV80LFRHPRcRJwBzgVOD4ZrM1W1bSYknrJa3fvXt3L5NpPeCAXlOumpcvIvYCa4HTgOmShvudzAF2tFhmeUTMi4h5M2fOnJyEWtc4oNfYRKrmLsmN3WS0wUuaKWl6/v4y4A3AJuA24II828XA13qeGJt0Dug1NpGquUtyfWsWcJukjcAdwOqIuBW4HLhU0gPA4cB1FabReqR2Xf/tpSJir6S1NFTNcym9ZdXc+lNEbARObjJ+C+mmbQVzCb2mXDU3K49L6PU1C7hB0hTSjX1lRNwq6XvACklXA3fjqrkVbhD6F3TKAb2mXDU3K4+bXMzMCtG2hC7paOCzwC8CzwPLI+Ljkg4DbgKGgK3ARRGxp3dJHZ+SqlNWPZ9P1s86KaHvA94bEceTfgXxh5JOAJYAa3KPwjV52MzMKtI2oEfEzoi4K39/gvRLiNnAeaSehOAehWZmlRtTG7qkIdKDtHXAkRGxE1LQB45osYx7FJqZTYKOA7qkg4EvA++JiB93upx7FJqZTY6OArqkaaRg/vmI+Eoe/aikWXn6LNL/AzEzs4q0DeiSROpcsikiPtow6RZST0Jwj0Izs8p10rHodOB3gXvzf+YDuAJYCqyUtAh4CLiwN0k0M+u+En+C2jagR8S3AbWYfEZ3k2NWpuHgMbRkFVuXLqw4NVYq9xQ1MyuEA7qZWSEc0M3MCuGAbtZGiQ/PrEwO6GZmhXBANzMrhAN6F7hKbmb9wAHdzKwQDuhmZoVwQDczK4QD+jgNLVnltvMacV7bIHBANzMrhAO6WUEkHS3pNkmbJN0v6ZI8/jBJqyVtzn8PrTqt1n0O6DXlC79Yfql7jTmg15cv/FEMapu5X+pebw7oNeULv3x+qXv9OKCbL/wC+aXu9eSAXnO+8Mvjl7rXlwN6jfnCL49f6t4dg9rPxAG9pnzhF2v4pe4LJG3In3NIL3U/U9Jm4Mw8bIVp+5JoK9bwhX+vpA153BWkC32lpEXAQ8CFFaXPxsEvda83B/Sa8oVvVh43uZiZFcIB3cysEA7oZmaFcEA3G8Ug/nTN6ssB3cysEA7oZmaFKCKgN+vV5aqymdVNEQHdzMwc0M3MiuGAbmZWCAd0M7NCOKCbteAH6zZoHNDNzArRNqBL+rSkXZLuaxjnN8ObWde4NtQdnZTQrwfOGjHOb4Y3M+szbQN6RPwT8PiI0X4zvBVpMkqKLo1ar4y3Db2jN8ND+W+H98VpZv2i5w9F/XZ4M7PJMd6A7jfDW7Fc67JBNd6A7jfDm5n1mU5+tvhF4HbgOEnb89vglwJnStoMnJmHzcysQlPbzRARb2sxyW+GNzPrI8X1FG32v9FL3q4NLp8v1m3FBXSzunPv7u4ZtIKaA3pN+aIv2vW4d3ctOaDX1/X4oi+Se3fXlwN6Tfmir52OeneX3rO7W/q1KcYB3Rr5XzrUnHt2D7aiAno/3jFL5Qt/4Lh3dw0UFdBtwnzRl8u9u2vAAd0a+aIvgHt3T55+axVo21PUypQv+vnADEnbgf9BushX5gDwEHBhdSm08XLv7voa6IA+tGQVW5curDoZA8kXvVl53ORilk1m9bnfquqTrV9/9tfKWNJb5b45oJuZFcIB3cysEA7oZmaFcEA3w23aVoaBD+j9eCEO2gMfMyvDwAd0MzNLBvp36GY2WAa95tos/cPj+qFPjEvoXTboJ6xNPp8z1i0O6GZmhXCTi5lVxrWT7hq4ErpPADOz5gYuoJuZWXNucjEz64KRrQcjhyfjVzADVUJ3c4uZWWsDFdDNzKw1B3Qzs0I4oPeQm4jMrJle/b8nB3Qzs0I4oFstDJeGmpWMXJOyUjigm1nP+V9KT07BYSAD+qCeHIOYZjMbHAMZ0M3M7KXcU9Rqo1nbedX/w7qxbX/r0oV9kSabPO1q7WM9FyZUQpd0lqQfSHpA0pKJrGtYs4dXg9rEAoPZzNKLfLX+4Lwt27gDuqQpwCeBs4ETgLdJOqFbCbNqOF/L5bwt30RK6KcCD0TEloh4BlgBnNedZFmFnK/lct4WbiJt6LOBbQ3D24FfHTmTpMXA4jz4tKT72q1Y1zT/3idmAI+NZYE+3J/jRpnWs3ztJ8P50JAfY87XbmuSprEaLV+hg7wd9HxtovJ8bTTWvG2Yv13eAhML6GoyLl4yImI5sBxA0vqImDeBbVaulH0YbXKTcc7XAdAmX6GDvHW+9qcO8haYWJPLduDohuE5wI4JrM/6g/O1XM7bwk0koN8BHCvpGEn7A78D3NKdZFmFnK/lct4WbtxNLhGxT9IfAd8EpgCfjoj72yy2fLzb6yNF74PzdaCNug/jyNvij8kA6Wg/FPGS5lEzMxtA7vpvZlYIB3Qzs0J0JaC3604s6ZWS1kjaKGmtpDkN0y6WtDl/Lu5GesZrgvvxnKQN+VPJgyZJn5a0q9Vvh5V8Iu/fRkmnNExrmg8l5O2g52tOR9fzts32BvpfBEg6WtJtkjZJul/SJVWnabwkTZF0t6Rb284cERP6kB6uPAi8CtgfuAc4YcQ8XwIuzt8XAJ/L3w8DtuS/h+bvh040TZO9H3n4ySrSPSJ9rwNOAe5rMf0c4O9Jv0c+DVg3Wj6UkLcl5Gsv8naix6zfP8As4JT8/RDgXwdtHxr25VLgC8Ct7ebtRgm9k+7EJwBr8vfbGqb/FrA6Ih6PiD3AauCsLqRpPCayH30hIv4JeHyUWc4DPhvJd4DpkmbROh9KyNuBz1foSd6OZuD/RUBE7IyIu/L3J4BNpJ6yAyXXFhcC13YyfzcCerPuxCMP3D3Ab+fvbwYOkXR4h8tOlonsB8CBktZL+o6k83ub1HFrtY9jHd+o3/O2DvkKE8vDTtc1kCQNAScD66pNybgsAy4Dnu9k5m4E9E66ir8PeL2ku4HXAw8D+zpcdrJMZD8A5kbqYvx2YJmkV/cspePXah/HOr5Rv+dtHfIVJpaHna5r4Eg6GPgy8J6I+HHV6RkLSW8EdkXEnZ0u042A3rY7cUTsiIi3RMTJwAfyuH/vZNlJNJH9ICJ25L9bgLWkEkG/abWPYx3/cwOQt3XIV5hAHo5hXQNF0jRSMP98RHyl6vSMw+nAuZK2kpq9Fki6cdQlutBgP5X0oOUYXniA8ssj5pkB7Je/fwi4Kl54YPNvpIc1h+bvh1X04GEi+3EocEDDPJup6AEMMETrB2cLefGDs++Olg8l5G0p+drtvJ3oMev3Tz4OnwWWVZ2WLu3PfDp4KNqtjZ1Deor8IPCBPO4q4Nz8/YJ8MfwrqXH/gIZl3w08kD/vqvigjWs/gF8D7s0n/r3AoorS/0VgJ/AsqZS1CPgD4A/ydJFecPBgTue8dvlQQt4Oer72Km/HeswG6QO8ltRMtBHYkD/nVJ2uCezPfDoI6O76b2ZWCPcUNTMrhAO6mVkhHNDNzArhgG5mVggHdDOzQjigm5kVwgHdzKwQ/x/6xVs9pHnlsQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "odds summary: 1.9014804902473739 0.2770084241734924 [0.79448631 1.51230717 1.83462062 2.19570031 3.17125921]\n"
     ]
    }
   ],
   "source": [
    "p0=beta(636,40)\n",
    "p1=beta(659,23)\n",
    "p0_sim=p0.rvs(1000)\n",
    "p1_sim=p1.rvs(1000)\n",
    "fig,ax=plt.subplots(1,3)\n",
    "p0hist=ax[0].hist(p0_sim,density=True,bins=50)\n",
    "ax[0].set_xlim(0.9,1)\n",
    "ax[1].set_xlim(0.9,1)\n",
    "c=ax[0].set_title('Control')\n",
    "p1hist=ax[1].hist(p1_sim,density=True,bins=50)\n",
    "t=ax[1].set_title('Treatment')\n",
    "odds_ratio=(p1_sim/(1-p1_sim)/(p0_sim/(1-p0_sim)))\n",
    "ax[2].set_xlim(0,4)\n",
    "ax[2].hist(odds_ratio,bins=50)\n",
    "o=ax[2].set_title('Odds Ratios')\n",
    "plt.show()\n",
    "print('odds summary:',np.mean(odds_ratio),np.var(odds_ratio),np.percentile(odds_ratio,[.025,25,50,75,97.5]))\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Here we use the beta(0,0) improper prior; the difference from the previous case is small"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEICAYAAABPgw/pAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAG6RJREFUeJzt3X+0XGV97/H3BxJAAQ0hgRsI8aBwKfQPfvSUcovVGIzll4BeQNFrUelKu/pLiy1E7A/lojfYLo32aruyQElFBYpSkbRoiqSWW4wECL+MNpBGEhJIgKSAKCHwvX88zyHD4Zwz+5yZOTPz7M9rrVkzs3/MfmY/s7/z7OfH3ooIzMys/+3W7QSYmVl7OKCbmRXCAd3MrBAO6GZmhXBANzMrhAO6mVkhHNC7SNJVki7rdjrMuk3SXEkbx5jfE8eKpEskXdHtdIymFgFd0nskrZL0jKTNkv5Z0htb/Mye+IGVIufN0ONFST9veP/eDmzvakkfb/fnVtz2YZKKHgAi6f2S7pP0rKRHJf2tpGldSEdI+ln+HT0i6TOSdq+47iv+ZCLiUxHx251JbeuKD+iSLgQWA58CDgTmAF8Ezuzwdqd08vNLExH7DD2Ah4G3N0z76vDlvX97l6SPAJcDfwq8FjgBeB2wXNIeXUjS0fl39WbgXcAHu5CGyRERxT5IP6ZngHNGmb8nKdhvyo/FwJ553lxgI/ARYAuwGfhAnrcAeB7YkT//23n6euBi4F7gOWAKcCSwAtgOPACc0bD9q4DLur2feu2R9+Nbh027DLgW+DrwNPB+UoHkEuAh4HHgGmC/vPxuwPXAo3nfrwCOzPN+b1j+3ZCnbwT+BLg/T19CKgR8B3gK+C4wrSFNJwI/yJ+/GnhTw7zbgE8A/57TezMwPc/bBETexjPAr3Z7n7cx716Tv9O5w6bvk4+jD+b3r8q//23Aj0jBf2PD8scCd+V9d23O28vyvBnATXm/Pwn8G7DbKOkJ4LCG99cBX2h4/wFgTd7OOuB38vS9gZ8DLzbk00HAx4GrG9Y/Ix/XL/uN5XkXA4/kz/4JcFLH93+3fwAd/nGdDOwEpowy/9J8QB4AzMwH3//O8+bmdS8FpgKnAs+yK2BcxbBgTApEq4FD8g92KvAgKejsAczLmXvEaJ/hx5gBfQfwdlKwfhUp+P4/4GBgL+BK4Ct5+d1IQX/fPO//AqsaPu9q4OPDtrEx/wYOAGYDTwCrgKPzZ/wr8LG87CF5/m/mbZ1M+lPZP8+/DVgLHA68mhR0hgLSYUB0ez93KO9GPeaApcDX8+tFeZ9Mz/vyfnJAz8fKT4E/zsfQ2aQ/4KH993+Av8vzpgK/AWiU9LwU0IFfIhXM/rhh/mnAGwCRSvDPAsfleXNp+JPJ0z5ODujAfwd+BszP6bgoH+97AEcAG4CD8rIDwBs6vf9Lr3LZH3g8InaOMv+9wKURsSUitpJKVO9rmP98nv98RPwT6V/6iCbb/HxEbIiIn5NONfcBFkXEjoj4HqlkcV4L36nObouIb0fEi3n//g5wSUQ8EhG/IB1s50raLS9zVUQ83TDvVyTt3WQbn8u/h42koHx7RNyTP+MfSSVHgN8CboyI7+Rt3QzcQwpoQ66MiLUR8SzwD8Ax7dkNPW0Gox9zm/N8gHOBT0bEkxGxAfh8w3InkALk4nzsXQ/c0TD/eWAW8Lo8/98iR81R3CXpZ6SS+ApSlSsAEbEsIh6K5F9JZ2G/UfG7vgtYFhHLI+J54K9JBY1fB14g1QAcJWlqRKyPiIcqfu6ElR7QnwBmjFHfehCpJDDkp3naS+sP+2E+SwrQY9kw7PM3RMSLw7ZxcJPPsJFtGPZ+DvBtSdslbQfuI5XIDpC0u6RPS1on6SlSyQl2BZTRPNbw+ucjvB/K/9cB5w1tO2//BF7++3m04XWV304JHmf0Y25Wng/52GiY13gcHgQ8MixIN87/K1J+fjfn78ImaTqOtO/fBfwaqToFAEmnSPqBpCdzHp5K899IYzpfSlc+zjcAB0fEg8CHSQWJLZKukXTQiJ/SRqUH9NuBXwBnjTJ/E+nAHDInT6titBJB4/RNwCGSGvfzHFK9mo3f8H2+EZgfEdMaHntFxKOkEvSppGqu15KqOSCdWo/0WeO1AfjysG3vHRF/NYHvUZLbSe1H72ycmM+MTgFuyZM2k6pahsxpeL0ZOFiSRpqfz7o+EhGvJ1XBXSjppLESlUvg1+X0/UVO057AN0gl6wMjYhrwT1T/jbwsfuT0HkI+viPiaxHxxrxMkBqKO6rogB4R/0XKvC9IOkvSqyVNzf/KnyY1sP2ZpJmSZuRlr6748Y8Br2+yzEpSHdtFebtzST/AaybyfewV/g74lKQ5AJIOkHRGnrcvKbA8QarD/uSwdavk31i+ArxD0vx8NrCXpLdULIVtAUJSK9vvSfmY+wTwN5JOzr/7AVKV00bSfoPUOPlRSftJmg38YcPH3E6qh/8jSVMkvRM4fmimpNNz10+RGqtfyI8qFgELJP03Ul33nsBWYKekU4C3NSz7GLC/pNeO8lnXAadJOknSVFIHiueAf5d0hKR5+U/jF6Szu6ppnLCiAzpARHwGuBD4M1LGbQD+gFQfehmp0ete0un6XXlaFVeS6se2S/rHUba9g9QKfgrpVPOLwG9FxI8n/IWs0WdIvUdukfQ0qUHzV/O8L7Or99IDeV6jK4CjJW2TdP14NxwR64F3AH9O+l09TDqgmx5TEfE0qWFvZf79DI53+70sIj5N6gjw16SAu5J03J0UEc/lxT5Bqq74T1K99Vca1t9BKuG/n9QL5l3ANxs2cTjwL6Q2rduBL0bEioppu4/UuP2nOR/+iBSYtwHvAW5sWPbHpELfupxPBw37rJ8A/wv4G9Lx/XZSd9sdpD+KRXn6o6SG9kuqpLEVGrstwczM+kXxJXQzs7pwQDczK4QDuplZIRzQzcwKMakXOJoxY0YMDAxM5iZtBHfeeefjETGzXZ/nfO0NztdyVc3bSQ3oAwMDrFq1ajI3aSOQ9NPmS1XnfO0NztdyVc1bV7nUVB74sLrh8ZSkD0uaLmm5pLX5eb9up9XMqnFAr6mI+ElEHBMRxwC/QrrWyA3AQuCWiDicNEy72XUyzKxHOKAbwEnAQxHxU9KNP5bm6UsZ/To4ZtZjHNAN4N2kIc6QLlK0GSA/HzDSCpIWKN3Wb9XWrVsnKZlmNhYH9JrLtwQ7g3TxpMoiYklEDEbE4MyZbetYYWYtcEC3U4C7ImLout+PSZoFkJ+3dC1lZjYuDuh2HruqWyBdbe78/Pp84FuTniIzmxAH9BqT9GrS/RAbL026CJgvaW2et6gbaTOz8as0sEjSetLNjV8AdkbEoKTppLtxD5Bu6ntuRGzrTDKtE/K9LvcfNu0JUq8XM+sz4ymhvyX3Wx66GH+x/ZUHFi7rdhKsz/k31BkDC5d5346hlSoX91c2M+shVQN6kO6wfaekBXma+yubmfWQqhfnOjEiNkk6AFguqfI9MSNiCbAEYHBw0Pe7MzPrkEol9IjYlJ+3kK73cTyF91d2PZ2Z9ZumAV3S3pL2HXoNvA24H/dXNjPrKVWqXA4EbpA0tPzXIuJmSXcA10m6AHgYOKdzyTQzs2aaBvSIWAccPcJ091c2M+shHilqZlYIB3Qzs0I4oJuZFcIB3cysEA7oZoWRNE3S9ZJ+LGmNpP/hm3/XgwO6WZv1wKC0zwE3R8QvkXqoraHgi+nZLg7oZgWR9BrgTcCVABGxIyK244vp1YIDullZXg9sBb4s6W5JV+QR3pUupmf9zQHdrCxTgOOAv42IY4GfMY7qFV8dtb85oJuVZSOwMSJW5vfXkwJ8pYvpRcSSiBiMiMGZM2dOSoKtfRzQzQoSEY8CGyQdkSedBPwIX0yvFqpeD93M+scfAl+VtAewDvgAqfDmi+kVzgE9G1i4jPWLTut2MsxaFhGrgcERZvlieoVzlYuZWSEc0GvMIwrNyuKAXm8eUdhhPTBq1GrEAb2mPKLQrDwO6PXV0ohCD0Ax6z0O6PXV0ohCD0Ax6z0O6PXV0ohCM+s9Dug15RGFVoqBhcvc+Jx5YFG9eUShWUEc0GvMIwrNyuIqFzOzQjigm1mR6li37oBu1iF1CybWfQ7oZmaFcEA3MyuEA/oYfMps4+Hfi3WbA7qZWSEc0M3MCuGAbmZWiMoBXdLu+TKrN+X3h0pame9sc20ePm5Wa65Ht24aTwn9Q6Q72gy5HPhsvrPNNuCCdibMzMzGp1JAlzQbOA24Ir8XMI90yVXwnW3MzLquagl9MXAR8GJ+vz+wPSJ25vcbgYNHWtF3trHSuZrFekXTgC7pdGBLRNzZOHmERWOk9Xv5zjY+EK1EktZLuk/Sakmr8rTpkpbnNq/lkvbrdjqt/aqU0E8EzpC0HriGVNWyGJgmaejyu7OBTR1JoZlNxFsi4piIGLo88kLgltzmdQvjuN2g9Y+mAT0iPhoRsyNiAHg38L2IeC9wK3B2Xsx3tjHrbWeS2rrAbV7FaqUf+sXAhZIeJNWpX9meJJlZiwL4rqQ7JS3I0w6MiM0A+fmArqXOOmZcAT0iVkTE6fn1uog4PiIOi4hzIuK5ziRx8oxUp151mlkPOTEijgNOAX5f0puqrtjLnRh83DXnkaJmhYmITfl5C3ADcDzwmKRZAPl5yyjr9mwnBmvOAd2sIJL2lrTv0GvgbcD9wI2kti5wm1exah/Qh5/GuYrF+tyBwG2S7gF+CCyLiJuBRcB8SWuB+fm9FWZK80WsVLkr6tPAC8DOiBiUNB24FhgA1gPnRsS2bqXRxici1gFHjzD9CeCkyU+RTabal9DN/ZUnymdu1mtcQrfhzgTm5tdLgRWkLqpmPWO0P9Oh6esXnTaZyekZLqGPU2Glsgn3V+7l7m1mdeUSer2dGBGbJB0ALJf046orRsQSYAnA4ODgiNfxqZPC/uitT7mEXmOt9Fc2s97jgF5T7q9sVh5XudTXgcAN6V4lTAG+FhE3S7oDuE7SBcDDwDldTKOZjYMDek25v7JZeVzlYmZWCAd0M7NCOKBX5G5pZtbrHNDNzArhgG5mVggHdDOzQtS22+JE6sRdj242+XzcVecSuplZIRzQzcwK4YBu1mGuMpgcjfu5rvvcAd3MrBAO6GZmhXBAN2tBXU/trTc5oJuZFcIBHZeyzKwMDuhmBZK0u6S7Jd2U3x8qaaWktZKulbRHt9M4WepUYHNANyvTh4A1De8vBz4bEYcD24ALupIq6ygHdLPCSJoNnAZckd8LmAdcnxdZCpzVndRZJzmgm5VnMXAR8GJ+vz+wPSJ25vcbgYNHWlHSAkmrJK3aunVr51NqbeWAblYQSacDWyLizsbJIywaI60fEUsiYjAiBmfOnNmRNFrnNL3aoqS9gO8De+blr4+Iv5R0KHANMB24C3hfROzoZGLNrKkTgTMknQrsBbyGVGKfJmlKLqXPBjZ1MY3WIVVK6M8B8yLiaOAY4GRJJ+BGliK4N0RZIuKjETE7IgaAdwPfi4j3ArcCZ+fFzge+1aUkWgc1DeiRPJPfTs2PwI0spXBviHq4GLhQ0oOkOvUru5yepgYWLqtVl8N2qFSHnktxq4EtwHLgIdzI0vfcG6JsEbEiIk7Pr9dFxPERcVhEnBMRz3U7fdZ+lQJ6RLwQEceQ6t6OB44cabFR1nUjS+9ybwizgoyrl0tEbAdWACeQG1nyLDey9Bn3hrA6Kr0ap2lAlzRT0rT8+lXAW0l1rn3byFJyho7DUG+I9aTeSvNo6A2Rl/Ef9TD+7Vgvq1JCnwXcKule4A5geUTcRB82stgu7g1hVp6m/dAj4l7g2BGmryPVp1tZLgaukXQZcDf+o7Yu8JnQxDQN6Fa+iFhBahvxH7VZH/PQfzOzQjigm5kVwgHdzKwQDuhmZoWoXUB367mZlRoHahfQzcxK5YBuZlYIB3Qzs0LUJqCXWmdmZjakNgHdzKx0DuhmZoVwQDczK4QDutkEjLdNpvQbK1hvcEA3MyuEA7qZWSEc0M3MClGrgO46TDMrWa0CulnpJO0l6YeS7pH0gKRP5OmHSlopaa2kayXt0e20Wvs5oJuV5TlgXkQcDRwDnCzpBOBy4LMRcTiwDbigi2m0DnFAN6ugsbqul6vuInkmv52aHwHMA67P05cCZ3UhedZhDugT1MsHtdWbpN0lrQa2AMuBh4DtEbEzL7IROHiUdRdIWiVp1datWycnwdY2DuhmhYmIFyLiGGA2cDxw5EiLjbLukogYjIjBmTNndjKZ1gEO6DXlxrPyRcR2YAVwAjBN0pQ8azawqVvpss5xQK8vN54VSNJMSdPy61cBbwXWALcCZ+fFzge+1Z0UWic5oNeUG8+KNQu4VdK9wB3A8oi4CbgYuFDSg8D+wJVdTKN1yJTmi1ipJO0O3AkcBnyBcTaeAQsA5syZ0/nEWiURcS9w7AjT15Hq061gLqHXmBvPzMrigG5uPDMrhAN6TbnxzKw8rkOvr1nA0lyPvhtwXUTcJOlHwDWSLgPuxo1nZn3DAb2m3HhmVp6mVS6SDpF0q6Q1eQDKh/L06ZKW5wEoyyXt1/nkmpnZaKrUoe8EPhIRR5IazX5f0lHAQuCWPADllvy+Vnw9FzPrJU2rXCJiM7A5v35a0hpS3+Qzgbl5saWkXhIXdySVZmYtqEvha1y9XCQNkOpdVwIH5mA/FPQPGGUdX73NzGwSVA7okvYBvgF8OCKeqrqeB6CYmU2OSgFd0lRSMP9qRHwzT35M0qw8fxbp2ss9qS6nW9ZZ/h2VqaR8rdLLRaS+yGsi4jMNs24kDTwBD0AxM+u6Kv3QTwTeB9yX74ICcAmwCLhO0gXAw8A5nUmimZlVUaWXy22ARpl9UnuTY2Y2OUqqahnia7mYmRXCAd3MrBAO6GZmhXBANzMrhAO6mVkhHNDNzArhgG5mNoKBhcv6rmujA7qZWSEc0NugH//Jzaw8Duhmk6jTf/z9focxF45a44BuVhbfYazGHNDNChIRmyPirvz6aaDxDmNL82JLgbO6k0LrJAf0mur3U3NrzncYqx8H9PryqXnBfIexenJArymfmper3+8wZhPngG4+NS+I7zBWbw7oNedT8+IM3WFsnqTV+XEq6Q5j8yWtBebn91aYKregs0KNdWoeEZt9at5/fIexenMJvaZ8am62y1gDmvppoJNL6PXlm39X1E8HtNVb8QG90wdjvx7sPjU3K4+rXMzMCuGAbmZWCAd0M7NCFBfQ+7VO28ysVcUFdDOzunJANxtDJ874fBZpneKAbmZWCAd0M+uqXrrtXC+lZSIc0M3MCuGAbmZWCAd0M7NCNA3okr4kaYuk+xum9fR9J/u9HszMbCKqlNCvAk4eNs33nTQz6zFNA3pEfB94cthk33fSzIo1/Ay/X876J1qHXum+k+B7T5qZTZaON4p2896Tk/2P2g//4Fad89P6zUQD+mP5fpP4vpNmZr1hogHd9500s5b5LKi9qnRb/DpwO3CEpI35XpOLgPmS1gLz83szM+uipvcUjYjzRpnVM/edHFi4jPWLTut2MszMusojRc0K04+DAftRL3ZldECvKR/03dfBgHAVHgxYSw7o9XUVPuiL5MGA9eWA3gG9dho2Eh/0tVNpMGA3BwL2w3HT6xzQrZFHANdcNwcCWusc0G1CfOD3HQ8GrIGm3RatVh6TNCsiNvugL87QYMBFeDBgS3q5aqivS+i9vGP7lEcAF8CDAevLJfSaygf9XGCGpI3AX5IO8utyAHgYOKd7KeyuThcWOvn5/TAY0DrDAb2mfNCblaevq1zMzGwXB/Q2azyVdh2/mU0mB3Qzs4p6vZDmgG5mVggHdDOzQvR9L5ehU6BePxUyM+u0vg/oZmbd1lig7ObNdlzlYmZWCAd0s2Emu/rO1YXWLn0b0PvlIOiXdJpZb5lI7OjbgG5mZi/ngG6159G9Vgr3cjGzSeE/y1ca2ift6hnTVyX0fupz7lKfmU22vgroZmY2Ogd0M7MWNDsDn8wzdAd0M7NC9E1A76f6czOzbnAvF6utgYXLXupd4IKCtUs3q2D6poRuZmZjc0A3MyuEA/okGli4zKf2PaAxH3olPxrT0ytpsv7TUh26pJOBzwG7A1dExKK2pMq6yvlarm7krf+gmmvXPppwCV3S7sAXgFOAo4DzJB3VllRZ1zhfy+W8LV8rVS7HAw9GxLqI2AFcA5zZnmRZFzlfy+W8LVwrVS4HAxsa3m8Efm34QpIWAAvy2+ck3d/CNnvBDODx8a6ky0d+3SVHjDHP+dolbfiNjJWvUCFvna+dMdFjvmG9ZnkLtBbQNcK0eMWEiCXAEgBJqyJisIVtdl0p32Gs2SNMc772gSb5ChXy1vnamyrkLdBalctG4JCG97OBTS18nvUG52u5nLeFayWg3wEcLulQSXsA7wZubE+yrIucr+Vy3hZuwlUuEbFT0h8A3yF1gfpSRDzQZLUlE91eDyn6Ozhf+9qY32ECeVv8Pukjlb6HIl5RPWpmZn3II0XNzArhgG5mVoi2BHRJJ0v6iaQHJS0cYf7rJN0i6V5JKyTNbph3vqS1+XF+O9IzUS1+jxckrc6PrjQ0SfqSpC2j9R1W8vn8/e6VdFzDvBHzoYS87fd8zeloe9422d6Y+6zXSTpE0q2S1kh6QNKHup2miZK0u6S7Jd3UdOGIaOlBalx5CHg9sAdwD3DUsGX+ATg/v54HfCW/ng6sy8/75df7tZqmyf4e+f0z3Uj3sPS9CTgOuH+U+acC/0zqj3wCsHKsfCghb0vI107kbav7rNcfwCzguPx6X+A/+u07NHyXC4GvATc1W7YdJfQqw4mPAm7Jr29tmP+bwPKIeDIitgHLgZPbkKaJaOV79ISI+D7w5BiLnAn8fSQ/AKZJmsXo+VBC3vZ9vkJH8nYsfX+JgIjYHBF35ddPA2tII2X7Sj5bPA24osry7QjoIw0nHr7j7gH+Z379DmBfSftXXHeytPI9APaStErSDySd1dmkTtho33G80xv1et7WIV+htTys+ll9SdIAcCywsrspmZDFwEXAi1UWbkdArzJU/E+AN0u6G3gz8Aiws+K6k6WV7wEwJ9IQ4/cAiyW9oWMpnbjRvuN4pzfq9bytQ75Ca3lY9bP6jqR9gG8AH46Ip7qdnvGQdDqwJSLurLpOOwJ60+HEEbEpIt4ZEccCH8vT/qvKupOole9BRGzKz+uAFaQSQa8Z7TuOd/pL+iBv65Cv0EIejuOz+oqkqaRg/tWI+Ga30zMBJwJnSFpPqvaaJ+nqMddoQ4X9FFJDy6HsakD55WHLzAB2y68/CVwauxps/pPUWLNffj29Sw0PrXyP/YA9G5ZZS5caYIABRm84O42XN5z9cKx8KCFvS8nXdudtq/us1x95P/w9sLjbaWnT95lLhUbRdm3sVFIr8kPAx/K0S4Ez8uuz88HwH6TK/T0b1v0g8GB+fKDLO21C3wP4deC+/MO/D7igS+n/OrAZeJ5UyroA+F3gd/N8kW5w8FBO52CzfCghb/s9XzuVt+PdZ/30AN5Iqia6F1idH6d2O10tfJ+5VAjoHvpvZlYIjxQ1MyuEA7qZWSEc0M3MCuGAbmZWCAd0M7NCOKCbmRXCAd3MrBD/Hw1LGMnk4ugzAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "odds summary: 1.906550111317809 0.2971471582695201 [0.75098136 1.51957027 1.82655737 2.20790153 3.1427537 ]\n"
     ]
    }
   ],
   "source": [
    "p0=beta(635,39)\n",
    "p1=beta(658,22)\n",
    "p0_sim=p0.rvs(1000)\n",
    "p1_sim=p1.rvs(1000)\n",
    "fig,ax=plt.subplots(1,3)\n",
    "p0hist=ax[0].hist(p0_sim,density=True,bins=50)\n",
    "ax[0].set_xlim(0.9,1)\n",
    "ax[1].set_xlim(0.9,1)\n",
    "c=ax[0].set_title('Control')\n",
    "p1hist=ax[1].hist(p1_sim,density=True,bins=50)\n",
    "t=ax[1].set_title('Treatment')\n",
    "odds_ratio=(p1_sim/(1-p1_sim)/(p0_sim/(1-p0_sim)))\n",
    "ax[2].set_xlim(0,4)\n",
    "ax[2].hist(odds_ratio,bins=50)\n",
    "o=ax[2].set_title('Odds Ratios')\n",
    "plt.show()\n",
    "print('odds summary:',np.mean(odds_ratio),np.var(odds_ratio),np.percentile(odds_ratio,[.025,25,50,75,97.5]))\n"
   ]
  },
  {
   "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.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}