pangolin-jags/Example-Simulation-Based-Calibration.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from pangolin import d, t, Given, I, IID, vmap, recurse, scan, sample, E, P, var, std, cov, corr, jags_code, makerv\n",
    "from matplotlib import pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "outputs": [
    {
     "data": {
      "text/plain": "<Figure size 640x480 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGdCAYAAADAAnMpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAeKUlEQVR4nO3df3DX9X3A8VcI5ht7JQFHSYClS3W12qnAoGTRetVdbE45Wv/YyWkPGOeP2bKeI7ebxB+k1pYwTx27GctJ5ezuqqCedV3hsC4b51mz4wzmzq2iR4HBrIlwrgnGNZHksz8841IC5huTvEl4PO4+f+ST9+fzeX/7LuZ5n++vgizLsgAASGRK6gkAAGc2MQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAElNTT2B4ejv749f//rXMW3atCgoKEg9HQBgGLIsi2PHjsWcOXNiypST3/+YEDHy61//OioqKlJPAwAYgcOHD8fv//7vn/T3EyJGpk2bFhEfPJiSkpLEswEAhqOrqysqKioG/o6fzISIkQ+fmikpKREjADDBfNxLLLyAFQBISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQ1NTUE0itcu321FMYFQc3LEk9BQAYEXdGAICkxAgAkNQZ/zQNp5fJ8rRZhKfOAIbLnREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASfnQM+CUJssH0fkQOjh9uTMCACSVd4y88MILsXTp0pgzZ04UFBTEs88+e8rxzzzzTFx11VXxmc98JkpKSqK6ujqee+65kc4XAJhk8o6R7u7umDdvXjQ1NQ1r/AsvvBBXXXVV7NixI1pbW+PKK6+MpUuXxiuvvJL3ZAGAySfv14xcffXVcfXVVw97/MaNGwf9vH79+vinf/qn+Od//udYsGBBvpcHACaZcX8Ba39/fxw7dizOOeeck47p6emJnp6egZ+7urrGY2oAQALj/gLW+++/P95999247rrrTjqmsbExSktLB7aKiopxnCEAMJ7GNUYef/zxuOeee+LJJ5+MWbNmnXRcfX19dHZ2DmyHDx8ex1kCAONp3J6m2bp1a9x0003x1FNPRU1NzSnH5nK5yOVy4zQzACClcbkz8sQTT8SqVaviiSeeiCVLfPAQAPCRvO+MvPvuu7Fv376Bnw8cOBBtbW1xzjnnxGc/+9mor6+PN998M/7xH/8xIj54amblypXx93//91FVVRXt7e0REXH22WdHaWnpKD0MgFPzSbJw+sr7zsjLL78cCxYsGHhbbl1dXSxYsCDWrVsXERFvvfVWHDp0aGD8I488EsePH4/Vq1fH7NmzB7bbbrttlB4CADCR5X1n5Iorrogsy076+8cee2zQz7t27cr3EgDAGcR30wAASYkRACApMQIAJDXuHwfP2Jgs7xQA4MzjzggAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAIKmpqScAABNV5drtqacwKg5uWJL0+u6MAABJiREAIKm8Y+SFF16IpUuXxpw5c6KgoCCeffbZjz1m165d8cd//MeRy+XiD//wD+Oxxx4bwVQBgMko7xjp7u6OefPmRVNT07DGHzhwIJYsWRJXXnlltLW1xV/91V/FTTfdFM8991zekwUAJp+8X8B69dVXx9VXXz3s8Zs2bYrPfe5z8cADD0RExIUXXhgvvvhi/N3f/V3U1tbme3kAYJIZ89eMtLS0RE1NzaB9tbW10dLSctJjenp6oqura9AGAExOYx4j7e3tUVZWNmhfWVlZdHV1xf/+7/8OeUxjY2OUlpYObBUVFWM9TQAgkdPy3TT19fXR2dk5sB0+fDj1lACAMTLmH3pWXl4eHR0dg/Z1dHRESUlJnH322UMek8vlIpfLjfXUAEhksnxYGKNjzO+MVFdXR3Nz86B9zz//fFRXV4/1pQGACSDvGHn33Xejra0t2traIuKDt+62tbXFoUOHIuKDp1hWrFgxMP7WW2+N/fv3x9/8zd/E3r174+GHH44nn3wy1qxZMzqPAACY0PKOkZdffjkWLFgQCxYsiIiIurq6WLBgQaxbty4iIt56662BMImI+NznPhfbt2+P559/PubNmxcPPPBA/PCHP/S2XgAgIkbwmpErrrgisiw76e+H+nTVK664Il555ZV8LwUAnAFOy3fTAABnDjECACQlRgCApMQIAJDUmH/oGQCjx4eFMRm5MwIAJCVGAICkxAgAkJQYAQCSEiMAQFLeTQNjxLseAIbHnREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACCpEcVIU1NTVFZWRnFxcVRVVcXu3btPOX7jxo3xhS98Ic4+++yoqKiINWvWxG9/+9sRTRgAmFzyjpFt27ZFXV1dNDQ0xJ49e2LevHlRW1sbb7/99pDjH3/88Vi7dm00NDTEa6+9Fo8++mhs27Yt7rjjjk88eQBg4ss7Rh588MG4+eabY9WqVfHFL34xNm3aFJ/61Kdiy5YtQ45/6aWX4rLLLosbbrghKisr46tf/Wpcf/31H3s3BQA4M+QVI729vdHa2ho1NTUfnWDKlKipqYmWlpYhj7n00kujtbV1ID72798fO3bsiGuuueak1+np6Ymurq5BGwAwOU3NZ/DRo0ejr68vysrKBu0vKyuLvXv3DnnMDTfcEEePHo0vf/nLkWVZHD9+PG699dZTPk3T2NgY99xzTz5TAwAmqDF/N82uXbti/fr18fDDD8eePXvimWeeie3bt8e999570mPq6+ujs7NzYDt8+PBYTxMASCSvOyMzZ86MwsLC6OjoGLS/o6MjysvLhzzm7rvvjuXLl8dNN90UEREXX3xxdHd3xy233BJ33nlnTJlyYg/lcrnI5XL5TA0AmKDyujNSVFQUCxcujObm5oF9/f390dzcHNXV1UMe8957750QHIWFhRERkWVZvvMFACaZvO6MRETU1dXFypUrY9GiRbF48eLYuHFjdHd3x6pVqyIiYsWKFTF37txobGyMiIilS5fGgw8+GAsWLIiqqqrYt29f3H333bF06dKBKAEAzlx5x8iyZcviyJEjsW7dumhvb4/58+fHzp07B17UeujQoUF3Qu66664oKCiIu+66K9588834zGc+E0uXLo3vf//7o/coAIAJqyCbAM+VdHV1RWlpaXR2dkZJScmonrty7fZRPR8ATDQHNywZk/MO9++376YBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAIKkRxUhTU1NUVlZGcXFxVFVVxe7du085/je/+U2sXr06Zs+eHblcLs4///zYsWPHiCYMAEwuU/M9YNu2bVFXVxebNm2Kqqqq2LhxY9TW1sbrr78es2bNOmF8b29vXHXVVTFr1qx4+umnY+7cufFf//VfMX369NGYPwAwweUdIw8++GDcfPPNsWrVqoiI2LRpU2zfvj22bNkSa9euPWH8li1b4p133omXXnopzjrrrIiIqKys/GSzBgAmjbyepunt7Y3W1taoqan56ARTpkRNTU20tLQMecxPf/rTqK6ujtWrV0dZWVlcdNFFsX79+ujr6zvpdXp6eqKrq2vQBgBMTnnFyNGjR6Ovry/KysoG7S8rK4v29vYhj9m/f388/fTT0dfXFzt27Ii77747Hnjggfje97530us0NjZGaWnpwFZRUZHPNAGACWTM303T398fs2bNikceeSQWLlwYy5YtizvvvDM2bdp00mPq6+ujs7NzYDt8+PBYTxMASCSv14zMnDkzCgsLo6OjY9D+jo6OKC8vH/KY2bNnx1lnnRWFhYUD+y688MJob2+P3t7eKCoqOuGYXC4XuVwun6kBABNUXndGioqKYuHChdHc3Dywr7+/P5qbm6O6unrIYy677LLYt29f9Pf3D+x74403Yvbs2UOGCABwZsn7aZq6urrYvHlz/OhHP4rXXnstvvnNb0Z3d/fAu2tWrFgR9fX1A+O/+c1vxjvvvBO33XZbvPHGG7F9+/ZYv359rF69evQeBQAwYeX91t5ly5bFkSNHYt26ddHe3h7z58+PnTt3Dryo9dChQzFlykeNU1FREc8991ysWbMmLrnkkpg7d27cdtttcfvtt4/eowAAJqyCLMuy1JP4OF1dXVFaWhqdnZ1RUlIyqueuXLt9VM8HABPNwQ1LxuS8w/377btpAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhqRDHS1NQUlZWVUVxcHFVVVbF79+5hHbd169YoKCiIa6+9diSXBQAmobxjZNu2bVFXVxcNDQ2xZ8+emDdvXtTW1sbbb799yuMOHjwYf/3Xfx2XX375iCcLAEw+ecfIgw8+GDfffHOsWrUqvvjFL8amTZviU5/6VGzZsuWkx/T19cU3vvGNuOeee+Lcc8/9RBMGACaXvGKkt7c3Wltbo6am5qMTTJkSNTU10dLSctLjvvvd78asWbPixhtvHNZ1enp6oqura9AGAExOecXI0aNHo6+vL8rKygbtLysri/b29iGPefHFF+PRRx+NzZs3D/s6jY2NUVpaOrBVVFTkM00AYAIZ03fTHDt2LJYvXx6bN2+OmTNnDvu4+vr66OzsHNgOHz48hrMEAFKams/gmTNnRmFhYXR0dAza39HREeXl5SeM/9WvfhUHDx6MpUuXDuzr7+//4MJTp8brr78e55133gnH5XK5yOVy+UwNAJig8rozUlRUFAsXLozm5uaBff39/dHc3BzV1dUnjL/gggvi1Vdfjba2toHta1/7Wlx55ZXR1tbm6RcAIL87IxERdXV1sXLlyli0aFEsXrw4Nm7cGN3d3bFq1aqIiFixYkXMnTs3Ghsbo7i4OC666KJBx0+fPj0i4oT9AMCZKe8YWbZsWRw5ciTWrVsX7e3tMX/+/Ni5c+fAi1oPHToUU6b4YFcAYHgKsizLUk/i43R1dUVpaWl0dnZGSUnJqJ67cu32UT0fAEw0BzcsGZPzDvfvt1sYAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgqRHFSFNTU1RWVkZxcXFUVVXF7t27Tzp28+bNcfnll8eMGTNixowZUVNTc8rxAMCZJe8Y2bZtW9TV1UVDQ0Ps2bMn5s2bF7W1tfH2228POX7Xrl1x/fXXx7/9279FS0tLVFRUxFe/+tV48803P/HkAYCJryDLsiyfA6qqquJLX/pSPPTQQxER0d/fHxUVFfHtb3871q5d+7HH9/X1xYwZM+Khhx6KFStWDOuaXV1dUVpaGp2dnVFSUpLPdD9W5drto3o+AJhoDm5YMibnHe7f77zujPT29kZra2vU1NR8dIIpU6KmpiZaWlqGdY733nsv3n///TjnnHPyuTQAMElNzWfw0aNHo6+vL8rKygbtLysri7179w7rHLfffnvMmTNnUND8rp6enujp6Rn4uaurK59pAgATyLi+m2bDhg2xdevW+MlPfhLFxcUnHdfY2BilpaUDW0VFxTjOEgAYT3nFyMyZM6OwsDA6OjoG7e/o6Ijy8vJTHnv//ffHhg0b4uc//3lccsklpxxbX18fnZ2dA9vhw4fzmSYAMIHkFSNFRUWxcOHCaG5uHtjX398fzc3NUV1dfdLj7rvvvrj33ntj586dsWjRoo+9Ti6Xi5KSkkEbADA55fWakYiIurq6WLlyZSxatCgWL14cGzdujO7u7li1alVERKxYsSLmzp0bjY2NERHxt3/7t7Fu3bp4/PHHo7KyMtrb2yMi4tOf/nR8+tOfHsWHAgBMRHnHyLJly+LIkSOxbt26aG9vj/nz58fOnTsHXtR66NChmDLloxsuP/jBD6K3tzf+7M/+bNB5Ghoa4jvf+c4nmz0AMOHl/TkjKficEQAYOxPqc0YAAEabGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmNKEaampqisrIyiouLo6qqKnbv3n3K8U899VRccMEFUVxcHBdffHHs2LFjRJMFACafvGNk27ZtUVdXFw0NDbFnz56YN29e1NbWxttvvz3k+Jdeeimuv/76uPHGG+OVV16Ja6+9Nq699tr4j//4j088eQBg4ivIsizL54Cqqqr40pe+FA899FBERPT390dFRUV8+9vfjrVr154wftmyZdHd3R0/+9nPBvb9yZ/8ScyfPz82bdo0rGt2dXVFaWlpdHZ2RklJST7T/ViVa7eP6vkAYKI5uGHJmJx3uH+/p+Zz0t7e3mhtbY36+vqBfVOmTImamppoaWkZ8piWlpaoq6sbtK+2tjaeffbZk16np6cnenp6Bn7u7OyMiA8e1Gjr73lv1M8JABPJWPx9/f/n/bj7HnnFyNGjR6Ovry/KysoG7S8rK4u9e/cOeUx7e/uQ49vb2096ncbGxrjnnntO2F9RUZHPdAGAYSjdOLbnP3bsWJSWlp7093nFyHipr68fdDelv78/3nnnnfi93/u9KCgoGLXrdHV1RUVFRRw+fHjUn/4hf9bj9GNNTi/W4/RiPT5elmVx7NixmDNnzinH5RUjM2fOjMLCwujo6Bi0v6OjI8rLy4c8pry8PK/xERG5XC5yudygfdOnT89nqnkpKSnxf6TTiPU4/ViT04v1OL1Yj1M71R2RD+X1bpqioqJYuHBhNDc3D+zr7++P5ubmqK6uHvKY6urqQeMjIp5//vmTjgcAzix5P01TV1cXK1eujEWLFsXixYtj48aN0d3dHatWrYqIiBUrVsTcuXOjsbExIiJuu+22+MpXvhIPPPBALFmyJLZu3Rovv/xyPPLII6P7SACACSnvGFm2bFkcOXIk1q1bF+3t7TF//vzYuXPnwItUDx06FFOmfHTD5dJLL43HH3887rrrrrjjjjvi85//fDz77LNx0UUXjd6jGKFcLhcNDQ0nPCVEGtbj9GNNTi/W4/RiPUZP3p8zAgAwmnw3DQCQlBgBAJISIwBAUmIEAEhq0sdIU1NTVFZWRnFxcVRVVcXu3btPOf6pp56KCy64IIqLi+Piiy+OHTt2jNNMzwz5rMfmzZvj8ssvjxkzZsSMGTOipqbmY9eP/OX7b+RDW7dujYKCgrj22mvHdoJnmHzX4ze/+U2sXr06Zs+eHblcLs4//3z/3RpF+a7Hxo0b4wtf+EKcffbZUVFREWvWrInf/va34zTbCSybxLZu3ZoVFRVlW7Zsyf7zP/8zu/nmm7Pp06dnHR0dQ47/xS9+kRUWFmb33Xdf9stf/jK76667srPOOit79dVXx3nmk1O+63HDDTdkTU1N2SuvvJK99tpr2Z//+Z9npaWl2X//93+P88wnr3zX5EMHDhzI5s6dm11++eXZ17/+9fGZ7Bkg3/Xo6enJFi1alF1zzTXZiy++mB04cCDbtWtX1tbWNs4zn5zyXY8f//jHWS6Xy3784x9nBw4cyJ577rls9uzZ2Zo1a8Z55hPPpI6RxYsXZ6tXrx74ua+vL5szZ07W2Ng45PjrrrsuW7JkyaB9VVVV2V/8xV+M6TzPFPmux+86fvx4Nm3atOxHP/rRWE3xjDOSNTl+/Hh26aWXZj/84Q+zlStXipFRlO96/OAHP8jOPffcrLe3d7ymeEbJdz1Wr16d/emf/umgfXV1ddlll102pvOcDCbt0zS9vb3R2toaNTU1A/umTJkSNTU10dLSMuQxLS0tg8ZHRNTW1p50PMM3kvX4Xe+99168//77cc4554zVNM8oI12T7373uzFr1qy48cYbx2OaZ4yRrMdPf/rTqK6ujtWrV0dZWVlcdNFFsX79+ujr6xuvaU9aI1mPSy+9NFpbWweeytm/f3/s2LEjrrnmmnGZ80R2Wn5r72g4evRo9PX1DXwy7IfKyspi7969Qx7T3t4+5Pj29vYxm+eZYiTr8btuv/32mDNnzgnByMiMZE1efPHFePTRR6OtrW0cZnhmGcl67N+/P/71X/81vvGNb8SOHTti37598a1vfSvef//9aGhoGI9pT1ojWY8bbrghjh49Gl/+8pcjy7I4fvx43HrrrXHHHXeMx5QntEl7Z4TJZcOGDbF169b4yU9+EsXFxamnc0Y6duxYLF++PDZv3hwzZ85MPR3igy8qnTVrVjzyyCOxcOHCWLZsWdx5552xadOm1FM7I+3atSvWr18fDz/8cOzZsyeeeeaZ2L59e9x7772pp3bam7R3RmbOnBmFhYXR0dExaH9HR0eUl5cPeUx5eXle4xm+kazHh+6///7YsGFD/Mu//EtccsklYznNM0q+a/KrX/0qDh48GEuXLh3Y19/fHxERU6dOjddffz3OO++8sZ30JDaSfyOzZ8+Os846KwoLCwf2XXjhhdHe3h69vb1RVFQ0pnOezEayHnfffXcsX748brrppoiIuPjii6O7uztuueWWuPPOOwd9bxuDTdr/ZYqKimLhwoXR3Nw8sK+/vz+am5ujurp6yGOqq6sHjY+IeP755086nuEbyXpERNx3331x7733xs6dO2PRokXjMdUzRr5rcsEFF8Srr74abW1tA9vXvva1uPLKK6OtrS0qKirGc/qTzkj+jVx22WWxb9++gSiMiHjjjTdi9uzZQuQTGsl6vPfeeycEx4ehmPkauFNL/QrasbR169Ysl8tljz32WPbLX/4yu+WWW7Lp06dn7e3tWZZl2fLly7O1a9cOjP/FL36RTZ06Nbv//vuz1157LWtoaPDW3lGU73ps2LAhKyoqyp5++unsrbfeGtiOHTuW6iFMOvmuye/ybprRle96HDp0KJs2bVr2l3/5l9nrr7+e/exnP8tmzZqVfe9730v1ECaVfNejoaEhmzZtWvbEE09k+/fvz37+859n5513XnbdddeleggTxqSOkSzLsn/4h3/IPvvZz2ZFRUXZ4sWLs3//938f+N1XvvKVbOXKlYPGP/nkk9n555+fFRUVZX/0R3+Ubd++fZxnPLnlsx5/8Ad/kEXECVtDQ8P4T3wSy/ffyP8nRkZfvuvx0ksvZVVVVVkul8vOPffc7Pvf/352/PjxcZ715JXPerz//vvZd77zney8887LiouLs4qKiuxb3/pW9j//8z/jP/EJpiDL3DsCANKZtK8ZAQAmBjECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQ1P8Bwovn8/EBFLwAAAAASUVORK5CYII=\n"
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "<Figure size 640x480 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGdCAYAAADAAnMpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAeHElEQVR4nO3dbWxd9X3A8Z/t4Guqxk5YGudh7lzYKHRAkiWNZygqTKYWROl4MREBSrKIh1GyisWaRsxD3JQ1zhDQTMM0IiWikwoJVJR1TRRG3UWI4inCwRJbCShNsmQUm0SsdjCrndhnL1BN3Tgh19j+x87nI90XPv6fc3+3p8FfnftUkGVZFgAAiRSmHgAAOLuJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASGpS6gFOR39/f/zyl7+MyZMnR0FBQepxAIDTkGVZHD16NGbNmhWFhSe//jEuYuSXv/xlVFRUpB4DABiGQ4cOxe///u+f9PfjIkYmT54cER8+mNLS0sTTAACno6urKyoqKgb+jp/MuIiR3zw1U1paKkYAYJz5uJdYeAErAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACCpSakHAIDxqnL1ttQjjIgD6xclvX9XRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASU1KPQD8tonyddwR6b+SG2C8cGUEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJn/YeeTZQP2fIBWwCMV66MAABJiREAICkxAgAkJUYAgKTyjpGXXnopFi9eHLNmzYqCgoJ4/vnnT7n+ueeei2uuuSY+85nPRGlpaVRXV8cLL7ww3HkBgAkm7xjp7u6OOXPmRFNT02mtf+mll+Kaa66J7du3R2tra1x99dWxePHieO211/IeFgCYePJ+a++1114b11577Wmv37Bhw6Cf161bF//yL/8S//qv/xrz5s3L9+4BgAlmzD9npL+/P44ePRrnnXfeSdf09PRET0/PwM9dXV1jMRoAkMCYx8hDDz0U77//ftxwww0nXdPY2Bhr164dw6kAxgcf1MhENKbvpnnqqadi7dq18cwzz8T06dNPuq6+vj46OzsHbocOHRrDKQGAsTRmV0a2bNkSt956azz77LNRU1NzyrW5XC5yudwYTQYApDQmV0aefvrpWLFiRTz99NOxaJFLcwDAR/K+MvL+++/H3r17B37ev39/tLW1xXnnnRef/exno76+Pt5+++3453/+54j48KmZ5cuXxz/+4z9GVVVVtLe3R0TEueeeG2VlZSP0MACA8SrvKyOvvvpqzJs3b+BtuXV1dTFv3rxYs2ZNRES88847cfDgwYH1jz/+eBw/fjxWrlwZM2fOHLjdddddI/QQAIDxLO8rI1dddVVkWXbS3z/55JODft65c2e+dwEAnEV8Nw0AkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACCpSakHABgLlau3pR4BOAlXRgCApMQIAJCUp2mAU/L0BjDaXBkBAJISIwBAUp6mmSBcSgdgvHJlBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTyjpGXXnopFi9eHLNmzYqCgoJ4/vnnP3afnTt3xp/8yZ9ELpeLP/zDP4wnn3xyGKMCABNR3jHS3d0dc+bMiaamptNav3///li0aFFcffXV0dbWFn/zN38Tt956a7zwwgt5DwsATDyT8t3h2muvjWuvvfa012/cuDE+97nPxcMPPxwRERdffHG8/PLL8e1vfztqa2vzvXsAYIIZ9deMtLS0RE1NzaBttbW10dLSctJ9enp6oqura9ANAJiYRj1G2tvbo7y8fNC28vLy6Orqiv/7v/8bcp/GxsYoKysbuFVUVIz2mABAImfku2nq6+ujs7Nz4Hbo0KHUIwEAoyTv14zka8aMGdHR0TFoW0dHR5SWlsa555475D65XC5yudxojwYAnAFG/cpIdXV1NDc3D9r24osvRnV19WjfNQAwDuQdI++//360tbVFW1tbRHz41t22trY4ePBgRHz4FMuyZcsG1t9xxx2xb9+++Lu/+7vYs2dPPPbYY/HMM8/EqlWrRuYRAADjWt4x8uqrr8a8efNi3rx5ERFRV1cX8+bNizVr1kRExDvvvDMQJhERn/vc52Lbtm3x4osvxpw5c+Lhhx+O7373u97WCwBExDBeM3LVVVdFlmUn/f1Qn6561VVXxWuvvZbvXQEAZ4Ez8t00AMDZQ4wAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJDUp9QAwUVWu3pZ6BIBxwZURACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkvKtvQCMOd9qzW9zZQQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJIaVow0NTVFZWVllJSURFVVVezateuU6zds2BCf//zn49xzz42KiopYtWpV/PrXvx7WwADAxJJ3jGzdujXq6uqioaEhdu/eHXPmzIna2tp49913h1z/1FNPxerVq6OhoSHeeOONeOKJJ2Lr1q1xzz33fOLhAYDxL+8YeeSRR+K2226LFStWxBe+8IXYuHFjfOpTn4rNmzcPuf6VV16JK664Im666aaorKyMr3zlK3HjjTd+7NUUAODskFeM9Pb2Rmtra9TU1Hx0gMLCqKmpiZaWliH3ufzyy6O1tXUgPvbt2xfbt2+P66677qT309PTE11dXYNuAMDElNd30xw5ciT6+vqivLx80Pby8vLYs2fPkPvcdNNNceTIkfjSl74UWZbF8ePH44477jjl0zSNjY2xdu3afEYDAMapUX83zc6dO2PdunXx2GOPxe7du+O5556Lbdu2xQMPPHDSferr66Ozs3PgdujQodEeEwBIJK8rI9OmTYuioqLo6OgYtL2joyNmzJgx5D73339/LF26NG699daIiLj00kuju7s7br/99rj33nujsPDEHsrlcpHL5fIZDQAYp/K6MlJcXBzz58+P5ubmgW39/f3R3Nwc1dXVQ+7zwQcfnBAcRUVFERGRZVm+8wIAE0xeV0YiIurq6mL58uWxYMGCWLhwYWzYsCG6u7tjxYoVERGxbNmymD17djQ2NkZExOLFi+ORRx6JefPmRVVVVezduzfuv//+WLx48UCUAABnr7xjZMmSJXH48OFYs2ZNtLe3x9y5c2PHjh0DL2o9ePDgoCsh9913XxQUFMR9990Xb7/9dnzmM5+JxYsXx7e+9a2RexQAwLhVkI2D50q6urqirKwsOjs7o7S0dESPXbl624geDwDGmwPrF43KcU/377fvpgEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFLDipGmpqaorKyMkpKSqKqqil27dp1y/a9+9atYuXJlzJw5M3K5XFx44YWxffv2YQ0MAEwsk/LdYevWrVFXVxcbN26Mqqqq2LBhQ9TW1sabb74Z06dPP2F9b29vXHPNNTF9+vT4wQ9+ELNnz47//u//jilTpozE/ADAOJd3jDzyyCNx2223xYoVKyIiYuPGjbFt27bYvHlzrF69+oT1mzdvjvfeey9eeeWVOOeccyIiorKy8pNNDQBMGHk9TdPb2xutra1RU1Pz0QEKC6OmpiZaWlqG3OdHP/pRVFdXx8qVK6O8vDwuueSSWLduXfT19X2yyQGACSGvKyNHjhyJvr6+KC8vH7S9vLw89uzZM+Q++/bti5/+9Kdx8803x/bt22Pv3r1x5513xrFjx6KhoWHIfXp6eqKnp2fg566urnzGBADGkVF/N01/f39Mnz49Hn/88Zg/f34sWbIk7r333ti4ceNJ92lsbIyysrKBW0VFxWiPCQAkkleMTJs2LYqKiqKjo2PQ9o6OjpgxY8aQ+8ycOTMuvPDCKCoqGth28cUXR3t7e/T29g65T319fXR2dg7cDh06lM+YAMA4kleMFBcXx/z586O5uXlgW39/fzQ3N0d1dfWQ+1xxxRWxd+/e6O/vH9j21ltvxcyZM6O4uHjIfXK5XJSWlg66AQATU95P09TV1cWmTZvie9/7Xrzxxhvxta99Lbq7uwfeXbNs2bKor68fWP+1r30t3nvvvbjrrrvirbfeim3btsW6deti5cqVI/coAIBxK++39i5ZsiQOHz4ca9asifb29pg7d27s2LFj4EWtBw8ejMLCjxqnoqIiXnjhhVi1alVcdtllMXv27Ljrrrvi7rvvHrlHAQCMWwVZlmWph/g4XV1dUVZWFp2dnSP+lE3l6m0jejwAGG8OrF80Ksc93b/fvpsGAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAIKlhxUhTU1NUVlZGSUlJVFVVxa5du05rvy1btkRBQUFcf/31w7lbAGACyjtGtm7dGnV1ddHQ0BC7d++OOXPmRG1tbbz77run3O/AgQPxt3/7t3HllVcOe1gAYOLJO0YeeeSRuO2222LFihXxhS98ITZu3Bif+tSnYvPmzSfdp6+vL26++eZYu3ZtnH/++Z9oYABgYskrRnp7e6O1tTVqamo+OkBhYdTU1ERLS8tJ9/vmN78Z06dPj1tuueW07qenpye6uroG3QCAiSmvGDly5Ej09fVFeXn5oO3l5eXR3t4+5D4vv/xyPPHEE7Fp06bTvp/GxsYoKysbuFVUVOQzJgAwjozqu2mOHj0aS5cujU2bNsW0adNOe7/6+vro7OwcuB06dGgUpwQAUpqUz+Jp06ZFUVFRdHR0DNre0dERM2bMOGH9L37xizhw4EAsXrx4YFt/f/+HdzxpUrz55ptxwQUXnLBfLpeLXC6Xz2gAwDiV15WR4uLimD9/fjQ3Nw9s6+/vj+bm5qiurj5h/UUXXRSvv/56tLW1Ddy++tWvxtVXXx1tbW2efgEA8rsyEhFRV1cXy5cvjwULFsTChQtjw4YN0d3dHStWrIiIiGXLlsXs2bOjsbExSkpK4pJLLhm0/5QpUyIiTtgOAJyd8o6RJUuWxOHDh2PNmjXR3t4ec+fOjR07dgy8qPXgwYNRWOiDXQGA01OQZVmWeoiP09XVFWVlZdHZ2RmlpaUjeuzK1dtG9HgAMN4cWL9oVI57un+/XcIAAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSw4qRpqamqKysjJKSkqiqqopdu3addO2mTZviyiuvjKlTp8bUqVOjpqbmlOsBgLNL3jGydevWqKuri4aGhti9e3fMmTMnamtr49133x1y/c6dO+PGG2+Mf//3f4+WlpaoqKiIr3zlK/H2229/4uEBgPGvIMuyLJ8dqqqq4otf/GI8+uijERHR398fFRUV8fWvfz1Wr179sfv39fXF1KlT49FHH41ly5ad1n12dXVFWVlZdHZ2RmlpaT7jfqzK1dtG9HgAMN4cWL9oVI57un+/87oy0tvbG62trVFTU/PRAQoLo6amJlpaWk7rGB988EEcO3YszjvvvJOu6enpia6urkE3AGBiyitGjhw5En19fVFeXj5oe3l5ebS3t5/WMe6+++6YNWvWoKD5XY2NjVFWVjZwq6ioyGdMAGAcGdN306xfvz62bNkSP/zhD6OkpOSk6+rr66Ozs3PgdujQoTGcEgAYS5PyWTxt2rQoKiqKjo6OQds7OjpixowZp9z3oYceivXr18dPfvKTuOyyy065NpfLRS6Xy2c0AGCcyuvKSHFxccyfPz+am5sHtvX390dzc3NUV1efdL8HH3wwHnjggdixY0csWLBg+NMCABNOXldGIiLq6upi+fLlsWDBgli4cGFs2LAhuru7Y8WKFRERsWzZspg9e3Y0NjZGRMQ//MM/xJo1a+Kpp56KysrKgdeWfPrTn45Pf/rTI/hQAIDxKO8YWbJkSRw+fDjWrFkT7e3tMXfu3NixY8fAi1oPHjwYhYUfXXD5zne+E729vfEXf/EXg47T0NAQ3/jGNz7Z9ADAuJf354yk4HNGAGD0jKvPGQEAGGliBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhqWDHS1NQUlZWVUVJSElVVVbFr165Trn/22WfjoosuipKSkrj00ktj+/btwxoWAJh48o6RrVu3Rl1dXTQ0NMTu3btjzpw5UVtbG+++++6Q61955ZW48cYb45ZbbonXXnstrr/++rj++uvjP//zPz/x8ADA+FeQZVmWzw5VVVXxxS9+MR599NGIiOjv74+Kior4+te/HqtXrz5h/ZIlS6K7uzt+/OMfD2z70z/905g7d25s3LjxtO6zq6srysrKorOzM0pLS/MZ92NVrt42oscDgPHmwPpFo3Lc0/37PSmfg/b29kZra2vU19cPbCssLIyamppoaWkZcp+Wlpaoq6sbtK22tjaef/75k95PT09P9PT0DPzc2dkZER8+qJHW3/PBiB8TAMaT0fj7+tvH/bjrHnnFyJEjR6Kvry/Ky8sHbS8vL489e/YMuU97e/uQ69vb2096P42NjbF27doTtldUVOQzLgBwGso2jO7xjx49GmVlZSf9fV4xMlbq6+sHXU3p7++P9957L37v934vCgoKRux+urq6oqKiIg4dOjTiT/+QP+fjzOOcnFmcjzOL8/HxsiyLo0ePxqxZs065Lq8YmTZtWhQVFUVHR8eg7R0dHTFjxowh95kxY0Ze6yMicrlc5HK5QdumTJmSz6h5KS0t9X+kM4jzceZxTs4szseZxfk4tVNdEfmNvN5NU1xcHPPnz4/m5uaBbf39/dHc3BzV1dVD7lNdXT1ofUTEiy++eNL1AMDZJe+naerq6mL58uWxYMGCWLhwYWzYsCG6u7tjxYoVERGxbNmymD17djQ2NkZExF133RVf/vKX4+GHH45FixbFli1b4tVXX43HH398ZB8JADAu5R0jS5YsicOHD8eaNWuivb095s6dGzt27Bh4kerBgwejsPCjCy6XX355PPXUU3HffffFPffcE3/0R38Uzz//fFxyySUj9yiGKZfLRUNDwwlPCZGG83HmcU7OLM7HmcX5GDl5f84IAMBI8t00AEBSYgQASEqMAABJiREAIKkJHyNNTU1RWVkZJSUlUVVVFbt27Trl+meffTYuuuiiKCkpiUsvvTS2b98+RpOeHfI5H5s2bYorr7wypk6dGlOnTo2ampqPPX/kL99/I7+xZcuWKCgoiOuvv350BzzL5Hs+fvWrX8XKlStj5syZkcvl4sILL/TfrRGU7/nYsGFDfP7zn49zzz03KioqYtWqVfHrX/96jKYdx7IJbMuWLVlxcXG2efPm7L/+67+y2267LZsyZUrW0dEx5Pqf/exnWVFRUfbggw9mP//5z7P77rsvO+ecc7LXX399jCefmPI9HzfddFPW1NSUvfbaa9kbb7yR/eVf/mVWVlaW/c///M8YTz5x5XtOfmP//v3Z7NmzsyuvvDL78z//87EZ9iyQ7/no6enJFixYkF133XXZyy+/nO3fvz/buXNn1tbWNsaTT0z5no/vf//7WS6Xy77//e9n+/fvz1544YVs5syZ2apVq8Z48vFnQsfIwoULs5UrVw783NfXl82aNStrbGwccv0NN9yQLVq0aNC2qqqq7K/+6q9Gdc6zRb7n43cdP348mzx5cva9731vtEY86wznnBw/fjy7/PLLs+9+97vZ8uXLxcgIyvd8fOc738nOP//8rLe3d6xGPKvkez5WrlyZ/dmf/dmgbXV1ddkVV1wxqnNOBBP2aZre3t5obW2NmpqagW2FhYVRU1MTLS0tQ+7T0tIyaH1ERG1t7UnXc/qGcz5+1wcffBDHjh2L8847b7TGPKsM95x885vfjOnTp8ctt9wyFmOeNYZzPn70ox9FdXV1rFy5MsrLy+OSSy6JdevWRV9f31iNPWEN53xcfvnl0draOvBUzr59+2L79u1x3XXXjcnM49kZ+a29I+HIkSPR19c38Mmwv1FeXh579uwZcp/29vYh17e3t4/anGeL4ZyP33X33XfHrFmzTghGhmc45+Tll1+OJ554Itra2sZgwrPLcM7Hvn374qc//WncfPPNsX379ti7d2/ceeedcezYsWhoaBiLsSes4ZyPm266KY4cORJf+tKXIsuyOH78eNxxxx1xzz33jMXI49qEvTLCxLJ+/frYsmVL/PCHP4ySkpLU45yVjh49GkuXLo1NmzbFtGnTUo9DfPhFpdOnT4/HH3885s+fH0uWLIl77703Nm7cmHq0s9LOnTtj3bp18dhjj8Xu3bvjueeei23btsUDDzyQerQz3oS9MjJt2rQoKiqKjo6OQds7OjpixowZQ+4zY8aMvNZz+oZzPn7joYceivXr18dPfvKTuOyyy0ZzzLNKvufkF7/4RRw4cCAWL148sK2/vz8iIiZNmhRvvvlmXHDBBaM79AQ2nH8jM2fOjHPOOSeKiooGtl188cXR3t4evb29UVxcPKozT2TDOR/3339/LF26NG699daIiLj00kuju7s7br/99rj33nsHfW8bg03Y/2WKi4tj/vz50dzcPLCtv78/mpubo7q6esh9qqurB62PiHjxxRdPup7TN5zzERHx4IMPxgMPPBA7duyIBQsWjMWoZ418z8lFF10Ur7/+erS1tQ3cvvrVr8bVV18dbW1tUVFRMZbjTzjD+TdyxRVXxN69eweiMCLirbfeipkzZwqRT2g45+ODDz44ITh+E4qZr4E7tdSvoB1NW7ZsyXK5XPbkk09mP//5z7Pbb789mzJlStbe3p5lWZYtXbo0W7169cD6n/3sZ9mkSZOyhx56KHvjjTeyhoYGb+0dQfmej/Xr12fFxcXZD37wg+ydd94ZuB09ejTVQ5hw8j0nv8u7aUZWvufj4MGD2eTJk7O//uu/zt58883sxz/+cTZ9+vTs7//+71M9hAkl3/PR0NCQTZ48OXv66aezffv2Zf/2b/+WXXDBBdkNN9yQ6iGMGxM6RrIsy/7pn/4p++xnP5sVFxdnCxcuzP7jP/5j4Hdf/vKXs+XLlw9a/8wzz2QXXnhhVlxcnP3xH/9xtm3btjGeeGLL53z8wR/8QRYRJ9waGhrGfvAJLN9/I79NjIy8fM/HK6+8klVVVWW5XC47//zzs29961vZ8ePHx3jqiSuf83Hs2LHsG9/4RnbBBRdkJSUlWUVFRXbnnXdm//u//zv2g48zBVnm2hEAkM6Efc0IADA+iBEAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICk/h88mukmS30MIgAAAABJRU5ErkJggg==\n"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def simulation_rank(prior,likelihood):\n",
    "    \"\"\"\n",
    "    1. Draw (θٰʹ,yʹ) ~ p(θ,y)\n",
    "    2. Use MCMC to generate samples approximately from p(θ|yʹ)\n",
    "    3. Find rank of θʹ among samples from θ\n",
    "    \"\"\"\n",
    "    latent_vars = prior()\n",
    "    observed_vars = likelihood(*latent_vars)\n",
    "\n",
    "    # draw latents and observed variables from prior\n",
    "    all_samps = sample(latent_vars + observed_vars, niter=1) # just draw one sample\n",
    "    all_samps = [samps[0] for samps in all_samps] # remove dummy first dimension\n",
    "    latent_samps = all_samps[:len(latent_vars)]\n",
    "    observed_samps = all_samps[len(latent_vars):]\n",
    "\n",
    "    # do inference on latents given observed variables\n",
    "    evidence = [observed_var == observed_samp for (observed_var, observed_samp) in zip(observed_vars, observed_samps)]\n",
    "    inference_samps = sample(latent_vars, Given(*evidence), niter=100)\n",
    "\n",
    "    # get rank of true latents compared to sampled latents\n",
    "    ranks = []\n",
    "    for latent_samp, inference_samp in zip(latent_samps, inference_samps):\n",
    "        ranks.append( np.mean(latent_samp < inference_samp, axis=0))\n",
    "\n",
    "    return ranks\n",
    "\n",
    "def SBC(prior,likelihood,nreps=1000):\n",
    "    \"\"\"\n",
    "    Call simulation_rank several times, make a histogram\n",
    "    \"\"\"\n",
    "    my_ranks = simulation_rank(prior,likelihood)\n",
    "    ranks = [np.zeros([nreps, *r.shape]) for r in my_ranks]\n",
    "    for i in range(nreps):\n",
    "        my_ranks = simulation_rank(prior,likelihood)\n",
    "        for r, my_r in zip(ranks, my_ranks):\n",
    "            r[i] = my_r\n",
    "\n",
    "    for r in ranks:\n",
    "        plt.figure()\n",
    "        plt.hist(r,np.arange(0,1,.1),density=True)\n",
    "        # annoying we can't do plt.title(nameof(r))...\n",
    "\n",
    "# Simple model with two latent variables and two observation variables, all scalar\n",
    "def prior():\n",
    "    u = d.norm(0,1)\n",
    "    s = d.exp(1)\n",
    "    return u,s\n",
    "\n",
    "def likelihood(u,s):\n",
    "    y = d.norm(u,s)\n",
    "    z = d.norm(u,s)\n",
    "    return y,z\n",
    "\n",
    "SBC(prior,likelihood)"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "outputs": [
    {
     "data": {
      "text/plain": "<Figure size 640x480 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGdCAYAAADAAnMpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAeLklEQVR4nO3df3DX9X3A8VcI5ht7JQFHCT+WLtXNaqcCg5JF61V3sZxydP6xk1MPGOePWVnPkdtN4g9S60qYp47djOWkcnZ3VdCedV3D4Ww6zrNmxwnmzq2iR4HBrIlwrgmGNYHksz88vy4lIN+Y5E3C43H3/SOfvD+f7+vbTzHP+3x/FWVZlgUAQCITUg8AAJzdxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQ1MfUAp6O/vz9+9atfxaRJk6KoqCj1OADAaciyLI4cORIzZ86MCRNOfv1jTMTIr371q6isrEw9BgAwBAcPHozf/d3fPenvx0SMTJo0KSI+fDBlZWWJpwEATkdXV1dUVlbm/46fzJiIkY+emikrKxMjADDGfNJLLLyAFQBISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQ1MTUAwDAWFW1ujn1CMNi/7pFSe/flREAICkxAgAkJUYAgKTECACQlBgBAJLybhrgrOBdD3DmcmUEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkVXCMvPzyy7F48eKYOXNmFBUVxQsvvHDK9c8//3xcc8018bnPfS7KysqipqYmXnzxxaHOCwCMMwXHSHd3d8yePTuamppOa/3LL78c11xzTWzdujV27twZV199dSxevDhef/31gocFAMafiYXucO2118a111572uvXr18/4Oe1a9fGP//zP8e//Mu/xNy5cwu9ewBgnCk4Rj6t/v7+OHLkSJx33nknXdPT0xM9PT35n7u6ukZjtDGtanVz6hGGxf51i1KPAMAoG/UXsD788MPxwQcfxA033HDSNY2NjVFeXp6/VVZWjuKEAMBoGtUYefrpp+OBBx6IZ599NqZNm3bSdfX19dHZ2Zm/HTx4cBSnBABG06g9TbN58+a49dZb47nnnova2tpTrs3lcpHL5UZpMhgZnjoDOD2jcmXkmWeeiRUrVsQzzzwTixb5DxsA8LGCr4x88MEHsWfPnvzP+/bti7a2tjjvvPPi85//fNTX18c777wT//RP/xQRHz41s3z58viHf/iHqK6ujvb29oiIOPfcc6O8vHyYHgYAMFYVfGXktddei7lz5+bflltXVxdz586NNWvWRETEu+++GwcOHMivf+KJJ+L48eOxcuXKmDFjRv521113DdNDAADGsoKvjFx11VWRZdlJf//UU08N+Hn79u2F3gUAcBbx3TQAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJDUqH03DQB8ZLx8dxPDw5URACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJDUx9QDw/1Wtbk49AgCjzJURACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABIquAYefnll2Px4sUxc+bMKCoqihdeeOET99m+fXv80R/9UeRyufj93//9eOqpp4YwKgAwHhUcI93d3TF79uxoamo6rfX79u2LRYsWxdVXXx1tbW3xV3/1V3HrrbfGiy++WPCwAMD4U/AnsF577bVx7bXXnvb6DRs2xBe+8IV45JFHIiLi4osvjldeeSX+/u//PhYuXFjo3QOjzKfiAiNtxF8z0traGrW1tQO2LVy4MFpbW0f6rgGAMWDEv5umvb09KioqBmyrqKiIrq6u+N///d8499xzT9inp6cnenp68j93dXWN9JgAQCJn5LtpGhsbo7y8PH+rrKxMPRIAMEJGPEamT58eHR0dA7Z1dHREWVnZoFdFIiLq6+ujs7Mzfzt48OBIjwkAJDLiT9PU1NTE1q1bB2x76aWXoqam5qT75HK5yOVyIz0aAHAGKPjKyAcffBBtbW3R1tYWER++dbetrS0OHDgQER9e1Vi2bFl+/R133BF79+6Nv/mbv4ndu3fH448/Hs8++2ysWrVqeB4BADCmFRwjr732WsydOzfmzp0bERF1dXUxd+7cWLNmTUREvPvuu/kwiYj4whe+EM3NzfHSSy/F7Nmz45FHHonvfe973tYLAETEEJ6mueqqqyLLspP+frBPV73qqqvi9ddfL/SuAICzwBn5bhoA4OwhRgCApMQIAJCUGAEAkhIjAEBSYgQASGrEP4EVgOFTtbo59Qgw7FwZAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgqYmpB0itanVz6hEA4KzmyggAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgqSHFSFNTU1RVVUVpaWlUV1fHjh07Trl+/fr18cUvfjHOPffcqKysjFWrVsVvfvObIQ0MAIwvBcfIli1boq6uLhoaGmLXrl0xe/bsWLhwYbz33nuDrn/66adj9erV0dDQEG+++WY8+eSTsWXLlrjnnns+9fAAwNhXcIw8+uijcdttt8WKFSviS1/6UmzYsCE+85nPxKZNmwZd/+qrr8YVV1wRN910U1RVVcXXvva1uPHGGz/xagoAcHYoKEZ6e3tj586dUVtb+/EBJkyI2traaG1tHXSfyy+/PHbu3JmPj71798bWrVvjuuuuO+n99PT0RFdX14AbADA+FfRx8IcPH46+vr6oqKgYsL2ioiJ279496D433XRTHD58OL7yla9ElmVx/PjxuOOOO075NE1jY2M88MADhYwGAIxRI/5umu3bt8fatWvj8ccfj127dsXzzz8fzc3N8eCDD550n/r6+ujs7MzfDh48ONJjAgCJFHRlZOrUqVFcXBwdHR0Dtnd0dMT06dMH3ef++++PpUuXxq233hoREZdeeml0d3fH7bffHvfee29MmHBiD+VyucjlcoWMBgCMUQVdGSkpKYl58+ZFS0tLflt/f3+0tLRETU3NoPscPXr0hOAoLi6OiIgsywqdFwAYZwq6MhIRUVdXF8uXL4/58+fHggULYv369dHd3R0rVqyIiIhly5bFrFmzorGxMSIiFi9eHI8++mjMnTs3qqurY8+ePXH//ffH4sWL81ECAJy9Co6RJUuWxKFDh2LNmjXR3t4ec+bMiW3btuVf1HrgwIEBV0Luu+++KCoqivvuuy/eeeed+NznPheLFy+O73znO8P3KACAMasoGwPPlXR1dUV5eXl0dnZGWVnZsB67anXzsB4PAMaa/esWjchxT/fvt++mAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUkOKkaampqiqqorS0tKorq6OHTt2nHL9r3/961i5cmXMmDEjcrlcXHjhhbF169YhDQwAjC8TC91hy5YtUVdXFxs2bIjq6upYv359LFy4MN56662YNm3aCet7e3vjmmuuiWnTpsUPf/jDmDVrVvzXf/1XTJ48eTjmBwDGuIJj5NFHH43bbrstVqxYERERGzZsiObm5ti0aVOsXr36hPWbNm2K999/P1599dU455xzIiKiqqrq000NAIwbBT1N09vbGzt37oza2tqPDzBhQtTW1kZra+ug+/z4xz+OmpqaWLlyZVRUVMQll1wSa9eujb6+vk83OQAwLhR0ZeTw4cPR19cXFRUVA7ZXVFTE7t27B91n79698bOf/Sxuvvnm2Lp1a+zZsyfuvPPOOHbsWDQ0NAy6T09PT/T09OR/7urqKmRMAGAMGfF30/T398e0adPiiSeeiHnz5sWSJUvi3nvvjQ0bNpx0n8bGxigvL8/fKisrR3pMACCRgmJk6tSpUVxcHB0dHQO2d3R0xPTp0wfdZ8aMGXHhhRdGcXFxftvFF18c7e3t0dvbO+g+9fX10dnZmb8dPHiwkDEBgDGkoBgpKSmJefPmRUtLS35bf39/tLS0RE1NzaD7XHHFFbFnz57o7+/Pb3v77bdjxowZUVJSMug+uVwuysrKBtwAgPGp4Kdp6urqYuPGjfH9738/3nzzzfjGN74R3d3d+XfXLFu2LOrr6/Prv/GNb8T7778fd911V7z99tvR3Nwca9eujZUrVw7fowAAxqyC39q7ZMmSOHToUKxZsyba29tjzpw5sW3btvyLWg8cOBATJnzcOJWVlfHiiy/GqlWr4rLLLotZs2bFXXfdFXfffffwPQoAYMwqyrIsSz3EJ+nq6ory8vLo7Owc9qdsqlY3D+vxAGCs2b9u0Ygc93T/fvtuGgAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkhhQjTU1NUVVVFaWlpVFdXR07duw4rf02b94cRUVFcf311w/lbgGAcajgGNmyZUvU1dVFQ0ND7Nq1K2bPnh0LFy6M995775T77d+/P/76r/86rrzyyiEPCwCMPwXHyKOPPhq33XZbrFixIr70pS/Fhg0b4jOf+Uxs2rTppPv09fXFzTffHA888ECcf/75n2pgAGB8KShGent7Y+fOnVFbW/vxASZMiNra2mhtbT3pft/+9rdj2rRpccstt5zW/fT09ERXV9eAGwAwPhUUI4cPH46+vr6oqKgYsL2ioiLa29sH3eeVV16JJ598MjZu3Hja99PY2Bjl5eX5W2VlZSFjAgBjyIi+m+bIkSOxdOnS2LhxY0ydOvW096uvr4/Ozs787eDBgyM4JQCQ0sRCFk+dOjWKi4ujo6NjwPaOjo6YPn36Cet/+ctfxv79+2Px4sX5bf39/R/e8cSJ8dZbb8UFF1xwwn65XC5yuVwhowEAY1RBV0ZKSkpi3rx50dLSkt/W398fLS0tUVNTc8L6iy66KN54441oa2vL377+9a/H1VdfHW1tbZ5+AQAKuzISEVFXVxfLly+P+fPnx4IFC2L9+vXR3d0dK1asiIiIZcuWxaxZs6KxsTFKS0vjkksuGbD/5MmTIyJO2A4AnJ0KjpElS5bEoUOHYs2aNdHe3h5z5syJbdu25V/UeuDAgZgwwQe7AgCnpyjLsiz1EJ+kq6srysvLo7OzM8rKyob12FWrm4f1eAAw1uxft2hEjnu6f79dwgAAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABIakgx0tTUFFVVVVFaWhrV1dWxY8eOk67duHFjXHnllTFlypSYMmVK1NbWnnI9AHB2KThGtmzZEnV1ddHQ0BC7du2K2bNnx8KFC+O9994bdP327dvjxhtvjH/7t3+L1tbWqKysjK997WvxzjvvfOrhAYCxryjLsqyQHaqrq+PLX/5yPPbYYxER0d/fH5WVlfHNb34zVq9e/Yn79/X1xZQpU+Kxxx6LZcuWndZ9dnV1RXl5eXR2dkZZWVkh436iqtXNw3o8ABhr9q9bNCLHPd2/3wVdGent7Y2dO3dGbW3txweYMCFqa2ujtbX1tI5x9OjROHbsWJx33nknXdPT0xNdXV0DbgDA+FRQjBw+fDj6+vqioqJiwPaKiopob28/rWPcfffdMXPmzAFB89saGxujvLw8f6usrCxkTABgDBnVd9OsW7cuNm/eHD/60Y+itLT0pOvq6+ujs7Mzfzt48OAoTgkAjKaJhSyeOnVqFBcXR0dHx4DtHR0dMX369FPu+/DDD8e6devipz/9aVx22WWnXJvL5SKXyxUyGgAwRhV0ZaSkpCTmzZsXLS0t+W39/f3R0tISNTU1J93voYceigcffDC2bdsW8+fPH/q0AMC4U9CVkYiIurq6WL58ecyfPz8WLFgQ69evj+7u7lixYkVERCxbtixmzZoVjY2NERHxd3/3d7FmzZp4+umno6qqKv/aks9+9rPx2c9+dhgfCgAwFhUcI0uWLIlDhw7FmjVror29PebMmRPbtm3Lv6j1wIEDMWHCxxdcvvvd70Zvb2/82Z/92YDjNDQ0xLe+9a1PNz0AMOYV/DkjKficEQAYOWPqc0YAAIabGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSGlKMNDU1RVVVVZSWlkZ1dXXs2LHjlOufe+65uOiii6K0tDQuvfTS2Lp165CGBQDGn4JjZMuWLVFXVxcNDQ2xa9eumD17dixcuDDee++9Qde/+uqrceONN8Ytt9wSr7/+elx//fVx/fXXx3/8x3986uEBgLGvKMuyrJAdqqur48tf/nI89thjERHR398flZWV8c1vfjNWr159wvolS5ZEd3d3/OQnP8lv++M//uOYM2dObNiw4bTus6urK8rLy6OzszPKysoKGfcTVa1uHtbjAcBYs3/dohE57un+/Z5YyEF7e3tj586dUV9fn982YcKEqK2tjdbW1kH3aW1tjbq6ugHbFi5cGC+88MJJ76enpyd6enryP3d2dkbEhw9quPX3HB32YwLAWDISf1///3E/6bpHQTFy+PDh6Ovri4qKigHbKyoqYvfu3YPu097ePuj69vb2k95PY2NjPPDAAydsr6ysLGRcAOA0lK8f2eMfOXIkysvLT/r7gmJktNTX1w+4mtLf3x/vv/9+/M7v/E4UFRUN2/10dXVFZWVlHDx4cNif/qFwzseZxzk5szgfZxbn45NlWRZHjhyJmTNnnnJdQTEyderUKC4ujo6OjgHbOzo6Yvr06YPuM3369ILWR0TkcrnI5XIDtk2ePLmQUQtSVlbm/0hnEOfjzOOcnFmcjzOL83Fqp7oi8pGC3k1TUlIS8+bNi5aWlvy2/v7+aGlpiZqamkH3qampGbA+IuKll1466XoA4OxS8NM0dXV1sXz58pg/f34sWLAg1q9fH93d3bFixYqIiFi2bFnMmjUrGhsbIyLirrvuiq9+9avxyCOPxKJFi2Lz5s3x2muvxRNPPDG8jwQAGJMKjpElS5bEoUOHYs2aNdHe3h5z5syJbdu25V+keuDAgZgw4eMLLpdffnk8/fTTcd9998U999wTf/AHfxAvvPBCXHLJJcP3KIYol8tFQ0PDCU8JkYbzceZxTs4szseZxfkYPgV/zggAwHDy3TQAQFJiBABISowAAEmJEQAgqXEfI01NTVFVVRWlpaVRXV0dO3bsOOX65557Li666KIoLS2NSy+9NLZu3TpKk54dCjkfGzdujCuvvDKmTJkSU6ZMidra2k88fxSu0H8jH9m8eXMUFRXF9ddfP7IDnmUKPR+//vWvY+XKlTFjxozI5XJx4YUX+u/WMCr0fKxfvz6++MUvxrnnnhuVlZWxatWq+M1vfjNK045h2Ti2efPmrKSkJNu0aVP2n//5n9ltt92WTZ48Oevo6Bh0/c9//vOsuLg4e+ihh7Jf/OIX2X333Zedc8452RtvvDHKk49PhZ6Pm266KWtqaspef/317M0338z+/M//PCsvL8/++7//e5QnH78KPScf2bdvXzZr1qzsyiuvzP70T/90dIY9CxR6Pnp6erL58+dn1113XfbKK69k+/bty7Zv3561tbWN8uTjU6Hn4wc/+EGWy+WyH/zgB9m+ffuyF198MZsxY0a2atWqUZ587BnXMbJgwYJs5cqV+Z/7+vqymTNnZo2NjYOuv+GGG7JFixYN2FZdXZ39xV/8xYjOebYo9Hz8tuPHj2eTJk3Kvv/974/UiGedoZyT48ePZ5dffnn2ve99L1u+fLkYGUaFno/vfve72fnnn5/19vaO1ohnlULPx8qVK7M/+ZM/GbCtrq4uu+KKK0Z0zvFg3D5N09vbGzt37oza2tr8tgkTJkRtbW20trYOuk9ra+uA9RERCxcuPOl6Tt9QzsdvO3r0aBw7dizOO++8kRrzrDLUc/Ltb387pk2bFrfccstojHnWGMr5+PGPfxw1NTWxcuXKqKioiEsuuSTWrl0bfX19ozX2uDWU83H55ZfHzp0780/l7N27N7Zu3RrXXXfdqMw8lp2R39o7HA4fPhx9fX35T4b9SEVFRezevXvQfdrb2wdd397ePmJzni2Gcj5+29133x0zZ848IRgZmqGck1deeSWefPLJaGtrG4UJzy5DOR979+6Nn/3sZ3HzzTfH1q1bY8+ePXHnnXfGsWPHoqGhYTTGHreGcj5uuummOHz4cHzlK1+JLMvi+PHjcccdd8Q999wzGiOPaeP2ygjjy7p162Lz5s3xox/9KEpLS1OPc1Y6cuRILF26NDZu3BhTp05NPQ7x4ReVTps2LZ544omYN29eLFmyJO69997YsGFD6tHOStu3b4+1a9fG448/Hrt27Yrnn38+mpub48EHH0w92hlv3F4ZmTp1ahQXF0dHR8eA7R0dHTF9+vRB95k+fXpB6zl9QzkfH3n44Ydj3bp18dOf/jQuu+yykRzzrFLoOfnlL38Z+/fvj8WLF+e39ff3R0TExIkT46233ooLLrhgZIcex4byb2TGjBlxzjnnRHFxcX7bxRdfHO3t7dHb2xslJSUjOvN4NpTzcf/998fSpUvj1ltvjYiISy+9NLq7u+P222+Pe++9d8D3tjHQuP1fpqSkJObNmxctLS35bf39/dHS0hI1NTWD7lNTUzNgfUTESy+9dNL1nL6hnI+IiIceeigefPDB2LZtW8yfP380Rj1rFHpOLrroonjjjTeira0tf/v6178eV199dbS1tUVlZeVojj/uDOXfyBVXXBF79uzJR2FExNtvvx0zZswQIp/SUM7H0aNHTwiOj0Ix8zVwp5b6FbQjafPmzVkul8ueeuqp7Be/+EV2++23Z5MnT87a29uzLMuypUuXZqtXr86v//nPf55NnDgxe/jhh7M333wza2ho8NbeYVTo+Vi3bl1WUlKS/fCHP8zefffd/O3IkSOpHsK4U+g5+W3eTTO8Cj0fBw4cyCZNmpT95V/+ZfbWW29lP/nJT7Jp06Zlf/u3f5vqIYwrhZ6PhoaGbNKkSdkzzzyT7d27N/vXf/3X7IILLshuuOGGVA9hzBjXMZJlWfaP//iP2ec///mspKQkW7BgQfbv//7v+d999atfzZYvXz5g/bPPPptdeOGFWUlJSfaHf/iHWXNz8yhPPL4Vcj5+7/d+L4uIE24NDQ2jP/g4Vui/kf9PjAy/Qs/Hq6++mlVXV2e5XC47//zzs+985zvZ8ePHR3nq8auQ83Hs2LHsW9/6VnbBBRdkpaWlWWVlZXbnnXdm//M//zP6g48xRVnm2hEAkM64fc0IADA2iBEAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICk/g+coeLskBlQnwAAAABJRU5ErkJggg==\n"
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "<Figure size 640x480 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGdCAYAAADAAnMpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAePElEQVR4nO3dfWyd1X3A8Z/j4Guqxk6yNM7LzFzYKHRAkiWNZyjqmEwtiNLxx0QEVZJlvIw2qxjWNGJe4lLaOEPAMg23ESkRrVRIAFHW1VEYdRchiqcIh0hs5UVpkiWj2CRitYNZbWI/+6PiMjdOyDW2T2x/PtL9w0/Oee65nAR/9Tz32kVZlmUBAJDIlNQLAAAmNzECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJTU29gNMxMDAQv/zlL2PatGlRVFSUejkAwGnIsiyOHTsW8+bNiylTTn79Y1zEyC9/+cuorKxMvQwAYBgOHz4cv/u7v3vSPx8XMTJt2rSI+M2LKSsrS7waAOB0dHd3R2VlZf77+MmMixj54NZMWVmZGAGAceaj3mLhDawAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKSmpl4AwFioWteSegkj4uDGZamXACPOlREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABIquAYef7552P58uUxb968KCoqimeeeeaU459++um48sor41Of+lSUlZVFTU1NPPvss8NdLwAwwRQcIz09PbFgwYJobm4+rfHPP/98XHnllbFjx45ob2+PK664IpYvXx4vv/xywYsFACaeqYVOuOqqq+Kqq6467fGbNm0a9PWGDRvin//5n+Nf/uVfYtGiRYU+PQAwwRQcIx/XwMBAHDt2LGbOnHnSMb29vdHb25v/uru7eyyWBgAkMOZvYL3//vvj3XffjWuvvfakY5qamqK8vDz/qKysHMMVAgBjaUxj5LHHHot77rknnnjiiZg9e/ZJxzU0NERXV1f+cfjw4TFcJQAwlsbsNs22bdvixhtvjCeffDJqa2tPOTaXy0UulxujlQEAKY3JlZHHH3881qxZE48//ngsW7ZsLJ4SABgnCr4y8u6778a+ffvyXx84cCD27t0bM2fOjHPOOScaGhrizTffjO9///sR8ZtbM6tXr45//Md/jOrq6ujo6IiIiLPPPjvKy8tH6GUAAONVwVdGXnrppVi0aFH+Y7n19fWxaNGiWL9+fUREvPXWW3Ho0KH8+IcffjiOHz8ea9eujblz5+Yft9566wi9BABgPCv4ysif/MmfRJZlJ/3zRx99dNDXu3btKvQpAIBJxO+mAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSBf84eGByqVrXknoJwATnyggAkJQYAQCScpuGM8pEuiVwcOOy1EsAGBdcGQEAkhIjAEBSbtMAwDBNlFvLqW8ruzICACQlRgCApCb9bRqX2AAgLVdGAICkxAgAkNSkv00Do2Wi3ALkzDJR/l65tcz/58oIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKT80LMJYqL8ICQAJh9XRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUgXHyPPPPx/Lly+PefPmRVFRUTzzzDMfOWfXrl3xR3/0R5HL5eL3f//349FHHx3GUgGAiajgGOnp6YkFCxZEc3PzaY0/cOBALFu2LK644orYu3dv/M3f/E3ceOON8eyzzxa8WABg4pla6ISrrroqrrrqqtMev3nz5vj0pz8dDzzwQEREXHjhhfHCCy/EP/zDP0RdXV2hTw8ATDCj/p6Rtra2qK2tHXSsrq4u2traTjqnt7c3uru7Bz0AgIlp1GOko6MjKioqBh2rqKiI7u7u+N///d8h5zQ1NUV5eXn+UVlZOdrLBAASOSM/TdPQ0BBdXV35x+HDh1MvCQAYJQW/Z6RQc+bMic7OzkHHOjs7o6ysLM4+++wh5+RyucjlcqO9NADgDDDqV0ZqamqitbV10LHnnnsuampqRvupAYBxoOArI++++27s27cv//WBAwdi7969MXPmzDjnnHOioaEh3nzzzfj+978fERG33HJLPPTQQ/F3f/d38Zd/+Zfx05/+NJ544oloaWkZuVcBwLhStc73AD5U8JWRl156KRYtWhSLFi2KiIj6+vpYtGhRrF+/PiIi3nrrrTh06FB+/Kc//eloaWmJ5557LhYsWBAPPPBAfPe73/WxXgAgIiKKsizLUi/io3R3d0d5eXl0dXVFWVnZiJ5bnQMw2R3cuGxUznu637/PyE/TAACThxgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASGpYMdLc3BxVVVVRWloa1dXVsXv37lOO37RpU3zmM5+Js88+OyorK+O2226LX//618NaMAAwsRQcI9u3b4/6+vpobGyMPXv2xIIFC6Kuri7efvvtIcc/9thjsW7dumhsbIxXX301Hnnkkdi+fXvccccdH3vxAMD4V3CMPPjgg3HTTTfFmjVr4rOf/Wxs3rw5PvGJT8TWrVuHHP/iiy/GZZddFtdff31UVVXFF7/4xbjuuus+8moKADA5FBQjfX190d7eHrW1tR+eYMqUqK2tjba2tiHnXHrppdHe3p6Pj/3798eOHTvi6quv/hjLBgAmiqmFDD569Gj09/dHRUXFoOMVFRXx2muvDTnn+uuvj6NHj8bnP//5yLIsjh8/Hrfccsspb9P09vZGb29v/uvu7u5ClgkAjCOj/mmaXbt2xYYNG+Lb3/527NmzJ55++uloaWmJe++996Rzmpqaory8PP+orKwc7WUCAIkUdGVk1qxZUVxcHJ2dnYOOd3Z2xpw5c4acc/fdd8fKlSvjxhtvjIiIiy++OHp6euLmm2+OO++8M6ZMObGHGhoaor6+Pv91d3e3IAGACaqgKyMlJSWxePHiaG1tzR8bGBiI1tbWqKmpGXLOe++9d0JwFBcXR0RElmVDzsnlclFWVjboAQBMTAVdGYmIqK+vj9WrV8eSJUti6dKlsWnTpujp6Yk1a9ZERMSqVati/vz50dTUFBERy5cvjwcffDAWLVoU1dXVsW/fvrj77rtj+fLl+SgBACavgmNkxYoVceTIkVi/fn10dHTEwoULY+fOnfk3tR46dGjQlZC77rorioqK4q677oo333wzPvWpT8Xy5cvjW9/61si9CgBg3CrKTnav5AzS3d0d5eXl0dXVNeK3bKrWtYzo+QBgvDm4cdmonPd0v3/73TQAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJDWsGGlubo6qqqooLS2N6urq2L179ynH/+pXv4q1a9fG3LlzI5fLxfnnnx87duwY1oIBgIllaqETtm/fHvX19bF58+aorq6OTZs2RV1dXbz++usxe/bsE8b39fXFlVdeGbNnz46nnnoq5s+fH//1X/8V06dPH4n1AwDjXMEx8uCDD8ZNN90Ua9asiYiIzZs3R0tLS2zdujXWrVt3wvitW7fGO++8Ey+++GKcddZZERFRVVX18VYNAEwYBd2m6evri/b29qitrf3wBFOmRG1tbbS1tQ0550c/+lHU1NTE2rVro6KiIi666KLYsGFD9Pf3n/R5ent7o7u7e9ADAJiYCoqRo0ePRn9/f1RUVAw6XlFRER0dHUPO2b9/fzz11FPR398fO3bsiLvvvjseeOCB+OY3v3nS52lqaory8vL8o7KyspBlAgDjyKh/mmZgYCBmz54dDz/8cCxevDhWrFgRd955Z2zevPmkcxoaGqKrqyv/OHz48GgvEwBIpKD3jMyaNSuKi4ujs7Nz0PHOzs6YM2fOkHPmzp0bZ511VhQXF+ePXXjhhdHR0RF9fX1RUlJywpxcLhe5XK6QpQEA41RBV0ZKSkpi8eLF0dramj82MDAQra2tUVNTM+Scyy67LPbt2xcDAwP5Y2+88UbMnTt3yBABACaXgm/T1NfXx5YtW+J73/tevPrqq/GVr3wlenp68p+uWbVqVTQ0NOTHf+UrX4l33nknbr311njjjTeipaUlNmzYEGvXrh25VwEAjFsFf7R3xYoVceTIkVi/fn10dHTEwoULY+fOnfk3tR46dCimTPmwcSorK+PZZ5+N2267LS655JKYP39+3HrrrXH77beP3KsAAMatoizLstSL+Cjd3d1RXl4eXV1dUVZWNqLnrlrXMqLnA4Dx5uDGZaNy3tP9/u130wAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKSGFSPNzc1RVVUVpaWlUV1dHbt37z6tedu2bYuioqK45pprhvO0AMAEVHCMbN++Perr66OxsTH27NkTCxYsiLq6unj77bdPOe/gwYPxt3/7t3H55ZcPe7EAwMRTcIw8+OCDcdNNN8WaNWvis5/9bGzevDk+8YlPxNatW086p7+/P7785S/HPffcE+eee+7HWjAAMLEUFCN9fX3R3t4etbW1H55gypSora2Ntra2k877xje+EbNnz44bbrjhtJ6nt7c3uru7Bz0AgImpoBg5evRo9Pf3R0VFxaDjFRUV0dHRMeScF154IR555JHYsmXLaT9PU1NTlJeX5x+VlZWFLBMAGEdG9dM0x44di5UrV8aWLVti1qxZpz2voaEhurq68o/Dhw+P4ioBgJSmFjJ41qxZUVxcHJ2dnYOOd3Z2xpw5c04Y/4tf/CIOHjwYy5cvzx8bGBj4zRNPnRqvv/56nHfeeSfMy+VykcvlClkaADBOFXRlpKSkJBYvXhytra35YwMDA9Ha2ho1NTUnjL/gggvilVdeib179+YfX/rSl+KKK66IvXv3uv0CABR2ZSQior6+PlavXh1LliyJpUuXxqZNm6KnpyfWrFkTERGrVq2K+fPnR1NTU5SWlsZFF100aP706dMjIk44DgBMTgXHyIoVK+LIkSOxfv366OjoiIULF8bOnTvzb2o9dOhQTJniB7sCAKenKMuyLPUiPkp3d3eUl5dHV1dXlJWVjei5q9a1jOj5AGC8Obhx2aic93S/f7uEAQAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSw4qR5ubmqKqqitLS0qiuro7du3efdOyWLVvi8ssvjxkzZsSMGTOitrb2lOMBgMml4BjZvn171NfXR2NjY+zZsycWLFgQdXV18fbbbw85fteuXXHdddfFv/3bv0VbW1tUVlbGF7/4xXjzzTc/9uIBgPGvKMuyrJAJ1dXV8bnPfS4eeuihiIgYGBiIysrK+NrXvhbr1q37yPn9/f0xY8aMeOihh2LVqlWn9Zzd3d1RXl4eXV1dUVZWVshyP1LVupYRPR8AjDcHNy4blfOe7vfvgq6M9PX1RXt7e9TW1n54gilTora2Ntra2k7rHO+99168//77MXPmzJOO6e3tje7u7kEPAGBiKihGjh49Gv39/VFRUTHoeEVFRXR0dJzWOW6//faYN2/eoKD5bU1NTVFeXp5/VFZWFrJMAGAcGdNP02zcuDG2bdsWP/zhD6O0tPSk4xoaGqKrqyv/OHz48BiuEgAYS1MLGTxr1qwoLi6Ozs7OQcc7Oztjzpw5p5x7//33x8aNG+MnP/lJXHLJJaccm8vlIpfLFbI0AGCcKujKSElJSSxevDhaW1vzxwYGBqK1tTVqampOOu++++6Le++9N3bu3BlLliwZ/moBgAmnoCsjERH19fWxevXqWLJkSSxdujQ2bdoUPT09sWbNmoiIWLVqVcyfPz+ampoiIuLv//7vY/369fHYY49FVVVV/r0ln/zkJ+OTn/zkCL4UAGA8KjhGVqxYEUeOHIn169dHR0dHLFy4MHbu3Jl/U+uhQ4diypQPL7h85zvfib6+vvjzP//zQedpbGyMr3/96x9v9QDAuFfwzxlJwc8ZAYDRM65+zggAwEgTIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJDStGmpubo6qqKkpLS6O6ujp27959yvFPPvlkXHDBBVFaWhoXX3xx7NixY1iLBQAmnoJjZPv27VFfXx+NjY2xZ8+eWLBgQdTV1cXbb7895PgXX3wxrrvuurjhhhvi5ZdfjmuuuSauueaa+I//+I+PvXgAYPwryrIsK2RCdXV1fO5zn4uHHnooIiIGBgaisrIyvva1r8W6detOGL9ixYro6emJH//4x/ljf/zHfxwLFy6MzZs3n9Zzdnd3R3l5eXR1dUVZWVkhy/1IVetaRvR8ADDeHNy4bFTOe7rfv6cWctK+vr5ob2+PhoaG/LEpU6ZEbW1ttLW1DTmnra0t6uvrBx2rq6uLZ5555qTP09vbG729vfmvu7q6IuI3L2qkDfS+N+LnBIDxZDS+v/7/837UdY+CYuTo0aPR398fFRUVg45XVFTEa6+9NuScjo6OIcd3dHSc9HmamprinnvuOeF4ZWVlIcsFAE5D+abRPf+xY8eivLz8pH9eUIyMlYaGhkFXUwYGBuKdd96J3/md34mioqIRe57u7u6orKyMw4cPj/jtHwpnP8489uTMYj/OLPbjo2VZFseOHYt58+adclxBMTJr1qwoLi6Ozs7OQcc7Oztjzpw5Q86ZM2dOQeMjInK5XORyuUHHpk+fXshSC1JWVuYv0hnEfpx57MmZxX6cWezHqZ3qisgHCvo0TUlJSSxevDhaW1vzxwYGBqK1tTVqamqGnFNTUzNofETEc889d9LxAMDkUvBtmvr6+li9enUsWbIkli5dGps2bYqenp5Ys2ZNRESsWrUq5s+fH01NTRERceutt8YXvvCFeOCBB2LZsmWxbdu2eOmll+Lhhx8e2VcCAIxLBcfIihUr4siRI7F+/fro6OiIhQsXxs6dO/NvUj106FBMmfLhBZdLL700HnvssbjrrrvijjvuiD/4gz+IZ555Ji666KKRexXDlMvlorGx8YRbQqRhP8489uTMYj/OLPZj5BT8c0YAAEaS300DACQlRgCApMQIAJCUGAEAkprwMdLc3BxVVVVRWloa1dXVsXv37lOOf/LJJ+OCCy6I0tLSuPjii2PHjh1jtNLJoZD92LJlS1x++eUxY8aMmDFjRtTW1n7k/lG4Qv+NfGDbtm1RVFQU11xzzegucJIpdD9+9atfxdq1a2Pu3LmRy+Xi/PPP9/+tEVTofmzatCk+85nPxNlnnx2VlZVx2223xa9//esxWu04lk1g27Zty0pKSrKtW7dm//mf/5nddNNN2fTp07POzs4hx//sZz/LiouLs/vuuy/7+c9/nt11113ZWWedlb3yyitjvPKJqdD9uP7667Pm5ubs5Zdfzl599dXsL/7iL7Ly8vLsv//7v8d45RNXoXvygQMHDmTz58/PLr/88uzP/uzPxmaxk0Ch+9Hb25stWbIku/rqq7MXXnghO3DgQLZr165s7969Y7zyianQ/fjBD36Q5XK57Ac/+EF24MCB7Nlnn83mzp2b3XbbbWO88vFnQsfI0qVLs7Vr1+a/7u/vz+bNm5c1NTUNOf7aa6/Nli1bNuhYdXV19ld/9Vejus7JotD9+G3Hjx/Ppk2bln3ve98brSVOOsPZk+PHj2eXXnpp9t3vfjdbvXq1GBlBhe7Hd77znezcc8/N+vr6xmqJk0qh+7F27drsT//0Twcdq6+vzy677LJRXedEMGFv0/T19UV7e3vU1tbmj02ZMiVqa2ujra1tyDltbW2DxkdE1NXVnXQ8p284+/Hb3nvvvXj//fdj5syZo7XMSWW4e/KNb3wjZs+eHTfccMNYLHPSGM5+/OhHP4qamppYu3ZtVFRUxEUXXRQbNmyI/v7+sVr2hDWc/bj00kujvb09fytn//79sWPHjrj66qvHZM3j2Rn5W3tHwtGjR6O/vz//k2E/UFFREa+99tqQczo6OoYc39HRMWrrnCyGsx+/7fbbb4958+adEIwMz3D25IUXXohHHnkk9u7dOwYrnFyGsx/79++Pn/70p/HlL385duzYEfv27YuvfvWr8f7770djY+NYLHvCGs5+XH/99XH06NH4/Oc/H1mWxfHjx+OWW26JO+64YyyWPK5N2CsjTCwbN26Mbdu2xQ9/+MMoLS1NvZxJ6dixY7Fy5crYsmVLzJo1K/VyiN/8otLZs2fHww8/HIsXL44VK1bEnXfeGZs3b069tElp165dsWHDhvj2t78de/bsiaeffjpaWlri3nvvTb20M96EvTIya9asKC4ujs7OzkHHOzs7Y86cOUPOmTNnTkHjOX3D2Y8P3H///bFx48b4yU9+EpdccsloLnNSKXRPfvGLX8TBgwdj+fLl+WMDAwMRETF16tR4/fXX47zzzhvdRU9gw/k3Mnfu3DjrrLOiuLg4f+zCCy+Mjo6O6Ovri5KSklFd80Q2nP24++67Y+XKlXHjjTdGRMTFF18cPT09cfPNN8edd9456Pe2MdiE/S9TUlISixcvjtbW1vyxgYGBaG1tjZqamiHn1NTUDBofEfHcc8+ddDynbzj7ERFx3333xb333hs7d+6MJUuWjMVSJ41C9+SCCy6IV155Jfbu3Zt/fOlLX4orrrgi9u7dG5WVlWO5/AlnOP9GLrvssti3b18+CiMi3njjjZg7d64Q+ZiGsx/vvffeCcHxQShmfg3cqaV+B+1o2rZtW5bL5bJHH300+/nPf57dfPPN2fTp07OOjo4sy7Js5cqV2bp16/Ljf/azn2VTp07N7r///uzVV1/NGhsbfbR3BBW6Hxs3bsxKSkqyp556Knvrrbfyj2PHjqV6CRNOoXvy23yaZmQVuh+HDh3Kpk2blv31X/919vrrr2c//vGPs9mzZ2ff/OY3U72ECaXQ/WhsbMymTZuWPf7449n+/fuzf/3Xf83OO++87Nprr031EsaNCR0jWZZl//RP/5Sdc845WUlJSbZ06dLs3//93/N/9oUvfCFbvXr1oPFPPPFEdv7552clJSXZH/7hH2YtLS1jvOKJrZD9+L3f+70sIk54NDY2jv3CJ7BC/438f2Jk5BW6Hy+++GJWXV2d5XK57Nxzz82+9a1vZcePHx/jVU9chezH+++/n33961/PzjvvvKy0tDSrrKzMvvrVr2b/8z//M/YLH2eKssy1IwAgnQn7nhEAYHwQIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEn9H6cW6TN2VdVNAAAAAElFTkSuQmCC\n"
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "<Figure size 640x480 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGdCAYAAADAAnMpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAfpUlEQVR4nO3df2xV5eHH8U9/0FuM9AKrvS3suioZogMLKdIVJIq72iCp44/FRgytjeDQahg3y2wFWhGlnQPWZRYbUIZmYlGizNmmyOoIQboQCk2YAg6BtUPuhcbZW8psofd8/zBev5UWemrbh3v7fiXnj54+55zn7hn2nXN/RVmWZQkAAMCQaNMTAAAAwxsxAgAAjCJGAACAUcQIAAAwihgBAABGESMAAMAoYgQAABhFjAAAAKNiTU+gL4LBoD7//HONGjVKUVFRpqcDAAD6wLIstbW1ady4cYqO7v3+R1jEyOeffy632216GgAAoB+am5v1wx/+sNffh0WMjBo1StLXDyYhIcHwbAAAQF8EAgG53e7Q3/HehEWMfPPUTEJCAjECAECYudpLLHgBKwAAMMp2jOzZs0fZ2dkaN26coqKitGPHjj4f+9FHHyk2NlZTp061e1kAABChbMdIe3u70tLSVFFRYeu4L7/8Urm5ufrZz35m95IAACCC2X7NyNy5czV37lzbF1qyZIkWLFigmJgYW3dTAABAZBuS14z86U9/0okTJ1RSUtKn8R0dHQoEAt02AAAQmQY9Rv71r3+psLBQf/7znxUb27cbMaWlpXI6naGNzxgBACByDWqMdHV1acGCBVq1apUmTpzY5+OKiorU2toa2pqbmwdxlgAAwKRB/ZyRtrY2HThwQIcOHdKTTz4p6euPdrcsS7Gxsfrggw90zz33XHacw+GQw+EYzKkBAIBrxKDGSEJCgg4fPtxt34YNG/Thhx9q+/btuummmwbz8gAAIAzYjpHz58/r+PHjoZ9PnjypxsZGjR07VjfeeKOKiop0+vRpvf7664qOjtbkyZO7HZ+UlKT4+PjL9gMAgOHJdowcOHBAc+bMCf3s9XolSXl5edqyZYvOnDmjpqamgZshAACIaFGWZVmmJ3E1gUBATqdTra2tfDcNAABhoq9/v/luGgAAYBQxAgAAjBrUd9MAABDJUgurTU9hQJwqm2f0+twZAQAARhEjAADAKGIEAAAYRYwAAACjiBEAAGAUMQIAAIwiRgAAgFHECAAAMIoYAQAARhEjAADAKGIEAAAYRYwAAACjiBEAAGAUMQIAAIwiRgAAgFHECAAAMIoYAQAARhEjAADAKGIEAAAYFWt6AgAwFFILq01PYUCcKptnegrAgOPOCAAAMIoYAQAARhEjAADAKGIEAAAYRYwAAACjiBEAAGAUMQIAAIwiRgAAgFHECAAAMIoYAQAARhEjAADAKGIEAAAYRYwAAACjiBEAAGAUMQIAAIwiRgAAgFHECAAAMMp2jOzZs0fZ2dkaN26coqKitGPHjiuOf+edd3TvvffqhhtuUEJCgjIzM7Vz587+zhcAAEQY2zHS3t6utLQ0VVRU9Gn8nj17dO+996qmpkYNDQ2aM2eOsrOzdejQIduTBQAAkSfW7gFz587V3Llz+zy+vLy8289r1qzRX/7yF/31r3/VtGnT7F4eAABEGNsx8n0Fg0G1tbVp7NixvY7p6OhQR0dH6OdAIDAUUwMAAAYM+QtY165dq/Pnz+vBBx/sdUxpaamcTmdoc7vdQzhDAAAwlIY0RrZu3apVq1bprbfeUlJSUq/jioqK1NraGtqam5uHcJYAAGAoDdnTNFVVVVq0aJHefvtteTyeK451OBxyOBxDNDMAAGDSkNwZefPNN5Wfn68333xT8+bNG4pLAgCAMGH7zsj58+d1/Pjx0M8nT55UY2Ojxo4dqxtvvFFFRUU6ffq0Xn/9dUlfPzWTl5enP/zhD8rIyJDP55MkjRw5Uk6nc4AeBgAACFe274wcOHBA06ZNC70t1+v1atq0aSouLpYknTlzRk1NTaHxGzdu1KVLl1RQUKCUlJTQtnTp0gF6CAAAIJzZvjNy9913y7KsXn+/ZcuWbj/v3r3b7iUAAMAwwnfTAAAAo4gRAABgFDECAACMIkYAAIBRxAgAADCKGAEAAEYRIwAAwChiBAAAGDVkX5QH9EVqYbXpKQyYU2V8DxMA9AV3RgAAgFHECAAAMIoYAQAARhEjAADAKGIEAAAYRYwAAACjiBEAAGAUMQIAAIwiRgAAgFHECAAAMIoYAQAARhEjAADAKGIEAAAYRYwAAACjiBEAAGAUMQIAAIwiRgAAgFHECAAAMCrW9ASASJVaWG16CgPiVNk801MAEOG4MwIAAIwiRgAAgFE8TRMhIuUpAQBXxr91RCLujAAAAKOIEQAAYBQxAgAAjCJGAACAUcQIAAAwihgBAABGESMAAMAoYgQAABhFjAAAAKOIEQAAYJTtj4Pfs2ePfve736mhoUFnzpzRu+++q/nz51/xmN27d8vr9erjjz+W2+3WihUr9Mgjj/RzygOLj1YGAMAs23dG2tvblZaWpoqKij6NP3nypObNm6c5c+aosbFRv/rVr7Ro0SLt3LnT9mQBAEDksX1nZO7cuZo7d26fx1dWVuqmm27SunXrJEm33nqr9u7dq9///vfKysqye3kAABBhBv01I/X19fJ4PN32ZWVlqb6+vtdjOjo6FAgEum0AACAyDXqM+Hw+uVyubvtcLpcCgYD+97//9XhMaWmpnE5naHO73YM9TQAAYMg1+W6aoqIitba2hrbm5mbTUwIAAIPE9mtG7EpOTpbf7++2z+/3KyEhQSNHjuzxGIfDIYfDMdhTAwAA14BBvzOSmZmpurq6bvt27dqlzMzMwb40AAAIA7Zj5Pz582psbFRjY6Okr9+629jYqKamJklfP8WSm5sbGr9kyRKdOHFCv/nNb3T06FFt2LBBb731lpYtWzYwjwAAAIQ12zFy4MABTZs2TdOmTZMkeb1eTZs2TcXFxZKkM2fOhMJEkm666SZVV1dr165dSktL07p16/TKK6/wtl4AACCpH68Zufvuu2VZVq+/37JlS4/HHDp0yO6lAADAMHBNvpsGAAAMH8QIAAAwihgBAABGESMAAMAoYgQAABhFjAAAAKMG/ePgAYS31MJq01MAEOG4MwIAAIwiRgAAgFHECAAAMIoYAQAARhEjAADAKGIEAAAYRYwAAACjiBEAAGAUMQIAAIwiRgAAgFHECAAAMIoYAQAARhEjAADAKGIEAAAYRYwAAACjiBEAAGAUMQIAAIwiRgAAgFHECAAAMIoYAQAARhEjAADAKGIEAAAYRYwAAACjiBEAAGAUMQIAAIwiRgAAgFHECAAAMIoYAQAARhEjAADAKGIEAAAYRYwAAACjiBEAAGAUMQIAAIwiRgAAgFH9ipGKigqlpqYqPj5eGRkZ2r9//xXHl5eX65ZbbtHIkSPldru1bNkyffXVV/2aMAAAiCy2Y2Tbtm3yer0qKSnRwYMHlZaWpqysLJ09e7bH8Vu3blVhYaFKSkp05MgRvfrqq9q2bZueeeaZ7z15AAAQ/mzHyPr167V48WLl5+frtttuU2Vlpa677jpt3ry5x/H79u3TrFmztGDBAqWmpuq+++7TQw89dNW7KQAAYHiwFSOdnZ1qaGiQx+P59gTR0fJ4PKqvr+/xmJkzZ6qhoSEUHydOnFBNTY3uv//+Xq/T0dGhQCDQbQMAAJEp1s7glpYWdXV1yeVyddvvcrl09OjRHo9ZsGCBWlpadOedd8qyLF26dElLliy54tM0paWlWrVqlZ2pAQCAMDXo76bZvXu31qxZow0bNujgwYN65513VF1drdWrV/d6TFFRkVpbW0Nbc3PzYE8TAAAYYuvOSGJiomJiYuT3+7vt9/v9Sk5O7vGYlStXauHChVq0aJEkacqUKWpvb9djjz2m5cuXKzr68h5yOBxyOBx2pgYAAMKUrTsjcXFxSk9PV11dXWhfMBhUXV2dMjMzezzmwoULlwVHTEyMJMmyLLvzBQAAEcbWnRFJ8nq9ysvL0/Tp0zVjxgyVl5ervb1d+fn5kqTc3FyNHz9epaWlkqTs7GytX79e06ZNU0ZGho4fP66VK1cqOzs7FCUAAGD4sh0jOTk5OnfunIqLi+Xz+TR16lTV1taGXtTa1NTU7U7IihUrFBUVpRUrVuj06dO64YYblJ2drRdeeGHgHgUAAAhbUVYYPFcSCATkdDrV2tqqhISEAT13amH1gJ4PAIBwc6ps3qCct69/v/luGgAAYBQxAgAAjCJGAACAUcQIAAAwihgBAABGESMAAMAoYgQAABhFjAAAAKOIEQAAYBQxAgAAjCJGAACAUcQIAAAwihgBAABGESMAAMAoYgQAABhFjAAAAKOIEQAAYBQxAgAAjCJGAACAUcQIAAAwihgBAABGESMAAMAoYgQAABhFjAAAAKOIEQAAYBQxAgAAjCJGAACAUcQIAAAwihgBAABGESMAAMAoYgQAABhFjAAAAKOIEQAAYBQxAgAAjCJGAACAUcQIAAAwihgBAABGESMAAMAoYgQAABhFjAAAAKP6FSMVFRVKTU1VfHy8MjIytH///iuO//LLL1VQUKCUlBQ5HA5NnDhRNTU1/ZowAACILLF2D9i2bZu8Xq8qKyuVkZGh8vJyZWVl6dixY0pKSrpsfGdnp+69914lJSVp+/btGj9+vP79739r9OjRAzF/AAAQ5mzHyPr167V48WLl5+dLkiorK1VdXa3NmzersLDwsvGbN2/WF198oX379mnEiBGSpNTU1O83awAAEDFsPU3T2dmphoYGeTyeb08QHS2Px6P6+voej3nvvfeUmZmpgoICuVwuTZ48WWvWrFFXV1ev1+no6FAgEOi2AQCAyGQrRlpaWtTV1SWXy9Vtv8vlks/n6/GYEydOaPv27erq6lJNTY1WrlypdevW6fnnn+/1OqWlpXI6naHN7XbbmSYAAAgjg/5ummAwqKSkJG3cuFHp6enKycnR8uXLVVlZ2esxRUVFam1tDW3Nzc2DPU0AAGCIrdeMJCYmKiYmRn6/v9t+v9+v5OTkHo9JSUnRiBEjFBMTE9p36623yufzqbOzU3FxcZcd43A45HA47EwNAACEKVt3RuLi4pSenq66urrQvmAwqLq6OmVmZvZ4zKxZs3T8+HEFg8HQvk8//VQpKSk9hggAABhebD9N4/V6tWnTJr322ms6cuSIHn/8cbW3t4feXZObm6uioqLQ+Mcff1xffPGFli5dqk8//VTV1dVas2aNCgoKBu5RAACAsGX7rb05OTk6d+6ciouL5fP5NHXqVNXW1oZe1NrU1KTo6G8bx+12a+fOnVq2bJluv/12jR8/XkuXLtXTTz89cI8CAACErSjLsizTk7iaQCAgp9Op1tZWJSQkDOi5UwurB/R8AACEm1Nl8wblvH39+8130wAAAKOIEQAAYBQxAgAAjCJGAACAUcQIAAAwihgBAABGESMAAMAoYgQAABhFjAAAAKOIEQAAYBQxAgAAjCJGAACAUcQIAAAwihgBAABGESMAAMAoYgQAABhFjAAAAKOIEQAAYBQxAgAAjCJGAACAUcQIAAAwihgBAABGESMAAMAoYgQAABhFjAAAAKOIEQAAYBQxAgAAjCJGAACAUcQIAAAwihgBAABGESMAAMAoYgQAABhFjAAAAKOIEQAAYBQxAgAAjCJGAACAUcQIAAAwihgBAABGESMAAMAoYgQAABhFjAAAAKP6FSMVFRVKTU1VfHy8MjIytH///j4dV1VVpaioKM2fP78/lwUAABHIdoxs27ZNXq9XJSUlOnjwoNLS0pSVlaWzZ89e8bhTp07p17/+tWbPnt3vyQIAgMhjO0bWr1+vxYsXKz8/X7fddpsqKyt13XXXafPmzb0e09XVpYcfflirVq3SzTff/L0mDAAAIoutGOns7FRDQ4M8Hs+3J4iOlsfjUX19fa/HPffcc0pKStKjjz7ap+t0dHQoEAh02wAAQGSyFSMtLS3q6uqSy+Xqtt/lcsnn8/V4zN69e/Xqq69q06ZNfb5OaWmpnE5naHO73XamCQAAwsigvpumra1NCxcu1KZNm5SYmNjn44qKitTa2hrampubB3GWAADApFg7gxMTExUTEyO/399tv9/vV3Jy8mXjP/vsM506dUrZ2dmhfcFg8OsLx8bq2LFjmjBhwmXHORwOORwOO1MDAABhytadkbi4OKWnp6uuri60LxgMqq6uTpmZmZeNnzRpkg4fPqzGxsbQ9sADD2jOnDlqbGzk6RcAAGDvzogkeb1e5eXlafr06ZoxY4bKy8vV3t6u/Px8SVJubq7Gjx+v0tJSxcfHa/Lkyd2OHz16tCRdth8AAAxPtmMkJydH586dU3FxsXw+n6ZOnara2trQi1qbmpoUHc0HuwIAgL6JsizLMj2JqwkEAnI6nWptbVVCQsKAnju1sHpAzwcAQLg5VTZvUM7b17/f3MIAAABGESMAAMAoYgQAABhFjAAAAKOIEQAAYBQxAgAAjCJGAACAUcQIAAAwihgBAABGESMAAMAoYgQAABhFjAAAAKOIEQAAYBQxAgAAjCJGAACAUcQIAAAwihgBAABGESMAAMAoYgQAABhFjAAAAKOIEQAAYBQxAgAAjCJGAACAUcQIAAAwihgBAABGESMAAMAoYgQAABhFjAAAAKOIEQAAYBQxAgAAjCJGAACAUcQIAAAwihgBAABGESMAAMAoYgQAABhFjAAAAKOIEQAAYBQxAgAAjCJGAACAUcQIAAAwql8xUlFRodTUVMXHxysjI0P79+/vdeymTZs0e/ZsjRkzRmPGjJHH47nieAAAMLzYjpFt27bJ6/WqpKREBw8eVFpamrKysnT27Nkex+/evVsPPfSQ/v73v6u+vl5ut1v33XefTp8+/b0nDwAAwl+UZVmWnQMyMjJ0xx136KWXXpIkBYNBud1uPfXUUyosLLzq8V1dXRozZoxeeukl5ebm9umagUBATqdTra2tSkhIsDPdq0otrB7Q8wEAEG5Olc0blPP29e+3rTsjnZ2damhokMfj+fYE0dHyeDyqr6/v0zkuXLigixcvauzYsXYuDQAAIlSsncEtLS3q6uqSy+Xqtt/lcuno0aN9OsfTTz+tcePGdQua7+ro6FBHR0fo50AgYGeaAAAgjAzpu2nKyspUVVWld999V/Hx8b2OKy0tldPpDG1ut3sIZwkAAIaSrRhJTExUTEyM/H5/t/1+v1/JyclXPHbt2rUqKyvTBx98oNtvv/2KY4uKitTa2hrampub7UwTAACEEVsxEhcXp/T0dNXV1YX2BYNB1dXVKTMzs9fjXnzxRa1evVq1tbWaPn36Va/jcDiUkJDQbQMAAJHJ1mtGJMnr9SovL0/Tp0/XjBkzVF5ervb2duXn50uScnNzNX78eJWWlkqSfvvb36q4uFhbt25VamqqfD6fJOn666/X9ddfP4APBQAAhCPbMZKTk6Nz586puLhYPp9PU6dOVW1tbehFrU1NTYqO/vaGy8svv6zOzk794he/6HaekpISPfvss99v9gAAIOzZ/pwRE/icEQAABk9Yfc4IAADAQCNGAACAUcQIAAAwihgBAABGESMAAMAoYgQAABhFjAAAAKOIEQAAYBQxAgAAjCJGAACAUcQIAAAwihgBAABGESMAAMAoYgQAABhFjAAAAKOIEQAAYBQxAgAAjCJGAACAUcQIAAAwihgBAABGESMAAMAoYgQAABhFjAAAAKOIEQAAYBQxAgAAjCJGAACAUcQIAAAwihgBAABGESMAAMAoYgQAABhFjAAAAKOIEQAAYBQxAgAAjCJGAACAUcQIAAAwihgBAABGESMAAMAoYgQAABhFjAAAAKOIEQAAYBQxAgAAjOpXjFRUVCg1NVXx8fHKyMjQ/v37rzj+7bff1qRJkxQfH68pU6aopqamX5MFAACRx3aMbNu2TV6vVyUlJTp48KDS0tKUlZWls2fP9jh+3759euihh/Too4/q0KFDmj9/vubPn69//vOf33vyAAAg/EVZlmXZOSAjI0N33HGHXnrpJUlSMBiU2+3WU089pcLCwsvG5+TkqL29Xe+//35o309/+lNNnTpVlZWVfbpmIBCQ0+lUa2urEhIS7Ez3qlILqwf0fAAAhJtTZfMG5bx9/fsda+eknZ2damhoUFFRUWhfdHS0PB6P6uvrezymvr5eXq+3276srCzt2LGj1+t0dHSoo6Mj9HNra6ukrx/UQAt2XBjwcwIAEE4G4+/r/z/v1e572IqRlpYWdXV1yeVyddvvcrl09OjRHo/x+Xw9jvf5fL1ep7S0VKtWrbpsv9vttjNdAADQB87ywT1/W1ubnE5nr7+3FSNDpaioqNvdlGAwqC+++EI/+MEPFBUVNWDXCQQCcrvdam5uHvCnf2Af63HtYU2uLazHtYX1uDrLstTW1qZx48ZdcZytGElMTFRMTIz8fn+3/X6/X8nJyT0ek5ycbGu8JDkcDjkcjm77Ro8ebWeqtiQkJPB/pGsI63HtYU2uLazHtYX1uLIr3RH5hq1308TFxSk9PV11dXWhfcFgUHV1dcrMzOzxmMzMzG7jJWnXrl29jgcAAMOL7adpvF6v8vLyNH36dM2YMUPl5eVqb29Xfn6+JCk3N1fjx49XaWmpJGnp0qW66667tG7dOs2bN09VVVU6cOCANm7cOLCPBAAAhCXbMZKTk6Nz586puLhYPp9PU6dOVW1tbehFqk1NTYqO/vaGy8yZM7V161atWLFCzzzzjH784x9rx44dmjx58sA9in5yOBwqKSm57CkhmMF6XHtYk2sL63FtYT0Gju3PGQEAABhIfDcNAAAwihgBAABGESMAAMAoYgQAABgV8TFSUVGh1NRUxcfHKyMjQ/v377/i+LfffluTJk1SfHy8pkyZopqamiGa6fBgZz02bdqk2bNna8yYMRozZow8Hs9V1w/22f038o2qqipFRUVp/vz5gzvBYcbuenz55ZcqKChQSkqKHA6HJk6cyH+3BpDd9SgvL9ctt9yikSNHyu12a9myZfrqq6+GaLZhzIpgVVVVVlxcnLV582br448/thYvXmyNHj3a8vv9PY7/6KOPrJiYGOvFF1+0PvnkE2vFihXWiBEjrMOHDw/xzCOT3fVYsGCBVVFRYR06dMg6cuSI9cgjj1hOp9P6z3/+M8Qzj1x21+QbJ0+etMaPH2/Nnj3b+vnPfz40kx0G7K5HR0eHNX36dOv++++39u7da508edLavXu31djYOMQzj0x21+ONN96wHA6H9cYbb1gnT560du7caaWkpFjLli0b4pmHn4iOkRkzZlgFBQWhn7u6uqxx48ZZpaWlPY5/8MEHrXnz5nXbl5GRYf3yl78c1HkOF3bX47suXbpkjRo1ynrttdcGa4rDTn/W5NKlS9bMmTOtV155xcrLyyNGBpDd9Xj55Zetm2++2ers7ByqKQ4rdtejoKDAuueee7rt83q91qxZswZ1npEgYp+m6ezsVENDgzweT2hfdHS0PB6P6uvrezymvr6+23hJysrK6nU8+q4/6/FdFy5c0MWLFzV27NjBmuaw0t81ee6555SUlKRHH310KKY5bPRnPd577z1lZmaqoKBALpdLkydP1po1a9TV1TVU045Y/VmPmTNnqqGhIfRUzokTJ1RTU6P7779/SOYczq7Jb+0dCC0tLerq6gp9Muw3XC6Xjh492uMxPp+vx/E+n2/Q5jlc9Gc9vuvpp5/WuHHjLgtG9E9/1mTv3r169dVX1djYOAQzHF76sx4nTpzQhx9+qIcfflg1NTU6fvy4nnjiCV28eFElJSVDMe2I1Z/1WLBggVpaWnTnnXfKsixdunRJS5Ys0TPPPDMUUw5rEXtnBJGlrKxMVVVVevfddxUfH296OsNSW1ubFi5cqE2bNikxMdH0dKCvv6g0KSlJGzduVHp6unJycrR8+XJVVlaantqwtHv3bq1Zs0YbNmzQwYMH9c4776i6ulqrV682PbVrXsTeGUlMTFRMTIz8fn+3/X6/X8nJyT0ek5ycbGs8+q4/6/GNtWvXqqysTH/72990++23D+Y0hxW7a/LZZ5/p1KlTys7ODu0LBoOSpNjYWB07dkwTJkwY3ElHsP78G0lJSdGIESMUExMT2nfrrbfK5/Ops7NTcXFxgzrnSNaf9Vi5cqUWLlyoRYsWSZKmTJmi9vZ2PfbYY1q+fHm3721DdxH7v0xcXJzS09NVV1cX2hcMBlVXV6fMzMwej8nMzOw2XpJ27drV63j0XX/WQ5JefPFFrV69WrW1tZo+ffpQTHXYsLsmkyZN0uHDh9XY2BjaHnjgAc2ZM0eNjY1yu91DOf2I059/I7NmzdLx48dDUShJn376qVJSUgiR76k/63HhwoXLguObULT4GrgrM/0K2sFUVVVlORwOa8uWLdYnn3xiPfbYY9bo0aMtn89nWZZlLVy40CosLAyN/+ijj6zY2Fhr7dq11pEjR6ySkhLe2juA7K5HWVmZFRcXZ23fvt06c+ZMaGtrazP1ECKO3TX5Lt5NM7DsrkdTU5M1atQo68knn7SOHTtmvf/++1ZSUpL1/PPPm3oIEcXuepSUlFijRo2y3nzzTevEiRPWBx98YE2YMMF68MEHTT2EsBHRMWJZlvXHP/7RuvHGG624uDhrxowZ1j/+8Y/Q7+666y4rLy+v2/i33nrLmjhxohUXF2f95Cc/saqrq4d4xpHNznr86Ec/siRdtpWUlAz9xCOY3X8j/x8xMvDsrse+ffusjIwMy+FwWDfffLP1wgsvWJcuXRriWUcuO+tx8eJF69lnn7UmTJhgxcfHW26323riiSes//73v0M/8TATZVncOwIAAOZE7GtGAABAeCBGAACAUcQIAAAwihgBAABGESMAAMAoYgQAABhFjAAAAKOIEQAAYBQxAgAAjCJGAACAUcQIAAAwihgBAABG/R+8v4EvNjwwkwAAAABJRU5ErkJggg==\n"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# proof that this works with different models\n",
    "\n",
    "def prior():\n",
    "    u1 = d.norm(0,1)\n",
    "    u2 = d.norm(0,1)\n",
    "    w = d.unif(0,1)\n",
    "    return u1,u2,w\n",
    "\n",
    "def likelihood(u1,u2,w):\n",
    "    x1 = d.norm(u1,1)\n",
    "    x2 = d.norm(x1 + w*u2, w)\n",
    "    x3 = d.norm(x2 + w*u2, w)\n",
    "    x4 = d.norm(x3 + w*u2, w)\n",
    "    return x1, x2, x3, x4\n",
    "\n",
    "SBC(prior,likelihood)"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Alternatively, I guess you might just pass the random variables directly?"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "outputs": [
    {
     "data": {
      "text/plain": "<Figure size 640x480 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGdCAYAAADAAnMpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAeMElEQVR4nO3dbWxd9X3A8Z/t4Guqxk5YGjvJ3LmwUeiAJEsaz1BUmEwjiNLxYiKCKskyHkabVSzWNGIe4lLaOENAM43QiJSITip1KKKsa6Iw6i1CFE8RDpbYyoPSJEtGsUnEagez2sQ+e4G4zI0Tco3tf+x8PtJ94eP/Ofd3OTj+6tx7fYuyLMsCACCR4tQDAABnNjECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJTUk9wKkYHByMX/3qVzF16tQoKipKPQ4AcAqyLIujR4/G7Nmzo7j4xNc/JkSM/OpXv4rq6urUYwAAI3Do0KH43d/93RN+f0LEyNSpUyPi/QdTXl6eeBoA4FT09PREdXV1/vf4iUyIGPngqZny8nIxAgATzEe9xMILWAGApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASU1JPQAAp65m7fbUI4yKAxuWpB6B04grIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAIKmCY+S5556LpUuXxuzZs6OoqCiefvrpk65/6qmn4qqrropPfepTUV5eHnV1dfHMM8+MdF4AYJIpOEZ6e3tj7ty5sWnTplNa/9xzz8VVV10VO3bsiPb29rjyyitj6dKl8dJLLxU8LAAw+UwpdIerr746rr766lNev3HjxiFfr1+/Pv7pn/4p/vmf/znmz59f6N0DAJPMuL9mZHBwMI4ePRrnnHPOeN81AHAaKvjKyMd1//33xzvvvBPXXXfdCdf09fVFX19f/uuenp7xGA0ASGBcr4w8/vjjcc8998QTTzwRM2fOPOG65ubmqKioyN+qq6vHcUoAYDyNW4y0tLTETTfdFE888UTU19efdG1jY2N0d3fnb4cOHRqnKQGA8TYuT9P88Ic/jL/4i7+IlpaWWLJkyUeuz+VykcvlxmEyTjc1a7enHmHUHNjw0f+vAxPbZPk3K/W/VwXHyDvvvBN79+7Nf71///7o6OiIc845Jz796U9HY2NjvPHGG/GP//iPEfH+UzMrV66Mv//7v4/a2tro7OyMiIizzz47KioqRulhAAATVcFP07z44osxf/78/NtyGxoaYv78+bFu3bqIiHjzzTfj4MGD+fWPPPJIHDt2LFavXh2zZs3K32677bZReggAwERW8JWRK664IrIsO+H3H3vssSFf79q1q9C7AADOID6bBgBISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBS4/LZNADw/02Wz3RhdLgyAgAkJUYAgKTECACQlBgBAJISIwBAUt5NM0l4ZTpjZbL8v3Vgw5LUIwAn4MoIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQ1JfUAqdWs3Z56BGAc+FmH05crIwBAUmIEAEhKjAAASYkRACApMQIAJFVwjDz33HOxdOnSmD17dhQVFcXTTz/9kfvs2rUr/uiP/ihyuVz8/u//fjz22GMjGBUAmIwKjpHe3t6YO3dubNq06ZTW79+/P5YsWRJXXnlldHR0xF//9V/HTTfdFM8880zBwwIAk0/Bf2fk6quvjquvvvqU12/evDk+85nPxAMPPBARERdeeGE8//zz8Z3vfCcWL15c6N0DAJPMmL9mpK2tLerr64dsW7x4cbS1tZ1wn76+vujp6RlyAwAmpzGPkc7OzqisrByyrbKyMnp6euJ///d/h92nubk5Kioq8rfq6uqxHhMASOS0fDdNY2NjdHd352+HDh1KPRIAMEbG/LNpqqqqoqura8i2rq6uKC8vj7PPPnvYfXK5XORyubEeDQA4DYz5lZG6urpobW0dsu3ZZ5+Nurq6sb5rAGACKDhG3nnnnejo6IiOjo6IeP+tux0dHXHw4MGIeP8plhUrVuTX33rrrbFv377427/923j11Vfj4YcfjieeeCLWrFkzOo8AAJjQCo6RF198MebPnx/z58+PiIiGhoaYP39+rFu3LiIi3nzzzXyYRER85jOfie3bt8ezzz4bc+fOjQceeCC+973veVsvABARI3jNyBVXXBFZlp3w+8P9ddUrrrgiXnrppULvCgA4A5yW76YBAM4cYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASU1JPQBMVjVrt6ceAWBCcGUEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSGlGMbNq0KWpqaqKsrCxqa2tj9+7dJ12/cePG+OxnPxtnn312VFdXx5o1a+I3v/nNiAYGACaXgmNk27Zt0dDQEE1NTbFnz56YO3duLF68ON56661h1z/++OOxdu3aaGpqildeeSUeffTR2LZtW9xxxx0fe3gAYOIrOEYefPDBuPnmm2PVqlXxuc99LjZv3hyf+MQnYuvWrcOuf+GFF+Kyyy6LG264IWpqauJLX/pSXH/99R95NQUAODMUFCP9/f3R3t4e9fX1Hx6guDjq6+ujra1t2H0uvfTSaG9vz8fHvn37YseOHXHNNdd8jLEBgMmioM+mOXLkSAwMDERlZeWQ7ZWVlfHqq68Ou88NN9wQR44ciS984QuRZVkcO3Ysbr311pM+TdPX1xd9fX35r3t6egoZEwCYQMb83TS7du2K9evXx8MPPxx79uyJp556KrZv3x733nvvCfdpbm6OioqK/K26unqsxwQAEinoysiMGTOipKQkurq6hmzv6uqKqqqqYfe5++67Y/ny5XHTTTdFRMTFF18cvb29ccstt8Sdd94ZxcXH91BjY2M0NDTkv+7p6REkADBJFXRlpLS0NBYsWBCtra35bYODg9Ha2hp1dXXD7vPuu+8eFxwlJSUREZFl2bD75HK5KC8vH3IDACangq6MREQ0NDTEypUrY+HChbFo0aLYuHFj9Pb2xqpVqyIiYsWKFTFnzpxobm6OiIilS5fGgw8+GPPnz4/a2trYu3dv3H333bF06dJ8lAAAZ66CY2TZsmVx+PDhWLduXXR2dsa8efNi586d+Re1Hjx4cMiVkLvuuiuKiorirrvuijfeeCM+9alPxdKlS+Pb3/726D0KAGDCKspO9FzJaaSnpycqKiqiu7t71J+yqVm7fVSPBwATzYENS8bkuKf6+9tn0wAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKRGFCObNm2KmpqaKCsri9ra2ti9e/dJ1//617+O1atXx6xZsyKXy8X5558fO3bsGNHAAMDkMqXQHbZt2xYNDQ2xefPmqK2tjY0bN8bixYvjtddei5kzZx63vr+/P6666qqYOXNmPPnkkzFnzpz4r//6r5g2bdpozA8ATHAFx8iDDz4YN998c6xatSoiIjZv3hzbt2+PrVu3xtq1a49bv3Xr1nj77bfjhRdeiLPOOisiImpqaj7e1ADApFHQ0zT9/f3R3t4e9fX1Hx6guDjq6+ujra1t2H1+8pOfRF1dXaxevToqKyvjoosuivXr18fAwMAJ76evry96enqG3ACAyamgGDly5EgMDAxEZWXlkO2VlZXR2dk57D779u2LJ598MgYGBmLHjh1x9913xwMPPBDf+ta3Tng/zc3NUVFRkb9VV1cXMiYAMIGM+btpBgcHY+bMmfHII4/EggULYtmyZXHnnXfG5s2bT7hPY2NjdHd352+HDh0a6zEBgEQKes3IjBkzoqSkJLq6uoZs7+rqiqqqqmH3mTVrVpx11llRUlKS33bhhRdGZ2dn9Pf3R2lp6XH75HK5yOVyhYwGAExQBV0ZKS0tjQULFkRra2t+2+DgYLS2tkZdXd2w+1x22WWxd+/eGBwczG97/fXXY9asWcOGCABwZin4aZqGhobYsmVLfP/7349XXnklvvrVr0Zvb2/+3TUrVqyIxsbG/PqvfvWr8fbbb8dtt90Wr7/+emzfvj3Wr18fq1evHr1HAQBMWAW/tXfZsmVx+PDhWLduXXR2dsa8efNi586d+Re1Hjx4MIqLP2yc6urqeOaZZ2LNmjVxySWXxJw5c+K2226L22+/ffQeBQAwYRVlWZalHuKj9PT0REVFRXR3d0d5efmoHrtm7fZRPR4ATDQHNiwZk+Oe6u9vn00DACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkNSIYmTTpk1RU1MTZWVlUVtbG7t37z6l/VpaWqKoqCiuvfbakdwtADAJFRwj27Zti4aGhmhqaoo9e/bE3LlzY/HixfHWW2+ddL8DBw7E3/zN38Tll18+4mEBgMmn4Bh58MEH4+abb45Vq1bF5z73udi8eXN84hOfiK1bt55wn4GBgfjKV74S99xzT5x77rkfa2AAYHIpKEb6+/ujvb096uvrPzxAcXHU19dHW1vbCff75je/GTNnzowbb7zxlO6nr68venp6htwAgMmpoBg5cuRIDAwMRGVl5ZDtlZWV0dnZOew+zz//fDz66KOxZcuWU76f5ubmqKioyN+qq6sLGRMAmEDG9N00R48ejeXLl8eWLVtixowZp7xfY2NjdHd352+HDh0awykBgJSmFLJ4xowZUVJSEl1dXUO2d3V1RVVV1XHrf/nLX8aBAwdi6dKl+W2Dg4Pv3/GUKfHaa6/Feeedd9x+uVwucrlcIaMBABNUQVdGSktLY8GCBdHa2prfNjg4GK2trVFXV3fc+gsuuCBefvnl6OjoyN++/OUvx5VXXhkdHR2efgEACrsyEhHR0NAQK1eujIULF8aiRYti48aN0dvbG6tWrYqIiBUrVsScOXOiubk5ysrK4qKLLhqy/7Rp0yIijtsOAJyZCo6RZcuWxeHDh2PdunXR2dkZ8+bNi507d+Zf1Hrw4MEoLvaHXQGAU1OUZVmWeoiP0tPTExUVFdHd3R3l5eWjeuyatdtH9XgAMNEc2LBkTI57qr+/XcIAAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJjShGNm3aFDU1NVFWVha1tbWxe/fuE67dsmVLXH755TF9+vSYPn161NfXn3Q9AHBmKThGtm3bFg0NDdHU1BR79uyJuXPnxuLFi+Ott94adv2uXbvi+uuvj3/7t3+Ltra2qK6uji996UvxxhtvfOzhAYCJryjLsqyQHWpra+Pzn/98PPTQQxERMTg4GNXV1fH1r3891q5d+5H7DwwMxPTp0+Ohhx6KFStWnNJ99vT0REVFRXR3d0d5eXkh436kmrXbR/V4ADDRHNiwZEyOe6q/vwu6MtLf3x/t7e1RX1//4QGKi6O+vj7a2tpO6RjvvvtuvPfee3HOOeeccE1fX1/09PQMuQEAk1NBMXLkyJEYGBiIysrKIdsrKyujs7PzlI5x++23x+zZs4cEzW9rbm6OioqK/K26urqQMQGACWRc302zYcOGaGlpiR//+MdRVlZ2wnWNjY3R3d2dvx06dGgcpwQAxtOUQhbPmDEjSkpKoqura8j2rq6uqKqqOum+999/f2zYsCF+9rOfxSWXXHLStblcLnK5XCGjAQATVEFXRkpLS2PBggXR2tqa3zY4OBitra1RV1d3wv3uu+++uPfee2Pnzp2xcOHCkU8LAEw6BV0ZiYhoaGiIlStXxsKFC2PRokWxcePG6O3tjVWrVkVExIoVK2LOnDnR3NwcERF/93d/F+vWrYvHH388ampq8q8t+eQnPxmf/OQnR/GhAAATUcExsmzZsjh8+HCsW7cuOjs7Y968ebFz5878i1oPHjwYxcUfXnD57ne/G/39/fFnf/ZnQ47T1NQU3/jGNz7e9ADAhFfw3xlJwd8ZAYCxM6H+zggAwGgTIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAICkxAgAkJUYAgKTECACQlBgBAJISIwBAUmIEAEhKjAAASYkRACApMQIAJCVGAICkxAgAkJQYAQCSEiMAQFJiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASGpEMbJp06aoqamJsrKyqK2tjd27d590/Y9+9KO44IILoqysLC6++OLYsWPHiIYFACafgmNk27Zt0dDQEE1NTbFnz56YO3duLF68ON56661h17/wwgtx/fXXx4033hgvvfRSXHvttXHttdfGf/zHf3zs4QGAia8oy7KskB1qa2vj85//fDz00EMRETE4OBjV1dXx9a9/PdauXXvc+mXLlkVvb2/89Kc/zW/74z/+45g3b15s3rz5lO6zp6cnKioqoru7O8rLywsZ9yPVrN0+qscDgInmwIYlY3LcU/39PaWQg/b390d7e3s0NjbmtxUXF0d9fX20tbUNu09bW1s0NDQM2bZ48eJ4+umnT3g/fX190dfXl/+6u7s7It5/UKNtsO/dUT8mAEwkY/H79f8f96OuexQUI0eOHImBgYGorKwcsr2ysjJeffXVYffp7Owcdn1nZ+cJ76e5uTnuueee47ZXV1cXMi4AcAoqNo7t8Y8ePRoVFRUn/H5BMTJeGhsbh1xNGRwcjLfffjt+53d+J4qKikbtfnp6eqK6ujoOHTo06k//UDjn4/TjnJxenI/Ti/Px0bIsi6NHj8bs2bNPuq6gGJkxY0aUlJREV1fXkO1dXV1RVVU17D5VVVUFrY+IyOVykcvlhmybNm1aIaMWpLy83P9IpxHn4/TjnJxenI/Ti/Nxcie7IvKBgt5NU1paGgsWLIjW1tb8tsHBwWhtbY26urph96mrqxuyPiLi2WefPeF6AODMUvDTNA0NDbFy5cpYuHBhLFq0KDZu3Bi9vb2xatWqiIhYsWJFzJkzJ5qbmyMi4rbbbosvfvGL8cADD8SSJUuipaUlXnzxxXjkkUdG95EAABNSwTGybNmyOHz4cKxbty46Oztj3rx5sXPnzvyLVA8ePBjFxR9ecLn00kvj8ccfj7vuuivuuOOO+IM/+IN4+umn46KLLhq9RzFCuVwumpqajntKiDScj9OPc3J6cT5OL87H6Cn474wAAIwmn00DACQlRgCApMQIAJCUGAEAkpr0MbJp06aoqamJsrKyqK2tjd27d590/Y9+9KO44IILoqysLC6++OLYsWPHOE16ZijkfGzZsiUuv/zymD59ekyfPj3q6+s/8vxRuEJ/Rj7Q0tISRUVFce21147tgGeYQs/Hr3/961i9enXMmjUrcrlcnH/++f7dGkWFno+NGzfGZz/72Tj77LOjuro61qxZE7/5zW/GadoJLJvEWlpastLS0mzr1q3Zf/7nf2Y333xzNm3atKyrq2vY9T//+c+zkpKS7L777st+8YtfZHfddVd21llnZS+//PI4Tz45FXo+brjhhmzTpk3ZSy+9lL3yyivZn//5n2cVFRXZf//3f4/z5JNXoefkA/v378/mzJmTXX755dmf/umfjs+wZ4BCz0dfX1+2cOHC7Jprrsmef/75bP/+/dmuXbuyjo6OcZ58cir0fPzgBz/Icrlc9oMf/CDbv39/9swzz2SzZs3K1qxZM86TTzyTOkYWLVqUrV69Ov/1wMBANnv27Ky5uXnY9dddd122ZMmSIdtqa2uzv/zLvxzTOc8UhZ6P33bs2LFs6tSp2fe///2xGvGMM5JzcuzYsezSSy/Nvve972UrV64UI6Oo0PPx3e9+Nzv33HOz/v7+8RrxjFLo+Vi9enX2J3/yJ0O2NTQ0ZJdddtmYzjkZTNqnafr7+6O9vT3q6+vz24qLi6O+vj7a2tqG3aetrW3I+oiIxYsXn3A9p24k5+O3vfvuu/Hee+/FOeecM1ZjnlFGek6++c1vxsyZM+PGG28cjzHPGCM5Hz/5yU+irq4uVq9eHZWVlXHRRRfF+vXrY2BgYLzGnrRGcj4uvfTSaG9vzz+Vs2/fvtixY0dcc8014zLzRHZafmrvaDhy5EgMDAzk/zLsByorK+PVV18ddp/Ozs5h13d2do7ZnGeKkZyP33b77bfH7NmzjwtGRmYk5+T555+PRx99NDo6OsZhwjPLSM7Hvn374l//9V/jK1/5SuzYsSP27t0bX/va1+K9996Lpqam8Rh70hrJ+bjhhhviyJEj8YUvfCGyLItjx47FrbfeGnfcccd4jDyhTdorI0wuGzZsiJaWlvjxj38cZWVlqcc5Ix09ejSWL18eW7ZsiRkzZqQeh3j/g0pnzpwZjzzySCxYsCCWLVsWd955Z2zevDn1aGekXbt2xfr16+Phhx+OPXv2xFNPPRXbt2+Pe++9N/Vop71Je2VkxowZUVJSEl1dXUO2d3V1RVVV1bD7VFVVFbSeUzeS8/GB+++/PzZs2BA/+9nP4pJLLhnLMc8ohZ6TX/7yl3HgwIFYunRpftvg4GBEREyZMiVee+21OO+888Z26ElsJD8js2bNirPOOitKSkry2y688MLo7OyM/v7+KC0tHdOZJ7ORnI+77747li9fHjfddFNERFx88cXR29sbt9xyS9x5551DPreNoSbtf5nS0tJYsGBBtLa25rcNDg5Ga2tr1NXVDbtPXV3dkPUREc8+++wJ13PqRnI+IiLuu+++uPfee2Pnzp2xcOHC8Rj1jFHoObngggvi5Zdfjo6Ojvzty1/+clx55ZXR0dER1dXV4zn+pDOSn5HLLrss9u7dm4/CiIjXX389Zs2aJUQ+ppGcj3ffffe44PggFDMfA3dyqV9BO5ZaWlqyXC6XPfbYY9kvfvGL7JZbbsmmTZuWdXZ2ZlmWZcuXL8/Wrl2bX//zn/88mzJlSnb//fdnr7zyStbU1OStvaOo0POxYcOGrLS0NHvyySezN998M387evRoqocw6RR6Tn6bd9OMrkLPx8GDB7OpU6dmf/VXf5W99tpr2U9/+tNs5syZ2be+9a1UD2FSKfR8NDU1ZVOnTs1++MMfZvv27cv+5V/+JTvvvPOy6667LtVDmDAmdYxkWZb9wz/8Q/bpT386Ky0tzRYtWpT9+7//e/57X/ziF7OVK1cOWf/EE09k559/flZaWpr94R/+YbZ9+/ZxnnhyK+R8/N7v/V4WEcfdmpqaxn/wSazQn5H/T4yMvkLPxwsvvJDV1tZmuVwuO/fcc7Nvf/vb2bFjx8Z56smrkPPx3nvvZd/4xjey8847LysrK8uqq6uzr33ta9n//M//jP/gE0xRlrl2BACkM2lfMwIATAxiBABISowAAEmJEQAgKTECACQlRgCApMQIAJCUGAEAkhIjAEBSYgQASEqMAABJiREAIKn/A1ga3zEAiiRxAAAAAElFTkSuQmCC\n"
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "<Figure size 640x480 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGdCAYAAADAAnMpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAftElEQVR4nO3df2xV9f3H8Vd/0FuM9AKrvaXsuiobogMLK9IVJIq72iDpRpbFRgztOsGh1ThultnKj4oo7QywmllsQBmaiUWJMmcbEOsIQboQCk3YBBwCa4fcC42zt5TZQu/5/mG8fist9NTefri3z0dy/+jp55z7vjtjfe7cXzGWZVkCAAAwJNb0AAAAYGgjRgAAgFHECAAAMIoYAQAARhEjAADAKGIEAAAYRYwAAACjiBEAAGBUvOkB+iIYDOrTTz/ViBEjFBMTY3ocAADQB5Zlqa2tTWlpaYqN7f36R0TEyKeffiq32216DAAA0A/Nzc367ne/2+vvIyJGRowYIenLB5OUlGR4GgAA0BeBQEButzv0d7w3EREjXz01k5SURIwAABBhrvQSC17ACgAAjCJGAACAUcQIAAAwihgBAABGESMAAMAoYgQAABhFjAAAAKOIEQAAYBQxAgAAjCJGAACAUcQIAAAwihgBAABGESMAAMAo2zGye/du5ebmKi0tTTExMdq2bVuf9/3www8VHx+vyZMn271bAAAQpeLt7tDe3q6MjAz96le/0s9//vM+7/f5558rPz9fP/nJT+T3++3eLYaI9OIa0yMMmJPlc0yPAAARwXaMzJ49W7Nnz7Z9R4sWLdK8efMUFxdn62oKAACIboPympE//elPOn78uEpLS/u0vqOjQ4FAoNsNAABEp7DHyL/+9S8VFxfrz3/+s+Lj+3YhpqysTE6nM3Rzu91hnhIAAJgS1hjp6urSvHnztGLFCo0fP77P+5WUlKi1tTV0a25uDuOUAADAJNuvGbGjra1N+/fv18GDB/Xoo49KkoLBoCzLUnx8vN577z3dddddl+zncDjkcDjCORoAALhKhDVGkpKSdOjQoW7b1q1bpw8++EBbt27VDTfcEM67BwAAEcB2jJw7d07Hjh0L/XzixAk1NjZq9OjRuv7661VSUqJTp07p1VdfVWxsrCZOnNht/5SUFCUmJl6yHQAADE22Y2T//v2aNWtW6Gev1ytJKigo0KZNm3T69Gk1NTUN3IQAACCqxViWZZke4koCgYCcTqdaW1uVlJRkehyEER96BgDRo69/v/luGgAAYBQxAgAAjCJGAACAUcQIAAAwihgBAABGESMAAMAoYgQAABhFjAAAAKOIEQAAYBQxAgAAjCJGAACAUcQIAAAwihgBAABGESMAAMAoYgQAABhFjAAAAKOIEQAAYBQxAgAAjCJGAACAUcQIAAAwihgBAABGESMAAMAoYgQAABhFjAAAAKOIEQAAYBQxAgAAjIo3PQAADIb04hrTIwyIk+VzTI8ADDiujAAAAKOIEQAAYBQxAgAAjCJGAACAUcQIAAAwihgBAABGESMAAMAoYgQAABhFjAAAAKOIEQAAYBQxAgAAjCJGAACAUbZjZPfu3crNzVVaWppiYmK0bdu2y65/6623dPfdd+u6665TUlKSsrOztWPHjv7OCwAAooztGGlvb1dGRoYqKyv7tH737t26++67VVtbq4aGBs2aNUu5ubk6ePCg7WEBAED0ibe7w+zZszV79uw+r6+oqOj286pVq/SXv/xFf/3rXzVlyhS7dw8AAKKM7Rj5toLBoNra2jR69Ohe13R0dKijoyP0cyAQGIzRAACAAYP+AtbVq1fr3Llzuu+++3pdU1ZWJqfTGbq53e5BnBAAAAymQY2RzZs3a8WKFXrjjTeUkpLS67qSkhK1traGbs3NzYM4JQAAGEyD9jRNdXW1FixYoDfffFMej+eyax0OhxwOxyBNBgAATBqUKyOvv/66CgsL9frrr2vOnDmDcZcAACBC2L4ycu7cOR07diz084kTJ9TY2KjRo0fr+uuvV0lJiU6dOqVXX31V0pdPzRQUFOj5559XVlaWfD6fJGn48OFyOp0D9DAAAECksh0j+/fv16xZs0I/e71eSVJBQYE2bdqk06dPq6mpKfT79evX6+LFiyoqKlJRUVFo+1frAQCIVOnFNaZHGBAny80+a2E7Ru68805ZltXr778ZGLt27bJ7FwAAYAjhu2kAAIBRxAgAADCKGAEAAEYRIwAAwChiBAAAGEWMAAAAo4gRAABgFDECAACMIkYAAIBRxAgAADCKGAEAAEYRIwAAwChiBAAAGEWMAAAAo4gRAABgVLzpAQAAfZdeXGN6hAFxsnyO6RFwFeHKCAAAMIorI1EiWv7fEgBg6OHKCAAAMIoYAQAARhEjAADAKGIEAAAYRYwAAACjiBEAAGAUMQIAAIwiRgAAgFHECAAAMIoYAQAARhEjAADAKGIEAAAYRYwAAACj+NZeIEyi5ZuUT5bPMT0CgCjHlREAAGAUMQIAAIwiRgAAgFHECAAAMIoYAQAARg35d9NEyzseAACIVFwZAQAARtmOkd27dys3N1dpaWmKiYnRtm3brrjPrl279KMf/UgOh0Pf//73tWnTpn6MCgAAopHtGGlvb1dGRoYqKyv7tP7EiROaM2eOZs2apcbGRv3mN7/RggULtGPHDtvDAgCA6GP7NSOzZ8/W7Nmz+7y+qqpKN9xwg9asWSNJuvnmm7Vnzx794Q9/UE5Ojt27BwAAUSbsrxmpr6+Xx+Ppti0nJ0f19fW97tPR0aFAINDtBgAAolPYY8Tn88nlcnXb5nK5FAgE9L///a/HfcrKyuR0OkM3t9sd7jEBAIAhV+W7aUpKStTa2hq6NTc3mx4JAACESdg/ZyQ1NVV+v7/bNr/fr6SkJA0fPrzHfRwOhxwOR7hHAwAAV4GwXxnJzs5WXV1dt207d+5UdnZ2uO8aAABEANsxcu7cOTU2NqqxsVHSl2/dbWxsVFNTk6Qvn2LJz88PrV+0aJGOHz+u3/3udzpy5IjWrVunN954Q4sXLx6YRwAAACKa7RjZv3+/pkyZoilTpkiSvF6vpkyZouXLl0uSTp8+HQoTSbrhhhtUU1OjnTt3KiMjQ2vWrNFLL73E23oBAICkfrxm5M4775RlWb3+vqdPV73zzjt18OBBu3cFAACGgKvy3TQAAGDoIEYAAIBRxAgAADCKGAEAAEYRIwAAwChiBAAAGEWMAAAAo4gRAABgFDECAACMCvu39gIA8E3pxTWmR8BVhCsjAADAKGIEAAAYRYwAAACjiBEAAGAUMQIAAIwiRgAAgFHECAAAMIrPGQFwWXweBIBw48oIAAAwihgBAABGESMAAMAoYgQAABhFjAAAAKOIEQAAYBQxAgAAjCJGAACAUcQIAAAwihgBAABGESMAAMAoYgQAABhFjAAAAKOIEQAAYBQxAgAAjCJGAACAUcQIAAAwihgBAABGESMAAMAoYgQAABhFjAAAAKP6FSOVlZVKT09XYmKisrKytG/fvsuur6io0E033aThw4fL7XZr8eLF+uKLL/o1MAAAiC62Y2TLli3yer0qLS3VgQMHlJGRoZycHJ05c6bH9Zs3b1ZxcbFKS0t1+PBhvfzyy9qyZYuefPLJbz08AACIfLZjZO3atVq4cKEKCwt1yy23qKqqStdcc402btzY4/q9e/dqxowZmjdvntLT03XPPffo/vvvv+LVFAAAMDTYipHOzk41NDTI4/F8fYDYWHk8HtXX1/e4z/Tp09XQ0BCKj+PHj6u2tlb33ntvr/fT0dGhQCDQ7QYAAKJTvJ3FLS0t6urqksvl6rbd5XLpyJEjPe4zb948tbS06Pbbb5dlWbp48aIWLVp02adpysrKtGLFCjujAQCACBX2d9Ps2rVLq1at0rp163TgwAG99dZbqqmp0cqVK3vdp6SkRK2traFbc3NzuMcEAACG2LoykpycrLi4OPn9/m7b/X6/UlNTe9xn2bJlmj9/vhYsWCBJmjRpktrb2/XQQw9pyZIlio29tIccDoccDoed0QAAQISydWUkISFBmZmZqqurC20LBoOqq6tTdnZ2j/ucP3/+kuCIi4uTJFmWZXdeAAAQZWxdGZEkr9ergoICTZ06VdOmTVNFRYXa29tVWFgoScrPz9fYsWNVVlYmScrNzdXatWs1ZcoUZWVl6dixY1q2bJlyc3NDUQIAAIYu2zGSl5ens2fPavny5fL5fJo8ebK2b98eelFrU1NTtyshS5cuVUxMjJYuXapTp07puuuuU25urp599tmBexQAACBixVgR8FxJIBCQ0+lUa2urkpKSBvTY6cU1A3o8AAAizcnyOWE5bl//fvPdNAAAwChiBAAAGEWMAAAAo4gRAABgFDECAACMIkYAAIBRxAgAADCKGAEAAEYRIwAAwChiBAAAGEWMAAAAo4gRAABgFDECAACMIkYAAIBRxAgAADCKGAEAAEYRIwAAwChiBAAAGEWMAAAAo4gRAABgFDECAACMIkYAAIBRxAgAADCKGAEAAEYRIwAAwChiBAAAGEWMAAAAo4gRAABgFDECAACMIkYAAIBRxAgAADCKGAEAAEYRIwAAwChiBAAAGEWMAAAAo4gRAABgFDECAACMIkYAAIBRxAgAADCqXzFSWVmp9PR0JSYmKisrS/v27bvs+s8//1xFRUUaM2aMHA6Hxo8fr9ra2n4NDAAAoku83R22bNkir9erqqoqZWVlqaKiQjk5OTp69KhSUlIuWd/Z2am7775bKSkp2rp1q8aOHat///vfGjly5EDMDwAAIpztGFm7dq0WLlyowsJCSVJVVZVqamq0ceNGFRcXX7J+48aN+uyzz7R3714NGzZMkpSenv7tpgYAAFHD1tM0nZ2damhokMfj+foAsbHyeDyqr6/vcZ933nlH2dnZKioqksvl0sSJE7Vq1Sp1dXX1ej8dHR0KBALdbgAAIDrZipGWlhZ1dXXJ5XJ12+5yueTz+Xrc5/jx49q6dau6urpUW1urZcuWac2aNXrmmWd6vZ+ysjI5nc7Qze122xkTAABEkLC/myYYDColJUXr169XZmam8vLytGTJElVVVfW6T0lJiVpbW0O35ubmcI8JAAAMsfWakeTkZMXFxcnv93fb7vf7lZqa2uM+Y8aM0bBhwxQXFxfadvPNN8vn86mzs1MJCQmX7ONwOORwOOyMBgAAIpStKyMJCQnKzMxUXV1daFswGFRdXZ2ys7N73GfGjBk6duyYgsFgaNvHH3+sMWPG9BgiAABgaLH9NI3X69WGDRv0yiuv6PDhw3r44YfV3t4eendNfn6+SkpKQusffvhhffbZZ3r88cf18ccfq6amRqtWrVJRUdHAPQoAABCxbL+1Ny8vT2fPntXy5cvl8/k0efJkbd++PfSi1qamJsXGft04brdbO3bs0OLFi3Xrrbdq7Nixevzxx/XEE08M3KMAAAARK8ayLMv0EFcSCATkdDrV2tqqpKSkAT12enHNgB4PAIBIc7J8TliO29e/33w3DQAAMIoYAQAARhEjAADAKGIEAAAYRYwAAACjiBEAAGAUMQIAAIwiRgAAgFHECAAAMIoYAQAARhEjAADAKGIEAAAYRYwAAACjiBEAAGAUMQIAAIwiRgAAgFHECAAAMIoYAQAARhEjAADAKGIEAAAYRYwAAACjiBEAAGAUMQIAAIwiRgAAgFHECAAAMIoYAQAARhEjAADAKGIEAAAYRYwAAACjiBEAAGAUMQIAAIwiRgAAgFHECAAAMIoYAQAARhEjAADAKGIEAAAYRYwAAACjiBEAAGBUv2KksrJS6enpSkxMVFZWlvbt29en/aqrqxUTE6O5c+f2524BAEAUsh0jW7ZskdfrVWlpqQ4cOKCMjAzl5OTozJkzl93v5MmT+u1vf6uZM2f2e1gAABB9bMfI2rVrtXDhQhUWFuqWW25RVVWVrrnmGm3cuLHXfbq6uvTAAw9oxYoVuvHGG7/VwAAAILrYipHOzk41NDTI4/F8fYDYWHk8HtXX1/e639NPP62UlBQ9+OCD/Z8UAABEpXg7i1taWtTV1SWXy9Vtu8vl0pEjR3rcZ8+ePXr55ZfV2NjY5/vp6OhQR0dH6OdAIGBnTAAAEEHC+m6atrY2zZ8/Xxs2bFBycnKf9ysrK5PT6Qzd3G53GKcEAAAm2boykpycrLi4OPn9/m7b/X6/UlNTL1n/ySef6OTJk8rNzQ1tCwaDX95xfLyOHj2qcePGXbJfSUmJvF5v6OdAIECQAAAQpWzFSEJCgjIzM1VXVxd6e24wGFRdXZ0effTRS9ZPmDBBhw4d6rZt6dKlamtr0/PPP99rYDgcDjkcDjujAQCACGUrRiTJ6/WqoKBAU6dO1bRp01RRUaH29nYVFhZKkvLz8zV27FiVlZUpMTFREydO7Lb/yJEjJemS7QAAYGiyHSN5eXk6e/asli9fLp/Pp8mTJ2v79u2hF7U2NTUpNpYPdgUAAH0TY1mWZXqIKwkEAnI6nWptbVVSUtKAHju9uGZAjwcAQKQ5WT4nLMft699vLmEAAACjiBEAAGAUMQIAAIwiRgAAgFHECAAAMIoYAQAARhEjAADAKGIEAAAYRYwAAACjiBEAAGAUMQIAAIwiRgAAgFHECAAAMIoYAQAARhEjAADAKGIEAAAYRYwAAACjiBEAAGAUMQIAAIwiRgAAgFHECAAAMIoYAQAARhEjAADAKGIEAAAYRYwAAACjiBEAAGAUMQIAAIwiRgAAgFHECAAAMIoYAQAARhEjAADAKGIEAAAYRYwAAACjiBEAAGAUMQIAAIwiRgAAgFHECAAAMIoYAQAARhEjAADAqH7FSGVlpdLT05WYmKisrCzt27ev17UbNmzQzJkzNWrUKI0aNUoej+ey6wEAwNBiO0a2bNkir9er0tJSHThwQBkZGcrJydGZM2d6XL9r1y7df//9+tvf/qb6+nq53W7dc889OnXq1LceHgAARL4Yy7IsOztkZWXptttu0wsvvCBJCgaDcrvdeuyxx1RcXHzF/bu6ujRq1Ci98MILys/P79N9BgIBOZ1Otba2Kikpyc64V5ReXDOgxwMAINKcLJ8TluP29e+3rSsjnZ2damhokMfj+foAsbHyeDyqr6/v0zHOnz+vCxcuaPTo0b2u6ejoUCAQ6HYDAADRyVaMtLS0qKurSy6Xq9t2l8sln8/Xp2M88cQTSktL6xY031RWVian0xm6ud1uO2MCAIAIMqjvpikvL1d1dbXefvttJSYm9rqupKREra2toVtzc/MgTgkAAAZTvJ3FycnJiouLk9/v77bd7/crNTX1svuuXr1a5eXlev/993Xrrbdedq3D4ZDD4bAzGgAAiFC2rowkJCQoMzNTdXV1oW3BYFB1dXXKzs7udb/nnntOK1eu1Pbt2zV16tT+TwsAAKKOrSsjkuT1elVQUKCpU6dq2rRpqqioUHt7uwoLCyVJ+fn5Gjt2rMrKyiRJv//977V8+XJt3rxZ6enpodeWXHvttbr22msH8KEAAIBIZDtG8vLydPbsWS1fvlw+n0+TJ0/W9u3bQy9qbWpqUmzs1xdcXnzxRXV2duoXv/hFt+OUlpbqqaee+nbTAwCAiGf7c0ZM4HNGAAAIn4j6nBEAAICBRowAAACjiBEAAGAUMQIAAIwiRgAAgFHECAAAMIoYAQAARhEjAADAKGIEAAAYRYwAAACjiBEAAGAUMQIAAIwiRgAAgFHECAAAMIoYAQAARhEjAADAKGIEAAAYRYwAAACjiBEAAGAUMQIAAIwiRgAAgFHECAAAMIoYAQAARhEjAADAKGIEAAAYRYwAAACjiBEAAGAUMQIAAIwiRgAAgFHECAAAMIoYAQAARhEjAADAKGIEAAAYRYwAAACjiBEAAGAUMQIAAIwiRgAAgFHECAAAMIoYAQAARvUrRiorK5Wenq7ExERlZWVp3759l13/5ptvasKECUpMTNSkSZNUW1vbr2EBAED0sR0jW7ZskdfrVWlpqQ4cOKCMjAzl5OTozJkzPa7fu3ev7r//fj344IM6ePCg5s6dq7lz5+of//jHtx4eAABEvhjLsiw7O2RlZem2227TCy+8IEkKBoNyu9167LHHVFxcfMn6vLw8tbe369133w1t+/GPf6zJkyerqqqqT/cZCATkdDrV2tqqpKQkO+NeUXpxzYAeDwCASHOyfE5YjtvXv9/xdg7a2dmphoYGlZSUhLbFxsbK4/Govr6+x33q6+vl9Xq7bcvJydG2bdt6vZ+Ojg51dHSEfm5tbZX05YMaaMGO8wN+TAAAIkk4/r7+/+Ne6bqHrRhpaWlRV1eXXC5Xt+0ul0tHjhzpcR+fz9fjep/P1+v9lJWVacWKFZdsd7vddsYFAAB94KwI7/Hb2trkdDp7/b2tGBksJSUl3a6mBINBffbZZ/rOd76jmJiYAbufQCAgt9ut5ubmAX/6B/ZxPq4+nJOrC+fj6sL5uDLLstTW1qa0tLTLrrMVI8nJyYqLi5Pf7++23e/3KzU1tcd9UlNTba2XJIfDIYfD0W3byJEj7YxqS1JSEv9FuopwPq4+nJOrC+fj6sL5uLzLXRH5iq130yQkJCgzM1N1dXWhbcFgUHV1dcrOzu5xn+zs7G7rJWnnzp29rgcAAEOL7adpvF6vCgoKNHXqVE2bNk0VFRVqb29XYWGhJCk/P19jx45VWVmZJOnxxx/XHXfcoTVr1mjOnDmqrq7W/v37tX79+oF9JAAAICLZjpG8vDydPXtWy5cvl8/n0+TJk7V9+/bQi1SbmpoUG/v1BZfp06dr8+bNWrp0qZ588kn94Ac/0LZt2zRx4sSBexT95HA4VFpaeslTQjCD83H14ZxcXTgfVxfOx8Cx/TkjAAAAA4nvpgEAAEYRIwAAwChiBAAAGEWMAAAAo6I+RiorK5Wenq7ExERlZWVp3759l13/5ptvasKECUpMTNSkSZNUW1s7SJMODXbOx4YNGzRz5kyNGjVKo0aNksfjueL5g312/418pbq6WjExMZo7d254Bxxi7J6Pzz//XEVFRRozZowcDofGjx/P/24NILvno6KiQjfddJOGDx8ut9utxYsX64svvhikaSOYFcWqq6uthIQEa+PGjdY///lPa+HChdbIkSMtv9/f4/oPP/zQiouLs5577jnro48+spYuXWoNGzbMOnTo0CBPHp3sno958+ZZlZWV1sGDB63Dhw9bv/zlLy2n02n95z//GeTJo5fdc/KVEydOWGPHjrVmzpxp/exnPxucYYcAu+ejo6PDmjp1qnXvvfdae/bssU6cOGHt2rXLamxsHOTJo5Pd8/Haa69ZDofDeu2116wTJ05YO3bssMaMGWMtXrx4kCePPFEdI9OmTbOKiopCP3d1dVlpaWlWWVlZj+vvu+8+a86cOd22ZWVlWb/+9a/DOudQYfd8fNPFixetESNGWK+88kq4Rhxy+nNOLl68aE2fPt166aWXrIKCAmJkANk9Hy+++KJ14403Wp2dnYM14pBi93wUFRVZd911V7dtXq/XmjFjRljnjAZR+zRNZ2enGhoa5PF4QttiY2Pl8XhUX1/f4z719fXd1ktSTk5Or+vRd/05H990/vx5XbhwQaNHjw7XmENKf8/J008/rZSUFD344IODMeaQ0Z/z8c477yg7O1tFRUVyuVyaOHGiVq1apa6ursEaO2r153xMnz5dDQ0Noadyjh8/rtraWt17772DMnMkuyq/tXcgtLS0qKurK/TJsF9xuVw6cuRIj/v4fL4e1/t8vrDNOVT053x80xNPPKG0tLRLghH9059zsmfPHr388stqbGwchAmHlv6cj+PHj+uDDz7QAw88oNraWh07dkyPPPKILly4oNLS0sEYO2r153zMmzdPLS0tuv3222VZli5evKhFixbpySefHIyRI1rUXhlBdCkvL1d1dbXefvttJSYmmh5nSGpra9P8+fO1YcMGJScnmx4H+vKLSlNSUrR+/XplZmYqLy9PS5YsUVVVlenRhqRdu3Zp1apVWrdunQ4cOKC33npLNTU1WrlypenRrnpRe2UkOTlZcXFx8vv93bb7/X6lpqb2uE9qaqqt9ei7/pyPr6xevVrl5eV6//33deutt4ZzzCHF7jn55JNPdPLkSeXm5oa2BYNBSVJ8fLyOHj2qcePGhXfoKNaffyNjxozRsGHDFBcXF9p28803y+fzqbOzUwkJCWGdOZr153wsW7ZM8+fP14IFCyRJkyZNUnt7ux566CEtWbKk2/e2obuo/U8mISFBmZmZqqurC20LBoOqq6tTdnZ2j/tkZ2d3Wy9JO3fu7HU9+q4/50OSnnvuOa1cuVLbt2/X1KlTB2PUIcPuOZkwYYIOHTqkxsbG0O2nP/2pZs2apcbGRrnd7sEcP+r059/IjBkzdOzYsVAUStLHH3+sMWPGECLfUn/Ox/nz5y8Jjq9C0eJr4C7P9Ctow6m6utpyOBzWpk2brI8++sh66KGHrJEjR1o+n8+yLMuaP3++VVxcHFr/4YcfWvHx8dbq1autw4cPW6Wlpby1dwDZPR/l5eVWQkKCtXXrVuv06dOhW1tbm6mHEHXsnpNv4t00A8vu+WhqarJGjBhhPfroo9bRo0etd99910pJSbGeeeYZUw8hqtg9H6WlpdaIESOs119/3Tp+/Lj13nvvWePGjbPuu+8+Uw8hYkR1jFiWZf3xj3+0rr/+eishIcGaNm2a9fe//z30uzvuuMMqKCjotv6NN96wxo8fbyUkJFg//OEPrZqamkGeOLrZOR/f+973LEmX3EpLSwd/8Chm99/I/0eMDDy752Pv3r1WVlaW5XA4rBtvvNF69tlnrYsXLw7y1NHLzvm4cOGC9dRTT1njxo2zEhMTLbfbbT3yyCPWf//738EfPMLEWBbXjgAAgDlR+5oRAAAQGYgRAABgFDECAACMIkYAAIBRxAgAADCKGAEAAEYRIwAAwChiBAAAGEWMAAAAo4gRAABgFDECAACMIkYAAIBR/wd+WHhAS/8R3gAAAABJRU5ErkJggg==\n"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Define a simple model with two latent variables and two observation variables, all scalar\n",
    "\n",
    "def simulation_rank(latent_vars,observed_vars):\n",
    "    \"\"\"\n",
    "    1. Draw (θٰʹ,yʹ) ~ p(θ,y)\n",
    "    2. Use MCMC to generate samples approximately from p(θ|yʹ)\n",
    "    3. Find rank of θʹ among samples from θ\n",
    "    \"\"\"\n",
    "    # draw latents and observed variables from prior\n",
    "    all_samps = sample(latent_vars + observed_vars, niter=1) # just draw one sample\n",
    "    all_samps = [samps[0] for samps in all_samps] # remove dummy first dimension\n",
    "    latent_samps = all_samps[:len(latent_vars)]\n",
    "    observed_samps = all_samps[len(latent_vars):]\n",
    "\n",
    "    # do inference on latents given observed variables\n",
    "    # intentionally screw up calibration by using a small number of Gibbs iterations\n",
    "    evidence = [observed_var == observed_samp for (observed_var, observed_samp) in zip(observed_vars, observed_samps)]\n",
    "    inference_samps = sample(latent_vars, Given(*evidence), niter=100)\n",
    "\n",
    "    # get rank of true latents compared to sampled latents\n",
    "    ranks = []\n",
    "    for latent_samp, inference_samp in zip(latent_samps, inference_samps):\n",
    "        ranks.append( np.mean(latent_samp < inference_samp, axis=0))\n",
    "\n",
    "    return ranks\n",
    "\n",
    "def SBC(prior,likelihood,nreps=1000):\n",
    "    \"\"\"\n",
    "    Call simulation_rank several times, make a histogram\n",
    "    \"\"\"\n",
    "    my_ranks = simulation_rank(prior,likelihood)\n",
    "    ranks = [np.zeros([nreps, *r.shape]) for r in my_ranks]\n",
    "    for i in range(nreps):\n",
    "        my_ranks = simulation_rank(prior,likelihood)\n",
    "        for r, my_r in zip(ranks, my_ranks):\n",
    "            r[i] = my_r\n",
    "\n",
    "    for r in ranks:\n",
    "        plt.figure()\n",
    "        plt.hist(r,np.arange(0,1,.1),density=True)\n",
    "\n",
    "u = d.norm(0,1)\n",
    "s = d.exp(1)\n",
    "y = d.norm(u,s)\n",
    "z = d.norm(u,s)\n",
    "\n",
    "latent_vars = [u,s]\n",
    "observed_vars = [y,z]\n",
    "\n",
    "SBC(latent_vars,observed_vars)"
   ],
   "metadata": {
    "collapsed": false
   }
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}