mutex_bayes.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# This notebook will compute the bayes factors of mutual exclusivity and co-occurrence under the model:\n",
    "\n",
    "\\begin{equation}\n",
    "\\begin{split}\n",
    "K = \\frac{P(D | M_{A,B})}{P(D) | H_0} = \\frac{\\int \\int \\int P(\\theta_A, \\theta_B, \\theta_{A,B}^M | M_{A,B}) P(D| \\theta_A, \\theta_B, \\theta_{A,B}^M , M_{A,B}) d\\theta_{A,B}^M d\\theta_A d\\theta_b }{\\int \\int \\int P(\\theta_A, \\theta_B, \\theta_{A,B}^M | H_0) P(D| \\theta_A, \\theta_B, \\theta_{A,B}^M , H_0) d\\theta_{A,B}^M d\\theta_A d\\theta_b }\n",
    "\\end{split}\n",
    "\\end{equation}\n",
    "\n",
    "where we assume uniform priors on $\\theta_A $ and $\\theta_B$:\n",
    "$ P(\\theta_A \\in dp | M_{A,B}) = dp$ (since  $\\int_0^1 P(\\theta_A \\in dp | M_{A,B}) dp = 1$)\n",
    "\n",
    "and\n",
    "$ P(\\theta_B \\in dp | M_{A,B}) = dp$ (since  $\\int_0^1 P(\\theta_B \\in dp | M_{A,B}) dp = 1$)\n",
    "\n",
    "Then, because under $M_{A,B}$, $\\theta_{A,B}^M$ must be greater than the value inferred by independence:\n",
    "\n",
    "\\[ P(\\theta_{A,B}^M \\in dp | \\theta_A, \\theta_B, M_{A,B}) = \\begin{cases}\n",
    "0 & p \\leq \\theta_A (1 - \\theta_B) + \\theta_B (1 - \\theta_A) \\\\\n",
    "\\frac{dp}{\\int_{\\theta_A (1 - \\theta_B) + \\theta_B (1 - \\theta_A)}^1 dp^*} & p > \\theta_A (1 - \\theta_B) + \\theta_B (1 - \\theta_A)\\\\\n",
    "\\end{cases}\n",
    "\\]\n",
    "\n",
    "Now, the likelihoods:\n",
    "\n",
    "\\begin{equation}\n",
    "P(D | \\theta_A, \\theta_B, \\theta_{A,B}^M, M_{A,B}) = P(X_A | \\theta_A, M_{A,B}) P(X_B | \\theta_B, M_{A,B}) P(X_{A,B}^M | \\theta_A, \\theta_B, \\theta_{A,B}^M, M_{A,B})\n",
    "\\end{equation}\n",
    "\n",
    "All the $X_E \\sim Binom(\\theta_E)$ for event $E$.\n",
    "\n",
    "\n",
    "So, $P(X_A | \\theta_A, M_{A,B}) = {n \\choose X_A} (\\theta_A)^{X_A} (1 - \\theta_A)^{n - X_A}$\n",
    "\n",
    "$P(X_B | \\theta_B, M_{A,B}) = {n \\choose X_B} (\\theta_B)^{X_B} (1 - \\theta_B)^{n - X_B}$\n",
    "\n",
    "$P(X_{A,B}^M | \\theta_{A,B}^M, M_{A,B}) = {n \\choose X_{A,B}^M} (\\theta_{A,B}^M)^{X_{A,B}^M} (1 - \\theta_{A,B}^M)^{n - X_{A,B}^M}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from scipy import special\n",
    "# Let's calculate the top integral\n",
    "\n",
    "\n",
    "\n",
    "n = 100\n",
    "x_A = 20\n",
    "x_B = 30\n",
    "x_AB = 15\n",
    "\n",
    "def make_model(n, x_A, x_B, x_AB):\n",
    "    theta_A = pymc.DiscreteUniform('theta_A', lower=0, upper=1)\n",
    "    theta_B = pymc.DiscreteUniform('theta_B', lower=0, upper=1)\n",
    "    \n",
    "    @pymc.stochastic(dtype=float)\n",
    "    def theta_AB(value = 0.3, theta_A=0.1, theta_B=0.2):\n",
    "        \"\"\"Return log-likelihood of it\"\"\"\n",
    "        if value < theta_A * theta_B or value > 1:\n",
    "            return -np.inf\n",
    "        else:\n",
    "            return -np.log(1 - theta_B*theta_A)\n",
    "\n",
    "    \n",
    "    @pymc.deterministic\n",
    "    def likelihood_A(theta=theta_A, n=n, x=x_A):\n",
    "        return special.binom(n, x) * np.power(theta,x) * np.power((1- theta), n-x)\n",
    "\n",
    "    @pymc.deterministic\n",
    "    def likelihood_B(theta=theta_B, n=n, x=x_B):\n",
    "        return special.binom(n, x) * np.power(theta,x) * np.power((1- theta), n-x)\n",
    "\n",
    "    @pymc.deterministic\n",
    "    def likelihood_AB(theta=theta_AB, n=n, x=x_AB):\n",
    "        return special.binom(n, x) * np.power(theta,x) * np.power((1- theta), n-x)\n",
    "\n",
    "# do you have to say the likelihood is a different variable?\n",
    "# It's conditional. Has it been observed?\n",
    "# partially, not completely. still needs theta_A, theta_B"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " [-----------------100%-----------------] 200 of 200 complete in 0.0 sec<pymc.MCMC.MCMC object at 0x10aabbb10>\n"
     ]
    },
    {
     "ename": "KeyError",
     "evalue": "'x_A'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyError\u001b[0m                                  Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-28-74fd808d90e5>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0;32mprint\u001b[0m \u001b[0mM\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mvalue\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      6\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mhist\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mM\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtrace\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'x'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 7\u001b[0;31m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mhist\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mM\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtrace\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'x_A'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/pymc/Model.pyc\u001b[0m in \u001b[0;36mtrace\u001b[0;34m(self, name, chain)\u001b[0m\n\u001b[1;32m    838\u001b[0m         \"\"\"\n\u001b[1;32m    839\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0misinstance\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mname\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstr\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 840\u001b[0;31m             \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdb\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtrace\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mname\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mchain\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    841\u001b[0m         \u001b[0;32melif\u001b[0m \u001b[0misinstance\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mname\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mVariable\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    842\u001b[0m             \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdb\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtrace\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mname\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m__name__\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mchain\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/pymc/database/base.pyc\u001b[0m in \u001b[0;36mtrace\u001b[0;34m(self, name, chain)\u001b[0m\n\u001b[1;32m    371\u001b[0m           \u001b[0mthe\u001b[0m \u001b[0mith\u001b[0m \u001b[0mcall\u001b[0m \u001b[0mto\u001b[0m \u001b[0;34m`\u001b[0m\u001b[0msample\u001b[0m\u001b[0;34m`\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    372\u001b[0m         \"\"\"\n\u001b[0;32m--> 373\u001b[0;31m         \u001b[0mtrace\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcopy\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcopy\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_traces\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mname\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    374\u001b[0m         \u001b[0mtrace\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_chain\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mchain\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    375\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0mtrace\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mKeyError\u001b[0m: 'x_A'"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAEACAYAAACuzv3DAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEK9JREFUeJzt3W+sZHV9x/H3B5AgKq4Es7sKBh5orAkpNNE00YRpDWZr\nBOFBrZpQYqghaURCohVt1Akm/qEFkqYpT0CztYaG1ErZNFHWLRPxQTWaBZGVUhs2kZa9VIupFEwE\nvn0wZ9e717v3zt879zf7fiUTzpw55/y+WWY+e/Y7c84vVYUkqV2nLLoASdJ0DHJJapxBLkmNM8gl\nqXEGuSQ1ziCXpMaNFORJTk1yMMm+7nk/yRPduoNJ9sy3TEnSiZw24nbXA4eAV3TPC7i1qm6dS1WS\npJFtekae5FzgncAdQI6uXrUsSVqgUVortwEfBV5cta6A65I8lOTOJDvmUp0kaVMbBnmSdwFPVdVB\njj8Dvx24ALgIeBK4ZW4VSpI2lI3utZLks8BVwPPAGcBZwFer6o9XbXM+sK+qLlxnf2/kIkkTqKqR\n29cbBvlxGyaXAB+pqsuS7K6qJ7v1NwBvrqr3r7NPeVOu2ej3+/T7/UWX0awkDDuCR/W7x9xH5mT4\nDPj+nK0kYwX5qL9agWFr5eg78uYkv909fxy4dozjSJJmaOQgr6oBMOiWr5pTPZKkMXllZyN6vd6i\nS1gyvUUXsFR8fy7WyD3yiQ5uj1zbxG/2yLds5JOiR67ZGrdH7hm5JDXOIJekxhnkktQ4g1ySGmeQ\nS1LjDHJJapxBLkmNM8glqXEGuSQ1ziCXpMYZ5JLUOINckhpnkEtS4wxySWrcSEGe5NQkB5Ps656f\nnWR/kseS3Jdkx3zLlCSdyKhn5NcDh/j1DZ1vBPZX1RuAA91zSdICbBrkSc4F3gncwXDeToDLgb3d\n8l7girlUJ0na1Chn5LcBHwVeXLVuZ1WtdMsrwM5ZFyZJGs2Gky8neRfwVFUdTNJbb5uqqiQnnMuq\n3+8fW+71es7tJ0lrDAYDBoPBxPtvOGdnks8CVwHPA2cAZwH/CLwZ6FXVkSS7gfur6o3r7O+cndoW\nnLNTLZnpnJ1V9YmqOq+qLgDeC/xLVV0F3Atc3W12NXDPpAVLkqYz7u/Ij55afB64NMljwO93zyVJ\nC7Bha2Xqg9ta0TZha0UtmWlrRZK0/RnkktQ4g1ySGmeQS1LjDHJJapxBLkmNM8glqXEGuSQ1ziCX\npMYZ5JLUOINckhpnkEtS4wxySWqcQS5JjTPIJalxG87ZKc3a8L7gkmZp0zPyJGck+U6SB5McSvK5\nbn0/yRNJDnaPPfMvV8uhFvCQltdIMwQlObOqnk1yGvBt4CPA24FfVNWtG+znDEE6ziJn6nGGILVi\nLjMEVdWz3eLpwKnA00fHG688SdKsjRTkSU5J8iCwAtxfVY90L12X5KEkdybZMbcqJUknNNKXnVX1\nInBRklcC30jSA24Hbuo2+QxwC3DN2n37/f6x5V6vR6/Xm6pgSVo2g8GAwWAw8f4j9ciP2yH5JPBc\nVf3lqnXnA/uq6sI129oj13HskUubm3mPPMk5R9smSV4KXAocTLJr1WZXAg+PW6wkaXqjtFZ2A3uT\nnMIw+L9cVQeS/G2Sixie5jwOXDvHOiVJJzB2a2Wsg9ta0Rq2VqTNzeXnh5Kk7csgl6TGGeSS1DiD\nXJIaZ5BLUuMMcklqnEEuSY0zyCWpcQa5JDXOIJekxhnkktQ4g1ySGmeQS1LjDHJJapxBLkmNM8gl\nqXEbBnmSM5J8J8mDSQ4l+Vy3/uwk+5M8luS+o1PBSZK23qYzBCU5s6qeTXIa8G3gI8DlwE+r6uYk\nHwNeVVU3rrOvMwTpOM4QJG1u5jMEVdWz3eLpwKnA0wyDfG+3fi9wxZh1SpJmZNMgT3JKkgeBFeD+\nqnoE2FlVK90mK8DOOdYoSdrAaZttUFUvAhcleSXwjSS/t+b1SnLCfzv2+/1jy71ej16vN3Gxmp1h\ni0PSdjAYDBgMBhPvv2mP/LiNk08CzwF/AvSq6kiS3QzP1N+4zvb2yLepk7FXbY9crZhpjzzJOUd/\nkZLkpcClwEHgXuDqbrOrgXsmK1eSNK3NWiu7gb1JTmEY+l+uqgNJDgJ3J7kGOAy8Z75lSpJOZKzW\nytgHt7Wybdla2bpx/QxoXDP/+aEkaXszyCWpcQa5JDXOIJekxhnkktQ4g1ySGmeQS1LjDHJJapxB\nLkmNM8glqXGb3sZWUnsWeZtib0mw9QxyaWkt6p422mq2ViSpcQa5JDXOIJekxo0y+fJ5Se5P8kiS\nHyb5cLe+n+SJJAe7x575lyu1J8mWP3Ry2XRiiSS7gF1V9WCSlwPfB65gOCvQL6rq1g32dWKJbcqJ\nJRx3XuP6mZ/euBNLbPqrlao6Ahzplp9J8iPgtUfHm6hKSdLMjNUjT3I+cDHwr92q65I8lOTOo5M0\nS5K21shB3rVV/gG4vqqeAW4HLgAuAp4EbplLhZKkDY10QVCSlwBfBf6uqu4BqKqnVr1+B7BvvX37\n/f6x5V6vR6/Xm7xaSVpCg8GAwWAw8f6jfNkZYC/ws6q6YdX63VX1ZLd8A/Dmqnr/mn39snOb8stO\nx53XuH7mpzful52jBPnbgG8BP+DX74xPAO9j2FYp4HHg2qpaWbOvQb5NGeSOO69x/cxPb+ZBPmUx\nBvk2ZZA77rzG9TM/vXGD3Cs7JalxBrkkNc4gl6TGGeSS1DiDXJIaZ5BLUuMMcklqnEEuSY0zyCWp\ncQa5JDXOIJekxhnkktQ4g1ySGmeQS1LjDHJJapxBLkmN2zTIk5yX5P4kjyT5YZIPd+vPTrI/yWNJ\n7kuyY/7lSpLWGmWqt13Arqp6MMnLge8DVwAfAH5aVTcn+Rjwqqq6cc2+zhC0TTlDkOPOa1w/89Ob\n+QxBVXWkqh7slp8BfgS8Fric4aTMdP+9YvxyJUnTGqtHnuR84GLgO8DOVZMtrwA7Z1qZJGkkIwd5\n11b5KnB9Vf1i9Wtd/8R/T0nSApw2ykZJXsIwxL9cVfd0q1eS7KqqI0l2A0+tt2+/3z+23Ov16PV6\nUxUsSctmMBgwGAwm3n+ULzvDsAf+s6q6YdX6m7t1X0hyI7DDLzvb4Zedjjuvcf3MT2/cLztHCfK3\nAd8CfsCv3xkfB74L3A28DjgMvKeqfr5mX4N8mzLIHXde4/qZn97Mg3zKYgzybcogd9x5jetnfnoz\n//mhJGl7M8glqXEGuSQ1ziCXpMYZ5JLUOINckho30pWdmo/hTwAlaToG+cIt6je3/iUiLQtbK5LU\nOINckhpnkEtS4wxySWqcQS5JjTPIJalxBrkkNc4gl6TGbRrkSb6YZCXJw6vW9ZM8keRg99gz3zIl\nSScyyhn5l4C1QV3ArVV1cff4+uxLkySNYtMgr6oHgKfXeclrvCVpG5imR35dkoeS3Jlkx8wqkiSN\nZdKbZt0O3NQtfwa4BbhmvQ37/f6x5V6vR6/Xm3BISVpOg8GAwWAw8f4ZZcbrJOcD+6rqwjFfK2fU\nPrHFzWQPJ+Ps7o67NeP6mZ9eEqpq5Pb1RK2VJLtXPb0SePhE20qS5mvT1kqSu4BLgHOS/AT4NNBL\nchHDv/IfB66da5WSpBMaqbUy8cFtrWzI1orjLuO4fuantyWtFUnS9mGQS1LjDHJJapxBLkmNM8gl\nqXEGuSQ1ziCXpMYZ5JLUOINckhpnkEtS4ya9ja0krWt464mtdzLfGsAglzRji7q3zMnL1ookNc4g\nl6TGGeSS1DiDXJIat2mQJ/likpUkD69ad3aS/UkeS3Jfkh3zLVOSdCKjnJF/CdizZt2NwP6qegNw\noHsuSVqATYO8qh4Anl6z+nJgb7e8F7hixnVJkkY0aY98Z1WtdMsrwM4Z1SNJGtPUFwRVVSU54RUA\n/X7/2HKv16PX6007pCQtlcFgwGAwmHj/jHJZa5LzgX1VdWH3/FGgV1VHkuwG7q+qN66zX53Ml81u\nZngp86L+fE6+2d0dd7nHXaasSUJVjXy56qStlXuBq7vlq4F7JjyOJGlKm56RJ7kLuAQ4h2E//FPA\nPwF3A68DDgPvqaqfr7OvZ+Qb8IzccR13duMuU9aMe0Y+UmtlimIM8g0Y5I7ruLMbd5myZqtaK5Kk\nbcIgl6TGGeSS1DiDXJIaZ5BLUuMMcklqnEEuSY0zyCWpcQa5JDXOIJekxhnkktQ4g1ySGmeQS1Lj\nDHJJapxBLkmNM8glqXFTTb6c5DDwv8ALwK+q6i2zKEqSNLqpgpzhVCC9qvqfWRQjSRrfLForI09H\nJEmavWmDvIBvJvlekg/OoiBJ0nimba28taqeTPJqYH+SR6vqgdUb9Pv9Y8u9Xo9erzflkJK0XAaD\nAYPBYOL9M6uZp5N8Gnimqm5Zta6WaWbrWUsWNeM4nIyzrDvuco+7TFmThKoauW09cWslyZlJXtEt\nvwx4B/DwpMeTJE1mmtbKTuBrw7NKTgO+UlX3zaQqSdLIZtZaWffgtlY2ZGvFcR13duMuU9ZsWWtF\nkrQ9GOSS1DiDXJIaZ5BLUuMMcklqnEEuSY0zyCWpcQa5JDXOIJekxk1798OZOXToEDfd9Be88MLW\nj33uua/mtttu3vqBJWkGts0l+gcOHOCyyz7Ec8/92dzqWd9TwI1bPOZqXqLvuI47m3EXYx4ZOu4l\n+tvmjBzg9NNfw3PPfWCLR/0PhkF+cr35pOVz8n6G7ZFLUuMMcklqnEEuSY2bKsiT7EnyaJJ/T/Kx\nWRUlSRrdNFO9nQr8NbAHeBPwviS/NavCtNZg0QUsmcGiC1gyg0UXcFKb5oz8LcCPq+pwVf0K+Hvg\n3bMpS79psOgClsxg0QUsmcGiCzipTRPkrwV+sur5E906SdIWmuZ35DP/0eYvf/kQZ5112awPu6EX\nX/w/nnlmS4eUpJma+MrOJL8L9KtqT/f848CLVfWFVdssz2yokrSFxrmyc5ogPw34N+DtwH8B3wXe\nV1U/muiAkqSJTNxaqarnk3wI+AZwKnCnIS5JW2+uN82SJM3fzK/sTPKHSR5J8kKS31nz2se7i4ce\nTfKOWY+97JL0kzyR5GD32LPomlrkhWyzleRwkh9078nvLrqeliT5YpKVJA+vWnd2kv1JHktyX5Id\nmx1nHpfoPwxcCXxr9cokbwL+iOHFQ3uAv0niLQLGU8CtVXVx9/j6ogtqjReyzUUBve49+ZZFF9OY\nLzF8L652I7C/qt4AHGCE+2zPPEir6tGqemydl94N3FVVv6qqw8CPGV5UpPFsj/tmtssL2ebD9+UE\nquoB4Ok1qy8H9nbLe4ErNjvOVp4Rv4bhRUNHeQHRZK5L8lCSO0f5J5d+gxeyzV4B30zyvSQfXHQx\nS2BnVa10yyvAzs12mOhXK0n2A7vWeekTVbVvjEP5TesaG/zZ/jlwO3BT9/wzwC3ANVtU2rLwPTd7\nb62qJ5O8Gtif5NHuTFNTqqoa5XqciYK8qi6dYLf/BM5b9fzcbp1WGfXPNskdwDh/aWpo7fvwPI7/\nl6LGVFVPdv/97yRfY9i+Msgnt5JkV1UdSbKb4XyUG5p3a2V13+xe4L1JTk9yAfB6hhcRaUTd/9Sj\nrmT4xbLG8z3g9UnOT3I6wy/g711wTc1KcmaSV3TLLwPege/Lad0LXN0tXw3cs9kOM5+zM8mVwF8B\n5wD/nORgVf1BVR1KcjdwCHge+NORZ2bWUV9IchHD9sDjwLULrqc5Xsg2czuBryWBYZ58paruW2xJ\n7UhyF3AJcE6SnwCfAj4P3J3kGuAw8J5Nj2OWSlLb/B23JDXOIJekxhnkktQ4g1ySGmeQS1LjDHJJ\napxBLkmNM8glqXH/D0dhcBSbNTQQAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10a74f110>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "M = pymc.MCMC(make_model(n, x_A, x_B, x_AB))\n",
    "\n",
    "\n",
    "M.sample(200)\n",
    "print M.value\n",
    "plt.hist(M.trace('x')[:])\n",
    "plt.hist(M.trace('x_A')[:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " [-----------------100%-----------------] 50000 of 50000 complete in 2.8 sec"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pymc\n",
    "import scipy.optimize as opt\n",
    "import scipy.stats as stats\n",
    "\n",
    "sigma = 1.0\n",
    "tau = 1/ sigma**2\n",
    "\n",
    "sigma0 = 3.0\n",
    "tau0 = 1/sigma0**2\n",
    "\n",
    "mu = pymc.Normal(\"mu\", mu=0, tau=tau0)\n",
    "x = pymc.Normal(\"x\", mu=mu, tau=tau) # a bunch of samples for x\n",
    "mcmc = pymc.MCMC([mu, x])\n",
    "mcmc.sample(50000, 10000, 1) #1 is thinning. Results in autocorrelation , serial dependent\n",
    "# First set burn=0, thin = 1, examine traces"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "x_samples = mcmc.trace(\"x\")[:]\n",
    "mu_samples = mcmc.trace(\"mu\")[:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/matplotlib/figure.py:387: UserWarning: matplotlib is currently using a non-GUI backend, so cannot show the figure\n",
      "  \"matplotlib is currently using a non-GUI backend, \"\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUMAAAE4CAYAAADB1pEFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xt01PWd//HnOwkhEO53hXARUECLhQJV0BK7aFl+Vmpd\nUap1reuWU8X668WyulrC6TnbrW23W49d63ZZ125VWrV4QamX/oxi5S53gYIQScIl3CTcScjn98dc\nHMYkM5PM5DOX1+OcOcx85/udeWeYvPL5fL+f7+drzjlERHJdnu8CRETSgcJQRASFoYgIoDAUEQEU\nhiIigMJQRASIIwzNbKqZbTGzbWY2p5HnR5jZUjM7ZWbfa+T5fDNbY2YvJ6toEZFkazYMzSwfeBSY\nCowCZprZyKjVDgL3AD9r4mXuBT4ANKBRRNJWrJbhBGC7c67COVcHLACmR67gnNvvnFsF1EVvbGYD\ngGnAfwGWnJJFRJIvVhj2ByojHlcFl8XrF8B9QEOCdYmItKlYYdjirq2ZXQvUOOfWoFahiKS5ghjP\nVwMlEY9LCLQO4zERuM7MpgFFQBcz+61z7rbIlcxM+xJFJCWcc3E3xGK1DFcBw81ssJkVAjcBLzWx\n7jlv6px7wDlX4pwbAtwM/L/oIIxYNyNvc+fO9V5DLtWdybVnat2ZXHuimm0ZOufqzWw28BqQD8x3\nzm02s1nB5x83s37ASqAL0GBm9wKjnHPHol8u4epERNpIrG4yzrnFwOKoZY9H3N/LuV3pxl7jbeDt\nFtYoIpJyOgOlFUpLS32X0CKZWjdkbu2ZWjdkdu2JsJb0rZNagJnzXYOIZB8zwyXxAIqISE5QGIqI\noDAUEQEUhiIigMJQRARQGIqIAApDERFAYSgiAigMRUQAhaGICKAwFBEBFIYiIoDCUEQEUBiKiAAK\nQxERQGEoIgIoDEVEAIWhiAigMBQRARSGIiJAHJcKFfHt6NGj/OEPf2Dt2rUUFxczZcoUvvjFL5KX\nF/hbbvbpa/7oImOSKF0dT9LawoUL+eY3v8mBAweaXa+srOyc+/pOia6OJ1nj17/+NV/96lc5cOAA\nl19+OT/72c8AKC4uBqBz587Mnj3bZ4mSRRSGkpYWLVrEt771rfDjpUuX8v3vfx+Au+++m0GDBnH0\n6FGeeuopXyVKltE+Q0k71dXVfP3rXwfgqquuYvLkyeHnysrK6NixI7fccgtPPPEEe/bsAQL7CBvb\ndygSL7UMJe3Mnj2bjz/+GIAvfOELja5TWFjITTfdRGFhIQCbNm1qs/okOykMJa0sWbKEF154Ibxf\nsLnWXrdu3bjmmmsAeP311zlz5kyb1CjZSWEoacM5xw9+8AOA8P7BWMaOHQtAbW0tK1euTFltkv0U\nhpI2Fi5cyLJly+jTpw/f+9734tomNNYQAgdZ6urqUlWeZLm4wtDMpprZFjPbZmZzGnl+hJktNbNT\nZva9iOUlZvaWmW0ys41m9u1kFi/Zw8y44YYbAKipqaFLly4Jbd+vXz+OHTvG2rVrU1Ge5ICYYWhm\n+cCjwFRgFDDTzEZGrXYQuAf4WdTyOuA7zrmLgcuAuxvZViSsQ4cO/PM///M5g6jjceWVVwLw7rvv\npqAqyQXxtAwnANudcxXOuTpgATA9cgXn3H7n3CoC4Re5fK9zbm3w/jFgM3B+UiqXrDR27FjatWuX\n8HYjR46kZ8+eHDlyJAVVSS6IJwz7A5URj6uCyxJiZoOBMcDyRLeV7BYaK2hmjB8/vkWvkZeXx7hx\n45JZluSYeMKw1Sd5mlkn4Dng3mALUSRs/vz5AIwYMYJu3bq1+HVGjx4dPqBSVVWVlNokd8RzBko1\nUBLxuIRA6zAuZtYOeB74nXPuhcbWidw/VFpaSmlpabwvLxnOOcfvfvc74JNhMi1VXFzMyJEj2bRp\nE//zP//Dgw8+mIwSJUOUl5dTXl7e4u3jCcNVwPBgN3c3cBMws4l1zxkha4ERs/OBD5xz/97UGyS6\ns1yyx+rVq9m6dSsAF1xwQatfb8yYMWzatIn58+fzwAMPnDP0RrJbdENq3rx5CW0f85vinKsHZgOv\nAR8Av3fObTazWWY2C8DM+plZJfAd4EEz2xXsGk8CbgWuMrM1wdvUhCqUrBY50UJ+fn6rXy8UqBUV\nFSxbtqzVrye5I66JGpxzi4HFUcsej7i/l3O70iHvooHd0oT6+nqeeeaZpL5mZEvw2WefZeLEiUl9\nfcleCirx5q233mLfvn0MHz48Ja//3HPP0dDQkJLXluyjMBRvFi5cCMBNN92U9NceMGAAVVVVrFix\nIumvLdlJYSheNDQ08OKLLwLwla98Jemv/3d/93dAoKssEg+FobQ5MyM/P5/du3cDpGSw9I033ggE\nusq6HorEQ2EoXoTOJR4/fnxKhlZNmjQJgF27dpGXl6dZsCUmhaF4ERpbeNFFF6Xk9cvKyvjsZz8L\nwJQpU1LyHpJdFIbiRU1NDe3bt2fw4MEpe4/QUept27al7D0keygMxZuhQ4dSUJC6a5INHTqUvLw8\ndu3albL3kOyhMBRvhg0bltLXLyoqYuDAgTqAInFRGEqbipyWf+jQoSl/v1QN6JbsozCUNrV06VIA\nevXqRdeuXVP+fpFhqBaiNEdhKG3qtddeA1LfRQ7p3bs3nTp1AmDz5s1t8p6SmRSG0qZef/11oG26\nyBAY4D1kyBAA3nzzzTZ5T8lMCkNpMwcPHmT16tUADBo0qM3eNzSt15///Oc2e0/JPApDaTPvvPNO\neL9dYWFhm71vqGVYXl5OfX19m72vZBaFobSZt956y8v7hq6rUltby6pVq7zUIOlPYShtpjXXp0gW\ndZWlKQpDaRMHDhxgw4YNFBUVea0jHQJZ0pPCUNrE22+/DcDll1/utY6lS5dqv6E0SmEobSLUIrvq\nqqu81TBs2DCOHz/O2rVrvdUg6UthKG0iFIY+r4kdmkNxyZIl3mqQ9KUwlJQ7dOgQGzdupH379kyY\nMMFbHQpDaY7CUFKuZ8+eAJw+fdrrAZQ77rgDCFyIysw0+7WcQ2EobWbSpEkpmeI/XnPnzg2fp3z3\n3Xd7q0PSk8JQ2kxJSYnX9zczBg4cCKAJX+VTFIaSUpHzF/oOQ0BhKE1SGEpKrVmzBgjsNywuLvZc\nzScTRHz00UeeK5F0ozCUlPrLX/4CpEerEKBv3760b9+ejz/+2HcpkmYUhpJSoTAMdU99y8vLS5tg\nlvSiMJSUcc6lXcsQ0qsWSR8KQ0mZiooK9u7dC3wy1jAdDBgwwHcJkoZihqGZTTWzLWa2zczmNPL8\nCDNbamanzOx7iWwr2S3UKoRA9zRd9O/fP3z/9OnTHiuRdNLsN9TM8oFHganAKGCmmY2MWu0gcA/w\nsxZsK1ksMgzTSVFREb169QLQpA0SFuvP9QRgu3OuwjlXBywApkeu4Jzb75xbBdQluq1kt2XLlvku\noUmhrnI61yhtK1YY9gcqIx5XBZfFozXbSoY7efIkGzZsSKvucaRQGC5fvtxzJZIuYn1TW3PVbV2x\nO4etWbOGs2fPcvHFF/supVFqGUq0ghjPVwOR4xBKCLTw4hH3tpEn75eWlnqd806SY8WKFQBMmDCB\nDRs2eK7m03r37g3Azp07qampoU+fPp4rktYqLy9v1WUdYoXhKmC4mQ0GdgM3ATObWDd6PqS4t/U5\nk4mkRigMx48fz/z58z1X82n5+fnh+8uXL+fLX/6yx2okGaIbUvPmzUto+2a7yc65emA28BrwAfB7\n59xmM5tlZrMAzKyfmVUC3wEeNLNdZtapqW0Tqk4y1sqVKwG8TuYaL3WVBWK3DHHOLQYWRy17POL+\nXs7tDje7rWS/Q4cOsX37doqKirjkkkt8lxOTDqII6AwUSYFQq3Ds2LG0a9fOczWxrVixgrNnz/ou\nQzxTGErShcJw/PjxniuJraSkhKNHj7J5s/bg5DqFoSRd5JHkdBeqcdWqVZ4rEd8UhpJUZsbLL78M\nwC233JL2F10aN24cAKtXr/ZcifimMJSUKCoqYu7cuWk/bOr+++8H4NFHH9UV83KcwlBSon///hkR\nLHPmBCZTKigo4KGHHvJcjfikMJSUiJwmK5116NCB7t27U19fz/79+32XIx4pDCUlMiUMAc4//3wA\ndu/e7bkS8UlhKEkTOVYvFDCZQGEooDCUJNq6dSsAXbp0oXPnzp6riZ/CUEBhKEkUGl+YSV1kgPPO\nOw+Affv2ea5EfFIYStK8//77QGZ1kSEwDKhnz546JS/HKQwlaUJhGGppZZJMC3BJPoWhJMXZs2fD\nF1dSGEomUhhKUmzbto3jx48DUFxc7LmaxCkMRWEoSRHqImeqfv36he+fOnXKYyXii8JQkiLTw7B9\n+/bhaymn4zVbJPUUhpIUmR6G8ElXWdN55SaFobSacy6rwlDTeeUmhaG02s6dOzly5Ah9+/b1XUqr\nhMIwNFO35BaFobRaqFU4duxYz5W0Tuggyvr168NzG2bCNGSSHApDabVsCcPCwsLw/TvvvDPtJ6aV\n5FIYSqtlSxhG2rNnj+8SpI0pDKVVIg+eZFMYagab3KMwlFaprq5m//79dO/enUGDBvkuJ2nUMsw9\nCkNplchWYTYdbKipqaG+vt53GdKGFIbSKtnYRe7VqxcNDQ3U1NT4LkXakMJQWiUbwzA06472G+YW\nhaG0SigMx4wZ47mS5AmFofYb5pYC3wVI5orcRzhixAiPlSRX6EwUhWFuURhKqw0cOJA77rgDICsG\nKofORNE1UXJLzG6ymU01sy1mts3M5jSxziPB59eZ2ZiI5feb2SYz22BmT5tZ+2QWL+khE2e2bk5R\nURE9evTQNVFyTLNhaGb5wKPAVGAUMNPMRkatMw0Y5pwbDnwTeCy4fDDwj8BY59xngHzg5iTXL2kg\n28IQNPN1LorVMpwAbHfOVTjn6oAFwPSoda4DngRwzi0HuplZX6AWqAM6mlkB0BGoTmbxkh6yMQyz\n8WeS5sUKw/5AZcTjquCymOs45w4BPwd2AbuBj51zb7auXEkXhw8fBqCgoCA8Q3Q2URjmnlgHUFyc\nr/OpUw/MbCjwf4HBwBHgWTO7xTn3VPS6kTvdS0tLKS0tjfNtxZfQlfD69u1Lfn6+52qSLzIM6+rq\naNeuncdqJB7l5eWUl5e3ePtYYVgNlEQ8LiHQ8mtunQHBZaXAe865gwBm9kdgItBsGEpmyORrJMej\nQ4cOdO/encOHD7N582ZGjx7tuySJIbohNW/evIS2j9VNXgUMN7PBZlYI3AS8FLXOS8BtAGZ2GYHu\n8D5gK3CZmXWwwIC0KcAHCVUnaSvbwxA++dmy4ZIGEluzYeicqwdmA68RCLLfO+c2m9ksM5sVXOdV\nYIeZbQceB+4KLl8L/JZAoK4PvuR/puSnkDaXC2Goa6LklpiDrp1zi4HFUcsej3o8u4ltHwYebk2B\nkn6OHTvG1q1bAejTp4/nalInFPQKw9ygc5MlYevWrcO5wLG1goLsPYkpFIZr167VAOwcoDCUhOXK\nPrSOHTsCcPLkSbZs2eK5Gkk1haEkLFfCMJK6ytlPYSgJy8UwzMWfOdcoDCUhp06dYtOmTeTl5dZX\nRy3D7Jdb32hptQ0bNnD27Nmsmr8wHmvWrKGhocF3GZJCCkNJSDZO8x9LSUkJx48f569//avvUiSF\nFIaSkFwMw9DPqq5ydlMYSlzMDDPjP/8zcBLRd7/7Xc8VtZ3Pfe5zgA6iZDuFocTtoYceCs9QM2dO\no5OeZ6VQGKplmN0UhhK3/fv3c/bsWbp3706HDh18l9NmQt3k999/XwdRspjCUOIWulpcNk/O0Jh+\n/fpx/vnnc/ToUT788EPf5UiKKAwlbrkahqCDKLlAYShxy+Uw1EGU7Je9U45I0u3duxfIvTAMzE0c\n8NOf/pSf/vSn4Vl7JHuoZShxq6uro2vXrhQXF/supU2VlZWFhxIVFRV5rkZSRWEoCcm1VmFI586d\nKS4u5tSpU75LkRRRGEpCcjUMzUwXls9yCkNJSC4HQq7+IcgVCkOJKXKgcS4HQi7/7LlAYSgxbdu2\nDQjsN+vUqZPnavyJDEMdTc4+CkOJKTTQOJe7yABdu3YNXxdl165dnquRZFMYSky5cI3keJiZLh+a\nxRSGElPoFz/XwxB0LeVspjCUZjU0NKhlGCG0q0Cn5WUfhaE0a8eOHdTW1gLQpUsXz9X4F9ky1EGU\n7KIwlGapBXSubt26AYG5HauqqjxXI8mkMJRmad/YuSInbdAfiuyiMJRmKQybps8muygMpUnOObV+\nmqHPJrsoDKVJFRUVHD58mN69e/suJS2pZZhdYoahmU01sy1mts3MGr0kmpk9Enx+nZmNiVjezcye\nM7PNZvaBmV2WzOIltUItn9Asz/KJrl27snfvXnbv3u27FEmSZsPQzPKBR4GpwChgppmNjFpnGjDM\nOTcc+CbwWMTTvwRedc6NBEYDm5NYu6RYqOWTSxeMj1fkFfMkO8RqGU4AtjvnKpxzdcACYHrUOtcB\nTwI455YD3cysr5l1Ba50zv138Ll659yR5JYvqRQKQ7UMP03XUs4+scKwP1AZ8bgquCzWOgOAIcB+\nM3vCzN43s9+YWcfWFixtI/LgiVqGn6ar5WWfWBeEineIvUU9dsHXHgvMds6tNLN/B/4J+GH0xmVl\nZeH7paWllJaWxvm2kiqVlZUcOHCAHj16MGjQIN/lpB1dLS/9lJeXU15e3uLtY4VhNVAS8biEQMuv\nuXUGBJcZUOWcWxlc/hyBMPyUyDCU9BB58CRyoLEEDBs2jM6dO1NdXc2+ffvo27ev75JyXnRDat68\neQltH6ubvAoYbmaDzawQuAl4KWqdl4DbAIJHiz92zu1zzu0FKs3swuB6U4BNCVUn3qxatQpQF7kp\neXl5jBkTGDih1mF2aLZl6JyrN7PZwGtAPjDfObfZzGYFn3/cOfeqmU0zs+3AceAbES9xD/BUMEg/\njHpO0lRkS/AnP/kJP/nJTzxWk54iP6Np06YBmv0608W8iLxzbjGwOGrZ41GPZzex7TpgfGsKFD+K\nioo4deoU3/nOd+jatat2ZUQpKytj3bp1LFy4kBEjRrBlyxbfJUkr6QwUadSpU6fo1KmTpu1qRmhu\nwz179niuRJJBYShN6t+/vw6eNKNnz560a9eOI0c0fDYbKAylSbl+AahY8vLy6Nevn+8yJEkUhtKk\n/v2jx9dLNP3ByB4KQzlHXV1d+L5+0WPTdWGyh8JQzrFpU2AoaI8ePcLXCJamKQyzh8JQzrFixQpA\nrcJ49erVi4KCwAi1Q4cOea5GWkNhKOcIhaH2F8YnPz8/fBBFZ6JkNoWhnGPlysCp5ArD+IW6ygrD\nzKYwlLDjx4+zceNGAA0ZSUBol4Km88psCkMJW7NmDQ0NDQAUFhZ6riZzhFqGocktJDMpDCUstL9Q\nEhO6YNaOHTvYv3+/52qkpRSGEhbaXyiJyc/PD99ftmyZx0qkNRSGEqaWYestXbrUdwnSQgpDAeDg\nwYPs2LGDDh06+C4lo/34xz/GzMI3yRwKQwFg+fLlgK6E11rt2rXjoYce0vyPGUhhKAC89957AEyc\nONFzJZmre/fu1NXVUVNT47sUaQGFoQCf7Ou6/PLLPVeSuQYMGABAVVX0NdMkEygMhfr6+vDBE4Vh\nyykMM5vCUNi4cSPHjh3jggsu0CUvW6GkJHDF3MrKSs+VSEsoDEVd5CTp27cvBQUFHDp0iBMnTvgu\nRxKkMBSFYZLk5+eHz1NWVznzKAwlfCRZYdh66ipnrpjXTZbsVlNTw4cffghojGEy6CBK5lIY5rjQ\nubSDBw/m9ttvB9CA4VYIhWF1dbXnSiRR6ibnuFAXOfRLLK3TuXNnunXrxpkzZ3yXIglSGOa40MGT\n0L4uaT39YclMCsMcdubMmfC0XfoFTh59lplJYZjDVq9ezcmTJwEoLi72XE32GDRokO8SpAUUhjls\nyZIlvkvISn379qV9+/aAjipnkphhaGZTzWyLmW0zszlNrPNI8Pl1ZjYm6rl8M1tjZi8nq2hJjnfe\necd3CVkpLy8vvA9Wf3AyR7NhaGb5wKPAVGAUMNPMRkatMw0Y5pwbDnwTeCzqZe4FPgBcsoqW1jt7\n9izvvvuu7zKyVqirrD84mSNWy3ACsN05V+GcqwMWANOj1rkOeBLAObcc6GZmfQHMbAAwDfgvQNP+\nppGNGzdy5MgR7d9KkdDnqpZh5ogVhv2ByPOKqoLL4l3nF8B9QEMrapQUCLVYrrzySs+VZKfQOcqb\nNm3i4MGDnquReMQKw3i7ttGtPjOza4Ea59yaRp4Xz0Itli984QueK8lOBQWfnNzVq1cvXRclA8Q6\nHa8aiByNW0Kg5dfcOgOCy24ArgvuUywCupjZb51zt0W/SeTpX6WlpZSWlsZZvrSEcy7cMlQYpt7l\nl1/Ol770JUCnOqZSeXk55eXlLd4+VhiuAoab2WBgN3ATMDNqnZeA2cACM7sM+Ng5txd4IHjDzCYD\n328sCEFfkLa2bds29u3bR58+fbjwwgt9l5P1PvroI98l5ITohtS8efMS2r7ZMHTO1ZvZbOA1IB+Y\n75zbbGazgs8/7px71cymmdl24DjwjaZeLqHKJGVCXeQrr7xS3bYUy8vLY8+ePZw+fTo89lDSU8xZ\na5xzi4HFUcsej3o8O8ZrvA283ZICJfl08KTtnHfeeVRXV1NZWcmwYcN8lyPN0BkoOcY5x1tvvQXA\n5MmTPVeT/UJDbCoqKvwWIjEpDHPM9u3bqayspFevXowePdp3OVnvggsuAGDHjh2eK5FYNLlrjgkd\nMDlw4AD5+fmeq8l+AwcODO83DE2KIelJYZijrr32WsaNGwfoaH4qFRYWUlJSwkcffaSucppTNzmH\nNDR8ciLQkCFDPFaSW0KftbrK6U1hmEPWrl0LQNeuXenRo4fnanJHaL/hzp07PVcizVEY5pA///nP\nQOCXU+ML207//v0pLCzkwIEDvkuRZigMc0goDNVFblv5+fmaHSgDKAxzxOnTp8NnnigM216oqyzp\nS2GYI5YtW8aJEyeAwOUspW1F/gFyTmempiOFYY4IdZHFjz59+tCxY0cAtm7d6rkaaYzCMEcsXrw4\n9kqSMnl5eQwdOhSAP/3pT56rkcYoDHPAvn37WLVqFUVFRb5LyWnDhw8H4NVXX/VciTRGYZgDQi2R\nq666ynMluS00a83bb7/NsWPHPFcj0RSGOeCVV14BYNq0aZ4ryW2hfYZnzpzRPtw0pDDMcnV1dbz+\n+uuAwjCdqKucfhSGWW7p0qUcOXKEiy66SGPd0sirr76qITZpRmGY5UItELUK00e/fv2oqqpi48aN\nvkuRCArDLKcwTD979+4FYPTo0TpHPI0oDLNYZWUlGzZsoLi4WNc7SSM33ngjgM5XTjMKwyy2aNEi\nAKZMmaIrs6WRoUOHkpeXx65du3yXIhEUhlns+eefB+D666/3XIlEKioqYsiQITqAkmYUhlnq4MGD\n4bFst99+O2am/VNpZOTIkb5LkCi6BkqWevHFF4HA1FG33XZbeLmud5IeRowYEd6NUVtbS5cuXTxX\nJGoZZqk//vGPAIwaNcpzJdKYTp06hQ+ghM4QEr8UhlmotraWN954Awi0QCQ9hbrKoX274pfCMAst\nWrSIM2fOAIEWiKSnUBguXrw4PPGu+KMwzEKhLrKkt65duwJw4sQJzXGYBhSGWebEiROayDUDqavs\nn8Iwy7z44oucOHGCz3/+875LkQS88MILmuPQM4VhlgiNI/za174GwPLlyz1XJPGaOHEiJ06cYOHC\nhb5LyWlxhaGZTTWzLWa2zczmNLHOI8Hn15nZmOCyEjN7y8w2mdlGM/t2MouXc913333k5eVhZnz/\n+9/3XY7EKTQO9H//9389V5LbYoahmeUDjwJTgVHATDMbGbXONGCYc2448E3gseBTdcB3nHMXA5cB\nd0dvK8mzadMmGhoaGDZsmI4iZ5AZM2ZQWFjIm2++SXV1te9yclY8LcMJwHbnXIVzrg5YAEyPWuc6\n4EkA59xyoJuZ9XXO7XXOrQ0uPwZsBs5PWvVyjvXr1wPwmc98xnMlkoju3btz7bXX4pzj6aef9l1O\nzoonDPsDlRGPq4LLYq0zIHIFMxsMjAG0MytFqqqqKCws1EDrDGNm4eFQP/jBD3QOuSfxnJsc79Qa\n0f+D4e3MrBPwHHBvsIUoKXLJJZdQWFjouwxJQFlZGfX19fz85z/n5MmTvsvJWfGEYTVQEvG4hEDL\nr7l1BgSXYWbtgOeB3znnXmjsDSInDygtLaW0tDSOsiQkdLYJwNixYz1WIi1VUFDAJZdcwsqVK32X\nkrHKy8spLy9v8fbxhOEqYHiwm7sbuAmYGbXOS8BsYIGZXQZ87JzbZ4H2/nzgA+fcvzf1BppJpXVC\ns5/07t2b/v2j92BIphgzZkw4DE+cOBG+tKjEJ7ohNW/evIS2j7nP0DlXTyDoXgM+AH7vnNtsZrPM\nbFZwnVeBHWa2HXgcuCu4+STgVuAqM1sTvE1NqEKJ6Te/+Q0QaBVqf1PmOv/888N/zP7whz94rib3\nxDWfoXNuMbA4atnjUY9nN7Ldu2hgd0pt27YtfF7r6NGjPVcjrTVu3Diqq6t57LHHuP32232Xk1MU\nVBnuV7/6Vfh+cXGxx0okGS6++GIAVqxYwfvvv++5mtyiMMxgx44d44knnvBdhiRR5EiARx55xGMl\nuUdhmMF++9vfUltbyxVXXOG7FEmyvLw8nn76aXbv3u27lJyhMMxQZ8+e5Re/+AUA99xzj+dqJNmu\nv/566urqztkNIqmlMMxQzz//PNu3b+eCCy7gq1/9qu9yJMlC8xv+y7/8i65s2EYUhhnIOcePf/xj\nIDBTTUGBLnKYbcrKyhgwIHBG65e+9CXP1eQGhWEGev3111m7di0A3/rWt9RqyFJXXnklAH/5y188\nV5Ib1KTIMM455s6dC8CUKVPCB090Fk/2ufDCCznvvPPYs2eP71JyglqGGWbRokXhWazHjx/vuRpJ\nJTNj8uTJ4ceaxCG1FIYZpKGhgQcffDD8uH379h6rkbZw0UUX0a9fP0DjDlNNYZhBFixYwPr168M7\n1iX7mRmPb1D6AAAJgklEQVRXX301EDiyfODAAc8VZS+FYYY4fvw4c+YELj+j/YO5ZejQoQDU1tbS\nu3fv8FAbHThLLoVhhnj44Yepqqpi7NixOoE/h+Xl5XHXXXfpD2IKKAwzwM6dO3n44YcB+OUvf0l+\nfr7nisSHz33uczQ0NPDKK6/gXLwT0Eu8FIZpzjnHrFmzOHXqFDNnztR5yDnsb/7mb+jYsSMfffRR\n+OJfkjwKwzT35JNP8sYbbwDwzDPPaF9RDuvYsSPXXHMNQHgOS0keDbpOY9XV1Xz3u98FAifuX3rp\npeHntM8oN1166aWsX7+eHTt2AIGeg/44Jodahmnq7Nmz3HrrrRw+fBjQLNYSYGZMnz49PMZ0/vz5\nnivKHgrDNPWv//qvlJeX07dvXwD99Zewrl27Mm3aNAC+/e1va/9hkigM09Brr73GD3/4QyAwgatI\ntFBP4eTJk9xwww0cOXLEc0WZT2GYZsyMqVOn0tDQAGj6JmlcqKdw6aWXsn37dm6++Wbq6+s9V5XZ\ndAAljezduzd8f8SIEcyYMYO8vDwdLJEmrVu3DggcXW7Xrh0NDQ3apdJCahmmiY8//jjcCjzvvPO4\n/vrrycvTf480r6ysjH/4h38IT/D74IMPakB2C+m3LQ0cPHiQa665Jrwj/NZbb9WMNBK3kpISbrjh\nBiAwmcPcuXMViC2gMPRs9+7dTJ48mZUrVzJkyBBA1z+WxI0cOTJ8/0c/+hF5eXkaoJ8ghaFHZkb/\n/v3ZtGkTEDgHWaQ1brzxxvC56yNGjPBcTWZRGHrgnDtnsOzAgQO57777dKBEWu3iiy8O72bZsmUL\nAB988IHnqjKDwrCN7du3jxkzZnDnnXcCMGbMGG677TZ1jSVphgwZwp133knv3r2BwGw3Dz/8sIbe\nxKAwbCP19fU88sgjXHjhhTz33HN07twZgOnTp+tSn5J0vXv3Dv/BPXXqFHPmzKFdu3baj9gMhWGK\n1dfX88wzz9CuXTvuvfdeamtrATh69KjnyiTbhUYk3HLLLXTt2vWc5zZs2OCjpLSmMEyRQ4cO8atf\n/YoRI0bwta99DYBu3bpx8803M3fuXO0flDYzfPhw7rrrLiZNmhTuhVx66aVce+21LFq0iLNnz3qu\nMD3EDEMzm2pmW8xsm5nNaWKdR4LPrzOzMYlsm01qa2t57rnnmDFjBj179mT27Nl8+OGH4ednz57N\niBEj1E2RNte+fXuuvvpq7r33XiBwEO+VV17hy1/+MgUFBZgZ7733Xvg00FzU7M4qM8sHHgWmANXA\nSjN7yTm3OWKdacAw59xwM/s88BhwWTzbZroXX3yRwsJCli1bxpIlS1iyZMk5O6mHDh3KmDFjGDly\nJD/60Y/SZt/gzp07w2MapW2ky2ce2ld93333sWbNGlavXh2eJm7SpEn06dOHyZMnU1paysSJExk1\nahTvvfcepaWlHqtuG7F+OycA251zFQBmtgCYDkQG2nXAkwDOueVm1s3M+gFD4tg2rTnnOHToELt3\n72bPnj3s3LmTrVu3snXrVrZs2RKeYLMpX//619uo0sRUVFSkxS9mLkm3z7y4uJgrrriCiRMnUllZ\nyRNPPAFATU0Nzz77LM8++ywABQUF9OjRg6uvvprhw4czaNCg8G3AgAEUFhb6/DGSKlYY9gcqIx5X\nAZ+PY53+wPlxbAvA8uXLcc6dc2toaIhrWax1GxoaOHPmDGfOnOH06dOcPn36nPunT5/m2LFjHDly\n5Jzb4cOH2bt3L2fOnInnc2TGjBkMHjyYjh07ApqJWjJDXl4egwYNAmDu3LkcPHiQiooKKioq2Lhx\nI/X19dTU1PDUU081un2XLl3o0aMHPXr0oGfPnvTo0YPOnTvToUMHioqK6NChQ/hWVFREUVER+fn5\n5OfnU1BQEL7f1C20S6ml/yYiVhjGe4Jjq3aCXXbZZa3ZvM3NmDGDXr16sXHjRr74xS9SVlbGqFGj\nfJcl0ipmRq9evejVqxfjxo1j48aNPPDAAyxevJgBAwbw8ssvf2qb2tpaamtrqaioaPuCk62xllbo\nBlwG/Cni8f3AnKh1fg3cHPF4C9A3nm2Dy51uuummWypuzeVb9C1Wy3AVMNzMBgO7gZuAmVHrvATM\nBhaY2WXAx865fWZ2MI5tcc7p0KqIeNdsGDrn6s1sNvAakA/Md85tNrNZwecfd869ambTzGw7cBz4\nRnPbpvKHERFpKdO8ZyIiHs9AMbMbzWyTmZ01s7ERyweb2UkzWxO8/YevGhvTVN3B5+4PDjDfYmbX\n+KoxHmZWZmZVEZ/zVN81NSeTB/CbWYWZrQ9+zit819MUM/tvM9tnZhsilvUwszfM7K9m9rqZdfNZ\nY1OaqD2h77jP0/E2ANcD7zTy3Hbn3Jjg7a42riuWRus2s1EE9ouOAqYC/2Fm6Xy6owP+LeJz/pPv\ngpoSMYB/KoHPd6aZjWx+q7TigNLg5zzBdzHNeILAZxzpn4A3nHMXAn8OPk5HjdWe0Hfc2y+rc26L\nc+6vvt6/pZqpezrwjHOuLjjQfDuBQevpLFMOXoUH/zvn6oDQAP5MkvaftXNuCXA4anH4pIrgv19p\n06Li1ETtkMDnnq4tlyHBZm25mV3hu5g4nU9gYHlIaPB5OrsneD75/HTt/gQ1NbA/UzjgTTNbZWb/\n6LuYBPV1zu0L3t9HYNhcJon7O57SMAzua9jQyO3LzWy2Gyhxzo0Bvgs8bWadU1lntBbW3RivR6ea\n+TmuI3AO+RDgs8Ae4Oc+a40h04/yTQp+n/8WuNvMrvRdUEu4wNHWTPq/SOg7ntKZA5xzV7dgmzPA\nmeD9983sQ2A48H6Sy2uuhoTrJjAZRUnE4wHBZd7E+3OY2X8Bnz69IH1Ef7YlnNsKT2vOuT3Bf/eb\n2UIC3f4lfquK2z4z6+ec22tm5wE1vguKl3MuXGs83/F06SaH+/Vm1iu4wxwzu4BAEDY/I4I/kfsj\nXgJuNrNCMxtCoO50PnJ4XsTD6wkcGEpX4cH/ZlZI4EDVS55riouZdQz1bMysGLiG9P6so70E/H3w\n/t8DL3isJSEJf8cTOV0lmbdgcZXASWAvsDi4/AZgI7AGWA38H181JlJ38LkHCBw42QJ8yXetMX6O\n3wLrgXUEvuB9fdcUo96/BbYGP9/7fdeTQN1DgLXB28Z0rh14hsBuqjPB7/g3gB7Am8BfgdeBbr7r\njLP2OxL9jmvQtYgI6dNNFhHxSmEoIoLCUEQEUBiKiAAKQxERQGEoIgIoDEVEAIWhiAgA/x8cuBbE\nuKCtbwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10aa4ec10>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "fig = plt.figure(figsize=(5,5))\n",
    "axes = fig.add_subplot(111)\n",
    "axes.hist(x_samples, bins=50, normed=True, color=\"gray\");\n",
    "x_range = np.arange(-15, 15, 0.1)\n",
    "x_precision = tau * tau0 / (tau + tau0)\n",
    "axes.plot(x_range, stats.norm.pdf(x_range, 0, 1 / np.sqrt(x_precision)), color='k', linewidth=2);\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "x_pdf = stats.kde.gaussian_kde(x_samples)\n",
    "\n",
    "def bayes_factor_sim(x_obs, x_pdf):\n",
    "    return x_pdf.evaluate(x_obs) / stats.norm.pdf(x_obs, 0, sigma)\n",
    "\n",
    "def bayes_factor_exact(x_obs):\n",
    "    return np.sqrt(tau0 / (tau + tau0)) * np.exp(0.5 * tau**2 / (tau + tau0) * x_obs**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAT8AAAE4CAYAAAAto/QTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XlcVXX+x/HXF8SNRVAEFVDTXNLcsjS3pKkmc0pNyyXN\npcXUNn/V5JSZVNNMVrZqSpaoLWou5W6NKe5ZTmpuaGoibpSCLAKyfX5/QHcA2blwLtzP8/G4D+/h\nfO+5H4/HN9/vPed+jxERlFLK2bhYXYBSSllBw08p5ZQ0/JRSTknDTynllDT8lFJOScNPKeWUigw/\nY8xcY0y0MWZ/Aet9jTHrjTF7jTEHjDGj7V6lUkrZWXF6fmFAn0LWPwHsEZGOQDAw3RhTzQ61KaVU\nuSky/ERkKxBbSJNzgFf2cy/gooik26E2pZQqN/booc0BNhpjzgKewGA7bFMppcqVPU54vAjsFZFG\nQEdgpjHG0w7bVUqpcmOPnl934HUAETlujPkNaAXsztnIGKNfIlZKlQsRMSV9jT16fhHA7QDGGH+y\ngu9Efg1FRB95HlOnTrW8Bkd76D7R/VKSR2kV2fMzxiwEegO+xpgoYCrglh1mocC/gDBjzD6ywvR5\nEYkpdUVKKVUBigw/ERlWxPoLwD12q0gppSqAfsPDYsHBwVaX4HB0n+RP94t9mbKMmUv0RsZIRb2X\nUsp5GGOQUpzwsPybGMaUuGalKhX9pe+YLA8/0INDVV36y91x6Wd+SimnpOGnlHJKGn5KKaek4VdK\nX3zxBXfeeWe5bHv06NFMmTKlXLad17x58+jVq1eFvFd5SE5O5p577sHb25shQ4ZYXY6qRDT8CrFt\n2za6d++Ot7c39erVo2fPnuzenfWV5eHDh/Ptt9+Wy/saY4r9QXlwcDCffvppudRx8uRJXFxc8PT0\nxNPTkwYNGvD444+Tnu44M5YtXbqU33//nZiYGBYvXlzq7VT2XwKq5DT8ChAfH8/dd9/N008/TWxs\nLGfOnGHq1KnUqFGjQt6/uGfAK+JsYlxcHAkJCezfv5+dO3cyc+bMcn/P4oqMjKRly5a4uFh7KGdk\nZFj6/qrkNPwKcPToUYwxDBkyBGMMNWvW5I477qBdu3bA1T0FFxcXZs2aRYsWLfDy8uLll1/m+PHj\ndOvWDW9vb4YOHUpaWlq+r/3z9SdOXD0fRGxsLHfffTd+fn7UrVuXe+65hzNnzgAwefJktm7dyhNP\nPIGnpydPPfUUABEREdxxxx3Uq1eP1q1bs2TJEtv2Ll68SL9+/ahTpw5du3bl+PHjxd4n9evX5447\n7uDQoUO2n73xxhtce+21eHl50bZtW7755hsAUlNTqVevHgcOHLC1/f3333F3d+fixYsArF69mo4d\nO+Lj40OPHj3Yv/9/d0qYNm0agYGBeHl50bp1azZu3HhVPVOnTuW1115j8eLFeHp6EhYWxokTJ/jL\nX/6Cr68v9evXZ8SIEcTFxdleExUVxcCBA/Hz88PX15cnn3ySiIgIxo0bx86dO/H09KRu3bpAVuiP\nHDkSPz8/mjZtyuuvv277pTRv3jx69OjBM888g6+vL6+88kqx96NyEBU484Lkp6CfZ6+036OE4uPj\npV69ejJq1ChZt26dxMTE5FofFhYmPXv2tC0bY2TAgAGSkJAgBw8elOrVq8utt94qv/32m8TFxUmb\nNm1k/vz5+b72z9cfP35cRERGjx4tL730koiIXLx4UZYvXy7JycmSkJAg999/vwwYMMD2uuDgYPn0\n009ty4mJiRIYGCjz5s2TjIwM2bNnj/j6+sqhQ4dERGTIkCEyZMgQSUpKkgMHDkhAQID06tUr333w\n22+/iTFG0tPTRUTkzJkz0qFDBwkLC7O1WbJkiZw7d05ERBYvXizu7u5y/vx5ERGZMGGCTJo0ydb2\nvffek379+omIyM8//yx+fn7y448/SmZmpsyfP1+aNm0qqampEhERIUFBQbbtRkZG2vZNXiEhIfLg\ngw/alo8dOyYbNmyQ1NRU+eOPP+SWW26RiRMniohIenq6tG/fXp555hlJSkqSlJQU2b59u4iIzJs3\n76p/kwcffFAGDBggiYmJcvLkSWnZsqVtX4eFhUm1atVkxowZkpGRIcnJyfnWV+jxrewiex+XPJNK\n86JSvVElCz8RkcOHD8vo0aMlMDBQqlWrJv369ZPo6GgRyT/8duzYYVvu3LmzvPnmm7blZ5991vaf\nsCThl9eePXvEx8fHthwcHCyffPKJbXnRokVXhdnYsWPllVdekfT0dHFzc5MjR47Y1r344otX1fKn\nP8PP29tbvL29xRgjPXr0kPj4+Hzbi4h07NhRVqxYISIiP/zwgzRu3DjXPlmyZImIiIwbN06mTJmS\n67WtWrWSzZs3y7Fjx8TPz88WYoWZOnWqjBgxosD1X3/9tXTq1ElERHbs2CH169eXjIyMq9rl/TdJ\nT0+X6tWry+HDh20/Cw0NleDgYFv7nH+3gmj4lb/Shp9jD3vtGX+l0Lp1a8LCwoiKiuLAgQOcPXuW\niRMnFtje39/f9rxWrVq5lmvWrEliYmKJa0hKSuKxxx6jadOm1KlTh969exMXF/fnLxQg9+d+kZGR\n7Nq1Cx8fH9vjyy+/JDo6mgsXLpCenk5QUJCtfePGjYus4eLFi8TGxpKUlET37t1zneVesGABnTp1\nsr3XgQMHbMParl27UqtWLcLDw4mIiOD48eP069fPVuf06dNz1Xn69GnOnTtH8+bNee+99wgJCcHf\n359hw4Zx7ty5Yu2v6Ohohg4dSmBgIHXq1OHBBx+01RMVFUWTJk2K9fnghQsXSEtLo0mTJrn21Z8f\nOQC59qOqfBw7/BxIq1atGDVqVK7PsEoiZ0C5u7uTlJRkWz5//nyB7adPn87Ro0f58ccfiYuLY/Pm\nzTl701ed8GjcuDG9e/cmNjbW9khISGDmzJn4+vpSrVo1Tp06ZWuf83lRatasyahRo/jhhx+IiYkh\nMjKSsWPHMnPmTGJiYoiNjeX666/PFcyjRo3i888/57PPPuP++++nevXqtjonT56cq87ExETb5SrD\nhg1j69atREZGYoxh0qRJRe5XgBdffBFXV1cOHDhAXFwcn332GZmZmUBWWJ06dSrfkxN5t+Pr64ub\nmxsnT57Mta8CAwMLfI2qXDT8CnDkyBHeeecd22/6qKgoFi5cSLdu3Yq9jZwhkPN5hw4dOHjwIPv2\n7SMlJYWQkJCrXvdn+8TERGrVqkWdOnWIiYm56oN1f3//XCct7r77bo4ePcrnn39OWloaaWlp/PTT\nT0RERODq6srAgQMJCQkhOTmZQ4cOMX/+/CL/E/9Zy5UrV/jss89o2LAhdevW5fLlyxhj8PX1JTMz\nk7CwsKt+OYwYMYLly5fzxRdfMHLkSNvPH330UWbPns2PP/6IiHD58mXWrFlDYmIiR48eZePGjVy5\ncoUaNWpQs2ZNXF1di9zHf+4vd3d3vLy8OHPmDG+99ZZtXZcuXWjYsCH/+Mc/SEpKIiUlhR07dtj2\n4+nTp20npVxdXRk8eDCTJ08mMTGRyMhI3n33XUaMGFHovlKVh4ZfATw9Pdm1axddu3bFw8ODbt26\n0b59e6ZPnw5cfS1efgGSd/2fyy1btuTll1/m9ttvp1WrVvTq1avAthMnTiQ5ORlfX1+6d+/OXXfd\nlavt008/zdKlS6lbty4TJ07Ew8OD7777jkWLFhEQEEDDhg154YUXSE1NBWDGjBkkJibSoEEDHnro\nIR566KEi94W3t7ftOr9du3axcuVKANq0acOzzz5Lt27daNCgAQcOHKBnz565XhsUFMQNN9yAi4tL\nrnWdO3dmzpw5PPHEE9StW5cWLVqwYMECICtkX3jhBerXr0/Dhg25cOEC//73v/OtLe+/w9SpU/n5\n55+pU6cO99xzD4MGDbKtd3V1ZdWqVRw7dozGjRsTFBTEV199BcBtt91G27ZtadCgAX5+fgB8+OGH\nuLu706xZM3r16sXw4cMZM2ZMvu+rKh/L5/PLnourQmpQ1nj44YcJCAjg1VdftbqUCqfHd/mrtPP5\nqart5MmTLF++nL1791pdilK56LBXlZspU6bQrl07nn/++VxnTZVyBDrsVaoc6fFd/ko77NWen1LK\nKWn4KaWckoafUsopafgppZyShp9SyikVGX7GmLnGmGhjzP5C2gQbY/YYYw4YY8LtWqGqdL7++muC\ngoLw9PRk3759VpejVL6K0/MLA/oUtNIY4w3MBO4RkeuB++xUm+WaNm1K7dq1bdO455wwtDyEh4cX\nOVPI6NGjqVGjBp6ennh5eXHjjTeyZcuWcqupNJ577jk++ugjEhIS6NChQ6m307Rp03wnMVXKHooM\nPxHZCsQW0uQBYJmInM5uf8FOtVnOGMPq1atJSEiwPT744APLa5o0aRIJCQnEx8czfvx4Bg4c6DDX\nkokIp06dok2bNmXeVlmukXOk+4wox2SPz/xaAHWNMZuMMbuNMQ/aYZsOb/z48dx33/86uZMmTeL2\n228HCp96HiAmJoYxY8YQEBBA3bp1GThwIElJSdx1112cPXvW1qvLb6qrvIYNG0ZMTAzR0dEAHD9+\nvMBp3N96661cNQM89dRTtjkK4+LiePjhh2nUqBGBgYFMmTLFNh3UsWPH6N27N97e3tSvX5+hQ4de\nVcuVK1fw9PQkIyODDh060KJFC6Dgqe7/NGfOHNq0aWNbv2fPHh588EFOnTrFPffcg6enJ2+//TYA\nK1eupG3btvj4+HDrrbcSERFh207Tpk158803ad++PZ6enrbaleObNw8+/7yC37Q4M54CTYH9Bayb\nAewAagH1gKNAi3zaFTYLq0Nq2rSpbNiwId91SUlJ0rJlS5k3b55s2bJFfH195cyZMyJS9NTzffv2\nlaFDh8qlS5ckLS1NtmzZIiIi4eHhEhgYWGhNOWd5Tk9Pl1mzZknz5s0lMzNTRAqfxv3cuXPi7u4u\nly5dEhGRtLQ08fPzk59//llERAYMGCDjxo2TpKQk+f3336VLly4SGhoqIiJDhw6Vf/3rXyIicuXK\nFdv07/nJOSu1SOFT3X/11VcSEBAgu3fvttUfGRlp2//ff/+9bTtHjhwRd3d32bBhg6Snp8ubb74p\n1157raSlpYmISJMmTaRTp05y+vRpSUlJKXQ/VhRHPr4dSc+eImvWlO61lOc09kWE3yQgJMfyJ8B9\n+bSTqVOn2h6bNm3KWXghfzH7PEqjSZMm4uHhYZvG3dvbO9eU8bt27RIfHx9p0qSJLFq0qMDt5Jx6\n/uzZs+Li4mILoJw2bdpUZPiNGjVKatasKd7e3lKrVi2pVauWfPnllwW2zzmNu4hInz59ZM6cOSIi\nsmrVKmnTpo2IiJw/f15q1KiR614UX375pdx6660iIjJy5EgZO3asnD59utD6RK4Ov7w6duwoK1eu\nFBGRv/71r/LBBx/k2y5v+L366qsyZMgQ23JmZqYEBATI5s2bbe1z3l/EEWj4Fe3sWRFvj1RJOXmu\nWO03bdqUK0usDL/WwAbAFagN7Afa5NMu37+IIx8cef/z5adz587SpEmTXPeFuHz5sowdO1aaNGki\nXl5e4uXlJS4uLpKZmSm7du0SX1/ffLdVnPAbPXp0rntfHDhwQBo1aiTr1q0TkawQGzJkiAQEBIiX\nl5d4eHjkutfEwoULpXfv3iKSdTOjN954Q0SygtzFxSVX0Ht5ecn1119v2+6jjz4qjRo1krZt28rc\nuXMLrDFv+M2fP186duxo2261atVsr2/Tpo2sKeBXft79P378ePn73/+eq83NN99sC//CeupWceTj\n21HMfD9NhrstFjFGJLvXXxKlDb/iXOqyMHtY28oYE2WMecgY85gx5rHsRIsA1gO/ALuAOSJyqOAt\nVh0zZ84kNTWVRo0a8eabb9p+XtjU80FBQcTExOS6neKfSjM5Ztu2benRowdr164FCp/GHaB///78\n8ssvHDhwgDVr1jB8+HAga9LRGjVq2O7XERsbS1xcnO12kv7+/nz88cecOXOG0NBQJkyYkO+tNvMq\naqr7oKAgjh07lu9r8+6PRo0aERkZaVsWEaKioggICCjwNcrxLZt7iUFpC+H666EY95Sxl+Kc7R0m\nIo1EpLqIBInIXBEJFZHQHG3eFpG2ItJORKw9HWpnf/4nzevo0aNMmTKFL774ggULFvDmm2/armkr\nbOr5hg0bctdddzFhwgQuXbpEWlqa7VIVf39/Ll68SHx8fKH15KwpIiKCbdu20bZtW9t7FzSNO2Td\nWGnQoEE88MADdO3a1XZPioYNG/LXv/6VZ555hoSEBDIzMzl+/LittiVLlnD69Gkga2ZnY0yxbgRU\n1FT3jzzyCG+//TY///wzIsKxY8ds9xXJO0X/4MGDWbNmDRs3biQtLY3p06dTs2ZNunfvXmQdyjFd\nuAD/PVybPqyH+++v2DcvTXexNA8q6bC3Vq1a4uHhYXsMHDhQ0tPTpUuXLjJt2jRb21mzZkm7du0k\nNTVVzp49K8HBweLh4SGtWrWS0NBQcXFxsQ2NY2JiZNSoUeLv7y8+Pj4yaNAg23YeeughqVevnvj4\n+NhOEuQ0evRoqV69unh4eIi7u7s0btxYJk+ebFt/8OBB6dy5s3h4eEinTp1k+vTpEhQUlGsbW7du\nFWOMzJs3L9fP4+LiZPz48RIYGCh16tSRTp06yeLFi0VE5Pnnn5eAgADx8PCQ5s2b2z43zI+Li0uu\nYe/kyZOlbt264uvrK88888xV9xqePXu2tGrVSjw8PKRdu3ayd+9eERFZsWKFNG7cWLy9vWX69Oki\nkvUZZps2baROnToSHBxsux/xn/9eRX1MUdEc+fh2BJ+Epsn9bl9nfQKX49+yJCjlsFfn83NCUVFR\ntG7dmujoaDw8PKwup0rT47twd930B6N3P8GQtgehDHdGFJ3PTxUlMzOT6dOnM2zYMA0+ZalLl2D7\nPg/6shbuq/gvhuk9PJzI5cuX8ff355prrmH9+vVWl6Oc3Kqv07mVcDxJrPjP+9Dwcyru7u4kJiZa\nXYZSACybE8N9aV/CdddB9gm7iqTDXqVUhUtIgE27PbmHVZYMeUHDTyllgbWrMujODryJs2TICxp+\nSikLLAu9kHVhc6tWWRc3W8AhPvPTq/KVch7JyfDtD17MZAXc9xhY9P/f8vDTa6CUci7frsugs+ym\nPhcsG/KCDnuVUhVs2aw/uC9tIbRoAe3bW1aHhp9SqsJcuQKrt3hxL19nneW18CMvDT+lVIX5/rsM\n2sp+GnLe0iEvaPgppSrQstm/MyhtETRrBh07WlqLhp9SqkKkp8PKjR4MZHlWr8/iqzw0/JRSFWLz\npkyuyThGE05ZPuQFDT+lVAVZ9tH5rCFv06Zwww1Wl2P9dX5KqaovIwOWf+fBNpY5xJAXtOenlKoA\nO7Zl4p96mms57hBDXtDwU0pVgGUzzzMofRE0aQI33mh1OYAOe5VS5UwElq+vxTqWWX5hc07a81NK\nlaufdmXinnyBNhxymCEvaPgppcrZ0hnnGJS+GBMUBF26WF2OjYafUqrciMCy1TUY5GBDXtDwU0qV\no317BUlMoiN7HWrICxp+SqlytOzDswzKWIwJCICuXa0uJ5ciw88YM9cYE22M2V9Eu5uMMenGmIH2\nK08pVZktW1GN+1iaNeR1cay+VnGqCQP6FNbAGOMKTAPWA44zqFdKWebQQSEhLoOb+MmyO7QVpsjw\nE5GtQGwRzZ4ElgJ/2KMopVTlt+zDMwzMWIJLo4bQvbvV5VylzP1QY0wA0B+Ylf0jvSmHUoplX7tk\nneUdNMjhhrxgnxMe7wH/kKw7ERl02KuU0zt+TDh/0Y0ebHfIIS/Y5+ttnYFF2bef9AXuMsakicjK\nvA1DQkJsz4ODgwkODrbD2yulHM2yD88yIGMVrg38oEcPu247PDyc8PDwMm/HFOfWkcaYpsAqEWlX\nRLuw7HbL81kneptKpZxDl0ZR/OvcGG5/vDXMmFGu72WMQURKPOIssudnjFkI9AZ8jTFRwFTADUBE\nQkv6hkqpqu1UpHAi2p3ebIb7XrK6nAIVGX4iMqy4GxORMWUrRylV2S2fcZZ+md/i5l8PevWyupwC\nOd4pGKVUpbZsUVrWWd6BA8HV1epyCqTz+Sml7ObcWeHAWR9uZwPct87qcgql4aeUsptlH56lb2Y4\nNerXgVtusbqcQumwVyllFyLwSZgLo5kH994L1Ry7b6Xhp5Syi10/CIkxqdzG9w43fVV+NPyUUnYR\n+mo0Y9Nm4tLAHyrBFxg0/JRSZRYbC99878EYwuCRRxx+yAt6wkMpZQeffRRPn/R11HeJgbFjrS6n\nWIr19Ta7vJF+vU2pKkkErm/wBx/9fh+9+/vAN99U6PuX29fblFKqMNs2Z5AZE8ctbIHx660up9j0\nMz+lVJmEvnKOsekzMc2bwx13WF1OsWn4KaVK7cIFWL3Nm1HMh8cec8hJSwuiw16lVKnNf/ci/dLX\nULdGEoypXPOaaPgppUpFBEJnZTKP2TB4MPj6Wl1SiVSePqpSyqFs+jaVGvF/0I2dMH681eWUmIaf\nUqpUQkPOMS5jJqZDB7j5ZqvLKTENP6VUiUVHw3e76zKCz7N6faby3bdMP/NTSpVY2L/PMzBjDXU8\nBYYPt7qcUtHwU0qVSGYmfDy3GouZDSNHgoeH1SWVioafUqpE/rMiCe/Lp7mR3TB+ntXllJp+5qeU\nKpHQV88zLvMjTK9e0Lat1eWUmoafUqrYzp4RwvfXYxgLK+XlLTnpsFcpVWyfhkQxJGMNnvVrZd2d\nrRLT8FNKFUtGBsz5sjYrCc2asLRGDatLKhMNP6VUsaxbFEfD5BN0NL/A2Iqds6886Gd+SqliCX39\nD8bJLOjbF5o2tbqcMisy/Iwxc40x0caY/QWsH26M2WeM+cUYs90Y097+ZSqlrHTqZCY7jtRjCIsr\n/YmOPxWn5xcG9Clk/QngFhFpD7wGfGyPwpRSjuOTF08wPPMzajfxgz6FxUHlUeRnfiKy1RjTtJD1\nO3Ms7gICy16WUspRpKXBp1/78B2hMG4cuLpaXZJd2PuEx8PAWjtvUyllodVhf3BNyhHauv0KDz1k\ndTl2Y7fwM8bcCjwE9CioTUhIiO15cHAwwZXgxsZKObvQabGMYzbcdx/4+VldDuHh4YSHh5d5O8W6\ndWX2sHeViLQrYH17YDnQR0SOFdBGb12pVCVzIiKVrm0SiJIAam75D/TqZXVJV7Hs1pXGmMZkBd+I\ngoJPKVU5zfnHMUbKOmq2vRZ69rS6HLsqMvyMMQuB3oCvMSYKmAq4AYhIKPAy4APMMlkTGqaJSJdy\nq1gpVSFSUyFsbQM28zFMeLpSTlhamGINe+3yRjrsVapS+eqd08x+9igb3fvB2bPg5WV1Sfkq7bBX\nv+GhlMpX6LuXeYxQGDHCYYOvLLTnp5S6ytG9SfS6IZEoCaT63p+gQwerSyqQZSc8lFJVz7sTT/KI\nfE31bjc6dPCVhYafUiqXU5HCV9sacYR3YcJ7VpdTbjT8lFK5/HvcScZmLMK3gVvWhc1VlIafUsrm\n1MlMvvqPD0eYDi9OhZo1rS6p3Gj4KaVs3hj3G2MzvsI3oCY8+qjV5ZQrDT+lFABRkZks3lAvu9f3\napXu9YGGn1Iq278fPcGjGUvxDaoNDz9sdTnlTsNPKUXUyQwWb/Qlgunw0uuV/uZExaHhp5Ti348e\n59GM5dRv6gGjR1tdToXQ8FPKyUX9ls7ijfWzen1TpkH16laXVCE0/JRycv9+5DiPZn5N/eZ14MEH\nrS6nwmj4KeXEok6ksTjcL6vX9/J0cHOzuqQKo+GnlBN745FjPJK5kvot68IDD1hdToXS8FPKSUUd\nT2XR5gZE8Da8/B5Uc644cK6/rVLK5o2Hf+WRzNXUb+0LQ4daXU6F0/BTyglFHbvCoq2Nsnp9ITOq\nzL14S0LDTykn9MbDR3kkcy312/rD/fdbXY4ldCZnpZzM6WMpdGiZTIS0pP7S2TBokNUllYnO5KyU\nKpY3HjrCI7Ke+h0C4N57rS7HMtrzU8qJnP41mQ6tkomQVtT/5hPo39/qkspMe35KqSK9MeYID8t3\n1L+hMfTrZ3U5ltKen1JO4vTRJNq3vkKEtMJvdRj87W9Wl2QX2vNTShXqjTERPCIb8OtyDfTta3U5\nltOen1JO4HREIu3bpGX1+tYtgD59rC7Jbkrb83MpxobnGmOijTH7C2nzgTHmV2PMPmNMp5IWoZQq\nX2+MOcIjMge/btfCnXdaXY5DKM6wNwz4EFiQ30pjTF/gWhFpYYzpCswCbrZfiUqpsjgTkcCXu5oR\nwV3w6pdgStxJqpKK7PmJyFYgtpAm/YD52W13Ad7GGH/7lKeUKqs3Rh/mYfkEv16t4bbbrC7HYdjj\nhEcAEJVj+TQQCETbYdtKqTI4cyiOL35sQQR3w6tfaa8vB3ud7c27R/M9sxESEmJ7HhwcTHBwsJ3e\nXimVn0lDT/KYrMXv1uuhivx/Cw8PJzw8vMzbKdbZXmNMU2CViLTLZ91sIFxEFmUvRwC9RSQ6Tzs9\n26tUBdq4MJoxw69wSK7Dfct66NXL6pLKRbmd7S2GlcDI7CJuBi7lDT6lVMW6kiJMeDSVD+RJ3Aff\nXWWDryyKHPYaYxYCvQFfY0wUMBVwAxCRUBFZa4zpa4w5BlwGxpRnwUqpor390EFaXj5Bf69weC/C\n6nIckl7krFQVc2JfAl1uSOO/mZ1oMnMSTJhgdUnlqrTDXg0/paoQEfjbNQfpHbmASV3CYceOKj9L\ns363VynF8jePERkJ/+fyAYTurPLBVxba81OqikiITaeN/wW+SBvMLc91hbfesrqkCqE9P6WcXMi9\n+7gtbT+3NI6EkHVWl+PwNPyUqgL2fXuezzY35iB9YGYYuLtbXZLD02GvUpVcZib09I1gdOw7jB0U\nA0uXWl1ShdJhr1JO6tMn9iCxKTzisRjeP2R1OZWGhp9SldgfvyUyOTSI/3AbLv/6JwQEWF1SpaHD\nXqUqsTFtd+FzaDvv3LgQfvjBKS9t0WGvUk5m66dH2HCoEYfMKxC6ySmDryzsMbGBUqqCpSZnMP6J\narzL/+GXWIpEAAAVB0lEQVQ58WG44QarS6p0tOenVCX07pAfCEqJZ1DALnh1ntXlVEoafkpVMpE/\nnOOt1a3ZRRfMzBng4WF1SZWShp9SlcxTA6OYKCtpPqA99O9vdTmVloafUpXIypd+5Mg5L75ynwUf\n7LW6nEpNw0+pSuLy75d56o2GzGUUNf45BYKCrC6pUtPwU6qSeO2eXfTIOMdfboiDJ56wupxKT8NP\nqUpg/9IjfPpjO/abERC6Eqrpf92y0uv8lHJwiRdSGDyiOm/zHA2evB9uvNHqkqoE/XqbUg5MBEa1\n/gHXo4cJu+Y12LsXvLysLsuh6NfblKqCwh7byX+PevGj2zOwZIMGnx1p+CnloPZ/c5xJc65lM71x\n//AN6NzZ6pKqFB32KuWAEs8nclOTaF5IfYWRwzPhs8/AlHhk5xR02KtUFSGZwvhue+mWepSRbf4L\ns3dp8JUDDT+lHMzcUZvZc9KXH2tPgqVb9Lu75UTDTykH8svCg/zj87Zs4RZqz3kfrrvO6pKqrCKv\n8zPG9DHGRBhjfjXGTMpnva8xZr0xZq8x5oAxZnS5VKpUFZdwKpbBI2vyDs9w3fhb4YEHrC6pSiv0\nhIcxxhU4AtwOnAF+AoaJyOEcbUKAGiLygjHGN7u9v4ik59mWnvBQqgCSkcmIJluodeYYn3SeDdu3\nQ40aVpdVKZTXCY8uwDEROZn9JouA/sDhHG3OAe2zn3sBF/MGn1KqcJ8M/o5fzgSwq84DsESDryIU\nFX4BQFSO5dNA1zxt5gAbjTFnAU9gsP3KU6rq2zf3v7y4vDNb6UXtz0LhmmusLskpFBV+xRmnvgjs\nFZFgY0xz4D/GmA4ikpC3YUhIiO15cHAwwcHBJShVqaon4fjvDB7rzXtMpPWkAXDPPVaX5PDCw8MJ\nDw8v83aK+szvZiBERPpkL78AZIrItBxt1gKvi8j27OXvgUkisjvPtvQzP6VykPQMhgeG4x59gjm3\nfA7ff6+ztZRCeX3mtxtoYYxpCpwFhgDD8rSJIOuEyHZjjD/QCjhR0kKUcjZz+q/mQPQ17PJ7GBbt\n1OCrYIXubRFJN8Y8AXwLuAKfishhY8xj2etDgX8BYcaYfWRdOvO8iMSUc91KVWr7Zm5j8trubDO3\nUGtRGDRsaHVJTke/26tUBYs/GMWN7a8QkvkyD7x+Pbz4otUlVWqlHfZq+ClVgeRKKsMabaZOzAlC\n+66EVavARecULgud2EApRyfCm71XExHTjJ2Bj8OCnRp8FtLwU6qCzO33DbN3dWKb21+otXQR1Ktn\ndUlOTcNPqQqw4qEVTF59M5td/kLAV+9C17zfFVAVTcNPqXK2ZdIaHg3rxlr60jLsBRgwwOqSFHr3\nNqXK1b7pG7j/zRtZyDBufH8kjBxpdUkqm4afUuXkxIJt9H3uOmbyOLeF3AJPPWV1SSoHvdRFqXJw\nfvVuevbz4Tl5i3FP1YD33tOp6MuJXuenlIOI236A3rcIgzK/YsqoKJg7Vy9pKUcafko5gJSDx7mz\n0+90SPuJ9/tvwixdot/ZLWcafkpZLD3yDPddd5DayRf4/Na5uKxdDTVrWl1Wlaff8FDKQnLhIo91\n/IGUZHe+unEmLivWa/A5OA0/pcoqIYEX263i4KXWfH/dk1Rfvx48Pa2uShVBw0+pskhJ4Z1On7Hi\nfDBbG4/AfcNq/dpaJaHhp1RppaezoPts3j9+L9vqD6TexiXQqJHVVali0vBTqjQyM1nz1/d5fs9w\nwr36E/T9PGje3OqqVAlo+ClVUunp/Kffh4zZ9CCra95P62/fh3btrK5KlZBeealUSSQm8sVN7zFi\n3QMsdxtKl5Uvwc03W12VKgXt+SlVXNHRTL/xS94/PZiNXvfSdt3b0L271VWpUtLwU6oYMg8f4e83\nb2V9/F/ZHjSMoP+EQatWVpelykDDT6kipG7azpg+ZzmV2optnZ7CZ91y8Pe3uixVRvqZn1KFSFjw\nNX+7LYWk1Gp81/d9fLau1OCrIjT8lCpA9KuhBI9qTHP5laWPbaDWikXg7m51WcpONPyUyisjg2Oj\n/0n3qbfTnxXM+tclXGfN0NlZqhj911Qqp+Rkdv9tKv02TSTE9TXGLugFDzxgdVWqHBTZ8zPG9DHG\nRBhjfjXGTCqgTbAxZo8x5oAxJtzuVSpVES5c4LvOL9B303PMqv0cYzcM0eCrwgqdz88Y4wocAW4H\nzgA/AcNE5HCONt7AduBOETltjPEVkQv5bEvn81OO6/hxPu8xi+ein2OZ3wR6bHwN2ra1uipVDOU1\nn18X4JiInMx+k0VAf+BwjjYPAMtE5DRAfsGnlCOTXT8y/ba1fHj5STa2mkCbjTN0ggInUNSwNwCI\nyrF8OvtnObUA6hpjNhljdhtjHrRngUqVp4zFS3i2xw/Mu3wf23u9QJsf52nwOYmien7FGae6ATcA\ntwG1gZ3GmB9E5NeyFqdUubl8mejHXmbEF31Ix5etD8zGZ958cHOzujJVQYoKvzNAUI7lILJ6fzlF\nARdEJBlINsZsAToAV4VfSEiI7XlwcDDBwcElr1ipstq7l+/7vc/IqH/ykOsCpr7tSbWnP9RbS1YS\n4eHhhIeHl3k7RZ3wqEbWCY/bgLPAj1x9wqM1MAO4E6gB7AKGiMihPNvSEx7KWiKkv/MBrz6fyCeZ\nY1jQeAq3r3oa2re3ujJVBuVywkNE0o0xTwDfAq7ApyJy2BjzWPb6UBGJMMasB34BMoE5eYNPKcv9\n/jtnhj7LA5sewY00fh75Pg1mfQi1a1tdmbKI3rpSVX3ffce6IfN46NJ0nqg1l38saIPrffdaXZWy\nE71vr1J5paaSNuklXnqvHl/yAF90eItbVj8PgYFWV6bsSO/bq1ROR48SOfD/GHZwMt7E8fOLy6j/\n6rvg6mp1ZcpB6MQGqmoRgbAwVrSfQpeDc7nXZzOrt/tQ//WJGnwqFx32qqrj0iVSH32c55fexDcM\nYOEdYXRb8gzUqWN1ZaoclXbYqz0/VTVs3Mjxtv3osXQika7N2fPRD3T7NkSDTxVIw09VbqdPk3H/\nUGbdtoRuZ5cyMiic5Yevw2f8UL1oWRVKT3ioyik1Fd55h59fWcW4lHep6ZLKxse/4fq3n4bq1a2u\nTlUC2vNTlc933xHXphtPvVCbvinLmHDjT2w+0ZjrPxirwaeKTcNPVR6Rkci9A1l856e0Ob6S5DoN\nOLg0gtE/PY5p0tjq6lQlo8Ne5fhSUuCtt/j19a94/Mp0zptGfDV+Mz3evU97eqrUtOenHNuaNaS0\nuYFXXk6n25VN3Nnxd/57woceMx/Q4FNloj0/5ZiOH4eJE/nP6hQmsIJ2npHs+eRXggaPsLoyVUVo\n+CnHcvkyTJvG2Wmf8Uzqv9llbmbGw3v520d/04lGlV3psFc5hsREmDaNtCbX8sFrl2if+hPNrq/N\nwWM1+ducARp8yu6056esFR8PM2aQPv19Povpy2ts51qvP9g6M5LrRvSzujpVhWn4KWvExcEHH5Dx\nzvssvNSHV9hGoFcC81++QK9nuui3M1S50/BTFSs2Ft5/n4x3P2BJ/F8JYRv1va7w8csJ3PpMJw09\nVWE0/FTFiImBd98l8/0PWZ5wOyFswdPT8OGUFG5/tgPGRUNPVSwNP1W+LlyAd99FPviQFYl/YSqb\nqe5Zg7deSqPP36/Xjp6yjIafKh9nzsCHHyIfzmBNUjAvE454evHPyanc/XxrDT1lOQ0/ZT8isHEj\nfPQR8s0K1mfewVQ2kuJRn1cmpzJg0rUaesphaPipsouNhfnzYdYs4o+eYz6j+Ij9uHnVZso/Uhk0\n6Rpc9IpS5WA0/FTp7d4Ns2bBwoUcTL6GmTzNIvMAd7Q5w8ev1afnAF/t6SmHpeGnSiYpCRYvhlmz\nSPtpD98wgJms5Wj1dowd8DsH3vKgUeO2VlepVJE0/FTxHD0Ks2dDWBjnLtVkDo/ysfmGZo2SeeLv\ntbl3Qj3c3OpZXaVSxabhpwoWHw8rVsCCBciGDWynBzOYxbeufRnS6xxrp9WlfZeaVlepVKkUGX7G\nmD7Ae4Ar8ImITCug3U3ATmCwiCy3a5Wq4iQlwZo1sGgRrFlD7JVafMVgPjL7SPHyY8Ijacx+yQtv\nby+rK1WqTAoNP2OMKzADuB04A/xkjFkpIofzaTcNWA/oR9yVzZUr8N13WYG3YgUJlw0r6ccilrDF\n9VbuaB/N2y/5cdsATz1rq6qMonp+XYBjInISwBizCOgPHM7T7klgKXCTvQtU5SQ9HTZtygq85ctJ\nvpTCOu5iEXP51rUvPVv+wbBxdfhitAdeXh5WV6uU3RUVfgFAVI7l00DXnA2MMQFkBeJfyAo/sWeB\nyo4yM2HHjqzAW7KEtN9j+A93sIj3WeU6gBuaxDD0YXdmPeZBvXoaeKpqKyr8ihNk7wH/EBExxhh0\n2OtYEhPh++9h3TpYs4aM02fZTG8W8RrLXe6nVaMEho50Y9oTXjRsqJ/jKedRVPidAYJyLAeR1fvL\nqTOwKCv38AXuMsakicjKvBsLCQmxPQ8ODiY4OLjkFavCicDhw1lht24dbNlCUlo1NtObtTzPUpch\nNKqfytChhv9O9KZJUx+rK1aqRMLDwwkPDy/zdoxIwZ07Y0w14AhwG3AW+BEYlveER472YcCq/M72\nGmOksPdSZZCYmPWd2uzAk8hIDnMd6+nDt9zJDpee3BB0gT5/c2XQk41o2VrPWqiqwxiDiJR4xFlo\nz09E0o0xTwDfknWpy6cictgY81j2+tBSVavKRgQiInL17i6l1uJ7bmM9k/nW5S5catWgz82XGDey\nAV/1d6dOHXerq1bKoRTa87PrG2nPr/T+DLtt22DrVtiyhczIU/zMDaynD+vpwz6XTvRsdo4+/apz\n58OBtLrORb9Xq5xCaXt+Gn6OKDUVfv75f2G3fTvpFy/xC+3ZQXe204Pvze34eqbSp0cCfcY0pNfd\ndahVy+rClap4Gn6VWXw87NyZFXbbtsGuXVxIrs1Outkeu81NNPGMpVvbeLrd4cHtowNpfI2r1ZUr\nZTkNv8riyhU4dAj27Mnq3e3YQcbe/RyS1uygOzvpxg66E+3SkK4No+h2YxrdB/jRtX8DvH10HKtU\nXhp+jig2Fvbtg717s8Ju717SDx7h14xrOMD1/EJ7dtGVXXSlQc1LdGv+B91uqUb3IY1p07Murtqx\nU6pIGn5WEoGoqFwhJ3v2cjYylf20y/U4Qisauf1BO79o2rVKo8vtXtw8vDm+jWtb/bdQqlLS8KsI\nly7Br79mPY4etT2POxrNwbiAq4KuGum084qkXeN42nVwod2tvrS9pxnufnrZiVL2ouFnLwkJ/wu4\n7EfakROcOprCb7F1OEEz2+M3ruEEzbhCDa5z/TW7N5dKu27utOsTiH/35lBNp0xUqjxp+BUlMxMu\nXsy6peLZs7n+TIq6yPnIK0SfyyTyktdV4XaWRjTkHM1cImlW5yLNGqVwTXMXmrVzp1nX+vjedA3G\n3w+9sE6piud84ZecnDUMze8RG0vymRiiT1zm/Ol0os8L52OqE53pSzT+nKdBrj/TqYY/0fgTTROX\nKJrVvUSzwFSuaeFGsw6eNO7aELfrroVGjTTglHIwlSL8MpctJyMljbTkdNKTs/9MSSctJYP0KxlZ\ny1cySLuSafszKSGD+EuZxMdDfIIhPsmVhGQ34jNqE4/XVY8EPInHizTc8CeaBpz/358142ngcwX/\n+pk0CHDFv0lN/Ft44dXMFxMYkBVuDRqgM3YqVXlUivADwZV0qpGOG2nF+rM2SflEXDyersl4uWfg\n5ZGJlxd41TF41q2GV73qePrWwDOwzv8CLSAAGjZEvwKhVNVTKcIvo/+9uNRwg+rV//eoUSP3cs6f\nu7mBhwf4+IC3d+5HzZo6BFVKVY7wqxRne5VSlUppw08/3FJKOSUNP6WUU9LwU0o5JQ0/pZRT0vBT\nSjklDT+llFPS8FNKOSUNP6WUU9LwU0o5JQ0/pZRT0vBTSjklDT+llFPS8FNKOaVihZ8xpo8xJsIY\n86sxZlI+64cbY/YZY34xxmw3xrS3f6lKKWU/RYafMcYVmAH0AdoAw4wx1+VpdgK4RUTaA68BH9u7\n0KoqPDzc6hIcju6T/Ol+sa/i9Py6AMdE5KSIpAGLgP45G4jIThGJy17cBQTat8yqSw/oq+k+yZ/u\nF/sqTvgFAFE5lk9n/6wgDwNry1KUUkqVt+LcVLbY0y8bY24FHgJ6lLoipZSqAEVOY2+MuRkIEZE+\n2csvAJkiMi1Pu/bAcqCPiBzLZzs6h71SqlyUZhr74vT8dgMtjDFNgbPAEGBYzgbGmMZkBd+I/IKv\ntMUppVR5KTL8RCTdGPME8C3gCnwqIoeNMY9lrw8FXgZ8gFkm645qaSLSpfzKVkqpsqmwu7cppZQj\nsfs3PIq6IDq7zQfZ6/cZYzrZuwZHU4yLxIONMXHGmD3Zj5esqLMiGWPmGmOijTH7C2njVMcJFL1f\nnPFYATDGBBljNhljDhpjDhhjniqgXfGPGRGx24OsYfExoCngBuwFrsvTpi+wNvt5V+AHe9bgaI9i\n7pNgYKXVtVbwfukFdAL2F7DeqY6TEuwXpztWsv/eDYCO2c89gCNlzRZ79/yKvCAa6AfMBxCRXYC3\nMcbfznU4kuLsEwCnOiEkIluB2EKaONtxAhRrv4CTHSsAInJeRPZmP08EDgON8jQr0TFj7/ArzgXR\n+bWpyt8IKc4+EaB7dld9rTGmTYVV57ic7TgpLqc/VrKvPOlE1rfJcirRMVOcS11KorhnT/L+5qrK\nZ12K83f7GQgSkSRjzF3AN0DL8i2rUnCm46S4nPpYMcZ4AEuBp7N7gFc1ybNc4DFj757fGSAox3IQ\nWelbWJvA7J9VVUXuExFJEJGk7OfrADdjTN2KK9EhOdtxUizOfKwYY9yAZcDnIvJNPk1KdMzYO/xs\nF0QbY6qTdUH0yjxtVgIjwfbtkUsiEm3nOhxJkfvEGONvsi+QNMZ0IesSpJiKL9WhONtxUizOeqxk\n/50/BQ6JyHsFNCvRMWPXYa8U44JoEVlrjOlrjDkGXAbG2LMGR1OcfQLcB4w3xqQDScBQywquIMaY\nhUBvwNcYEwVMJetsuFMeJ38qar/ghMdKth7ACOAXY8ye7J+9CDSG0h0zepGzUsop6TT2SimnpOGn\nlHJKGn5KKaek4aeUckoafkopp6Thp5RyShp+SimnpOGnlHJK/w+n5EfNcREGUQAAAABJRU5ErkJg\ngg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10b11fc90>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize=(5,5))\n",
    "axes = fig.add_subplot(111)\n",
    "x_range = np.arange(0, 2, 0.1)\n",
    "axes.plot(x_range, bayes_factor_sim(x_range, x_pdf), color=\"red\", label=\"Simulated Bayes factor\", linewidth=2)\n",
    "axes.plot(x_range, bayes_factor_exact(x_range), color=\"blue\", label=\"Exact Bayes factor\")\n",
    "axes.legend(loc=2)\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}