CBayes - Nonlinear Stochastic Example.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 154,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Imports necessary Packages\n",
    "import matplotlib\n",
    "import numpy as np\n",
    "import scipy as sp\n",
    "\n",
    "# Specific Plotting Packages\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import cm\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "\n",
    "# Specific Stats Packages\n",
    "from scipy.stats import gaussian_kde\n",
    "from scipy.stats import truncexpon\n",
    "from scipy.stats import uniform\n",
    "from scipy.stats import norm\n",
    "from scipy.stats import chi2\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "# for Q map\n",
    "from scipy import optimize\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Thing Exampl\n",
    "\n",
    "Stuff\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define Map and Simulate the Data\n",
    "Here we've embedded the stochasticity inside the model map Q. In addition, we create a simulated data set by drawing samples from the \"true\" distribution of lambda 1 and lambda 2 and passing them through the model Q. Since the stochasticity is embedded in the model Q, the error $\\eta$ is not accessible (or saved).\n",
    "\n",
    "In addition, we don't have access to a prior distribution to sample the error term!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 227,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Defines the Q Map (including stochastic parameter)\n",
    "def Q(lam1,lam2):\n",
    "    ## input checks\n",
    "    lam1_m = lam1\n",
    "    lam2_m = lam2\n",
    "    \n",
    "    # checks lam1 is numeric and converts to np.ndarray if possible\n",
    "    if lam1!=np.ndarray:\n",
    "        try: lam1_m = np.asarray(lam1)\n",
    "        except ValueError: (\"Couldn't convert to nupmy array! Check that lam1 is a numeric array.\")        \n",
    "    \n",
    "    # checks lam2 is numeric and converts to np.ndarray if possible\n",
    "    if lam2_m!=np.ndarray:\n",
    "        try: lam2_m = np.asarray(lam2)\n",
    "        except ValueError: (\"Couldn't convert to nupmy array! Check that lam1 is a numeric array.\")\n",
    "    \n",
    "    \n",
    "    \n",
    "    # number of values to evaluate QoI\n",
    "    lam1_size = np.size(lam1_m)\n",
    "    lam2_size = np.size(lam2_m)\n",
    "    \n",
    "    ### RANDOM STOCHASTICISTY INTRODUCED\n",
    "    lam1_m = lam1_m + np.random.normal(0,0.5,lam1_size)**2\n",
    "    lam2_m = lam2_m + np.random.normal(0,0.5,lam2_size)**2\n",
    "    \n",
    "    \n",
    "    # Empty QoI vector\n",
    "    q = np.zeros(lam1_size)\n",
    "    \n",
    "    \n",
    "    ### Calculates the QoI for the given lam1 and lam2 values\n",
    "    for i in np.arange(lam1_size):\n",
    "        # sets up nonlinear system for the solver\n",
    "        def G(x):    \n",
    "            # create coefficient matrix of nonlinear system\n",
    "            A = np.reshape([lam1_m.item(i),1,1,-lam2_m.item(i)],(2,2))\n",
    "\n",
    "            # compute matrix-vector multiplication, then subtract 1\n",
    "            return A @ x**2 - 1\n",
    "\n",
    "        # finds the x-coordinates in Q1 that satisfy the model\n",
    "        x = optimize.fsolve(G,[1,1])\n",
    "        \n",
    "        # save QoI\n",
    "        q[i] = x[1] # quantity of interest is second coordinate \n",
    "        \n",
    "    return q"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generating the simulated data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 228,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "data_sample_size = 500"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 229,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\ProgramData\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:8: DeprecationWarning: elementwise != comparison failed; this will raise an error in the future.\n",
      "  \n",
      "C:\\ProgramData\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:13: DeprecationWarning: elementwise != comparison failed; this will raise an error in the future.\n",
      "  del sys.path[0]\n",
      "C:\\ProgramData\\Anaconda3\\lib\\site-packages\\scipy\\optimize\\minpack.py:161: RuntimeWarning: The iteration is not making good progress, as measured by the \n",
      "  improvement from the last ten iterations.\n",
      "  warnings.warn(msg, RuntimeWarning)\n",
      "C:\\ProgramData\\Anaconda3\\lib\\site-packages\\scipy\\optimize\\minpack.py:161: RuntimeWarning: The iteration is not making good progress, as measured by the \n",
      "  improvement from the last five Jacobian evaluations.\n",
      "  warnings.warn(msg, RuntimeWarning)\n"
     ]
    }
   ],
   "source": [
    "# generate data\n",
    "# simulate data by sampling uniforms \n",
    "lam1_data = uniform.rvs(0.5,0.5,data_sample_size)\n",
    "lam2_data = uniform.rvs(0.75,0.5,data_sample_size)\n",
    "\n",
    "# simulated QoI data from data variable samples\n",
    "q_data = Q(lam1_data,lam2_data)\n",
    "\n",
    "# calculate a gaussian kde for data\n",
    "data_kde_pdf = gaussian_kde(q_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Reasonability check:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 230,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x227126a46a0>]"
      ]
     },
     "execution_count": 230,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAzYAAAEICAYAAACNoGG1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3Xe8HFX9//HX+97chJbQEiGkXXoH\nAxGCINKlCBFECSoIlihfsJcfogKKBRuKIipNBBVRQAyE3lsoISRACIEAgYQECKEkoSS5uZ/fHzMX\nls0te3N3d3b2vp+Pxz7u7s7Zmc/Mzpy7nzlnzigiMDMzMzMzy7OGrAMwMzMzMzPrKSc2ZmZmZmaW\ne05szMzMzMws95zYmJmZmZlZ7jmxMTMzMzOz3HNiY2ZmZmZmuderEhtJ0yTtkXUc9UjSqZL+3o3y\nF0r6SSVjKljWIEkzJK1SjeX1VCnbRtIVkvavVkz2Xq5LKsd1Sfm4LikfH/Mg6VBJsyUtljQy63g6\nIunuSsUn6RBJ/6rEvDtYXrfqwyxJWlXSVZJel/SfrOKom8RG0ixJ+xS9d4yku9peR8TWEXFbF/Np\nlhSS+lQo1Koq3gb1QNIJkiZJWiLpwhI+ciLw14h4O/38bZL2SCuMU4vmvaGkVklndyOeilU8aZy3\npc8Lbzp1OvDTSiyzt3Nd0r56q0sk9ZN0vqRnJS2S9JCkA7r4mOuSOuRjvmS/Bk6IiDUi4qHiiem6\nv5EmPgsk3SzpiFJnnu6jc3oSoKSDgUWF8UnaStL49Af3Ikm3SBrdjXm+s39ExHhgG0nb9STOcije\nXpL6picq7pY0IK1PlqXrvEjSE5LOkjS4aB6t6XdW+NhlJUI6HFgPWDciPtFBzGX7LjpSN4lNXtRx\nhVdNc4GfABd0VVBSP+CzQKk/Fo4GXgXGpp+tSRFxPzBA0qisY7FsuC7psT7AbODDwJrAD4F/S2pu\nr7DrEstaDRzzI4BpXZTZPiLWADYHLgTOknRKpQMr8GXg4rYXkjYG7gYeATYENgCuBG6UtNNKLuMS\nYFwP4yyrtI65AlgL2C8iFqaTLo2I/sA6wKHA+sCDhckNMDdNVgsfE1cijBHAExHR0kGMlfguVtCr\nEpvCTE/STulZ/4WSXpR0RlrsjvTva21Zq6QGST9Iz+y9JOkiSWsWzPfodNoCST8sWs6pki6T9HdJ\nC4Fj0mVPlPSapHlpBt23YH4h6f8kPZlmtKdJ2jj9zEJJ/y4s34Ptcayk6ekynpb0pYJpe0iaI+m7\n6TrPk/QxSQemWf8rkk4qmuUqki5N5zdZ0vYF8xuZvrdI0qXAKgXT1pZ0taT5kl5Nnw/tKO6IuCIi\nrgQWlLCaOwOvRUSpZ4GOBn4ALAMOLpwgaWtJN6br/qKkk5R04TgJOCLdX6amZd9zVkFFZ2Il/UfS\nC+lZizskbV1ifIVuAw5aic9ZD7kuWWF75K4uiYg3IuLUiJgVEa0RcTXwDLBjB6vpuqQX6w3HfEex\nKmndXAw0AlMlPdXV9oqIlyPiYuA44HuS1k2X0W5dIWl14FpgA73barBBV+tbFH9fYC/g9oK3TwUm\nRsT3I+KViFgUEb8nOUHxi4LPHqKku+FrSlpit+xk9W6jk+NF0omSnkrX8TFJhxZMO0bSXZJ+ndZR\nz6igpVhJS+/t6WdvBAZ2EkfbZ1YDrgKagIMi4o3iMhGxLCKmAUcA84FvdTXfDpa1Zbp9Xku31yHp\n+z8CTubd+uvz7Xz8VMr/XawoIuriAcwC9il67xjgrvbKABOBo9LnawCj0+fNQAB9Cj73OWAmsFFa\n9grg4nTaVsBiYDegL0lT7bKC5Zyavv4YSSK5Ksk/ztEkZwybgenA1wuWF8B4YACwNbAEuDld/prA\nY8BnC8q/BuzWwXZ5zzYomnYQsDEgkrOWbwI7pNP2AFpIdtQm4IskB8M/gf5pXG8DGxWt5+Fp+W+T\n/EhoSrfLs8A30teHp2V/kn52XeDjwGrpvP8DXFnCd/4T4MIuyhwPTChxH/pQuq3XBv4AjC+Y1h+Y\nR1IZrJK+3rlg3f/e2f5YXCbdp/oD/YDfAVMKpl3Ytm26iPebwBVZH3v19ij+7tL33nMc4bqkeFqu\n65L0s+ulcWzRwXTXJXX6KN7G7e3v1OkxX7TOHcZaMO9NOtmOK0wnOU5bgAPS113VFXOKPt/p+haV\n3Rp4o+i9F4Bj2ym7ZxrXKsBmwBvAvmm83023Q98OjsF10nUd0EEcnyBpjWggSSTeAAYX7FfLSOrB\nRpLEby6ggn3rDJLjeXdgEUV1QsFy9iCpS29Pv/N+RdNPbe+zwI+B+zra5p18v03pdjmJZH/dK41v\n886WV8nvor1HvbXYXJlmeK9Jeg3orG/zMmATSQMjYnFE3NtJ2U8DZ0TE0xGxGPgeSfeCPiT/WK+K\niLsiYinJP+8o+vzEiLgykrOCb0XEgxFxb0S0RMQs4C8kB3ihX0TEwkgy7EeBG9Llv05yVuOdC+Mi\nYq2I6Hbf94iYEBFPReJ24AaSf8htlgE/jYhlwL9IzhycGUmWPY2kSbqwn+mDEXFZWv4Mkp10dPpo\nAn4XyVmDy4AHCuJYEBGXR8SbEbGIpL938fZYWWuRHHil+CxwbUS8SvKj6wBJ70unfRR4ISJ+ExFv\np9vgvpUNKiIuSOexhKQy2L7wLF6JFpGsn5Wf65JuyHtdIqkJ+Afwt4h4vINirkvqW6885rsR60pJ\nj+GXSZKBUuqK4s+Xsr5t2jtGB5KcSCg2jySxWIck+ZgQETem8f6aJIn8YAfLaVtGu8dMRPwnIuam\n39mlwJNAYVerZyPi3IhYDvwNGAysJ2k48AHghxGxJCLuIGmJ6Ux/YBeSumtJF2XbzCX9PlIbFO77\n6WP1dj43miThPT0ilkbELcDVwJElLrcS38UK6i2x+Vj6j3mtiFgL+L9Oyn6eJDN8XNIDkj7aSdkN\nSM4StnmW5OzBeum02W0TIuJNVuwiNbvwhaTNlHSReCFtXv4ZKzY3vljw/K12Xq/RSbwlkXSApHuV\ndId4DTiwKI4F6YHXtsz24iqMo3A7tAJzSLbPBsDzEVFYYb+zPSWtJukvafP3QpLm/LUkNfZwFSHp\n496/q0KSViU5y/KPNP6JwHPAp9Iiw4Aum99LIalR0ulpU/VCkjMQUEKTc5H+JGfYrfxcl3RDnusS\nSQ0kffKXAid0spquS+qbj/nOY10p6UmDQcAr6euu6oriz5eyvm3aO0ZfJkkcig0mSSJfpmi90zpn\nNjCkg+W0LaPdY0ZJ98IpBUnyNkUxv1CwrDfTp2ukcbwa7+1KVvh9tOdlYCzwN0kf6aJsmyGk30dq\nbuG+nz5W6M6Wxjc73T6F8XW0ndqLtdzfxQrqLbEpWUQ8GRFHAu8j6dt3WZqhFp8tgSS7HVHwejhJ\ns9mLJJnmO324039q6xYvruj1n4DHgU0jYgBJs55Wfm26T8mFZpeTZMPrpRX5NT2MY1jB/BtItstc\nkm00RFLhvIcXPP8WyYWGO6fbY/e22fQgljYPk/wD6sqhJE33Z6cV6AskB9LR6fTZJM3n7Wlvn3mD\npDtMm/ULnn8KGAPsQ9I1oDl9v7vruyUwtZufsTJzXZLfuiSdz/kkP9w+np4h7IjrEgPq+pjvLNaV\nNSadx/0l1BXtbb/urO+TJId14Y/gm0hONBT7JHBv2lL2nvVO64VhwPMdLGdLYFa8e4H+OySNAM4l\nOUmybrqOj3YSc6F5wNpFrSXDOyrcJiKuIOnadpmkPTsrm9anBwN3lhBPsbnAsHQehfF1tJ2KVeK7\nWEGvTWwkfUbSoDQbbMu6l5P0V2wl6WPa5hLgG0ou6lqD5IzBpZGM/HAZcLCkDyq5cO1HdL0D9wcW\nAoslbUHSx7KSJGmVwgdJ/8h+JOvbouTitf16uJwdJR2WNlt/naRv770kfUZbgK9K6iPpMN7bLNuf\n5CzSa5LWAU7pYmX6pOvQCDSm69RRU/n9JGdsu8r2P0syytq2wPvTx67A+yVtS9Lcur6kryu5kLK/\npJ3Tz74INBcd7FNImvCblIw2dHjR+i4hOTO3Gsn+tDI+TNKtwDLkuiS/dQnJj6YtgYMj4q0uyrou\nMaCuj/nOYu0WSetI+jTwR5KucQvouq54EVhX7+1KWfL6picmbuK9XdV+BHxQ0k/TmPpL+gpwLO/W\nD/8GDpK0d9rC9C2S4+qeDhbV2fHSluDOT7fDsSQtNl2KiGeBScCPlAzdvBtFA4908tlLSJKp/0na\ntXh6Wn9sSfIdr0/Sxbe77iM50fLddH57pPGVel+fSnwXK+i1iQ2wPzBNyUgfZwJj0/7Ob5L0y75b\nSTPiaJJ/UheTdGt4huQC068ApH1Yv0Lyxc4j6Xv5EskX0ZFvk5xpW0SS2V/akxVRMgJFh31USfom\nvtXO46skO9GraTzjexIH8D+S/pGvAkcBh0XSD34pcBjJRXOvpmWuKPjc70j6UL5M8uPlui6W84M0\n/hOBz6TPf9BewXTZF6bl2pX+UNmbpN/+CwWPB9NYPhtJf/19SQ7iF0jODLWdGWm7EdUCSZPT5z8k\nOSv7KsnB/M+CRV5E0tT6PMmFnJ31z+4o5g+QXCR5f3c/a2XnuiSHdUl6ZvVLJInHC3p3JKZPt1fe\ndYkVqJtjvkiHsXbD1HS7zAS+AHwjIk4GSPf9DuuKSK5vuwR4Ot1+G9D99f0LSZ3RNs8nSQZn2J6k\nq+ZrwGnAoRFxY1pmBslx/QeSuuNgkpMdSztYxpHpclYQEY8BvyE5CfMiyQmOu7uIudCnSEZgfIXk\nx/5FpX4wIv5GkghM0LvDJx+Rfh+vkWzrBcCOETG34KOFI9G1PT7ezvyXAocAB5Bsp7OBo6Pj6xKL\nP1+J72IFbaMwWJmkZzleI2k2fSbreAwkDSJpdh1ZwlnZXJB0OXB+RFyTdSxWGa5Lao/rEqskH/Pl\noeTGql+J9m8iOpTkBMApEXH+Ssz7YJIR8T7Z80h7t55+Fx3O14lNz6U7+s0kTci/Icm2dwhvXDPr\nBtclZr2Lj/nqS7uEjiFpWV2cdTy9WSW+i6zvYlsvxpA034qkf+RYV0pmthJcl5j1Lj7mqywiHgEe\nyToOq8x34RYbMzMzMzPLvd48eICZVZmkYZJulTRd0jRJX2unzB6SXldyH4Apkk7OIlYzMzPLl8y6\nog0cODCam5uzWryZtePBBx98OSIGVXARLcC3ImKypP7Ag5JuTEeSKXRnRHR207v3cH1iVnuqUJ9U\nhOsTs9rSnboks8SmubmZSZMmZbV4M2uHpK7uctwjETGPZFhTImKRpOkkNy8sTmy6xfWJWe2pdH1S\nKa5PzGpLd+oSd0Uzs0xIagZGktz0q9gukqZKulbS1h18fpykSZImzZ8/v4KRmpmZWR44sTGzqkvv\n13A58PWIWFg0eTIwIiK2J7lJ15XtzSMizomIURExatCg3PV2MTMzszJzYmNmVSWpiSSp+UdEXFE8\nPSIWto1nn940sEnSwCqHaWZmZjnjxMbMqkaSgPOB6RFxRgdl1k/LIWknknpqQfWiNDMzszzyDTrN\nrJp2BY4CHpE0JX3vJGA4QET8GTgcOE5SC/AWvmGdmZmZlcCJjZlVTUTcRXKH7c7KnAWcVZ2IzKye\nSdocuLTgrY2AkyPidxmFZGYV5MTGzMzM6lJEzADeDyCpEXge+G+mQZlZxfgaGzMzM+sN9gaeiohc\n3l/HzLrmxMbMzMx6g7HAJe1N8H2xzOqDu6KZ1ZjmEyeUfZ6zTj+o7PO0ysnLPpCXOM0k9QUOAb7X\n3vSIOAc4B2DUqFE1MVhJT44vH0fWW7nFxszMzOrdAcDkiHgx60DMrHKc2JiZmVm9O5IOuqGZWf1w\nYmNmZmZ1S9JqwL7AFVnHYmaV5WtszMzMrG5FxJvAulnHYWaV5xYbMzMzMzPLPSc2ZmZmZmaWe05s\nzMzMzMws95zYmJmZmZlZ7jmxMTMzMzOz3HNiY2ZmZmZmuefExszMzMzMcs+JjZmZmZmZ5Z4TGzMz\nMzMzy70+WQdgZmbWpvnECWWf56zTDyr7PM3MrPa4xcbMzMzMzHLPiY2ZmZmZmeWeu6KZWa/hbk5m\nZmb1yy02ZmZmZmaWe05szMzMzMws95zYmJmZmZlZ7jmxMTMzMzOz3OsysZE0TNKtkqZLmibpa+2U\n2UPS65KmpI+TKxOumZmZmZnZikoZFa0F+FZETJbUH3hQ0o0R8VhRuTsj4qPlD9HMzMzMzKxzXSY2\nETEPmJc+XyRpOjAEKE5szMx6nUoMIW1mZmbd161rbCQ1AyOB+9qZvIukqZKulbR1B58fJ2mSpEnz\n58/vdrBmZmZm3SFpLUmXSXo87Va/S9YxmVlllJzYSFoDuBz4ekQsLJo8GRgREdsDfwCubG8eEXFO\nRIyKiFGDBg1a2ZjNzMzMSnUmcF1EbAFsD0zPOB4zq5CSEhtJTSRJzT8i4ori6RGxMCIWp8+vAZok\nDSxrpGZmZmbdIGkAsDtwPkBELI2I17KNyswqpZRR0URSIUyPiDM6KLN+Wg5JO6XzXVDOQM3MzMy6\naSNgPvBXSQ9JOk/S6sWF3FXerD6U0mKzK3AUsFfBcM4HSvqypC+nZQ4HHpU0Ffg9MDYiokIxm5mZ\nmZWiD7AD8KeIGAm8AZxYXMhd5c3qQymjot0FqIsyZwFnlSsoMzMzszKYA8yJiLZBjy6jncTGzOpD\nKfexMTMrC0nDgIuA9YFW4JyIOLOojEgu9j0QeBM4JiImVztWqx+VGJJ71ukHlX2eVn4R8YKk2ZI2\nj4gZwN74dhVmdcuJjZlVUyk3/D0A2DR97Az8Kf1rZrYyvgL8Q1Jf4Gng2IzjMbMKcWJjZlVT4g1/\nxwAXpdfp3Zveg2Jw+lkzs26JiCnAqKzjMLPK69YNOs3MyqWTG/4OAWYXvJ6TvmdmZmbWISc2ZlZ1\nXdzwt73BSlYYZdHDs5qZmVkhJzZmVlVd3fCXpIVmWMHrocDc4kIentXMzMwKObExs6op5Ya/wHjg\naCVGA6/7+hozMzPrigcPMLNqarvh7yOSpqTvnQQMB4iIPwPXkAz1PJNkuGePYGRmZmZdcmJjZlVT\n4g1/Azi+OhGZmZlZvXBXNDMzMzMzyz232JiZ9QLNJ07IOgQzM7OKcouNmZmZmZnlnltsrNeoxBnr\nWacfVPZ5mpmZmVn3ucXGzMzMzMxyzy02ZmZm3VTuFmC3/pqZ9ZxbbMzMzMzMLPec2JiZmZmZWe45\nsTEzMzMzs9xzYmNmZmZmZrnnxMbMzMzMzHLPiY2ZmZmZmeWeExszMzMzM8s9JzZmZmZmZpZ7TmzM\nzMzMzCz3+mQdgJmZmVmlSJoFLAKWAy0RMSrbiMysUpzYmJmZWb3bMyJezjoIM6ssd0UzMzMzM7Pc\n67LFRtIw4CJgfaAVOCciziwqI+BM4EDgTeCYiJhc/nCtFjWfOKHs85x1+kFln6eZmfVKAdwgKYC/\nRMQ5xQUkjQPGAQwfPrzK4ZlZuZTSYtMCfCsitgRGA8dL2qqozAHApuljHPCnskZpZmZmtnJ2jYgd\nSH6rHC9p9+ICEXFORIyKiFGDBg2qfoRmVhZdJjYRMa+t9SUiFgHTgSFFxcYAF0XiXmAtSYPLHq2Z\nmZlZN0TE3PTvS8B/gZ2yjcjMKqVbgwdIagZGAvcVTRoCzC54PSd9b17R593Um7FKdBurhLzEaWZm\ntUvS6kBDRCxKn+8H/DjjsMysQkoePEDSGsDlwNcjYmHx5HY+Eiu84aZeMzMzq571gLskTQXuByZE\nxHUZx2RmFVJSi42kJpKk5h8RcUU7ReYAwwpeDwXm9jw8MzMzs5UTEU8D22cdh5lVR5ctNumIZ+cD\n0yPijA6KjQeOVmI08HpEzOugrJmZmZmZWVmV0mKzK3AU8IikKel7JwHDASLiz8A1JEM9zyQZ7vnY\n8odqZmZmZmbWvi4Tm4i4i/avoSksE8Dx5QrKzMzMzMysO0oePMDMzMzMzKxWObExMzMzM7Pcc2Jj\nZmZmZma558TGzMzMzMxyz4mNmZmZmZnlXkk36DQzM7N8aT5xQtnnOev0g8o+TzOzcslFYuPK2czM\nzMzMOuOuaGZWNZIukPSSpEc7mL6HpNclTUkfJ1c7RjMzM8unXLTYmFnduBA4C7iokzJ3RsRHqxOO\nmZmZ1QsnNmVSie5yZvUmIu6Q1Jx1HGZmZlZ/3BXNzGrNLpKmSrpW0tZZB2NmZmb54BYbM6slk4ER\nEbFY0oHAlcCm7RWUNA4YBzB8+PDqRWhmZmY1yS02ZlYzImJhRCxOn18DNEka2EHZcyJiVESMGjRo\nUFXjNDMzs9rjxMbMaoak9SUpfb4TSR21INuozMzMLA96bVc0X+xv5eD9qHskXQLsAQyUNAc4BWgC\niIg/A4cDx0lqAd4CxkZEZBSumZmZ5UivTWzMrPoi4sgupp9FMhy0mZmZWbe4K5qZmZmZmeWeExsz\nMzMzM8s9JzZmZmZW1yQ1SnpI0tVZx2JmlePExszMzOrd14DpWQdhZpXlxMbMzMzqlqShwEHAeVnH\nYmaV5cTGzMzM6tnvgO8CrR0VkDRO0iRJk+bPn1+9yMysrJzYmJmZWV2S9FHgpYh4sLNyEXFORIyK\niFGDBg2qUnRmVm5ObMzMzKxe7QocImkW8C9gL0l/zzYkM6sUJzZmZmZWlyLiexExNCKagbHALRHx\nmYzDMrMKcWJjZmZmZma51yfrAMzMzMwqLSJuA27LOAwzqyC32JiZmZmZWe51mdhIukDSS5Ie7WD6\nHpJelzQlfZxc/jDNzMzMzMw6VkpXtAuBs4CLOilzZ0R8tCwRmZmZmZmZdVOXLTYRcQfwShViMTMz\nMzMzWynlusZmF0lTJV0raeuOCvnOvmZmZmZmVgnlSGwmAyMiYnvgD8CVHRX0nX3NzMzMzKwSejzc\nc0QsLHh+jaSzJQ2MiJd7Om8zM7PeoPnECVmHYGaWez1usZG0viSlz3dK57mgp/M1MzMzMzMrVZct\nNpIuAfYABkqaA5wCNAFExJ+Bw4HjJLUAbwFjIyIqFrGZmZmZmVmRLhObiDiyi+lnkQwHbWZmZmZm\nlolyjYpmZmZmZmaWGSc2ZmZmZmaWe05szMzMzMws95zYmJmZmZlZ7jmxMTMzMzOz3HNiY2ZmZmZm\nuefExszMzKzGDOAN9m+4nyHMzzoUs9zo8j42ZmZmZlZFcx5kQt+TGNaQJDXTW4dxwrKv8lQMyTgw\ns9rmFhszMzOzWjHjWrhgP6Tgi0u/yWnLPsNAvc45TWfQnzezjs6sprnFxszMzKwWLH0TJnwbBm3J\ngc9+hYWsAcAjrRvyj74/44ymsxm37JuEz0ubtctHhpmZmdUlSatIul/SVEnTJP0o65g6NfEsWDgH\nDvjFO0kNwP2xJae1fIZ9GydzRONt2cVnVuOc2JiZmVm9WgLsFRHbA+8H9pc0OuOY2rdwLtz1W9hq\nDDTvusLki5bvx5TWjTmucTyNLM8gQLPa58TGzMzM6lIkFqcvm9JHZBhSx+48A1pbYN8fd1BA/LFl\nDCMaXuKjDROrGppZXjixMbOqkXSBpJckPdrBdEn6vaSZkh6WtEO1YzSz+iKpUdIU4CXgxoi4L+uY\nVrDsbXjk37DVx2Dt5g6L3dS6A4+3DuP4Pv9DtFYvPrOccGJjZtV0IbB/J9MPADZNH+OAP1UhJjOr\nYxGxPCLeDwwFdpK0TXEZSeMkTZI0af78DO4bM2MCvP06vP9TnRYLGji75RA2a3iefRserFJwZvnh\nxMbMqiYi7gBe6aTIGOCitPvIvcBakgZXJzozq2cR8RpwG+2cXImIcyJiVESMGjRoUNVjY8o/YcBQ\n2HD3LotOaB3N3FjHgwiYtcOJjZnVkiHA7ILXc9L3zMy6TdIgSWulz1cF9gEezzaqIgvnwlO3wPZj\noaGxy+LLaeTK5bvx4YapDOT1KgRolh9ObMyslqid99q90DfzriNmlgeDgVslPQw8QHKNzdUZx/Re\nD18K0dplN7RCVyzfjT5q5eDGeyoYmFn+OLExs1oyBxhW8HooMLe9gpl3HTGzmhcRD0fEyIjYLiK2\niYiOhhzLzrT/wtAPwLobl/yRmTGUh1s35LDGOysYmFn+OLExs1oyHjg6HR1tNPB6RMzLOigzs4pY\n9ALMmwqbdTamSvuuWP4htm2YxaaaU4HAzPLJiY2ZVY2kS4CJwOaS5kj6vKQvS/pyWuQa4GlgJnAu\n8H8ZhWpmVnlP3pD83ewj3f7oVct3YVk0cmjjXWUOyiy/+mQdgJn1HhFxZBfTAzi+SuGYmWXrieth\nwBBYb4URqLu0gDWZ2LoV+zfczy85gvYvUTTrXdxiY2ZmZlZtLUvg6dtg031BK5eUXN/6ATZqeIFN\n9Xx5YzPLKSc2ZmZmZtX27D2wdDFs2v1uaG1uWL4jrSE+0vBAGQMzyy8nNmZmZmbV9uQN0NgPNvrw\nSs9iPmvzUGzCRxqd2JiBExszMzOz6pt5EzTvCn1X79Fsrl8+im0bZjFUvp+XmRMbMzMzs2pa9AK8\n/ARsuPKtNW2ub/0AAPs1TOrxvMzyrsvERtIFkl6S9GgH0yXp95JmSnpY0g7lD9PMzMysTsxKh2je\n8EM9ntWzsT6Ptw5jv0YnNmaltNhcCHR256gDgE3TxzjgTz0Py8zMzKxOPXMH9BsA629fltnd3DqS\nUZrBABaXZX5medVlYhMRdwCvdFJkDHBRJO4F1pI0uFwBmpmZmdWVWXfBiA9CY3luJ3jz8h3oo1Y+\n3PBwWeZnllfluMZmCDC74PWc9L0VSBonaZKkSfPn+yI3MzMz62UWzoVXnoLmnndDazMlNmFB9Gev\nxofKNk+zPCpHYtPeXaWivYIRcU5EjIqIUYMGDSrDos3MzMxy5Jk7k79luL6mTSsN3No6kj0bptDI\n8rLN1yxvypHYzAGGFbweCswtw3zNzMzM6susO2CVtWC9bcs625uXj2QtvcEOerKs8zXLk3IkNuOB\no9PR0UYDr0fEvDLM18zMzKy+zLoLRuwKDeW948adrduyNBrZu3FyWedrliddXrUm6RJgD2CgpDnA\nKUATQET8GbgGOBCYCbwJHFupYM3MzMxya+FceHUW7DSu7LNezGrc17ol+zQ4sbHeq8vEJiKO7GJ6\nAMeXLSIzMzOzevTcvcnf4aPpZm6pAAAXjElEQVQrMvtbWkdyStPF8MrTsM5GFVmGWS0rbzuomZmZ\nmbXvuYnQtHrZ7l9T7ObW9B7pT1xfkfmb1TonNmZmZmbV8NxEGDqqbPevWWH2sR5Ptg6BGddWZP5m\ntc6JjZmZmVmlvf06vDgNhu9S0cXc3LoDPHs3vL2wossxq0VObMzMzMwqbc4DEK0Vu76mzc3LR0Jr\nCzx1c0WXY1aLnNiYmZlZXZI0TNKtkqZLmibpa5kF8+xEUCMM/UBFFzM5NoVV1/Z1NtYrVaaTp5mZ\nmVn2WoBvRcRkSf2BByXdGBGPVT2S5+6FwdtBvzUqupjlNMIm+yaJzfKWil3PY1aL3GJjZmZmdSki\n5kXE5PT5ImA6MKTqgbQshecnVfz6mndscRC89UoyWIFZL+LExszMzOqepGZgJHBfO9PGSZokadL8\n+fPLv/B5U6Hl7YpfX/OOTfaBPqvA9KuqszyzGuHExszMzOqapDWAy4GvR8QKw4VFxDkRMSoiRg0a\nNKj8AcxOb8w5rEqJTb81YOO94PEJEFGdZZrVACc2ZmZmVrckNZEkNf+IiCsyCeK5e2HtZui/XvWW\nueXBsHAOzH2oess0y5gTGzMzM6tLkgScD0yPiDMyCSIiSWyqdX1Nm832T0Zhc3c060Wc2JiZmVm9\n2hU4CthL0pT0cWBVI3jlaXjzZRi2c1UXy2rrQPNuTmysV/EYgGZmZlaXIuIuQJkG8Vx6fU21Bg4o\ntNUhMOFb8MKjsP421V++WZW5xcbMzMysUp6bCKusCQM3r/6ytzoUGvrAI/+u/rLNMuDExszMzKxS\nZt+XjIbWkMFPrtXXTUZHe+RyaG2t/vLNqsxd0czMzMwq4c1X4OUnYPux2cWw7SfhyS8kLUfNu2YX\nRw1pPnHCSn921ukHlTESKze32JiZmZlVwnMTk7/VHhGt0BYHQtPq7o5mvYITGzMzM7NKePYeaOwL\nG+yQXQx9V4ctDoJpV0LLkuziMKsCJzZmZmZmlfDsPTBkFDStkm0c2x8Bb7/moZ+t7jmxMTMzMyu3\nJYth3lQYkWE3tDYb7QVrbwgPnJ91JGYV5cTGzMzMrNzmPACxHIZ/MOtIkhHZPvB5eO4eeHFa1tGY\nVYxHRTOzqpK0P3Am0AicFxGnF00/BvgV8Hz61lkRcV5VgzQz66nnJoIaYNhOWUeSeP+n4ZafwAPn\nwUd/m3U01RMBr86C+TNgwUx4Yz4/6/MwECyhLwtiALNifR6PYcyMIWR9P1frGSc2ZlY1khqBPwL7\nAnOABySNj4jHiopeGhEnVD1AM7NyefYeWH9bWGVA1pEkVlsHtjkcpl4K+5ya3DS0HkXAy0/CU7fA\n07fBnPvhzQXvTm9oYt/GVQkaWIWlDNCb70yaH2tyR+t2XLF8Nya2bk2rOzbljhMbM6umnYCZEfE0\ngKR/AWOA4sTGzCy/WpYmXdF2PDbrSN5rpy/ClL/D/efC7t/OOpoea7sfjWhlpGayf+MD7NcwieaG\nFwF4pnU9HmjdhodiE2a0DuPpGMxrrEFhq8wqLKFZL7Jtw9Ps2vAo+zQ8yMcb72RODOS8lgP51/I9\neZt+WayerQQnNmZWTUOA2QWv5wA7t1Pu45J2B54AvhERs4sLSBoHjAMYPnx4BUI1M1tJcx+Clrdr\nY+CAQhu8HzY7AO75Q5Lk5LnVprWVHTWDgxrv44DG+xmsV1gSfZjYujXnLTuQ21q3Z068r8vZvE0/\nHo/hPL58OP9Zvgf9WMo+DZM5us8NnNp0Ecf3uZLftnyCfy3f0y04OeDExsyqqb3Oy1H0+irgkohY\nIunLwN+AvVb4UMQ5wDkAo0aNKp6HmVl2Zt2R/B2xW7ZxtGfP78FfdoeJZyfP86R1Ocy+Hx77Hzz2\nPy7vN5cl0cTtrdtx+vKx3NK6A4tYrUeLWEJfJrSOZsLS0XxAj/Ptpn/zs6bzOarxBr6/7PPAQeVZ\nF6sIJzZmVk1zgGEFr4cCcwsLRERBZ2jOBX5RhbjMzMrnmTtgvW1h9XWzjmRFg7eHLQ+Ge8+Gnb+U\nXHtTy5a9BU/fDjOuSR5vzIfGfrDJPnz9kY9xU+sOLO5hMtORB2ILjlj6Q/ZveIAfNP2dy/r+CK6Z\nB3ufDP3WqMgyrWdKalOTtL+kGZJmSjqxnenHSJovaUr6+EL5QzWzOvAAsKmkDSX1BcYC4wsLSBpc\n8PIQYHoV4zMz65llb8Nz98GGu2cdScf2OAmWLoabTs06khW1Xfx/31/g74fDL5rhkiPg0StgxK5w\n+AXwnZlw5D+5snW3iiU17xLXte7ER5b8gr8t3w/uPwf+vFvScmQ1p8sWG49iZGblEhEtkk4AricZ\n7vmCiJgm6cfApIgYD3xV0iFAC/AKcExmAZuZddec+2H5ktpObNbbCnY5Ae75PWw1BjbZO9t4Fs6F\nZ+5MWrqeuR1eTy+rXGcj2OGzsNl+0Pwh6JPdRfxvsCo/avksx37x63DFl+CCj8Du34XdvwON7gBV\nK0r5JjyKkZmVTURcA1xT9N7JBc+/B+Ss47eZWeqZO0CNMKIGbszZmT1PghnXwlVfg+Puqe6w1Itf\nSrbTrDuThOaVp5L3V10bmneD3b4BG++ZJDa1ZsQH4bi74drvwu2nJ4nYYefCWsO6/qxVXCmJjUcx\nMjMzMyvFM3fABiNr5/41HWlaFT52Npy/H/z3S/DJiyvX8rDsLZh1d3pvmVvhpfTceL8BSaIw6nOw\n4YeS65IacjDy2CoD4NA/w8Z7wdXfSLqmjTkruXbJMlXKHuxRjMzMzMy6smQRPP8gfPCrWUdSmmE7\nwYG/gmu+nbTcjDkL1N7PvpXw2nMw4zp48nqYdVcy/HVjPxg+GvY+BTb8MAzenubvXw8PQ3LefE55\nll0t230Sho6Cyz4Hl34GRn0ePvLTJGm0TJSS2HgUIzMzM7OuPDsRWltq+/qaYjt9Ed54OelWFa1w\n0K+h7+rdnk3ziVezjZ5h38YH2a/hQbZseA5IbpJ5a+se3N66Pfe1bsHb0/ulQ8K8kD5ybp2N4HM3\nwM0/golnwbP3wOHnw3pbZx1Zr1RKYvPOKEbA8ySjGH2qsICkwRExL33pUYzMzMys95l5I/RZNWmV\nyJM9TgQCbv9lMvjBwb9Puoh11XqzZFFyjczMG5nY70oG6xWWh5gUm/OTZZ/m5tYdeCYGdz6PetCn\nb9JSs/Ge8N/j4Jw9YZ9TYecv56NrXR3pMrHxKEZmZmaWR5IuAD4KvBQR21R0YRHw5A1Ja03euiJJ\nyWACI3ZNrre58EBYf1vY5uPwvq1hzSFJa86SRbDgKXhpOsy+F+ZNTVqomlZnSutW/Kb1E9y8fCSv\nUuPXF1XKJvskAzGM/wpc/z144joY80cPLFBFJV0l5lGMzMzMLIcuBM4CLqr4khbMhFdnwQe/UvFF\nVcxGH4YTJsEj/4EHzu/4PjeN/WDIjsm6brwXDNuZ435wU1VDrVlrDIIjL4HJf4PrToKzd4H9fwYj\njyrf9UvWIQ+8bWZmZnUpIu6Q1FyVhT15Q/J3k32rsriK6bcGjDo2ebz5Csx/HBa/CA19kpaodTaC\nNYf73i2dkWDHY2CjPeB/JyQtOI9eDgefCWs3ZxtbnfNeaWZmZr1aWW5H8eQNMGgLWHtEGSPL2Grr\n1P79eGrZ2s1w9Hh48AIWX/19Gn63E79p+QQXLv8Iy2kseTazTj+ocjHWGSc2ZmZm1qv1+HYUSxYn\n92kZ/eVyh1ZVzSdOyDqEmrdy22gwG/BLTmv6Kz9s+jtjGu/me8u+yLRoLnd4vZ6HajAzMzPriWdu\nh9ZlsOl+WUdiNWouA/n8sm9z/NKvMlivML7v9/l+n7+zKm9nHVpdcWJjZmZm1hPTr4J+a8KwnA3z\nbFUmJrSOZu8lv+LS5XvyxT7XcEPf/8duDY9kHVjdcFc0MzMzq0uSLgH2AAZKmgOcEhHnl3UhLUvg\n8Qmw5cHJ/UxqgLuU1baFrMFJLV/gv8t34/Smc/l735/z75YPc1rLUSxitazDyzW32JiZmVldiogj\nI2JwRDRFxNCyJzUAM2+GJQth68PKPmurbw/EFhy49Of8seUQPt54B9f3+y4fang467ByzYmNmZmZ\n2cqa9l9Yde3kHjBm3bSEvvyqZSyHLf0Rb8SqXNz3dE7ucxH9WJp1aLnkxMbMzMxsZSx7C2Zck3RD\na2zKOhrLsamxCR9d+lMubNmPz/W5jiv7nsxGmpt1WLnjxMbMzMxsZcy8GZYuhq0PzToSqwNL6Mup\nLcdwzNLvsJ5e4aq+3+eQhruzDitXnNiYmZmZrYypl8BqA6F596wjsTpyW+tIDlzyc6ZFM7/v+0e4\n9kRYvizrsHLBiY2ZmZlZdy2cCzOuhZGfhkYPMmvl9QLr8qml3+f8lgPgvj/BRR+DN17OOqya58TG\nzMzMrLsmXwyxHHY8JutIrE610IfTWo6Cw86F5yfBuXvCi9OyDqumObExMzMz647lLTD5b7DxXrDO\nRllHY/Vuu0/CsddAy1I4fz944vqsI6pZTmzMzMzMuuPJG2Dh8zDqc1lHYr3FkB1h3K2w7sZwyVi4\n908QkXVUNceJjZmZmVl3TLoA+g+GzfbPOhLrTQZsAMdeC5sfCNedCNd8J2k9tHf4ajczMzOz7jjg\nF/DqLN+7xqqv7+rwyYvhplPgnt/Dq8/A4X+FVQZkHVlNcIuNmZmZWXesuzFssnfWUVhv1dAA+50G\nB58JT90KF3wEXn0266hqgltszMzMzMxqVPOJEzqYMohdG77L2S+eybLf7caXln6DB2Pzd6bOOv2g\n6gRYQ9xiY2ZmZmaWQ3e3bsuhS3/MwliNS/r+hCMbb846pEw5sTEzMzMzy6mnYwM+tvQ07mndhp83\nnc/P+5xLP5ZmHVYm3BXNzMzMzCzHFrI6n1v2Hb4Z/+GEPv9ju4anYcF2yfVgvYhbbMzMzMzMcq6V\nBn7dcgTHLv0OQ/Qy/Hk3mPTXXnW/Gyc2ZmZmZmZ14tbWkey/5HQY+gG4+uvwj0/0mlHTnNiYmZmZ\nmdWRF1gXjroSDvglPDcR/rgz3PFrWPZW1qFVlBMbMzMzM7N609AAO38Jjr8/ue/SLafBH3aEyRdB\nS30OLuDExszMzMysXq05BMb+A46ZAP3Xh/FfgTO3h7t/D28syDq6snJiY2ZmZmZW75p3gy/cDJ++\nPBkt7cYfwm82h0uPgkcvh7cXZh1hj5WU2EjaX9IMSTMlndjO9H6SLk2n3yepudyBmll9cH1iZtXU\nVZ1j1qtIsOk+cMzVcNxE2OmLyTU4l30OfrkRXHAA3PJTmHEdLJyXuxHVuryPjaRG4I/AvsAc4AFJ\n4yPisYJinwdejYhNJI0FfgEcUYmAzSy/XJ+YWTWVWOeY9U7rbQX7/xz2+wnMvh9mXAOz7oI7fw3R\nmpTpt2bSurP2CBgwBNZYD1YfCKuuA/36Q781oGl1aFoFGvtBn77Q0ASNTaDG5DqfKirlBp07ATMj\n4mkASf8CxgCFlcIY4NT0+WXAWZIUkbM0z8wqzfWJmVVTKXWOWe/W0AgjdkkeAEsWwQuPwrypsOBJ\nWDAT5j2ctOK0rMSoampIHih9ruT9TfZJrv0po1ISmyHA7ILXc4CdOyoTES2SXgfWBV4uLCRpHDAu\nfblY0owSlj+weD41ynGWl+MsI/2i5DhHVDiUrOuTSsjFPoDjLLc8xFn2GPWLbhWvdH1SilLqnFqq\nT4rlYT/ritchI0XHaw2uwz/hyH+WUrDkuqSUxEbtvFd85rSUMkTEOcA5JSzz3RlLkyJiVHc+kwXH\nWV6Os7xqKM5M65NKqKFt2ynHWV55iDMPMVZBruqTYvXwHXodakM9rEMpSun4NgcYVvB6KDC3ozKS\n+gBrAq+UI0AzqyuuT8ysmkqpc8ysTpSS2DwAbCppQ0l9gbHA+KIy44HPps8PB25xf3gza4frEzOr\nplLqHDOrE112RUv7uJ8AXA80AhdExDRJPwYmRcR44HzgYkkzSc6sji1jjDXXNNwBx1lejrO8aiLO\nGqhPKqEmtm0JHGd55SHOPMRYUR3VORmH1R318B16HWpDPaxDl+QToWZmZmZmlnfVHVzazMzMzMys\nApzYmJmZmZlZ7tVMYiNpf0kzJM2UdGI704+RNF/SlPTxhVqMMy3zSUmPSZomqaQBusuthO3524Jt\n+YSk12o0zuGSbpX0kKSHJR1YgzGOkHRzGt9tkoZWO8Y0jgskvSTp0Q6mS9Lv0/V4WNIO1Y4xz/Kw\nr5YYZ+b7a1721RLi3ELSRElLJH272vEVxNFVnJ9Ot+PDku6RtH21Y7SOlXDM9pN0aTr9PknN1Y+y\ncyWswzfT30UPp/VPLdzjaAWl/MZLyx0uKSTV3PDJefmdWjERkfmD5IK+p4CNgL7AVGCrojLHAGfl\nIM5NgYeAtdPX76vFOIvKf4Xkgsqai5PkYrfj0udbAbNqMMb/AJ9Nn+8FXFztbZkue3dgB+DRDqYf\nCFxLcl+H0cB9WcSZx0ce9tVuxJn5/pqXfbWEON8HfAD4KfDtLGIsMc4PFvxPOsDHfu08Sjxm/w/4\nc/p8LHBp1nGvxDrsCayWPj+u1tah1PVIy/UH7gDuBUZlHfdKfBeZ/06t5KNWWmx2AmZGxNMRsRT4\nFzAm45jaU0qcXwT+GBGvAkTES1WOEbq/PY8ELqlKZO9VSpwBDEifr0n17z9QSoxbATenz29tZ3pV\nRMQddH6/lzHARZG4F1hL0uDqRJd7edhXISf7a1721a7ijIiXIuIBYFn1omo3jq7ivKftfxLJj7FM\nWpWtXaUcs2OAv6XPLwP2ltTejUez0uU6RMStEfFm+rJW98FSfzudBvwSeLuawZUoL79TK6ZWEpsh\nwOyC13PS94p9PG3GvEzSsHamV1opcW4GbCbpbkn3Stq/atG9q9TtSdocvCFwSxXiKlZKnKcCn5E0\nB7iGpHWpmkqJcSrw8fT5oUB/SetWIbbuKnm/sBXkYV+F+tlfva9WzudJWsOsNpSyr79TJiJagNeB\nWjpmu3u81uo+2OV6SBoJDIuIq6sZWDfk5XdqxdRKYtPemYficaivApojYjvgJt49e1FNpcTZh6SZ\nbw+SlpDzJK1V4biKlRJnm7HAZRGxvILxdKSUOI8ELoyIoSTdUy6WVM39tpQYvw18WNJDwIeB54GW\nSge2ErqzX9h75WFfhfrZX72vVoCkPUl+VP6/rGOxd5Syr9f68VByfJI+A4wCflXRiFZOp+uR1ue/\nBb5VtYi6Ly+/UyumVhKbOUBhC8xQirpxRMSCiFiSvjwX2LFKsRXqMs60zP8iYllEPAPMINmBqqmU\nONuMJZtuaFBanJ8H/g0QEROBVYCBVYkuUcq+OTciDouIkcD30/der16IJevOfmHvlYd9Fepnf/W+\nWmaStgPOA8ZExIKs47F3lPq7YhiApD4kXV0768pZbSUdr5L2IalzDin4PVdLulqP/sA2wG2SZpFc\n/ze+xgYQyMvv1IqplcTmAWBTSRtK6kvyY3t8YYGi/tWHANOrGF+bLuMEriS5SA5JA0ma/J6uapSl\nxYmkzYG1gYlVjq9NKXE+B+wNIGlLkh+L82spRkkDC87Mfw+4oIrxdcd44Oh0xKnRwOsRMS/roHIi\nD/sq1M/+6n21jCQNB64AjoqIJ7KOx96jlLplPPDZ9PnhwC2RXvVdI0qpd0YCfyFJamr1mo5O1yMi\nXo+IgRHRHBHNJNcKHRIRk7IJt115+Z1aOVmPXtD2IOm68QTJaA7fT9/7MclOA/BzYBpJ//BbgS1q\nNE4BZwCPAY8AY2sxzvT1qcDpNf69bwXcnX7vU4D9ajDGw4En0zLnAf0y2paXAPNILmSeQ9KC8GXg\nywX75h/T9XiEGhvNpdYfedhXS4wz8/01L/tqCXGun76/EHgtfT6gBuM8D3g13S+nAJOy2J5+dPj9\ndXXMrkIymuFM4H5go6xjXol1uAl4sWAfHJ91zCuzHkVlb6vF/6MlfBc18Tu1Ug+lK2lmZmZmZpZb\ntdIVzczMzMzMbKU5sTEzMzMzs9xzYmNmZmZmZrnnxMbMzMzMzHLPiY2ZmZmZmeWeExszMzMzM8s9\nJzZmZmZmZpZ7/x+SSGAd5s6fdAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x2271275e550>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# check that the KDE is reasonable for the q_data\n",
    "q_values = np.linspace(-0.1,0.75,100) # plot pdf for these \"x-values\"\n",
    "\n",
    "# check that the simulation makes sense by plotting!\n",
    "f, (ax1, ax2, ax3) = plt.subplots(1, 3)\n",
    "f.set_figwidth(14)\n",
    "\n",
    "# Comparison of exact pdf vs. density\n",
    "ax1.set_title('Histogram: Lambda 1 (\"Actual\")')\n",
    "ax1.hist(lam1_data, density=True)\n",
    "ax2.set_title('Histogram: Lambda 2 (\"Actual\")')\n",
    "ax2.hist(lam2_data, density=True)\n",
    "\n",
    "# histogram of data\n",
    "ax3.set_title(\"Histogram of Data (QoI) and KDE of QoI\")\n",
    "ax3.hist(q_data, density=True) \n",
    "ax3.plot(q_values,data_kde_pdf(q_values))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solve for $\\lambda_1$ and $\\lambda_2$ Using Consistent Bayes\n",
    "Now that we have the dataset and a model for our simulation, we attempt to construct a solution distribution for the consistent Bayes method using the Accept-Reject paradigm.\n",
    "\n",
    "The purpose of using the paradigm in this example is to show how we **may not even need to assume a prior on the stochastic term** to calculate the distribution of $\\lambda_1$ and $\\lambda_2$!\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Construct the Push-Forward KDE\n",
    "We construct the push-forward KDE for the denominator of the consistent baye's ratio by sampling the prior. Note that we don't need to sample a prior for our stochastic error term because it is embedded in the calculation of the map $Q$!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 231,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# defines prior class objects using scipy package (imported at the top)\n",
    "# use \"rvs\" method to generate samples from these objects\n",
    "# use \"pdf\" method to evaluate the pdf at a given value\n",
    "lam1_prior = uniform(0.5,0.5)\n",
    "lam2_prior = uniform(0.75,0.5)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 232,
   "metadata": {},
   "outputs": [],
   "source": [
    "# generate a sample from our PRIOR pdfs\n",
    "lam1_sample_prior = lam1_prior.rvs(1000)\n",
    "lam2_sample_prior = lam2_prior.rvs(1000)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 233,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\ProgramData\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:8: DeprecationWarning: elementwise != comparison failed; this will raise an error in the future.\n",
      "  \n",
      "C:\\ProgramData\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:13: DeprecationWarning: elementwise != comparison failed; this will raise an error in the future.\n",
      "  del sys.path[0]\n",
      "C:\\ProgramData\\Anaconda3\\lib\\site-packages\\scipy\\optimize\\minpack.py:161: RuntimeWarning: The iteration is not making good progress, as measured by the \n",
      "  improvement from the last five Jacobian evaluations.\n",
      "  warnings.warn(msg, RuntimeWarning)\n",
      "C:\\ProgramData\\Anaconda3\\lib\\site-packages\\scipy\\optimize\\minpack.py:161: RuntimeWarning: The iteration is not making good progress, as measured by the \n",
      "  improvement from the last ten iterations.\n",
      "  warnings.warn(msg, RuntimeWarning)\n"
     ]
    }
   ],
   "source": [
    "q_sample_prior = Q(lam1_sample_prior,lam2_sample_prior)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 234,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x2271279d3c8>]"
      ]
     },
     "execution_count": 234,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAEICAYAAAB/Dx7IAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3Xl8VNX9//HXZ9YkBLJA2JewiMq+\nBFREFAQBN3BDrFqhVrQu1da2trX+amvV2m+tXdRaV1xxQxRckDWyq6Cgsq+BhCUhbNlnO78/7kRH\nTMgAmUxu8nk+HvNg5s6dez9zZvLmzrnn3ivGGJRSStmHI94FKKWUOj4a3EopZTMa3EopZTMa3Eop\nZTMa3EopZTMa3EopZTMa3HVIRLJF5KfxriOWRCRTRIyIuOJdC4BYXhCRgyLyWbzriafw59Ktrl9b\nxbKeEpH7amNZjZUGdzVEZIeIlIlIsYjsC//xJ9dxDdkiUh6uofJ2Vl3WUNvC7Toy4vHEcKieGxH6\nxRHt/r6IjKpiGWVHtcvj1axyKDAKaG+MGRzDt1YlEfGIyKMikhuuc7uIPFbXdcTSUd/T/SLyjoi0\nqW5+Y8wtxpgH6rLGhkaD+9guMcYkAwOAQcAf4lDD7caY5Ijb8uN5cXiLMyaf88luVYvIDcATwEXG\nmE8inkoNt3tfYC4wQ0QmHfXyS45ql9urWU0nYIcxpiQW7yEKvwOygMFAU2A48GWM1xkPt4c/s+5A\nKlDlf04i4jzRFdSXX3H1gQZ3FIwxecBHQC+ocqvxfhF5JXw/QUReEZFCETkkIp+LSKuIxXUSkaUi\nUiQic0SkxYnUJCJDwss+HP53SMRz2SLyoIgsBUqBG0Xk64jn50V2G4jIEhEZH77/WxHZGq5vnYhc\nFjHfpHDtj4nIAeB+EXGKyN/DW1rbgIuirH8K8Cgw2hizrKp5jDF7jTH/Au4HHjne/4BE5EbgWeCs\n8Nbgn0TkvPDW7z0ishd4ITzvTSKyRUQOiMhMEWkbsRwjIreKyOZwuzwgIl1FZLmIHBGRN0XEU00Z\ng4AZxpjdxrLDGPNSxLKjbe9DIrIt/LlPEpFdIpIf/s+vcv6pYnVDzA0v7xMR6VRN23jDn9vO8C+b\np0QkMeL5X4vIHhHZLSI/ibbNjTEHgOl897cyVUT+KyIfikgJMDw87S8R66qp7W8Tkc3A5vCGyGPh\n935YRL4SkV7R1tdgGGP0VsUN2AGMDN/vAKwFHjj6ufDj+4FXwvdvBmYBSYATGAg0Cz+XDWzF2ipJ\nDD/+6zFqyAZ+WsX0dOAgcD3gAq4JP24e8bqdQM/w802BMqBF+PFeYHd4emL4ucrXXgW0xfpP/Wqg\nBGgTfm4SEADuCC8nEbgF2BBuo3RgIWAA1zHadTqwD+h71HOZVb0W6BKefnpV7V/D5zgJWBLx+Lzw\ne3gE8IbfwwhgP9YvKy/wH2BRxGsMMBNoFm7TCmB+uK4UYB1wQzXr/0P4s7gV6A3IUc9H096Tw9+l\nv4SX9US4zguAIiA5PP/U8ONh4ef/ddR7N0C38P1/ht9Tevh7MAt4OPzcmPDn0wtoArwW+dpjfU+x\nvmMLgJcjajoMnB1+jwnhaX8JPx9N288N15kIjAZWYW3VC3B6ZXs1plvcC6ivt3A4FAOHgBzgSSAx\n4rnqgvsnwDKgTxXLzAb+EPH4VmD2MWrIxtpiPhS+fRGefj3w2VHzLgcmRbzuz0c9vxi4HDgTmAO8\nGf4DHQ58dYwaVgPjwvcnATuPen4BcEvE4wuoObiPAO8BjqOey6zqteE/dgOcXcVnU3m7qZr1TeKH\nwe0DEiKmPQf8LeJxMuAHMsOPv113+PEq4J6Ix48C/6xm/U7gNmApVuDvppqQr6a9N0c81ztcS6uI\naYVAv/D9qcDrR72PINAh4n10wwq8EqBrxLxnAdvD958nYoMCa0OjpuCu/J7mAa8CGRE1vXTU/FP5\nLrijafsREc+PADZhfY8dVdXTGG7aVXJs440xqcaYTsaYW40xZVG85mXgY+D18M/Mv4mIO+L5vRH3\nS7G+qJV72it3tP0+Yp6fh2tINcYMCE9ri/WfSaQcoF3E411HPf8JVmgNC9/PBs4N377tXxaRH4vI\n6vBP80NYW12R3TlHL7ftUdOOrqsqt2CFwbMiIlHMX/m+DkRMGx/RLqnGmGeiWE6lAmNMecTj77Wn\nMaYYKxAj23NfxP2yKh5XuePaGBM0xjxhjDkbayvxQeB5ETkdomrvo9eDMeZY6/72swi/jwPh9xcp\nA+sX4aqI9c4OT4cT+0wrv6ftjDHXGmMKqqqpCtG0feR7WgA8jvWrY5+IPC0izaKor0HR4D4xJVhf\n/EqtK+8YY/zGmD8ZY3oAQ4CLgR/XtEBj7Wmv3NH2UA2z78ba6RapI9bWzreLPOr5o4P7E44K7nB/\n6DPA7VhdJ6nAN1hbaNUtdw9WN0lkHTXJB84HzsH6JVOTy8Kv2RjFvNE4+j18rz1FpAnQnO+358mv\n1JgyY8wTWN1aPaJs7+P17Wch1iiodKz3F2k/VuD3jPiPL8VYOxfhxD7TYznWKUijafvvvd4Y829j\nzECsbqvuwK9Psj7b0eA+MauBiSLiFpEs4MrKJ0RkuIj0Fmvv+RGsn33BWl7/h0B3EfmRiLhE5Gqg\nB/D+MV6zDDgVa3TDZ8aYtVh/MGcAi8LzNMH6IykIv5fJhHcyHcObwM9FpL2IpAG/jeYNGGN2Y/3s\nHSPVDI8TkVYicjvwR+B3xphQNMs+Aa8Bk0Wkn4h4gYeAT40xO052wSJyV3iHaGL4s7oBq0/5S06s\nvWtyoYgMDe8sfQDrfXxvizfcjs8Aj4lIy/C624nI6PAsbwKTRKSHiCRhtX+sHFfbi8ggETkj/Cu2\nBCin9v++6j0N7hNzH9AVa8vpT1hfvkqtgbexQns91tbsK7W5cmNMIdaW/N1YPyt/A1xsjNl/jNeU\nAF8Aa40xvvDk5UCOMSY/PM86rP7a5Vg/0Xtj9c0eyzNYXUNrwst/5zjexy6s8L5SRB6OeOpQeATC\n18CFwFXGmOePevks+f447hnRrreKOuZjfabTsbY2uwITT3R5RynDatO9WFu6twFXGGO2nWB71+Q1\nrKA9gLVj/Npq5rsH2AKsEJEjwDys/9gxxnyEtfNyQXieBSdZU7VOoO2bYX3nDmJ1sRQCf49VffWV\nhDv8lVI2JyJTgVxjTDyON1B1SLe4lVLKZjS4lVLKZrSrRCmlbEa3uJVSymZictKWFi1amMzMzFgs\nWimlGqRVq1btN8Zk1DxnjII7MzOTlStXxmLRSinVIIlINEeoAtpVopRStqPBrZRSNqPBrZRSNqPB\nrZRSNqPBrZRSNqPBrZRSNqPBrZRSNqPBrZRSNqPBrZRSNhOTIycbgszfflDry9zx14tqfZlKqcZH\nt7iVUspmNLiVUspmNLiVUspmNLiVUspmNLiVUspmNLiVUspmNLiVUspmNLiVUspmNLiVUspmNLiV\nUspmNLiVUspmNLiVUspmogpuEUkVkbdFZIOIrBeRs2JdmFJKqapFe3bAfwGzjTFXiogHSIphTUop\npY6hxuAWkWbAMGASgDHGB/hiW5ZSSqnqRNNV0gUoAF4QkS9F5FkRaXL0TCIyRURWisjKgoKCWi9U\nKaWUJZrgdgEDgP8aY/oDJcBvj57JGPO0MSbLGJOVkZFRy2UqpZSqFE1w5wK5xphPw4/fxgpypZRS\ncVBjcBtj9gK7ROTU8KTzgXUxrUoppVS1oh1VcgfwanhEyTZgcuxKUkopdSxRBbcxZjWQFeNalFJK\nRUGPnFRKKZvR4FZKKZvR4FZKKZvR4FZKKZvR4FZKKZvR4FZKKZvR4FZKKZvR4FZKKZvR4FZKKZvR\n4FZKKZvR4FZKKZvR4FZKKZvR4FZKKZvR4FZKKZvR4FZKKZvR4FZKKZvR4FZKKZvR4FZKKZvR4FZK\nKZvR4FZKKZvR4FZKKZvR4FZKKZvR4FZKKZtxRTOTiOwAioAgEDDGZMWyKKWUUtWLKrjDhhtj9ses\nEqWUUlHRrhKllLKZaIPbAHNEZJWITKlqBhGZIiIrRWRlQUFB7VWolFLqe6IN7rONMQOAscBtIjLs\n6BmMMU8bY7KMMVkZGRm1WqRSSqnvRBXcxpjd4X/zgRnA4FgWpZRSqno1BreINBGRppX3gQuAb2Jd\nmFJKqapFM6qkFTBDRCrnf80YMzumVSmllKpWjcFtjNkG9K2DWpRSSkVBhwMqpZTNaHArpZTNaHAr\npZTNaHArpZTNaHArpZTNaHArpZTNaHArpZTNaHArpZTNaHArpZTNaHArpZTNaHArpZTNaHArpZTN\naHArpZTNaHArpZTNaHArpZTNaHArpZTNaHArpZTNaHArpZTNaHArpZTNaHArpZTNaHArpZTNaHAr\npZTNaHArpZTNRB3cIuIUkS9F5P1YFqSUUurYjmeL+05gfawKUUopFZ2ogltE2gMXAc/GthyllFI1\niXaL+5/Ab4BQdTOIyBQRWSkiKwsKCmqlOKWUUj9UY3CLyMVAvjFm1bHmM8Y8bYzJMsZkZWRk1FqB\nSimlvi+aLe6zgUtFZAfwOjBCRF6JaVVKKaWqVWNwG2N+Z4xpb4zJBCYCC4wx18W8MqWUUlXScdxK\nKWUzruOZ2RiTDWTHpBKllFJR0S1upZSyGQ1upZSyGQ1upZSyGQ1upZSyGQ1upZSyGQ1upZSyGQ1u\npZSyGQ1upZSyGQ1upZSyGQ1upZSyGQ1upZSyGQ3uahkudSxjtuce/u5+ivaSH++ClFIK0OCu2oFt\nvOT+K//2PI6TEJc4lrPAczd/dr1ACsXxrk4p1cgd19kBG4X9m+F/w+jvMPw//w28EhxFSw5yh+td\nrnEuoIPkM9n/G0DiXalSqpHS4D7ax/eCw8WYigfIw7oE216ac2/gRraYtvzR/TKXBpczMzQkzoUq\npRor7SqJtGUebP4Yhv3629CO9GJwNKtDXfl/7pdIpSgOBSqllAb3d4IBa2s7rTOccXOVs4RwcI//\nJlIo4Q/uV+u4QKWUsmhwV1r1AhRsgNEPgstb7WwbTUeeCl7Clc5FDHV8XYcFKqWURfu4AcoOwcIH\nofMwOPVCa5r48KQvxZ2+GEKJBMs6ECxrT+BIfx4PjGecYyl3uqazxNc7vrUrpRod3eIG+PotKDsI\nI/+EAd7d8i5Nuv4db8uPCZV3IFjeGmfSNhJav09S5pP43EW8GLyAQY5NnC458a5eKdXIaHADrHkd\nWvWCdgN4Y+Mb3Lf0Pow/ldIdN1O2azLleddTsuX3lO64BXGWkpT5FG+5TqPcuLneOTfe1SulGhkN\n7v2bIW8l9J3ImoI1PPL5IwxrP4zSnFsIlnX+3qzBskxKd9wCRgh1epknXQMZ71xKU0rjVLxSqjGq\nMbhFJEFEPhORNSKyVkT+VBeF1Zk1r4M4ONB9FHdn302rpFY8NPQhqmuakK8VpTm3YIKJvNS6CBw+\nrnAuqtualVKNWjRb3BXACGNMX6AfMEZEzoxtWXUkFIKv3iDYZTi/+eLvHKo4xGPnPUaKN+WYLzP+\ndMp3TyDoKeK+1E7h7hJTNzUrpRq9GoPbWCpP0OEO3xpGSuUshcO7mNnhdD7d8ym/P+P3nN789Khe\nGizrjO/AWcxJCVGUWMgQx9oYF6uUUpao+rhFxCkiq4F8YK4x5tMq5pkiIitFZGVBQUFt1xkba17H\n72nK//Z/Ts/mPbms22XH9fKKgjGYQAr3tsjgKte8GBWplFLfF1VwG2OCxph+QHtgsIj0qmKep40x\nWcaYrIyMHx4uXu/4SmHdu8zoOoi8kj3c1u82RI7zxFEhL2V7rmCnx0lO+haSKI9NrUopFeG4RpUY\nYw4B2cCYmFRTlzbPocJfzNOh/fTN6MvQdkNPaDHBku40O5LJaylJnONeUctFKqXUD0UzqiRDRFLD\n9xOBkcCGWBcWc5vn8nZaBvt8h7i9/+3Hv7UdYW/BeMpF8DbX0SVKqdiLZou7DbBQRL4CPsfq434/\ntmXFWChE2Za5PJvSjKxWWZzR+oyTWlzQ15ruxSmsSCmmibOwlopUSqmqRTOq5CtjTH9jTB9jTC9j\nzJ/rorCY2vsVs6SE/fhPrG+7CuUF51MuQufmM2qhQKWUql7jPHJy81zebppM95QuDGw1sFYWubZi\nEMNKAuxN2wrOklpZplJKVaVRBvfazR+w3uvhqtOuqZWtbQCDgw6Fp+GTEE3TF9TKMpVSqiqNL7hL\nCnmrLIcEcXJRl4tqddHLy4YzsrQMV9qnIL5aXbZSSlVqdMFdsukjPkxOYkzrITT1NK3VZa82XRl9\nyEHAGcCd8mWtLlsppSo1uuD+cOOblDkcXNl3SgyWLuwsHcBpFX686UtoKGcGUErVL40ruENB3ire\nQndHEn1a9o3JKuYGB3H9kSOItwBn0paYrEMp1bg1quBeu/4d1rudXNnu3FrbKXm0VaY7g4uEJgEH\nnvSlMVmHUqpxa1TB/f6mN3Ebw0UDbovZOoI4WRwcwISiYpzJGxH3/pitSynVODWa4A6Ggnx8eDPn\nhDw0S+0U03XNDQ3k+qIDOBA86ctjui6lVOPTaIL7i93LKZAgY9P7xHxdi0J9SA64OLUoGXfKKh0a\nqJSqVY0muD9c+wqJoRDDTr085usqx8viUG9+WrQPcZbjarYm5utUSjUejSK4/SE/c/M/57zScpI6\nn1cn65wbGsgFvgLcFWl40n5w3QmllDphjSK4V+xeweGQj7HeNpDQrE7WOT84gJARsg4n40zMxZGQ\nWyfrVUo1fI0iuD/aOpOmwRBndxpRZ+s8QDNWmlO5uWQnJuTGnapb3Uqp2tHgg7s8UM6CXQsZWVqK\np8vwOl33nOBABrIL75FTcaespshXVKfrV0o1TA0+uJfkLaEkWMGY0grocHIXTDhec0JZAAw77EIc\nft7fZu/rTyil6ocGH9xzcuaQZoTBGX3Bk1Sn695lWrE+1JGJgU0Ey9rx5sY3MUbPX6KUOjkNOrh9\nQR+Ldn3C8OIiXJ3PjUsNc0JZZMlGPIf6sOXQFlYXrI5LHUqphqNBB/eKPSsoCZRyfkkpZJ4Tlxo+\nDmbhFMOIYh/J7mTe2PhGXOpQSjUcDTq45++cTxNxcaYfaJ8VlxrWmU7kmhZcKKu5pOslzNkxh4Pl\nB+NSi1KqYWiwwR0IBVi4cyHD/IKnw2BweeNUiTAnmMU5jm+Y0Pli/CE/7255N061KKUaggYb3F/m\nf8nBioOMKtwTt26SSh8HB+EVP90KtjGg5QDe2vQWIROKa01KKftqsME9N2cuXnFxdlkZZA6Nay0r\nTXcOmGTY8AETTp3ArqJdrNizIq41KaXsq8bgFpEOIrJQRNaLyFoRubMuCjsZIRNi/s75nO1pTpLT\nC+0GxrWeIE7mBQfCptmMajeMNG8ab258M641KaXsK5ot7gBwtzHmdOBM4DYR6RHbsk7ON/u/Ib80\nn5FFRdB+UBz7t7/zQehMqDiCZ/tixp8ynuxd2ewp3hPvspRSNlRjcBtj9hhjvgjfLwLWA+1iXdjJ\nmLdzHi5xMmzPprj3b1daGuoJiemw9h0mnjoRg+H1ja/HuyyllA0dVx+3iGQC/YEfnDFJRKaIyEoR\nWVlQUFA71Z0AYwwLdi5gUNMupIRCkHl23GqJFMAFp18CGz6krSeV8zuez9ub3qYsUBbv0pRSNhN1\ncItIMjAduMsYc+To540xTxtjsowxWRkZGbVZ43HZemgrOUdyOJ8EcHqhXXzGb1ep1xXgL4HNc7j2\n9Gs54jui5y9RSh23qIJbRNxYof2qMead2JZ0chbsWgDA8Pyd0GEwuBPiXFGEzKHQpCWsfYcBLQdw\nevrpvLb+NT1/iVLquEQzqkSA54D1xph/xL6kkzN/53z6NO9Jy71roVP96Cb5lsMJPcbBpo8RXzE/\nOv1HbDm0RYcGKqWOSzRb3GcD1wMjRGR1+HZhjOs6IXuK97CucB0jkjqBCcV9/HaVel0BgXLYOJux\nnceSnpDOa+tfi3dVSikbiWZUyRJjjBhj+hhj+oVvH9ZFccerspvk/LJyq387TucnOaYOZ0DTtrD2\nHbxOLxNOncAnuZ+w7dC2eFemlLKJBnXk5Pyd8+ma0pXMXV+E+7cT413SDzkc0Oty2DwXSvZzzWnX\nkOBK4Nmvn413ZUopm2gwwX2o/BCr9q1iRNshsPdriNP5t6PS/zoI+WHNNNIT0rmq+1V8uP1DdhXt\nindlSikbaDDBnZ2bTciEOF+SAQOdh8W7pOq1PN3qMln1IhjDDT1vwClOnv/m+XhXppSygQYT3PNz\n5tO6SWt67NsK7ibQbkC8Szq2gZOgcDPkLKNlUksuO+Uy3t3yLntL9sa7MqVUPdcggrvYV8yy3csY\n2XEksmMxdBoCTne8yzq2HuPBmwKrpgLwk14/AQNT106Na1lKqfrPFe8CasOi3EX4Qj5GZQyE/Q9Z\nfcj1nScJ+kyAL16CsY/QNrktl3S9hLc3vc3knpNp1aRVvCusNZm//aBWl7fjrxfV6vKUspsGscU9\nb+c8MhIz6HfkgDWhPvdvRxp4AwQr4CvrOpRT+kwhZEL8+8t/x7kwpVR9ZvvgLvWXsjh3MSM6jsCx\nYxEkpELr3vEuKzqte1vnCl81FYyhfdP2XNfjOmZuncnawrXxrk4pVU/Zvqtk6e6llAfLuaDTBbDi\nRutoSYcz3mVFL+tGeO9Wa1x39wuY0nsK7215j7999jemjpmKdcaBhkkIkcFh2sl+Wshh0qSIdIpI\nljIS8JGADwehb+f24aIcL2R/DZ4mkJgWvqVDckvr5kmGBtxmSkEDCO65O+aS5k1jgKc5HNoJZ90R\n75KOT58JkP1X+OQROGUUyZ5kbut3Gw+seIB5O+cxqtOoeFd40oQQnWUvfWQb3R25nCJ5dJU82ksB\nHgn+YP6AcVCOh3I8BHFQGcNuAiTgg+xZ1a/MnQRNW1tHpzZtbd2atYWmbcLh3gqaZEBCCkX+YrYf\n3s6ekj3sK9lHfmk+Rf4iSv2llAXKEAS3043H6SHVm0qrpFa0TGpJZrNMuqV1w+uM/wU6VONk6+Cu\nCFbwSe4njO08FlfOUmuiXfq3KzndcM4v4f27YOt86DaSy0+5nGkbpvHoykcZ1n6Y/QIi4IO8VbBj\nCeQsYY33U5qJdd5xv3Gy3bRmg+nInNAgck0L8kwL8k0qB01TDpJMGcc+o+OOh8aCrwjKDkHZQSgt\nhJICKN4HRfugaA8U7YW8lVC0FxMoJ9flZK3XyzqPh/VeN9vcHvJd3/9lloDQDCdJ4iARJwbwE6LC\nhDhIkBK+u8Cz00DnEPTzBRniCzC43E+KMeBwg8NlHbWbmAZJ6dYZIdM7Q/Ou0KI7ZJxmr1+Fqt6x\ndXAvy1tGaaDU2ipd/pz1B5JxarzLOn79roVFf4fsR6Dr+bgcLu4ZfA83zbmJ/3zxH3416FfxrrBm\nRftg88ewcTZsWwj+Umt6q97MDA5hjenKV6EubDVtrYtKnITM339UxdSm4Vs3IIAjMQ9n4g6cqTm4\nE3PAVQKAhByk+pLIKPVwuk/o6A/Rxh+kbTBA81BFeBvfj5MQBgjiIIQDPy4Oi5d9Tic7PUKON0Su\nN8jMhAreTnQizby0K0ukV1EyPQ4n0cL4SaGENNlEK/mMlnLou6YyiXwZ6sbK0KnMDQ1kvekIRNe9\noyNqFNg8uOfmzKWppymDWw6ErddC97H27N90eeCcX8AHd8O2bOg6nDPbnMmE7hN4ad1LDGs/jMFt\nBse7yh8qKYT178HX0yFnKWCgWXvoOxG6DLf2NySl84daHg74A+LHmbgTZ9J265a4E3H4AQhVtMBf\nfCrBsk4Ey9sTKm/FEVzsPJn1lUY+COJM3IUzeRM5Tb8mt1UBH7V0ECg+Df/BsQRLugIOkiink+zj\nNNnJAMdmBjo2c5drOr+Ut9kRasXs0GCmBYeTY1qfTGWqkbBtcJcHypm/cz5jOo/Bvfcr6yfzKSPj\nXdaJ6389LHrU6uvuch6IcHfW3azYs4J7l97LO5e+Q1NP05iXUdOYawchhjnWMNGZzfmOL3BLkC2h\ntswKXs6cUBbr8ztCvsASgOWxKVIqcCaFgzpxO87EXYgjgDFCqKI1/kODCJZ2JliaiQnGus2cBMsy\nCZZl4isYhcO7B3fKalwpq3A3XUeoogW+g2dRemgQ600n1ptOzAhZ10FtwWFGOlcxxvE5Nzo/5BbX\nLLKDfXkxeAELQ/2IditcNT62De5Pcj+hNFDKhZ0vhPVzQRzWVp5dubxWX/eHv4K1M6DX5SS5k3j4\nnIf58Uc/5qFPH+Lhcx6OW3npHOEa5wKudc2jrRxgv2nG88ExvBc8m3WmE7EMGXEdCm9R5+BM3IEj\nYQ8iIYxxECpvi//gWQTCQU0oKWZ1RFEpoYq2VOS3paLgAlxNv8aTvoyE1rPwtFiA/8AQfAfP+rbG\n/aTwenAErwdHkMFBrnEu5FrXPF5w/h9fhzJ5NDCB7FBfNMDV0Wwb3B9u+5CMxAyyWmXBrHug/SBr\nR5CdDZwMq1+FD39t7WRt0oI+GX24qc9NPLXmKfpl9OPq066u05JOkVxudH7IZc6leMXPomBv/hz8\nMfNDA/DH4uvjqMCZkIszIRdH4i6ciTtxuK1LnJqQm2BZB3z7zwtv5XaCUD3dcWtcBI70J3CkP87E\nHXhaLMTbci6e5ovwHRiK78A5EPpuJ2wBafw7eDlPBi9lnGMZd7mmM9XzNz4Lncqf/dfzjekSxzej\n6htbBvfhisMszlvMxNMm4iw9ALu/hOF/iHdZJ8/pgnFPwv+GwUe/gSutswXe3Odm1hWu48FPH6R5\nYnNGdop9l9AA2cTPXLMY5VxFmfHwdnAYLwRHs8W0r72VOEpxJuzGkbAbZ0IezoQ8xFOIiHUNzpAv\nnWBpZ3xlHQmWdSRU3gY7fmWDZZmU7ZqMw7sbT4v5eDPm40lfhq9wGL4DZ4PxfDtvABfTQ8OY6RvC\n1c6F3Ol6h/c89/FS8AIeDVwVx3eh6hP7/RVgXTDBH/JzUeeLrCF0YO/+7UitesC5v4GFD0LPy+H0\ni3E5XPz93L/z0zk/5Z5F9/DUqKcY1HpQ7a/bGM5yrOXnzhmc5VzHQZPMY/4reDF4AYc4mb5ig7gP\n4EzYg8O7G0fCHuu++7uRFiFdj4MYAAAS30lEQVR/CsHydoSO9CdY1p5QeXtMsMnJv6d6JFTRlvK8\n6/EV5uFtMRdvy49xpy/Ft/98/AcHEfnn6MfFK8FRzAwO4W7Xm9zgnMOFzk9hUzPoPjp+b0LVCxKL\nK4xnZWWZlStX1vpyK/3045+yt3Qvs8bPQqb/FLZ/Andvsq4uU0tq+8RIcBxDuYJ+eGa4NcTulsXW\nQSRYF4u4YfYNFJQW8OTIJ+nXsl/tFGaM9R/gJ/8Hu1aw16TxdOBipgWH1zimuoricXjzcSTkfbc1\n7d2DOCvCqxJCvgxC5W0IlbclWNGWUHkbTDC5dt6LjTgSc/C2/AhX0g5CvnQqCkYTONKbqs5E0Ue2\n8oj7GU537LSOtr3gL9aJylSDISKrjDFRXW/Rdlvc+aX5fLb3M27pewtiQlbgnDK6VkM77pxuGP9f\neG40vDQOJn0ATVqQmpDKUyOf4sY5NzL548n8bvDvmHDqhBNfjzGweY41kiVvFTRrzx/8k3kreC4V\neGp+PSEcnv1WX3RCLs7EXBzePYgjYC0+5CZU3gb/4f6EKtoSLG9DqKLV97oGGrNQWSfKcm7G2WQj\n3pazSWw3jWD6IiryxxAsPeV7835lujLO9wCbRqyEZf+B7YusrrQ2feJUvYon2wX37O2zMRjGdh4L\neV+EhwHa/7DwH2jdG370Brx6Jbw8Hm6YBYlptEluw7SLpnHPont4YMUDrCtcx28H/5YE13FsGYeC\nsO49WPIY7P0KUjvCJf+Cvj/ilT/Mrf51jlJrzHJiTvjfXYizHAAT9BAsb4f/4Jnh8dJtCfla0ADO\nYxZjQrDkNEq3d8fVbA3ejI9J6vQcgZJuVOSPJlTe4ds5fbitLe1uo2DGzfDcKLjoH9D/2jjWr+LB\nVsFtjGHWtln0aN6Dzimd4YvXrWGAXUfEu7TY6HwOTHwVpl0Dr1wJ10yD5JakeFN44vwneHz14zz7\n9bMsyVvCz/r+jHHdxuFyHOMjrSi2TiG7/HE4sA2ad4NLH7cOmPnBhScqt6Z3WkGdlIPTmw/w3Xjp\nI30JlnUgVNaBkC8DDemT4bBGoRT1xp26Ak+LhTTp/AT+Iz3xFVxAyBdxfvYu58LNi+HtydYJynI/\nh7GPWENKVaNgq+BeV7iODQc2cO8Z91oTNnzQMIYBHku3kXDVVHhrEjyeBaP+DP1/jNPh5M4BdzKk\n7RD++cU/uX/5/UxdO5Wrul/F8A7D6dDsuy018tfDyhdgzTSoOAJt+8OEl+C0i8HhxBjD3uI9bDiw\nAU/GnG+H44nLOkTQBBMIlnWi4nA/gmUdCZZ1AKMhERPGhf/gUPyHB+FJX4wnfTGupusIHOlLRUHE\nDvjkDLj+XVjwACz9J+Svg6tftaarBq/GnZMi8jxwMZBvjOkVzUJjtXPyT8v/xPtb32fBhAU0LcqH\n/wyA0Q/DWbfW+rriunOyKgWb4P1fQM4S60LDWTdavzSSMzDGsHDXQp5a8xTrD6wHoGuTtvQ0Htoe\n3EXbI/kk4SDUfhDBrsM5ktyC/NJ88kvzySnKYeuhrZT4rXN5GOMgVNGSYFkHawheWSft8ognZwme\n9EV40peBBBnX7WKm9JlCp2advptn7QyY8TMrtK95wxqZpGzneHZORhPcw4Bi4KV4Bnepv5Thbw5n\nZKeRPDj0QeukTAsegF+sg5R2tbouqIfBDdbOxDXTYN6foDh8UeHWfazTlTpcgLDr4BY+qdjHoiQv\n29xu8p0uTBUH3rnERUZSBu2S29EttRvdUrvRPb07lz22Q3ce1kPiLMLTPJumLVfiD/m5sPOF3NT7\nJrqkhg/MyfvC6lLzlcBVLzTM/T4NXK2OKjHGLBKRzJMt6mTN3jGb0kApV3a/0pqw7l1oPzgmoV1v\niUC/H0GfibB3DWyZB9sXQ+kBCAXAhOiQ1o3rWo/nula9IHMo/oQU9pbspSJYgUMcOMRBsieZ9IR0\nHFLFVrTZXffvS9XIBJtSkX8JvsLz8DRfzKzgx8za+gGBoh74CocTKm9Pa+7lWc+jnPbKBO4N3Mgb\nwZpPAaFnG7SnWuvjFpEpwBSAjh071tZivzV903S6pHShX0Y/KNwKe7+G0Q/V+npsweGw+qnb9odh\nvz7mrG74fn+3sjUrwC/EVzgMd/oyPGnLcDdbS6C4GwWF53F16R94wv0fHnE/Q1sp5LHAFei5Thqe\nWuu4NMY8bYzJMsZkZWTU7g6STQc38dX+r7jilCusS3mte9d6ose4Wl2PUnZhgsn4Ci6geMtvqcgf\ngyNhH0mdnsVkPseUhNFMC5zLna53+JvraZz88CpDyt5sMapk+qbpuB1uLul6iTVh7bvWaJKUWjxv\nhg3Foh9e2UwoAV/hefgOnI075Us86Yvwtn+dv/jSyT40lH+VLiJdjnC7/+fW9TpVg1DvhwoU+Yp4\nb+t7jOw0krSEtHA3yVfQY3y8S1Oq/jBu/IcGU7Ltl5TlXocJNmFZy52c3b4r29O28kzCg6RSFO8q\nVS2pcYtbRKYB5wEtRCQX+KMx5rlYF1bp7U1vU+Iv4YaeN1gTbNxNolvIKvYcBIp6ESjqiTNpG4Hm\nn/Dv5j6apJUz/vADzC64g71B3edhd9GMKrmmLgqpij/o55V1r3BG6zPo2bynNRzuq7esbpJU/fIp\nVT0hWNqVstKuOBLyaNZ8Ju+k7sCb8jitDg4kv3BsHVwdSMVKve4q+WD7B+SX5TO512Rrws4VULAe\nBvw4voUpZSOh8nbk5f2MVtsnMrTET3n6Kpp2ewRPxmwOVxyOd3nqBNTb4A6ZEFO/mUr3tO4MaTvE\nmrjyefCmQK8r4lucUja0uaI/K/Pu4YldfkaVFuFtns2Y6WP435r/UeovrXkBqt6ot8G9JG8JWw9v\nZVLPSdYQwJJCq3+770TwNKwT7CtVV3aZVtxZ8kdu2JfMm7l7GZzYmsdXP87FMy5m+qbpBEM6dNAO\n6mVwG2N47uvnaN2kNWM6j7Emrn4Vgj7Imhzf4pSyuUJSmOi7j/yK0/jXmgW83Ho0bZPbcv/y+7ly\n1pWs2LMi3iWqGtTL4F6ct5gv8r/gJ71+gtvhhlDI6ibpOARanh7v8pSyvVISuNH/K+h7Df2WP8PL\nvmb8Y9j/URYo46Y5N3HXwrvILcqNd5mqGvUuuAOhAI+ufJROzTp9d16S7dlwcDtk/SSutSnVkARw\nWVdaOvceZPWrjFr8X94b8xI/7/9zlu1exrh3x/Hf1f+lIlgR71LVUepdcL+z+R22Hd7GLwb+wtra\nBvj8OUhqDj0ujW9xSjU0IjD89zDuSdixBO+L47ipwyhmjZ/F+R3P58k1T3LZe5exOHdxvCtVEepV\ncJf4S3hi9RMMaDmAER3CV7XZ8xVseN86/7Re4UOp2Oh/LVw3HY7kwTMjaFWwmb+d+zeeueAZnOLk\n1vm3cnf23eSX5se7UkU9C+7nv3meA+UH+FXWr6yRJAAL/gIJqXDWbfEtTqmGrst5cNNCSEy3LlK9\n8gXObHMm71z6Dnf0v4PsXdmMe3cc0zZM09EncVZvgvuI7wgvr3uZsZ3H0jujtzVx5wrY/DEMvQsS\nU+NboFKNQfOu8NN5Voi/fxe8dzvuUIApfaYwY9wMerfozUOfPsT1H13PxgMb411to1VvgruZpxkv\njH6BXwz4hTXBGJj/Z0huBYOnxLc4pRqTxFT40Ztwzq/gy5fhuQvgwHY6NuvI/0b9j4fPeZi84jyu\nfv9q/rHyH3rwThzUm+AG6NmiJ22S21gPts6HnKXWhQL0gBul6pbDCeffZwX4oRx4+lxYOwMR4eIu\nFzNz/EzGdxvPC2tfYPx748nelR3vihuVehXc3wr6Yd79kNoRBtwQ72qUary6j4abF0HzU+CtSdZF\niSuKSPGmcP+Q+3lxzIs0cTfhjgV3cOeCO9ldrJe+qwv180IK2X+1Lk024WVw6YVrlYqVaE817OLn\n3OGawe2rp5H75TzuCUxhRajyavKT8DRfwvzQfObnLOaOATczqdckvE4dBRYr9W+LO2c5LPkH9LtO\nx20rVU8EcPFY4Cqu9t0HwOuev/CQ6xmaUQK48BWeR8nWuwkUn8bjqx9n/LvjmZ8zH2NMfAtvoOpX\ncJcfhhlTrC6SsX+NdzVKqaOsNKcx2vcI/wtcxNXObOZ4f8OljqWAwQRSKc+7lqdHPY3H6eGu7LuY\n/PFk1u5fG++yG5z6Fdwf/gYO58Hlz4BXT/KuVH1UjpeHA9cy3vcABSaFf3ueYLrnfvrIVgDOansW\n0y+dzn1n3sf2w9uZ+MFEfvXJr9h2aFucK2846k9wlx6AncutUSQdBse7GqVUDb42XRjn+wu/9k+h\no+xjpvc+Hnf/C/LX43K4mHDqBD647ANu6n0Ti3IXcdnMy/jd4t+x/fD2eJduexKLPqisrCyzcuXK\n439hRRG4EsDprvWajpdeH1Kp6CVTys2u95nsnE2yVFgXOxn6C2jdC4AD5QeY+s1Upm2YRkWwgvM6\nnMfkXpPp37J/nCuvP0RklTEmK6p561Vw1yMa3EodvzSO8OX5G+Czp8FfCp2HwZm3wikXgMNJYVkh\n0zZM4/WNr3O44jA9m/fkqu5XMbbzWJLcSfEuP640uGuBBrdSJ2bHXy+yuj6/eBE+e8Y6cVWzdtDn\nauj3I2hxCqX+Ut7b+h5vbnyTLYe2kORKYkznMYzOHM3g1oNxOernSOVY0uCuBRrcSp2YHX+96LsH\nQT9s+MC6gtWW+WCC0Lo3nHoRnHYhplVv1uz/irc3vc28nfMo8ZeQ5k1jRMcRDG03lDPbnEmyJzl+\nb6YOaXDXAg1upU7M94I7UtE++Pot6zTNO1cAxjoXUeZQyDyH8nYDWFqRz+ycOSzJW0KxvxiXuOiT\n0Yd+LfvRv2V/+mb0JS0hrU7fT12p9eAWkTHAvwAn8Kwx5piDrDW4lWq8qg3uSMUF1pk/t2XD9sVQ\nvNea7m4Cbfrgb9mDNcnNWBI8wuclO1l3eCuBUACAlokt6Z7ene5p3enUrBMdm3akQ9MOtEhsgdPh\njN0bi7HjCe4aO5JExAk8AYwCcoHPRWSmMWbdyZWplGq0kjOg/3XWzRgo3AK5K2HPatj9Je6v3yKr\n4giVKVbucPBNalu+SU5hY4WPjXtXsyJvGQFC3y7SiYMMbwotE5qT7k0hzZtGakIazTwpJHubkuxJ\nIcnThERXExLdTfC6EvC4vHgcHtxONy6HC7fDjVOcOMSBy+HCIQ4c4kCQ764RUA9EswdgMLDFGLMN\nQEReB8YBGtxKqZMnAi1OsW79rrGmGQPF+2D/JjiwjYTDuWQdziXrSB4U7YfifAKlhexxOdjldrPL\n5WKfy8k+5xHyXbvZ43CyzungoNOJvxYD12EMAlQuUSJvBpobYfZPvq619VUnmuBuB+yKeJwLnHH0\nTCIyBag8cXaxiDSUs6y3APbHu4h6Stvm2Bpl+8gjUc3WYNtGbjzh/yg6RTtjNMFdVRU/6Bg3xjwN\nPB3tiu1CRFZG2+/U2GjbHJu2T/W0bU5ONIe85wIdIh63B/Sku0opFSfRBPfnwCki0llEPMBEYGZs\ny1JKKVWdGrtKjDEBEbkd+BhrOODzxpjGdJ7GBtf9U4u0bY5N26d62jYnISYH4CillIqd+nNaV6WU\nUlHR4FZKKZvR4MY6pF9ENorIFhH5bRXPe0XkjfDzn4pIZt1XGT9RtM8vRWSdiHwlIvNFJOrxqA1B\nTe0TMd+VImJEpNEMg4umbURkQvj7s1ZEXqvrGm3JGNOob1g7XLcCXQAPsAbocdQ8twJPhe9PBN6I\nd931rH2GA0nh+z/T9vl++4TnawosAlYAWfGuu760DXAK8CWQFn7cMt512+GmW9wRh/QbY3xA5SH9\nkcYBL4bvvw2cL/XpxAWxVWP7GGMWGmNKww9XYI31byyi+f4APAD8DSivy+LiLJq2uQl4whhzEMAY\nk1/HNdqSBnfVh/S3q24eY0wAOAw0r5Pq4i+a9ol0I/BRTCuqX2psHxHpD3Qwxrxfl4XVA9F8d7oD\n3UVkqYisCJ+JVNWg8V1m4oeiOaQ/qsP+G6io37uIXAdkAefGtKL65ZjtIyIO4DFgUl0VVI9E891x\nYXWXnIf1S22xiPQyxhyKcW22plvc0R3S/+08IuICUoADdVJd/EV1ygMRGQncC1xqjKmoo9rqg5ra\npynQC8gWkR3AmcDMRrKDMtq/rfeMMX5jzHZgI1aQq2PQ4I7ukP6ZwA3h+1cCC0x4T0ojUGP7hLsC\n/ocV2o2tj/KY7WOMOWyMaWGMyTTGZGLtA7jUGGPvK41EJ5q/rXexdm4jIi2wuk621WmVNtTogzvc\nZ115SP964E1jzFoR+bOIXBqe7TmguYhsAX4JVDvkq6GJsn3+D0gG3hKR1SLSaM5lE2X7NEpRts3H\nQKGIrAMWAr82xhTGp2L70EPelVLKZhr9FrdSStmNBrdSStmMBrdSStmMBrdSStmMBrdSStmMBrdS\nStmMBrdSStnM/wd39FnE/iqt1gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x22712615ef0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# calculate a gaussian kde for the push-forward of the sampled priors\n",
    "push_forward_kde_pdf = gaussian_kde(q_sample_prior)\n",
    "\n",
    "# sanity check: Is push-forward KDE close to what we expect?\n",
    "q_values = np.linspace(-0.1,0.75,100) # plot pdf for these \"x-values\"\n",
    "\n",
    "plt.title(\"Push-Forward KDE from Sampled Priors\")\n",
    "plt.hist(q_sample_prior, density=True) # histogram of data\n",
    "plt.plot(q_values,push_forward_kde_pdf(q_values)) # plot of KDE\n",
    "\n",
    "plt.plot(q_values,data_kde_pdf(q_values)) # this is the data KDE compared to the push forward"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Use Accept-Reject to Generate Posterior Samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 235,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(598,)\n"
     ]
    }
   ],
   "source": [
    "# calculate maximum of the ratio\n",
    "M = np.max(data_kde_pdf(q_sample_prior)/push_forward_kde_pdf(q_sample_prior))\n",
    "\n",
    "# generate random numbers from uniform for accept-reject for each sample value\n",
    "test_value = np.random.uniform(0,1,np.shape(q_sample_prior))\n",
    "\n",
    "# calculate the ratio for accept reject: data_kde/push_kde/M and compare to test sample\n",
    "# is the kde ratio > test value?\n",
    "accept_or_reject_samples = np.greater(data_kde_pdf(q_sample_prior)/push_forward_kde_pdf(q_sample_prior)/M,\n",
    "                            test_value)\n",
    "\n",
    "# accepted values of posterior sample\n",
    "posterior_lam1_sample = lam1_sample_prior[accept_or_reject_samples]\n",
    "posterior_lam2_sample = lam2_sample_prior[accept_or_reject_samples]\n",
    "\n",
    "print(np.shape(posterior_lam1_sample))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We check to make sure the posterior samples are reasonable by plotting the posterior samples of $\\lambda_1$ and $\\lambda_2$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 236,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\ProgramData\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:8: DeprecationWarning: elementwise != comparison failed; this will raise an error in the future.\n",
      "  \n",
      "C:\\ProgramData\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:13: DeprecationWarning: elementwise != comparison failed; this will raise an error in the future.\n",
      "  del sys.path[0]\n",
      "C:\\ProgramData\\Anaconda3\\lib\\site-packages\\scipy\\optimize\\minpack.py:161: RuntimeWarning: The iteration is not making good progress, as measured by the \n",
      "  improvement from the last ten iterations.\n",
      "  warnings.warn(msg, RuntimeWarning)\n",
      "C:\\ProgramData\\Anaconda3\\lib\\site-packages\\scipy\\optimize\\minpack.py:161: RuntimeWarning: The iteration is not making good progress, as measured by the \n",
      "  improvement from the last five Jacobian evaluations.\n",
      "  warnings.warn(msg, RuntimeWarning)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(array([ 6.96784648,  0.87943693,  0.98091043,  1.18385741,  1.21768191,\n",
       "         1.48827789,  1.75887387,  0.98091043,  1.01473492,  0.43971847]),\n",
       " array([-0.00256269,  0.05656605,  0.11569479,  0.17482354,  0.23395228,\n",
       "         0.29308102,  0.35220976,  0.41133851,  0.47046725,  0.52959599,\n",
       "         0.58872473]),\n",
       " <a list of 10 Patch objects>)"
      ]
     },
     "execution_count": 236,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA00AAAEICAYAAACDCR/aAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3X24JGV95//3RxhEBcE4Y0RgHB9Q\ng14RcAIYoyFqsjwouMqVhUQU1zgb4xMb/Cm6iUE3mzXZ/Zmo+AuZgEswiriAZIyYxAcIEoU44IjC\nqBlxDCMoI8IMiE+j398fVYM9PX369Jnp0w/nvF/X1dep6rq7+lvVd92nvn3fVZ2qQpIkSZLU2wPG\nHYAkSZIkTTKTJkmSJEnqw6RJkiRJkvowaZIkSZKkPkyaJEmSJKkPkyZJkiRJ6mNRJk1JbkpyzLjj\nWIiSnJ3kb+dQ/oIkfzyfMU2yJG9Oct6449Cusz2ZP7Ync2N74vE4SaahPibZmOS5445jPiU5Jsmm\ncccxSZJcleR35vq6BZc09ToAkpye5Jrt81X15Kq6apb1rEhSSfacp1BHqnsfLARJXp1kbZIfJrlg\nlrKnJ/lJknuTbE2yLsnzdvP9d/sErar+pKrmfOBqNGxPelto7UmSByY5P8k3ktyT5PNJjutT3vZk\nDDwee1ss9bH9zL7Xbuc3k7wjyR7DWv+AMWxM8v02hu2PR40yhmnQ7qcVbT06vWvZMe1n+YY5rG/e\nvhBrj58L2ng39iu74JKmabFQGusxuw34Y+C9A5b/bFXtA+wPnA98KMnPzVdws9mdOmD9USfrw27b\nE7gV+FVgP+APadqHFX1eY3uinsa0PxdLfXxqu53PAX4LeMWuvs9ueH5V7dPxuG0uL57P+jElx/JL\nge+2f6fKokyaOr+tSnJk21uxNcm3k7yjLXZ1+/fu9puEpyd5QJI/aL+NvCPJhUn261jvS9pldyb5\nw673OTvJJUn+NslW4PT2vT+b5O4ktyc5J8leHeurJL+X5N/abz//e5LHta/ZmuRDneV3Y3+8LMn6\n9j1uSfJfOpYdk2RTkje023x7khckOT7JV5N8N8mbu1a5d5KL2/XdkOSpHes7vH3uniQXA3t3LHtY\nkr9PsjnJXe30QTPFXVWXVdXlwJ1z2d6q+ilNovUg4LHte78iyYZ2e9Zs/+YojT9vt31LkhuTPCXJ\nKuC3gTe09eMjbflHJbm03YavJ3ltx/b1qgM7DD9KcmKa4SV3p+k+/oWOZRuTvDHJjcD3pqRxXPBs\nT3baH1PXnlTV96rq7KraWFU/raq/B74OPG227bU9mSyL/XhcLPWxqr4MfBp4Ssf+fHzH+u7vmUiy\ntD3+7273waeTdJ7/HtZu+5Y0bc3e7II5bt8rtu/XdvmGJB/qmL81yWHt9Dvb+a1Jrk/yzI5yvfb7\ng9rtvyvJzcAv9Yn53CT/u+u5v0vy++30G9P06t2T5CtJnrMr+6Zj3Q8GTgZeBRySZGXX8l9J8pl2\nH96aphdopvrY7zOf0/nkwKpqQT2AjcBzu547HbimVxngs8Bp7fQ+wNHt9AqggD07XvefgQ00DdE+\nwGXA+9plhwL3Ar8C7AX8b+DHHe9zdjv/Appk9UE0/5CPpvmWcwWwHjij4/0KWAM8FHgy8EPgk+37\n7wfcDLy0o/zdwK/MsF922Addy04AHgeE5pvW+4Aj2mXHANuAtwBLaL7V2Qx8ANi3jesHwGO7tvPk\ntvzraU4+lrT75RvAf23nT27L/nH72ocDLwIe3K77/wKXD/CZ/zFwwSxl7t/+dn+/Drin3Y/PBr4D\nHAE8EHg3cHVb9j8A19N8exfgF4AD2mUXbI+9nX9AW/Yt7bY+FrgF+A996sDZwN+2y58AfA/49Xb/\nvIGmvu3VUW/XAQcDDxr3sbYYHtieLLr2pH3tz7dxPGm27cf2xONxgo7HhVwf2332+I7P5FvAy7uX\ndccP/E/g3PZ9lwDPBNLxvv8KPAr4ufZz+t251MFd2b52/93d7qsDaNqyb7ZlHwvcBTygnX8xTXu2\nJ3Bmu91799nvb6dJKH+ufb8vAZtm2J5n0fSyb98fDwO+3+6PJ7bLHtVx3DxuN4/h04DbgT2AjwDv\n6li2nKbentruw4cDh/WqjwN85n3bf+Aq4HfmHP+oG735frQV8962Mm5/3MfMjerVwFuBpV3rWcHO\njeongd/rmH9iW1n3pGlILupY9mDgR+zYqF49S+xnAB/uqhDP6Ji/Hnhjx/z/C/zFgPvldGY4yelR\n9nLgde30Me0BtEc7v28b11Fdcb2gYzuv7Vj2gPYAeWZ7cN62/eBsl3+m+0DoWHYYcNcA8Q6aNG1r\n68N3gGs7PpvzgT/rKLtP+7muoPmH81Waf34P6FrnBez4T+Uo4N+7yrwJ+D8z1QF2/Kfyh8CHuvbd\nN4FjOurtfx73MbaYHtie9DueFmp7sgT4BPBXs2y/7cmIH3g8Lur62O6zrTQJxddo/vc/oGPZTCfQ\nbwP+rnN5V315ccf8nwHnzqEOXr6r20eTkBwBnAKspknengS8DFjTJ4a7aIYpzrTfbwGO7ZhfxcxJ\nU4B/B57Vzr8C+FQ7/XjgDuC5wJIhHcOfoK3XNMnR5u3rbuvTh2d43Q71cbbPvMfrd2j/2cWkaaEO\nz3tBVe2//QH8Xp+yL6f5huDLST6X/hdPPorm24DtvkHToP58u+zW7Quq6j52HjZ2a+dMkie0XYbf\nartV/wRY2vWab3dMf7/H/D594h1IkuOSXNt2W98NHN8Vx51V9ZOO9+wVV2ccnfvhp8Ammv3zKJpv\nUqqj7P37M8mDk/xVmiEQW2n+4e2f4V3oeW1bJ5ZW1dFV9Yn2+R0+16q6l+azO7CqPgWcA7wH+HaS\n1UkeOsP6Hw08qu1Wvrvdl2+mqR/b3dr7pT3j+Glb/sABX6/5YXsyB9PcnrRDdt5Hc0L86lk21fZk\nPDwee1ss9fGIqnpYVT2uqv6gXc9s/hdNr88/pRkyfFbX8m91TN9Hu9+TfCw/u9nDb3eU6ayDL9iN\n7ftnmi+SntVOX0XTO/+r7TxtHGemGfK8pd3v+7FjXepe76O6nvsGM2jbzw/SJDDQXCf2/nbZBppk\n/2zgjiQfzG7c9CLJwcCvbV8/TSK7N83oBGh6xb62q+vveq95OZ9cqEnTwKrq36rqVOARwJ8ClyR5\nCE0G2+02moZju+U03+58m+bbz/vHSyZ5EE334A5v1zX/l8CXgUOq6qE0DVB2fWvmLskDgUtphhv8\nfPtP6IrdjOPgjvU/gGa/3Eazjw5M0rnu5R3TZ9J8u3dUuz+etX01uxHLIHb4XNvP/+E03xJRVe+q\nqqfRDKF4AvD/tEW7P89bga93/kOvqn2r6viOMr3q1UxxhGZffnPA12vMbE+mtz1p13M+zUngi6rq\nx7sYr+3JhFjsx2NrsdTH+2h6ALd75P0rrbqnqs6sqscCzwd+f5Brc6rquPrZzR7eP0vxXdm+7UnT\nM9vpf6YraUpz/dIbgd8EHta2qVvYsS51r/d2OtpNdmwXe7kIODnJo2l6FC+9f8VVH6iqX2m3rWiO\no111Gk3e8ZEk36LpEdsbeEm7/Faaod299KobM37mzNP55KJPmpK8OMmy9luBu9unf0LTZfhT2gsp\nWxcB/zXJY5LsQ/PN0cVVtQ24BHh+kl9Oc/HmW5n9w9mXpqv53iRPAl45tA3rLUn27nzQjE1+IM32\nbktzm93f2M33eVqSF6a5kPMMmrHa19KML98GvDbJnkleCBzZ8bp9ab5duzvNXX/+aJaN2bPdhj2A\nPdpt2pWLmT8AvCzJYe1J358A11XVxiS/lOSoJEtoxiv/gKZ+QPPPtLN+/CuwNc2Fkw9Kskeai2pn\nvAizy4eAE5I8p32/M2n23Wd2YZs0BrYn09ue0Jzk/gLNnbG+P0vZfmxPJsQiOx5nsljq4zrgt9q4\njqVJPABI8rwkj28Tma002/iTGdazq3Zl+/6ZpuflQVW1ieY6pGNpktrPt2X2pWnnNgN7JnkLzXVw\ns8XypjQ3QzgIeE2/wlX1+Xb95wH/WFV3AyR5YpJnt/XmBzTt6e7st5fQHDuHdTxeRLPfHk7TA/Xc\nJL/ZtukPT3szDHauj9DnM2fu7f9AFn3SRFNBb0pyL/BO4JSq+kHb/f4/gH9J0xV9NM0dad5H0833\ndZpK9BqAqrqpnf4gTZZ/D81Y0B/2ee/X03SF3gP8NXDx7mxImi7kZ/Yp8ss0laj78Vqag+yuNp41\nuxMHTZfrf2rXdxrwwqr6cVX9CHghzfjru9oyl3W87i9oLmLcPib7H2Z5nz9o4z+L5kLJ77fPzUlV\nfZJmPPKlNJ/d42jGGEPTOP11G+83aIY1bL/TzPnAoW39uLwdcvR8mobg6+12nEfTlT5IHF9pt+Pd\n7WufT3MC96O5bpPGxvZkCtuT9hvW/0Jz7H4rvYfkDMT2ZKIspuOxp0VUH1/XrvNumjutXd6x7BCa\na2nupfmy5f+rWX7La652Zfuq6qttTJ9u57fS9L78S8cQ5n8EPkZz7dk3aOrlbMMY39qW/TrwTzT1\nejYX0Vy79IGO5x5Ic1OJ79AMYXwETY8pSX47yU0DrJe2/NE019G9p6q+1fFYQzN08tSq+nea4dxn\n0tySfB2w/W6pO9TH9rl+n/lczycH244dh4NrWNpvqu6m6Zr/+rjjkTS9bE+kyeHxKC1O9jQNUZLn\np7n47CE039Z8keaOKZI0J7Yn0uTweJRk0jRcJ9FcDHgbTXfwKWVXnqRdY3siTQ6PR2mRc3ieJEmS\nJPVhT5MkSZIk9bErt2ceiqVLl9aKFSvG9faSerj++uu/U1XLxh3HXNmeSJPH9kTSMExKWzK2pGnF\nihWsXbt2XG8vqYckM/5y+CSzPZEmj+2JpGGYlLbE4XmSJEmS1IdJkyRJkiT1YdIkSZIkSX2YNEka\nmSQHJ7kyyfokNyV5XY8yxyTZkmRd+3jLOGKVJEnabmw3gpC0KG0DzqyqG5LsC1yf5ONVdXNXuU9X\n1fPGEJ8kSdJO7GmSNDJVdXtV3dBO3wOsBw4cb1SSFqokT+zotV6XZGuSM8Ydl6TpY9IkaSySrAAO\nB67rsfjpSb6Q5GNJnjzD61clWZtk7ebNm+cxUknTqqq+UlWHVdVhwNOA+4APjzksSVPIpEnSyCXZ\nB7gUOKOqtnYtvgF4dFU9FXg3cHmvdVTV6qpaWVUrly0b+2/eSZp8zwG+VlUT8ZsvkqaLSZOkkUqy\nhCZhen9VXda9vKq2VtW97fQVwJIkS0ccpqSF5xTgol4L7LmWNBtvBCHthhVnfXTo69z49hOGvs5J\nkSTA+cD6qnrHDGUeCXy7qirJkTRf7tw5jPf389Kksm7OryR7AScCb+q1vKpWA6sBVq5cWQOt9Oz9\nhhVexzq3DH+dkobCpEnSKD0DOA34YpJ17XNvBpYDVNW5wMnAK5NsA74PnFJVg53ESFJvxwE3VNW3\nxx2IpOlk0iRpZKrqGiCzlDkHOGc0EUlaJE5lhqF5kjQIr2mSJEkLVpIHA78O7HQNpSQNyp4mSZK0\nYFXVfcDDxx2HpOlmT5MkSZIk9WHSJEmSJEl9mDRJkiRJUh8mTZIkSZLUx6xJU5KDk1yZZH2Sm5K8\nrkeZY5JsSbKufbxlfsKVJEmSpNEa5O5524Azq+qGJPsC1yf5eFXd3FXu01X1vOGHKEmSJEnjM2tP\nU1XdXlU3tNP3AOuBA+c7MEmSJEmaBHO6pinJCuBw4Loei5+e5AtJPpbkyTO8flWStUnWbt68ec7B\nSpIkSdKoDZw0JdkHuBQ4o6q2di2+AXh0VT0VeDdwea91VNXqqlpZVSuXLVu2qzFLkiRJ0sgMlDQl\nWUKTML2/qi7rXl5VW6vq3nb6CmBJkqVDjVSSJEmSxmCQu+cFOB9YX1XvmKHMI9tyJDmyXe+dwwxU\nkiRJksZhkLvnPQM4DfhiknXtc28GlgNU1bnAycArk2wDvg+cUlU1D/FKkiRJ0kjNmjRV1TVAZilz\nDnDOsIKSJGlYVpz10aGvc+PbTxj6OiVJk2tOd8+TJEmSpMXGpEmSJEmS+hjkmiZJkjRlHJYoScNj\nT5MkSZIk9WHSJEmSJEl9mDRJkiRJUh8mTZIkSZLUh0mTJElasJLsn+SSJF9Osj7J08cdk6Tp493z\nJEnSQvZO4B+q6uQkewEPHndAkqaPSZMkSVqQkjwUeBZwOkBV/Qj40ThjkjSdHJ4nSZIWqscCm4H/\nk+TzSc5L8pDuQklWJVmbZO3mzZtHH6WkiWdPkzRh/EFKSRqaPYEjgNdU1XVJ3gmcBfxhZ6GqWg2s\nBli5cmWNPEpJE8+eJkmStFBtAjZV1XXt/CU0SZQkzYlJkyRJWpCq6lvArUme2D71HODmMYYkaUo5\nPE+SJC1krwHe39457xbgZWOOR9IUMmmSJEkLVlWtA1aOOw5J082kSZI0MebjRiiSJO0ur2mSJEmS\npD7saRoSbxMtzS7JwcCFwCOBnwKrq+qdXWUCvBM4HrgPOL2qbhh1rJIkSduZNEkapW3AmVV1Q5J9\ngeuTfLyqOu9mdRxwSPs4CvjL9q8kSdJYODxP0shU1e3be42q6h5gPXBgV7GTgAurcS2wf5IDRhyq\nJEnS/expkjQWSVYAhwPXdS06ELi1Y35T+9ztXa9fBawCWL58+XyFuWA4hFiSpF1nT5OkkUuyD3Ap\ncEZVbe1e3OMltdMTVauramVVrVy2bNl8hClJkgSYNEkasSRLaBKm91fVZT2KbAIO7pg/CLhtFLFJ\nkiT1YtIkaWTaO+OdD6yvqnfMUGwN8JI0jga2VNXtM5SVJEmad17TJGmUngGcBnwxybr2uTcDywGq\n6lzgCprbjW+gueX4y8YQpyRJ0v1MmiSNTFVdQ+9rljrLFPCq0UQkSZI0O4fnSZIkSVIfJk2SJEmS\n1IfD87RozMfv1EiSJGnhm7WnKcnBSa5Msj7JTUle16NMkrwryYYkNyY5Yn7ClSRJkqTRGqSnaRtw\nZlXdkGRf4PokH6+qmzvKHAcc0j6OAv6y/StJkiRJU23Wnqaqur2qbmin7wHWAwd2FTsJuLAa1wL7\nJzlg6NFKkiRJ0ojN6ZqmJCuAw4HruhYdCNzaMb+pfW6HH6RMsgpYBbB8+fK5RSpJmiheJyhJWiwG\nTpqS7ANcCpxRVVu7F/d4Se30RNVqYDXAypUrd1oubefJmCRJkibFQLccT7KEJmF6f1Vd1qPIJuDg\njvmDgNt2PzxJkiRJGq9B7p4X4HxgfVW9Y4Zia4CXtHfROxrYUlW3z1BWkiRJkqbGIMPzngGcBnwx\nybr2uTcDywGq6lzgCuB4YANwH/Cy4YcqSZI0N0k2AvcAPwG2VdXK8UYkaRrNmjRV1TX0vmaps0wB\nrxpWUJIkSUP0a1X1nXEHIWl6zenueQvJNNxoYD5i3Pj2E4a+TkmSJGkhG+hGEJIkSVOqgH9Kcn37\n0yc7SbIqydokazdv3jzi8CRNg0Xb0yRJ0q6ahtEKut8zquq2JI8APp7ky1V1dWcBfxJF0mymImny\nn5MkSdoVVXVb+/eOJB8GjgSu7v8qSdqRw/MkSdKClOQhSfbdPg38BvCl8UYlaRpNRU+TJEnSLvh5\n4MPNT06yJ/CBqvqH8YYkaRqZNC0y3pFPkrRYVNUtwFPHHYek6efwPEmSJEnqw6RJkiRJkvowaZIk\nSZKkPkyaJEmSJKkPkyZJkiRJ6sOkSZIkSZL6MGmSJEmSpD78nSZJmjDz8XtqkiRp19nTJEmSJEl9\nmDRJkiRJUh8mTZIkSZLUh0mTJEmSJPVh0iRJkiRJfZg0SZIkSVIfJk2SRibJe5PckeRLMyw/JsmW\nJOvax1tGHaMkSVI3f6dJ0ihdAJwDXNinzKer6nmjCUeSJGl29jRJGpmquhr47rjjkCRJmguTJkmT\n5ulJvpDkY0mePFOhJKuSrE2ydvPmzaOMT5IkLTIOz5M0SW4AHl1V9yY5HrgcOKRXwapaDawGWLly\nZY0uxB2tOOuj43prSZI0IvY0SZoYVbW1qu5tp68AliRZOuawJEnSImfSJGliJHlkkrTTR9K0UXeO\nNypJkrTYOTxP0sgkuQg4BliaZBPwR8ASgKo6FzgZeGWSbcD3gVOqamxD7yQtDEn2ANYC3/TunJJ2\nhUmTpJGpqlNnWX4OzS3JJWmYXgesBx467kAkTadZh+f5Y5SSJGlaJTkIOAE4b9yxSJpeg1zTdAFw\n7CxlPl1Vh7WPt+1+WJIkSUPxF8AbgJ/OVMCfMJA0m1mTJn+MUpIkTaMkzwPuqKrr+5WrqtVVtbKq\nVi5btmxE0UmaJsO6e54/RilJkibNM4ATk2wEPgg8O8nfjjckSdNoGEnT9h+jfCrwbpofo+zJb3Ik\nSdKoVNWbquqgqloBnAJ8qqpePOawJE2h3U6a/DFKSZIkSQvZbidN/hilJEmadFV1lb/RJGlXzfo7\nTf4YpSRJkqTFbNakyR+jlCRJkrSYDevueZIkSZK0IJk0SZIkSVIfJk2SJEmS1Mes1zRJs1lx1kfH\nHYIkSZI0b+xpkiRJkqQ+TJokSZIkqQ+TJkmSJEnqw6RJkiRJkvowaZIkSZKkPkyaJEmSJKkPkyZJ\nkiRJ6sOkSZIkSZL6MGmSJEmSpD5MmiRJkiSpD5MmSZIkSerDpEmSJEmS+jBpkiRJkqQ+TJokSZIk\nqQ+TJkmStCAl2TvJvyb5QpKbkrx13DFJmk57jjsASZKkefJD4NlVdW+SJcA1ST5WVdeOOzBJ08Wk\nSZIkLUhVVcC97eyS9lHji0jStHJ4niRJWrCS7JFkHXAH8PGqum7cMUmaPiZNkiRpwaqqn1TVYcBB\nwJFJntJdJsmqJGuTrN28efPog5Q08UyaJEnSgldVdwNXAcf2WLa6qlZW1cply5aNPDZJk8+kSdLI\nJHlvkjuSfGmG5UnyriQbktyY5IhRxyhp4UiyLMn+7fSDgOcCXx5vVJKmkUmTpFG6gB7f8nY4Djik\nfawC/nIEMUlauA4ArkxyI/A5mmua/n7MMUmaQt49T9LIVNXVSVb0KXIScGF7x6trk+yf5ICqun0k\nAUpaUKrqRuDwccchafrZ0yRpkhwI3Noxv6l9bideuC1JkkbFpEnSJEmP53r+pooXbkuSpFExaZI0\nSTYBB3fMHwTcNqZYJEmSgAGSJu92JWmE1gAvaduVo4EtXs8kSZLGbZAbQVwAnANcOMPyzrtdHUVz\nt6ujhhGcpIUlyUXAMcDSJJuAPwKWAFTVucAVwPHABuA+4GXjiVSSJOlnZk2avNuVpGGpqlNnWV7A\nq0YUjiRJ0kCGcU2Td7uSJEmStGANI2nybleSJEmSFqxhJE3e7UqSJEnSgjWMpMm7XUmSJElasGa9\nEYR3u5IkSZK0mA1y9zzvdiVJkiRp0RrG8DxJkiRJWrBMmiRJkiSpD5MmSZIkSerDpEmSJEmS+jBp\nkiRJkqQ+TJokSZIkqQ+TJkmSJEnqw6RJkiRJkvowaZIkSZKkPkyaJEnSgpTk4CRXJlmf5KYkrxt3\nTJKm057jDkCSJGmebAPOrKobkuwLXJ/k41V187gDkzRd7GmSJEkLUlXdXlU3tNP3AOuBA8cblaRp\nZNIkSZIWvCQrgMOB63osW5VkbZK1mzdvHnVokqaASZMkSVrQkuwDXAqcUVVbu5dX1eqqWllVK5ct\nWzb6ACVNPJMmSZK0YCVZQpMwvb+qLht3PJKmk0mTJElakJIEOB9YX1XvGHc8kqaXSZMkSVqongGc\nBjw7ybr2cfy4g5I0fbzluCRJWpCq6hog445D0vSzp0mSJEmS+jBpkiRJkqQ+HJ4nSZIkDeLs/Ya8\nvi3DXZ/mjUmTJEmSZrXirI8OdX0b337CUNc3bL22d+PeYwhEE8HheZIkSZLUh0mTJEmSJPVh0iRJ\nkiRJfZg0SZIkSVIfJk2SJEmS1IdJk6SRSnJskq8k2ZDkrB7LT0+yOcm69vE744hTkiRpO285Lmlk\nkuwBvAf4dWAT8Lkka6rq5q6iF1fVq0ceoCRpanXfInzj3r+1W+vz9uLqZE+TpFE6EthQVbdU1Y+A\nDwInjTkmSZKkvgZKmhxOI2lIDgRu7Zjf1D7X7UVJbkxySZKDe60oyaoka5Os3bx583zEKkmSBAyQ\nNHUMpzkOOBQ4NcmhPYpeXFWHtY/zhhynpIUhPZ6rrvmPACuq6heBTwB/02tFVbW6qlZW1cply5YN\nOUxJkqSfGeSapvuH0wAk2T6cpvsaBEmazSags+foIOC2zgJVdWfH7F8DfzqCuCRJI9Z9DZI0yQYZ\nnudwGknD8jngkCSPSbIXcAqwprNAkgM6Zk8E1o8wPkmSpJ0MkjQ5nEbSUFTVNuDVwD/SJEMfqqqb\nkrwtyYltsdcmuSnJF4DXAqePJ1pJkqTGIMPzHE4jaWiq6grgiq7n3tIx/SbgTaOOS5IkaSaDJE33\nD6cBvkkznGaHG98nOaCqbm9nHU4jSZKkge3ubypJ823WpKmqtiXZPpxmD+C924fTAGurag3NcJoT\ngW3Ad3E4jSRJ0oJkgqPFaJCeJofTSJIkSVq0BkqaJEmSpk2S9wLPA+6oqqeMOx6p2+7edn3j208Y\nUiSazSB3z5MkSZpGFwDHjjsISdPPpEmSJC1IVXU1zbXWkrRbHJ4nSZIWtSSrgFUAy5cvH3M00m46\ne795WOeW4a9zypg0SZKkRa2qVgOrAVauXFljDmdoVpz10aHe6W7j3kNblVq7/fmcPZQwNACH50mS\nJElSHyZNkiRJktSHw/MkSdKClOQi4BhgaZJNwB9V1fnjjaqPIV+L4nA6aXhMmiRJ0oJUVaeOOwZJ\nC4PD8yRJkiSpD5MmSZIkSerDpEmSJEmS+vCaJkmSJEkzWnHWR3fpdRvffsKQIxkfe5okSZIkqQ+T\nJkmSJEnqw6RJkiRJkvowaZIkSZKkPkyaJEmSJKkPkyZJkiRJ6sOkSZIkSZL68HeaJEmSJM1o496/\ntWsvPLvfsi27ts4xsadJkiRJkvowaZIkSZKkPkyaJEmSJKkPkyZJkiRJ6sOkSZIkSZL6MGmSJEmS\npD5MmiRJkiSpD5MmSZIkSerDpEmSJEmS+hgoaUpybJKvJNmQ5Kweyx+Y5OJ2+XVJVgw7UEkLg+2J\npFGarc2RpEHMmjQl2QN4D3A2cTqfAAAINklEQVQccChwapJDu4q9HLirqh4P/Dnwp8MOVNL0sz2R\nNEoDtjmSNKtBepqOBDZU1S1V9SPgg8BJXWVOAv6mnb4EeE6SDC9MSQuE7YmkURqkzZGkWe05QJkD\ngVs75jcBR81Upqq2JdkCPBz4TmehJKuAVe3svUm+MsD7L+1ez4QyzuEyziHKnw4c56PnOZRxtyfz\nYSrqAMY5bNMSJwwx1syt33e+25NBDNLm7E57Mmn1YJLiMZbeJikWGGc8b93p+9CZYpmEtmSgpKnX\nN7y1C2WoqtXA6gHe82crTtZW1cq5vGYcjHO4jHO4JijOsbYn82GC9m1fxjlc0xInTFes82Be25NJ\n27eTFI+x9DZJscBkxTNJsfQyyPC8TcDBHfMHAbfNVCbJnsB+wHeHEaCkBcX2RNIoDdLmSNKsBkma\nPgcckuQxSfYCTgHWdJVZA7y0nT4Z+FRV7fRNjqRFz/ZE0igN0uZI0qxmHZ7XXlPwauAfgT2A91bV\nTUneBqytqjXA+cD7kmyg+Ub4lCHGOPbhNwMyzuEyzuGaiDgnoD2ZDxOxbwdgnMM1LXHCdMU6VDO1\nOUN8i0nbt5MUj7H0NkmxwGTFM0mx7CR+gStJkiRJMxvox20lSZIkabEyaZIkSZKkPiYiaUpybJKv\nJNmQ5Kwey09PsjnJuvbxO5MYZ1vmN5PcnOSmJB8YdYxtDLPtzz/v2JdfTXL3OOJsY5kt1uVJrkzy\n+SQ3Jjl+QuN8dJJPtjFeleSgMcT43iR3JPnSDMuT5F3tNtyY5IhRxzjNrKtDj3Mq6usAcT4pyWeT\n/DDJ60cdX0ccs8X52+1+vDHJZ5I8ddQxTrMBjqsHJrm4XX5dkhVjjOVZSW5Isi3JyfMVx4Cx/H57\nTnRj2+7M6+/tDBDP7yb5Ynv+c02SQ8cVS0e5k5NUknm71fYA+2Wk59mD7JtMwPn0TqpqrA+aCzO/\nBjwW2Av4AnBoV5nTgXOmIM5DgM8DD2vnHzGJcXaVfw3NhbGTuk9XA69spw8FNk5onP8XeGk7/Wzg\nfWOI81nAEcCXZlh+PPAxmt8tORq4bhyf+zQ+rKvzEutU1NcB4nwE8EvA/wBeP44YB4zzlzv+Nx3n\n8T+nfTvIcfV7wLnt9CnAxWOMZQXwi8CFwMlj3i+/Bjy4nX7lfO2XOcTz0I7pE4F/GFcsbbl9gauB\na4GVY9wvpzOi8+wB4xn7+XSvxyT0NB0JbKiqW6rqR8AHgZPGHFMvg8T5CuA9VXUXQFXdMeIYYe77\n81TgopFEtrNBYi3goe30fozn9zUGifNQ4JPt9JU9ls+7qrqa/r9ndBJwYTWuBfZPcsBoopt61tUh\nm5b6OlucVXVHVX0O+PHoouoZx2xxfmb7/yaaE7Sx9DBOqUGOq5OAv2mnLwGek6TXD+vOeyxVtbGq\nbgR+Og/vP9dYrqyq+9rZ+a53g8SztWP2IfT4oeNRxdL678CfAT+YpzjmEsuoTMv59E4mIWk6ELi1\nY35T+1y3F7Xdu5ckObjH8vk2SJxPAJ6Q5F+SXJvk2JFF9zOD7k/abvLHAJ8aQVy9DBLr2cCLk2wC\nrqDpGRu1QeL8AvCidvo/AvsmefgIYpuLgeuGdmJdHT3r6/x5OU0vngYzSF28v0xVbQO2APNxXE3S\ncTHXWOa73g0UT5JXJfkaTbLy2nHFkuRw4OCq+vt5imHgWFqjOs+elvPpnUxC0tTrm5juzP8jwIqq\n+kXgE/zs25xRGiTOPWm6FI+h6cE5L8n+8xxXt0Hi3O4U4JKq+sk8xtPPILGeClxQVQfRDNd5X5JR\n19tB4nw98KtJPg/8KvBNYNt8BzZHc6kb2pF1dfSsr/Mgya/RnLy+cdyxTJFB6uKo6uskHRcDx5Lk\nxcBK4H+NO56qek9VPY7mGPiDccTS/m/4c+DMeXr/gWNpjfI8e1rOp3cyCUnTJqAzoz2IrmEtVXVn\nVf2wnf1r4Gkjiq3TrHG2Zf6uqn5cVV8HvkLzoY/SIHFudwrjG5oHg8X6cuBDAFX1WWBvYOlIovuZ\nQerobVX1wqo6HPhv7XNbRhfiQOZSN7Qj6+roWV+HLMkvAucBJ1XVneOOZ4oM+v//YIAke9IM0e03\n/HQ+YxmVgWJJ8lyatubEjnO5scXT4YPAC8YUy77AU4CrkmykuW5zzTzdDGLSzrOn5Xx6J5OQNH0O\nOCTJY5LsRXMiv6azQNc49hOB9SOMb7tZ4wQup7nokSRLaboXbxlplIPFSZInAg8DPjvi+DoNEuu/\nA88BSPILNCeim0ca5WB1dGlHr8KbgPeOOMZBrAFeksbRwJaqun3cQU0J6+roWV+HKMly4DLgtKr6\n6rjjmTKDHP9rgJe20ycDn6r2CvYxxDIqg7Q3hwN/RZMwzfd1KYPE03nifQLwb+OIpaq2VNXSqlpR\nVStorvc6sarWjjoWGPl59rScT+9s3HeiaNuU44Gv0txN47+1z72NpgIB/E/gJpqx+FcCT5rQOAO8\nA7gZ+CJwyiTG2c6fDbx9Cj77Q4F/aT/7dcBvTGicJ9M0vl+l+Sb3gWOI8SLgdpoL0jfR9Hz8LvC7\nHfXzPe02fJF5ulPPQn1YV4ce51TU1wHifGT7/Fbg7nb6oRMY53nAXW3dXAesHcf+nNbHAMfV3jR3\nptwA/Cvw2DHG8kttHfgecCdw0xhj+QTw7Y56t2bMn9M7ac4n19GcTz55XLF0lb1qPtu4AfbLSM+z\nB4hnIs6nux9pg5MkSZIk9TAJw/MkSZIkaWKZNEmSJElSHyZNkiRJktSHSZMkSZIk9WHSJEmSJEl9\nmDRJkiRJUh8mTZIkSZLUx/8Pxv9WzGiTn0QAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x227128deeb8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# this is the push forward of the posterior samples\n",
    "push_forward_posterior_samples = Q(posterior_lam1_sample,posterior_lam2_sample)\n",
    "\n",
    "# check that the simulation makes sense by plotting!\n",
    "f, (ax1, ax2, ax3) = plt.subplots(1, 3)\n",
    "f.set_figwidth(14)\n",
    "# Comparison of exact pdf vs. density\n",
    "ax1.set_title('Histogram: Lambda 1 Posterior')\n",
    "ax1.hist(posterior_lam1_sample, density=True)\n",
    "ax2.set_title('Histogram: Lambda 2 Posterior')\n",
    "ax2.hist(posterior_lam2_sample, density=True)\n",
    "\n",
    "ax3.set_title('Histogram: Posterior Push-Forward vs. \"Actual\"')\n",
    "ax3.hist(push_forward_posterior_samples, density=True)\n",
    "ax3.hist(q_data, density=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The graphs above show that our posterior sample is doing a decent job at reproducing our original \"true\" distributions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## How sensitive is this kind of modeling to different priors?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}