From 770e7f432c0567f029e56c820b6302eab0c60928 Mon Sep 17 00:00:00 2001 From: Jeremy Teitelbaum Date: Fri, 27 Apr 2018 09:15:32 -0400 Subject: [PATCH 1/8] added 4.7.5 --- BDA 4.7.5.ipynb | 166 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 166 insertions(+) create mode 100644 BDA 4.7.5.ipynb diff --git a/BDA 4.7.5.ipynb b/BDA 4.7.5.ipynb new file mode 100644 index 0000000..1320e4b --- /dev/null +++ b/BDA 4.7.5.ipynb @@ -0,0 +1,166 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Problem 4.7.5.\n", + "\n", + "Approximate mean and variance.\n", + "1. Suppose x and y are independent normally distributed random variables, where x~N(4,1) and y~N(3,2). What are the mean and standard deviations of y/x? Compute this using simulation.\n", + "\n", + "2. Do the same computation without simulation.\n", + "\n", + "3. What assumptions do you need for part (2)?" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "from scipy.stats import norm" + ] + }, + { + "cell_type": "code", + "execution_count": 48, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "mean: 0.8079667283774233 sd: 0.6574088802979349\n" + ] + } + ], + "source": [ + "x=norm(4.0,1.0)\n", + "y=norm(3.0,2.0)\n", + "\n", + "samples_x=x.rvs(10000)\n", + "samples_y=y.rvs(10000)\n", + "z=samples_y/samples_x\n", + "print('mean:',np.mean(z),'sd:',np.sqrt(np.var(z)))" + ] + }, + { + "cell_type": "code", + "execution_count": 57, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD8CAYAAACb4nSYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XuYXFWZ7/Hvr5skQEAIJiBCQgcNl3ANNAFkJujhYkAn0VGPiccRBE/OzHCZGceZE2RGMzBqvMxB58goyOQIoxIRBaJGMdy8kpAGAiGBkBBCCIkkGAgCIaHT7/lj74ZKpS+7u6v2rur6fZ6nnq7al1pvde16a9Xaa6+liMDMzBpDU9EBmJlZfpz0zcwaiJO+mVkDcdI3M2sgTvpmZg3ESd/MrIE46ZuZNRAnfTOzBuKkb2bWQHYrOoByI0eOjJaWlqLDsEHs/vvvfy4iRuVdro9tq6asx3XNJf2Wlhba2tqKDsMGMUlPFVGuj22rpqzHtZt3zLogabKkFZJWSZrZw3YflBSSWvOMz6y/nPTNykhqBq4GzgHGA9Mlje9iu72BS4FF+UZo1n9O+ma7mgisiojVEbEdmAtM7WK7K4EvAa/mGZzZQGRq05c0Gfga0AxcFxGzy9ZfBbwrfbgnsH9E7Juu2wEsTdetjYgplQjcsmuZ+dPX76+Z/Z4CI6kbBwFPlzxeB5xcuoGkCcDoiPiJpE/lGZwlfFz3T69Jv+Sn7lkkB/9iSfMiYnnnNhHxdyXbXwJMKHmKrRFxfOVCNqs6dbHs9YknJDUBVwHn9/pE0gxgBsCYMWMqFJ5Z/2Vp3sn6U7fTdODGSgRnVpB1wOiSxwcD60se7w0cDdwjaQ1wCjCvq5O5EXFtRLRGROuoUbn3EjXbRZbmnV5/6naSdAgwFrirZPHuktqAdmB2RNzaz1jN8rIYGCdpLPAMMA34SOfKiNgCjOx8LOke4FMR4f6YOSht1rG+y5L0e/ypW2YacHNE7ChZNiYi1ks6FLhL0tKIeGKnAvwT2GpIRLRLuhi4neQ81pyIWCbpCqAtIuYVG2FjcHKvjixJv7efuqWmAReVLoiI9enf1WmNaALwRNk21wLXArS2tnrSXitcRMwH5pct+0w3274zj5jMKiFLm/7rP3UlDSVJ7LvUdCQdDowA7i1ZNkLSsPT+SOA0YHn5vmZmlo9ea/p9+Kk7HZgbEaU19SOBayR1kHzBzC7t9WNmZvnK1E8/y0/diJjVxX6/A44ZQHxmZr1yn/3sfEWumVkDqblRNs3MBsK1/p65pm9m1kBc028w5X2fXRMyayxO+oOUL2wxs6446ZvZoOX2/V25Td/MrIE46ZuZNRAnfTOzBuI2fTOrGe6AUH2u6ZuZNRDX9BvI4VrLpKaH2cowbt1xGi+xZ9EhmVnOnPQbxFFaw9yhV7K3tgLw3uaFnLf9fxcclVnCzTr5cfNOAxjFC3x76Gy2MJxJ267ib7f/Nac0PcrsId8qOjQzy5lr+g3gE7v9lP34I5O3f5G1cQBr4wAObV/PpbvdChseggOPKzpEM8uJa/qD3Jt4if/RfCc/7jiVlXHw68uva38Pf4w94DdfLTA6M8ubk/4g97HmBeylV/lm+5Sdlr/IcL674wxYfitsXl1QdGaWNyf9QS2Yttvd/HLHsTwWY3ZZO6f9HEDwwA35h2ZmhXDSH8SO0xMcrOe4bcc7uly/kREwdhIsvw12mtrYzAYrJ/1B7NzmRWyPZu7oOLH7jY56X9K88+wj+QVmZoVx0h+0gvc0L+I3HcfwIsO73eqEHwxjR4j/+/Wv5BibNaqWmT99/WbFyJT0JU2WtELSKkkzu1h/vqRNkpakt0+UrDtP0sr0dl4lg7fuHavVHKznmN9xco/bbeZNLOwYz7lNi3KKzCzhL4Bi9Jr0JTUDVwPnAOOB6ZLGd7Hp9yPi+PR2XbrvfsBngZOBicBnJY2oWPTWrdObHqIjxJ07JvS67S86Wnlb0wbY/GQOkZlZkbJcnDURWBURqwEkzQWmAssz7PtuYEFEbE73XQBMBm7sX7jWk9Ia0/eHPsKyOITneVOv+/2m4+jkzup7YL+xVYrOzGpBluadg4CnSx6vS5eV+4CkhyXdLGl0H/e1ChrOVk7QSn7TcUym7Z+It7Ih9oPVd1c5MrPiuCkpkSXpq4tl5f37fgy0RMSxwB3A9X3YF0kzJLVJatu0aVOGkKwnJzc9yhDt4NcZkz6I3+w4Glb/Ejp2VDU2MytWlqS/Dhhd8vhgYH3pBhHxh4jYlj78FnBi1n3T/a+NiNaIaB01alTW2K0bf9q0lFdjCPd3HJZ5n193HAOvvgAbllQxMjMrWpY2/cXAOEljgWeAacBHSjeQdGBEbEgfTgEeTe/fDny+5OTt2cBlA47aenRa0yPc13EE2xiaeZ/flrbrH9RDv36zOlfaxLNm9nsKjKQYvdb0I6IduJgkgT8K3BQRyyRdIalzQJdLJS2T9BBwKXB+uu9m4EqSL47FwBWdJ3WtOkbwIoc1PcPCjq46WHXvD+wDo46Ap35XpcjMrBZkGlo5IuYD88uWfabk/mV0U4OPiDnAnAHEaH1wUtMKABZ1HNH3ncecCktvTtr1m5orHJmZ1QJfkTvITGx6jG0xhKVxaN93PuQ02P5H+P3SygdmZjXBSX+QOalpBUvibWxnSN93PuTU5O/aeysblJnVDCf9QWQ4WzlaT/avaQdgn4NhnzFu1zcbxJz0B5ETmlbSrGBxf5M+wCHvSGr6HmrZbFBy0h9EWpseZ0eIBzrG9f9JDjkVXt7k2bTMBilPjD6InKDHeSzG8DJ79P9JRqejcj69CN78tsoEVockTQa+BjQD10XE7LL1nwQ+AbQDm4ALIuKp3AOtA7U89EEj9tl3TX+w6NjBhKZVfboKt0sjD4dh+yRJv0FlHFn2QaA1HXrkZuBL+UZp1j9O+oPFxuXspVdpG2jSb2qC0SfB0/dVJq769PrIshGxHegcWfZ1EXF3RLySPlxIMsSIWc1z0h8s0pr5AzHApA9JE8/GR2HrCwN/rvrU19FhLwR+1tUKDyZotcZJf7B4+j6ejX1ZFyMH/lyjTwYC1rUN/LnqU6bRYQEkfRRoBb7c1XoPJmi1xidyB4u1C9NeO13lq2w6T2rtyass36MJnl4I486sUIB1JdPosJLOBC4HTi8ZZdasprmmPxj88ffwwlMDb89PvcLucMDRjdyu//rIspKGkowsO690A0kTgGuAKRGxsYAYzfrFNf3BIE3OD1Qo6QNJE89DN8KOdmhurMMkItoldY4s2wzM6RxZFmiLiHkkzTl7AT+QBLA2IqZ0+6QNqJa7ajayxvo0D1ZPL4LmYTwSFZzfdvTJsPhbsHEZHHhc5Z63TmQYWbYh272s/rl5ZzB4ehEcdAKvVfI7fPTE9LkbtonHbFBy0q93r70K65e8kaQrZd8xsPeBDX2Rltlg5KRf7zYsgY7X3hg+oVKk5IvESd9sUHHSr3NfuObbAJz47S2Vf/LRJ8MLa+HFDb1va2Z1wSdy69xJTY/xRMeByRy3lTbmlOTv2nvh6D+v/POb1ZBGGXzNNf161tFBa9PjLO44vDrP/5ZjYcienknLbBBx0q9nmx5jX708sElTetI8BA4+CZ5y0jcbLDIlfUmTJa2QtErSzC7Wf1LSckkPS7pT0iEl63ZIWpLe5pXvawOwNpnWcHFUqaYPyUxazz7SyIOvmQ0qvSb9CowtvjUijk9vvmKxktYu5NnYl7Wxf/XKGHMqyeBri6tXhpnlJktN32OL16IIeOretGmn/4Os9ergVmjazZOlW0NpmfnTQTuMRJbeO12NLd5Tp/DyscV3l9RGMq3c7Ii4tc9R2q6eXwMvrmNRR5VHAxg6HN46AZ76bXXLsUFhsCbKwSRL0u/P2OKnlyweExHrJR0K3CVpaUQ8UbbfDGAGwJgxYzIF3vCe/BUAv+s4qipPv1P3tbMmwW+/Btv+CMP2rkp5ZpaPLM07fR1bfErp2OIRsT79uxq4B5hQvq8nmuiHJ38Fex3AE/HW6pc1dhJ0tMPahdUvy8yqKkvS7/fY4pJGSBqW3h8JnAYsr1TwDSsiSfpjJ1HV9vxOo0+G5qHw5C+rX5aZVVWvST8i2oHOscUfBW7qHFtcUmdvnNKxxUu7Zh4JtEl6CLibpE3fSX+gnnscXt6YJv0cDNkjSfxpk5KZ1a9MwzD0d2zxiPgdcMxAArQurE5r3GMnAcvyKXPsJLj78/DKZthzv3zKNCvYYByawVfk1qMn7oQRLcktL4e+EwhYfXd+ZZpZxXnAtXrz2takpn/Cx/It96ATYY8RsHIBHP2BfMu2muZumvXFNf16s+a30L4Vxp2db7lNzfD2M5Ok39GRb9lmVjFO+vVm5S9gtz2g5bT8yx73bnjlOVj/YP5lm1lFOOnXkwhYeXtyUnXIHvmX//YzACUxmFldctKvJxsfTYZfOCznpp1Oe+6XTKH42PzetzWzmuQTufVk2S2gJlp/uCfP/bCgk2fjp8Ltn4bnVsHItxcTgxWuEU/eDpbum67p14sIWPYjaPkTnqvG1Ig96BxxsGXmT2H8+5KFy27JNQYzqwzX9OvFs4/AH1bBqRcl10UXZZ+DkjH2l/0ITv+HAgOxvDVi7X4wck2/XjzyQ1AzHFkD89Ac9eewcXlyjsHM6oqTfj3Y8Ros+R6MOwuGjyw6Gjjq/dA0BB64oehIzKyP3LxTD1bMh5eehdYLio4ksdcoOPLPYMl34b/9Mwzds+iIzHJVzyd1XdOvB21zYJ/RyRWxteKkC+HVLT6h2wAG89SBjcg1/Vq38TFYfQ+86/JkKISCvfHhD9YcfBjcdw0c/xFQDuP6m9mAuaZf6375RRi6F7ReWHQkZQTvuAQ2PASP/7zoYMwsIyf9Wrbx0aT5ZOIMGP7moqPZ1XHTk+Gd7/lCch2BmdU8J/1aFQF3/AsMHZ7UqGtR8xCY9I9Jbd9t+9ag6u2ch9v0a9Xy2+Dxn8FZV8Ce+9XuQXXsh5N2/Z/PhLe9Kxlz38xqlpN+LXr5OZj/D3DgcXDKRUVH063OL6Lx+jDzd/9n+Pmn4f3fKDgqq4SarWTUsHrpxunmnVrTvg2+/1HY9iJM+To01/738vJogT/9e3joe7DomqLDMbMe1H5GaSTt2+BHM2DtvfDBOXDgsUVHlN07Z8Kzy5Jmnj3fDMd8sOiIzApTy7V+J/1a8dJGuPkCWPNrOPtf624e2pZP/5w9+ADXD13FxB9eCC88Baf9bU1cW2DdczNO48mU9CVNBr4GNAPXRcTssvXDgBuAE4E/AB+OiDXpusuAC4EdwKUR4WmXSr22FR78Dtx1ZXL//dfCcR+uyw/jVnbnL7ZfxleGfJM/u/MKHlzwPT7/2kf4wec/WXcXbw3kmDcr1dVnucjaf69JX1IzcDVwFrAOWCxpXkQsL9nsQuD5iHi7pGnAF4EPSxoPTAOOAt4K3CHpsIjYUekXUjciYMvTyTyzq+6ER38MWzfDIafBe78Kow4rOsIB2cZQLnntEu7YcQL/NOQ7/GDYFfDNW+DI90LLn8BbjoHd850PoK8GcsznH2129ViRsMrLUtOfCKyKiNUAkuYCU4HSD8BUYFZ6/2bg65KULp8bEduAJyWtSp/v3m5Le/GZpBdIXUgvSIpI7kdAdEDsgI522NEO7Vth+8vw6ovw8qbk9bW/muw3dO9k5MyTPgGHvIOWy+YDK4t6MRUkbuv4E27fdhJ/3vwbPrDhV0z4/WyalPy/NsdePBsjeD725iX2YCvDmHrioclJazWnTUJKfx2o5FdCbr8W+n3MR/gqNetdd1/AefwCyJL0DwKeLnm8Dji5u20iol3SFuDN6fKFZfseVF6ApBnAjPThSzrnCysyRd+1kcBzA9g/x3JfBL6d3vIuuyJ6LfcL6W1nLwLrq1puLw7pZf1Ajvmd4trl2Jb6e2wX9R6Xcxw7q2gc+uKA4ujtuAayJf2uqlfltZnutsmyLxFxLXBthlh6JaktIlor8Vz1UG6RZQ/icgdyzO+8oELHdpHHl+OomzhasmybpZ/+OmB0yeOD2bWa9vo2knYD9gE2Z9zXrNYM5Jg3q2lZkv5iYJyksZKGkpyYnVe2zTzgvPT+B4G70rbNecA0ScMkjQXGAfdVJnSzqhnIMW9W03pt3knbKy8GbifpvjYnIpZJugJoi4h5wH8C/5WeqN1M8iEh3e4mkhNg7cBFOfTcqUgzUR2VW2TZg7LcgRzzVVTk8VXKceys7uKQKydmZo3DY++YmTUQJ30zswYyaJO+pE9JCkkjcyzzy5Iek/SwpFsk7Vvl8iZLWiFplaSZ1SyrrNzRku6W9KikZZL+Jq+y0/KbJT0o6Sd5lluEot7jshgKfb+7iKcm3n9J+0q6Of3MPyrp1ILi+Lv0fXlE0o2Sdu9p+0GZ9CWNJrmEfm3ORS8Ajo6IY4HHgcuqVVDJUAHnAOOB6emwF3loB/4+Io4ETgEuyrFsgL8BHs2xvEIU/B6XKvr9Llcr7//XgJ9HxBHAcRQQk6SDgEuB1og4mqTjQY+dCgZl0geuAv6RLi6WqaaI+EVEtKcPF5L0766W14cKiIjtQOdQAVUXERsi4oH0/h9JDvZdrrSuBkkHA+8BrsujvIIV9h6XKvL9Llcr77+kNwGTSHpxERHbI+KFgsLZDdgjvV5kT3q5FmrQJX1JU4BnIuKhgkO5APhZFZ+/q6ECcv8gSmoBJgCLciryqyRf6B05lVekmniPSxXwfperlff/UGAT8P/SpqbrJA3PO4iIeAb4CkmrxgZgS0T8oqd96jLpS7ojbb8qv00FLgc+U1DZndtcTvKT+LvVioOMwwBUk6S9gB8CfxsRL+ZQ3nuBjRFxf7XLqhGFv8el8n6/uyi/lt7/3YATgG9ExATgZSD3cy6SRpD8+htLMpLxcEkf7WmfupxEJSLO7Gq5pGNIXvxDySCfHAw8IGliRPy+mmWXxHAe8F7gjCpfoVnoEBeShpAkgO9GxI9yKvY0YIqkc4HdgTdJ+k5E9HiQ17GaGcakoPe7XC29/+uAdRHR+YvnZgpI+sCZwJMRsQlA0o+AdwDf6XaPiBi0N2ANMDLH8iaTXH08KoeydgNWk3zJDQUeAo7K6XWKZAKRrxb43r4T+ElR5ef0Ggt7j2vt/a7F9x/4NXB4en8W8OUCYjgZWEbSli/geuCSnvapy5p+Dfs6MAxYkP7SWBgRf1mNgqKboQKqUVYXTgP+AlgqaUm67NMRMT+n8htCwe9xKb/fXbsE+G46PtNq4ON5BxARiyTdDDxA0qT8IL0MydDrMAyS5pA0V2yMpEtQ+XqRdF06F3gFOD/SM/1pU8c/pZv+a0Rc36dXZGZmFZXlRO63SZotunMOyeiZ40gmi/gGgKT9gM+S/PyYCHw2PelgZmYF6TXpR8Sv6Hmc8KnADZFYCOwr6UDg3cCCiNgcEc+TXLjU05eHmZlVWSW6bHbXl7jm+hibmTW6SpzIHdBUibDzPKLDhw8/8YgjjqhAWGZdu//++5+LiFF5lzty5MhoaWnJu1hrEFmP60ok/e76Eq8j6VZVuvyerp4gSuYRbW1tjba2tgqEZdY1SU8VUW5LSws+tq1ash7XlWjemQd8TIlTSC4D3kDSzexsSSPSE7hnp8vMzKwgvdb0Jd1IUmMfKWkdSY+cIQAR8U1gPkl3zVUkXTY/nq7bLOlKkvlGAa6ICE8cbWZWoCxz5E7vZX0AF3Wzbg4wp3+hmZlZpdXlgGtmNkjN2ie5WdU46ZuZNRAnfTOzBuKkbw3lggsuADhO0iOdyyTNkvSMpCXp7dySdZel89OukPTukuWFz11r1h9O+tZQzj//fICVXay6KiKOT2/zAdJ5YKcBR5EMIfIf6aTctTJ3rVmfeWhlayiTJk2CZAjaLBWeqcDciNgGPClpFcnggZDOXQsgqXPu2uWVj9isslzTN0tcLOlhSXNKRoMd8LhSkmZIapPUtmnTpmrEbdYnTvpmyXDgbwOOJ5lc+t/S5QMeVyoiro2I1ohoHTUq9+F+zHbh5h1reBHxbOd9Sd8CfpI+7GmO2pqYu9asr1zTt4aXzv/Q6f1AZ8+eecA0ScMkjSWZKOg+kqFFxkkam06VNy3d1qzmuaZvDWX69OkAR5DM9Nk5ltQ7JR1P0kSzBvhfABGxTNJNJCdo24GLImIHyc61MHetWZ856VtDufHGG5k7d+7DEdFasvg/u9s+Ij4HfK6L5fNJBhs0qytu3jGz2lA65o7H36kaJ30zswbipG9m1kCc9M3MGkimpN/b4FKSrioZrOpxSS+UrNtRss7d2szMCpRlusTOwaXOIrlYZbGkeRHx+jgjEfF3JdtfAkwoeYqtEXF85UI2M7P+ylLTn0g6uFREbAc6B5fqznTgxkoEZ2ZmlZUl6fdlcKlDgLHAXSWLd08HnFoo6X39jtTMzAYsy8VZmQeXIrkc/ebOqxZTYyJivaRDgbskLY2IJ3YqQJoBzAAYM2ZMhpDMzKw/stT0exp0qtw0ypp2ImJ9+nc1cA87t/d3buORCM3McpAl6WcaXErS4cAI4N6SZSMkDUvvjwROwxNNmJkVptfmnYho72pwKUlXAG0R0fkFMJ1klqHSpp8jgWskdZB8wcwu7fVjZmb5yjTgWleDS0XEZ8oez+piv98BxwwgPjMzqyBfkWtm1kCc9M3MGoiTvplZA3HSN7PaNGsfj6tfBU76ZmYNxEnfzKyBOOmbmTUQJ31rKBdccAHAcZIe6VwmaT9JCyStTP+OSJdL0r+n80g8LOmEkn3OS7dfKem8/F/JIOO2+9w46VtDOf/88wFWli2eCdwZEeOAO9PHAOcA49LbDOAbkHxJAJ8FTiYZevyznV8UZrXOSd8ayqRJkwDayxZPBa5P718PvK9k+Q2RWAjsK+lA4N3AgojYHBHPAwuAyVUP3qwCnPTN4ICI2ACQ/t0/Xd7dXBKZ55gwqzVO+mbd624uicxzTEiakU4i1LZp06aKBmfWH076ZvBs2mxD+ndjury7uSQyzzHhuSKs1jjpmyXzQ3T2wDkPuK1k+cfSXjynAFvS5p/bgbPT+SJGAGeny8xqXqahlc0Gi+nTpwMcQdIjcx1JL5zZwE2SLgTWAh9KN58PnAusAl4BPg4QEZslXUkywRDAFRGxObcXYTYATvrWUG688Ubmzp37cES0lq06o3zbdEKgi7p6noiYA8ypQohmVZWpeUfSZEkr0otUZnax/nxJmyQtSW+fKFnni1jMzGpErzV9Sc3A1cBZJCewFkua18W0h9+PiIvL9u28iKWVpHfD/em+z1ckejMz65MsNf2JwKqIWB0R24G5JBetZOGLWMzMakiWpJ/1QpQPpOOT3Cypsztbpn3dl9nMLB9Zkn6WC1F+DLRExLHAHbxxSXumi1jcl9nMLB9Zkn6vF6JExB8iYlv68FvAiVn3NTOz/GRJ+ouBcZLGShoKTCO5aOV1nVczpqYAj6b3fRGLmVkN6bX3TkS0S7qYJFk3A3MiYpmkK4C2iJgHXCppCsnohZuB89N9fRGLmVkNyXRxVkTMJ7k6sXTZZ0ruXwZc1s2+vojFzPpv1j4wa0vRUQwaHnvHzKyBOOmbmTUQJ30zswbipG9m1kCc9M3MGoiTvplZA/F4+mZWnFn7FB1Bw3FN38ysgTjpm5k1ECd9M7MG4qRvZtZAnPTNUpLWSFqazvPcli7bT9KCdI7nBelosSjx7+m80Q9LOqHY6M2ycdI329m7IuL4iGhNH88E7oyIccCd6WOAc4Bx6W0G8I3cIzXrByd9s55N5Y2Z4K4H3ley/IZILAT2LZtXwqwmOembvSGAX0i6X9KMdNkBEbEBIP27f7o869zRVgnuz18xmZK+pMmSVqTtlzO7WP9JScvTts07JR1Ssm5H2ka6RNK88n3NashpEXECSdPNRZIm9bBtpvmfJc2Q1CapbdOmTZWK06zfek36kpqBq0k+COOB6ZLGl232INCaTox+M/ClknVb0zbS4yNiSoXiNqu4iFif/t0I3AJMBJ7tbLZJ/25MN880/3NEXBsRrRHROmrUqGqGb5ZJlpr+RGBVRKyOiO3AXJL2zNdFxN0R8Ur6cCHJB8CsbkgaLmnvzvsk8zk/QjIf9HnpZucBt6X35wEfS3vxnAJs6WwGMqtlWcbe6art8uQetr8Q+FnJ493T7m/twOyIuLXPUZpV3wHALZIg+Vx8LyJ+LmkxcJOkC4G1wIfS7ecD5wKrgFeAj+cfslnfZUn6mdouASR9FGgFTi9ZPCYi1ks6FLhL0tKIeKJsvxkk3d4YM2ZMpsDNKikiVgPHdbH8D8AZXSwP4KIcQjOrqCzNO5naLiWdCVwOTImIbZ3LS9pJVwP3ABPK93W7p5lZPrIk/cXAOEljJQ0FppG0Z75O0gTgGpKEv7Fk+QhJw9L7I4HTgOWVCt7MzPqm1+adiGiXdDFwO9AMzImIZZKuANoiYh7wZWAv4Adpm+jatKfOkcA1kjpIvmBmR4STvplZQTJNohIR80lOXJUu+0zJ/TO72e93wDEDCdDMzCrHV+SamTUQJ30zK4aHViiEk76Z1YdZ+/iLogKc9M3MGoiTvplZA3HSNzNrIE76ZmYNxEnfzOqLT+YOiJO+mVkDcdI3M2sgTvpmZg3ESd/MrIE46ZtZvipxZa1P5vabk76ZWQNx0jczayBO+mZmDSRT0pc0WdIKSaskzexi/TBJ30/XL5LUUrLusnT5CknvrlzoZsXq7XNhXahkW3znuQG37/dJr0lfUjNwNXAOMB6YLml82WYXAs9HxNuBq4AvpvuOJ5lT9yhgMvAf6fOZ1bWMnwuzmpOlpj8RWBURqyNiOzAXmFq2zVTg+vT+zcAZSibLnQrMjYhtEfEksCp9PrN6l+VzYXlxrT+zLHPkHgQ8XfJ4HXByd9ukE6lvAd6cLl9Ytu9B/Y7WrHZk+VxYEUm4tMxZW954PGtL/rHUoCxJX10si4zbZNkXSTOAGenDbZIeyRBXNYwEnmugcossu8jXfHgFnqM/x/ZLklak94t8/Z1qIQaoZhz/oq7v5xlD3wwkjkOybJQl6a8DRpc8PhhY38026yTtBuwDbM64LxFxLXAtgKS2iGjNEnxuYyE+AAAFKUlEQVSlFVW2X3P+ZVfgafp8bJfHUNTrr6UYaiWOWoghrziytOkvBsZJGitpKMmJ2Xll28wDzkvvfxC4KyIiXT4t7d0zFhgH3FeZ0M0KleVzYVZzeq3pp230FwO3A83AnIhYJukKoC0i5gH/CfyXpFUkNfxp6b7LJN0ELAfagYsiYkeVXotZbrr7XBQcllmvsjTvEBHzgfllyz5Tcv9V4EPd7Ps54HN9iGmXn8I5Kqpsv+Y6LLurz0XeMQxQLcQAtRFHLcQAOcShpBXGzMwagYdhMDNrIIUl/YEM7ZBD2Z+UtFzSw5LulJSpK9RAyy3Z7oOSQlJFzuJnKVfSf09f8zJJ36tEuVnKljRG0t2SHkz/3+dWqNw5kjZ21/1XiX9P43pY0gmVKLeXmD6U/n87yt/booYrkTRL0jOSlqS3ivz/M5ZdE8NYSFojaWn6+ivRsytLmbscn5L2k7RA0sr074iqFB4Rud9ITnw9ARwKDAUeAsaXbfPXwDfT+9OA7+dY9ruAPdP7f1WJsrOUm263N/ArkovaWnN6veOAB4ER6eP9c/xfXwv8VXp/PLCmQmVPAk4AHulm/bnAz0j6258CLKr0cd5FmUeSXCNwT+l7m77uh4BhwNj0f9Zc7XjSsmcBn8qjrL4eGznGsgYYmXOZuxyfwJeAmen9mcAXq1F2UTX9gQztUPWyI+LuiHglfbiQpA921ctNXUny5r9agTKzlvs/gasj4nmAiNiYY9kBvCm9vw9d9HXvj4j4FUlPsu5MBW6IxEJgX0kHVqLsHmJ6NCJWdLGqEYcraehhLLo5Pktz3vXA+6pRdlFJv6tL2MuHZ9hpaAegc2iHPMoudSFJjbDq5UqaAIyOiJ9UoLzM5QKHAYdJ+q2khZIm51j2LOCjktaR9IS5pEJl96avx8FgjuXitIlrTtWaFHZV9GsuFcAvJN2fXkFdlAMiYgNA+nf/ahSSqctmFQxkaIc8yk42lD4KtAKnV7tcSU0kI5SeX4GyMpeb2o2kieedJL9qfi3p6Ih4IYeypwPfjoh/k3QqyfUeR0dExwDLrkRsfX9S6Q7gLV2sujwibsszliwxAd8g+XUZ6d9/Ay6oVNk9hdXFsqK6Ep4WEesl7Q8skPRYWhMflIpK+gMZ2iGPspF0JsmH4vSI2JZDuXsDRwP3pK1YbwHmSZoSEQM5uZT1f70wIl4DnlQyPsw4kqtOByJL2ReSDLtNRNwraXeS8Ucq1cQ0kNj6LCLOrJVYOmWNSdK3gEr+yuxJVV9zX0TE+vTvRkm3kDQ9FZH0n5V0YERsSJsaq/IZKKp5ZyBDO1S97LSZ5RpgSgXbt3ssNyK2RMTIiGiJiBaScwkDTfi9lpu6leTkNZJGkjT3rB5guVnLXguckZZ9JLA7sKkCZfdmHvCxtBfPKcCWzp/WBShsuJKy8xjvB/Ia7LAmhrGQNFzS3p33gbPJ739QrjTnnQd098twYPI8Y1129vpc4HGSM/iXp8uuIEl0kHz4f0ByUus+4NAcy74DeBZYkt7m5VFu2bb3UIHeOxlfr4D/QzJcxlJgWo7/6/HAb0l6bywBzq5QuTcCG4DXSGqVFwJ/CfxlyWu+Oo1raaX+173E9P40lm3p8XV7ybrL01hWAOdUO5aScv8rff0PkySdA3Mse5djI+8bSe+hh9Lbsrzi6Ob4fDNwJ7Ay/btfNcr2FblmZg3EV+SamTUQJ30zswbipG9m1kCc9M3MGoiTvplZA3HSNzNrIE76ZmYNxEnfzKyB/H/GxpPhCZ8T7QAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "h=norm(.807,.657)\n", + "fig,ax=plt.subplots(2,2)\n", + "ax[0,0].set_xlim(-5,5)\n", + "j=ax[0,0].hist(z,bins=100,density=True)\n", + "j=ax[0,1].hist(samples_x,bins=50,density=True)\n", + "j=ax[1,1].hist(samples_y,bins=50,density=True)\n", + "ax[0,0].plot(np.linspace(-5,5,100),h.pdf(np.linspace(-5,5,100)))\n", + "\n", + "b=plt.hist(z,bins=100)" + ] + }, + { + "cell_type": "code", + "execution_count": 58, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(array([1.000e+00, 0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00,\n", + " 0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00,\n", + " 0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00,\n", + " 0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00,\n", + " 0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00,\n", + " 0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00,\n", + " 0.000e+00, 1.000e+00, 1.000e+00, 6.000e+00, 7.000e+00, 1.800e+01,\n", + " 5.800e+01, 1.330e+02, 3.230e+02, 6.530e+02, 1.059e+03, 1.401e+03,\n", + " 1.615e+03, 1.410e+03, 1.124e+03, 7.740e+02, 5.150e+02, 3.100e+02,\n", + " 2.150e+02, 1.320e+02, 7.000e+01, 5.800e+01, 2.700e+01, 2.600e+01,\n", + " 2.200e+01, 7.000e+00, 5.000e+00, 7.000e+00, 4.000e+00, 2.000e+00,\n", + " 3.000e+00, 2.000e+00, 1.000e+00, 0.000e+00, 1.000e+00, 1.000e+00,\n", + " 0.000e+00, 0.000e+00, 2.000e+00, 0.000e+00, 0.000e+00, 1.000e+00,\n", + " 0.000e+00, 1.000e+00, 1.000e+00, 0.000e+00, 0.000e+00, 0.000e+00,\n", + " 0.000e+00, 0.000e+00, 1.000e+00, 0.000e+00, 1.000e+00, 0.000e+00,\n", + " 0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00,\n", + " 0.000e+00, 0.000e+00, 0.000e+00, 1.000e+00]),\n", + " array([-9.34383882, -9.13726736, -8.9306959 , -8.72412443, -8.51755297,\n", + " -8.31098151, -8.10441005, -7.89783858, -7.69126712, -7.48469566,\n", + " -7.2781242 , -7.07155274, -6.86498127, -6.65840981, -6.45183835,\n", + " -6.24526689, -6.03869543, -5.83212396, -5.6255525 , -5.41898104,\n", + " -5.21240958, -5.00583812, -4.79926665, -4.59269519, -4.38612373,\n", + " -4.17955227, -3.97298081, -3.76640934, -3.55983788, -3.35326642,\n", + " -3.14669496, -2.9401235 , -2.73355203, -2.52698057, -2.32040911,\n", + " -2.11383765, -1.90726619, -1.70069472, -1.49412326, -1.2875518 ,\n", + " -1.08098034, -0.87440888, -0.66783741, -0.46126595, -0.25469449,\n", + " -0.04812303, 0.15844844, 0.3650199 , 0.57159136, 0.77816282,\n", + " 0.98473428, 1.19130575, 1.39787721, 1.60444867, 1.81102013,\n", + " 2.01759159, 2.22416306, 2.43073452, 2.63730598, 2.84387744,\n", + " 3.0504489 , 3.25702037, 3.46359183, 3.67016329, 3.87673475,\n", + " 4.08330621, 4.28987768, 4.49644914, 4.7030206 , 4.90959206,\n", + " 5.11616352, 5.32273499, 5.52930645, 5.73587791, 5.94244937,\n", + " 6.14902083, 6.3555923 , 6.56216376, 6.76873522, 6.97530668,\n", + " 7.18187815, 7.38844961, 7.59502107, 7.80159253, 8.00816399,\n", + " 8.21473546, 8.42130692, 8.62787838, 8.83444984, 9.0410213 ,\n", + " 9.24759277, 9.45416423, 9.66073569, 9.86730715, 10.07387861,\n", + " 10.28045008, 10.48702154, 10.693593 , 10.90016446, 11.10673592,\n", + " 11.31330739]),\n", + " )" + ] + }, + "execution_count": 58, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "b\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 +} From aa114e22fee44b01a869476c89774f97e12ca604 Mon Sep 17 00:00:00 2001 From: Jeremy Teitelbaum Date: Fri, 27 Apr 2018 09:33:26 -0400 Subject: [PATCH 2/8] Added a useful formula page --- Useful Formulae.ipynb | 61 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 61 insertions(+) create mode 100644 Useful Formulae.ipynb diff --git a/Useful Formulae.ipynb b/Useful Formulae.ipynb new file mode 100644 index 0000000..fd0c133 --- /dev/null +++ b/Useful Formulae.ipynb @@ -0,0 +1,61 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Conjugate Normal Distributions (known variance)\n", + "\n", + "We are trying to learn about the unknown mean of a normal distribution with known variance. \n", + "We choose a prior distribution is normal with mean $\\mu_{0}$ and variance $\\tau_{0}^2$. \n", + "We draw $n$ values $y_1,\\ldots, y_n$ from the distribution with known variance $\\sigma^2$. The posterior distribution\n", + "$p(\\mu|y_1,\\ldots,y_n)=p(y_1,\\ldots,y_n|\\mu)p(\\mu)$ is again normal. Let \n", + "$$\n", + "\\overline{y}=\\frac{1}{n}\\sum_{i=1}^{n} y_i\n", + "$$\n", + "be the sample mean. \n", + "\n", + "The posterior variance\n", + "is\n", + "$$\\frac{1}{\\tau_1^2}=\\frac{1}{\\tau_0^2}+\\frac{n}{\\sigma^2}$$\n", + "and the posterior mean is\n", + "$$\n", + "\\mu_1=\\frac{\\frac{\\mu_0}{\\tau_0^2}+\\frac{n\\overline{y}}{\\sigma^2}}{\\frac{1}{\\tau_{1}^2}}\n", + "$$\n" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [], + "source": [ + "def posterior(prior_mean,prior_variance,sample_mean,pop_variance,n):\n", + " post_var=1/((1/prior_variance) + n/pop_variance)\n", + " post_mean=(prior_mean/prior_variance+sample_mean*n/pop_variance)/(1/post_var)\n", + " return post_mean, post_var" + ] + } + ], + "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 1715a6d08308697b3c286fc2368abd85d934f85d Mon Sep 17 00:00:00 2001 From: Jeremy Teitelbaum Date: Fri, 27 Apr 2018 12:01:51 -0400 Subject: [PATCH 3/8] updated useful formulae --- Useful Formulae.ipynb | 89 +++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 86 insertions(+), 3 deletions(-) diff --git a/Useful Formulae.ipynb b/Useful Formulae.ipynb index fd0c133..3fc65f4 100644 --- a/Useful Formulae.ipynb +++ b/Useful Formulae.ipynb @@ -21,20 +21,103 @@ "and the posterior mean is\n", "$$\n", "\\mu_1=\\frac{\\frac{\\mu_0}{\\tau_0^2}+\\frac{n\\overline{y}}{\\sigma^2}}{\\frac{1}{\\tau_{1}^2}}\n", - "$$\n" + "$$\n", + "\n", + "The posterior sampling distribution $\\theta$ is\n", + "$$\n", + "p( z |y)=\\int_{\\theta} p(z|\\theta) d\\theta\n", + "$$\n", + "is a normal distribution with mean equal to the posterior mean $\\mu_1$ and variance equal to $\\sigma^2+\\tau_1^2$\n", + "where $\\tau_1$ is the posterior variance.\n", + "\n", + "See Pages 39-42 of BDA (Section 2.5) for more information." ] }, { "cell_type": "code", - "execution_count": 8, + "execution_count": 24, "metadata": {}, "outputs": [], "source": [ + "from scipy.stats import norm\n", + "import numpy as np\n", "def posterior(prior_mean,prior_variance,sample_mean,pop_variance,n):\n", " post_var=1/((1/prior_variance) + n/pop_variance)\n", " post_mean=(prior_mean/prior_variance+sample_mean*n/pop_variance)/(1/post_var)\n", - " return post_mean, post_var" + " return post_mean, post_var\n", + "\n", + "def post_sample(y,prior_mean,prior_variance,sample_mean,pop_variance,n):\n", + " post_mean,post_var=posterior(prior_mean,prior_variance,sample_mean,pop_variance,n)\n", + " return norm.pdf(y,post_mean,np.sqrt(pop_variance+post_var))\n", + " " ] + }, + { + "cell_type": "code", + "execution_count": 26, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "0.7403867575800461" + ] + }, + "execution_count": 26, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "post_sample(-.25,1,.25,-.25,1,10)+post_sample(-.25,-1,.25,-.25,1,10)" + ] + }, + { + "cell_type": "code", + "execution_count": 29, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "0.05095226579074726" + ] + }, + "execution_count": 29, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + ".1*post_sample(-.25,-1,.25,-.25,1,10)/(post_sample(-.25,1,.25,-.25,1,10)+post_sample(-.25,-1,.25,-.25,1,10))" + ] + }, + { + "cell_type": "code", + "execution_count": 30, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "0.4414296078832747" + ] + }, + "execution_count": 30, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + ".9*post_sample(-.25,1,.25,-.25,1,10)/(post_sample(-.25,1,.25,-.25,1,10)+post_sample(-.25,-1,.25,-.25,1,10))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { From 8e76a6b3b1770c250e246d97bcfbce6de7f51bae Mon Sep 17 00:00:00 2001 From: Jeremy Teitelbaum Date: Fri, 27 Apr 2018 14:42:03 -0400 Subject: [PATCH 4/8] 5.9.8 added (in progress) --- BDA 5.9.8.ipynb | 72 +++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 72 insertions(+) create mode 100644 BDA 5.9.8.ipynb diff --git a/BDA 5.9.8.ipynb b/BDA 5.9.8.ipynb new file mode 100644 index 0000000..9a60fc0 --- /dev/null +++ b/BDA 5.9.8.ipynb @@ -0,0 +1,72 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Discrete Mixture Models\n", + "\n", + "Discrete mixture models: if $p_m(\\theta)$, for $m=1,\\ldots,M$ are conjugate prior densities for the sampling model $y|\\theta$, show that the class of finite mixture prior densities given by \n", + "$$\n", + "p(\\theta)=\\sum_{1}^{M} \\lambda_m p_m(\\theta)\n", + "$$\n", + "is also a conjugate class, where the $\\lambda_m$’s are nonnegative weights that sum to 1. This can provide a useful extension of the natural conjugate prior family to more flexible distributional forms. As an example, use the mixture form to create a bimodal prior density for a normal mean, that is thought to be near $1$, with a standard deviation of $0.5$, but has a small probability of being near $−1$, with the same standard deviation. If the variance of each observation $y_1,\\ldots,y_{10}$ is known to be $1$, and their observed mean is $y =−0.25$, derive your posterior distribution for the mean, making a sketch of both prior and posterior densities. Be careful: the prior and posterior mixture proportions are different.\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let's skip the theory part and look at the example.\n", + "\n", + "We have\n", + "$$\n", + "p(\\theta|y_1,\\ldots,y_{10})\\propto p(y_1,\\ldots,y_10|\\theta)p(\\theta)$$\n", + "so\n", + "$$\n", + "p(\\theta|\\{y_{i}\\})\\propto \\sum \\lambda_{m}p(\\{y_{i}\\}|\\theta)p_{m}(\\theta)\n", + "$$\n", + "\n", + "Each of the terms $p_{m}(\\theta)p(\\{y_{i}\\}|\\theta)$\n", + "is equal to $p_{m}(\\theta|\\{y_{i}\\})p_{m}(\\{y_{i}\\})$.\n", + "\n", + "Therefore the total posterior density is a weighted sum\n", + "of the individual posteriors:\n", + "\n", + "$$p(\\theta|\\{y_{i}\\})=\\sum c_{m}p_{m}(\\theta|\\{y_{i}\\})$$\n", + "where \n", + "$$\n", + "c_{m}=\\frac{\\lambda_m p_{m}(\\{y_{i}\\})}{\\sum_{m} \\lambda_m p_{m}(\\{y_{i}\\}}\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 +} From 681a084f8f5be0932cfaface8f589aec74dd7f9c Mon Sep 17 00:00:00 2001 From: Jeremy Teitelbaum Date: Sun, 29 Apr 2018 09:15:40 -0400 Subject: [PATCH 5/8] still trying to figure out BDA 5.9.8 --- BDA 5.9.8.ipynb | 74 ++++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 73 insertions(+), 1 deletion(-) diff --git a/BDA 5.9.8.ipynb b/BDA 5.9.8.ipynb index 9a60fc0..d3ea631 100644 --- a/BDA 5.9.8.ipynb +++ b/BDA 5.9.8.ipynb @@ -37,7 +37,79 @@ "where \n", "$$\n", "c_{m}=\\frac{\\lambda_m p_{m}(\\{y_{i}\\})}{\\sum_{m} \\lambda_m p_{m}(\\{y_{i}\\}}\n", - "$$" + "$$\n", + "\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In the special case under consideration, $p_1$ is normal with mean $-1$ and $\\sigma=.5$, $p_2$ is normal with mean $1$ and $\\sigma=.5$ and we can set $\\lambda_1=.1$ and $\\lambda_2=.9$. The $p_m(\\{y_{i}\\})$ can be calculated from the $t$ distribution. Drawing a sample of size $10$ from $p_1$ and getting a sample mean of $-.25$ and a sample variance of $1$ gives a $t$-statistics of $\\sqrt{10}(-.25+1)$ in the first case and $\\sqrt{10}(-.25-1)$ in the second. " + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "0.041664931082753924\n", + "0.0035119750957915393\n" + ] + } + ], + "source": [ + "import numpy as np\n", + "from scipy.stats import norm, t\n", + "t_1=np.sqrt(9)*.75\n", + "t_2=np.sqrt(9)*1.25\n", + "print(t.pdf(t_1,df=9))\n", + "print(t.pdf(t_2,df=9))" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "0.5655172413793104" + ] + }, + "execution_count": 18, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + ".1*.041/(.1*.041+.9*.0035)" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "0.43448275862068964" + ] + }, + "execution_count": 17, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + ".9*.0035/(.1*.041+.9*.0035)" ] }, { From 9c03e4d65c419d81cbac1dec49539b2dfdf21ff6 Mon Sep 17 00:00:00 2001 From: Jeremy Teitelbaum Date: Sun, 29 Apr 2018 13:00:00 -0400 Subject: [PATCH 6/8] reconciling home and web --- BDA 5.9.3.ipynb | 96 ++++++++++++++++++++++++------------------------- 1 file changed, 48 insertions(+), 48 deletions(-) diff --git a/BDA 5.9.3.ipynb b/BDA 5.9.3.ipynb index abfbe93..850878d 100644 --- a/BDA 5.9.3.ipynb +++ b/BDA 5.9.3.ipynb @@ -128,16 +128,16 @@ }, { "cell_type": "code", - "execution_count": 350, + "execution_count": 355, "metadata": {}, "outputs": [], "source": [ - "answers=sm.sampling(data=dict({'means':p55['effect'],'se':p55['se']}),iter=00)" + "answers=sm.sampling(data=dict({'means':p55['effect'],'se':p55['se']}),iter=5000)" ] }, { "cell_type": "code", - "execution_count": 351, + "execution_count": 356, "metadata": {}, "outputs": [ { @@ -145,31 +145,31 @@ "output_type": "stream", "text": [ "Inference for Stan model: anon_model_1c0a010b4129370aa04f0b4b9f729b4d.\n", - "4 chains, each with iter=500; warmup=250; thin=1; \n", - "post-warmup draws per chain=250, total post-warmup draws=1000.\n", + "4 chains, each with iter=5000; warmup=2500; thin=1; \n", + "post-warmup draws per chain=2500, total post-warmup draws=10000.\n", "\n", " mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat\n", - "theta[0] 12.07 0.66 9.38 -3.89 5.93 11.21 17.23 34.12 200 1.01\n", - "theta[1] 7.54 0.39 6.89 -6.07 3.02 7.48 12.15 21.17 312 1.01\n", - "theta[2] 5.45 0.46 8.58 -12.34 0.2 6.07 10.92 20.96 342 1.0\n", - "theta[3] 7.16 0.37 7.21 -7.68 2.72 7.31 11.83 20.44 383 1.01\n", - "theta[4] 4.51 0.43 6.58 -10.0 0.23 5.07 9.18 15.99 235 1.02\n", - "theta[5] 5.53 0.41 7.27 -10.59 1.01 6.16 10.74 18.24 313 1.01\n", - "theta[6] 11.03 0.53 7.5 -2.83 5.64 11.03 15.77 27.15 199 1.01\n", - "theta[7] 8.3 0.39 7.93 -6.88 3.08 8.28 13.26 24.5 414 1.0\n", - "mu 7.66 0.41 5.62 -3.36 4.06 7.92 11.47 18.73 192 1.02\n", - "tau 7.74 0.71 6.11 0.4 3.91 6.29 9.9 22.49 75 1.05\n", - "results[0] 11.81 0.57 12.82 -9.18 3.71 10.54 17.23 43.38 513 1.01\n", - "results[1] 7.26 0.52 12.33 -16.66 0.94 7.04 13.57 33.52 563 1.01\n", - "results[2] 5.43 0.42 12.64 -22.88 -1.07 5.99 12.55 29.09 902 1.0\n", - "results[3] 7.01 0.41 11.92 -16.05 1.03 7.18 13.39 29.57 861 1.0\n", - "results[4] 4.62 0.46 11.55 -23.93 -1.48 5.45 11.65 26.17 638 1.0\n", - "results[5] 5.83 0.44 11.71 -20.83 0.05 6.61 11.82 28.91 724 1.0\n", - "results[6] 10.76 0.56 11.87 -11.69 3.67 10.23 17.43 36.29 443 1.0\n", - "results[7] 8.01 0.54 12.75 -16.27 1.62 7.82 14.62 31.63 561 1.0\n", - "lp__ -18.65 1.21 5.8 -27.82 -22.04 -19.13 -16.31 0.06 23 1.18\n", + "theta[0] 11.66 0.17 8.25 -2.16 6.27 10.51 16.07 30.71 2223 1.0\n", + "theta[1] 8.02 0.1 6.27 -4.44 4.18 7.82 11.76 21.02 3674 1.0\n", + "theta[2] 6.36 0.12 7.68 -10.59 2.17 6.78 11.15 20.84 3806 1.0\n", + "theta[3] 7.79 0.11 6.43 -5.29 4.02 7.62 11.79 20.79 3689 1.0\n", + "theta[4] 5.19 0.13 6.38 -8.66 1.46 5.72 9.44 16.84 2357 1.0\n", + "theta[5] 6.3 0.12 6.62 -7.68 2.38 6.54 10.57 18.98 3277 1.0\n", + "theta[6] 10.99 0.15 6.78 -0.97 6.41 10.4 14.99 25.86 2181 1.0\n", + "theta[7] 8.61 0.13 7.97 -7.29 4.13 8.21 12.91 25.66 3946 1.0\n", + "mu 8.2 0.11 5.3 -1.85 5.01 7.96 11.32 18.88 2354 1.0\n", + "tau 6.84 0.2 5.58 0.93 2.99 5.46 9.11 20.69 751 1.01\n", + "results[0] 11.6 0.18 12.19 -8.79 4.97 9.96 16.83 40.26 4692 1.0\n", + "results[1] 7.97 0.13 10.69 -14.26 2.73 7.74 13.36 29.84 6547 1.0\n", + "results[2] 6.4 0.15 11.84 -18.96 1.18 6.95 12.45 28.63 5970 1.0\n", + "results[3] 7.79 0.14 11.1 -14.73 2.66 7.66 13.29 30.0 6438 1.0\n", + "results[4] 5.24 0.16 11.09 -19.44 0.35 6.13 11.03 24.67 4830 1.0\n", + "results[5] 6.36 0.15 10.98 -17.1 1.36 6.66 11.93 28.06 5695 1.0\n", + "results[6] 10.93 0.16 11.19 -8.59 4.94 9.77 16.05 35.59 4800 1.0\n", + "results[7] 8.64 0.14 11.63 -13.75 2.97 8.25 14.03 33.19 6700 1.0\n", + "lp__ -17.68 0.41 5.52 -27.62 -21.48 -18.16 -14.18 -5.46 179 1.04\n", "\n", - "Samples were drawn using NUTS at Sat Apr 21 16:45:21 2018.\n", + "Samples were drawn using NUTS at Sat Apr 21 16:45:42 2018.\n", "For each parameter, n_eff is a crude measure of effective sample size,\n", "and Rhat is the potential scale reduction factor on split chains (at \n", "convergence, Rhat=1).\n" @@ -182,12 +182,12 @@ }, { "cell_type": "code", - "execution_count": 352, + "execution_count": 357, "metadata": {}, "outputs": [ { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAEICAYAAACzliQjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAGC1JREFUeJzt3XuUHGWdxvHvswkEBQyQhBWSyAQTkbAiSgyoiEoUgyhBN0gianRZcc/Cqkc9mngBxBt4dkFW2UuOoIiXgPGWhWhkQV1FjRkuAmOIDDGQIQoTEsCIXBJ++0e9o0WnZ6Z6pjPTM+/zOWfOdL31VvWvq2uefuftnhpFBGZmloe/Ge4CzMxs6Dj0zcwy4tA3M8uIQ9/MLCMOfTOzjDj0zcwy4tC3ppK0TdLBQ3h/L5O0ron7+76kRen22yX9rIn7Pk3SD5u1vwbu96WS7kzPzclDff/WWhz6I5SkDZL+nH6Q75P0JUl7DWJ/bZJC0tjB1BURe0XE+sHso1TTuZKekPTH9PVbSV+QdEDp/n4aEYdU3NdX++sXESdExOVNqH2n4xkRX4uI4we77wE4D/hCem6+W7uy2edSbyTNk3SLpIclbZZ0naS2tK7f59qaw6E/sr0+IvYCXgi8CPjocBUy2BeLPra/MiL2BvYD3gA8E7ix2WGgwmj9eTgI6Oinzy49lyRNB74CvB8YD0wD/gN4stRtSJ7r3I3WkzwrEXEv8H3g7wAkHShphaQtkjolvbOnr6TZktrTaOs+SRemVf+Xvj+YRnwvTv3/QdJaSVslrZJ0UGlfIelMSXcCd5bapqfb4yV9RVK3pLslfbQnWNPUyQ2SLpK0BTi3n8f4RER0AKcC3RThgaRXSOoq1fQhSfem0eI6SXMkzQU+DJyaHtuvU98fS/qUpBuAR4CDU9s/lu5akj4v6SFJd0iaU1qxQdKrSsvl3yZ2Op6100WSXiJpTdr3GkkvKa37saRPpGP0R0k/lDSxt+Mj6Z3pud6SnvsDU/tdwMHA/6Q6xvVznBs5l86VtFzSlanGmyQ9v5ddHwH8LiKui8IfI+JbEXFPnRrqPtfWHA79UUDSVOC1wM2p6RtAF3AgMB/4dCmsLgYujohnAM8Grkrtx6bv+6RpgF+omP/9MPBGYBLw07TvspOBo4CZdUr7PMWo7mDg5cDbgHeU1h8FrAf2Bz5V5bFGxA7ge8DLatdJOgQ4C3hRGjG+BtgQET8APk0xktwrIsrB9FbgDGBv4O46d9lT40TgHODbkvarUOpOx7Om1v2Aa4B/ByYAFwLXSJpQ6vZmiuO1P7A78IF6dyTpOOAzwJuAA9LjWAYQEc8G7iGN5CPisb6KbvBcApgHfJNidP514LuSdquz65uA56YX+VeqwvRRX8+1DZxDf2T7rqQHgZ8BP6H4gZwKHAN8KCIejYhbgC9ShBvAE8B0SRMjYltE/LKP/b8L+ExErI2I7RTBeUR5tJ/Wb4mIP5c3lDSGYqS2JI3qNgD/VqoDYFNEfD4ittdu349NFCFTawcwDpgpabeI2BARd/Wzry9HREeq4Yk66+8HPpdGn1cC64ATG6i1NycCd0bEFem+vwHcAby+1OdLEfHbdGyuohgt13MacFlE3JRCfQnwYqX58ooGci4B3BgRy9OxuxDYAzi6dufpfZ5XAJPTY9ks6csVwr+359oGyKE/sp0cEftExEER8c8pHA4EtkTEH0v97qb4YQM4HXgOcEeaUnhdH/s/CLhY0oMpELYAKu0LYGMv206kGJ2WR8/lOvratj+TUy1PERGdwHspporul7SsZ5qjD/3VcG889aqEd1Mc48E6kJ1/s6g9Pn8o3X4E6C0gn7KviNgGPFCzr/4M5FyC0vGLiCf5628FO4mIX0bEmyJiEsXo/VjgI/3UVfe5toFz6I8+m4D9JO1dansWcC9ARNwZEQsppgwuAJZL2hOod7nVjcC7Uhj0fD0tIn5e6tPbZVo3U/xWUf6t4C919LNtr9J7Aq+nmGraSUR8PSKOSfcbFI+xr/vqr4bJklRafhbFMQb4E/D00rpnNrDfTTz12PTs+946ffvzlH2l53PCAPdVu99ez6Vkaul+/waYwl+PT68iYg3wbdJ7B/X091zbwDj0R5mI2Aj8HPiMpD0kHU4xuv8agKS3SJqURmUPps12ULxh9iTF/HuP/wKWSDosbTte0ikV69hB8Wv8pyTtnaaE3gf0+7HJeiTtJulQijnmZ1JMJdT2OUTScenNykeBP6fHBnAf0KbGP6GzP/DudP+nAIcCK9O6W4AFad0sijnvHvWOZ9lK4DmS3ixprKRTKd4XubrB+qCYS3+HpCPSY/80sDpNqQ1Yf+dScqSkN6r49NV7gceAnaYMJR2T3mzePy0/Fzipl779Ptc2cA790Wkh0EYx4voOcE5EXJvWzQU6JG2jeFN3QZqvfYTizdQb0nTO0RHxHYqR8jJJDwO3Ayc0UMe/UIyG11PMFX8duKzBx3JqqvVBYAXFtMWREVFvNDkOOJ/it4w/UAT2h9O6b6bvD0i6qYH7Xw3MSPv8FDA/Ih5I6z5G8Wb4VuDjFI8PgHrHs7zTtI/XUXwy5QHgg8DrImJzA7X17Ou6VMu3gN+nmhY0up9e9HUuQfFG66kUx+CtwBt7eW/kQYqQvy09nz9I+/tsqU8jz7UNkPxPVMxsICSdC0yPiLcMdy1WnUf6ZmYZceibmWXE0ztmZhnxSN/MLCODukjWrjBx4sRoa2sb7jLMzEaUG2+8cXP6w7c+tVzot7W10d7ePtxlmJmNKJLqXTtqJ57eMTPLiEPfzCwjDn0zs4w49M3MMuLQNzPLiEPfzCwjDn0zs4w49M3MMuLQNzPLSMv9Re5Qa1t8Td32Dec3439fm5m1Fo/0zcwy4tA3M8uIQ9/MLCMOfTOzjDj0zcwy4tA3M8uIQ9/MLCMOfTOzjDj0zcwy4tA3M8uIQ9/MLCMOfTOzjDj0zcwy4tA3M8tIpdCXNFfSOkmdkhbXWX+spJskbZc0v2bdIkl3pq9FzSrczMwa12/oSxoDXAKcAMwEFkqaWdPtHuDtwNdrtt0POAc4CpgNnCNp38GXbWZmA1FlpD8b6IyI9RHxOLAMmFfuEBEbIuJW4MmabV8DXBsRWyJiK3AtMLcJdZuZ2QBUCf3JwMbScldqq6LStpLOkNQuqb27u7virs3MrFFV/l2i6rRFxf1X2jYilgJLAWbNmlV137tUb/9GsTf+94pmNhJUGel3AVNLy1OATRX3P5htzcysyaqE/hpghqRpknYHFgArKu5/FXC8pH3TG7jHpzYzMxsG/YZ+RGwHzqII67XAVRHRIek8SScBSHqRpC7gFOC/JXWkbbcAn6B44VgDnJfazMxsGFSZ0yciVgIra9rOLt1eQzF1U2/by4DLBlGjmZk1if8i18wsIw59M7OMOPTNzDLi0Dczy4hD38wsIw59M7OMOPTNzDLi0Dczy4hD38wsIw59M7OMOPTNzDLi0Dczy4hD38wsIw59M7OMOPTNzDLi0Dczy4hD38wsIw59M7OMOPTNzDLi0Dczy4hD38wsIw59M7OMOPTNzDLi0Dczy4hD38wsIw59M7OMOPTNzDLi0Dczy0il0Jc0V9I6SZ2SFtdZP07SlWn9akltqX03SZdLuk3SWklLmlu+mZk1ot/QlzQGuAQ4AZgJLJQ0s6bb6cDWiJgOXARckNpPAcZFxPOAI4F39bwgmJnZ0Ksy0p8NdEbE+oh4HFgGzKvpMw+4PN1eDsyRJCCAPSWNBZ4GPA483JTKzcysYVVCfzKwsbTcldrq9omI7cBDwASKF4A/Ab8H7gH+NSK2DLJmMzMboCqhrzptUbHPbGAHcCAwDXi/pIN3ugPpDEntktq7u7srlGRmZgNRJfS7gKml5SnApt76pKmc8cAW4M3ADyLiiYi4H7gBmFV7BxGxNCJmRcSsSZMmNf4ozMyskiqhvwaYIWmapN2BBcCKmj4rgEXp9nzg+ogIiimd41TYEzgauKM5pZuZWaP6Df00R38WsApYC1wVER2SzpN0Uup2KTBBUifwPqDnY52XAHsBt1O8eHwpIm5t8mMwM7OKxlbpFBErgZU1bWeXbj9K8fHM2u221Ws3M7Ph4b/INTPLiEPfzCwjDn0zs4w49M3MMlLpjdzRoG3xNcNdgpnZsPNI38wsIw59M7OMOPTNzDLi0Dczy4hD38wsIw59M7OMOPTNzDLi0Dczy4hD38wsIw59M7OMOPTNzDLi0Dczy4hD38wsIw59M7OMOPTNzDLi0Dczy4hD38wsIw59M7OMZPPvEne13v4d44bzTxziSszMeueRvplZRhz6ZmYZceibmWXEoW9mlhGHvplZRiqFvqS5ktZJ6pS0uM76cZKuTOtXS2orrTtc0i8kdUi6TdIezSvfzMwa0W/oSxoDXAKcAMwEFkqaWdPtdGBrREwHLgIuSNuOBb4K/FNEHAa8AniiadWbmVlDqoz0ZwOdEbE+Ih4HlgHzavrMAy5Pt5cDcyQJOB64NSJ+DRARD0TEjuaUbmZmjaoS+pOBjaXlrtRWt09EbAceAiYAzwFC0ipJN0n6YL07kHSGpHZJ7d3d3Y0+BjMzq6hK6KtOW1TsMxY4BjgtfX+DpDk7dYxYGhGzImLWpEmTKpRkZmYDUSX0u4CppeUpwKbe+qR5/PHAltT+k4jYHBGPACuBFw62aDMzG5gqob8GmCFpmqTdgQXAipo+K4BF6fZ84PqICGAVcLikp6cXg5cDv2lO6WZm1qh+L7gWEdslnUUR4GOAyyKiQ9J5QHtErAAuBa6Q1Ekxwl+Qtt0q6UKKF44AVkZE/SuTmZnZLlfpKpsRsZJiaqbcdnbp9qPAKb1s+1WKj22amdkw81/kmpllxKFvZpYRh76ZWUYc+mZmGfG/S9zF/G8UzayVeKRvZpYRh76ZWUYc+mZmGXHom5llxKFvZpYRh76ZWUYc+mZmGXHom5llxKFvZpYRh76ZWUYc+mZmGXHom5llxKFvZpYRh76ZWUYc+mZmGXHom5llxKFvZpYRh76ZWUYc+mZmGXHom5llxKFvZpYRh76ZWUbGDncBuWpbfE3d9g3nnzjElZhZTjzSNzPLSKXQlzRX0jpJnZIW11k/TtKVaf1qSW01658laZukDzSnbDMzG4h+Q1/SGOAS4ARgJrBQ0syabqcDWyNiOnARcEHN+ouA7w++XDMzG4wqI/3ZQGdErI+Ix4FlwLyaPvOAy9Pt5cAcSQKQdDKwHuhoTslmZjZQVUJ/MrCxtNyV2ur2iYjtwEPABEl7Ah8CPt7XHUg6Q1K7pPbu7u6qtZuZWYOqhL7qtEXFPh8HLoqIbX3dQUQsjYhZETFr0qRJFUoyM7OBqPKRzS5gaml5CrCplz5dksYC44EtwFHAfEmfBfYBnpT0aER8YdCVm5lZw6qE/hpghqRpwL3AAuDNNX1WAIuAXwDzgesjIoCX9XSQdC6wzYFvZjZ8+g39iNgu6SxgFTAGuCwiOiSdB7RHxArgUuAKSZ0UI/wFu7JoMzMbmEp/kRsRK4GVNW1nl24/CpzSzz7OHUB9ZmbWRP6LXDOzjIy6a+/0dk0bMzPzSN/MLCsOfTOzjDj0zcwyMurm9Ec6X2ffzHYlj/TNzDLi0Dczy4hD38wsIw59M7OMOPTNzDLi0Dczy4hD38wsIw59M7OMOPTNzDLi0Dczy4hD38wsIw59M7OMOPTNzDLi0Dczy4hD38wsIw59M7OMOPTNzDLi0Dczy4hD38wsIw59M7OMOPTNzDLi0Dczy4hD38wsI5VCX9JcSeskdUpaXGf9OElXpvWrJbWl9ldLulHSben7cc0t38zMGtFv6EsaA1wCnADMBBZKmlnT7XRga0RMBy4CLkjtm4HXR8TzgEXAFc0q3MzMGldlpD8b6IyI9RHxOLAMmFfTZx5webq9HJgjSRFxc0RsSu0dwB6SxjWjcDMza1yV0J8MbCwtd6W2un0iYjvwEDChps/fAzdHxGO1dyDpDEntktq7u7ur1m5mZg0aW6GP6rRFI30kHUYx5XN8vTuIiKXAUoBZs2bV7tuAtsXX1G3fcP6JQ1yJmY1kVUb6XcDU0vIUYFNvfSSNBcYDW9LyFOA7wNsi4q7BFmxmZgNXJfTXADMkTZO0O7AAWFHTZwXFG7UA84HrIyIk7QNcAyyJiBuaVbSZmQ1Mv6Gf5ujPAlYBa4GrIqJD0nmSTkrdLgUmSOoE3gf0fKzzLGA68DFJt6Sv/Zv+KMzMrJIqc/pExEpgZU3b2aXbjwKn1Nnuk8AnB1mjmZk1if8i18wsIw59M7OMVJresdbV20c5wR/nNLOdeaRvZpYRh76ZWUYc+mZmGXHom5llxKFvZpYRh76ZWUb8kU37C1/J02z0c+hnqK/P9pvZ6ObQH8Uc7mZWy3P6ZmYZceibmWXEoW9mlhGHvplZRhz6ZmYZceibmWXEH9m0AfMfc5mNPB7pm5llxKFvZpYRT+/YkPF0kNnw80jfzCwjDn0zs4x4esf65Qu3mY0eHumbmWXEoW9mlhGHvplZRjynb03X6HsA/iin2dCpFPqS5gIXA2OAL0bE+TXrxwFfAY4EHgBOjYgNad0S4HRgB/DuiFjVtOrNKujrRWikv7D4BdMa1W/oSxoDXAK8GugC1khaERG/KXU7HdgaEdMlLQAuAE6VNBNYABwGHAj8r6TnRMSOZj8Qy0czP03UrNB0+A5cjsduOB9zlTn92UBnRKyPiMeBZcC8mj7zgMvT7eXAHElK7csi4rGI+B3QmfZnZmbDoMr0zmRgY2m5Cziqtz4RsV3SQ8CE1P7Lmm0n196BpDOAM9LiNknrKlW/s4nA5gFuO9RGUq0wDPXqggFvOuhaB3HfA9lP049ts+rvxZCcC016DCPq50wXDKreg6p0qhL6qtMWFftU2ZaIWAosrVBLnyS1R8Sswe5nKIykWmFk1TuSagXXuyuNpFphaOqtMr3TBUwtLU8BNvXWR9JYYDywpeK2ZmY2RKqE/hpghqRpknaneGN2RU2fFcCidHs+cH1ERGpfIGmcpGnADOBXzSndzMwa1e/0TpqjPwtYRfGRzcsiokPSeUB7RKwALgWukNRJMcJfkLbtkHQV8BtgO3DmLv7kzqCniIbQSKoVRla9I6lWcL270kiqFYagXhUDcjMzy4Evw2BmlhGHvplZRkZF6EuaK2mdpE5Ji4e7nlqSLpN0v6TbS237SbpW0p3p+77DWWMPSVMl/UjSWkkdkt6T2lu13j0k/UrSr1O9H0/t0yStTvVemT6E0BIkjZF0s6Sr03Ir17pB0m2SbpHUntpa8lwAkLSPpOWS7kjn8ItbsV5Jh6Rj2vP1sKT3DkWtIz70S5eJOAGYCSxMl39oJV8G5ta0LQaui4gZwHVpuRVsB94fEYcCRwNnpuPZqvU+BhwXEc8HjgDmSjqa4lIgF6V6t1JcKqRVvAdYW1pu5VoBXhkRR5Q+P96q5wIU1wj7QUQ8F3g+xXFuuXojYl06pkdQXLPsEeA7DEWtETGiv4AXA6tKy0uAJcNdV50624DbS8vrgAPS7QOAdcNdYy91f4/iukstXy/wdOAmir8Y3wyMrXeODHONU9IP83HA1RR/wNiStaZ6NgATa9pa8lwAngH8jvQBlVavt1Tf8cANQ1XriB/pU/8yETtd6qEF/W1E/B4gfd9/mOvZiaQ24AXAalq43jRdcgtwP3AtcBfwYERsT11a6Zz4HPBB4Mm0PIHWrRWKv6D/oaQb0+VSoHXPhYOBbuBLafrsi5L2pHXr7bEA+Ea6vctrHQ2hX+lSD9YYSXsB3wLeGxEPD3c9fYmIHVH8mjyF4oJ+h9brNrRV7UzS64D7I+LGcnOdrsNea8lLI+KFFNOnZ0o6drgL6sNY4IXAf0bEC4A/0QJTOX1J79+cBHxzqO5zNIT+SL3Uw32SDgBI3+8f5nr+QtJuFIH/tYj4dmpu2Xp7RMSDwI8p3ovYJ10SBFrnnHgpcJKkDRRXqz2OYuTfirUCEBGb0vf7KeacZ9O650IX0BURq9PycooXgVatF4oX05si4r60vMtrHQ2hX+UyEa2ofOmKRRRz58MuXRL7UmBtRFxYWtWq9U6StE+6/TTgVRRv3v2I4pIg0CL1RsSSiJgSEW0U5+n1EXEaLVgrgKQ9Je3dc5ti7vl2WvRciIg/ABslHZKa5lBcDaAl600W8tepHRiKWof7TYwmvRHyWuC3FHO5HxnueurU9w3g98ATFKOR0ynmcq8D7kzf9xvuOlOtx1BML9wK3JK+XtvC9R4O3JzqvR04O7UfTHGdp06KX53HDXetNXW/Ari6lWtNdf06fXX0/Gy16rmQajsCaE/nw3eBfVu1XooPHjwAjC+17fJafRkGM7OMjIbpHTMzq8ihb2aWEYe+mVlGHPpmZhlx6JuZZcShb2aWEYe+mVlG/h9vN194yrGc8gAAAABJRU5ErkJggg==\n", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAEICAYAAACzliQjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAFzBJREFUeJzt3X+0XWV95/H3p4lEBQsSwigBvKHgj9gfViOoVetAa6GosQ6UoLbUocVZI2MddWlwWkVaFbo6UkeZdliCxZ9A8UczJUodaTvWKpMAVk2REmOEEH8EAigqQuA7f+yderw5N/fcm0su9z7v11p33bOf/ey9n7PPvp/znGfvs2+qCklSG35qthsgSdp7DH1JaoihL0kNMfQlqSGGviQ1xNCXpIYY+ppRSe5OcsRe3N5zktw4g+v7ZJLT+se/k+QfZ3DdL0vytzO1vils95eS3NS/Ni/e29vXQ4uhP0cl2Zzkh/0f8reTvC/JfnuwvrEklWThnrSrqvarqk17so6BNp2d5L4k3+t//jXJe5I8dmB7n62qJ4y4rg9OVq+qTqiqS2ag7bvsz6r6UFU9f0/XPQ3nAO/pX5tPjJ8508fSRJKsTPLFJN9NcluSzyQZ6+dN+lprZhj6c9sLq2o/4KnA04E/mK2G7OmbxW6Wv6yqHgUcCPwG8Bjg2pkOg3Tm69/D44ANk9R5UI+lJEcC7wdeB+wPLAP+J/DAQLW98lq3br4e5E2pqluBTwI/C5DkkCRrkmxPsjHJ7+2sm+ToJOv73ta3k7yzn/V/+9939j2+Z/b1/2OSG5LckeSqJI8bWFcleVWSm4CbBsqO7B/vn+T9SbYl+UaSP9gZrP3QyeeSnJ9kO3D2JM/xvqraAJwCbKMLD5I8L8mWgTa9McmtfW/xxiTHJTkeeBNwSv/c/rmv+/dJ3pbkc8APgCP6st8d2HSSvDvJXUm+muS4gRmbk/zKwPTgp4ld9uf44aIkz0qyrl/3uiTPGpj390n+qN9H30vyt0kOmmj/JPm9/rXe3r/2h/TlXwOOAP53345Fk+znqRxLZye5IsllfRuvS/ILE6z6KcDXq+oz1fleVX20qm4e0oahr7VmhqE/DyQ5DPh14Pq+6CPAFuAQ4CTg7QNh9S7gXVX108DPAJf35c/tfx/QDwN8Pt3475uAlwBLgM/26x70YuAYYPmQpr2brld3BPDLwG8DrxiYfwywCTgYeNsoz7Wq7gf+GnjO+HlJngCcCTy97zH+GrC5qj4FvJ2uJ7lfVQ0G028BZwCPAr4xZJM723gQ8BbgY0kOHKGpu+zPcW09ELgS+B/AYuCdwJVJFg9Ueynd/joY2Ad4/bANJTkWeAfwm8Bj++dxKUBV/QxwM31Pvqp+tLtGT/FYAlgJ/BVd7/zDwCeSPGzIqq8Dnti/yf/7jDB8tLvXWtNn6M9tn0hyJ/CPwD/Q/UEeBjwbeGNV3VNVXwTeSxduAPcBRyY5qKrurqov7Gb9rwTeUVU3VNUOuuB8ymBvv5+/vap+OLhgkgV0PbWz+l7dZuC/D7QDYGtVvbuqdoxffhJb6UJmvPuBRcDyJA+rqs1V9bVJ1vWXVbWhb8N9Q+Z/B/izvvd5GXAjcOIU2jqRE4GbquoD/bY/AnwVeOFAnfdV1b/2++Zyut7yMC8DLq6q6/pQPwt4Zvrx8hFN51gCuLaqruj33TuBhwPPGL/y/jzP84Cl/XO5LclfjhD+E73WmiZDf257cVUdUFWPq6r/3IfDIcD2qvreQL1v0P2xAZwOPB74aj+k8ILdrP9xwLuS3NkHwnYgA+sCuGWCZQ+i650O9p4H27G7ZSeztG/LT6iqjcBr6IaKvpPk0p3DHLsxWRturZ+8K+E36PbxnjqEXT9ZjN8/3xp4/ANgooD8iXVV1d3A7ePWNZnpHEswsP+q6gF+/KlgF1X1har6zapaQtd7fy7w3yZp19DXWtNn6M8/W4EDkzxqoOxw4FaAqrqpqk6lGzI4D7giyb7AsNut3gK8sg+DnT+PqKp/Gqgz0W1ab6P7VDH4qeDf2jHJshPqzwm8kG6oaRdV9eGqena/3aJ7jrvb1mRtWJokA9OH0+1jgO8DjxyY95gprHcrP7lvdq771iF1J/MT6+pfz8XTXNf49U54LPUOG9juTwGH8uP9M6GqWgd8jP7cwTCTvdaaHkN/nqmqW4B/At6R5OFJfp6ud/8hgCQvT7Kk75Xd2S92P90Jswfoxt93+gvgrCRP7pfdP8nJI7bjfrqP8W9L8qh+SOi1wKSXTQ6T5GFJnkQ3xvwYuqGE8XWekOTY/mTlPcAP++cG8G1gLFO/Qudg4NX99k8GngSs7ed9EVjVz1tBN+a907D9OWgt8PgkL02yMMkpdOdF/maK7YNuLP0VSZ7SP/e3A9f0Q2rTNtmx1Htakpeku/rqNcCPgF2GDJM8uz/ZfHA//UTgRRPUnfS11vQZ+vPTqcAYXY/r48BbqurT/bzjgQ1J7qY7qbuqH6/9Ad3J1M/1wznPqKqP0/WUL03yXeArwAlTaMd/oesNb6IbK/4wcPEUn8spfVvvBNbQDVs8raqG9SYXAefSfcr4Fl1gv6mf91f979uTXDeF7V8DHNWv823ASVV1ez/vD+lOht8BvJXu+QEwbH8OrrRfxwvorky5HXgD8IKqum0Kbdu5rs/0bfko8M2+Taumup4J7O5Ygu5E6yl0++C3gJdMcG7kTrqQ/3L/en6qX9+fDNSZymutaYr/REXSdCQ5Gziyql4+223R6OzpS1JDDH1JaojDO5LUEHv6ktSQPbpJ1oPhoIMOqrGxsdluhiTNKddee+1t/RffdushF/pjY2OsX79+tpshSXNKkmH3jtqFwzuS1BBDX5IaYuhLUkMMfUlqiKEvSQ0x9CWpIYa+JDXE0Jekhhj6ktSQh9w3ch8sY6uvHFq++dyZ+B/XkjQ3zLvQnyjcJUkO70hSUwx9SWqIoS9JDTH0Jakhhr4kNcTQl6SGGPqS1BBDX5IaYuhLUkMMfUlqyEihn+T4JDcm2Zhk9ZD5z01yXZIdSU4aN++0JDf1P6fNVMMlSVM3aegnWQBcAJwALAdOTbJ8XLWbgd8BPjxu2QOBtwDHAEcDb0ny6D1vtiRpOkbp6R8NbKyqTVV1L3ApsHKwQlVtrqovAQ+MW/bXgE9X1faqugP4NHD8DLRbkjQNo4T+UuCWgektfdkoRlo2yRlJ1idZv23bthFXLUmaqlFCP0PKasT1j7RsVV1YVSuqasWSJUtGXLUkaapGCf0twGED04cCW0dc/54sK0maYaOE/jrgqCTLkuwDrALWjLj+q4DnJ3l0fwL3+X2ZJGkWTBr6VbUDOJMurG8ALq+qDUnOSfIigCRPT7IFOBn4X0k29MtuB/6I7o1jHXBOXyZJmgUj/bvEqloLrB1X9uaBx+vohm6GLXsxcPEetFGSNEP8Rq4kNcTQl6SGGPqS1BBDX5IaYuhLUkNGunpnPhtbfeXQ8s3nnriXWyJJDz57+pLUEENfkhpi6EtSQwx9SWqIoS9JDTH0Jakhhr4kNcTQl6SGGPqS1BBDX5IaYuhLUkMMfUlqiKEvSQ0x9CWpIYa+JDXE0Jekhhj6ktQQQ1+SGmLoS1JDDH1JaoihL0kNMfQlqSGGviQ1ZKTQT3J8khuTbEyyesj8RUku6+dfk2SsL39YkkuSfDnJDUnOmtnmS5KmYtLQT7IAuAA4AVgOnJpk+bhqpwN3VNWRwPnAeX35ycCiqvo54GnAK3e+IUiS9r5RevpHAxuralNV3QtcCqwcV2clcEn/+ArguCQBCtg3yULgEcC9wHdnpOWSpCkbJfSXArcMTG/py4bWqaodwF3AYro3gO8D3wRuBv60qrbvYZslSdM0SuhnSFmNWOdo4H7gEGAZ8LokR+yygeSMJOuTrN+2bdsITZIkTccoob8FOGxg+lBg60R1+qGc/YHtwEuBT1XVfVX1HeBzwIrxG6iqC6tqRVWtWLJkydSfhSRpJKOE/jrgqCTLkuwDrALWjKuzBjitf3wScHVVFd2QzrHp7As8A/jqzDRdkjRVk4Z+P0Z/JnAVcANweVVtSHJOkhf11S4CFifZCLwW2HlZ5wXAfsBX6N483ldVX5rh5yBJGtHCUSpV1Vpg7biyNw88vofu8szxy909rFySNDv8Rq4kNcTQl6SGGPqS1BBDX5IaYuhLUkMMfUlqiKEvSQ0x9CWpIYa+JDXE0Jekhhj6ktQQQ1+SGmLoS1JDRrrLZovGVl85tHzzuSfu5ZZI0syxpy9JDTH0Jakhhr4kNcTQl6SGGPqS1BBDX5IaYuhLUkMMfUlqiKEvSQ0x9CWpIYa+JDXE0Jekhhj6ktQQQ1+SGmLoS1JDDH1JashIoZ/k+CQ3JtmYZPWQ+YuSXNbPvybJ2MC8n0/y+SQbknw5ycNnrvmSpKmYNPSTLAAuAE4AlgOnJlk+rtrpwB1VdSRwPnBev+xC4IPAf6qqJwPPA+6bsdZLkqZklJ7+0cDGqtpUVfcClwIrx9VZCVzSP74COC5JgOcDX6qqfwaoqtur6v6ZabokaapGCf2lwC0D01v6sqF1qmoHcBewGHg8UEmuSnJdkjcM20CSM5KsT7J+27ZtU30OkqQRjRL6GVJWI9ZZCDwbeFn/+zeSHLdLxaoLq2pFVa1YsmTJCE2SJE3HKKG/BThsYPpQYOtEdfpx/P2B7X35P1TVbVX1A2At8NQ9bbQkaXoWjlBnHXBUkmXArcAq4KXj6qwBTgM+D5wEXF1VleQq4A1JHgncC/wy3YneOWts9ZVDyzefe+JebokkTd2koV9VO5KcCVwFLAAurqoNSc4B1lfVGuAi4ANJNtL18Ff1y96R5J10bxwFrK2q4akpSXrQjdLTp6rW0g3NDJa9eeDxPcDJEyz7QbrLNiVJs8xv5EpSQwx9SWqIoS9JDTH0Jakhhr4kNcTQl6SGGPqS1BBDX5IaYuhLUkMMfUlqiKEvSQ0x9CWpIYa+JDXE0Jekhhj6ktQQQ1+SGmLoS1JDDH1JaoihL0kNMfQlqSGGviQ1xNCXpIYY+pLUEENfkhpi6EtSQwx9SWqIoS9JDTH0Jakhhr4kNWSk0E9yfJIbk2xMsnrI/EVJLuvnX5NkbNz8w5PcneT1M9NsSdJ0TBr6SRYAFwAnAMuBU5MsH1ftdOCOqjoSOB84b9z884FP7nlzJUl7YuEIdY4GNlbVJoAklwIrgX8ZqLMSOLt/fAXwniSpqkryYmAT8P0Za/VD0NjqK4eWbz73xL3cEkma2CjDO0uBWwamt/RlQ+tU1Q7gLmBxkn2BNwJv3d0GkpyRZH2S9du2bRu17ZKkKRol9DOkrEas81bg/Kq6e3cbqKoLq2pFVa1YsmTJCE2SJE3HKMM7W4DDBqYPBbZOUGdLkoXA/sB24BjgpCR/AhwAPJDknqp6zx63XJI0ZaOE/jrgqCTLgFuBVcBLx9VZA5wGfB44Cbi6qgp4zs4KSc4G7jbwJWn2TBr6VbUjyZnAVcAC4OKq2pDkHGB9Va0BLgI+kGQjXQ9/1YPZaEnS9IzS06eq1gJrx5W9eeDxPcDJk6zj7Gm0T5I0g/xGriQ1xNCXpIYY+pLUEENfkhpi6EtSQwx9SWqIoS9JDTH0Jakhhr4kNcTQl6SGGPqS1BBDX5IaYuhLUkNGusumps//nSvpocSeviQ1xNCXpIYY+pLUEENfkhpi6EtSQwx9SWqIoS9JDTH0Jakhhr4kNcTQl6SGGPqS1BBDX5IaYuhLUkMMfUlqiKEvSQ0ZKfSTHJ/kxiQbk6weMn9Rksv6+dckGevLfzXJtUm+3P8+dmabL0maiklDP8kC4ALgBGA5cGqS5eOqnQ7cUVVHAucD5/XltwEvrKqfA04DPjBTDZckTd0oPf2jgY1Vtamq7gUuBVaOq7MSuKR/fAVwXJJU1fVVtbUv3wA8PMmimWi4JGnqRvl3iUuBWwamtwDHTFSnqnYkuQtYTNfT3+k/ANdX1Y+m39z5w3+jKGk2jBL6GVJWU6mT5Ml0Qz7PH7qB5AzgDIDDDz98hCZJkqZjlOGdLcBhA9OHAlsnqpNkIbA/sL2fPhT4OPDbVfW1YRuoqgurakVVrViyZMnUnoEkaWSjhP464Kgky5LsA6wC1oyrs4buRC3AScDVVVVJDgCuBM6qqs/NVKMlSdMzaehX1Q7gTOAq4Abg8qrakOScJC/qq10ELE6yEXgtsPOyzjOBI4E/TPLF/ufgGX8WkqSRjDKmT1WtBdaOK3vzwON7gJOHLPfHwB/vYRslSTPEb+RKUkMMfUlqiKEvSQ0x9CWpIYa+JDXE0Jekhhj6ktSQka7T197jjdgkPZjs6UtSQwx9SWqIoS9JDTH0Jakhhr4kNcTQl6SGGPqS1BBDX5IaYuhLUkP8Ru4c4Td1Jc0Ee/qS1BBDX5IaYuhLUkMMfUlqiCdy57iJTvCCJ3kl7cqeviQ1xNCXpIYY+pLUEMf057HdjfcP4zkAaf6zpy9JDTH0JakhDu/o33h/H2n+Gyn0kxwPvAtYALy3qs4dN38R8H7gacDtwClVtbmfdxZwOnA/8OqqumrGWq+9wnMD0vwxaegnWQBcAPwqsAVYl2RNVf3LQLXTgTuq6sgkq4DzgFOSLAdWAU8GDgH+T5LHV9X9M/1E9NDhJwbpoWuUnv7RwMaq2gSQ5FJgJTAY+iuBs/vHVwDvSZK+/NKq+hHw9SQb+/V9fmaar7lkqp8YZsru3mxm6g3KNzrNFaOE/lLgloHpLcAxE9Wpqh1J7gIW9+VfGLfs0vEbSHIGcEY/eXeSG0dq/Y8dBNw2xWVa4H4Bct4uRZPulyHLzNS2H6o8VoabS/vlcaNUGiX0M6SsRqwzyrJU1YXAhSO0Zagk66tqxXSXn6/cL8O5X3blPhluPu6XUS7Z3AIcNjB9KLB1ojpJFgL7A9tHXFaStJeMEvrrgKOSLEuyD92J2TXj6qwBTusfnwRcXVXVl69KsijJMuAo4P/NTNMlSVM16fBOP0Z/JnAV3SWbF1fVhiTnAOurag1wEfCB/kTtdro3Bvp6l9Od9N0BvOpBunJn2kND85z7ZTj3y67cJ8PNu/2SrkMuSWqBt2GQpIYY+pLUkDkf+kmOT3Jjko1JVs92e2ZDksOS/F2SG5JsSPL7ffmBST6d5Kb+96Nnu62zIcmCJNcn+Zt+elmSa/r9cll/gUJTkhyQ5IokX+2Pm2e2frwk+a/9389XknwkycPn47Eyp0N/4BYRJwDLgVP7Wz+0Zgfwuqp6EvAM4FX9flgNfKaqjgI+00+36PeBGwamzwPO7/fLHXS3EWnNu4BPVdUTgV+g2z/NHi9JlgKvBlZU1c/SXbSy85Yy8+pYmdOhz8AtIqrqXmDnLSKaUlXfrKrr+sffo/sDXkq3Ly7pq10CvHh2Wjh7khwKnAi8t58OcCzd7UKgwf2S5KeB59JddUdV3VtVd+LxshB4RP9do0cC32QeHitzPfSH3SJil9s8tCTJGPCLwDXAv6uqb0L3xgAcPHstmzV/BrwBeKCfXgzcWVU7+ukWj5kjgG3A+/phr/cm2ZeGj5equhX4U+BmurC/C7iWeXiszPXQH+k2D61Ish/wUeA1VfXd2W7PbEvyAuA7VXXtYPGQqq0dMwuBpwJ/XlW/CHyfhoZyhunPX6wEltHdEXhfumHj8eb8sTLXQ9/bPPSSPIwu8D9UVR/ri7+d5LH9/McC35mt9s2SXwJelGQz3dDfsXQ9/wP6j/DQ5jGzBdhSVdf001fQvQm0fLz8CvD1qtpWVfcBHwOexTw8VuZ66I9yi4h5rx+nvgi4oareOTBr8PYYpwF/vbfbNpuq6qyqOrSqxuiOjaur6mXA39HdLgTa3C/fAm5J8oS+6Di6b823fLzcDDwjySP7v6ed+2TeHStz/hu5SX6drve28xYRb5vlJu11SZ4NfBb4Mj8eu34T3bj+5cDhdAf1yVW1fVYaOcuSPA94fVW9IMkRdD3/A4HrgZf3//OhGUmeQndyex9gE/AKuk5gs8dLkrcCp9BdDXc98Lt0Y/jz6liZ86EvSRrdXB/ekSRNgaEvSQ0x9CWpIYa+JDXE0Jekhhj6ktQQQ1+SGvL/ARuZcyj6Jn9hAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] @@ -204,7 +204,7 @@ }, { "cell_type": "code", - "execution_count": 353, + "execution_count": 358, "metadata": {}, "outputs": [], "source": [ @@ -217,7 +217,7 @@ }, { "cell_type": "code", - "execution_count": 354, + "execution_count": 359, "metadata": {}, "outputs": [ { @@ -227,14 +227,14 @@ "Chance is empirical probability that given school is best\n", " effect se Chance\n", "school \n", - "A 28 15 0.191\n", - "B 8 10 0.109\n", - "C -3 16 0.105\n", - "D 7 11 0.124\n", - "E -1 9 0.077\n", - "F 1 11 0.088\n", - "G 18 10 0.186\n", - "H 12 18 0.120\n", + "A 28 15 0.2047\n", + "B 8 10 0.1160\n", + "C -3 16 0.1010\n", + "D 7 11 0.1075\n", + "E -1 9 0.0766\n", + "F 1 11 0.0854\n", + "G 18 10 0.1761\n", + "H 12 18 0.1327\n", "Only thing that worries me is that Gelman has A best with prob=10%\n" ] } @@ -248,7 +248,7 @@ }, { "cell_type": "code", - "execution_count": 332, + "execution_count": 360, "metadata": {}, "outputs": [ { @@ -259,14 +259,14 @@ " as good as or better than corresponding column\n", " \n", " A B C D E F G H\n", - "A 1.00 0.61 0.64 0.60 0.68 0.66 0.53 0.58 \n", - "B 0.39 1.00 0.53 0.50 0.59 0.56 0.41 0.47 \n", - "C 0.36 0.47 1.00 0.47 0.55 0.52 0.37 0.44 \n", - "D 0.40 0.50 0.53 1.00 0.58 0.56 0.42 0.47 \n", - "E 0.32 0.41 0.45 0.42 1.00 0.47 0.33 0.38 \n", - "F 0.34 0.44 0.48 0.44 0.53 1.00 0.35 0.42 \n", - "G 0.47 0.59 0.63 0.58 0.67 0.65 1.00 0.57 \n", - "H 0.42 0.53 0.56 0.53 0.62 0.58 0.43 1.00 \n" + "A 1.00 0.59 0.62 0.60 0.66 0.64 0.52 0.58 \n", + "B 0.41 1.00 0.54 0.51 0.57 0.56 0.42 0.49 \n", + "C 0.38 0.46 1.00 0.47 0.54 0.51 0.39 0.45 \n", + "D 0.40 0.49 0.53 1.00 0.56 0.54 0.41 0.48 \n", + "E 0.34 0.43 0.46 0.44 1.00 0.48 0.35 0.41 \n", + "F 0.36 0.44 0.49 0.46 0.52 1.00 0.38 0.43 \n", + "G 0.48 0.58 0.61 0.59 0.65 0.62 1.00 0.57 \n", + "H 0.42 0.51 0.55 0.52 0.59 0.57 0.43 1.00 \n" ] } ], @@ -336,7 +336,7 @@ " }\n", "}\n", "'''\n", - "sm=pystan.StanModel(model_code=stan_code_2)" + "sm2=pystan.StanModel(model_code=stan_code_2)" ] }, { @@ -345,7 +345,7 @@ "metadata": {}, "outputs": [], "source": [ - "answers=sm.sampling(data=dict({'means':p55['effect'],'se':p55['se']}),iter=5000)\n" + "answers=sm2.sampling(data=dict({'means':p55['effect'],'se':p55['se']}),iter=5000)\n" ] }, { From 92c245e86749a901a24377bb323d9961935f82df Mon Sep 17 00:00:00 2001 From: Jeremy Teitelbaum Date: Sun, 29 Apr 2018 13:22:49 -0400 Subject: [PATCH 7/8] More on 5.9.8 --- BDA 5.9.8.ipynb | 97 ++++++++++++++++++++++++++++++------------- Useful Formulae.ipynb | 2 +- 2 files changed, 68 insertions(+), 31 deletions(-) diff --git a/BDA 5.9.8.ipynb b/BDA 5.9.8.ipynb index d3ea631..5228ce3 100644 --- a/BDA 5.9.8.ipynb +++ b/BDA 5.9.8.ipynb @@ -46,70 +46,107 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "In the special case under consideration, $p_1$ is normal with mean $-1$ and $\\sigma=.5$, $p_2$ is normal with mean $1$ and $\\sigma=.5$ and we can set $\\lambda_1=.1$ and $\\lambda_2=.9$. The $p_m(\\{y_{i}\\})$ can be calculated from the $t$ distribution. Drawing a sample of size $10$ from $p_1$ and getting a sample mean of $-.25$ and a sample variance of $1$ gives a $t$-statistics of $\\sqrt{10}(-.25+1)$ in the first case and $\\sqrt{10}(-.25-1)$ in the second. " + "In the special case under consideration, $p_1$ is normal with mean $-1$ and $\\sigma=.5$, $p_2$ is normal with mean $1$ and $\\sigma=.5$ and we can set $\\lambda_1=.1$ and $\\lambda_2=.9$. The $p_m(\\{y_{i}\\})$ can be calculated from the $t$ distribution. Drawing a sample of size $10$ from $p_1$ and getting a sample mean of $-.25$ and a sample variance of $1$ gives a $t$-statistics of \n", + "$$\\frac{(\\overline{y}-\\mu)}{s/\\sqrt{N}}=\\frac{(-.25+1)}{1/\\sqrt{10}}=\\sqrt{10}(-.25+1)$$ in the first case and $\\sqrt{10}(-.25-1)$ in the second. " ] }, { "cell_type": "code", - "execution_count": 14, + "execution_count": 23, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "from scipy.stats import norm, t\n", + "import matplotlib.pyplot as plt\n", + "t_1=np.sqrt(10)*.75\n", + "t_2=np.sqrt(10)*(-1.25)\n", + "p1y=t.pdf(t_1,df=9)\n", + "p2y=t.pdf(t_2,df=9)\n", + "lambda1,lambda2=(.1,.9)\n" + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "metadata": {}, + "outputs": [], + "source": [ + "wt1=lambda1*p1y/(lambda1*p1y+lambda2*p2y)\n", + "wt2=lambda2*p2y/(lambda1*p1y+lambda2*p2y)" + ] + }, + { + "cell_type": "code", + "execution_count": 25, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "0.041664931082753924\n", - "0.0035119750957915393\n" + "0.600590082806086 0.3994099171939141\n" ] } ], "source": [ - "import numpy as np\n", - "from scipy.stats import norm, t\n", - "t_1=np.sqrt(9)*.75\n", - "t_2=np.sqrt(9)*1.25\n", - "print(t.pdf(t_1,df=9))\n", - "print(t.pdf(t_2,df=9))" + "print(wt1, wt2)" ] }, { "cell_type": "code", - "execution_count": 18, + "execution_count": 26, "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "0.5655172413793104" - ] - }, - "execution_count": 18, - "metadata": {}, - "output_type": "execute_result" - } - ], + "outputs": [], + "source": [ + "def posterior(prior_mean,prior_variance,sample_mean,pop_variance,n):\n", + " post_var=1/((1/prior_variance) + n/pop_variance)\n", + " post_mean=(prior_mean/prior_variance+sample_mean*n/pop_variance)/(1/post_var)\n", + " return post_mean, post_var" + ] + }, + { + "cell_type": "code", + "execution_count": 27, + "metadata": {}, + "outputs": [], "source": [ - ".1*.041/(.1*.041+.9*.0035)" + "post_mean1,post_var1=posterior(-1,.25,-.25,1,10)\n", + "post_mean2,post_var2=posterior(1,.25,-.25,1,10)" ] }, { "cell_type": "code", - "execution_count": 17, + "execution_count": 28, "metadata": {}, "outputs": [ { "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl8VOX1+PHPSQDD6gYiAhpUUMBElriwWFFcICCIIoJVwd1atOpXf6Kt2tL67dd9qUu1VbSWgiCCKCAILkgBJW6ALMoSMIKK7IgsCc/vj5OBEBIyJPfOvXfmvF+veU1mcufeM4Gc3Dn3ec4jzjmMMcYkl7SgAzDGGOM9S+7GGJOELLkbY0wSsuRujDFJyJK7McYkIUvuxhiThCy5G2NMErLkbowxSciSuzHGJKFqFW0gIi8BPYEfnXMnlfF9AZ4EcoGtwCDn3GcV7bd+/fouMzPzgAM2xphU9umnn/7knGtQ0XYVJnfgZeBp4F/lfL870Lz4dhrwXPH9fmVmZpKXlxfH4Y0xxsSIyIp4tquwLOOcmw6s288mvYF/OTUbOEREGsUXpjHGGD94UXNvDHxb4nFB8XPGGGMC4kVylzKeK7PVpIhcLyJ5IpK3Zs0aDw5tjDGmLPHU3CtSADQt8bgJsKqsDZ1zLwAvAOTk5FivYWOS2M6dOykoKGDbtm1BhxJJGRkZNGnShOrVq1fq9V4k9/HAYBEZiV5I3eicW+3Bfo0xEVZQUEDdunXJzMxEB9WZeDnnWLt2LQUFBTRr1qxS+4hnKOQIoAtQX0QKgPuB6sUB/B2YiA6DXIIOhbyqUpEYY5LKtm3bLLFXkohw+OGHU5XydYXJ3Tk3oILvO+C3lY7AGJO0LLFXXlV/dl6UZYwJxNy5MGkStGgBvXtDms23NmY3+3UwkfTKK9CuHQwZAhddBBdcANu3Bx2Viar77ruPqVOnBh2Gpyy5m8hZtAiuuw7OOgu+/x6eegomToTbbgs6MhNFRUVFDB06lHPOOeeAXhN2ltxN5Pzud1C7NgwfDg0bws03w513wnPPwXvvBR2dCZP8/HxOPPFEBg4cSHZ2Nn379mXr1q1kZmYydOhQOnfuzOjRoxk0aBCvv/46ANOmTaNt27ZkZWVx9dVXs734I2Hp14Sd1dxNpHzxBUyZAg89BEccsef5oUNh9Gi49VbdxurvIRP7h/FSmzbwxBMVbrZ48WJefPFFOnXqxNVXX82zzz4L6DjyGTNmAPDOO+8AOsJn0KBBTJs2jRYtWnDllVfy3HPPceutt+7zmrCzXwETKX/7G9SqBddeu/fzGRnwwAMwbx6MHx9MbCacmjZtSqdOnQC4/PLLdyfnSy+9dJ9tFy9eTLNmzWjRogUAAwcOZPr06bu/X9ZrwsrO3E1kbNumZ+eXXgqHHrrv9/v1gz/8AR58UEfP2Ci8EInjDNsvpYcUxh7Xrl17n211ZHf5ynpNWNmZu4mMKVNg82ZN7mWpVg3+539g9my9GQOwcuVKZs2aBcCIESPo3LlzudueeOKJ5Ofns2TJEgBeffVVzjzzzITE6TVL7iYyxozRM/azzy5/m4EDoU4d+Oc/ExeXCbeWLVvyyiuvkJ2dzbp16/jNb35T7rYZGRkMGzaMSy65hKysLNLS0rjxxhsTGK13rCxjIsE5PXPv1g3210epTh09sx8xAh5/HOrVS1yMJpzS0tL4+9//vtdz+fn5ez1++eWXd3/dtWtXPv/88332U/o1YWdn7iYSFi7UMe1du1a87bXXwtat8Npr/sdlTFhZcjeRMG2a3seT3E87DVq3hmHD/I3JhF9mZibz588POoxAWHI3kfDee9CsGcSzproIXH45zJoFEfskbYxnLLmb0Nu1Cz78cP8XUkvr31/vR470JyZjws6Suwm9JUtg/Xro2DH+12RmQocOemHVmFRkyd2E3ief6P2ppx7Y6wYM0LbACxZ4H5MxYWfJ3YTeJ59oo7CWLQ/sdf36aY8ZO3s3VTFu3DgWVOIMYfz48fzf//2fDxHFx5K7Cb1PPoGcHEhPP7DXNWyodfoRI3ScvDGVUZnkXlhYSK9evRgyZMgBvcZLltxNqO3YAZ9/fuAlmZjLLoOlS2HOHG/jMtFQXsvf8tr6DhkyhFatWpGdnc0dd9zBzJkzGT9+PHfeeSdt2rRh6dKlLF26lG7dutG+fXvOOOMMFi1aBMCgQYO4/fbbOeuss7jrrrt4+eWXGTx4MAArVqyga9euZGdn07VrV1auXFnma7xkM1RNqM2frwk+J6dyr+/TB268Uc/eK/sHwlRdgB1/92n5+9hjj/H888/v09b3yiuvZOzYsSxatAgRYcOGDRxyyCH06tWLnj170rdvX0BnsP7973+nefPmfPzxx9x00028V7yQwNdff83UqVNJT0/fa9br4MGDufLKKxk4cCAvvfQSt9xyC+PGjdvnNV6yM3cTavPm6f3JJ1fu9YccAj166JDICCyeY3xQuuXvtGnTymzrW69ePTIyMrj22mt54403qFWr1j772rJlCzNnzuSSSy6hTZs23HDDDaxevXr39y+55JIyk/SsWbO47LLLALjiiiv26glf3muqys7cTajNm6e92o8/vvL7uOwyGDsWPvggvhmuxnsBdvzdp+VveapVq8Ynn3zCtGnTGDlyJE8//fTuM/KYXbt2ccghh/BFOR9D4m0JXDImv9oI25m7CbW5c6FVqwO/mFpSjx5Qty785z/exWWio3TL33POOafMtr5btmxh48aN5Obm8sQTT+xO4HXr1mXz5s0A1KtXj2bNmu1eZs85x5dffllhDB07dmRk8Yy64cOH77ftsFcsuZtQmzcPsrKqto+aNbX2PmYMFF83MymkdMvf2267rcy2vps3b6Znz55kZ2dz5pln8vjjjwPQv39/Hn74Ydq2bcvSpUsZPnw4L774IieffDKtW7fmzTffrDCGp556imHDhpGdnc2rr77Kk08+6ffbRipaecQvOTk5Li8vL5Bjm2j46Sdo0AAeeUQX4aiKyZO1XfDYsXDhhd7EZ/Zv4cKFtDzQyQkey8/Pp2fPnpFtHlbWz1BEPnXOVTjEwM7cTWjFLqZW9cwdtNbeoIGVZkzqsORuQiuW3LOzq76vatV0EY+33oJNm6q+PxMN1vLXmBCaOxfq19eZpl4YMEAX2R471pv9mYoFVfZNBlX92VlyN6G1aJGOlIlzJFuFOnSAY4+FEnNLjI8yMjJYu3atJfhKcM6xdu1aMjIyKr0PG+duQmvxYh3l4hURuO46uPtu+PprKJ7DYnzSpEkTCgoKWLNmTdChRFJGRgZNmjSp9OstuZtQWrdOR8t4nYAHDYJ774V//AMeftjbfZu9Va9enWbNmgUdRsqysowJpa+/1vsTTvB2v0ceCb17a2nGxrybZBZXcheRbiKyWESWiMg+PSxF5GgReV9EPheRuSKS632oJpXEkrsfpZPrr9dPBcV9m4xJShUmdxFJB54BugOtgAEi0qrUZn8ARjnn2gL9gWe9DtSklsWLteXAscd6v+9zztHFtp9+2vt9GxMW8Zy5nwoscc4tc87tAEYCvUtt44B6xV8fDKzyLkSTir7+WhN79ere7zstDW65BWbMgI8/9n7/xoRBPMm9MfBticcFxc+V9EfgchEpACYCN3sSnUlZixd7X28v6ZprtB3wI4/4dwxjghRPci9rlHHpgasDgJedc02AXOBVEdln3yJyvYjkiUieDY8y5dm1C775xt/kXreuLuLxxhu6UpMxySae5F4ANC3xuAn7ll2uAUYBOOdmARlA/dI7cs694JzLcc7lNGjQoHIRm6RXUKAzSf0eh37LLdqW4MEH/T2OMUGIJ7nPAZqLSDMRqYFeMB1fapuVQFcAEWmJJnc7NTeVsnix3vt55g7QqJGOnBk2zM7eTfKpMLk75wqBwcBkYCE6KuYrERkqIr2KN/sf4DoR+RIYAQxyNufYVJKfwyBLu+cePXsfOtT/YxmTSHHNUHXOTUQvlJZ87r4SXy8AOnkbmklVS5ZA7do64chvjRrBb38Ljz+ubQlOPNH/YxqTCDZD1YTOihWQmeldw7CK3HWXrtZ0//2JOZ4xiWDJ3YROfj4cc0zijtegAdx2G4waBZ99lrjjGuMnS+4mdPLz9cw9ke64Aw47TGvwxiQDS+4mVDZtgvXrE5/cDz5Ya+6TJ8MHHyT22Mb4wZK7CZUVK/Q+kWWZmN/+Fho31iRvY71M1FlyN6GSn6/3iT5zhz0XVWfP1rVWjYkyS+4mVGJn7kEkd4CrrtLx9ffcA0VFwcRgjBcsuZtQyc/XM+igulNUqwZ/+Qt89RX85z/BxGCMFyy5m1CJDYNM1Bj3slx8MbRrB/fdZ6s1meiy5G5CJYhhkKWlpcEDD2gsI0YEG4sxlWXJ3YTKihXBjJQp7fzzISsLHnssoiNnnNO/TgsWwC+/BB2NCYAldxMaW7bo2qZBn7mDloVuvx3mzYOpU4OO5gDs3AkPPQRNm+pagq1b66okAwbsGYpkUoIldxMaQY+UKW3AAGjYEJ54IuhI4lRQAB06aLOcVq3guedg+HBdleStt/S5114LOkqTIHF1hTQmEYKcwFSWgw6Ca6+Fv/4VvvtOJziF1ooVcNZZ+tFnzBi46KI937vsMu2vcNll0L+/TgO+7rrgYjUJYWfuJjSCnMBUnquu0mX/Xnkl6Ej2Y/NmuOACWLcOpk3bO7HHNG0K774L3bvrmfykSYmP0ySUJXcTGvn5erbcsGHQkexx3HHQpQu89FKIL6zecINeOB09Gk45pfztMjK09WV2Nvz611rGMUnLkrsJjdhImbSQ/a+85hpdhu+//w06kjK88YaO17z/fjj33Iq3r1NHE/yOHTBoUIj/YpmqCtmvkUllie7jHq8LL9ST3tGjg46klHXr4De/0RlXQ4bE/7rmzeHRR7WEYwP5k5YldxMaYZjAVJY6dbRU/frrWn8Pjb/8RS+gvvQSVK9+YK+97jot4dxxh9bsTdKx5G5C4Zdf4Mcfw5ncAfr2hVWrtGNkKCxdCk8/rVd8Tz75wF+flqavX70aHnzQ+/hM4Cy5m1AI2zDI0nr2hBo19Ow9FO69V8/Whw6t/D5OPRX69YMnn9RPACapWHI3oRDGYZAl1aun1yvHjQvBNchvvtHJSIMHw1FHVW1f998PP/+sNXiTVCy5m1AIe3IHyM2F5cs1twbqoYf0rP2226q+r1atdGLT3/6m6xuapGHJ3YTCihWarxo1CjqS8nXvrveBzv/57judUXX11XDkkd7s86679Oz9H//wZn8mFCy5m1DIz4ejjw7fGPeSmjWDE06Ad94JMIinntIhO3fe6d0+Tz4Zzj5bz9537vRuvyZQIf5VMqkkrMMgS+vWDT74IKAuutu2wYsvQu/e+pfGS7fdpjNWx4zxdr8mMJbcTSiEpY97Rbp31xw7fXoAB3/9dVi7VicueS03Vyc3Pfmk9/s2gbDkbgK3bZsOt47CmXvnzrrO6gcfBHDw557TBHz22d7vOy1N/2jMng3z53u/f5NwltxN4Fau1PsoJPfatXV4eMKT+9y5MHOmdnT068LEFVfoVe0XX/Rn/yahLLmbwIV9AlNpZ54JeXm6clTCPP+8NrgZNMi/Y9Svr/X8V1+1lcGTgCV3E7gojHEv6cwzobBQT6QTYscObfDVpw8cdpi/x7rmGq3rjx/v73GM7yy5m8Dl52sdu6qTLROlUydIT4cPP0zQASdO1AlGV1zh/7HOPVcX9rDSTOTFldxFpJuILBaRJSJSZm9REeknIgtE5CsR+Y+3YZpklp+v+aRaRBZ9rFMHcnISmNxffRWOOCK+fu1VlZ4OAwfClCl6ldtEVoXJXUTSgWeA7kArYICItCq1TXPgbqCTc641cKsPsZokFZVhkCV17qx19x07fD7Q+vXw9tu6Wnei/vpddpk20Bk1KjHHM76I58z9VGCJc26Zc24HMBLoXWqb64BnnHPrAZxzP3obpklmUZnAVNLpp+s1xy+/9PlAo0frX5BElGRiWrbUWasjRybumMZz8ST3xsC3JR4XFD9XUgughYj8V0Rmi0g3rwI0yW3HDu2THsXkDgno7/7vf8OJJ+pqS4nUv7++ueXLE3tc45l4kruU8VzppqfVgOZAF2AA8E8ROWSfHYlcLyJ5IpK3Zs2aA43VJKFvv9UKQNTKMk2aQOPGPif3lSvho490MWsp69fQR/376/1rryX2uMYz8ST3AqBpicdNgFVlbPOmc26nc245sBhN9ntxzr3gnMtxzuU0aNCgsjGbJBK1YZAlnX66z8n9jTf0/tJLfTxIOTIzoUMHW2M1wuJJ7nOA5iLSTERqAP2B0oNgxwFnAYhIfbRMs8zLQE1yinpyX7ZMlwf0xeuvQ3a2thwIQv/+OjN2wYJgjm+qpMLk7pwrBAYDk4GFwCjn3FciMlREehVvNhlYKyILgPeBO51za/0K2iSPFSt0Nn3j0ldxIiBWd//4Yx92vnq1zpK6+GIfdh6nfv20HDR6dHAxmEqLa5y7c26ic66Fc+4459wDxc/d55wbX/y1c87d7pxr5ZzLcs7ZZXYTl/x8rV9Xrx50JAeuXTsdnThrlg87HztWL0YEmdyPPBI6dtRYTOTYDFUTqCgOg4ypVUurJnPm+LDz11/XUTKtWlW8rZ/69NHxnjZqJnIsuZtARXECU0nt2sFnn3m8aPaaNTr99eKLEz9KprQ+ffR+3Lhg4zAHzJK7CczOnbr4T1TP3EGT+7p1e9oWe+LNN3Upvb59PdxpJR17rH48sdJM5FhyN4EpKNAcFuXk3r693n/2mYc7HTNGk+rJJ3u40yro0wdmzPBxWJDxgyV3E5jYMMgol2WysrTXlmfJfcsWeO89uPDC4EsyMX36aN3J2gBHiiV3E5jYIh1RPnOvWVOveXqW3KdO1Z4MPXt6tEMPZGfrgtxWmokUS+4mMPn5enLatGmFm4Za7KKqJyZMgHr1tO1kWIjo2fvUqbBpU9DRmDhZcjeByc/XBTpq1Ag6kqpp1w6+/96D9ufO6cIc550XvoH/vXvrJ4p33w06EhMnS+4mMCtWRLskExNr2Fjls/cvvtAWmT16VDkmz3XsCIccop8sTCRYcjeBifIEppLatNHKxaefVnFHb7+tO+re3ZO4PFWtGnTrpsl9166gozFxsORuAlFYqO1+kyG516kDJ5zgQXKfMAFOOQUaNvQkLs/17KnDIfPygo7ExMGSuwnEqlVQVBTtYZAltWmjDRQrbc0a+OSTcJZkYrp10y5vVpqJBEvuJhBRbvVbluxsfU8bN1ZyB5Mm6QXVMCf3ww/XHu9vvx10JCYOltxNIJIxuQPMn1/JHUyYoF0Y27b1LCZf9OypV45XlV6vx4SNJXcTiNgEpqiPcY+JJfdKlWZ27oTJkyE3V8seYRb7ZDFxYrBxmAqF/H+SSVb5+dCoEWRkBB2JN5o00ZGClUruM2dqPSfMJZmYk06Co4+20kwEWHI3gUiWYZAxInr2XqnkPmGCTlo691zP4/KciP4Revdd2LYt6GjMflhyN4HIz0+ekTIx2dkwb14leru//TaceSbUretLXJ7r2RO2btWe8ya0LLmbhCsqSp4x7iVlZcHmzXuuJ8Rl+XJYuDAaJZmYs87SjmlWmgk1S+4m4Vav1muIyZbcK3VRNTZmPErJvWZN6NpVk7unS1AZL1lyNwmXDH3cy3LSSXp/wMm9eXO9RUmPHvoPuWhR0JGYclhyNwmXDH3cy1KnDhx33AEk959/hvffj9ZZe0xurt7bbNXQsuRuEm75cr0/+uhg4/DDAY2Yee892L49msn96KP1o4qNdw8tS+4m4fLzdTJmrVpBR+K97Gz45hsdTFKhCRP0dP9Xv/I9Ll/06AEffVSFngvGT5bcTcItX558JZmY7GztiLtgQQUbOqfJ/bzzortaSW6utvecOjXoSEwZLLmbhFu+XJfkTEZZWXpfYWlm7lwoKIhmSSamY0c4+GCru4eUJXeTULE+7sma3I89VstNFSb3WEKMXZiMomrV4Pzzte5uC3iEjiV3k1DffacJPlnLMunp0Lq1zlTdrwkToH17vfgQZT16wA8/wOefBx2JKcWSu0mo2EiZZD1zhz0jZsqd37N2LcyeHe2STEy3btpvxkozoWPJ3SRUKiT3rCz46Sc9oS3TO+9oGSMZkvsRR+jSgDYkMnQsuZuEys/XE71k6eNellgbgnJLMxMmaFLMyUlYTL7q0UOXCFyzJuhITAmW3E1CLV+uvc+jOvovHrERM2Um98JCXVIvCgtzxCs3V2tQkyYFHYkpIa7/XSLSTUQWi8gSERmyn+36iogTkSQ5JTFeS+ZhkDH16+t10jJHzMyaBRs2JEdJJqZdO2jY0EozIVNhcheRdOAZoDvQChggIq3K2K4ucAvwsddBmuSRbIt0lCfW230fEyboEMIoLMwRr7Q0PXufPFk/mZhQiOfM/VRgiXNumXNuBzAS6F3Gdn8GHgJseRZTpu3bdShksp+5g5ZmFiwoI9dNmABnnKGTf5JJbq5+Ipk1K+hITLF4kntj4NsSjwuKn9tNRNoCTZ1z++3eLyLXi0ieiOStsYsvKWflSi3Npkpy37YNliwp8eTKlTB/fnKVZGLOPVc/kdiQyNCIJ7lLGc/tHsErImnA48D/VLQj59wLzrkc51xOgwYN4o/SJIXYMMhUKctAqdJMFBfmiNfBB+snEqu7h0Y8yb0AKDlwrQmwqsTjusBJwAcikg+cDoy3i6qmtNgiHalw5t6ypc5W3Su5v/02HH88nHBCYHH5KjdX3/DKlUFHYogvuc8BmotIMxGpAfQHxse+6Zzb6Jyr75zLdM5lArOBXs65PF8iNpG1fLl+cm/cuOJtoy4jQxdX2j1i5uefYdo0XVxayvownARin0js7D0UKkzuzrlCYDAwGVgIjHLOfSUiQ0Wkl98BmuSxfLmu8ZCeHnQkibHXiJnYwhw9ewYak69OPFFrbpbcQ6FaPBs55yYCE0s9d18523apelgmGeXnp0ZJJiYrC0aNgi1boM7bb0PdulqXTlYievY+bJheTc7ICDqilJYkU+RMFCxblnrJHWD+PKf19vPPT+6puaDJfetW+PDDoCNJeZbcTUJs3KitR5o3DzqSxNk9YmbCSli1KrlLMjFdukDNmjYkMgQsuZuEWLpU748/Ptg4EumYY3SJ1Lnv/qAli+7dgw7JfzVrwtlna3Ivt+exSQRL7iYhvvlG71MpuaelaWlm3sJ0OO007QSZCnJztQb39ddBR5LSLLmbhIjN1DzuuGDjSLSs47Yyb3MmrkcKlGRiYkMirTQTKEvuJiGWLIGjjoLatYOOJLGydn3JOg5n1akXBh1K4hxzjK41aEMiA2XJ3STEkiWpVZKJyf5Wz17nFe3TSDW55ebC9OmweXPQkaQsS+4mIVIyuW/fTtanLwMwb36SzkotT48esHMnTJ0adCQpy5K78d2WLfD99ymY3N9/n0O3fkfjw38pe+GOZNaxozYTs7p7YCy5G9+l4jBIAMaOhTp1yM6pUf56qsmqenU47zytu9uQyEBYcje+i42USankXlQEb74J3buTdXI6CxdqlSKl9OgBq1fDF18EHUlKsuRufJeKY9yZPRt++AH69CErC3bsSMFh39266f3b+13Dx/jEkrvx3ZIlun5y3bpBR5JAY8dqaSI3t+yFO1JBw4Zw+ukwblzQkaQkS+7Gdyk3UsY5Te5du8LBB3PiidrHPuWSO8BFF8Fnn8GKFUFHknIsuRvfpVxynzdPp9/36QNoI8gTTiD1RszA7p8BY8cGG0cKsuRufPXzz/DddymW3MeO1UZhvXvvfmqvhTtSyfHH65t/442gI0k5ltyNrxYv1vuWLYONI6HGjoVOnbTmXCwrSysTGzcGGFdQLroIZszQC8wmYSy5G18tXKj3J54YbBwJs2wZfPnlnnJEsdhF1fnzA4gpaH366HWIN98MOpKUYsnd+GrhQl0zNWUW6Xj9db0vldxjqzKlZGkmK0vbgVrdPaEsuRtfLVqkv9fJvrrcbiNHau/2UusJNm2qs/FT8qKqiJZmpk2DDRuCjiZlWHI3vlq4MIVKMl9/DZ9/Dpdeus+3RKBNGx0VmJIuukin6FqvmYSx5G58U1ios1NT5mLqa69pFu/Xr8xvt2+v5fjCwgTHFQannqoN/W3UTMJYcje+WbpUT9ZSJrmPHAlnnAGNG5f57XbtYNs2WLAgwXGFQVqaXoeYNMl6vCeIJXfjm0WL9D4lyjLz52vW7t+/3E3at9f7lC3N9O8Pv/wC48cHHUlKsORufJNSwyBHjtSz04svLneTFi2gTh349NMExhUmHTvqleURI4KOJCVYcje+WbhQy6wHHxx0JD5zTpN7165wxBHlbpaWBm3bpnByT0vTs/fJk2Ht2qCjSXqW3I1vUmakzMcf6wWGAQMq3LR9e21vnpIXVUF/RoWFe+YDGN9Ycje+2LULvvoKTjop6EgS4OWXoVYt6Nu3wk3btdOyc+x6RMpp00a7qFlpxneW3I0vli2DrVv3TLtPWr/8oiWZiy+Oq2F97KJqypZmRPTsffp07ShnfGPJ3fgiNhPz5JODjcN3b76p3cAGDYpr8xNOgNq1U3jEDGhyd07nBRjfWHI3vpg7V6+ftWoVdCQ+e+UVOPpo6NIlrs3T07UyMWeOv2GFWosWkJOjPztbPNs3cSV3EekmIotFZImIDCnj+7eLyAIRmSsi00TkGO9DNVEyd642C6tVK+hIfPTddzBlClx5pf4li9Ppp+uZ+44dPsYWdlddpf9JPv886EiSVoX/I0UkHXgG6A60AgaISOnzsc+BHOdcNvA68JDXgZpomTs3Bert//qXXjkeOPCAXnb66bB9u46aSVkDBsBBB8GwYUFHkrTiOd04FVjinFvmnNsBjAR6l9zAOfe+c25r8cPZQBNvwzRRsmWLjgxM6uReVATPPw9nn33Ay0x16KD3s2b5EFdUHHqoNhMbPlx7MhjPxZPcGwPflnhcUPxcea4BJlUlKBNtsQUpkjq5T5qkSyvddNMBv7RxY2jSJMWTO2hpZv16W8TDJ/EkdynjuTKvgojI5UAO8HA5379eRPJEJG/NmjVNjywxAAAR/UlEQVTxR2kiJVZGTeqRMs8+C40aQa9elXp5hw4we7bHMUXN2WfrxWgrzfginuReADQt8bgJsKr0RiJyDvB7oJdzbntZO3LOveCcy3HO5TRo0KAy8ZoImDMH6tfX39uktGwZvPMOXH89VK9eqV106KAn/qtXexxblKSn6xDSKVPg228r3NwcmHiS+xyguYg0E5EaQH9gr7ZuItIWeB5N7D96H6aJkrw8OOUUna+SlJ5/XkfHXHddpXdx+ul6n/Klmdj8gH/+M9AwklGFyd05VwgMBiYDC4FRzrmvRGSoiMQ+kz4M1AFGi8gXImI9PVPUzz9r24GcnKAj8cnPP2si6t273L7t8WjXTpceTPnk3qwZdO8OL7yQ4mNDvVctno2ccxOBiaWeu6/E1+d4HJeJqC++0NGBp5wSdCQ+GTYM1q2D22+v0m4OOkgXJ/rwQ4/iirKbb9YEP2ZMXM3XTHxshqrxVGzmZVKeuRcWwqOPal/yTp2qvLsuXbTHzKZNVQ8t0s47T2e8/e1vQUeSVCy5G0/NmaPVikaNgo7EB2PGQH4+/L//58nuunTRTzkzZniyu+hKS4Pf/lZrVCnbUc17ltyNp+bMSdKzdufgoYe089cFF3iyyw4dtO7+wQee7C7aBg3SjmpPPx10JEnDkrvxzA8/wDffaNUi6UyYoA1h7rzzgPrI7E+tWnDaaZbcAV2ua9AgnbFqrYA9YcndeCZWXjjjjGDj8NyuXXDvvXDccdokzEOxuvvGjZ7uNpruuEN/1o8+GnQkScGSu/HMRx9BzZp7FqRIGmPH6jCg+++v9KSl8sTq7h995OluoykzU0fLvPCCrbHqAUvuxjMffaRlhho1go7EQ0VFcN99uhjsZZd5vvsOHfQP4pQpnu86moYM0bkENnKmyiy5G09s2qQnt0lXknnpJViwAIYO1enyHqtZU1usTJhg61YA0Lq19ut56inYvDnoaCLNkrvxxKxZWl5IquS+fj3cc4++qTgWv66s3FxtV/PNN74dIlruuUd/9nb2XiWW3I0npk3TcnSsV3lS+NOfdDbqU0/52ignN1fvJ07c/3Yp47TT9Oz9wQet9l4FltyNJyZPhs6doU6doCPxyLx5Oub6hht00VMfZWbqWrMTJvh6mGh54AEtyzz4YNCRRJYld1Nlq1bpsnrdugUdiUd27tSFJA47DP7854QcMjdX+8xs2ZKQw4XfSSfBFVdoaaagIOhoIsmSu6myyZP1PmmS+0MP6eDz556Dww9PyCF79tS/KZNsDbM9/vQnvZAzZEjQkUSSJXdTZZMnay+ZrKygI/HAp59qUunXDy6+OGGH7dwZjjwSXnstYYcMv8xM7eMzfDhMnx50NJFjyd1Uyc6dOkb7/POTYHGO9et1VMyRR+oyegmUnq6HnjDBRgDu5e67dUmvwYO1K6eJmyV3UyXvv685sU+foCOpoqIiGDhQ+5qMGpWwckxJl14K27bBW28l/NDhVasWPP64XuC2oZEHxJK7qZLXX9cRMuedF3QkVeCcLr7x1lvw2GN71sBLsI4dtV3yiBGBHD68+vSBHj3g97+Hr78OOprIsORuKq2wUNuuXHABZGQEHU0VPPKIjmW/7Tb9+B+QtDQdIDJpkjVG3IuI9pvJyNDOkUVFQUcUCZbcTaVNeWcXP/0E/c7bAMuX6zTL1athw4bo1Ef/+le9aHfppZrkA3bttZq7hg0LOpKQOeooLcvMmgUPPxx0NJEgLqCGFjk5OS4vLy+QY5sDtHatNo758ksd0L50KXz7LX1XPsp0dwYFNKEGO/d+TVoaNGyodYbYLJ2TTtLbCSd41hO90oqKdIjdI49oQ7BXXoFqcS0p7Ltzz9VWBEuX+tLOJrqc0z/Cb7wB770Hv/pV0BEFQkQ+dc5VuCSOJXezr3Xr9Epp7LZgwZ7vNWoELVqwpkErGr/xFDd3/oJH+8/RDlgi8Msvelu/XmsL332nWWrp0j2dsQ4+WOvaHTvqrUMHXYUnUdasgV//Gt59V5d3e/LJUGXR0aN1JOb48Z4t+pQ8Nm3S1dc3bdLFU5JyPcf9s+RuDkxBAYwbp0X0Dz/UM9vatXUAdpcuunbeySdDgwaAzgofMkQHMZx0Uhz7/+UXWLRIz/5nz4aZM2H+fE341appE/gzz9SzsU6d4JBDvH+PzulA8ptv1uTw7LNwzTXeH6eKdu7U9aKPOgr++98kGGLqtfnztf9MmzYwdaqeWKSQeJM7zrlAbu3bt3cmYGvXOvfMM86ddppzmvqca9nSuXvuce6//3Vux44yX7Z9u3NHHeVc165VPP6GDc5NmuTc3Xc716mTc9WrawwizrVt69zvfufcmDHOrVlTteMUFTn3zjvOnX667v+UU5ybO7eKwfvrmWc01A8+CDqSkBo9Wv+f9OnjXGFh0NEkFJDn4sixltxTzc6dzk2c6Fy/fs7VqKH/BbKznfvf/3Vu4cK4dvHKK/qyiRM9jm3rVufee8+5P/7RubPPdq5mzT1/dFq1cu7GG50bMcK5ggLndu3a/7527HBuxgzn7r3XuWOP1X00bercCy9EIhls3ercEUfoH9CK3mrKeuIJ/Xf9zW9S6ocUb3K3skyqWLQIXn4ZXn1VO30dfrjWnQcNgrZt497Nzp16bbRmTa2w+Foy2LED8vK0TDR9ui7SGuusVaeOXqht0kS/rlkTtm/Xckt+PixZoiN2RLSsdPXVcMklcNBBPgbsrSefhFtv1eH3PXsGHU1I3XWX9gK68UZ45pngL9QngNXcjQ5JfO01TeqzZ+tFw9xcTeg9e1ZqPbznnoObboK339Z5JQlVWKijdmbO1GGXy5frH6qtW/V20EFQty40barL4rVrp8scHXZYggP1xs6d2q9n1y4tMyfV8oVecU5bFDz4oM4w/sc/PF/nNmys5p6qCgudmzzZuQEDnMvI0I+trVs798gjzq1eXaVdr1njXP36zv3qVyn1KThQEyfqP+H99wcdSYjt2uXc0KH6g+rSperXaEKOOMsy4RjYa6pu8WIdq/2vf+nww0MP1ZEggwbpSBQP6ie33QYbN+ogExvBkRjdu8Pll8Nf/qItlQPqjBBuInDvvdCsmc4CO/VU7eFw2mlBRxao5C9QJbPvvtNp8x06aBnioYd0eNjo0TpT9OmndQijB5n41Vfh3//WT8CtW3sQu4nb00/rpYX+/eGHH4KOJsQuv3zPMN5OneD++/W6TYqymnvUrFoFY8ZoAp8xQ2uOWVlw5ZV6gdSHSR0ffwxnnaUnQu++G5qJnCklL0+nAbRsCR98kETLGfph40a45Rb9FNu8OTz6qF5jSpKPm/HW3O3MPeyKirSfxr33anmlcWP9j7thAwwdqqNg5s6FO+7wJbF/+qmWBho1gpEjLbEHJSdHr41//jl07Qo//RR0RCF28MFaopw4UUfP9OqlF9anTNkzSzoF2Jl72BQW6hjDGTPgo490+v+6dTrSpUMHzbR9+ugpnM9GjdKSfYMG+mk3M9P3Q5oKjBsHAwbo7NXhw60GX6GdO/Ui0YMPaqmyTRuty/fvH0jPfi94OloG6AYsBpYAQ8r4/kHAa8Xf/xjIrGifNlrG6ciWhQt1Ys5ddzl3zjnO1amzZ+JOZqZzAwc699przq1bl7CwVq7UOU7gXIcOzn3/fcIObeIwc6ZzxxzjXHq6czff7NyPPwYdUQRs2+bciy/qhD3Q2dA9ezr37LPOLVsWdHQHBK8mMYlIOvA1cC5QAMwBBjjnFpTY5iYg2zl3o4j0B/o45y7d335T5sx92zb4/ntYuVIn1ixdqvdLlsDChdpzBXRsbuvW2kjrjDO0p0uTJgkLs7BQPyT8+9860CAtDf7wB+2Ga+Orw2fjRp2/889/6n+dfv20F3znzhHvrZ8IX36pIwTGjNEJbwDHHKMNyU45RedHtGihv38hnBTl2SQmEekA/NE5d37x47sBnHN/LbHN5OJtZolINeB7oIHbz85DndyLinS24/5umzbpb9iGDXofu23YoMn8++/1Y+D69Xvvu1o1rW8cf7yOcGnbVhtytWzpaxbdtQt+/lnX59y4UfuE5efrXKC8PL1ounmzzgEaOFBL+Mcc41s4xiOLFulM1uHD9d+vZk0t1WRl6blC06a6JGzDhvpvW6tWqBpgBss57a38zjvaoW3OHJ0YF5ORAccdp2u4Nmy453bEEVCvnl7VrltX7+vU0UZ71avvufn0h8HL5N4X6Oacu7b48RXAac65wSW2mV+8TUHx46XF25R72aeyyf2lq6bzyH+OAsA5Ady+XyPgYo8o8X39buz7UNY2Jfax+xVS7nN7fZ2WhpM0/e1JS8cV35OehktLh2rV9D62vxIHr+zX8WxbVKSJvax/6mrVtKtjhw7aR7x7dzvzi6Kff9ZPXu++q9ffFyzQ58pSo4b+EcjI0Pyzv5tI0gwyiU9RIWzbrkMoY7fCnfrRtrCIvTNGHESA4h+isPuHef8NP3Dp02dUKsR4k3s8Yx/K+qct/Q7j2QYRuR64HuDoo4+O49D7qn/UQZzUcM3uQ8Z+diBlfi27NyqxrUjZzyNIGruTs6QXJ+r0NP06LR2ppvfUqI7UqKG/KTVqINXS997f7ve879cVff9Av47n+3Xr6q1ePb1v3FjnfBx1lI2ASQa1a+tov1gPml274NtvdeTsDz/oB8ktW7QKGOvWsG2b/sHftWv/t9RSrfhW1voCDnbshO3bNNnvLCxO+iVubhfsKv6hlvV18RW1Q4+pl5B3UpECoGmJx02AVeVsU1BcljkYWFd6R865F4AXQM/cKxNwrwdOo9cDlXmlMakjLU3LalZa85IANYpv4RdPUWgO0FxEmolIDaA/ML7UNuOBgcVf9wXe21+93RhjjL8qPHN3zhWKyGBgMpAOvOSc+0pEhqJDcsYDLwKvisgS9Iy9v59BG2OM2b+4qq3OuYnAxFLP3Vfi623AJd6GZowxprLCN4jTGGNMlVlyN8aYJGTJ3RhjkpAld2OMSUKW3I0xJgkF1vJXRNYAKyr58vpAsnS0tvcSPsnyPsDeS1hV5b0c45xrUNFGgSX3qhCRvHh6K0SBvZfwSZb3AfZewioR78XKMsYYk4QsuRtjTBKKanJ/IegAPGTvJXyS5X2AvZew8v29RLLmbowxZv+ieuZujDFmPyKb3EXkzyIyV0S+EJEpInJU0DFVlog8LCKLit/PWBE5JOiYKkNELhGRr0Rkl4hEclSDiHQTkcUiskREhgQdT2WJyEsi8mPxKmmRJSJNReR9EVlY/H/rd0HHVFkikiEin4jIl8Xv5U++Hi+qZRkRqeec21T89S1AK+fcjQGHVSkich7aA79QRB4EcM7dFXBYB0xEWgK7gOeBO5xzIV0kt2zxLAYfFSLyK2AL8C/n3ElBx1NZItIIaOSc+0xE6gKfAhdG9N9EgNrOuS0iUh2YAfzOOTfbj+NF9sw9ltiL1eaAFzcMD+fcFOdcYfHD2ehqV5HjnFvonFscdBxVcCqwxDm3zDm3AxgJ9A44pkpxzk2njNXQosY5t9o591nx15uBhUDjYKOqHKe2FD+sXnzzLW9FNrkDiMgDIvIt8Gvgvoq2j4irgUlBB5GiGgPflnhcQEQTSTISkUygLfBxsJFUnoiki8gXwI/Au845395LqJO7iEwVkfll3HoDOOd+75xrCgwHBgcb7f5V9F6Kt/k9UIi+n1CK531EWFwLvZvEE5E6wBjg1lKf2iPFOVfknGuDfjo/VUR8K5mFet1759w5cW76H2ACcL+P4VRJRe9FRAYCPYGuYV5/9gD+TaIonsXgTYIV16fHAMOdc28EHY8XnHMbROQDoBvgy0XvUJ+574+INC/xsBewKKhYqkpEugF3Ab2cc1uDjieFxbMYvEmg4ouQLwILnXOPBR1PVYhIg9hIOBGpCZyDj3kryqNlxgAnoKMzVgA3Oue+CzaqyileWPwgYG3xU7OjOPJHRPoAfwMaABuAL5xz5wcb1YERkVzgCfYsBv9AwCFVioiMALqg3Qd/AO53zr0YaFCVICKdgY+AeejvOsA9xes6R4qIZAOvoP+30oBRzrmhvh0vqsndGGNM+SJbljHGGFM+S+7GGJOELLkbY0wSsuRujDFJyJK7McYkIUvuxhiThCy5G2NMErLkbowxSej/A2QntWHSphnaAAAAAElFTkSuQmCC\n", "text/plain": [ - "0.43448275862068964" + "
" ] }, - "execution_count": 17, "metadata": {}, - "output_type": "execute_result" + "output_type": "display_data" } ], "source": [ - ".9*.0035/(.1*.041+.9*.0035)" + "x=np.linspace(-3,3,1000)\n", + "y1=lambda1*norm.pdf(x,-1,.5)+lambda2*norm.pdf(x,1,.5)\n", + "y2=wt1*norm.pdf(x,post_mean1,np.sqrt(post_var1))+wt2*norm.pdf(x,post_mean2,np.sqrt(post_var2))\n", + "fig,ax=plt.subplots(1)\n", + "ax.plot(x,y1,color='red',label='prior')\n", + "ax.plot(x,y2,color='blue',label='posterior')\n", + "ax.legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This isn't consistent with the solutions, which I don't understand. The key issue is the meaning of the statement that \"the variance of each observation is known to be 1\". How does one compute $p_{m}(\\{y_{i}\\})$?" ] }, { @@ -136,7 +173,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.6.4" + "version": "3.6.5" } }, "nbformat": 4, diff --git a/Useful Formulae.ipynb b/Useful Formulae.ipynb index 3fc65f4..5a3d891 100644 --- a/Useful Formulae.ipynb +++ b/Useful Formulae.ipynb @@ -136,7 +136,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.6.4" + "version": "3.6.5" } }, "nbformat": 4, From 025db3a873daa6111df00dc165ef2861e80e8fe7 Mon Sep 17 00:00:00 2001 From: Jeremy Teitelbaum Date: Sun, 29 Apr 2018 13:26:52 -0400 Subject: [PATCH 8/8] more on 5.9.8 --- BDA 5.9.8.ipynb | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/BDA 5.9.8.ipynb b/BDA 5.9.8.ipynb index 5228ce3..8ad5732 100644 --- a/BDA 5.9.8.ipynb +++ b/BDA 5.9.8.ipynb @@ -6,6 +6,10 @@ "source": [ "# Discrete Mixture Models\n", "\n", + " This solution differs from the one published here:\n", + "http://www.stat.columbia.edu/~gelman/book/solutions3.pdf\n", + "\n", + "\n", "Discrete mixture models: if $p_m(\\theta)$, for $m=1,\\ldots,M$ are conjugate prior densities for the sampling model $y|\\theta$, show that the class of finite mixture prior densities given by \n", "$$\n", "p(\\theta)=\\sum_{1}^{M} \\lambda_m p_m(\\theta)\n", @@ -52,7 +56,7 @@ }, { "cell_type": "code", - "execution_count": 23, + "execution_count": 29, "metadata": {}, "outputs": [], "source": [ @@ -68,7 +72,7 @@ }, { "cell_type": "code", - "execution_count": 24, + "execution_count": 30, "metadata": {}, "outputs": [], "source": [