From f3398184a9ad87be08f4e2be8556f1504ea76665 Mon Sep 17 00:00:00 2001 From: jeremyteitelbaum Date: Thu, 26 Apr 2018 11:50:35 -0400 Subject: [PATCH 1/6] Problem 5.9.6 added, but is the solution right? --- BDA 5.9.6.ipynb | 144 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 144 insertions(+) create mode 100644 BDA 5.9.6.ipynb 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 +} From 770e7f432c0567f029e56c820b6302eab0c60928 Mon Sep 17 00:00:00 2001 From: Jeremy Teitelbaum Date: Fri, 27 Apr 2018 09:15:32 -0400 Subject: [PATCH 2/6] 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 3/6] 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 4/6] 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 5/6] 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 6/6] 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)" ] }, {