{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Unsupervised Learning\n",
    "\n",
    "In unsupervised learning, the task is to infer hidden structure from\n",
    "unlabeled data, comprised of training examples $\\{x_n\\}$.\n",
    "\n",
    "We demonstrate with an example in Edward. A webpage version is available at\n",
    "http://edwardlib.org/tutorials/unsupervised."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from __future__ import absolute_import\n",
    "from __future__ import division\n",
    "from __future__ import print_function\n",
    "\n",
    "import edward as ed\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.cm as cm\n",
    "import numpy as np\n",
    "import six\n",
    "import tensorflow as tf\n",
    "\n",
    "from edward.models import (\n",
    "    Categorical, Dirichlet, Empirical, InverseGamma,\n",
    "    MultivariateNormalDiag, Normal, ParamMixture)\n",
    "\n",
    "plt.style.use('ggplot')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data\n",
    "\n",
    "Use a simulated data set of 2-dimensional data points\n",
    "$\\mathbf{x}_n\\in\\mathbb{R}^2$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def build_toy_dataset(N):\n",
    "  pi = np.array([0.4, 0.6])\n",
    "  mus = [[1, 1], [-1, -1]]\n",
    "  stds = [[0.1, 0.1], [0.1, 0.1]]\n",
    "  x = np.zeros((N, 2), dtype=np.float32)\n",
    "  for n in range(N):\n",
    "    k = np.argmax(np.random.multinomial(1, pi))\n",
    "    x[n, :] = np.random.multivariate_normal(mus[k], np.diag(stds[k]))\n",
    "\n",
    "  return x\n",
    "\n",
    "\n",
    "N = 500  # number of data points\n",
    "K = 2  # number of components\n",
    "D = 2  # dimensionality of data\n",
    "ed.set_seed(42)\n",
    "\n",
    "x_train = build_toy_dataset(N)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We visualize the generated data points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAEICAYAAABCnX+uAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl4VOX5N/DvmUkwZCOYQmhQEiNJWGSrmgK2Llilv/Jq\nXXvhghLFhUUQUJDtVbZi+QHKZlBBEClaKUYrtF5W8RXaS0XFYlkTCAQIICEw2ZhAZs55/zhMSDJn\nmTlzZjkz3891eQGZM2eeE/A+T+5zP/cjSJIkgYiILMsW7gEQEVFgGMiJiCyOgZyIyOIYyImILI6B\nnIjI4hjIiYgsjoGcdLlcLrzxxhsYPHgwrrnmGgwYMAATJkxAeXl50zFLly7FPffcE5TP/+CDD/DL\nX/7S5+O3b9+OvXv3Gv68devWYdCgQUE7PtDx+WPfvn34+uuvQ/JZFD4M5KRr0aJF+Otf/4qpU6fi\nk08+weuvv476+no89NBDqKmpAQA89thjWLVqVZhHKhs2bBhOnjwZ7mGoCuX4Ro0ahQMHDoTksyh8\nGMhJ18aNGzF69GjcdNNNuOKKK9C7d2+8+uqrqK+vx2effQYASEpKQvv27cM8UqLYxEBOugRBwNdf\nfw2Xy9X0tbZt2+LDDz/EbbfdBqBlauWbb77BDTfcgM2bN+PGG29Ev379MGvWLJSXl+Phhx9G7969\ncd9996GsrAyAcurkhRdewNixYxXHs23bNvzhD39A79690adPHzzyyCM4dOgQADSlOJ5++mm88MIL\nAIBDhw7h8ccfR58+fXDLLbdg/vz5uHDhQtP5fvzxx6bzPfDAA7qzZb3j/R2f1vEAsHz5ctx0003o\n1asX7rrrLnz55ZdNr1VWVmLcuHHo168ffvWrX2HatGmora0FIM/8KyoqMHv2bAwbNkzzmsjaGMhJ\n12OPPYYPPvgAN910E1544QUUFxfj9OnTyMrKQkpKiuJ7HA4HNmzYgFWrVmHWrFn485//jMLCQgwf\nPhzvv/8+Lly4gEWLFvk9lmPHjmHkyJH47W9/i82bN+Ptt99GdXU15s+fDwD461//CgCYP38+pk2b\nhvPnz+Pxxx9HdnY2iouLMX/+fGzbtg1z5sxpGucTTzyBbt26obi4GPfccw/efvtt1c/XO97f8ekd\n/9lnn2HVqlWYN28e/vGPf+Cmm27Cs88+i7q6OgDAM888AwD4y1/+gqKiIhw5cgTjx48HIN9cO3Xq\nhAkTJmDp0qV+f6/JOuLCPQCKfE8++SSysrLw7rvvYtOmTSguLobdbscDDzyAqVOnwm63e73H5XLh\n+eefR25uLnJzczFv3jzceuut+M1vfgMAuPPOO7Fhwwa/x+JyuTB58uSmGeaVV16J3//+91i/fj0A\n4PLLLwcApKamIiUlBRs3bkR8fDxmzJgBAMjJycHMmTPx0EMPYdKkSfj73//e9Hp8fDyuvvpq/Pe/\n/8W//vUvxc/XO97f8VVVVWkef+zYMcTHxyMzMxNXXHEFRo8ejeuvvx5xcXH4+uuvsX//fqxduxZt\n2rQBACxYsAA33ngjSktLkZubC7vdjqSkJKSlpfn9vSbrYCAnnwwePBiDBw9GfX09tm/fjg8//BDr\n1q1DRkYGnnzyScX3XHnllU2/T0hI8PpzY2Oj3+PIzs5G27Zt8eabb6K0tBSHDh3C3r170bFjR8Xj\nDxw4gKNHj6Jfv35NX5MkCaIo4vDhwygtLUV+fj7i4+ObXu/du7dqINc73t/x6R1/xx134L333sPg\nwYPRo0cP3Hzzzbj33nuRkJCAAwcOwOl0Klb0lJWVITc3V+e7SdGCgZw07du3Dxs2bGia0SYlJeGW\nW27BLbfcgvHjx2Pbtm2qgbz1TN1mU87kCYLg9bXm+fjm9u/fjwceeAADBw7E9ddfj/vuuw87d+7E\nu+++q3i8y+VC3759MW/ePK/XMjIyIAgCWjcAbR6klcaqdby/49M7Pj09HZs3b8Y333yDL7/8Eh9/\n/DHeeecdrFu3Di6XC5mZmVi9erXXedPT01WvgaIPc+SkSRRFrFu3Dtu3b/d6LTk52ZRKlfj4eDid\nzhYB8tixY4rHvv/+++jevTuWLVuGRx99FAUFBaioqPAKrh5XX301ysvL0alTJ2RlZSErKwtnz57F\n/Pnz0djYiLy8POzbt6/Fw889e/aojlXveH/Hp3f8l19+iXfeeQcDBw7ElClT8MknnyAlJQVbt27F\n1VdfjVOnTiEpKanp2uLi4jBv3jycOXNG9Roo+jCQk6YePXrg9ttvx7hx47BhwwYcOXIEe/fuxapV\nq/Dxxx+jsLAw4M+45pprcP78ebz55ps4evQo3njjDdVgmpGRgbKyMnz33Xc4evQo1qxZgw0bNrQI\nrImJiSgtLYXD4cCdd94Jm82GyZMno6SkBDt27MCUKVPQ2NiIlJQUDBkyBDabDdOnT8fBgwfx8ccf\na+bu9Y73d3x6x0uShIULF2LTpk2oqKjAp59+isrKSvTq1Qs33HADcnNzMX78eOzatQt79+7FxIkT\nUVFRgc6dOwOQf4I6ePAgqqqqAvo7osjGQE66Fi5ciGHDhmHNmjW444478OCDD2Lbtm1YuXJli9yz\nUdnZ2ZgyZQrWrl2LO++8EwcPHsSjjz6qeOywYcMwcOBAPP3007jnnnvwxRdf4KWXXkJVVVVTGeDj\njz+O5cuXY9q0aUhMTMRbb72Fmpoa3H///Rg1ahT69u2LBQsWAABSUlKwZs0aHD9+HHfffTfeeust\nDB8+XHWsesf7Oz6942+++WZMnjwZixcvxuDBg/HKK69gxowZGDBgAGw2G4qKipCWloZHHnkEw4YN\nQ4cOHfDmm282pbUeeeQR/O1vf8Pjjz8e8N8TRS6BOwQREVkbZ+RERBbHQE5EZHGGyw9FUcSKFStw\n4sQJAMATTzyBLl26mDYwIiLyjeEZ+XfffQcAmD17NoYOHYr33nvPtEEREZHvDM/ICwoKcO211wKQ\nG/ckJiaaNigiIvJdQCs77XY7li1bhm+//RYTJkzQPf748eOBfFxEy8zMjNrri+ZrA3h9VhcL16cn\n4IedY8aMweLFi/H666+joaEh0NMREZGfDAfyrVu3ori4GADQpk0bCIKg2kuDiIiCJ6Ac+WuvvYYX\nX3wRLpcLw4cPb2qlSUREoWM4kCckJPiUFyciouBiLoSIyOIYyImILI6BnIjI4hjIiYgsjoGciMji\nGMiJiCyOgZyIyOIYyImILI6BnIjI4hjIiYgsjoGciMjiGMiJiCyOgZyIyOIYyImILI6BnIjI4hjI\niYgsjoGciMjiGMiJiCyOgZyIyOIYyIlinOSoglSyG5KjKtxDIYMMb75MRNYmNTghrlwIHC4FahxA\nahqQnQvbiIkQEtqGe3jkB87IiWKUuHIhsHM7UH0WkCT5153b5a+TpTCQE8UgyVElz8SVHC5lmsVi\nGMiJYtGpk3I6RUltNVD5U2jHQwFhICeKRR07yTlxJSntgA4ZoR0PBYSBnCgGCWnpQHau8ovZufLr\nZBkM5EQxyjZiItCnAGjXHrDZ5F/7FMhfJ0th+SFRjBIS2sI+Zrr8YLPyJ6BDBmfiFmUokLtcLhQV\nFaGyshKNjY249957cd1115k9NiIKASEtHWAAtzRDgXzbtm1ISUnBM888g7q6Ojz//PMM5EREYWIo\nkA8YMAD9+/cHAEiSBLvdbuqgiIjId4IkSZLRNzudTsyfPx+33norfvWrX5k5LiIi8pHhh52nT5/G\nggULcPvtt/scxI8fP2704yJeZmZm1F5fNF8bwOuzuli4Pj2GArnD4cDcuXPx2GOPoVevXkZOQURE\nJjEUyIuLi1FXV4eNGzdi48aNAICpU6eiTZs2pg6OiIj0GQrkhYWFKCwsNHssRERkAFd2EhFZHAM5\nkYVwNx9SwiX6RBZgZDcfyVEFnDoJ92Xx+ue/eCw6duIyfQtiICeygKbdfDya7eZjHzO9xbGtg/7J\n9ukQr8xpCvrNgzYSErndWxRgICcKISMzX19282l+rtZBXzxzGjhzGuIb8wGbvWXQlqSWG0xo3CAo\ncjGQE4VAQBsd+7Kbz8VALjmqgEMlysfu2Qm4XZf+XH1W/TMVbhAUufiwkygEAtro2J/dfLSCfvMg\nrofbvVkKAzlRkAW60bE/u/lIl10GCILRoV4SpO3eWHUTHEytEAWbH6kRNbYREy+lZmqr5UB7MTXT\nnHD+PALog3eJydu9NaWWyvZfGn9OPh+qmoSBnCjYPKkRpZy0jzNfn3fz6dgJsMcD7kaFkwhyWkfz\ngwQgpR2EYaN1x+QP8fX5wK7vL32htlpOLb0+H/ZxL5r6WbGIqRWiIDNzo2MhLR1Cbg/V90jVZwGI\nym+224EefYE4jbryi1Us0jvLfR6THslRBez7UfnFfT8yzWICzsiJQkC4vxBSXS1QeQKoq1FNjehR\nK1+UGpwQi+YBZSWA2638ZpcL6PtLYP8u/Q86XAqx/ACE8+cDXiQklZUALoWfEADA1QiprBTCL1gd\nEwgGcqIg8io7TE4FrsqHUDgWtozOxs/TrHwRAMRpT6nn4Zt7701AVJmxN1d9FtIrL0I6Vxf8RUKC\nCTn9GMfUClEQeZUd1lYDB/dC2rA6sPM0K18Ui+b5FsQB34K4R31ty88qmmeo4kTIyVNP58TFQ7gq\nz6/zkTfOyImCRLPssGQ33Ns+ha3XtbppC83zlO2XUyahsHcnxL07/Z6hC2npQLfeLR92enTrzUVH\nJuCMnChYtMoOnfXA2mUQpzwJ9+KZkBqcxs5TWy2fKxQkyf/FTBcJQ0cAV3eXU0uCIN8M+hTA9tSk\nIA44dnBGThQsWmWHHq5GYNf32r1NOnaSH46qBXNfygqD4VAJxB1fQcjJU6+iUXpGkNPN72cEpI0z\ncqIg0Sw7bO1QiWruWTf1EI4gDshlikXzIM6ZAPeyOYo/VZj1jIC0MZATBZFtxESgTwGQlKJ9YI1D\ntbdJxNdZq6RaAm1NQL5jICcKIs+KTGH8TKBtksaBAsTidyDu2uEd4E6dlGeyka51cPalNQGZgjly\nohCwZXWFO69ny80hmpMkoHQ3pMW7IbVr37IqpGMnIDFZLgeMZK37xpjQmoB8wxk5UYjYRkwErrkW\nsOvMnzypigXT5BluQiJgQkND06iVHLYKzkZbE7BDov84IycKlYZzsP3PfRB/MQBYu0z/+PIDECeP\nABITgboImo03XlD++sXg3LyNgF7XRm47Zw4GcqIgUyzBi4tX7z/SnOiOrCAOePdySU0DrsgGBt4K\n9yv/F6go924j0HCuRddGqcEJ97I53HbOJAzkREHmtXGyFR5c+kMUgaOHgKJ5Lb/eOhBr7CvKbecC\nwxw5kYla53c199CMi1Pfws1K6mq0b06tqlk0yxKVsMJFF2fkRCZQ604o3Pw79RI8lwvCQ6MgNZ4H\n/rERqDgc0jGHTOtqFq2yRCWscNHFQE5kAsVUwc7tkGqqNZfQS59+AJz+SQ5s9jj/Nki2itaB2JfW\nBc2ZvO1cNAootVJaWoqXXnrJpKEQWZNmquDQfu0l9Af3XVq+Ho1BHPAKxJplialpQLv2gM0m/9qn\nwO/NN2KR4Rn5Rx99hK1btyIhIcHM8RBZj7+pgliRmgZclacYiDXLEltVuJA+w4E8IyMDzz33HJYt\n86Eeliia+ZsqCJgAIMiNsi5LAM43GH//VfmwjXpBNRBrbiad0LZFhQvpMxzI+/fvj1OnTvn1nszM\nTKMfZwnRfH3RfG1AgNeXmYnKbr3Q8M1W8wakKQTdDo0GcbsdCdfdgPTn58DWNlH/+MxMAL2MfVaL\n00T3v089IX3Yefz48VB+XEhlZmZG7fVF87UB5lyf9PBooLYG2LszfG1lw81uBybOQWNuT5w86wDO\nhibdFAv/PvWwjpzIBEJCW9jHzwK69zF4Ahvgyww2kuX3gj23Z7hHEZMYyIlMZBs5Re4/3q69XHbo\nI+GBJ4H7HwviyELg2GHVDSa0sElW4AJKrXTs2BFz5841ayxElud5iCf+VAGp6GW574gPpA4ZEC5c\nCEX2O3hqHD71RvE0ypLapck7BbFJVsC4IIgoCKQNq30O4gCAlYsgnasz58ONVpwEWqniodIbxWv1\nq83WsgEXm2QZxtQKxYRQ/vjudy8RQN40wqyHpLf9HrDZ/X9fu8thSuNzld4oXvt3tu6i6MFt4PzG\nGTlFNbUeKMH68V1yVEHa8VV4Fwht+otf+fkmp46b0yYgORVSbTXQbFbu182tdW8W0sVATlFNrQeK\n2T++t7hhVJ81FkjNZHR2b0abgHP1kFa8DCk1DeicBeG2u4AL532/ubFJlt8YyClq+bKLu1lLwL1u\nGFauJbfHyaWQdTWaDb+aCAIAQa4jdzVe2jCj+ixQfRbSnv/Iwbl1TlwNm2T5jTlyil4h2sVd84Yh\nCHKNeLv2QI++wBMTgaQUUz43aNwuOfDe/Dvgiqv0j+/eB8LIydp18LXV+kE8NY1NsgzijJyiV6h2\ncde6YUgScNfDsN0wqGmW6d7yd+DgXt/Pr5W3TkqRH5Sa7cRR+T8tKe2AnHw58B4pg1RX48OJVfrE\nJCZDmDQPtozORkYb8zgjp6hldBd3v3luGGo+/xvEdUVNC2WEwrHynp2+Usu3t0kA7njAj4G2OmeX\nHGPvlU8A3D0M9jHT5YfGet+D5p+rNHM/VyeXbJIhDOQU1WwjJl5aaRmkHteaNwxATitcfMAKQJ51\n9uynfnxc/KWxXt1dfZPmCw3A+yuNDToxGUhMMvZeAIAEFL/TtJJT93vgkZyifhNj2aFhTK1QVNNs\nl2qipv7aZfvV969s9oBVuL8Q0u4flIN020QIw0YBmV0grX9d+4GjKBobcGYX4MQxY+/1aHaDso+Z\n3rLHuFpL346ZQNk+9fOx7NAQzsgpJghp6RByewStGsJzwxAeHgXVRTXNHrAK1Q71vHd9LYTkdnKq\nYc9/zK+ASUyG8Lv7zat199ygLn4PbNMXQXh2pvxwt9VPQkLhWPUUDMsODeOMnMhEQk4epHYqD1jb\nJkJq00YO8zoPYqU2bfxfHeoLexyEqf8L4bIEuc7b180w2iYCznPKr7WaSQtp6fINs2c/xZ+E3Nm5\nLUs1PVh2aBhn5EQm0swV19dBWjob7mVzgIRE7Qex5/1YQKPErjJHu+YXsGV01h6nPe7SA1ZBkG84\nz82VZ9VKNGbSSj8JheK5RazhjJzIZJq54mYrS3X3rQxk+zjRDXTtDlSe9D630jhbHSP9VAGUlQA5\nebBldQVg3kw6VM8tYgkDOZHJmlrZlh+A9MqLynXeh0uBhnPa+1aqBU5fpKbB9tQkAEB7x2mcOXMG\nQk5ei/4yWgFVyOoKXAzgHpo3HgOEtHQ+2DQJAzlRkAjnz6u3pm2WV1YLaGqBE4PvBpbOAZz16h+e\nnQskJEJcuRBnj5TJzbxUGob5GlA5k45cDOREwaL1QDMpxatDYGtqgVNyVEHM7KK8OjQuHujZ79JN\nYOd2NBUomtQwjDPpyMNAThQkTQ8UldIjznOXOgTqtNX1BE6pwSk/KPW05PUsrHG7gORUoGMmhMKx\nsGV0NqVhmGcnH3TsxJl3hGMgJwoir/SITaFDoI+zZK8Oi55zdO0O21OTWgZbXxqGqQTnUPdwp8Cx\n/JAoiFosknnqBfUOgTrL0zVn2JUnvb+m1ftEZ+GN104+zW42FJkYyIlCQEhLh5CcIvf4VqLXVtfP\nlrxGG4b5kpKhyMNAThQqAcySjbzXs/DGdvnPfF94E6Ie7mQu5siJQkTz4afOohrN90qSvFK09Xsu\npnUyLovHyV07fSsXDFUPdzIVZ+REIRTI8nTbiInKs/Iah2b+2p7eweeGYSHr4U6m4oycKIQCWlTT\ncE59kwkT9yBVXIjUOQvCzb8zdZ9TMg8DOVEYGFpUE0BJoV9ja3azkSqOQPq0GKgoh7Rkpk917xR6\nTK0QWUUgD0sNENLSIX2xWe6JzlLEiMZATmQRoc5fsxTROgynVkRRxMqVK1FeXo74+Hg8/fTT6NSp\nk5ljI6JWzO5AqClEqRwKnOFA/u2336KxsRFz585FSUkJ1q5di0mTJpk5NiJqJaQdCFmKaBmGUyv7\n9u1D3759AQB5eXk4ePCgaYMiIm3B3oPU8xksRbQGwzNyp9OJxMRLixBsNhvcbjfsdrvqezIzM41+\nnCVE8/VF87UBvD414owFqPrf6bhQugei4wxsaZejTW4PpD8/Bza1vjFhEO1/f3oMB/K2bdvC6XQ2\n/VmSJM0gDgDHjx83+nERLzMzM2qvL5qvDeD16RrxHOCogu1iKqcxLR0nzzqAswHsKWqiWPj702M4\ntZKfn48ffvgBAFBSUoIuXboYPRUFkeSoglSy25QKAzPPRdYSilQOGWd4Rl5QUIAff/wR06dPhyRJ\nGDVqlJnjogCZ2VNadJ5ruaEBF4UQRRRBkiQpVB8W7T/+RNL1uZfNUW6w1KfA722+4lcuQMM3W005\nVySKtL87s/H6rC2oqRWKXGYu5JAcVbiwf5cp5yKi4GAgj0Ym9ZSWGpwQV8yH6DijeS7mzonCi02z\nopFJCznElQuVd2r3SE6FuOk9oKKcuXOiMOKMPArpLeQAoDuD1kzPNMeGSkRhxxl5lFLsyXFlDuB2\nQ5wzQX8GrZWeAYDsrsCZ08qvmdgbm4j0cUYepZrv3m577o+wTV8E2O3Aru99m0HrtEwV/ud++Qah\nhHs7EoUUA3kU0HrY6FnIAcCnShbPuQCop2dy8oH0DkBisvLrYWioxAeuFMuYWokgkqNKTml07KSY\nlmj9eotFP9VngaRkIKsrbCOn+JcqqT4L8fAB4F9FLRf9XJkDXHMtbMfL5cqVZukZaelsoL5W+Xwh\nbKhk5sInIqtiII8AesFI7XW43XKqxKO+DtjzH4jTnoJt7ustApnULg2wxwGuRuVBvL0UqKu59Ofq\ns0D190CfAnR69Z2mXdjFdUXKC40AeSPhYPXGViGuXNhyPM3SRdGwWInIF0ytRICmYKSSu1Z7Hft+\nVD5hjQNi0cstviStf109iAPqs+vdP0BscOqnZ5JSIDwzA/Yx032aCZuRCuEONkQyzsjDTC8YieUH\ngLL9yq9rBebyA/K5ExIhLpkFlO7RGYhKpwZXI868OgtS4TiI32xVrk0HAGc9hAsXtD/CUQXpWDmk\nf35oTu05d7AhAsBAHn5awajGAWn1YvXqEC3n6uUd0FcubJkyMeBCyW5g1rPa40hOVX3AKf5UAWn1\nEuDUce9zBJIK4Q42RAAYyMNPKxjZ4+SZqxp7HOB2Kb+WmASp6GXgvFP5da/jk4FzdcqvuRr1byai\nKP+y4ysAgJCTJ/80sHIhsPsH7Z8eAOBQCcQdX0HIyfPpQS/QbOGTUs6eO9hQDGEgDzPNYKSnex9g\n307ApRDM1XLercXFAz37AXc+CMx7Xj/gqqmrhTh5BCC6AQBSXDxwWYLv46hxQCqaB6nZA1OtB72e\n10O6GTFRhGIgjwCKwehnnYCyfepv6todwtARkOZP0V6BqeXnXSCMngJbRmc5YCYmGT8XpKYgDkC+\nIRi5KbRKtehVpfiyGbFeWSeR1TGQRwClYAQA4uzxyoE1ORW2pyZBKisBagzkzz2qz0CqOAL3F3+X\nH6gaDuJBcLgU7tLdclpG5fXmbQCEtHSvB5usMadYwUAeQZSCkfKBglzPXbYfQAD7gpyrA4rmGX9/\nMNVWA+tWqM/qW1WlKM26WWNOsYKBPEJp1kDX1RjLqUeStonAQyOBv6xUfpCanApUq/RBB4CkFKBD\nhuqsW7i/ULfGnGkWihYM5CHid5721En1SpHQ7c4XPPFtIGR3hSQIyq93+Ln2M4KOmRDS0r23tLs4\n65bqalljTjGDgTzIvGaMyalyECocC1tGZ/U3duwkH2ukhtwKamvk2nK1YHviqPr12+3AvY9oL6aq\nPKH+ftaYU5RhIA8yrzxtbTVQWw3ppbFw9+yn+uBNSEsHOmZGbyCHpL1xhVpNOyD3mCl6GWJyivpK\n07oa4Kp85e8fa8wpyrDXShBpzhhdjbq76QiFY+VFP1YmqPwTkyT1xUzNJSbLzbhaq60GThxTf19K\nO/n716dAfr/NJv/ap4A15hR1LB4lIpzeLjuA6oM3qcEJacNqQC2HnNIOqK1BQFUroZCY5PuiICWi\nGxg+EVi5yL/zdM6SU1c6NeZE0YCBPJi0lt97qDx480rJeHhWYjrrrJF2aWgI8P1O4Puv/L4ZCLff\nfen3vpZ1ElkUUytBpLkJsofCgzfNlExiEoQ7hgIHVToiRhq3zupOX1JH//rUv89s1x5C5y7+vYfI\nwhjIg8w2YqKcp42LVz5AaVd7zfasNZD+3z/kB37RIDsXmPGKnMM28ZxMoVAsYWolyDzL78WfKiCt\nWQKcOiFXVGjsai/cX6iekpFE4F//DP2FBIPdDtwxFPjbu3KpYKAtAsKwQxFRJGAgDxFbRmdg8p9a\nPHjz2jbNs5gFMN4R0UrcbmDZHP+ba3XqLG9rV18r3xA7Z0G4/W4InbtwJk4xiYE8xDwP3vR2BhKe\n/6Mc0A+XXmyMJUXHis7WjHRIrK8DrrzKK3h7to9jl0OKNQEF8u3bt+Orr77CuHHjzBpP7NDZpkyo\nqYbtYumctONrSO++HtrxRbLaamDPfyDFt5G/R+xySDHO8BOm1atXY/369ZCicZYYCp7SRCXNKlmE\ntHQIv+ivvCgm1l2swdfbvJoo2hmekefn5+P666/HZ5995vN7MjMzjX6cJfh1fZmZqOzWCw3fbPV6\nKaFbL3To0cunY2NabTXaO6pwtvwARIWXbUfLkHFZPOzpHXRPxX+b1hbt16dHN5Bv2bIFmzdvbvG1\nkSNHYuDAgdi9e7dfH3b8+HH/RmchmZmZfl+f9PBoecFMq23KLjw82utcTceW7tHuQxJWAoKy0jQp\nRXlBUHIqqtatABzK7W7FM6dxYsd22Hr20zy9kb87K+H1WZsvNyndQD5o0CAMGjTIlAFRS75sU9b6\nWPeu74HFM0M8Ul8FIYinpgFXZAN7/qP8evkB7RF9+qG8EpYoinFBUAQQ0tIh5PZQDeKeagzJUQXh\nwoUQjy4cgwNYAAALeUlEQVTMrsqDbeQU7+ZXPfr69v6Kw9qbdBBFAZYfRjCvagx7XHSWIDZntwOi\n6FV54rWn6amTEBdM1T8fN5GgGBBQIO/Zsyd69uxp1lioFa/GWb7UXLdNBJzngjeoYEpKgTB+pvxT\nh0KaqXnzKwnQb0gGcBMJiglMrUQozQVDSuLi5aAVyUHcHgdkdVV/vWt32LK6aqaZPHxqSAaw7wrF\nBAbySOVLL3OPu4ZBeOFP6r3LI4XbJT+cjFP5QbCxEVKD0+fTNTUka9devva4ePk/gZtIUGxhjjwA\nzTdUhtl1rL70MgeAdu1hu2EQcOokpECbToWKS2VnoD0/QFy5EPYx0306jWLuHOAmEhRzGMgNUFoS\nXtmtF6SHR5u2JFxIS5e7I1Z/r33gxdSBWH1WnpVa/WHooRLFHZNaa34T9do4ggGcYgwDuQFeDyGr\nz8qrLhsafJ5NBqxVy1bh/PnoaJegU2XCvipE3pgj95Ne10IzapYlRxXEHV+pL3ZJSoHwzAzYx0xv\nCl5Su7TI2ag5kE0idKpM2FeFyFuE/J9vITpdCwOpWW4x29TKjTvrIVy40CK9IG1Y7duu9KEgSUBi\nsvFWAgmJyqf14SbKvDjFIgZyf2k9hAywZll1w+XWklMhbnoPqCiXbyrJqcC5esOfa7rUNKBzlvqy\nei01DvUHnkG8iRJZGVMrftKsXw6gZtnvuvE9/7mUXqitjpzZOCA/pLXZW6V6/CiNVEtR+dj6lyjW\nMJAb0KJ++WLvj4Rf3hhYzbJe3bgg+NdjJFA2u2916Xk9WvZA6VMgL7Hf9X2rm4sk/+Sgtgl1c57Z\ndSvBuokSWR1TKwYo1S936NErsFaaOikbYdgoCFfl+d5jJFCJScDPOgKnT8mbRatpmwzb9EUt6rjF\nKU8qH6t1nuY0Zte2ERMvPUdo1vqXC38oljGQB8CrfjnQc6ltuJyTD1u/AQD86DHiL3tcyxl0XY38\nn1oqw+NiOkjI7QEAcrWNkX04m9OYXfvT+pcoVjC1EkGUUjatl5kLaenyg0TDVNIlGmmU+ItBWpFK\nGsSw1DSfZtfNW/82b/NLFIs4I48gvs42hdvugmSkIgQArvkFILqB8oNypUtqO+BnnYCDe5WPr6tB\n6jPTULVkjhy0W7PZIaW2a7o9CDl5kOLijc/KBQFoOAf4sLiHi4OIZAzkEUgvZSNckQWpXXv/0itx\n8UDPfk1BrnV/EnHOBNX8/GX51wA5+cppH1ejXMN+sVxQSEsHuvWWH3Ya4UcZodIKW8/ioJCtsCWK\nAEytWJBm9UZOPpDbU97nUhDkWWrX7hBeWtJiJWjz1IReNYg9vQOE+wvVK05alQvanpokp4hS2vl/\ncT6WEYZihS2RVXBGblFa1RtKM26cOgnpsgTVB4N61SBCtQOSWq16q1l08xSRVFYK6Y35vte5+1pG\nyMVBRE0YyC1KL58upKVDSkj0OYcsJLSF7eGRkMpKAAEQrsprGVANrGgV0tIh/CId7u599FMtggB0\n7+N7GWEQV9gSWQ0DucVp5dN9ySFLjipIx8oh/fPDS0v+U9MgtQr4muWROrNo21OT5LGU7Aac6q0E\nbP9nqM8PKQMZD1G0YY48SunlkMWfKuBeNgfinAmQFr/Ucsm/SkdBX8ojlXh+ehAmzpZz90pS0/ye\nRRsdD1G04Yw8WunkkKXVS9RLDj08Dw0v7n4U6GIcW1ZXuLt2N20WzcVBRDLOyKOVVoOp5FTglA/t\nBDR6nviyQbISxVl0j74Qbv6d4UqTQMZDFA04I49SmjnkDj8Hyvbpn0TloaHXNmv+jKt5RUvFEUif\nFgMV5ZCWzITEBT1EhnBGHsXUcshC4Vj9HiqAV7pDanA25dXFBVMhzpkA97I5kBqcfo9NSEuH9MVm\nn3LzRKSNM/IoppVDdqvN1gGv/UA9zFxJyd1+iMzDQB4DlEoUFRcAdc6CcPvdEDp38QqipgdeLugh\nMg0DeYzyu+LD7MDLBT1EpmGOPMb5XPFh8jZr3O2HyDyGZuTnzp3DkiVL4HQ64XK58OijjyIvL8/s\nsVEECcZKSu72Q2QOQ4F806ZN6NWrF4YMGYLjx49j8eLF+NOf/mT22CjCmB14uaCHyByGAvmQIUMQ\nHy+3NHW73U2/p+gWrMBr5pZ5RLFIkCRJ0jpgy5Yt2Lx5c4uvjRw5El27doXD4cAf//hHDB8+HD16\naGwHRkREQaMbyNUcOXIEr776KoYNG4Z+/fr59J6AdpmPcJmZmVF7fdF8bQCvz+pi4fr0GEqtHDt2\nDIsWLcKzzz6L7OxsI6cgIiKTGArk69evR2NjI9asWQMASExMxKRJk8wcFxER+chQIGfQJiKKHFwQ\nRERkcQzkREQWx0BORGRxDORERBbHQE5EZHEM5EREFsdATkRkcQzkREQWx0BORGRxDORERBbHQE5E\nZHEM5EREFsdATkRkcQzkREQWx0BORGRxDORERBbHQE5EZHEM5EREFsdATkRkcQzkREQWx0BORGRx\nDORERBbHQE5EZHEM5EREFsdATkRkcQzkREQWx0BORGRxDORERBYXZ+RNDQ0NWLJkCerr6xEXF4fR\no0fj8ssvN3tsRETkA0Mz8s8//xw5OTmYOXMmfv3rX+Ojjz4ye1xEROQjQzPyIUOGQBRFAMDp06eR\nlJRk6qCIiMh3giRJktYBW7ZswebNm1t8beTIkejatStmzpyJI0eOYMaMGcjOzg7mOImISIVuINdT\nUVGBl19+GUuXLjVrTERE5AdDOfLi4mJs3boVAJCQkACbjcUvREThYmhG7nA4sHz5cjQ2NkIURTz4\n4IPo1q1bMMZHREQ6Ak6tEBFReDEnQkRkcQzkREQWx0BORGRxhhYE+Sval/SfO3cOS5YsgdPphMvl\nwqOPPoq8vLxwD8t027dvx1dffYVx48aFeyimEEURK1euRHl5OeLj4/H000+jU6dO4R6WqUpLS/Hn\nP/8ZL730UriHYiqXy4WioiJUVlaisbER9957L6677rpwD8s0oihixYoVOHHiBADgiSeeQJcuXVSP\nD8mMPNqX9G/atAm9evXCzJkzMXr0aKxatSrcQzLd6tWrsX79ekTTs/Fvv/0WjY2NmDt3Lh588EGs\nXbs23EMy1UcffYQVK1agsbEx3EMx3bZt25CSkoJZs2Zh2rRpUff/3HfffQcAmD17NoYOHYr33ntP\n8/iQzMijfUn/kCFDEB8fDwBwu91Nv48m+fn5uP766/HZZ5+Feyim2bdvH/r27QsAyMvLw8GDB8M8\nInNlZGTgueeew7Jly8I9FNMNGDAA/fv3BwBIkgS73R7mEZmroKAA1157LQCgsrISiYmJmsebHsh9\nXdJvVVrX53A4sHTpUgwfPjw8gzOB2vUNHDgQu3fvDtOogsPpdLb4H8Rms8HtdkdNUOjfvz9OnToV\n7mEERUJCAgD573DRokUYOnRomEdkPrvdjmXLluHbb7/FhAkTtA+WQuzYsWPSmDFjQv2xQVdeXi6N\nHz9e2rFjR7iHEjS7du2SXnnllXAPwzRr1qyR/v3vfzf9+amnngrjaILjp59+kqZOnRruYQRFZWWl\nNHnyZOnzzz8P91CC6uzZs9LIkSMlp9OpekxIcuTRvqT/2LFjWLRoEcaOHYt+/fqFezjko/z8fPzw\nww8AgJKSEs2HSRRZHA4H5s6di4ceegiDBg0K93BMt3XrVhQXFwMA2rRpA0EQNONmSHLkt9xyC5Yv\nX44tW7ZAFEWMHDkyFB8bMuvXr0djYyPWrFkDAEhMTMSkSZPCOyjSVVBQgB9//BHTp0+HJEkYNWpU\nuIdEPiouLkZdXR02btyIjRs3AgCmTp2KNm3ahHlk5igoKMBrr72GF198ES6XC8OHD9e8Ni7RJyKy\nuOjKcRARxSAGciIii2MgJyKyOAZyIiKLYyAnIrI4BnIiIotjICcisrj/D4OmCIgmy6fvAAAAAElF\nTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x118c557d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(x_train[:, 0], x_train[:, 1])\n",
    "plt.axis([-3, 3, -3, 3])\n",
    "plt.title(\"Simulated dataset\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model\n",
    "\n",
    "A mixture model is a model typically used for clustering.\n",
    "It assigns a mixture component to each data point, and this mixture component\n",
    "determines the distribution that the data point is generated from. A\n",
    "mixture of Gaussians uses Gaussian distributions to generate this data\n",
    "(Bishop, 2006).\n",
    "\n",
    "For a set of $N$ data points,\n",
    "the likelihood of each observation $\\mathbf{x}_n$ is\n",
    "\n",
    "\\begin{align*}\n",
    "  p(\\mathbf{x}_n \\mid \\pi, \\mu, \\sigma)\n",
    "  &=\n",
    "  \\sum_{k=1}^K \\pi_k \\, \\text{Normal}(\\mathbf{x}_n \\mid \\mu_k, \\sigma_k).\n",
    "\\end{align*}\n",
    "\n",
    "The latent variable $\\pi$ is a $K$-dimensional probability vector\n",
    "which mixes individual Gaussian distributions, each\n",
    "characterized by a mean $\\mu_k$ and standard deviation $\\sigma_k$.\n",
    "\n",
    "Define the prior on $\\pi\\in[0,1]$ such that $\\sum_{k=1}^K\\pi_k=1$ to be\n",
    "\n",
    "\\begin{align*}\n",
    "  p(\\pi)\n",
    "  &=\n",
    "  \\text{Dirichlet}(\\pi \\mid \\alpha \\mathbf{1}_{K})\n",
    "\\end{align*}\n",
    "\n",
    "for fixed $\\alpha=1$. Define the prior on each component $\\mathbf{\\mu}_k\\in\\mathbb{R}^D$ to be\n",
    "\n",
    "\\begin{align*}\n",
    "  p(\\mathbf{\\mu}_k)\n",
    "  &=\n",
    "  \\text{Normal}(\\mathbf{\\mu}_k \\mid \\mathbf{0}, \\mathbf{I}).\n",
    "\\end{align*}\n",
    "\n",
    "Define the prior on each component $\\mathbf{\\sigma}_k^2\\in\\mathbb{R}^D$ to be\n",
    "\n",
    "\\begin{align*}\n",
    "  p(\\mathbf{\\sigma}_k^2)\n",
    "  &=\n",
    "  \\text{InverseGamma}(\\mathbf{\\sigma}_k^2 \\mid a, b).\n",
    "\\end{align*}\n",
    "\n",
    "We build two versions of the model in Edward: one jointly with the\n",
    "mixture assignments $c_n\\in\\{0,\\ldots,K-1\\}$ as latent variables,\n",
    "and another with them summed out.\n",
    "\n",
    "The joint version includes an explicit latent variable for the mixture\n",
    "assignments. We implement this with the `ParamMixture` random\n",
    "variable; it takes as input the mixing probabilities, the components'\n",
    "parameters, and the distribution of the components. It is the\n",
    "distribution of the mixture conditional on mixture assignments. (Note\n",
    "we can also write this separately by first building a `Categorical`\n",
    "random variable for `z` and then building `x`; `ParamMixture` avoids\n",
    "requiring `tf.gather` which is slightly more efficient.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "pi = Dirichlet(tf.ones(K))\n",
    "mu = Normal(tf.zeros(D), tf.ones(D), sample_shape=K)\n",
    "sigmasq = InverseGamma(tf.ones(D), tf.ones(D), sample_shape=K)\n",
    "x = ParamMixture(pi, {'loc': mu, 'scale_diag': tf.sqrt(sigmasq)},\n",
    "                 MultivariateNormalDiag,\n",
    "                 sample_shape=N)\n",
    "z = x.cat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The collapsed version marginalizes out the mixture assignments. We\n",
    "implement this with the `Mixture` random variable; it takes as\n",
    "input a Categorical distribution and a list of individual distribution\n",
    "components. It is the distribution of the mixture summing out the\n",
    "mixture assignments. Gibbs sampling does not work with `Mixture` random\n",
    "variables, please try an alternative."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'\\npi = Dirichlet(tf.ones(K))\\nmu = Normal(tf.zeros(D), tf.ones(D), sample_shape=K)\\nsigmasq = InverseGamma(tf.ones(D), tf.ones(D), sample_shape=K)\\ncat = Categorical(probs=pi, sample_shape=N)\\ncomponents = [\\n    MultivariateNormalDiag(mu[k], tf.sqrt(sigmasq[k]), sample_shape=N)\\n    for k in range(K)]\\nx = Mixture(cat=cat, components=components)\\n'"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"\n",
    "pi = Dirichlet(tf.ones(K))\n",
    "mu = Normal(tf.zeros(D), tf.ones(D), sample_shape=K)\n",
    "sigmasq = InverseGamma(tf.ones(D), tf.ones(D), sample_shape=K)\n",
    "cat = Categorical(probs=pi, sample_shape=N)\n",
    "components = [\n",
    "    MultivariateNormalDiag(mu[k], tf.sqrt(sigmasq[k]), sample_shape=N)\n",
    "    for k in range(K)]\n",
    "x = Mixture(cat=cat, components=components)\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will use the joint version in this analysis."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Inference\n",
    "\n",
    "Each distribution in the model is written with conjugate priors, so we\n",
    "can use Gibbs sampling. It performs Markov chain Monte Carlo by\n",
    "iterating over draws from the complete conditionals of each\n",
    "distribution, i.e., each distribution conditional on a previously\n",
    "drawn value. First we set up Empirical random variables which will\n",
    "approximate the posteriors using the collection of samples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "T = 500  # number of MCMC samples\n",
    "qpi = Empirical(tf.get_variable(\n",
    "    \"qpi/params\", [T, K],\n",
    "    initializer=tf.constant_initializer(1.0 / K)))\n",
    "qmu = Empirical(tf.get_variable(\n",
    "    \"qmu/params\", [T, K, D],\n",
    "    initializer=tf.zeros_initializer()))\n",
    "qsigmasq = Empirical(tf.get_variable(\n",
    "    \"qsigmasq/params\", [T, K, D],\n",
    "    initializer=tf.ones_initializer()))\n",
    "qz = Empirical(tf.get_variable(\n",
    "    \"qz/params\", [T, N],\n",
    "    initializer=tf.zeros_initializer(),\n",
    "    dtype=tf.int32))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run Gibbs sampling. We write the training loop explicitly, so that we can track\n",
    "the cluster means as the sampler progresses."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " 50/500 [ 10%] ███                            ETA: 6s | Acceptance Rate: 1.000\n",
      "Inferred cluster means:\n",
      "[[-0.99587333 -0.95255095]\n",
      " [ 0.95759684  1.00283527]]\n",
      "100/500 [ 20%] ██████                         ETA: 4s | Acceptance Rate: 1.000\n",
      "Inferred cluster means:\n",
      "[[-1.00205481 -0.96149278]\n",
      " [ 0.97067398  1.01528645]]\n",
      "150/500 [ 30%] █████████                      ETA: 3s | Acceptance Rate: 1.000\n",
      "Inferred cluster means:\n",
      "[[-1.0040288  -0.96395862]\n",
      " [ 0.97497743  1.01705074]]\n",
      "200/500 [ 40%] ████████████                   ETA: 2s | Acceptance Rate: 1.000\n",
      "Inferred cluster means:\n",
      "[[-1.00589085 -0.96495295]\n",
      " [ 0.97453934  1.01985526]]\n",
      "250/500 [ 50%] ███████████████                ETA: 2s | Acceptance Rate: 1.004\n",
      "Inferred cluster means:\n",
      "[[-1.00728667 -0.96597838]\n",
      " [ 0.97563273  1.02098477]]\n",
      "300/500 [ 60%] ██████████████████             ETA: 1s | Acceptance Rate: 1.000\n",
      "Inferred cluster means:\n",
      "[[-1.00836658 -0.9667052 ]\n",
      " [ 0.97756976  1.02122355]]\n",
      "350/500 [ 70%] █████████████████████          ETA: 1s | Acceptance Rate: 1.000\n",
      "Inferred cluster means:\n",
      "[[-1.0085218  -0.9666754 ]\n",
      " [ 0.97716254  1.02096808]]\n",
      "400/500 [ 80%] ████████████████████████       ETA: 0s | Acceptance Rate: 1.000\n",
      "Inferred cluster means:\n",
      "[[-1.00917315 -0.9672277 ]\n",
      " [ 0.97824979  1.02170002]]\n",
      "450/500 [ 90%] ███████████████████████████    ETA: 0s | Acceptance Rate: 1.000\n",
      "Inferred cluster means:\n",
      "[[-1.00925195 -0.96741939]\n",
      " [ 0.97901499  1.02244294]]\n",
      "500/500 [100%] ██████████████████████████████ Elapsed: 4s | Acceptance Rate: 1.000\n",
      "\n",
      "Inferred cluster means:\n",
      "[[-1.00926054 -0.96760303]\n",
      " [ 0.9795807   1.02280807]]\n"
     ]
    }
   ],
   "source": [
    "inference = ed.Gibbs({pi: qpi, mu: qmu, sigmasq: qsigmasq, z: qz},\n",
    "                     data={x: x_train})\n",
    "inference.initialize()\n",
    "\n",
    "sess = ed.get_session()\n",
    "tf.global_variables_initializer().run()\n",
    "\n",
    "t_ph = tf.placeholder(tf.int32, [])\n",
    "running_cluster_means = tf.reduce_mean(qmu.params[:t_ph], 0)\n",
    "\n",
    "for _ in range(inference.n_iter):\n",
    "  info_dict = inference.update()\n",
    "  inference.print_progress(info_dict)\n",
    "  t = info_dict['t']\n",
    "  if t % inference.n_print == 0:\n",
    "    print(\"\\nInferred cluster means:\")\n",
    "    print(sess.run(running_cluster_means, {t_ph: t - 1}))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Criticism\n",
    "\n",
    "We visualize the predicted memberships of each data point. We pick\n",
    "the cluster assignment which produces the highest posterior predictive\n",
    "density for each data point.\n",
    "\n",
    "To do this, we first draw a sample from the posterior and calculate a\n",
    "a $N\\times K$ matrix of log-likelihoods, one for each data point\n",
    "$\\mathbf{x}_n$ and cluster assignment $k$.\n",
    "We perform this averaged over 100 posterior samples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Calculate likelihood for each data point and cluster assignment,\n",
    "# averaged over many posterior samples. ``x_post`` has shape (N, 100, K, D).\n",
    "mu_sample = qmu.sample(100)\n",
    "sigmasq_sample = qsigmasq.sample(100)\n",
    "x_post = Normal(loc=tf.ones([N, 1, 1, 1]) * mu_sample,\n",
    "                scale=tf.ones([N, 1, 1, 1]) * tf.sqrt(sigmasq_sample))\n",
    "x_broadcasted = tf.tile(tf.reshape(x_train, [N, 1, 1, D]), [1, 100, K, 1])\n",
    "\n",
    "# Sum over latent dimension, then average over posterior samples.\n",
    "# ``log_liks`` ends up with shape (N, K).\n",
    "log_liks = x_post.log_prob(x_broadcasted)\n",
    "log_liks = tf.reduce_sum(log_liks, 3)\n",
    "log_liks = tf.reduce_mean(log_liks, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We then take the $\\arg\\max$ along the columns (cluster assignments)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Choose the cluster with the highest likelihood for each data point.\n",
    "clusters = tf.argmax(log_liks, 1).eval()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the data points, colored by their predicted membership."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAEICAYAAABCnX+uAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xd0FNXbB/DvzPZseqMjEEyASJMeQUMRUIoiIohUQQWR\n+iIdxUJXpBfpRUGRqggIIqH8kKJ0pARISKjpbUt2Zu77x0KSJbvJZrMh2eT5nMM57NQ7u5tn7j5z\nC8cYYyCEEOKy+OIuACGEkMKhQE4IIS6OAjkhhLg4CuSEEOLiKJATQoiLo0BOCCEujgJ5CdGmTRuE\nhIRk/atTpw5efvllzJo1CxkZGU49V7NmzbB9+3YAwIQJEzBixAi79jtw4ADu37/v8Hlnz56Nvn37\nFni/kydPIiQkxCnvQ2xsLA4ePFjo4zwLbdq0waZNm4q7GA7JyMjA1q1bi7sYZYa8uAtAso0ZMwZv\nvfUWAECSJNy8eROffvop0tLSMH369CI55+TJk2FPV4K7d+/ik08+wa+//ooKFSoUSVmehYkTJ+KF\nF15Au3btirso+frll1+g0WiKuxgOWbt2LQ4dOoQePXoUd1HKBKqRlyBarRYBAQEICAhAuXLlEBYW\nhn79+mH//v1Fdk4PDw94enrmux31G3v2fH19XTaQ0/fl2aJAXsLJ5XIolUoAwKJFizB48GAMGjQI\njRo1wo4dOwAAq1evRuvWrdGwYUO8++67OHfuXNb+oihi7ty5aN68OZo1a4b169dbHP/p1Mq+ffvQ\ntWtX1KtXD507d8ahQ4cAAG3btgUAdOnSBYsWLQIAnD9/Hr169ULdunXRvn17rFy5EpIkZR3ryJEj\n6NKlC+rVq4chQ4YgLS0tz2s9efIkevbsifr166Ndu3b45Zdfcm0TGxuLkJAQXL9+PWvZ9u3b0axZ\ns6zXP/30E1599VW88MIL6NixI3bu3Jl1radOncKaNWvQpk0bAEB6ejqmTp2Kpk2bolmzZhgxYgQe\nPnyYdayQkBDMnz8fYWFh6Nq1K0RRzFWmo0eP4p133kG9evVQv3599OvXD7dv385av2TJErzyyiuo\nW7cu3nzzTURERORbVsAytcIYw4IFCxAWFoZGjRphxowZ6Nu3r0WKbNq0aZgwYQIaNmyIFi1aYMmS\nJVnHmjBhAmbOnImJEyeiQYMGaNOmDSIiIrBz506Eh4ejcePGmDp1qkUAzut7ldf5tm/fjsWLF+Py\n5csICQlBbGwsIiMj0bdvXzRs2BDNmzfH5MmTodPprH8RSMExUiK0bt2abdy4Meu1KIrs/PnzLDw8\nnE2dOpUxxtjChQtZcHAwW758OYuMjGQJCQls8+bNLDw8nB0+fJjdvn2bLVu2jNWrV4/FxMRk7RMW\nFsYiIiLYf//9x/r27cuCg4PZtm3bGGOMjR8/ng0fPpwxxtiJEydYrVq12OrVq1lUVBRbs2YNCw0N\nZZGRkez8+fMsODiYnTx5kqWnp7P4+HjWqFEjtmzZMhYVFcUOHz7MwsPD2YoVKxhjjN26dYuFhoay\nRYsWsVu3brEVK1awkJAQ1qdPH6vXf/PmTRYaGspmz57Nbt26xXbt2sVCQ0PZ//73P/b333+z4OBg\nlp6ezmJiYlhwcDC7du1a1r7btm1jTZs2ZYwxdvnyZVarVi22d+9eFhsby3744QcWEhLCbt++zVJT\nU1nPnj3Z559/zhISEhhjjI0ePZr16dOHXbhwgV27do2NGDGCde7cmZlMJsYYY8HBwaxdu3bsxo0b\n7MqVK7nKHRMTw0JDQ9nq1avZnTt32NmzZ1nXrl3ZkCFDGGOMHThwgDVs2JAdP36cxcTEsHnz5rEG\nDRqwtLS0PMv69HdixYoVrFmzZuzPP/9k165dY4MHD2YhISEWn2NoaChbsGABi46OZitWrGDBwcHs\n8uXLFutXrlzJoqOj2YgRI1ijRo1Y37592dWrV9lvv/3GateuzQ4ePMgYY/l+r/I6n16vZ7NmzWJd\nu3Zljx49YoIgsK5du7Jx48ax6Ohodu7cOda6dWv23Xff2fOnQexAgbyEaN26NQsNDWUNGjRgDRo0\nYHXq1GGhoaFs+PDhLCUlhTFmDsr169dnkiRl7RceHs52795tcayBAweyWbNmMUmSWFhYGNu0aVPW\nuvv377PatWtbDeTDhw9nw4YNszjW4sWL2YULF3IF0AULFrCBAwdabLt7927WrFkzxhhjc+bMYW+9\n9ZbF+vfff99mIJ81axZ78803LZZt3LiRHTt2rECB/I8//mC1a9dm586dy1p/7NixrPewT58+bNas\nWYwxxu7cucOCg4PZgwcPsrY1Go2sQYMG7K+//mKMmQP5k5uTNbdv32YbNmywWLZ69WrWtm1bxhhj\na9euZU2bNs0KzkajkR09epTp9fp8y5ozkLds2ZKtX78+a7vExERWv359i8/xtddesyhH06ZN2ZYt\nW7LWd+7cOWvd4cOHWXBwsMXN6fXXX2fLly9njOX9vbLnfAsXLmTdunXLWvfiiy+yGTNmZN0gr169\nym7evGntLSUOoIedJchHH32Erl27AgAUCgX8/f2z0ipPVKpUCRzHATC3DLh37x6mTJmCzz77LGub\nzMxMKJVKJCUlIT4+HnXq1MlaV758eQQGBlo9/82bN7PO/8SwYcMAmFMaOUVGRuLkyZNo2LBh1jJJ\nkmAwGJCUlIQbN24gNDTUYp969erhzJkzNs9dt25di2V9+vQBYE652KtVq1aoV68e3nnnHdSoUQPh\n4eHo1q2b1ecAkZGRAICOHTtaLNfr9bh16xbCw8MBAFWqVLF5vmrVqkGj0WDlypW4ceMGbt++jf/+\n+y/rPe7SpQu2bNmCDh06oE6dOggPD0f37t2hVqvtLmtiYiIePXpk8f74+PigWrVqFttVrVrV4rVW\nq4UgCFmvc16HWq3OtUylUiEzMzPf75W958tp+PDhmD17NrZv346WLVvi1VdfxWuvvWZ1W1JwFMhL\nEB8fHzz33HN5bqNSqbL+/yRfO2vWLItgDWT/oQK5HzwpFAqrx1YoFHY/pBIEAe3bt8eoUaNyrfPw\n8ADHcXaftyDnfnITyyln3lqtVmPz5s04e/YsIiIicOjQIWzatAnLly/HSy+9lGs/hUJhkZd+wsvL\ny+KYtly7dg3vvvsuwsLC0KRJE7z99ts4f/48Nm/eDADw8/PDnj17cPLkSURERODXX3/Fxo0bsWnT\nJtSqVcuusj5533I+f7Dm6Zs+YPnZy+W5/9zzej/z+17ld76cBgwYgI4dO+LPP//E0aNHMW7cOBw9\nehQzZ860cTWkIOhhpwvz9PREQEAAHj58iOeeey7r3/r163H06FH4+PggICAAFy5cyNonMTHRZlvw\natWq4cqVKxbLBg0ahPXr1+f6gw8KCsLt27ctznvjxg0sWrQIPM8jODjY4rwAch376XNfvnzZYtmU\nKVMwa9Ysi2VPglrONuUxMTFZ/z979iwWLVqEF198EaNHj8avv/6K0NBQ/PHHH7nOWaNGDZhMJuh0\nuqxr8Pf3x8yZMxEVFWWzrDn9/PPPqF27NhYvXoz+/fujadOmuHv3blZAi4iIwMaNGxEWFoaJEydi\n37598PDwwJEjR+wuq4eHB8qXL2/x/qSlpSE6OtquMhZUft8re+T8vhiNRkyfPh0mkwnvvfceli9f\njs8++wx79uwpkvKXRRTIXdzgwYOxdOlS/P7777hz5w4WL16Mn376CTVq1ADHcRgwYACWLVuGgwcP\n4saNG5g4caLVlhcA0L9/fxw8eBCbNm3CnTt3sG7dOpw+fRotW7aEm5sbAOC///5DWloa3nvvPURF\nReHrr7/GrVu3cOzYMXz++efw8PAAz/Po2bMnoqOjMWfOHNy+fRubNm3C4cOHbV7Hu+++ixs3bmD+\n/PmIiorC7t27sWvXLrzyyisW2/n7+6NChQpYu3Yt7ty5gwMHDmS13AAAjUaDFStWYMOGDYiNjcWx\nY8dw8+ZN1KtXD4D55390dDQePnyIGjVqoE2bNhg3bhzOnDmT1W7/0qVLCAoKsuv9L1euHG7duoUz\nZ84gJiYG69atw9atW5GZmQnAXEP99ttv8dtvv+Hu3bv4448/EBcXh7p16+Zb1pwGDBiA5cuX4/Dh\nw4iMjMTEiROh0+ms1qidIa/vlT3c3NwQHx+PmJgYyGQy/Pvvv/jyyy9x/fp13Lx5EwcOHLB6ncQx\nlFpxcf369YPBYMDcuXMRHx+P6tWrY+HChXjxxRcBmGvUmZmZ+Pzzz2EwGNC3b1+bNbmGDRti9uzZ\nWLJkCWbPno2goCAsWbIkK6j16NEDU6ZMQa9evTB58mSsWrUK33zzDd544w14eXmha9euGD16NACg\ncuXKWLVqFWbMmIGNGzeiQYMG6NWrF65du2b13JUrV8by5cvxzTffYPXq1ahYsSKmT5+OFi1aWOTI\neZ7HzJkz8fXXX+P1119H/fr1MXr0aMyZMwcAUKtWLcydOxdLly7F3Llz4evri4EDB6J79+4AzDeM\n8ePHo2vXrjhx4gRmz56NmTNnYtiwYcjMzETDhg2xbt06eHh42PX+9+3bF9euXcOQIUPAcRzq1KmD\nadOmYerUqXjw4AHCw8Mxfvx4LFiwAPfv30eFChUwdepUtGjRAgDyLGtO/fv3R1xcHMaPHw9BENCr\nVy9UqlQpz3RVYeT3vcpPhw4dsHXrVrz++uv44YcfMH/+fHz99dfo3bs3RFFEWFgYvvrqqyIpe1nE\nMXuTooSQYhMREYHQ0FD4+/sDMD+jaNasGZYvX44mTZoUc+lIcaPUCiEuYOvWrfj0009x/fr1rJSW\nl5cX6tevX9xFIyUA1cgJcQGPHj3CV199hb///huCIODFF1/ElClTUL169eIuGikBHA7kkiRh+fLl\nWS0gPvjgg1ztSgkhhBQ9h1MrTzp2fPXVV+jVqxe2bNnitEIRQgixn8OtVpo2bYpGjRoBAOLi4rKa\npxFCCHm2CtX8UCaTYfHixTh9+jTGjBmT7/b37t0rzOlKtIoVK5ba6yvN1wbQ9bm6snB9+Sl0q5VP\nPvkECxYswIoVK2AwGAp7OEIIIQXkcCA/cuRI1njYSqUSHMeB56k1IyGEPGuFypEvXboUn3/+OQRB\nwIABA6wOokMIIaRoORzI1Wq1XXlxQgghRYtyIYQQ4uIokBNCiIujQE4IIS6OAjkhhLg4CuSEEOLi\nKJATQoiLo0BOCCEujgI5IYS4OArkhBDi4iiQE0KIi6NATgghLo4COSGEuDgK5IQQ4uIokBNCiIuj\nQE4IIS6OAjkhhLg4CuSEEOLiKJATQoiLo0BOSFlnMJj/EZfl8JydhBDXJr94EV5ffglZdDTAcRCq\nVUPqtGkQatcu7qKRAqJATkgZxN+9C98PP4T8zp2sZfLYWMgGDULC9u2QypcvxtKRgqLUCiFlkMc3\n31gE8ScU0dHw+O67YigRKQwK5ISUQfKYGNvroqKeXUGIU1AgJ6QMYhqNzXVSHutIyUSBnJAyKKNv\nX4ju7rmWS56e0A0cWAwlIoVBgZyQMsjYvj10ffpALFcua5lQvjwy+vWD8ZVXirFkxBHUaoWQMipt\n6lRkfPAB3H76CeA46Hr1ghQYWNzFIg6gQE5IGSaVL4/0kSOLuxikkBwK5IIgYNmyZYiLi4PJZEL3\n7t3RuHFjZ5eNEEKIHRwK5EePHoWHhweGDx+O9PR0fPrppxTICSGkmDgUyFu0aIHmzZsDABhjkMlk\nTi0UIYQQ+3GMMeboznq9HnPmzEHbtm3RsmVLZ5aLEEKInRx+2BkfH49vvvkG7du3tzuI37t3z9HT\nlXgVK1YstddXmq8NoOtzdWXh+vLjUCBPTk7G9OnT8f7776Nu3bqOHIIQQoiTOBTId+zYgfT0dGzb\ntg3btm0DAEyaNAlKpdKphSOEEJI/hwL5wIEDMZC68RJCSIlAXfQJIcTFUc9OQlwAl5gI90WLoPjv\nPzC1Gro+fWBs1865J2EM4DjnHpM8ExTICSnh+Lt34de7NxSRkVnLlP/7H/Q9eyL1q6+s7yQIUO/a\nBcXFi8CrrwLNmwNW+ntwOh08J0+G8t9/wRmNEMuVQ/qQITC+9lpRXQ4pAhTICSnhvKZMsQjiACDL\nyIBm505kDBwIsUYNy3WRkfAZMgSKyEhwJhOwYQP8n38eiStXQnX8ODR79wKMwdC+PTQ7d0L1999Z\n+8pjYiCfMAHJHAdjx47P5PpI4VEgJ+RZYAxumzZBs3MnOL0eYmAg0kaNgtCgQb67yp8K4k/IEhOh\nXbvWslbOGHxGjoTyv/+ylxmNUF66hICOHcHrdObgDkAVEWFOpzx93Ph4uC9fToHchVAgJ+QZ8Jw0\nCW4//wzeYMhaprh4EcnffovM8PA89+Xy6nwtihYv5f/9B/mlS1Y35VNSkDMDzj21b06yhw8pZ+5C\nqNUKIUWMv3cPmv37LYI4AMgfPIDnvHn57i9Ur251uejjA13//pbHvHYNvCBY3b4gIZmp1UUTxBnL\ndfMhhUeBnJAiptm501zDtYK/exdcSkqe+6dMm5YrmEtqNQwdOkAICbFYLrtzp3CFBcAAZDZqVOjj\n5MTpdPAaNQoBr7yCwLAw+HftCvWvvzr1HGUZpVYIKWJMq7W9UiYDFIo89xeDghC/bRs8vvkG8ps3\nwVQq6Lp3h6F791zbSv7+YLBe+7a1XFIowD/Om0sqFUyNGyPFVmsYRzAG3379oDpxIntZbCy8oqLA\nOA7Gzp2dd64yigI5IUVM360b3Jcvh9xKbVmoUQPMzS3fY0jlyiFl7tx8t2N5DZPBcVYfbj55+AkA\nkMkgVKsGqNX5nsteqogIKE6dyrVclpAA95UrKZA7AaVWCClizNMT6R99BDEgwGK5qWZNpMyc6ZRz\ncElJ8B46FN6ff24zFy76+UH09c29b47/8zodNDt2QPXXX04pFwC4LV0K3kZeXBYVZfXmQgqGauSE\nFDHFP/9AFh+PtMGDobx2DVx6OoRatZA+ZAiYl5fdx+Hj4qBduxZcYiL03brB1LQpwHHgHjxAYIcO\nkMXH57m/LDERzI5JYHidDm7r10MWHQ3ZgwfQv/46hPr17S7n0xR379o+l15PLWOcgAI5IUWE0+ng\n8/77UJ47Bz4tDUwmg1CjBpLnzoWpSZMCHUu7ejW0y5ZBfv8+APMDVFOTJkiZOBH+PXtClpiY/0Ek\nCZwk2XU+1bFjUB88CA6AduNGZDZqhMTVqwEHRjiVfH2BqCjr6zw8Cnw8khulVggpIl7jxkF99Cj4\ntDQA5nbbihs34P3pp4DBYHdKQRYbC/clS7KCOADI0tKgOnwYfv372xfEYU6h2FP3ZQB4gyFrWz4l\nBapDh+AzcCAU//xT4FSIUK2azXX6Tp0KdCxiHQVyQoqCyQTl2bNWV8lv3ED52rUR2KgRvIcOzbf5\noXbZMqvNFzlJAv/okVOKa3FcG8s0hw/Dr2dP+L/2GuSXL9t9vLQxYyBUqpRrualmTaSPHu14QUkW\nCuSEFAFOpwOn11tfB4DPzIT84UO47d4Nv169ABudeACAT0qyfZ489isKvF4P5cWL8PnkE/BRUVD8\n8w+45GTrGzMG93nz4DtoELj0dEgqFSSNBkKVKjC0aYPEjRvBfHyeaflLKwrkhBQB5ukJ0d/frm0V\nFy9CvXOnzfWGl18Gs/FAsLgeE8qvX0dAhw7w794dAR06wGvUqFw3I/f58+G+dCkU165BlpIC3mgE\nr9dDeO45JG7cCLFq1WIqfelDgZyQosBx0PfsCdGOh3kcY1Dn0dyPE8US10SPAyBLTwdnMkEeGwu3\n7dvhNXVq9gaiCM2vv5pbpTxFefEi5BcuPLvClgEUyAkpIhmDBiF97FhkvvACpLx6dwJgGo3NdU9a\nj5RknChCdfSo+SEuAD4xEbyNh7B8SgpUx48/y+KVehTICSlCGYMHI37fPiSuXAkpj6746j/+gG/P\nnnD7/vvcOWc72n6XBFxKCmQJCQDMzQptDU0gKZUQgoKeZdFKPQrkhBQ1jkPmK68gMywMthIksoQE\nqI8dg/cXXyCgQwe4rVmTtU73zjtgfMn4U2UcZ7MszNMzu+eoWo3Mhg2tbifUrAlj27ZFVcQyiToE\nEfIsSBJSJ0yA1scHmr17wRuNNjeVx8bC6/PPoV23DplNmpi7sdvZkaeoSR4eEJ57DqqLFy2WM56H\nMSwMyJEiSpk7F3xSEhTnzkGWnAxJpYIQFITkhQut/spQnjgB94ULIXv0CEyjgSE83Nw80UV+kRQn\nCuSEFDH17t1wX7TIPGgWx9nVJZ2TJChu3oTi5k2boxYWB1lqKnD/PgRfX3CiCD4lBWLFishs1gwp\nM2ZAfukSPJYsAZecDLFiRSTPmQM+JQXK48ch1qgBY+vWwOMavfzCBfOQA4IAU3Aw3NesgSxHu3jF\nxYuQR0Yiefny4rpcl0GBnJAiJL9wAV6ffQZZXJzDxygpQfyJJ2O6iN7eMLz+OgwtW0Kzbx8CwsMh\nu3/f4teG6tgxJC9cCN0HH2QfgDF4TZhgbtXyuDMU4/lcwwdwggDVsWOQX7uWa9x1YqlkJN4IKS1E\nEcgRyDwWLSpUEC/JZMnJUO/fD6/PPoP6yBEooqJypYzksbHwfGpsc9WhQ9Bs354VxAHYHANGlpQE\nzY4dzi98KUM1ckKcgE9IgNe4cVBcvQqYTBArVkTaiBHgbcwMBACSVgtJq4UsPt7uwaxKGk4U8/3F\nILtzB/zdu5Aed9N3++EH8Dqd3eeQqPdnviiQE1JYRiN8e/eGMsekx/K7dyEfO9Zmj0zg8Ww8deuC\nO38efGKiywbzfDFmMdEzl5lp965ClSrQ9epVFKUqVSi1QkghuW3ZYq6JP0X28CF4W+OQwDw+uObP\nP126Rm4PsXJliFWqZL3OfPFFm80wc974hPLlkfbJJwUas72sKlQgv3HjBqZNm+akohDimlQnTtgc\nvIp73NPR6joby/PqjF+yOurnTyhXDmmjR1u01Mn46COYXngh17am6tWROn489F26IH3AAMTv3g19\nnz7Psrguy+HUyq5du3DkyBGonTi3HyGuyNr0aU840uIkr3042J5E2VmccXxJrYaxbVukjRoFoU4d\ny+NrtUjYvBle06ZBceUKIEkQgoKQ+tlnEKtUQUYhz10WORzIy5Urh7Fjx2Lx4sXOLA8hLidj2DCo\n9+2DPI8Hm84iqVTg8uhM5AyFvVlIWi3Sxo9HxqBBNrdhvr7mjkHEKTjGHB9W7dGjR1iwYAGmT5/u\nzDIR4nq+/x4YOTJr0KgyieOAOnWAvn2B8eOLuzRlyjNttXLv3r1nebpnqmLFiqX2+krztQFOur7O\nneF+8SI8Fi8ucE32Se23JPXgLCjR0xNJy5Yhs2VLQC4HnuH3pSx8P/NDrVYIcZKMkSNhqlu3wPsx\ntRrpLtrEjgEQ/P1h6NgRpsaNzUGcPHMUyAlxEubmhoSffkJG797IrFsXUh5jjFvsJ5cjbfJkMJWq\niEtYBHge8vh4uP38MwLat4d25cr898mRzeXv34dq3z7IL10qcZNnuJJC5cgLqrT//Cmt11earw0o\nmutznzsX7osWgc/REcYWJpNBqF4d8lu3iq09ubPSOqKvLxI2bIDw9BC2jEG7eDE0v/0GPjkZzM0N\nkCTwqamQPXoEyd0dpuefR/LixRCrVSvQOcvC9zM/VCMnxMnkZ89Cu26dXUEcMHdzV0RGOiWIM8Du\nXwJPMxUwgFojS0yEh5XRCj2//hoeCxZAeekS5LGxUFy/DnlkZNZoh3x6OlRnz8Lno49KzJC9roQC\nOSn13Natg1/XrggID4ffG29A8+OPRXo+95UrIcujR2eR4nkYWrYscMchDsizF2qBjpWaavlap4N6\n375c83da+wUgj4yE6vBhp5SjLKEnE6RU85g5E9q1a8FnZHczkV+7BllcHNJHjnTqubjkZKj/+APy\nW7ecetwClUGS4HbggEP75hyNsFBlMBqh2boVhtdeA3N3h/zaNfB2pj54gwGKS5dgbNPGKWUpKyiQ\nk1KLS0+H+rffLII4AMjS0qDZtg3pH30EOKlnsvvcuXDbtg3ymBgwuGZTQi6Px2VP1uR3TUwmg+rk\nSShPnoQ4dy7ESpXAFAq7JtMAAFGrhbFxY/sKTLJQICelluL8ecijo62uk8XEQHH1KkwNGhT6PKrf\nf4d21SrI0tMB5DGGCs+DcZzduXN7FfSmkd/2TCaDWL48wPPg0tPB6XTgjUab6RpJoQBzcwOflpY1\nyiGHxyNA3r1bgJIBYu3aMLVoUaB9COXISSnGPD3BbD3402ggubs75TzajRuzgvjTJI0GQpUqMLz8\nMhI2bULqlCkFzl8/PRQuy/FP0mgKXPPnrBzTYr0oQnb3LvgHD8AnJ2dNFvFkjyfllzQaZNavj/hf\nf4WpTh27HtZau3YGc2sXQ3g4Eteutbv2TrJRICellumFFyDUqGF1nVCjBsSaNZ1yHt5GEAcATq8H\nRBGSry8yW7SAvlcvmzcQW0Hu6ZQHB0Dy8kLyd99BrFDBoTKbQkIg5ZFW4gDwJpPVdAsHwBQcjKTv\nv0f8nj0Q6tYFZzLZdV7GcZCsBGqpXDkkLVsGKY8ByIhtFMhJ6cVxSJ4+HUL16llBknEcTEFBSJ49\n22mnEQMCbBcBgPzePWh27YLXuHFgnp7IDA+3uW2usGmjp6QsJQVeU6eCfzx/ZkGJlSqZ23I7ShAg\nVK+eVXuWypWzbz+5HFAoLBZxABT//QePmTMdL08ZR4GclGpC48aI27sXaWPHQvfGG0ibMAHxe/dC\nCA112jnSRo2CkE8g4xiD6tQpcOnpSP7uO5iCgqzWwDkAopcXjM2bI23oUIje3jaPyaenm2e1LyAm\nk8HQoUOhelIqbt2C/5tvwnvIEEAUkTZ6NITy5fPfUaUCb2OGIOX58w6Xp6yjQE5KPebhgfTRo5G8\ndCnSP/kETKt16vGFevWQMmsWMuvXh5jHsfmkJPAPHoC5uUEICbGZ25YqV0bC1q2QxcZClpTk1LIy\nAPqOHcG02kIfWxYfD83evfCYORNC7dpI/vZbZDZsCNHfH5K7O5hMZrG96O0NoWrVPApHXfQdRa1W\nCHECY/sV4imlAAAgAElEQVT2ML76KpQREfAdMgR8WlqubSRvb0iPa615TSgsqdVw27gRmv37Lea6\ndAbJywvJy5dDfvUqJC8vu9uO22rpwgkCVBERSJsyBZnh4YgPDzdPOC2XQ3HmjLkNf2oqJC8vZAwe\nDM5ohPcnn2Q9QM3J5MRfSWUN1cgJcRaOQ2Z4uNWAxGB+8OkxYwb4hASkDxsGMTAw13aSUglDhw5Q\n79lToEmKc53PygNFBiBj4ECA5yHUqQNTrVpWt3m6XswUChibN7d5Ll6ns3gtlSsHyc8Pxg4dkLhl\nC+J//x2JmzfD2LYtDB07wtiyJdhTuf/M2rWROmmSvZdHnkKBnBAnS1y5EoZWrSD6+1t0DpLFxcF9\n/Xr4desG8DxSx4yBkGNSYtHPD/ouXZDx8cd2twKxhanVFmkeSamEsVUri96sSStXwhAWlpWHFwMD\nYWjfHhn9+kEoVw6Slxcy69TBo4MHoevdG0yptHouyd/f/oLxPJLWrkXK1KkwvPSS+VnAhx8iYds2\nMGqx4jAa/dBJSvMIbKX52oCiuz7Vn3/C5+OPrTZP1HfogKQ1a8ClpkLz88/gU1Oh79YNYvXqAACv\nsWOh3bzZ4XOLfn5Inj4dmt9/h5tcjsS2bWHo0gV4Km8NALKbNyGPioKpVi1IlSpZP6AgwL9LFygv\nXLBYLHl7I+Wrr6B/6y2Hy1pYZeH7mR/KkRNSRNT799tsYy6PjARg7rSkGzw41/q0ceOgPHkSiqfG\nbTFVrGhXaxWxYkUY27WD+sgR4Nw5eJ4+Dfc1a5D+wQfmgJ5z26AgiEFBeV+MXI7ETZvgPWYM5Neu\ngTMYIJYrh4y+fYs1iBMzCuSEFId8ei9KgYFI/PFHeE6bZh6Ei+NgCglB6ldfQXXoELzHj7eZQxf9\n/JDx/vvwHTQIqogIAI//0GNi4BUVBUgSDG+8UeAiS35+SFy/HlxGBjidzpxSoV6YJQIFckKKSMaA\nAVDv2WN1SFsuLQ2+774LoVo1pI8eDcnKg0+xShUkrV6d9ZpPTIT73LlQXL0KSau16HnJYM6Lmxo0\nQPrQoZA8PaH4559cLU1kCQlwX7XKoUD+BNNqnd6EkxQOBXJCiohQpw4Mr78OzY4dFmNxM56H/OFD\nyB8+BI4cgSoiAknffw/hhRdsHot/9Ah+PXtCcf26xXFEPz+Y6tSB8eWXzS1SHne79/jyS5vjv/AP\nHwImU64eljlxOh2UR4+CKZXIfOklwMaDTlIyUCAnpAilzJkDY6tWcNu6FVxGBhSXL+fKmyuio+H1\n+edI2LbN5nE8v/jCIogD5rHHuYwMZHz4IYytW1usEytUsNn2m6lUeU6SrF28GNrNmyGLigIeT0OX\nPmIE9N2753e5pJhQ80NCihLHwdC1KxI3boS+a1ebDz9l0dHgrHQiekJ+44bV5bxeD83WrbmW6999\nF+Jzz1ndx9Swoc3cturPP+GxdCnkUVHgkD0NnefXX0P2+AEtKXkokBPyrAiC7XWMOT5XJZ/7z5i5\nuyNl0iQIOYK55OYGY/PmSMljcCrtmjVWe3vKHj2Cx8KFjpWPFDlKrRDyjBjefBPCihWQW2nzLFau\nDOblZXNf0wsvQHn5cq7ljOehb9/e6j7Gzp0R9/LLqLBzJ9KvXIGxQwcYw8PzbGlibWiBrHWJiTbX\nkeJFNXJCnhHJ39/c6eep8ciFSpWQlk/39NTPPoMQGJir+zwnSfBctMhmbZ95egITJiB11ixzHj2f\n5oKin5/tdbY6C5FiR4GckGcobdIkJC9YAEN4OIxNm0L3xhtI2LwZmc2a5bkf8/aG5O9vc+Z59a+/\nOqd8w4dbHV9dqFwZaaNGOeUcxPkotULIM2bs2BHGjh0LtpMogrfRm5PLzITyn39g6Nat0GUTXnwR\nKdOmwX3ZMvNE0pIEyGSQfH3hOX060kaOhPj884U+D3EuqpET4gpkMkienlZXMbkcprp1nXYqw5tv\nIn7fPqSMHQuoVJAlJ0N54QLcduyAX8+eUO/e7bRzEeegQE6IizC2a2d1BEIhKAh6J9TGLTAG9x9/\nhOypqeTkDx/CfcECwMnjpJPCoUBOiItI+/RT6N55B8Ljh46SVovM+vWRuHy503teyq9cgSw62vq6\nW7egeGoURFK8HM6RS5KEVatWITo6GgqFAkOGDEF5e+bsI4Q4hueRMns2uMREKC5cgOTvb557tAgG\nruJE0Wa7do4xqpGXMA7XyE+fPg2TyYTp06ejd+/e2LBhgzPLRQixgfn6IjM83Dw2SxGNPmgKDbXZ\nM1SoXh2m+vWL5LzEMQ4H8qtXr6JBgwYAgODgYNy8edNphSKEFDO5HOlDhkB8avYf0c8P6R98kOeA\nW+TZczi1otfr4ebmlvWa53mIogiZlRlInrBnpgtXVpqvrzRfG0DXZ9WoUUBYGDBzJpCQAPj6QjZ+\nPHxatIDtqaOLR2n//PLjcCDXaDTQ5xyak7E8gzhAU725qtJ8bQBdX54qVwaWLLFcVsLeq7Lw+eXH\n4UAeEhKCf/75B2FhYbh+/TqqVq3q6KFIEUhN5fD991pcv65A1aoChg7NgJ+fg4MyAbh2TY4FC9yR\nkCBDQICIUaPSULMmPfAipCRwOJA3bdoUFy5cwJQpU8AYw8cff+zMcpFCuHRJjqFDfXHrVvbHu2eP\nBrNmJeOVV6xPD5aXLVuAkSP98OhR9i+uEyeU+OqrFLz+utEpZSaEOI5jjD09Dk+RKe0/f0rK9XXp\n4o9//83drrhWLRMOHIizNuqpTYIAdOpUEZcu5V5Xp44J+/cX7HglUUn67IoCXZ9rsye14uJ/guRp\n0dEy3L5t/VnF7dsynDlTsNYG584pcO2a9XVRUTJcuULD9RBS3OivsJTJyOBgY3J1GI0c0tLsv3df\nvy7D6NHeMJmsr5ck4NQpJebP94AgAK1bG/Huuzqa3pGQZ4wCeSkTHCygUiUJ16/nrpU/95yAZs3s\ny5GLIvDxxz64dct2DV6pZJg50xM6nfnm8NdfauzYocHmzYnQaJ5Zxo6QMo9SK6WMXA706aODl5dl\nCxWtVkKXLgakp3OYMsUTgwb5YPp0DyQmWu8ZuH+/GpGRtu/znp4iDAYuK4gDgCBwOH1ahTlz3G3u\nRwhxPqqRl0KDBmWgXDkB69a5IzmZh4eHhLff1iEgQELXrv64ezf7Y//9dzWWLElCgwaWM8xcvSqH\nyWT9Pq9SMdSpY8Lff6utrj99WgXA9pRhhBDnokBeSnXubETnztlNAwUBaN8+wCKIA0BUlAJTp3rj\n118thytt0cIIrVZERkbuFE2NGiZ4edlOnTg6hzAhxDEUyF1cXByP7ds1YAzo3l2PgADrUfTUKSVu\n3bLemiUqSoZr1+T44Qc3nDmjhCgC1asLqFlTxPnzlvuoVBK8vKTHrVUYYGXyseeft/F0tIikpnJY\nvVqLyEg56tY1oX9/HeXoSZlCgbyEuHhRjiVLPJCczKFCBRFjxqSjShVzz0mjEdiwwQ3/+58KSiXQ\nr186XnrJhDlzPPDTTxo8eGD+GOfNc0evXnp8+WXuKcGSkjiYTNbz4RkZHEaM8MalS9nNTS5dUuK5\n50xo2dKA6Gg1UlMl+PuL4HngzBkVBMH6sZ5/3oTJk59dWuX0aQVGj/bG7dvmh7K7djFs2eKG779P\nRHAw9TwlZQMF8hJgwwYNvv3WE/Hx2bXf48dVmD8/GXXrmvDuu364cEEBUTQHz4gIFZo3N+Lvv1UW\nzQkzMmRYs0YLhYJh6tTsYMoYsGqV7QeQJhOHK1dyt06JjlYgMFDEuXPAlStxMBqBnj39rQZxNzcJ\nnTvrMWlSms1fBTkJAsDzKFRnIkkCJk3yygriAMAYhxs3FPj0U2/s2pXg+MEJcSHUaqWY6fXAypXu\nFkEcAO7elePrrz3x5ZeeOHtWmRXEASAtjcdff6mstglnjMP69Vqkppq3P3pUgbCwQJw6pYS1NAgA\nSBIHSbK+7t9/VTh0CNBoGNav1yIhwXp6xtNTwhdfpOYZxNPSOBw5osDbb/vhpZcC0bJlIPr398G9\ne459Dc+eVVgMQ5DTrVty3LmT9yBuhJQWVCMvZhERKkRFWf8Y7tyRIT7eepATBNvBT6/n8OOPboiM\nlGHrVq3NNIg9RJHD4MGAQhHwuCzW8+IKhblduTUnT8oxc6YXbt+WITFRZnHTiI6WIyZGjt274+Hu\nXrC8dnIyD4PB+vtgMHBZNzNCSjsK5MVMFDmbrTzS0ngkJTlyVA6//qrG+fNKMJZ/MPP0FKHTcTZv\nDikpAJB37TYkRIBCATx8yMPDg8HNjSE6WoYxY7xw+rTK4hfF065dk2P2bA906WJAcLAJ3t72BfTG\njTNRpYqAmJjcX+NKlUQEBwtW9iKk9KFAXszCw42oVk20WisXBMBWOkSlkiAInI0AyewO4jzP0KWL\nAZcuyXH+vKpghc9xjNu3edSsWQGCACgUDPXrZyI+XpZnz9BsHDZu1GL9ei0qVBDRsqURc+akQCYD\nDh5UYdUqLVJTeXh6Snj//Qy0b29uVunlxfDaawZs3OgGvT77JuThIaJHDxoqgJQdFMiLmVbL8N57\nGViyxB3Jydm1XoVCstkhB2Bo0SITAQECtm7VInew55D/mJYMHh4S/u//0jB4sA4JCRxatSqH1NSC\n56slCbh5MztqGo0cTp2y3lnIlictamJj5fjlFxnc3BiqVxfw7bceFu/LhQtKjByZho8+ygAAfPZZ\nKipWFLF7txqpqTx8fCS8+64OPXtmT3qi13NIS+Pg5ychn7lPCHFJNIytkxR2KM3jx5X4/ntzzTMw\nUEJkpBxXr1qvzdaoYcKhQ3H4v//zxrZtbla3yQ/HMWzenIBWrTKRnMzhvffMLWNy5q/NNxMOtn4V\nFKWgIBM4DoiMzP0eBAWZh+NV5fMDIi2Nw9ix3rhwQQG9noO/v4hu3fQYNizDYruyMAwqXZ/rKtIZ\ngohzvfRSJl56KXtAq+HDvW0EcoZevcyBqDAP8xjj0KuXH+RyZrPVilYrIT2dL9TDUkclJPBISbF+\n3uhoOS5cUKJJk+z3SxAAkwnQaMyvGQMGDvTFiRPZ0T4uToY7d+TgeWDo0IynD0uIy6LmhyVUv34Z\n4HnrP5b++kuN9u0D8NdfBUtf5GZ+wGmr6WFysjzP1jGOk9CokRE9emRAo7H+pFcmg80cP88zqFTm\n9yYxkcOgQT54+eVAtGoViK5d/bF3rwqnTilw8WLuekpGRnZPWEJKC6qRl1A7dmhsBFjucS2zuJvW\nWW+GaA83N2DNmiQkJnLYty/3zUguZ4+bMlo/vkbD8MILJggC0KePH86fz87P378PTJzojbAwI9LT\nrSfEExN5pKVx8PSkaE5KBwrkRUyv5zB3rjtOn1ZCkjjUqmXCpElp+U6EnPdDx4IGUMeDbt5lcOy4\nOh2PY8eUmDvXA2lpTwdbBjc383gu9+9b3796dfNQAb/8orHaIzUuTobLl+WQy5mNXqjm5pGElBYU\nyIuQ0Qj06uWLM2ey87Tnzilx7JgK27YloHJl22OBdOxowG+/aWyOj1JwRRXMHdiLY5gwwctKEDcf\nMzVVlueNLDWVw5kzCuzdq87j/eFQo4aA69efDvQMjRplQk7ffFKKUI68CP3wgxvOns3dmDk2Vo52\n7QLw3Xe2xz95/XWDA6MI2qplFk/LE1sY42wE8Zye1Phzu3VLgW7d/HHggO1mKxoNw7x5SQgJMUGh\nMB/H21tEmzZGzJqV4mDJCSmZqF5ShCIi1DZ7NKal8Vi+XIuqVUV0767PtZ7ngR9/TMBLLwVaHRPc\nupITrAFAo5HAmLm7vGNl4+DmJiEzk8uVIjE/P7B1TIaXXzaiYUMBf/wRh4MHVYiOlqNVKyPq1KHe\nnqT0oUBehGSyvPOw6ekybN7sZjWQr1ypxaZNT3osWqZF3N1FpKfzKGmB2xKDj4+Ee/cK9xWrUME8\nfMD9+/YfR6lk6N/f3LxQLgc6djQCMOa9EyEujFIrRejtt3VQqwv+UPPAARXmz3dHZKTCouapVDK0\nbauH81MlDICzp/XhcP++Pb8k8r7Z3bols/M42Tw9GfXgJGUKBfIi9NprRnToYIBMZjtIenpKePiQ\nR1JSdmBevVpr0S39icxMDu7uEnS6onho6fyvgj1jvdSqlQm53Pb7w1jBf3lUrSqifHmab46UHRTI\nixDHAUuWJGPhwmS4ueUOLBqNhHv3eHToEIBXXw1Ajx5+uHlThpQU2x/LwYOaUjMnpp+fiJUrk/Dm\nm/rHv1xYjn+OqVhRwIQJuWdIIqQ0oxx5EeM44M03DahSJR5Tp3ojKkoOoxHw95eQlMQjOjq7edz9\n+3K8/74vAgNtN0vMyCgd916OY9BoGLp0CUBycs5ad97NJOVyBoWCITMTEEUeHh4iVCogKEhAxYoi\nxoxJQ40aNMUbKVsokD8jjRoJ2LMnHjdvyqHTcZgxwwOxsbl7NUZGyhEamglvb9FqeqW0YIxDbKy1\nr1/eaRSOA3r3zkCLFpl49EiGli2NCAoSYTQCMTFyu8cyJ6Q0oUD+DHEcULOmuflbUpKtmjWHzEwe\no0al44cf3BAVJXs8+URJbqHy7JhMHFav9sC//2Zi69Z4qNXArFke2LtXjUePZNBqJdSqJWDhwmT4\n+paSHBQh+SjU7/RTp05hwYIFzipLmZJXF3F/fwkffJCBP/6Iw7ZtCQgNLWjHoNLv3DkFVq/WYsEC\nd6xapUVkpAKpqTzu35fjr7/UGDDAlwbGImWGw4F87dq1+PHHH/EMhzMvVd54Q2915L8KFQR88kk6\nAECpBBo1MqF1a1dpA10U3wXrx2SMw/HjKuzZo7aYHeiJ//6TIyLCsRmPCHE1DqdWQkJC0KRJExw8\neNDufewZIN2VFeT6Jk40j9S3cycQG2vuyRkUBEydKkfTpuUstp05Ezh9Gjhxwtkldi6etz3/qKO8\nvTkkJ1tfp1Sqba7T6XicO+eH3r3tOw99N11bab++/OQ7Q9ChQ4ewZ88ei2VDhw5FzZo1cfnyZRw4\ncACjRo2y62SlfRYPR64vLo7H3r1qeHhI6NjRkDUxwtP0emDVKncsW6ZFSkrpfQj6tJo1Tbh1S57r\nGYFSydC6tR6HDtkeWCwoyITDh+PA5/O7syzMMEPX57qcMkNQmzZt0KZNG6cUiOQWECChXz9dntsk\nJPCIi+MxeHAGIiJUOHGibARymYxhwIAM7N+vxt9/q7ICtlIpITjYhP/9T5Xn6JCxsTLs3atGp06G\nZ1VkQooFtVopwZKSOIwY4YMrV+RIT+fh7S0hObmstF5haNIkE3366NC3rw47dmjw++/qx+3y9diw\nQZvvCIpGI499+yiQk9KPAnkJxRjw/vu+OHUq+4GdeaCsvHEcs6trfPHIu7OPt7eI554ToFQCLVpk\nYuTINCge95fq0UOPHj2yBxdbvNjDrjP6+VHnIFL6FSqQh4aGIjQ01FllITn8848CV67k9/FkB0aO\nYyhfXkRqKo+MDKCkjYzIcQyenpLN/H7lygL274+zu0OPl1f+T1XLlxfw0Uc0yTIp/UpHf+9S6PRp\npc05J7Nx0GolvPqqHm+/rUNKCve4C3/JCuKAubmgKAIq1dMBmEGjkTB5ckqBemUOGJCeZzCvUEHA\n8OHpqFCBOgWR0o9SKyVUrVomqNUSDIa877VVqohYty4Jb77pB52uZD8EfXJjMqd/gCfD8er1wJgx\nPkhJSUHfvrnHZremUycj7t5Nw48/uiEmRgaFAvD3FxEUJKBWLQGDB2cgIICCOCkbKJA7gDFgzx41\nfvlFA0Hg0KqVERMnOvccr7ySieefF3DxYu6p4nKqXNnc5f/mTdf5KHPn8Dno9RymTfNC48Ym1K5t\n3yw+H36YgQEDMnD1qgIaDUPNmgK4kvdjhJAi5zp//SUEY8DHH3tj/341jEZzbTkiQoUDB4D16zlo\ntc7p3chxQNWqAi5fVtgcZ6VyZQETJqQBQKkY2tZg4LFokTuWLrXRy+cxvZ7DTz9pEBUlR+vWRrz8\nspECOCnTKEdeQH/+qcKBA9lBHDDPH3niBDBnjn0tKexx4IAKf/2lthLEGfz8BLRpY8CaNYlZtdca\nNUrSXJSO38wSE/P+Sh4/rkT79gGYOtULK1e6Y/BgH7z1lh9SUymSk7KLauQF9NNPblbH9gCAf/7J\nOw1ij9u3ZVi5UosDB9TQ6aydh8Pzz4vYuDExa4kxayiWvJv3PSscB4cHrCpf3nZzQaMRmDjRC7du\nZX9tdToep06p8Omn3lixIsmxkxLi4iiQF1BeKYzCpjdmzPDAli1uSEjI+6GlIJhrrmfOKODvL2H9\nei3+/VeJkhDEAcDbW0JSUsEfvPI8y3NSjT171IiKsv6VvXhRAaMRUNE4WaQMotRKAb3+uh5KpfXq\nZp06jg83+++/CmzcqM03iANAcjKPDh38MXCgH955xw+7dqlRPEE89/vg4SFi+PB0VKxoT6rHcn9J\n4rB5sxYXL1oP1rGxcoii9es0GDhkZJSMGxkhzxoF8gJ6800DmjUzguMsg9ALLwCTJqU5fNxVq7RI\nTc3/4/DxEREdLcO9e+Zgp9fzMJmK52NUqSQ8/7y5mSTPM9SoIeDDD81jozx69OSGZJ6DU6m09nMl\nd+BNTJRh2TJ3q+dr184Ab2/rNfaAABE+PjSkMimbKLVSQDIZsGFDIpYudcfRoyoIAhAaasK337rD\naHQ8t6LT2a5NajQSqlcXULOmgP/+kyMpqfC5+Lx4eIgoV07C7du2a8BmHHbujEd0tBwZGRwaNjTh\no4+8cfJkzjRPds/TgAARcXH5d1hKS7N+Y6pTR0CjRiYcOsRbNGH09BTx7rs6arlCyiwK5A5QKoFR\no9IxalR61jI/P3cUZiTNpk0zcfCg2uo4Ka++asCyZeYmec2bBzp+kjwxyGSAKHLIzOSgUjGEhJhw\n5Yrtm4bRyGHDBi1GjMh+H44fV8FaoDYaOaSmwuq6p4WE2E5RrVyZiKlTvXDypBJ6PQc/PwnvvZeB\nPn3s60hESGlEgbyEGDhQh127NLh0yTJwVq0qYNy47JSNj4+EmJiiKcOT2rfRyOPyZSUqVjShUiUB\nd+/a+ppwOH9ekfXKYEAew8pyMBrzD+JBQSYMG5Zuc71KBcyZkwLGzA99FQrzQ2ZRNP9aIqQsohx5\nCaHRMGzenIA33tAhONiEoCAT2rQxYO3aBFSvnp0X7tZNnys/by+el6DVinj6IaNcLsFaTfnBAzne\nfluHgADbx/T0zE4nKZWFC6ZeXiLWrk20K9fNcUBKCo/Bg33QqlUgwsIC0b27H/73v6JNOxFSElGN\nvATx9WX59mr88MMMLFjgjuTkgkVMnmd45x095sxJwfbtavz0kxaSZJ7o+Z9/FHjwIPc9XZI4PHwo\nw8SJwIQJEjIz+VzHbN3akOM1EBxsyvWrwl4ymf03AoMBeO89X4tzxcbKMWKEDCtXJqJhw5LUQYqQ\nokU1chdUt671HDLPM6jV1lp1MLRrZ8C336ZAJgN69DDgl18SsH17Ar7/Pgnly1t/SCuXmyd3GDXK\n+rCxksRh0SIPi/bzP/+cgEqVBDjSuzMxkce9e/ZF8k2b3HDliiLX8vv35Zg3z7PA5ybElVEgd0Gf\nfpqaq502xzG89JIRFy8+QqdOelSsKECtllCxooAuXfRYutR2r8fu3XXQanMH6uefF9C9ux4PHgCc\njSYht27JH7dSMfPyYjhx4hEWL05CixYGFCSgly8vomZN+2rSJ0+qbI5B8+ABJctJ2UKpFRfUqJGA\nlSsT8c03nrh3TwalkqF580xMnJgKlQr4/vskJCbyiI2Vwd1dxNmzShw5okLr1kYorWQ93n9fh5QU\nHjt3anD/vgxubgxBQQLmz0+GQgEkJ5sHqrLGYOARH29ZH5DJgG7dDOjSxYCXXw5AdHTumnNuDE2b\nZiIw0L4mnHmNRa5Wl4IRxAgpAArkLqpBAwGbNiXaXO/jI2HRInfs2aPG3btyyGQM1asL+PTTVHTu\nbLTYVhSR1RlJo5Hg7S2hVSsjKlc2p2mCgsy1ZWvtuytWFNCiRabVMsjlwCefpGPaNE9kZORdS9Zq\nGebNy/v5QE4ff5yOgwfViIuzPK5cztC2rdHGXoSUTpRaKaU2bnTDpk1uWU0HRZFDZKQCX3zhhfv3\nzR/7tWtyjBjhjWbNArFypRaRkQrEx8tx44YSixa5Y8YM82iOSiXQrZsO7u6W+XelkuGVV4zw97dd\nA+7dW48NGxLRurUevr62x1Hx8ZGy5ue0R40aIoYNsxwKwNtbRMeOeot27YSUBVQjL6V27NBYHT3x\n3j05lixxx3PPiVi82B3x8dZrygaDeQb6MWPMbdhHjsyApyfDL7+4ITGRh4eHhHbtjBg7Nv9hCZo3\nN6F58yTcvClDt27+VseTqVpVhLyA38YPPsjAW2/psWmTG5KTOfTooUedOtRahZQ9FMhLKfPcndbF\nxMhw4IDaZhB/4sEDGW7fliMoyPx64EAdBg7UOVymoCARbdsasGuXxmI890qVBEycmOLQMf38JIwc\nSTVwUrZRIC+lzGkMa7kKBpPJ3OY6P+7uDL6+udMmjJn/8Q4k5ubNS0HDhibs3q1BQgKPxEQePM8w\nfrw3XnzRhGnTUqHR0OBXhBQE5chLqSFD0uHnlzsnHRQkoF49+4bbDQ42WbQxz8jgMGKEN155JQAt\nWgTijTf8sH9/wQYA5zigXz8dJk1KQVoah/h4GWJiFLhyRYlNm7To08e3VExbR8izRIG8lAoPz8Tk\nyakIDc2Ej4+IwEARzZoZsWpVIt57Tw9/f9sPHtVqCfXrZ2LRouxWJJIE9O3ri23b3HDzpgKxsXKc\nOaPC+PHeOHiw4LM5fPONJ+7fz/2r4Nw5RYFvDoSUdZRaKcV69tSjRw89YmNlUKuZRRvt9u0N2L5d\nA4Mh+17u7S3i5ZcN6NbNgHbtjBapk4MHVRYDZD0RFyfD0qXuaNeuYE3+7t61/ZD19981eO01akJI\niPkAhoIAAAamSURBVL0okJdyPG9uEfK0OXNSULeuCbt2aWAwcAgMFDF2bBpCQ623+ti3T20R9HN6\n+LDgPSltzbIE5N3ZhxCSGwXyMupJrrpfP/taoVSokFcqpuAPJxs3znw8Voplj9HAQBFDhmQU+HiE\nlGWUIyd2GTRIh8qVrdXWGZo0KXga5LPPUtGiRaZFd/rAQBEffpie1aOUEGIfh2rkOp0OCxcuhF6v\nhyAI6N+/P4KDg51dNlKC+PpKGDcuDXPneiAmxvy1cXcX8eKLJnzxRWqBj6fRmEdKPHhQhd9+08DL\nS8KQIRmoVImCOCEF5VAg/+2331C3bl106tQJ9+7dw4IFCzB79mxnl42UMN2769G2rQEbN2rx8CGP\nzp31aNbM5PBcmTwPtG9vRPv29GCTkMJwKJB36tQJiscDY4iimPV/Uvp5ezMMH049KQkpSTjGWJ5P\nqg4dOoQ9e/ZYLBs6dChq1qyJ5ORkzJgxAwMGDECdOnWKtKCEEEKsyzeQ23Lnzh3Mnz8fffv2RcOG\nDe3a515hppkv4SpWrFhqr680XxtA1+fqysL15ceh1EpsbCzmzZuHUaNGoVq1ao4cghBCiJM4FMh/\n/PFHmEwmrFu3DgDg5uaGcePGObNchBBC7ORQIKegTQghJQd1CCKEEBdHgZwQQlwcBXJCCHFxFMgJ\nIcTFUSAnhBAXR4GcEEJcHAVyQghxcRTICSHExVEgJ4QQF0eBnBBCXBwFckIIcXEUyAkhxMVRICeE\nEBdHgZwQQlwcBXJCCHFxFMgJIcTFUSAnhBAXR4GcEEJcHAVyQghxcRTICSHExVEgJ4QQF0eBnBBC\nXBwFckIIcXEUyAkhxMVRICeEEBdHgZwQQlwcBXJCCHFxFMgJIcTFyR3ZyWAwYOHChcjIyIBcLsew\nYcPg6+vr7LIRQgixg0M18j///BM1atTAF198gVatWmHXrl3OLhchhBA7OVQj79SpEyRJAgDEx8dD\nq9U6tVCEEELsxzHGWF4bHDp0CHv27LFYNnToUNSsWRNffPEF7ty5g6lTp6JatWpFWU5CCCE25BvI\n83P37l3MmjULixYtclaZCCGEFIBDOfIdO3bgyJEjAAC1Wg2ep8YvhBBSXByqkScnJ2PJkiUwmUyQ\nJAm9e/dGrVq1iqJ8hBBC8lHo1AohhJDiRTkRQghxcRTICSHExVEgJ4QQF+dQh6CCKu1d+nU6HRYu\nXAi9Xg9BENC/f38EBwcXd7Gc7tSpUzhx4gRGjhxZ3EVxCkmSsGrVKkRHR0OhUGDIkCEoX758cRfL\nqW7cuIEffvgB06ZNK+6iOJUgCFi2bBni4uJgMpnQvXt3NG7cuLiL5TSSJGH58uW4f/8+AOCDDz5A\n1apVbW7/TGrkpb1L/2+//Ya6deviiy++wLBhw7B69eriLpLTrV27Fj/++CNK07Px06dPw2QyYfr0\n6ejduzc2bNhQ3EVyql27dmH58uUwmUzFXRSnO3r0KDw8PPDll19i8uTJpe5v7syZMwCAr776Cr16\n9cKWLVvy3P6Z1MhLe5f+Tp06QaFQAABEUcz6f2kSEhKCJk2a4ODBg8VdFKe5evUqGjRoAAAIDg7G\nzZs3i7lEzlWuXDmMHTsWixcvLu6iOF2LFi3QvHlzAABjDDKZrJhL5FxNmzZFo0aNAABxcXFwc3PL\nc3unB3J7u/S7qryuLzk5GYsWLcKAAQOKp3BOYOv6wsLCcPny5WIqVdHQ6/UWfyA8z0MUxVITFJo3\nb45Hjx4VdzGKhFqtBmD+DOfNm4devXoVc4mcTyaTYfHixTh9+jTGjBmT98bsGYuNjWWffPLJsz5t\nkYuOjmajR49m//77b3EXpchcunSJfffdd8VdDKdZt24dO378eNbrjz76qBhLUzQePnzIJk2aVNzF\nKBJxcXFs/Pjx7M8//yzuohSppKQkNnToUKbX621u80xy5KW9S39sbCzmzZuHESNGoGHDhsVdHGKn\nkJAQnD17FgBw/f/bt0MjC2EojMI/JoKC8Lw2MFgoAMEwg6IE0ChcJh0gKAH5qkDgYGA7YFe8gc3u\n+RrIjTniZvJ+Xz4m4XdZlkVN0yhJEr1er6fH+bhpmuSckyQZYxQEwWU3b9mRx3Gstm01jqOO41CW\nZXcce5thGLRtm/q+lySFYaiiKJ4dCt+KokjzPKuqKp3nqTzPnx4JP+Sc07qustbKWitJKstSxpiH\nJ/uMKIrUdZ3quta+70rT9PJufNEHAM/9rR0HAPxDhBwAPEfIAcBzhBwAPEfIAcBzhBwAPEfIAcBz\nX/xUQK2x2MpwAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x118ccecd0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(x_train[:, 0], x_train[:, 1], c=clusters, cmap=cm.bwr)\n",
    "plt.axis([-3, 3, -3, 3])\n",
    "plt.title(\"Predicted cluster assignments\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The model has correctly clustered the data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Remarks: The log-sum-exp trick\n",
    "\n",
    "For a collapsed mixture model, implementing the log density can be tricky.\n",
    "In general, the log density is\n",
    "\n",
    "\\begin{align*}\n",
    "  \\log p(\\pi) +\n",
    "  \\Big[ \\sum_{k=1}^K \\log p(\\mathbf{\\mu}_k) + \\log\n",
    "  p(\\mathbf{\\sigma}_k) \\Big] +\n",
    "  \\sum_{n=1}^N \\log p(\\mathbf{x}_n \\mid \\pi, \\mu, \\sigma),\n",
    "\\end{align*}\n",
    "\n",
    "where the likelihood is\n",
    "\n",
    "\\begin{align*}\n",
    "  \\sum_{n=1}^N \\log p(\\mathbf{x}_n \\mid \\pi, \\mu, \\sigma)\n",
    "  &=\n",
    "  \\sum_{n=1}^N \\log \\sum_{k=1}^K \\pi_k \\, \\text{Normal}(\\mathbf{x}_n \\mid\n",
    "  \\mu_k, \\sigma_k).\n",
    "\\end{align*}\n",
    "\n",
    "To prevent numerical instability, we'd like to work on the log-scale,\n",
    "\n",
    "\\begin{align*}\n",
    "  \\sum_{n=1}^N \\log p(\\mathbf{x}_n \\mid \\pi, \\mu, \\sigma)\n",
    "  &=\n",
    "  \\sum_{n=1}^N \\log \\sum_{k=1}^K \\exp\\Big(\n",
    "  \\log \\pi_k + \\log \\text{Normal}(\\mathbf{x}_n \\mid \\mu_k, \\sigma_k)\\Big).\n",
    "\\end{align*}\n",
    "\n",
    "This expression involves a log sum exp operation, which is\n",
    "numerically unstable as exponentiation will often lead to one value\n",
    "dominating the rest. Therefore we use the log-sum-exp trick.\n",
    "It is based on the identity\n",
    "\n",
    "\\begin{align*}\n",
    "  \\mathbf{x}_{\\mathrm{max}}\n",
    "  &=\n",
    "  \\arg\\max \\mathbf{x},\n",
    "  \\\\\n",
    "  \\log \\sum_i \\exp(\\mathbf{x}_i)\n",
    "  &=\n",
    "  \\log \\Big(\\exp(\\mathbf{x}_{\\mathrm{max}}) \\sum_i \\exp(\\mathbf{x}_i -\n",
    "  \\mathbf{x}_{\\mathrm{max}})\\Big)\n",
    "  \\\\\n",
    "  &=\n",
    "  \\mathbf{x}_{\\mathrm{max}} + \\log \\sum_i \\exp(\\mathbf{x}_i -\n",
    "  \\mathbf{x}_{\\mathrm{max}}).\n",
    "\\end{align*}\n",
    "\n",
    "Subtracting the maximum value before taking the log-sum-exp leads to\n",
    "more numerically stable output. The $\\texttt{Mixture}$ random variable\n",
    "implements this trick for calculating the log-density."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}