03 - SCMs.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "eeea2eaa",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.stats as stats\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "plt.style.use('fivethirtyeight')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8baca376",
   "metadata": {},
   "source": [
    "# SCMs\n",
    "\n",
    "This notbook contains solutions to some problems presented by Peters et al. in **Elements of Causal Inference** (2017)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0fd1c54f",
   "metadata": {},
   "source": [
    "## Joint sample from SCM\n",
    "\n",
    "Consider the SCM:\n",
    "\n",
    "$$\\large X := Y^2 + N_X$$\n",
    "$$\\large Y := N_Y$$\n",
    "\n",
    "where:\n",
    "\n",
    "$N_X, N_Y \\sim  N(0, 1)$ \n",
    "\n",
    "Generate an iid sample of 200 observation from joint distr. $(X, Y)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "96534b6e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Sample from Y\n",
    "y = stats.norm(0, 1).rvs(200)\n",
    "\n",
    "# Compute X\n",
    "x = y**2 + stats.norm(0, 1).rvs(200)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "31d3831c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ8AAAEJCAYAAABL3SrKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAAAaYUlEQVR4nO3dfVDVZf7/8RfCFzBScVj26CbYaCToLJlMYnZDsVu75a4aSabN5KIkCpnmLbqzOqBmire1LLsOMZnpKO3QrNmNTcZseIM1lqvl2h5Xc2wVWFCUm1iG4/n94cSvEyAHzjnXOQefj7861+f6fK73FcrL63M+5zoBtbW1dgEAYFAvbxcAALj5ED4AAOMIHwCAcYQPAMA4wgcAYBzhAwAwjvABABhH+AAAjCN83MRqtXq7BLdiPr6rJ81FYj6+zlPzIXwAAMYRPgAA4wgfAIBxhA8AwLggbxcAAD1RS0uLGhoavF2Gy0JDQ3XlypUOj4eFhSkoqOtRQvgAgJu1tLSorq5O4eHhCggI8HY5LgkJCVFoaGi7x+x2u2pra9WnT58uBxC33QDAzRoaGnpE8HQmICBA4eHh3VrhET4A4AE9PXi+1915ctutBwvdtNSpfk0vrvFwJQDgiJUPAMA4wgcAYBy33QDAgMkf1Rgdb/cvI7rU/9q1axo3bpz69eunXbt2tbY3Njbq/vvv14MPPqiNGze6rT5WPgAA9erVSwUFBSorK9P27dtb21etWqWWlhatXLnSreOx8gEASJJuv/12rVy5UsuWLVNSUpLOnj2rbdu2ae/evQoLC3PrWIQPAKDV9OnTtXfvXmVkZOj8+fPKyMjQvffe6/ZxuO0GAHCwceNGlZeXKzg4WEuWLPHIGIQPAMDBm2++qd69e+vChQs6d+6cR8YgfAAArT7//HNt3rxZ27Zt00MPPaS5c+fKZrO5fRzCBwAgSWpqatKsWbM0depUPfLII9qyZYu++eYbbdmyxe1jET4AAElSTk6OmpqatHr1akmSxWLRmjVr9PLLL+vkyZNuHYvwAQDo4MGD2rp1q/Lz89WnT5/W9okTJ+rxxx9XZmamWlpa3Dae049aFxYW6pVXXlFlZaViY2O1Zs0ajR07ttPz/v3vfyspKUl2u13/+c9/XCoWAPxVV3ccMO2+++5TTU37uzC8/vrrbh/PqZVPSUmJsrOztWDBAn3yyScaPXq0UlNTdf78+Rue19zcrOnTpzsVUgCAm4dT4ZOfn6+pU6dq2rRpGjZsmPLy8mSxWFRUVHTD81asWKERI0ZowoQJbikWANAzdBo+zc3NOnbsmJKTkx3ak5OTdeTIkQ7P27dvn/bt26e1a9e6XiUAoEfp9D2fmpoa2Ww2RUZGOrRHRkaqqqqq3XMqKio0d+5cbd++3eGNq85YrVan+/oiX6t/iJNfbXumg7p9bT6u6knz6UlzkXrefK5evaqQkBBvl+E2TU1NNzx+9erVdvMgJiamw3OcfuDgx1+VarfbO/z61JkzZ2r69Om65557nL28pBsX6uusVqvP1R/q5EaA7dXti/NxRU+aT0+ai9Qz59O3b1+FhoZ6uxS3aGpq6nQuffv2VVRUVJeu2+ltt4iICAUGBrZJterq6jaroe998sknWrt2rSIiIhQREaE5c+aooaFBERERHnlqAgDgXzpd+QQHB2vkyJEqLS3VxIkTW9tLS0s1fvz4ds85dOiQw+v33ntPGzZs0P79+/Wzn/3MtYoBAH7PqdtuWVlZysjIUEJCghITE1VUVKSKigqlpaVJuv6p2KNHj2rPnj2SpOHDhzuc/8UXX6hXr15t2gEANyenwiclJUWXLl1SXl6eKisrFRcXp+LiYkVHR0u6/oDB2bNnPVooAKDncPqBg/T0dKWnp7d7rKCg4IbnPvPMM3rmmWe6VhkAoMfim0wBwIDQTUuNjtf04pou9Z85c6ZOnjypjz/+WMHBwa3tZWVlmjJlit59910lJia6rT42FgUAaN26dbp8+bLDxgBXr17VvHnz9MILL7g1eCTCBwAgKTw8XPn5+dqyZYuOHj0qSVq2bJnCw8OVnZ3t9vEIHwCAJOmhhx7S9OnTNWvWLP3tb3/TW2+9pT/+8Y8Ot+HchfABALTKycmRJKWlpen3v/+94uLiPDIO4QMAaNW7d2/NmTNHISEhev755z02DuEDAHAQGBioXr16qVcvz0UE4QMAMI7wAQAYR/gAAIxjhwMAMKCrOw54k4kt0Vj5AACMI3wAAMYRPgAA4wgfAIBxhA8AwDjCBwA8wG63e7sEI7o7Tx61RrtfcjWkoUGhYWHdup4/PVIKeEJYWJhqa2sVHh6ugIAAb5fjMXa7XbW1terTp0+XzyV8AMDNgoKC1KdPH129etXbpbjs6tWr6tu3b4fH+/Tpo6CgrkcJ4QMAHhAUFKR+/fp5uwyXVVVVKSoqyu3X5T0fAIBxhA8AwDjCBwBgHOEDADCO8AEAGEf4AACMI3wAAMYRPgAA4wgfAIBxhA8AwDjCBwBgHOEDADCO8AEAGEf4AACMI3wAAMYRPgAA4wgfAIBxhA8AwDjCBwBgHOEDADDO6fApLCxUfHy8LBaLkpKSdOjQoQ77njp1Sr/5zW8UExMji8Wiu+66S7m5uWpubnZL0QAA/xbkTKeSkhJlZ2drw4YNGjNmjAoLC5Wamqry8nJFRUW16R8cHKwpU6YoPj5e/fr105dffqm5c+eqpaVFubm5bp8EAMC/OBU++fn5mjp1qqZNmyZJysvL0/79+1VUVKQVK1a06T9kyBANGTKk9XV0dLQOHDigw4cPu6lsAIA/6zR8mpubdezYMc2ZM8ehPTk5WUeOHHFqkDNnzmj//v167LHHulclHIRuWurtEgDAJZ2GT01NjWw2myIjIx3aIyMjVVVVdcNzH330Uf3jH//Q//73P02bNk3Lly+/YX+r1epEyb7LVP1DGhqMjNPQzXHO+OjP0d//fP1QT5qLxHx8XXfnExMT0+Exp267SVJAQIDDa7vd3qbtx4qKilRfX68vv/xSy5cv1+bNmzV//vxuFerrrFarsfpDw8I8PkZDQ4PCujmOL/4cTf58PK0nzUViPr7OU/PpNHwiIiIUGBjYZpVTXV3dZjX0Y4MGDZIkxcbGymaz6YUXXtALL7ygoCCnMw8A0AN1+qh1cHCwRo4cqdLSUof20tJSJSYmOj3QtWvX1NLSIpvN1vUqAQA9ilNLkKysLGVkZCghIUGJiYkqKipSRUWF0tLSJEk5OTk6evSo9uzZI0natWuXQkNDNXz4cAUHB+uLL75Qbm6uJkyYoJCQEM/NBgDgF5wKn5SUFF26dEl5eXmqrKxUXFyciouLFR0dLUmqqKjQ2bNn//9Fg4K0ceNGnTlzRna7XVFRUUpPT1dmZqZnZgEA8CtOv/mSnp6u9PT0do8VFBQ4vJ40aZImTZrkWmUAgB6Lvd0AAMYRPgAA4wgfAIBxhA8AwDjCBwBgHOEDADCO8AEAGEf4AACMI3wAAMYRPgAA4wgfAIBxfLEOXHb0v80Or3M/qnH5mrt/GeHyNQD4LlY+AADjCB8AgHGEDwDAOMIHAGAc4QMAMI7wAQAYR/gAAIzjcz5wu+Ufv+RUv9zkZR6uBICvYuUDADCOlc9N6Mc7ErTHZuulwMbO+wFAd7DyAQAYR/gAAIwjfAAAxhE+AADjCB8AgHGEDwDAOMIHAGAc4QMAMI7wAQAYR/gAAIwjfAAAxhE+AADjCB8AgHGEDwDAOMIHAGAc4QMAMI7wAQAYR/gAAIxzOnwKCwsVHx8vi8WipKQkHTp0qMO+ZWVlmjJlioYNG6aBAwdq7Nix2r59u1sKBgD4P6fCp6SkRNnZ2VqwYIE++eQTjR49WqmpqTp//ny7/T/99FONGDFC27Zt0+HDhzVjxgzNmzdPb731lluLBwD4pyBnOuXn52vq1KmaNm2aJCkvL0/79+9XUVGRVqxY0ab/ggULHF7PmDFDZWVl2rNnj1JTU91QNgDAn3W68mlubtaxY8eUnJzs0J6cnKwjR444PVBdXZ3Cw8O7XCAAoOfpdOVTU1Mjm82myMhIh/bIyEhVVVU5NcgHH3ygv//979q3b1/3qgQA9ChO3XaTpICAAIfXdru9TVt7ysvL9dxzz2nt2rVKSEi4YV+r1epsOT7JVP1DGhpcOt9mc+45E5vN5tI4nWmo73geVuslt4/n73++fqgnzUViPr6uu/OJiYnp8Fin4RMREaHAwMA2q5zq6uo2q6EfO3z4sJ566iktXbpUM2bMcKlQX2e1Wl2uP3TTUuc6hoW5NE5gY3OnfWw2mwIDA10apzNht3Y8j5iYaLeO5Y6fj6/oSXORmI+v89R8Ov0ncHBwsEaOHKnS0lKH9tLSUiUmJnZ43sGDB5WamqrFixcrMzPT9UoBAD2GU/dfsrKytHPnTr3xxhv6+uuvtWTJElVUVCgtLU2SlJOTo/Hjx7f2LysrU2pqqtLS0vTUU0+psrJSlZWVqq6u9swsAAB+xan3fFJSUnTp0iXl5eWpsrJScXFxKi4uVnT09VsjFRUVOnv2bGv/nTt3qrGxUa+++qpeffXV1vaoqCidOHHCzVMAAPgbpx84SE9PV3p6ervHCgoK2rz+cRvQFZM/qnHr9VYNduvlALiIvd0AAMYRPgAA4wgfAIBxhA8AwDjCBwBgHOEDADCO8AEAGEf4AACMI3wAAMYRPgAA4wgfAIBxhA8AwDjCBwBgHOEDADCO8AEAGEf4AACMI3wAAMYRPgAA4wgfAIBxhA8AwDjCBwBgHOEDADAuyNsF3CxCNy31dgkA4DNY+QAAjCN8AADGcdvNxx39b7O3SwAAt2PlAwAwjvABABhH+AAAjOM9H/i85R+/1Gmf3ORlBioB4C6sfAAAxhE+AADjCB8AgHGEDwDAOMIHAGAcT7vhpvDiVyEKO1fj1mvu/mWEW68H3EwIH/QInT2ObbPZFBgYyCPZgI/gthsAwDjCBwBgHOEDADCO8AEAGOd0+BQWFio+Pl4Wi0VJSUk6dOhQh32bmpo0e/ZsjR07Vj/5yU80btw4txQLAOgZnHraraSkRNnZ2dqwYYPGjBmjwsJCpaamqry8XFFRUW3622w2hYaGaubMmfrwww915coVtxcO/+fMhqEAeianVj75+fmaOnWqpk2bpmHDhikvL08Wi0VFRUXt9g8LC9OmTZv0u9/9TrfddptbCwYA+L9Ow6e5uVnHjh1TcnKyQ3tycrKOHDniscIAAD1Xp7fdampqZLPZFBkZ6dAeGRmpqqoqtxZjtVrdej3TblT/kIaGbl3TZvPeMyE2m81rY3uCzWZTQ333fg7tsVovue1aXR/bv/+u/Bjz8W3dnU9MTEyHx5ze4SAgIMDhtd1ub9PmqhsV6uusVusN6w8NC+vWdQMbm7tbkku+3xGgp/h+PmG3du/n0J6YmGi3XasrOvuz5m+Yj2/z1Hw6/Wd1RESEAgMD26xyqqur26yGAABwRqfhExwcrJEjR6q0tNShvbS0VImJiR4rDADQczl12y0rK0sZGRlKSEhQYmKiioqKVFFRobS0NElSTk6Ojh49qj179rSec+rUKTU3N6umpkYNDQ06fvy4JCk+Pt4D0wAA+BOnwiclJUWXLl1SXl6eKisrFRcXp+LiYkVHX7/nXVFRobNnzzqck5qaqvPnz7e+fvDBByVJtbW1biodAOCvnH7gID09Xenp6e0eKygoaNN24sSJ7lflpyZ/1PH3xSz/r3ceHAAAX8TebgAA4wgfAIBxfJMp0E03us3aHXwtN24mrHwAAMax8nFR6Kalkq5vn7O88f+8XA064+xO2rnJyzxcCXBzY+UDADCO8AEAGEf4AACMI3wAAMYRPgAA4wgfAIBxPGrdge8foQYAuB8rHwCAcYQPAMA4wgcAYBzhAwAwjvABABhH+AAAjONRa8BHOPv9QA31IQo751xfviMIvoqVDwDAOMIHAGAc4QMAMI7wAQAYR/gAAIwjfAAAxhE+AADjCB8AgHGEDwDAOHY4AHowZ3dNcBY7JsBdWPkAAIwjfAAAxnHbDYDT3H0bT5JWDXb7JeEHWPkAAIwjfAAAxhE+AADjCB8AgHGEDwDAOMIHAGBcj3nUOnTTUqf6Nb24xsOVAAA6w8oHAGCc0+FTWFio+Ph4WSwWJSUl6dChQzfs/9VXX+nxxx/XgAEDFBcXp7Vr18put7tcMADA/zl1262kpETZ2dnasGGDxowZo8LCQqWmpqq8vFxRUVFt+l+9elVPPPGExo4dq48//lhWq1VZWVm65ZZbNGfOHLdPAoD/evGrEIWdc//OCd7SUO/cfNy9Sasndp/w5EayTq188vPzNXXqVE2bNk3Dhg1TXl6eLBaLioqK2u3/1ltv6bvvvlNBQYGGDx+uCRMmaO7cufrTn/7E6gcAoIDa2tobpkFzc7MGDhyo1157TRMnTmxtX7hwoU6ePKn33nuvzTkZGRm6fPmyiouLW9s+//xzJScn69ixY7r99tvdNgEAgP/pdOVTU1Mjm82myMhIh/bIyEhVVVW1e05VVVW7/b8/BgC4uTn9wEFAQIDDa7vd3qats/7ttQMAbj6dhk9ERIQCAwPbrFiqq6vbrG6+99Of/rTd/pI6PAcAcPPoNHyCg4M1cuRIlZaWOrSXlpYqMTGx3XNGjx6tw4cPq6mpyaH/wIEDNXgwX94BADc7p267ZWVlaefOnXrjjTf09ddfa8mSJaqoqFBaWpokKScnR+PHj2/tP2nSJPXu3VuZmZk6efKk9uzZo82bNyszM5PbbgAA58InJSVFa9asUV5enh544AGVl5eruLhY0dHRkqSKigqdPXu2tX+/fv309ttv6+LFi3r44Ye1aNEiZWVl6fnnn/fMLHzE5cuXtWjRIt1zzz0aMGCARowYofnz5+vSpUveLs1pXf0wsa/auHGjHn74YUVFRWno0KGaPHmyTp486e2y3GbDhg0KDw/XokWLvF1Kt1VUVGjWrFkaOnSoLBaLEhMTdeDAAW+X1S02m02rVq1q/bsTHx+vVatWqaWlxdulOeXgwYN6+umnFRcXp/DwcO3YscPhuN1u15o1axQbG6sBAwZo3Lhx+uc//+nSmE7v7Zaenq709PR2jxUUFLRpGzFihN5///3uV+aHLl68qIsXLyonJ0exsbG6cOGCFi5cqBkzZujtt9/2dnmd6uqHiX3ZgQMHNGPGDI0aNUp2u10vvfSSJk6cqCNHjqh///7eLs8ln332mbZt26YRI0Z4u5Ruq62t1a9+9SuNGTNGxcXFioiI0Llz5/z2PeHNmzersLCw9bONX331lWbPnq3g4GAtXrzY2+V1qqGhQcOHD9eUKVM0a9asNse3bNmi/Px85efnKyYmRuvWrdMTTzyhzz77TH369OnWmJ1+zgeu+fDDDzV58mSdO3dOffv29XY5N/SLX/xCI0aM0CuvvNLaNmrUKE2YMEErVqzwYmWuq6+vV3R0tHbs2KHHHnvM2+V025UrV5SUlKQtW7Zo3bp1Gj58uPLy8rxdVpfl5ubq4MGD2rdvn7dLcYvJkyerf//++vOf/9zaNmvWLF2+fFm7d+/2YmVdd9ttt2ndunV65plnJF1f9cTGxuq5557TwoULJUnfffedYmJitHLlyta3X7qKjUU9rK6uTiEhIbrlllu8XcoNNTc369ixY0pOTnZoT05O1pEjR7xUlfvU19fr2rVrCg8P93YpLpk3b54mTJigpKQkb5fiknfffVcJCQlKS0vTHXfcofvvv19bt2712x1QxowZowMHDuhf//qXJOnUqVMqKyvTI4884uXKXHfu3DlVVlY6/G7o3bu3xo4d69Lvhh7zlQq+qLa2VqtXr9azzz6roCDf/l/dnQ8T+5Ps7Gz9/Oc/1+jRo71dSrdt27ZNZ86c0V/+8hdvl+Kyb775Rq+99poyMzM1b948nThxQkuWLJEkzZw508vVdd28efNUX1+vxMREBQYGqqWlRQsXLuzwrQp/UllZKantx2QiIyN18eLFbl/Xt38j+ohVq1Zp/fr1N+zzzjvv6IEHHmh93dDQoClTpmjgwIHKzc31dIlu09UPE/uDZcuWqby8XB988IECAwO9XU63WK1W5ebm6v3331dwcLC3y3HZtWvXdPfdd7fezr3rrrt05swZFRYW+mX4lJSUaNeuXSosLFRsbKxOnDih7OxsRUdH69lnn/V2eW7h7t8NhI8TZs+eraeeeuqGfQYNGtT63/X19UpNTZUk7d69W6GhoR6tzx2682Fif7B06VKVlJTonXfe8es9BT/99FPV1NTo3nvvbW2z2Ww6dOiQioqKdOHCBYWEhHixwq6xWCwaNmyYQ9udd96pb7/91ksVuWb58uV6/vnn9eSTT0q6/sDV+fPntWnTJr8PH4vFIun61mg//D3n6u8GwscJERERiohwbmvxuro6paamym63669//atuvfVWD1fnHj/8MPEPN5AtLS11+AyXP1myZIlKSkq0d+9e3Xnnnd4uxyXjxo3T3Xff7dCWlZWloUOHav78+X63GhozZoxOnz7t0Hb69Gm/e6rye42NjW1W1YGBgbp27ZqXKnKfwYMHy2KxqLS0VKNGjZIkNTU16fDhwy7d1SF83Kiurk4pKSmqq6vTjh071NjYqMbGRklS//79ff4XRFZWljIyMpSQkKDExEQVFRU5fJjYnyxcuFC7d+/Wm2++qfDw8Nb71mFhYX7zD4IfCg8Pb/OwxC233KL+/ftr+PDh3inKBZmZmXr00Ue1fv16paSk6Pjx49q6dav+8Ic/eLu0bvn1r3+tzZs3a/DgwYqNjdXx48eVn5+vp59+2tulOaW+vl5nzpyRdP2W6Lfffqvjx4+rf//+ioqK0uzZs7VhwwbFxMTojjvu0Pr16xUWFqZJkyZ1e0wetXajsrIy/fa3v2332I/fE/JVhYWF2rJliyorKxUXF6eXXnpJ9913n7fL6rKOnmpbsmSJli5darYYDxk3bpzfPmotSfv27VNubq5Onz6tQYMG6bnnnlNGRoZfvsdYV1en1atXa+/evaqurpbFYtGTTz6pxYsX+8Vt945+d02ZMkUFBQWy2+16+eWX9frrr6u2tlYJCQlav369S//wIXwAAMbxOR8AgHGEDwDAOMIHAGAc4QMAMI7wAQAYR/gAAIwjfAAAxhE+AADjCB8AgHH/D17kyNpXtUxtAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot histograms\n",
    "plt.hist(x, bins=15, density=True, label='X', alpha=.7)\n",
    "plt.hist(y, bins=15, density=True, label='Y', alpha=.7)\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "429214ec",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Joint sample\n",
    "xy = np.hstack([x[:, np.newaxis], y[:, np.newaxis]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "bbc1b3ec",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\aleks\\AppData\\Local\\Temp/ipykernel_7100/973994305.py:1: MatplotlibDeprecationWarning: Auto-removal of grids by pcolor() and pcolormesh() is deprecated since 3.5 and will be removed two minor releases later; please call grid(False) first.\n",
      "  plt.hist2d(x, y, density=True, bins=15)\n",
      "C:\\Users\\aleks\\AppData\\Local\\Temp/ipykernel_7100/973994305.py:4: MatplotlibDeprecationWarning: Auto-removal of grids by pcolor() and pcolormesh() is deprecated since 3.5 and will be removed two minor releases later; please call grid(False) first.\n",
      "  plt.colorbar()\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbMAAAEfCAYAAADLMygSAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAAAr9klEQVR4nO3df1hUZd4/8Pc0hqCtDI0wWIr0sISgEYbOKAYkXLk+uo8kibb2QzET+dFjT5JC3y3DNRVJdFeBaMfZspWKDAPEanuULlAE1M0wXDF/rvsowwUKisqiw3z/2HVyHJzhx+CZc3i/rutcl3Of+9znc67KT/d97vs+submZiOIiIhE7D6hAyAiIuotJjMiIhI9JjMiIhI9JjMiIhI9JjMiIhI9JjMiIhI9JjMiIhI9JjMiIhK9AUIH0JdmuS0UOgTqp1pemmjX9pT76u3anr3dPHnaru0N8HnEru3ZOz572n5J2yftOl8ba/V826Dv++S+QpF0MiMi6q860L82d2IyIyKSoA50CB3CPcVkRkQkQTeM1pOZ7B7Fca8wmRERSZDBxjCj1P7yl9rzEBER+M6MiIgkwGBkMiMiIpHrX9M/mMyIiCTJ1jszqWEyIyKSIEP/ymVMZkREUsRhRiIiEj2D5FaSWcdkRkQkQR0cZiQiIrFjz8xBZWZmori4GCdOnICTkxPGjRuHFStWICAgQOjQJMneu5bXJXnatT3vne12bW/gmYt2bc916367tnfTrq05Pkfe5V4sbhj71xe+RPO0e/fuxcsvv4xvvvkGRUVFGDBgAJ555hlcunRJ6NCIiByOATKrh9SIJpkVFBTghRdeQEBAAEaPHo3c3Fw0NjaisrJS6NCIiByOAfdZPWzRarUIDAyESqVCeHg4Kioq7lq3ra0N8fHxCAkJwdChQzF9+nSLOuXl5VAoFBbH8ePHzeoVFhZCo9HAw8MDGo0GxcXFXXpe0SSzO7W2tqKjowMKhULoUIiIHE6HUWb1sKagoAApKSlYunQpysrKoFarERMTg3PnznVa32AwwNnZGYsWLcKUKVOstl1ZWYm6ujrT4ePjYzpXXV2NBQsWICYmBuXl5YiJicH8+fNx8OBBm88rmndmd0pJScFjjz0GtVotdChERA6nN0OJWVlZmDt3LubNmwcAyMjIwO7du6HT6bBixQqL+oMHD8aGDRsAALW1tWhpablr2+7u7lAqlZ2ey8nJQWhoKJKTkwEAfn5+KC8vR05ODrZs2WI1ZlH2zN58801UVlbi448/hlwuFzocIiKHYzDeZ/W4m/b2dhw+fBgRERFm5REREaiqqup1XE899RT8/PwwY8YMlJWVmZ07cOCAxX0jIyO7dF/R9cxSU1NRUFCA4uJieHt7Cx0OEZFD6uhhX6WpqQkGgwHu7u5m5e7u7mhoaOhxPJ6ensjMzMQTTzyB9vZ2fPbZZ4iKisLOnTsxadIkAIBer+/xfUWVzJYvX46CggLs3LkTjz76qNDhEBE5rN7OWJTJzK83Go0WZd3h6+sLX19f02+1Wo2///3v2LRpkymZ9ea+ohlmTE5ORl5eHrRaLRQKBfR6PfR6PVpbW4UOjYjI4fR0mFGpVEIul1v0hhobGy16Tb0VHByMU6dOmX6rVKoe31c0yUyr1eLKlSuIioqCn5+f6di0aZPQoREROZwOyKwed+Pk5ISgoCCUlpaalZeWlkKj0dg1xiNHjkClUpl+jx8/vsf3Fc0wY3Nzs9AhEBGJRrux53+9JyYmIi4uDsHBwdBoNNDpdKivr0dsbCwAIC0tDYcOHUJRUZHpmmPHjqG9vR1NTU24evUqampqAACBgYEAgOzsbHh5ecHf3x/t7e3Iz89HSUkJtm7dampj8eLFmDZtGjIzM/HrX/8aO3fuRHl5Ob7++mubMYsmmRERUdf1dAIIAERHR+PixYvIyMiAXq+Hv78/8vPz4eXlBQCor6/H6dPmW47duQ4tLCwMwM8dkRs3buCtt97ChQsX4OzsbGrz9nVptxLnqlWrsGbNGjzyyCPQ6XQYN26czZhlzc3Nkt1beZbbQqFDEC3uzdg73FuQumr7JW2ftLu78Smr5yOHftcn9xUKe2ZERBLUlS2rpITJjIhIgjr62a75TGYCsPcQXl9ommTfYUG/zfV2bc/e8Q08Y9fmiATHnhkREYmewcZmwlLDZEZEJEG9mc0oRkxmREQSZG2XDyliMiMikiBru3xIEZMZEZEE9WYHEDHqX09LRNRP2PqatNQwmRERSRCn5hMRkehx0TQREYlebz/OKTZMZkREEsSeGRERiR57ZkREJHrsmRERkehxBxAiIhI97gBCRESid6NDLnQI95Skk5m9vht28+Rpu7TTV+21vDTRru2JgevW/XZt76ZdWyMSHhdNExGR6HE7KyIiEj1+z4yIiESPX5omIiLR4zAjERGJXn9bNN2/npaIqJ8wQGb1sEWr1SIwMBAqlQrh4eGoqKi4a922tjbEx8cjJCQEQ4cOxfTp0y3qFBUVYebMmfDx8cHw4cMRGRmJXbt2mdXZtm0bFAqFxdHW1mYzXiYzIiIJ6jDKrB7WFBQUICUlBUuXLkVZWRnUajViYmJw7ty5TusbDAY4Oztj0aJFmDJlSqd19u3bh7CwMOTn56OsrAxPP/00XnjhBYskOWjQINTV1Zkdzs7ONp+Xw4xERBLUm2HGrKwszJ07F/PmzQMAZGRkYPfu3dDpdFixYoVF/cGDB2PDhg0AgNraWrS0tFjUSU9PN/udkpKCv/zlLygpKUFISIipXCaTQaVSdTtmUfXM9u3bh+eeew7+/v5QKBTYtm2b0CERETmkDsisHnfT3t6Ow4cPIyIiwqw8IiICVVVVdo2xtbUVCoXCrOz69esYM2YMAgICMGfOHPzwww9daktUyezq1asICAjA2rVr4eLiInQ4REQO60aH3OpxN01NTTAYDHB3dzcrd3d3R0NDg93i++Mf/4jz589jzpw5pjJfX19s3rwZeXl50Gq1GDhwIKZOnYqTJ0/abE9Uw4xTpkwxjccmJCQIHA0RkePq7dR8mcz8eqPRaFHWU4WFhXj77bexZcsWeHl5mcrVajXUarXpt0ajQWhoKHJzc7Fu3TqrbYoqmRERUdf0dNd8pVIJuVxu0QtrbGy06K31RGFhIRYvXoz3338f06ZNs1pXLpcjKCgIp06dstmuqIYZiYioa3o6m9HJyQlBQUEoLS01Ky8tLYVGo+lVTDt27EBcXByys7MRFRVls77RaERtbW2XJoSwZ0ZEJEG9mc2YmJiIuLg4BAcHQ6PRQKfTob6+HrGxsQCAtLQ0HDp0CEVFRaZrjh07hvb2djQ1NeHq1auoqakBAAQGBgIAvvjiC8TFxeF3v/sdQkJCoNfrAfwrebq5uQEA1q5di/Hjx8PHxweXL19Gbm4uamtrkZmZaTNmJrMuMEQG27U9/biBdm1vU9z7dm0PAFa9Emv3Nono3unNO7Po6GhcvHgRGRkZ0Ov18Pf3R35+vun9Vn19PU6fNv+U1Z3r0MLCwgAAzc3NAACdToebN28iNTUVqamppnqTJk1CSUkJAKClpQVLlixBQ0MDhgwZgsDAQOzatQvBwbb/DmYyIyKSoN5+aXrhwoVYuHBhp+dycnIsyo4cOWK1vVsJy5o1a9ZgzZo1XQvwDqJKZq2traYXgR0dHfjHP/6BmpoauLm5YcSIEQJHR0TkOPrbRsOimgDy/fffIywsDGFhYbh+/TrWrFmDsLAwrF69WujQiIgcSm+2sxIjUfXMQkNDTeOvRER0d1JMWNaIKpkREVHX3Oxnn4BhMiMikiD2zIiISPSYzIiISPSYzIiISPSYzIiISPSMTGZERCR2vd0BRGyYzIiIJIjDjEREJHocZiQiItFjz4yIiETP0MEdQCSjaZKnXdpxnX/OdqVueMiurQGv5i62c4vA4IeNdm2v9WH7/LO4xfXkaduViPoxo33/E3Z4kk5mRET9FWczEhGR6HECCBERiR4ngBARkejxnRkREYkehxmJiEj0mMyIiEj0+M6MiIhEj+/MiIhI9PrbMGP/2u+EiKif6DDKrB62aLVaBAYGQqVSITw8HBUVFXet29bWhvj4eISEhGDo0KGYPn16p/X27t2L8PBwqFQqPP7449DpdBZ1CgsLodFo4OHhAY1Gg+Li4i49L5MZEZEEGY0yq4c1BQUFSElJwdKlS1FWVga1Wo2YmBicO9f51n4GgwHOzs5YtGgRpkyZ0mmdM2fOYPbs2VCr1SgrK8Prr7+OZcuWobCw0FSnuroaCxYsQExMDMrLyxETE4P58+fj4MGDNp+Xw4xERFLUi3dmWVlZmDt3LubNmwcAyMjIwO7du6HT6bBixQqL+oMHD8aGDRsAALW1tWhpabGo86c//Qmenp7IyMgAAPj5+eHgwYPYvHkzoqKiAAA5OTkIDQ1FcnKyqU55eTlycnKwZcsWqzGzZ0ZEJEE97Zm1t7fj8OHDiIiIMCuPiIhAVVVVj+Oprq62aDMyMhLff/89bty4AQA4cOBAp3W6cl8mMyIiCTIarR9309TUBIPBAHd3d7Nyd3d3NDQ09DiehoaGTtu8efMmmpqaAAB6vb7H95X0MOPVYfaZzdOkH2qXdvpM0LU+aHSQXVvz2n7Bru3dtGtrRNLT29mMMpn59Uaj0aLMHm3eWd7T+0o6mRER9Vs9TGZKpRJyudyiN9TY2GjRa+oODw+PTtscMGAAHnzwQQCASqXq8X05zEhEJEE9HWZ0cnJCUFAQSktLzcpLS0uh0Wh6HI9arcZ3331n0ebYsWNx//33AwDGjx/f4/symRERSZHRxmFFYmIi8vLysHXrVtTV1WH58uWor69HbGwsACAtLQ0zZswwu+bYsWOoqalBU1MTrl69ipqaGtTU1JjOx8bG4vz580hJSUFdXR22bt2KvLw8JCUlmeosXrwYZWVlyMzMxPHjx5GZmYny8nLEx8fbfFzRDTNqtVr84Q9/gF6vx6hRo7BmzRqEhIQIHRYRkUPpzTuz6OhoXLx4ERkZGdDr9fD390d+fj68vLwAAPX19Th9+rTZNXeuQwsLCwMANDc3AwC8vb2Rn5+PN998EzqdDp6enkhPTzdNywcAjUYDnU6HVatWYc2aNXjkkUeg0+kwbtw4mzHLmpubRbODV0FBARYtWoT169djwoQJ0Gq1yMvLQ2VlJUaMGGFRP2St5erynmjrkwkWjs35sINPADl52nYlIhHYfknbJ+0+XpRj9fwPM2z3dsREVMOMty/k8/PzQ0ZGBlQqVadbohAR9Wu9GGYUI9Eks75ayEdEJE0yG4e0iOadWV8t5CMikiQJ9r6sEU0yu6UvFvIREUkOk5lj6quFfEREksTvmTmmvlrIR0QkRT1dNC1WoumZAf9ayBcXF4fg4GDTeoTbF/IREdG/STBhWSOqZGZrIR8REf1bPxtmFFUyA4CFCxdi4cKFQodBROTQZOyZERGR6HWwZyYZ14d12KUdFztv7XRkifVtZrrruzb7/0v76uHFdm3P3ttPDfB5xK7tcXsskhz2zIiISPSYzIiISPT6WTKzus7s0qVL9yoOIiKyJ6PM+iExVpPZ5MmTcezYsXsVCxER2YnMaP2QGqvJ7L777sOUKVOwa9euexUPERHZAz8B87M9e/Zg3LhxePHFF/Hee+/dq5iIiIi6xWoyUygU2L59O+Li4vDuu+9iwYIFuH79+r2KjYiIeojDjHdWuO8+rF69GllZWdi1axemTp2K//u//7sXsRERUU/1swkgXZ6aP3fuXPj5+SEmJgahoaEIDw/HmDFjMHr0aIwZMwbDhw/vyziJiKg77LNnhGh0KZl1dHSguLgYubm5uHTpEgYPHoyvvvoKX375penDmEOGDMHo0aNRUlLSpwETEZFtUhxKtMZqMrt48SI+/PBD6HQ6nD9/Hg899BBWrlyJ+fPnY9CgQfjpp59QW1uLH3/8ET/++CNqa2vvVdxERGQNk9nPAgIC8M9//hOPPfYY3n77bTz77LOQy+Wm86NGjcKoUaPw7LPP9nmgRETUDUxmPwsNDUVSUhLCw8PvVTxERGQH/W2Y0epsxs8//5yJjIhIjHo5m1Gr1SIwMBAqlQrh4eGoqKiwWr+2thbTpk2Dp6cn/P39kZ6eDqPx54waHx8PhUJhcTz00EOmOtu2beu0Tltbm814Jb3RsM9rlXZpxxAZbJd2bnmrYYxd2/vyk1C7tgcAbUHX7NqevT/ZYm/8pAxJTi96ZgUFBUhJScH69esxYcIEaLVaxMTEoLKyEiNGjLCof/nyZcycORMhISHYs2cPfvrpJyQmJmLQoEF49dVXAQBr167FO++8Y3bdr371K4SEhJiVDRo0CN9//71ZmbOzs82Yba4zIyIi8enNoumsrCzMnTsX8+bNg5+fHzIyMqBSqaDT6Tqt//nnn+P69evIyclBQEAAoqKisGTJEmRnZ5t6Z66urlCpVKbj9OnTOHPmDObNm2cet0xmVk+lUnXpeZnMiIikqId7M7a3t+Pw4cOIiIgwK4+IiEBVVVWn11RXV2PixIlwcXExlUVGRuLChQs4e/Zsp9d89NFH8Pf3h0ajMSu/fv06xowZg4CAAMyZMwc//PCDrScFwGRGRCRJPe2ZNTU1wWAwwN3d3azc3d0dDQ0NnV7T0NDQaf1b5+7U0tKCwsJCvPTSS2blvr6+2Lx5M/Ly8qDVajFw4EBMnToVJ0+etPm8kn5nRkTUb/VyNuOtDTFMzRmNFmW26ndWDgD5+fkwGAx47rnnzMrVajXUarXpt0ajQWhoKHJzc7Fu3Tqr8TKZERFJkKyH21kplUrI5XKLHlVjY6NF7+sWDw+PTusD6PSajz76CDNmzICbm5vVWORyOYKCgnDq1CmbcXOYkYiITJycnBAUFITS0lKz8tLSUov3W7eo1Wrs37/fbAp9aWkphg0bhpEjR5rVPXjwIH788UeLIcbOGI1G1NbWdmkSCJMZEZEU9eLjnImJicjLy8PWrVtRV1eH5cuXo76+HrGxsQCAtLQ0zJgxw1R/1qxZcHFxQUJCAo4ePYqioiJs3LgRCQkJFsOMH330EXx8fPDkk09a3Hft2rXYvXs3zpw5g5qaGiQlJaG2thYLFiyw+bgcZiQikqDe7AASHR2NixcvIiMjA3q9Hv7+/sjPz4eXlxcAoL6+HqdP/7yW0tXVFTt27EBycjImT54MhUKBxMREJCUlmbV75coVFBQUYNmyZZ2+S2tpacGSJUvQ0NCAIUOGIDAwELt27UJwsO21vrLm5mbJbnoyy22hXdqx96Lpieur7dqeGBZNP/rbFru25+i4aJq6avslbZ+0q9nQ+ZqwW6r+x3ZvR0zYMyMikiLJdlM6x2RGRCRB/W2jYSYzIiIpYjIjIiKx6289M9FMzf/www/x61//Gl5eXlAoFHfd74uIiNCrqfliJJpkdu3aNURERCAlJUXoUIiIHJ6sw/ohNaIZZkxISAAAi+/c3AsDz1y0a3vfZlouFuyVYfZtDgCUX7nYrtQN//SW27U9+e5Ddm3P3ssv5JyaT0KTYO/LGtEkMyIi6gYmMyIiEjtOALmHVq1aBYVCYfUoLy8XMkQiInHqZxNABO2ZxcfHY/bs2VbrDB8+/B5FQ0QkHf2tZyZoMlMqlVAqlUKGQEQkTUxmjkmv10Ov1+PEiRMAgLq6OrS0tGDEiBE2P/BGRNTv9LNkJpp1ZjqdDmFhYXjllVcAALNnz0ZYWBh27dolcGRERI5HZuOQGtH0zFJTU5Gamip0GERE4tDPemaiSWZERNR1UtzlwxomMyIiKWLPjIiIxI5T84mISPyYzIiISOzYMyMiIvFjMiMiIrFjz4ws3LTzt6lc7dzeA3b+Fhdg/2+4NU3ytGt7V5eH2LU9r+0X7NreTbu2RtQDTGZERCR6TGZERCR2/W2YUTR7MxIRUdfJOoxWD1u0Wi0CAwOhUqkQHh6OiooKq/Vra2sxbdo0eHp6wt/fH+np6TAaf75PeXl5p9+sPH78uFk7hYWF0Gg08PDwgEajQXFxcZeel8mMiEiKevFxzoKCAqSkpGDp0qUoKyuDWq1GTEwMzp0712n9y5cvY+bMmfDw8MCePXuwdu1abNq0CZs3b7aoW1lZibq6OtPh4+NjOlddXY0FCxYgJiYG5eXliImJwfz583Hw4EGbj8tkRkQkQTKj9cOarKwszJ07F/PmzYOfnx8yMjKgUqmg0+k6rf/555/j+vXryMnJQUBAAKKiorBkyRJkZ2eb9c4AwN3dHSqVynTI5XLTuZycHISGhiI5ORl+fn5ITk7Gk08+iZycHJvPy2RGRCRFPeyZtbe34/Dhw4iIiDArj4iIQFVVVafXVFdXY+LEiXBxcTGVRUZG4sKFCzh79qxZ3aeeegp+fn6YMWMGysrKzM4dOHDA4r6RkZF3ve/tmMyIiCSopz2zpqYmGAwGuLu7m5W7u7ujoaGh02saGho6rX/rHAB4enoiMzMTH3/8MT7++GP4+voiKioK+/btM12j1+u7dd/bcTYjEZEU9XI2o0xm/glPo9FoUWar/u3lvr6+8PX1NZ1Xq9X4+9//jk2bNmHSpEk9vu8t7JkREUlQT3tmSqUScrncojfU2Nho0Wu6xcPDo9P6AO56DQAEBwfj1KlTpt8qlapb970dkxkRkRT18J2Zk5MTgoKCUFpaalZeWloKjUbT6TVqtRr79+9HW1ubWf1hw4Zh5MiRd73XkSNHoFKpTL/Hjx/frfvejsmMiEiCejObMTExEXl5edi6dSvq6uqwfPly1NfXIzY2FgCQlpaGGTNmmOrPmjULLi4uSEhIwNGjR1FUVISNGzciISHBNESYnZ2NnTt34uTJk/jb3/6GtLQ0lJSU4JVXXjG1s3jxYpSVlSEzMxPHjx9HZmYmysvLER8fb/N5+c6MiEiKjD1/aRYdHY2LFy8iIyMDer0e/v7+yM/Ph5eXFwCgvr4ep0//vMesq6srduzYgeTkZEyePBkKhQKJiYlISkoy1blx4wbeeustXLhwAc7OzqY2p0yZYqqj0Wig0+mwatUqrFmzBo888gh0Oh3GjRtnM2ZZc3OzZDc9meW2UOgQ7glDf9xoeJjtF8LdYfeNhu28mTRJ1/ZL2j5pd2q89Xa/zpHW34/smRERSZDMIHQE9xaTmQTIdx+ye5v2/oSJ0s7tPeD9oF3b+6ed25PbuWc2wOcRu7bHnmM/INkxt84xmRERSVB/2zWfyYyISIp6MQFEjJjMiIgkiD0zIiISPyYzIiISO/bMiIhI/PjOjIiIxK6/9cxEsTfjpUuX8MYbb2D8+PHw9PTE6NGj8frrr+PiRfvuUkFEJBk93GhYrESRzC5cuIALFy4gLS0NFRUVyM3NRUVFBV5++WWhQyMickgyg9HqITWiGGYMCAjAn//8Z9Pv//iP/8DKlSsxZ84cXL58GUOGDBEwOiIiByS9fGWVKJJZZ65cuYKBAwdi0KBBQodCRORw+ts7M1Ems+bmZrz77rt46aWXMGCAKB+BiKhv9bPZjIK+M1u1ahUUCoXVo7y83Oyaq1ev4je/+Q2GDRuGlStXChQ5EZFj683HOcVI0G5NfHw8Zs+ebbXO8OHDTX9ubW1FTEwMAOCzzz6Ds7Nzn8ZHRCRaEkxY1giazJRKJZTKrn0c5MqVK4iJiYHRaMT27dvxwAMP9HF0RETiJetnw4yieOF05coVREdH48qVK9i2bRuuXbuGa9euAQDc3Nzg5OQkcIRki72/n2Xv74XZm72//m3v7yzye2v9QIfQAdxbokhmhw8fxoEDBwAAwcHmf0kUFxcjNDRUiLCIiBwWe2YOKDQ0FM3NzUKHQUQkHv0rl4kjmRERUffIOvpXNmMyIyKSIg4zEhGR2Mk4AYSIiESvn/XMRLFrPhERdVMvPwGj1WoRGBgIlUqF8PBwVFRUWK1fW1uLadOmwdPTE/7+/khPT4fxtoRaVFSEmTNnwsfHB8OHD0dkZCR27dpl1sa2bds63Qmqra3NZrxMZkREEiQzGq0e1hQUFCAlJQVLly5FWVkZ1Go1YmJicO7cuU7rX758GTNnzoSHhwf27NmDtWvXYtOmTdi8ebOpzr59+xAWFob8/HyUlZXh6aefxgsvvGCRJAcNGoS6ujqzoyu7PXGYkYhIinoxzJiVlYW5c+di3rx5AICMjAzs3r0bOp0OK1assKj/+eef4/r168jJyYGLiwsCAgJw/PhxZGdnIykpCTKZDOnp6WbXpKSk4C9/+QtKSkoQEhJiKpfJZFCpVN2OmT0zIiIp6rBx3EV7ezsOHz6MiIgIs/KIiAhUVVV1ek11dTUmTpwIFxcXU1lkZCQuXLiAs2fP3vVera2tUCgUZmXXr1/HmDFjEBAQgDlz5uCHH36w8pA/YzIjIpKgng4zNjU1wWAwwN3d3azc3d0dDQ0NnV7T0NDQaf1b5zrzxz/+EefPn8ecOXNMZb6+vti8eTPy8vKg1WoxcOBATJ06FSdPnrT5vBxmJCKSol7OZpTJZHc0Z7Qos1W/s3IAKCwsxNtvv40tW7bAy8vLVK5Wq6FWq02/NRoNQkNDkZubi3Xr1lmNl8mMiEiKDD1LZkqlEnK53KJH1djYaNH7usXDw6PT+gAsriksLMTixYvx/vvvY9q0aVZjkcvlCAoKwqlTp2zGzWFGIiIJ6ukwo5OTE4KCglBaWmpWXlpaCo1G0+k1arUa+/fvN5tCX1paimHDhmHkyJGmsh07diAuLg7Z2dmIioqy+QxGoxG1tbVdmhDCnhlRH5DvPiR0CPcUP9nigHoxzJiYmIi4uDgEBwdDo9FAp9Ohvr4esbGxAIC0tDQcOnQIRUVFAIBZs2YhPT0dCQkJSE5OxokTJ7Bx40YsW7bMNMz4xRdfIC4uDr/73e8QEhICvV4P4F/J083NDQCwdu1ajB8/Hj4+Prh8+TJyc3NRW1uLzMxMmzEzmRERSVEvkll0dDQuXryIjIwM6PV6+Pv7Iz8/3/R+q76+HqdP//w/MK6urtixYweSk5MxefJkKBQKJCYmIikpyVRHp9Ph5s2bSE1NRWpqqql80qRJKCkpAQC0tLRgyZIlaGhowJAhQxAYGIhdu3ZZfPqrM7Lm5mbJ7nkyy22h0CEQEVm1/ZK2T9qd86T1CROf7V3WJ/cVCntmRERSxI2GiYhI7PilaSIiEj8mMyIiEj1+aZqIiESPPTMiIhI9JjMiIhI9Q/+azshkRkQkRUYmMyIiEjsOMxIRkehxNiMREYkee2ZERCR6TGZERCR6TGZERGTNAJ9HhA7Btg7OZiQiIrHrZz2z+4QOoKv++7//G0FBQfD09ISPjw9+85vfoK6uTuiwiIgck9Fo/ZAY0SSzsWPHIjs7G1VVVfjiiy9gNBrxzDPP4MaNG0KHRkTkcIwGg9VDakQzzBgbG2v688iRI/Hb3/4WTz75JM6cOQNfX18BIyMickBcZ+b4rl69im3btmH48OHw8vISOhwiIscjwaFEa0QzzAgAWq0WDz/8MB5++GH87//+L4qKijBw4EChwyIicjwdHdYPiRE0ma1atQoKhcLqUV5ebqofExODsrIylJSUwMfHB/PmzcO1a9cEfAIiIgfVzyaACDrMGB8fj9mzZ1utM3z4cNOfXV1d4erqCh8fH4wfPx7e3t4oKirCc88919ehEhGJilGCvS9rBE1mSqUSSqWyR9cajUYYjUa0t7fbOSoiIgmQYO/LGlFMADl16hSKiorw1FNPQalU4vz589iwYQOcnJzwq1/9SujwiIgcTz+bzSiKCSBOTk7Yu3cvZs2ahSeeeAILFizAL37xC3z77bdQqVRCh0dE5HiMHdYPG7RaLQIDA6FSqRAeHo6Kigqr9WtrazFt2jR4enrC398f6enpMN7RO9y7dy/Cw8OhUqnw+OOPQ6fTWbRTWFgIjUYDDw8PaDQaFBcXd+lxRdEzGz58OLZv3y50GEREomHsRc+soKAAKSkpWL9+PSZMmACtVouYmBhUVlZixIgRFvUvX76MmTNnIiQkBHv27MFPP/2ExMREDBo0CK+++ioA4MyZM5g9ezaef/55fPDBB6isrMTSpUuhVCoRFRUFAKiursaCBQuQmpqK//qv/0JxcTHmz5+Pb775BuPGjbMas6y5uVmyfdFZbguFDoGIJMieGw1/evD/2a2t2816cJHV89svfnDXc5GRkRg9ejT+8Ic/mMqeeOIJREVFYcWKFRb1t2zZgnfeeQfHjx+Hi4sLACAjIwM6nQ5Hjx6FTCbDihUrUFxcjL/+9a+m61599VUcO3YM3377LYB/bY5x6dIlfPnll6Y6UVFRGDp0KLZs2WL1eUTRM+up7Ze0QodARCQIa8nKmvb2dhw+fNjUo7olIiICVVVVnV5TXV2NiRMnmhIZ8K+E+O677+Ls2bPw9vZGdXU1IiIizK6LjIzEJ598ghs3buD+++/HgQMHsGjRIos6H3xg+1lE8c6MiIjujaamJhgMBri7u5uVu7u7o6GhodNrGhoaOq1/65y1Ojdv3kRTUxMAQK/Xd+u+t2MyIyIiCzKZzOy30Wi0KLNV/87yntaxdt9bmMyIiMhEqVRCLpdb9IYaGxstek23eHh4dFof+LmHdrc6AwYMwIMPPggAUKlU3brv7ZjMiIjIxMnJCUFBQSgtLTUrLy0thUaj6fQatVqN/fv3o62tzaz+sGHDMHLkSFOd7777zqLNsWPH4v777wcAjB8/vlv3vR2TGRERmUlMTEReXh62bt2Kuro6LF++HPX19aZPcaWlpWHGjBmm+rNmzYKLiwsSEhJw9OhRFBUVYePGjUhISDANEcbGxuL8+fNISUlBXV0dtm7diry8PCQlJZnaWbx4McrKypCZmYnjx48jMzMT5eXliI+Ptxkzk5kD6u5iRUeUmZmJyZMnY8SIEfDx8cGcOXNw9OhRocPqtfXr10OhUOCNN94QOpQeqa+vx+LFi+Hj4wOVSgWNRoO9e/cKHVa3GAwGrFq1yvTfSGBgIFatWoWbN28KHZpkREdHY82aNcjIyEBoaCgqKyuRn59v+uRWfX09Tp8+barv6uqKHTt24MKFC5g8eTLeeOMNJCYmmiUqb29v5Ofno6KiAqGhoXjvvfeQnp5uWmMGABqNBjqdDp988gkmTZqETz/9FDqdzuYaM0Di68zEqKCgAIsWLTJbrJiXl3fXxYqOKjo6GtHR0XjiiSdgNBqxevVqHDhwAFVVVXBzcxM6vB45cOAAXn75ZfziF79ASEgIMjIyhA6pW5qbmxEeHo4JEyZg0aJFUCqVOHv2LDw9PeHn5yd0eF22fv16bNq0CTk5OQgICEBtbS3i4+ORmJiIZcuWCR0eCYTJzMF0d7GiWLS2tsLLywvbtm3Df/7nfwodTre1tLQgPDwcv//977Fu3ToEBASILpmtXLkS+/btwzfffCN0KL0yZ84cuLm54f333zeVLV68GJcuXcJnn30mYGQkJA4zOpBbixXvXFhobbGiWLS2tqKjowMKhULoUHrktddeQ1RUFMLDw4UOpcdKSkoQHByM2NhY/PKXv8STTz6JDz74wGL/PEc3YcIE7N27F8ePHwcAHDt2DOXl5Xj66acFjoyEJOkdQMSmJ4sVxSIlJQWPPfYY1Gq10KF020cffYRTp04hNzdX6FB65cyZM9iyZQsSEhLw2muv4ciRI1i+fDkAWOy64Mhee+01tLa2QqPRQC6X4+bNm0hOTsbChdy+rj9jMnNAPV006KjefPNNVFZW4uuvv4ZcLhc6nG756aefsHLlSnz11VdwcnISOpxe6ejowNixY03D1Y8//jhOnToFrVYrqmRWUFCATz/9FFqtFqNGjcKRI0eQkpICLy8vvPTSS0KHRwJhMnMgPVms6OhSU1NRUFCA4uJieHt7Cx1Ot1VXV6OpqQkTJ040lRkMBlRUVECn0+H8+fMYOHCggBF2nUqlspjo8eijj+If//iHQBH1zNtvv42kpCQ8++yzAIDRo0fj3Llz2LBhA5NZP8Zk5kBuX6z4zDPPmMpLS0vN1nSIxfLly1FQUICdO3fi0UcfFTqcHpk+fTrGjh1rVpaYmAgfHx+8/vrrouqtTZgwASdOnDArO3HihKhmyQLAtWvXLHr4crkcHR22v9FF0sVk5mASExMRFxeH4OBg05qL2xcrikVycjI+++wz/PnPf4ZCoYBerwcADB48GA888IDA0XWdQqGwmLQyaNAguLm5ISAgQJigeighIQFTpkzBe++9h+joaNTU1OCDDz7AW2+9JXRo3TJ16lRs3LgRI0eOxKhRo1BTU4OsrCw899xzQodGAuLUfAek1Wrx+9//Hnq9Hv7+/li9ejUmTZokdFjdcrdZi8uXL0dqauq9DcbOpk+fLsqp+QDwzTffYOXKlThx4gSGDx+OV155BXFxcaJ6J3vlyhW8++672LlzJxobG6FSqfDss89i2bJlcHZ2Fjo8EgiTGRERiR7XmRERkegxmRERkegxmRERkegxmRERkegxmRERkegxmRERkegxmRERkegxmRHdxcaNG6FQKPDXv/7V4tw777wDhUKBP/3pTwJERkR3YjIjuov58+dj8ODByMrKMisvKCjAxo0bsWDBAtFtM0YkVUxmRHehUCjw/PPPo7Cw0LSzfE1NDZKSkjBx4kSkp6cLHCER3cJkRmRFQkICOjo68P7776OxsRHPP/88HnzwQXz88ce4//77hQ6PiP6Nu+YTWeHt7Y1p06Zh69atOHToEJqamvD1119j6NChQodGRLfhRsNENlRWVmLq1KkAAJ1Oh+joaIEjIqI7cZiRyIa6ujoAgIeHhyg/kkrUHzCZEVmxf/9+JCcn45e//CUaGhpQWFgodEhE1AkOMxLdxblz5xAREQF3d3d89dVXCA0NhVKpRGlpqdChEdEd2DMj6sS1a9cwd+5cGAwG5OXlwdXVFfHx8fj+++9RUVEhdHhEdAcmM6I7GI1GxMfH429/+xs+/PBDeHt7AwBefPFFuLq6YvPmzcIGSEQWmMyI7rBu3ToUFhZi9erVCAsLM5U/8MADmD9/Pr7++mucPHlSwAiJ6E5MZkS32blzJ9auXYsXX3wRixYtsjgfFxcHuVyO7OxsAaIjorvhBBAiIhI99syIiEj0mMyIiEj0mMyIiEj0mMyIiEj0mMyIiEj0mMyIiEj0mMyIiEj0mMyIiEj0mMyIiEj0mMyIiEj0/j9PvMZI4oxtsgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist2d(x, y, density=True, bins=15)\n",
    "plt.xlabel('$X$')\n",
    "plt.ylabel('$Y$')\n",
    "plt.colorbar()\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:econml-dowhy-py38]",
   "language": "python",
   "name": "conda-env-econml-dowhy-py38-py"
  },
  "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}