Part_2_Rebuilding_the_Transportation_model_Pyro.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "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.8.5+"
    },
    "colab": {
      "name": "Part 2: Rebuilding the Transportation model Pyro.ipynb",
      "provenance": [],
      "collapsed_sections": []
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "CgpdSs0fu5bx"
      },
      "source": [
        "# Rebuilding the Transporation Model in Pyro\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "zXkGLCrnw8og",
        "outputId": "8e990a06-81bf-482e-e36b-424c5523df27"
      },
      "source": [
        "!pip install torch torchvision torchaudio\n",
        "!pip3 install pyro-ppl"
      ],
      "execution_count": 3,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Requirement already satisfied: torch in /usr/local/lib/python3.7/dist-packages (1.8.1+cu101)\n",
            "Requirement already satisfied: torchvision in /usr/local/lib/python3.7/dist-packages (0.9.1+cu101)\n",
            "Collecting torchaudio\n",
            "\u001b[?25l  Downloading https://files.pythonhosted.org/packages/aa/55/01ad9244bcd595e39cea5ce30726a7fe02fd963d07daeb136bfe7e23f0a5/torchaudio-0.8.1-cp37-cp37m-manylinux1_x86_64.whl (1.9MB)\n",
            "\u001b[K     |████████████████████████████████| 1.9MB 4.8MB/s \n",
            "\u001b[?25hRequirement already satisfied: numpy in /usr/local/lib/python3.7/dist-packages (from torch) (1.19.5)\n",
            "Requirement already satisfied: typing-extensions in /usr/local/lib/python3.7/dist-packages (from torch) (3.7.4.3)\n",
            "Requirement already satisfied: pillow>=4.1.1 in /usr/local/lib/python3.7/dist-packages (from torchvision) (7.1.2)\n",
            "Installing collected packages: torchaudio\n",
            "Successfully installed torchaudio-0.8.1\n",
            "Collecting pyro-ppl\n",
            "\u001b[?25l  Downloading https://files.pythonhosted.org/packages/aa/7a/fbab572fd385154a0c07b0fa138683aa52e14603bb83d37b198e5f9269b1/pyro_ppl-1.6.0-py3-none-any.whl (634kB)\n",
            "\u001b[K     |████████████████████████████████| 634kB 7.0MB/s \n",
            "\u001b[?25hRequirement already satisfied: torch>=1.8.0 in /usr/local/lib/python3.7/dist-packages (from pyro-ppl) (1.8.1+cu101)\n",
            "Collecting pyro-api>=0.1.1\n",
            "  Downloading https://files.pythonhosted.org/packages/fc/81/957ae78e6398460a7230b0eb9b8f1cb954c5e913e868e48d89324c68cec7/pyro_api-0.1.2-py3-none-any.whl\n",
            "Requirement already satisfied: tqdm>=4.36 in /usr/local/lib/python3.7/dist-packages (from pyro-ppl) (4.41.1)\n",
            "Requirement already satisfied: numpy>=1.7 in /usr/local/lib/python3.7/dist-packages (from pyro-ppl) (1.19.5)\n",
            "Requirement already satisfied: opt-einsum>=2.3.2 in /usr/local/lib/python3.7/dist-packages (from pyro-ppl) (3.3.0)\n",
            "Requirement already satisfied: typing-extensions in /usr/local/lib/python3.7/dist-packages (from torch>=1.8.0->pyro-ppl) (3.7.4.3)\n",
            "Installing collected packages: pyro-api, pyro-ppl\n",
            "Successfully installed pyro-api-0.1.2 pyro-ppl-1.6.0\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "sTDSnDGUu5b3"
      },
      "source": [
        "\n",
        "\n",
        "from statistics import mean\n",
        "import torch\n",
        "import numpy as np\n",
        "import pyro\n",
        "import pyro.distributions as dist\n",
        "from pyro.infer import Importance, EmpiricalMarginal\n",
        "import matplotlib.pyplot as plt\n",
        "%matplotlib inline\n",
        "\n",
        "pyro.set_rng_seed(101)"
      ],
      "execution_count": 4,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "1BIiYV7vu5b4"
      },
      "source": [
        "Now that we have validated our model, let's rebuild the model with Pyro.\n",
        "\n",
        "Usually we would take advantage of Pyro's ability to employ arbitrarily complicated parametric assumptions.  But here, we're going to simply replicate the original model to establish proof of concept."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "YNiQg0aNu5b6"
      },
      "source": [
        "- Age (A): Recorded as young (young) for individuals below 30 years, adult (adult) for individuals between 30 and 60 years old, and old (old) for people older than 60.\n",
        "- Sex (S): The biological sex of an individual, recorded as male (M) or female (F).\n",
        "- Education (E): The highest level of education or training completed by the individual, recorded either high school (high) or university degree (uni).\n",
        "- Occupation (O): Employee (emp) or a self employed (self) worker.\n",
        "- Residence (R): The size of the city the individual lives in, recorded as small (small) or big (big).\n",
        "- Travel (T): The means of transport favored by the individual, recorded as car (car), train (train) or other (other)\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "YBPO4mC3u5b7"
      },
      "source": [
        "![dag.png](attachment:dag.png)"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "tbecTBMBu5b_",
        "outputId": "c18d1de6-2002-47a2-be53-a9a85acaae69"
      },
      "source": [
        "A_alias = ['young','adult','old']\n",
        "S_alias = ['M','F']\n",
        "E_alias = ['high','uni']\n",
        "O_alias = ['emp','self']\n",
        "R_alias = ['small','big']\n",
        "T_alias = ['car','train','other']\n",
        "\n",
        "A_prob = torch.tensor([0.3,0.5,0.2])\n",
        "S_prob = torch.tensor([0.6,0.4])\n",
        "E_prob = torch.tensor([[[0.75,0.25], [0.72,0.28], [0.88,0.12]],\n",
        "                     [[0.64,0.36], [0.7,0.3], [0.9,0.1]]])\n",
        "O_prob = torch.tensor([[0.96,0.04], [0.92,0.08]])\n",
        "R_prob = torch.tensor([[0.25,0.75], [0.2,0.8]])\n",
        "T_prob = torch.tensor([[[0.48,0.42,0.1], [0.56,0.36,0.08]],\n",
        "                     [[0.58,0.24,0.18], [0.7,0.21,0.09]]])\n",
        "\n",
        "def model():\n",
        "    A = pyro.sample(\"A\", dist.Categorical(probs=A_prob))\n",
        "    S = pyro.sample(\"S\", dist.Categorical(probs=S_prob))\n",
        "    E = pyro.sample(\"E\", dist.Categorical(probs=E_prob[S][A]))\n",
        "    O = pyro.sample(\"O\", dist.Categorical(probs=O_prob[E]))\n",
        "    R = pyro.sample(\"R\", dist.Categorical(probs=R_prob[E]))\n",
        "    T = pyro.sample(\"T\", dist.Categorical(probs=T_prob[R][O]))\n",
        "    return{'A': A,'S': S,'E': E,'O': O,'R': R,'T': T}\n",
        "\n",
        "print(model())"
      ],
      "execution_count": 5,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "{'A': tensor(1), 'S': tensor(1), 'E': tensor(0), 'O': tensor(0), 'R': tensor(1), 'T': tensor(2)}\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "hLI04CTpu5cA"
      },
      "source": [
        "def prob_dict(alias, prob):\n",
        "    dict = {}\n",
        "    for i, _ in enumerate(alias):\n",
        "        dict[alias[i]] = prob[i]\n",
        "    return dict"
      ],
      "execution_count": 6,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "UKuYB1C9u5cC"
      },
      "source": [
        "A_dict = prob_dict(A_alias, A_prob)\n",
        "S_dict = prob_dict(S_alias, S_prob)\n",
        "E_dict = prob_dict(E_alias, E_prob)\n",
        "O_dict = prob_dict(O_alias, O_prob)\n",
        "R_dict = prob_dict(R_alias, R_prob)\n",
        "T_dict = prob_dict(T_alias, T_prob.T)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "EezE_q06u5cC"
      },
      "source": [
        "### Inference\n",
        "\n",
        "There are two broad categories of causal inference questions, one is to predict effects from causes, referred to as forward causal inference, and the other is to predict causes from effects, referred to as reverse causal inference. Forward causal inference is straightforward to do through forward sampling, without the need of any inference algorithm. Forward sampling doesn't work for reverse causal inference, since there is no guarantee the sample generated matches the evidence we observed. An easy solution is rejection sampling, which rejects all the samples that do not match the evidence. But it can be very inefficient when the probability of the observed evidence is very small. An improvement is likelihood weighting, a special case of importance sampling, which samples unobserved variables like in forward sampling, but sets the observed variables to the same value as evidence, and weights each sample by its likelihood. An caveat of importance sampling is its accuracy depends on how close the proposal distribution is to the target distribution. In the case of Bayesian network, the proposal distribution is often the CPTs, and the target distribution is the posterior distribution we want to estimate. (For details of the theory of different sampling methods, see \"Probabilistic Graphical Models\" Chapter 12)\n",
        "\n",
        "Pyro has implemented many inference algorithms including importance sampling, and we encourage you to experiment with them. Detailed understanding about how these methods work in Pyro is not a requirement. For the purpose of the homework, it's sufficient to have a high level understanding of the general idea and be able to use at least one of the inference algorithms in Pyro in your model. \n",
        "\n",
        "Most of the inference algorithms implemented in Pyro are approximate inference algorithms either through sampling or variational inference. Exact inference using variable elimination is NP hard. It becomes intractable when the model is large. However, most models in homeworks are fairly small. Although inference can be done manually without much efforts in small networks, inference using Pyro is required in this question."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "6FoXl0_cu5cC"
      },
      "source": [
        "### Forward inference\n",
        "\n",
        "Implement forward inference using `pyro.condition`. Assuming you observed a person with a university degree. What is your prediction of this person's means of travel? Provide a histogram of the marginal distribution on the variable \"T\". (6 points)\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "TVMbr-mJu5cC"
      },
      "source": [
        "conditioned_model = pyro.condition(model, data={'E':torch.tensor(1),'O':torch.tensor(1), 'R':torch.tensor(0)})"
      ],
      "execution_count": 7,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 312
        },
        "id": "-K3nWHz8u5cC",
        "outputId": "752bdd1d-7fad-4369-8cff-143e5eca7cba"
      },
      "source": [
        "T_posterior = pyro.infer.Importance(conditioned_model, num_samples=5000).run()\n",
        "T_marginal = EmpiricalMarginal(T_posterior,\"T\")\n",
        "T_samples = [T_marginal().item() for _ in range(1000)]\n",
        "T_unique, T_counts = np.unique(T_samples, return_counts=True)\n",
        "\n",
        "plt.bar(T_unique, T_counts, align='center', alpha=0.5)\n",
        "plt.xticks(T_unique, T_alias)\n",
        "plt.ylabel('count')\n",
        "plt.xlabel('T')\n",
        "plt.title('P(T | E = uni, O = emp, R = big) - Importance Sampling')"
      ],
      "execution_count": 9,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "Text(0.5, 1.0, 'P(T | E = uni, O = emp, R = big) - Importance Sampling')"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 9
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEWCAYAAACJ0YulAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAb2UlEQVR4nO3deZRcZZ3/8fcHwh4ghMQASSAMBhHGQTGDgHhkQESYYYIOqygh4kSP/BDFBdzG+Pu5oMOuIzM5oARhgIAiEUGJYVGQxSCRfclgIAmBrOwgE/z+/niefnJTqU5XN327Op3P65w6fbe69b11b9XnPnepVkRgZmYGsF67CzAzs/7DoWBmZoVDwczMCoeCmZkVDgUzMyscCmZmVgz4UJD0CUnn9PC5kyVN7uWS+oykL0u6oN11GEgKSW/uZNyxkm5ocT6HSrqid6uz/ip/B12Su7eX9KKk9et8zbU+FCTNlfRKfrOekXSRpMF53IbAV4F/l/SePM2Lkl7KH9IXK4/te/DakedVnc8Xe3sZeyoivh0RH291ekm7Spou6TlJL0i6SdI+ddbYbg3rcIGks+r+0DWKiEsj4v0tTvsLYDdJf9fT18ufkW/29Pm9qfql159I2lfS7/NnYZmk2yT9fTtriognI2JwRLxe5+us9aGQHRoRg4E9gHGkIAAYDzwcEQsi4nf5DR0M7JbHD+kYFhFP9vC1d6/MY3BEfO+NLUp7SNoJuA24D9gR2A64GrhB0t7trK0P7J63i/cCRwEfa3M9XbkMmNTuIt4oSYPaXUMzkrYArgW+DwwFRgLfAP7Szrr6ykAJBQAiYgFwPfC3edDBwC3tq+iNyy2h91X6q83JMXlPd4KkJyUtkfSVZtO2YDJwe0R8JSKWRcQLEXEe8BPgu724SB21fUzSQ5KWS/q1pB0q40LSpyQ9llss/0/STnnP7XlJ03IrEEn7SZqfD5Utye/XsT2pKSLmkILx7b2zlKs5RNLjuc5/l7QegKTjJd3aMZGk90t6JO+l/lDSLZKqLb6bgX/sjYIq29BESfPy+vikpL+XdK+kZyX9oDL98Xmv+Qe5voclHVAZv11ubS6TNEfSv1bGTZZ0laRLJD0PfBL4MnBUbqn9KU83MW8bL+T36xOVeXSs789JWiRpoaSJlfGbSDpT0hO5vlslbZLH7ZW3oWcl/UnSfp28LTsDRMRlEfF6RLwSETdExL15PjtJulHS0rwuL5U0pFLDXElfyO/fS5IulDRC0vV5mX4jaauG93+SpKfy8ny+i3U1KPffnD8bt+X53iBpWGX64/L7sFTS19TwXdKZARUKkkYDhwD35EFvAx5pX0UrSTotb4xNH29w9vsCbwEOAP5N0lt7MI8DgSubDJ8GvLvjg1WldIyz02WS9OFmLyRpPOnL4EPAcOB3pL3fqoOAdwJ7AV8EpgAfAUaTQv+YyrTbAMNIe3QTgCmS3tLqglfq2gV4DzBnDdNcu4blvbaLl/ggqSW7B6kVu1qLJH+orwK+BGxN2n4bD+E9BIxR2qPtLe8CxpJaSucAXwHeR2pVHynpvQ3T/g/pPf868DNJQ/O4y4H5pJbm4cC3Je1fee540vINAS4Evg1ckVvZu+dpFgH/BGwBTATOlrRHZR7bAFuS1vcJwH90fMkCZ5C2m31Ie/lfBP4qaSTwS+CbefjngZ9KGt7kvXgUeF3SVEkHV+bdQcB38jK+lbRNTm6Y5l9In6mdgUNJO6tfJm3v6wGfbpj+H0jv//uBU1v58s4+THqP3gRsmJcLSbsCPwSOBbZl5fvVtYhYqx/AXOBF4FngifxGbJLHPQZ8oMlzxgABDOpi3pOByWsYH8Dz+bU7HgfVsHzva6jpkoblGFUZfxdwdOO0LbzOik7eq13ya4zsxWW6Hjih0r8e8DKwQ+V9fXdl/N3AqZX+M4Fzcvd+ufbNKuOnAV9rsZaOdfhS7r4M2KiG7TSq7y/wKWBm7j4euDV3H0dqsXVMJ2Ae8PHKsA3y/LbvYS0XAd9s2IZGVsYvBY6q9P8U+Eyl1qcANWxzHyV9Ob4ObF4Z9x3gosr2+NuGWrrcRoGfAydX1vcrVD67pBDZK29Hr5AOBzbO41TgJw3Dfg1M6OQ135rfp/l5+5oOjOhk2sOAeyr9c4FjG96/8yv9JwE/b3j/d6mM/x5wYeP7Q8P3FqnF+NWGbepXufvfgMsq4zYFXqPyXdLZY6C0FA6LiCERsUNEfCoiXsnDlwOb1/zae+TX7nj8uubXa+bpSvfLwOAezGMJaY+i0bbAX0nvZW/ZATi30kpaRvryq+7JPFPpfqVJf3UZl0fES5X+J0h7ca3aI8/vKNJe8GbdeG53zKt0d1bjdtXpIn2i5zdM07FNr9bCVLqSqeOih+u7UVt33u8Fua4OHcuyHbAsIl5oGFddr9X3oKm8d35HPgT1LKn1P6wyydKIWFHp79jmhwEbk1oxjXYAjmhone9L822eiHgoIo6PiFGklul2pBYU+VDQ5UoXJjwPXNJQH3Tv/YTWto1mOvvsN25HL5PCvksDJRQ6cy/5+GC7KR3zfrGzxxqe+hIp5TtsU1OJvwGOaDL8SNKe68uNI7TyErnOHp0d258HfKIhTDeJiN/3sPatJFW/yLcn7c22LJJpwO2kvaym8nHhzpa3qy/h0S3UuBAYVXk9VfuztwJzI+L5Jstxaay86OHgLurpqZG5rg4dy/IUMFTS5g3jFlRLbJjXKv2SNiLtWZ9B2jMfAlxH2mnoyhLgVWCnJuPmkVoK1W1us4g4vauZRsTDpFZDx7nKb+e63xYRW5AOa7ZS35q0sm10R+N2tAnpcGSXBnooXEe6oqTtIl0eOrizxxqeOhs4WtIGksaRjtP2SD7RdHwno78B7CPpW5KGStpc0kmkwxmndrJMT65pmSLi0k5e6z+BL0naLde1paRmgdQd35C0oaT3kI5HX5nnfbykud2Yz+nAv0pqGr4RcfAalrerL+EvSNoqn/s6GWh2v8EvgbdJOkzphOKJrL4j8F7SIbh2eRPw6bxNHkEKqesiYh7we+A7kjZWumz2BNKedGeeIZ0f6fgu2hDYCFgMrJB0MOk4e5ci4q/Aj4CzlE54ry9p7xw0lwCHSjooD99Y6aR1Y+AiaRelE9mjcv9o0jmsO/Ikm5MOWT+Xz1V8oZX6uvA1SZvmz8REmm8b3XEVaXn3UbooYzItBtdAD4VfALtI6s6hhO76U8PeYo9ulFuDr5H2fJaTvrj/uyczyRvG1qzcsFcREY+RmtO7k46JLiSdLDsoIm7ryWt2JiKuJl3RdHluft9PulKsp54mvT9PAZcCn8x7d5D2wFquPyLuA35L73zQG11DOj8ym/Tlf2GT119CarF9j9Tc3xWYxaqXQx4D/FcN9bXqTtJJ0SXAt4DDI6Lj0MQxpGPfT5Euaf56RPxmDfPquLhhqaQ/5kNPnyadF1pOOpE6vRu1fZ50WfUfSIclvwuslwOr4wKHxaSWwxdo/h34Aukw4p2SXiJ9Zu4HPpfHf4N0yPE50nr8WTfq68wtpAscZgJnRERLNzN2JiIeIJ27uJz0WX6RdO6ly8tqteqhwYFH0iRg14j4TA+eOxkgIib3cll9TtK+wIkRcUyXE69FlC4rvCQf+202/gbSScqH+rSwXpL3oOeTTlzeJOlQ4KMRcWSb6jmedNJ733a8/kAjaQzwZ2CDhvMkvf06g0nnoMZGxJ/XNG2/vHmkN0XElHbX0B9ExK3ArV1OOMBEi3cK9yeSDiLtjb9C2psVuYUX6Y7mX7SvOltb5B2ImaTt5wxSC2puV8+r9fCRpCFKN6s8rHQzyt75ePUMpRuTZmjlTRySdJ7SDS/3atXrktvl5vww60t7k66gWUK6xv2wyhV1Zq0az8oLAMaSLlXv8tBQrYePJE0FfhcRF+Rj2puSjukti4jTJZ0GbBURp0o6hHQM7BDS8bxzI+JdtRVnZmarqS0UJG1JOqH2N9V0kvQIsF9ELJS0LXBzRLxF0n/l7ssap6ulQDMzW02d5xR2JJ3l/7Gk3UlXXZxMuva444v+aWBE7h7JqjdwzM/DOg2FYcOGxZgxY3q5bDOzge3uu+9eEhHNfuKj1lAYRLps66SIuFPSucBp1QkiIiR1q6mSryaaBLD99tsza9as3qrXzGydIOmJzsbVeaJ5PjA/Iu7M/VeRQuKZfNiI/HdRHr+AVe/qG8Wqd0IC6WqiiBgXEeOGD28adGZm1kO1hUJEPA3M08pfqzwAeJB0I8qEPGwC6YYe8vDj8lVIewHP+XyCmVnfqvs+hZOAS/OVR4+Tbt9eD5gm6QTSDz913IRzHenKozmkH3aauPrszMysTrWGQkTMJv1+fKMDmkwbpN95MTOzNhnov31kZmbd4FAwM7PCoWBmZoVDwczMCoeCmZkVA/6nsztz9oxH213CgPXZA/vFf0A1sx5wS8HMzAqHgpmZFQ4FMzMrHApmZlY4FMzMrHAomJlZ4VAwM7PCoWBmZoVDwczMCoeCmZkVDgUzMyscCmZmVjgUzMyscCiYmVnhUDAzs8KhYGZmhUPBzMwKh4KZmRUOBTMzKxwKZmZWOBTMzKxwKJiZWVFrKEiaK+k+SbMlzcrDhkqaIemx/HerPFySzpM0R9K9kvaoszYzM1tdX7QU/iEi3h4R43L/acDMiBgLzMz9AAcDY/NjEnB+H9RmZmYV7Th8NB6YmrunAodVhl8cyR3AEEnbtqE+M7N1Vt2hEMANku6WNCkPGxERC3P308CI3D0SmFd57vw8bBWSJkmaJWnW4sWL66rbzGydNKjm+e8bEQskvQmYIenh6siICEnRnRlGxBRgCsC4ceO69VwzM1uzWlsKEbEg/10EXA3sCTzTcVgo/12UJ18AjK48fVQeZmZmfaS2UJC0maTNO7qB9wP3A9OBCXmyCcA1uXs6cFy+Cmkv4LnKYSYzM+sDdR4+GgFcLanjdf47In4l6Q/ANEknAE8AR+bprwMOAeYALwMTa6zNzMyaqC0UIuJxYPcmw5cCBzQZHsCJddVjZmZd8x3NZmZWOBTMzKxwKJiZWeFQMDOzwqFgZmaFQ8HMzAqHgpmZFQ4FMzMrHApmZlY4FMzMrHAomJlZ4VAwM7PCoWBmZoVDwczMCoeCmZkVDgUzMyscCmZmVjgUzMyscCiYmVnhUDAzs8KhYGZmhUPBzMwKh4KZmRUOBTMzKxwKZmZWOBTMzKxwKJiZWVF7KEhaX9I9kq7N/TtKulPSHElXSNowD98o98/J48fUXZuZma2qL1oKJwMPVfq/C5wdEW8GlgMn5OEnAMvz8LPzdGZm1odqDQVJo4B/BC7I/QL2B67Kk0wFDsvd43M/efwBeXozM+sjdbcUzgG+CPw1928NPBsRK3L/fGBk7h4JzAPI45/L069C0iRJsyTNWrx4cZ21m5mtc2oLBUn/BCyKiLt7c74RMSUixkXEuOHDh/fmrM3M1nmDapz3u4F/lnQIsDGwBXAuMETSoNwaGAUsyNMvAEYD8yUNArYEltZYn5mZNaitpRARX4qIURExBjgauDEijgVuAg7Pk00Arsnd03M/efyNERF11WdmZqtrx30KpwKnSJpDOmdwYR5+IbB1Hn4KcFobajMzW6fVefioiIibgZtz9+PAnk2meRU4oi/qMTOz5nxHs5mZFQ4FMzMrHApmZlY4FMzMrHAomJlZ4VAwM7PCoWBmZoVDwczMCoeCmZkVDgUzMyscCmZmVjgUzMyscCiYmVnhUDAzs6JPfjrbrDecPePRdpcwYH32wJ3bXYL1E24pmJlZ4VAwM7PCoWBmZoVDwczMipZCQdLMVoaZmdnabY1XH0naGNgUGCZpK0B51BbAyJprMzOzPtbVJamfAD4DbAfczcpQeB74QY11mZlZG6wxFCLiXOBcSSdFxPf7qCYzM2uTlm5ei4jvS9oHGFN9TkRcXFNdZmbWBi2FgqSfADsBs4HX8+AAHApmZgNIqz9zMQ7YNSKizmLMzKy9Wr1P4X5gmzoLMTOz9mu1pTAMeFDSXcBfOgZGxD939oR8OetvgY3y61wVEV+XtCNwObA16Yqmj0bEa5I2Ih2OeiewFDgqIuZ2f5HMzKynWg2FyT2Y91+A/SPiRUkbALdKuh44BTg7Ii6X9J/ACcD5+e/yiHizpKOB7wJH9eB1zcysh1q9+uiW7s44n394MfdukB8B7A98OA+fSgqc84HxrAyfq4AfSJLPY5iZ9Z1Wf+biBUnP58erkl6X9HwLz1tf0mxgETAD+B/g2YhYkSeZz8o7o0cC8wDy+OdIh5ga5zlJ0ixJsxYvXtxK+WZm1qKWQiEiNo+ILSJiC2AT4F+AH7bwvNcj4u3AKGBPYJc3Umye55SIGBcR44YPH/5GZ2dmZhXd/pXUSH4OHNSN5zwL3ATsDQyR1HHYahSwIHcvAEYD5PFbkk44m5lZH2n15rUPVXrXI9238GoXzxkO/G9EPCtpE+BA0snjm4DDSVcgTQCuyU+Znvtvz+Nv9PkEM7O+1erVR4dWulcAc0knhtdkW2CqpPVJQTItIq6V9CBwuaRvAvcAF+bpLwR+ImkOsAw4usXazMysl7R69dHE7s44Iu4F3tFk+OOk8wuNw18Fjuju65iZWe9p9eqjUZKulrQoP34qaVTdxZmZWd9q9UTzj0nH/LfLj1/kYWZmNoC0GgrDI+LHEbEiPy4CfD2omdkA02ooLJX0kXwz2vqSPoIvFzUzG3BaDYWPAUcCTwMLSZeMHl9TTWZm1iatXpL6f4EJEbEcQNJQ4AxSWJiZ2QDRakvh7zoCASAiltHkclMzM1u7tRoK60naqqMntxRabWWYmdlaotUv9jOB2yVdmfuPAL5VT0lmZtYurd7RfLGkWaT/hQDwoYh4sL6yzMysHVo+BJRDwEFgZjaAdfuns83MbOByKJiZWeFQMDOzwqFgZmaFQ8HMzAqHgpmZFQ4FMzMrHApmZlY4FMzMrHAomJlZ4VAwM7PCoWBmZoVDwczMCoeCmZkVDgUzMyscCmZmVtQWCpJGS7pJ0oOSHpB0ch4+VNIMSY/lv1vl4ZJ0nqQ5ku6VtEddtZmZWXN1thRWAJ+LiF2BvYATJe0KnAbMjIixwMzcD3AwMDY/JgHn11ibmZk1UVsoRMTCiPhj7n4BeAgYCYwHpubJpgKH5e7xwMWR3AEMkbRtXfWZmdnq+uScgqQxwDuAO4EREbEwj3oaGJG7RwLzKk+bn4c1zmuSpFmSZi1evLi2ms3M1kW1h4KkwcBPgc9ExPPVcRERQHRnfhExJSLGRcS44cOH92KlZmZWayhI2oAUCJdGxM/y4Gc6Dgvlv4vy8AXA6MrTR+VhZmbWR+q8+kjAhcBDEXFWZdR0YELungBcUxl+XL4KaS/gucphJjMz6wODapz3u4GPAvdJmp2HfRk4HZgm6QTgCeDIPO464BBgDvAyMLHG2szMrInaQiEibgXUyegDmkwfwIl11WNmZl3zHc1mZlY4FMzMrHAomJlZ4VAwM7PCoWBmZoVDwczMCoeCmZkVDgUzMyscCmZmVjgUzMyscCiYmVnhUDAzs8KhYGZmhUPBzMwKh4KZmRUOBTMzKxwKZmZWOBTMzKxwKJiZWeFQMDOzwqFgZmaFQ8HMzAqHgpmZFQ4FMzMrHApmZlY4FMzMrHAomJlZUVsoSPqRpEWS7q8MGypphqTH8t+t8nBJOk/SHEn3StqjrrrMzKxzdbYULgI+0DDsNGBmRIwFZuZ+gIOBsfkxCTi/xrrMzKwTtYVCRPwWWNYweDwwNXdPBQ6rDL84kjuAIZK2ras2MzNrrq/PKYyIiIW5+2lgRO4eCcyrTDc/DzMzsz7UthPNERFAdPd5kiZJmiVp1uLFi2uozMxs3dXXofBMx2Gh/HdRHr4AGF2ZblQetpqImBIR4yJi3PDhw2st1sxsXdPXoTAdmJC7JwDXVIYfl69C2gt4rnKYyczM+sigumYs6TJgP2CYpPnA14HTgWmSTgCeAI7Mk18HHALMAV4GJtZVl5mZda62UIiIYzoZdUCTaQM4sa5azMysNb6j2czMitpaCmZmZ894tN0lDFifPXDnWubrloKZmRUOBTMzKxwKZmZWOBTMzKxwKJiZWeFQMDOzwqFgZmaFQ8HMzAqHgpmZFQ4FMzMrHApmZlY4FMzMrHAomJlZ4VAwM7PCoWBmZoVDwczMCoeCmZkVDgUzMyscCmZmVjgUzMyscCiYmVnhUDAzs8KhYGZmhUPBzMwKh4KZmRUOBTMzK/pVKEj6gKRHJM2RdFq76zEzW9f0m1CQtD7wH8DBwK7AMZJ2bW9VZmbrln4TCsCewJyIeDwiXgMuB8a3uSYzs3XKoHYXUDESmFfpnw+8q3EiSZOASbn3RUmP9EFt/cEwYEm7i2jFKe0uoH9Ya9YXeJ1l69I626GzEf0pFFoSEVOAKe2uo69JmhUR49pdh7XG62vt43WW9KfDRwuA0ZX+UXmYmZn1kf4UCn8AxkraUdKGwNHA9DbXZGa2Tuk3h48iYoWk/wP8Glgf+FFEPNDmsvqTde6Q2VrO62vt43UGKCLaXYOZmfUT/enwkZmZtZlDwczMCoeCWYskDZH0qR487zpJQ+qoyVrXuP4k7Sfp2nbW1B85FNZykvrNxQLrgCHAaqHQ1TqIiEMi4tnaqrJWNV1/PTVQP3sDcqHWVpKOAz4PBHAvMA34KrAhsBQ4NiKekTQZ2An4G+BJ4Ji2FLzuOR3YSdJs4H+BV4HlwC7AzpJ+TrrXZmPg3HyjJZLmAuOAwcD1wK3APqT7cMZHxCt9vBzrBEmnAB/LvRcAe7Fy/c0AfgkMlnQV8LfA3cBHIiIkvRM4i7TOlgDHR8RCSTcDs4F9gcuAM/twkfpGRPjRDx7AbsCjwLDcPxTYipVXiH0cODN3TyZtwJu0u+516QGMAe7P3fsBLwE7VsYPzX83Ae4Hts79c0k/oTAGWAG8PQ+fRvoSavuyDbQH8E7gPmAz0hf7A8A7OtZfZR0+R7pRdj3gdtKX/QbA74HhebqjSJfIA9wM/LDdy1fnwy2F/mN/4MqIWAIQEcskvQ24QtK2pNbCnyvTTw/vYbbbXRFRXSeflvTB3D0aGEtq4VX9OSJm5+67SUFhvW9f4OqIeAlA0s+A9zSZ7q6ImJ+nmU1aH8+SWg4zJEG6b2ph5TlX1Fd2+zkU+rfvA2dFxHRJ+5FaCB1eaktFVlXWQV4/7wP2joiX82GGjZs85y+V7tdJrQprn8b1MQgQ8EBE7N3Jcwb0Z88nmvuPG4EjJG0NIGkosCUrf/9pQrsKs+IFYPNOxm0JLM+BsAvp+LW1z++AwyRtKmkz4IPAbXS+/qoeAYZL2htA0gaSdquv1P7FLYV+IiIekPQt4BZJrwP3kFoGV0paTgqNHdtY4jovIpZKuk3S/cArwDOV0b8CPinpIdKXyh3tqNGSiPijpIuAu/KgCyLi7sr6u550ornZc1+TdDhwnqQtSd+T55DOSwx4/pkLMzMrfPjIzMwKh4KZmRUOBTMzKxwKZmZWOBTMzKzwJalmvSTfYzIz925Duhlqce7fMyJea0thZt3gS1LNapB/tPDFiDij3bWYdYcPH5mZWeFQMDOzwqFgZmaFQ8HMzAqHgpmZFQ4FMzMrfEmqmZkVbimYmVnhUDAzs8KhYGZmhUPBzMwKh4KZmRUOBTMzKxwKZmZW/H/2CP9QoqnTXgAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "odb4KqIcu5cE"
      },
      "source": [
        "We can implement reverse inference using `pyro.condition` and an inference algorithm in Pyro. Assuming you observed a self-employed person who lives in a big city. What is your prediction of this person's age?  Let's generate a histogram of the marginal on the variable \"A\". "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "XQX3xUh8u5cE"
      },
      "source": [
        "conditioned_model = pyro.condition(model, data={'O':torch.tensor(1),'R':torch.tensor(1)})"
      ],
      "execution_count": 10,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 312
        },
        "id": "W8TO1iNxu5cE",
        "outputId": "a9575f47-4b70-4fdf-fd89-a8b6dfee9ed1"
      },
      "source": [
        "A_posterior = pyro.infer.Importance(conditioned_model, num_samples=5000).run()\n",
        "A_marginal = EmpiricalMarginal(A_posterior,\"A\")\n",
        "A_samples = [A_marginal().item() for _ in range(1000)]\n",
        "A_unique, A_counts = np.unique(A_samples, return_counts=True)\n",
        "\n",
        "plt.bar(A_unique, A_counts, align='center', alpha=0.5)\n",
        "plt.xticks(A_unique, A_alias)\n",
        "plt.ylabel('count')\n",
        "plt.xlabel('A')\n",
        "plt.title('P(A | O = self, R = big) - Importance Sampling')"
      ],
      "execution_count": 11,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "Text(0.5, 1.0, 'P(A | O = self, R = big) - Importance Sampling')"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 11
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEWCAYAAACJ0YulAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAapUlEQVR4nO3debhddX3v8fdH5goakIiQROJV1GJbEKOCxYqgKLQKrYCilUFq2qf2qrXOtRW9jlcramnppWLB4Yo4IEgtFRkcWkGDIqIopAqSyBAio0O94Pf+sX5nsTmeQ3aSs88+Sd6v59nPWeu31l77u8fP+v3W2vukqpAkCeB+4y5AkjR3GAqSpJ6hIEnqGQqSpJ6hIEnqGQqSpJ6hMMOS/GmS967H9U9NcswMljSrkuyXZMXA/KOSXJbkjiQvHWdtg5JUkkdMs+wFST4/5HaeleTjM1ud5qokxyf5SJt+aJI7k2w27rpmkqGwlpJck+Tn7cVwY/sQ37Yt2xJ4A/CuSdfZtq3/bzNw+1sleXuSH7U6rk7yqiRZ322PyKuBC6tqu6p6/5pWbh/WP22P18ok75ntN11VfbSqDhxy3c8Cj0nyO+t6e+019JZ1vf5MGvzQm0uS7JvkP5PcluQnSf4jyePHWVNV/aiqtq2qu8dZx0wzFNbNs6pqW2AvYAldEAAcAnyvqlZOWv85wH8DT0/ykPW87U8ABwAHA9sBLwSWAu9bz+2Oyq7Ad9byOnu0x/cpwHOBF814VTPrY3TPwQYtyebjrmEqSR4AnAP8PbADsAB4E917SjPMUFgP7cP/34Dfak0HAV+cYtWjgX8CLgf+eF1vL8kBwIHAc6rqiqq6q6oubtt8yXTDIet4W0lyQpKbktye5NtJfqst2yrJu1tv5cYk/5Rkmym2cQHwVODEtuf/yLWpoaqWA/8B7DkT92kKByf5QZKbk7wryf0AkhyT5CsTKyU5MMn3217qPyb5YpI/GdjORcDvz0RBSRa33tKxSa5LckuSP0vy+CSXJ7k1yYkD6x/T9ppPbPV9r71OJpbvkuTstne9PMmLB5Ydn+STST6S5Hbgz4DXA89tz9e32nrHJrmyDQH+IMmfDmxjvyQrkvxVe61cn+TYgeXbJPm7JNe2+r4y8VpJsnfb+781ybeS7DfNw/JIgKr6WFXdXVU/r6rPV9XlbTsPT3JBktXtufxoknkDNVzTetOXt17oKUl2SvJv7T59Icn2kx7/pUl+3O7PK9fwXG3e5i9K8r/a83FHks8n2XFg/aPa47A6yd+0up62xhfFLDMU1kOSRXR77N9sTb8NfH/SOrsC+wEfbZej1uMmnw5cUlXXDTZW1SXACroexFR1ntPeeFNdzpnmtg4Efo/uDflA4AhgdVv2jta+J/AIuj23v528garaH/gy8Betm33V2tzZJI8Gngwsv4911uW+TfhDup7eXnS9vF/rkbQ39SeB1wEPont+nzRptSuBxen2aGfKE4Hd6HpK7wX+Gnga8BjgiCRPmbTufwE7Am8EPp1kh7bsdLrXxi7AYcDbkuw/cN1D6O7fPOAU4G3Ax9vztUdb5ybgD4AHAMcCJyTZa2AbD6F7jSwAjgP+YeJDFng38Di6x2wHuuHEXyVZAPwr8JbW/krgU0nmT/FYXAXcneS0JAcNbHtCgLe3+/ibwCLg+EnrPIfu/fNI4Fl0O3OvB+bTfQ5OPt71VLrH/0DgNWvx4f18usfowcCW7X6RZHfgH4EXADtzz+M191SVl7W4ANcAdwK3AtfSPdHbtGVXA8+ctP4bgMva9ALgbuCx97H9U4Fjpln2AeD0aZZdDPz1DN7P/enejHsD9xtoD/BT4OEDbfsAP2zT+wErBpZdBPzJWtxuAbe32yi6oZmtRvA81uBzBfw5cH6bPgb4Sps+CvjqpPt/3eB9ArZo23voOtZyKvCWNr24bWvBwPLVwHMH5j8FvHyg1h8DGVj+NbphxUXt9bbdwLK3A6e26eOBL02q5XjgI2uo9zPAywae758Dmw8sv2niddOW7THFNl4DfHhS278DR09zm7/ZHqcVwF3A2cBO06x7KPDNgflrgBdMevxOGpj/n8BnJj3+jx5Y/r+BUyY/PgPrbj7wWn/DpNfUuW36b4GPDSz7DeCXwNNm+rW9vhd7Cuvm0KqaV1W7VtWfV9XPW/stdOP8g46i6yFQ3XDTF+mGk9bFzXR7GVPZuS2fEVV1AXAi8A/ATUlObnvC8+le0JdO7JED57b2mbIXsC3dXvITgfvP4LYHDfa4rqXb05xsl8H1qntHr5i0zsRzfuvkK6c7k+nOrP2JBjcOTP98ivltB+ZXtromTNyXXYCfVNUdk5YN7qHeq9c5lbZ3fnEbgrqVrne848Aqq6vqroH5n7X6dgS2puvFTLYrcPhgzw7Yl2le31V1ZVUdU1UL6YZrd6HrQdGGgk5Pd2LC7cBHJtUHa/d4wnCvjancMDA98TjAr7+OfsY9Pe85xVCYWZfTxj8BkjyJrgv6uiQ3JLmB7kPu+Vm3g3pfAJ7Yhq16SZ5It1d4wVRXamOnd05zmfaDqqreX1WPA3Zv9+tVdMHzc+AxLRjnVdUDqzswPGOqcwbwVaYYmlrf+9YMPo4Ppdvjnux6YOHA7WVwvvlN4Jqqun2K+/HR6oZitq2qg9ZQz7pa0OqaMHFffgzskGS7ScsGT4SY/DPJ95pPshXdnvW76fbM5wGfo+sxrcnNwC+Ah0+x7Dq6nsK8gcv9q+oda9poVX2PrtcwcSzvba3u366qB9AdY1vfs/GGeW2sjcmvo23ohiPnHENhZn2O7oyZCUcD59F9qO7ZLr8FbEN3UHqtVNUXgPPpxl4fk2SzJHvT7RmdVFVXT3O9gwY+mCZfpqwj3YHNJybZgm4o5xfAr6rqV8A/040rP7ituyDJM4a5D+kOjF6zFnf7HcCLM81ZW+ty3wa8Ksn2LWRfBkz1fYN/BX47yaEtyF9CN4Y+6Cl0Y9Tj8mDgpUm2SHI4XUh9rrpjT/8JvD3J1ulOmz2O7vUynRvpjo9MfDZsCWwFrALuSnIQ3Tj7GrXXygeB96Q74L1Zkn1a0HwEeFaSZ7T2rdMdtJ4cuCR5dLoD2Qvb/CLgSLohU+h6ancCt7VjFa8apr41+Jskv5HkMXTHCNb3uyifpLu/T0p36vrxrH9wjYShMLM+Czy6vQG2pjs4+/dVdcPA5YfAh1n3IaTnABfSDdncSffmOoVuXHQmPYDuw/8Wuu7zau75/sVr6A7+Xty6618AHjXkdhfRnVE0lKr6NvAlZuaNPtlZwKXAZXQf/qdMcfs3A4fTjSuvpgv4Zdz7dMgjgf8zgvqGdQldj/Rm4K3AYVU1MTRxJN3Y94+BM4E3tp2L6Xyi/V2d5Btt6OmlwBl0r4Xn043nD+uVwLeBrwM/Ad5Jd4zqOrqD3K+nC5zr6J7jqT6T7qDrYV+S5Kd0YXAF8Fdt+Zvohhxvo3seP70W9U3ni3Sv8fOBd1fVUF9mnE5VfYfuPXo6Xa/hTrpjL3PutNrceyhS6yvJUmD3qnr5Ol7/VOCiqjp1JuuaK9J9U/hlVXXluGtZF20PegXdgcsLkzwLeGFVHTGmeo6hO+i97zhuf2OTZDHwQ2CLScdJZvp2tqU7BrVb21GcM+bkl1U2ZFV18rhrmMtqyG8KzyVtaOwSumMpr6Lr9l8M/TeaPzu+6rShaDsQ59O9ft5N14O6Zpw1TcXho7nnM3TDGZo79qE7g+ZmunPcDx0440wa1iHccwLAbsDzag4O1Th8JEnq2VOQJPU26GMKO+64Yy1evHjcZUjSBuXSSy+9uaqm/MLpBh0KixcvZtmyZeMuQ5I2KEmunW6Zw0eSpJ6hIEnqGQqSpJ6hIEnqGQqSpJ6hIEnqGQqSpJ6hIEnqGQqSpN4G/Y1mbVpOOO+qcZew0frLpz9yzStpk2BPQZLUMxQkST1DQZLUMxQkST1DQZLUMxQkST1DQZLUMxQkSb2RhkKSa5J8O8llSZa1th2SnJfk6vZ3+9aeJO9PsjzJ5Un2GmVtkqRfNxs9hadW1Z5VtaTNvxY4v6p2A85v8wAHAbu1y1LgpFmoTZI0YBzDR4cAp7Xp04BDB9o/VJ2LgXlJdh5DfZK0yRp1KBTw+SSXJlna2naqquvb9A3ATm16AXDdwHVXtLZ7SbI0ybIky1atWjWquiVpkzTqH8Tbt6pWJnkwcF6S7w0urKpKUmuzwao6GTgZYMmSJWt1XUnSfRtpT6GqVra/NwFnAk8AbpwYFmp/b2qrrwQWDVx9YWuTJM2SkYVCkvsn2W5iGjgQuAI4Gzi6rXY0cFabPhs4qp2FtDdw28AwkyRpFoxy+Ggn4MwkE7fzf6vq3CRfB85IchxwLXBEW/9zwMHAcuBnwLEjrE2SNIWRhUJV/QDYY4r21cABU7QX8JJR1SNJWjO/0SxJ6hkKkqSeoSBJ6hkKkqSeoSBJ6hkKkqSeoSBJ6hkKkqSeoSBJ6hkKkqSeoSBJ6hkKkqSeoSBJ6hkKkqSeoSBJ6hkKkqSeoSBJ6hkKkqSeoSBJ6hkKkqSeoSBJ6hkKkqSeoSBJ6hkKkqSeoSBJ6hkKkqSeoSBJ6hkKkqSeoSBJ6hkKkqTeyEMhyWZJvpnknDb/sCSXJFme5ONJtmztW7X55W354lHXJkm6t9noKbwMuHJg/p3ACVX1COAW4LjWfhxwS2s/oa0nSZpFIw2FJAuB3wc+0OYD7A98sq1yGnBomz6kzdOWH9DWlyTNklH3FN4LvBr4VZt/EHBrVd3V5lcAC9r0AuA6gLb8trb+vSRZmmRZkmWrVq0aZe2StMkZWSgk+QPgpqq6dCa3W1UnV9WSqloyf/78mdy0JG3yNh/htn8XeHaSg4GtgQcA7wPmJdm89QYWAivb+iuBRcCKJJsDDwRWj7A+SdIkI+spVNXrqmphVS0GngdcUFUvAC4EDmurHQ2c1abPbvO05RdUVY2qPknSrxvH9xReA7wiyXK6YwantPZTgAe19lcArx1DbZK0SRvl8FGvqi4CLmrTPwCeMMU6vwAOn416JElT8xvNkqSeoSBJ6hkKkqSeoSBJ6hkKkqSeoSBJ6hkKkqSeoSBJ6hkKkqSeoSBJ6hkKkqSeoSBJ6hkKkqSeoSBJ6hkKkqSeoSBJ6hkKkqSeoSBJ6hkKkqSeoSBJ6hkKkqSeoSBJ6hkKkqSeoSBJ6hkKkqSeoSBJ6hkKkqSeoSBJ6hkKkqSeoSBJ6g0VCknOH6Zt0vKtk3wtybeSfCfJm1r7w5JckmR5ko8n2bK1b9Xml7fli9f+7kiS1sd9hkL7YN8B2DHJ9kl2aJfFwII1bPu/gf2rag9gT+CZSfYG3gmcUFWPAG4BjmvrHwfc0tpPaOtJkmbR5mtY/qfAy4FdgEuBtPbbgRPv64pVVcCdbXaLdilgf+D5rf004HjgJOCQNg3wSeDEJGnbmXEnnHfVKDYr4C+f/shxlyBpHd1nT6Gq3ldVDwNeWVX/o6oe1i57VNV9hgJAks2SXAbcBJwH/Bdwa1Xd1VZZwT09jgXAde127wJuAx40xTaXJlmWZNmqVauGvJuSpGGsqacAQFX9fZInAYsHr1NVH1rD9e4G9kwyDzgTePS6l9pv82TgZIAlS5aMpBchSZuqoUIhyYeBhwOXAXe35gLuMxQmVNWtSS4E9gHmJdm89QYWAivbaiuBRcCKJJsDDwRWD3tHJEnrb6hQAJYAu6/N+H6S+cD/a4GwDfB0uoPHFwKHAacDRwNntauc3ea/2pZfMKrjCZKkqQ0bClcADwGuX4tt7wyclmQzumMXZ1TVOUm+C5ye5C3AN4FT2vqnAB9Oshz4CfC8tbgtSdIMGDYUdgS+m+RrdKeaAlBVz57uClV1OfDYKdp/ADxhivZfAIcPWY8kaQSGDYXjR1mEJGluGPbsoy+OuhBJ0vgNe/bRHXRnGwFsSfdFtJ9W1QNGVZgkafYN21PYbmI6Sei+fbz3qIqSJI3HWv9KanU+AzxjBPVIksZo2OGjPxqYvR/d9xZ+MZKKJEljM+zZR88amL4LuIZuCEmStBEZ9pjCsaMuRJI0fsP+k52FSc5MclO7fCrJwlEXJ0maXcMeaP4Xut8m2qVdPtvaJEkbkWFDYX5V/UtV3dUupwLzR1iXJGkMhg2F1Un+uP3TnM2S/DH+rLUkbXSGDYUXAUcAN9D9UuphwDEjqkmSNCbDnpL6ZuDoqroFIMkOwLvpwkKStJEYtqfwOxOBAFBVP2GKn8WWJG3Yhg2F+yXZfmKm9RSG7WVIkjYQw36w/x3w1SSfaPOHA28dTUmSpHEZ9hvNH0qyDNi/Nf1RVX13dGVJksZh6CGgFgIGgSRtxNb6p7MlSRsvQ0GS1DMUJEk9Q0GS1DMUJEk9Q0GS1DMUJEk9Q0GS1DMUJEk9Q0GS1DMUJEm9kYVCkkVJLkzy3STfSfKy1r5DkvOSXN3+bt/ak+T9SZYnuTzJXqOqTZI0tVH2FO4C/qqqdgf2Bl6SZHfgtcD5VbUbcH6bBzgI2K1dlgInjbA2SdIURhYKVXV9VX2jTd8BXAksAA4BTmurnQYc2qYPAT5UnYuBeUl2HlV9kqRfNyvHFJIspvv3nZcAO1XV9W3RDcBObXoBcN3A1Va0tsnbWppkWZJlq1atGlnNkrQpGnkoJNkW+BTw8qq6fXBZVRVQa7O9qjq5qpZU1ZL58+fPYKWSpJGGQpIt6ALho1X16dZ848SwUPt7U2tfCSwauPrC1iZJmiWjPPsowCnAlVX1noFFZwNHt+mjgbMG2o9qZyHtDdw2MMwkSZoFQ/87znXwu8ALgW8nuay1vR54B3BGkuOAa4Ej2rLPAQcDy4GfAceOsDZJ0hRGFgpV9RUg0yw+YIr1C3jJqOqRJK2Z32iWJPUMBUlSz1CQJPUMBUlSz1CQJPUMBUlSz1CQJPVG+eU1SZu4E867atwlbLT+8umPHMl27SlIknqGgiSpZyhIknqGgiSpZyhIknqGgiSpZyhIknqGgiSpZyhIknqGgiSpZyhIknqGgiSpZyhIknqGgiSpZyhIknqGgiSpZyhIknqGgiSpZyhIknqGgiSpZyhIknqGgiSpN7JQSPLBJDcluWKgbYck5yW5uv3dvrUnyfuTLE9yeZK9RlWXJGl6o+wpnAo8c1Lba4Hzq2o34Pw2D3AQsFu7LAVOGmFdkqRpjCwUqupLwE8mNR8CnNamTwMOHWj/UHUuBuYl2XlUtUmSpjbbxxR2qqrr2/QNwE5tegFw3cB6K1rbr0myNMmyJMtWrVo1ukolaRM0tgPNVVVArcP1Tq6qJVW1ZP78+SOoTJI2XbMdCjdODAu1vze19pXAooH1FrY2SdIsmu1QOBs4uk0fDZw10H5UOwtpb+C2gWEmSdIs2XxUG07yMWA/YMckK4A3Au8AzkhyHHAtcERb/XPAwcBy4GfAsaOqS5I0vZGFQlUdOc2iA6ZYt4CXjKoWSdJw/EazJKlnKEiSeoaCJKlnKEiSeoaCJKlnKEiSeoaCJKlnKEiSeoaCJKlnKEiSeoaCJKlnKEiSeoaCJKlnKEiSeoaCJKlnKEiSeoaCJKlnKEiSeoaCJKlnKEiSeoaCJKlnKEiSeoaCJKlnKEiSeoaCJKlnKEiSeoaCJKlnKEiSeoaCJKlnKEiSenMqFJI8M8n3kyxP8tpx1yNJm5o5EwpJNgP+ATgI2B04Msnu461KkjYtcyYUgCcAy6vqB1X1S+B04JAx1yRJm5TNx13AgAXAdQPzK4AnTl4pyVJgaZu9M8n3Z6G2uWBH4OZxFzGMV4y7gLlhg3m+wOes2ZSes12nWzCXQmEoVXUycPK465htSZZV1ZJx16Hh+HxteHzOOnNp+GglsGhgfmFrkyTNkrkUCl8HdkvysCRbAs8Dzh5zTZK0SZkzw0dVdVeSvwD+HdgM+GBVfWfMZc0lm9yQ2QbO52vD43MGpKrGXYMkaY6YS8NHkqQxMxQkST1DQZohSY5JcuIa1lmc5Io2vWeSg2enOt2XJHdO035qksNmu55xMhSk8dkTMBQ0pxgKsyjJm5O8fGD+rUleluRdSa5I8u0kz23L9ktyzsC6JyY5pk1fk+RNSb7RrvPo1j4/yXlJvpPkA0muTbLjLN/NjVaSzyS5tD2+S1vbsUmuSvI14HcH1r3XHubkPdF22vWbgecmuWziedfoJXlFe79dMfh+bMvS3mvfT/IF4MFjKnNsDIXZ9UHgKIAk96P7LsYKuj3GPYCnAe9KsvMQ27q5qvYCTgJe2dreCFxQVY8BPgk8dGbL3+S9qKoeBywBXppkAfAmujDYl+6HHIfSft/rb4GPV9WeVfXxURSse0vyOOBYup/Q2Rt4cZLHDqzyh8Cj6J7Lo4AnzXqRY2YozKKqugZY3V6EBwLfpPsw+VhV3V1VNwJfBB4/xOY+3f5eCixu0/vS/ZAgVXUucMuMFS/oguBbwMV0375/IXBRVa1qH/J+sM99+wJnVtVPq+pOuvfRkweW/x73vB9/DFwwjiLHyVCYfR8AjqHbW/ngfax3F/d+fraetPy/29+7mUNfQtxYJdmPrie3T1XtQRfo37uPq/TPX+sVbjnqGqWZYCjMvjOBZ9L1Bv4d+DLduPJmSebT7al8DbgW2D3JVknmAQcMse3/AI4ASHIgsP0I6t9UPRC4pap+1o7h7A1sAzwlyYOSbAEcPrD+NcDj2vSzgS2m2OYdwHajK1lT+DJwaJLfSHJ/uuGiLw8s/xL3vB93Bp46jiLHyT3MWVZVv0xyIXBrVd2d5ExgH+BbQAGvrqobAJKcAVwB/JBuz3RN3gR8LMkLga8CN9B98Gj9nQv8WZIrge/TDSFdDxxP91jfClw2sP4/A2e14aZzgZ9Osc0LgdcmuQx4u8cVRq+qvpHkVLodL4APVNU3k0ysciawP/Bd4Ed0z+0mxZ+5mGVtKOEbwOFVdfUMb3sr4O72O1L7ACdV1Z4zeRuSNm72FGZR+/ei59Ad6JrRQGgeCpzRgueXwItHcBuSNmL2FCRJPQ80S5J6hoIkqWcoSJJ6hoI0w5IcmqQmfpNK2pAYCtLMOxL4SvsrbVA8+0iaQUm2pfty21OBz1bVo8ZckrRW7ClIM+sQ4Nyquoruxw8ft6YrSHOJoSDNrCNpv1Tb/jqEpA2Kw0fSDEmyA93/x1hF9ztWm7W/u5ZvNG0g7ClIM+cw4MNVtWtVLa6qRXQ/ZvjkNVxPmjMMBWnmHEn3K5uDPoVDSNqAOHwkSerZU5Ak9QwFSVLPUJAk9QwFSVLPUJAk9QwFSVLPUJAk9f4/xc6P649Tlu4AAAAASUVORK5CYII=\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "0YrBPZunvm_-"
      },
      "source": [
        "### Intervention\n",
        "\n",
        "Finally, we'll use the `do` operator to simulate intervention.  We'll set residence to \"small\"."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 312
        },
        "id": "TmiZD4okvlLG",
        "outputId": "4ebc1738-d5c3-4d21-8d00-a739f76ba468"
      },
      "source": [
        "intervention_model = pyro.do(model, data={'R':torch.tensor(0)})\n",
        "\n",
        "T_posterior = pyro.infer.Importance(intervention_model, num_samples=5000).run()\n",
        "T_marginal = EmpiricalMarginal(T_posterior,\"T\")\n",
        "T_samples = [T_marginal().item() for _ in range(1000)]\n",
        "T_unique, T_counts = np.unique(T_samples, return_counts=True)\n",
        "\n",
        "plt.bar(T_unique, T_counts, align='center', alpha=0.5)\n",
        "plt.xticks(T_unique, T_alias)\n",
        "plt.ylabel('count')\n",
        "plt.xlabel('T')\n",
        "plt.title('P(T | do(R = small)) - Importance Sampling')"
      ],
      "execution_count": 12,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "Text(0.5, 1.0, 'P(T | do(R = small)) - Importance Sampling')"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 12
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEWCAYAAACJ0YulAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAZm0lEQVR4nO3debgddZ3n8feHTVSUAImoSTAoOIrtnkFstcVtWnABfZTWAUEaOzIu3W7j2I49E7tbR20FUVtHFBXUVnBHRVsGxXZDDUori0tUaIgIAQICigp+54/63eIkuZecQOqem+T9ep7z3Kpf1an7vafOrc/51XIqVYUkSQDbTLoASdLcYShIknqGgiSpZyhIknqGgiSpZyhIknqGwoQkeX6St97K5y5Psnwj5j8zyfM2Yv5vJHnwraltrktyYZLHt+HlST7UhndPckGS243M++Ikb5xUrZpdST6Q5B/b8KOS/HjSNU2CoTCQtvH5bZLrklzW3nA7tWk7AK8B/qm9+a5rj+uT1Mj4dUn2mOW6nwJcW1Xfb+PLk/yh1XJ1km8mefhs1jQbquoy4CvAspHm9wCHJrnLrV3uxgbykEY3enNJkoOSnJPk10muSPLlJHtOsqaq+lpV/adJ1jAphsKwnlJVOwEPAZbSBQHAQcCPqmpVe/Pt1Oa7X5s+b6qtqv5jlms+GvjgOm0nt/rm0204PzbLNc2WDwPPnxqpqhuALwCHT6yiTSTJtpOuYTpJ9gJOAl4O7AzsCfwzcNMk69qaGQqzoKpW0W1c/qQ1HQB8dajfl+QJSX6U5Jok7wAyMm2bJK9JclGSy5OclGTnNm0H4LEz1VZVN9JtOBcmWbCJa943yYr2afGyJMe09iWt93RkkouTrElydJL/nOQHrffyjpHl3Kt90ryyfer8cJJ5Y5bxbeCeSe4x0nYm8KRN9Dfun+SSJK9sr/2lSQ5OcmCSnyS5KsmrR+ZfnuTjSU5Ocm2S7yV54Mj0+7aeyNVJzkvy1JFpH0jyriSnJbkeOAo4FHhl6/V9ts33qiQ/a8s/P8nTRpbx3CRfT/Lm9rr/IskBI9N3TfL+JL9s0z89Mu3J7dP/VO/yATO8LA8CflFVZ1Tn2qr6xNSHofa++FZbzqVJ3tHep1O/p5K8IMlP29/wD+098M32Xjplav6R1//V7b1xYZJDb2ldjYxfmOQV7T13TVsnO45Mf2Wr75dJntfq2usW3g5zlqEwC5IsBg4Evt+a7g8Msr8yyXzgk3S9kvnAz4BHjMzy3PZ4DHBPYCdgaqO6N/DHqrqEabR/rsOBK4E1M8zzyPYPPNPjkTOUfhxwXFXdGbgXcMo60x/W6vsL4K3A/wQeT9e7OiTJo6dKAP4PcHfgvsBiYPkMv3MtLfRWAg8cab5gnfHb6q7AjsBC4H/R7aI6DHgo8Cjg77L2rpOD6HpmuwL/Anw6yfZJtgc+C3wJuAvwYuDDSUZ3efxX4HXAneg+jX8YeFPrgT6lzfOz9nt3Bl4LfCjJ3UaW8TC69+p84E3ACUmmPmR8ELgD3Tq4C3AsQLrjUe+j63XtBrwbODUjx2tGfA+4T5JjkzwmbRfriJuAl7bf/3DgccAL1pnnz9vrtx/wSuD49poupvsg9uyRee/alrUQOAI4fp3X7JYcAjyRrjfzALr/I5I8EXgZ3ftxL2D/MZc3JxkKw/p0kquBr9N9+n59a58HXDvQ7zwQOK+qPl5Vf6DbgP5qZPqhwDFV9fOqug74W+BZSba7hboOaX/Hb4G/Ap7RNqDrqaqvV9W8W3h8fYa6/wDslWR+VV1XVWetM/0fquqGqvoScD3wkaq6vPXCvgY8uP3+lVV1elX9rqpWA8cAj2Z817bXYXR85414/ob8AXhdWzcfpdtAHdc+IZ8HnM/aIXT2yLo8hi5Q9muPnYA3VNXvq+rLwOdYewP4mar6RlX9se0KW09VfayqftnmORn4KbDvyCwXVdV7quom4ETgbsDuLTgOAI6uqjVV9YeqmuphLgPeXVXfrqqbqupE4Het5nV//8/pNqIL6T4IXJGR429VdXZVnVVVN1bVhXQBs+76fFNV/bq9fucCX2rv72voeujrnjTxd+398VXg83Qb+3G8rb1WV9EF8oNa+yHA+6vqvKr6DWN+CJmrDIVhHdw2hPeoqhdU1W9b+xq6T29DuDtw8dRIdd94ePE60y8aGb8I2A7Y/RbqOqWq5rV5zqX7VLapHQXcG/hRku8mefI60y8bGf7tNONTB/F3T/LRJKuS/Br4EN2Gd1x3Aq5eZ/ya6WZsuyGmTgj4v2Mu/8q2gZ2qG2b4W5rRdflH4BK6dXh34OLWNuUiuo3res+dSZLDR3bzXE33yXr09eo/ULQNHq2+xcBVVTVdj/EewMtHe4ht/rtPV0Pb6B9SVQvoei1/RtcTJMm9k3wuya/a+nw966/Psd4bzZqqun5k/KKZ6prG6Ier34wsd63/OcZ43ecyQ2EyfkC3ARzCpXT/gAC0rv7ikem/pPunnbIHcCPdP9LK9pTRDUuvqq6g+xS4fJ1dDL2sfTbVdI9HzbDsn1bVs+l2Q7wR+HiSO477R494PVDA/duuqMMYOaZyS1pvaS/g30ea77vO+GjNrx85IeDoW1HrOEbX5TbAIrp1+EtgcWubsgewarTEdUseHUl37OQ9wIuA3Vrwn8t4r9fFwK6Z/njNxXS9odEe4h2q6iMbWmhVfZdu9+fU8bd3AT8C9m7r89Vj1jeTXdZ5X+1B91reFpfSrZcpi2eacXNgKEzGaWzcLo2N8Xngfkme3jZyf023H3XKR4CXJtmzddFfT3d20Y1V9Xvg/91SbVX1Y+Bf6fbdTjf9ayMbyukeX5vueUkOS7KgffKd+qT+x+nm3YA7AdcB17Rw++8b8dx9gQurarQn9Wi6XRCT8tCRdfkSut0wZ9EdFP8N3YHj7ZPsDzyFbpfUTC6jO4405Y50QbEaIMmR3LwxvkVVdSnd6/LOJLu0Gv6sTX4PcHSSh6VzxyRPSrJeL7Qdg/qrtNN+k9wHeGr7G6Fbn78GrmvT/ts49W3Aa5Ps0D6gPJnbfjbdKcCR6Q783wH4u9tc4QQZCpPxWbqDa+N2W8fWPs0/E3gD3QHhvYFvjMzyProDhP8G/AK4ge4g5ZR3A8/ZwK/5J2BZbsP5+9N4InBekuvoDjo/a2R328Z4Ld0pwNfQBeQnN+K5hwL9bqB2dsmBdPvSJ+UzdAfX19Ctl6e3/fe/pwuBA4ArgHcCh1fVj25hWScA+7RdOp+uqvOBtwDfoguM+7P2e2VDnkN3jORHwOV0oUVVraA79vSOVvdK2kHZaVxNFwI/bOv+i8Cn6A5qA7yC7oD5tXRhc/JG1DedX7Wafkl34P3oDbxmG1RVXwDeRne69kpuDrTf3ZblTkrKm+xMRJJlwD5V9ZJb8dzlAFW1fBOXNbX8bwAvqnYB29agBdxXgQdPHZRN8mJgcVVN2yuahZqWA3tV1WGT+P1bmtab+lBVLdrQvLfx99yXbjfc7WY6IWMu227SBWytqur4Sdcwk6p6xIbn2rJU1eV0xw9G294+oXK0mUl3fcdpdKfovhH47OYYCODuo83Vme0haW54Pt0utJ/RXVuxKY59TIS7jyRJPXsKkqTeZn1MYf78+bVkyZJJlyFJm5Wzzz77inax4Ho261BYsmQJK1asmHQZkrRZSXLRTNPcfSRJ6hkKkqSeoSBJ6g0aCuluTPHD9i2MK1rbrklOT3dTjNOT7NLak+RtSVamu5HFQ4asTZK0vtnoKTymqh5UVUvb+KuAM6pqb+CMNg7dd7js3R7L6L4dUZI0iyax++ggbv6CsROBg0faT6rOWcC8mb6eWZI0jKFDoYAvJTm7fQEcwO7ta3eh+8bC3dvwQta+OcUlrH3DEEnSwIa+TuGRVbWqfQPl6UnW+oraqqokG/U9Gy1clgHssccem65SSdKwPYV2/9ypb6D8FN1NTC6b2i3Ufl7eZl/F2ncsWsTad5GaWubxVbW0qpYuWDDtBXmSpFtpsJ5Cu+XdNlV1bRv+L8DfA6cCR9DdBOYIupuI0NpflOSjwMOAa0Z2M21yx57+k6EWvdV76ROGutOopKENuftod+BT3S2C2Q74l6r6YpLvAqckOYruptmHtPlPo7vL1Uq62wweOWBtkqRpDBYKVfVz4IHTtF8JPG6a9gJeOFQ9kqQN84pmSVLPUJAk9QwFSVLPUJAk9QwFSVLPUJAk9QwFSVLPUJAk9QwFSVLPUJAk9QwFSVLPUJAk9QwFSVLPUJAk9QwFSVLPUJAk9QwFSVLPUJAk9QwFSVLPUJAk9QwFSVLPUJAk9QwFSVLPUJAk9QwFSVLPUJAk9QwFSVLPUJAk9QwFSVLPUJAk9babdAHSuI49/SeTLmGL9dIn3HvSJWiOsKcgSeoZCpKk3uChkGTbJN9P8rk2vmeSbydZmeTkJDu09tu18ZVt+pKha5MkrW02egp/A1wwMv5G4Niq2gtYAxzV2o8C1rT2Y9t8kqRZNGgoJFkEPAl4bxsP8Fjg422WE4GD2/BBbZw2/XFtfknSLBm6p/BW4JXAH9v4bsDVVXVjG78EWNiGFwIXA7Tp17T515JkWZIVSVasXr16yNolaaszWCgkeTJweVWdvSmXW1XHV9XSqlq6YMGCTbloSdrqDXmdwiOApyY5ENgRuDNwHDAvyXatN7AIWNXmXwUsBi5Jsh2wM3DlgPVJktYxWE+hqv62qhZV1RLgWcCXq+pQ4CvAM9psRwCfacOntnHa9C9XVQ1VnyRpfZO4TuF/AC9LspLumMEJrf0EYLfW/jLgVROoTZK2arPyNRdVdSZwZhv+ObDvNPPcADxzNuqRJE3PK5olST1DQZLUMxQkST1DQZLUMxQkST1DQZLUMxQkST1DQZLUMxQkST1DQZLUMxQkST1DQZLUMxQkST1DQZLUMxQkST1DQZLUMxQkST1DQZLUMxQkST1DQZLUMxQkST1DQZLUMxQkST1DQZLUMxQkST1DQZLUMxQkST1DQZLUMxQkST1DQZLUMxQkST1DQZLUMxQkSb3BQiHJjkm+k+Tfk5yX5LWtfc8k306yMsnJSXZo7bdr4yvb9CVD1SZJmt6QPYXfAY+tqgcCDwKemGQ/4I3AsVW1F7AGOKrNfxSwprUf2+aTJM2iwUKhOte10e3bo4DHAh9v7ScCB7fhg9o4bfrjkmSo+iRJ6xv0mEKSbZOcA1wOnA78DLi6qm5ss1wCLGzDC4GLAdr0a4DdplnmsiQrkqxYvXr1kOVL0lZn0FCoqpuq6kHAImBf4D6bYJnHV9XSqlq6YMGC21yjJOlms3L2UVVdDXwFeDgwL8l2bdIiYFUbXgUsBmjTdwaunI36JEmdIc8+WpBkXhu+PfAE4AK6cHhGm+0I4DNt+NQ2Tpv+5aqqoeqTJK1vuw3PcqvdDTgxybZ04XNKVX0uyfnAR5P8I/B94IQ2/wnAB5OsBK4CnjVgbZKkaQwWClX1A+DB07T/nO74wrrtNwDPHKoeSdKGjbX7KMkZ47RJkjZvt9hTSLIjcAdgfpJdgKnrBu7MzaeSSpK2EBvaffR84CXA3YGzuTkUfg28Y8C6JEkTcIuhUFXHAccleXFVvX2WapIkTchYB5qr6u1J/hRYMvqcqjppoLokSRMwVigk+SBwL+Ac4KbWXIChIElbkHFPSV0K7OPFZJK0ZRv3iuZzgbsOWYgkafLG7SnMB85P8h26+yQAUFVPHaQqSdJEjBsKy4csQpI0N4x79tFXhy5EkjR54559dC3d2UYAO9DdRe36qrrzUIVJkmbfuD2FO00Nt1tkHgTsN1RRkqTJ2Oj7KbR7L38a+PMB6pEkTdC4u4+ePjK6Dd11CzcMUpEkaWLGPfvoKSPDNwIX0u1CkiRtQcY9pnDk0IVIkiZv3JvsLEryqSSXt8cnkiwaujhJ0uwa90Dz+4FT6e6rcHfgs61NkrQFGTcUFlTV+6vqxvb4ALBgwLokSRMwbihcmeSwJNu2x2HAlUMWJkmafeOGwl8ChwC/Ai4FngE8d6CaJEkTMu4pqX8PHFFVawCS7Aq8mS4sJElbiHF7Cg+YCgSAqroKePAwJUmSJmXcUNgmyS5TI62nMG4vQ5K0mRh3w/4W4FtJPtbGnwm8bpiSJEmTMu4VzSclWQE8tjU9varOH64sSdIkjL0LqIWAQSBJW7CN/upsSdKWy1CQJPUMBUlSz1CQJPUMBUlSb7BQSLI4yVeSnJ/kvCR/09p3TXJ6kp+2n7u09iR5W5KVSX6Q5CFD1SZJmt6QPYUbgZdX1T7AfsALk+wDvAo4o6r2Bs5o4wAHAHu3xzLgXQPWJkmaxmChUFWXVtX32vC1wAXAQrp7O5/YZjsROLgNHwScVJ2zgHlJ7jZUfZKk9c3KMYUkS+i+QO/bwO5VdWmb9Ctg9za8ELh45GmXtLZ1l7UsyYokK1avXj1YzZK0NRo8FJLsBHwCeElV/Xp0WlUVUBuzvKo6vqqWVtXSBQu8+ZskbUqDhkKS7ekC4cNV9cnWfNnUbqH28/LWvgpYPPL0Ra1NkjRLhjz7KMAJwAVVdczIpFOBI9rwEcBnRtoPb2ch7QdcM7KbSZI0C4a8J8IjgOcAP0xyTmt7NfAG4JQkRwEX0d3mE+A04EBgJfAb4MgBa5MkTWOwUKiqrwOZYfLjppm/gBcOVY8kacO8olmS1DMUJEk9Q0GS1DMUJEk9Q0GS1DMUJEk9Q0GS1DMUJEk9Q0GS1DMUJEk9Q0GS1DMUJEk9Q0GS1DMUJEk9Q0GS1DMUJEk9Q0GS1DMUJEk9Q0GS1DMUJEk9Q0GS1DMUJEk9Q0GS1DMUJEk9Q0GS1DMUJEk9Q0GS1DMUJEk9Q0GS1DMUJEk9Q0GS1DMUJEk9Q0GS1BssFJK8L8nlSc4dads1yelJftp+7tLak+RtSVYm+UGShwxVlyRpZkP2FD4APHGdtlcBZ1TV3sAZbRzgAGDv9lgGvGvAuiRJMxgsFKrq34Cr1mk+CDixDZ8IHDzSflJ1zgLmJbnbULVJkqY328cUdq+qS9vwr4Dd2/BC4OKR+S5pbetJsizJiiQrVq9ePVylkrQVmtiB5qoqoG7F846vqqVVtXTBggUDVCZJW6/ZDoXLpnYLtZ+Xt/ZVwOKR+Ra1NknSLJrtUDgVOKINHwF8ZqT98HYW0n7ANSO7mSRJs2S7oRac5CPA/sD8JJcA/xt4A3BKkqOAi4BD2uynAQcCK4HfAEcOVZckaWaDhUJVPXuGSY+bZt4CXjhULZIm49jTfzLpErZYL33CvQdZrlc0S5J6hoIkqWcoSJJ6hoIkqWcoSJJ6hoIkqWcoSJJ6hoIkqWcoSJJ6hoIkqWcoSJJ6hoIkqWcoSJJ6hoIkqWcoSJJ6hoIkqWcoSJJ6hoIkqWcoSJJ6hoIkqWcoSJJ6hoIkqWcoSJJ6hoIkqWcoSJJ6hoIkqWcoSJJ6hoIkqWcoSJJ6hoIkqWcoSJJ6hoIkqWcoSJJ6cyoUkjwxyY+TrEzyqknXI0lbmzkTCkm2Bf4ZOADYB3h2kn0mW5UkbV3mTCgA+wIrq+rnVfV74KPAQROuSZK2KttNuoARC4GLR8YvAR627kxJlgHL2uh1SX48C7XNBfOBKyZdxDheNukC5obNZn2B66zZmtbZPWaaMJdCYSxVdTxw/KTrmG1JVlTV0knXofG4vjY/rrPOXNp9tApYPDK+qLVJkmbJXAqF7wJ7J9kzyQ7As4BTJ1yTJG1V5szuo6q6McmLgH8FtgXeV1XnTbisuWSr22W2mXN9bX5cZ0CqatI1SJLmiLm0+0iSNGGGgiSpZyhIY0oyL8kLbsXzTksyb4iaNL5111+S/ZN8bpI1zUWGwmYuyZw5WWArMA9YLxQ2tA6q6sCqunqwqjSuadffrbWl/u9tkX/U5irJ4cArgAJ+AJwCvAbYAbgSOLSqLkuyHLgXcE/gP4BnT6Tgrc8bgHslOQf4A3ADsAa4D3DvJJ+mu9ZmR+C4dqElSS4ElgI7AV8Avg78Kd11OAdV1W9n+e/YKiR5GfCXbfS9wH7cvP5OBz4P7JTk48CfAGcDh1VVJXkocAzdOrsCeG5VXZrkTOAc4JHAR4C3zOKfNDuqyscceAD3A34CzG/juwK7cPMZYs8D3tKGl9O9gW8/6bq3pgewBDi3De8PXA/sOTJ91/bz9sC5wG5t/EK6r1BYAtwIPKi1n0K3EZr437alPYCHAj8E7ki3YT8PePDU+htZh9fQXSi7DfAtuo399sA3gQVtvr+gO0Ue4EzgnZP++4Z82FOYOx4LfKyqrgCoqquS3B84Ocnd6HoLvxiZ/9TyE+akfaeqRtfJXyd5WhteDOxN18Mb9YuqOqcNn00XFNr0Hgl8qqquB0jySeBR08z3naq6pM1zDt36uJqu53B6Euium7p05DknD1f25BkKc9vbgWOq6tQk+9P1EKZcP5GKNKpfB239PB54eFX9pu1m2HGa5/xuZPgmul6FJmfd9bEdEOC8qnr4DM/Zov/3PNA8d3wZeGaS3QCS7ArszM3f/3TEpApT71rgTjNM2xlY0wLhPnT7rzU5XwMOTnKHJHcEngZ8g5nX36gfAwuSPBwgyfZJ7jdcqXOLPYU5oqrOS/I64KtJbgK+T9cz+FiSNXShsecES9zqVdWVSb6R5Fzgt8BlI5O/CByd5AK6jcpZk6hRnar6XpIPAN9pTe+tqrNH1t8X6A40T/fc3yd5BvC2JDvTbSffSndcYovn11xIknruPpIk9QwFSVLPUJAk9QwFSVLPUJAk9TwlVdpE2jUmZ7TRu9JdDLW6je9bVb+fSGHSRvCUVGkA7UsLr6uqN0+6FmljuPtIktQzFCRJPUNBktQzFCRJPUNBktQzFCRJPU9JlST17ClIknqGgiSpZyhIknqGgiSpZyhIknqGgiSpZyhIknr/H8S+oMwfdx5SAAAAAElFTkSuQmCC\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    }
  ]
}