Using tensorflow_probability's HMC sampler with pymc4

hmc.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import tensorflow as tf\n",
    "import tensorflow_probability as tfp\n",
    "import numpy as np\n",
    "tfd = tf.contrib.distributions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "class HMC(object):\n",
    "    def __init__(self, model,  target_accept_rate=0.651, num_leapfrog_steps=3):\n",
    "#         tf.reset_default_graph()\n",
    "        self.target_accept_rate = target_accept_rate\n",
    "        # Tuning acceptance rates:\n",
    "        self.dtype = np.float32\n",
    "        \n",
    "        self.x = {}\n",
    "        self.step_size = {}        \n",
    "        for i in model.named_vars:\n",
    "            self.x[i] = tf.get_variable(name='x', initializer=self.dtype(1))#, reuse=tf.AUTO_REUSE)\n",
    "            self.step_size[i] = tf.get_variable(name='step_size', initializer=self.dtype(1))\n",
    "    \n",
    "    def sample(self, model, draws, tune):\n",
    "        print(\"Intializing HMC Sampler...\")\n",
    "        dictionary = {}\n",
    "        for i in model.named_vars:\n",
    "            print(model.named_vars[i].distribution.log_prob)\n",
    "            dictionary[i] = tfp.mcmc.HamiltonianMonteCarlo(\n",
    "              target_log_prob_fn=model.named_vars[i].distribution.log_prob,\n",
    "              step_size=self.step_size[i],\n",
    "              num_leapfrog_steps=3)\n",
    "        # One iteration of the HMC\n",
    "            next_x = {}\n",
    "            other_results = {}\n",
    "            \n",
    "            next_x[i], other_results[i] = dictionary[i].one_step(\n",
    "              current_state=self.x[i],\n",
    "              previous_kernel_results=dictionary[i].bootstrap_results(self.x[i]))\n",
    "            x_update = {}\n",
    "            x_update[i] = self.x[i].assign(next_x[i])\n",
    "\n",
    "            step_size_update = {}\n",
    "            step_size_update[i] = self.step_size[i].assign_add(\n",
    "              self.step_size[i] * tf.where(\n",
    "                  tf.exp(tf.minimum(other_results[i].log_accept_ratio, 0.)) >\n",
    "                      self.target_accept_rate,\n",
    "                  0.01, -0.01))\n",
    "        # Note, the adaptations are performed during warmup only.\n",
    "            warmup = {}\n",
    "            warmup[i] = tf.group([x_update[i], step_size_update[i]])\n",
    "        init = tf.global_variables_initializer()\n",
    "        with tf.Session() as sess:\n",
    "            sess.run(init)\n",
    "            # Warm up the sampler and adapt the step size\n",
    "            for _ in range(tune):\n",
    "                sess.run(warmup)\n",
    "            # Collect samples without adapting step size\n",
    "#             for i in model.named_vars:\n",
    "#                 samples[i] = np.zeros([draws])\n",
    "            samples = []\n",
    "            for j in range(draws):\n",
    "                _, x_,= sess.run([x_update, self.x])\n",
    "                samples.append(x_)\n",
    "#                 print(x_)\n",
    "        return samples\n",
    "#     print(samples.me~~an(), samples.std())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pymc4 as pm\n",
    "from tensorflow_probability import edward2 as ed\n",
    "from tensorflow_probability import distributions as tfd\n",
    "from pymc4 import RandomVariable as RV"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "with pm.Model() as model:\n",
    "    x = RV(\"x\", tfd.Normal(0.,1.))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "hmc = HMC(model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<bound method Distribution.log_prob of <tf.distributions.Normal 'Normal' batch_shape=() event_shape=() dtype=float32>>\n"
     ]
    }
   ],
   "source": [
    "for i in model.named_vars:\n",
    "    print(model.named_vars[i].distribution.log_prob)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "trace = hmc.sample(model, draws=1000, tune=500)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/sharan/anaconda3/envs/pymc3/lib/python3.6/site-packages/matplotlib/axes/_axes.py:6462: UserWarning: The 'normed' kwarg is deprecated, and has been replaced by the 'density' kwarg.\n",
      "  warnings.warn(\"The 'normed' kwarg is deprecated, and has been \"\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xd4XFed//H3d0a9V0u2iiW3xHKJi2wnMTGBFOIADoF0lrZZsssmu+xvl4e6y283LD9YeB7KLgE2sAkLJCQhCek9pDqJbdlxr3KVZFu993J+f0g2inEZyzO6M6PP68EPmtG15oPRfHR07rnnmnMOERGJLj6vA4iISPCp3EVEopDKXUQkCqncRUSikMpdRCQKqdxFRKKQyl1EJAqp3EVEopDKXUQkCsV49cI5OTmupKTEq5cXEYlI69evb3DO5Z7pOM/KvaSkhIqKCq9eXkQkIpnZwUCO07SMiEgUUrmLiEQhlbuISBRSuYuIRCGVu4hIFFK5i4hEIZW7iEgUUrmLiEQhlbuISBTy7ApVkVC5f82hgI+9ZVlxCJOIeEcjdxGRKKRyFxGJQip3EZEopHIXEYlCOqEqE5pOvkq00shdRCQKqdxFRKKQyl1EJAqp3EVEolBA5W5mV5nZLjOrNLOvnua4T5iZM7Py4EUUEZGzdcZyNzM/cBewEigDbjazspMclwp8EVgT7JAiInJ2Ahm5LwUqnXP7nHN9wAPANSc57lvAfwA9QcwnIiJjEEi5FwBVox5Xjzx3nJktAoqcc08HMZuIiIzROV/EZGY+4AfAZwM49jbgNoDiYl0QIt4bco6+gSGGhhyDzpEcH4PPzOtYIucskHKvAYpGPS4cee6YVGAu8KoNvynygSfMbJVzrmL0F3LO3Q3cDVBeXu7OIbfIOTvS2s1v3zlIc1f/8edyU+O5YXERBZmJHiYTOXeBlPs6YKaZlTJc6jcBtxz7pHOuFcg59tjMXgW+dGKxi4STnUfaeKCiioQYHyvn5uP3GUNDjjcrG/jZa5V88Pw83j8rF79Po3iJTGcsd+fcgJndATwP+IF7nHPbzOxOoMI590SoQ4oE01t7G3h68xEmZyTwqQtLSE+MPf65xVOzeGJTDS/tqOVgYyefvqhEBS8RKaA5d+fcM8AzJzz3zVMce+m5xxIJjYONnTy1+Qiz81O5cUkxcTHvXVOQGOfnxiXFlOY08djGGh7fWMO1CwswzcNLhNGukDJhDA45Ht94mPTEWG5YUvRnxT7a0tIsWrv7eGVXPTkp8ayYlTuOSUXOncpdJoy39jZwtK2Hv1hWTHyM/4zHXzY7j4aOPp7fdpTslLhxSCgSPNpbRiaElq4+Xt5Rx/n5qcyenBbQ3/GZcd3iQgozE3l4fTXVzV0hTikSPCp3mRCe2nwEh+OjF0w5q/nzWL+PG5cU44CvPLKZoSGt4JXIoHKXqHeoqYvtR9r4wHmTyEw6++mVrOQ4Vs7NZ3VlI/etORiChCLBpzl3iXqv764nMdbPRdOzx/w1lpZksf1wG3c+tZ2Wrn6yU+JPeaxuxyfhQCN3iWp17T3sONLGhdOyAzqJeipmxrULC/CZ8ciGGpzT9IyEN5W7RLU39jQQ47dzGrUfk5EUx9VzJ3OgsZPN1a1BSCcSOip3iVqt3f1sPNTC4qmZpMQHZwZycUkmUzISeG7bUfoGhoLyNUVCQeUuUeutygYcjvfNCN4FSD4zPjp/Cq3d/by2uz5oX1ck2FTuEpV6+gdZe6CJuQXpZCUH9wKkqdnJzC9M54099TR39gX1a4sEi8pdotL6g830DgxxyczQbBuwcu5kzODZrUdC8vVFzpWWQkrQ3b/mUMDHhmLZ4JBzvL2vkalZSRRkhGZf9vTEWFbMzOXlnXUcbulmSoheR2SsNHKXqLP7aDtNnX1BWSFzOstn5JAQ6+OVXXUhfR2RsVC5S9R5e18jaQkxzJmSHtLXSYj1c9G0HLYdbqO2TfeFl/CicpeoUlnXzp66DpZNyx6Xm2wsn55NXIyPVzV6lzCjcpeo8r9vHSTGZywpyRqX10uKj+HC0iw2V7fS0NE7Lq8pEgiVu0SN1u5+HtlQzfzCjKBdtBSI5TNy8PuM13Zp3buED5W7RI3frT1EV98gF4f4ROqJUhNiWVKaxbtVzbR194/ra4ucispdokLfwBD3rt7P8hnZnixLvHhaNkMO1h9qHvfXFjkZlbtEhSc2Haa2rZfbVkz35PWzU+KZlpNMxYEm3dBDwoLKXSKec45fvL6P8/NTWTEzx7McS0qyaO7qZ/XeBs8yiByjcpeI9+ruenbVtvP5S6ad1S30gq1sShpJcX4eWFvlWQaRY1TuEvHufm0f+WkJfPSCKZ7miPX7WFiUwQvbj9KoZZHiMZW7RLSNVS28va+Rzy0vIS7G+2/nJSVZ9A86HtlQ7XUUmeC8fzeIjJFzju88s4Ps5LiwuW/ppLQEyqdm8sC6Kt2KTzylcpeI9fKOOtbsb+IfLp9JakKs13GOu3FJEfvqO9lY1eJ1FJnAVO4SkQYGh/jOszuYlpPMTUvDY9R+zJVz8on1G89tPep1FJnAVO4SkR6sqGJvfSdfWXk+sf7w+jZOT4xl+Ywcnt16VFMz4pnweleIBKCzd4AfvriHJSWZXFmW53Wck7pqTj6HmrrYfqTN6ygyQancJeL8+9M7aOjo5etXz/Z0XfvpXFGWh8/Q1Ix4RuUuEeXJTYf53dpDfOHS6SwszvQ6zillp8SzrDSbZ1Xu4hGVu0SMg42dfO3RLSwqzuAfr5jldZwzWjkvn8q6Dirr2r2OIhOQyl0iQt/AEH/3u3fxGfznzQvD7iTqyXxoTj4Az27R6F3GX/i/Q2TC6xsY4osPvMvm6la+d90FFGYmeR0pIHlpCSyemqmpGfGEyl3CWk//IH/9mwqe3XqUf/7wbK6am+91pLOycm4+24+0caixy+soMsGo3CVsdfQO8Ll71/Hq7nq+8/F5/NUl07yOdNYunz28VPPV3bqBtoyv8bvRpMhZ2F3bzt/et4H9DZ388IYFfGxhAfevOeR1rLNWkpPM1OwkXttVz6cvKvE6jkwgGrlL2Hl4fTXX/GQ1LV19/Povl/KxhQVeRzon75+Vy1t7G+kdGPQ6ikwgKncJGwODQ/zLY1v50u83Mb8wnWf+/hKWz/DuzkrB8v5ZuXT3D7L+gO6vKuMnoGkZM7sK+DHgB37pnPvuCZ//G+B2YBDoAG5zzm0PclaJYl19A/zd/e/y8s46blsxjS9/6DxiImC548mcOH3UOzCI34yfvbaXAyecWA2XrYol+pzx3WNmfuAuYCVQBtxsZmUnHHa/c26ec24B8D3gB0FPKlGrrr2HG//7HV7ZVce3PjaXr189O2KL/WTiY/xMzUliT22H11FkAgnkHbQUqHTO7XPO9QEPANeMPsA5N3p3pGRAW+FJQHoHBrn1VxVU1nXwi0+X86kLp3odKSRmTUrlaFsPrd39XkeRCSKQci8ARt/xt3rkufcws9vNbC/DI/e/P9kXMrPbzKzCzCrq6+vHkleizHef3cmWmlZ+dNMCLpsdnjs8BsOsvFQA9tRqKwIZH0H73dc5d5dzbjrwFeCfT3HM3c65cudceW5ubrBeWiLU9sNt3Lv6AJ+9uOT4pfrRKi8tnrSEGHbXaWpGxkcg5V4DFI16XDjy3Kk8AHzsXEJJ9Gvp6uORDdXMLUjja1ef73WckDMzZk5KpbKuncEhzVpK6AVS7uuAmWZWamZxwE3AE6MPMLOZox5+GNgTvIgSjR7bWMOQc/zk5kXEx/i9jjMuZual0NM/RHWztiKQ0DvjUkjn3ICZ3QE8z/BSyHucc9vM7E6gwjn3BHCHmV0O9APNwGdCGVoiW01LN7trO7iyLI+SnGSv44ybGbkpGFBZ38HU7Inzv1u8EdA6d+fcM8AzJzz3zVEffzHIuSSKvba7nvgYHxdOy/Y6yrhKio9hckYCe+s6uSz6Z6LEY9GzmFgiQkNHL9tqWrlwWjYJsRNjOma06bkpVDV10Tcw5HUUiXIqdxlXr++ux+8zLp4+sUbtx8zITWHQOQ40dnodRaKcyl3GTWt3P+8eamHx1ExSE2K9juOJqdnJ+H3GXi2JlBBTucu4WV3ZgMNxycyJe41DXIyP4qwkKutV7hJaKncZF4NDjg2Hmimbkk5WcpzXcTw1PTeFI609dPYOeB1FopjKXcbFvvoOuvoGWVCY7nUUz82YlALAXo3eJYRU7jIuttS0EhfjY+bIHisTWUFGIvExPvbW66SqhI7KXUJucMix7XAbs/NTiY2irXzHyu8zpuUka+QuIaV3moTcvvoOuvsHmVeQ4XWUsDF9UgpNnX1UNWkrAgkNlbuE3JaaVuJjfMzMS/E6StiYnjv8b7G6ssHjJBKtVO4SUsenZCanaUpmlEmpw1sAv6FylxDRu01Cau/xKRmtkhnNzJgxKZXVlQ3aAlhCQuUuIbV1ZErm2PI/+ZMZk1Jo6epna02r11EkCqncJWSGnGPHkTbO1yqZkzr2A+9NTc1ICOgdJyFztLWHzr7B4/cPlfdKiY9hzpQ0Xt+t+wlL8KncJWQqRzbHmq4pmVN638wcNhxqpkNbEUiQqdwlZCrrOkZuDD0xd4AMxIqZufQPOtbsa/Q6ikQZlbuERP/gEAcaO5mRq1H76SyemklCrI839mjeXYJL5S4hcaCxk4Ehp1UyZ5AQ62dpaTZv7NG8uwRXQPdQFTlbe+s68JtRmnP6cr9/zaFxShS+VszM4d+f3sHhlm6mZCR6HUeihEbuEhKVdR0UZycRF6NvsTM5dvMSrZqRYNI7T4Kuo3eAw609mpIJ0Ky8FAoyEnlpR53XUSSKqNwl6I5tZauTqYExM64oy+ONPfV09WlJpASHyl2CrrKug8RYPwWZmj8O1BVlefQODGnVjASNyl2CyjnH3roOpuUm4zPzOk7EWFqaRVpCDC9sq/U6ikQJlbsEVXVzNy3d/UzTlMxZifX7+OD5k/jjzloGBoe8jiNRQOUuQfXOyJWWpTnJHieJPFfOyae5q5+Kg81eR5EooHKXoFqzv4mkOD+TUuO9jhJxVszKJc7v48XtmpqRc6dyl6Bas7+RkmzNt49FSnwMy2dk88L2ozinG3jIuVG5S9AcbummqqlbUzLn4IqyfKqautlV2+51FIlwKncJmrX7mwDNt5+Ly8smYQbPbDnqdRSJcCp3CZo1+xtJTYghPz3B6ygRa1JqAsun5/DI+mqGdG9VOQcqdwmaNfuaWFKSpfn2c3TDkiJqWrpZvVcXNMnYqdwlKOrae9jX0Mmy0iyvo0S8K8vySE+M5cF1VV5HkQimcpegODbfvmxatsdJIl9CrJ9rFxbwwrZamjv7vI4jEUrlLkGxZt/w+vY5U9K8jhIVblxSRN/gEH94t8brKBKhVO4SFGv2N7J4aiaxfn1LBcPsyWlcUJjOQxVVWvMuY6J3opyzxo5edtd2cKGmZILqhiVF7DzazubqVq+jSARSucs5OzbffuE0nUwNplUXTCEx1s//vn3A6ygSgQIqdzO7ysx2mVmlmX31JJ//RzPbbmabzexlM5sa/KgSrt7Z10hirJ95BRleR4kqqQmx3LKsmMferTl+AxSRQJ2x3M3MD9wFrATKgJvNrOyEw94Fyp1z84GHge8FO6iEr3f2NVFekqn7pYbAFy6dTnyMnx+/tMfrKBJhAnk3LgUqnXP7nHN9wAPANaMPcM694pzrGnn4DlAY3JgSrho7etlV26759hDJSYnnc8tLeHLzYXYebfM6jkSQQMq9ABh9NUX1yHOncivw7Mk+YWa3mVmFmVXU1+tO79HgT/PtKvdQuW3FNFLiYvjhi7u9jiIRJKi/R5vZXwDlwPdP9nnn3N3OuXLnXHlubm4wX1o8cmy+fX5hutdRolZGUhy3XlLK89tq2aKVMxKgQMq9Biga9bhw5Ln3MLPLgW8Aq5xzvcGJJ+Hu2Hy71reH1l++r5SMpFi+9fR2bSgmAQnkHbkOmGlmpWYWB9wEPDH6ADNbCPw3w8VeF/yYEo403z5+0hJi+frK2azd38Sv3z7gdRyJAGcsd+fcAHAH8DywA3jIObfNzO40s1Ujh30fSAF+b2YbzeyJU3w5iSJrRubbL5quch8P15cX8oHzcvnuczs50NDpdRwJcwH9Lu2ce8Y5N8s5N9059+2R577pnHti5OPLnXN5zrkFI39Wnf4rSjR4Z18jSXF+5hVovn08mBnf+fh84vw+vvT7TQxqekZOQxOlMmbv7GukvCRL8+3jKD89gX9dNYeKg83c8+Z+r+NIGIvxOoBEpvr24f1kPrbwdKti5UzuX3Mo4GNvWVYMwLULC3hu61G+9/xOFk3NZPHUzFDFkwimIZeMyeu7h69TWDFTS1rHm5nx/esvYHJ6Irfft4GGDi1Okz+nkbuMyWu768lJiadssvZvHy8njvJXXTCFn7+2lxt+/jafW16K3zd8e8NjI3yZ2DRyl7M2OOR4fU89K2bl4PPpfqlemZKRyDULCtjX0MlLO2q9jiNhRuUuZ21zdQstXf1cet4kr6NMeIunZrKkJJPXdtez62i713EkjKjc5ay9uqsen8ElM3K8jiLAR+ZPIT8tgd+vr6K1u9/rOBImVO5y1l7bXc8FRRlkJsd5HUWAWL+Pm5cWMzDoeHDdIQYGh7yOJGFA5S5npamzj03VLbx/llbJhJPc1HiuWTCFA41d/Phl7f0uKnc5S2/sqcc5NN8ehhYWZ7KoOJO7XqnU7pGicpez89quejKTYrXlQJj68LzJ5KTE89VHN2t6ZoJTuUvAhoYcr+2uZ8Ws3ONrqiW8JMb5+bdVc9h2uI17Vmt7golMFzFJwDZWt9DY2af59jDX1NnH7PxUvv/8LvoGHFmnOfGtC56il0buErAnNx0mzu/j8rI8r6PIaZgZqxYU4DPj8Y01OKfdIycilbsEZHDI8fTmI1x6Xi5pCbFex5EzSE+M5fLZeeyp66CyrsPrOOIBlbsEZM3+Rurae1m1YIrXUSRAy0qzyEyK5fntRxnS6H3CUblLQJ7cdISkOD+Xna8pmUgR4/dx2ew8Drf0sO1wm9dxZJzphKqcUf/gEM9uPcLMSSn84d0/uze6hLEFRRm8vrueF7cfpWxymlY5TSAaucsZvbmngZaufuYXZngdRc6Sz4wry/Jp6Ohjw6Fmr+PIOFK5yxk9uekwaQkxzMxL8TqKjMHsyakUZSbyx511urBpAlG5y2n19A/ywvZaVs6dTIxP3y6RyMy4bHYerd39bNK2BBOG3q1yWs9uPUJH7wDXaJVMRJs5KYX8tISRvYG0cmYi0AnVCe50N2h2zvHTV/eSmxLP/oZOzHQyLlKZGZfMzOH366vZU9fBrLxUryNJiGnkLqdU1dRFTUs3F03PVrFHgXmF6aQlxPDGnnqvo8g4ULnLKb21r5GEWB8Li7VKJhrE+HxcPD2HvfWd1LR0ex1HQkzlLifV2t3P1ppWyqdmER/j9zqOBMnS0iziY3y8qdF71FO5y0mt2d+Ic3DhtGyvo0gQJcT6WVKSxZaaVt1vNcqp3OXP9A8OsXZ/E+dPTjvtdrESmS6alo1zsO5Ak9dRJIRU7vJnNhxqpqtvkIuna9QejTKT45iVl8q6A03066KmqKVyl/foHxzilZ11TM1KYlpOstdxJESWlWbR3jPAS9trvY4iIaJyl/dYu7+Jtp4BrijL0/LHKDYrP5WMpFh+u+ag11EkRFTuclzvwCCv7q5nem4y03K1j0w085mxtCSL1ZWN7KvXzTyikcpdjntnbyOdvQNcMVt7tk8Ei6dmEus37jvNVcoSuVTuAgxvEPb6ngbOy0ulOFtz7RNBakIsH5qTz8Prq+nuG/Q6jgSZyl0AeH13Pd39g7r59QTzyWVTae3u55ktR7yOIkGmchdauvp4s7KBCwrTKchI9DqOjKMLp2VRmpPMg+uqvI4iQaZyF14cWQ535Zx8j5PIeDMzblxSxNoDTVTW6cRqNFG5T3A1zd28W9XC8hk5ZCbpatSJ6BOLConxGQ+u04nVaKJyn8Ccczyz9QjJcX7ePyvX6zjikdzUeK4oy+ORDTX0DujEarQIqNzN7Coz22VmlWb21ZN8foWZbTCzATO7LvgxJRRe3F7L/oZOLpudR0Ksdn6cyG5aWkxTZ9/xKTqJfGcsdzPzA3cBK4Ey4GYzKzvhsEPAZ4H7gx1QQqNvYIjvPLuT3JR4lpRkeR1HPHbJjBwKMhJ5YK1OrEaLQEbuS4FK59w+51wf8ABwzegDnHMHnHObAe1CFCF+885B9jd0cvW8fPw+bTMw0fl8wydW36xs4FBjl9dxJAgCKfcCYPSP8+qR5yRCtXT18Z8v7+GSmTm6l6Ycd315IT6DB3RiNSqM6wlVM7vNzCrMrKK+XneC8cqPXtpDe08///zhMm0OJsdNTk/kg+fn8VBFtbYCjgKBlHsNUDTqceHIc2fNOXe3c67cOVeem6vVGV7YW9/Bb985yE1LizkvX6N2ea9PLiumoaNXWwFHgUDKfR0w08xKzSwOuAl4IrSxJFS+88wOEmL9/J/LZ3kdRcLQilm5FGQkcv9aTc1EujOWu3NuALgDeB7YATzknNtmZnea2SoAM1tiZtXA9cB/m9m2UIaWsXlzTwMv7ajjjg/OIDc13us4Eob8IydW39jTwMHGTq/jyDkIaM7dOfeMc26Wc266c+7bI8990zn3xMjH65xzhc65ZOdctnNuTihDy9kbGBziW09tpygrkc8tL/E6joSxG5cU4fcZv9OyyIimK1QniAcrqthV287XV84mPkYXLMmp5aUlcNn5k3h4fRV9AzqxGqlU7hNAW08/P3hhN0tLs7hqrjYHkzO7ZVkxDR19PLftqNdRZIxU7hPAT/5YSVNXH/+ipY8SoBUzcynJTuJXq/d7HUXGSOUe5XbXtnPPm/u5YXER8wrTvY4jEcLnMz5zcQkbDrWwqarF6zgyBir3KOac458f20pKQgxfWXm+13Ekwly3uJCU+Bh+9dYBr6PIGKjco9hjG2tYu7+Jr1x1PlnJ2qtdzk5qQizXLS7kqc2HqWvr8TqOnKUYrwNIaLR29/Ptp3ewoCiDG8uLzvwXZEK6f83pL1bKTo5jYNDx1Ue3cM9nl4xTKgkGjdyj1Pee20lTZx///rG5+LTro4xRdko8s/JSWbO/STfyiDAauUeh//v4Nu5bc4jl07PZXN3K5upWryNJBLt4Rjb3rj7A4xsPc4N+C4wYGrlHmcaOXh7dUE1eWrxueC1BMSM3hcnpCfz81b0MDjmv40iAVO5RxDnH1x7dQlf/IDeUFxHr1/+9cu7MjEvPm8S+hk6e3nLE6zgSIL37o8hDFVW8sL2WK8vymJye6HUciSJzpqQxY1IKd/2xkiGN3iOCyj1KrNnXyL88to2Lp2ezfEaO13EkyvjMuOMDM9hV284L2us9IuiEahTYU9vO539dQVFWIj/95CKe2aL9QCT42nsGyE6O486nttHY0XvKrSxuWVY8zsnkZDRyj3C1bT189t51xMf6+dXnlpKRpIuVJDT8PuP9s3I53NLDrtp2r+PIGajcI1h1cxef+p81NHf1ce9nl1CUleR1JIlyC4szyUqO47mtR7VyJsyp3CPUpqoWPnbXWxxp7eGXnylnboE2BZPQ8/uMq+bkU9feS8XBJq/jyGmo3CPQ05uPcOPdb5MQ6+PRL1zMxdN1AlXGz5wpaZTmJPPS9lp6+nXVarjSCVWPnWlvj9GunpfPvz25nT+8W8PC4gx+8elyclJ0L1QZX2bG1fMm89NXKnllVx0r5072OpKchMo9Ajjn2H6kjR+8uJuWrj7+/rKZ3P6B6bpdnnimICORhcWZvLW3kWWl2dp1NAxpWibM1bX18Ku3DnDfmkNMSo3n8TuW849XzFKxi+euLMvDb8bjG2twTidXw41G7mGqp3+Ql3fU8va+RuJifHx43mR+dNMCbSkgYSMtMZYPzc3nyU2HqTjQzJLSLK8jySgq9zDjnGNzTSvPbDlCR88A5SVZXFGWR0p8jIpdws6y0iy2H27l6a1HmD4pRdMzYURtEUaau/q4d/UBHlxXRVpCLF+4dDrXLiwgJV4/gyU8+cz4xKJCDHhkQzVDmp4JGyr3MOCcY8OhZv7z5T1UNXfx0Qum8IVLp1OYqYuSJPxlJMXxkfmT2d/QyerKBq/jyAgNCT3W3TfIo+9Ws+1wGyXZSVy/uIhM/WorEWZRcSY7j7bz3NajvLKzjg+cP8nrSBOeRu4e2lrTyk9e2cPOI+1cNSefv7pkmopdIpKZcf3iIianJ/B3v3uXnUfbvI404ancPfLQuio+/rO3GHLw+RXTWDErF98pdtkTiQRxMT4+dVEJSXF+bv1VBfXtvV5HmtBU7uOsd2CQrz26hS8/spmlJVnc/oEZFGvDL4kS6Ymx/PIz5TR29vLpe9ZS197jdaQJS3Pu46i2rYe/+e163j3Uwhcunc6XrjyPB9dVBfz3z2arAhGvzC/M4O5PlfPXv1nPdT97m9/cupSp2clex5pwNHIfJ+/sa+Qj//Umu46287NPLuIrV52P36dpGIlOK2blcv/nl9HW088nfvY22w63eh1pwlG5h9jgkONHL+3mll+8Q2p8DI/dvpyV87TRkkS/hcWZPPw3FxHrN6796Vv88o19uv/qODKv9oQoLy93FRUVnrx2qB2bPmnt7uf3FVXsa+hkYVEGqxZM0Z4wEvVOvM1eXXsPX390Cy/tqGNpaRbfv27+8Wmas5lq1O37hpnZeudc+ZmO05x7CAw5xzv7Gnlhey3OOa5bVMiiqZlexxLxxKTUBH7x6XIeXl/NnU9u5/IfvMaNS4r420tneB0tqqncg2z9wSZ+/tpeqpu7mZWXwqoLCrTfhkx4Zsb15UW8b2YO//XHSh5cV8WD66pYUJTBkpIsCjIST3nDbRkblXuQrN3fxI9f3s3qykZS4mO4cUkR8wvS9Q0rMsrk9ET+37XzuP0DM/jpK5U8VFHFugPN5KXFs6g4kzlT0jUYChKV+1k4cX6wp3+QTdUtrD/YTHVzN8nxMaycm8+y0mziYnSuWuRUCjIS+fa185iem8Lm6lbWH2zi2a1HeXbrUaYGnocuAAAHqklEQVSkJ1A2JY3z8tOYkp6gAdIYqdzPUnffILtq29h+pJ1dR9voH3TkpcXzkfmTKZ+apVIXOQsJsX6WlmaxtDSLps4+th1uZdvhNl7aUcdLO+pITYhhVl4qMyal8KE5eWTrtpIBU7mfweCQY3N1C6srG3hkQw0HGzsZcpAaH8Oi4kwWT83UfKFIEGQlx3HJzFwumZlLe08/e2o72FXbzrbDraw/2MyD66qYMyWNpaVZLCnJorwkk0mpCV7HDlsBlbuZXQX8GPADv3TOffeEz8cDvwYWA43Ajc65A8GNOj56+gfZdriVtfubWbu/kYoDzbT3DgAwJT2BS2bmMntyGoWZidoLRiREUhNiWTQ1k0VTMxkcchxu6SY53s/qykZ+t/YQ964+AEB+WgJzC9Iom5LOjEkpTMtJpiQnWfdAIIByNzM/cBdwBVANrDOzJ5xz20cddivQ7JybYWY3Af8B3BiKwMHU0z/I7tp2dhxpY8eRdt6tamH74Vb6B4fX/k/PTeajC6Zw4bRslk/P5vlttR4nFpl4/D6jKCuJW5YVc8cHZ9I3MHR8NL+1ppWth9v44846Rl8flZ4YS35aAvnpCeSkxJOdEkdWchyZSbGkJ8aRnhhLemIsGUnD/50U54+6374D+fG2FKh0zu0DMLMHgGuA0eV+DfCvIx8/DPzEzMyN0xVSQ0OO/qEhBgYdvQND9A4M0t03SHvPAO09A7T19FPf3kt9ey+1bT1UNXdxsLGLI61/2tQoKc7PvIJ0bn3fNBYWZ7CoOJPcVM3viYSLExc0JMXFsLQ0m6Wl2fQPDtHY0UdDRy+NHb20dPfT1jNAZV0HG6ta6OwdYOA0V8fG+Ox40Z/4Jy0xlrSEWFITYkhJiCE5PoaU+BgSY/0kxPpJiPURF+Mjzu8j1u8jxm/4zfD7zNMfGIGUewEweneramDZqY5xzg2YWSuQDQT9tiy/fGMf33t+F845htzwBUOB/gjxGeSkxFOUlcRF07OZmpXMrLwUZk9OozgrCZ/2ehGJSLF+H/npwyP1k3HO0Tc4RHffIF0jf7r7B5lfmE5rd/+f/nT109LdR31HL5X1HbR29dPeOxBwx5zIbPhWhD4Dwxj5D/+6ag43Lw3tFbfjOjFlZrcBt4087DCzXUAOIfghcCr7gXWBHTquuc6Sso2Nso3NWWX7ZAiDnERE/rvd8m24Zexfd2ogBwVS7jVA0ajHhSPPneyYajOLAdIZPrH6Hs65u4G7Rz9nZhWB7JMw3sI1FyjbWCnb2Cjb2HidLZBF2euAmWZWamZxwE3AEycc8wTwmZGPrwP+OF7z7SIi8ufOOHIfmUO/A3ie4aWQ9zjntpnZnUCFc+4J4H+A35hZJdDE8A8AERHxSEBz7s65Z4BnTnjum6M+7gGuH2OGu898iCfCNRco21gp29go29h4ms2z/dxFRCR0tBGKiEgUCptyN7N/MjNnZjleZznGzL5lZpvNbKOZvWBmU7zOdIyZfd/Mdo7k+4OZZXid6Rgzu97MtpnZkJmFxUoGM7vKzHaZWaWZfdXrPMeY2T1mVmdmW73OMpqZFZnZK2a2feT/yy96nekYM0sws7Vmtmkk2795nelEZuY3s3fN7CmvMoRFuZtZEXAlEPg9t8bH951z851zC4CngG+e6S+MoxeBuc65+cBu4Gse5xltK/Bx4HWvg8B7ttBYCZQBN5tZmbepjvsVcJXXIU5iAPgn51wZcCFwexj9m/UCH3TOXQAsAK4ysws9znSiLwI7vAwQFuUO/BD4MhBWJwCcc22jHiYTRvmccy845wZGHr7D8PUHYcE5t8M5t8vrHKMc30LDOdcHHNtCw3POudcZXmEWVpxzR5xzG0Y+bme4qAq8TTXMDesYeRg78ids3ptmVgh8GPillzk8L3czuwaocc5t8jrLyZjZt82siuEL78Jp5D7aXwLPeh0ijJ1sC42wKKpIYGYlwEJgjbdJ/mRk2mMjUAe86JwLm2zAjxgerA55GWJcth8ws5eA/JN86hvA1xmekvHE6bI55x53zn0D+IaZfQ24A/i/4ZJt5JhvMPwr9H3jlSvQbBL5zCwFeAT4hxN+k/WUc24QWDByrukPZjbXOef5eQsz+whQ55xbb2aXepllXMrdOXf5yZ43s3lAKbBpZPe0QmCDmS11zh31MttJ3MfwWv9xK/czZTOzzwIfAS4b7yuCz+LfLRwEsoWGnMDMYhku9vucc496nedknHMtZvYKw+ctPC93YDmwysyuBhKANDP7rXPuL8Y7iKfTMs65Lc65Sc65EudcCcO/Li8ar2I/EzObOerhNcBOr7KcaOQGKl8GVjnnurzOE+YC2UJDRrHh0db/ADuccz/wOs9oZpZ7bHWYmSUyfK+JsHhvOue+5pwrHOmzmxjeimXcix3CYM49zH3XzLaa2WaGp47CZjkY8BMgFXhxZKnmz70OdIyZXWtm1cBFwNNm9ryXeUZOPB/bQmMH8JBzbpuXmY4xs98BbwPnmVm1md3qdaYRy4FPAR8c+f7aODIaDQeTgVdG3pfrGJ5z92zJYbjSFaoiIlFII3cRkSikchcRiUIqdxGRKKRyFxGJQip3EZEopHIXEYlCKncRkSikchcRiUL/H/4T7TOc8YwlAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "sns.distplot([i['x'] for i in trace])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python (pymc3)",
   "language": "python",
   "name": "pymc3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}