Skip to content
Permalink
master
Switch branches/tags

Name already in use

A tag already exists with the provided branch name. Many Git commands accept both tag and branch names, so creating this branch may cause unexpected behavior. Are you sure you want to create this branch?

BDA/BDA 5.9.1-2-4-5 -Exchangeability and DeFinetti's Theorem.ipynb

Go to file
 
 
Cannot retrieve contributors at this time
267 lines (267 sloc) 34.8 KB
Raw Blame
Open in GitHub Desktop
{
"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": [
"<Figure size 432x288 with 2 Axes>"
]
},
"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",
"After a small amount of cheating (by looking at some solutions by Gelman) the suggestion is to look at the covariance of $y_1$ and $y_2$. Informally, they should have negative covariance because if $y_1$ is large, it suggests that it came from $N(1,1)$; but then $y_2$ comes from $N(-1,1)$ so it should be small.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"from scipy.stats import norm\n",
"y1=norm(loc=1,scale=1)\n",
"y2=norm(loc=-1,scale=1)\n",
"y1_sample=y_1.rvs(500)\n",
"y2_sample=y_2.rvs(500)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"A sample of our situation is $(y_1,y_2)$ or $(y_2,y_1)$ with equal probability. So the mean of each variable is zero. The covariance is $-1=E(y_1y_2)=E(y_1)E(y_2)$. The next problem (5.9.5) shows that mixtures of iid variables have positive covariances."
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Covariance= -1.0023058663082174\n"
]
}
],
"source": [
"cov=sum(y1_sample*y2_sample)/500\n",
"print(\"Covariance=\",cov)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In the general case, the joint probablity distribution can be written\n",
"$$\n",
"P(y_1,\\ldots,y_{2J})=(\\binom{2J}{J})^{-1}\\sum_{{S\\subset [2J]}\\atop{|S|=J}} P_{S}(y_1,\\ldots,y_{2J})\n",
"$$\n",
"where $$P_{S}(y_1,\\ldots,y_{2J})=\\prod_{i\\in S} P(y_i,N(1,1))\\prod_{j\\not\\in S} P(y_j,N(-1,1)).$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To understand the covariance of, say $y_1$ and $y_2$, we need to know how often they are chosen from the same distribution and how often they are chosen from different ones. That raises the combinatorial question of how many of the partitions of $[2J]$ have $y_1$ and $y_2$ together, and how many of them separate $y_1$ and $y_2$. To have them together, we first pick $y_1$ and $y_2$, and then choose $J-2$ additional elements from the remaining $2J-2$. So there are $\\binom{2J-2}{J-2}$ subsets of size $J$ that contain both $y_1$ and $y_2$. To split them, we pick $J-1$ elements from the $2J-2$ elements other than $y_1$ and $y_2$ and combine those with $y_1$ (for example) so there are $\\binom{2J-2}{J-1}$ sets that split them. \n",
"\n",
"When computing the covariance, the cases where $y_1$ and $y_2$ are together contribute $+1$, and the cases where they are split contribute $-1$. This gives the following:\n",
"$$\n",
"\\mathrm{cov}(y_1,y_2)=\\frac{2(\\binom{2J-2}{J-2}-\\binom{2J-2}{J-1})}{\\binom{2J}{J}}\n",
"$$\n",
"The two in the numerator comes from the fact that the number of partitions is $1/2$ of $\\binom{2J}{J}$. \n",
"\n",
"Some trial computations gives the explicit formula that the covariance is $-\\frac{1}{(2J-1)}$ and this goes to zero as $J\\to\\infty$.\n",
"\n",
"The next problem (5.9.5) shows that in a mixture, the correlations are non-negative, so this shows we don't have a mixture of iid variables."
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"1 -1.0\n",
"2 -3.0\n",
"3 -5.0\n",
"4 -7.0\n",
"5 -9.0\n",
"6 -11.0\n",
"7 -13.0\n",
"8 -15.0\n",
"9 -17.0\n"
]
}
],
"source": [
"# an illustration of the last point from the discussion above\n",
"from scipy.special import binom\n",
"for i in range(1,10):\n",
" print(i,(binom(2*i,i)/2/(binom(2*i-2,i-2)-binom(2*i-2,i-1))))\n",
" "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For the last point, I think the way to think of it is that $y_1$ is a combination of $\\binom{2J-1}{J-1}$ copies of $N(1,1)$ -- corresponding to the partitions in which $y_1$ is in the first half -- and $\\binom{2J-1}{J}$ -- corresponding to the partitions in which $y_1$ is in the second half. (Note that since $2J-1$ is odd, these numbers are actually equal). In the limit the correlation between different $y_i$'s drops to zero and so they become independent, and there's no contradiction to deFinetti's theorem."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Problem 5.9.5\n",
"\n",
"Suppose that the distribution of $\\theta=(\\theta_1,\\ldots,\\theta_{J})$ can be written as a mixture of independent and identically distributed components:\n",
"$$\n",
"p(\\theta)=\\int \\prod_{j=1}^{J} p(\\theta_{j}|\\phi)p(\\phi)d\\phi.\n",
"$$\n",
"Prove that the covariances $\\mathrm{cov}(\\theta_{i},\\theta_{j})$ are all non-negative.\n",
"\n",
"Here we apply the formula:\n",
"$$\n",
"\\mathrm{cov}(y_1,y_2)=E_{\\phi}(cov(y_1,y_2|\\phi))+\\mathrm{cov}_{\\phi}(E(y_1|\\phi),E(y_2|\\phi))\n",
"$$\n",
"The first term is zero (since $y_1$ and $y_2$ are independent, conditional on $\\phi$), and the second term is positive since $E(y_1|\\phi)=E(y_2|\\phi)=\\mu(\\phi)$ since the $y_1$ are identically distributed given $\\phi$;\n",
"thus this term is $\\mathrm{var}(\\mu(\\phi))\\ge 0$.\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.5"
}
},
"nbformat": 4,
"nbformat_minor": 2
}