{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Estimate probabilities with the Probability Estimator"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `ProbabilityEstimator` is an objective class in Halerium, with which one can estimate the probability of certain events."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ingredients"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To estimate probabilities you only need a graph and data. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Imports"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the examples below, we need to import the following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# for handling data:\n",
    "import numpy as np\n",
    "\n",
    "# for plotting:\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# for some Halerium specific functions:\n",
    "import halerium.core as hal\n",
    "\n",
    "# for building graphs:\n",
    "from halerium.core import Graph, Entity, Variable, StaticVariable, show\n",
    "\n",
    "# for estimating probabilities:\n",
    "from halerium import ProbabilityEstimator "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Basic example"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's create a graph first."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "g = Graph(\"g\")\n",
    "with g:\n",
    "    Variable(\"z\", mean=0, variance=1)\n",
    "    Variable(\"x\", mean=z, variance=0.1)\n",
    "    Variable(\"y\", mean=z, variance=0.1)\n",
    "    \n",
    "show(g)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this graph, both x and y depend on the variable z and follow it fairly closely. The means of x and y are set to z, and the variances of x and y are much smaller than that of z."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's look at a certain set of data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = {g.z: [-1, 0, 1., 2., 3.]}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We evaluate the probability that the graph we defined would generate such data points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "probability_estimator = ProbabilityEstimator(graph=g, data=data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can ask the probability estimator for the probability of any scopetor in the graph."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-1.41893853, -0.91893853, -1.41893853, -2.91893853, -5.41893853])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display(probability_estimator(g)) # get the probability of the whole graph"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We get the logarithmic probability densities for each data point individually. The higher this value, the more likely the event (=data point). We see that the data point with g.z=0 has the highest (least negative) value. It is the most probable event."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Detailed example"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's get some data for g.x and g.y instead of g.z."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = {g.x: [1., -1, -1.],\n",
    "        g.y: [1., -1,  1.]}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We created three cases here:\n",
    " - one where both x and y are 1,\n",
    " - one where both x and y are -1, and\n",
    " - one where x is -1 and y is 1.\n",
    " \n",
    "Now we ask the probabilty estimator for the likelihood of these events."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ -1.6222337 ,  -1.62543181, -11.00601139])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "probability_estimator = ProbabilityEstimator(g, data)\n",
    "probability_estimator(g)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that the examples where both x and y are +1 or both are -1 have the same likelihood. But it is very unlikely that one of them is -1 and the other +1.\n",
    "\n",
    "This reflects that both x and y follow z closely. So they are expexted to be close to each other."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Probabilities of individual variables\n",
    "So far we have always asked for the probabilities for the whole graph. We can also ask for probabilities of some sub-structure of the graph, e.g. a specific variable. We just specify it in the call.\n",
    "\n",
    "For example, we get the probabilities for g.x by:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-1.47250095, -1.49544156, -1.41192295])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "probability_estimator(g.x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "or for g.y by:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-1.47250095, -1.49544156, -1.45960488])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "probability_estimator(g.y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that the probabilites of the three cases are much closer now. This is because neither g.y = +1 nor g.y = -1 is particularly unlikely. It is only the combination of g.x = -1 and g.y = +1 that was unlikely.\n",
    "\n",
    "Since g.x and g.y are not causally connected in the sense that one leads to the other you get their independent probabilities if you ask for them seperately."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Missing data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If some of your data are missing, the probabilities are calculated by marginalizing over the missing values.\n",
    "You can indicate missing value by `np.nan`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'g': array([-11.01642622,  -1.45289939]),\n",
       " 'g/z': array([0., 0.]),\n",
       " 'g/x': array([-1.4624225,  0.       ]),\n",
       " 'g/y': array([-1.41741181, -1.45289939])}"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data = {g.x: [-1., np.nan],\n",
    "        g.y: [ 1.,     1.]}\n",
    "probability_estimator = ProbabilityEstimator(g, data)\n",
    "probability_estimator() # the call without any arguments returns a dict of all possible answers."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that the missing value in g.x yields a logarithmic probability of 0 and the total probability of the graph is equal to the of of g.y."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Background: Logarithmic probability densities"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the class outputs logarithmic probability densities. If we exponentiate these, we get the density."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEGCAYAAAB/+QKOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAxvklEQVR4nO3dd3xc5ZXw8d/RqFrFsoqbLCFblm1scEPu2GAIxFTTkkAKoew6bCA9m/Am75vsbrLZZDekERLWYUlIYR0SIDSHZoqxZYPlgnuR5CYXVXdZbea8f8wVjOWRNLJ1dUfS+X4+85mZe+9z54ytmTNPuc8jqooxxhjTVozXARhjjIlOliCMMcaEZQnCGGNMWJYgjDHGhGUJwhhjTFixXgfQnbKysjQ/P9/rMIwxptdYu3Ztjapmh9vXpxJEfn4+JSUlXodhjDG9hojsbW+fNTEZY4wJyxKEMcaYsCxBGGOMCcsShDHGmLAsQRhjjAnL1QQhIgtEZIeIlIrIgx0cN01E/CJyW1fLGmOMcYdrCUJEfMAjwDXAeOAOERnfznE/Al7palljjDHucfM6iOlAqaqWA4jIEmAhsLXNcV8AngamnUNZY6JWIKCs2VPH+xVHafYrBdkpXDYmm6R4n9ehGRMRNxNEDrA/5HkFMCP0ABHJAW4GruDMBNFp2ZBzLAIWAeTl5Z130MZ0h+LSGr77/BZ2VZ08Y/vApDi+dvUYPj3jAmJixKPojImMmwki3F9/29WJfgZ8U1X9ImccHknZ4EbVxcBigKKiIlv9yHhKVXnkzVIeem0n+ZnJPPSxSVx54WASYn2s33eER94q5TvPbWHFrhp+dvtkBsT3qckMTB/j5l9nBZAb8nwEcLDNMUXAEic5ZAHXikhLhGWNiTo/eW0nD79Rys1TcvjBzRef0Zw0e3QWswoy+e3KPXz/pa384+9LePyuaSTEWpOTiU5ujmJaAxSKyEgRiQduB54PPUBVR6pqvqrmA38FPq+qf4ukrDHR5i8l+3n4jVJun5bLQx+bFLavQUS459KR/Ndtk1hZWsuDT2/Clv010cq1GoSqtojIAwRHJ/mAx1V1i4jc5+x/tKtl3YrVmPO1s/IE3352M5eOzuL7N13Uaf/CrZeMoOLIaX76+k6mj8zgjunWf2aij/SlXy9FRUVqs7mantbsD3Dzr1Zy6GgDr3xlHlkpCRGV8weUu377HiV7jvDqV+aRmzHA5UiNOZuIrFXVonD77EpqY87T4uXlbD5wnH+/+aKIkwOAL0b44a0TEYHvPLfZmppM1LEEYcx5qDrewCNvlnL1+CEsuGhYl8vnpCfx1avG8OaOal7ZUulChMacO0sQxpyHh17dSbM/wLeuvfCcz3HX7HxGD07hP1/ZTos/0I3RGXN+LEEYc45Kq07y1Nr93Dkrn/ys5HM+T6wvhq9fPZby6lM8s+5AN0ZozPmxBGHMOfr1W2UkxMbw+csLzvtcH50whEm56fzs9Z00tVgtwkQHSxDGnIP9dfX8bcMB7pieR2YXOqbbIyJ8+cpCDh5r4MWNdk2oiQ6WIIw5B/+9vIwYgUXzRnXbOS8fm82YISksXl5uI5pMVLAEYUwXHa1v4i8lFdw6dQTDBiZ123lFhH+cO4rth0+wfFdNt53XmHNlCcKYLnqqZD+NLQHumpPf7edeODmHIWkJPPZOebef25iusgRhTBcEAsofV+9jen4G44amdfv542Nj+NSMC3hnVw17a091+/mN6QpLEMZ0wds7q9lXV89nZl3g2mt8vCgXX4zwv+/t7/xgY1xkCcKYLvjj6r1kpybw0QlDXXuNoQMTuWLcYP66dr8NeTWesgRhTISqTzTy1s5qbp06gvhYdz86n5yRR83JJl7batNvGO9YgjAmQs9tOIA/oNw6Ncf115pXmE1OehJPlVgzk/GOJQhjIvTMugNMHDGQwiGprr+WL0ZYOHk4K0prqDnZ6PrrGROOqwlCRBaIyA4RKRWRB8PsXygiG0Vkg4iUiMilIfv2iMim1n1uxmlMZ7YdOs7WQ8e5deqIHnvNm6bk4A8oL75vV1Ybb7iWIETEBzwCXAOMB+4QkfFtDlsGTFLVycA9wGNt9s9X1cntLWZhTE95Zl0FcT7hhknDe+w1xwxJZfywNJ7dYAnCeMPNGsR0oFRVy1W1CVgCLAw9QFVP6odzCiQDNr+AiTqBgPLchoNcPnYwGcnxPfraN0/J4f39R9ldY9dEmJ7nZoLIAUJ72CqcbWcQkZtFZDvwEsFaRCsFXhWRtSKyqL0XEZFFTvNUSXV1dTeFbsyHSvYeoepEY4/WHlrdMGk4IvC39TYNuOl5biaIcKu2n1VDUNVnVXUccBPwvZBdc1R1KsEmqvtFZF64F1HVxapapKpF2dnZ3RC2MWdauukQCbExXDlucI+/9tCBicwcmWkzvBpPuJkgKoDckOcjgHb/ylV1OVAgIlnO84POfRXwLMEmK2N6VCCg/H3zIS4fm01yQqwnMVx78VDKqk+xq/KEJ69v+i83E8QaoFBERopIPHA78HzoASIyWkTEeTwViAdqRSRZRFKd7cnA1cBmF2M1Jqx1+45QebyRay/u+nrT3eVq56rtlzcf9iwG0z+5liBUtQV4AHgF2AY8papbROQ+EbnPOexWYLOIbCA44ukTTqf1EGCFiLwPvAe8pKovuxWrMe1Zuukw8bExXOFB81KrIWmJXHLBIP5uCcL0MFfrzKq6FFjaZtujIY9/BPwoTLlyYJKbsRnTmdbmpXmF2aQmxnkay4IJQ/n3pdvYV1tPXuYAT2Mx/YddSW1MOzZUHOXQsQaum+jexHyRWnCR08y05ZDHkZj+xBKEMe14bWslsTHCFeOGeB0KuRkDmDA8zfohTI+yBGFMO97YVsW0/AwGJnnbvNTqmouGsm7fUSqPN3gdiuknLEEYE8b+unp2VJ7gygu965xu66rxwWamN7ZXeRyJ6S8sQRgTxrJtwXUYPnKh981LrcYMSSEnPckShOkxliCMCWPZ9ioKspPJz0r2OpQPiAhXjBvMytIaGpr9Xodj+gFLEMa0caKhmdXltVwZRbWHVleMG0x9k593d9d5HYrpByxBGNPGO7tqaParJ3MvdWZWQSaJcTG8ac1MpgdYgjCmjde3VTIwKY5LLhjkdShnSYzzMacgi2XbK/lwpnxj3GEJwpgQ/oDy1o5q5o/NJtYXnR+P+eMGs7/uNGXVJ70OxfRx0fkJMMYjmw4co+5UE/OjsHmpVeu8UDaaybjNEoQxIZbvrEYE5hZG79oiw9OTGDc0lWXbLEEYd1mCMCbE8p3VXJwzsMeXFu2q+eMGs3bvEU40NHsdiunDLEEY4zje0Mz6/UeZF8W1h1bzCrNpCSiry224q3GPJQhjHMWltfgDytzCLK9D6dTUC9IZEO/jnV22DrtxjyUIYxzLd1WTkhDL1Cgc3tpWQqyPmaMyeWdXjdehmD7M1QQhIgtEZIeIlIrIg2H2LxSRjSKyQURKROTSSMsa051UleU7q5lVkElclA5vbWtuYRa7a06xv67e61BMH+XaJ0FEfASXEb0GGA/cISLj2xy2DJikqpOBe4DHulDWmG6zu+YUFUdOM29M9Pc/tGodaWW1COMWN38qTQdKVbVcVZuAJcDC0ANU9aR+eDloMqCRljWmO7V+yV7WCzqoWxVkJzN8YKL1QxjXuJkgcoD9Ic8rnG1nEJGbRWQ78BLBWkTEZZ3yi5zmqZLqavugmHOzfGc1+ZkDetV6zyLC3MJsVpbW0OIPeB2O6YPcTBASZttZk8eo6rOqOg64CfheV8o65RerapGqFmVn955ffyZ6NLUEWFVe26ual1rNHZPF8YYWNh445nUopg9yM0FUALkhz0cAB9s7WFWXAwUiktXVssacj5K9ddQ3+aP66un2zCnIQgTe2Wn9EKb7uZkg1gCFIjJSROKB24HnQw8QkdEiIs7jqUA8UBtJWWO6S3FpLbExwqyCTK9D6bJByfFMzBlo/RDGFbFunVhVW0TkAeAVwAc8rqpbROQ+Z/+jwK3AnSLSDJwGPuF0Woct61aspn9bWVbDpNx0UhJc+zi4am5hNr9+u4zjDc2kJcZ5HY7pQzqtQYjIj0VkwrmcXFWXquoYVS1Q1X93tj3qJAdU9UeqOkFVJ6vqLFVd0VFZY7rbiYZmNlYcY3YvrD20urQwC39AWV1W63Uopo+JpIlpO7BYRN4VkftEZKDbQRnTU9bsqcMf0F7ZvNRqSl46iXExFFuCMN2s0wShqo+p6hzgTiAf2CgiT4rIfLeDM8ZtxaW1JMTGMDUv+qfXaE9CrI9p+RmssgRhullEndTOlc3jnFsN8D7wVRFZ4mJsxrhuZVktRfmDSIzzeR3KeZlVkMmOyhNUn2j0OhTTh0TSB/ETgs1M1wI/UNVLnL6DG4ApbgdojFvqTjWx7dBxZhdE/+ytnZnjvIdV5VaLMN0nkhrEZoLzJX1OVd9rs2+6CzEZ0yNWO1+mvbn/odWE4WmkJsayqsyuhzDdJ5IE8SlVPWO6SBFZBqCqdvmm6bVWltaQkhDLxJzeP+4i1hfDjJGZ1lFtulW7CUJEEkUkA8gSkUEikuHc8oHhPRahMS5ZVVbL9JEZxPaS6b07M7sgk7219VQcsem/Tffo6JPxOWAtwY7pdc7jtcBzBKfiNqbXOnTsNOU1p3r19Q9tzRkd7IewWoTpLu0mCFX9uaqOBL6uqiNDbpNU9Zc9GKMx3a51SGhf6KBuNWZICpnJ8Tbc1XSbducWEJErVPUN4ICI3NJ2v6o+42pkxrhoZWktgwbEMW5oqtehdBuR4HxSxWU1qCrONGfGnLOOJp+5DHgDuCHMPgUsQZheSVVZVVbDrIJMYmL61pfo7IIsXtx4iPKaUxRkp3gdjunl2k0Qqvpd5/7ungvHGPftra3n4LEG/qkPNS+1mjM62KdSXFpjCcKct0gulPuSiKRJ0GMisk5Eru6J4IxxQ2sn7pw+1EHdKi9jADnpSdZRbbpFJOP77lHV48DVwGDgbuCHrkZljItWltUwNC2RkVnJXofS7Vr7IVaV1xIIhF2E0ZiIRZIgWhtprwV+q6rvE35JUGOiXsCZFnt2QWaf7cSdXZDJ0fpmth0+7nUoppeLJEGsFZFXCSaIV0QkFbAV0k2vtLPqBLWnmvrE9BrtaR26W1xqzUzm/ESSIO4FHgSmOVNuxBNsZuqUiCwQkR0iUioiD4bZ/ykR2ejcikVkUsi+PSKySUQ2iEhJhO/HmA61fmnOHt33OqhbDR2YyKjsZFbavEzmPHW6xqKqBkSkEhgvIhGvyehMEf4IcBVQAawRkedVdWvIYbuBy1T1iIhcAywGZoTsn6+q9lduuk1xWQ35mcGO3L5sdkEmz647QLM/QFwfmUrE9LxOv/BF5EfAJ4CtgN/ZrMDyTopOB0pVtdw5zxJgoXOe4ElUi0OOXw2MiDhyY7qoxR/g3fI6rp/U96cSm12QxR9X72NjxTEuuaD3LoZkvBVJjeAmYKyqdnUlkhxgf8jzCs6sHbR1L/D3kOcKvCoiCvy3qi4OV0hEFgGLAPLy8roYoulPNh88zonGlj41/1J7Zo4KvsdVZTWWIMw5i6TuWQ7EncO5ww0RCTvuzlm+9F7gmyGb56jqVOAa4H4RmReurKouVtUiVS3Kzs4+hzBNf1HstMn35Q7qVhnJ8Vw4LM2uhzDnJZIaRD2wwVkD4oNahKp+sZNyFUBuyPMRwMG2B4nIROAx4BpV/eCvWVUPOvdVIvIswSarzpq1jGnXqrJaxg1NJSslwetQesTsgkz+sHovDc3+Xr+kqvFGJDWI54HvAcV8OOX32gjKrQEKRWSkiMQDtzvn+oCI5BGc0+kzqrozZHuyM5wWEUkmeJHe5ghe05iwGlv8rNlT90HTS38wuyCTppYA6/Yd8ToU00tFMorpCRFJAvJUdUekJ1bVFhF5AHgF8AGPq+oWEbnP2f8o8B0gE/iVc9FSi6oWAUOAZ51tscCTqvpy196aMR/asO8oDc2BftH/0Gr6yAx8McKqsto+Na256TmRjGK6AfgxwesfRorIZODfVPXGzsqq6lJgaZttj4Y8/gfgH8KUKwcmtd1uzLkqLqslRmBGP6pBpCbGcXHOQIrLavma18GYXimSJqZ/Idj+fxRAVTcAI12LyBgXrCqv5aKcgQxMOpfxFr3X7IJM3t9/lJONLV6HYnqhSBJEi6oea7PNZgEzvcbpJj/r9x1hVj+qPbSaXZBFS0BZs6fO61BMLxRJgtgsIp8EfCJSKCIPE+ywNqZXKNlbR7Nf+8Xw1rYuuWAQ8b4YW4bUnJNIEsQXgAkEh7j+L3Ac+LKLMRnTrYrLaomNEablZ3gdSo9LivcxJS/9g2tAjOmKThOEqtar6rdVdZpzQdq3VbWhJ4IzpjusKqtlUm46yQkRTyXWp8wuyGLLweMcrW/yOhTTy7T7iRGRF+igryGSUUzGeO14QzMbK45y//zRXofimdmjM/np67C6vI4FFw31OhzTi3RUg/gx8BDBGVdPA79xbiexi9ZML7Fmdx0B7R/Ta7Rn0oh0kuJ8rLJmJtNF7dYgVPVtABH5nqqGzoP0gojYlBemV1hVVkt8bAxT8/rvhHXxsTFMG5lh8zKZLoukkzpbREa1PhGRkYDNimd6heKyWi7JG9Tv5yKaU5DJrqqTVJ2w7kMTuUgSxFeAt0TkLRF5C3gT+JKrURnTDY6camLroeP9anqN9rROtWHDXU1XRDIX08siUgiMczZtP4e1IYzpcavLg1+G/bn/odX44WmkJcayqqyWhZNzvA7H9BIRjftzEsL7LsdiTLdaVV7LgHgfE0ekex2K53wxwsxRmdYPYbrEFqs1fVZxWS3T8jOIj7U/cwjOy7Svrp79dfVeh2J6CfvkmD6p6ngDpVUnrf8hxOzRTj9EudUiTGQ6TRAi8rSIXCcilkxMr7HK+h/OUjg4hayUeOuoNhGL5Ev/18AngV0i8kMRGddZAWO8tqqsltTEWCYMH+h1KFFDRJhVkEVxWQ2qNiGz6VwkczG9rqqfAqYCe4DXRKRYRO4WkQ4n1xeRBSKyQ0RKReTBMPs/JSIbnVuxiEyKtKwxHSkuq2XmqEx8MeJ1KFFldkEmlccbKa855XUopheIqNlIRDKBuwiu/rYe+DnBhPFaB2V8wCPANcB44A4RGd/msN3AZao6keC614u7UNaYsCqO1LOvrr5frv/QmdY+GRvNZCIRSR/EM8A7wADgBlW9UVX/rKpfAFI6KDodKFXVclVtApYAC0MPUNViVW1dUX01MCLSssa0p7WNffZoSxBt5WUMICc9yeZlMhGJ5DqIx5y1pT8gIgmq2qiqRR2UywH2hzyvAGZ0cPy9wN+7WlZEFgGLAPLy8jo4vekvistqyUyOZ8zgVK9DiTrBfohMlm2rJBBQYqwJznQgkiam74fZtiqCcuH+8sL2jInIfIIJ4ptdLauqi511Koqys22KqP5OVVlRWsPs0Vn25deO2QWZHKlvZvvhE16HYqJcR+tBDCX4Sz5JRKbw4Zd2GsHmps5UALkhz0cAB8O8zkTgMeAaVa3tSllj2tpReYLqE43Mdcb8m7PN+qAfoobxw9M8jsZEs46amD5KsGN6BPCTkO0ngG9FcO41QKEz++sB4HaCw2U/ICJ5wDPAZ1R1Z1fKGhPOil3BtvVLCy1BtGfYwCRGZSVTXFbLP8wd1XkB0291tB7EE8ATInKrqj7d1ROraouIPAC8AviAx1V1i4jc5+x/FPgOkAn8SkQAWpzmorBluxqD6X/e2VXDqOxkhqcneR1KVJtVkMnf1h+g2R8gzmfXwJrwOmpi+rSq/hHIF5Gvtt2vqj8JU6ztMUuBpW22PRry+B8IDp2NqKwxHWls8fPe7jo+XjSi84P7udkFWfzp3X1sOnCsXy+mZDrW0U+HZOc+BUgNczMmqqzbe5TTzX7mWP9Dp2aOygBsfQjTsY6amP7buf/XngvHmHO3orQ6OK21zb/UqcyUBC4clsaKXTXcP3+01+GYKNVRE9MvOiqoql/s/nCMOXcrdtUwOTedtMQOZ4AxjnmFWTy+cjenGltITohoaRjTz3TUxLS2k5sxUeNofRMbDxzjUmteitjcwmya/frBynvGtNXZKCZjeoXislpUYa4Nb41YUf4gEuNieGdXDVdeOMTrcEwU6qiJ6Weq+mUReYEwVzGr6o2uRmZMF6worSElIZZJueleh9JrJMb5mDkqk+U7q70OxUSpjhoe/+Dc/7gnAjHmfKzYVcPMUZk2pr+L5hZm870dW9lfV09uRiQTJJj+pN1Pk6qude7fJjj30hGgDljlbDMmKuyrDU7vbc1LXXfZmOC/2fJdVoswZ4tkuu/rgDLgF8AvgVIRucbtwIyJ1DulwS83u/6h6wqyUxg+MJF3dtr03+ZskYxtewiYr6qlACJSALzEh1NzG+Op5TurGT4wkYLs5M4PNmcQEeaNyealTYdo8QeItSY6EyKSv4aq1uTgKAeqXIrHmC5pagmwYlcNl48bjDOfl+miuYXZnGhoYcP+o16HYqJMR6OYbnEebhGRpcBTBEczfYzgbKvGeK5kbx2nmvxcPsbWAjlXl47OIkaCNbGi/AyvwzFRpKMaxA3OLRGoBC4DLgeqAZvdy0SFt3ZUE+cT6384DwMHxDEpN523d1k/hDlTRxfK3d2TgRhzLt7aUcX0kRk2VcR5mleYzS/e2MXR+ibSB8R7HY6JEpGMYkoUkftF5Fci8njrrSeCM6YjB46eZmflSeaPHex1KL3evDHZqMJyq0WYEJF0Uv8BGEpwhbm3Ca4wF9FitiKyQER2iEipiDwYZv84EVklIo0i8vU2+/aIyCYR2SAiJZG8nulf3toRHCtx+Vjrfzhfk3PTGTQgjje32/gT86FI6uWjVfVjIrJQVZ8QkScJrvTWIRHxAY8AVxFcY3qNiDyvqltDDqsDvgjc1M5p5quq/aQxYb21o5qc9CQKslO8DqXX88UIl48dzFs7qvAHFF+MjQgzkdUgmp37oyJyETAQyI+g3HSgVFXLVbUJWAIsDD1AVatUdU3IaxgTkcYWPytLa5g/LtuGt3aTK8YN5kh9Mxv2H/E6FBMlIkkQi0VkEPD/gOeBrcCPIiiXA+wPeV7hbIuUAq+KyFoRWdSFcqYfKNlzhPomP5ePsf6H7jJvTDa+GGHZNmtmMkGdJghVfUxVj6jq26o6SlUHt64214lwP+vOmhW2A3NUdSpwDXC/iMwL+yIii0SkRERKqqttPpn+4s3tVcT7Ypg92laP6y4Dk+KYlj+IN6wfwjgiGcWUKSIPi8g659f8z0Qkkk9lBZAb8nwEcDDSwFT1oHNfBTxLsMkq3HGLVbVIVYuys62zsr94c0cVM0ZlMCDehrd2pyvHDWH74RNUHKn3OhQTBSJpYlpCcGqNW4HbgBrgzxGUWwMUishIEYkHbifYRNUpEUkWkdTWx8DVwOZIypq+b3fNKcqqT3HFOGte6m5XXBj8N7XRTAYiG8WUoarfC3n+fRG5qbNCqtoiIg8QHPHkAx5X1S0icp+z/1ERGQqUAGlAQES+DIwHsoBnnc7HWOBJVX058rdl+rLXth4G4KrxtgpadxuVlUx+5gCWba/iM7PyvQ7HeCySBPGmiNxOcC4mCNYiXork5Kq6FFjaZtujIY8PE2x6aus4MCmS1zD9z2tbK7lwWBojBtkCN91NRLhi3BD++O5e6ptarAmvn2u3iUlETojIceBzwJNAk3NbAnylZ8Iz5ky1JxtZu/eI1R5cdOWFg2lqCVBcWut1KMZjHa0ol6qqac59jKrGOrcYVU3rySCNabVsexUBhastQbhmWn4GKQmxLNte6XUoxmMRrQ4iIjeKyI+d2/VuB2VMe17bWsnwgYlMGG6/UdwSHxvDZWOzeW1rJf5AV0amm74mkmGuPwS+RPACua3Al5xtxvSo001+3tlVzVXjh9jV0y675qKh1JxsomRPndehGA9FUoO4FrhKVR9X1ceBBc42Y3rUitIaGpoDXDV+qNeh9HmXjx1MfGwML2857HUoxkORLkCbHvJ4oAtxGNOp17YeJjUxlhmjbNUzt6UkxDKvMJtXNh9G1ZqZ+qtIEsQPgPUi8jsReQJY62wzpse0+AMs21bF/LGDifNF+rvGnI8FFw3l4LEGNlYc8zoU45EOBzmLSAwQAGYC0wjOr/RN5/oFY3rMe7vrqD3VxDUXWfNST/nIhYOJjRH+vvkwk3LTvQ7HeKDDn2KqGgAeUNVDqvq8qj5nycF44cVNhxgQ7+NyWz2ux6QPiGdWQSYvbz5kzUz9VCR19ddE5OsikisiGa031yMzxtHiD/Dy5sNceeEQkuJ9XofTryy4aCh7auvZURnRIpKmj4kkQdwD3A8sJ9j/sJbg/EnG9IjV5XXUnWriuouHeR1KvxMcUgwvb7aGg/4okvUgRoa5jeqJ4IwBeGnTQZLjfbb2tAcGpyYyLT+DFzdaM1N/FMmFcoki8lUReUZEnhaRL4tIYk8EZ0xzSPNSYpw1L3nhxknDKa06ybZD1szU30TSxPR7YALwMPBLgtNx/8HNoIxptaqsliP1zVw30ZqXvHLtxcOIjRGee/+A16GYHhZJghirqveq6pvObREwxu3AjAF4aeMhkuN9XDbGmpe8kpEcz7wx2byw4SABm5upX4kkQawXkZmtT0RkBrDSvZCMCWpo9rN08yE+OmGoNS95bOHk4Rw81kDJ3iNeh2J6UCQJYgZQLCJ7RGQPsAq4TEQ2icjGjgqKyAIR2SEipSLyYJj940RklYg0isjXu1LW9H2vb6vkREMLt0wNt6aU6UkfuXAISXE+nttgzUz9SSTLRS04lxOLiA94BLgKqADWiMjzqro15LA64IvATedQ1vRxz647wJC0BGYVZHodSr+XnBDLVeOH8NKmQ3z3hgnEx9p0J/1BJMNc93Z066DodKBUVctVtXUluoVtzl2lqmuA5q6WNX1b7clG3t5ZzU1TcvDF2NTe0WDh5OEcrW/mnV3VXodieoibPwNygP0hzyucbd1aVkQWiUiJiJRUV9sfbl/xwvsHaQkot0yx5qVoMbcwm4zkeP66tsLrUEwPcTNBhPvZF+kQiIjLqupiVS1S1aLsbBvp0lc8s/4A44elMXZoqtehGEd8bAw3T8nh9W2V1J5s9Doc0wPcTBAVQG7I8xHAwR4oa3q5nZUn2FhxjFumRlrhND3lE9NyafYrz663zur+wM0EsQYoFJGRIhIP3A483wNlTS/3v+/tI84n3DzFEkS0GTMklcm56fx5zX6beqMfcC1BqGoL8ADwCrANeEpVt4jIfSJyH4CIDBWRCuCrwP8VkQoRSWuvrFuxmujR0OznmXUH+OiEoWSmJHgdjgnjE9Ny2VV1kvX7j3odinFZJMNcz5mqLgWWttn2aMjjwwSbjyIqa/q+lzcf5tjpZj45Pc/rUEw7rp84jH97YStPrdnP1LxBXodjXGSDmU1UefK9feRnDmDmKLv2IVqlJsZx3cRhvPD+QU41tngdjnGRJQgTNUqrTvLe7jo+MS2PGLv2IardMT2XU01+/mZXVvdpliBM1PjTu3uJjRFuu8SufYh2U/MGMWF4Gk8U77HO6j7MEoSJCicamvlLSQXXTRxGdqp1Tkc7EeGu2fnsrDzJqrJar8MxLrEEYaLCX0oqONnYwt1zRnodionQDZOGk5Ecz++K93gdinGJJQjjOX9A+V3xHi65YBCTc9O9DsdEKDHOx+3Tcnl9WyX76+q9Dse4wBKE8dyybZXsq6vnHqs99DqfnnkBIsIfVnc0b6fprSxBGM/9z4rd5KQn8dEJQ7wOxXTR8PQkrr14GE++u49j9W0nZTa9nSUI46l1+47w7u467pqdT6zP/hx7o3+6rICTjS38ftUer0Mx3cw+kcZTj7xRyqABcXxyhl053VuNH57G/LHZ/LZ4D6eb/F6HY7qRJQjjmc0HjrFsexX3XjqS5ARXZ30xLvv8/NHUnWpiyZp9XodiupElCOOZR94sJTUxljtn53sdijlP0/IzmJY/iN8sL6epJeB1OKabWIIwnthZeYK/bz7M3bPzSUuM8zoc0w2+cEUhB481WC2iD7EEYTzx0Ks7SEmItQvj+pC5hVnMGJnBL5aVUt9kk/j1BZYgTI9bu/cIr2yp5HPzRjEoOd7rcEw3ERG+sWAcNScb+e3KPV6HY7qBJQjTo1SVH/59G9mpCdw712oPfc0lFwziIxcO4dG3yjhyqsnrcMx5cjVBiMgCEdkhIqUi8mCY/SIiv3D2bxSRqSH79ojIJhHZICIlbsZpes7r26pYs+cIX/5IIQPibeRSX/TPHx3LyaYWHn6j1OtQzHlyLUGIiA94BLgGGA/cISLj2xx2DVDo3BYBv26zf76qTlbVIrfiND2nodnP91/aSkF2Mh8vyvU6HOOSsUNTuX1aHk+s2sPOyhNeh2POg5s1iOlAqaqWq2oTsARY2OaYhcDvNWg1kC4iw1yMyXjoN8vL2Vtbz7/eeBFxdtV0n/bPHx1LamIs33lus60X0Yu5+SnNAfaHPK9wtkV6jAKvishaEVnU3ouIyCIRKRGRkurq6m4I27hhf109v3yzlOsuHsalhVleh2NclpEcz9evHsvq8jpe3HjI63DMOXIzQYRbM7LtT4mOjpmjqlMJNkPdLyLzwr2Iqi5W1SJVLcrOzj73aI1rVJXvPr8FX4zwf6+/0OtwTA+5Y3oeF+Wk8f2XtnLstE3k1xu5mSAqgNCG5hHAwUiPUdXW+yrgWYJNVqYXenrdAd7YXsXXrh7LsIFJXodjeogvRvjBzRdTc7KJ77241etwzDlwM0GsAQpFZKSIxAO3A8+3OeZ54E5nNNNM4JiqHhKRZBFJBRCRZOBqYLOLsRqXHD7WwL++sIXp+RncbVNq9DsTR6TzT5cV8Ne1FbyxvdLrcEwXuZYgVLUFeAB4BdgGPKWqW0TkPhG5zzlsKVAOlAK/AT7vbB8CrBCR94H3gJdU9WW3YjXuCASUbzy9kWZ/gP+8bSIxMeFaFE1f94UrRzN2SCoPPr2J2pONXodjukD60giDoqIiLSmxSyaixSNvlvJfr+zgezddxGdmXuB1OMZDWw4e4+ZfFTNzVCa/u2ua/ViIIiKytr1LCWysoXHFe7vreOjVHVw/cRiftrUe+r0JwwfynevHs3xnNb9+u8zrcEyELEGYbnfo2GkeeHIdF2Qm8x+3XIyI/Vo08KkZedwwaTgPvbqD5TttSHpvYAnCdKtTjS3c+7sS6pv8/PrTU0m1qbyNQ0T4j1suZsyQVO7/0zq7yroXsARhuk2LP8AX/3c9OypP8MtPTmHc0DSvQzJRJiUhlsfvmkZivI+7f7uGqhMNXodkOmAJwnQLf0D52l/eZ9n2Kv7lxglcPnaw1yGZKDU8PYn/+WwRdaea+PRj71Jns75GLUsQ5rwFAsqDT2/kuQ0H+caCsTZiyXRq4oh0/uezReytrefTj73L0XpLEtHIEoQ5L40tfr6wZD1/WVvBl64s5POXj/Y6JNNLzB6dxeI7iyitOsnH/3sVB4+e9jok04YlCHPOjjc0c/dv1/DSxkN869pxfOWqMV6HZHqZy8Zk87u7p3HoaAO3/KqY7YePex2SCWEJwpyTbYeOc+PDK3hvdx0/+fgkFs0r8Dok00vNHp3Fnz83i4Aqt/16FS/Z7K9RwxKE6RJV5amS/dz8q5WcbvazZNFMbpk6wuuwTC83fngaf7t/DoVDUrj/yXV897nNNLb4vQ6r37MEYSJ28Ohp7v7dGr7x141Mzk3nxS/MpSg/w+uwTB8xPD2JPy+axT1zRvLEqr1c+/N3WLOnzuuw+jWbi8l0qqHZz29X7uGRN0vxB5RvLhjLnbPybT4d45q3d1bzrWc2ceDoaW6flstXrxrD4LREr8Pqkzqai8kShGlXQ7OfZ9cf4OFluzh4rIErxw3muzdMIC9zgNehmX7gVGMLP31tJ0+s2oMvRrhnzkjuvXQkmSkJXofWp1iCMF1SdbyBJWv28/tVe6g52cSkEQN58JoLmVWQ6XVoph/aW3uKh17dyfPvHyQ+NoZbpuRw15x8u1K/m1iCMJ06Vt/MGzsqeXb9QVbsqiagMH9sNv84bxSzRmXahHvGc6VVJ3h85R6eXltBY0uA8cPSuGnKcK69eBgjBlmt9lxZgjBnqW9qYWPFMUr21PH2zmrW7j1CQCEnPYmbp+Rw89QcCrJTvA7TmLPUnWriuQ0H+Nv6A7xfcQyA0YNTuGxMNrMLMpmcm27NUF3gWYIQkQXAzwEf8Jiq/rDNfnH2XwvUA3ep6rpIyoZjCeJsjS1+Dh5toLTqJGXVJymtOsm2Q8fZfvgE/kDw//7inIHMH5vNZWMHMyU33TqfTa+xu+YUy7ZV8vbOat7dXUdTSwCA3IwkJo5IZ3R2CgWDUxiVlczIrGSSE2I9jjj6eJIgRMQH7ASuAioIrlF9h6puDTnmWuALBBPEDODnqjojkrLh9KUEEQgoTf4Azf4ATS0Bmpz7Zn+AxpYAp5v8HG9o5kRDC8dPN3O8oYUTDS0cO93E4WMNHD7eSOXxhrMmQhucmsCYIalMyUtnat4gJuemMyg53qN3aUz3Od3kZ2PFUTbsD942HzxGxZHThH7FpSbGMiQtkaFpiQxOSyAzOZ60xDjSkuJITYz94HFSnI+EuBjifTEkxMWQEOsjITaGhNgYYn196+qAjhKEm+l0OlCqquVOEEuAhUDol/xC4PcazFKrRSRdRIYB+RGU7TbXP/wOp5v8KICCErwgzHmKKigavA/5Y1PVs/aHlqG1TLjjncecUSa4rdkfoNnf9cQd74shLSmOIWkJDB+YyJS8dIakJjI8PZGCwSkUZKcwMMnWZzB9U1K8jxmjMpkx6sPBFA3NfvbUnqK8+hR7ak9ReayByuONHD7eQFnZSY7WN3O6uWsX5ImAT4SYGAneC8TECDEi+Jz7GOHDxzEgyBnlP3h81rkl/L52yrQenzEgnqfum9Wl9xEJNxNEDrA/5HkFwVpCZ8fkRFgWABFZBCwCyMs7t6UtCwenBqumEvzHFxHn3nmNkG3BY8R5bT44TpDg/Qf/e3LW/jPLhD4/c198bAxxvhjinV8srY/jfTHEOfdJ8T7SEmPP+OWTEBtjncnGhEiM8zFuaFqHI56aWgKcaAjWwo+fDtbKG5r9NLYEaGxx7p3nrbV5f0Dxa/AHoD+g+AOKanBbQIMtAP6A8zjkV2Voi03bn4Bn/Pg8Y3s7ZdrUjNzgZoII903V9t+kvWMiKRvcqLoYWAzBJqauBNjqp5+YfC7FjDF9QHxsDJkpCdaxHYabCaICyA15PgI4GOEx8RGUNcYY4yI3e1vWAIUiMlJE4oHbgefbHPM8cKcEzQSOqeqhCMsaY4xxkWs1CFVtEZEHgFcIDlV9XFW3iMh9zv5HgaUERzCVEhzmendHZd2K1RhjzNnsQjljjOnHOhrm2rcG9BpjjOk2liCMMcaEZQnCGGNMWJYgjDHGhNWnOqlFpBrY63UcXZQF1HgdRA+z99w/2HvuHS5Q1exwO/pUguiNRKSkvREEfZW95/7B3nPvZ01MxhhjwrIEYYwxJixLEN5b7HUAHrD33D/Ye+7lrA/CGGNMWFaDMMYYE5YlCGOMMWFZgogiIvJ1EVERyfI6FreJyH+JyHYR2Sgiz4pIutcxuUVEFojIDhEpFZEHvY7HbSKSKyJvisg2EdkiIl/yOqaeICI+EVkvIi96HUt3sQQRJUQkF7gK2Od1LD3kNeAiVZ0I7AT+j8fxuEJEfMAjwDXAeOAOERnvbVSuawG+pqoXAjOB+/vBewb4ErDN6yC6kyWI6PFT4Bu0s7RqX6Oqr6pqi/N0NcFVA/ui6UCpqparahOwBFjocUyuUtVDqrrOeXyC4JdmjrdRuUtERgDXAY95HUt3sgQRBUTkRuCAqr7vdSweuQf4u9dBuCQH2B/yvII+/mUZSkTygSnAux6H4rafEfyBF/A4jm7l5prUJoSIvA4MDbPr28C3gKt7NiL3dfSeVfU555hvE2yS+FNPxtaDJMy2flFLFJEU4Gngy6p63Ot43CIi1wNVqrpWRC73OJxuZQmih6jqR8JtF5GLgZHA+yICwaaWdSIyXVUP92CI3a6999xKRD4LXA9cqX33gpwKIDfk+QjgoEex9BgRiSOYHP6kqs94HY/L5gA3isi1QCKQJiJ/VNVPexzXebML5aKMiOwBilS1t80I2SUisgD4CXCZqlZ7HY9bRCSWYCf8lcABYA3wyb68xroEf+k8AdSp6pc9DqdHOTWIr6vq9R6H0i2sD8J45ZdAKvCaiGwQkUe9DsgNTkf8A8ArBDtrn+rLycExB/gMcIXzf7vB+XVtehmrQRhjjAnLahDGGGPCsgRhjDEmLEsQxhhjwrIEYYwxJixLEMYYY8KyBGHMORKRky6fP1tE3nVmCJ3r5msZE45dSW1M9LoS2K6qn/U6ENM/WQ3CGEBEfiQinw95/i8i8jURSRGRZSKyTkQ2ichZM7GKyOWhawCIyC9F5C7n8SUi8raIrBWRV0RkWJjyFzivsdG5zxORycB/Atc6F5olhRxfFHIB2iYRsYuZjCssQRgTtAT4RMjzjwN/ARqAm1V1KjAfeMiZSqJTznxEDwO3qeolwOPAv4c59JfA7521Mf4E/EJVNwDfAf6sqpNV9XTrwapa4mybDLwM/LhL79SYCFkTkzGAqq4XkcEiMhzIBo6o6j7nS/4HIjKP4FTOOcAQIJKJFMcCFxGcTgTABxwKc9ws4Bbn8R8I1hw6JSIfB6bSB2cCNtHBEoQxH/orcBvBKcqXONs+RTBhXKKqzc5kioltyrVwZm28db8AW1R1Vhfj6LTJSEQmAP8KzFNVfxfPb0xErInJmA8tAW4nmCT+6mwbSHCu/2YRmQ9cEKbcXmC8iCSIyECCncsAO4BsEZkFwSYn54u9rWLndSGYkFZ0FKTzGkuAO/vyTLjGe1aDMMahqltEJJXg6n6tTUF/Al4QkRJgA7A9TLn9IvIUsBHYBax3tjeJyG3AL5wv9ViCK4+1nc31i8DjIvLPQDVwdyeh3kQwUf2mtTvE6Y8wplvZbK7GGGPCsiYmY4wxYVmCMMYYE5YlCGOMMWFZgjDGGBOWJQhjjDFhWYIwxhgTliUIY4wxYf1/uy+h9sPKXRAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "delta_z = 0.01\n",
    "z_vals = np.arange(-5, 5, delta_z)\n",
    "log_probs = ProbabilityEstimator(g, {g.z: z_vals})(g)\n",
    "probs = np.exp(log_probs)\n",
    "plt.plot(z_vals, np.exp(log_probs));\n",
    "plt.ylabel(\"probability density\");\n",
    "plt.xlabel(\"value of z\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we integrate that, we get 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9999994265729888"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.sum(probs*delta_z)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Background: skewed probability densities"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Graphs are typically constructed out of normal distributions, but the chain of conditional probabilities can produce very different probabilities if non-linear functions are connecting the variables.\n",
    "\n",
    "Let's create an example graph."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEJCAYAAACZjSCSAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAA+uUlEQVR4nO3deXzU9Z348dd7ct93AkkIBAggh6AgqPWud1fpYeu1tbW11rb0bnftdrvtb+3ttrtttUup62pPdVvbosXiUcUTBeW+QwgkBMh935n374/vJISQY4B858i8n4/HkDm+M/OeTJj3fK73R1QVY4wxkcsT7ACMMcYElyUCY4yJcJYIjDEmwlkiMMaYCGeJwBhjIpwlAmOMiXCuJgIRuVZE9ohIqYjcO8ztaSLylIhsEZEdInKnm/EYY4w5mbi1jkBEooC9wFVAJbABuFVVdw465l+ANFX9ZxHJAfYAk1S125WgjDHGnCTaxcdeCpSqahmAiDwGLAd2DjpGgRQRESAZqAd6R3vQ7OxsnTZtmisBG2PMRPX222/XqmrOcLe5mQgKgIpBlyuBZUOOeQBYDVQBKcDNquod7UGnTZvGxo0bxzNOY4yZ8ETk4Ei3uTlGIMNcN7Qf6hpgM5APLAIeEJHUkx5I5G4R2SgiG2tqasY7TmOMiWhuJoJKYMqgy4U43/wHuxN4Uh2lwAFgztAHUtVVqrpEVZfk5AzbsjHGGHOa3EwEG4ASESkWkVjgFpxuoMEOAe8GEJE8YDZQ5mJMxhhjhnBtjEBVe0VkBbAWiAIeVtUdInKP7/aVwH3AIyKyDacr6Z9VtdatmIwxxpzMzcFiVHUNsGbIdSsHna8CrnYzBmOMMaOzlcXGGBPhLBEYY0yEs0RgjAk7T22pora1K9hhTBiWCIwxYaW2tYvP/n4Tj75eHuxQJgxLBMaYsHKovh2AXUdaAOjt81L8tb/y0Cs28/x0WSIwxoSVCl8i2HOsGYC3yutRhZ88vw+A2x9az6/eKA9WeGHJEsEE09DWzffW7KKjuy/YoRjjikN1TiKoqO+gtauXZ3ccA2B6bjLNnT28VlrHm2X1wQwx7FgimGB+s/4gv3i5jIdfO8CB2jaWfPu5gf84xkwE/V1DANsPN7F2x1EAGtu7Ka1uBeBoc2dQYgtXlggmmNZup4r3q/tq2X64idrWbvYcawlyVMaMn4qGdialxgPwpcc3c6SpkzmTUjjS2Mk+39/60SZLBKfCEsEEc6CmDYA3D9SxvaoJgIZ22+fHTBwV9R2cPz2TGxfmU9fWzeeumMnty4ro7vPyxv46AKpbOvF63dl0ayJytcSECbwDtW1kJMbQ0N7Di7urAafJbMxE0N3rpaqpg6KsQr501SxUFRHhhV3OOMG6vU6Z+p4+pb69m+zkuGCGGzasRTCB9HmVg3XtXD47F4C9x5z+0ob2nmCGZcy4qWntQhUmpzldQ87mhpCfngA4f+uJsVGAdQ+dCksEE8jhhg66+7wsLc4kLvr4W2stAjNR1LY4q4mHftPvTwQAty4tAuCYDRj7zRLBBLGxvJ5L7n8RgOk5yRRlJg7c1tBmLQIzMdS19SeC2BOuT0uIIT8tnivPyuOui4sBmzl0KmyMYIJ4aU8NUR7hrouKOaconaLMRPZV93cNWYvATAy1Lc7f8nB9/+v+6XKiPUKfV/EIHLOuIb9ZIpggtlQ2Mjsvha9dfxYAUwa1CBptjMBMANsqm6hq6gAga0iLACAmyungiI4SclPiqWzsCGh84cwSwQSgqmytbOL6BZMGruvvGkqKjbIWgQl7pdWt3PDAq6TER5MYG0Vi7OgfXdNzkijzTaU2Y3N1jEBErhWRPSJSKiL3DnP7V0Vks++0XUT6RCTTzZgmooN17TR19HB2YfrAddNzkgCYl59GY3sPqjan2oSv53Y600NbOnv9mhI6IyeZ/TWt9nfvJ9cSgYhEAQ8C1wFzgVtFZO7gY1T1flVdpKqLgK8B61TVioScom2HnYVjCwrSBq67pCSHX398KVeclUt3n5d2qz1kwthzO48OnB+uW2ioGTlJtHT2UttqrWF/uNkiWAqUqmqZqnYDjwHLRzn+VuD3LsYzYe071oJHoCQveeA6j0e4uCSHjMQYwAaMTfhqbO9mU0Uj0R5nzYBfLYJc5//C/ppWV2ObKNxMBAVAxaDLlb7rTiIiicC1wB9HuP1uEdkoIhtramrGPdBwt7+mjaLMROKio066LT3R+fZkA8YmXJXXtaMKF8zIAk6eOjqcGTmWCE6Fm4lAhrlupA67G4DXRuoWUtVVqrpEVZfk5OSMW4ATRWl168Af/lA5Kc63p8EVG40JJ0d8s38uneX83/enRTApNZ7E2KiBaqRmdG4mgkpgyqDLhUDVCMfegnULnZY+r3Kgto2ZucMnggUFaWQkxvDsjqPD3m5MqDvsSwSXz8klJkpOmBo9Eo9HKEhP4EijrSXwh5uJYANQIiLFIhKL82G/euhBIpIGXAr8xcVYJqyK+na6+7wjtghiojxcNTeP53dV85fNh+np8wY4QmPOTFVjJ0mxUUzPTuL5L13K+84Ztof5JLmpcRxrsUTgD9cSgar2AiuAtcAu4AlV3SEi94jIPYMOfR/wrKrapN/T0N8HOiM3acRj3nN2Pq1dvXz+sc0875uGZ0y4ONLUweT0BESEqVlJAwvHxpKbEk91c5fL0U0Mrq4jUNU1qjpLVWeo6nd8161U1ZWDjnlEVW9xM46J7KBv97FpWSMngktKsnn0Y0sBbJMaE3aqGjtOKCrnr9zUOGpaumwtgR+s6FyYq2rsID7GQ2bSyDMpRIRLZ+WcUH/ImHBxuLGTgvT4U75fbko83X1emjpsxtxYLBGEuaom59tSf1320czMTWa/JQITRrp6+6ht7WJy2mm0CHwz5o5Z99CYLBGEOefbkn//SUpykymrbaPXBoxNmOjv+vT3b3yw/kRQbQPGY7JEEOaqGjvI9/Pb0ozcZLp7vXzrqR20dvW6HJkxZ+7VfbUALC0+9RJkeb4N7m3AeGyWCMJYV28fNS1dfg+kzZ2cCsBv1h/iqS0jLekwJnS8sq+G4uwkv9YODJWb2t8isEQwFksEYax/T9Z8PwfS5hek8btPLANg95Fm1+IyZjx09faxvqyeS0qyT+v+ibHRJMdF25aVfrBEEMb6V1yeSv/phTOyObcond1HbRqpCW0Hatvo6Olj8bTTr0w/KS2eKtugZkyWCMJYZYPzBz75FAfSZk9KZffRFptfbUJana+EdP+g7+kozk7iQK2tVR2LJYIwtrmikeS46BM2qvfHnEkpNHX0cP/aPTZobEJWXZuTCLJGWSMzluk5SRysa6fPa196RmOJIIxtOFDP4qkZRHnGXkMw2JxJKQD8/KX9/GnTYTdCM+aM1bc6g7yjLZYcy/TsJLr7vBxusO6h0VgiCFP1bd3sq249rWl15xRlcOe7pgFQ2WDlqU1oqm/rRuT4nhqnY7qvGGNZrS2kHI0lgjC1odzZuuF0EkFstIdv3jCPaVmJ9k3JhKy6tm4yEmNPucU7WHG2U4PLNrIfnSWCMFXuGwDr7+Y5HYUZiQMDzsaEAlXlj29XcqSpg7rW7jPqFgJnfCElPtoGjMcQHewAzOlpaO8hJkpIjjv9t7AwI4Hnd1WPY1TGnJnNFY18+f+2AJCWEMPsM/iiA07Bxek5ydY1NAZrEYSpxvZu0hNj/So2N5KC9ARqW7vo7Okbx8iMOX21vimjAE0dPWc0Y6jf9OwkDljX0KgsEYSphvZuMhJjzugxCjOd9QeHbcGNCRENvimjsdHOR9OZdg2Bkwiqmjpp77ap0iOxRBCmGtp6yDiD2RTgjBEANk5gQkb/2oFLSpyN6sejRVCc4wwYl9faDLmRuJoIRORaEdkjIqUicu8Ix1wmIptFZIeIrHMznonEaRGc2X+S/oVotkeBCRUN7d3ERXuYO9kZG+juO/OFYNOzbQrpWFxLBCISBTwIXAfMBW4VkblDjkkHfg7cqKrzgA+6Fc9E09DeQ0bSmXUN5aXGU5CeMDAV1Zhgq29zZgot8dUXOpPyEv2mZTtfeGycYGRuzhpaCpSqahmAiDwGLAd2DjrmNuBJVT0EoKo2hcUPqjowWHymlhVnsm5vDap6RgPPxoyHBt/agUtm5fCHey7gnKKMM37MxNho8tPibZvWUbjZNVQAVAy6XOm7brBZQIaIvCQib4vIHcM9kIjcLSIbRWRjTU2NS+GGj9auXnq9esaDxeAsSKtr62a/fVsyIaC+/fjagSXTMs9oMdlg507N4PX9tXit5tCw3EwEw72DQ9+FaGAx8B7gGuAbIjLrpDuprlLVJaq6JCcnZ/wjDSPHmjtZ9XIZwBmPEQAsm54FwOv7a8/4sYw5Uw1t3WSMwwDxUFeelUdtazdbKhvH/bEnAjcTQSUwZdDlQmDotliVwN9UtU1Va4GXgYUuxhT2/u0v2/nZ30uB8UkE07ISmZ6dxLM7jp3xYxlzpurbuskch5buUJfNziHKI7xgCyiH5WYi2ACUiEixiMQCtwCrhxzzF+BiEYkWkURgGbDLxZjCXnTU8bfsTAeLwVl5ee38SbxRVjcwh9uYYOjp89Lc2Utm0pkPEA+VnhjL7LwUtlc1jftjTwSuJQJV7QVWAGtxPtyfUNUdInKPiNzjO2YX8DdgK/AW8JCqbncrpokgLvr4WzYeg8UA182fTJ9XeWG3fVsywdPY3gNA5jh8wRlOZlIszR09rjx2uHO11pCqrgHWDLlu5ZDL9wP3uxnHRNL/nwXGZ7ENwPyCVLKSYnl9fy03LS4cl8c05lQ1tDstUjfGCABSE6I50mSLJ4djRefCTH1bN4unZvDVa2aPW4tARFhanMmbZbaewATPH9+pBI6Xjh5vaQkxNHVYmYnhWImJMNPY3k1+egLn+2b7jJdlxZkcbuywjWpMUJRWt7Dq5TJuXVrEvPw0V54jNSGG5k7rGhqOJYIw09DeMy7rB4ZaWuwkFmsVmGB4o6weVfj0ZTNce47U+Bi6e71WbXcYlgjCSJ9Xae7sGbcuocFmT0ohNsrD3uqWcX9sY8aytaKRrKRYCjMSXHuOtATnC5QNGJ/MEkEYaeroQRVXWgRRHqEoK9HqsZig2FrZxILCNFfLnKT6EkGTJYKTjJkIROQ/RGReIIIxo6v3zfMfj4Vkw5mWlUR5nSUCE1jt3b3sq27h7MJ0V59noEVg4wQn8adFsBtYJSJv+tYAuDOSY8bU6Jtel+5CiwCgODuRg3XtVo/FBNTWyia8CgsL3f1oSY13Jklai+BkYyYCVX1IVd8F3AFMA7aKyO9E5HK3gzMnavCtIXCtRZCdRFevlyPNna48vjHDeWVfDVEeGSg97ZbjYwQ2hXQov8YIfHsLzPGdaoEtwJd8paVNgPQvuBmP7fuGU5zlzN9++2ADfdYqMAHy0p4aFhdlDHxQu8XGCEbmzxjBj3G6h64Hvquqi1X1B6p6A3CO2wGa49zuGprmW8jzud9v4uFXD7jyHMb06+zp45HXDrCjqplLZ7tfVdhmDY3MnxbBdmChqn5SVd8acttSF2IyIzja1EVibBTJce4sCJ+cFs+Ky2cCsGb7EVeew5h+L+6u5ltP7SQ22sM18/Jcf76YKA+JsVHWIhiGP4ngdlU9YbmpiLwAoKpWyi+Aqho7mJwW79oUOxHhK9fM5otXzmJzRSO1rV2uPI8xAAd8M9Q2/MuVzMxNCchzpsbHWCIYxoiJQETiRSQTyBaRDBHJ9J2mAfkBi9AMONLUQX66ewtu+r37rFxUnb5bY9xSUd9OVlIsaS51dQ7HqTdkiWCo0VoEnwTexhkgfsd3/m2cPQQedD80M9Thxk7y09xPBPPyU0lPjGGjbWpvXHSwrp2irMSAPmdGUszApAtz3IiJQFV/oqrFwFdUtXjQaaGqPhDAGA3Q1dtHbWtXQFoEIsKCgjS2VlrPn3HPwbp2ijIDmwiykuOoa7VEMNRoXUNX+M4eFpH3Dz0FKD7jc7TJmds/OT0+IM93dmEae461WIEu44ruXi9HmjqYGuBEkJMcR42NfZ1ktOknlwJ/B24Y5jYFnnQlIjOsqkYnERQEoEUAsKAgnT6v8uPn9vLxi4rJSw1MAjKR4XBjB16Foix39h4YSVZSLC2dvXT29BEfExXQ5w5lIyYCVf2m7+edp/vgInIt8BMgCmcbyu8Puf0ynDGH/knrT6rqv5/u801kVY3OzkqT0wLXIgBY9XIZda3d/OhDCwPyvCYy9Ne0mhrgMYLsFGc/5Pq27oB0s4YLfxaUfV5EUsXxkIi8IyJX+3G/KJxB5euAucCtIjJ3mENfUdVFvpMlgREc9ZV9mByAwWLneeK5813TANh+2MYKzPh6u7yBKI9w1uTUgD5v//auNjX6RP6sI/iYqjYDVwO5wJ3A90e/C+AsNitV1TJV7QYeA5afdqQR7lhzJ6nx0STEBqY5KyJ884Z5fOHKEvZWt1jFRjOu3iirY0FBmmuLI0fS3yKwAeMT+ZMI+lcvXQ/8r6puGXTdaAqAikGXK33XDXWBiGwRkWdGKnctIneLyEYR2VhTE5lz26ubu8gNQj/9uUUZqMLWCmsVmPHR1tXLlopGLpgxvtut+iM7yUkENmB8In8Swdsi8ixOIlgrIimA14/7DZcshlYyeweYqqoLgZ8Bfx7ugVR1laouUdUlOTnu1yQJRTWtXeT6vs0E0qKidERg40FbU2BOX0+fl91HmwH47ZsH6fUqF4zzvtv+yE5xuoasRXAifxLBx4F7gfN8pSZicbqHxlIJTBl0uRCoGnyAqjaraqvv/BogRkSy/Qk80lS3dAYlEaTGx7CwMJ0XdlUH/LnNxPH4hgqu/a9XuOvRDXx3zW4um50TlBZBYmw0CTFRNkYwhD/7EXiBY8BcEbkEmAek+/HYG4ASESkWkVjgFmD14ANEZJL4CueIyFJfPHWn9AoigKpS3dxFThASAcD1Cyax7XATFfXtYx9szDD2HXP2wn5+VzXXzpvE/3zkPGKigrNTbnZKrCWCIcYcqRGRHwA3AzuB/tVFCrw82v1UtVdEVgBrcaaPPqyqO0TkHt/tK4GbgE+JSC/QAdyiqlYIf4jmzl66er3kpgRnLv918yfz3TW7eWb7Ea5fMJn0xNiAD/KZ8Haovp2UuGg+cuE0Pn35DKI87u1NPJbs5DhqWiwRDObP/+b3ArNV9ZR/c77unjVDrls56PwDgJWrGEP/H21uanBaBFMyE5mRk8SbZfV8d81uUuKi2fiNK4mLtgU5xj9ltW1cMjuHr1wzO9ihkJ+WwM4jzcEOI6T40zYrAwJXHtCcpLrFWUOQkxycRACwsDCdV/bVAtDS1cuDL+4PWiwmvHT19lFR386M7MCuIh5Jfno8hxs7sM6H4/xJBO3AZhH5hYj8tP/kdmDmuGC3CMBZadzd50wWy0mJ46U9NnhsRqaq9Pr+Xg7VteNVmJ6THOSoHPnpCXT3eqlrs5lD/fxJBKuB+4DXOV6K+m03gzIn6k8EOUEaIwA4e0o6ABmJMXxoSSE7qppp77ZNwM3wfv9WBed953k6uvsGumGKQ6ZF4KzO7y/bYvwYI1DVR0UkAShS1T0BiMkMUd3SRWy0h9T44A3Qzp2cSrRHWDw1kyVTM+nz7mdrZRPnB2EuuAl9v3qjnIb2Hv7v7Qq+t2Y3BekJzMoLzC5kYykYSASdnF0Y5GBChD+1hm4ANgN/811eJCKrR72TGVfVzc4aAre2qPRHfEwU337vfFZcMZNzitIBePtgQ9DiMaFrZ1Uzu48600V/9OxeAP706QsDVh5lLAXWIjiJP18xv4VTN+glAFXdLCLFLsZkhgjWquKhbllaNHB+enYSWysbgxeMCVl/236EKI8QG+WhqaOHd83MCkp5lJGkJ8aQEBNliWAQf8YIeofZpN6G2wMomIvJRjIrL4V91a1sqWhkf01rsMMxIWTjwQbmTEphbr5TWXTptNDqPhQR8tPjqWqyRNDPn0SwXURuA6JEpEREfoYzcGwCpLqlK2iLyUZSkpfMwbp2lj/4Gu/+0bpgh2NCRJ9X2VLRyOKpGZw12RkTOK84I8hRnSwvNZ5jzbaorJ8/ieCzOGUluoDfA83AF1yMyQzS1dtHU0dPSHQNDTYzN5k+rzUMzYn2HG2hrbuPc4syuGJOLnMmpXDOlNBLBBmJsbaJ/SD+zBpqB77uO5kAC4U1BMMpyQ2NGSAmdKzeUsXnfr8JcMqXF2UlcsWcvCBHNbz0xBga222PjX4jJgIReYpRxgJU9UZXIjInqB5YQxBaiWB6ThIegf5Gge0Ba/6+6xgAH71wGlMyQ3sbyIzEWBrbu/F6FU8Q6x6FitG6hv4D+BHOfsIdwC99p1Zgu/uhGXAGioGQGyOIj4k6YV74Md9WmiZy7atu5ZJZOXzrxnlBnersj/TEGLzqlEsxoyQCVV2nquuAc1T1ZlV9yne6DbgocCFGtv6dlEJtjADgVx9fysp/PBeAI02WCCJZn1cprW5lVm5olJEYS0ais0FNo40TAP4NFueIyPT+C741BJG5TVgQVDd3IgKZvk23Q0luSjwlvlbB0aZO+rzK7948xF82Hw5yZCbQKurb6er1hszq4bGkJzp1NBtsnADwb0HZF4GXRKTMd3kacLdrEZkTlNW0UZSZSHSQNvEYyyTfQqGqpg4+9Zu3eXan00/c3NnLh8+fGszQTADtq3bWkszMC48WQbqvRWAzhxz+zBr6m4iUAHN8V+0+nb0JzOnZV91CSQg3t5PioomL9vDDvzllqL56zWxe2lPNQ6+UWSKIENXNnfx6/UGAkP5bHSzD1yKwriGHX18zVbVLVbf4Tn4nARG5VkT2iEipiNw7ynHniUifiNzk72NHgp4+Lwdq2wa6X0JVie9b4MUl2Xzq0hlcMSePg3XtNFiZ34jwtSe3sb6sjs9eMZOU+PDYuqR/jKChzbqGwL+uodMiIlHAg8BVOBvZbxCR1aq6c5jjfoCzpaUZpLy2jZ4+DflvWQ/cei51bV0sLEzH4xEW+UpWb6ls5LLZucENzrjmqS1V/GXzYV7YXc0Xr5zF568sCXZIfktNiEHEWgT93KxrvBQoVdUyABF5DFiOs/fxYJ8F/gic52IsYam/3zXUB+CmZScxbVCt+QWFaYjAloomSwQTSHNnD509fXR2e9lf28q3/7qTY81dJMVGcccF4dUNGOUR0hJibLDYx5/N6/8IPAw8o6reU3jsAqBi0OVKYNmQxy4A3gdcgSWCk+w71ooIzAiRnZ38lRwXTUluMn/dVsWNi/JDZkMSc2b+5cltbCxvICcljm2HnTqUD9x2DgsL08kIwVltY8lIjKWxwxIB+DdG8N/AbcA+Efm+iMwZ6w4+w60oGbpS+b+Af1bVvlEfSORuEdkoIhtramr8fPrwd7S5g6ykuJCp434qPnP5TKoaO/nCY5uCHYoZB16v8mppLUebO9l2uInzpmWwfFE+71kwmSmZicEO77SkJcTYOJaPP7OGngeeF5E04FbgORGpwFll/BtVHSmlVgJTBl0uBKqGHLMEeMy3CjEbuF5EelX1z0NiWAWsAliyZEnEVDqrbe0mKwy/aQEsX1TAjqpmHn29nD6vEmXL+MPa7qMtNLb3kBATRXyMh0c/tpTE2ODtmDceJqfFs8e3gU6k8+udFJEs4B+BDwObgN/irC7+CHDZCHfbAJT4FqAdBm7BaVkMUNWBDW5E5BHg6aFJIJLVt3WH5EIyf83ISaKr18uql50lKJ+6bAZ7jrbwxv5a3jUzO+RnQxlHc2fPwCLBJz55Acnx0WGfBMCpl/XczmP09HmJCdF1OoHizxjBkzhrCH4N3KCqR3w3PS4iG0e6n6r2isgKnNlAUcDDqrpDRO7x3b7yjKOf4Orbupnn29wjHE33jW3cv3Y3AMXZiXz1/7bS0tVLXmocb/7LlSPeV1XxKtaSCAErfreJl/fWUJydxILCtGCHM26Ks5Pp9SqH6tvDbhxuvPmT1h9S1TWDrxCRON/agiWj3dF3vzVDrhs2AajqR/2IJaLUtnaFbdcQHB/k7q9Qes9v3qEgPYEbFuXzuzcPjVqx9K/bjvCvf97OK/90edjMTZ+IjjV38sq+Gj64uJAVV8wMdjjjanqOM4nhQE1bxCcCf9pD3x7mujfGOxBzou5eLy2dvWQlh16xOX9lJsUOrOA8a3Iqealx/PauZSyZ6mxUMlqhuk2HGmls72HXEevDDYTOnj66e70D53v6nPNPbalC1enWm5o1sWZ/TffNZiurta1WR9uPYBLOFNAEETmH47OAUoHwnCYQRup9sxnCeYwAnFbB7qMt/OnTFyICcdFRAwmgqrHjhKml68vq+MPblfzwA2dTUd8OwK4jzSwtzgxK7JHktl+upzAjkZ/csog7/uctjrV08r33LWDlujLOLUof6OabSNITY8lMiqWspi3YoQTdaF1D1wAfxZnt8+NB17cA/+JiTAaoa3MqeWQnh3ci+PhFxdS0dp3QBZSf7itU19jBvmMtbDrUyIfOm8KDL5byyr5a7rl0OocGJQLjriNNHbxzqJEdVc28WlrLW+X1eARue+hNEmKi+P4Hzg52iK6Znp1EWa0lghETgao+CjwqIh9Q1T8GMCbD4BZB+HYNAVy3YPJJ101K608EnWwoL+OJjZXMyE3mtdJaADaWN5zQIjDuWrfHWZvT1evlw//zFilx0fzpMxeyuaKJefmpIb+y/UzkpcWzq8r+xkbrGvpHVf0NME1EvjT0dlX98TB3M+OkrnVidA0NJy46iuzkWI40dbDT90H/hcc34VWIj/Hw7M5jtHX3kRATxZ5jLbR19ZIUF23rEVzy/K5jTE6LH+iy+8o1s5mZm8LMCNiXOispltpWK6Y8WtdQf+ftxOscDAN1vhZBuHcNjSQ/PYFD9e3sPeYM1FXUd3Djwnw6evp4zrenwc3nTeHRN8pZ/uBrfOXq2Xzpic2sXvGuiPiActv9a3czLz+N8ro2nt9Vzacvm8Hlc3JRJaLGZLKS4mju7KW710tsdOSuJRita+gXvp//L3DhmH51rV1EeYTUCTp1cnJaPGt3OB/4t5w3Ba8q/758Pk9srBhIBLctK2LhlDS++PgWfvTsHtq7+/jP5/bx4O3O9piqGvJ74wZKb5+Xr/5hK3ddXMy8/NHn+pfXtvHgi/spSE+gprWLa+dN4stXz47I1lZW8vENavJSQ2tf8EAarWvop6PdUVU/N/7hmH51rc6qYs8E/c85Oy9lIBF87KLigX7o25dN5fXSOl4traUoM5H89AQ8soV91U4Bvr9uO8KnDjdR3dLJxx7ZyMtfvZwer5fPP7aJ/3fjfBb7pqZGmvK6dv606TBPb61i33eup6O7j1tWvcEdF0zjA4sLae3q5ddvHCQrOZbdvim5hxs7APj8lSURmQTgeIu7trXLEsEI3g5YFOYk1S2dIblh/Xj55KUz+OnfS4Hj87nBWUm88sOLT2iqzy9IY2tlEx+5YBpPvlPJj5/by9ZKp/rl+rI6NlU0sP1wM19+YjMrrijhxoX5rjTzVRVVQjI5V/k+1Hv6lMvuf5GclDi2VDbxnTW7uPKsPG5a+fpAWXOA86ZlsO1wE3MmpXLW5PBdvX6m+idj1Ed48bmxZg2ZIKlp7SJnAieCpLho3vr6u6mobx92P+bBH+TLijPZWtnEZbNzyEmJ4/61ewZu21zZyJ82HWbRlHT2HG3hK/+3hdhoDxfNzCYjMeakrqMzqSvz0f/dQEN7N6tXXHRa9x9JRX07v3nzIF+9evawv4va1i6yx1hY2J8IACobOiiva2dGThL7a9r44hOb2Vfdynfft2DguotLsimvaxvYczpS9XcN9U/OiFSjdQ39l6p+QUSe4uTy0ajqja5GFuFqWrqYO8G/qeWmxJObMvYH0XvPKWDPsVaWFmdy0cxsspNj2XeslWe2H+UPGyvp7vNy3/L5FOckseTbz/HkO5V86fHN/PKOJVw+5/jGOB3dfZx733PcurSIf7th7inHu26vM83y7YP1LJ56+gOqrV297Kw6vlBu1ctl/Hr9Qf5hQf5JtXzWl9Vx6y/X89u7lnHhjOyTHktVeXlfLftrWvEIfPu9C1g2PZM1W49w/dmT+fendvL33dVEeYTrF0wiPTGWZdOzAMK2fPR4yva1COqsRTCiX/t+/kcgAjHHeb1KbWv3hG4RnIp5+Wn86mNLBy7ffF4RAFVNHazZdpSclDjmF6QiIiwoSOMl37z4N8rqWLe3ho9fVMyUzERKq1vp6Onj4dcOcMWcXL7+52186tIZPPJ6OTctLuSui6cP+/wv7q4mcdCeECvXlfHLO8ZOBKrKqpfLWLPtCBeX5PClq2bR2NHD/Wt389iGCtZ/7d1kJcWyZptTx3HXkeaTEsGzO46hCt9ds4uYKA/ffd8CZuel8OPn9vK7tw5x4Ywsnt7q3D8/LZ7bljm/m8++29k28vZlRazbW8PiqRmkJ07MGWhnIjUhmmiPUBfhU0hH6xp62/dznYjE4lQgVWCPqkZ2+nRZQ3s3fV4lJ4zrDAWCM430KBeXZA90AS0sTGdDeQMAj2+ooKmjh7gYD1+77ixKa47XLfrs79+hob2HX7xcxoHaNlauK+OOC6bx5oE68lLjT1hE9a9/3o6q0yiOjfawvqzOrzUND79Wzvee2c20rEQeeLGUNw/UsaXCGdtQhddKa8lKjhv4Ntq/pqLPq7xWWsuM3GRe3uckte2Hndvue3onK66YyQMvlpKdHDuQBMCZkjvUFXNyuXBGFjefN+Wk2wyICJlJsdY1NNYBIvIeYCWwH6feULGIfFJVn3E7uEhV3eJ8O8nxo9skks3Kc5a4XDorZ+C6hVPSB843+bYhfG7nMWKjPBysayfKI3zi4umsXLcfgAO+8gK1rV38edNh/umPWwH406cv5JyiDNq7ewdm1wC8d1E+T2ysZGtlI15V8lLjKcw43sXS0NZNU0cP07KTWLvjKPPyU/nDPRdy8Q//zobyBhJioujq7SM5LppXS2sRhJT4aKZmJQ4kgme2H2HF747v7PbJS6ZT3dJFemIM//taOV5VROCPn7qQB18sZffRFrZWNlGQcXIiiI7y8LtPnH9Gv+eJbnAyjlT+jJr9CLhcVS9T1UuBy4H/dDesyFbjSwS5qdYiGM275+Tx1Wtmc828SQPXXTgji/kFqbz/3IKB68pq2vjZ30tZvaWKqVmJ3HHBVKI9Qr6v1EVGYgwF6Qn86Lnjg9D9iWJoQbL3nuM87vt+/jof+O83uPT+l/iPtXvo7Omjqb2Hz/zuHW5a+QadPX1sqWhkWXEWCbFRfPnq2cwvSOWlr17G375wCZfNzuGFXdWs3nKYa+dNYmFhOm8dqGfd3hq2VDQCTp2mqVmJ/OP5U/nPmxfxz9fOISU+mvVl9Zw1KZWpWUn88KaF3O7rDhquRWDGlp1sq4v92Y+gWlVLB10uA6pdisdwPBFY19DoEmKj+MzlJ9bIz0qO4+nPXszLe2t48p3DXDd/En/bcZT46Cg6evqYmZNMfnoCaz5/MV09Xm544FUWT82gMCORR14vB+Dikmw2VzTyvWd2sXNQHRqPwHnTjo8N3Pfe+bx1oJ4HXizlqa1V1LV209rVC8Bjbx2iq9fLkmnOuoZblxZx61LnAzsvNZ7liwoGunVuXJRPVWMHv30TPvLwW2Qnx7KwMI1v/MNcvvEPxwe142OieM+CyTy2oeKE1b/nT89CBIonWJnoQMlLjWfvscjZC304o80aer/v7A4RWQM8gTNG8EGcbSiNS2pa+7uGLBGcrnOnZnDRzGw+f2UJX756FqXVbdzzm7eZmet0J83KS8HrVZYVZ3LDwnwyk2J55PVyZuYmc8WcXF7ZV8sv1pWd8JiT0xKIifLwzRvm0t7dx4fPn8r7zyngrQN1VNS3k50cR2y0h8b2br67xtmVbaQFblfNzeNXH1vKa6W1XDA9iz5VspLi+ORv3qa2tZur5k4a9n4fXFLIYxsquGjm8RlEU7OSeGrFRRO6OJybCtITqG7piugyE6O1CG4YdP4YcKnvfA3g1/JNEbkW+AnOVpUPqer3h9y+HLgP8AK9wBdU9VX/Qp+4qpu7SIyNIiku/PeFDZbkuGh+c9eygctFmUm8/5wCrh9UDdXjER7/5AWAsxFQRmIMF83MPmGcoV9eatxAH/yd7xrYapukuGgeuXMpDW3dzC9Mo72rjzsf2cCuI80sLEwbdbXqJbNyuMQ3vhENXDk3j6XTMnmjrI75BcNPHV48NZMXv3IZ07JOnPo5v2DibCEZaAUZCag65bgn2uY7/hpt1tCdZ/LAIhIFPAhcBVQCG0RktaruHHTYC8BqVVURORun1THnTJ53Ijjc2E6B9feOq9hoDz++edGotz/9uYtJT4ghyiPERAmFGYkcqG0j2iN85erZpCUMX/dp8Mrc1PgYvn79Wby2v5ZPXzbjlOO8el4eb5TVsWCUD/bBm/mYM1fo+792uMESwYhEJB74ODAPGPh6o6ofG+OuS4FSVS3zPc5jwHJgIBGo6uA94pIYZuFaJDpU30GRLfYJuMHJd8XlJczKSyYtMYaspDhmT/K/2+WikmwuKjl58Zc/bltWxOS0hFETgRlf/S29ykGzwyKNP30PvwZ24+xY9u/A7cAuP+5XAFQMulwJLBt6kIi8D/gekAu8x4/HndBUlYr6dpZFUCngUPT5K0uC8rxx0VFcO3/48QHjjslpCYg4LYJI5c/IyExV/QbQ5qs/9B5ggR/3G261zXClKv6kqnOA9+KMF5z8QCJ3i8hGEdlYUzOxR/fr25yZJ9YiMCYwYqM95KbEnbBeJNL4kwh6fD8bRWQ+kAZM8+N+lcDg5YyFQNVIB6vqy8AMETmpTa2qq1R1iaouycnJGebeE0f/Xr2WCIwJnIL0BGsRjGGViGQA3wBW4/Tx/8CP+20ASkSk2Fei4hbf/QeIyEzx1QYQkXOBWKDuFOKfcAYSQZYlAmMCpSgzkfK6yN3EfswxAlV9yHd2HTB8Va7h79crIiuAtTjTRx9W1R0ico/v9pXAB4A7RKQH6ABu1v6iLhGqf9P2KRmWCIwJlJK8FP68uYrWrl6SI3Datj+zhrKAbwHvwunjfwW4T1XH/OauqmuANUOuWzno/A/wr3URMQ7Vt5OTEkfCoGqXxhh39S803F/dOuw6konOn66hx3BKSnwAuAmoBR53M6hIdqi+3cYHjAmw/kQweBe3SOJPIshU1ftU9YDv9G0g3eW4IlaFrSEwJuCmZiYSEyXsq24Z++AJyJ9E8KKI3CIiHt/pQ8Bf3Q4sEnX3eqlq6rCdo4wJsOgoD9Ozk9kfoS2C0YrOteCMCQjwJeA3vps8QCvwTdejizCHGztQtamjxgTD9Jwk9hy1FsEJVDVFVVN9Pz2qGu07eVR1Ym+mGyS2hsCY4MlLjR8oAR9p/JonJSI3Apf4Lr6kqk+7F1LkskRgTPDkpMTR0tVLR3dfxM3aG3OMQES+D3weZyHZTuDzvuvMOKuobx9Y7m6MCaz+/3eR2Crwp0VwPbBIVb0AIvIosAm4183AItGhOmfqqGeMTdGNMeOvfyOomtbOiFvZ7+92POmDzlt9XJfYGgJjgifHWgSj+i6wSURexJlBdAnwNVejikD95aeXWvlpY4IiN8XZbqXaEsGJRMSDs43k+cB5OIngn1X1aABiiyiN7T20dPXaGgJjgiQzKRaPWIvgJKrqFZEVqvoEQyqHmvFlM4aMCa4oj5CdHBeRicCfMYLnROQrIjJFRDL7T65HFmEsERgTfDkpcdY1NIL+vYk/M+g65RRKUpux9SeCKZm2ab0xweIkgs5ghxFw/uxHUByIQCLdobp2spNjSYyNvFroxoSKosxENpY34PVqRE3j9mdBWbyIfElEnhSRP4rIF0QkPhDBRZKy2lamZycHOwxjItqCgjRau3opq42s4nP+jBH8CpgH/Ax4AJgL/NrNoCLR/po2ZuQmBTsMYyJa/6Y0WyqaghtIgPmTCGar6sdV9UXf6W5glj8PLiLXisgeESkVkZNWIovI7SKy1Xd6XUQWnuoLmAjq27qpb+tmRo61CIwJphk5ySTGRrG1sjHYoQSUP4lgk4ic339BRJYBr411JxGJAh4ErsNpRdwqInOHHHYAuFRVzwbuA1b5G/hEsr/GaYbOyLVEYEwwRXmE+QVpbKm0FsFQy4DXRaRcRMqBN4BLRWSbiGwd5X5LgVJVLVPVbpwtL5cPPkBVX1fVBt/F9UDhKb+CCaB/M4yZ1iIwJuhm5SVTXtcW7DACyp8pKtee5mMXABWDLlfiJJWRfBx4ZrgbRORu4G6AoqKi0wwndJVWtxIX7aEg3aaOGhNsk9MSaGzviahy1P5MHz14mo893NwrHfZAkctxEsFFI8SwCl+30ZIlS4Z9jHB2oLaN4uykiJquZkyompTqTIo82txJcXZkTODwt/ro6agEpgy6XAhUDT1IRM4GHgKWq2qdi/GErIP17UyNsLK3xoSqyWlOIjjS1BHkSALHzUSwASgRkWIRiQVuYUi9IhEpAp4EPqyqe12MJWR5vcqh+namZkXGNw9jQt0kXyI42hQ5K4xdW8aqqr0isgJYC0QBD6vqDhG5x3f7SuDfgCzg5yIC0KuqS9yKKRQda+mku9drNYaMCRGTBloElgjGhaquAdYMuW7loPN3AXe5GUOoO1jn1BiyriFjQkNibDRpCTER1SJws2vI+OFQfyLItK4hY0LF5LT4iGoRWCIIsvK6NqI9Qn66lW8yJlRMSovnWLMlAhMgB+vbKchIIDrK3gpjQkVBegIH69pQnXCz1Ydlnz5Btr+61WoMGRNi5uan0tzZS2VDZEwhtUQQRL19Xspq2yixGkPGhJQFBWkAbDscGTWHLBEEUUVDB929Xis2Z0yImZWXQrRH2G6JwLit1FdszloExoSW+JgoSvJS2F7VHOxQAsISQRDtq24BYKYlAmNCzvz8VHZaIjBuK61uZVJqPCnxMcEOxRgzxPScZGpbu2jp7Al2KK6zRBBE/VVHjTGhZ5pvtX//6v+JzBJBEB2qa2eaJQJjQlJ/IUhLBMY1zZ091LV1D3zrMMaElv76X5GwW5klgiAZqDFk5aeNCUlJcdHkpMRx0BKBcUv/t4xp2dYiMCZUTctKpNy6hoxb+vsdbR8CY0LX1Kwkymomfs0hSwRBUl7bRl5qHImxrm4JYYw5AwsL06ht7aKifmLXHLJEECR7j7UwPdsWkhkTypZNzwJg/YGJvZ26q4lARK4VkT0iUioi9w5z+xwReUNEukTkK27GEkq6evvYeaSZhVPSgx2KMWYUJbnJZCbF8mZZfbBDcZVr/RIiEgU8CFwFVAIbRGS1qu4cdFg98DngvW7FEYp2HWmhp09ZNCUt2KEYY0YhIiydlsmb1iI4bUuBUlUtU9Vu4DFg+eADVLVaVTcAE38N9yBbKhoBrEVgTBhYMi2DyoYOqifwjmVuJoICoGLQ5UrfdadMRO4WkY0isrGmpmZcggumzRWN5KbEMSnVtqc0JtSdU5QBwDuHGoMbiIvcTAQyzHWnNQdLVVep6hJVXZKTk3OGYQWX16u8WlrLecWZiAz3KzLGhJL5BanERnnYdKgh2KG4xs1EUAlMGXS5EKhy8fnCwrbDTdS0dHHlWbnBDsUY44e46Cjm5qeyyVoEp2UDUCIixSISC9wCrHbx+cLCC7uO4RG4bJYlAmPCxeKpGWypbKSzpy/YobjCtUSgqr3ACmAtsAt4QlV3iMg9InIPgIhMEpFK4EvAv4pIpYikuhVTKFi3r5ZzijLISIoNdijGGD9dMD2Lrl4v70zQ7iFXl7Wq6hpgzZDrVg46fxSnyygidHT3seNwE3dfMj3YoRhjTsGy6ZlEeYTXSmu5cEZ2sMMZd7ayOIA2VzTS61XOm5YZ7FCMMacgJT6GRVPSebV0Yq4nsEQQQBvLndWJ5/qmoxljwsdFM7PZWtlIbWtXsEMZd5YIAuit8npm5SWTlmh7FBsTbq6ZNwlVeH7nsWCHMu4sEQRIU3sP68vquGy2zRYyJhydNTmFosxE/rbjaLBDGXeWCALk+V3H6OlTrl8wOdihGGNOg4hw7fxJvFZaS3PnxKqKY4kgQNZsO0JBegILC63QnDHh6pp5k+jpU17cXR3sUMaVJYIAaO7s4ZV9tVw3f5KVlTAmjJ0zJZ3clDj+tn1idQ9ZIgiAv++qprvPy3XWLWRMWPN4hGvmTeLFPdU0dUyc7iFLBAHw9NYqJqXGc46VnTYm7N183hQ6e7w8vuFQsEMZN5YIXFZa3coLu6t5/7kFeDzWLWRMuJtfkMYF07N45LVyevu8wQ5nXFgicNnPXywlLtrDxy4qDnYoxphx8pELp1LV1Mlr+yfGSmNLBC5aX1bHk5sO89ELi8lOjgt2OMaYcXL5nFzSEmJ48p3KYIcyLiwRuKSzp497/7iVqVmJfP7dJcEOxxgzjuKio7hh4WTW7jhKZUN7sMM5Y5YIXPLj5/ZSXtfO996/gITYqGCHY4wZZ5+8ZAbRHg9ffHxz2O9TYInABU9vrWLVy2XctqxoQpasNcbAlMxEvvO++Wwob+ATv9pId2/4DhxbIhhnWyoa+fITW1g8NYNv3jA32OEYY1y0fFEBP7zpbF7ZV8t9T+9E9bS2ZQ86VxOBiFwrIntEpFRE7h3mdhGRn/pu3yoi57oZj5tUlbU7jvKR/32LnJQ4fvHhxcRFW5eQMRPdh5ZM4RMXF/Pr9QdZ8btNYbnQzLUdykQkCngQuApnI/sNIrJaVXcOOuw6oMR3Wgb8t+9n2FBVtlY2cf/aPbxaWsu8/FR+fvu5NkvImAjytevOIjMpjh89u4dNhxr44lWzeM/Zk0mMdXUTyHHjZpRLgVJVLQMQkceA5cDgRLAc+JU67an1IpIuIpNV9YiLcZ0WVaWtu4/G9m4ON3RQXtfG2wcbeKOsjor6DtITY/jmDXO5fdlUYqOtx82YSOLxCJ+6bAYXzMji63/axlf/sJVv/GU7l5TksHBKOnMmpTA5LYHc1DgyE2NDbnGpm4mgAKgYdLmSk7/tD3dMATDuiWDd3hrue3onXlVUnQ92rzL8ZU687PU6SaDPe2L/X1pCDEuLM7nn0hn8w9n5pCXYhjPGRLJFU9J5asVFvHmgnme2H+HFPdU8O2QjG49AbLSH2CgPsdEeYqKcU39uEBEG0oQc/yEi3HLeFO66ePz3PHczEQyX8oaOpPhzDCJyN3A3QFFR0WkFkxwXzey8FEScX6hHwOP7hZ9weejtvp/JcdGkJkSTGh/D5PQEpmYmUpSZGHKZ3RgTXB6PcMGMLC6YkQU41Yf3HWulurmT6pYu6lq76Or10tXrpafPS7fvp/MF9PgHYP/Asw78g2tdzm4mgkpgyqDLhUDVaRyDqq4CVgEsWbLktIblF0/NYPFU2yvYGBNYqfExIf/Z42Zn9gagRESKRSQWuAVYPeSY1cAdvtlD5wNNoTg+YIwxE5lrLQJV7RWRFcBaIAp4WFV3iMg9vttXAmuA64FSoB240614jDHGDM/VuU2qugbnw37wdSsHnVfgM27GYIwxZnQ2z9EYYyKcJQJjjIlwlgiMMSbCWSIwxpgIZ4nAGGMinIRb2VQRqQEOnubds4HacQwnmOy1hCZ7LaHJXgtMVdWc4W4Iu0RwJkRko6ouCXYc48FeS2iy1xKa7LWMzrqGjDEmwlkiMMaYCBdpiWBVsAMYR/ZaQpO9ltBkr2UUETVGYIwx5mSR1iIwxhgzxIROBCLyQRHZISJeERlxlF1EykVkm4hsFpGNgYzRX6fwWq4VkT0iUioi9wYyRn+JSKaIPCci+3w/hy3WHsrvy1i/Z19p9Z/6bt8qIucGI05/+PFaLhORJt/7sFlE/i0YcY5FRB4WkWoR2T7C7eH0noz1Wsb3PVHVCXsCzgJmAy8BS0Y5rhzIDna8Z/pacMp97wemA7HAFmBusGMfJs4fAvf6zt8L/CCc3hd/fs845dWfwdmF73zgzWDHfQav5TLg6WDH6sdruQQ4F9g+wu1h8Z74+VrG9T2Z0C0CVd2lqnuCHcd48PO1LAVKVbVMVbuBx4Dl7kd3ypYDj/rOPwq8N3ihnBZ/fs/LgV+pYz2QLiKTAx2oH8Llb2ZMqvoyUD/KIeHynvjzWsbVhE4Ep0CBZ0Xkbd/+yOGqAKgYdLnSd12oyVPfTnS+n7kjHBeq74s/v+dweS/8jfMCEdkiIs+IyLzAhDbuwuU98de4vSeubkwTCCLyPDBpmJu+rqp/8fNh3qWqVSKSCzwnIrt9GTmgxuG1yDDXBWVa2Giv5RQeJiTel2H483sOmfdiDP7E+Q5OeYJWEbke+DNQ4nZgLgiX98Qf4/qehH0iUNUrx+Exqnw/q0XkTzjN5YB/4IzDa6kEpgy6XAhUneFjnpbRXouIHBORyap6xNc0rx7hMULifRmGP7/nkHkvxjBmnKraPOj8GhH5uYhkq2q41e4Jl/dkTOP9nkR815CIJIlISv954Gpg2JH6MLABKBGRYhGJBW4BVgc5puGsBj7iO/8R4KTWToi/L/78nlcDd/hmqpwPNPV3h4WYMV+LiEwSEfGdX4rzuVEX8EjPXLi8J2Ma9/ck2KPjLo+8vw/nW0AXcAxY67s+H1jjOz8dZ6bEFmAHTjdM0GM/ndfiu3w9sBdnJkiovpYs4AVgn+9nZri9L8P9noF7gHt85wV40Hf7NkaZtRbskx+vZYXvPdgCrAcuDHbMI7yO3wNHgB7f/5WPh/F7MtZrGdf3xFYWG2NMhIv4riFjjIl0lgiMMSbCWSIwxpgIZ4nAGGMinCUCY4yJcJYIjBmFiLS6/Pg5IvKmiGwSkYvdfC5jRhL2K4uNCXPvBnar6kfGPNIYl1iLwEQMEfmBiHx60OVviciXRSRZRF4QkXd8+x+cVH3TV//96UGXHxCRj/rOLxaRdb7ieGuHq2gpIlN9z7HV97NIRBbhlOS+3ldTPmHQ8Wm+PQJm+y7/XkQ+MZ6/D2P6WSIwkeQx4OZBlz8E/B/QCbxPVc8FLgd+1L98fywiEgP8DLhJVRcDDwPfGebQB3BKIJ8N/Bb4qapuBv4NeFxVF6lqR//BqtqEs3r0ERG5BchQ1V+e0qs1xk/WNWQihqpuEpFcEckHcoAGVT3k+zD/rohcAnhxShPnAUf9eNjZwHyc6qjgbPQyXP2aC4D3+87/GqclMFa8z4nIB3HKIiz0IxZjToslAhNp/gDchFMi+zHfdbfjJIbFqtojIuVA/JD79XJiC7r/dgF2qOoFpxjHmLVdRMSDszNdB5CJU3PGmHFnXUMm0jyGU2HzJpykAJAGVPuSwOXA1GHudxCYKyJxIpKGM8gLsAfIEZELwOkqGmGTkNd9zwtO4nnVj1i/COwCbgUe9rVcjBl31iIwEUVVd/jKWx/W4yWIfws8JSIbgc3A7mHuVyEiTwBbcaqmbvJd3y0iNwE/9SWIaOC/cCpDDvY5nA/zrwI1wJ2jxSkis4C7gKWq2iIiLwP/Cnzz1F+1MaOz6qPGGBPhrGvIGGMinCUCY4yJcJYIjDEmwlkiMMaYCGeJwBhjIpwlAmOMiXCWCIwxJsJZIjDGmAj3/wF3FaKzzX/WAQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "g2 = Graph(\"g2\")\n",
    "with g2:\n",
    "    Variable(\"z\", mean=0., variance=1.)\n",
    "    mu_x = hal.cos(z*3.)\n",
    "    Variable(\"x\", mean=mu_x, variance=0.01)\n",
    "\n",
    "delta_x = 0.01\n",
    "x_vals = np.arange(-1.5, 1.5, delta_x)\n",
    "\n",
    "log_probs = ProbabilityEstimator(g2, {g2.x: x_vals}, n_samples=10000)(g2)\n",
    "\n",
    "plt.plot(x_vals, np.exp(log_probs));\n",
    "plt.xlabel(\"value of x\");\n",
    "plt.ylabel(\"probability density\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This probability density is clearly not a normal distribution, but it still integrates to 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9984480589771683"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.sum(np.exp(log_probs) * delta_x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}