Emcee + Edward.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.stats import cauchy, gamma\n",
    "\n",
    "import edward as ed\n",
    "import numpy as np\n",
    "import tensorflow as tf\n",
    "from edward.models import Normal, Gamma, StudentT, Empirical\n",
    "from edward.models import NormalWithSoftplusScale, GammaWithSoftplusConcentrationRate\n",
    "\n",
    "import itertools\n",
    "import time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We generate 1 million samples from the Gaussian distribution with the model parameters:\n",
    "\n",
    "1. <b>delta</b> ~ cauchy(0, 1)\n",
    "2. <b>mu</b> ~ cauchy(0, 1)\n",
    "3. <b>sigma</b> ~ gamma(2, 2)\n",
    "4. <b>alpha</b> = delta * sigma\n",
    "5. <b>x</b> ~ normal(mu+alpha, sigma)\n",
    "6. <b>y</b> ~ normal(mu, sigma)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "num_points = 1000000 # [100, 10000, 1000000]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def generate_data_normal(N):\n",
    "    delta = 1.5\n",
    "    mu = 0.0\n",
    "    sigma = 1.0\n",
    "    alpha = delta * sigma\n",
    "    \n",
    "    xdata = np.random.normal(mu+alpha, sigma, N)\n",
    "    ydata = np.random.normal(mu, sigma, N)\n",
    "    return (xdata, ydata)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def plot_points(x, y):\n",
    "    plt.figure(figsize=(10,5))\n",
    "    plt.plot(x, y, 'ok', markersize=3)\n",
    "    plt.grid()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def plot_data_hist(x, y):\n",
    "    plt.figure(figsize=(10,5))\n",
    "    plt.hist(x, alpha=0.5, bins=500, label=\"X data\")\n",
    "    plt.hist(y, alpha=0.5, bins=500, label=\"Y data\")\n",
    "    plt.legend(loc='lower right', fontsize=11)\n",
    "    plt.show()\n",
    "    \n",
    "def plot_trace_hist(trace):\n",
    "    plt.hist(trace, alpha=0.5, bins=500, label=\"delta\")\n",
    "    plt.legend(loc='lower right', fontsize=11)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlkAAAEuCAYAAAC9NwejAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XtQ3PW9//HXsgsoy2KCs71m0KDZk7EOYyBiMqdE06rU\n2h45GU0CR9RqHcOJF9KeSLzk4kki8rPBaHrw1liVuhBsNKdT55xWc5ygJmKk2mhGrILHNsbqivTI\nbmAhyff3h4WEcFku+93vd3efj5lOw3c/7H4+H8jXV96fz37WYRiGIQAAAERVitUdAAAASESELAAA\nABMQsgAAAExAyAIAADABIQsAAMAEhCwAAAATuCI16O/v16pVq/TRRx8pJSVF69evl8vl0qpVq+Rw\nODRr1iytXbtWKSkpampqUmNjo1wulyoqKrRw4UL19vZq5cqV6uzslNvtVk1NjbKzs8d8zUCgO2oD\ntJvp0zPU1XXI6m4kHebdGsx77DHn1mDerWGXefd6PSNej1jJ2rVrlw4fPqzGxkYtX75cmzdvVnV1\ntSorK+X3+2UYhnbu3KlAIKD6+no1NjZq69atqq2tVV9fnxoaGuTz+eT3+1VSUqK6urqoDy6euFxO\nq7uQlJh3azDvscecW4N5t4bd5z1iyJo5c6aOHDmio0ePKhgMyuVyaf/+/SosLJQkLViwQLt379a+\nffs0Z84cpaWlyePxKCcnR21tbWptbVVRUdFg2z179pg7IgAAABuIuFyYkZGhjz76SJdccom6urr0\n0EMPae/evXI4HJIkt9ut7u5uBYNBeTzHymVut1vBYHDI9YG2kUyfnmH7dDoVo5UVYS7m3RrMe+wx\n59Zg3q1h53mPGLIef/xxffvb39ZPf/pTffzxx7r66qvV398/+HgoFFJWVpYyMzMVCoWGXPd4PEOu\nD7SNxA7rq2bxej0JvefMrph3azDvscecW4N5t4Zd5n3Se7KysrIGK1GnnHKKDh8+rLPOOkstLS2S\npObmZs2dO1d5eXlqbW1VOBxWd3e32tvb5fP5lJ+fr127dg22LSgoiNaYAAAAbMsR6QOiQ6GQbr/9\ndgUCAfX39+uqq67S2WefrdWrV6u/v1+5ubnasGGDnE6nmpqatG3bNhmGoRtuuEHFxcXq6elRVVWV\nAoGAUlNTtWnTJnm93jE7ZYdUaha7pO5kw7xbg3mPPebcGsy7Newy76NVsiKGLCvYYcLMYpdfiGTD\nvFuDeY895twazLs17DLvk14uBAAAwMQRsgAAAExAyAIAADABIQsAAMAEhCwAAAATELIARN1zHb+3\nugsAYDlCFgAAgAkIWQBihgrX2Ha81GF1FwBEESELgCkIVEMdH6DGE6YIXED8I2QBQAwRnoDk4bK6\nAwASH1WtL41UzSopyh32GIDEQCULQNScGKaa3v5txDaJZiAsnfj/AJIPIQtAVEw0PCV62JLYewUk\nO0IWgAk5MRwlQ1gy22hBiwAGxDdCFgBEQbQCEcEKSByELABRNdllQypixxC0gMRAyAIwaeMJRmO1\nGemx5zp+H3eB68RQFOlrAMmBkAVgwqa6LytR93URpgAcj3OyAIxLLIJRooQtAJCoZAHAlMSiesWZ\nW0B8ImQBwAQRegCMByELQEQs4w0Xq4BFkAPiFyELwIgIVvZE6ALiByELAOIA4QqIP4QsAGOysqJl\np2oaIQfARBGyAAyyU6ixU1/shsAHxIeI52Q988wzevbZZyVJ4XBY77zzjvx+v+6++245HA7NmjVL\na9euVUpKipqamtTY2CiXy6WKigotXLhQvb29WrlypTo7O+V2u1VTU6Ps7GzTBwYA0bbjpQ6VFOVa\n3Q0AcSJiJWvRokWqr69XfX29vvWtb+nOO+/Uf/zHf6iyslJ+v1+GYWjnzp0KBAKqr69XY2Ojtm7d\nqtraWvX19amhoUE+n09+v18lJSWqq6uLxbgATNLxFSQ7VZPs1Bc72PFSBxUtwObGvVz41ltv6f33\n39eSJUu0f/9+FRYWSpIWLFig3bt3a9++fZozZ47S0tLk8XiUk5OjtrY2tba2qqioaLDtnj17zBkJ\nAACAjYz7Y3UefvhhLV++XJJkGIYcDockye12q7u7W8FgUB6PZ7C92+1WMBgccn2gbSTTp2fI5XJO\naCDxxOv1RG6EqGPeI8v4JF2S9OInu5ThTo/Oc07xeQZ+bhmfpFv6M3T/fRxer2fwz3Yw0pzwu24N\n5t0adp73cYWsL774Qh988IHmzZsnSUpJOVYAC4VCysrKUmZmpkKh0JDrHo9nyPWBtpF0dR2a0CDi\nidfrUSAQOWgiupj38TkUCkf1+TLc6VN+zsdbtuvS3It1KBS29GcY+vs4AoHuwT/bwYlzwu+6NZh3\na9hl3kcLeuNaLty7d6/mz58/+PVZZ52llpYWSVJzc7Pmzp2rvLw8tba2KhwOq7u7W+3t7fL5fMrP\nz9euXbsG2xYUFEx1LAAQU+x9AjAZ4wpZH3zwgWbMmDH4dVVVlbZs2aIlS5aov79fxcXF8nq9Ki8v\nV1lZma6++mqtWLFC6enpKi0t1XvvvafS0lJt27ZNN954o2mDATB58bCxPB76GGvHB0DCIGAvDsMw\nDKs7cSI7lP7MYpfSZrJh3iMzI8BEY7nwRJfmXhzV54skHoLLwLESO17qkNudrovyv2lxj5IP9xhr\n2GXep7RcCAAAgIkhZAFJLt6W4OKtvwCSFyELQNyJVdCKh6VCiYNJAbsiZAEAAJiAkAUgLpfg4rHP\nAJILIQsAAMAEhCwgiVENAgDzELIAAABMMO4PiAaAZMC79ABEC5UsAHGL5c7RERYB6xGygCRFQAEA\ncxGyAODvqP4AiCZCFoC4RkUOgF0RsoAkQyhJfFTkAHsgZAGAEjeYJOq4gHhAyAIAADABIQtIEscv\nE7JkOFSiVnsSdVxAvCBkAUCCI2wB1iBkAUgIA9W58VbpBoIHAQSAWQhZQIJLhmXCRB0XgPhGyAIA\nADABIQtA0mKpEICZCFlAEmA5DQBij5AFAEmAqh0Qe4QsAAmDih0AO3GNp9HDDz+s//mf/1F/f79K\nS0tVWFioVatWyeFwaNasWVq7dq1SUlLU1NSkxsZGuVwuVVRUaOHChert7dXKlSvV2dkpt9utmpoa\nZWdnmz0uAAAAS0WsZLW0tOiNN95QQ0OD6uvr9de//lXV1dWqrKyU3++XYRjauXOnAoGA6uvr1djY\nqK1bt6q2tlZ9fX1qaGiQz+eT3+9XSUmJ6urqYjEuAEmOqhYAq0UMWS+//LJ8Pp+WL1+uZcuW6YIL\nLtD+/ftVWFgoSVqwYIF2796tffv2ac6cOUpLS5PH41FOTo7a2trU2tqqoqKiwbZ79uwxd0QAkh4B\na2TsywJiK+JyYVdXlw4ePKiHHnpIBw4cUEVFhQzDkMPhkCS53W51d3crGAzK4/EMfp/b7VYwGBxy\nfaAtAABAoosYsqZNm6bc3FylpaUpNzdX6enp+utf/zr4eCgUUlZWljIzMxUKhYZc93g8Q64PtI1k\n+vQMuVzOyYwnLni9nsiNEHXJOu8Zn6Rb+/pu615/pJ+5/3dtKiueLbeF/TLbWGN7/g8fSZLKimfH\nqjtJI1nvMVaz87xHDFkFBQV68skn9aMf/Uiffvqpenp6NH/+fLW0tOi8885Tc3Oz5s2bp7y8PG3e\nvFnhcFh9fX1qb2+Xz+dTfn6+du3apby8PDU3N6ugoCBip7q6DkVlcHbk9XoUCFDNi7VknvdDobBl\nr53hTrf09Uf6mYdCYT36zB8t6E1suN3pCo1jzpP174NZkvkeYyW7zPtoQS9iyFq4cKH27t2ryy+/\nXIZhaM2aNZoxY4ZWr16t2tpa5ebmqri4WE6nU+Xl5SorK5NhGFqxYoXS09NVWlqqqqoqlZaWKjU1\nVZs2bYr64AAMNbAn6dLciy3uCQAkL4dhGIbVnTiRHVKpWeySupNNss27XTZ+W13JGilkJvrm7/FW\nskqKcmPQm+SRbPcYu7DLvE+6kgUAiSDRwxUA++HEdwAAABMQsgAkrOc6fm+bpVM7oroHmIuQBQAA\nYAJCFoCER8VmOOYEMB8hCwAAwASELAAJr/3IXqu7YGtUtQBzELIAIIkRsADzELKABMO76Y5p+3OX\n1V0AkMQIWUACIFiNbiBosWQ4NipaQPQRsgAkJKpYAKxGyAISCBWt8aGqBSAWCFkAAAAmIGQBSBpU\nsADEEiELSBAsFQKAvRCyAAAATEDIAuIY1auJY8kQQKwQsgAAgzgvC4geQhYAAIAJXFZ3AMDUsGQ4\nFIeQArALKlkAkhJ7s4ZjqRCILkIWAACACQhZQJximRAA7I2QBQAAYAJCFgBgiB0vdbA/C4gCQhaA\nhND2565JvbOQDfAAzDKuIxz++Z//WZmZmZKkGTNmaNmyZVq1apUcDodmzZqltWvXKiUlRU1NTWps\nbJTL5VJFRYUWLlyo3t5erVy5Up2dnXK73aqpqVF2drapgwISHfuxAMD+IlaywuGwDMNQfX296uvr\nVV1drerqalVWVsrv98swDO3cuVOBQED19fVqbGzU1q1bVVtbq76+PjU0NMjn88nv96ukpER1dXWx\nGBcAREQVa2wsGQJTEzFktbW1qaenR9dee62uuuoqvfnmm9q/f78KCwslSQsWLNDu3bu1b98+zZkz\nR2lpafJ4PMrJyVFbW5taW1tVVFQ02HbPnj3mjggAAMAGIi4XnnTSSbruuut0xRVX6H//9391/fXX\nyzAMORwOSZLb7VZ3d7eCwaA8Hs/g97ndbgWDwSHXB9pGMn16hlwu52THZHterydyI0Rdosx709u/\nVYY73epujFus+pqaOvkPsHCfHD/zOR7uKM7583/4SGXFs6P2fIksUe4x8cbO8x7xrjRz5kyddtpp\ncjgcmjlzpqZNm6b9+/cPPh4KhZSVlaXMzEyFQqEh1z0ez5DrA20j6eo6NJmxxAWv16NAIHLQRHQl\n0rwfCoWt7sK4ZbjTY9bf/v7Dk/7eff0v6wznuVHsjXXc7nSFojznifJ3x0yJdI+JJ3aZ99GCXsTl\nwl//+te65557JEmffPKJgsGg/vEf/1EtLS2SpObmZs2dO1d5eXlqbW1VOBxWd3e32tvb5fP5lJ+f\nr127dg22LSgoiNaYAAAAbCtiJevyyy/XbbfdptLSUjkcDt19992aPn26Vq9erdraWuXm5qq4uFhO\np1Pl5eUqKyuTYRhasWKF0tPTVVpaqqqqKpWWlio1NVWbNm2KxbgAJJFofCh0+5G9CVPNMsOOlzpU\nUpRrdTeAuOIwDMOwuhMnskPpzyx2KW0mm0SZ93g7uiFWy4XRCFmSEiJkmbFcWFKUS8iKIFHuMfHG\nLvM+6eVCAAAATBwhC4gD8VbBQmLhvCxgcghZAAAAJiBkAXGCatbIorUfCwCijZAF2BzhamSEKwB2\nR8gCgL/jswwBRBMhC0DcMqOaRdACEC2ELAAAABMQsgDgBFSzAEQDIQsAAMAEhCzAxnhnIQDEL0IW\nAIyAJUMAU0XIAgCMGx+xA4wfIQtA3OEgUgDxwGV1BwAMx16skbX9uUuzc6Zb3Y2kRRULmBgqWQAw\nCvZljYywBYwPIQtAXGGpEEC8IGQBAACYgJAF2Ah7sQAgcRCyAGAM7MsaGfuygMgIWYDNUM0CgMRA\nyAIQF9jwbj9Us4CxEbIAIIKBJUOWDocjaAGjI2QBsD2qWADiESELgK0RsADEq3GFrM7OTp1//vlq\nb2/Xhx9+qNLSUpWVlWnt2rU6evSoJKmpqUmLFi3S4sWL9eKLL0qSent7ddNNN6msrEzXX3+9Pv/8\nc/NGAsQ5NrzbG0uFo2PJEBhZxJDV39+vNWvW6KSTTpIkVVdXq7KyUn6/X4ZhaOfOnQoEAqqvr1dj\nY6O2bt2q2tpa9fX1qaGhQT6fT36/XyUlJaqrqzN9QAAAAHYQMWTV1NRo6dKl+spXviJJ2r9/vwoL\nCyVJCxYs0O7du7Vv3z7NmTNHaWlp8ng8ysnJUVtbm1pbW1VUVDTYds+ePSYOBYhfVLEAIPGMGbKe\neeYZZWdnDwYlSTIMQw6HQ5LkdrvV3d2tYDAoj8cz2MbtdisYDA65PtAWAJB4WDIEhnON9eD27dvl\ncDi0Z88evfPOO6qqqhqyryoUCikrK0uZmZkKhUJDrns8niHXB9qOx/TpGXK5nJMZT1zwej2RGyHq\n7DrvTW//VhnudKu7YZqpjO2t9z9TauqYt6mYc59s/5+V26LfJ7v+HYuVZB+/Vew872PevZ566qnB\nP5eXl2vdunW699571dLSovPOO0/Nzc2aN2+e8vLytHnzZoXDYfX19am9vV0+n0/5+fnatWuX8vLy\n1NzcrIKCgnF1qqvr0NRGZWNer0eBABW9WLPzvB8Kha3ugmky3OlTGl9//+Eo9iY6Qkft/fNyu9MV\nsuh3yq5/x2LBzveYRGaXeR8t6E34n4hVVVVavXq1amtrlZubq+LiYjmdTpWXl6usrEyGYWjFihVK\nT09XaWmpqqqqVFpaqtTUVG3atGnKAwEAAIgHDsMwDKs7cSI7pFKz2CV1Jxs7z3sib3qfaiXLrmdk\nneE81+oujMrKSpYklRTlWvbaVrLzPSaR2WXeR6tkcRgpAACACQhZAGzJrlUsABgvQhYA27F7wOL0\ndwDjQcgCYCt2D1gDCFoj47ws4BhCFgAAgAkIWYCFEvmdhQCQ7AhZAGwjXpYKAWA8CFkAMEnsywIw\nFkIWAACACQhZgEXYj4VExTsMgS8RsgALELAAIPERsgAAAExAyAKAKWDz++hYNkSyI2QBAKKOgAUQ\nsgDYRDyfkUU1a3SELSQzQhYQY2x6B4DkQMgCYLl4rmINoJoF4ESELCCGqGINlwgBCwBGQsgCAAAw\nASELiBGqWEhWbH5HsiJkATFAwAKA5EPIAgCYjmoWkhEhCwCiZOAdhrzTEIBEyAJMx1JhciFgARhA\nyAJgGY5vAJDICFkAgJhgXxaSjStSgyNHjujOO+/UBx98IIfDobvuukvp6elatWqVHA6HZs2apbVr\n1yolJUVNTU1qbGyUy+VSRUWFFi5cqN7eXq1cuVKdnZ1yu92qqalRdnZ2LMYGwMYSvYrVfmSvznCe\na3U3AFgoYiXrxRdflCQ1NjaqsrJS9913n6qrq1VZWSm/3y/DMLRz504FAgHV19ersbFRW7duVW1t\nrfr6+tTQ0CCfzye/36+SkhLV1dWZPigAAACrRaxkXXjhhbrgggskSQcPHlRWVpZ2796twsJCSdKC\nBQv0yiuvKCUlRXPmzFFaWprS0tKUk5OjtrY2tba26sc//vFgW0IWgESvYgGANI6QJUkul0tVVVV6\n/vnn9cADD+iVV16Rw+GQJLndbnV3dysYDMrj8Qx+j9vtVjAYHHJ9oG0k06dnyOVyTmY8ccHr9URu\nhKizat4zPkm35HXtIsM9fPypqeO69cQ998nW/OzdI8y5XSTy/S+Rx2Zndp73cd/pampq9G//9m9a\nvHixwuHw4PVQKKSsrCxlZmYqFAoNue7xeIZcH2gbSVfXoYmMIa54vR4FApGDJqLLinnn6IYvA9ah\nUHjItWSqYoWOhiM3ijK3O12hUOxfd7wefeaPkqSSolyLexJd3NutYZd5Hy3oRdyTtWPHDj388MOS\npJNPPlkOh0Nnn322WlpaJEnNzc2aO3eu8vLy1NraqnA4rO7ubrW3t8vn8yk/P1+7du0abFtQUBCt\nMQEAANhWxErWxRdfrNtuu03/8i//osOHD+v222/XGWecodWrV6u2tla5ubkqLi6W0+lUeXm5ysrK\nZBiGVqxYofT0dJWWlqqqqkqlpaVKTU3Vpk2bYjEuAAAASzkMwzCs7sSJ7FD6M4tdSpvJhuVCa5y4\nXJhMS4UDYn2Mg92XCweUFOVqx0sdCbNsyL3dGnaZ90kvFwKYOAIWAICQBQAm4rMMgeRFyAIAADAB\nIQuIMpYKgcj4HEMkA0IWgJhIxk3vAJIbIQuA6QhY7M0CkhEhC4gilgoxEgLW6Fg2RCIjZAEwFVUs\nRELQQqIiZAFRQhULAHA8QhYQBQQsYGqoZiEREbIAIEbYmwUkF0IWAMQQQQtIHoQsAKZ56/3PrO4C\nAFiGkAUAsA32ZiGRELKAKWLT+8g4umF0LBmOjICFREPIAgAAMAEhCwAsQDULSHyELGCSWCYcHUuF\nAEDIAqaEoAUAGA0hCwAAwASELAAAABMQsoBJYJlwdOzHGj82v4+MoxyQKAhZAADbIWghERCygAmi\nijU6qlgAcAwhC0BUELAAYCjXWA/29/fr9ttv10cffaS+vj5VVFTozDPP1KpVq+RwODRr1iytXbtW\nKSkpampqUmNjo1wulyoqKrRw4UL19vZq5cqV6uzslNvtVk1NjbKzs2M1NgAAAMuMWcn6zW9+o2nT\npsnv9+sXv/iF1q9fr+rqalVWVsrv98swDO3cuVOBQED19fVqbGzU1q1bVVtbq76+PjU0NMjn88nv\n96ukpER1dXWxGhdgCpYKh6OCNTVsfh/bjpc62J+FuDVmyPre976nW265RZJkGIacTqf279+vwsJC\nSdKCBQu0e/du7du3T3PmzFFaWpo8Ho9ycnLU1tam1tZWFRUVDbbds2ePycMBYAWC1tQRtoYjXCHe\njRmy3G63MjMzFQwGdfPNN6uyslKGYcjhcAw+3t3drWAwKI/HM+T7gsHgkOsDbQEAQxGwgMQ05p4s\nSfr444+1fPlylZWV6Yc//KHuvffewcdCoZCysrKUmZmpUCg05LrH4xlyfaDteEyfniGXyznRscQN\nr9cTuRGibqrz3vT2b5XhTo9SbxLDW+9/ptTUsW8jkR7HMe6To/P75U7A39N4uG/GQx8TkZ3nfcy7\n32effaZrr71Wa9as0fz58yVJZ511llpaWnTeeeepublZ8+bNU15enjZv3qxwOKy+vj61t7fL5/Mp\nPz9fu3btUl5enpqbm1VQUDCuTnV1HZr6yGzK6/UoEKCiF2tTnXf2Yo2sv//wmI+nproitsEx+/pf\n1hnOc6f0HG53ukKhcJR6ZB92v29yb7eGXeZ9tKDnMAzDGO2bNmzYoP/6r/9Sbm7u4LU77rhDGzZs\nUH9/v3Jzc7VhwwY5nU41NTVp27ZtMgxDN9xwg4qLi9XT06OqqioFAgGlpqZq06ZN8nq9ETtrhwkz\ni11+IZINIcsckfZiEbImZypBK1FDliSVFOVGbmQR7u3WsMu8TypkWcUOE2YWu/xCJJupzDsBa3SE\nLHMQskZGyMKJ7DLvo4UsDiMFMCm8oxCxxrsNEW8IWQBgM+1H9vKOQyABELIAAABMQMgCRsFeLADA\nVBCygDEQtEbGfixYhX1ZiCeELACwKfZljYyghXhByAIwblSwAGD8CFkAxmUgYBG0YotqFhC/CFnA\nCNiLBdgbS4aIB4Qs4AQELNgN1SwgPhGygOMQsIZjeRB2teOlDipasDVCFvB3BKzREbSsRzULiD+E\nLABAXKOaBbsiZAGiigXEO4IW7IiQBQBxgiVDIL4QspD0qGIhnhC0gPhByAIwKja82xNBC4gPhCwA\nIyJg2RtBC7A/QhaSGkuFQOJg8zvshpCFpEXAAhIPQQt2QshCUiJgja7tz10sFSKuEbRgF4QsAIhT\n7MsC7I2QBWAQFSwkCqpZsANCFpIKy4RINFSzAPsiZCHpELRGRhUrfrUf2UvYGgHVLFiNkAWAgIWE\nRtiCVcYVsv74xz+qvLxckvThhx+qtLRUZWVlWrt2rY4ePSpJampq0qJFi7R48WK9+OKLkqTe3l7d\ndNNNKisr0/XXX6/PP//cpGEAkVHBQqKjmjUcAQtWihiyHn30Ud15550Kh8OSpOrqalVWVsrv98sw\nDO3cuVOBQED19fVqbGzU1q1bVVtbq76+PjU0NMjn88nv96ukpER1dXWmDwgYSdPbv7W6C7ZFFQvJ\ngLAFK7giNcjJydGWLVt06623SpL279+vwsJCSdKCBQv0yiuvKCUlRXPmzFFaWprS0tKUk5OjtrY2\ntba26sc//vFgW0IWYC8ErMTzTs+r6j9yWGc4z7W6K0gQZgXUkqLciG3+8IfXddddd+jxxxs0fXq2\nJMnvr9f+/fu0ceO9Y37fY489op///JFR27z8crMOHPizli69cuKdH6eIIau4uFgHDhwY/NowDDkc\nDkmS2+1Wd3e3gsGgPB7PYBu3261gMDjk+kDb8Zg+PUMul3NCA4knXq8nciNEzUAVK8OdbnFP7Cc1\nNeItIC5eA0OlprrkPpnf9xN5vR75f9emsuLZpj1/InKbdO8cz3wVFy/Um29epvvuu0cPPvig3nzz\nTT333A79+te/VlaWZ9TnmTYtQ2lprjFf48CBjnH3Y7ImfPdLSTm2whgKhZSVlaXMzEyFQqEh1z0e\nz5DrA23Ho6vr0ES7FTe8Xo8CgfGFTUzdwD6sDHe6DoXCFvfGXmJRxUpNdam//7Dpr4NjBuY8dJTf\n9xM9+swfJcmUe3Ai39tDJt07xztf5eXX6/rrr1Zd3aPavn2b7rzzLoXDDgUC3UPm/bXXXtUDD9Qq\nLS1Np512uvr6DisQ6NYbb7TqkUfqFA73qru7WxUVN2vmzFz5/Q2SJI8nW4WF81RdvV7BYLc6Oz/T\nhRcWq6LipnGPZbSgNuF3F5511llqaWmRJDU3N2vu3LnKy8tTa2urwuGwuru71d7eLp/Pp/z8fO3a\ntWuwbUFBwURfDoAJWCYEEC9SU1O1Zs16/fzn9+m7371YZ5+dN6xNX1+fNm5cqw0bavTYY79Sevqx\n6tv27du0atVqPfbYU1q1arUef/xRzZyZq8suW6TLLlukSy/9Jz3//O900UXFeuSRx/XEE4169tlf\n629/+9uU+z7hkFVVVaUtW7ZoyZIl6u/vV3Fxsbxer8rLy1VWVqarr75aK1asUHp6ukpLS/Xee++p\ntLRU27Zt04033jjlDgOYPD6XMHnwTkMkkrfe+qNOOWWaXn/9NR0+PLw63tHxvk491avTT58pSbrk\nkh8MPrZ69Xp1dLyvxx//hRobf6Wenp5h319WVq6vfvVr8vvrdf/9P9Phw/3q7R3ebqLGtVw4Y8YM\nNTU1SZI8DdwLAAAK/UlEQVRmzpypX/3qV8PaLF68WIsXLx5y7eSTT9YDDzww5U4Ck8GRDUh27Uf2\nsgF+BAMbucez8RrW++CDDj322MN68MGtqq7+dz3xxFZdd90NJ7RyyDCODn7ldB7b1718+fXKzy/Q\nnDkFKig4V3fddeew19iy5T4dPPiRLrroe1qw4AK9/vprMgxjyn3nMFIkJALWUFSwkhcVLcSzcDis\ntWtv07/+6y365jdn6M4779L27U16++23hrQ788xZ6urq0nvv/UmS9MILv5MkffHF/+kvf/lQ1123\nTPPnf1uvvfbq4PmeTqdTR44ckSS9/nqLysrK9Z3vXKhPP/1EgcCng+2mgrf9AECCo6KFqbCy4rdl\nS61yc89UcfH3JUlf+9rXdfPNP9H69av1y1/6JX254dzlcmnduo3asGGNnE6nfL4v30GalXWKfvCD\nEpWXL5bb7da3vpWn3t5e9fT06Jxz8rVx4zplZ2fryiuv0fr1a5SZ6VF2drZmzz5LBw9+pG9+c8aU\n+u8wolEPi7JEfYeGlNjvQLGD0SpYyfzuQisrWLy7MPZGm3NC1uiiESK4t1vDLvMetXcXAgDiD8uG\nY+NEeJiB5UIkBPZgjYx9WDgRS4dA7FDJAoAkMVDNoqo11EAVi2oWoo1KFuIaFayRUcFCJFS0APNR\nyQISCOEKE0FFa7gdL3VQ0ULUELIQt6hiAQDsjOVCxCUC1uioZmGiWDrEWMy6316ae3HENrW1Nfr8\n805t2PD/Bq+99tqruvfeaj3xxLFzsk70hz+8rscee0Q///kjoz73yy8368CBP2vp0isn3PfxopKF\nuEPAGmrgNHfCFSaDJcORsWRoD8uW3aR3323Tyy83S5J6enr0s59V67bbVisjwz2l53733XcUCoWi\n0c1RUcmC7T3X8Xtdmnsx4QowEdWs4Xa81KGSolw+69BCGRkZuvXWO1Rd/e8qKDhXv/jFQ/r2txco\nP3/usLavvfaqHnigVmlpaTrttNMHr7/xRqseeaRO4XCvuru7VVFxs2bOzNV//uczkr48Rb6wcJ6q\nq9crGOxWZ+dnuvDCYlVU3DTl/hOyEBcIWCOjeoVoImgNd3xFayB0IbbOPfc8nXfefN1991368MMP\n9OijTwxr09fXp40b1+r++x/S6afP1D33rB98bPv2bVq1arVOO+10tbbu1f33/0xPPrlNl122SJJ0\n6aX/JL+/XhddVKxLLvmBgsGgFi26VKWl5Zo2bdqU+s5yIRCHWB6EWVg+hB3deGOl9u59VZWVK5We\nftKwxzs63tepp3p1+ukzJUmXXPKDwcdWr16vjo739fjjv1Bj46/U09Mz7PvLysr11a9+TX5/ve6/\n/2c6fLhfvb3D200UIQu2RgULiL32I3sJW6Ngr5Y13O5MZWZ69PWvf2OUFg4ZxtHBr5xO5+Cfly+/\nXu+8s1//8A+zddVV12qkj2zesuU+Pf10o772ta/r6quv0ymnTBux3UQRsmA7z3X8fvB/GIoKFmKJ\nsDUyToi3nzPPnKWuri69996fJEkvvPA7SdIXX/yf/vKXD3Xddcs0f/639dprr+ro0S/DmNPp1JEj\nRyRJr7/eorKycn3nOxfq008/USDw6WC7qWBPFmyFYDVU25+7NDtn+uCfASsMBC32ax1DwLIXl8ul\ndes2asOGNXI6nfL5ZkuSsrJO0Q9+UKLy8sVyu9361rfy1Nvbq56eHp1zTr42blyn7OxsXXnlNVq/\nfo0yMz3Kzs7W7Nln6eDBj/TNb86YUr8cRjTqYVEWCHRb3QXTeL2ehB7fVJgZsDLc6ToUCpv2/GaJ\n92CVmupSf/9hq7uRVMyec4LWcCVFudzbLWKXefd6Rz6vi+VCWI6lwZHFe8BCYmIJcbjjlw+pcOF4\nhCzAhghYQHzx/67N6i7AhtiThZijajXcwN4rwhXiyYkVLZYSj+EAU0hUshBjBKyhjg9VBCzEk5GW\nDFlG/NKJB5iyhJi8qGTBdASr4QhXSFS8E3FkxwctqlvJg5CFqCFMjY7lQCSb46taZzjP5SN7jjNS\nZYvglZgIWZgywtVQowUpAhaS1fGBi7CFZELIwrg91/F7XZp7MaHqBMdXqQYODgUw3EDYOrHKNXAt\nmcPXaPu2qHDFN9MPIz169KjWrVund999V2lpadqwYYNOO+20Mb/HDgeLmcUuB6dFkmiBKtqHkVKV\nGh8OI429RJvzeAlebne6QjE88Jjw9SW7/Dd1tMNITa9kvfDCC+rr69O2bdv05ptv6p577tGDDz5o\n9stiHCIFqEQJWJNBiALsIdI7FuMlhEXbWO9YJIDZh+khq7W1VUVFRZKkc845R2+//bbZL5lwjq8q\nJVJ1aaIGgs/xG8iP/1y/EzeWH/91ov3rHsCXon1sxPHLl8d/HU+memRESVHukHO+OPNr8kxfLrzj\njjt08cUX6/zzz5ckXXDBBXrhhRfkcrEdDAAAJC7TDyPNzMxUKBQa/Pro0aMELAAAkPBMD1n5+flq\nbm6WJL355pvy+XxmvyQAAIDlYvbuwj/96U8yDEN33323zjjjDDNfEgAAwHKmhywAAIBkxAdEAwAA\nmICQBQAAYAJClkXa29tVUFCgcDh2JwQns+7ubi1btkxXXnmllixZojfeeMPqLiWso0ePas2aNVqy\nZInKy8v14YcfWt2lpNDf36+VK1eqrKxMl19+uXbu3Gl1l5JKZ2enzj//fLW3t1vdlaTx8MMPa8mS\nJVq0aJGefvppq7szIs5SsEAwGFRNTY3S0tKs7krS+OUvf6l58+bpmmuuUUdHh37605/q2Weftbpb\nCYlPebDGb37zG02bNk333nuv/va3v6mkpETf/e53re5WUujv79eaNWt00kknWd2VpNHS0qI33nhD\nDQ0N6unp0WOPPWZ1l0ZEJSvGDMPQ6tWr9ZOf/EQnn3yy1d1JGtdcc42WLl0qSTpy5IjS09Mt7lHi\n4lMerPG9731Pt9xyi6Qv7zNOp9PiHiWPmpoaLV26VF/5yles7krSePnll+Xz+bR8+XItW7ZMF1xw\ngdVdGhGVLBM9/fTTeuKJJ4Zc+8Y3vqHvf//7mj17tkW9Snwjzfvdd9+tvLw8BQIBrVy5UrfffrtF\nvUt8wWBQmZmZg187nU4dPnyYQ4hN5na7JX05/zfffLMqKyst7lFyeOaZZ5Sdna2ioiI98sgjVncn\naXR1dengwYN66KGHdODAAVVUVOi///u/5XA4rO7aENz1THTFFVfoiiuuGHLtoosu0vbt27V9+3YF\nAgFde+21euqppyzqYWIaad4l6d1339VPfvIT3XrrrSosLLSgZ8mBT3mwzscff6zly5errKxMP/zh\nD63uTlLYvn27HA6H9uzZo3feeUdVVVV68MEH5fV6re5aQps2bZpyc3OVlpam3Nxcpaen6/PPP9ep\np55qddeG4M4XY88///zgn7/zne/Ydh050bz//vu65ZZbtHnzZqqIJsvPz9eLL76o73//+3zKQwx9\n9tlnuvbaa7VmzRrNnz/f6u4kjeP/kVxeXq5169YRsGKgoKBATz75pH70ox/p008/VU9Pj6ZNm2Z1\nt4YhZCEpbNq0SX19fdq4caOkL6stbMY2x0UXXaRXXnlFS5cuHfyUB5jvoYce0hdffKG6ujrV1dVJ\nkh599FE2YyMhLVy4UHv37tXll18uwzC0Zs0aW+5D5MR3AAAAE/DuQgAAABMQsgAAAExAyAIAADAB\nIQsAAMAEhCwAAAATELIAAABMQMgCAAAwASELAADABP8f346dS0jdKagAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1313dc198>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "np.random.seed(42)\n",
    "xdata, ydata = generate_data_normal(num_points)\n",
    "plot_data_hist(xdata, ydata)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Further we will investigate sampling package [Emcee](https://github.com/dfm/emcee) and variational inference using [Edward](https://github.com/blei-lab/edward) package."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Emcee"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.2.1\n"
     ]
    }
   ],
   "source": [
    "import emcee\n",
    "print(emcee.__version__)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "_xdata = np.array(xdata, dtype=float)\n",
    "_ydata = np.array(ydata, dtype=float)\n",
    "\n",
    "assert(all(np.isnan(_xdata)) == False)\n",
    "assert(all(np.isnan(_ydata)) == False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def lnprior(theta):\n",
    "    delta, mu, sigma = theta\n",
    "    \n",
    "    prob_delta = cauchy.pdf(delta, 0, 1)\n",
    "    prob_mu = cauchy.pdf(mu, 0, 1)\n",
    "    prob_sigma = gamma.pdf(sigma, 2, 2)\n",
    "    \n",
    "    if prob_sigma <= 0:\n",
    "        return -np.inf\n",
    "    else:\n",
    "        return np.log(prob_delta) + np.log(prob_mu) + np.log(prob_sigma)\n",
    "    \n",
    "    \n",
    "def lnlike(theta, x, y):\n",
    "    delta, mu, sigma = theta\n",
    "    \n",
    "    log_coeff = np.log( 1 / np.sqrt(2 * np.pi * sigma**2) )\n",
    "    denom = 2 * (sigma ** 2)\n",
    "    x_g = ( (x - (mu + delta/sigma))**2 ) / denom\n",
    "    y_g = ( (y - mu)**2 ) / denom\n",
    "                \n",
    "    return (2 * log_coeff) - x_g - y_g \n",
    "\n",
    "def lnposterior(theta, x, y):\n",
    "    return lnprior(theta) + lnlike(theta, x, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ndim = 3  # number of parameters in the model\n",
    "nwalkers = 100  # number of MCMC walkers\n",
    "nburn = 1000  # \"burn-in\" period to let chains stabilize\n",
    "nsteps = 2000  # number of MCMC steps to take\n",
    "\n",
    "# set theta near the maximum likelihood, with \n",
    "starting_guesses = np.random.random((nwalkers, ndim))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/ddedik/anaconda2/envs/py3k/lib/python3.6/site-packages/emcee/ensemble.py:335: RuntimeWarning: invalid value encountered in subtract\n",
      "  lnpdiff = (self.dim - 1.) * np.log(zz) + newlnprob - lnprob0\n",
      "/Users/ddedik/anaconda2/envs/py3k/lib/python3.6/site-packages/emcee/ensemble.py:336: RuntimeWarning: invalid value encountered in greater\n",
      "  accept = (lnpdiff > np.log(self._random.rand(len(lnpdiff))))\n"
     ]
    }
   ],
   "source": [
    "start_time = time.time()\n",
    "\n",
    "sampler = emcee.EnsembleSampler(nwalkers, ndim, lnposterior, args=[xdata, ydata])\n",
    "sampler.run_mcmc(starting_guesses, nsteps)\n",
    "\n",
    "end_time = time.time()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeoAAAFJCAYAAABU5W56AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAH8xJREFUeJzt3X9sVfX9x/HXbW9vLbe3A0JnTFxRGHdMSGN/BJeQVtjE\nGqcLgly8F7tFEJU4XDslRSygqRMb17pBqL9GtqSsvWvEGKbbjCLSb7Ajs3M0sHVL+kWSqXEVSrj3\nIrdiz/ePb7jjXuD+aO9tP708H//Q8zmfc8778+6599V7aW5tlmVZAgAARsqZ6AIAAMDlEdQAABiM\noAYAwGAENQAABiOoAQAwGEENAIDB7BNdwKUMDgbSdq5p06ZoaOhM2s53paKP6UEf04M+jh09TI90\n9bG42HXZfVn/itpuz53oErICfUwP+pge9HHs6GF6jEcfsz6oAQCYzAhqAAAMRlADAGAwghoAAIMR\n1AAAGIygBgDAYAQ1AAAGI6gBADAYQQ0AgMGSCuoTJ07o5ptv1sDAgI4fPy6v1yufz6etW7dqZGRE\nktTV1aVly5bJ4/Fo//79kqSzZ89q/fr18vl8Wrt2rU6ePJm5lQAAkIUSBvWXX36pLVu26KqrrpIk\nbdu2TXV1dero6JBlWdq3b58GBwfV3t4uv9+vXbt2qbW1VcPDw+rs7JTb7VZHR4eWLl2qtra2jC8I\nAIBskjCom5ubdc899+jrX/+6JOno0aNasGCBJKm6ulrvv/+++vr6VFZWJofDIZfLpZKSEvX396u3\nt1dVVVWRuT09PRlcCgAA2SfuX8967bXXNH36dFVVVenll1+WJFmWJZvNJklyOp0KBAIKBoNyuf77\nlz+cTqeCwWDU+Pm5yZg2bUpaP+g83l8lkaSOt/qjtn01c9N27WySqI9IDn1MD/o4dvQwPTLdx7hB\nvWfPHtlsNvX09Ogf//iHGhoaov6fORQKqaioSIWFhQqFQlHjLpcravz83GSk80+vFRe7Ev7ZzFAo\nHLWdzj+zmS2S6SMSo4/pQR/Hjh6mR7r6OOo/c/nb3/5Wu3fvVnt7u7797W+rublZ1dXVOnTokCSp\nu7tblZWVKi0tVW9vr8LhsAKBgAYGBuR2u1VeXq4DBw5E5lZUVIx5MQAAXEnivqK+lIaGBm3evFmt\nra2aNWuWampqlJubq9raWvl8PlmWpfr6euXn58vr9aqhoUFer1d5eXlqaWnJxBoAAMhaNsuyrIku\nIlY6345J5m2J1//nf6O2l1bNStv1swVvk6UHfUwP+jh29DA9JvytbwAAMLEIagAADEZQAwBgMIIa\nAACDEdQAABiMoAYAwGAENQAABiOoAQAwGEENAIDBCGoAAAxGUAMAYDCCGgAAgxHUAAAYjKAGAMBg\nBDUAAAYjqAEAMBhBDQCAwQhqAAAMRlADAGAwghoAAIMR1AAAGIygBgDAYAQ1AAAGI6gBADAYQQ0A\ngMEIagAADEZQAwBgMIIaAACD2RNN+Oqrr9TY2Khjx47JZrPpqaee0rlz5/Tggw/quuuukyR5vV7d\nfvvt6urqkt/vl91u17p167R48WKdPXtWGzZs0IkTJ+R0OtXc3Kzp06dnel0AAGSFhEG9f/9+SZLf\n79ehQ4f0/PPP67vf/a7uu+8+rV69OjJvcHBQ7e3t2rNnj8LhsHw+nxYuXKjOzk653W6tX79eb775\nptra2tTY2Ji5FQEAkEUSBvUtt9yiRYsWSZI++eQTFRUV6ciRIzp27Jj27dunmTNnatOmTerr61NZ\nWZkcDoccDodKSkrU39+v3t5e3X///ZKk6upqtbW1ZXRBAABkk4RBLUl2u10NDQ16++23tX37dn32\n2WdasWKF5s+frxdeeEE7d+7U3Llz5XK5Isc4nU4Fg0EFg8HIuNPpVCAQSHi9adOmyG7PHeWSLlZc\n7Iq73+nMT2n+lYq+pAd9TA/6OHb0MD0y3cekglqSmpub9dhjj8nj8cjv9+vqq6+WJC1ZskRNTU2q\nrKxUKBSKzA+FQnK5XCosLIyMh0IhFRUVJbzW0NCZVNdxWcXFLg0Oxv/hIBQKR20nmn8lSqaPSIw+\npgd9HDt6mB7p6mO8sE/4W9+vv/66XnrpJUlSQUGBbDabfvzjH6uvr0+S1NPTo3nz5qm0tFS9vb0K\nh8MKBAIaGBiQ2+1WeXm5Dhw4IEnq7u5WRUXFmBcEAMCVIuEr6ltvvVWPP/64Vq1apXPnzmnTpk26\n5ppr1NTUpLy8PM2YMUNNTU0qLCxUbW2tfD6fLMtSfX298vPz5fV61dDQIK/Xq7y8PLW0tIzHugAA\nyAo2y7KsiS4iVjrfjknmbYnX/+d/o7aXVs1K2/WzBW+TpQd9TA/6OHb0MD2MeOsbAABMHIIaAACD\nEdQAABiMoAYAwGAENQAABiOoAQAwGEENAIDBCGoAAAxGUAMAYDCCGgAAgxHUAAAYjKAGAMBgBDUA\nAAYjqAEAMBhBDQCAwQhqAAAMRlADAGAwghoAAIMR1AAAGIygBgDAYAQ1AAAGI6gBADAYQQ0AgMEI\nagAADEZQAwBgMIIaAACDEdQAABjMnmjCV199pcbGRh07dkw2m01PPfWU8vPztXHjRtlsNs2ZM0db\nt25VTk6Ourq65Pf7ZbfbtW7dOi1evFhnz57Vhg0bdOLECTmdTjU3N2v69OnjsTYAACa9hK+o9+/f\nL0ny+/2qq6vT888/r23btqmurk4dHR2yLEv79u3T4OCg2tvb5ff7tWvXLrW2tmp4eFidnZ1yu93q\n6OjQ0qVL1dbWlvFFAQCQLRK+or7lllu0aNEiSdInn3yioqIivf/++1qwYIEkqbq6WgcPHlROTo7K\nysrkcDjkcDhUUlKi/v5+9fb26v7774/MJagBAEhewqCWJLvdroaGBr399tvavn27Dh48KJvNJkly\nOp0KBAIKBoNyuVyRY5xOp4LBYNT4+bmJTJs2RXZ77mjWc0nFxa64+53O/JTmT0Ydb/VHbftq5qZ8\njgv7ko7zZdp415js9bLx/poI6erjZLiXL5TOerkX0yPTfUwqqCWpublZjz32mDwej8LhcGQ8FAqp\nqKhIhYWFCoVCUeMulytq/PzcRIaGzqSyhriKi10aHIz/w0EoFI7aTjR/MhrrGmP7OBl6Nt41JnO9\nZO5HJJbOPk6Ge/lC6aqXezE90tXHeGGf8P+oX3/9db300kuSpIKCAtlsNs2fP1+HDh2SJHV3d6uy\nslKlpaXq7e1VOBxWIBDQwMCA3G63ysvLdeDAgcjcioqKMS8IAIArRcJX1Lfeeqsef/xxrVq1SufO\nndOmTZs0e/Zsbd68Wa2trZo1a5ZqamqUm5ur2tpa+Xw+WZal+vp65efny+v1qqGhQV6vV3l5eWpp\naRmPdQEAkBUSBvWUKVP0y1/+8qLx3bt3XzTm8Xjk8XiixgoKCrR9+/YxlAgAwJWLDzwBAMBgBDUA\nAAYjqAEAMBhBDQCAwQhqAAAMRlADAGAwghoAAIMR1AAAGIygBgDAYAQ1AAAGI6gBADAYQQ0AgMEI\nagAADEZQAwBgMIIaAACDEdQAABiMoAYAwGAENQAABiOoAQAwGEENAIDBCGoAAAxGUAMAYDCCGgAA\ngxHUAAAYjKAGAMBgBDUAAAYjqAEAMBhBDQCAwezxdn755ZfatGmTPv74Yw0PD2vdunW65ppr9OCD\nD+q6666TJHm9Xt1+++3q6uqS3++X3W7XunXrtHjxYp09e1YbNmzQiRMn5HQ61dzcrOnTp4/HugAA\nyApxg3rv3r2aOnWqnnvuOZ06dUpLly7Vww8/rPvuu0+rV6+OzBscHFR7e7v27NmjcDgsn8+nhQsX\nqrOzU263W+vXr9ebb76ptrY2NTY2ZnxRAABki7hvfd922236yU9+IkmyLEu5ubk6cuSI3nvvPa1a\ntUqbNm1SMBhUX1+fysrK5HA45HK5VFJSov7+fvX29qqqqkqSVF1drZ6ensyvCACALBL3FbXT6ZQk\nBYNBPfLII6qrq9Pw8LBWrFih+fPn64UXXtDOnTs1d+5cuVyuqOOCwaCCwWBk3Ol0KhAIJFXUtGlT\nZLfnjnZNFykudsXd73TmpzR/MkrHGi88ZjL0bLxrTPZ6JvZqMkpXHyfDvXyhdNZr+loni0z3MW5Q\nS9Knn36qhx9+WD6fT3feeadOnz6toqIiSdKSJUvU1NSkyspKhUKhyDGhUEgul0uFhYWR8VAoFDku\nkaGhM6NZyyUVF7s0OBj/B4RQKBy1nWj+ZDTWNcb2cTL0bLxrTOZ6ydyPSCydfZwM9/KF0lUv92J6\npKuP8cI+7lvfn3/+uVavXq0NGzbo7rvvliStWbNGfX19kqSenh7NmzdPpaWl6u3tVTgcViAQ0MDA\ngNxut8rLy3XgwAFJUnd3tyoqKsa8GAAAriRxX1G/+OKLOn36tNra2tTW1iZJ2rhxo5555hnl5eVp\nxowZampqUmFhoWpra+Xz+WRZlurr65Wfny+v16uGhgZ5vV7l5eWppaVlXBYFAEC2iBvUjY2Nl/wt\nbb/ff9GYx+ORx+OJGisoKND27dvHWCIAAFcuPvAEAACDEdQAABiMoAYAwGAENQAABiOoAQAwGEEN\nAIDBCGoAAAxGUAMAYDCCGgAAgxHUAAAYjKAGAMBgBDUAAAYjqAEAMBhBDQCAwQhqAAAMRlADAGAw\nghoAAIMR1AAAGIygBgDAYAQ1AAAGI6gBADAYQQ0AgMEIagAADEZQAwBgMIIaAACDEdQAABiMoAYA\nwGD2eDu//PJLbdq0SR9//LGGh4e1bt06ffOb39TGjRtls9k0Z84cbd26VTk5Oerq6pLf75fdbte6\ndeu0ePFinT17Vhs2bNCJEyfkdDrV3Nys6dOnj9faAACY9OK+ot67d6+mTp2qjo4O/epXv1JTU5O2\nbdumuro6dXR0yLIs7du3T4ODg2pvb5ff79euXbvU2tqq4eFhdXZ2yu12q6OjQ0uXLlVbW9t4rQsA\ngKwQ9xX1bbfdppqaGkmSZVnKzc3V0aNHtWDBAklSdXW1Dh48qJycHJWVlcnhcMjhcKikpET9/f3q\n7e3V/fffH5lLUAMAkJq4Qe10OiVJwWBQjzzyiOrq6tTc3CybzRbZHwgEFAwG5XK5oo4LBoNR4+fn\nJmPatCmy23NHtaBLKS52xd3vdOZHbb/914+jtn01c6O2O97qj7t/NGLPmUiqNcWuMVFPLuXCY9Jx\nvliJepBqn1P9vo5Vsj1JtleX6ke6a55s0tGT2HOk+17OxPPDhZKpN9kazh+b7sfeWGW6h6MRr6Z0\nPP/FEzeoJenTTz/Vww8/LJ/PpzvvvFPPPfdcZF8oFFJRUZEKCwsVCoWixl0uV9T4+bnJGBo6k+o6\nLqu42KXBwfg/IIRC4bj7Y4+PnZ/o/MlIVMNYaxprzbF9nIgepHqNdJ8v1etd6vzJ3I+XO9/lznkl\nOd8TpzM/8rVp90UmHhupnj/Ve3G8HyuJZLqHo3G5mlJ5TMcTL+zj/h/1559/rtWrV2vDhg26++67\nJUk33HCDDh06JEnq7u5WZWWlSktL1dvbq3A4rEAgoIGBAbndbpWXl+vAgQORuRUVFWNeDAAAV5K4\nr6hffPFFnT59Wm1tbZH/X37iiSf09NNPq7W1VbNmzVJNTY1yc3NVW1srn88ny7JUX1+v/Px8eb1e\nNTQ0yOv1Ki8vTy0tLeOyKAAAskXcoG5sbFRjY+NF47t3775ozOPxyOPxRI0VFBRo+/btYywRAIAr\nFx94AgCAwQhqAAAMRlADAGAwghoAAIMR1AAAGIygBgDAYAQ1AAAGI6gBADAYQQ0AgMEIagAADEZQ\nAwBgMIIaAACDEdQAABiMoAYAwGAENQAABiOoAQAwGEENAIDBCGoAAAxGUAMAYDCCGgAAgxHUAAAY\njKAGAMBgBDUAAAYjqAEAMBhBDQCAwQhqAAAMRlADAGAwghoAAIMlFdSHDx9WbW2tJOnvf/+7qqqq\nVFtbq9raWv3hD3+QJHV1dWnZsmXyeDzav3+/JOns2bNav369fD6f1q5dq5MnT2ZoGQAAZCd7ogmv\nvPKK9u7dq4KCAknS0aNHdd9992n16tWROYODg2pvb9eePXsUDofl8/m0cOFCdXZ2yu12a/369Xrz\nzTfV1tamxsbGzK0GAIAsk/AVdUlJiXbs2BHZPnLkiN577z2tWrVKmzZtUjAYVF9fn8rKyuRwOORy\nuVRSUqL+/n719vaqqqpKklRdXa2enp7MrQQAgCyU8BV1TU2N/v3vf0e2S0tLtWLFCs2fP18vvPCC\ndu7cqblz58rlckXmOJ1OBYNBBYPByLjT6VQgEEiqqGnTpshuz011LZdVXOyKu9/pzE/p+Nj5ic6f\njEQ1jLWmdNR84TET0YNUr5Hu86V6vcudP9nrXqr+dNc82VzYk/Nfm3ZfZOKxker5U70Xx/uxkkim\nezga8WrKdH0JgzrWkiVLVFRUFPm6qalJlZWVCoVCkTmhUEgul0uFhYWR8VAoFDkukaGhM6mWdVnF\nxS4NDsb/ASEUCsfdH3t87PxE509GohrGWtNYa47t40T0INVrpPt8qV7vUudP5n683Pkud84ryfme\nOJ35ka9Nuy8y8dhI9fyp3ovj/VhJJNM9HI3L1ZTKYzqeeGGf8m99r1mzRn19fZKknp4ezZs3T6Wl\npert7VU4HFYgENDAwIDcbrfKy8t14MABSVJ3d7cqKipGuQQAAK5MKb+ifvLJJ9XU1KS8vDzNmDFD\nTU1NKiwsVG1trXw+nyzLUn19vfLz8+X1etXQ0CCv16u8vDy1tLRkYg0AAGStpIL62muvVVdXlyRp\n3rx58vv9F83xeDzyeDxRYwUFBdq+fXsaygQA4MrEB54AAGAwghoAAIMR1AAAGIygBgDAYAQ1AAAG\nI6gBADAYQQ0AgMEIagAADEZQAwBgMIIaAACDEdQAABiMoAYAwGAENQAABiOoAQAwGEENAIDBCGoA\nAAxGUAMAYDCCGgAAgxHUAAAYjKAGAMBgBDUAAAYjqAEAMBhBDQCAwQhqAAAMRlADAGAwghoAAIMR\n1AAAGCypoD58+LBqa2slScePH5fX65XP59PWrVs1MjIiSerq6tKyZcvk8Xi0f/9+SdLZs2e1fv16\n+Xw+rV27VidPnszQMgAAyE4Jg/qVV15RY2OjwuGwJGnbtm2qq6tTR0eHLMvSvn37NDg4qPb2dvn9\nfu3atUutra0aHh5WZ2en3G63Ojo6tHTpUrW1tWV8QQAAZJOEQV1SUqIdO3ZEto8ePaoFCxZIkqqr\nq/X++++rr69PZWVlcjgccrlcKikpUX9/v3p7e1VVVRWZ29PTk6FlAACQnRIGdU1Njex2e2TbsizZ\nbDZJktPpVCAQUDAYlMvlisxxOp0KBoNR4+fnAgCA5NkTT4mWk/PfbA+FQioqKlJhYaFCoVDUuMvl\niho/PzcZ06ZNkd2em2ppl/X2Xz+O2vbVzI3adjrz4x5fXOyK2o6dH7u/463+VEtMWEOs2DUlqil2\nf+zxsWJ7FHvORNeLlUxPEvUg1ZpT/b7G1nipHsSbn2qPYyVTf6rnTHSNWIl6kOj7ONbzJ3JhT85/\nnejxHe8cl5LqvZzo+5aovlR7msx9lmwN59ea6mMv1TWkKt2PJSn9j+8L75NE98xYpRzUN9xwgw4d\nOqSbbrpJ3d3d+s53vqPS0lL94he/UDgc1vDwsAYGBuR2u1VeXq4DBw6otLRU3d3dqqioSOoaQ0Nn\nUl7I5RQXuxQKhaPGBgejX9nH7o+VaH6q5xsPY60p9vjiYlfUWKIexBqPnmT6+xprNGtyOvMve9x4\n3EeprmmsPU31/ImcPz6VPiaqYazHp/u+G+v8ZAwOBqIe02N9fjDhOe9Cl/oepvvxff742OfG0YoX\n9ikHdUNDgzZv3qzW1lbNmjVLNTU1ys3NVW1trXw+nyzLUn19vfLz8+X1etXQ0CCv16u8vDy1tLSM\naSEAAFxpkgrqa6+9Vl1dXZKk66+/Xrt3775ojsfjkcfjiRorKCjQ9u3b01AmAABXJj7wBAAAgxHU\nAAAYjKAGAMBgBDUAAAYjqAEAMBhBDQCAwQhqAAAMRlADAGAwghoAAIMR1AAAGIygBgDAYAQ1AAAG\nI6gBADAYQQ0AgMEIagAADEZQAwBgMIIaAACDEdQAABiMoAYAwGAENQAABiOoAQAwGEENAIDBCGoA\nAAxGUAMAYDCCGgAAgxHUAAAYjKAGAMBgBDUAAAazj/bAu+66S4WFhZKka6+9Vg899JA2btwom82m\nOXPmaOvWrcrJyVFXV5f8fr/sdrvWrVunxYsXp614AACy3aiCOhwOy7Istbe3R8Yeeugh1dXV6aab\nbtKWLVu0b98+3XjjjWpvb9eePXsUDofl8/m0cOFCORyOtC0AAIBsNqqg7u/v1xdffKHVq1fr3Llz\n+ulPf6qjR49qwYIFkqTq6modPHhQOTk5Kisrk8PhkMPhUElJifr7+1VaWprWRQAAkK1GFdRXXXWV\n1qxZoxUrVuijjz7S2rVrZVmWbDabJMnpdCoQCCgYDMrlckWOczqdCgaDCc8/bdoU2e25oyntkpzO\n/Kjt4mJX3P2xEs1P9XzjYaw1xR4fO5aoB7HGoyeZ/r7GGu2aLnfceNxHqa5prD1N9fyJXHh8sn1M\nVMNYj0/3fTfW+ck4f87z/471+cGE57wLXep7mO7H94XHp3ofp2pUQX399ddr5syZstlsuv766zV1\n6lQdPXo0sj8UCqmoqEiFhYUKhUJR4xcG9+UMDZ0ZTVmXVFzsUigUjhobHAxEbcfuj5VofqrnGw9j\nrSn2+OJiV9RYoh7EGo+eZPr7Gms0a3I68y973HjcR6muaaw9TfX8iZw/PpU+JqphrMen+74b6/xk\nDA4Goh7TY31+MOE570KX+h6m+/F9/vjY58bRihf2o/qt71dffVXPPvusJOmzzz5TMBjUwoULdejQ\nIUlSd3e3KisrVVpaqt7eXoXDYQUCAQ0MDMjtdo/mkgAAXJFG9Yr67rvv1uOPPy6v1yubzaZnnnlG\n06ZN0+bNm9Xa2qpZs2appqZGubm5qq2tlc/nk2VZqq+vV36+WW+RAABgslEFtcPhUEtLy0Xju3fv\nvmjM4/HI4/GM5jIAAFzx+MATAAAMRlADAGAwghoAAIMR1AAAGIygBgDAYAQ1AAAGI6gBADAYQQ0A\ngMEIagAADEZQAwBgMIIaAACDEdQAABiMoAYAwGAENQAABiOoAQAwGEENAIDBCGoAAAxGUAMAYDCC\nGgAAgxHUAAAYjKAGAMBgBDUAAAYjqAEAMBhBDQCAwQhqAAAMRlADAGAwghoAAIMR1AAAGMye6QuM\njIzoySef1D//+U85HA49/fTTmjlzZqYvCwBAVsj4K+p33nlHw8PD+t3vfqdHH31Uzz77bKYvCQBA\n1sh4UPf29qqqqkqSdOONN+rIkSOZviQAAFnDZlmWlckLPPHEE7r11lt18803S5IWLVqkd955R3Z7\nxt91BwBg0sv4K+rCwkKFQqHI9sjICCENAECSMh7U5eXl6u7uliT97W9/k9vtzvQlAQDIGhl/6/v8\nb33/61//kmVZeuaZZzR79uxMXhIAgKyR8aAGAACjxweeAABgMIIaAACDZU1Qj4yMaMuWLVq5cqVq\na2t1/PjxqP3vvvuuli9frpUrV6qrq2uCqjRfoj6+8cYbWrFihe655x5t2bJFIyMjE1SpuRL18LzN\nmzfr5z//+ThXN3kk6mNfX598Pp+8Xq8eeeQRhcPhCarUbIn6uHfvXt11111avny5Ojo6JqjKyeHw\n4cOqra29aDzj+WJlibfeestqaGiwLMuyPvzwQ+uhhx6K7BseHrZuueUW69SpU1Y4HLaWLVtmDQ4O\nTlSpRovXxy+++ML63ve+Z505c8ayLMuqr6+33nnnnQmp02TxenheZ2en5fF4rOeee268y5s04vVx\nZGTE+sEPfmB99NFHlmVZVldXlzUwMDAhdZou0f24cOFCa2hoyAqHw5HnSVzs5Zdftu644w5rxYoV\nUePjkS9Z84o63iegDQwMqKSkRF/72tfkcDhUUVGhv/zlLxNVqtHi9dHhcMjv96ugoECSdO7cOeXn\n509InSZL9Gl8f/3rX3X48GGtXLlyIsqbNOL18dixY5o6dap+85vf6N5779WpU6c0a9asiSrVaInu\nx29961sKBAIaHh6WZVmy2WwTUabxSkpKtGPHjovGxyNfsiaog8GgCgsLI9u5ubk6d+5cZJ/L5Yrs\nczqdCgaD417jZBCvjzk5OZoxY4Ykqb29XWfOnNHChQsnpE6Txevhf/7zH+3cuVNbtmyZqPImjXh9\nHBoa0ocffqh7771Xv/71r/XnP/9ZPT09E1Wq0eL1UZLmzJmj5cuX6/vf/74WLVqkoqKiiSjTeDU1\nNZf8sK7xyJesCep4n4AWuy8UCkU1Fv+V6JPkRkZG1NzcrIMHD2rHjh389H0J8Xr4pz/9SUNDQ3rg\ngQf08ssv64033tBrr702UaUaLV4fp06dqpkzZ2r27NnKy8tTVVUVf0fgMuL1sb+/X++995727dun\nd999VydPntQf//jHiSp1UhqPfMmaoI73CWizZ8/W8ePHderUKQ0PD+uDDz5QWVnZRJVqtESfJLdl\nyxaFw2G1tbVF3gJHtHg9/OEPf6jXXntN7e3teuCBB3THHXdo2bJlE1Wq0eL18Rvf+IZCoVDkF6M+\n+OADzZkzZ0LqNF28PrpcLl111VXKz89Xbm6upk+frtOnT09UqZPSeORL1nzo9pIlS3Tw4EHdc889\nkU9A+/3vf68zZ85o5cqV2rhxo9asWSPLsrR8+XJdffXVE12ykeL1cf78+Xr11VdVWVmpH/3oR5L+\nP3iWLFkywVWbJdG9iOQk6uPPfvYzPfroo7IsS2VlZVq0aNFEl2ykRH1cuXKlfD6f8vLyVFJSorvu\numuiS54UxjNf+GQyAAAMljVvfQMAkI0IagAADEZQAwBgMIIaAACDEdQAABiMoAYAwGAENQAABiOo\nAQAw2P8BoGIGEjXF9yYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1240f4ba8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--- emcee sampling vi 69.027498960495 seconds ---\n"
     ]
    }
   ],
   "source": [
    "emcee_trace = sampler.chain[:, nburn:, :].reshape(-1, ndim).T\n",
    "plt.hist(emcee_trace[0], alpha=0.5, bins=100, label=\"delta\")\n",
    "plt.show()\n",
    "print(\"--- emcee sampling vi %s seconds ---\" % (end_time - start_time))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Edward"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "sess = ed.get_session()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# FORWARD MODEL, Prior\n",
    "mu = StudentT(1.0,[0.0],[1.0])\n",
    "delta = StudentT(1.0,[0.0],[1.0])\n",
    "sigma = Gamma([2.0],[2.0])\n",
    "alpha = delta * sigma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "x = Normal(tf.tile(mu+delta,[num_points]), tf.tile(sigma,[num_points]))\n",
    "y = Normal(tf.tile(mu,[num_points]), tf.tile(sigma,[num_points]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Mean and delta are best approximated by the <b>NormalWithSoftplusScale</b> distribution with the softplus function on the scale (sigma)(variance should be positive) parameter since Cauchy and Normal are both defined on positive and negative scales. Sigma as a variance should be always positive: we approximate sigma with <b>GammaWithSoftplusConcentrationRate</b> distribution ensuring the positive concentration and rate parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# BACKWARD MODEL\n",
    "q_mu = NormalWithSoftplusScale(loc=tf.Variable([0.0]), scale=tf.Variable([1.0]))\n",
    "q_sigma = GammaWithSoftplusConcentrationRate(tf.nn.softplus(tf.Variable([1.0])), tf.nn.softplus(tf.Variable([1.0])))\n",
    "q_delta = NormalWithSoftplusScale(loc=tf.Variable([0.0]), scale=tf.Variable([1.0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# INFERENCE\n",
    "inference = ed.KLqp({delta: q_delta, mu: q_mu, sigma: q_sigma}, data={x: xdata, y: ydata})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10000/10000 [100%] ██████████████████████████████ Elapsed: 1537s | Loss: 3539549.500\n",
      "--- Edward vi 1537.8405618667603 seconds ---\n"
     ]
    }
   ],
   "source": [
    "start_time = time.time()\n",
    "inference.run(n_iter=10000, n_print=200, n_samples=10)\n",
    "end_time = time.time()\n",
    "print(\"--- Edward vi %s seconds ---\" % (end_time - start_time))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "T = 10000\n",
    "delta_sample = sess.run(delta.sample(sample_shape=T))\n",
    "q_delta_sample = sess.run(q_delta.sample(sample_shape=T))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfUAAAFJCAYAAACCQLQfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XtcVHX+x/HXwADqMAgkll2oMMnbkqCh5i27Ua22bqkF\nLVm6pqypmBGkKZrmZUvyUnjJS66uEmbbr3uauZKXJWMzk7ILpamVoWHMoNzk/P5wHcUbqIMMx/fz\n8ehRc+bMme/ne0685/s9Z85YDMMwEBERkTrPq7YbICIiIu6hUBcRETEJhbqIiIhJKNRFRERMQqEu\nIiJiEgp1ERERk7DWdgNOpbz8MAUFB2u7GTUmKKiBaeszc22g+uo61Vd3mbk2gJAQu1u245EjdavV\nu7abUKPMXJ+ZawPVV9epvrrLzLW5k0eGuoiIiJw9hbqIiIhJKNRFRERMQqEuIiJiEh559buIXHj9\np3zk1u0tTLnFrdsTkapppC4iImISCnURERGTUKiLiIiYhEJdRETEJBTqIiIiJqFQFxERMQmFuoiI\niElU+T31srIyUlJS2LNnD15eXkyYMAGr1UpKSgoWi4VmzZqRmpqKl5cXmZmZZGRkYLVaSUhIoHv3\n7hQXF5OUlMT+/fux2WxMnTqV4ODgC1GbiIi4gWEYWCyW2m6GVEOVob5u3TrKy8vJyMhgw4YNTJ8+\nnbKyMhITE2nfvj1jx45lzZo1tGnThiVLlrBy5UpKSkqIi4ujU6dOLF++nPDwcIYOHco777xDeno6\nTz/99IWoTUTkrLn7Jjxny5037XnssUdp0KABf//79Gqt/9//fsqwYYOZP/8fNG/ekl9/3cuUKRMZ\nO3YCgYGBbmuXJ1mwYC4ZGUtZvfrj2m6KW1Q5/X7ttddy+PBhKioqcDqdWK1WcnNziY6OBqBr165s\n3LiRrVu3EhkZia+vL3a7ndDQULZv305OTg5dunRxrbtp06aarUhERNzi008/4ZNP9De7LqlypN6g\nQQP27NnDXXfdRUFBAXPmzGHz5s2uqRibzYbD4cDpdGK3H/uRd5vNhtPprLT86LrV4a4fjPdUZq7P\nzLWB+etzF0/tJ09t11Hn277jX+/ra8XX11rtbQYGNnD9OyTEjt1eD4BLLrERHFz7/VYT+85m88Ni\nsXj8cVFdVYb6K6+8QufOnRk5ciQ///wz/fr1o6yszPV8UVERAQEB+Pv7U1RUVGm53W6vtPzoutWR\nn1+98K+LQkLspq3PzLWB+etzJ0/sp7qw/861fQcPHuTll2fxwQersFggNjae0tJyrNZy1zbLy8t5\n5ZX5vPvuWxw4UMC11zYlIWEo7dodmXk9cOCg69+LFy9j0qTxAHTs2JFHHhnIgAGD2LdvH/PmvUR2\n9iYOHCggMDCIW265nYSEofj6+p6ybYcOHWLGjOfZtGk9DoeTa665hn79BtCt27FTDatWvUdm5nJ+\n+CEPgGbNrmfw4Mdo0yYKgGnTnqWg4HdatvwDK1Ysx+EopEOHTjz11BgyM5fz+usrqKg4zO2338mw\nYSPx8vLi3Xff4oUX/s6ECVOZPv159u37lebNWzJ8+EiaNbsegKKiEgzDqNTvq1e/z5Ili9i160dC\nQhrTt28svXs/4Ho+N3cb6ekz+Oabr7FarbRteyOPPZbIZZc1Oad9B+77wFLl9HtAQIBrpN2wYUPK\ny8tp2bIl2dnZAGRlZdGuXTsiIiLIycmhpKQEh8NBXl4e4eHhREVFsW7dOte6bdu2dUvDRUTkmHHj\nRrF69Wr+9rdhpKSM4cMPV7Ft29ZK6/z978+SkbGUPn1imTTpea6++hqeeGIYX3zx+Unb69ixM/36\nDQBg2rRZ9OzZi4qKCkaOHMo333zN448nk5b2IjExd7NixXLefPP107Ztxoznycn5lOHDk3j++Rlc\nc00YY8aksGPHDwCsXfshEyaMpWPHTjz33AxGjUrF6XSQmvpUpUHkJ59kk5W1liefHM3f/jacrKy1\nDBjwELm5X/D00+Pp2fPPvPbaq6xZs8r1mtLSUiZOHMt99/Vl3LhJlJQUM3ToYAoKfjtlW997723G\nj3+aNm2imDr1Be66qwezZr3AsmX/AMDpdJKUNJxGjUKYMmUaycmj+eab7aSmjqrmnqpZVY7UH374\nYUaNGkVcXBxlZWWMGDGC1q1bM2bMGNLS0ggLCyMmJgZvb2/i4+OJi4vDMAxGjBiBn58fsbGxJCcn\nExsbi4+PD9OmTbsQdYmIXDS+++5bNm5czwsvvMCNNx65hqlFi9b07XuPa52dO3fw7rtvkZz8ND17\n9gKgQ4eb2L9/Hy+/PJuZM+dU2mZQUBBXXHElANdf34LAwED27v0Fu91OYmIS113XDIC2bW8kO3sT\nW7b8t9Jo9nhbt27hxhujueWW2wD4wx9uICgomMOHDwOwe/du7r23DwMGDHK9xmr1YfToJHbt2klY\n2HUAHDp0kIkT/06jRo0A+OCDd9mx43sWLlxCgwY22rfvyKpV7/Hll7ncfvudABw+fJi//nUwvXr1\nBqBVqz/Qp09P/vWv1+jf/9FK7ayoqGDu3Je44467ePzxZACioztgsVh45ZUF/PnPfdix4wcKC3+n\nT58HaN06AoCGDQP5738/paKiAi+v2v2meJWhbrPZmDFjxknLly5detKyvn370rdv30rL6tevz8yZ\nM8+jiSIiciZHR9pdu3bl0CEDgEaNGtGq1R9c63z22acAdOzYifLyctfyDh06MXfui5VGxKdz6aWX\n8eKL86ioqGDXrh/ZtetHvvvuGwoKfuPSSy877etuuCGSt956g/3793HTTV3o1KkLQ4eOcD0fH/8w\nAA6Hg507d7Br107Wrz8yw1taeqxdjRtf6gp0gODgYCoqDtOggc21LCCgIU5n5VMYt94a4/rvoKAg\nWrWK4PPPt5zUzl27fmTfvvxT9NFNzJ8/h6++yqV58xYEBDQkOXkEt956Bx07dqZt2xuJjPSMWWj9\nnrqISB3ncBRitVrx9/fn0KFjgRYcfAmHDh05T/77778D0KvXXafcxu+/H6jWe7399hvMmzeb337b\nzyWXNKJly9b4+fkBxmlfk5j4BI0ahfD++++yYcPHTJvmRYcOnRg1KpXAwED279/HlCkT+M9/NuLj\n48M114TRpMnR89PHttugQYOTtu3nV++M7fX19at0ETdAYGAgu3btPGndo30wfvzTjB9/8lev9+3b\nR4MGNl58cR6vvDKf9957h9dfX4G/v534+Id58MF+Z2zLhaBQFxGp4wICjlzvVFhYCBy7SUxh4e/4\n+PgA4O/vj8ViYfbshVit3idto2HDQODkoDveZ5/lMHXqszz88F+5996+BAUFATBw4ENnfJ2fXz0G\nDBjEgAGD+PHHHaxdu4bFixcwf/5snnjiKcaPf5pff/2VOXMWcf31zbFarWzatJ5169aeXUecQmlp\nCcXFxdSrdyz8j1zgd/JN0Pz9/QF4/PFkWrZsddLzTZpcDkBYWFOeeWYyZWVlfP75Z6xYsZzZs2fR\npk1bWrVqfd5tPh+6TayISB0XFXVk6nfVqmMXiBUWFpKbu831OCKiDYZhcPBgEc2bt3T9s3nzJ7z6\n6jK8vU8O+hPPD+fmfoHFYqFfvwGuQN+3L5+8vDwM49Qj9cOHDxMf35fMzGUAhIYeufK9Vas/sHfv\nLwBs2/YFt912B61atcZqPTLWzM4+8v340233bGzYcOzGMgUFv7Ft21ZXnx0vNPQaGjZsSH7+r5X6\n6Pfff+fll+fgdDr5z3820qPH7RQUFODj40O7dtGMGPEkgKue2qSRuohIHRcaeg0xMXcxadIk9u8v\n5NJLL2Pp0kWVzgs3a3Y9N998C888M4b+/R/lmmuu5bPPcli8eAFxcQ+d8gIvf/8j09br1n1EdHQH\nWrRoRUVFBTNmTKN791vZu/cX/vGPhZSVlVJcXHzKtnl7e9OyZWsWLnwZX19fQkOvITd3G1u3biEp\n6SkAWrRoybvvvkXTptdhtweQlbWWN95YCXDa7Z6NF16YyqFDRQQGBrFo0XwCAhq6Lpw7ntVqpX//\nR5k16wXgyEWAP//8E3PnvsiVV4Zy+eVX/G8q32D06CQefLAfPj4+ZGYuw9/fTlRUu/Nu6/lSqIuI\nHMedt2m9kFJSxrJ48VwWLpxLWVkZPXr8iUaNGlNSciwUx46dyPz5c1i69JX/XdzWhMGDHyM2Nv6U\n22zXLpro6I5Mn/4cPXv24vHHkxk6dAQrVmTwzjtv0rhxY7p3vw2r1cqKFcspLS095XfVExOTqF+/\nPv/4xyLX+w4ZkkiPHkeuwh81KpVp06YwadIz+Pn50rRpM2bOnENS0nC2bfvivC9Ce+yxESxYMI8D\nB36jbdsbmThxqmuq/UT33Xc/fn71ePXVf5KRsZSAgIbcfPNtPPro37BYLAQENOT552cxZ84sJkwY\nS3l5GS1btmb69HSPuJWuxXDH3EYN8PQbRJyPunADjHNl5trA3PW5+57nnhiOZt5/YO76zqW2d999\ni0mTxvP22x96ROCeyQW7+YyIiIjUDQp1ERERk1Coi4iIKd19d0/Wr//U46fe3UmhLiIiYhIKdRER\nEZNQqIuIiJiEQl1ERMQkFOoiIiImoVAXERExCYW6iIiISSjURURETEKhLiIiYhIKdREREZNQqIuI\niJiEQl1ERMQkrNVZae7cuXz00UeUlZURGxtLdHQ0KSkpWCwWmjVrRmpqKl5eXmRmZpKRkYHVaiUh\nIYHu3btTXFxMUlIS+/fvx2azMXXqVIKDg2u6LhERkYtOlSP17OxsPvvsM5YvX86SJUv45ZdfmDx5\nMomJiSxbtgzDMFizZg35+fksWbKEjIwMFixYQFpaGqWlpSxfvpzw8HCWLVtGr169SE9PvxB1iYiI\nXHSqDPX169cTHh7OkCFDGDx4MDfffDO5ublER0cD0LVrVzZu3MjWrVuJjIzE19cXu91OaGgo27dv\nJycnhy5durjW3bRpU81WJCIicpGqcvq9oKCAn376iTlz5rB7924SEhIwDAOLxQKAzWbD4XDgdDqx\n2+2u19lsNpxOZ6XlR9etjpAQe9Ur1WFmrs/MtYH563MXT+0nT22Xu5i5PjPX5i5VhnpgYCBhYWH4\n+voSFhaGn58fv/zyi+v5oqIiAgIC8Pf3p6ioqNJyu91eafnRdasjP7964V8XhYTYTVufmWsD89fn\nTp7YT2bff2auz8y1gfs+sFQ5/d62bVs+/vhjDMNg7969HDp0iI4dO5KdnQ1AVlYW7dq1IyIigpyc\nHEpKSnA4HOTl5REeHk5UVBTr1q1zrdu2bVu3NFxEREQqq3Kk3r17dzZv3kzv3r0xDIOxY8dy5ZVX\nMmbMGNLS0ggLCyMmJgZvb2/i4+OJi4vDMAxGjBiBn58fsbGxJCcnExsbi4+PD9OmTbsQdYmIiFx0\nLIZhGLXdiFMx+zSLWeszc21g7vr6T/nIrdtbmHKLW7fnDmbef2Du+sxcG1zA6XcRERGpGxTqIiIi\nJqFQFxERMQmFuoiIiEko1EVERExCoS4iImISCnURERGTUKiLiIiYhEJdRETEJBTqIiIiJqFQFxER\nMQmFuoiIiEko1EVERExCoS4iImISCnURERGTUKiLiIiYhEJdRETEJBTqIiIiJqFQFxERMQmFuoiI\niEko1EVERExCoS4iImIS1Qr1/fv3061bN/Ly8ti5cyexsbHExcWRmppKRUUFAJmZmdx777307duX\ntWvXAlBcXMzQoUOJi4tj4MCB/PbbbzVXiYiIyEWuylAvKytj7Nix1KtXD4DJkyeTmJjIsmXLMAyD\nNWvWkJ+fz5IlS8jIyGDBggWkpaVRWlrK8uXLCQ8PZ9myZfTq1Yv09PQaL0hERORiVWWoT506lQce\neIDGjRsDkJubS3R0NABdu3Zl48aNbN26lcjISHx9fbHb7YSGhrJ9+3ZycnLo0qWLa91NmzbVYCki\nIiIXN+uZnnz99dcJDg6mS5cuzJs3DwDDMLBYLADYbDYcDgdOpxO73e56nc1mw+l0Vlp+dN3qCgmx\nV71SHWbm+sxcG5i/Pnfx1H7y1Ha5i5nrM3Nt7nLGUF+5ciUWi4VNmzbx1VdfkZycXOm8eFFREQEB\nAfj7+1NUVFRpud1ur7T86LrVlZ9f/Q8AdU1IiN209Zm5NjB/fe7kif1k9v1n5vrMXBu47wPLGaff\n//nPf7J06VKWLFlCixYtmDp1Kl27diU7OxuArKws2rVrR0REBDk5OZSUlOBwOMjLyyM8PJyoqCjW\nrVvnWrdt27ZuabSIiIic7Iwj9VNJTk5mzJgxpKWlERYWRkxMDN7e3sTHxxMXF4dhGIwYMQI/Pz9i\nY2NJTk4mNjYWHx8fpk2bVhM1iIiICGAxDMOo7UacitmnWcxan5lrA3PX13/KR27d3sKUW9y6PXcw\n8/4Dc9dn5trgAk2/i4iISN2hUBcRETEJhbqIiIhJKNRFRERMQqEuIiJiEgp1ERERk1Coi4iImIRC\nXURExCQU6iIiIiahUBcRETEJhbqIiIhJKNRFRERMQqEuIiJiEgp1ERERk1Coi4iImIRCXURExCQU\n6iIiIiahUBcRETEJhbqIiIhJKNRFRERMQqEuIiJiEgp1ERERk7Ce6cmysjJGjRrFnj17KC0tJSEh\ngeuuu46UlBQsFgvNmjUjNTUVLy8vMjMzycjIwGq1kpCQQPfu3SkuLiYpKYn9+/djs9mYOnUqwcHB\nF6o2ERGRi8oZR+pvvvkmgYGBLFu2jPnz5zNhwgQmT55MYmIiy5YtwzAM1qxZQ35+PkuWLCEjI4MF\nCxaQlpZGaWkpy5cvJzw8nGXLltGrVy/S09MvVF0iIiIXnTOO1O+8805iYmIAMAwDb29vcnNziY6O\nBqBr165s2LABLy8vIiMj8fX1xdfXl9DQULZv305OTg5//etfXesq1EVERGrOGUPdZrMB4HQ6GTZs\nGImJiUydOhWLxeJ63uFw4HQ6sdvtlV7ndDorLT+6bnWFhNirXqkOM3N9Zq4NzF+fu3hqP3lqu9zF\nzPWZuTZ3OWOoA/z8888MGTKEuLg4evbsyXPPPed6rqioiICAAPz9/SkqKqq03G63V1p+dN3qys+v\n/geAuiYkxG7a+sxcG5i/PnfyxH4y+/4zc31mrg3c94HljOfU9+3bR//+/UlKSqJ3794AtGzZkuzs\nbACysrJo164dERER5OTkUFJSgsPhIC8vj/DwcKKioli3bp1r3bZt27ql0SIiInKyM47U58yZQ2Fh\nIenp6a7z4aNHj2bixImkpaURFhZGTEwM3t7exMfHExcXh2EYjBgxAj8/P2JjY0lOTiY2NhYfHx+m\nTZt2QYoSERG5GFkMwzBquxGnYvZpFrPWZ+bawNz19Z/ykVu3tzDlFrduzx3MvP/A3PWZuTa4QNPv\nIiIiUnco1EVERExCoS4iImISCnURERGTUKiLiIiYhEJdRETEJKq8o5yIeCZ3fwVNROo+jdRFRERM\nQqEuIiJiEgp1ERERk1Coi4iImIRCXURExCQU6iIiIiahUBcRETEJhbqIiIhJKNRFRERMQqEuIiJi\nEgp1ERERk1Coi4iImIRCXURExCQU6iIiIiahUBcRETGJGv899YqKCsaNG8fXX3+Nr68vEydO5Oqr\nr67ptxUREbno1PhI/cMPP6S0tJRXX32VkSNHMmXKlJp+SxERkYtSjY/Uc3Jy6NKlCwBt2rRh27Zt\nNf2WIh6r/5SParsJF4y7a12YcotbtydiRjUe6k6nE39/f9djb29vysvLsVrP/NYhIfaablqtMnN9\nZq4Nzq++t6b9yY0tkXOh47PuMnNt7lLj0+/+/v4UFRW5HldUVFQZ6CIiInL2ajzUo6KiyMrKAmDL\nli2Eh4fX9FuKiIhclCyGYRg1+QZHr37/5ptvMAyDSZMm0bRp05p8SxERkYtSjYe6iIiIXBi6+YyI\niIhJKNRFRERMotYuQ1+9ejXvv/8+06ZNA45cRPfss8/i7e1N586deeyxxwB48cUX+fe//43VamXU\nqFFERETw22+/8cQTT1BcXEzjxo2ZPHky9evXr61STmvevHl8/PHHABQWFrJv3z42bNjA6tWrmTp1\nKk2aNAFg6NChREdHn7JWT2UYBl27duWaa64BjtyDYOTIkWe1Hz2Zw+EgKSkJp9NJWVkZKSkpREZG\nmmLfnchMd30sKytj1KhR7Nmzh9LSUhISEmjSpAmDBg1yHauxsbHcfffdZGZmkpGRgdVqJSEhge7d\nu9du46vpz3/+s+trwldeeSWDBw8mJSUFi8VCs2bNSE1NxcvLq87V9/rrr/Ovf/0LgJKSEr766ite\nffVVU+y7zz//nOeff54lS5awc+fOau+v4uJikpKS2L9/PzabjalTpxIcHHzmNzNqwYQJE4yYmBgj\nMTHRteyee+4xdu7caVRUVBh//etfjdzcXGPbtm1GfHy8UVFRYezZs8e49957Xa9fuXKlYRiGMXfu\nXGPRokW1UcZZefTRR42PP/7YMAzDSEtLM95///1Kz5+uVk+1Y8cOY9CgQSctP5v96MlmzJjhOq7y\n8vKMXr16GYZhjn13og8++MBITk42DMMwPvvsM2Pw4MG13KJz99prrxkTJ040DMMwCgoKjG7duhmZ\nmZnGggULKq3366+/Gj169DBKSkqMwsJC1397uuLiYuNPf/pTpWWDBg0y/vOf/xiGYRhjxowxVq1a\nVWfrO2rcuHFGRkaGKfbdvHnzjB49ehh9+vQxDOPs9tfChQuNmTNnGoZhGG+//bYxYcKEKt+vVqbf\no6KiGDdunOux0+mktLSU0NBQLBYLnTt3ZuPGjeTk5NC5c2csFguXX345hw8f5rfffqt0l7quXbuy\ncePG2iij2latWkVAQACdO3cGIDc3l5UrVxIXF8eUKVMoLy8/ba2eKjc3l7179xIfH8/AgQP5/vvv\nz3o/erKHH36YBx54AIDDhw/j5+cHmGPfnchMd3288847GT58OHBkNsnb25tt27bx73//mwcffJBR\no0bhdDrZunUrkZGR+Pr6YrfbCQ0NZfv27bXc+qpt376dQ4cO0b9/fx566CG2bNlCbm4u0dHRwLG/\nh3W1PoAvvviC7777jvvvv98U+y40NJRZs2a5Hp/N/jox6zZt2lTl+9Xo9PuKFStYvHhxpWWTJk3i\n7rvvJjs727XsxLvO2Ww2du3ahZ+fH4GBgZWWOxwOnE4ndru90rLadrpaIyIimDt3Lmlpaa7lnTp1\n4rbbbuPKK68kNTWVjIwMnE7nKWutcqrlAjhVbWPHjuXRRx/lrrvu4tNPPyUpKYmXXnrprPajJ9QG\nZ953+fn5JCUlMWrUKKDu7bvqONe7Pnoim80GHKlp2LBhJCYmUlpaSp8+fWjdujWzZ8/mpZdeonnz\n5q6/IUdf53Q6a6vZ1VavXj0GDBhAnz592LFjBwMHDsQwDCwWC3Dqv5FHl9eF+gDmzp3LkCFDAIiI\niKjz+y4mJobdu3e7Hp/N/jqXrKvR/2v79OlDnz59qlzvxLvOFRUVERAQgI+Pz0nL7Xa7a/169eq5\n1q1tp6v1u+++IyAgoNI5yvvuu8/V5ltvvZUPPviA5s2bn7JWT3Cq2g4dOoS3tzcA7dq149dff8Vm\ns53VfvQUp9t3X3/9NY8//jhPPvmk65N1Xdt31WG2uz7+/PPPDBkyhLi4OHr27ElhYaFrn91+++1M\nmDCBdu3a1cl9du2113L11VdjsVi49tprCQwMJDc31/X80f/nTvU3tS7UV1hYyA8//ECHDh2AI/vL\nLPvuKC+vYxPkVe2v45dXN+s84up3f39/fHx8+PHHHzEMg/Xr19OuXTuioqJYv349FRUV/PTTT1RU\nVBAcHExUVBTr1q0DICsri7Zt29ZyBae3ceNGunbt6npsGAb33HMPv/zyCwCbNm2iVatWp63VU734\n4ouu0e327dtp0qQJdrv9rPajJ/vuu+8YPnw406ZNo1u3boB59t2JzHTXx3379tG/f3+SkpLo3bs3\nAAMGDGDr1q3AsX0WERFBTk4OJSUlOBwO8vLy6kTdr732muuXLvfu3YvT6aRTp06umc+srCzatWtX\nZ+vbvHkzHTt2dD020747qmXLltXeX+eSdR7zcXz8+PE88cQTHD58mM6dO3PDDTcAR0aB999/PxUV\nFYwdOxaAhIQEkpOTyczMJCgoyHUFvSf64Ycf6NSpk+uxxWJh4sSJPPbYY9SrV4+mTZvSt29ffHx8\nTlmrp3r00UdJSkpi3bp1eHt7M3nyZODs9qMnmzZtGqWlpTz77LPAkQ+es2fPNsW+O9Htt9/Ohg0b\neOCBB1x3fayr5syZQ2FhIenp6aSnpwOQkpLCpEmT8PHxoVGjRkyYMAF/f3/i4+OJi4vDMAxGjBjh\num7Ck/Xu3ZunnnqK2NhYLBYLkyZNIigoiDFjxpCWlkZYWBgxMTF4e3vXyfp++OEHrrzyStfjcePG\nMWHCBFPsu6OSk5Orvb9iY2NJTk4mNjYWHx+famWd7ignIiJiEh4x/S4iIiLnT6EuIiJiEgp1ERER\nk1Coi4iImIRCXURExCQU6iIiIiahUBcRETEJhbqIiIhJKNRFRERMQqEuIiJiEgp1ERERk1Coi4iI\nmITH/Erb8fLzq/4h+LoqKKgBBQUHa7sZpqY+rlnq35ql/q15ntjHISHu+U14jdQvMKvVu7abYHrq\n45ql/q1Z6t+aZ+Y+VqiLiIiYhEJdRETEJBTqIiIiJqFQFxERMQmPvPpdRGpH/ykfVXq8MOWWWmqJ\niJwLhbqIVJtCX8SzafpdRETEJBTqIiIiJqFQFxERMQmFuoiIiEnoQjmRi9iJF76JSN2mUBcxMV2t\nLnJxOadQP3z4ME8//TQ//PADFouF8ePH4+fnR0pKChaLhWbNmpGamoqXlxeZmZlkZGRgtVpJSEig\ne/fu7q5BRGqIRvJiGAYWi6W2myHVdE6hvnbtWgAyMjLIzs7mhRdewDAMEhMTad++PWPHjmXNmjW0\nadOGJUuWsHLlSkpKSoiLi6NTp074+vq6tQgRkdr+AOKJsyC9e/fkpps68/jjydVa/91332LSpPG8\n/faHBAZrd2ATAAAW0klEQVQG8v33eUyf/hwzZ86p4ZbWnmefHcf27V+yZElmbTfFLc4p1G+77TZu\nvvlmAH766ScCAgLYuHEj0dHRAHTt2pUNGzbg5eVFZGQkvr6++Pr6Ehoayvbt24mIiHBbASJSfe4O\nPk3vm9vatR/y1Ve5td0MOQvnfE7darWSnJzM6tWrmTlzJhs2bHBN0dhsNhwOB06nE7v92A+/22w2\nnE5nldsOCmpg6t+7DQmxV72SnBf1ce24mPvdnbW7a1ve3l7Ur+9b7e3Z7fUAuOQSG8HBdmw2PywW\niyn369Ga6tXzwWr1Nk2N53Wh3NSpU3niiSfo27cvJSUlruVFRUUEBATg7+9PUVFRpeXHh/zpFBQc\nPJ9mebSQEDv5+Y7aboapqY9rz8Xc7+dae1GRk1mzXmDdurV4eVn4y1/+ws6du/ntt/38/e/Tq72d\n/fv38cILz7F583+oX78BgwYN4fDhCg4dKnW17dChQ8yZM4uPPvqQoqIiWrZsxbBhjxMe3hwAh6P4\nf9sqYt68hSxa9DIA119/PaNGpXL33T358ccdzJs3my1bcnA6nTRqFEKPHn+iX78Bpz33vn//PqZP\nf57//nczxcXFXH99CwYOTCAysi1w5Lz9ihUZvP32G+zevQtvbyutWrVm6NDHadr0OgAee+xRmjW7\nHh8fH959901KS8u44447GTIkkblzX+T999/Bz8+P3r0fID7+EQAWLJjLv/+9hri4h3j55dk4HA4i\nI6OYMGE8fn4NASguLqO8/HCl/bdiRQYrV77K3r2/cMUVV/HII3/l1lvvcD2/adN65s+fy44d31O/\nfgNuuqkzjz2WSEBAw2rvrxO560PFOYX6G2+8wd69exk0aBD169fHYrHQunVrsrOzad++PVlZWXTo\n0IGIiAimT59OSUkJpaWl5OXlER4e7paGi4imv83g6aeT2b79K4YMGYbd3pBFi+byyy97iYi4odrb\nOHz4MI8/PpSDB4t48snRGIbB3LkvkZ//q2sdwzBISXmcvLxvefTRIVxySSNWrsxk6NBBLFz4T664\n4spK2+zZsxf5+b+yevX7zJgxhyuuuJKDBw8ydOggrr76WkaPHofVamX16g+YP38O110XTufOXU/Z\nvmeeGYvD8TtPPZWKn58vy5cv5cknE1m58m0CAhqyfPlS5s+fzeDBQ2nWLJyff/6JefPSefbZcSxc\nuNS1nXfeeZM2baJITX2WL774nIUL5/Hpp5tp3foPTJgwldWr32fu3JeIjGxL69ZHTvP+8ssvzJ49\ni0GDhtCgQQPmzHmRhx9+mMWLXz3l9V0LF85j8eIF/OUvD3PDDZFs2rSBceNGY7F4ccstt7F79y5G\nj36Se+65lyFDhvPrr3t58cUXKCkpYfz4SdXeZzXlnEL9jjvu4KmnnuLBBx+kvLycUaNG0bRpU8aM\nGUNaWhphYWHExMTg7e1NfHw8cXFxGIbBiBEj8PPzc3cNIiJ10tdfb2fz5mzGjp3IHXfcCcBNN7Xj\n9ttvP6vtbNq0nry8b5kzZxGtW/8BgKuuupoBA/7iWueTT/5DTs5mXnjhJW68sT0A7dt3JD6+L4sX\nL2DUqNRK22zc+FJCQhrj5eXl2ub27V9xxRVXMX78ZIKCggBo2zaarKx/s2XLf08b6l98sYVHHhno\nev7aa6/j1Vf/yaFDhwgIaMivv+6lX78B9O0bC0BkZFscjkJmzXqBgwcP0qBBAwC8vb2ZMGEyfn71\nuPHG9vzf/72OYVQwatSRb1vdcEMkH3zwLl9+uc0V6ocOHWTChCl06HATAFdffQ39+sWyZs0q7rqr\nR6V2OhwOli5dzIMP9mPgwAQAoqM7cPBgEXPmzOKWW25j+/YvKS0t5S9/eZhGjRoBUL9+A/bu/fms\n9llNOadQb9CgATNmzDhp+dKlS09a1rdvX/r27XsubyMiZ6m2rwDXzMHZ2bZtKwA33dTZtezyyy8n\nIqLNWW1n69bPsdsDXOELcP31zWnS5HLX4//+91Pq1atHZGRbysvLXcujozuwfn1Wtd6nefMWpKfP\np7y8nB9++J5du37k22+/5vDhcsrKSk/7uoiINsyfP4fvvvuWm27qTMeOnRgyZLjr+cTEJwAoKCjg\nxx93sHPnDjZs+Bjgf9s9EupNm16Hn1891+uCgoK56qpQvLyO3BzV19eX+vXr43Acm0r39/d3BTpA\nWNh1XHXVVXz++WcnhXpu7heUlpZw002dK/VRhw438c47b/LTT3to0aIVvr6+PPpoP2699Q46duxE\n585d8fb2jOvAdPMZEZFa4nQ6sFqt+Pv7V1reqFEITmf1z9E7HA4aNgw8aXlw8CWu/y4s/J3i4mJu\nvrnDSetZrdWPgsWLF7B8+RKcTieXXdaE1q0j8Pa2YhjGaV/zzDOTWbRoPh99tJo1a1ZhtVq59dY7\nePLJUfj51WPnzh1MnTqRrVu3UK9ePa67LpwGDWwAHL/ZoyP249WrV++kZce75JJGJy0LDg6msLDw\npOWFhb8DMHhw/1Nua//+ffzhDzcwfXo6S5e+wsqVr7J8+RKCgy8hIWHoSR8SaoNCXUSkljRsGEh5\neTkOh6PSRcSFhb+7Rp/V205DDhz47aTlv//+u+u/bTZ/goKCee656l98d6L33nub+fPnMHJkMrfd\ndqfrw0iPHmc+XRAQ0JDhw0cyfPhIvv32a1atep+MjKVce20YcXEPkZw8goCAhvzjHxlcc00YXl5e\nvP76Cj75ZNM5t/Wo4/vgqP3793P11WEnLbfZjtQzadLzNG7c+KTnQ0OvBo7MPPz979MpLi7m008/\nYdmyfzBlygTatYsmJOTk111I+kEXEZFaEhV15OrvdeuOnbZwOByuafnqb6cdTqeTnJzNrmU//riD\nn37a7XocEdGGAwcKqF+/Ac2bt3T9s2rVe3zwwXun3O6JU8q5uV8QEtKYXr16uwL966+3c+BAAacb\nqB84cIB77/2jq8Zmza5nyJDhXHZZE/bu3cuBAwXs3r2Le+75M2Fh17k+zGRnbwQ44wxAdRw4UMCX\nX25zPf7uu2/ZtWsXUVHtTlq3ZcvWWK1WCgp+q9RH33+fx6JFL2MYBm+99QZ9+txDeXk59erVo3Pn\nrgwcmMDhw4fZty//vNrqDhqpi4jUktDQa7j77p7MnJlGSUkJTZpcTmbm0kpfBa6OG29sT5s2UTzz\nzNMkJAyjfv36vPzybKxWH9c6nTp1oXnzljzxxHD69x/IpZdextq1a/jXv1bwxBNPnXK7/v7+FBcX\n8/HH/6ZFi9Y0b96SN95YyaJFL9OmTRQ7d/7AokUvY7FYKCkpPuU2AgMDueqqUGbMmMahQ4do3PhS\nNm3awC+//EzXrjcTFBTMpZdexooVywkKCsbb25v33nubjRvXA5x2u9VlsVhITR3F4MGPARbmzXuJ\nFi1a0K3bydd7BAUF0bv3A7z44nQcjkJatmzNt99+zbx56XTu3A2bzZ8bbohk//79jBmTzJ//3Iey\nsjIWL15AkyZX0KzZ9efVVndQqIvIBVOTF9LV1YvynnxyNIGBgSxaNI/y8nL69u1LWdnhs9qGxWJh\nypQ0Zs6cxvTpz2O1WomN/Qvr1q11rePt7U1a2ovMnj2T9PSZFBUVcdVVV7m+f34qt94aw/vvv8uY\nMSkMHJhAbGw8u3b9yBtvrOSf/1zMZZddTmxsPDt37mDr1i2nbd+4cc/y0kszmD17JoWFhYSGXs3Y\nsRNcV+E/++xzTJ/+HKmpT2Gz2WjRohXTp6czfHgC27Zt5bLLmpxVfxyvXr16PPLIQGbOnEZJSQmd\nOnVh/PhUystPHX9/+9swgoKCePPNf7FgwVwuuaQRffvG8cgjA4EjU/BTp6axYMFcnn76yO13o6La\nkpo68ayuTagpFuN85zZqgJlvYKEbo9S8i6mPa/tq96qcGLS6Or5qISF27r8/lgYNGpzVzWfkZAsW\nzCUjYymrV39cabkn/o2o1ZvPiEjt8PQQP1Fda6+nKSgoYM+e3VWud/xX2eTiplAXEfFQmzatZ9Kk\n8VWut379pxegNVIXaPr9AvPEaR+zMXMfm23kq+n3k5n5+PUUntjH7pp+11faRERETELT7yK1pKqL\nxsw2KheRmqeRuoiIiElopC5ygWjkLSI1TaEuIh5D32MXOT8KdREPoZG8iJwvnVMXERExCYW6iIiI\nSSjURURETEKhLiIiYhIKdREREZNQqIuIiJiEQl1ERMQk9D11kRqi752LyIWmUBeRWqMPPiLupel3\nERERk1Coi4iImISm30WkzjjVdL1+9EXkGI3URURETEIjdRE30UVfIlLbNFIXERExCYW6iIiISZzT\n9HtZWRmjRo1iz549lJaWkpCQwHXXXUdKSgoWi4VmzZqRmpqKl5cXmZmZZGRkYLVaSUhIoHv37u6u\nQURERDjHUH/zzTcJDAzkueee48CBA/Tq1YvmzZuTmJhI+/btGTt2LGvWrKFNmzYsWbKElStXUlJS\nQlxcHJ06dcLX19fddYiICek6BZGzc06hfueddxITEwOAYRh4e3uTm5tLdHQ0AF27dmXDhg14eXkR\nGRmJr68vvr6+hIaGsn37diIiItxXgYhc1E4Mfn3FTS5m5xTqNpsNAKfTybBhw0hMTGTq1KlYLBbX\n8w6HA6fTid1ur/Q6p9NZ5faDghpgtXqfS9PqhJAQe9UryXlRH1+8zLDvzVCDpzNrH5/zV9p+/vln\nhgwZQlxcHD179uS5555zPVdUVERAQAD+/v4UFRVVWn58yJ9OQcHBc22WxwsJsZOf76jtZpia+vji\nVtf3vY7fmueJfeyuDxnndPX7vn376N+/P0lJSfTu3RuAli1bkp2dDUBWVhbt2rUjIiKCnJwcSkpK\ncDgc5OXlER4e7paGi4iISGXnNFKfM2cOhYWFpKenk56eDsDo0aOZOHEiaWlphIWFERMTg7e3N/Hx\n8cTFxWEYBiNGjMDPz8+tBYiIiMgRFsMwjNpuxIk8bVrEnTxx2sdsaquPdaW2Z6jrF8rpb0TN88Q+\nrtXpdxEREfE8uve7SDXpq1Mi4ukU6iJiavowJhcThbrIaegced2k/SYXM51TFxERMQmFuoiIiEko\n1EVERExC59RFzpHO3YqIp9FIXURExCQ0UheRi5q+8iZmolAXkYuKTpuImWn6XURExCQ0Uhf5H43g\nRKSu00hdRETEJBTqIiIiJqFQFxERMQmFuoiIiEko1EVERExCoS4iImISCnURERGTUKiLiIiYhEJd\nRETEJHRHOblo6Q5yImI2GqmLiIiYhEbqIiLH0U+xSl2mkbqIiIhJKNRFRERMQqEuIiJiEjqnLqal\nc6MicrHRSF1ERMQkNFIXETkPmhEST6JQlzpLf0xFRCo7r1D//PPPef7551myZAk7d+4kJSUFi8VC\ns2bNSE1NxcvLi8zMTDIyMrBarSQkJNC9e3d3tV3krOgOciJidud8Tv3ll1/m6aefpqSkBIDJkyeT\nmJjIsmXLMAyDNWvWkJ+fz5IlS8jIyGDBggWkpaVRWlrqtsaLiIjIMec8Ug8NDWXWrFk8+eSTAOTm\n5hIdHQ1A165d2bBhA15eXkRGRuLr64uvry+hoaFs376diIgI97ReROQC04yPeLJzDvWYmBh2797t\nemwYBhaLBQCbzYbD4cDpdGK3213r2Gw2nE5nldsOCmqA1ep9rk3zeCEh9qpXkrN2fL+qj8VdzvZY\ncsexp+O35pm1j912oZyX17GZ/KKiIgICAvD396eoqKjS8uND/nQKCg66q1keJyTETn6+o7abYUpH\n+1V9LO7Uc+T/ndX653vs6fiteZ7Yx+76kOG2UG/ZsiXZ2dm0b9+erKwsOnToQEREBNOnT6ekpITS\n0lLy8vIIDw9311uKVKJpURG52Lkt1JOTkxkzZgxpaWmEhYURExODt7c38fHxxMXFYRgGI0aMwM/P\nz11vKSIiIsexGIZh1HYjTuRp0yLu5InTPnWVRubiic73fgn6G1HzPLGP3TX9rtvEioiImIRCXURE\nxCR0m1gRkRqk2xnLhaRQlzpD59ClLtBxKrVJ0+8iIiImoZG6iEgtOnFk/9a0P9VSS8QMNFIXEREx\nCYW6iIiISSjURURETELn1MVj6Ks/IiLnR6EuHktfDRIROTsKdakVCmwREffTOXURERGT0EhdROQC\n0iyV1CSN1EVERExCoS4iImISCnURERGT0Dl1EREP0nPk/1V6fOL9GnQ/BzkThbqck7P9w6KLg0RE\nap5CXWqEQlxE5MJTqIuIeDB9QJazoVAXETEZnXe/eOnqdxEREZPQSF3cQlOEIiK1T6EuIlKH6QO1\nHE/T7yIiIiahkbpUi0YDIiKeT6Eup6QQFzGPqv5/1tXx5qFQF0AhLnIx01fgzEOhfpFSiIvIudKH\nAM+lC+VERERMQiN1ERGpRCPxuqvGQ72iooJx48bx9ddf4+vry8SJE7n66qtr+m1Nr6r/6TS9LiLu\ncrZ/T/ShoPbUeKh/+OGHlJaW8uqrr7JlyxamTJnC7Nmza/ptPZ67D3qFuIjUVfoQ4D4WwzCMmnyD\nyZMnExERwR//+EcAunTpwscff3zG1+TnO9zahrP9Osf5ru+OUbNG3iIiF05VHyRq+oNHSIjdLdup\n8VAfPXo0d9xxB926dQPg5ptv5sMPP8Rq1el8ERERd6rxq9/9/f0pKipyPa6oqFCgi4iI1IAaD/Wo\nqCiysrIA2LJlC+Hh4TX9liIiIhelGp9+P3r1+zfffINhGEyaNImmTZvW5FuKiIhclGo81EVEROTC\n0B3lRERETEKhLiIiYhIK9RpSUVHB2LFjuf/++4mPj2fnzp2Vnn/llVf44x//SHx8PPHx8Xz//fe1\n1NK67fPPPyc+Pv6k5R999BH33Xcf999/P5mZmbXQMnM4Xf/q+D1/ZWVlJCUlERcXR+/evVmzZk2l\n53UMn5+q+te0x7AhNeKDDz4wkpOTDcMwjM8++8wYPHhwpedHjhxpfPHFF7XRNNOYN2+e0aNHD6NP\nnz6VlpeWlhq33XabceDAAaOkpMS49957jfz8/FpqZd11uv41DB2/7vDaa68ZEydONAzDMAoKCoxu\n3bq5ntMxfP7O1L+GYd5jWCP1GpKTk0OXLl0AaNOmDdu2bav0fG5uLvPmzSM2Npa5c+fWRhPrvNDQ\nUGbNmnXS8ry8PEJDQ2nYsCG+vr60bduWzZs310IL67bT9S/o+HWHO++8k+HDhwNgGAbe3t6u53QM\nn78z9S+Y9xhWqNcQp9OJv7+/67G3tzfl5eWux3/84x8ZN24cixcvJicnh7Vr19ZGM+u0mJiYU97I\nyOl0Yrcfu+WizWbD6XReyKaZwun6F3T8uoPNZsPf3x+n08mwYcNITEx0Padj+PydqX/BvMewQr2G\nnOlOeoZh0K9fP4KDg/H19aVbt258+eWXtdVU0zmx74uKiir9gZTzo+PXfX7++Wceeugh/vSnP9Gz\nZ0/Xch3D7nG6/jXzMaxQryFnupOe0+mkR48eFBUVYRgG2dnZtG7duraaajpNmzZl586dHDhwgNLS\nUj799FMiIyNru1mmoePXPfbt20f//v1JSkqid+/elZ7TMXz+ztS/Zj6GdRP2GnL77bezYcMGHnjg\nAded9N566y0OHjzI/fffz4gRI3jooYfw9fWlY8eOrh+8kXN3fP+mpKQwYMAADMPgvvvu49JLL63t\n5tV5On7da86cORQWFpKenk56ejoAffr04dChQzqG3aCq/jXrMaw7yomIiJiEpt9FRERMQqEuIiJi\nEgp1ERERk1Coi4iImIRCXURExCQU6iIiIiahUBcRETEJhbqIiIhJ/D/jezLLxpWh9wAAAABJRU5E\nrkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x137a0c240>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.subplot(211)\n",
    "plt.hist(delta_sample, bins=100, label='delta samples')\n",
    "plt.legend(loc='upper right', fontsize=16)\n",
    "plt.grid(True)\n",
    "plt.xlim(-1000, 1000)\n",
    "\n",
    "plt.subplot(212)\n",
    "plt.hist(q_delta_sample, bins=100, label='q_delta samples')\n",
    "plt.legend(loc='upper right', fontsize=16)\n",
    "plt.grid(True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mean: 1.6855731010437012\n",
      "Variance: 0.0\n"
     ]
    }
   ],
   "source": [
    "print(('Mean: {}').format(q_delta.eval().mean()))\n",
    "print(('Variance: {}').format(q_delta.eval().std()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "delta mu: [ 1.43979084], delta scale: [-1.06268966]\n"
     ]
    }
   ],
   "source": [
    "print(('delta mu: {}, delta scale: {}').format(q_delta.parameters['loc'].eval(), q_delta.parameters['scale'].eval()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Comparison results for different <b>number of samples</b> for variation inference with Edward."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "|Num Samples| Loss    | Time   |\n",
    "|-----------|--------:| ------:|\n",
    "| 200       |55458.379| 2296s  |\n",
    "| 100       |33341.516| 1182s  |\n",
    "| 50        |49572.426| 600s   |\n",
    "| 20        |44003.395| 227s   |\n",
    "| 15        |46201.785| 169s   |\n",
    "| 10        |47803.379| 122s   |\n",
    "| 5         |56676.551| 68s    |\n",
    "| 1         |63954.184| 24s    |"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Comparison results for different <b>optimization algorithms</b> for KL divergence. Number of samples:10, number of iterations: 10 000"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "|Optimization method              | Description | Loss    | Time |\n",
    "|---------------------------------|--------:| ------:|------:|\n",
    "| KLqp                            |minimizes the objective by automatically selecting from a variety of black box inference techniques|35386.125|128s|\n",
    "| ReparameterizationKLqp          |minimizes the objective using the reparameterization gradient|43876.039|240s|\n",
    "| ReparameterizationEntropyKLqp   |minimizes the objective using the reparameterization gradient and an analytic entropy term|54193.547|184s|\n",
    "| ReparameterizationKLKLqp        |minimizes the objective using the reparameterization gradient and an analytic KL term|-|-|\n",
    "| ScoreEntropyKLqp                |minimizes the objective using the score function gradient and an analytic entropy term|38887.348|127s|\n",
    "| ScoreKLKLqp                     |minimizes the objective using the score function gradient and an analytic KL term|-|-|\n",
    "| ScoreKLqp                       |minimizes the objective using the score function gradient|43540.645|123s|\n",
    "\n",
    "'-' means 'No KL(dist_a || dist_b) registered for dist_a type NormalWithSoftplusScale and dist_b type StudentT' (Is automatically computed when p(z) and q(z; lambda) are Normal)"
   ]
  }
 ],
 "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.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}