From eafa025d180003c3fc36d9a5949742d0b7cef835 Mon Sep 17 00:00:00 2001 From: jeremyteitelbaum Date: Thu, 12 Apr 2018 12:28:14 -0400 Subject: [PATCH 1/3] Did another problem --- BDA 3.10.11.ipynb | 310 +++++----------------------------------------- BDA 3.10.4.ipynb | 168 +++++++++++++++++++++++++ Untitled.ipynb | 82 ------------ gelman.r | 3 +- 4 files changed, 203 insertions(+), 360 deletions(-) create mode 100644 BDA 3.10.4.ipynb delete mode 100644 Untitled.ipynb diff --git a/BDA 3.10.11.ipynb b/BDA 3.10.11.ipynb index 4c37fc8..a7e4f65 100644 --- a/BDA 3.10.11.ipynb +++ b/BDA 3.10.11.ipynb @@ -179,7 +179,7 @@ }, { "cell_type": "code", - "execution_count": 15, + "execution_count": 9, "metadata": {}, "outputs": [], "source": [ @@ -205,7 +205,7 @@ }, { "cell_type": "code", - "execution_count": 16, + "execution_count": 10, "metadata": {}, "outputs": [ { @@ -240,7 +240,7 @@ }, { "cell_type": "code", - "execution_count": 17, + "execution_count": 11, "metadata": { "scrolled": true }, @@ -282,301 +282,57 @@ }, { "cell_type": "code", - "execution_count": 225, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "array([[-2. , -2. , -2. , ..., -2. ,\n", - " -2. , -2. ],\n", - " [-1.9798995 , -1.9798995 , -1.9798995 , ..., -1.9798995 ,\n", - " -1.9798995 , -1.9798995 ],\n", - " [-1.95979899, -1.95979899, -1.95979899, ..., -1.95979899,\n", - " -1.95979899, -1.95979899],\n", - " ...,\n", - " [ 1.95979899, 1.95979899, 1.95979899, ..., 1.95979899,\n", - " 1.95979899, 1.95979899],\n", - " [ 1.9798995 , 1.9798995 , 1.9798995 , ..., 1.9798995 ,\n", - " 1.9798995 , 1.9798995 ],\n", - " [ 2. , 2. , 2. , ..., 2. ,\n", - " 2. , 2. ]])" - ] - }, - "execution_count": 225, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "Y\n", - "\n" - ] - }, - { - "cell_type": "code", - "execution_count": 231, + "execution_count": 58, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "[0.82616995 1.33837647 1.80934486 2.49570902 5.68990452]\n" + "With Rounding\n", + "mu: 10.40618592964824 0.49106811542637807 [ 9.10552764 10.01005025 10.38693467 10.7638191 11.74371859]\n", + "sigma: 1.3555419607805625 0.5063610652735683 [0.61111977 0.90438284 1.18631786 1.55614414 3.30679651]\n" ] } ], "source": [ - "Psigma=list(map(sum,np.exp(Z2-np.max(Z2))))\n", - "p=Psigma/sum(Psigma)\n", - "s=np.random.choice(np.exp(np.linspace(-2,2,200)),p=p,size=1000)\n", - "print(np.percentile(s,q=[2.5,25,50,75,97.5]))" - ] - }, - { - "cell_type": "code", - "execution_count": 178, - "metadata": {}, - "outputs": [], - "source": [ - "p=[1/3,1/3].append([1/200.0]*198)" - ] - }, - { - "cell_type": "code", - "execution_count": 180, - "metadata": {}, - "outputs": [], - "source": [ - "p=[1/3,1/3]\n", - "p.append([1/200.0]*198)\n" + "print('With Rounding')\n", + "Pmu=list(map(np.sum,np.exp(Z1-np.max(Z1)).T))\n", + "p=Pmu/np.sum(Pmu)\n", + "mu_samples=np.random.choice(np.linspace(3,18,200),p=p,size=5000)\n", + "print('mu:',np.mean(mu_samples),np.var(mu_samples),np.percentile(mu_samples,q=[2.5,25,50,75,97.5]))\n", + "Psigma=list(map(np.sum,np.exp(Z1-np.max(Z1))))\n", + "p=Psigma/np.sum(Psigma)\n", + "sigma_samples=np.random.choice(np.exp(np.linspace(-2,4,200)),p=p,size=5000)\n", + "print('sigma:',np.mean(sigma_samples),np.var(sigma_samples),np.percentile(sigma_samples,q=[2.5,25,50,75,97.5]))\n" ] }, { "cell_type": "code", - "execution_count": 181, + "execution_count": 55, "metadata": {}, "outputs": [ { - "data": { - "text/plain": [ - "[0.3333333333333333,\n", - " 0.3333333333333333,\n", - " [0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005,\n", - " 0.005]]" - ] - }, - "execution_count": 181, - "metadata": {}, - "output_type": "execute_result" + "name": "stdout", + "output_type": "stream", + "text": [ + "Ignoring Rounding\n", + "mu: 10.3954824120603 0.49760868662407504 [ 9.03015075 10.01005025 10.38693467 10.7638191 11.81909548]\n", + "sigma: 1.4352129437474488 0.5752559640752803 [0.68945124 0.99000008 1.26005551 1.65286899 3.30679651]\n" + ] } ], "source": [ - "p" + "print('Ignoring Rounding')\n", + "Pmu=list(map(np.sum,np.exp(Z2-np.max(Z2)).T))\n", + "p=Pmu/np.sum(Pmu)\n", + "s=np.random.choice(np.linspace(3,18,200),p=p,size=5000)\n", + "print('mu:',np.mean(s),np.var(s),np.percentile(s,q=[2.5,25,50,75,97.5]))\n", + "Psigma=list(map(np.sum,np.exp(Z2-np.max(Z2))))\n", + "p=Psigma/np.sum(Psigma)\n", + "s=np.random.choice(np.exp(np.linspace(-2,4,200)),p=p,size=5000)\n", + "print('sigma:',np.mean(s),np.var(s),np.percentile(s,q=[2.5,25,50,75,97.5]))" ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] } ], "metadata": { diff --git a/BDA 3.10.4.ipynb b/BDA 3.10.4.ipynb new file mode 100644 index 0000000..3cdfa04 --- /dev/null +++ b/BDA 3.10.4.ipynb @@ -0,0 +1,168 @@ +{ + "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": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#(a). Since p0 and p1 are binomially distributed and independent,\n", + "# the noninformative prior is the beta(1,1) distribution and the posterior \n", + "# distributions are p0~beta(636,40) and p1~(659,23)." + ] + }, + { + "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": [ + "
" + ] + }, + "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", + "\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 effect is pretty minimal" + ] + }, + { + "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": [ + "
" + ] + }, + "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 +} diff --git a/Untitled.ipynb b/Untitled.ipynb deleted file mode 100644 index 6617670..0000000 --- a/Untitled.ipynb +++ /dev/null @@ -1,82 +0,0 @@ -{ - "cells": [ - { - "cell_type": "code", - "execution_count": 6, - "metadata": { - "collapsed": false - }, - "outputs": [ - { - "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAADXBJREFUeJzt3X+s3Xddx/Hni8rUTAkJu4lL29lGG7HRoct1+Af+IDjp\n5o+yiMmmgaAsTY0V+cO4xSX6ByFs0ShqJrWBocTF/jNmGigZETBAyLR3yxx0c6SpmnYBVwaIkz9G\n2ds/7hk7XHt7vqf33Hvufff5SJqc7/f76fm+29776ud+z/v7+aaqkCT18pJ5FyBJmj3DXZIaMtwl\nqSHDXZIaMtwlqSHDXZIaMtwlqSHDXZIaMtwlqaHvmNeJr7rqqtq1a9e8Ti9JW9LDDz/8papamDRu\nbuG+a9culpaW5nV6SdqSkvznkHFelpGkhgx3SWrIcJekhgx3SWrIcJekhgx3SWrIcJekhgx3SWrI\ncJekhuZ2h6q02ey648Pfev0fd/3CxP3rcS5pVpy5S1JDhrskNWS4S1JDhrskNWS4S1JDhrskNWQr\npDQFWxi1VThzl6SGDHdJashwl6SGDHdJashwl6SGDHdJashWSF3WxlsbZ/U+07ZI2l6p9eDMXZIa\nMtwlqSHDXZIaMtwlqSHDXZIasltGW8JGd5TMqotmvd9TWs2gmXuSfUmeTHIqyR0XGfcTSc4neePs\nSpQkTWtiuCfZBtwD3AjsBW5NsneVcXcDH511kZKk6QyZuV8PnKqq01X1HHAU2H+Bcb8D3A88PcP6\nJEmXYEi4bwfOjG2fHe37liTbgZuB91zsjZIcSLKUZOncuXPT1ipJGmhW3TLvBm6vqucvNqiqjlTV\nYlUtLiwszOjUkqSVhnTLPAXsHNveMdo3bhE4mgTgKuCmJOer6h9mUqUkaSpDwv0EsCfJbpZD/Rbg\n18YHVNXuF14n+RvgQwa7hlrvNsf1akG0tVGb2cRwr6rzSQ4BDwLbgHur6mSSg6Pjh9e5RknSlAbd\nxFRVx4HjK/ZdMNSr6i1rL0uStBYuPyBJDRnuktSQ4S5JDRnuktSQ4S5JDRnuktSQ4S5JDRnuktSQ\n4S5JDRnuktSQ4S5JDfmAbG1pG/3gbGmrcOYuSQ0Z7pLUkOEuSQ0Z7pLUkOEuSQ0Z7pLUkK2Q0iZl\nm6fWwpm7JDVkuEtSQ4a7JDVkuEtSQ4a7JDVkuEtSQ7ZCSutovJ1xPcZLq3HmLkkNGe6S1JDhLkkN\nGe6S1JDhLkkNGe6S1JCtkNKMrXc7o6tFaghn7pLUkOEuSQ0Z7pLU0KBwT7IvyZNJTiW54wLH9yd5\nLMmjSR5J8rrZlypJGmriB6pJtgH3ADcAZ4ETSY5V1eNjwz4GHKuqSnIt8ADwA+tRsCRpsiEz9+uB\nU1V1uqqeA44C+8cHVNWzVVWjzSuBZ2ZbpiRpGkNaIbcDZ8a2zwKvXjkoyc3Au4CrgdfPpDrpAlZr\nNbwcV1S0LVKrmdkHqlX1QFW9Evgl4ANJ/t97JzmQZCnJ0rlz52Z1aknSCkPC/Slg59j2jtG+C6qq\nT7L8E8ErLnDsSFUtVtXiwsLCtLVKkgYaEu4ngD1Jdie5ArgFODY+IMkPJsno9XVAqsqpuSTNycRr\n7lV1Pskh4EFgG3BvVZ1McnB0/DDwK8Cbk3wD+F+W/wOQJM3JoLVlquo4cHzFvsNjr+8G7p5taZKk\nS+UdqpLUkOEuSQ0Z7pLUkOEuSQ0Z7pLUkOEuSQ0Z7pLUkOEuSQ35gGytK1ctlObDmbskNWS4S1JD\nhrskNWS4S1JDhrskNWS3jC7ZWjphfA6qtL6cuUtSQ4a7JDVkuEtSQ4a7JDVkuEtSQ4a7JDVkK6TU\nkAu2yZm7JDVkuEtSQ4a7JDVkuEtSQ4a7JDVkuEtSQ7ZCSlvAWlbLtC3y8uTMXZIaMtwlqSHDXZIa\nMtwlqSHDXZIaMtwlqSFbIaUmhrRLDm2LtH1y63PmLkkNDQr3JPuSPJnkVJI7LnD815M8luSzST6T\n5FWzL1WSNNTEcE+yDbgHuBHYC9yaZO+KYf8O/ExV/SjwDuDIrAuVJA03ZOZ+PXCqqk5X1XPAUWD/\n+ICq+kxVfWW0+RCwY7ZlSpKmMSTctwNnxrbPjvat5q3AR9ZSlCRpbWbaLZPktSyH+2tWOX4AOABw\nzTXXzPLUkqQxQ8L9KWDn2PaO0b5vk+Ra4L3AjVX1zIXeqKqOMLoev7i4WFNXK2nD2Ra5NQ25LHMC\n2JNkd5IrgFuAY+MDklwDfBB4U1V9fvZlSpKmMXHmXlXnkxwCHgS2AfdW1ckkB0fHDwN/CLwC+Ksk\nAOeranH9ypYkXcyga+5VdRw4vmLf4bHXtwG3zbY0SdKl8g5VSWrIcJekhgx3SWrIcJekhgx3SWrI\ncJekhgx3SWrIcJekhgx3SWrIcJekhnxAtjbMkAc4S5oNZ+6S1JDhLkkNGe6S1JDhLkkNGe6S1JDh\nLkkN2QqpC/KhyP1txL+xX0fz48xdkhoy3CWpIcNdkhoy3CWpIcNdkhoy3CWpIVshNXOu/rj1bMS/\nmW2RG8uZuyQ1ZLhLUkOGuyQ1ZLhLUkOGuyQ1ZLhLUkO2Ql6GVmt7W609zRY2aetx5i5JDRnuktSQ\n4S5JDRnuktSQ4S5JDQ3qlkmyD/hzYBvw3qq6a8XxVwLvB64D7qyqP5l1oVp/a1k8ysXCLg92Tm0d\nE8M9yTbgHuAG4CxwIsmxqnp8bNiXgbcBb1iXKiVJUxlyWeZ64FRVna6q54CjwP7xAVX1dFWdAL6x\nDjVKkqY0JNy3A2fGts+O9k0tyYEkS0mWzp07dylvIUkaYEM/UK2qI1W1WFWLCwsLG3lqSbqsDAn3\np4CdY9s7RvskSZvUkHA/AexJsjvJFcAtwLH1LUuStBYTu2Wq6nySQ8CDLLdC3ltVJ5McHB0/nOT7\ngCXgZcDzSd4O7K2qr61j7WL6RcCkWbH9dXMb1OdeVceB4yv2HR57/UWWL9dIkjYB71CVpIYMd0lq\nyHCXpIYMd0lqyHCXpIZ8huomMGSlvWlX47NNTZuZLbzrz5m7JDVkuEtSQ4a7JDVkuEtSQ4a7JDVk\nuEtSQ7ZCSto0fAD37Dhzl6SGDHdJashwl6SGDHdJashwl6SGDHdJashWyE1m2hUiN5qrTWoe1mPl\n1O6cuUtSQ4a7JDVkuEtSQ4a7JDVkuEtSQ4a7JDVkK+SYtbRSrfZ71/IgYNsOpYvze2R1ztwlqSHD\nXZIaMtwlqSHDXZIaMtwlqSHDXZIa2pKtkBu9+tu0bY6S1m69v7/W0qa88vdP29q8EbnlzF2SGjLc\nJamhQeGeZF+SJ5OcSnLHBY4nyV+Mjj+W5LrZlypJGmpiuCfZBtwD3AjsBW5NsnfFsBuBPaNfB4D3\nzLhOSdIUhszcrwdOVdXpqnoOOArsXzFmP/CBWvYQ8PIkV8+4VknSQEPCfTtwZmz77GjftGMkSRsk\nVXXxAckbgX1Vddto+03Aq6vq0NiYDwF3VdWnR9sfA26vqqUV73WA5cs2AD8EPDmrP8glugr40pxr\nuBTWvbGse2NZ98V9f1UtTBo0pM/9KWDn2PaO0b5px1BVR4AjA865IZIsVdXivOuYlnVvLOveWNY9\nG0Muy5wA9iTZneQK4Bbg2Ioxx4A3j7pmfhL476r6woxrlSQNNHHmXlXnkxwCHgS2AfdW1ckkB0fH\nDwPHgZuAU8DXgd9Yv5IlSZMMWn6gqo6zHODj+w6PvS7gt2db2obYNJeIpmTdG8u6N5Z1z8DED1Ql\nSVuPyw9IUkOXfbgnecdoyYR/TfLxJNfMu6Yhkvxxkn8b1f5AkpfPu6YhkvxqkpNJnk+yaToLVjNp\n6Y3NKMm9SZ5O8rl51zKNJDuTfCLJ46Ovkd+dd01DJPmuJP8yypAnktw175rAyzIkeVlVfW30+m3A\nq6rqrXMua6IkPw98fPSB990AVXX7nMuaKMkPA88Dfw383sp7ITaT0dIbnwduYPnGvBPArVX1+FwL\nmyDJTwPPsnzX+I/Mu56hRne1X11VjyT5XuBh4A1b4O87wJVV9WySlwKfZvlr+1PzrOuyn7m/EOwj\nVwLPzKuWaVTVR6vq/GjzIZbvLdj0quqJqpr3zWtDDVl6Y9Opqk8CX553HdOqqi9U1SOj1/8DPMEW\nuNN9tOzKs6PNl7LcVfiVOZYEGO4AJHlnkjMst3C+a971XILfBD4y7yIaclmNOUmyC/hx4J/nW8kw\nSbYleRR4Gvinqpr7JbHLItyT/GOSz13g136AqrqzqnYC7wf+bL7VvmhS3aMxdwLngfvmV+m3G1K3\ntJok3wPcD7x9xU/Wm1ZVfbOqfozln6B/Kslr513TlnzM3rSq6ucGDr2PTTQDnlR3krcAvwi8rjbR\nhydT/H1vdoOW1dDsjK5Z3w/cV1UfnHc906qqryb5MLAIfGKetVwWM/eLSbJnbHM/8Oi8aplGkn3A\n7wO/XFVfn3c9TQ1ZekMzMvpg8n3AE1X1p/OuZ6gkCy90qyX5bpY/gJ97jtgtk9zP8gqV3wROA79V\nVV+cb1WTJTkFfCcvfgD8UFUdnGNJgyS5GfhLYAH4KvBoVb1+vlWtLslNwLt5cemNd865pImS/D3w\nsyyvUvhfwB9V1fvmWtQASV4DfAr4LMsdVQB/MLpDftNKci3wtyxPll8C/F1V3T3fqgx3SWrpsr8s\nI0kdGe6S1JDhLkkNGe6S1JDhLkkNGe6S1JDhLkkNGe6S1ND/Afs8xd13koSRAAAAAElFTkSuQmCC\n", - "text/plain": [ - "" - ] - }, - "metadata": {}, - "output_type": "display_data" - } - ], - "source": [ - "import matplotlib.pyplot as plt\n", - "from scipy.stats import norm\n", - "x=norm.rvs(size=1000)\n", - "plt.hist(x,bins=100,normed=True)\n", - "plt.show()" - ] - }, - { - "cell_type": "code", - "execution_count": 4, - "metadata": { - "collapsed": false - }, - "outputs": [ - { - "data": { - "text/plain": [ - "998.39023336044352" - ] - }, - "execution_count": 4, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "x" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": true - }, - "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.5.2" - } - }, - "nbformat": 4, - "nbformat_minor": 2 -} diff --git a/gelman.r b/gelman.r index f39d03f..87928f5 100644 --- a/gelman.r +++ b/gelman.r @@ -33,4 +33,5 @@ ldens} mu <- mugrid[muindex] sd <- rep (NA,nsim) for (i in (1:nsim)) sd[i] <- exp (sample(logsdgrid, 1, prob=dens[muindex[i],])) - kprint (rbind (summ(mu),summ(sd))) \ No newline at end of file + + print (rbind (summ(mu),summ(sd))) \ No newline at end of file From 8b473947f84445bdc6c45b9cb9e5ab427ef17001 Mon Sep 17 00:00:00 2001 From: jeremyteitelbaum Date: Thu, 12 Apr 2018 14:06:43 -0400 Subject: [PATCH 2/3] added 3.10.7 --- BDA 3.10.4.ipynb | 11 ++------ BDA 3.10.7.ipynb | 71 ++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 74 insertions(+), 8 deletions(-) create mode 100644 BDA 3.10.7.ipynb diff --git a/BDA 3.10.4.ipynb b/BDA 3.10.4.ipynb index 3cdfa04..d888741 100644 --- a/BDA 3.10.4.ipynb +++ b/BDA 3.10.4.ipynb @@ -31,14 +31,10 @@ ] }, { - "cell_type": "code", - "execution_count": null, + "cell_type": "markdown", "metadata": {}, - "outputs": [], "source": [ - "#(a). Since p0 and p1 are binomially distributed and independent,\n", - "# the noninformative prior is the beta(1,1) distribution and the posterior \n", - "# distributions are p0~beta(636,40) and p1~(659,23)." + "# Here we use the beta(1,1) uniform prior" ] }, { @@ -81,7 +77,6 @@ "ax[2].hist(odds_ratio,bins=50)\n", "o=ax[2].set_title('Odds Ratios')\n", "plt.show()\n", - "\n", "print('odds summary:',np.mean(odds_ratio),np.var(odds_ratio),np.percentile(odds_ratio,[.025,25,50,75,97.5]))\n", "\n" ] @@ -90,7 +85,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Here we use the beta(0,0) improper prior; the effect is pretty minimal" + "# Here we use the beta(0,0) improper prior; the difference from the previous case is small" ] }, { diff --git a/BDA 3.10.7.ipynb b/BDA 3.10.7.ipynb new file mode 100644 index 0000000..545fe3f --- /dev/null +++ b/BDA 3.10.7.ipynb @@ -0,0 +1,71 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Problem 3.10.7\n", + "\n", + "Poisson and binomial distributions: a student sits on a street corner for an hour and records the number of bicycles $ b $ and the number of other vehicles $v$ that go by. Two models are considered: \n", + "\n", + "* The outcomes $b$ and $v$ have independent Poisson distributions, with unknown means $\\theta_b$ and $\\theta_v$ . \n", + "\n", + "* The outcome $b$ has a binomial distribution, with unknown probability $p$ and sample size $b + v$. \n", + "\n", + "Show that the two models have the same likelihood if we define $p = \\theta_b/( \\theta_b +\\theta_v)$.\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 81). CRC Press. Kindle Edition. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This problem has no computational element, it's a fact about poisson distributions. In the first case we have\n", + "$$\n", + "P(b=b_0)=\\frac{\\theta_b^{b_0}e^{-\\theta_b)}}{b_0!}\n", + "$$\n", + "and\n", + "$$\n", + "P(v=v_0)=\\frac{\\theta_v^{v_0}e^{-\\theta_v)}}{v_0!}\n", + "$$\n", + "It's also a fact that the sum of two poisson variables with rates $\\theta_v$ and $\\theta_b$ is poisson with \n", + "rate $\\theta_v+\\theta_b$. \n", + "\n", + "A direct calculation gives \n", + "$$\n", + "P(b=b_0,v=v_0|b_0+v_0=N)=\\binom{N}{b_0}\\frac{\\theta_b^{b_0}\\theta_v^{N-b_0}}{(\\theta_b+\\theta_v)^{N}}\n", + "$$\n", + "which is what we're supposed to show." + ] + }, + { + "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 +} From f8632e30fe8b89896ed0fbf258653883205028de Mon Sep 17 00:00:00 2001 From: jeremyteitelbaum Date: Thu, 12 Apr 2018 16:11:29 -0400 Subject: [PATCH 3/3] trying to do 3.10.8 and 3.10.6 --- BDA 3.10.8.ipynb | 67 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 67 insertions(+) create mode 100644 BDA 3.10.8.ipynb diff --git a/BDA 3.10.8.ipynb b/BDA 3.10.8.ipynb new file mode 100644 index 0000000..b87de8c --- /dev/null +++ b/BDA 3.10.8.ipynb @@ -0,0 +1,67 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Problem 3.10.8\n", + "\n", + "Analysis of proportions: a survey was done of bicycle and other vehicular traffic in the neighborhood of the campus of the University of California, Berkeley, in the spring of 1993. Sixty city blocks were selected at random; each block was observed for one hour, and the numbers of bicycles and other vehicles traveling along that block were recorded. The sampling was stratified into six types of city blocks: busy, fairly busy, and residential streets, with and without bike routes, with ten blocks measured in each stratum. Table 3.3 displays the number of bicycles and other vehicles recorded in the study. For this problem, restrict your attention to the first four rows of the table: the data on residential streets. \n", + "\n", + "(a) Let $y_1$ , . . . , $y_{10}$ and $z_1$ , . . . , $z_8$ be the observed proportion of traffic that was on bicycles in the residential streets with bike lanes and with no bike lanes, respectively (so $y_1 = 16/(16 + 58)$ and $z_1 = 12/(12 + 113)$, for example). Set up a model so that the $y_i$ ’s are independent and identically distributed given parameters $\\theta_y$ and the $z_i$ ’s are independent and identically distributed given parameters $\\theta_z$ . \n", + "\n", + "(b) Set up a prior distribution that is independent in $\\theta_y$ and $\\theta_z$ . \n", + "\n", + "(c) Determine the posterior distribution for the parameters in your model and draw 1000 simulations from the posterior distribution. (Hint: $\\theta_y$ and $\\theta_z$ are independent in the posterior distribution, so they can be simulated independently.) \n", + "\n", + "(d) Let $\\mu_y = E(y_i |\\theta_y )$ be the mean of the distribution of the $y_i$ ’s; $\\mu_y$ will be a function of $\\theta_y$. Similarly, define $\\mu_z$ . Using your posterior simulations from (c), plot a histogram of the posterior simulations of $\\mu_y-\\mu_z$, the expected difference in proportions in bicycle traffic on residential streets with and without bike lanes. We return to this example in Exercise 5.13.\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 81). CRC Press. Kindle Edition. \n", + "\n", + "#### Data\n", + "|Type |Bike lane? |Counts of Bikes/others|\n", + "|--- |----------|----|\n", + "|Residential |yes |16/58, 9/90, 10/48, 13/57, 19/103, 20/57, 18/86, 17/112, 35/273, 55/64 |\n", + "|Residential |no |12/113, 1/18, 2/14, 4/44, 9/208, 7/67, 9/29, 8/154|\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 81). CRC Press. Kindle Edition. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Probably best to first do 3.10.6\n", + "For that problem see the reference [Raftery, 1988](https://www.stat.washington.edu/raftery/Research/PDF/bka1988.pdf)" + ] + }, + { + "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 +}