pymc3-svgd-averaging.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pymc3 as pm\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import theano.tensor as tt\n",
    "import theano\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [],
   "source": [
    "with pm.Model() as model:\n",
    "    dummy = pm.Normal('norm', 0, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 10000/10000 [00:22<00:00, 447.28it/s]\n"
     ]
    }
   ],
   "source": [
    "# the below code is model agnostic\n",
    "with model:\n",
    "    approx = pm.fit(method='svgd', inf_kwargs={'n_particles': 100})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<pymc3.variational.approximations.Empirical at 0x117666d30>"
      ]
     },
     "execution_count": 63,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# this is an empirical approximatoin\n",
    "approx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "histogram"
      ]
     },
     "execution_count": 64,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# let's introspect what is inside\n",
    "# every empirical has the corresponding histogram storage\n",
    "# the histogram is of size (n_samples, model_dim)\n",
    "approx.histogram # histogram storage"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "theano.tensor.sharedvar.TensorSharedVariable"
      ]
     },
     "execution_count": 65,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# this is regular theano variable\n",
    "type(approx.histogram)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [],
   "source": [
    "# this is a link to all plotting functions and those who accept trace object as input\n",
    "trace = approx.sample(100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x1175afe10>"
      ]
     },
     "execution_count": 62,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbMAAAEYCAYAAADWNhiqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAGzNJREFUeJzt3Xt8FPW9//H3bHaTAE2UMBhciuAFvFTwAgptELkJVqoghhFrKcLR/rxQpajVIyJoreg5NSoe+aGCKK0FB7R4xAqiWCmg1IJYLygqAgpBTQMkhFw2yZw/iHmwBshOyGb3m7yej0cej+zsZ775zJedeTOzk43leZ4AADBZINENAABwpAgzAIDxCDMAgPEIMwCA8QgzAIDxCDMAgPEIMwCA8QgzIAlYlhVMdA+AySx+aRpoHJZldZH0iqRVkn4iabuk4ZJOljRLUmtJn0sa73neLsuy/iZpg6S+kuZL6i6pVNJZko6RNF7SLyX9WNJaz/OuarKNAQzDmRnQuLpKeszzvB9J2i3pMknzJN3meV4PSe9LmnpAfarneb08z3uw5nFb7Q+v30j6X0kPSfqRpO6WZZ3ZRNsAGIcwAxrXF57nbaj5fp2kEyUd7XnemzXLnpHU74D65763/kve/ssl70v62vO89z3Pq5b0oaQu8WsbMBthBjSu8gO+r5J0dD31JYdYv/p7Y1VL4n014BAIMyC+9kjaZVnWeTWPx0h68zD1ABqA/+kB8TdW0izLslpL2ixpXIL7AZod7mYEABiPy4wAAOMRZgAA4xFmAADjEWYAAOPF5W7GHTt2NJu7SrKyslRYWJjoNpIKcxKN+YjGfNTFnETzMx/hcNiKpY4zs3oEAkzR9zEn0ZiPaMxHXcxJtHjMBzMMADAeYQYAMB5hBgAwHmEGADAeYQYAMB5hBgAwHmEGADAeYQYAMB5hBgAwHmEGADBeXP44Z1lZmddcPr4lGAyqsrIy0W0klWAwqF1llSqujP21kxG0lGHw3zUvrtQht9eyLH1/P0q27T1c/wdzJP2zz9TFnETzMx+pqakxfTZjXHa35vSBmrZtq6CgINFtJBXbtlVYWqHnt5XFvM5lx6WrvLwojl3FV0laptHb25T9s8/UxZxE8zMf4XA4prrmcfoEAGjRCDMAgPEIMwCA8QgzAIDxCDMAgPEIMwCA8QgzAIDxCDMAgPEIMwCA8QgzAIDxCDMAgPEIMwCA8QgzAIDxCDMAgPEIMwCA8QgzAIDxCDMAgPEIMwCA8QgzAIDxCDMAgPEIMwCA8QgzAIDxCDMAgPEIMwCA8QgzAIDxCDMAgPGC8Rg0KytLgUDzyMlgMCjbthPdRlIJBoMKhTxJZTGvEwqFZGfEPo/FlVJxpRdzfUbQUkZcXs37VZQ1bHs9z9OkSZO0dOlStW7dWrNnz9ZZZ51Vp379+vW6+uqrVVpaqgsvvFB5eXmyLEvTpk3TSy+9pEAgoPbt22v27NkKh8N68803lZubqy5dukiSRowYocmTJzda/2mpqaoItYu5/sD5Z5+pizmJFo/5iMvuX1hYGI9hE8K2bRUUFCS6jaRi27YikYivdSKRiAr2FsVcX5KWqee3xX7wvey4dJWXxz6+X5G0TH/1Ndv7+uuva+PGjVq5cqXWr1+v6667TkuWLKlTf91112n69Ok6++yzNWbMGC1cuFADBw7U2LFjNWHCBEnSnDlzNGXKFD3wwAPas2ePzjnnHM2bN692jMO9Tv32vy9SpcVbS2KuP3D+2WfqYk6i+ZmPcDgcU13zOH0Canz55Zfq16+fJk6cqL59+2rChAlauXKlhg8frpycHL377rvat2+fJk2apGHDhmnIkCFatmxZ7bqXXnqphg4dqqFDh+qdd96RJK1Zs0a/uOSnevbWccob+WMtmHytPC+2s8Zly5YpNzdXlmWpZ8+e2rNnj77++uuomq+//lrFxcXq2bOnLMtSbm6uli5dKknKyMiordu3b58sy2qMaQKanThemAESY8uWLXr88ceVl5eniy66SIsXL9bixYv16quv6tFHH1XXrl2Vk5OjvLw87dmzR8OGDdN5550n27Y1f/58paena/Pmzbrhhhv0yiuvSJI+ev9fmuD+XRntO2jWuGHaumGtupzVR0v+cKc2/3NVnR4KLh+lSddeo507d0b9z/LYY4/Vzp07lZ2dXbts586dOvbYY+vUfOf+++/XokWLlJmZqYULF9YuX7dunQYPHqwOHTpoypQpOvnkkxt1HgGTEGZodjp16qRTTz1VktStWzf17dtXlmXplFNO0Zdffqn8/HwtX75cs2bNkiSVl5dr+/btys7O1uTJk/XRRx8pEAho8+bNtWP2OLunjsreH0rhk0/Xrh1fqstZffSzW+49aA+XHZcuNdJlz9tvv1233367Hn30Uc2dO1e33HKLunfvrn/84x9q06aNXn/9dY0fP16rV69ulJ8HmIgwQ7OTlpZW+30gEFBqamrt91VVVUpJSdETTzyhk046KWq9Bx98UO3bt9fy5ctVXV2tE044ofa578aQJCsQUHVVpSTVe2bWoUMH7dixo3Z5fn6+OnToEFXboUMH5efnH7ZGkkaOHKkxY8bolltuibr8OGjQIN1xxx0qLCxUVlbW4ScHaKYIM7Q4559/vubOnat7771XlmXpgw8+0Omnn66ioiIde+yxCgQCWrhwoaqqquodq74zsyFDhujpp5/W8OHDtX79emVmZkZdYpSk7OxsZWRkaN26dTr77LO1aNEijRs3TpK0efPm2lBdtmyZTjzxREnSN998o/bt28uyLL377ruqrq5W27Ztj2RaAKMRZmhxJk6cqKlTp2rw4MGqrq5Wp06dNG/ePI0dO1a/+tWvtGjRIg0YMECtW7c+4p81aNAgrVixQjk5OWrVqpXy8vJqn7vgggu0fPlySdJ9992n3/zmNyorK9OAAQM0cOBASdL06dP1+eefKxAIqGPHjrr//vslSS+//LLmzZunlJQUpaena+bMmdwcghbNivWuLD927NjR+IMmCLfU1mXbtrYWV/i+db6Nj/eQGnJrvp/x/Uq2fvzy2/+Izm1835rfhlvzD4k5iebz1vyY/pfGrfkAAOMRZgAA4xFmAADjEWYAAOMRZgAA4xFmAADjEWYAAOMRZgAA4xFmAADjEWYAAOMRZgAA4xFmAADjEWYAAOMRZgAA4xFmAADjEWYAAOMRZgAA4xFmAADjEWYAAOMF4zFoVlaWAoHmkZPBYFC2bSe6jaQSDAYVCnmSymJeJxQKyc6IfR4ryuI7vl/J1k9xpVRc6cVcH/AkP/1bPvffA7eXfaYu5iRaPOYjLmFWWFgYj2ETwrZtFRQUJLqNpGLbtiKRiK91IpGICvYWxV6flhnX8f1Ktn5K0jL1/LbYw2lE5za+xveqq33VH7i97DN1MSfR/MxHOByOqa55nD4BAFo0wgwAYDzCDABgPMIMAGA8wgwAYDzCDABgPMIMAGA8wgwAYDzCDABgPMIMAGA8wgwAYDzCDABgPMIMAGA8wgwAYDzCDABgPMIMAGA8wgwAYDzCDABgPMIMAGA8wgwAYDzCDABgPMIMAGA8wgwAYDzCDABgPMIMAGA8wgwAYLxgPAbNyspSINA8cjIYDMq27US3kVSCwaBCIU9SWczrpKWmqiLULub6gKe4jp+eYqmsykuafjKCljJ87I0VZf7m3/K5P/qtP3B788sl7weH33a/22s6jiPR4jEfcXk5FRYWxmPYhLBtWwUFBYluI6nYtq1IJOJrnX2RKi3eWhJz/YjObeI+fjL1c9lx6SovL4q5PpKW6asfr7o6rvXx3l7TcRyJ5mc+wuFwTHXN4/QJANCiEWYAAOMRZgAA4xFmAADjEWYAAOMRZgAA4xFmAADjEWYAAOMRZgAA4xFmAADjEWYAAOMRZgAA4xFmAADjEWYAAOMRZgAA4xFmAADjEWYAAOMRZgAA4xFmAADjEWYAAOMRZgAA4xFmAADjEWYAAOMRZgAA4xFmAADjEWYAAOMF4zFoVlaWAoHmkZPBYFC2bddbV1wpFVd6MY+bEbSUEZfZj79gMKhQyJNUFvM6ls/XQ0urD4VCsjPqf519p6LM7Pn3u71+9y9JSk+xVFaVHPtkrMeRliIe8xGXf7rCwsJ4DJsQtm2roKCg3rqStEw9vy32g8tlx6WrvLzoSFpLGNu2FYlEfK3jVVdTfxiRSEQFe2N/PUTSMuPaT7Jtr9/9S5JGdG6jxVtLYq6P5z4Z63GkpfAzH+FwOKa65nH6BABo0QgzAIDxCDMAgPEIMwCA8QgzAIDxCDMAgPEIMwCA8QgzAIDxCDMAgPEIMwCA8QgzAIDxCDMAgPEIMwCA8QgzAIDxCDMAgPEIMwCA8QgzAIDxCDMAgPEIMwCA8QgzAIDxCDMAgPEIMwCA8QgzAIDxCDMAgPEIMwCA8QgzAIDxgvEYNCsrS4FA88jJYDAo27brraso8ySVxTxuKBSSnVH/uMkoGAwqFPK3vZbP10NLq09LTVVFqF3M9QFPMnn+4729DekpnvtkrMeRliIe8xGXMCssLIzHsAlh27YKCgrqrYukZfoaNxKJqGBvUUPbSijbthWJRHyt41VXU38Y+yJVWry1JOb6EZ3bxLUf07dX8t9TPPfJWI8jLYWf+QiHwzHVNY/TJwBAi0aYAQCMR5gBAIxHmAEAjEeYAQCMR5gBAIxHmAEAjEeYAQCMR5gBAIxHmAEAjEeYAQCMR5gBAIxHmAEAjEeYAQCMR5gBAIxHmAEAjEeYAQCMR5gBAIxHmAEAjEeYAQCMR5gBAIxHmAEAjEeYAQCMR5gBAIxHmAEAjEeYAQCMF4zHoFlZWQoE4peTxZVScaUXc31G0FJGA7c0GAzKtu166yrKPEllMY8bCoVkZ9Q/7neacpvrEwwGFQr5217L5+uBeuqPpL4h6/jdJ/2I9TiSzBrzGBSP+YjL4a6wsDAew9YqScvU89tiP5Bedly6ysuLGvSzbNtWQUFBvXWRtExf40YiERXsjb2nptzm+ti2rUgk4msdr7qaeuqbrL4h6/jdJ/2I9TiSzBrzGORnPsLhcEx1XGZEUlj958f18Kjz9FBuX616dlbt8tdm/ZemD+2uGaP7a8bo/vp41XJJ0pYNa/WIc77+58rBKtj2uSSptHiP5lw/StWHOIg9cc1wffXRhtrHu3Zs08OjzpMkbf7nak3rd4JmjO6vvJE/0WuP/3fU8pH9++jBS/vo8f+4WBtXvhqXOQDQcHG6EAXEbudnG/XOX/6k6+ctU0ooVXMnXK6tm0dIKR0kSTlXXqt+v7whap1Vf/z/uurR+dq1Y5vWLnpGV5/3oN6YnacB4yc2+BJ3lzP76KoZf1ZFaYlmjB6gU/sNqV3+wosvavHWEu345H39cdJYhdLSdVLvfke24QAaDWdmSLhvv9ikTqefrdRWrZUSDOr4nj/R8iUvHnadQDCoirJSVZSVKiUY1LYvNmv3zu06oVfOEfeT2qqNOp56hv795Rd1nguf3F2DrrlZb7lzjvjnAGg8nJkh4bJPPFXLHrtPJbsLFUpL1yerXlN271767u3ht56bo3eXuOp42hkaNuketco8Wv3H36SFU25QKD1dzu9m6uHfT9OQG/6z3p/13ORrFUpL1zOpARWWlB/0JoGS3YXa9v46DbzmZpXs+ned58On9tDKeY8d6WYDaESEGRLumBO66fyrfq2nrh+l1FatFT75dAVSUiRJvUddpYHX3CxZlpbPnK6X8+5S7rQZCp/cXdfPWypJ+mLdGrXP7iDPk/5829VKCYZ00aS7ldHumDo/6/Lfz9IPTztTIzq30dy3NuqZm66sfW7Lhrc144oBsqyA+o+7UdknnqLN/1xdt2Ev9ju6ADQNwgxJ4ZwRv9A5I34hSVr26L3qcmIXSYoKpHNHjokKH0nyPE8r5jykP/7pT7r21zfppzdN1a78bVoz/0kNnTDZVw/fvWdWnx0fv69jju/ma2wA8cV7ZkgKewu/lSTtzv9KH77xsn6We7kkqejbnbU1H674q7JPPCVqvfVLntPJOYN0dNssRcpKZQUCsqyAImWlcekzf9OHWjE7T32c8XEZH0DDcGaGpPDsLeO0b88uBYIhXXLbA8o86mhpd4leeeQe5W/6QJYstQ130ojJf6hdp6J0n9a/tEDjH1soSep75bV6+sYrlBIMafR9jzdab1s2vK2R/fvo26IS/aCtrYtvvY87GYEkQ5ghKfy/p5YcdPnl98485DqprVrrmicW1z4+/uwfa6K78pD1v3oy+g7JtuHjNHHh3yVJJ/TKOeidkCf0ytG0lZs1onMbLd5acthtAJA4XGYEABiPMAMAGI8wAwAYjzADABiPMAMAGI8wAwAYjzADABiPMAMAGI8wAwAYjzADABiPMAMAGI8wAwAYjzADABiPMAMAGI8wAwAYjzADABiPMAMAGC8uf2k6KytLgUD8crKizJNUFnN9KBSSnWHHXF9cKRVXepKk/HLJ+0G7etcJePLVU1pqqipC9Y/bVOOnp1gqq/Jiqs0vlwLBkK9+LJ+vB+qpP5L6hqxzJMeJ+uSXS2lH2THvY5KUEbSU4eMI7aefhozfmMfdYDAo2459rmMRlzArLCyMx7C1ImmZ/uojERXsLYq5viQtU89vi/0fTZJGdG7jq35fpEqLt5Yk1fjx7Merrqae+iarb8g68T5O+N3HLjsuXeXl8evH7/iNedy1bVsFBQUxjRMOh2Oq4zIjAMB4hBkAwHiEGQDAeL7CzLKsKy3L+pdlWe9blrXGsqwzDlc/ZcoUde3atfbx9u3blZubqyFDhmjw4MF6/fXXJe1/jy03N1ddu3bV5MmTG7IdABBXkyZNUo8ePTRw4MCDPr/ymf/RjNH9NWN0fz086jzd0Stb+/bs0rdbPtOM0f116fm9NWN0f00773itenaWJOm1Wf+l6UO716738arlTblJzYrfG0C+kHS+53m7LMv6qaQnJPU+WOF7772n3bt3Ry175JFHdPHFF2vs2LHatGmTxowZo7Vr1yo9PV2//e1v9fHHH+uTTz5p0IYAQDw5jqNx48bppptuOujz/cZOUL+xEyRJG99cplXPzlLro9qq9VFtdeOCv2lE5zZ6YXORpl/YXT8aMKx2vZwrr1W/X97QJNvQnPk6M/M8b43nebtqHr4t6YcHq6uqqtLvfvc73XnnnXWe27t3rySpqKhI2dnZkqTWrVvr3HPPVVpamp92AKDJ9OnTR0cffXRMte8te0FnXDiyzvLP/rFS7X7YRW3DnRq7vRbvSG7N/w9Jrxzsiblz52rIkCG1YfWdm2++WT//+c/11FNPqbS0VAsWLDiCHw8AyaeidJ82rVmhS267v85z/1r2F/UYGh1ybz03R+8ucdXxtDM0bNI9apUZW2AiWoNuALEsa4D2h9ltB3kuvGTJEo0fP77OeosXL9aoUaO0bt06zZs3TzfeeKOqG/D7IwCQrD5euUydzzhXrY9qG7W8oqJCG1cuU/cLLqld1nvUVbr1f9/Rrxe8oQw7Wy/n3dXU7TYb9YaZZVk3WJa1oeYrbFlWD0mzJQ33PO/fB1nlrC1btignJ0e9e/dWaWmpcnJyJEkLFizQxRdfLEnq1auXysvL4/4L1gDQlN57dfFBLzH+/bVlCp/SQxntjqldltHuGAVSUhQIBHTuyDH66sN3m7LVZqXeMPM87zHP8870PO9M7b8s+YKkMZ7nbTpE/csbNmzQ2rVrtXbtWrVq1UqrV6+WJHXs2FGrVq2SJH366acqLy9Xu3axf+QSACSzsuIifbFujU7rf2Gd515+YaHOGHpp1LKib3fWfv/hir8q+8RT4t5jc+X3PbO7JLWTNNOyLEmq9DyvlyRZlvVXSVd7nrfjkCvfdZduvfVWPfnkk7IsSw899JBqxlHv3r21d+9eVVRUaOnSpZo/f766devWoI0CgMZ2/fXX66233lJhYaF69uypX99+p975plSS1Dv3KknSh2+8rK59+iu1VfTHvVWUlmjN31Zo4sQHopa/8sg9yt/0gSxZahvupBGT/9Ak29Ic+Qozz/OulnT1IZ676GDLP/3009rvu3XrphdffPGgY69du9ZPKwDQpGbOnBn1uCQtU4HvfRZiz0uuUM9LrqizbmqrNnr7s6/qfDbj5ffOrFOLhuETQAAAxiPMAADGI8wAAMYjzAAAxiPMAADGI8wAAMYjzAAAxiPMAADGszzPa/RBd+zYcUSDduzYsbFaAYBGs3379trvS9Iy9fz3fmn6cEZ0blPnl6YP57Lj0tWmvCjmer/9JHJ827ZVUFAQ0zjhcNiKpY4zMwCA8QgzAIDxCDMAgPGO5C9Nx82B16UPJtmuDUv+r4dTTz318atvyDp+jxNILpyZAQCMF5e7Ge++++7GHxQA0CJNnTq1/jsaPc/j6zBf06ZN8xLdQ7J9MSfMB/PBnCTbfHCZEQBgPMKsfncnuoEkxJxEYz6iMR91MSfRGn0+4vKeGQAATYkzMwCA8QgzAIDxCDMAgPGS8hNAko3jOP8t6WJJFZI+lzTOdd3die0qcRzHGSVpmqRTJZ3ruu4/E9tRYjiOc6GkRySlSJrtuu79CW4poRzHeUrSzyR947ru6YnuJ9Ecx+kkaZ6kbEmepCdc130ksV0ljuM46ZJWSkrT/uxZ5Lru1MYanzOz2CyXdLrruj0kbZL0nwnuJ9E+kDRS+1+YLZLjOCmSHpP0U0mnSbrCcZzTEttVwj0t6cJEN5FEKiXd7LruaZL6SLqhhb9GyiUNdF33DElnSrrQcZw+jTU4Z2YxcF331QMevi0pN1G9JAPXdTdKkuM4iW4lkc6V9JnrupslyXGcBZKGS/oooV0lkOu6Kx3H6ZLoPpKF67r5kvJrvi92HGejpI5qoa8R13U9SXtrHoZqvhrtdnrCzL/xkp5LdBNIuI6Svjzg8VeSeieoFyS5mpA/S9LaBLeSUDVXNNZJOknSY67rNtp8EGY1HMd5TVKHgzw12XXdF2tqJmv/pYNnm7K3RIhlPgDUz3GcH0h6XtJE13Vb9Mfyu65bJelMx3GOlvQXx3FOd133g8YYmzCr4bru4MM97zjOVdr/5vagmtPlZq2++YC2S+p0wOMf1iwDajmOE9L+IHvWdd0XEt1PsnBdd7fjOG9o/3ushFlTqblr7beSznddd1+i+0FSeEdSV8dxjtf+EBst6eeJbQnJxHEcS9IcSRtd181LdD+J5jhOe0mRmiBrJekCSQ801vh8nFUMHMf5TPtvJ/13zaK3Xde9NoEtJZTjOJdKelRSe0m7JW1wXXdoYrtqeo7jXCTpYe2/Nf8p13V/n+CWEspxnPmS+kuyJX0taarrunMS2lQCOY7TV9LfJb0vqbpm8R2u6/41cV0ljuM4PSQ9o/37S0CS67ruPY01PmEGADAev2cGADAeYQYAMB5hBgAwHmEGADAeYQYAMB5hBgAwHmEGADDe/wHo+8zHfiOV+AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pm.plot_posterior(trace)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Trace object can be converted to Empirical approx that can manipulate the graph. This is useful to wrap over NUTS output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [],
   "source": [
    "approx2 = pm.Empirical(trace, model=model)\n",
    "# in our case it is just an example"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So the question was about model averaging"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [],
   "source": [
    "# suppose you have a graph\n",
    "out = tt.sin(dummy * 100.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Basicaly there are two ways to do do model averaging\n",
    "\n",
    "* $\\bar{f} = \\mathbb{E}{f(\\xi)}$ -- model averaging\n",
    "* $\\bar{f} = f(\\mathbb{E}{\\xi})$ -- mean propagation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# model averaging"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [],
   "source": [
    "mout = approx.sample_node(out, 100).mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [],
   "source": [
    "samples1 = [mout.eval() for _ in range(1000)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAADbZJREFUeJzt3W+MXOdVx/Hv4m0jBJSQDhhvbOFKNahuVYwwSRGlCkoKSRXZQSqnLlLrlMjLi0QI0TemfdEKhOSKQhsBCthtVRtB40OgzaqxWqgRlLxIS1v1T0qAmOISezc221ohKFKrmOHFXFsTZ3ZnvDt/1me+H2m19z73mbknJzO/vXv3zvVMu91GklTX90y6AEnSaBn0klScQS9JxRn0klScQS9JxRn0klScQS9JxRn0klScQS9Jxc1OuoCGH8+VpLWZ6TdhowQ9i4uLPcdbrRbLy8tjrmZjsycvZD9ezJ68WMWezM3NDTTPUzeSVJxBL0nFGfSSVJxBL0nFGfSSVJxBL0nFGfSSVJxBL0nFGfSSVNyG+WSstFFdPLBnIvvddGRhIvtVPR7RS1JxBr0kFWfQS1JxBr0kFWfQS1JxBr0kFWfQS1JxBr0kFWfQS1JxBr0kFectEHRNWO02BOfGWId0LfKIXpKKM+glqTiDXpKKM+glqTiDXpKKM+glqTiDXpKKM+glqTiDXpKKM+glqbi+t0CIiG3AMWAz0AYOZ+b9EXEDcBzYDpwGIjMvRMQMcD/wJuA54O7M/NJoypck9TPIEf3zwDszcyfwOuDeiNgJHAROZuYO4GSzDnAHsKP5mgceGHrVkqSB9Q36zFy6dESemc8CTwA3AnuBo820o8BdzfJe4FhmtjPzMeD6iNgy9MolSQO5qrtXRsR24KeAzwGbM3Op2fQ0nVM70Pkh8FTXw840Y0tIGthqd+zsZz139Nx0ZGEdj9ZGNHDQR8T3A38N/GZm/k9EXN6Wme2IaF/NjiNins6pHTKTVqvVu8DZ2RW3Tatp7Im3Ih6fqq+taXzfXDJQ0EfES+iE/F9k5t80w+ciYktmLjWnZs4342eBbV0P39qMvUBmHgYON6vt5eXlnvtutVqstG1a2RONUtXXVsX3zdzc3EDzBrnqZgb4MPBEZv5h16YFYD9wqPn+cNf4fRHxIHAz8EzXKR5J0pgNckT/c8DbgK9FxJebsXfRCfiMiHuAbwKXzuWcoHNp5Sk6l1e+Y6gVS5KuSt+gz8xHgZkVNt/aY34buHeddUmShsRPxkpScQa9JBVn0EtScQa9JBVn0EtScQa9JBVn0EtScQa9JBVn0EtScQa9JBVn0EtScQa9JBVn0EtScQa9JBVn0EtScQa9JBVn0EtScQa9JBVn0EtScQa9JBVn0EtScQa9JBVn0EtScQa9JBVn0EtScQa9JBVn0EtScQa9JBVn0EtScQa9JBVn0EtScQa9JBVn0EtScQa9JBVn0EtScQa9JBU3229CRHwEuBM4n5mvacbeCxwA/ruZ9q7MPNFs+23gHuAi8BuZ+ekR1C1JGlDfoAc+CvwxcOyK8Q9k5vu7ByJiJ7APeDUwB3wmIn48My8OoVZJ0hr0PXWTmZ8Fvj3g8+0FHszM72TmfwKngJvWUZ8kaZ0GOaJfyX0R8XbgC8A7M/MCcCPwWNecM82YJGlC1hr0DwC/C7Sb738A/NrVPEFEzAPzAJlJq9XqXeDs7IrbptU09uTcpAuYIlVfW9P4vrlkTUGfmZffdxFxBPhks3oW2NY1dWsz1us5DgOHm9X28vJyz321Wi1W2jat7IlGqeprq+L7Zm5ubqB5a7q8MiK2dK3+MvB4s7wA7IuI6yLiFcAO4PNr2YckaTgGubzyY8AtQCsizgDvAW6JiF10Tt2cBn4dIDO/HhEJ/AvwPHCvV9xI0mTNtNvtSdcA0F5cXOy5oeKvW+s1jT25eGDPpEuYGpuOLEy6hJGo+L5pTt3M9JvnJ2MlqTiDXpKKM+glqTiDXpKKM+glqTiDXpKKM+glqTiDXpKKM+glqTiDXpKKM+glqTiDXpKKM+glqbj1/FOCmkLeRVK69nhEL0nFGfSSVJxBL0nFGfSSVJxBL0nFGfSSVJxBL0nFGfSSVJxBL0nFGfSSVJxBL0nFGfSSVJxBL0nFGfSSVJxBL0nFGfSSVJxBL0nFGfSSVJxBL0nFGfSSVJxBL0nFGfSSVJxBL0nFzfabEBEfAe4Ezmfma5qxG4DjwHbgNBCZeSEiZoD7gTcBzwF3Z+aXRlO6pFG4eGDPxPa96cjCxPZd2SBH9B8Fbr9i7CBwMjN3ACebdYA7gB3N1zzwwHDKlCStVd+gz8zPAt++YngvcLRZPgrc1TV+LDPbmfkYcH1EbBlWsZKkq9f31M0KNmfmUrP8NLC5Wb4ReKpr3plmbIkrRMQ8naN+MpNWq9W7wNnZFbdNq0n25NxE9qppMcrX9TRnyVqD/rLMbEdEew2POwwcblbby8vLPee1Wi1W2jat7ImqGuXruuL7Zm5ubqB5a73q5tylUzLN9/PN+FlgW9e8rc2YJGlC1npEvwDsBw413x/uGr8vIh4Ebgae6TrFI0magEEur/wYcAvQiogzwHvoBHxGxD3AN4Fopp+gc2nlKTqXV75jBDVLkq5C36DPzLeusOnWHnPbwL3rLUqSNDx+MlaSijPoJak4g16SijPoJak4g16SijPoJak4g16SijPoJak4g16SijPoJak4g16SijPoJak4g16SijPoJak4g16SijPoJak4g16SijPoJak4g16SijPoJak4g16SijPoJak4g16SijPoJak4g16SijPoJak4g16SijPoJak4g16SijPoJak4g16SijPoJak4g16SijPoJak4g16SijPoJam42fU8OCJOA88CF4HnM3N3RNwAHAe2A6eByMwL6ytTkrRWwzii/4XM3JWZu5v1g8DJzNwBnGzWJUkTMopTN3uBo83yUeCuEexDkjSg9QZ9G/jbiPhiRMw3Y5szc6lZfhrYvM59SJLWYV3n6IHXZ+bZiPgR4O8i4l+7N2ZmOyLavR7Y/GCYb+bRarV6Fzg7u+K2aTXJnpybyF41LUb5up7mLJlpt3vm8FWLiPcC/wscAG7JzKWI2AL8Q2b+RJ+HtxcXF3tuaLVaLC8vD6XGKi4e2DPpEqSR2HRkYWTPXTFL5ubmAGb6zVvzqZuI+L6I+IFLy8AvAo8DC8D+Ztp+4OG17kOStH7rOUe/GXg0Ir4CfB54JDM/BRwC3hgRTwK3NeuSpAlZ8zn6zPwG8JM9xr8F3LqeoiRJw7PeP8ZK0tCM8u9Pq11IMMq/DWwE3gJBkooz6CWpOINekooz6CWpOINekooz6CWpOINekooz6CWpOINekooz6CWpOINekooz6CWpOINekooz6CWpOINekooz6CWpOINekooz6CWpOINekooz6CWpOINekooz6CWpOINekooz6CWpOINekooz6CWpOINekooz6CWpuNlJFyBJk3bxwJ6J7XvTkYWR78OgX4dJvjgkaVCeupGk4gx6SSrOoJek4gx6SSrOoJek4kZ21U1E3A7cD2wCPpSZh0axH698kaTVjeSIPiI2AX8C3AHsBN4aETtHsS9J0upGdermJuBUZn4jM78LPAjsHdG+JEmrGFXQ3wg81bV+phmTJI3ZxD4ZGxHzwDxAZjI3N7fi3NW28cgXhl2aJJUyqiP6s8C2rvWtzdhlmXk4M3dn5m5gZqWviPjiatun8cue2A97Yk+6vvoa1RH9PwM7IuIVdAJ+H/CrI9qXJGkVIzmiz8zngfuATwNPdIby66PYlyRpdSM7R5+ZJ4ATQ3iqw0N4jmrsyQvZjxezJy82tT2Zabfbk65BkjRC3gJBkorbcP/wSETcABwHtgOngcjMC1fM+THg43R+UL0E+KPM/NPxVjo+A/ZkF/AA8DLgIvB7mXl8vJWOxyD9aOZ9Cngd8Ghm3jnOGsel361GIuI64Bjw08C3gLdk5ulx1zkuA/TjDcAHgdcC+zLzofFXOX4b8Yj+IHAyM3cAJ5v1Ky0BP5uZu4CbgYMRscrF9te8QXryHPD2zHw1cDvwwYi4fow1jtMg/QD4feBtY6tqzAa81cg9wIXMfCXwAeB9461yfAbsx38BdwN/Od7qJmsjBv1e4GizfBS468oJmfndzPxOs3odG/O/Y5gG6cm/Z+aTzfIicB744bFVOF59+wGQmSeBZ8dV1AQMcquR7l49BNwaEQNde30N6tuPzDydmV8F/m8SBU7KRgzIzZm51Cw/DWzuNSkitkXEV+ncauF9TbhVNVBPLomIm4CXAv8x6sIm5Kr6Udggtxq5PKe57PkZ4OVjqW78vPXKCiZyjj4iPgP8aI9N7+5eycx2RPS8LCgznwJe25yy+UREPJSZ54Zf7XgMoyfN82wB/hzYn5nX7FHLsPohaUJBn5m3rbQtIs5FxJbMXGpC63yf51qMiMeBn6fzq+k1aRg9iYiXAY8A787Mx0ZU6lgM8zVSWN9bjXTNORMRs8AP0vmjbEWD9GMqbcRTNwvA/mZ5P/DwlRMiYmtEfG+z/EPA64F/G1uF4zdIT15K50qkY1NwJUHffkyJy7caaf7/76PTm27dvXoz8PeZWfU3oEH6MZU2YtAfAt4YEU8CtzXrRMTuiPhQM+dVwOci4ivAPwLvz8yvTaTa8RikJwG8Abg7Ir7cfO2aTLkjN0g/iIh/Av6Kzh8gz0TEL02k2hFZ6VYjEfE7EXHpn177MPDyiDgF/BYrX6F0zRukHxHxMxFxBvgV4M8iYipuzeInYyWpuI14RC9JGiKDXpKKM+glqTiDXpKKM+glqTiDXpKKM+glqTiDXpKK+3/Gl7TYxOP4tAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(samples1);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above still has some distribution as we take few MC samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {},
   "outputs": [],
   "source": [
    "mout1 = approx.sample_node(out, 10000).mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {},
   "outputs": [],
   "source": [
    "samples2 = [mout1.eval() for _ in range(1000)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAADlRJREFUeJzt3W+sZPVdx/H3utPWpJAgGVn3siRbzfYBkoDJBjT6YBVKoZJdNPEbsKnbluz6AGIqPJAWExppdB+0GpIa9G6tXWIFvkoN1xaLsAZrY7D/UrUtJsW6hP3DrreQFoN/wu34YM7CAJe9c+fMuXPvd9+v5GbmnPObc777u2c+89vfnJm7aTAYIEmq64dmXYAkqVsGvSQVZ9BLUnEGvSQVZ9BLUnEGvSQVZ9BLUnEGvSQVZ9BLUnG9WRfQ8OO5kjSZTSs1WC9Bz/Hjx1s9vt/vs7i4OKVqNjb74hX2xZD98IpKfTE3NzdWO6duJKk4g16SijPoJak4g16SijPoJak4g16SijPoJak4g16SijPoJam4dfPJWOlMlvbtnuhxJ6dw7M0HF6awF2l2DHppBZO+yLTlC4ymxakbSSrOoJek4gx6SSrOoJek4gx6SSrOoJek4gx6SSrOoJek4gx6SSrOoJek4gx6SSrOoJek4gx6SSrOoJek4gx6SSrOoJek4gx6SSrOoJek4gx6SSrOoJek4gx6SSrOoJek4gx6SSqut1KDiLgIuBfYAgyA+cy8OyLOBx4AtgNHgMjM5yNiE3A38C7gReC9mfm1bsqXJK1knBH9S8BtmXkx8NPAzRFxMXA7cDgzdwCHm2WAa4Edzc9+4J6pVy1JGtuKQZ+ZJ06PyDPzBeBJ4EJgD3CoaXYIuL65vwe4NzMHmfkEcF5EbJ165ZKksaw4dTMqIrYDPwX8E7AlM080m55lOLUDwxeBZ0YedrRZd2JkHRGxn+GIn8yk3++vtvZX6fV6rfdRRcW+ODnrAmZgmr/DiufEpM7Gvhg76CPiHOBB4AOZ+f2IeHlbZg4iYrCaA2fmPDDfLA4WFxdX8/DX6ff7tN1HFfZFDdP8HXpOvKJSX8zNzY3VbqyrbiLiTQxD/tOZ+Zlm9cnTUzLN7alm/THgopGHb2vWSZJmYJyrbjYBfwI8mZm/P7JpAdgLHGhuHxpZf0tE3A9cAXxvZIpHkrTGxpm6+VngPcC/RsTXm3UfYhjwGRE3AU8Dp+dyHmZ4aeVTDC+vfN9UK5YkrcqKQZ+ZXwQ2vcHmK5dpPwBublmXJGlK/GSsJBVn0EtScQa9JBVn0EtScQa9JBVn0EtScQa9JBVn0EtScQa9JBVn0EtScQa9JBVn0EtScQa9JBVn0EtScQa9JBVn0EtScQa9JBVn0EtScQa9JBVn0EtScQa9JBXXm3UBkpa3tG/31PZ1chVtNx9cmNpxtT44opek4gx6SSrOoJek4gx6SSrOoJek4gx6SSrOoJek4gx6SSrOoJek4gx6SSrOoJek4gx6SSrOoJek4lb89sqI+CRwHXAqMy9p1n0Y2Af8Z9PsQ5n5cLPtg8BNwBLwG5n5SAd1S5LGNM7XFH8K+Dhw72vW/0FmfnR0RURcDNwA/CQwBzwWEW/PzKUp1CpJmsCKUzeZ+QXguTH3twe4PzP/NzP/A3gKuLxFfZKkltr84ZFbIuLXgK8At2Xm88CFwBMjbY426yRJMzJp0N8D3AUMmtuPAe9fzQ4iYj+wHyAz6ff7E5Yy1Ov1Wu+jiop9sZq/kKR2qp07r1Xx+bGSiYI+M19+3kXEQeCzzeIx4KKRptuadcvtYx6YbxYHi4uLk5Tysn6/T9t9VGFfqI3q506l58fc3NxY7Sa6vDIito4s/hLwjeb+AnBDRLwlIt4G7AC+NMkxJEnTMc7llfcBu4B+RBwF7gR2RcRlDKdujgC/DpCZ34yIBL4FvATc7BU3kjRbmwaDwaxrABgcP3681Q4q/XesrYp9sbRv96xLOGtsPrgw6xI6Ven50UzdbFqpnZ+MlaTi2lxeqbOQI2tp43FEL0nFGfSSVJxBL0nFGfSSVJxBL0nFGfSSVJxBL0nFGfSSVJxBL0nFGfSSVJxBL0nFGfSSVJxBL0nFGfSSVJxBL0nFGfSSVJxBL0nFGfSSVJxBL0nFGfSSVJxBL0nFGfSSVJxBL0nFGfSSVJxBL0nFGfSSVJxBL0nFGfSSVJxBL0nFGfSSVJxBL0nFGfSSVJxBL0nFGfSSVFxvpQYR8UngOuBUZl7SrDsfeADYDhwBIjOfj4hNwN3Au4AXgfdm5te6KV2SNI5xRvSfAq55zbrbgcOZuQM43CwDXAvsaH72A/dMp0xJ0qRWDPrM/ALw3GtW7wEONfcPAdePrL83MweZ+QRwXkRsnVaxkqTVm3SOfktmnmjuPwtsae5fCDwz0u5os06SNCMrztGvJDMHETFY7eMiYj/D6R0yk36/36qOXq/Xeh9VdNkXJzvZq9aT6s+jszErJg36kxGxNTNPNFMzp5r1x4CLRtpta9a9TmbOA/PN4mBxcXHCUob6/T5t91GFfaE2qp87lZ4fc3NzY7WbNOgXgL3Ageb2oZH1t0TE/cAVwPdGpngkSTMwzuWV9wG7gH5EHAXuZBjwGRE3AU8D0TR/mOGllU8xvLzyfR3ULElahRWDPjNvfINNVy7TdgDc3LYoSdL0+MlYSSrOoJek4lpfXimplqV9u2d27M0HF2Z27Moc0UtScQa9JBVn0EtScQa9JBVn0EtScQa9JBVn0EtScQa9JBVn0EtScQa9JBVn0EtScQa9JBVn0EtScQa9JBVn0EtScQa9JBVn0EtScQa9JBVn0EtScQa9JBVn0EtScQa9JBVn0EtScQa9JBVn0EtScQa9JBVn0EtScQa9JBXXm3UBWr2lfbvPuP3kGtUhaWNwRC9JxRn0klScQS9JxRn0klRcqzdjI+II8AKwBLyUmTsj4nzgAWA7cASIzHy+XZmSpElNY0T/85l5WWbubJZvBw5n5g7gcLMsSZqRLqZu9gCHmvuHgOs7OIYkaUxtr6MfAH8bEQPgjzNzHtiSmSea7c8CW5Z7YETsB/YDZCb9fr9VIb1er/U+Ngqvk1dVa/EcPpuy4rS2Qf9zmXksIi4AHo2IfxvdmJmD5kXgdZoXhflmcbC4uNiqkH6/T9t9SJqttXgOV8qKubm5sdq1mrrJzGPN7Sngr4DLgZMRsRWguT3V5hiSpHYmDvqIeGtEnHv6PnA18A1gAdjbNNsLPNS2SEnS5NqM6LcAX4yIfwa+BHwuMz8PHADeERHfBq5qliVJMzLxHH1mfge4dJn13wWubFOUJGl6/GSsJBVn0EtScQa9JBVn0EtScQa9JBVn0EtScQa9JBVn0EtScW2/1EySpmZp3+7Oj7Hct79uPrjQ+XFnyRG9JBVn0EtScQa9JBVn0EtScQa9JBVn0EtScQa9JBVn0EtScQa9JBXnJ2NbWItP8UlSW47oJak4g16SijPoJak4g16SijPoJak4g16SijPoJak4g16SijPoJak4g16SijPoJak4g16SijPoJam4Df/tlae/QfLkjOuQpPVqwwe9JLU1y68c33xwofNjOHUjScV1NqKPiGuAu4HNwCcy80BXx5IkvbFORvQRsRn4Q+Ba4GLgxoi4uItjSZLOrKupm8uBpzLzO5n5f8D9wJ6OjiVJOoOugv5C4JmR5aPNOknSGpvZVTcRsR/YD5CZzM3NTbajz31lilVJUj1djeiPAReNLG9r1r0sM+czc2dm7gQ2tf2JiK9OYz8VfuwL+8J+OKv6YkVdjei/DOyIiLcxDPgbgF/t6FiSpDPoZESfmS8BtwCPAE8OV+U3uziWJOnMOpujz8yHgYe72v8y5tfwWOudffEK+2LIfnjFWdcXmwaDwaxrkCR1yK9AkKTiNtSXmkXE+cADwHbgCBCZ+fwy7fYCv90sfiQzD0XEucA/jDTbBvxZZn6g06I70qYvmvVvBj4O7AJ+ANyRmQ92XngHptAXjwNbgf9utl2dmae6rXr62vbDyPYF4Mcz85JOC+7QFM6JzzM8J3oMc+PmzFzqvvJubLQR/e3A4czcARxull+l+QXfCVzB8BO6d0bEj2TmC5l52ekf4GngM2tY+7RN3BfN5juAU5n5doZfU/H3a1J1N9r2BcC7R86PDRfyjdb9EBG/DPzX2pTbqbZ9EZl5KXAJ8KPAr6xJ1R3ZaEG/Bzg9+jgEXL9Mm3cCj2bmc80r+KPANaMNIuLtwAW8eoS/0bTti/cDvweQmT/IzMWO6+3SVM6LAlr1Q0ScA9wKfGQNau1aq77IzO83bXrAm4EN/Wbmhpq6AbZk5onm/rPAlmXajPP1CzcAD2TmRv7lTdwXEXFes3xXROwC/h24JTM36t9vmcZ58acRsQQ8yPC/8Bvx3GjbD3cBHwNe7KzCtdP6nIiIRxiO9P8G+MuO6lwT6y7oI+Ix4MeW2XTH6EJmDiJi0ifjDcB7JnzsmumwL3oM36P4x8y8NSJuBT7KOu6Tjs+Ld2fmseZ9nAcZ9sO9k1Xara76ISIuA34iM38zIra3q3JtdJ0VmfnOiPhh4NPALzAc8W9I6y7oM/OqN9oWEScjYmtmnoiIrcByc6nHGL7BeNo24PGRfVwK9DLzq9OpuDsd9sV3GY7aTr9H8RfATdOouStdnheZeay5fSEi/pzhKG5dBn2H/fAzwM6IOMIwFy6IiMczcxfrVNdZ0RzjfyLiIYZTQQb9GlkA9gIHmtuHlmnzCPC7I2+qXA18cGT7jcB9XRa5Ribui2aE89cMT/K/A64EvtV5xd2ZuC8iogecl5mLEfEm4DrgsTWouQttzonngHsAmhH9Z9dzyI+hzTlxDnBu8yLRA36Rjf1+3oZ7M/YA8I6I+DZwVbNMROyMiE8ANCfsXQy/b+fLwO80604LagR92774LeDDEfEvDKcqblvj+qepTV+8BXik6YevMxzlHVz7f8JUTOP5UUWbvngrsDByTpwC/mjt/wnT4ydjJam4jTailyStkkEvScUZ9JJUnEEvScUZ9JJUnEEvScUZ9JJUnEEvScX9Py86h8/Sta1QAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(samples2);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This distribution is much more concentrated. CLT is in work"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# mean propagation\n",
    "This is a bit more tricky. The procedure is not mathematically correct in sence that it may give unexpected results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [],
   "source": [
    "mean_theta = approx.histogram.mean(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {},
   "outputs": [],
   "source": [
    "out_prop = theano.clone(approx.to_flat_input(out), {approx.input: mean_theta})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(0.4845134, dtype=float32)"
      ]
     },
     "execution_count": 96,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "out_prop.eval()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x114ad7da0>"
      ]
     },
     "execution_count": 98,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAGBlJREFUeJzt3XtwlPW9x/F3LpRghXLZSklgxlZgWsACioDDTSRHLqGmjvSnHM8p3kjPDJZSqRSOFk+rzNhLQAZbJfU+Q5HvQYVUUQqBFh2goNWi0qOgTQ3hopEICIQmsOePXbZRQ/JkL1n2x+c1k8k+zz77PN8vGz558ttnf5sVDocRERF/Zae7ABERSS0FvYiI5xT0IiKeU9CLiHhOQS8i4jkFvYiI5xT0IiKeU9CLiHhOQS8i4rncdBcQpbfniojEJ6ulDc6WoGfv3r3pLiHpQqEQNTU16S4jJdRbZvK1t0zta8qUbgCsXPnRGbdprrf8/PxAx9HQjYiI5xT0IiKeU9CLiHhOQS8i4jkFvYiI5xT0IiKeU9CLiHhOQS8i4jkFvYiI586ad8amWn5BQez23urqNFYiItK2zskz+vyCgk8Fv4iIz87JoBcROZco6EVEPKegFxHxnIJeRMRzCnoREc8p6EVEPKegFxHxnIJeRMRzCnoREc8p6EVEPKegFxHxnIJeRMRzCnoREc8p6EVEPKegFxHxnIJeRMRzCnoREc8p6EVEPKegFxHxnIJeRMRzCnoREc8p6EVEPKegFxHxXG5LGzjnegFPAt2BMFBmZoudc12BFcCFQCXgzKzWOZcFLAYmAceAG83sL6kpX0REWhLkjL4BmG1m/YDhwAznXD9gLlBhZn2AiugywESgT/SrBHgw6VWLiEhgLQa9me07fUZuZkeAvwEFQDHwRHSzJ4BvR28XA0+aWdjMtgKdnXM9kl65iIgE0uLQTWPOuQuBwcCfge5mti96134iQzsQ+SVQ1ehhe6Lr9jVah3OuhMgZP2ZGKBRqbe0JS/Uxc3Nz09JXW1BvmcnX3jK1r3btIhHcXO3J6C1w0DvnzgeeBmaZ2WHnXOw+Mws758KtObCZlQFl0cVwTU1Nax7eavlNrEv1MUOhUMqPkS7qLTP52lum9lVf3w2AmpqPzrhNc73l5zeVbJ8X6Kob51w7IiG/zMyeia4+cHpIJvr9g+j6aqBXo4f3jK4TEZE0CHLVTRbwCPA3M1vY6K5yYBpwX/T76kbrb3POPQUMAw41GuIREZE2FmToZgTwn8AbzrnXo+v+m0jAm3PuFuAfwOmxnDVELq3cTeTyypuSWrGIiLRKi0FvZi8DWWe4e1wT24eBGQnWJSIiSaJ3xoqIeE5BLyLiOQW9iIjnFPQiIp5T0IuIeE5BLyLiOQW9iIjnFPQiIp5T0IuIeE5BLyLiOQW9iIjnFPQiIp5T0IuIeE5BLyLiOQW9iIjnFPQiIp5T0IuIeE5BLyLiOQW9iIjnFPQiIp5T0IuIeC433QWkU35BQez23urqNFYiIpI6OqMXEfGcgl5ExHMKehERzynoRUQ8p6AXEfGcgl5ExHMKehERzynoRUQ8p6AXEfGcgl5ExHMKehERzynoRUQ8p6AXEfFci7NXOuceBSYDH5jZgOi6/wGmAx9GN/tvM1sTvW8ecAtwEphpZmtTULeIiAQUZJrix4EHgCc/s36Rmf2q8QrnXD/geqA/kA+sd871NbOTSahVRETi0OLQjZltAg4G3F8x8JSZnTCzvwO7gaEJ1CciIglK5INHbnPOfRd4BZhtZrVAAbC10TZ7outERCRN4g36B4F7gHD0eylwc2t24JwrAUoAzIxQKBRnKcmRiuPn5uamva9UUW+ZydfeMrWvdu0iEdxc7cnoLa6gN7MDp287534LPBddrAZ6Ndq0Z3RdU/soA8qii+Gampp4Sgksv4X7U3H8UCiUkv2eDdRbZvK1t0ztq76+GwA1NR+dcZvmesvPbynZIuK6vNI516PR4jXAm9Hb5cD1zrn2zrmvAn2AbfEcQ0REkiPI5ZXLgSuAkHNuD3A3cIVzbhCRoZtK4HsAZvaWc86AnUADMENX3IiIpFeLQW9mU5tY/Ugz2y8AFiRSlIiIJI/eGSsi4rlELq/MCPkFurpTRM5tOqMXEfGcgl5ExHMKehERzynoRUQ8p6AXEfGcgl5ExHMKehERzynoRUQ8p6AXEfGcgl5ExHMKehERzynoRUQ8p6AXEfGcgl5ExHMKehERzynoRUQ8p6AXEfGcgl5ExHMKehERzynoRUQ8p6AXEfGcgl5ExHMKehERzynoRUQ8p6AXEfGcgl5ExHMKehERzynoRUQ8p6AXEfGcgl5ExHMKehERzynoRUQ8p6AXEfGcgl5ExHO5LW3gnHsUmAx8YGYDouu6AiuAC4FKwJlZrXMuC1gMTAKOATea2V9SU7qIiAQR5Iz+cWDCZ9bNBSrMrA9QEV0GmAj0iX6VAA8mp0wREYlXi0FvZpuAg59ZXQw8Eb39BPDtRuufNLOwmW0FOjvneiSrWBERab14x+i7m9m+6O39QPfo7QKgqtF2e6LrREQkTVoco2+JmYWdc+HWPs45V0JkeAczIxQKJVpKQlJx/Nzc3LT3lSrqLTP52lum9tWuXSSCm6s9Gb3FG/QHnHM9zGxfdGjmg+j6aqBXo+16Rtd9jpmVAWXRxXBNTU2cpTQvP+B2qTh+KBRKyX7PBuotM/naW6b2VV/fDYCamo/OuE1zveXnB0u4eIO+HJgG3Bf9vrrR+tucc08Bw4BDjYZ4REQkDYJcXrkcuAIIOef2AHcTCXhzzt0C/ANw0c3XELm0cjeRyytvSkHNIiLSCi0GvZlNPcNd45rYNgzMSLQoERFJHr0zVkTEcwp6ERHPKeij8gsKyC/QJf8i4h8FvYiI5xT0IiKeU9CLiHhOQS8i4jkFvYiI5xT0IiKeU9CLiHhOQS8i4jkFvYiI5xT0IiKeU9CLiHhOQS8i4jkFvYiI5xT0IiKeU9CLiHhOQS8i4jkFvYiI5xT0IiKeU9CLiHhOQS8i4jkFvYiI5xT0IiKeU9CLiHhOQS8i4jkFvYiI5xT0IiKeU9CLiHhOQS8i4rncdBeQCvkFBekuQUTkrKEzehERzynoRUQ8p6AXEfGcgl5ExHMJvRjrnKsEjgAngQYzG+Kc6wqsAC4EKgFnZrWJlSkiIvFKxhn9WDMbZGZDostzgQoz6wNURJdFRCRNUjF0Uww8Eb39BPDtFBxDREQCSvQ6+jDwB+dcGFhqZmVAdzPbF71/P9C9qQc650qAEgAzIxQKJVhKciSzjtzc3LOmr2RTb5nJ194yta927SIR3Fztyegt0aAfaWbVzrkLgHXOuf9rfKeZhaO/BD4n+kuhLLoYrqmpSbCUf8lP4LHJrCMUCiV1f2cT9ZaZfO0tU/uqr+8GQE3NR2fcprne8vODpV1CQzdmVh39/gHwLDAUOOCc6wEQ/f5BIscQEZHExB30zrkvOuc6nr4NXAW8CZQD06KbTQNWJ1qkiIjEL5Ez+u7Ay865vwLbgOfN7EXgPuDfnHO7gMLosoiIpEncY/Rm9h4wsIn1HwHjEilKRESSR++MFRHxnJfTFEvqhMNh6urqOHDgACdOnEh3OSlxrvYWDofJzs4mLy+PrKysNq5MUklBL61SV1dHu3btaN++PTk5OekuJyVyc3PP2d4aGhqoq6ujQ4cObViVpJqGbqRVTp06RW6uzg98lZuby6lTp9JdhiSZgl5aRX/S+0/PsX8U9CIinlPQiyRZQ0NDuksQ+RQNtkpGqaqq4oYbbuCSSy7hlVdeYdCgQTjnKC0tpaamhgceeIDBgwdz7Ngx7rrrLt5++23q6+uZPXs248ePp6qqipkzZ3Ls2DEA7r33Xi677DI2b97MwoUL6dKlC++88w4XX3wxS5Ys+dwwxpQpU+jXrx9btmzh5MmTlJaWMnjwYEpLS6msrOT999+noKCAhQsXMm/ePHbs2EFOTg533303I0aMYMWKFbz44oscOXKEffv2ce2113L77ben459SziEK+s/ILyiI3d5bXZ3GSs5+8+d3YufOdkndZ79+9fzsZ4eb3aayspKlS5eycOFCJk2axKpVq1i1ahV/+MMfWLJkCY8++iiLFy9mxIgRLFy4kEOHDlFUVMSoUaMIhUIsX76cvLw83nvvPWbMmMELL7wAwJtvvsmGDRvo2bMnRUVFbN++naFDh37u+MePH2fdunVs3bqV2bNns2HDBgB27drFs88+S4cOHXjooYfIysqioqKC3bt3M3XqVF566SUAXn/9dSoqKujQoQNFRUWMGzeOgQM/995DkaRR0EvG6dWrF9/4xjcA6Nu3LyNHjiQrK4uvf/3rVFVVAbBp0ybWrVvHQw89BMCJEyeorq6me/fu3HnnnezcuZPs7Gzee++92H4HDRpEfn4+2dnZ9O/fn6qqqiaDvri4GIDhw4dz5MgRDh06BMBVV10Vuyxx+/bt3HTTTQD07t2bnj17xo41atQounbtCsDEiRPZtm2bgl5SSkEvcWvpzDtV2rdvH7udnZ3NF77whdjtkydPApE3/5SVldG7d+9PPba0tJQvf/nLrFu3jlOnTvG1r30tdt/p/QDk5OSccaz9s8M5p5fPO++8QPWf6fEiqaIXY8VLY8aM4bHHHiMcjnwcwptvvgnA4cOHueCCC8jOzubpp5+O/WJojfLycgC2bdtGp06d6NSp0+e2GTp0KM8++ywA7777LtXV1Vx00UUAvPTSS9TW1nL8+HHWrl3LZZddFlePIkHpjF68NGvWLO6++24KCws5deoUvXr14sknn2TatGmUlJSwcuVKxo4dG/gsvLH27dtz1VVX0dDQQGlpaZPbTJs2jXnz5jFu3DhycnJYtGhR7C+RQYMGMX369NiLsRq2kVTLOn3Gk2bhvXv3Jm1njV9QTUSiL8Zm6qfeNOfYsWOcd9555ObmensZYXO9TZkyhZ/85Cdxh/OKFSvYsWMHCxYsSKTEuAV53k4/x5kkU/+vTZkS+YSplSsT+oSpFsf+NHQjIuI5Dd2ItMLKlSsTevx1113Hddddl6RqRILRGb2IiOcU9CIinvNq6CZZL8KKiPhEZ/QiIp5T0Ms5bdiwYRw8eDDhbdJp//79TJ8+Pd1lyFnMq6EbaXvJHi471yeSO3nyZKs/xvArX/kKv/3tb1NUkfhAZ/SSUaqqqhg9ejSzZs1i5MiR3HbbbWzatIni4mJGjBjBa6+9BkBtbS0333wzhYWFTJ48mZ07dwJw8OBBpk6dytixY/nRj35E4zcMPv300xQVFXHllVcyZ86cFqdHmDt3LhMnTmTs2LH86le/AmDjxo2UlJTEttm8eTPf/e53AfjTn/7Et771LcaPH09JSQlHjx4FIn8xLFiwgPHjx/Pcc8+xbNkyJk2aRGFhIdOnT+f48eNAZNbOyZMnM27cOH7+85/Tp0+f2L/JlVdeCUTekHXrrbdyww03MGLECO69995YLcuXL+fyyy+nqKiIO+64gzvvvDP+J0IyioJeMk5lZSXf+9732LRpE7t3745NUzx//nyWLFkCRCYvGzBgAOvXr2fu3Ln84Ac/AGDRokUMHTqUjRs3MmHCBKqjf0Hs2rWL8vJyVq1axYYNG8jJyeGZZ55pto4f//jHvPDCC6xfv56tW7eyc+dORo0axWuvvRab7768vJzi4mIOHjzI4sWLWbFiBWvXrmXgwIGUlZXF9tWlSxfWrl1LcXExEydOZM2aNaxfv57evXuzfPlyAObPn8+tt95KRUUFPXr0OGNdb731Fg8++CAVFRWUl5dTXV3N/v37uf/++1mzZg2rVq1i9+7d8T8BknE0dCMZJ8g0xdu2bYsNZ4wcOZLa2lqOHDnC1q1befjhhwEoLCykc+fOALz88su88cYbTJo0iaysLI4fP04oFGq2jt///vcsW7aMkydPcuDAAXbt2kW/fv0YO3Ys69ato6ioiIqKCu666y62bNnCO++8E5viuL6+nksvvTS2r6uvvjp2++233+YXv/gFhw8f5ujRo4wZMwaAV199lUcffRSAa665hnvuuafJukaOHBmbaK1v375UV1dz8OBBhg8fTpcuXWhoaGDy5MmfmqJZ/Kagl4wTZJri1gqHw3znO99h3rx5geaDef/991m6dCnPP/88nTt3ZtasWdTV1QGR0H788cfp3LkzAwcO5PzzzyccDjN69Gh+85vfNLm/xnPL/PCHP+SRRx6hf//+rFixgi1btrSql8bTLWdnZ3s7J5EEp6Eb8dKwYcNiQy+bN2+ma9eudOzYkeHDh8emD96wYQMff/wxEDkLfu6552KTR9XW1rJnz54z7v/IkSN06NCBTp068eGHH7Jx48bYfZdffjlvvPEGy5Yti52pX3rppWzfvp2///3vQGTisHfffbfJfX/yySd0796d+vr6WK0Al1xyCc8//zwAq1evbtW/x8CBA9m6dSsff/wxDQ0NrFmzplWPl8ymM3rx0u23387s2bMpLCwkLy+P+++/H4icLc+YMYOxY8cyZMgQCqJXDfXt25c5c+YwdepUwuEwubm5LFiwgJ49eza5//79+zNgwABGjx5Nfn7+p+aUz8nJobCwEDNj8eLFAHTr1o1FixYxY8YM/vnPfwIwZ86c2Bz1jd1xxx1MnjyZbt26MXjwYD755BMAfvrTnzJz5kyWLFnCFVdc0eQ8+GfSo0cPvv/97zNhwgQ6d+7MRRddRMeOHQM/XjKbV9MUn22X+mXq1KnNOdenKU6n48ePk5eXR1ZWFqtXr2bVqlU89thjgR9/9OhRvvSlL1FXV8ctt9zC9ddfz8SJEz+3naYpbjttNU2xzuhFMsSOHTtil0R26tTpjB96cialpaW8/PLL1NXVMWbMGCZMmJCKMuUspKAXyRDDhg1j/fr1cT9+/vz5Z+1fK5JaejFWRMRzGX9Grxkr29ZZ8pqOpJCeY//ojL4Z+QUF+kXyGbou228NDQ1kZysWfJPxZ/TStvLy8qirqyMrK4sTJ06ku5yUaN++/TnZWzgcJjs7m7y8vDauSlItZUHvnJsALAZygIfN7L5UHUvaTlZWFh06dMjYy9mCUG/im5T8jeacywF+DUwE+gFTnXP9UnEsERFpXqoG44YCu83sPTP7J/AUUJyiY4mISDNSFfQFQFWj5T3RdSIi0sbS9mKsc64EKAEws9Nv5W29NrgULM7KIo+Nt68MoN4yk6+9ZWJfmzefvtV87Yn2lqoz+mqgV6PlntF1MWZWZmZDzGwIkbkavPtyzr2a7hrUm3o7F3rzta+AvbUoVWf024E+zrmvEgn464F/T9GxRESkGSk5ozezBuA2YC3wt8gqeysVxxIRkealbIzezNYA5/qnG5S1vEnGUm+ZydfefO0LktDb2TIfvYiIpIgmtRAR8ZzmukmClqZ7cM61B54ELgU+Aq4zs8q2rjMeAXobDdwPfBO43sxWtn2VrRegr9uBW4EG4EPgZjP7R5sXGocAvf0XMAM4CXwClJjZzjYvNA5Bp1Zxzl0LrAQuM7NX2rDEuAV43m4Efsm/rmB8wMweDrJvndEnKOB0D7cAtWbWG1gE/Lxtq4xPwN7eB24Efte21cUvYF+vAUPM7JtEAuMXbVtlfAL29jszu9jMBhHpa2EblxmXoFOrOOc6Aj8A/ty2FcavFdPGrDCzQdGvQCEPCvpkCDLdQzHwRPT2SmCccy7Q9a9p1mJvZlZpZjuAU+koME5B+tpoZseii1uJvBckEwTp7XCjxS8CmfJCXdCpVe4hcjJV15bFJSil08Yo6BMXZLqH2DbRS08PAd3apLrE+DqVRWv7ugV4IaUVJU+g3pxzM5xz7xI5o5/ZRrUlqsXenHOXAL3M7Pm2LCwJgv5MXuuc2+GcW+mc69XE/U1S0Is0wzn3H8AQImOj3jCzX5vZRcCPgbvSXU8yOOeyiQxDzU53LSnye+DC6HDiOv41StAiBX3iWpzuofE2zrlc4EtEXpQ92wXpLRMF6ss5VwjcCVxtZpnySSStfc6eAr6d0oqSp6XeOgIDgD865yqB4UC5c25Im1UYvyDTxnzU6OfwYSIXdwSiq24SF2S6h3JgGrAFmAJsMLNMGBf1dSqLFvtyzg0GlgITzOyDti8xbkF662Nmu6KLRcAuMkOzvZnZISB0etk590fgRxly1U2Q562Hme2LLl5NZNaBQHRGn6AzTffgnPuZc+7q6GaPAN2cc7uB24G56am2dYL05py7zDm3B/gOsNQ5d9ZPdRHwOfslcD7wv865151z5Wkqt1UC9nabc+4t59zrRH4ep6Wp3FYJ2FtGCtjbzOjz9lcir6vcGHT/emesiIjndEYvIuI5Bb2IiOcU9CIinlPQi4h4TkEvIuI5Bb2IiOcU9CIinlPQi4h47v8BXPjrB5X59QAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(samples2, label='model averaging', color='r');\n",
    "plt.axvline(out_prop.eval(), label='mean prop', color='b');\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are done, hope I gave the full answer"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "pymc4",
   "language": "python",
   "name": "pymc4"
  },
  "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}