gibbs sampler

gistfile1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Gibbs Sampler for misture of Gaussians"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We implement a Gibbs sampler for detecting the centers of a specified number of clusters based on a Gaussian mixture model. Work done by Rachel Levanger for Fall 2016 Foundation of Graphical Models at Columbia University, HW #2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate the sample dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We specify K clusters and randomly select their means from a Gaussian distribution centered at zero with high variance (to spread them out). We then randomly choose a membership according to a uniform distribution for each datapoint, then generate that datpoint's position according to a multivariate normal distribution centered at the selected mean with some pre-specified variance common to all clusters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Generate the means of the Gaussian mixture"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Set the number of clusters\n",
    "K = 3\n",
    "\n",
    "# Set the variance of the distribution of means\n",
    "p_lambda = 1000\n",
    "\n",
    "# Randomly choose K means according to N((0,0),p_lambda) distribution\n",
    "mu_x, mu_y = np.random.multivariate_normal([0,0], p_lambda*np.identity(2) , 3).T\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Generate the sample data according to the generated means of the Gaussian mixture"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Set the number of datapoints in the sample\n",
    "N = 100\n",
    "\n",
    "# Set the variance of each of the clusters\n",
    "sigma = 50.\n",
    "\n",
    "# Initialize the dataset\n",
    "x=np.zeros(N)\n",
    "y=np.zeros(N)\n",
    "\n",
    "# Loop through each datapoint and assign it to a cluster, then generate its position\n",
    "for n in range(0,N):\n",
    "    k = np.random.choice(K)\n",
    "    x[n], y[n] = np.random.multivariate_normal([mu_x[k],mu_y[k]], sigma*np.identity(2) , 1).T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot the sample data and the means"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-30.0, 70.0, -50.0, 40.0)"
      ]
     },
     "execution_count": 108,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEACAYAAACj0I2EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X+wXHWZ5/H3YyKJIzNe2a0JCJGkopTEuOCPAZaMscfh\nJkhY0NqSMEumUFMrtcGEcidcA0FzKZJaxIyjRLCchdmiRCCsuhTcK4Sr0qQ2rqODGIiBRah7A1HB\n2SVNSc0m/Hr2j+/p2337dvftH+d0n3Pu51XVdfucPuf0OTedp7/3+T7n+zV3R0REsu9N/T4BERGJ\nhwK6iEhOKKCLiOSEArqISE4ooIuI5IQCuohITsQS0M1sjpk9amb3RcvHmdmYmT1lZg+a2UAc7yMi\nIo3F1UK/AjgAlIvaNwNj7n4K8KNoWUREEtR1QDezk4DzgFsAi1ZfANwWPb8N+Hi37yMiIs3F0UL/\nO+BK4I2qdQvc/YXo+QvAghjeR0REmugqoJvZ+cDv3f1RKq3zKTyMLaDxBUREEja3y/3PBi4ws/OA\n+cCfmNm3gRfM7Hh3f97MTgB+X7ujmSnIi4h0wN3rNqC7aqG7+9XuvtDdFwMXAz92978G7gUujTa7\nFLinwf6xPLZu3Rrbsfr5yMN15OEa8nIdebgGXcf0RzNx16GX3+16YNDMngI+Gi2LiEiCuk25THL3\nh4GHo+cvAufEdWwREZlZLu4ULRQK/T6FWOThOvJwDZCP68jDNYCuox02U04msTc28369t4hIVpkZ\nnkSnqIiIpIcCuohITiigi4jkhAK6xGJ0FEqlqetKpbBeRHpDAV1isXw5bNlSCeqlUlhevry/5yUy\nm6jKRWJTDuJXXglf+Qps3w4DGglfJFbNqlwU0CVWExOweDGMj8OiRf0+G5H8Udmi9ESpFFrm4+Ph\nZ21OXUSSpYAusSinW7ZvDy3z7dun5tRFJHlKuUgsRkdDB2h1zrxUgr17YfXq/p2XSN4ohy4ikhPK\noUvqqG5dJH4K6NIXqlsXiZ9SLtI3qlsXaV9iOXQzm0+Y1GIeYbKM77r7sJkdB+wCTgYmgIvcvVSz\nrwK6qG5dpE2J5dDd/QjwF+5+OnA6cK6ZnQlsBsbc/RTgR9GyyBSqWxeJV9c5dHf/l+jpMcCbCfOK\nXgDcFq2/Dfh4t+8j+aK6dZH4dZ1DN7M3Ab8AlgDfcPerzOywu789et2AF8vLVfsp5TKLqW5dpDPN\nUi5dTxLt7m8Ap5vZ24D/YWbLal53M6sbuYeHhyefFwqF3MwdmEW9DrD1jjkwoGAuUqtYLFIsFlva\nNtYqFzP7IvAvwH8ECu7+vJmdADzk7u+p2VYt9BSpToEMDExfFpF0SKxT1Mz+tZkNRM/fAgwCTwD3\nApdGm10K3NPN+0jyBgYqeeyJifaDuW4UEum/bssW30fo9JxD+HLY5e7borLFu4F3orLFTOm0jFAt\nfJHeSLJs8XF3/4C7n+bu73P3bdH6F939HHc/xd1X1gZzSaduygi7beHXo1a/SHt0p6gA8bWw47xR\nSK1+kek0OJfMaO/eSqAst4C3bw/robWWcdw3CiXR6hfJM7XQZZpOWsZJtqY1PIBIhVro0pZOWsbV\nLfzqY5Rb+J3S8AAirVMLXaYp32RUKlVaxgMDvb+LUzl0kek0Y5G0pVSCTZvC82uugW3bwvMdOzoL\npJ3eharhAUSmU8pF+qrTySxWr57+BaLhAUQaUwtdpkki5aLJLETioZSLtC2JAKxqFZHuKeUibUli\nrHJVq4gkTy10mSbuzkhVq4jERykXSVztl8DoKCxbBvv3h+VyB2j5S0HVKiKdUcpFEldbybJsGZx/\nfvhZDtrVlS2qVhGJnwK6xKL27tIbboCRkfBT47CI9IZSLhKr2koWVbaIxEspF+mJ2kqWgwdV2SLS\nS91OQbfQzB4ys1+Z2X4z2xitP87MxszsKTN7sDxNnaRLnBNI1JY6Dg2FHPrQUHyljyLSXLct9FeB\nz7v7e4GzgMvN7FRgMzDm7qcAP4qWJWU6vSW/ntrRFvfvDzn0cpVLXKMvikhjsebQzewe4BvR4yPu\n/oKZHQ8U3f09Ndsqh54CuiVfJFt6UoduZouAh4FlwLPu/vZovQEvlpertldATwl1XIpkR7OAPjem\nNzgW+B5whbv/IcTwwN3dzOpG7uHh4cnnhUKBQqEQx+lIG2o7MuNuoe8ZHeXBG29k7tGjvDZvHis3\nbmSFCtBFWlYsFikWiy1t23UL3czeDIwA97v716J1TwIFd3/ezE4AHlLKJX2SviV/z+gou6+4gu3P\nPDO5bsuSJaz6+tcV1EU6lFjZYpROuRU4UA7mkXuBS6PnlwL3dPM+koykpo0re/DGG6cEc4DtzzzD\n2M6d8byBiEzRbcplObAWeMzMHo3WXQVcD9xtZuuACeCiLt9HElCvkRznLflzjx6tu37OkSPxvIGI\nTNFVQHf3/0njVv453RxbsqHZyIyvzZtXd5/X58/v0dmJzC66U1S60qyWfeXGjWxZsmTK9lcvWcLg\nhg19OFOR/NNYLtK1ZrXse0ZHGdu5kzlHjvD6/PkMbtgwY4eoJocWaUzjoUviOq1lHx2Fl1+GVasq\nAfzgQfjsZ+Ezn4E1azQhhkg1Dc4liepmernly2FsDDZtCvuVSnDddfCnfxrWa+hdkdYpoEtdrQzc\nNToaWtO1g3KtXdt6UB8YgMFBOHoUNmwID4Bt2+Dss0Or/8orFcxFWqGALnW1MnDX8uWwfj186EOV\nbW64AW66CXbvbn3UxlWrws/bbw+PjRtDK/0nP9HQuyLtUA5dGppp4K7yvKHXXRda2PPmwRe/CD/9\nKezZ03qapFQKLfOf/QyWLoUDB+CMM2DnTk0qLVJLnaLSsWadneVAe9llcNppcNFFcMwxIbDv2NF6\nMN+0CZ59Fr70Jfjwh+HP/zy817ZtYfhdTSotUqGALh1pZWjdgwfDRBbf/CasWwdPPdVepUu5yuXU\nU+GSS+A734Ebb4TTT4cHHggpGLXKRSoU0KVtrQzcFUcLvfo4Q0MhB3/ZZSG4j4zAyScnd40iWaSA\nLm1r5eaeuHLo1e9VTvHs2wfPPdebFItuZJIsaRbQcfe+PMJbS5YdPuy+fr37XXeF5+Xl8vORkc6O\nNz5eOU4vVJ93vWWRNIliZ/242uiFpB8K6P01MjI9YLUbhLdudZ+YmLpuYiKsb1e/g2q/vkxE2qWA\nnlFxBN1G4gigcQbhJK+1VePj4X/E+Hjv3lOkXQroGZV0qzWOVmleWrZ5uQ7Jv2YBXZ2iKddK6WA3\n4pggOuuTTCc9FZ9InBIdnMvM/sHMXjCzx6vWHWdmY2b2lJk9aGb6b9GhgYEQzJMY06SbQbU6Ocbo\nKOzaNXWbUimsa3WYgCQkPRWfSM80arq3+gA+DLwfeLxq3Q3AUPT8C8D1dfZL+i+TXEgqFZB0Dr1e\nTnxiwn3lSvd16yqVMOvWVZZFZGYknUMHFtUE9CeBBdHz44En6+zTg0vPtiRz6ElXuTQ694mJEMDX\nrg0PBXOR9vQjoB+uem7Vy66A3rI0VH40M9MXTqO/LsrVJKooEWlfs4AeS6eomS0C7nP390XLh939\n7VWvv+jux9Xs41u3bp1cLhQKFAqFrs9FemumTtvaDtPyYFxHj4bX2x0mQGS2KRaLFIvFyeVrr722\nYadoUgH9SaDg7s+b2QnAQ+7+npp9PI73lv5rVOVSG+yHhsIwARCCOITgXl5WUBeZWT+moLsXuDR6\nfilwT0LvI33WqMqluvRv0aLwc/36MGZKOXjv3RvGfhkcrFSU1M6KJCJtaJSLafUB3An8FngFeA74\nNHAc8EPgKeBBYKDOfsklmaQn2q1yqc3/9/t2f5EsQjcWSRLiGKUw6RunRPJGw+dKqmX9TlORXupH\nDl2kJXHcrSoigQK69E2pBGvXhuqXcsfpli1hWjt1jIq0TwFd+mbvXrjppjDtXKkUcudDQ5VqGBFp\nj3Lo0nfqGBVpnTpFpSe6qXpRx6hIa9QpKj2xfHloadfeXDRT+kQdoyLxUAtdYtVu+kSTS4i0RykX\n6al20idx3JzUz+PPBvodpotSLtIz7aZPVq+e3hIfGIgvUHSaBpIK/Q4zpNGYAEk/0FguuZPWsVk0\nAXT39DtMDzSWi/RCmv80VxVN9/Q7TAelXKQnkk6fdEpVNN3T7zAbFNAl1+qNy16dD5aZ6XeYHUq5\nSK6lOQ2UFfodpktfyhbN7Fzga8Ac4BZ3/3LN6wroIiJt6nlAN7M5wP8GzgF+A/wc+Ct3f6JqGwV0\nEZE29aNT9AzgaXefcPdXgbuACxN6LxERIbmAfiJhftGyQ9E6ERFJSFIBXbkUEZEem5vQcX8DLKxa\nXkhopU8xPDw8+bxQKFAoFBI6HRGpRxUs6VcsFikWiy1tm1Sn6FxCp+hfAr8FfoY6RUVSp5PRLvUl\n0F897xR199eAzwG7gQPArupgLiLpMDBQuVFoYqK1oYs1WFd66cYikT5IWyu33XFaNG1g/2gsF5GU\nSVMrt91xWkZHw88rrwxfAldeOXW99I9a6CJ9koZWbic59FIJNm0Kz6+5Bj77WViwAHburOyjnHpy\n1EIXSaGBgamt3H6kLPbunRq8yzn1vXtbP8aCBfDII/DSS2FZOfX+UQtdpE/S0ELvRDn/XypV8u5m\ncPnl8I1vZOtaskhzioqkTNYnx673ZVQd4DUBRnIU0EVSJm1VLu2o92VUnVNXCz1ZCugiEpvaL6Ny\nQB8chDVrsvfXRtYooItIYrL810YWqcpFRBIz01yyo6PTa9tLpc7q1uM8Vh4poItIouK8iSpNN2Sl\nkVIuIpK4OEs0s1ruGRfl0EWk79odL6ZXx8oa5dBFpK/aHS+mV8fKGwV0kZyKuwOx0+NVlzEuWlQZ\nrreTQBznsfJIAV0kp+LuQOz0eHGMF5PEsfJIOXSRHIu7A3G2d0imQSKdomb2SWAYeA/wZ+7+i6rX\nrgI+A7wObHT3B+vsr4Au0gNxdyDO5g7JNEiqU/Rx4BPAnpo3WwqsAZYC5wI3m5lSOyJ9EHcHojok\n063jQOvuT7r7U3VeuhC4091fdfcJ4GngjE7fR0Q6E3cHojok0y+JlvM7gENVy4eAExN4HxFpIu4O\nRHVIpl/TgG5mY2b2eJ3Hv2vzfZQsF+mxmcZYqdZKSWI7x5P+mNvsRXcf7OCYvwEWVi2fFK2bZnh4\nePJ5oVCgUCh08HYi0q1ySWK9CTc0mmJ/FYtFisViaxu7e1cP4CHgg1XLS4FfAscAi4FniKppavZz\nEUmPw4fd1693Hx8PPw8fnrq+0XI7Rkam73f4cFgvrYliZ9143HEO3cw+YWbPAWcBo2Z2fxSlDwB3\nAweA+4H10UmISAP9GBZ2z+go16xaxXChwDWrVvHY3tG6k1aXc+VbtoSSxXLLfe/e9s9ZoyUmrFGk\nT/qBWugik+JsBbfi4ZERv3rJEneYfHxh8RJfc97ItBZ62fh42HR8vLtzbvSXgLSGJi10BXSRlOhl\noNuycuWUYD4Z1D+6asq51AbrRumYds+59stBWtcsoOuGH5EeaCWlMjBA3ZRHEuYePVp3/fzXj0ye\nS3VapVH9eSfnrJuTkqOALtIDreSOexnoXps3r+761+fPn3xeLklsVn/e7jnr5qSENWq6J/1AKReZ\nZZqlJ3qVQy9XmdTLoV+1ZIk/3Ea5SSfn3EmViypjpqJJykWjLYr0UKOBrZKq9a49bqkEmzbB4CCc\ncOwoP/jqTp594gjvPHU+5/3nDaxo4816VZ9e3aqvrZGfjSM9NhucSy10kR5ptwMxjpZpvVb0unXh\nkaUqE1XGVKAqF5H+6iQ9EVcapl4wzGKVSRbPOQkK6CJ91mlrO66WaXUwzGJrN4vnnJRmAV05dJGU\n63ZCiepZhrZtC+t27MhOPlo59KmSmuBCRBLWbSljbZngYM1we/WGwK1XM79rV3jUHjvJoQnKNGxv\nGxo13ZN+oJSLSFNx5NA7SfU060jt1dAE0hhKuYhkTz+Hra03GTRogug0SGSS6G4poCdH41dLHOrl\n7mfK5+uzlzzl0GcZDVGab70Yarde7r6VfL4+e33WKBeT9APl0BOlMq/8SnrCiYkJ99WrO8+h67OX\nLJRDn526LXeT9KqX4+4kn12vJPCSS+Dmm+HkkyvblStc1qyZum+jVIo+e8lJ5NZ/4CvAE8A+4PvA\n26peuwr4NfAksLLB/sl/lc1iaiXFL22DRMV152TcnxV99pJFEneKAoPAm6Ln1wPXR8/Lc4q+GVgE\nPF3ezhXQe6LXs9/MFmn6vcYdNOP+ckjD7yivEgnoPjU4fwK43Sut8y9UvfYAcFadfXpw6bNT2lqS\neZKG1mejoHnXXeFROyzvXXe1VncexzXps5e8ZgE9lhy6md0H3Onud5jZTuCn7v6d6LVbgPvd/Xs1\n+3gc7y3Sa53mh+Mq6Wt0nN27YWwsLO/YEX5u2lRZrpdj12312dMshz53hh3HgOPrvHS1u98XbbMF\neMXd72hyqLqRe3h4ePJ5oVCgUCg0Ox2Rvqst3Wsn8JVL+uoFz3bUC/4DA6HDctWqEMQ3bAjr581r\nHMyh+W31qhtPh2KxSLFYbG3jRk33Vh7Ap4C9wPyqdZuBzT415XJmnX0T/KNEJH5x5Id7kbIp58M1\n1Gw+kUTKxczOBf4W+Ii7/5+q9UuBO4AzgBOBHwLv8po3UspFsiaulEmSJX3lGYnKc0DP1EKX7Enq\nTtGdwLHAmJk9amY3A7j7AeBu4ABwP7BekVvyYPXq6YGxPJFyq5KcCLoczAF27gwPCOs0CfPsoBuL\nUkjjYeRT0h2Qo6Pw8sshj149h+ju3XDssfrs5IUG58qYNFUe6MslPvpdShwU0DMorlu74zqPNHy5\niIgCemalZTyMtHy5SGP1Wv/tjr8i2aDhczMoyc6zdg0MhGC+eHH4qWCePtXD1o6OwsGD4SajsbHK\n0Le7dmko27xTCz2F0pbmUAs9G8r/TpddFkZMHBmBt72tUsaoEsZ8UMolY9LUeZa2Lxdprpym27cP\nvvWt8CX8xS/C7bf3P3Un8VDKJWNarXduNnNNXLPaaMb17KhO033rW6GlvnhxeC0NqTtJngJ6hjWb\n7qvRay+/3F6gj+NmGkle9V9OixbB0BBcfDFcdFFItZS/iKs/E5JDjcYESPqBxnKJRbOxQeq9pvGq\n86l62Nryv+mtt1aG063+99dQttmGpqDLt2bljfVeUydnvqWpD0bipxx6jjUrb2z0msoQ801pstlL\nAT3DavOm1TnSmV5LS427iMRHKZcMa/anNTSe1WbPHpUhimSV6tBTpN/5zX6/v4h0Rzn0FGlWagjx\n1Y+LyOyjgN5j1fXAExPT0x0zBfxuJX18Eemfbqaguw64AHgD+D3wKXf/XfTaVcBngNeBje7+YJ39\nZ2XKpaxZqWHSZYUqWxTJrkRy6Gb2x+7+h+j5BmCpu/+nqjlF/4zKnKKnuPsbNfvP2oA+U0AdHYWF\nC+G00yoBP+48d1qG5hWR9iSSQy8H88ixhJY6wIXAne7+qrtPAE8TJowWmpcTli1bFkbL27cvBPyD\nB+NNi6hsUSSfusqhm9l2M3sW+A/Al6LV7wAOVW12iNBSF2Ye7KpUghtuCEOflgdYOv/8MDbH3r2d\nd5iWO1tryxRXrND4HtI+dd6nU9OUi5mNAcfXeelqd7+varvNwHx3HzazncBP3f070Wu3AD9w9+/X\nHNu3bt06uVwoFCgUCt1cSy5UlxVWD4X63HOVDs1OasjL265YESYRhsq+oLJFaY+GVe6dYrFIsVic\nXL722msbplziGmjrncDj0fPNwOaq1x4AzqyzT9eD1ORZo0G3mg3G1ekxRTqhz1N/kMTgXGb2bnf/\ndfR8A/Bhd7+oqlP0DCqdou/ymjeazZ2iM5mp9dNNh6Y6QyVO+jz1XlI3Fv0XM3vczPYB5wBXALj7\nAeBu4ABwP7Bekbs9zfLs3XRoqjNU4qTPUwo1aron/UApl7Z1M5a5xkGXOOnz1D9oPPR86GYcFo3h\nInHS56l/NDiXiEhOaHCuDGlU3zs8rLpfEWlOAT1lGg2e9elPtz+olm7+EJldlHJJoV27YGwMzj4b\nfvIT2LEjrN+9u7L+kUfqjwFTndcslWDTJhgchDVrdPOHSB4o5ZIxq1bB0aOwbl34+dJLIRCfdVZl\nfb25QGtb92VjY/WH6hWRfFFAT1inaY9582DtWnjlFfjYx8KYLtddF9Y3qvutN9b6jh1wzTWaEFpk\nVmhUz5j0g1lSh95uvW716+Pj7uB+1lnh59q1lf0mJtz/4kMjPvTRlb71Ix/xLStX+sMjI+5e2W98\nXLdni+QNTerQ5/b5+yT3qlvNrUwoUb5LFMK2+/bBJz8J731vZZtSCb6wfpQP/fMVfPmfnplcv+WZ\nZ3j5ZRjds5rxcdi2LazfsWPqeSjtIpJP6hTtkXbGvCh3Xg4NhaF0h4ZCuuWVV8A9pF3+1fgqvvzj\naRNB8e8XruLWxx5gYKDSuVoO6OVj6+YPkexSp2iftTvmRbmVvn9/+HnyySEoFwpw++2hyuUtrx+t\nu++p7zwyGbzXrAn7lcdahxDYFcxF8kkBPWGtzFBUa/XqSuCtTo088kj4UnjkEfh/c+bV3deOnT9l\nudsArlp2kexQQE/YTDMUtaLel8LB+RvZvHjJlO2uXrKEwQ0b4jt5Gt/oFNd0eCISH+XQM6DRQEi3\nfG2Ul366kzlHjvD6/PkMbtjAigTyKTNNai0ivaPBuVKo2Wh1kL6R7DSRgUg6qFM0hZqlMnqZ5mgl\nR66JDEQyolGBeqsP4G+AN4DjqtZdBfwaeBJY2WC/pOruM6PZTT+9uiFophufNJGBSLrQ5MaiboP5\nQsIk0OPlgA4sBX4JvBlYBDwNvKnOvj26/P4aGZke/A4fDuvdp97VWavZa3Fq9uUx0/mLSG81C+jd\nply+CgzVrLsQuNPdX3X3iSign9Hl+2RWOX2ya1dIVVSnTw4ehM99rn4qo5dpjoGB0OFZb7yX2tLJ\n8vaqZRdJn44DupldCBxy98dqXnoHcKhq+RBwYqfvk3XlMsWxMdiwIQxnu2IFPPssnH8+3HRT6GQc\nGgqDcVUH/XZq11vRKF++a1fly+Pyy8MXTfU+Bw9Oz6mrDl0kfZqO5WJmY8DxdV7aQsiTr6zevMmh\n6pazDA8PTz4vFAoUCoVmp5NZAwOVEQ/XroVTT4ULLoB77w13gZZK4Rb/m26qVLk0ql3vpmVc/muh\nfOzyeOlQGR7g5pvDF83ISDi3ZcsqyzD1y0ZEklcsFikWi61t3CgX0+wBLANeIOTOx4FXgQlgAbAZ\n2Fy17QPAmXWO0ZN8UxpU56jXrg158X37+jMKYm2+/K67pr/3xIT76tWVbSYmNGKjSFrQJIceSx26\nmY0DH3T3F81sKXAHIW9+IvBD4F1e80azpQ69tkW7aVOYpGLePNi4EU47Lb7a7lZnYm+lprx2G9Wh\ni6RDL+rQJyOzux8A7gYOAPcD62dF5G6gejjc8oQTO3eGwHvJJWF43Lg6PVupX2+ls7V2m4MHVYcu\nkgmNmu5JP5hFKRf3qeV/5bTHxERlfVypjFZq25vVlNeum5hwX7Ys/Gy0j4j0DkmnXDoxW1Iu9bSa\nGulUo/RIK+9bu83oaOgY3b+/sk2/hyEQmc00lsssooG0RPJNY7nMEknVr4tINqiFniNJp3JEpP+U\ncskpBXCR2Ucpl5zSbEIiUk0t9IxTJ6jI7KKUS87pLk6R2UMplxzTbEIiUqaAniKtTAdX+5rKFEWk\nTAE9Rdrt5CyPE1NvmF0RmX2UQ08ZdXKKSDPqFM0YdXKKSCPqFM0QdXKKSKcU0FOknG5ZsaKSDy/n\n1DWPp4jMpJtJoofN7JCZPRo9Plb12lVm9msze9LMVjY7jlSUOzlXrQqBHMLy7t26A1REZtZNC92B\nr7r7+6PH/QDRFHRrgKXAucDNZpboXwItT6Cacm99a5GBgemt8z17stM5mpd/izxcRx6uAXQd7eg2\n0NZLzF8I3Onur7r7BPA0YX7RxOTxH3xgIFS6LF4cfmYhmEM+/y2yKg/XALqOdnQb0DeY2T4zu9XM\nyiHnHcChqm0OESaLljaoc1RE2tU0oJvZmJk9XudxAfBNYDFwOvA74G+bHEr1iW3QHaAi0olY6tDN\nbBFwn7u/z8w2A7j79dFrDwBb3f0fa/ZRkBcR6UCjOvS5nR7QzE5w999Fi58AHo+e3wvcYWZfJaRa\n3g38rNUTEhGRznQc0IEvm9nphHTKOHAZgLsfMLO7gQPAa8B63RIqIpK8vt36LyIi8cr8naJm9jdm\n9oaZHVe1LjM3NpnZV8zsiaha6Ptm9raq1zJzHQBmdm50rr82sy/0+3xaYWYLzewhM/uVme03s43R\n+uOiooCnzOzBqiquVDOzOdGNfvdFy5m6DjMbMLPvRv8nDpjZmVm7BgAz+3z0eXrczO4ws3m9uI5M\nB3QzWwgMAger1vX8xqYuPQi8191PA54CroLsXYeZzQG+QTjXpcBfmdmp/T2rlrwKfN7d3wucBVwe\nnfdmYMzdTwF+FC1nwRWEdGf5T++sXcfXgR+4+6nAvwGeJGPXYGYnAhuAD7r7+4A5wMX04DpSGyBa\n9FVgqGZdz29s6oa7j7n7G9HiPwInRc8zdR2Ec3va3Sfc/VXgLsI1pJq7P+/uv4yevww8QejMvwC4\nLdrsNuDj/TnD1pnZScB5wC1UbvrLzHVEf51+2N3/AcDdX3P3l8jQNVSZC/yRmc0F/gj4LT24jswG\ndDO7EDjk7o/VvJTlG5s+A/wgep616zgReK5qOe3nO01Ufvt+whfrAnd/IXrpBWBBn06rHX8HXAm8\nUbUuS9exGPhnM/tvZvYLM/uvZvZWsnUNuPtvCPflPEsI5CV3H6MH19FNlUvizGwMOL7OS1sIqYnq\nvHKzMsi+9vw2uY6r3b2c69wCvOLudzQ5VJp7sNN8bjMys2OB7wFXuPsfzCofJ3f3tN83YWbnA793\n90fNrFBvmwxcx1zgA8Dn3P3nZvY1atISGbgGzOzthNb4IuAl4L+b2drqbZK6jlQHdHcfrLfezJYR\nvs33Rf+1hOwJAAABr0lEQVTxTgIeMbMzgd8AC6s2Pyla1zeNrqPMzD5F+FP5L6tWp+46ZlB7vguZ\n+hdGapnZmwnB/Nvufk+0+gUzO97dnzezE4Df9+8MW3I2cIGZnQfMB/7EzL5Ntq7jEOGv7p9Hy98l\nNNyez9A1AJwDjLv7/wUws+8D/5YeXEcmUy7uvt/dF7j7YndfTPggfCD6c+Ze4GIzO8bMFtPgxqa0\nMLNzCX8mX+juR6peytR1AP8EvNvMFpnZMYQO3Xv7fE4zstAiuBU44O5fq3rpXuDS6PmlwD21+6aJ\nu1/t7guj/w8XAz92978mQ9fh7s8Dz5nZKdGqc4BfAfeRkWuIHATOMrO3RJ+vcwgd1YlfR6pb6G2Y\n/NMlgzc27QSOAcaivzb+l7uvz9p1uPtrZvY5YDehV/9Wd3+iz6fViuXAWuAxM3s0WncVcD1wt5mt\nAyaAi/pzeh0rf1aydh0bgO9EjYJngE8TPk+ZuQZ3/5mZfRf4BeH/7i+Avwf+mISvQzcWiYjkRCZT\nLiIiMp0CuohITiigi4jkhAK6iEhOKKCLiOSEArqISE4ooIuI5IQCuohITvx/Ete3744APtoAAAAA\nSUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x106567150>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot the means\n",
    "plt.plot(x, y, 'x', markerfacecolor='blue')\n",
    "plt.plot(mu_x, mu_y, 'x', marker='o', markerfacecolor='red')\n",
    "plt.axis('equal')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Implement the Gibbs sampler"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Forthcoming..."
   ]
  }
 ],
 "metadata": {
  "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}