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 +} diff --git a/BDA 5.9.6.ipynb b/BDA 5.9.6.ipynb new file mode 100644 index 0000000..e75cef4 --- /dev/null +++ b/BDA 5.9.6.ipynb @@ -0,0 +1,144 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Problem 5.9.6. Exchangeable Models\n", + "\n", + "1. In the divorce rate example of Section 5.2, set up a prior distribution for the values $y_1 \\ldots, y_8$ that allows for one low value (Utah) and one high value (Nevada), with independent and identical distributions for the other six values. This prior distribution should be exchangeable, because it is not known which of the eight states correspond to Utah and Nevada. \n", + "\n", + "2. Determine the posterior distribution for $y_8$ under this model given the observed values of $y_1, \\ldots, y_7$ given in the example. This posterior distribution should probably have two or three modes, corresponding to the possibilities that the missing state is Utah, Nevada, or one of the other six. \n", + "\n", + "3. Now consider the entire set of eight data points, including the value for $y_8$ given at the end of the example. Are these data consistent with the prior distribution you gave in part (1) above? In particular, did your prior distribution allow for the possibility that the actual data have an outlier (Nevada) at the high end, but no outlier at the low end?\n", + "\n", + "The states are Arizona, Colorado, Idaho, Montana, Nevada, New Mexico, Utah, and Wyoming.\n", + "\n", + "The rates from seven of these states are 5.8, 6.6,7.8,5.6,7.0,7.1,5.4 divorces per 1000 population per year. \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 135). CRC Press. Kindle Edition. " + ] + }, + { + "cell_type": "code", + "execution_count": 73, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "mean: 6.471428571428571 variance: 0.6877551020408161\n" + ] + } + ], + "source": [ + "from scipy.stats import norm\n", + "from itertools import permutations\n", + "import matplotlib.pyplot as plt\n", + "\n", + "rates=np.array([5.8,6.6,7.8,5.6,7.0,7.1,5.4])\n", + "print('mean:',rates.mean(),'variance:',rates.var())\n" + ] + }, + { + "cell_type": "code", + "execution_count": 74, + "metadata": {}, + "outputs": [], + "source": [ + "s=range(8)\n", + "def prior(x):\n", + " L=0\n", + " for i,j in permutations(s,2):\n", + " L0=1\n", + " L0=L0*norm.pdf(x[i],loc=8.7,scale=.8)\n", + " L0=L0*norm.pdf(x[j],loc=4.3,scale=.8)\n", + " for k in s:\n", + " if k!=i and k!=j:\n", + " L0=L0*norm.pdf(x[k],loc=6.5,scale=.8)\n", + " L=L+L0\n", + " return(L)\n", + "\n", + "def like(x):\n", + " return prior(np.append(x,rates))" + ] + }, + { + "cell_type": "code", + "execution_count": 77, + "metadata": {}, + "outputs": [], + "source": [ + "x=np.linspace(0.0,15.0,1000)\n", + "y=np.array([like(i) for i in x])" + ] + }, + { + "cell_type": "code", + "execution_count": 78, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[]" + ] + }, + "execution_count": 78, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAD8CAYAAABZ/vJZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl8lNd56PHfo9GGhBaQBNoAiSADYrWtYLzEiYMdg12D0zgJTus4qR3SxNS58b1t7dt708at27i9rdsmdmPHOCGuE+wSJ1EcsOMtXtnEIkBsFhIgsUlISGLT/tw/5hUdxIxmtL6jmef7+ejDzHnPOe8zgPToPee87xFVxRhjjOmvGLcDMMYYMzpZAjHGGDMglkCMMcYMiCUQY4wxA2IJxBhjzIBYAjHGGDMglkCMMcYMiCUQY4wxA2IJxBhjzIDEuh3AcMrMzNSCggK3wzDGmFFl69atp1Q1K1i9iE4gBQUFlJWVuR2GMcaMKiJyOJR6NoRljDFmQCyBGGOMGRBLIMYYYwbEEogxxpgBsQRijDFmQCyBGGOMGRBLIMYYYwYkou8DMWawLrR3UVp+lJYLnSyZk03+uCS3QzImbFgCMSaAU2fbWP7MRirrzgLwL68f4Ol7rubGK4LeoGtMVLAhLGP86OpWvv78Vmoaz/Pjr36c9/7iJqZkJLHyZ9uoaTzvdnjGhAVLIMb4sWbLEbYePs3ff3YON02fwKTxSTxzTwndCo++ssft8IwJC5ZAjOmltaOLf33jIxYUjucPr8q7WD45I4kVN07l9T0n2X7ktIsRGhMeLIEY08svtx+l/kwb31pUhIhccuy+GwpJGxPHqverXYrOmPBhCcSYXp7fcJjinFSu+1jGZceSE2L5/NX5vLr7BHVnWl2IzpjwYQnEGB/7TrSw53gLXyjJv+zqo8eXrplMZ7fyq+1HRzg6Y8KLJRBjfPxy21FiY4Q75uUGrDM1ayyz81L57a4TIxiZMeHHEogxDlWltPwYn5qeRcbYhD7r3jYnh/KaJlvSa6KaJRBjHLuPtnC8uZXFs3OC1r19jrfO+t3HhzssY8JWSAlERBaLyH4RqRSRh/0cTxCRF53jm0SkwOfYI075fhG5NVifIrJKRMpFZKeIrBWRsU75V0SkXkR2OF/3D+aDG9PbG3tPIgI3TQ9+p/mUjGRmZKfw9r76EYjMmPAUNIGIiAd4ElgCFAN3i0hxr2r3AadVdRrwBPC407YYWA7MAhYDT4mIJ0if31bVeao6FzgCrPQ5z4uqOt/5enZgH9kY/97cd5KrJo8LOnzV45NXZFF2uJFzbZ3DHJkx4SmUK5AFQKWqVqlqO7AGWNarzjJgtfN6LbBIvEtYlgFrVLVNVauBSqe/gH2qaguA034MoIP5gMaEoq6lld1HW1g0c0LIbT5RlEVHl7KxqmEYIzMmfIWSQPKAGp/3tU6Z3zqq2gk0Axl9tO2zTxH5MXACmAF836fe53yGtiaFELsxIfnwoDcJ3FgU+oMSSwrGkRgXw7sHbBjLRKdQEoi/xfC9rwoC1elvufeF6leBXGAv8EWn+DdAgTO09Qb/fcVzaSAiK0SkTETK6uvtG9uE5sODp0gbE8fMnNSQ2yTGeVg4NYP3K08NY2TGhK9QEkgt4Pvbfj5wLFAdEYkF0oDGPtoG7VNVu4AXgc857xtUtc05/CPgan/BquozqlqiqiVZWfbYbROaDw82sHDqeDwx/m8eDGRB4XgO1p+j4Wxb8MrGRJhQEsgWoEhECkUkHu+keGmvOqXAvc7ru4C3VFWd8uXOKq1CoAjYHKhP8ZoGF+dA7gD2Oe9911YuxXt1Ysyg1TSep/b0Ba77WGa/2y4oGA9A2WF7uKKJPkE3lFLVThFZCbwGeIDnVLVCRB4FylS1FFgFPC8ilXivPJY7bStE5CVgD9AJPOBcWRCgzxhgtYik4h3mKge+4YTyoIgsdfppBL4yJH8DJuptcOY//D37Kpg5+WnEx8awpbqRW2dlD3VoxoS1kHYkVNV1wLpeZd/xed0KfD5A28eAx0Lssxu4PkA/jwCPhBKvMf2xsaqBzLHxTJswtt9tE2I9zM9PZ8uhxmGIzJjwZneim6i3o6aJKyePC/jwxGA+XjiO3cda7H4QE3UsgZio1nyhg6pT55iXnzbgPj5eMJ6ubmX7kaYhjMyY8GcJxES13UebAZg3KX3AfVw5aRwA5bWWQEx0sQRiotqOGu8P/bl5A08gaUlxFGYmU15jCcREF0sgJqrtrG2iMDOZtKS4QfUzLz/NrkBM1LEEYqJaeU0zcwcx/9Fjbn46J1vaONFs29ya6GEJxEStky2tnGhpZW7+wIevesyb5E1CdhViooklEBO1euYs5k8a/BXIrNw0PDHCTksgJopYAjFRa2dtM54YoThn8AkkMc7D9IkplNc0D0FkxowOlkBM1CqvbWL6xBTGxHuGpL95k9LZWduE9zFwxkQ+SyAmKqkqO2ubL85dDIV5+Wm0tHZyqOH8kPVpTDizBGKi0uGG8zRf6BiSCfQes/O8yWjPsZYh69OYcGYJxESlntVS84YwgRRNHEtsjFBxzOZBTHSwBGKiUnlNM4lxMVwxsf9P4A0kIdbDtAlj2XPcrkBMdLAEYqJSeW0Ts3PTiPUM7bdAcW4qFTaEZaKEJRATdTq7uqk41jyk8x89ZuWmUX+mjbozdke6iXyWQEzUOXDyLK0d3UO6AqtHcU4qYBPpJjqElEBEZLGI7BeRShF52M/xBBF50Tm+SUQKfI494pTvF5Fbg/UpIqtEpFxEdorIWhEZG+wcxvTHcEyg9yjOdRKIzYOYKBA0gYiIB3gSWAIUA3eLSHGvavcBp1V1GvAE8LjTthjv/uizgMXAUyLiCdLnt1V1nqrOBY4AK/s6hzH9tbO2ibQxcUzJSBryvtPGxDFp/BibBzFRIZQrkAVApapWqWo7sAZY1qvOMmC183otsEi8+4MuA9aoapuqVgOVTn8B+1TVFgCn/RhAg5zDmH7Z4TyBd7j++xTnpLLXEoiJAqEkkDygxud9rVPmt46qdgLNQEYfbfvsU0R+DJwAZgDfD3IOY0J2ob2LAyfPDMvwVY9ZuWlUN5yzPdJNxAslgfj7Na33w34C1elvufeF6leBXGAv8MV+xIGIrBCRMhEpq6+v99PERLOKY810deugtrANpjgnFVXYd8KuQkxkCyWB1AKTfN7nA8cC1RGRWCANaOyjbdA+VbULeBH4XJBz0KvdM6paoqolWVlZIXw8E03Ka5090IdgE6lAZuV5J9JtHsREulASyBagSEQKRSQe76R4aa86pcC9zuu7gLfU+0jSUmC5s4KqECgCNgfqU7ymwcU5kDuAfUHOYUzIymuayElLZEJq4rCdIzs1kXFJcVQctQRiIltssAqq2ikiK4HXAA/wnKpWiMijQJmqlgKrgOdFpBLvVcFyp22FiLwE7AE6gQecKwsC9BkDrBaRVLxDVuXAN5xQ/J7DmP7YWds0JFvY9kVEKM5NZa8NYZkIFzSBAKjqOmBdr7Lv+LxuBT4foO1jwGMh9tkNXB+gn4DnMCYUTefbOdRwns+XTApeeZCKc1JZveEwnV3dQ/64FGPChf3PNlFjpzP/MX8YJ9B7FOem0t7ZTdWpc8N+LmPcYgnERI2e/cp79u0YTj3b5O61O9JNBLMEYqLGjppmpmYlkzYmbtjPNTUrmfjYGHsmlololkBMVFBVymubhvUGQl9xnhimT0yxZ2KZiGYJxESFEy2t1J9pG9b7P3qbmZPCnmMt2GpzE6ksgZioUF7jnUCfOwIT6D2Kc1JpONdO3Zm2ETunMSPJEoiJCjtrm4iNkYv7dYyE4lzv1Y7Ng5hIZQnERIXy2iZm5KSQGOcZsXPOyEkBbG8QE7ksgZiI192t7Kwdni1s+5KaGMfk8UmWQEzEsgRiIl7VqbOcae0ckRsIe7O9QUwkswRiIt72I94bCK90IYHMzEm1vUFMxLIEYiLejpomUhJi+VjW2BE/d3Fuz94gZ0b83MYMN0sgJuJtP9LEvEnpxMSM/A7IxbneVV82D2IikSUQE9EutHex/+QZV+Y/AHLTEkkbE2dLeU1EsgRiItquo94tbK+c7E4CEfHee2JXICYSWQIxEW1HzWlgZB7hHkhxbir7T7TQ1W2PNDGRxRKIiWjbjzQxeXwSGWMTXIthZk4qrR3dVNveICbCWAIxEW1HTZOrVx/Axcen2DCWiTQhJRARWSwi+0WkUkQe9nM8QURedI5vEpECn2OPOOX7ReTWYH2KyAtO+W4ReU5E4pzyT4lIs4jscL6+gzF9ON58gePNra4nkGkTxhLnEZtINxEnaAIREQ/wJLAEKAbuFpHiXtXuA06r6jTgCeBxp20xsByYBSwGnhIRT5A+XwBmAHOAMcD9Pud5T1XnO1+PDuQDm+ixuboRgAWF412NIz42hqIJtjeIiTyhXIEsACpVtUpV24E1wLJedZYBq53Xa4FFIiJO+RpVbVPVaqDS6S9gn6q6Th3AZiB/cB/RRKvN1Y2kJMQycwSfwBtIcW6qXYGYiBNKAskDanze1zplfuuoaifQDGT00TZon87Q1T3Aqz7F14pIuYisF5FZ/oIVkRUiUiYiZfX19SF8PBOpNlc3UlIwDo8LNxD2NjMnlVNn26g70+p2KMYMmVASiL/vvt7rEQPV6W+5r6eAd1X1Pef9NmCKqs4Dvg/8yl+wqvqMqpaoaklWVpa/KiYKNJxt46O6sywozHA7FOC/J9L3HrdHmpjIEUoCqQUm+bzPB44FqiMisUAa0NhH2z77FJG/BrKAh3rKVLVFVc86r9cBcSKSGUL8JgptOdQz/zHO5Ui8Lq7EsmEsE0FCSSBbgCIRKRSReLyT4qW96pQC9zqv7wLecuYwSoHlziqtQqAI77xGwD5F5H7gVuBuVe3uOYGIZDvzKojIAif2hoF8aBP5NlU3khAbw5w8d1dg9UhLiiMvfYxNpJuIEhusgqp2ishK4DXAAzynqhUi8ihQpqqlwCrgeRGpxHvlsdxpWyEiLwF7gE7gAVXtAvDXp3PKHwKHgQ1OvnjZWXF1F/ANEekELgDLnSRlzGU2Vzdy1eRxxMeGz61O3on0ZrfDMGbIBE0gcHHIaF2vsu/4vG4FPh+g7WPAY6H06ZT7jUlVfwD8IJR4TXRrONvGnuMtfPvmK9wO5RLFOam8sfck59s7SYoP6VvPmLAWPr+eGTNE3q88hSrceEV4LaKYk5eGKuw+asNYJjJYAjER550D9YxLimNOXprboVxinnNHfHlNk8uRGDM0LIGYiNLdrbx74BQ3FGWFxf0fvrJSEshLH8MOSyAmQlgCMRFl74kWTp1t45NhNnzVY/7kdEsgJmJYAjER5e19dQDcWBSetwjNz0/naNMFuyPdRARLICairN99gisnpzMhNdHtUPyaP7lnHsSW85rRzxKIiRhHGs5TcayF22bnuB1KQLNz0/DEiE2km4hgCcREjPW7jwOweHa2y5EENibew/SJKTYPYiKCJRATEVSVX24/ytz8NCaNT3I7nD7Nn5xOeW0T3bZHuhnlLIGYiLD7aAv7Tpzh8yWTgld22fz8dM60dlJ16qzboRgzKJZATER4sewICbExLJ2X63YoQV1d4H1C8Obq0y5HYszg2AN5zIg7dbaNH/7+IG/uq6PxXDuTxyexeHY291w7hdTEuH7313y+g19tP8Ztc3JIG9P/9iNtamYymWMT2FTdwJeumex2OMYMmF2BmBG1ubqRW/7lHX7y4SGmZiZzx7wcEmJj+KfX9nPjP77NKzt7bzUT3OoNhzjb1smKG6cOfcDDQES4Zup4NlU1Yg+UNqOZXYGYEbOrtpkvP7eJ3PQxvPT1aymamHLJsf/z692s/Nl23t5Xz9/eOSukJ9aeae3gxx9Us2jGhLDY+zxUCwvH89udx6lpvMDkjPCe9DcmELsCMSOi8Vw7X3++jIzkhMuSB8Cc/DR+8afX8q1FRby8vZY7n/yAyrrgk8xPvP4RTRc6+NbNRcMV+rC4Zqp3q92N1bYnmhm9LIGYEfHYb/dSf7aNH/7x1WSOTfBbJ9YTw7dvuYKf/skCTp1tZ9kP3qe0PPCQ1oaDDfzkw2ruXjCZufnhsfNgqKZljWVcUhybqhrdDsWYAbMEYobd5upGfrGtlhU3TmVOfvBHrH+iKIvfPngDM3JSefDn2/m/v9pNW2fXJXUqjjXzwM+2UZiZzP++beZwhT5sYmKEBYXj2XzIrkDM6BVSAhGRxSKyX0QqReRhP8cTRORF5/gmESnwOfaIU75fRG4N1qeIvOCU7xaR50QkzikXEfl3p/5OEblqMB/cjAxV5Xvr95KTlsjKm0IfZspJG8OaFQv52icKeX7jYZb863s8/c5BXt19gn9Yv5c/fOpDEmNjWHXvxxmbMDqn8q6dmkFN4wWONJx3OxRjBiRoAhERD/AksAQoBu4WkeJe1e4DTqvqNOAJ4HGnbTHe/dFnAYuBp0TEE6TPF4AZwBxgDHC/U74EKHK+VgD/MZAPbEbWB5UNbDvSxDdvmsaYeE+/2sZ5Yvir24tZdW8JKWPi+If1+/jT/9zKs+9Vc/PMifzygespyEwepsiH3yenTwDg9wfqXI7EmIEJ5Ve3BUClqlYBiMgaYBmwx6fOMuBvnNdrgR+IiDjla1S1DagWkUqnPwL16eyVjlO+Gcj3OcdP1bvucaOIpItIjqoe7++HNiPnybcryU5N5Asl+cErB7Bo5kQWzZxIXUsrdWfayB83hvSk+CGM0h2FmckUZCTx9r46vnxtgdvhGNNvoQxh5QE1Pu9rnTK/dVS1E2gGMvpoG7RPZ+jqHuDVfsRhwshHJ8+woaqBe68rICG2f1cf/kxITWR2XlpEJI8en5o+gQ1VDbR2dAWvbEyYCSWB+NsXtPfdT4Hq9Lfc11PAu6r6Xj/iQERWiEiZiJTV19f7aWJGygubjhDviRnU1Uek+9T0LFo7utlYZZPpZvQJJYHUAr5PqMsHeq+tvFhHRGKBNKCxj7Z99ikifw1kAQ/1Mw5U9RlVLVHVkqys8NzWNBqcb+/kF1trWTInm4wAy3YNLJyaQWJcDG/sPel2KMb0WygJZAtQJCKFIhKPd1K8tFedUuBe5/VdwFvOXEUpsNxZpVWIdwJ8c199isj9wK3A3ara3escX3ZWYy0Emm3+I3yt33WCM22d/NE1U9wOJawlxnn49IwJvLr7BJ1d3cEbGBNGgiYQZ05jJfAasBd4SVUrRORREVnqVFsFZDiT5A8BDzttK4CX8E64vwo8oKpdgfp0+vohMBHYICI7ROQ7Tvk6oAqoBH4EfHNwH90Mp9/sPEZe+hg+7jx51gR2x9xcTp1tZ6PdVGhGmZAW0Dsro9b1KvuOz+tW4PMB2j4GPBZKn06535icK5oHQonXuKvxXDvvf3SK+z8xFe9iPNOXm2ZMIDnew2/Kj3FDUabb4RgTMrsT3Qy5V3efoLNbuWNe+O5NHk4S4zx8ZlY2r1acuOyOe2PCmSUQM+R+U36MqVnJFI+ip+O67bNX5tF8oYP1u064HYoxIbMEYoZU3ZlWNlY38Adzc234qh9umJZJQUYSP91wyO1QjAmZJRAzpH6/rx5VWDI72+1QRpWYGOGPF05h25EmdtU2ux2OMSGxBGKG1Jv7TpKblsiM7JTglc0lvvDxSaQmxvJvbx5wOxRjQmIJxAyZts4u3vvoFJ+eOcGGrwYgNTGOFTdO5Y29dWw/ctrtcIwJyhKIGTKbqho5397FohkT3Q5l1PrK9YVkjo3nO7+usBsLTdizBGKGzFv76kiMi+Haj2W4HcqoNTYhlu8unc2uo8388J2DbodjTJ8sgZghoaq8ue8kN0zLJDFu8E/ejWa3z83hjnm5/PPrB1i/y57WY8KXJRAzJCrrzlLTeIFP2/DVkPinu+Zy5aR0/uzn23nu/Wq6ui978LQxrhude4GasPNB5SkAPmGP4hgSiXEefvInC3joxXIefWUPL2w6zO1zci7uwNh4rp36s23Ut7RRd6aN8+2djE9OYHZeKrfPyaFooq2CM8PPEogZEhurGskfN4ZJ45PcDiVipCbG8cw9V/PbXcd57oNqfvB2Jb4XIvGeGLJSEshKSSA5wcORxnO8ue8k//rGR9w+J4dHl82yR+mbYWUJxAxad7eysbqBW2ba8NVQi4kR7piXyx3zcjnX1kn9mTa6VclITiB1TOxly6VPnW3jpx8e4ofvVrH5UCOrv7qA4lx7pIwZHjYHYgZt/8kzNJ3vYOFUW301nJITYinITGZq1ljSkuL83muTOTaBhz4znV9983riYoS7f7SRvcdbXIjWRANLIGbQNhz0bse60Jbvho3i3FRe/Pq1jInz8LWflnH6XLvbIZkIZAnEDNrGqgYmj08iL32M26EYH5PGJ/H0PVdTd6aNP19bjndLHWOGjiUQMyjd3cqm6kYWTh3vdijGj3mT0vmLW6fzxt46fr3jmNvhmAgTUgIRkcUisl9EKkXkYT/HE0TkRef4JhEp8Dn2iFO+X0RuDdaniKx0ylREMn3KPyUizc42t75b3RoX7T3RQvOFDrv7PIx99fpCrpqczqOv7OFMa4fb4ZgIEjSBiIgHeBJYAhQDd4tIca9q9wGnVXUa8ATwuNO2GFgOzAIWA0+JiCdInx8ANwOH/YTznqrOd74e7d9HNcOhZx9vm0APX54Y4btLZ9N4rp2n36lyOxwTQUK5AlkAVKpqlaq2A2uAZb3qLANWO6/XAovEu0RkGbBGVdtUtRqodPoL2KeqblfVQ4P8XGaEbDjYQEFGEjlpNv8Rzubkp7Fsfi7Pvl/F8eYLbodjIkQoCSQPqPF5X+uU+a2jqp1AM5DRR9tQ+vTnWhEpF5H1IjIrhPpmGHV1K5urG+zqY5T4X5+ZTmeX2lWIGTKhJBB/Gzv0Xs4RqE5/y/uyDZiiqvOA7wO/8ldJRFaISJmIlNXX1wfp0gzG3uMttLR22vzHKDFpfBJ3XpnHmi1HaLRlvWYIhJJAaoFJPu/zgd7LOS7WEZFYIA1o7KNtKH1eQlVbVPWs83odEOc7ye5T7xlVLVHVkqysrOCfzgzYxirv/R/XFFoCGS3+9JNTae3oZvWHh9wOxUSAUBLIFqBIRApFJB7vpHhprzqlwL3O67uAt9S76LwUWO6s0ioEioDNIfZ5CRHJduZVEJEFTuwNoXxIMzw2HGygMDOZ7LREt0MxIZo2IYVbiify0w2HaO3ocjscM8oFTSDOnMZK4DVgL/CSqlaIyKMistSptgrIEJFK4CHgYadtBfASsAd4FXhAVbsC9QkgIg+KSC3eq5KdIvKsc467gN0iUg78O7Bc7c4o13jnPxpt/mMU+sp1BZw+38H63bbXiBkcieSfwSUlJVpWVuZ2GBFpV20zd/zgff5t+XyWzQ9l/YMJF6rKon9+h3HJ8fziG9e5HY4JQyKyVVVLgtWzO9HNgGyo8u7/ca1dgYw6IsKXrpnM1sOn7UGLZlAsgZgB2VjVyNSsZCak2vzHaHTX1fkkxMbwwiZ/9+saExpLIKbfOru6bf5jlEtPimfJ7Gx+U36ctk6bTDcDYwnE9FvFsRbOtnXa8NUo99mr8mm+0MHb++rcDsWMUpZATL9dvP/DnsA7ql3/sQyyUhJ4edtRt0Mxo5QlENNvG6oamDZhLBNSbP5jNIv1xHDn/Fze3l9nG06ZAbEEYvqls6ubLbb/R8T47JX5dHQpr+yye0JM/1kCMf2y62gz59q7bAI9QhTnpjIjO4Vfbqt1OxQzClkCMf1i+39Enjvm5bLtSBPHmuwx76Z/LIGYftlQ1UDRhLFkjk1wOxQzRG6fkwPAOhvGMv1kCcSErKOrm7JDjfb49ghTkJnMrNxUSyCm3yyBmJDtOtrMeZv/iEi3zcmxYSzTb5ZATMg2HOzZ/8NWYEUaG8YyA2EJxIRsY1UD0yemkGHzHxGnIDOZ4hwbxjL9YwnEhKS9s5uyQ6ft/o8IdvtcG8Yy/WMJxIRk19EmLnR02QR6BLvNGcZav/uEy5GY0cISiAlJz/zHAtv/PGIVOsNY620Yy4TIEogJyYaqBmbmpDI+Od7tUMwwWjI7m7LDpznR3Op2KGYUCCmBiMhiEdkvIpUi8rCf4wki8qJzfJOIFPgce8Qp3y8itwbrU0RWOmUqIpk+5SIi/+4c2ykiVw30Q5v+aevsouzQaXt8exRY4gxjvVZhw1gmuKAJREQ8wJPAEqAYuFtEintVuw84rarTgCeAx522xcByYBawGHhKRDxB+vwAuBnovVXaEqDI+VoB/Ef/PqoZqO1Hmmjr7OY6m/+IeNMmjOWKiWNZv9uGsUxwoVyBLAAqVbVKVduBNcCyXnWWAaud12uBRSIiTvkaVW1T1Wqg0ukvYJ+qul1VD/mJYxnwU/XaCKSLSE5/PqwZmA8PNhAjsMBWYEWFJbNz2FzdSP2ZNrdDMWEulASSB9T4vK91yvzWUdVOoBnI6KNtKH0OJA5EZIWIlIlIWX19fZAuTSg2HDzFnPx0UhPj3A7FjIAlc7LpVvjdHhvGMn0LJYGInzINsU5/ywcbB6r6jKqWqGpJVlZWkC5NMOfbO9lR02TzH1Fk+sQUpmYms36XJRDTt1ASSC0wyed9PnAsUB0RiQXSgMY+2obS50DiMEOs7NBpOrrU5j+iiIiwZE42G6oabKdC06dQEsgWoEhECkUkHu+keGmvOqXAvc7ru4C3VFWd8uXOKq1CvBPgm0Pss7dS4MvOaqyFQLOq2kzfMPvwYANxHqGkYJzboZgRtGR2Dl3dyut7TrodigljQROIM6exEngN2Au8pKoVIvKoiCx1qq0CMkSkEngIeNhpWwG8BOwBXgUeUNWuQH0CiMiDIlKL9wpjp4g865xjHVCFdyL+R8A3B/3pTVAbqhqYPymdpPhYt0MxI2hWbiqTxo9hna3GMn0Q74VCZCopKdGysjK3wxi1Wlo7mP/d37Hy00U8dMsVbodjRtg/rNvLcx9UU/Z/biFtjC2giCYislVVS4LVszvRTUAbDzbQrdj8R5RaMieHji7lzb02jGX8swRiAnrnQD3J8R6ummzzH9FoXn4auWmJrLPVWCYASyDGL1XlnQP1XDcQdaVpAAAP40lEQVQtk/hY+28SjUSExbNzePejes62dbodjglD9pPB+FV96hy1py9w4xV2L000u21ONu2d3by1r87tUEwYsgRi/Hr3gPcu/k8WWQKJZldNHseElAR7xLvxyxKI8eudA/UUZiYzOSPJ7VCMi2JihMWzs3l7fx3n220Yy1zKEoi5TGtHFxurGvmkDV8ZvDcVtnZ0885+e7acuZQlEHOZskOnudDRxY1XZAavbCLegsLxZCTHs862ujW9WAIxl3nnQB3xnhgW2gMUDeCJET4zK5u39p6ktaPL7XBMGLEEYi7z5t46rpk63h5fYi66bU4259q7Li6uMAYsgZheKuvOUnXqHJ8pnuh2KCaMLJyaQXpSHK/aMJbxYQnEXKJnE6GbLYEYH3GeGG6ZOZHX956krdOGsYyXJRBzidf3nGRufho5aWPcDsWEmdvm5nCmtdNWY5mLLIGYi+paWtl+pMmGr4xfn5iWSebYBH6xrdbtUEyYsARiLnpjr/dxFbcUZ7sciQlHsZ4Y7pyfy1v76mynQgNYAjE+Xqs4wZSMJK6YONbtUEyY+tzV+XR0KaXltpu0sQRiHKfPtfNB5SkWz85GRNwOx4SpmTmpFOek2jCWAUJMICKyWET2i0iliDzs53iCiLzoHN8kIgU+xx5xyveLyK3B+nT2Sd8kIh85fcY75V8RkXoR2eF83T+YD24utW73cTq7lWXz8twOxYS5P7wqj521zXx08ozboRiXBU0gIuIBngSWAMXA3SJS3KvafcBpVZ0GPAE87rQtBpYDs4DFwFMi4gnS5+PAE6paBJx2+u7xoqrOd76exQyZX+84xrQJY5mZk+J2KCbMLZufhydGWLvVrkKiXShXIAuASlWtUtV2YA2wrFedZcBq5/VaYJF4x0GWAWtUtU1Vq4FKpz+/fTptPu30gdPnnQP/eCYUx5ousOVQI0vn5drwlQkqKyWBRTMm8F9ba+2ekCgXSgLJA2p83tc6ZX7rqGon0Axk9NE2UHkG0OT04e9cnxORnSKyVkQmhRC7CcErO4+hCkvn5bodihkl7rl2Co3n2lln+4REtVASiL9fSTXEOkNVDvAboEBV5wJv8N9XPJcGIrJCRMpEpKy+3m54CkZV+a+yWuZPSqcgM9ntcMwocf3HMpmamczzGw67HYpxUSgJpBbw/W0/H+i9hu9iHRGJBdKAxj7aBio/BaQ7fVxyLlVtUNU2p/xHwNX+glXVZ1S1RFVLsrJsP4tgth05zUd1Z1n+cbugM6GLiRH+aOEUth1pYvfRZrfDMS4JJYFsAYqc1VHxeCfFS3vVKQXudV7fBbylquqUL3dWaRUCRcDmQH06bd52+sDp89cAIpLjc76lwN7+fVTjz88315Ac7+EOG74y/XTXVfkkxsXwnxvtKiRaBU0gznzESuA1vD+0X1LVChF5VESWOtVWARkiUgk8BDzstK0AXgL2AK8CD6hqV6A+nb7+EnjI6SvD6RvgQRGpEJFy4EHgK4P76Kb5Qgev7DzG0vl5JCfYo9tN/6QlxfHZK/N5eftR6s60uh2OcYF4f+mPTCUlJVpWVuZ2GGFr9YeH+OvSCkpXXs/c/HS3wzGj0KFT5/j0P/+er904lUeWzHQ7HDNERGSrqpYEq2d3okeprm5l1fvVzJ+Uzpy8NLfDMaNUQWYyt83J4YWNR2i+0OF2OGaEWQKJUr+rOMGRxvN8/capdu+HGZRvfmoaZ9s6eX7DIbdDMSPMEkgUUlWefreKKRlJfGaWPXnXDE5xbio3Tc9i1fvVnGm1q5BoYgkkCm2qbmRHTRP331CIJ8auPszgPXTLdE6f7+Dpd6rcDsWMIEsgUUZV+afX9jMxNYHPl9i9H2ZozMlPY+m8XJ59v4oTzbYiK1pYAokyb+2rY+vh0zy4qIjEOI/b4ZgI8ue3TqerW/nn3+13OxQzQiyBRJGOrm4ef3UfUzKS+IJdfZghNml8En9yfSH/tbWWzdWNbodjRoAlkCiy+sNDHDh5lv9920ziPPZPb4bet24uIi99DI+8vNOe1BsF7KdIlDjR3MoTrx/g0zMm8JniiW6HYyJUUnwsf/fZ2RysP8f336x0OxwzzCyBRIHubuXP15bTrfA3d8yy+z7MsLpp+gTuujqfJ39fyQeVp9wOxwwjSyBR4LkPqnnvo1P83z8oZnJGktvhmCjw6LJZTM1M5ltrdthzsiKYJZAIt+FgA99bv4+bZ07k7gU2cW5GRlJ8LE/90dWcbevg/tVlnGvrDN7IjDqWQCJYVf1ZvvHCVgoyk/mXL86zoSszoqZnp/Dkl65i99FmHvjZNjq6ut0OyQwxSyARqvrUOe7+0UY8Iqy6t4TUxDi3QzJRaNHMifzdnXP4/f56vvbTMi6028qsSGIJJAKV1zTxhac30NGl/OxrC5mSYVvVGvd86ZrJ/P1n5/DOgXr+eNUm6lpsTiRSWAKJIKrKms1H+MLTG4j3xLBmxUKmZ6e4HZYxfOmayTz5pauoONbMbf/+Hu8eqHc7JDMELIFEiKr6s3z1J1t4+OVdXDV5HL9eeT1XTLTkYcLHbXNyKF15A+lJ8Xz5uc08+PPtHG++4HZYZhBCSiAislhE9otIpYg87Od4goi86BzfJCIFPsceccr3i8itwfp09knfJCIfOX3GBztHNNt9tJm/WFvOLU+8y+bqRr67dBYv3H8NmWMT3A7NmMtcMTGFV/7sBh5cVMSrFSf45D/+nkde3sne4y1uh2YGIOiWtiLiAQ4AtwC1wBbgblXd41Pnm8BcVf1TEVkOfFZVvygixcDPgQVALvAGcIXTzG+fIvIS8LKqrhGRHwLlqvofgc7RV+yRuKVtZ1c3u4428+6BU7y+9wS7j7aQEBvD3Qsm88BN08hKscRhRoeaxvM8/e5BXiqrpb2zm+kTU7h1djbXfSyDKyenkxBrD/t0S6hb2oaSQK4F/kZVb3XePwKgqv/gU+c1p84GEYkFTgBZwMO+dXvqOc0u6xP4HlAPZKtqp++5A51D+/gAoymBdHcrFzq6ON/exfn2Ts61dVF/to2Tza2cbGml5vR59p04w/4TZ2jr7EYE5uan84dX5nHn/DzSkmyVlRmdGs+189udx/jVjmNsP3KaboV4TwxTs5IpmphCYWYy2amJTExNICslgbEJsSQ7X0lxHmJsT5shF2oCiQ2hrzygxud9LXBNoDrOD/5mIMMp39irbZ7z2l+fGUCTqnb6qR/oHEP+rITf76/jb1/ZgwIoqPeczp+gqPdPJ3X5PUbPcd/3PvWc1yh0OsmjL5lj45mencI9C6cwd1I6N0zLZHxy/FB/dGNG3PjkeO65toB7ri2g+UIHm6sb2Xr4NAdOnmHb4dO8svMYff2eGx8bQ2yM4IkR50+f9x6hd3rpfT/UZelH+nwbvH2Y+OLHJ3H/J6YO6zlCSSD+/n56/3MGqhOo3N/cS1/1Q40DEVkBrACYPHmynybBpSTGMSM7FcR7UhFx/rz0vfe4+JT7vHcq+D3Gf/8nFIHYGCEpPpakeI/z5X2dmZJAdmoiE1IT7HLeRIW0MXHcUjyRW3we+NnR1c2ps22cbGmj/kwb59o6Odfeybm2Ts62ddHW2UV3t9LZrXT1/Nnl/bOz+9KbF3snot4/QHoPaFz2A+ay9n2P4LhpJOZBQ0kgtYDvMzDygWMB6tQ6w0tpQGOQtv7KTwHpIhLrXIX41g90jkuo6jPAM+Adwgrh813m6injuHrKuIE0NcYMsThPDDlpY8hJG+N2KKaXUFZhbQGKnNVR8cByoLRXnVLgXuf1XcBbztxEKbDcWUFVCBQBmwP16bR52+kDp89fBzmHMcYYFwS9AnHmG1YCrwEe4DlVrRCRR4EyVS0FVgHPi0gl3quC5U7bCmdV1R6gE3hAVbsA/PXpnPIvgTUi8nfAdqdvAp3DGGOMO4KuwhrNRtMqLGOMCRehrsKyO9GNMcYMiCUQY4wxA2IJxBhjzIBYAjHGGDMglkCMMcYMSESvwhKReuDwAJtnMgyPSRliFuPghXt8EP4xhnt8EP4xhlt8U1Q1K1iliE4ggyEiZaEsY3OTxTh44R4fhH+M4R4fhH+M4R5fIDaEZYwxZkAsgRhjjBkQSyCBPeN2ACGwGAcv3OOD8I8x3OOD8I8x3OPzy+ZAjDHGDIhdgRhjjBkQSyB+iMhiEdkvIpUi8rDb8fQmIpNE5G0R2SsiFSLyLbdj8kdEPCKyXURecTsWf0QkXUTWisg+5+/yWrdj8iUi33b+fXeLyM9FJDEMYnpOROpEZLdP2XgReV1EPnL+dHUznQAx/pPz77xTRH4pIunhFJ/Psf8lIioimW7E1l+WQHoREQ/wJLAEKAbuFpFid6O6TCfwP1V1JrAQeCAMYwT4FrDX7SD68G/Aq6o6A5hHGMUqInnAg0CJqs7Gu+1BOGxh8BNgca+yh4E3VbUIeNN576afcHmMrwOzVXUucAB4ZKSD8vETLo8PEZkE3AIcGemABsoSyOUWAJWqWqWq7cAaYJnLMV1CVY+r6jbn9Rm8P/jy+m41skQkH7gdeNbtWPwRkVTgRpz9ZlS1XVWb3I3qMrHAGGcHziQu3wl0xKnqu1y+E+gyYLXzejVw54gG1Yu/GFX1d84upwAb8e526ooAf4cATwB/gZ+ddMOVJZDL5QE1Pu9rCbMfzr5EpAC4EtjkbiSX+Ve83wzdwSq6ZCpQD/zYGWZ7VkSS3Q6qh6oeBf4f3t9GjwPNqvo7d6MKaKKqHgfvLzfABJfjCeZPgPVuB+FLRJYCR1W13O1Y+sMSyOXET1lY/kYgImOBXwD/Q1Vb3I6nh4j8AVCnqlvdjqUPscBVwH+o6pXAOdwfernImUdYBhQCuUCyiPyxu1GNfiLyV3iHgF9wO5YeIpIE/BXwHbdj6S9LIJerBSb5vM8nDIYOehOROLzJ4wVVfdnteHq5HlgqIofwDgF+WkT+092QLlML1Kpqz5XbWrwJJVzcDFSrar2qdgAvA9e5HFMgJ0UkB8D5s87lePwSkXuBPwD+SMPr/oWP4f1Fodz5nskHtolItqtRhcASyOW2AEUiUigi8XgnLktdjukSIiJ4x+73quq/uB1Pb6r6iKrmq2oB3r+/t1Q1rH57VtUTQI2ITHeKFgF7XAyptyPAQhFJcv69FxFGk/y9lAL3Oq/vBX7tYix+ichi4C+Bpap63u14fKnqLlWdoKoFzvdMLXCV8380rFkC6cWZaFsJvIb3G/YlVa1wN6rLXA/cg/c3+x3O121uBzUK/RnwgojsBOYDf+9yPBc5V0ZrgW3ALrzfq67frSwiPwc2ANNFpFZE7gO+B9wiIh/hXUX0vTCM8QdACvC68/3ywzCLb1SyO9GNMcYMiF2BGGOMGRBLIMYYYwbEEogxxpgBsQRijDFmQCyBGGOMGRBLIMYYYwbEEogxxpgBsQRijDFmQP4/PoM6VU2KnswAAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "plt.plot(x,y)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "I have a lot of questions about this problem, starting with: have I computed the right thing? I am not sure.\n", + "As far as the last part, No, my prior distribution did not allow for a high value without a low value; it expected one of each." + ] + }, + { + "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 5.9.8.ipynb b/BDA 5.9.8.ipynb new file mode 100644 index 0000000..d3ea631 --- /dev/null +++ b/BDA 5.9.8.ipynb @@ -0,0 +1,144 @@ +{ + "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", + "$$\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)" + ] + }, + { + "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/Useful Formulae.ipynb b/Useful Formulae.ipynb new file mode 100644 index 0000000..3fc65f4 --- /dev/null +++ b/Useful Formulae.ipynb @@ -0,0 +1,144 @@ +{ + "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", + "\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": 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\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": { + "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 +}