diff --git a/BDA 5.9.1-2-4.ipynb b/BDA 5.9.1-2-4.ipynb new file mode 100644 index 0000000..67737be --- /dev/null +++ b/BDA 5.9.1-2-4.ipynb @@ -0,0 +1,138 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Problem 5.9.1\n", + "\n", + "### Exchangeability with known model parameters\n", + "\n", + "For each of the following three examples, answer: \n", + "1. Are observations $y_1$ and $y_2$ exchangeable? \n", + "2. Are observations $y_1$ and $y_2$ independent? \n", + "3. Can we act as if the two observations are independent? \n", + "\n", + "Examples:\n", + "1. A box has one black ball and one white ball. We pick a ball $y_1$ at random, put it back, and pick another ball $y_2$ at random. \n", + "\n", + "Here the events are clearly independent and exchangeable. \n", + "\n", + "\n", + "2. A box has one black ball and one white ball. We pick a ball $y_1$ at random, we do not put it back, then we pick ball $y_2$.\n", + "\n", + "In this case there are four outcomes: (BB), (BW), (WB), (WW) and of these four only (WB) and (BW) have non-zero probability (1/2). Since the likelihood is symmetric, the observations are exchangeable. Clearly, though,\n", + "the events aren't independent; for example P(B|B)=0 and P(B|W)=1. And you clearly can't act as if they're independent since the second observation is determined by the first.\n", + "\n", + "3. A box has a million black balls and a million white balls. We pick a ball $y_1$ at random, we do not put it back, then we pick ball $y_2$ at random.\n", + "\n", + "These are exchangeable observations since P(BW)=(1/2)(1000000/1999999) =P(WB) and P(BB)=P(WW). They act independent since 1/2(1000000)/1999999) is just about 1/4 and so is 1/2(999999)/(1999999).\n", + "\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 134). CRC Press. Kindle Edition. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Problem 5.9.2\n", + "\n", + "We ask the same questions as in the preceeding problem but under the conditions:\n", + "\n", + "1. A box has $n$ balls colored black and white, but we don't know how many of each. We pick a ball, put it back, then pick another.\n", + "2. Same except we pick a ball, don't put it back, then pick another.\n", + "3. Suppose we know that there are a lot of balls of each color.\n", + "\n", + "In the first case, let $\\theta$ be the proportion of white balls in the urn. Then $P(BW)=(1-\\theta)\\theta$\n", + "and $P(WB)=\\theta(1-\\theta)$. So the events are exchangeable. Also $P(BB)=(1-\\theta)^2=P(B)^2$, $P(BW)=P(WB)=\n", + "P(W)P(B)$, and $P(WW)=P(W)^2$. So they are independent. \n", + "\n", + "In the second case, we have $P(WB)=(\\theta)(n(1-\\theta)/(n-1))$ and $P(BW)=(1-\\theta)(n\\theta)/(n-1)$\n", + "and these are the same, so it's exchangeable. But they are not independent since $P(WB)$ isn't $P(W)P(B)$.\n", + "Also $P(WW)=(\\theta)(n\\theta-1)/(n-1)$ and $P(BB)=(1-\\theta)(n-n\\theta-1)/(n-1)$.\n", + "\n", + "They do get close to independent if $n$ is large.\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEICAYAAABbOlNNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsvXl8FNeV6P892hFoayEJbSAWsW82wuDYsXFAhmxgTzYnmYQJ8fMkk5k3y3tJyOQ5fnbGL+Tl5ZfE45nEDomjzJLNGRvi2AaBjR1vGGH2VWLVvqMFof38/uhqaEktJOilernfz6c+3X3rVtXp7lt17rnn3HNFVTEYDAaDwUWU3QIYDAaDIbgwisFgMBgMQzCKwWAwGAxDMIrBYDAYDEMwisFgMBgMQzCKwWAwGAxDMIrBYDAYDEMwiiECEZECEVER6XTbHrZbLoNhLEQkTkSeFZHzVhteNWy/iMh3RaTZ2v6viIjb/qUisl9EuqzXpQH/EiGAUQyRTaqqTrK2b9stjMEwTt4A/hyo87DvIeA+YAmwGPgI8JfgVCrANuDfgTSgBNhmlRvcMIohBLF6S/9TRA6LSJuI/EZEEuyWy2AYC2/brqr2quoPVfUNYMBDlY3A91W1SlWrge8Df2HtWwXEAD9U1R5VfQIQ4APefKdwxCiG0OWTwDpgOs6e0V+IyFQRuXSd7TPDznFBRKpE5BkRmRzwb2CIVHzRdkdjAXDI7fMhq8y177AOzQN02G2/wSLGbgEMN80TqloDICJ/AJaq6k+A1HEc2wQsBw4C6cC/AP8BrPWTrAaDO9603bGYBLS5fW4DJll+huH7XPuTfHDdsMIohtDFfXy1C8gZ74Gq2gmUWR/rReSvgVoRSVbVdh/KaDB44qbb7jjoBJLdPicDnaqqIjJ8n2t/hw+vHxaYoaQwwjLHO6+zfXaUQ12mtYyy32DwK1603eEcw+l4drHEKnPtW+wepYRzKOsYhiEYiyGMUNWLOM3l6yIiK4BLQDnO6IwngD2qOtzMNhgCwnjbLoCIxHOtExNnOa97LN/BL4F/EJEXcXZ4/gfwz1bdPTgd1v9dRH4C/Der/BWffIkwwlgMkckM4GWcJvRRoAf4tK0SGQzj5xRwBcgFdljvp1n7ngL+ABzB2bb/aJWhqr04Q1k/j7NjtAm4zyo3uCFmoR6DwWAwuGMsBoPBYDAMwSgGg8FgMAzBKAaDwWAwDMEoBoPBYDAMISTDVSdPnqwFBQV2i2EIU/bv39+kqhmBvq5p1wZ/sn///nbgbVVdN1bdkFQMBQUFlJWVjV3RYLgJROSCHdc17drgT0SkfDxKAcxQksFgMBiGYRSDwWAwGIZgFIPBYDAYhmAUg8FgMBiGYBSDwWAwGIbgE8UgIutE5JSIVIjIZg/7460l/CpEZK+IFLjt+4ZVfkpEzEIxBr/x8ssvM2fOHGbNmsWWLVtG7O/p6eFTn/oUwMLxttPR2r6ITLfOUW61fbOusCFk8FoxiEg0zhXAPgjMBz4tIvOHVfsi0Kqqs4AfAN+1jp0PPIBzab11wL9a5zMYfMrAwABf+cpXeOmllzh+/Di/+tWvOH78+JA6P/vZz0hLSwNnVs4x2+kYbf+7wA9UtRBoxXkPGAwhgS8shtuAClU9a6Wv/TWwYVidDUCJ9f5ZYLW1WMYG4NfWwtzngArrfDdMU2cPPyg9zak6sxiTYSTvvvsus2bNYsaMGfzuQC2L7/oQ27ZtG1Jn27ZtbNy40fVxPO3UY9u3jvmAdQ5wtv37blb2tyqa+EHp6Zs93BAhdPb0852XTvjkGegLxZALVLp9rrLKPNZR1X6c66ymj/NYAETkIREpE5GyxsbGEfsHBpUf7S7n7TNNN/s9DGFMdXU1+fn5APzsT+eo6ptAdXX1qHXG2U5HK08HLlnncC8fwVjtGmD/hVZ+tLuczp5+j/sNBoBTde089dpZKlu6vD6XLxSDp+Ughy/yMFqd8RzrLFR9WlWLVLUoI2NktoLMpHhSE2M5VW8sBsNIXOuOdPcNcL75MtkpCQxd4fFaneGHcuPt12ftGmDOFOda9adN2zZch5OWpTA3O8nrc/lCMVQB+W6f84Ca0eqISAyQArSM89hxISLMnZJ09ccxGNzJy8ujsrKS8vpOBhXiui+Rk5PjsQ6Mu52OVt4EpFrncC+/KeZlO9evN8OkhutxsraDpPgYclMneH0uXyiGfUChFYURh9NJt31Yne2Aa/D248Ar1vqs24EHrKil6UAh8O7NCjJ3SjKn6zoYHDSr0hmGsnz5csrLy3l9/1F0oI+9u/7A+vXrh9RZv349JSUuV9i42qnHtm8d86p1DnC2/aEOjRsgN3UCE+OijWIwXJdTdR3MnpI0whK+GbxWDNY46l/jXHv1BPBbVT0mIo+JiOvO+xmQLiIVwD8Am61jjwG/BY7jXIP4K6o6cLOyzJmSxOXeAapar9z8FzKEJTExMTz55JM88lefpXbrl/nMpz/FggUL+Na3vsX27c5+zBe/+EWam5sBFjKOdjpa27cu+XWci9JX4PQ5/OxmZY+KEmZPSeJEbfvNnsIQ5qgqJ+vamTvF+2Ek8FF2VVV9EXhxWNm33N53A58Y5djHgcd9IYdrLPZkXTtT0xN9cUpDGPGhD32IdY/+hrYrfTz8N3cC8Nhjj13dn5CQwO9+9ztE5KiqDomOG62demr7VvlZbjLCzhNzpyTx0tE6VNUnPUJDeFHb1k17d7/PFENYzXyek+X8UYzJbRiNk3UdPrt5AsncKclc6uqjoaPHblEMQYjrmTdnSrJPzhdWimFifAxTHYmcNNEbBg80dfbQ1Nlz1bIMJa5Zw6ZtG0Zy8qpiMBaDR+aasVjDKLh6VXN91KsKJC4rx7RtgydO1rWTk5JAyoRYn5wv/BRDdjLnmy7T3XfTPmxDmOJ6qM7zQZx3oElNjCM7JcEMkxo8crK242pYsy8IO8Uwb0oSgwrl9Z12i2IIMk7WdZCRFE/6pHi7RbkpjDVs8ERP/wBnGjt9MrHNRdgphrmW1jxRZ24gw1B8Gc5nB3OzkznT2Elv/6DdohiCiDMNl+kfVJ8OkYadYpjqSGRCbDQna43JbbhG/8Agp+s7fWpuB5q5U5LoG1DONBpr2HCNk3W+HyINO8UQbU0GOmksBoMb55sv09s/GNIWg0upmbZtcOdkXQdxMVEUpE/02TnDTjEAzM92jsWOkhTNEIGcsCzIULYYZkyeSFx0lLGGDUM4UdvO7KxJxET77nEeloph7pRkWs1kIIMbJ2rbiYkSZmZMsluUmyYmOorCrEkcNw5ogxsnajt8HoIdporBOVxgbiCDixO17czKnERcTGg3+blTkq9aPwZDY4dz0qavh0hD+y4ZhauRSUYxGCxO+DjO2y7mZSfR1NlDo7GGDVx7xs33cdsOS8WQMiGWvLQJpmdlAKD1ci917d0+v3nsYH6O6fQYrnFt0qZRDONiXnayuXkMgP9uHjuYb6xhgxsnatvJTkkgbWKcT88b1orhbGOnSY1huOprCsVUGMNxpcYwisEAzrbtjw5P2CqG+dnO1Bgmt4zheG07mSGcCmM4TmvYtOtIp7tvgDONl/3S4QlbxTDPmNwGi3BxPLuYl53EmcZOevqNNRzJVDR0MjCoxmK4EfLTEpkUH2NCViOc3v5BKhrCSzHMz06hf1BNosgI57gffWdhqxiiooR52SYbZaRT0dBJ34CyICeMFIP1XUynJ7I5XtNOYly0T1NhuAhbxQDOCI7jNe0MDprUGJHKsZo24NrDNByY5kgkMS6a4zVGMUQyx2ud2YKjo3y/Bnh4K4acZC73DnCxpctuUQw2cby2nQmx/ulV2YXTGk42iiGCUVVO1LT7rcMT3oohOwUwJnckc7ymnbnZ/ulV2cmCnGSO1xprOFKpbLlCR08/C3JS/HL+sFYMhVmTiI4S07OKUFSV47XtYTHjeTjzs5Pp7OmnqvWK3aIYbOB4rTVE6qe27ZViEBGHiJSKSLn1mjZKvY1WnXIR2ehW/riIVIqIX8IrEmKjmZUx6eo4syGyqGq9Qkd3f1j5F1y4vpNp25HJ8Zp2ogTm+Gl9EW8ths3AblUtBHZbn4cgIg7gEWAFcBvwiJsC+YNV5jfm5yRzzFgMEYnrfw9Hi2F2lnN4zLTtyOR4bTszMyaREBvtl/N7qxg2ACXW+xLgPg911gKlqtqiqq1AKbAOQFXfUdVaL2W4LgtykmmwUtMaIovjNW1ECT7PVR8MuKxh4z+LTI750fEM3iuGLNeD3XrN9FAnF6h0+1xlld0QIvKQiJSJSFljY+O4j7tmcpsbKNI4VuPsVU2I80+vym4W5CSboaQIpLmzh9q2bhb6yfEM41AMIrJLRI562DaM8xqewkFuOJRCVZ9W1SJVLcrIyBj3cQusyKSj1eYGijSO1rSF1cS24czPSaa+3azNEGm4Orn+bNsxY1VQ1TWj7RORehHJVtVaEckGGjxUqwJWuX3OA/bcoJw3TUqic20GE5kUWTR19lDf3uO3cL5gwPXdjtW0sWqOJ2PdEI5c9Z0F8VDSdsAVZbQR2Oahzg7gXhFJs5zO91plAWNhTooxuSOMq72q3PC2GMAMk0YaR2vayEubQGqib9dgcMdbxbAFKBaRcqDY+oyIFInIVgBVbQG+DeyztsesMkTk/4pIFZAoIlUi8r+9lMcjC3KSOd/cRUd3nz9ObwhCXB2BBdkptLS0UFxcTGFhIcXFxbS2tno8pqSkhMLCQoCFw8Kql4nIERGpEJEnRESsco/h2iKySkTaROSgtX3LH98xZUIsUx2JxhqOMI7XtPt9iNQrxaCqzaq6WlULrdcWq7xMVR90q/dzVZ1lbc+4lX9NVfNUNcp6/d/eyDMarl6jyWEfORyraScvbQIpibFs2bKF1atXU15ezurVq9myZcuI+i0tLTz66KPs3bsX4ARDw6p/DDwEFFrbOqv8euHaf1LVpdb2mH++pXFARxod3X2ca7rs9yHSsJ757MLlvT9iHNARw7Hqa47nbdu2sXGj0wDYuHEjzz///Ij6O3bsoLi4GIfDATCAFVZt+c6SVfVtVVXgl1wLyx5PuLZfWZibwvnmLtqNNRwRuDq3QW0xhAqZyQlkJMWbnlWE0N7dx/nmLhblOjsE9fX1ZGdnA5CdnU1Dw8gYierqavLz892LXGHVudb74eVw/XDt20XkkIi8JCILRpP1ZsOwXbgeEGY4KTJwRVe62ra/GDMqKVxYlJtiQlbDmDVr1lBXVwdAV+8ANa1dfP+/JpD/ve+O63inMTCymJsLt34PmKaqnSLyIeB5nENQnq77NPA0QFFR0Q2HcS/MvRaOvXJG+o0ebggxjla3kZkUT2Zygl+vEzGKYWFOMntONXCldyBsJzxFMrt27br6fuufzvJPfzzBvv+1hsmT4snKyqK2tpbs7Gxqa2vJzBwZ2pmXl8eePXuGFOEMq66y3ruX11jvPYZrq+rV7ruqvigi/yoik1W1ySdf1o3Jk+LJTkkwnZ4I4WhN29XOgD+JiKEkcPasBtWk4I4EjlS3kZ2SwORJ8QCsX7+ekhKnK6CkpIQNG0bOzVy7di07d+50RSxFY4VVW0NEHSKy0opG+jzXwrI9hmuLyBS3yKXbcN5nzf74ruCcz2D8Z+FPV28/FQ2dRjH4EtePafwM4c/R6rYhURubN2+mtLSUwsJCSktL2bzZGTxUVlbGgw86g+ccDgcPP/wwy5cvB5iHW1g18GVgK1ABnAFesso9hmsDHweOisgh4AngAR1lrMoXLMxN5mzTZS739PvrEoYg4ERtB4PqHP3wNxEzlJSdkoBjYhxHqoxiCGc6e/o523SZjy7JuVqWnp7O7t27R9QtKipi69atVz9v2rSJTZs2ISJHh4VVlwELhx+vqs3Aag/lTwJPevtdxsui3BTUsoaXFzgCdVlDgHENFxqLwYeICAtzjckd7hyvaUfV/1EbwYTrQWE6PeHNkeo20ifGkZ3iX8czRJBiAFiUm0x5QyfdfQN2i2LwE0cCFM4XTGRZ4djGAR3eHK12Op4t95VfiTDFkMrAoHLCOKDDlqPVbWQl+z+cL9hYnJvCYaMYwpYrvQOUN3SyOC8wHZ6IUgyuH9UMJ4Uvh6susSg31W4xAs6ivBTONHYaB3SYcry2nYFBDZglHFGKITslgXTjgA5bXI7nSBpGcuFyQJtMq+HJ1RnPxmLwPSLCojzjgA5XjlW3oUrAzO1gwqUMTdsOTw5XtTF5UjxTAjREGlGKAZw3UHlDJ1d6jQM63DgSwHC+YCMzOYGs5HiOVF2yWxSDHzha3cai3OSAOJ4hQhXDwKByvNb0rMIN14znjKR4u0WxhUW5qRw2w6RhR1dvP+UNHSzKC5zvLOIUw5J8549rbqDw43BVW0QOI7lYkpfC2abLJgV3mHGspp1Bdf6/gSLiFEOWZXIbxRBetF1xLmCyOIC9qmBjsdXpOWradlhxqNI5PBgoxzNEoGIAp8l9yIzFhhWuSLNIthhcDmgznyG8OFzlHCLNTArc3JyIVAxL8lI422hM7nDCpegXR+AcBheOiXHkOyZw2HR6worDVZcC3uGJSMVgTO7w43DVJQrSE0lJjLVbFFtZnJfKoUrTrsOFti7naoSBHiKNTMVgmdyHjGIIG5yO58i1FlwsyUuh+tIVmjt77BbF4AMOVzutvyVGMfiftIlxTHUkGpM7TGho76a2rTui/QsuXMrRBFeEB67/MdCz+b1SDCLiEJFSESm3XtNGqbfRqlMuIhutskQR+aOInBSRYyKyxdOx/mJJfupVb78htHFZfkvzjcWwKDeFKIGDpm2HBQcrLzFj8sSAD5F6azFsBnaraiGw2/o8BBFxAI8AK4DbgEfcFMj/U9W5wC3AHSLyQS/lGTdL8lKoaeumob07UJc0+IlDlZeIjpIhq7ZFKhPjYyjMTDJRd2GAqnKw8tLVuVeBxFvFsAEosd6XAPd5qLMWKFXVFlVtBUqBdarapaqvAqhqL/AeQxdd9yu3THX+2MbPEPocqrrE3ClJTIiLtluUoGCpZQ37cTVRQwCoa++msaPHFkvYW8WQZS2WjvWa6aFOLlDp9rnKKruKiKQCH8VpdQSEBTkpREeJGU4KcQYHlUM29aqClSX5qbR29VHZcsVuUQxe4Ho22dG2x1zzWUR2AVM87PrmOK/hKevT1a6MiMQAvwKeUNWz15HjIeAhgKlTp47z0qOTEBvN3ClJZiw2xDnXfJn27n6WmoikqyzJdw6pHahsZWp6os3SGG6Wg5VtxEYL87KTAn7tMS0GVV2jqgs9bNuAehHJBrBeGzycogrId/ucB9S4fX4aKFfVH44hx9OqWqSqRRkZGWOJPS6W5DtnQA8OGpM7VHH1qhbnG/+Ci9lZSSTERplOT4hzqPIS87KTiY8J/BCpt0NJ24GN1vuNwDYPdXYA94pImuV0vtcqQ0T+CUgB/s5LOW6KpfmpdHQ7F3cxhCYHKy8xMS6awszA96qCldjoKBblpphh0hBmYFA5XHWJW2waIvVWMWwBikWkHCi2PiMiRSKyFUBVW4BvA/us7TFVbRGRPJzDUfOB90TkoIg86KU8N4TrRz9wsTWQlzX4kAMXL7E4L5XoqMDkqQ8VluancrSmnd7+QbtFMdwE5Q0dXO4dYOnUEFQMqtqsqqtVtdB6bbHKy1T1Qbd6P1fVWdb2jFVWpaqiqvNUdam1bfXu69wYMzMmkRQfY0zuEKW7b4ATte223TzBzNL8NHr7BzlRa5b6DEUOXHQ+k5bme5wa5ncicuazi6goYUl+qlEMIcqxmjb6B9U2czuYcYVjm7Ydmhy8eInUxFgKbAoeiGjFAE6T+2Rdh1nqMwS52qsyFsMInGma480waYhyoLKVpfmpAVvKczgRrxhumZp61dFjCC0OVF4iN3VCQPPUhwoiwi1TUzlgLIaQo6O7j/KGTm6xaRgJjGK4OqvQ3EChx4ELrcZauA63TE3jQnOXybQaYhyuakPVXks44hVD+qR4pqUn8t4FY3KHEnVt3dS0dXPrVPt6VcHOtag70+kJJVzPIjuTQka8YgC4dWoaB0xumZDCNXZ+q7EYRmVxXioxUcKBStPpCSXeu9hKYeYkUibYt+iUUQw4Hy6NHT1UtZrcMqHCexdbiYuJGldG1ZaWFoqLiyksLKS4uJjWVs8PypKSEgoLCwEWutLDA4jIMhE5IiIVIvKEWB5BEfmElTJ+UESK3M8lIt+w6p8SkbXefNebZUJcNPOyk3nvgrEYQgVV5UDlJdstYaMYcI7FgvNhYwgN3rt4iUW5KcTFjN2Et2zZwurVqykvL2f16tVs2TJy6Y+WlhYeffRR9u7dC3CCoenhf4wzT1ehta2zyo8Cfwa87n4uEZkPPAAssOr+q4jYkvr11qnOtC/9A2aiWyhwtukyl7r6uHWavZawUQzA3ClJJMZFm7HYEKG3f5Aj1W3jHkbatm0bGzc6DYCNGzfy/PPPj6izY8cOiouLcTgcAANY6eGtHGDJqvq2Oscaf4mVXl5VT6jqKQ+X3AD8WlV7VPUcUIFzLZKAc+u0NLp6BzhV32HH5Q03iOsZZCyGICAmOorFeSnsNw7okOBYTRu9/YPjvnnq6+vJzs4GIDs7m4aGkbkeq6uryc93z/V4NT18rvV+ePn1GDPVfKBw/UYmuCI0eO9iK0kJMczMmGSrHGOm3Y4Ulk1L4yevnaWrt5/EOPOzBDMuBb5s2jXFsGbNGurq6kbUffzxx8d1zlECD5Qx0saPwriP8XU6+eHkpU0gIyme/Rda+dztBT4/v8G37D/fyq1T04iyOfeXeQJaFE1zMDB4hkOVbdw+M91ucQzXYf+FVvIdE8hMvjaxbdeuXaPWz8rKora2luzsbGpra8nMHLmeVF5eHnv27BlSBOzB2dvPG1bunjbeE2Olmr+Kqj6NM/U8RUVFPg+LExGKpqWx3/jPgp62K32cbujgI4uz7RbFDCW5cOWWMQ7o4EZVKbvQyrIbGINdv349JSXOFWhLSkrYsGHDiDpr165l586droilaKz08NbKhB0istKKRvo8ntPLu7MdeEBE4kVkOk6H9bvjFtjHLJuWRmXLFbO+eZBz4GIrqkMtYbswisEiNTGOWZmTjJ8hyKlqvUJjR88N3TybN2+mtLSUwsJCSktL2bx5MwBlZWU8+KAzCbDD4eDhhx9m+fLlAPOw0sNbp/gysBWnE/kM8BKAiNwvIlXA7cAfRWQHgKoeA34LHAdeBr6iqrYl47p1mom6CwXeu9BKtJXY027MUJIbRdPSeOloHYODavsYn8Ez1/wLjnEfk56ezu7dI5cTLyoqYuvWa5neN23axKZNmxCRo6708OBMIw8sHH68qj4HPOfpmqr6ODA+B4efWZjjDOstO9/KuoX2D1MYPLP/YivzspOYGG//Y9lYDG7cOi2Ntit9nGnstFsUwyjsv9DKxLho5kwxK7aNl7iYKJbkpVBmrOGgpX9gkAMX7Z/Y5sIoBjeWFzh7oeYGCl72nW/h1mlpZsW2G6SowMGxmjaTXj5IOVHbQVfvwNVnkN0YxeBGQXoikyfFse98y9iVDQGn7Uofp+o7gubmCSWWF6TRN6AcMunlgxLXM6eowFgMQYeIsGxaGmXnjcUQjLxnRW0UBUHURqjhGqIoM52eoKTsQgu5qRPITplgtyiAUQwjWF7g4GJLF/UmtC/oKDvfQnSUmDUYboLUxDhmZ01in+n0BB2qStn5VpYHibUARjGMoMjlZzA3UNCx73wrC3OSzcz0m6SowMF7F1oZGDTp5YOJypYrNHT0XH32BANGMQxjQU4yE2KjjZ8hyOjpH+BQ5aWgunlCjeUFaXT09HOyrt1uUQxuvGs9a4LJd2YUwzBio6O4ZWoq754ziiGYOFLVRk//YFDdPKGG67fbZ9p2ULHvXAspE2IpzLQ3cZ47XikGEXGISKmIlFuvHgfJRGSjVad82AIoL4vIIWuxk5/YlbN+OLdNd3Cirp22K312i2Kw2HvO1asKnnHYUCMvLZHc1AlXe6iG4ODd8y0sL3AE1aRaby2GzcBuVS0EdlufhyAiDuARYAXOnPTuC6B8UlWX4JxVmgF8wkt5fMJt0x2omlTFwcS+8y0UZk4ifVK83aKENLdNd/DuuVazjG2Q0NDRzbmmy6yYHlyWsLeKYQNQYr0vwVrAZBhrgVJVbVHVVqwFUABU1TXYGQPEMXY644BwS34asdFytZdqsJeBQWX/+VaWB9nNE4osL3DQ1NnDuabLdotiAPadc3Y+g61te6sYsqzsk1ivI/MZj7FoiZV4rAHoAJ4d7UIi8pCIlIlIWWNjo5diX58JcdEsyk0xDugg4URtOx09/UHXqwpFbrN+Q9O2g4N3zzWTGBfNgpxku0UZwpiKQUR2ichRD9vI3MWjnMJD2VXLQFXXAtlAPPCB0U6iqk+rapGqFmVkZIzz0jfPbdPTOVx1yaQQCAKu+ReMYvCWmRkTSZ8Yx96zRjEEA3vPtXDr1DRio4MrDmhMaVR1jaou9LBtA+qtNXGxXkeumTiORUtUtRtnDvvxKhu/s2KGg74BNamKg4C9Z5uZ6kgkJzU4ZoWGMiLCihkOM0waBLRe7uVkXQcrZwRfh8dbNbUdcEUZbcTzAiY7gHtFJM1yOt8L7BCRSW5KJQb4EHDSS3l8RtG0NKLE+VAy2MfgoPLu+RYzjORDVkxPp/rSFSpbuuwWJaJxRYetmBF8K0Z6qxi2AMUiUg4UW58RkSIR2QpgLXbybWCftbkWQJkIbBeRw8AhnNbGT7yUx2ckJcSyMDeFd0zPylZO1XdwqasvKG+eUGWF1UN9x3R6bGXv2RbiY6JYnJditygj8Cq3gKo2A6s9lJcBD7p9/jnw82F16oHl3lzf36yckc4v3jxPd98ACbFBMcUi4nBZbMZi8B2zM5NIS4xl77kWPlGUP/YBBr/wztlmlk1LIz4m+J4tweXxCDJWTHfQay2gYbCHd846s07mOxLtFiVsiIoSbpvuMBaDjbRd6eNEXTsrpgenJWwUw3VYPt1BlMDb5gayhcFB5Z1zzaw0w0g+Z+WMdKpajZ/BLvaebUaVoHQ8g1EM1yXZ5Wc4YxSDHZysc/oXbp9EleZ6AAAgAElEQVRpFIOvcf2mptNjD2+fbSY+JipoU8gbxTAGt89I50Blq5nPYAOuh5ZRDL5ndmYSjolxptNjE2+faaaoIDj9C2AUw5jcPjOdvgGl7IKJTgo0b59pZlq6M/GbwbdERQm3z0jn7bPNJm9SgGnu7OFkXQfvmznZblFGxSiGMVhe4CAmSnjb9KwCysCgsvdcM7cb/4LfWDkzndq2bi40Gz9DIHFNLgxm35lRDGMwMT6GxXkpvGUUQ0A5VtNGR3e/GUbyIy6la9p2YHnrTBOJcdFBOX/BhVEM4+COWZM5XHWJ9m6zPkOgeLPC+Bf8zcyMiWQlx/PmmSa7RYko3qpoZsV0R9DlR3IneCULIt43czKDikk8FkDeOtPE7KxJZCYl2C1K2CIi3DFzMm+faWbQrAMdEGrbrnC26TJ3zApe/wIYxTAubp2WSkJsFG9WmJ5VIOjpH2Df+Zagds6FC++bNZmWy72cMOtABwSXJRzsbdsohnEQHxPN8gIHbxmTOyC8d+ES3X2DQd+rCgfumGX5GSqMnyEQvFXRhGNiHHOnJNktynUximGc3DFrMqfrO2lo77ZblLDnrTNNREfJ1WRvBv+RnTKBGRkTecNYw35HVXnzTBO3z0wPqvWdPWEUwzi50+q9mhvI/7xe3sSSvBSSE2LtFiUiuHPWZN4910JPv5nE6U8qGjqpb+/h/SFgCRvFME7mZyfjmBjHG+VGMfiTtq4+jlRd4s5C/6/SZ3Dy/sIMrvQNsP+CWZTKn7xuPTvuLDSKIWyIihLumDWZP1U0mZmifuStM00MKtzlw5unpaWF4uJiCgsLKS4uprXV8wOwpKSEwsJCgIUi4lqAChFZJiJHRKRCRJ4QEbHKPyEix0RkUESK3OoXiMgVETlobUGzzognVs5wEB0lptPjZ94ob2TG5InkpQV/pmCjGG6A98+aTGNHD6fqO+wWJWx5vbyJSfExLMn3XXKxLVu2sHr1asrLy1m9ejVbtmwZUaelpYVHH32UvXv3ApwAHrFWHAT4MfAQUGht66zyo8CfAa97uOwZVV1qbV/y2ZfxA0kJsdw6NdUMk/qR3v5B9p5rCQlrAYxiuCFcf+qfTpsbyB+oKm9UNLJyRrpPJ/9s27aNjRudBsDGjRt5/vnnR9TZsWMHxcXFOBwOgAGgFFhnLT+brKpvq9NU/CVwnyXvCVU95TNBbeTOWRkcqW6j9XKv3aKEJe9dbKWrd+CqrzLYMYrhBshJncCszEm8Xt5otyhhyfnmLipbrnDXbN/ePPX19WRnZwOQnZ1NQ0PDiDrV1dXk5w9ZzawKyLW2Kg/lYzFdRA6IyGsi8v7RKonIQyJSJiJljY32tav3z56Mqgmu8Bevn24kJkpCZia/V0t7RiJ3FWbw73svcKV3gAlxwZkyN1R57ZTzgX337Bt3PK9Zs4a6uroR5Y8//vi4jh/Fb6SAp7jCsZxMtcBUVW0WkWXA8yKyQFVHzCJT1aeBpwGKiopsc14tyUslZUIsr51u5KNLcuwSI2x57XQjt05LIylEIu2MYrhB7p6Twc/fPMfec82smpNptzhhxevlTRSkJzItfeINH7tr165R92VlZVFbW0t2dja1tbVkZo783/Ly8tizZ8+QImAPTgshb1h5zfVkUdUeoMd6v19EzgCzgbLxfBc7iI4S7iyczOunG1FVLP+6wQc0dvRwrKadr66dY7co48YMJd0gK6Y7iI+J4rXTZjjJl/T0D/D2meabshbGYv369ZSUlADOyKMNGzaMqLN27Vp27tzpiliKBu4FdqhqLdAhIiutaKTPA9uudz0RyRCRaOv9DJwO67M+/Ep+4e7ZGTR0ONcKMPiOP1lDz/5o2/7CKIYbJCE2mhUz0o1i8DFl51u50jfAXX64eTZv3kxpaSmFhYWUlpayefNm5zXLynjwwQcBcDgcPPzwwyxfvhxgHvCYqrqyJn4Z2ApUAGeAlwBE5H4RqQJuB/4oIjus+ncBh0XkEPAs8CW3cwUtd1lzR0zb9i2vnW4kfWIc87OT7RZl3JihpJtg1ewMHnvhOJUtXeQ7gj8mORTYc6qBuOgovzjn0tPT2b1794jyoqIitm7devXzpk2b2LRpEyJyVFWfcZWrahmwcPjxqvoc8JyH8t8Dv/eR+AFjSkoCc6ck8dqpRr5090y7xQkLBgaV1083smpOZtCnwXDHK4tBRBwiUioi5dZr2ij1Nlp1yt0nDrnt3y4iR72RJZDcM9c5Rr3n1MjoFsPN8eqpRlbMcJAYZ/oqdnLP3Ez2nW+hw6w94hMOVV2itavv6jMjVPB2KGkzsFtVC4Hd1uchiIgDeARYAdzG0IlDiMifAZ1eyhFQpk+eSEF6Iq+eMia3L6hs6aKioZN7jDPfdu6Zk0n/oJoU8z5iz8kGosS3M/kDgbeKYQNQYr0vwZr4M4y1QKmqtqhqK9bEIQARmQT8A/BPXsoRcFbNyeStM01095nEY97isrxCrVcVjtw6NZWkhBhePWk6Pb7g1VON3Do1jdTEOLtFuSG8VQxZVtQG1qunOzsXqHT77D5B6NvA94ExVyMPlolALlbNyaC7b5C3z5o89t7y6qlGpqUnMn3yjYepGnxLTHQUdxVm8OqpBpMTzEsaOro5Ut3GqjmhE43kYkzFICK7ROSoh21kzN8op/BQpiKyFJhlOfDGRFWfVtUiVS3KyLD/h145I50JsdG8etL4GbzhSu8Ab1Y08QFjLQQNH5ibSUNHD0erzapu3rDHsrpWz8uyWZIbZ0zFoKprVHWhh20bUG/lksF69fSUrALccw24JgjdDiwTkfPAG8BsEdnj3dcJHAmx0dxZOJndJ0zPyhverGiip3+Q1XND7+YJV1bNyUAEdp+st1uUkGbXiXpyrEivUMPboaTtgCvKaCOeJ/7sAO4VkTTL6eyaOPRjVc1R1QLgTuC0qq7yUp6AsnpuJtWXrphsq16w+2QDk+JjuG26Wa0tWEifFM8t+ansPmGs4Zulu2+ANyqa+MC8zJCcRe6tYtgCFItIOVBsfUZEikRkK4A1sefbwD5reywUJvuMB9fwh7mBbg5V5ZWT9dw1ezJxMWauZTCxel4WR6rbqDdL2d4Ue8+10NU7EJLDSOClYlDVZlVdraqF1muLVV6mqg+61fu5qs6ytmc8nOe8qo6YQBTsZCYnsDgvhV0njMl9Mxytbqe+vYcPmGGkoGP1PNPp8YZdx+uZEBvN7TNCI5vqcEw3zUuK52VxsPISDaZndcOUHq8jSjCO5yBkTlYS+Y4JlB4fmbHWcH1UlV0nnJZwQmxoZmA2isFLihdkoQq7TM/qhtl5vJ6iAgeOiaEV4x0JiAjF86bw5plmLvf02y1OSHG0up3atm6K50+xW5SbxigGLzE9q5ujsqWLk3Ud3DvfDCMFK/cuyKK3f5DXTVK9GyIcLGGjGLzEvWfVaXpW42bncadfptgohqClaFoaqYmxV/8rw/gIB0vYKAYfsNbqWZmkeuNnx9E65mQl3dSiPIbAEBMdxZp5Wew+UU9v/6Dd4oQE55suh4UlbBSDDygqcJA+MY6Xj5rhpPHQ2NHDvgstrFsYumOwkcK6BVNo7+7nHZP6ZVzsOOZ8BoR62zaKwQdERwn3Lsji1ZMNJqneOCg9Xo9q6N88kcCdhZNJjIvm5WOm0zMeXjpax6LcFPLSQnudFqMYfMTaBVO4bOX9MVyfl4/VMS09MSRTBUQaCbHR3DM3k53H6hkYNKlfrkdt2xUOVl4Kiw6PUQw+4n0zJ5OUEMOLR0zP6npc6urlrYom1i2YEpKpAiKRdQum0NTZQ9n5sEhY4DdcQ8lrFxjFYLCIi4mieH4WpcfrjKPuOuw8Xk//oPKhRdl2i2IYJx+Ym0l8TBQvHqm1W5Sg5sUjtczJSmJW5iS7RfEaoxh8yIcXZdPe3W+Gk67Di0dqyUubwOK8FLtFMYyTifExrJqTwUtH68xw0ijUtXVTdqE1bDo8RjH4kDsLncNJLxw2PStPXOrq5Y3yJj68KNsMI4UYH16cQ0OHGU4ajZeO1qIKH14c+sNIYBSDT4mPiaZ4fhY7j9fR02+ik4az85gZRgpVVlvDSX80w0ke+eNh1zBSeARUGMXgYz66JIeO7n5eP22Gk4bzh8M1TEtPNMNIIcjE+BhWz8vkxSO19A8YH5o71ZeuUHahlY8uCZ8Oj1EMPubOWZNJS4xl+6Eau0UJKho7enizoomPLs4xw0ghyvolOTR19pp1zofxgnWvf2Rxjs2S+A6jGHxMbHQUH1yUza7j9XT1mtxJLl48UsugOi0qQ2iyak4mk+Jj2H7QdHrc2X6ohiV5KRRMDp/0LkYx+IH1S3K40jdAqUk+dpXth2qYk5XEHDOpLWRJiI3m3gVZvHyszszwtzjT2Mmxmvaw6/AYxeAHbitwkJOSwPMHqu0WJSi42NzF/gut3HdLrt2iGLzk/lty6eju59WTJmEkwPMHqomS8LOEjWLwA1FRwoZbcnm9vInGjh67xbGd5ywFuWFpeN08kcj7Zk4mMyme/zKdHlSV5w5Uc8esyWQlJ9gtjk8xisFP3H9LLgODyguHI3s8VlV5/mA1K2c4yEmdYLc4Bi+JjhI2LM1hz6kGWi/32i2Orey/0EpV6xXuD0NL2CgGPzE7K4n52cn8/r0qu0WxlQOVlzjXdDksb55I5b5bcukbUP4Q4Z2e379XzYTY6LDIjTQcoxj8yMeX5XG0up2Tde12i2Ibz+6vIiE2ytZJbS0tLRQXF1NYWEhxcTGtra0e65WUlFBYWAiwUEQ2uspFZJmIHBGRChF5Qqx4WxH5noicFJHDIvKciKS6HfMNq/4pEVnr568YUBbkpDAvO5ln90dup+dK7wAvHKrhg4umMDE+xm5xfI5XikFEHCJSKiLl1mvaKPU2WnXKh91we6wb56C1he4iqR6475ZcYqOFZ8si8wbq7hvgD4dq+NDCbJISYm2TY8uWLaxevZry8nJWr17Nli1bRtRpaWnh0UcfZe/evQAngEfc2vOPgYeAQmtbZ5WXAgtVdTFwGvgGgIjMBx4AFlh1/1VEov32BW3g48vyOFzVxqm6DrtFsYWdx+vo6OnnE8vy7RbFL3hrMWwGdqtqIbDb+jwEEXEAjwArgNsYesMBfFZVl1pbWIU6OCbG8YG5mTx/sJq+CJwtuuNYHR3d/Xx8WZ6tcmzbto2NG539kY0bN/L888+PqLNjxw6Ki4txOBwAAzgf+utEJBtIVtW3VVWBXwL3AajqTlV1TVZ5B3B90Q3Ar1W1R1XPARU4237YcN/SHGKihGf3V9otii38rqyKvLQJrJjusFsUv+CtYtgAlFjvS7BumGGsBUpVtUVVW7FuOC+vGzJ8siifps5edp8IK503Ln5bVkle2gRWzki3VY76+nqys51DWdnZ2TQ0jPwvqquryc8f0vurAnKtrcpD+XA2AS9Z73MB9yfmaMcgIg+JSJmIlDU2No7r+wQD6ZPiWT0vk/96rzri0sxXtnTxRkUTH1+WR1RUeM7i91YxZKlqLYD16mkoaKyb5BlrGOlh19itJ0L1Brp7dgZZyfH8Zt9Fu0UJKBeaL/NmRTOfKsoPyM2zZs0aFi5cOGLbtm3buI53GgMjiwFPwg+pLCLfBPqB/3AVjXWM23WfVtUiVS3KyMgYl6zBwgPLp9J8uZddJyJrIudvyyqJEmenL1wZ02siIrsAT273b47zGte7ST6rqtUikgT8HvgcTlN95AGqTwNPAxQVFYVMUviY6Cg+WZTPv7xaQc2lKxETsvmbfc6b5xMBunl27do16r6srCxqa2vJzs6mtraWzMyR/Ze8vDz27NkzpAjYg7Mjkzes/Go4juUz+wiwWq9plyogf7RjwoW7ZmeQk5LAr/dVRkzG3P6BQX5bVsndszPC+l4e02JQ1TWqutDDtg2ot8ZgsV49jZeMepOoarX12gH8J2E2Duvik0X5KM6HZSTQNzDI7/ZXcc+cTKak2D/xZ/369ZSUOEc8S0pK2LBhw4g6a9euZefOna6IpWjgXmCHZQl3iMhKy6L9PLANQETWAV8H1qtql9vptgMPiEi8iEzH6bB+129f0Caio4RPFOXzp/JGKlu6xj4gDHj1VCP17T18avlUu0XxK94OJW0HXFFGG7FumGHsAO4VkTTL6XwvsENEYkRkMoCIxOLsdR31Up6gJN+RyF2FGfx638WIcELvPFZPY0cPn10ZHDfP5s2bKS0tpbCwkNLSUjZvdsZIlJWV8eCDDwLgcDh4+OGHWb58OcA84DFVda1K82VgK04n8hmu+RKeBJKAUms49CcAqnoM+C1wHHgZ+IqqhmVyoU8tz0eA/3w3MoZK//2dC2Qlx7NmXlgFUI5EVW96A9JxRiOVW68Oq7wI2OpWbxPOm6oC+IJVNhHYDxwGjgE/AqLHc91ly5ZpqFF6rE6nff0FffFwjd2i+J1PPfWW3rFlt/YPDNotyk0BlKkX98XNbqHYrlVV/1vJPr3lsZ3a3ddvtyh+5XxTpxZsfkF/UHrKblFuihtp115ZDKrarKqrVbXQem2xystU9UG3ej9X1VnW9oxVdllVl6nqYlVdoKp/q2HaqwK4Z24muakT+Ld3Ltgtil+paOjgnbMtfHbFNKLDNGLDMJTP3T6Nlsu9vHSkzm5R/Mp/7r1IlAifvi04LGF/YmY+B4joKOGzK6fy1plmTteH76SgX7x1nriYKD5ZZO/cBUPguGPmZGZMnsgv3jpvtyh+o6u3n1/vq2TtgqywS5jnCaMYAsgDy6cSHxPFM2+et1sUv9DW1cfv91ezYUkO6ZPi7RbHECCiooSN7yvgYOUlDlz0nG4k1HnuQDVtV/r4wh3T7RYlIBjFEEAcE+O4b2kuzx2oCsvMlL8pu8iVvgH+4o4Cu0UxBJiPLcsjKT4mLDs9qsov3jzPgpxkiqZ5zPoTdhjFEGC+cGcB3X2DYRfF0TcwyC/ePM9t0x0syEmxWxxDgJkUH8Mnl+fz4pFaai5dsVscn/La6UbKGzr5wh3TI2a9cqMYAszcKcncNTuDZ948T09/+PjaXzxSS01bN3951wy7RTHYxBfuKECBZ948Z7coPuWnfzpLVnI868NslbbrYRSDDTz0/hk0dfaEzdKfqspTr51lZsZE7pkT5vHdhlHJS0vkI4uz+dW7lbR399ktjk84Wt3GmxXNfOGO6cTFRM7jMnK+aRBxx6x0FuQk89RrZxkYDJnsHqPyenkTx2vbeeiuGWGbVMwwPh66awadPf3829vhEZb949fOMCk+hs+sCP8QVXeMYrABEeGvVs3ibNNlXjpaa7c4XvMvr1SQnZLA/beYENVIZ0FOCqvmZPDzN85xpTe0h0rPNHby4pFaPnf7NJJtXE/EDoxisIl1C6cwM2MiT75SMVpmz5Dg3XMtvHu+hb+8a0ZEmdqG0fnre2bRfLk35AMsfrznDPExUXzxzsgIUXXH3Mk2ER0lfOWeWZys62DHsdBNW/zDXaeZPCmOByJgNqhhfBQVOFg5w8FTr52huy80rYYLzZd57kA1n75tKpMjcE6OUQw2sn5JDjMmT+SHu04zGIK+hrfPNPPWmWa+vGoWCbFhtXKlwUv+bs1sGjp6+PcQTQHzo93lxEYLX141025RbMEoBhuJiY7ib9cUcrKugz8eCS1fg6ryg9LTZCXH89kIc8wZxmbljHTumJXOT147w+We/rEPCCLONHby/IFqPrdyGplJ4Z/+whNGMdjMRxbnMCcrie/vPBVSKblfPdXAu+db+Ot7jLVg8Mz/uHcOTZ29/OyN0JrX8L2XTzEhNpq/vDsyrQUwisF2oqOEr39wDuebu/h1iDjrBgaV7750ioL0RONbMIzKrVPTWLdgCk+9doamzh67xRkX711s5eVjdTx018yI9C24MIohCLhnTia3TXfww13lITEx6Nn9lZyq7+Cra+cSG22akGF0vrpuDt39g/xw12m7RRkTVeXxP55g8qR4Hnx/5EUiuWPu6iBARHj4w/Np6erlyVcq7BbnunR09/G9HacompbGhxZ5WgrcYLjGzIxJ/PmKqfzn3oucrGu3W5zrsv1QDfsvtPLVtbOZGB9jtzi2YhRDkLAoL4VPLMvjmTfPcbax025xRuWfX6mg+XIvj3x0QcQkFDN4x98XzyZ5QiyPbj8etHN2unr72fLSSRbmJvPxZfljHxDmGMUQRPzPtXNIiI3m4W1Hg/IGOlnXzs/eOMcnl+WzKM9kUDWMj9TEOP7HvXN4+2wz2w/V2C2OR360q5zatm4e+egCs/IgRjEEFZlJCXxt7RzerGhm28HguoEGB5VvPneUlAmxbP7gXLvFMYQYn7ltKkvyU/n2C8dp6wouP9rJuna2vnGOTxXls7zAYbc4QYFRDEHGZ1ZMY2l+Ko+9cDyoIjl++fZ59l9o5R8/NI+0iXF2i2MIMaKjhP9z/0Jau/p47IXjdotzlf6BQb727GHT4RmGUQxBRnSU8L2PL6azu5+Hnw+OIaULzZf57sunWDUng4/dmmu3OIYQZUFOCl++eya/f6+KV04GRxqYp14/y+GqNr69YaHp8LhhFEMQUpiVxN8Xz+alo3X8/j1712zoGxjk735zkJho4Tt/tsg4nA1e8TerZzEnK4mvPXvEdov4SFUbP9x1mg8vyubDi7NtlSXYMIohSHnorhmsmO7gW9uOcq7psm1y/GhXOQcuXuL/3L+I7JQJtslhCA/iY6L50aeX0t7dx1d/d8i2HGGXe/r5778+wORJ8Tx+/0JbZAhmvFIMIuIQkVIRKbdePa6ULSIbrTrlIrLRrTxORJ4WkdMiclJEPuaNPOFEdJTwg08tJS4mii//+366egOfb+bVkw38y54KPrEsj49G0LKGBv8yd0oy/+vD83j1VCM/fu1MwK+vqnzt2cNcaL7MDz61lNREM4Q0HG8ths3AblUtBHZbn4cgIg7gEWAFcBvwiJsC+SbQoKqzgfnAa17KE1bkpE7gRw/cwqn6Dr7++yMB9Teca7rM3/76APOmJPPYBtOjMviWz62cxvolOfy/nad49VRDQK/99Otn+eORWr62bi4rZ6QH9NqhgreKYQNQYr0vAe7zUGctUKqqLaraCpQC66x9m4DvAKjqoKo2eSlP2HH37Ay+unYOfzhUww9KA5NWoOVyL1945l2io4SnPreMCXEmSZ7Bt4gIWz62iLlTkvmb/zzA8ZrAzIp+6UgtW14+yYcXZfOXd80IyDVDEW8VQ5aq1gJYr55Wgs8FKt0+VwG5IpJqff62iLwnIr8TkazRLiQiD4lImYiUNTY2eil2aPHlu2fyyaI8nnilgn97+7xfr9XZ08+mX+yjpq2brRuLyHck+vV6hsglMS6Gn/9FEZPiY/jCL97lQrN/fWlvn2nm735zkKX5qXz/k0tMIMV1GFMxiMguETnqYdswzmt4+vUViAHygDdV9VbgbeD/jXYSVX1aVYtUtSgjI2Oclw4PRITH71/EmnmZPLztGL/Z558srC6lcKS6jX/+9C0sm2Ym+xj8S3bKBEo23UZP/yCf+eleKlu6/HKdfedb+GLJPqY6Etn6+SKTKn4MxlQMqrpGVRd62LYB9SKSDWC9ehosrALck4/kATVAM9AFPGeV/w641YvvEtbERkfx5Gdu5a7ZGXz990f46etnfepzaLncy2d/+g77L7Tyg08tZe0CkyDPEBjmTEni3zatoKO7j4//5C1O1XX49Pyvnmzgcz/by5TkBP7jwRWkR3A67fHi7VDSdsAVZbQR2Oahzg7gXhFJs5zO9wI71PlU+wOwyqq3GgieKZFBSEJsND/9/DI+vCibx188wTf+6wg9/d6vqXu8pp31T77ByboOnvrzZawPswiklpYWiouLKSwspLi4mNbWVo/1SkpKKCwsBFg4LHpumYgcEZEKEXlCrDEIEfmeFU13WESecw2PikiBiFwRkYPW9pMAfM2QZlFeCr/90u2owsd+/Balx72fAKeq/PT1s3yxZB+zMifx2y/dTmZyZK7IdsOo6k1vQDrOaKRy69VhlRcBW93qbQIqrO0LbuXTgNeBw9bxU8dz3WXLlmkkMzAwqN97+aRO+/oLuvYHr+mJ2rabOk//wKBu/dNZLfzmi7ri8V168GKrjyUNDr761a/qd77zHVVV/c53vqNf+9rXRtRpbm7W6dOna3NzswIHgLNAmjrb6bvA7TiHRV8CPmiV3wvEWO+/C3zXel8AHNUbvJ8ivV2rqtZc6tKP/vOfdNrXX9BvPX9Eu3r6b+o89e1XdNMz7+q0r7+gX/q3Mu3s7vOxpKEHUKbjbIuiQZBy4UYpKirSsrIyu8Wwnd0n6vn67w/T2tXH51ZO46/umTmuNWpVlbfONPOdl05wtLqd1XMz2fKxxWQkhaeJPWfOHPbs2UN2dja1tbWsWrWKU6dODanzq1/9ij179vDUU08hIvuB/cAea3tVVecCiMingVWq+pfux4vI/cDHVfWzIlIAvKCqNxTna9q1k+6+Ab778kmeefM8OSkJ/M+1c1i/JIeYcSwKdbmnn39/5wJPvlJB78AgX183ly/cUWAczYCI7FfVovHUjezVKEKc1fOy2Pn3d/P/lZ7il2+f5z/3XuSDi6bw4UXZrJieTkpi7NW6qsqZxk72nGrkuQPVHKtpJzslgR89sJT1S3LC+sapr68nO9uZ8iA7O5uGhpGusOrqavLzh+Thr8IZUZdrvR9ePpxNwG/cPk8XkQNAO/C/VPVPnmQTkYeAhwCmTjXLpIJzyPSRjy7gQ4uy+d/bj/EPvz3E93ee5mO35rJ6Xhbzc5KHrBx4pXeAA5Wt7Dhax/MHa2i70sc9czJ4+CPzmZExycZvEroYxRDiOCbG8U/3LeLBO2fwszfOsf1QzdWU3ZMnxZGaGMfAoFLX1s2VPqc/YkFOMt/5s0Xcf0tu2ERnrFmzhrq6uhHljz/++LiOH8VyVkaPqruKiCOWENwAAAatSURBVHwT6Af+wyqqxTks2iwiy4DnRWSBqo4I1lfVp4GnwWkxjEvYCGF5gYM//PWd7D7ZwC/eOseTr1bwxCsVxEYLU1ISSIiJprOnn/r2bgYV4mKiKJ6XxaY7p7NsmsckDIZxYhRDmFAweSLfvm8hD39kPvsvtHKgspXKli7arvQRJcIH5mYyJyuJ26Y7KJg80W5xfc6uXbtG3ZeVlUVtbe3VoaTMzJHTbfLy8tizZ8+QIpzDSFXWe/fyq4tlWE7qjwCrrXFcVLUH6LHe7xeRM8BswIwT3SBRUULx/CyK52fR2NHD3nPNHK1up67tCr0Dg0yIjSE3bQJL81NYXuAgKSF27JMaxsQohjAjLiaK22emc/tMM9Xfxfr16ykpKWHz5s2UlJSwYcPIKThr167lH//xH10RS9E4HcvfUNUWEekQkZXAXuDzwD8DiMg64OvA3ap6NQBfRDKAFlUdEJEZQCFOZ7bBCzKS4vnI4hw+sji8ouaCEZNd1RD2bN68mdLSUgoLCyktLWXzZmdKr7KyMh588EEAHA4HDz/8MMuXLweYBzymqi3WKb4MbMUZVXcGZ2QSwJNAElA6LCz1LuCwiBwCngW+5HYugyHoMVFJBsMwbiR6w5eYdm3wJzfSro3FYDAYDIYhGMVgMBgMhiEYxWAwGAyGIRjFYDAYDIYhGMVgMBgMhiEYxWAwGAyGIYRkuKqINAIXRtk9GQiGJUKDRQ4IHlmCRQ64vizTVDXgq0GFSLuG4JElWOSA0JClEHhbVdd52DeEkFQM10NEyuyIQQ9WOSB4ZAkWOSC4ZBkPwSRvsMgSLHJA+MlihpIMBoPBMASjGAwGg8EwhHBUDE/bLYBFsMgBwSNLsMgBwSXLeAgmeYNFlmCRA8JMlrDzMRgMBoPBO8LRYjAYDAaDFxjFYDAYDIYhhIxiEJF1InJKRCpEZLOH/fEi8htr/15rQXbXvm9Y5adEZG0AZPkHETkuIodFZLeITHPbN2Dl7j8oItv9LMdfiEij2/UedNu3UUTKrW2jN3KMU5YfuMlxWkQuue3z5W/ycxFpEJGjo+wXEXnCkvOwiNzqts+nv8kNyBwUbTtY2vU4ZQlI2w6Wdm2dL3BtW1WDfsO5otYZYAYQBxwC5g+r81fAT6z3DwC/sd7Pt+rHA9Ot80T7WZZ7gETr/ZddslifOwP4m/wF8KSHYx04VxRzAGnW+zR/yjKs/t8AP/f1b2Kd6y7gVuDoKPs/hHOhHQFWAnv98ZuEWtsOlnYdTG07mNp1oNt2qFgMtwEVqnpWVXuBXwPD12fcAJRY758FVouIWOW/VtUeVT2HcxWu2/wpi6q+qteWenyHoWsG+4rx/CajsRYoVdUWVW0FSoExZ0P6UJZPA7/y4nqjoqqvA9dbLW0D8Et18g6QKiLZ+P43GS/B0raDpV2PS5br4Mv/MWjaNQS2bYeKYsgFKt0+V1llHuuoaj/QBqSP81hfy+LOF7m2FCRAgoiUicg7InJfAOT4mGVWPisi+Td4rK9lwRp+mA684lbsq99kPIwmq69/E2/l8VjHj207WNr1jcji77YdSu0afNi2Y3wumn8QD2XD42xHqzOeY30ti7OiyJ8DRcDdbsVTVbVGnIvEvyIiR1T1jJ/k+APwK1XtEZEv4ex1fmCcx/paFhcPAM+q6oBbma9+k/EQqHYyXoKlbQdLux6vLIFo26HUrsGH7SRULIYqIN/tcx5QM1odEYkBUnCaXeM51teyICJrgG8C61W1x1WuqjXW61lgD3CLv+RQ1Wa3a/8UWHYj38GXsrjxAMPMbR/+JuNhNFl9/Zt4K4/HOn5s28HSrsclS4Dadii1a/Bl2/alc8RfG07L5ixOU83lBFowrM5XGOqg+631fgFDHXRn8c75PB5ZbsHptCocVp4GxFvvJwPlXMeZ5QM5st3e3w+8o9ecUecsedKs9w5//iZWvTnAeayJlb7+TdzOWcDoDroPM9RB964/fpNQa9vB0q6DqW0HW7sOZNv2a6P35YbT437aapjftMoew9lzAUgAfofTAfcuMMPt2G9ax50CPhgAWXYB9cBBa9tulb8POGI1sCPAF/0sx3eAY9b1Xv3/27ljE4ShKArDv9s4gAupTcawdAFncCBbe4ewsTCFz1olwvd1IVwIlwPnwYNU65fZ3byra7X99k7m50N1fJv79E7O1a269zwp7aupmub3q+o0f+el2nxrJ/+W7aXkeknZXkquf51tv8QAYPAvdwwA/IhiAGCgGAAYKAYABooBgIFiAGCgGAAYPACVdz2wO7B8QAAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "# These graphs show the difference Pnr(WW)-Pr(WW) where Pnr means no replacement and Pr means with replacement.\n", + "x=np.linspace(0,1,100)\n", + "n=5\n", + "fig,ax=plt.subplots(1,2)\n", + "ax[0].plot(x,x*(n*x-1)/(n-1)-x*x)\n", + "ax[0].set_title('n=5')\n", + "n=100\n", + "ax[1].plot(x,x*(n*x-1)/(n-1)-x*x)\n", + "#ax[1].plot(x,x*x)\n", + "ax[1].set_title('n=100')\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Problem 5.9.4\n", + "\n", + "Exchangeable prior distributions: suppose it is known a priori that the $2J$ parameters $\\theta_1,\\ldots,\\theta_{2J}$ are clustered into two groups, with exactly half being drawn from a $N(1, 1)$ distribution, and the other half being drawn from a $N(−1 , 1)$ distribution, but we have not observed which parameters come from which distribution. \n", + "\n", + "1. Are $\\theta_1,\\ldots,\\theta_{2J}$ exchangeable under this prior distribution? \n", + "2. Show that this distribution cannot be written as a mixture of independent and identically distributed components.\n", + "3. Why can we not simply take the limit as $J\\to\\infty$ and get a counterexample to de Finetti’s theorem?\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 134). CRC Press. Kindle Edition. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let's look at the case J=1. This is very much like the earlier problems, but with a continuous distribution.\n", + "So for example $P(x,y)=(1/2)P(x,N(1,1))P(y,N(-1,1))+(1/2)P(x,N(-1,1))P(y,N(1,1))$ and\n", + "$P(x,x)=P(x,N(1,1))P(x,N(-1,1))$ so it's exchangeable.\n", + "\n", + "However, in a mixture distribution there's no way to make sense of the requirement that the parameters are clustered into two groups. If we had a mixture then \n", + "$$P(x,y)=\\int P((x,y)|\\theta)p(\\theta) d\\theta=\\int P(x|\\theta)P(y|\\theta)p(\\theta)d\\theta$$" + ] + } + ], + "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 +}