diff --git a/BDA 3.10.11.ipynb b/BDA 3.10.11.ipynb index 62dec8a..d0e49b9 100644 --- a/BDA 3.10.11.ipynb +++ b/BDA 3.10.11.ipynb @@ -179,7 +179,11 @@ }, { "cell_type": "code", +<<<<<<< HEAD "execution_count": 5, +======= + "execution_count": 9, +>>>>>>> f8632e30fe8b89896ed0fbf258653883205028de "metadata": {}, "outputs": [], "source": [ @@ -205,7 +209,11 @@ }, { "cell_type": "code", +<<<<<<< HEAD "execution_count": 6, +======= + "execution_count": 10, +>>>>>>> f8632e30fe8b89896ed0fbf258653883205028de "metadata": {}, "outputs": [ { @@ -240,7 +248,11 @@ }, { "cell_type": "code", +<<<<<<< HEAD "execution_count": 7, +======= + "execution_count": 11, +>>>>>>> f8632e30fe8b89896ed0fbf258653883205028de "metadata": { "scrolled": true }, @@ -282,6 +294,7 @@ }, { "cell_type": "code", +<<<<<<< HEAD "execution_count": 9, "metadata": {}, "outputs": [], @@ -293,17 +306,27 @@ { "cell_type": "code", "execution_count": 16, +======= + "execution_count": 58, +>>>>>>> f8632e30fe8b89896ed0fbf258653883205028de "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ +<<<<<<< HEAD "1.437043700373185 [0.68945124 0.99000008 1.26005551 1.65286899 3.30679651]\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" +>>>>>>> f8632e30fe8b89896ed0fbf258653883205028de ] } ], "source": [ +<<<<<<< HEAD "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,4,200)),p=p,size=10000)\n", @@ -338,233 +361,45 @@ "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" +>>>>>>> f8632e30fe8b89896ed0fbf258653883205028de ] }, { "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..d888741 --- /dev/null +++ b/BDA 3.10.4.ipynb @@ -0,0 +1,163 @@ +{ + "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": [ + "
" + ] + }, + "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": [ + "
" + ] + }, + "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/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 +} 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 +} 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 dcfcd05..054752c 100644 --- a/gelman.r +++ b/gelman.r @@ -34,5 +34,3 @@ ldens} sd <- rep (NA,nsim) for (i in (1:nsim)) sd[i] <- exp (sample(logsdgrid, 1, prob=dens[muindex[i],])) print (rbind (summ(mu),summ(sd))) - -