BDA 5.9.8.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Discrete Mixture Models\n",
    "\n",
    "<font color='red'> This solution differs from the one published here:\n",
    "http://www.stat.columbia.edu/~gelman/book/solutions3.pdf so it's probably wrong.\n",
    "</font>\n",
    "\n",
    "Discrete mixture models: if $p_m(\\theta)$, for $m=1,\\ldots,M$ are conjugate prior densities for the sampling model $y|\\theta$, show that the class of finite mixture prior densities given by \n",
    "$$\n",
    "p(\\theta)=\\sum_{1}^{M} \\lambda_m p_m(\\theta)\n",
    "$$\n",
    "is also a conjugate class, where the $\\lambda_m$’s are nonnegative weights that sum to 1. This can provide a useful extension of the natural conjugate prior family to more flexible distributional forms. As an example, use the mixture form to create a bimodal prior density for a normal mean, that is thought to be near $1$, with a standard deviation of $0.5$, but has a small probability of being near $−1$, with the same standard deviation. If the variance of each observation $y_1,\\ldots,y_{10}$ is known to be $1$, and their observed mean is $y =−0.25$, derive your posterior distribution for the mean, making a sketch of both prior and posterior densities. Be careful: the prior and posterior mixture proportions are different.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's skip the theory part and look at the example.\n",
    "\n",
    "We have\n",
    "$$\n",
    "p(\\theta|y_1,\\ldots,y_{10})\\propto p(y_1,\\ldots,y_10|\\theta)p(\\theta)$$\n",
    "so\n",
    "$$\n",
    "p(\\theta|\\{y_{i}\\})\\propto \\sum \\lambda_{m}p(\\{y_{i}\\}|\\theta)p_{m}(\\theta)\n",
    "$$\n",
    "\n",
    "Each of the terms $p_{m}(\\theta)p(\\{y_{i}\\}|\\theta)$\n",
    "is equal to  $p_{m}(\\theta|\\{y_{i}\\})p_{m}(\\{y_{i}\\})$.\n",
    "\n",
    "Therefore the total posterior density is a weighted sum\n",
    "of the individual posteriors:\n",
    "\n",
    "$$p(\\theta|\\{y_{i}\\})=\\sum c_{m}p_{m}(\\theta|\\{y_{i}\\})$$\n",
    "where \n",
    "$$\n",
    "c_{m}=\\frac{\\lambda_m p_{m}(\\{y_{i}\\})}{\\sum_{m} \\lambda_m p_{m}(\\{y_{i}\\}}\n",
    "$$\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<font color=\"red\"> This is the part that was wrong -- I leave it here for historical purposes, but the correct part follows. </font>\n",
    "<p>\n",
    "In the special case under consideration, $p_1$ is normal with mean $-1$ and $\\sigma=.5$, $p_2$ is normal with mean $1$ and $\\sigma=.5$ and we can set $\\lambda_1=.1$ and $\\lambda_2=.9$.  The $p_m(\\{y_{i}\\})$ can be calculated from the $t$ distribution.  Drawing a sample of size $10$ from $p_1$ and getting a sample mean of $-.25$ and a sample variance of $1$ gives a $t$-statistics of \n",
    "$$\\frac{(\\overline{y}-\\mu)}{s/\\sqrt{N}}=\\frac{(-.25+1)}{1/\\sqrt{10}}=\\sqrt{10}(-.25+1)$$ in the first case and $\\sqrt{10}(-.25-1)$ in the second. \n",
    "<p>\n",
    "<font color=\"red\"> Now we continue with what is correct </font>\n",
    "\n",
    "We need to properly interpret $p_m(\\{y_{i}\\})$ and for that we should remember where it comes from.  We\n",
    "rewrote\n",
    "$$\n",
    "p(\\{y_{i}\\}|\\theta)p_{m}(\\theta)=p_m(\\{y_{i}\\},\\theta)=p_{m}(\\{y_{i}\\})p_{m}(\\theta|\\{y_{i}\\})\n",
    "$$\n",
    "We're dealing here with normal distributions.  On the left, the quadratic form that is the log-likelihood of the\n",
    "relevant bivariate normal is\n",
    "$$Q=\\frac{(\\overline{y}-\\theta)^2}{\\sigma^2}+\\frac{(\\theta-\\mu)^2}{\\tau^2}\n",
    "$$\n",
    "where, more specifically, $\\overline{y}=-.25$, $\\mu=\\pm 1$, $\\sigma^2=1/10$, and $\\tau^2=0.25=0.5^2$. \n",
    "The first term comes from $p_{m}(\\overline{y}|\\theta)$ and the second from $p_{m}(\\theta)$.\n",
    "\n",
    "Pure algebra (by expanding, writing $Q$ as a quadratic in $\\theta$, and completing the square) gives us\n",
    "the expression\n",
    "$$\n",
    "Q=\\frac{(\\theta-\\mu_{1})^2}{\\tau_1^2}+\\frac{(\\overline{y}-\\mu)^2}{\\sigma^2+\\tau^2}\n",
    "$$\n",
    "where \n",
    "$$\n",
    "\\mu_{1}=\\frac{\\mu/\\tau^2+\\overline{y}/\\sigma^2}{\\tau_1^2}\n",
    "$$\n",
    "and \n",
    "$$\n",
    "\\frac{1}{\\tau_1^2}=\\frac{1}{\\sigma^2}+\\frac{1}{\\tau^2}\n",
    "$$\n",
    "\n",
    "In the context of the problem under discussion, the two terms tell us that\n",
    "$$\n",
    "\\theta\\sim N(\\mu_1,\\tau_1^2)\n",
    "$$\n",
    "where $\\tau_1^2=1/(1/.1+1/.25)=1/14$ and $\\mu_1=(-.25/.1\\pm 1/.25)/14$ giving $-.46$ and $.11$. \n",
    "and $p_{m}(\\{y_{i}\\})=N(-.25,\\pm 1, .1+.25)$.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.30191827840729224 0.07235502834102417\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "from scipy.stats import norm, t\n",
    "import matplotlib.pyplot as plt\n",
    "p2y=norm.pdf(-.25,1,np.sqrt(.35))\n",
    "p1y=norm.pdf(-.25,-1,np.sqrt(.35))\n",
    "lambda1,lambda2=(.1,.9)\n",
    "print(p1y,p2y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "wt1=lambda1*p1y/(lambda1*p1y+lambda2*p2y)\n",
    "wt2=lambda2*p2y/(lambda1*p1y+lambda2*p2y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.3167705292212528 0.6832294707787472\n"
     ]
    }
   ],
   "source": [
    "print(wt1, wt2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-0.4642857142857143 0.10714285714285714\n"
     ]
    }
   ],
   "source": [
    "post_mean1,post_var1=posterior(-1,.25,-.25,1,10)\n",
    "post_mean2,post_var2=posterior(1,.25,-.25,1,10)\n",
    "print(post_mean1,post_mean2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl8VOX1+PHPIQEiCEUhbmzBCAY1LBpFBbcviKiIsoobWPdat9b2V2y/ar+0WNu6oNVaF3DBBQUj4lYEVAJIUBBUZFGUKGGRCLKJEJI8vz9OBkIIZJLcmXvnznm/XvOaZHLn3jNZTp4597nnEeccxhhjwqWe3wEYY4zxniV3Y4wJIUvuxhgTQpbcjTEmhCy5G2NMCFlyN8aYELLkbowxIWTJ3RhjQsiSuzHGhFCqXwdu0aKFy8jI8OvwxhiTkObPn/+Dcy69uu18S+4ZGRnMmzfPr8MbY0xCEpFvo9nOyjLGGBNCltyNMSaELLkbY0wI+VZzN8aE286dOyksLGT79u1+h5KQ0tLSaNWqFfXr16/V8y25G2NiorCwkCZNmpCRkYGI+B1OQnHOsX79egoLC2nXrl2t9mFlGWNMTGzfvp3mzZtbYq8FEaF58+Z1etdjyd0YEzOW2Guvrt87S+7G7ENpKUycCA88AKtW+R2NMTVjyd2YKpSVwSWXwODBcPvtcNxx8NlnfkdlYuWuu+5i2rRpfofhKUvuxlRh9GiYMAHuuQcWLYJGjeCyy6C42O/IjNdKS0sZOXIkvXr1qtFzgs6SuzGVbNgAI0fCeefBiBFw7LHwyCOa5MeP9zs6UxMFBQVkZWUxfPhwOnXqxKBBg9i2bRsZGRmMHDmSHj16MGHCBK688komTpwIwPTp0+natSvZ2dlcddVV7NixA2Cv5wSdTYU0ppKHHoLNm+HeeyFyTuuii+CYY+D+++GKK3Y/bqJ0222wcKG3++zSRd9iVWPZsmWMGTOG7t27c9VVV/Hvf/8b0Hnks2bNAuC///0voDN8rrzySqZPn06HDh0YNmwYjz32GLfddttezwk6G7kbU0FJCYwZA+eeC9nZux8X0fz02Wdg/e4SS+vWrenevTsAl19++a7kfPHFF++17bJly2jXrh0dOnQAYPjw4eTl5e36elXPCSobuRtTwZQpOjPmX//a+2sDB8KNN2ot/sQT4x9bQotihB0rlacURj5v3LjxXts65/a7r6qeE1Q2cjemgrFj4ZBDoG/fvb928MFw9tma3KvJASZAvvvuO+bMmQPASy+9RI8ePfa5bVZWFgUFBSxfvhyAcePGccYZZ8QlTq9Zcjem3LZt8M47Ov1xX+08LroICgpg2bK4hmbqoGPHjjz77LN06tSJDRs28Ktf/Wqf26alpfH0008zePBgsrOzqVevHjfccEMco/WOlWWMKTdtGvz8sybwfTn77N3bZmXFJy5TN/Xq1eM///nPHo8VFBTs8fkzzzyz6+OePXuyYMGCvfZT+TlBZyN3Y8pNmgS/+AXs7114u3Zw5JEwdWr84jKmNqpN7iIyVkTWiciifXxdRORhEVkuIp+JyPHeh2lMbJWVwZtv6tz26jqs9uwJM2boc0ywZWRksGhRlakr9KIZuT8D9NnP188F2pffrgMeq3tYxsTXp59CUZFOgazOqafCpk1WdzfBVm1yd87lARv2s8mFwHNO5QPNRORwrwI0Jh6mT9f7nj2r3/bkk/U+Pz928RhTV17U3FsCKyt8Xlj+mDEJY/p0PUF6xBHVb9uhg9bm586NfVzG1JYXyb2qC7GrnAUsIteJyDwRmVdUVOTBoY2pu+JiyMuLbtQOUK8edOtmI3cTbF4k90KgdYXPWwGrq9rQOfeEcy7HOZeTnp7uwaGNqbu5c3WOew2aAtKtG3z+OWzdGru4TDBMmjSJxYsX1/h5kydP5t57741BRNHxIrlPBoaVz5o5GdjknFvjwX6NiYv33tPeMTW5EDEnR2fLWI/38KtNci8pKaFfv36MGDGiRs/xUrUXMYnIS8CZQAsRKQTuBuoDOOf+A7wNnAcsB7YBv/Q0QmNi7MMPtUnYQQdF/5zOnfX+s8909owJpoKCAvr06UO3bt1YsGABHTp04LnnnmPOnDn87ne/o6SkhBNPPJHHHnuMhg0bMmLECCZPnkxqaiq9e/dmwIABTJ48mRkzZvDXv/6VV199FYBf//rXFBUV0ahRI5588kmysrK48sorOfjgg1mwYAHHH3882dnZzJs3j0ceeYRvv/2Wq666iqKiItLT03n66adp06bNXs+5//77PXvt1SZ359wl1XzdAb/2LCJj4qi0FObM0YU4aqJNGz2p+umnsYkrbHzs+LtXy98HHniAxx9/fK+2vsOGDeO1115j6dKliAgbN26kWbNm9OvXj759+zJo0CBAr2D9z3/+Q/v27Zk7dy433ngj7733HgBffvkl06ZNIyUlZY+rXm+66SaGDRvG8OHDGTt2LLfccguTJk3a6zlesitUTVL74gvYsgXKO8JGTQQ6dbKyTCKo3PJ3+vTpVbb1bdq0KWlpaVxzzTXk5ubSqFGjvfa1detWPvzwQwYPHkyXLl24/vrrWbNmdxV68ODBVSbpOXPmcOmllwJwxRVX7NETfl/PqSvrLWOS2uzZel+b0kqnTvDcc1p7r2fDpP3ysePvXi1/9yU1NZWPPvqI6dOnM378eB555JFdI/KIsrIymjVrxsJ9vA2JtiVwxZhi1UbYfiVNUvvwQzjsMO0ZU1OdO+uoP8H6SSWdyi1/e/XqVWVb361bt7Jp0ybOO+88Ro8evSuBN2nShC1btgDQtGlT2rVrt2uZPeccn0ZRmzv11FMZX75G4wsvvLDftsNeseRuktrs2Tpqr82yeZ066b2VZoKtcsvf3/zmN1W29d2yZQt9+/alU6dOnHHGGTz44IMADB06lH/+85907dqVr7/+mhdeeIExY8bQuXNnjj32WF5//fVqY3j44Yd5+umn6dSpE+PGjeOhhx6K9ctGqlt5JFZycnLcPFuvzPhozRq9IvW+++D222v+/C1boGlTuOceuOMO7+NLdEuWLKFjx46+xlBQUEDfvn0TtnlYVd9DEZnvnMup7rk2cjdJq/ydeo1PpkY0aaL/HJYu9S4mY7xiyd0krdmzoWFD6Nq19vvIyrLukEFmLX+NSUIffqgLXTdsWPt9ZGXpyN3WVK2aX2XfMKjr986Su0lKP/8M8+fX/erSrCzt7f79997EFSZpaWmsX7/eEnwtOOdYv349aWlptd6HzXM3SWn+fNi505vkDjp6P+ywuscVJq1ataKwsBDrAFs7aWlptGrVqtbPt+RuklKkXW9k4Y3aiiT3ZcvgzDPrtq+wqV+/Pu1qcwGB8YSVZUxSys/XC5cOPbRu+2nZEho1shkzJngsuZuklJ9f91E7aNuBo4+25G6Cx5K7STqFhbBqFZxyijf7s+RugsiSu0k6XtXbI446Cr77TpfrMyYoLLmbpJOfr3PbIwtu1FVmpnaG/PZbb/ZnjBcsuZukk58PJ5wADRp4s7/MTL3/+mtv9meMFyy5m6RSXKxz3L0qyYAldxNMltxNUvnsM9i+3dvkfvjhcMABltxNsFhyN0nF65OpoL3gjzzSkrsJFkvuJqnk52ub3jpc1V2lzExL7iZYLLmbpBK5eKk2Ky/tT2YmfPONdYc0wWHJ3SSNdet0dO1lSSYiM1M7Ta5d6/2+jakNS+4macydq/deXZlakc2YMUFjyd0kjfx8SE2F44/3ft+W3E3QWHI3SSM/X69KbdTI+323batNxCy5m6Cw5G6SQkkJfPRRbOrtoFe7tmljyd0EhyV3kxQWLoStW+G002J3jHbtYMWK2O3fmJqIKrmLSB8RWSYiy0VkRBVfbyMi74vIAhH5TETO8z5UY2pv1iy979EjdsfIyLDmYSY4qk3uIpICPAqcCxwDXCIix1Ta7H+BV5xzXYGhwL+9DtSYupg5U0fWLVvG7hht28Lq1bBjR+yOYUy0ohm5nwQsd85945wrBsYDF1baxgFNyz/+BbDauxCNqRvndOQey1E76MgdYOXK2B7HmGhEk9xbAhV/XQvLH6voz8DlIlIIvA3c7El0xnjgq6/0AqZY1ttBR+4ABQWxPY4x0YgmuVd1oXbli6wvAZ5xzrUCzgPGiche+xaR60RknojMKyoqqnm0xtTCzJl6H6+Ru9XdTRBEk9wLgdYVPm/F3mWXq4FXAJxzc4A0oEXlHTnnnnDO5TjnctLT02sXsTE1NGsWNG8OWVmxPU7LljrX3UbuJgiiSe4fA+1FpJ2INEBPmE6utM13QE8AEemIJncbmptAmDlTR+1eNwurrH597TZpI3cTBNUmd+dcCXATMAVYgs6K+UJERopIv/LNbgeuFZFPgZeAK52z/njGf6tW6YVFsa63R7RtayN3Ewyp0WzknHsbPVFa8bG7Kny8GOjubWjG1N306Xrfs2d8jpeRATNmxOdYxuyPXaFqQm36dGjRAjp1is/x2raFwkLYuTM+xzNmXyy5m9ByTpP7WWfpic54yMiAsjItBxnjJ0vuJrS+/FKTbLxKMmBz3U1wWHI3oRXvejvYXHcTHJbcTWhNn65teCMLacRD69Y65dJG7sZvltxNKJWWwvvv66g91vPbK2rYEA4/3Ebuxn+W3E0ozZ8PP/4IvXrF/9gZGTZyN/6z5G5C6a23dIbMOefE/9ht29rI3fjPkrsJpbfe0iX1mjeP/7EzMuC777Q0ZIxfLLmb0FmzRssyffv6c/y2bXXN1tV+rGrgHLzxBlx4IRxyCDRpAscfD/feq3UqkzQsuZvQebu8Ucb55/tz/Mhc9+++i/OB162Dc8+Ffv3g44/1/uqroXFjuOMObYv55ptxDsr4xZK7CZ233tLujNnZ/hw/ktzjWnf/6ivo1g3y8uDhh/U/y1NPwejR2hZz3jydxtOvHzz2WBwDM36JqnGYMYlixw6YOhUuvzy+UyAratNG7+OW3Neu1TPHW7dq17ITT9x7mxNOgDlzYMgQuPFGaNQIhg+PU4DGDzZyN6HywQea4/wqyYBWQZo3j1Ny37FDR+Pff6/1qKoSe8QBB0Burk7+v+YamD07DgEav1hyN6HyyivQtKk/89srats2TjX3O+/U+vrzz+8/sUfUrw8TJ2qAl10GmzbFPkbjC0vuJjSKi3VgeuGFkJbmbyxxmev+wQdw331www3Qv3/0z2vWTP8ZFBbCTTfFLDzjL0vuJjTefRc2boSLL/Y7kt3JPWbrkRUXw69+BUceCfffX/Pnn3wy/PGPmuTff9/7+IzvLLmb0Hj5ZTjoIDj7bL8j0ZOqP/0EGzbE6AAPPQRLl+rMmEaNarePO+6Adu109G6ri4SOJXcTCtu3w+uvw4AB0KCB39HEeK57URGMHAkXXADnnVf7/RxwgE6VXLwYHn/cu/hMIFhyN6HwzjuwZUswSjIQ47nu//gHbNum93V1wQVw+ukwahT8/HPd92cCw5K7CYWXX4b0dF1SLwhiltzXroVHH9WZLllZdd+fiL4LWLvWLm4KGUvuJuH99JO2Uxk4EFIDclle8+Za9fA8uf/tb3oy9a67vNvnGWfo3Pd779WLBEwoWHI3Ce/NN7VKEZSSDOiA2PO57j/8AE8+CcOGwVFHebhj4P/+T2v5Y8d6u1/jG0vuJuG9/LK2TTntNL8j2ZPnc90ff1zr4rff7uFOy3XvDqecoidYrVdxKFhyNwlt82a96n7wYEhJ8TuaPXma3HfsgEcegT594NhjPdppJbffDitWwKRJsdm/iStL7iahTZ6seS9IJZmINm200rFtmwc7Gz9eT3r+9rce7GwfLrqo9hdFmcCx5G4S2vjx0Lq1XnAZNJEZMytX1nFHzsGDD8Jxx8W2aU5KCtx6q3aPnD8/dscxcWHJ3SSsH3/UlgNDhuh6qUHj2XTI/Hz49FO45ZbY9zEeNkyn+dhFTQkvqj8JEekjIstEZLmIjNjHNkNEZLGIfCEiL3obpjF7e+01vWo+iCUZ8DC5P/mk9hEeOrTOMVWrWTM9zosv6gkNk7CqTe4ikgI8CpwLHANcIiLHVNqmPXAH0N05dyxwWwxiNWYPL7+sJeKcHL8jqdoRR2ilo07JffNmfaGXXKLrocbD9dfrxQMv2hgtkUUzcj8JWO6c+8Y5VwyMBy6stM21wKPOuR8BnHPrvA3TmD0VFcH06Tpq92vFpeqkpkLLlnVM7i++qGdkr73Ws7iqddJJ0LmzlmZi1tbSxFo0yb0lUPGUUGH5YxV1ADqIyGwRyReRPl4FaExVcnN1OnZQSzIRdb6Q6cknoVOn6Bbi8IqIjt4XLrQTqwksmuRe1bio8r/zVKA9cCZwCfCUiDTba0ci14nIPBGZV1RUVNNYjdllwgTo0EHzXpDVaa77J5/o7dpr4//2ZOhQaNgQxo2L73GNZ6JJ7oVA6wqftwJWV7HN6865nc65FcAyNNnvwTn3hHMuxzmXk56eXtuYTZIrKtL1JQYPDm5JJqJNG13wqKSkFk9+9lntX3zZZZ7HVa2DDtKOkS++aL3eE1Q0yf1joL2ItBORBsBQYHKlbSYBZwGISAu0TPONl4EaE5GbC2VlmtyDrm1bLR+trjwcqk5JiU7iv+ACTbR+GDZM+9n897/+HN/USbXJ3TlXAtwETAGWAK84574QkZEi0q98synAehFZDLwP/N45tz5WQZvkNmECtG8f/JIM1GHRjmnTYN06f0btEX36QIsWVppJUFE1SHXOvQ28Xemxuyp87IDflt+MiZlISWbEiOCXZGDPue49etTgiS+8oHPO67LSUl3Vr69TMJ94Qq8Y8+sdhKmVAF7XZ8y+vfZa4pRkQFsjQA1Pqv70k77QwYP1pKafrrhCm/dMmOBvHKbGLLmbhDJhgrYy79zZ70ii07ixVjZqlNxff10T/OWXxyyuqOXkwNFH2wVNCciSu0kY69cnziyZimo81/3553XIX6M6ToyI6MUEeXmwZo3f0ZgasORuEsbbb+vMk/79/Y6kZmo0172oSLuhXXZZcLqhXXyxXqk6caLfkZgaCMhvjzHVe+MNOOwwOOEEvyOpmTZtNLlHdSX/a6/pf7B4NAmL1jHHaLvhV17xOxJTA5bcTUIoLtbp1n37BmdAG622bbU9zPpoJgcHdZ7nkCEwa5ZekWUSQoL9mZhklZcHW7boNT2JJuq57j/8ENyTCpEmPlaaSRiW3E1CeOMNSEuL7UJEsRJ1X/dJk7QkM2hQzGOqsQ4doEsXbT9sEoIldxN4zmly/5//gUaN/I6m5tq00ftqk/uECZCZqUk0iIYM0VWhPFv128SSJXcTeEuWwIoViVmSAWjeXP8p7Tcnrl+vDeqDWJKJiJRm7IKmhGDJ3QTeG2/ofd++/sZRWyJRzHWPlGSCfOntkUfqVCVL7gnBkrsJvClTIDsbWrXyO5Laq3au+8SJ0K4ddO0at5hqZeBA+OgjmzWTACy5m0Dbtg1mz4bevf2OpG4ic92rtGGDdoEMckkmYsAAvZ80yd84TLUsuZtAy8vTOe5nn+13JHXTtq3OdPzppyq++Prr2r89yCWZiKOP1ouacnP9jsRUw5K7CbSpU3UxotNO8zuSuolMh1y5soovTpgAGRmJc+ntgAEwY4b+tzKBZcndBNrUqZrYE3EKZEX7nOu+caOWZAYNCn5JJmLAAO27PLnygmwmSCy5m8BaswY+/zzxSzKwn7nub76pa5QOHBj3mGqtSxd9p2GlmUCz5G4Ca9o0vQ9Dcj/iCEhJqSK55+ZCy5Zw0km+xFUrIjp6nzoVNm/2OxqzD5bcTWBNnaoLXQT1gs2aSE3VqZx7JPefftJuaP37J143tP799Uz3O+/4HYnZhwT7jTLJwjkduffqlXh5b1/2upBpyhT4+efd0wsTySmnwKGHWmkmwELyZ2PC5quvtOZ+1ll+R+Kdvea65+Zqb4JEnAqUkgIXXQRvvQXbt/sdjamCJXcTSDNn6v3pp/sbh5fatoVVq3RKO8XF2lfhwgu1ZpOIBgzQ0tLUqX5HYqpgyd0EUl4epKfrNTNh0batto9ZtQptErZ5c2KWZCLOPBOaNbPSTEBZcjeBNHOmVisSZep3NPZYtCM3F5o0gZ49fY2pTho00FadkyfrdE4TKJbcTeAUFmqL30QsRe/PrrnuK8q0N8v55+sKJIlswADtjZOX53ckphJL7iZwwlhvhwrJPe9bvXQ/kUsyEb176+XDr77qdySmEkvuJnBmztSKRefOfkfirUaNdPbgijlroWFDOPdcv0Oqu0aN9HVMmqQtCUxgWHI3gZOXB6eeqrPtwiYz0/H11w7OOQcOPNDvcLwxcKDOW83P9zsSU0FUyV1E+ojIMhFZLiIj9rPdIBFxIpLjXYgmmaxfD198Eb6STMSRv1jP1ztahaMkE3H++Xpy1UozgVJtcheRFOBR4FzgGOASETmmiu2aALcAc70O0iSP2bP1PmwnUyMyNy+kkFbs6J2gC8JWpWlTbQCUm6uXFptAiGbkfhKw3Dn3jXOuGBgPXFjFdn8B/gHY5Wqm1vLydBB44ol+RxIDzpG5fAqOeqzYdLDf0XhrwAAoKIAFC/yOxJSLJrm3BCouMVBY/tguItIVaO2ce3N/OxKR60RknojMKyoqqnGwJvxmzoRu3RJ/hmCVFi8m83t9a/L11z7H4rV+/fQkiV3QFBjRJPeqLiPZ9d5LROoBDwK3V7cj59wTzrkc51xOenp69FGapLB1K8yfH96SDLm5ZPINEMLk3qIFnHGG1d0DJJrkXgi0rvB5K2B1hc+bAMcBH4hIAXAyMNlOqpqays/Xy/PDejKV3FwOOfUoGjcOYXIHnTWzdCksWeJ3JIbokvvHQHsRaSciDYChwK71tZxzm5xzLZxzGc65DCAf6OecmxeTiE1ozZyp7X1POcXvSGLgm29g4UJk4AAyM0Oa3C+6SO9t9B4I1SZ351wJcBMwBVgCvOKc+0JERopIv1gHaJJHXp4uzNG0qd+RxMBrr+l9//5kZmquD50jjtALFKzuHghRzXN3zr3tnOvgnMt0zo0qf+wu59xeK+Q65860UbupqeJiLcuEtiTz6qv6n6tdu13JPZQXdA4YoDNmQvnfK7HYFaomEObP1zUfQnky9bvvYM4cGDwYgMxM2LEDVq+u5nmJKHJxVuSdivGNJXcTCJGmgqFM7hMn6n2F5A4hrbu3awddu1rdPQAsuZtAmDkTsrJ0gY7QeeUVTXjt2wMhT+6gs2bmzAnpW5PEYcnd+K60FGbNCumovaAA5s6FIUN2PdSmja6sF9rkbqWZQLDkbny3aBFs2hTSk6mVSjKgib1tW1i+3KeYYq1jR73ZrBlfWXI3vosszhHKkfsrr8AJJ+yuxZTr0AG+/NKnmOJhwACYMUMXJTG+sORufJeXB61b715jNDRWrICPP96jJBNx9NGa3EM5HRK07l5aaqUZH1lyN75yTkfuoSzJTJig9xVKMhFZWbBtm64XG0pduugJ5Jdf9juSpGXJ3fjq669h7doQl2ROOkmnB1Zy9NF6v2xZnGOKFxEYOhTef19/wCbuLLkbX0Xmt4du5P7113plVhUlGUiC5A6a3MvKdr+DMXFlyd34auZM7RableV3JB6LlCMGDaryy4cdpj10li6NY0zxdswx0KkTjB/vdyRJyZK78VVeHvTooe/iQ8M5GDdOX9g+zhKL6Og91CN30NH7hx/Ct9/6HUnSseRufLN6tfaXCl1J5pNPdEh++eX73SwpkvvFF+u9nViNO0vuxjeh7Sfz/PO6EOw+6u0RWVmwciX89FOc4vLDkUfqSWUrzcSdJXfjm7w8aNJEZ82FRkkJvPQS9O0LBx20300jJ1VDfTETaGlmwYIkeJsSLJbcjW9mzIDu3fVy/NCYNg2+/77akgwkyYwZ0HcwIjZ6jzNL7sYXRUWweLGuqRwq48bpiP2886rdtH17XVYw9EuOtmypP+gXXtCTzSYuLLkbX0T6yYQquW/ZopfbDxkCDRtWu3lamracWbQoDrH5bfhw+OorbQVs4sKSu/HFjBlwwAHaUys0XnsNfv4Zrrgi6qdkZ8Pnn8cwpqAYOBAaNYJnnvE7kqRhyd34YsYMXUu5QQO/I/HQmDE6FD/11Kifkp2trX+3bYthXEHQpIle0PXyy/oP0MScJXcTdz/+CJ99FrL57cuW6fSfa66p0RVZ2dlahl68OIaxBcXw4bB5M0ya5HckScGSu4m7WbM0oYWq3v7UUzrt58ora/S07Gy9T4rSzJln6jJUVpqJC0vuJu7y8rQc062b35F4pLgYnn0WLrhAm8bUQGamnlhNiuRerx4MG6bTRVet8jua0LPkbuJuxgxN7Glpfkfikddf17md115b46empGh/raRI7qClmbIynTJqYsqSu4mrzZu19UqoSjJPPqnlht69a/X07OwkmQ4JcNRR2lBt7NgQL0MVDJbcTVzNmKGrr/Xs6XckHlmxAqZOhauu0mF4LWRn63oWSbPc6PXX65z3997zO5JQs+Ru4mraNJ3ffsopfkfikUcf1aR+9dW13kXkpOpnn3kUU9ANGgTNm8Njj/kdSahFldxFpI+ILBOR5SIyooqv/1ZEFovIZyIyXUTCttSx8ci0adoFMooLOINv61adJTNoELRqVevdRBqnLVjgUVxBl5am73Ref91OrMZQtcldRFKAR4FzgWOAS0TkmEqbLQBynHOdgInAP7wO1CS+1at1PnevXn5H4pFnn4VNm+DWW+u0m0MOgdatYd48j+JKBNdfr/W5p57yO5LQimbkfhKw3Dn3jXOuGBgPXFhxA+fc+865yDV2+UDthzEmtKZP1/tQJPeyMnj4YTjxRDj55DrvLidHl1xNGpmZcM45ejK6pMTvaEIpmuTeElhZ4fPC8sf25WrgnboEZcJp2jQttXbu7HckHpgyRRux33abJ2sE5uToOcaNGz2ILVH86ldalpk82e9IQima5F7Vb26VfTtF5HIgB/jnPr5+nYjME5F5RUVF0UdpEp5zmtx79tRrWRLeAw/A4YfvcwHsmsrJ0ftPPvFkd4nh/PN1jdnRo/2OJJSi+TMrBFpX+LwVsLryRiLSC/gT0M85t6OqHTnnnnDO5Thv0bnzAAAQ3klEQVTnctLT02sTr0lQS5ZozT0UUyA/+kj/U/3mN551Pot0x0yquntqqp6vmDlTv6fGU9Ek94+B9iLSTkQaAEOBPd5HiUhX4HE0sa/zPkyT6N5+W+/79PE3Dk+MGqULctxwg2e7bN4c2rVLsuQO2mjtF7+A++/3O5LQqTa5O+dKgJuAKcAS4BXn3BciMlJE+pVv9k/gQGCCiCwUESuimT289ZbO527Txu9I6ujzz7VGfOut2sbWQzk5MHeup7sMviZNdObMxIlQUOB3NKESVfXTOfe2c66Dcy7TOTeq/LG7nHOTyz/u5Zw71DnXpfzWb/97NMlk0ybtBHn++X5H4oF77oEDD4Sbb/Z81927w3ffwcqV1W8bKjffrCdiHnrI70hCJQyntkzAvfuuznZL+OT+6ae62MRNN8HBB3u++x499H72bM93HWytWsEll+i0SJto4RlL7ibm3npLS9QeTAf31x13aH34//2/mOy+c2d9UzBrVkx2H2x//KOu0HTffX5HEhqW3E1MlZbCO+/oidTUVL+jqYMZM/SF3HGH/qeKgdRU7bmTlMk9K0tH7488AutsToYXLLmbmJo9W/9W+yXyWRjnYMQIOOKImNTaK+rRQxuIbdoU08ME0513wvbt8M8qL5MxNWTJ3cTUxInaJ6pvX78jqYNx4yA/H/7yF21pGUM9euj/kg8/jOlhgunoo+HSS7XT5vff+x1NwrPkbmKmrAxefRXOPVdryQlp40b4/e/1hEEN10etjZNP1o6Z06bF/FDBdOedumzhyJF+R5LwLLmbmJkzR69KHTzY70jq4K67dBWNRx+NS9+ERo20JfKUKTE/VDB16KAXhz3+uLYQNbVmyd3EzIQJOgpN2JJMfr4m9RtugOOPj9the/eGL76AwsK4HTJY/vxnfat3++1+R5LQLLmbmCgt1eTep4/nF3LGx7ZtMGyYzsG+5564Hvqcc/T+3XfjetjgaNFCyzP//a/eTK1YcjcxMXWqlmSGDfM7klr6wx+0B+8zz+jc9jjKztaGk0lbmgG9UCwzU1sq76iyD6GphiV3ExPPPqsXcSbkValvvKHzrX/zGzjrrLgfXkRH7+++q+cWk1LDhvozWLYs7u+cwsKSu/Hcxo3w2ms6qy3h1kr96iu4/HLtwetjUhkwQL+PkdWrklKfPnDZZfC3v+lJCFMjltyN555/Xt9JDx/udyQ1tGkT9O8P9evrHM60NN9C6d0bmjbV8xZJ7cEH9RtxzTW2HF8NWXI3nior03fTJ564e3WhhLB9O1x4oS6d9/LLukKQjxo21Kt6J02CnTt9DcVf6en6C5WfD3/9q9/RJBRL7sZT06ZpmfSWW/yOpAZ27IChQ7V/zLPPBma5qMGD4ccfk7w0A/qzGTZMrxCeOdPvaBKGJXfjqYcegkMPTaALl7Zt0xH766/Dv/6lzasC4pxz9KT02LF+RxIAjzyiS1Vddpm1BY5SIvfpM/HmHKxfDytW6Ko569fD5s1aq96xgwXr2/D227cw8qz3aThmib6lPuQQLXG0aRO8lbELC3WB648/hqeegquv9juiPTRsqOct/vUvbbVy6KF+R+SjJk20XNajBwwcqG8RPVq/NqzEOefLgXNycty8pFswMoE4p6taf/wxfPIJzJ+v7Qq3bNl725QUaNiQ/ttf5P2yMyggg2ZUamuYlgbt22tzqOOP312Ub9YsPq+nsqlTdVbMtm3aGOyii/yJoxpLl0LHjjpx5447/I4mAF58UUfv116rLQpE/I4o7kRkvnOu2jNaNnI3u333nSa96dPhvfd2d+Zr3Bi6dNG651FH6dvjjAwdlTdtCo0asfBTYVJXuPtuaPa/P+ioft063ceKFVqI//JLWLBAW0VGHHUUdOumHbNOOQU6ddLZKrHyww96gdLYsdpD/IMPNHsGVFYWnHmm5rHf/z7Be+J74dJLYdEinR7Zpg387//6HVFwOed8uZ1wwgnOBMC33zp3333OnXSSczped+6ww5y79FLnnnrKucWLnSsp2e8uysqcO+ss5w4+2LkNG6I45oYNzk2d6tw99zjXv79zhx+++9hpac716OHc73/v3KuvOrdqlTevc80a5+6807kDD3QuJcW5P/zBuW3bvNl3jE2apN+aceP8jiQgSkudGzZMvymjR/sdTdwB81wUOdbKMslo5UqdQD1hgk4xA71oZ8gQ7fLVsWON3u6++qqWrh99FG68sRbxOKf17zlzNJ78fC0DRS7PbN1aR/Unn6y3rl2rn4Me2ee778Kbb+qtpEQDHTky0KP1ysrK9I3Tzp06aE1J8TuiACgpgYsvhtxcGD0abr3V74jiJtqyjCX3ZFFYqOWQV17RJApa+x48WG+ZmbXa7aZN2gvloIM0H3tWNtixAxYu3J3s58yBb7/Vr6WmakOvVq2gZUstGzVsqNnvxx+1FLRkiZaGQLcZOhSuu05byiagiRP1xxTA877+2bFDyzS5ubpS1j33JEUN3pK7gVWrNCtMmKDr3YEOAYcM0Uxx1FF1PsSwYXqOa/ZsLZ3H1Jo1MHcuzJunib6wUF/jzz/rH3pqqv6XadFCi9XHHacF6+OOS/g/eufg9NN3n7rw6zx04JSW6tvFJ57Q3+kxYxK0DWn0ok3uVnMPm1WrnHv4Ya1bR+rYnTs7N2qUc19+6emhxo3T3d99t6e7Nfswf75z9eo598tf+h1JwJSVOff3v+s3p2NH5xYt8juimCLKmrsl9zAoKHDu/vudO/XU3Qk9O9u5v/zFuaVLY3LIWbOca9DAudNPd664OCaHMFX405/0xzthgt+RBND06c6lp+sv5qhRof3FjDa5W1kmETmnXfLefFPPZka+j1266AnDgQO1LBEjn3yija0OOkjL4c2bx+xQppKdO3UZvs8/1yvx47hAVGL4/nu4+WYtRR57LPzjH7qIb4KX5SqKtiwTsEsGzT5t3qwnjq69Vuf3ZmfrVS0i8Pe/w/LlOof8T3+KaWJ//31tcX7ggbpIjiX2+KpfX9spN2+uHXE/+cTviALm0EN10kBurjaDO/987RU0daoOipKIjdyDat06PUs5c6beFizQk0dNm0KvXjoaOeccnSYYBzt36iDorrt0wsm778bt0KYKX34JZ5+tk4PGjdP2OKaS4mI90TpqFKxdqyfWr79eJxQccojf0dWap7NlRKQP8BCQAjzlnLu30tcbAs8BJwDrgYudcwX726cl9wrWrtVpf5HbggX61ws6n7tbN+2p0bu3zveO5RWclZSVafXnj3/UStDQofr3EvIJCQmhsFDbAi9YAFdeqR1xW7b0O6oA2rFD+9KMHq3frJQUHc2ff74OlGp4XYffPEvuIpICfAmcDRQCHwOXOOcWV9jmRqCTc+4GERkK9HfOXby//SZVci8u1gS+ejV8842WUCK3r77SS+IjMjK0dn7KKZrQTzgh7ssZOadT7nJz4bnn9ON27fRv44ILEurvIPSKi+HPf4b77tOcNXSoJvru3a1VQZUWLYKXXtKa/Fdf6WOHH64DqOOP19uxx+o1FAH9BnqZ3E8B/uycO6f88zsAnHN/q7DNlPJt5ohIKrAWSHf72Xngk3tpqf7H37FDa3eRjyvetmzRq3git40bd3+8bp3Oy16zZvfFNBW1bq3zzDMz9e1ily7QuXNcJjCXlmroW7ZoKX/NGm0rU1Cgbxw+/lj/F4FeEHrzzTqFOI5vGEwNrVih7VZeegm2btV3Vqedpq16jj5axwwtWmijzoMPtp8loL/w06ZpH6X583e/WwZN7K1bw5FH6tuh9HT9Bka+ieU9lWjcWG+Rj9PS9L9sDEdAXib3QUAf59w15Z9fAXRzzt1UYZtF5dsUln/+dfk2P1S1T6h9ch/7yzzue/EIAJwTwFX4GMDhkMjDVHx1e2y/a5vy7fe1za69yj4f2+PjevVwUg/qpUBKPVy9VEipBykpuHop+oNPSa10TOL28c6d2gixKiKaCE48Ud84XHCBDmBM4vjpJ3jrLT3xPXOm5quqVnKqV0/fEEZuDRroY5GcJLLnbX+PhUZZKWzfAcU7oHgn7CzW+5ISKC2p+QnZPb5JArL747uvX8vFj5xWqzC97ApZ1Y+v8quMZhtE5DrgOoA2bdpEcei9tTiiIccdWrTrkJHvG8juXzQRRByRb6iUfz3y/dXfzgofA7Jr3pAg9UQTcr0UJDWSqFOQlHp7Pl6/PtKggf5lNKiPpKbu3l+F70iQPk5N1VFdxdvhh2vL9ZYtbUSX6Bo31vOFQ4bo5yUlWgksLNTq3w8/wIYNVb8RLSvT5+y+WGJ3PtvXY+GSAjQqv1WhtGT3N2tnecIvKd3zvrR09zeprOI3rWyPzw9qE/uTVtEk90Kg4ryIVsDqfWxTWF6W+QWwofKOnHNPAE+AjtxrE3C/Ud3oN6o2zzQm+aSm6uymBG2pEzCp5bfGfgcSlWjmuX8MtBeRdiLSABgKTK60zWQgstb9IOC9/dXbjTHGxFa1I3fnXImI3ARMQd+3jHXOfSEiI9HLYCcDY4BxIrIcHbEPjWXQxhhj9i+quT7OubeBtys9dleFj7cDibIksjHGhJ61HzDGmBCy5G6MMSFkyd0YY0LIkrsxxoSQJXdjjAkh31r+ikgR8G0tn94C2GdrgwRjryV4wvI6wF5LUNXltbR1zqVXt5Fvyb0uRGReNL0VEoG9luAJy+sAey1BFY/XYmUZY4wJIUvuxhgTQoma3J/wOwAP2WsJnrC8DrDXElQxfy0JWXM3xhizf4k6cjfGGLMfCZvcReQvIvKZiCwUkXdF5Ai/Y6otEfmniCwtfz2viUjs19qLAREZLCJfiEiZiCTkrAYR6SMiy0RkuYiM8Due2hKRsSKyrnyVtIQlIq1F5H0RWVL+u3Wr3zHVloikichHIvJp+Wv5v5geL1HLMiLS1Dm3ufzjW4BjnHM3+BxWrYhIb7QHfomI/B3AOfcHn8OqMRHpCJQBjwO/c84FeJHcvUWzGHyiEJHTga3Ac8654/yOp7ZE5HDgcOfcJyLSBJgPXJSgPxMBGjvntopIfWAWcKtzLj8Wx0vYkXsksZdrTBXL+iUK59y7zrmS8k/z0dWuEo5zbolzbpnfcdTBScBy59w3zrliYDxwoc8x1YpzLo8qVkNLNM65Nc65T8o/3gIsAVr6G1XtOLW1/NP65beY5a2ETe4AIjJKRFYClwF3Vbd9grgKeMfvIJJUS2Blhc8LSdBEEkYikgF0Beb6G0ntiUiKiCwE1gFTnXMxey2BTu4iMk1EFlVxuxDAOfcn51xr4AXgJn+j3b/qXkv5Nn8CStDXE0jRvI4EFtVC7yb+RORA4FXgtkrv2hOKc67UOdcFfXd+kojErGQW1UpMfnHO9Ypy0xeBt4C7YxhOnVT3WkRkONAX6Bnk9Wdr8DNJRNEsBm/irLw+/SrwgnMu1+94vOCc2ygiHwB9gJic9A70yH1/RKR9hU/7AUv9iqWuRKQP8Aegn3Num9/xJLFoFoM3cVR+EnIMsMQ594Df8dSFiKRHZsKJyAFAL2KYtxJ5tsyrwNHo7IxvgRucc6v8jap2yhcWbwisL38oPxFn/ohIf+BfQDqwEVjonDvH36hqRkTOA0azezH4UT6HVCsi8hJwJtp98HvgbufcGF+DqgUR6QHMBD5H/9YB/li+rnNCEZFOwLPo71Y94BXn3MiYHS9Rk7sxxph9S9iyjDHGmH2z5G6MMSFkyd0YY0LIkrsxxoSQJXdjjAkhS+7GGBNCltyNMSaELLkbY0wI/X/nkFfdAZY+dwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x=np.linspace(-3,3,1000)\n",
    "y1=lambda1*norm.pdf(x,-1,.5)+lambda2*norm.pdf(x,1,.5)\n",
    "y2=wt1*norm.pdf(x,post_mean1,np.sqrt(post_var1))+wt2*norm.pdf(x,post_mean2,np.sqrt(post_var2))\n",
    "fig,ax=plt.subplots(1)\n",
    "ax.plot(x,y1,color='red',label='prior')\n",
    "ax.plot(x,y2,color='blue',label='posterior')\n",
    "ax.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This now matches the solution published on the web!"
   ]
  },
  {
   "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
}