5_2_meta_learners.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "colab": {
      "name": "5_2_meta_learners.ipynb",
      "provenance": [],
      "collapsed_sections": []
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "aoxI3DOK9vm2",
        "colab_type": "text"
      },
      "source": [
        "# 5.2 Meta-learners（T-Learner、S-Learner、X-Learner）の実装\n",
        "\n",
        "本ファイルは、5.2節の実装です。\n",
        "\n",
        "第1章から取り上げている、人事研修の効果について因果推論を実施します。"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "2XdIDbdlejUk",
        "colab_type": "text"
      },
      "source": [
        "## プログラム実行前の設定など"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "XZFKJwcu-_Oj",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# 乱数のシードを設定\n",
        "import random\n",
        "import numpy as np\n",
        "\n",
        "np.random.seed(1234)\n",
        "random.seed(1234)\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "hx1idArc_F15",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# 使用するパッケージ（ライブラリと関数）を定義\n",
        "# 標準正規分布の生成用\n",
        "from numpy.random import *\n",
        "\n",
        "# グラフの描画用\n",
        "import matplotlib.pyplot as plt\n",
        "\n",
        "# SciPy 平均0、分散1に正規化（標準化）関数\n",
        "import scipy.stats\n",
        "\n",
        "# シグモイド関数をimport\n",
        "from scipy.special import expit\n",
        "\n",
        "# その他\n",
        "import pandas as pd\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "AWqP6yeQlI_t",
        "colab_type": "text"
      },
      "source": [
        "## データの作成"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "DpnGB2KZ_L8x",
        "colab_type": "code",
        "outputId": "dc238d6f-1c8d-4ed9-984a-a90a14c021e9",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 282
        }
      },
      "source": [
        "# データ数\n",
        "num_data = 500\n",
        "\n",
        "# 部下育成への熱心さ\n",
        "x = np.random.uniform(low=-1, high=1, size=num_data)  # -1から1の一様乱数\n",
        "\n",
        "# 上司が「上司向け：部下とのキャリア面談のポイント研修」に参加したかどうか\n",
        "e_z = randn(num_data)  # ノイズの生成\n",
        "z_prob = expit(-5.0*x+5*e_z)\n",
        "Z = np.array([])\n",
        "\n",
        "# 上司が「上司向け：部下とのキャリア面談のポイント研修」に参加したかどうか\n",
        "for i in range(num_data):\n",
        "    Z_i = np.random.choice(2, size=1, p=[1-z_prob[i], z_prob[i]])[0]\n",
        "    Z = np.append(Z, Z_i)\n",
        "\n",
        "# 介入効果の非線形性：部下育成の熱心さxの値に応じて段階的に変化\n",
        "t = np.zeros(num_data)\n",
        "for i in range(num_data):\n",
        "    if x[i] < 0:\n",
        "        t[i] = 0.5\n",
        "    elif x[i] >= 0 and x[i] < 0.5:\n",
        "        t[i] = 0.7\n",
        "    elif x[i] >= 0.5:\n",
        "        t[i] = 1.0\n",
        "\n",
        "e_y = randn(num_data)\n",
        "Y = 2.0 + t*Z + 0.3*x + 0.1*e_y \n",
        "\n",
        "# 介入効果を図で確認\n",
        "plt.scatter(x, t, label=\"treatment-effect\")\n"
      ],
      "execution_count": 3,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "<matplotlib.collections.PathCollection at 0x7f6b8ebd4160>"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 3
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAASc0lEQVR4nO3df5BdZX3H8fc3mwQjVkPIaiUJJDgRTRsFvYO0zBT8SWA6SfzZZMoIlppqxc60yhRGxjIMjra049BKa5GhqNggUqXrFCejAuOMQ5ClCBiYQIgjJFBZQZxpoZAf3/5xT+R0c3fv3ey5e7MP79fMzt7zPM8955vn3nz27Dnn7onMRJI0+80ZdAGSpGYY6JJUCANdkgphoEtSIQx0SSrE3EFtePHixbl8+fJBbV6SZqW77rrrF5k53KlvYIG+fPlyRkdHB7V5SZqVIuJnE/V5yEWSCmGgS1IhDHRJKoSBLkmFMNAlqRBdr3KJiGuA3weeyMzf7tAfwBXAWcAzwLmZ+Z9NFyrp8HTxTffxtTseYbp/52/unGDv/slXEkACEUx7e4O2cME8Lln7W6w/aUlj6+xlD/1aYM0k/WcCK6uvTcA/Tb8sSbPBxTfdx3Vbpx/mQNcwh3aYw+wPc4Cnn93DBd+4h5vu3t3YOrsGemb+AHhqkiHrgK9k21ZgYUS8uqkCJR2+Nt/x6KBLmNX27E8u37K9sfU1cQx9CVB/VXdVbQeJiE0RMRoRo2NjYw1sWtIg7SthV3nAHnv62cbWNaMnRTPzqsxsZWZreLjjJ1clzSJDEYMuYdY7ZuGCxtbVRKDvBpbVlpdWbZIKt/Ety7oP0oTmzQkuOOOExtbXRKCPAB+MtlOAX2Xm4w2sV9Jh7rL1qzn7lGNpYkd97pzuKzkwooRfDBYumMfl739jo1e5RLd7ikbEZuB0YDHwc+CvgHkAmfnF6rLFL9C+EuYZ4EOZ2fWvbrVarfSPc0nS1ETEXZnZ6tTX9Tr0zNzYpT+Bjx1ibZKkhvhJUUkqhIEuSYUw0CWpEAa6JBXCQJekQhjoklQIA12SCmGgS1IhDHRJKoSBLkmFMNAlqRAGuiQVwkCXpEIY6JJUCANdkgphoEtSIQx0SSqEgS5JhTDQJakQBrokFcJAl6RCGOiSVAgDXZIKYaBLUiEMdEkqhIEuSYUw0CWpEAa6JBXCQJekQhjoklQIA12SCtFToEfEmojYHhE7IuLCDv3HRcT3I+LeiLgtIpY2X6okaTJdAz0ihoArgTOBVcDGiFg1btjfAl/JzDcAlwKfbbpQSdLketlDPxnYkZk7M/N54Hpg3bgxq4Bbqse3duiXJPVZL4G+BHi0tryraqu7B3hP9fjdwG9ExNHjVxQRmyJiNCJGx8bGDqVeSdIEmjop+kngtIi4GzgN2A3sGz8oM6/KzFZmtoaHhxvatCQJYG4PY3YDy2rLS6u2X8vMx6j20CPiZcB7M/PppoqUJHXXyx76ncDKiFgREfOBDcBIfUBELI6IA+u6CLim2TIlSd10DfTM3AucD2wBHgBuyMxtEXFpRKythp0ObI+IB4FXAZ/pU72SpAlEZg5kw61WK0dHRweybUmarSLirsxsderzk6KSVAgDXZIKYaBLUiEMdEkqhIEuSYUw0CWpEAa6JBXCQJekQhjoklQIA12SCmGgS1IhDHRJKoSBLkmFMNAlqRAGuiQVwkCXpEIY6JJUCANdkgphoEtSIQx0SSqEgS5JhTDQJakQBrokFcJAl6RCGOiSVAgDXZIKYaBLUiEMdEkqhIEuSYUw0CWpED0FekSsiYjtEbEjIi7s0H9sRNwaEXdHxL0RcVbzpUqSJtM10CNiCLgSOBNYBWyMiFXjhl0M3JCZJwEbgH9sulBJ0uR62UM/GdiRmTsz83ngemDduDEJvLx6/ArgseZKlCT1opdAXwI8WlveVbXVXQKcHRG7gJuBj3daUURsiojRiBgdGxs7hHIlSRNp6qToRuDazFwKnAV8NSIOWndmXpWZrcxsDQ8PN7RpSRL0Fui7gWW15aVVW915wA0AmXk78BJgcRMFSpJ600ug3wmsjIgVETGf9knPkXFjHgHeDhARr6cd6B5TkaQZ1DXQM3MvcD6wBXiA9tUs2yLi0ohYWw37BPDhiLgH2Aycm5nZr6IlSQeb28ugzLyZ9snOetuna4/vB05ttjRJ0lT4SVFJKoSBLkmFMNAlqRAGuiQVwkCXpEIY6JJUCANdkgphoEtSIQx0SSqEgS5JhTDQJakQBrokFcJAl6RCGOiSVAgDXZIKYaBLUiEMdEkqhIEuSYUw0CWpEAa6JBXCQJekQhjoklQIA12SCmGgS1IhDHRJKoSBLkmFMNAlqRAGuiQVwkCXpEIY6JJUCANdkgoxt5dBEbEGuAIYAq7OzM+N6/888NZq8aXAKzNzYZOFSv128U33cd3WRwZdxqx06msW8bUP/86gy3jR67qHHhFDwJXAmcAqYGNErKqPycw/z8wTM/NE4B+Ab/ajWKlfDPPp+eHDT/GHX7p90GW86PVyyOVkYEdm7szM54HrgXWTjN8IbG6iOGmmbL7j0UGXMOv98OGnBl3Ci14vgb4EqL/bd1VtB4mI44AVwC0T9G+KiNGIGB0bG5tqrVLf7MscdAnStDV9UnQDcGNm7uvUmZlXZWYrM1vDw8MNb1o6dEMRgy5BmrZeAn03sKy2vLRq62QDHm7RLLTxLcu6D9KkTn3NokGX8KLXS6DfCayMiBURMZ92aI+MHxQRrwOOAjwzolnnsvWrOfuUYwddxqzlVS6Hh66XLWbm3og4H9hC+7LFazJzW0RcCoxm5oFw3wBcn+nBSM1Ol61fzWXrVw+6DOmQ9XQdembeDNw8ru3T45Yvaa4sSdJU+UlRSSqEgS5JhTDQJakQBrokFcJAl6RCGOiSVAgDXZIKYaBLUiEMdEkqhIEuSYUw0CWpEAa6JBXCQJekQhjoklQIA12SCmGgS1IhDHRJKoSBLkmFMNAlqRAGuiQVwkCXpEIY6JJUCANdkgphoEtSIQx0SSqEgS5JhTDQJakQBrokFcJAl6RCGOiSVIieAj0i1kTE9ojYEREXTjDmAxFxf0Rsi4h/bbZMSVI3c7sNiIgh4ErgncAu4M6IGMnM+2tjVgIXAadm5i8j4pX9KliS1Fkve+gnAzsyc2dmPg9cD6wbN+bDwJWZ+UuAzHyi2TIlSd30EuhLgEdry7uqtrrXAq+NiB9GxNaIWNNpRRGxKSJGI2J0bGzs0CqWJHXU1EnRucBK4HRgI/CliFg4flBmXpWZrcxsDQ8PN7RpSRL0Fui7gWW15aVVW90uYCQz92TmT4EHaQe8JGmG9BLodwIrI2JFRMwHNgAj48bcRHvvnIhYTPsQzM4G65QkddE10DNzL3A+sAV4ALghM7dFxKURsbYatgV4MiLuB24FLsjMJ/tVtCTpYJGZA9lwq9XK0dHRgWxbkmariLgrM1ud+vykqCQVwkCXpEIY6JJUCANdkgphoEtSIQx0SSqEgS5JhTDQJakQBrokFcJAl6RCGOiSVAgDXZIKYaBLUiEMdEkqhIEuSYUw0CWpEAa6JBXCQJekQhjoklQIA12SCmGgS1IhDHRJKoSBLkmFMNAlqRAGuiQVwkCXpEIY6JJUCANdkgphoEtSIQx0SSqEgS5JhZjby6CIWANcAQwBV2fm58b1nwtcDuyumr6QmVc3WCcAN929m8u3bGf30882vepJHTl/iMzkmT37Z3S707Vg3hze++alfP1Hj9BE6XOArL6atmThApYfvYDbdz7F/n5sADhi7hye29t9Io6cP8Rn3r2a9Sct6U8hUp90DfSIGAKuBN4J7ALujIiRzLx/3NCvZ+b5fagRaIf5Rd+8j2f37OvXJib0P8/P/Dab8Oye/Vy39ZHG1tfPH2e7n3627z+oewlzaL/en/jGPQCGumaVXg65nAzsyMydmfk8cD2wrr9lHezyLdsHEuZ6cdq3P7l8y/ZBlyFNSS+BvgR4tLa8q2ob770RcW9E3BgRyzqtKCI2RcRoRIyOjY1NqdDHZvgwi+R7TrNNUydFvw0sz8w3AN8FvtxpUGZelZmtzGwNDw9PaQPHLFww/SqlKfA9p9mml0DfDdT3uJfywslPADLzycx8rlq8GnhzM+W94IIzTmDBvKGmVyt1NDQnuOCMEwZdhjQlvQT6ncDKiFgREfOBDcBIfUBEvLq2uBZ4oLkS29aftITPvmc1Swaw13Tk/CFeOm/2XeG5YN4czj7lWJoqfQ4QzazqIEsWLuDU1yxiTr82QPsql14cOX+Iv3v/Gz0hqlmn61Uumbk3Is4HttC+bPGazNwWEZcCo5k5AvxZRKwF9gJPAef2o9j1Jy3xP9khuGz96kGXIGkGRGafLvrtotVq5ejo6EC2LUmzVUTclZmtTn2z7ziCJKkjA12SCmGgS1IhDHRJKsTATopGxBjws0N8+mLgFw2W0xTrmhrrmhrrmrrDtbbp1HVcZnb8ZObAAn06ImJ0orO8g2RdU2NdU2NdU3e41tavujzkIkmFMNAlqRCzNdCvGnQBE7CuqbGuqbGuqTtca+tLXbPyGLok6WCzdQ9dkjSOgS5JhTgsAz0i3h8R2yJif0RMeGlPRKyJiO0RsSMiLqy1r4iIO6r2r1d/9rep2hZFxHcj4qHq+1Edxrw1In5c+/rfiFhf9V0bET+t9Z04U3VV4/bVtj1Sa+/LnPU4XydGxO3Va35vRPxBra/R+ZroPVPrP6L69++o5mN5re+iqn17RJwxnToOoa6/iIj7q/n5fkQcV+vr+JrOUF3nRsRYbft/XOs7p3rdH4qIc2a4rs/XanowIp6u9fVzvq6JiCci4icT9EdE/H1V970R8aZa3/TnKzMPuy/g9cAJwG1Aa4IxQ8DDwPHAfOAeYFXVdwOwoXr8ReCjDdb2N8CF1eMLgb/uMn4R7T8p/NJq+VrgfX2Ys57qAv57gva+zFkvdQGvBVZWj48BHgcWNj1fk71namP+FPhi9XgD7ZufA6yqxh8BrKjWMzSDdb219h766IG6JntNZ6iuc4EvdHjuImBn9f2o6vFRM1XXuPEfp/1nv/s6X9W6fw94E/CTCfrPAr5D+9YCpwB3NDlfh+UeemY+kJnd7tDb8ebVERHA24Abq3FfBtY3WN46XrjFXi/rfh/wncx8psEaOplqXb/W5znrWldmPpiZD1WPHwOeAKZ2j8Le9HLD83q9NwJvr+ZnHXB9Zj6XmT8FdlTrm5G6MvPW2ntoK+07h/XbdG4Qfwbw3cx8KjN/SfvWlGsGVNdGYHND255UZv6A9g7cRNYBX8m2rcDCaN8gqJH5OiwDvUcT3bz6aODpzNw7rr0pr8rMx6vH/wW8qsv4DRz8ZvpM9evW5yPiiBmu6yXRvlH31gOHgejvnE1pviLiZNp7XQ/Xmpuar15ueP7rMdV8/Ir2/PR6s/R+1VV3Hu29vAM6vaYzWVenG8QfFvNVHZpaAdxSa+7XfPViotobma+udyzql4j4HvCbHbo+lZn/PtP11E1WW30hMzMiJrzus/rJu5r23Z4OuIh2sM2nfS3qXwKXzmBdx2Xm7og4HrglIu6jHVqHrOH5+ipwTmbur5oPeb5KFBFnAy3gtFrzQa9pZj7ceQ2N+zawOTOfi4g/of3bzdtmaNu92ADcmJn7am2DnK++GligZ+Y7prmKiW5e/STtX2PmVntYB93Uejq1RcTPI+LVmfl4FUBPTLKqDwDfysw9tXUf2Ft9LiL+BfjkTNaVmbur7zsj4jbgJODfmMacNVFXRLwc+A/aP9C31tZ9yPPVQdcbntfG7IqIucAraL+nenluP+siIt5B+4fkafnCTdknek2bCKiebhBfW7ya9jmTA889fdxzb2ugpp7qqtkAfKze0Mf56sVEtTcyX7P5kEvHm1dn+wzDrbSPXQOcAzS5xz9SrbOXdR907K4KtQPHrdcDHc+G96OuiDjqwCGLiFgMnArc3+c566Wu+cC3aB9bvHFcX5Pz1fWG5+PqfR9wSzU/I8CGaF8FswJYCfxoGrVMqa6IOAn4Z2BtZj5Ra+/4ms5gXRPdIH4L8K6qvqOAd/H/f1Pta11Vba+jfYLx9lpbP+erFyPAB6urXU4BflXttDQzX/062zudL+DdtI8hPQf8HNhStR8D3FwbdxbwIO2frp+qtR9P+z/bDuAbwBEN1nY08H3gIeB7wKKqvQVcXRu3nPZP3Tnjnn8LcB/tYLoOeNlM1QX8brXte6rv5/V7znqs62xgD/Dj2teJ/ZivTu8Z2odw1laPX1L9+3dU83F87bmfqp63HTiz4fd8t7q+V/1fODA/I91e0xmq67PAtmr7twKvqz33j6p53AF8aCbrqpYvAT437nn9nq/NtK/S2kM7w84DPgJ8pOoP4Mqq7vuoXcXXxHz50X9JKsRsPuQiSaox0CWpEAa6JBXCQJekQhjoklQIA12SCmGgS1Ih/g9Sp439On0bGgAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "BHcdUlW9koTa",
        "colab_type": "text"
      },
      "source": [
        "## データをまとめた表を作成し、可視化する"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "1EMwdGIIIPrK",
        "colab_type": "code",
        "outputId": "f402ab06-bf25-441a-e94d-46a46b6ae342",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 196
        }
      },
      "source": [
        "df = pd.DataFrame({'x': x,\n",
        "                   'Z': Z,\n",
        "                   't': t,\n",
        "                   'Y': Y,\n",
        "                   })\n",
        "\n",
        "df.head()  # 先頭を表示\n"
      ],
      "execution_count": 4,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/html": [
              "<div>\n",
              "<style scoped>\n",
              "    .dataframe tbody tr th:only-of-type {\n",
              "        vertical-align: middle;\n",
              "    }\n",
              "\n",
              "    .dataframe tbody tr th {\n",
              "        vertical-align: top;\n",
              "    }\n",
              "\n",
              "    .dataframe thead th {\n",
              "        text-align: right;\n",
              "    }\n",
              "</style>\n",
              "<table border=\"1\" class=\"dataframe\">\n",
              "  <thead>\n",
              "    <tr style=\"text-align: right;\">\n",
              "      <th></th>\n",
              "      <th>x</th>\n",
              "      <th>Z</th>\n",
              "      <th>t</th>\n",
              "      <th>Y</th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>0</th>\n",
              "      <td>-0.616961</td>\n",
              "      <td>1.0</td>\n",
              "      <td>0.5</td>\n",
              "      <td>2.303183</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>1</th>\n",
              "      <td>0.244218</td>\n",
              "      <td>1.0</td>\n",
              "      <td>0.7</td>\n",
              "      <td>2.668873</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>2</th>\n",
              "      <td>-0.124545</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.5</td>\n",
              "      <td>2.193123</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>3</th>\n",
              "      <td>0.570717</td>\n",
              "      <td>0.0</td>\n",
              "      <td>1.0</td>\n",
              "      <td>2.245229</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>4</th>\n",
              "      <td>0.559952</td>\n",
              "      <td>0.0</td>\n",
              "      <td>1.0</td>\n",
              "      <td>2.139868</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "</div>"
            ],
            "text/plain": [
              "          x    Z    t         Y\n",
              "0 -0.616961  1.0  0.5  2.303183\n",
              "1  0.244218  1.0  0.7  2.668873\n",
              "2 -0.124545  0.0  0.5  2.193123\n",
              "3  0.570717  0.0  1.0  2.245229\n",
              "4  0.559952  0.0  1.0  2.139868"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 4
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "L6Tb2Hjk9vno",
        "colab_type": "code",
        "outputId": "bfad21f6-aa31-4c4a-c6d9-ca0159033148",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 282
        }
      },
      "source": [
        "plt.scatter(x, Y)\n"
      ],
      "execution_count": 5,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "<matplotlib.collections.PathCollection at 0x7f6b8be691d0>"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 5
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO2df5Qc1XXnv3d6WlKPyGqEURwYEJK9jghYQRO0QFY5a0ubSNgEmAAOkODgxD6svU42EKITYbNGUsiiRBuz3nV2bR3bazsmWBjBRBiIkCN5HbMWtsYjWR6QbGzMj4Y1Y9DISBqknpm7f1TV6HX1e69e/equ7rmfc+bMTNXrqlevq+57dX8SM0MQBEHoXLpa3QFBEAQhX0TQC4IgdDgi6AVBEDocEfSCIAgdjgh6QRCEDqe71R3QccYZZ/CiRYta3Q1BEIS2YWho6GfMvEC3r5CCftGiRdi7d2+ruyEIgtA2ENFzpn2iuhEEQehwRNALgiB0OCLoBUEQOhwR9IIgCB2OCHpBEIQOp5BeN4IgCO3G4HAVm3ccwktj4zirt4K1a5ZgoL+v1d0CIIJeEAQhNYPDVdz+4AGM1yYBANWxcdz+4AEAKISwj1TdENEcIvo2Ee0nohEi2qBp8z4iGiWiff7PB5R9NxHRD/2fm7K+AEEQhFazecehaSEfMF6bxOYdh1rUo3pcVvQnAKxi5qNEVAbwTSJ6jJn3hNptZeY/VjcQ0ekA7gSwHAADGCKi7cx8OIvOC4IgFIGXxsZjbW82kSt69jjq/1v2f1yrlawBsJOZX/OF+04AlyXqqSAIQkE5q7cSa3uzcfK6IaISEe0D8Ao8wf2kptk1RPQ9InqAiM7xt/UBeEFp86K/TXeOm4loLxHtHR0djXEJgiAIrWXtmiWolEt12yrlEtauWdKiHtXjJOiZeZKZlwE4G8DFRPT2UJOHASxi5l+Ft2r/QtyOMPMWZl7OzMsXLNDm5REEQSgkA/19uPvqpejrrYAA9PVWcPfVSwthiAViet0w8xgR7Yanfvm+sv1VpdlnAPyN/3cVwDuVfWcD+HqSjgqCIBSZgf6+wgj2MC5eNwuIqNf/uwLgtwAcDLU5U/n3SgBP+3/vALCaiOYT0XwAq/1tgiAIQpNwWdGfCeALRFSCNzHcz8xfJaKNAPYy83YA/4mIrgQwAeA1AO8DAGZ+jYj+EsB3/GNtZObXsr4IQRAEwQwxuzrQNI/ly5ez5KMXBEFwh4iGmHm5bp/kuhEEQehwRNALgiB0OJLrRhCEGU2Rk5FlhQh6QRBmLEVPRpYVoroRBGHGUvRkZFkhgl4QhBlL0ZORZYUIekEQZixFT0aWFSLoBUGYsRQ9GVlWiDFWEIQZS2BwFa8bQRCEDsaUjCxPt8tmu3SKoBcEQQiRp9tlK1w6RUcvCIIQIk+3y1a4dIqgFwRBCJGn22UrXDpF0AuCIITI0+2yFS6dIugFQRBCZO12OThcxYpNu7B43SM4dmIC5RJldmwXxBgrCIIQIku3y7DxdWy8hnIXYX5PGWPHa+J1IwiCEJesXBezqgGrM77Wphg9s7ox/LHVqY/vQqSgJ6I5AL4BYLbf/gFmvjPU5s8AfABeKcFRAH/EzM/5+yYBHPCbPs/MV2bXfUEQhFOkdV3Mw7+9CPl0XFb0JwCsYuajRFQG8E0ieoyZ9yhthgEsZ+bjRPQhAH8D4Dp/3zgzL8u224IgCI3YXBejBHZek8RZvRVUNUI9ML42I3gqUtCzV1T2qP9v2f/hUJvdyr97ANyYVQcFQeg88hJuplVydWwcKzbtsp4vr0li7ZoldfsAz/i68rwF6N/4OA4fr9X1M4/gKSevGyIqEdE+AK8A2MnMT1qavx/AY8r/c4hoLxHtIaIByzlu9tvtHR0ddeq8IAjtRyAUq2PjYJwSboPD1dTHNrkokn8e2/nSqFiiJom7r16Kvt4KCEBfbwXXXNSHbUPVOiEf/lyWOAl6Zp701S9nA7iYiN6ua0dENwJYDmCzsvlcvzL57wH4b0T0VsM5tjDzcmZevmDBglgXIQhC+5BnZKjOLZIQUkEYzpfGvz1qkhjo78MT61bh2U2X44l1q7D74GjDGLgcLymx/OiZeQzAbgCXhfcR0W8C+CiAK5n5hPKZqv/7xwC+DqA/RX8FQWhz8jRO6lbPYSFvOl8a3/m4k0TUtWYdPBUp6IloARH1+n9XAPwWgIOhNv0APg1PyL+ibJ9PRLP9v88AsALAU9l1XxCEdiPvyNDw6rnP8Xy6SeLuq5c66crjThK2a80jeMrF6+ZMAF8gohK8ieF+Zv4qEW0EsJeZt8NT1ZwG4CtEBJxyo/wVAJ8moin/s5uYWQS9IMxgTMbJvCJD45wvqe983AArXZ8AoLdSxvorL8jc64Y8p5pisXz5ct67d2+ruyEIQk40Ox972vPl0d+sj0lEQ749tHGfCHpBEAQzYddJwHsjcFXrNAuboJcUCIIgzDjirKbT+NcXBRH0giDMKOJGwBYhhUFaJE2xIAgzirh+/K3IH581IugFQWgL1JzuKzbtShxJq8s7Y9uedW76ViCqG0EQCk+WBbVLRJjUOKGUiDSts81N3ypE0AuCUDjCxtLjJycyM4jqhLxtO5A+N32z3UnDiKAXBKFQ6FbvJpIYRPsMaYNNEbRpyfJtJCmioxcEoVDojKUmkhhEm61zzzOJmyuyohcEoVC4rtKTCudm69yL4J4pgl4QhEJhqsjUWylj7uzuTIRzVvVgXYiqMNUMRNALgpA7cYyRpiRkrsm+Wm34DNPsJG46RNALgpArcY2RaVQrRTB8himCe6YkNRMEIVdWbNpl9HJ5Yt2qQp2raG8DcZCkZoIgZEocgdhMY2SacxXxbSArxL1SEBzJKgS/3Ylb3LuZuWLSnKsIbpB54VJKcA4RfZuI9hPRCBFt0LSZTURbiegZInqSiBYp+273tx8iojXZdl8QmkNc4VYU8pic4grEZvqtpzlXEdwg88JFdXMCwCpmPkpEZQDfJKLHmHmP0ub9AA4z878mousB/DWA64jofADXA7gAwFkAvkZEv8zMbtEQglAQ2jEneV6qiLgCMQ9jpEl1lOZcRXCDzItIQc+etfao/2/Z/wlbcK8CsN7/+wEAnySveOxVAL7MzCcAPEtEzwC4GMC30nddEJpHO6728pqckgjELP3WoyawpOcqghtkXjjp6ImoRET7ALwCYCczPxlq0gfgBQBg5gkARwC8Sd3u86K/TXeOm4loLxHtHR0djXcVgpAz7ZiT3DQJVcfGU6lw0qhHslAl5aVLH+jvw91XL0VfbwUEz1OnaOUCk+LkdeOrWpYRUS+Ah4jo7cz8/Sw7wsxbAGwBPPfKLI8tCGlpx9WeaeUNILYKJ6wqueaiPuw+OBpLPZKVKinPt6tmRsw2k1julcw8RkS7AVwGQBX0VQDnAHiRiLoBzAPwqrI94Gx/myC0FUUIeomLbnIKcFHhBMK9OjYOwil9bXVsHNuGqrFXu1mpkjpZl54XkYKeiBYAqPlCvgLgt+AZW1W2A7gJnu79WgC7mJmJaDuAfyCij8Mzxr4NwLezvABBaBbtttoL+nrL1n3a/bYVcHj1HX7FTiKgs1qJt+PbVatxWdGfCeALRFSCp9O/n5m/SkQbAexl5u0APgvg731j62vwPG3AzCNEdD+ApwBMAPiweNwIQja4BC0N9PdNr8rD2FbALqmCX/J1/a5vOVmtxNvx7arVSAoEQWhDwituwFvVhtUpg8NVbHh4BIeP1+o+r2ursnjdIw2r+DC9lTJOTExF9iFun4VkSAoEQSgYppWw6wrZRd+tE6yAJ6CjMkHaDLmAJ6CJEEvnLivx1iGCXhCazOBwFWu/sh+1KW/NXB0bx9qv7Mfe517DtqGqk1eKi77bpH6ZO7s7Urjq9OCBQbbPF9C3JtD9t5udo1MQQS8ITWb99pFpIR9Qm2Lcu+d5rdHztvv349at++pWwKYV97xKGSs27cJLfqoGHS7GT5fVdxLdv9AaRNALQpMZG69pt5sE8ySfWvkHK3zdirvcRTh2csJ4/ABXQRy1+hbvl/ZBslcKQhuh6sDDUZynzelGbdJuQs1SEHdyJGmnISt6QWgy83vKDV4wADC7uwtdRE5ujUCjesUm4gnIxfgpOvf2QAS9IDSZO6+4ALd9ZT8mQ3r6qSnGey4+ezq1QBfRtNpGJVC9mLxqwhCAe65blmlSMfGcaS9E0AtCkxno78P67SMNuvTaFGP3wdHpkncmv/NA9eIS1AR4uv8g4VdaAd3JVZg6GdHRC0ILOGIwmKoeMQP9fbjmoj6UiAAAJSJcc9EpVUmc1AGBQLYVTnHJLNnJVZg6GRH0gtACXNIeDw5XsW2oOq2+mWTGtqHqtACO48ZY0uj+VQHtWkGrHfPyCyLoBSE3bCtkXU53ArDyvAXT/9tWz4PDVRw/OeHUj0q5pNX1A6cEtOlc67eP1G1rx7z8ggh6QXAmTtGMqBVyoJYh5TMM1K3YbYVDbn/wgNZzR0V1eeyLENCmc42N1yInKAA4dmKi8PVzZzIi6AXBgbjFwV102bsPjhrT/wLmVbJODROmr7eCZzddjifWrcJAf19kVaioTJYBge/8/J5yXZux8VpbFEufqYigFwQH4hohXXTZUW1MwtmkhlHbhIOiooKbbEFU4X4O9PdB1wUxyhYXca8UBAfiGiFdcq9HtTHlmzHlmAFOJRwzZY80uUAO9Pdp0xkDXv4clcHhqjHNghhli4ms6AXBgbhGSJcC2i5tBvr78MS6VbjnumUAgFu37sOxExMol6jucwTgxksXTqtqknDnFReg3EUN24+drNe/21btYpQtJpGCnojOIaLdRPQUEY0Q0Z9q2qwlon3+z/eJaJKITvf3/YSIDvj7pJqI0Ja4CGUVlzwwrrliwvaBsfFaQ06bsCE3CQP9fThtTuNLfm2S64S7bdUuCc2KiYvqZgLAbcz8XSL6BQBDRLSTmZ8KGjDzZgCbAYCIrgBwKzO/phxjJTP/LMuOC0IzSVI0wyUPjEsb1wjYJHVcw4wZPHlU4W5SOc3vKUt0bEGJFPTM/DKAl/2/XyeipwH0wasDq+MGAPdl1kNhRtAO+VNcE3hlfS1x9N4ubW39c7EtmNIT33nFBc79FJpLLB09ES0C0A/gScP+HgCXAdimbGYAjxPREBHdbDn2zUS0l4j2jo6OxumW0ObEdV0sMnlcSxy9d1RbXf9u2boPyzY87lW+crQbSHri9sK5ODgRnQbg/wD4K2Z+0NDmOgA3MvMVyrY+Zq4S0S8C2AngT5j5G7ZzSXHwmcWKTbu0q8i+3sp0gq+iEl4dHzuhL/yR5lpcs1QC0fVgTWMNnCrUDUhd13YkdXFwIirDW6XfaxLyPtcjpLZh5qr/+xUiegjAxQCsgl6YWbRD/hSdugNAQyZHE2muJRCyJvdHlSBwSf2ciq2PgY4/jeeOUEwiBT0REYDPAniamT9uaTcPwDsA3Khsmwugy9ftzwWwGsDG1L0WOgoXvXAr0aXmvXXrPmuhjzBhX/S4TAc2KUXFTZjqzA4OV6cLfJsIJqR2sJkI7ris6FcAeC+AA0QUlH3/CICFAMDMn/K3/Q6Ax5n5mPLZNwN4yJsr0A3gH5j5n7LouNA5FL32qM7rJY6QBwBqdE/XYhOwm3ccihTyAbo6s5t3HIrsNwNYtuFxHDs5Me3CKTnn2x8Xr5tvAoi8TZn58wA+H9r2YwAXJuybMENI4rrYTLJQIY0dr0WukqOKeiTtR6CScf28zsaQheum0DokBYJQCIpce9SkWorDvEo5sjKTLZ/OQH9fqn4Ek0ua6yiSzUSIh6RAEIQITKl5XamUSyBCZFK0JEnOXAneIHQ58OMcQ2hPRNALQgSq3zgQTzgG5f9cI051qEnObLnlbQRFSnQ58F2up0g2EyE+zn70zUT86IUiE9a1Hz85Een22EWAzo6q+tebioHrgpEWrXvEej6dd02lXMLs7i6tDj7cvtxFOG1ON8aO1wpnMxH0pPajF4SZjs2Q6hLQpBPywSpZPXZvTxmzu7twZNwuYEtE1rz0uj3jtUljHxnepFNEY7iQHhH0ghBBlDeMakyNMnaWiDDFbAy6Ony8hkq5hHuuW2YVtDdccg6+tOf51NcW0A5RyEJyREcvCBG4VJcK8sZH6c+nmOtK/MWtXBVw18BS3HjpQq1+vVIuNZT6C5jfU46VblnoDGRFLwgKqhplXqUMIhj17zovGV3wl0rY4Jom/cNdA0tx18BSp/QMQH2GyaLGLAj5IIJeEHzCKhpTubwAnZdMIDDXbx9p+Lxu5ZxF+gdbDIJJoItgn1mIoBcEH9cCH0B0dakgt4xupb1i067pbSvPW4BtQ9Vc0j8UOQhNaC4i6AXBxzXy01aAWyUsaHVG3W1DVVxzUR92HxwVVYqQGyLohY4hbcZFlxQBabxTTIbX3QdHxeNFyBXxuhE6giwqO0WlGAirVAaHq1ixaRcWr3sEKzbtijxXO+TdFzoTWdELuZJVXvOo40QlBHMhnEUz8LrRRYdG+dbrKHrefaFzkRQIbUgS4ZlU4KYR1LqI0SDU3lXPbTpOODXA4nWPaKNBCcCzmy536m8ckpQ/jJPiQBDiYkuBIKqbNiOJiiKpWiOtOsRWsCPOsVyCiqISgmVNEjWMFNUWWoWobtqMJCqKpGoN0+d0Zep0ROmeXVUrLkJVF6hULhGOnZjA4nWPRPY17ptLUjWMuDwKrSByRU9E5xDRbiJ6iohGiOhPNW3eSURHiGif//MxZd9lRHSIiJ4honVZX8BMI8lK0rSvOjZuNSaaPjfJ7LTCd1lNuxgiXVbr4dXy/J4ywF7QU9DXW7buQ//Gxxv6q3tzMbUN0BluJZWAUFRcVDcTAG5j5vMBXArgw0R0vqbdvzDzMv9nIwAQUQnA3wF4F4DzAdxg+KzgSBIVhWkfAVa1TK8hX4qKLS+LS6GMs3orkd4rrkI1yDfz7KbL0TOrW1tf9fDxWsN1bnh4RBsopWurnkvUMEK74FIz9mUAL/t/v05ETwPoA/CUw/EvBvCMXzsWRPRlAFc5flbQkKSQtu4zunzlqiplcLiKo29MOPXJtCoPZ3UMn7NSLmHleQsivVdMNWWB+ihTVd1i84cPX6ctl7xNVSVqGKFdiKWjJ6JFAPoBPKnZ/etEtB/ASwD+nJlH4E0ILyhtXgRwieHYNwO4GQAWLlwYp1uFI5xfnBmR+cVdSVJIW/cZkyAMhPbmHYe0K2IdtrcJVRjq9OCu9gOXKNNgggD0E5npOqMI8r67uFAKQhFxFvREdBqAbQBuYeafh3Z/F8C5zHyUiN4NYBDA2+J0hJm3ANgCeO6VcT5bJMICSF0tZiEokro7hgWlyT0wENquQTy6ICJbIq1wX2/duk97XNP5g+Pr+q6qkaJuoLjXGT6HCHqhnXByrySiMjwhfy8zPxjez8w/Z+aj/t+PAigT0RkAqgDOUZqe7W/rWKISY7nkGjeRRfRnQJTe27RKr5S7UCIvC3pQDzUcRKT279at+3DH4AHtsWzn0W1Xj2/ipbHxyDQGLtdpQyJZhXbDxeuGAHwWwNPM/HFDm1/y24GILvaP+yqA7wB4GxEtJqJZAK4HsD2rzmdF3FB2Gy5CIKmgSFqkQkeUMVE3EZS7CBNTPK3KmGTGtqHq9HiZ/Obv3fO8cUxXnregoXgG+dvDuGSX7O0pW4tdu1xn8PlgQgsjkaxCu+GiulkB4L0ADhBR8J79EQALAYCZPwXgWgAfIqIJAOMArmcv5HaCiP4YwA4AJQCf83X3LcGlQENa9YpLYqywoHBVx2SdK8VmTAwbUktEWp29OtGYrpsB3Hb//rrjAt51bxuqNqhZGMC2oSqWn3t6Xfuo66yUS2DWq20I0Jbns9k8TJGs4kLZGWSVnqMdmDEpEEwP7Zxyl9brImmWwqhC0eGQ98HhKtY+sB+1yVPfQ7lE2HzthQ03nUmv3lspY+7s7lxuWJfC14B3XS5t1Gs3XU9A+DuwtQ9SKty6dZ9RP/+TBKkQZpIwaHfifFe2dBRAe1bgkhQIMKs94pSJc0EXuNNb8dQJvZUy5pS7cOvWfdMqog0Pj9QJeQCoTTI2PNz44mNSpxw7OZGJ3l6HazEO1zaqmilqjMP7TX75vZXy9MNoUqtE1XIVmkeWqlL1mHHsVyZ5sH77SGZ2sCIxY1IgxBXccfWwttXE4HAVH33oQF1pueAGMgnIw8drDaH7OjXD8ZMTDZOVKlDTrkyi1FBxCb6HweEquoim9f06wt9B0PcND4/UXfPYeG1a3ZYkzsBEkgyVgp28xjRumg+TPNCVjwxiKdL2sZXMGEFv0p1Xyl0AKJVgiPLpDqtmAqJWweqKAjilU1dvtsXrHtF+NvB4CScRC47jel1R/uhxCSJhb3/wgFXIA8CxExMYHK42+NNv3nHIOLkFqp4sXr2zSH0s1JPXmMa1X7nY0lQmmTOb5O8YPID7nnwBk8woEeGGS87BXQNLUx0zihkj6NeuWYK1X9nfYFCcmGJc92/OTlXKbf32xhB6dVWtE/JxCEdyqkKst6dsVD/ZIl9d2LzjUKZCvlwirDxvAW67f3+kkAfqV+ouRtmk6jbT25gUCsmevMY0bpI505ufyWYH6N+UbTULdNwxeABf2vP89P+TzNP/5ynsZ4yOfqC/D6fNaZzXapM8Xcrt2U2X44l1q2IJ+cHhqvZ1D/Bu3qgbeH5PGeUum0PgqWPp9JBH35hAuRT9efU4ebR1YXKSsfU7LzgJ+QCd+6jN9z6urtbWvtmpj2cCeY1p3CRzJvfiO6+4wJqfKbg/gvtlbLyGw8drzvr8+558Idb2rJgxgh7wZl0dgRBNYiCy+bCf1VuJvIGZgesuPmf6hrP5butee2tTjNokGz/XcD7A+friPnw9ZfvtNIVkbzcuRtngoY4ba2BrLxkqsyevMU2SZE5Nghcs8ILjmJ4mgl3lGhXXYlrkxFn8JKHjBb0qwLsMwnBepZzY0m5b9a5dswRr1yyxrrjHxmvYNlTF2jVL8Oymy/G3v3uh8UGwnWuS2RoopOJ6fSYvH51AL3cR/svVv4reSnTGy7jojLKmh9qWklk3idtUCZKhMnvyHFOd4E56HFPmVhdxrLunAjlkI0svpDAdraMPG0l1s2alXAJR4yw9XpvEhodH6nS3K89b0KDLN+kG5/eU6260sKdI+FyB7lznWTK72xOsUQYkRmMyL5Mx1UVfHxVMZDJ4uvjdu2Ja7ZmCvWxjpDNuR+l2JUPlKbKKKSjSmJquyfT270JvT7kuq+rK8xZg21A1kfNFVnR0wJQpwKZEhCnm6S/WFmRjo1Iu4ZqL+hq+RFMd0MHhKm4xJPFSa5uagjl059LR11upu3FN15dXPVXbdYYplwjlLsLx2hQAzwtqTrnkbNzSndtlogmCsaSOqx01iZwuzXQ7j5Puuy+XCHNndRvtblGUSwQw6pw+kniuJQnYtAVMdfSK3vRaPsVcJ+BM2RCjGK9NYvfBUdx99VKnlU6UPl9tp3vDuO/JF3DDJedg98FRY39LRA39MF1fUgOY+vCXfF/4OMW+A5J8JorwW4jpAQvujSRpn4tIHhG8YUGY1ouraGhtXpOcWMjP91OShz9vugfJsj/r+JWOFvSuLldr1ywx+rpHEehyXW72KH1+VLsgiVgQpq1buepyp+cZRKQ73/rt0emMCEiUYsIF9fuISsccbt+ONDMIKUyUZ1aRU0hk7VXWM6s71jGDe1B3fxLQEEOSho41xg4OV3H8ZGOFJKOAS6jBCr6sOwYP4K23P4pF6x7BW29/VJua17SCDuvzbSttteLR7O4uzPezNeq8btQVV1YGMNvDH4SQu6yImuWiOBM8Z7LMaqqSpp4vkG1a7TzI+h4MJjMXyiWadtbQOVEw3IriuNKRK3qTnra3Usb6Ky9oEHBxqimpBALDJQhicLiKYyf0E8+dV1xQt80U3KUeH/BeESvlEu65bplRJ66qKJIIdp2axoaLkG+moE2iminyKlRHs4OQAqK+xzhRsHETkiVVHarnmVcpo1yi1AGNAUG/XWxEE5M8XZ4ySr2YBR0p6E3FnufO7o6V9yKKOb6boS0I4q6BpcaJp4vqV151fXP0lQxW0SaDj2mFEb7hddF9Ll5LcTE9kHkK1ziTXDvmt4kbFeqKrdawi2B1nYB0Y37L1n1Yv32kYWHmojrU9SlIIBjOkVTuIszvKWPMD3py4c2/MAs/ff1k3bZg0nO1EampSeI+u0noOEFvK/as3mCqYDGNdJS1/PDxmjVnS7DdpO4IFuzhm3TzjkOxVhmmVTQB2hVX+GHRJVuz9duELYR8fk8Zwx9brf1ckYRrO+a3ydIGo5LWUG2agIKgPdVZQHef6VJgRKkOdalCenvKOPrGhPYNuTbFOHy8hr4YuW9+drSGGy9dOO0UUSJqWLAF/e3f+Li1+HwwHjqPpizfejtOR+/i2RLWHZoWqj2zSpHpCWyCMNCbu7wxqDdKVq9sDL2gdC136NKP4BptIeQ69VRUf7LQMSehHfPbFDWwy5RWGjg1md8xaC8Nqd4Hg8PVSGGsSxVy+HgtUjUbN8HZtqEqVp63AJVyqeGtIrBBDA5XccTRHz94S8rr+4tc0RPROQC+CODNfn+2MPMnQm1+H8BfwJuYXgfwIWbe7+/7ib9tEsCEyc8zK1w8W1xXqsdOTsbKIxPmhku8crmumfKCvsfJrBdVPMV2HhsuOnmbr2+cVWCRhGteahAdWaqrotRTSc4VlZU1KphQXYWbirnfq9i2TATCe+0D+yPbmlKFZE3g7hx+PlRnCZC7j0fSQkeuuKhuJgDcxszfJaJfADBERDuZ+SmlzbMA3sHMh4noXQC2ALhE2b+SmX+WXbfNmB7U3sopzxZXAVIiN0NNn3+jh1OPLj/39Gn3PpegiSDsOsqgUyl34Y3a1PQDtfe513DvnucbXv1WnregLkJPLc7hMpHYhHxUwqg4AquZwjWKvNQgYZqprkp6LlNW1vXbR3BiYqrueKozgi619uJ1j2jvfxdB2EWEjzz4vchnMfiebnUM1ktLZN6aGCat6th4nTorayIFPTO/DOBl/+/XiehpAH0AnlLa/L0d/cIAAB3lSURBVF/lI3sAnJ1xP50xPajrrzylOrCl9g0guBkfVSOMmmbUFmxiEvpv+G2jVkIATdc/1dVdJQC/tnBeXRRtlF+9jWBlnyY4ysbK8xbUCQp1e7NpVgBVlC0gvAK3rZjTnkuHLSuri2dV+Phx87+rTDLjeM3+LKr3ZNIAyFaT52QfyxhLRIsA9AN40tLs/QAeU/5nAI8TEQP4NDNvMRz7ZgA3A8DChQvjdKsOlwf1DQcBpzOQBIRTKMTRg9tW9uO1qekgieBn2YbHGx4sVW+py+3OAPb8+LD2tVJXnCPwujFNflPMudZb3X1wVPt50/a8aUYAlU1dpVuBm1bMQPSklEQ1loV9JDj+HYMH8PKR/ARv+N40LWSIPHtcb6WMkxOT02k3dKh1mG3G3KzJy/DvLOiJ6DQA2wDcwsw/N7RZCU/Q/4ay+TeYuUpEvwhgJxEdZOZvhD/rTwBbAC/XTYxraCDqQR23fMF1fdJsi8rvofr4mo5p032rX7JtVRU87HHTntr86l2iSF2Joy4oko6+WcyrlLXfrauOOUi690ZtKnKMk6jGshj7s3orDTEmWaPLlhpcezh4L3gkot5Iyl2kdeuMcpnMijzeRpy8boioDE/I38vMDxra/CqAzwC4iplfDbYzc9X//QqAhwBcnLbTaUgTlUcArrnIPIncMXgAt27dF/lF2VRC6gNmSyUQuHTZ9uuwPdxZRpHG8aSZaQU+BoerOKaJ2i53UWQ6apXDx2tOYxz3ew3q+aYhOH7eBTV++8IztdsH+vswd3Z87/HeShmb33NhwzM+0H8qBbJr7Yek5HH8SEFPRATgswCeZuaPG9osBPAggPcy8w+U7XN9Ay6IaC6A1QC+n0XHw7gWDknzSsrwgqB05xgcrjYYRE309VYw35DvWhVutpVHlKH0hkvOiS20s3TTi7NKzytNQdJiMnljipM4bU73tLE8DeExjvO9utbzVenrreDGSxdqj593QQ2bei/OWwnBUwHtu3N15P0eeNPlRR5j5jLlrQDwXgAHiCgwZ38EwEIAYOZPAfgYgDcB+J/evDDtRvlmAA/527oB/AMz/1OmVwC9mmDtV/Zjw8MjDdGeaV9JTZF4rvVV1SCmNN4dJvVPiWj6IVt+7umxjYpZ6afjqAvyMIAWKQgrjOmNL7CRuBjLK+USZnd3GdU/YVy/1ySuidWxcew+ONqSVBG25zmOATgqZ0/43tz67efhqAGOjcktOg0uXjffRERAPjN/AMAHNNt/DODCxL1zxFRiL3hw1Ifc9uXPnVVCb88svDQ2ji6HvC6q4cR1AgkHMdmE23yLd5CpiEog5FudryWum2LWBtAiR7jaJmlAP/HpvG4A+2IhyT2QdCHUqonUluJDl1tKhymCPDiObsFQLnWhNhVP0pe6CJMRBt288kB1RAqEOJGna9csMSYAO35yEiMbPY8U1wIWcYOc1NlaFW7BQxkkOlq7ZgnuvOIC5/TJ6kq+CKvZZrkpmmiGgTfpZOpiQA9PfIPD1QY1hW2Mk94DadwgdROpKbVAX28Fx09ORLo52zAJRddnN+D3L11oHBPTgiEJJiHv4sWXlo4Q9K43Z3Vs3Kqjn1epLwF2zUV906so0wo/WFG4vm673JjBQ3n31Uux+doL6x5k03VOMVvzgbRiNdsMN0UTprHqIsLidY+kfqjSTKY24ace35SrRReQFCaun34wFqZEZl0EuKRf0hVyN711RAU2ke+LrL7BuBb4sQUbBl5383vKuPOKxmy2tuvJg3AhpDzoiFKCcUrXmSh3EUCoWz2HVSFRJeeSBrmY3Bp1YdEubU1RiHmVDiwiLqu6NKXw4nxnLn2Lutd0qOUQwwLQVj7ynuuWGUvoHRmvYY4iDOOi+p9HCeiohF8ugliH6f4PjhkUCDFlbVUxfc9EQHdXNimOs0p/YCsl2BGCHoA2sMgV22uk+iXkpfeOI5hdJpw0QqiTUL8v0xtZ0jFJMpmGV+nMwJHxRiFj+v5059IJbZf8R3lFjpoCAsNBSnNnlXDsZLQKxFZ/2bSogqEPSc5lm3TLXYTT5nRPTxSHj52wBmHp6slmWXfXJug7Jnvl+isbsya68sS6VcaVRZCDYvG6R6Z1/PdctwwAcOvWfZm47cXxI3dxlZsJVZVcUH2fpyICyOIS1/dfl1HxxMQU7rluGZ5Yt6ru+3Ptkym4arw2CWYY74E81REmATs2XqsTgi5CHtDHBegqV31pz/OnstEm67r2XMHzpvNtr00xemZ149lNl+OJdausQr6vt+KpYd9zYUuyjHaEjh5wyQ9jx5aaIDhe4LapqniyMHRm7aHSakNoEck6cVrc7yyO3cTF5hSl5z4yXsM91y3T3gPtlgsmPDHlmZ1SNwkO9PcZx9l10lTfGlvxHHaMoFchmEPMwwQh1K6rAF2+i7SGzjwEcxJDaBrVVKvdOaPIOitl3O8sbgCZLVeLSwKvs3orxnsgblI7wBsrF+eEPAhPxnm+kdje1KIWCr0GmaNL09BsOkbQ26ommQhyWmRB2puvlR4qQDovkiK4c0bRzMlUN+klCSDT5WoJl6xLMoGFx0KXtCusfw6PVVwXxqSEryVIz5DHJGNL7e0yzuuvvKCh1nOWMiYNHSPo477OhQuFu+SLt9HueVnSuGQWxZ0zimZMpqZJ75qL+urSRgPR6rnNOw4ZM5cG15F0AtP56asFtwP9c1AZLBzjEXxWlz3VlahnTo0NCfromp5BnahsrXt8D6PAoGtK7e06znNnd09/Z0m9hvKgYwS964qa4AVIqLnj4W/TZdkLR7OZ3DBdVQBFVXGkCTAqUvbJVo+vadLbfXDUi4vIofJWWAiptUtdCdrqinSr6ISf0ZWTgHlzyjgyXsM8jdfNyYkpa+pfNTYkuD6XxVy4XoLJi6m3Usa+O0/VMV6xaZd1wRL8qMGNgYMG0Bil/EZeORIS0DGC3lZZKuzXq3sAAsGvqxIVfjiBZCqAIqs40hgri1IhKsvxTTphmIRzdWx8ekUcFI2Jwhb0FdQtCPqaRu0W5YIaJiz8TBXObB4lLi6kcXXzpnO6FCOyHV/tp2msZ3d3FfqttmP86F38y1tNkf3b04xfUcY+q/HVXY8aUJQk+E2FAPzbt56On7w6bp1IbHpwdXyTXncaPXs4XiDuxGgLagL0949tbKOqnrn0zxgcBUxPzq4xDupnXYIUs3gTnREBU0DrX9ujKHrEapG9bnTHB+rfrEwPYByhNDhcddY763SwdwwecE5XrWILDDL1JxDkSe+ruEJLZX5PGcMfW63dl0aoAqdqMLsmcUuyoDDdTyY1VNRYm8giUtoVm6DvGNUN0HrPlSiKouJQyUpAZzH2pr7oXpd1umOTcU8d38Hhap1nxHRshE+cXOyHj9fqVCS6+r2uBEW3w2Po4sNtqoEcdV+lsaGYhshVjWRSp9x9tadCNeV+imvniNPHu69eavzuopIXzu8p11X7Cq7HxXbXDGeGjhL0AUVd2Wfty52WItkMbH1xNcLp6vw2uMBtH2kwANamGOu3j2Du7O7Yagz1gUwbyDM2XkP/xscb3hJsC4TB4SqOvqGpVlWiBrdEV5fPLgKiyqMeMbgvm4TW+u0jDec3CW2bUTQcRZwEm2A1JZyzJS+slEt13klZ2XWydGboOEFfJOEVpmgRq0Vyi7T1Jc4Nz/Bel03ja4qvGBuvJc6VFKTJyCLaNPyWANgXCJt3HNJ6rsyd1W011gYqCt3EGKyqbRG0prcF03eljq+6gtapNfIWfLbjm3IHBZNm1DOc5Llpxpt+pKAnonMAfBFetSgGsIWZPxFqQwA+AeDdAI4DeB8zf9ffdxOAO/ymdzHzFzLrvYYiCS8dRVIvFckt0taXODnS0xi208RS2FRHcVFzrgQCZV6ljDnlroYAJpNaZ2y8hhWbdk1PBuFngpXfQb/DBk2d2gywv4W6fle2ZzJvwWc7vstiLOtnuBlv+i4r+gkAtzHzd/36r0NEtJOZn1LavAvA2/yfSwD8LwCXENHpAO4EsBzevTRERNuZ+XBmVxCiSMKr6DTTZhClTrP1JU7I/rETE3Wuh2FMVbuC9AJpsH08cPN1nbCCVa8a6V0plxpcM22CNXwMW79NE2Tct9A435Xpmcxb8EUdv9mLsWa86buUEnwZwMv+368T0dMA+gCogv4qAF9kz4VnDxH1EtGZAN4JYCczvwYARLQTwGUA7svsCkIU0eBZVJplM3BRp9n6oj4IUSvnsfFG1YeKrmpXuRQvr7ipFKAJ1Wc7nNbABFFjJSPdKjhKsI7XJp36a5uA4gg+ndAypQA3PZN5C76iqVCDPuV5/ljulUS0CMA3ALydmX+ubP8qgE1+fVkQ0T8D+At4gn4OM9/lb//PAMaZ+b9qjn0zgJsBYOHChRc999xziS6oKD7d7ULSYilxcPXzdjWip80zrzuPSR9t0mG7ZoAM1CFAo2ugCVttUVONgqj+VMol67lLRPjR3e+O7FsS5JlsDpm4VxLRaQC2AbhFFfJZwcxbAGwBPD/6pMcp4mxdZNSVRF6GbBd1WhxPKbXPi9c9kujYuklAJ4zUjI3hfkUJbsKp9LQ6TxIdUQU5TEnQogKnbJMZYK5jmwXyTLYeJ0FPRGV4Qv5eZn5Q06QK4Bzl/7P9bVV4q3p1+9eTdDQORTJ4thM217g04xmlTkszwWR5bK+EntcunPROh9re1LcAVxtRVH4Um1otSv0VNRnkiTyTrSWywpTvUfNZAE8z88cNzbYD+APyuBTAEV+3vwPAaiKaT0TzAaz2twkFxOYal6aKVlTFK5unVNpjb3h4JPLYwWSg6pFPTJgFrq59mLCtw9VGZFtZz+8pT3vCBFXP1ApnA/1SfUzQ47KiXwHgvQAOEFHgy/URAAsBgJk/BeBReK6Vz8Bzr/xDf99rRPSXAL7jf25jYJgViofNgyPPwippPKVMxwbsdYTVY8d1yTUFRgVGzxJR3WQy0K/PZx4HgmdIjnpDiVo5ixplZuLidfNNePeZrQ0D+LBh3+cAfC5R74TMcNGBr12zpCG1QEDaYCCbAErjKWXKWRIlVF3UKnG3TzLXGT1t+cx1xT5MBKm1o6JG43jGiGCfWXRcZKzQiKueeqC/DxseHtH7mfvHyUNARLl5xsmBY0oZqztnQNyJxtQ+WMmr6PKZB7h4y4SDmCROREhCpI5eaE9UPe5t9+930oHb9PAMOOnMTfpjG4FuWa2tOafcNX282x88gOrYOBinhHkgJHXX5eKnvnnHoem+xdVbm9qb9OsmITzQ34cn1q0yGkIDF9FwgJQOiRMRbMiKPidamVgtvNJ1EUAuucmjVo1p3TNVA2iQ80Xn2RJMUmnUSUnKxAWE28+rlLUBTgFRQjhO4FrREuMJ7YEI+hxodWI11yyKqgBy+UyUwMqj7qypTy+NjRsjPrsImN1tDxAK9y2u3jpoHzVBugjhOBONGFOFJIigz4FWJ1Zz0deGBZBLmbYogZVH3VkTNg+hKUZDClxT27S6bdsEqat6ZHrTi5tmQAS7EAcR9CkwPbStMpgF/TH5cZSIMMUcK6kYEF2mLeoYujeB8NiZimcA5hzzJvVNn5+FUI34NVVpSqvbNn2nalRsQKvf9ISZS8cZY5MYA5Oex2QkbIXBTO2Pjkq5hL/93Qvx7KbLtcUbdAZGArDiracD8EqsRY2n7hjAqYySur4GY3f0jQmUS3ov3iCVLlAfBORiRA3OpRPyWei243zXaQLDBCENHbWib+aKyfTQ3rJ1H3or5YaMiEmEShyDblwVQpiB/j7sfe61unqnDOCJH52Kb4saz2Bb2EUzyCi597nXsPvgqHYyqk2xNZWvLpWui77aFtyURVKtOMbRNOqjolZNE9qDjhL0zdSN2x7OsfEayl2E+T3lhkIRrpjqpG54eKSh1JytPzoVgondB0cjC2dEjedAv1dSL6yGGa9NRhbNPjJew747VxsLMOuuMUpfbRqXKX+Fv2LTrlTC09U4Ojhcdappq0NUPkJaOkrQN1M3HlVJpzbF6JnVjeGPrTa2sa3STCtRXak5W3/iqItcxymqnWl/1CQS9DXLmgKmY82rlDMTni7GUZPthGBPVBZ8tshV04Ti01E6+qx04y56fpM+WsUmEG06/qjP6vS6SZJVha+zt6dsbKsSNZ5JBLLa1ywTb5mOZSvskQe2yS9KWEs0rJCWjhL0WQiIKAEcoGYKNNFFZJwsogxzUcIyvEp1yVwYdZ02g2iAy3iaDLsmwn2Ney02TMcaM3j45CU8Td+nS3pgiYYV0hKrwlSzWL58Oe/duzfRZ9MarVyrIYXPGRVVGq6oY9JDB+eKKpmXtiKQ6ToDg2jailO6ylXbhqqFqTKU5HtOQ5oqS1KhSXAhkwpT7ULaYJIkr8lhg5yuvF1Yp2rSHRNOrdZtU3DaikCm6wkMomnRfQ/Lzz29MJ4jzU4lkCaiVaJhhbR0nKBPS1JDoCrYXMrb6QSNbQUfJm1FoFYUUW92RKft7a4VwjPN9Us0rJCGjhf0cVU5Waz0XISoTtC4JunKYuXZ6cmxXFwSRXgKM4WOFvRJ/I+zWOm5CtGwoDHpjefOKuGN2tR09aJrLkovoDpdHSAuiYJwikhBT0SfA/DbAF5h5rdr9q8F8PvK8X4FwAK/jOBPALwOYBLAhMlQkBdJH/a0K72kQlQ3QZRLhJMTU9M6+UlmbBuqYvm5p2ci7DtV6IlLoiCcwmVF/3kAnwTwRd1OZt4MYDMAENEVAG4N1YVdycw/S9nPRLTyYbcJ0TsGD+C+J1+YXqHfcMk5uGtgad0EUfXT8KppFAJkZRqtkmuFDUIQikqkHz0zfwOAa0HvGwDcl6pHGeLif9ysJGgBdwwewJf2PF+3Qv/SnuexyD8/cMoP3eZZExWM1cxrajYusQ5ZBl0JQruTWcAUEfUAuAzANmUzA3iciIaI6OaIz99MRHuJaO/o6GgmfYp62F2Do7LkvidfMO4Lzr/h4ZHERUBacU3NxiULZJZBV4LQ7mRpjL0CwBMhtc1vMHOViH4RwE4iOui/ITTAzFsAbAG8gKksOhSlK2+FwS7K/91WVSnAtjKdCUZIV5VcJ9sgBCEOWQr66xFS2zBz1f/9ChE9BOBiAFpBnxe2h73ZOvwsVtVRKYdnghFS9O+CEI9MBD0RzQPwDgA3KtvmAuhi5tf9v1cD2JjF+ZKgM941W2C4JszqrZRxYmIqUch7q4RgM/Old3oMgCBkTaSOnojuA/AtAEuI6EUiej8RfZCIPqg0+x0AjzPzMWXbmwF8k4j2A/g2gEeY+Z+y7LwrJr31yvMWNNVg51rLdf2VFyTWL7fCCNlsu4Do3wUhHh2X1EyHLYFVUHu0GStRUz9stVyT0OxqRM1OECYIQiMzKqmZDpveupkGO5PKIevVaLONkDPBLiAI7UxH5aM3UZR83p2qcijK+AqCoGdGrOiLZLzrRJe/Io2vIAiNzAhB3+kJvFqNjK8gFJsZYYyNotnGS0EQhKyZ8cZYG0lSGQuCILQTM8IYa8Mlb4ogCEI7M+MFvbgGCoLQ6cx41U3eKQNE/y8IQquZ8Sv6tWuWoFyium3lEmXiGjgTUgYLglB8ZrygB+Blzbf9nxDR/wuCUARmvKDfvOMQalP1kr02xZkIY9H/C4JQBGa8oM9TGEtqAEEQisCMF/R5CmOpWyoIQhGY8YI+T2HcqUnMBEFoL2a8e2XeeVo6MYmZIAjtRaSgJ6LPAfhtAK8w89s1+98J4B8BPOtvepCZN/r7LgPwCQAlAJ9h5k0Z9TtTRBgLgtDJuKhuPg/gsog2/8LMy/yfQMiXAPwdgHcBOB/ADUR0fprOCoIgCPGJFPTM/A0AryU49sUAnmHmHzPzSQBfBnBVguMIgiAIKcjKGPvrRLSfiB4jogv8bX0AXlDavOhv00JENxPRXiLaOzo6mlG3BEEQhCwE/XcBnMvMFwL4HwAGkxyEmbcw83JmXr5gwYIMuiUIgiAAGQh6Zv45Mx/1/34UQJmIzgBQBXCO0vRsf5sgCILQRFK7VxLRLwH4KTMzEV0Mb/J4FcAYgLcR0WJ4Av56AL/ncsyhoaGfEdFzCbpzBoCfJfhc3ki/4lPUvkm/4iH9ikeafp1r2uHiXnkfgHcCOIOIXgRwJ4AyADDzpwBcC+BDRDQBYBzA9ezVJ5wgoj8GsAOee+XnmHnEpbfMnEh3Q0R7TaW0Won0Kz5F7Zv0Kx7Sr3jk1a9IQc/MN0Ts/ySATxr2PQrg0WRdEwRBELJgxqdAEARB6HQ6TdBvaXUHDEi/4lPUvkm/4iH9ikcu/SJPnS4IgiB0Kp22ohcEQRBCiKAXBEHocNpO0BPRe4hohIimiMjohkRElxHRISJ6hojWKdsXE9GT/vatRDQro36dTkQ7ieiH/u/5mjYriWif8vMGEQ34+z5PRM8q+5Y1q19+u0nl3NuV7a0cr2VE9C3/+/4eEV2n7Mt0vEz3i7J/tn/9z/jjsUjZd7u//RARrUnTjwT9+jMiesofn38monOVfdrvtIl9ex8RjSp9+ICy7yb/u/8hEd3U5H7do/TpB0Q0puzLZcyI6HNE9AoRfd+wn4jov/t9/h4R/ZqyL/1YMXNb/QD4FQBLAHwdwHJDmxKAHwF4C4BZAPYDON/fdz88X38A+BSAD2XUr78BsM7/ex2Av45ofzq8ZHE9/v+fB3BtDuPl1C8ARw3bWzZeAH4ZwNv8v88C8DKA3qzHy3a/KG3+I4BP+X9fD2Cr//f5fvvZABb7xyk1sV8rlXvoQ0G/bN9pE/v2PgCf1Hz2dAA/9n/P9/+e36x+hdr/CbwYn1zHDMC/A/BrAL5v2P9uAI8BIACXAngyy7FquxU9Mz/NzFGVu7WZM4mIAKwC8IDf7gsABjLq2lX+8VyPey2Ax5j5eEbnNxG3X9O0eryY+QfM/EP/75cAvAIgj0RILplW1f4+AODf++NzFYAvM/MJZn4WwDP+8ZrSL2berdxDe+ClGmkGabLTrgGwk5lfY+bDAHYiOhV6Xv26AcB9GZ3bCEdnAb4KwBfZYw+AXiI6ExmNVdsJekdMmTPfBGCMmSdC27Pgzcz8sv/3/wPw5oj216PxBvsr/7XtHiKa3eR+zSEve+ieQJ2EAo0Xeek1ZsFbrQVkNV4umVan2/jjcQTe+MTK0ppDv1TeD29VGKD7TrPCtW/X+N/RA0QU5L4qxJj5aq7FAHYpm/McMxumfmcyVoUsJUhEXwPwS5pdH2Xmf2x2fwJs/VL/YWYmIqPfqj9TL4WXHiLgdngCbxY8X9q/ALCxif06l5mrRPQWALuI6AA8YZaYjMfr7wHcxMxT/ubE49WJENGNAJYDeIeyueE7ZeYf6Y+QCw8DuI+ZTxDRf4D3RrSqieeP4noADzDzpLKt1WOWC4UU9Mz8mykPYcqc+Sq8V6Juf1UWK6OmrV9E9FMiOpOZX/YF0yuWQ/0ugIeYuaYcO1jdniCi/w3gz5vZL2au+r9/TERfB9APYBtaPF5E9K8APAJvkt+jHDvxeGlwybQatHmRiLoBzIN3P+WZpdXp2ET0m/Amz3cw84lgu+E7zUpoRfaNmV9V/v0MPLtM8Nl3hj779Wb1S+F6AB9WN+Q8ZjZM/c5krDpVdfMd+JkzyfMSuR7AdvasG7vh6ccB4CZ49W6zYLt/PJfjNugFfWEX6MUHAGit83n0i4jmB6oP8lJMrwDwVKvHy//uHoKnu3wgtC/L8dLeL5b+Xgtglz8+2wFcT55XzmIAbwPw7RR9idUvIuoH8GkAVzLzK8p27XeaUb9c+3am8u+VAJ72/94BYLXfx/kAVqP+7TbXfvl9Ow+ecfNbyra8x8zGdgB/4HvfXArgiL+YyWas8rAw5/kD4Hfg6alOAPgpgB3+9rMAPKq0ezeAH8CbjT+qbH8LvAfxGQBfATA7o369CcA/A/ghgK8BON3fvhxeYfSg3SJ4s3RX6PO7AByAJ7C+BOC0ZvULwL/1z73f//3+IowXgBsB1ADsU36W5TFeuvsFniroSv/vOf71P+OPx1uUz37U/9whAO/K+H6P6tfX/OcgGJ/tUd9pE/t2N4ARvw+7AZynfPaP/LF8BsAfNrNf/v/rAWwKfS63MYO3sHvZv59fhGdP+SCAD/r7CV6N7R/5516ufDb1WEkKBEEQhA6nU1U3giAIgo8IekEQhA5HBL0gCEKHI4JeEAShwxFBLwiC0OGIoBcEQehwRNALgiB0OP8fZCI1W25tGuYAAAAASUVORK5CYII=\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "AeC7Uv29KsXC"
      },
      "source": [
        "## T-Learner"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "EZZFqyFB9vn1",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# 集団を2つに分ける\n",
        "df_0 = df[df.Z == 0.0]  # 介入を受けていない集団\n",
        "df_1 = df[df.Z == 1.0]  # 介入を受けた集団\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "xp2P-IDT9vql",
        "colab_type": "code",
        "outputId": "de2ff6c7-58be-4219-cf6e-541066df7ccc",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 229
        }
      },
      "source": [
        "# ランダムフォレストモデルを作成\n",
        "from sklearn.ensemble import RandomForestRegressor\n",
        "\n",
        "# 介入を受けていないモデル\n",
        "reg_0 = RandomForestRegressor(max_depth=3)\n",
        "reg_0.fit(df_0[[\"x\"]], df_0[[\"Y\"]])\n",
        "\n",
        "# 介入を受けたモデル\n",
        "reg_1 = RandomForestRegressor(max_depth=3)\n",
        "reg_1.fit(df_1[[\"x\"]], df_1[[\"Y\"]])\n"
      ],
      "execution_count": 7,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "/usr/local/lib/python3.6/dist-packages/ipykernel_launcher.py:6: DataConversionWarning: A column-vector y was passed when a 1d array was expected. Please change the shape of y to (n_samples,), for example using ravel().\n",
            "  \n",
            "/usr/local/lib/python3.6/dist-packages/ipykernel_launcher.py:10: DataConversionWarning: A column-vector y was passed when a 1d array was expected. Please change the shape of y to (n_samples,), for example using ravel().\n",
            "  # Remove the CWD from sys.path while we load stuff.\n"
          ],
          "name": "stderr"
        },
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "RandomForestRegressor(bootstrap=True, ccp_alpha=0.0, criterion='mse',\n",
              "                      max_depth=3, max_features='auto', max_leaf_nodes=None,\n",
              "                      max_samples=None, min_impurity_decrease=0.0,\n",
              "                      min_impurity_split=None, min_samples_leaf=1,\n",
              "                      min_samples_split=2, min_weight_fraction_leaf=0.0,\n",
              "                      n_estimators=100, n_jobs=None, oob_score=False,\n",
              "                      random_state=None, verbose=0, warm_start=False)"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 7
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Rx2_b0CD9vqw",
        "colab_type": "code",
        "outputId": "21d57dea-1111-48ac-abc5-641b20874d57",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 34
        }
      },
      "source": [
        "# ATEを求める\n",
        "mu_0 = reg_0.predict(df[[\"x\"]])\n",
        "mu_1 = reg_1.predict(df[[\"x\"]])\n",
        "\n",
        "ATE = (mu_1-mu_0).mean()\n",
        "print(\"ATE：\", ATE)\n"
      ],
      "execution_count": 8,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "ATE： 0.6703869421982507\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "XvOWVBt99vq7",
        "colab_type": "code",
        "outputId": "e823b038-93c9-4dff-af7c-929be1e86cc0",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 52
        }
      },
      "source": [
        "# 処置群における平均処置効果ATTと、対照群における平均処置効果ATU\n",
        "ATT = df_1[\"Y\"] - reg_0.predict(df_1[[\"x\"]])\n",
        "ATU = reg_1.predict(df_0[[\"x\"]]) - df_0[\"Y\"]\n",
        "\n",
        "print(\"ATT：\", ATT.mean())\n",
        "print(\"ATU：\", ATU.mean())\n"
      ],
      "execution_count": 9,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "ATT： 0.5884244766831146\n",
            "ATU： 0.7470417076496992\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "8WXoUe6E9vrE",
        "colab_type": "code",
        "outputId": "aa6bc811-5ab4-46d8-f0b2-9011d57ccd9d",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 282
        }
      },
      "source": [
        "# 推定された治療効果を各人ごとに求めます\n",
        "t_estimated = reg_1.predict(\n",
        "    df[[\"x\"]]) - reg_0.predict(df[[\"x\"]])\n",
        "plt.scatter(df[[\"x\"]], t_estimated,\n",
        "            label=\"estimated_treatment-effect\")\n",
        "\n",
        "# 正解のグラフを作成\n",
        "x_index = np.arange(-1, 1, 0.01)\n",
        "t_ans = np.zeros(len(x_index))\n",
        "for i in range(len(x_index)):\n",
        "    if x_index[i] < 0:\n",
        "        t_ans[i] = 0.5\n",
        "    elif x_index[i] >= 0 and x_index[i] < 0.5:\n",
        "        t_ans[i] = 0.7\n",
        "    elif x_index[i] >= 0.5:\n",
        "        t_ans[i] = 1.0\n",
        "\n",
        "\n",
        "# 正解を描画\n",
        "plt.plot(x_index, t_ans, color='black', ls='--', label='Baseline')\n"
      ],
      "execution_count": 10,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "[<matplotlib.lines.Line2D at 0x7f6b86b9bda0>]"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 10
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3de3xU1bnw8d+TySQkIITITQIhEVOsijV9UwHpp/XQo6hUiZcKFHu0b9W21mMpyhFeOZVarPRQtdp6ar219oiIok5RRGpV1BcIEhowgnKHwABCDeGSBJJM1vlj9sTJ3DOZO8/388mHmb3W7P1kT3hmzVprry3GGJRSSqW/rGQHoJRSKjY0oSulVIbQhK6UUhlCE7pSSmUITehKKZUhspN14H79+pmSkpJkHV4ppdLSunXr/mmM6R+oLGkJvaSkhOrq6mQdXiml0pKI7A5Wpl0uSimVITShK6VUhtCErpRSGUITulJKZQhN6EoplSGSNstFKaVONY4aJ/OXb2ZfQzODC/KYMX4EleVFMdu/JnSllEoAR42TWa/U0tzqAsDZ0MysV2oBYpbUtctFKaUSYP7yzR3J3KO51cX85ZtjdgxN6EoplQD7Gpq7tD0amtCVUioBBhfkdWl7NDShK6VUAswYP4I8u63Ttjy7jRnjR8TsGDooqpRSCeAZ+NRZLkoplQEqy4timsB9aZeLUkplCE3oSimVITShK6VUhtCErpRSGUITulJKZQid5aKUUgkQ74W5QBO6UkrFXSIW5oIIErqIPAN8GzhojDkvQLkAjwBXAE3ATcaYf8QsQqWUSiFTn1zNyu31QctzbMJ/XfeVTok61MJcib6w6M/A74G/BCm/HCizfkYBf7D+VUqptOfpKnFGuIhWi8swbdF6pi1aH7ZuLBfmgggSujHmfREpCVFlIvAXY4wBqkSkQETOMMbsj1GMSqk00NjYyNKlS/22l5eXU1ZWRkNDA3/729/8yr/2ta9RWlrKP//5T9555x2/8osuuoghQ4Zw4MAB3n//fb/yb3zjGwwaNIi9e/eyatUqv/Jx48bRr18/du7cydq1a/3KL730UgoKCti6dSs1NTWdyl5YW0e1q5SsnB60/nMPLYd2+b0+v2w0km2n5eBOWj/f618+4iIky8bJA9vI7t0fW36fjrJYLswFselDLwL2eD3fa23zS+gicitwK0BxcXEMDq2USqbjx4/z7LPPMm7cOPLy8pg0aZJfnUcffZSysjL27NkTsPxPf/oTpaWlbNmyJWD54sWLGTJkCBs2bAhYvnz5cgYNGsSaNWsClq9atYp+/frx3nvv8f3vf9+v/KOPPqKgoIA333yTO+64w6+86EdPk5XTg6atq2l437+jYsgdz2PLttP4yXscrVrsV15816uAjcbav5N/9tc7ErpATBfmAhB3wzpMJXcL/fUgfeivA/OMMf/fev42cLcxpjrUPisqKkx1dcgqSqkUt3v3bkpKSnjmmWeYOnUq27Zt86szaNAgCgsLOXHiBDt27PArHzx4MAUFBTQ1NbFr1y6/8iFDhtC7d2+OHz9OXV2dX3lxcTG9evXi6NGj7N3r30IuKSkhPz+fhoYG9u3b51d+5pln0qNHD+rr6zlw4EDH9qlPVvHZsZPY+56B2Oy4mo7gajri93p7YRGSZcPVeBhX8zH/8tOHIiK0Ha8nKyefrJweHWW75k3wqx+OiKwzxlQEKotFC90JDPV6PsTappTKcK2trQBkZ2eTk5PDOeecE7Rujx49Qpbn5+eHLO/Vq1fI8t69e4csLygooKCgIGh5YWEhhYWFHc8bcneSk/tFuS2/T6fuEl+2nn2x9ewbtDy7V2HQsliJxYVFS4B/E7fRwBHtP1fq1NDW1ga4E3qmiXX/diJEMm1xIXAx0E9E9gL3AnYAY8zjwBu4pyxuwz1t0b+TSimVkTwJ3W63JzmS2MvPCd/eLbIuEAK6NBMGYOzw2LfYI5nlMiVMuQF+ErOIlFJpI1Nb6Jc8tIKtBxsDlt0wupi5lSP9tvvOJ5/tqGVBVR2BRinHDi9kwS1jYhFqJxENisaDDooqlf5OnjzJvn37GDBgAD179kx2ODEx21HLc1X+g68e0QxkxlK8B0WVUqeo3NxcSktLkx1GTDhqnMx4aT2t7cmOJHq62qJSKmp79uxh3rx5AacbppOpT65m2qL0TuagLXSlVDfs2LGDWbNmMWrUKEpKSpIWR6DL820iTBk1NGB/t/frIrlEP11oC10pFbVUGBT1rGToO8PEZQzPVdVxyUMrgr6uq8k8HjNTYkkTulIqaqmQ0AOtZOht68FGpj65OuDruiJeM1NiSbtclFJRS2ZCDzcbxdvK7fXMdtR26n6JdKXDdEjkHtpCV0pFLVkJvSvJ3OO5qjpmO2o7nkdyJegNo4vTJpmDJnSlVDdcfvnlHDx4kPPO81u3L64WrtkTvlIAz1XV4ahxLzUVaqXDgaflsGvehJADqqlIE7pSKmo5OTn0798/4Zf+u7pxQaSn77yyvIjfTroAu08WHDu8kDX3XNKd8JJG+9CVUlGrqalh8eLFTJ8+ndNPPz1hx7WJRJ3UnQ3NlMzsfCOO3rk2PvrFZbEILam0ha6UitqGDRv41a9+xdGjRxN63Cmjhoav1AVHT7o4/943Y7rPZNCErpSKWrIGRedWjuSG0f53PRs7vNDqRpEu7/PoyeBTH9OFdrkopaKWzGmLcytHhhy0nLNkIw3NrV3ap6PG6bdqYjrRFrpSKmoul7tVm2rL51aWF7H+3kvpaju9qxcbpRpN6EqpqKXClaKhdPWuQ5FebJSqNKErpaJ2++2309zcHPJenckUaq55IOl42zlvmtCVUlGz2Wz06NEDka4PQiZCZXkRA0/LCVjmO26aZ7d1+QMg1WhCV0pFbdmyZdxxxx20t6fmQuKzHbUcPNbit71sQE8euv4CigryENz3Bn3gmpFpPSAKOstFKdUNa9as4Xe/+x2PPPJIskPx46hxBl3vpamlncryorRP4L60ha6UilpbWxs2my0lu1zuebU2aFm6D34GowldKRW1tra2lJ3h0tgS/EKhdB/8DCY13wmlVFpIVkL33HJuX0MzgwvymDF+RJe6T9J98DMYTehKqW7Jzc2Ny34dNU7+3ysf0WTduVmAPHtWx3MPZ0Mzs15xd694J3UBgi3flWl95x5iIlixTEQuAx4BbMBTxph5PuXDgGeA/kA9cIMxZm+ofVZUVJjq6upo41ZKpbHZjloWVNUFTbjRKCrIY+XMcZ2OEWhQ9IbRxWm3zrk3EVlnjKkIVBa2hS4iNuAx4BJgL7BWRJYYYzZ5VfsN8BdjzLMiMg54APhe90NXSmUC942cP6K5NX7TG30HOj1Je+GaPbiMwSbClFFD0zqZhxNJl8uFwDZjzA4AEXkBmAh4J/RzgOnW43cBRyyDVEqlpqeffppNmzbx4IMPBq3jqHEyfdF64j1TPdBAZ7gFvDJNJLNcigDv+z3ttbZ52wBcYz2+GjhNRPxWuxeRW0WkWkSqDx06FE28SqkU8sEHH/Dyyy8HLJvtqKV01lKmJSCZC5k70NkVsZq2eBfwTRGpAb4JOAG/OUPGmCeMMRXGmIr+/fvH6NBKqWQJNsvF03/djTvFRUyAqaOLM3agsysi6XJxAt63BxlibetgjNmH1UIXkV7AtcaYhlgFqZRKTcESerQ3ce6qoiimLGaySBL6WqBMREpxJ/LJwHe9K4hIP6DeGNMOzMI940UpleGCJfTu3MQ5nDy7LSPWXYmHsAndGNMmIrcDy3FPW3zGGLNRRO4Dqo0xS4CLgQdExADvAz+JY8xKqRSRn59P3759/bZHehPnUFMIPRcPORuaO/anLfLQIpqHHg86D12pzBVsDrivXfMmJCCazNKteehKKdVVnlb3gjXBB0YD3eRZdY+20JVSUbv33ntpb2/nl7/8Zch6sx21p9QFPvGkLXSlVFy89957EdU71S7wSRZdPlcpFbVUXj73VKQJXSkVNU3oqUUTulIqaprQU4u+E0qpqA0YMICBAwd22jb1ydWs3F7f8Xzs8EIW3DIm0aGdknSWi1IqZi55aAVbDzb6bdekHjuhZrlol4tSKiYcNc6AyRzo1GJX8aNdLkqpqN14441QMIRtZ3wLp88NJlTiaUJXSkVt2VvvcKLfCAqvuCjZoSi0y0Up1Q0NjScwYgtbb+zwwgREozShK6Wi5nK1QVboNFI2oKcOiCaIJnSlVPTa25Gs0C30t6ZfnJhYlCZ0pVT0evQbgq2X3+2DO0gCY1Ga0JVS3bDgr3+j75jrg5ZP1SVyE0pnuSiloua5c9CcJRtpaG7t2C4CU0cFvxuRig9N6EqpqI0dO5Ybb7yR9ffemuxQFNrlopSKkjGGVatWsW/fvmSHoiya0JVSUXG5XAC62mIK0YSulIpKW1sboAk9lWhCV0pFRRN66tGErpSK2pgxYygqKkp2GMoS0UeriFwGPALYgKeMMfN8youBZ4ECq85MY8wbMY5VKZVCevXqxapVq5IdhvIStoUuIjbgMeBy4Bxgioic41NtNvCiMaYcmAz8d6wDVUopFVokXS4XAtuMMTuMMS3AC8BEnzoG6G097gPoPCalMtyBAwc477zzeOWVV5IdirJEktCLgD1ez/da27zNAW4Qkb3AG8C/B9qRiNwqItUiUn3o0KEowlVKpYoTJ06wceNGjh49muxQlCVWg6JTgD8bY4YAVwD/IyJ++zbGPGGMqTDGVPTv3z9Gh1ZKJYPOckk9kSR0JzDU6/kQa5u3HwAvAhhjVgM9gH6xCFAplZo0oaeeSBL6WqBMREpFJAf3oOcSnzp1wLcAROTLuBO69qkolcFaW92LcWlCTx1hE7oxpg24HVgOfIJ7NstGEblPRK6yqt0J3CIiG4CFwE3GGBOvoJVSyZefn8+ll17KoEGDkh2Kskiy8m5FRYWprq5OyrGVUipdicg6Y0xFoDK9UlQppTKEJnSlVFQ+/PBDhg0bxsqVK5MdirJoQldKRaWpqYm6urqOwVGVfJrQlVJR0WmLqUcTulIqKprQU48mdKVUVDShpx5N6EqpqAwcOJBrrrmGwsLCZIeiLPrRqpSKyte+9jVefvnlZIehvGgLXSmlMoQmdKVUVBYtWsTpp5/Ozp07kx2KsmiXi1JxMttRy8I1e3AZg02EKaOGMrdyZLLDipnm5mbq6+sRkWSHoiya0JWKg9mOWp6rqut47jKm43mmJHWd5ZJ6tMtFqThYuGZPl7anm6lPrubul2oAGDNvBVOfXJ3kiBRoQlcqLlxBVjENtj2dTH1yNSu312PaXe4NWVms3F6vST0FaEJXKg5sIfqVHTW+N/xKLyu31wNgLxxCz/O+hWTndtqukkcTulJxMGXU0KBl85dvTmAk8ZNXWk6/CT8jK6dHskNRFk3oSsVBqIHPfQ3NCYxEnUp0eFqpOCkqyMMZIHkboGTmUsR6DNA33869V55L9e76lJ/qOHZ4ISu313Nk1SIaVi6k+M6XkSwbY4frEgDJpi10peLk82MnQpZ7D48ebmpl2qL1PFdV1zFw6pnqONtRG8cou27BLWMYO7wQ094G7W0gWYwdXsiCW8YkO7RTniZ0peLkhCs2M1pScarjglvGcPvFZ2Kz2dj9629rMk8R2uWiVBiOGiezXvmI5tZ2v7KeOTbuv3okleVFHXX/Y/EGWmKUzCF1pzq2tbXpRUUpRt8NpUJw1DiZvmg9/qncrbHFxbRF65m2aH1C40oFmtBTj74bSoUwf/nmoMn8VDd27Nhkh6B8RNSHLiKXichmEdkmIjMDlD8sIuutny0i0hD7UJVKvECzVJLhy/+5LOUuSLr66qv5zW9+k+wwlJewLXQRsQGPAZcAe4G1IrLEGLPJU8cY8zOv+v8OlMchVqXiznuFxFTS3NrOdKtbx9Nfn2wtLS0A5OTkJDkS5RFJl8uFwDZjzA4AEXkBmAhsClJ/CnBvbMJTKnE8KyS2Ha/HtHRumWf37o9k59B+shFXo/8X0Ow+AxCbHdeJ47Q3HQlQPhCxZeNqPkp78zFybNJp4DS7YBCSZcPVdIT2E8f9X993MO0i3L94Nef0PL9TmYhQVlYGwP79+zl27FincpvNxvDhwwFwOp00NjZ2Krfb7ZSWlgKwZ88emps7/+65ubkMGzYMgN27d3Py5EkA7rnnHlavXs3evXv94lXJEUlCLwK8503tBUYFqigiw4BS4J3uh6ZUYj2/po7Wz/ey76kf+ZUN+reHyT2jjMZPPqB++e/9ygff/Dj204fQ+NFbHH73ab/yIbf9Gdtp/ZCNy9n39rN+5UOnvYjk5nO0ajFH177qV178H68B8OnSpxnxqzc7lfXs2ZPjx90fAnfddRfPP/9859gHDWL//v0A/PjHP+a1117rVN7j9CIG3fxHBhfkcfTl2Xxc3XmRrfLycv7xj38AcO2117Ju3bqOsi996Ut+sarkifWg6GRgsTHGFahQRG4FbgUoLi6O8aGV6p52A2LPpVf5BLJ7FZJdMLCjzPO4x7Dz6XflXeRlZ9Hc9sVwqa2X+yrJvDMrsPXqC0C+3ca1/2cIFSWFVFZWkp+fT21tMbW1lwJQvaue1z/az+GmFor69aG5DU6e801yBg0PGmOv8y/hj//xvU7bvGea3HbbbUyYMKFTeY8eX6y1cueddzJ58uSO47+wdg9ttlwM7vGC9i9dyc+u/i4VJV9c9dm3b9+Ox3PnzqW+/otFuEaOTK2rWE91YsL0FYrIGGCOMWa89XwWgDHmgQB1a4CfGGNWhTtwRUWFqa6ujipopeKhZObSiOp5ror0LCPruz1a7vnutTS3BmwPddg1r3PCdtQ4+cVrGznc1BpR3B3P570TcNC3qCCPlTPHBdxHpt+FKR2IyDpjTEWgskha6GuBMhEpBZy4W+HfDXCQs4G+gC6KrNJWe8sJ2puPYOtViNjsfuXeSTHWV0d6BjvnL98cdHaN77K8jhonMxZvoDWCC5k8a5Z74g62SFiw7afCXZjSXdiEboxpE5HbgeWADXjGGLNRRO4Dqo0xS6yqk4EXTLgmv1Ip7MSuGg69ej9n3PQoOQPP7Nh+w+jihCStyvIiKsuL/JKnh8uYjm8SZQN60tTSHlEy9/D+RjE4yOJhgwvyOj131DiZs2QjDc2BvwEsXLNHE3qKiKgP3RjzBvCGz7af+zyfE7uwlEo8mwjGZSUtW3an7YlOWHMrR7Lz0PGQN43YerAxaFkojhonleVFzBg/wq+LJ89uY8b4ER3PL3loRdjjpNoUz1OZXimqlGXKqKE8Xuu+8bF4JfRQN6uIp12fx+eiJt/57N5dPM2tLu58cQPVu+vZeeh4RB8aoe7OpBJLV1tUyjK3ciSjS3oDIDY7NpGEdbUEEq8bYbQD0xat55z/XAbAv5zdv1O5p2880lvKJesDT/nTFrpSXiacO4ClwLp7L2fgwIFh68dTsD7uWGlqbWf6i+vpTo9Jvj1L+89TiLbQlfLy9a9/nUcffZTevXsnOxRmjB9BtJ0ZRQV5FPkMbgbSbjrfaKMr7FnCr645P3xFlTDaQlfKy8iRI1PmYhlPH3dX11e3Z0nHwGa8lvW1Z8H873wlZdaVUW6a0JXycuDAAQ4dOsR5552HpMBgn2caI0DpzKVhW9MFeXbmXHVux2uqd9cHnP7YHckcV1ChaUJXyssf/vAH7rvvPtrbU28V9FB96r+ddEHA1vLcypFUDCvknldraWwJfQVqOAI8HOQ4KjVoH7pSXlpaWsjJyUmJ1rmvGeNHkGe3ddqWZ7cFTeYeleVFbLzvMm4YXUx3fq2po4s1mac4TehKefEk9FRUWV7EA9eMpKggD8E98PnANSMjTrJzK0ey84EJ3DC6awvj2bPc3wC0myX1aZeLUl5SOaFD5z71aM2tHMmCqrqQ/fHaT56eNKEr5SXVE3qshErmvqs5qvShCV0pLzfddBPjxgVeOjaT2EQCrsGil/GnN+1DV8rLmDFjmDRpUrLDiLtgl+vrZfzpTRO6Ul42bdrEhg0bkh1G3M2tHMkNo4s7WuTJXrdGxUbYOxbFi96xSKWiiRMnUldXR01NTbJDUSqgUHcs0ha6Ul5OlUFRlZk0oSvlRRO6Smea0JXyogldpTNN6Ep50YSu0pnOQ1fKy69//Wvsdnuyw1AqKprQlfJy8cUXJzsEpaKmXS5KeVmxYgUff/xxssNQKiqa0JXycuONN/Lggw8mOwyloqIJXSkvOiiq0llECV1ELhORzSKyTURmBqlzvYhsEpGNIvJ8bMNUKjE0oat0FnZQVERswGPAJcBeYK2ILDHGbPKqUwbMAsYaYw6LyIB4BaxUPGlCV+kskhb6hcA2Y8wOY0wL8AIw0afOLcBjxpjDAMaYg7ENU6nE0ISu0lkk0xaLgD1ez/cCo3zqfAlARFYCNmCOMeZN3x2JyK3ArQDFxV27DZZSibBs2TKGDtUlZFV6itU89GygDLgYGAK8LyIjjTEN3pWMMU8AT4B7tcUYHVupmDkVbm6hMlckXS5OwLvJMsTa5m0vsMQY02qM2QlswZ3glUobra2tLFq0iC1btiQ7FKWiEklCXwuUiUipiOQAk4ElPnUcuFvniEg/3F0wO2IYp1Jxd+zYMSZPnsyyZcuSHYpSUQmb0I0xbcDtwHLgE+BFY8xGEblPRK6yqi0HPheRTcC7wAxjzOfxClqpeGhpaQHQQVGVtiLqQzfGvAG84bPt516PDTDd+lEqLbW2tgKa0FX60itFlbJoC12lO03oSlk0oat0p8vnKmUpKSmhqqqK4cOHJzsUpaKiCV0pS15eHqNG+V4zp1T60C4XpSz79+/n6aefZv/+/ckORamoaEJXyvLpp59y8803s3Xr1mSHolRUNKErZdFBUZXuNKErZdGErtKdJnSlLJ6EbrfbkxyJUtHRhK6U5eTJk4C20FX60oSulOXb3/42GzdupKSkJNmhKBUVnYeulKV3796cc845yQ5DqahpQo8xR42T+cs3s6+hmcEFecwYP4LK8qJkh6UiUFNTw/vvv8+tt95KXl5essNRqsvEvVBi4lVUVJjq6uqkHDsepj65mpXb60PWyc3O4tfXnq8JPkU9/PDDTJ8+nYaGBvr06ZPscJQKSETWGWMqApVpCz0GLnloBVsPNoatd7KtnWmL1jNt0XpsIkwZNZS5lSMTEKGKhA6KqnSnCb2bHDXOiJK5L5cxPFdVx3NVdRRp10xK0HnoKt1pQu8GR42Tny1a3+39OBuaO1ruY4cXsuCWMcx21LKgqg7vDrEbRhdriz6OWlpayMrKwmazJTsUpaKiCT1Kjhon019cT6xHIFZur2fU/W/x2bEWv7LnquoAMjapz3bUsnDNHlzGJKVLqqWlhdzc3IQdT6lY00HRKI2d9w7OhuakHPu3ky7IuO6Z2Y7ajg+sYERg6qj4fUs5cuQIR48eZejQoXHZv1KxoIOicbAvSckcYMbiDQAZldQXrtkTto4x8f2W0qdPH53dotKaXinaBY4aJ8NnLaVk5tKwXS1FBXnsmjeBscMLYx5Hq8swbdF6Sma6Y7ngF3/DUeOM+XESydWFb4qRJP9o/PWvf+WRRx6Jy76VSgRtoUfIUeNkWoQDoFnAjPEjAL4Y4FxThydnCcS0772huZVpi9ZTvbs+5fvXfccHBHh40gXYRCJO6l1J/l3x8ssv88EHH/DTn/40LvtXKt40oUdo/vLNEdXLs2fxwDWdLx6aWzkyYKLtyodEJBZU1VExrDBlu2ICDfYa6JjdE+7CLG+zHbUx//BqbW3VKYsqrUWU0EXkMuARwAY8ZYyZ51N+EzAf8Hzv/70x5qkYxtltjhonM15aT2t75+2eaYLhRNJnLsAnv7w84pg8iTdWSd3g/uBJ1YQeaOaOx8rt9diz8Ht/golHX3pLS4smdJXWwiZ0EbEBjwGXAHuBtSKyxBizyafqImPM7XGIMSqOGie/eG0jh5taQ9Zbub2es+95g0/vvyJkvcEFeWFntQwu6Pr6H57kO33ReiLMZSHFerB2tqOW59fU0W71ctizoGeunSPNrTFfq6a1Hew2ITtLaI4gsy9cs0cTulJeImmhXwhsM8bsABCRF4CJgG9CT5h3332Xffv2ddrWu3dvrrzySgDm/GEhj7/5D9rav+hrtfU4jbzh7pk+zdvX4jpxvKPsOPDVH33CPx6/E4AlS5Zw7NixTvu/vF8uTzW4E3bj5pWYts6tzR59BzFj0hQAXnjhBVwuV6fys846i1GjRmGM4fnnn/f7nRZf/2W++tWv8tKHO7nnwSc53NRC3/wcJpx/BgAv7sjC3r+E9tYTNG1Z7ff63IFnYe83lAE9XCxYsMCvfNSoUZx11ll8/vnnvPnmm37lY8eOpaSkhAMHDvD222/z2Ltb2fLZF+eoR/FIsk/rR3PDIQ7v+RiAzcBtVW/y3teGMvPm6xk4cCA7d+5k1apVfvt/9VA/AFrrnZzcv8WvPP+sUWTl5tP42W56HdvD/VedC8BPX6hxl3/pIrLsubQc3EHLod0dr1uw4DAAkyZNIjs7m3Xr1vHpp5/67X/q1KkAVFVVsX379k5ldrud66+/XhO6Sn/GmJA/wHW4u1k8z7+Hu0vFu85NwH7gI2AxMDTIvm4FqoHq4uJiE60rrrjC4O5h6PgZMWJER3l+8Xl+5TmDysywu183w+5+3dgHnOlXnlt8fsfrzzrrLL/yq666yrz6j73mzJmvm6yeBX7l37i88ovj5+f7lf/whz80xhjjcrn8ygBz1113GWOMaWhoCFg++Ud3mtKZr5ui2/4csLzvuJvN2bOXmd+9vCJg+RNPPGGMMebDDz8MWL5w4UJjjDFvv/12wPL+191rht39uul/zeyA5e+9954xxpiJ0x4IWH7GTY+YYXe/bgovvS1g+eBb/miG3f26Kbj4+wHLh/zkf8ywu183fS6aHLD8+PHjxhhjpk2b5lcmIh3vzc033+xX3rt3b2OMMY2NjebIkSNR/10qlQhAtQmSr8NeWCQi1wGXGWNutp5/DxhlvLpXROR04Lgx5qSI/BCYZIwZF2q/3bmwyOl00tTU1Gmb3W6npKSE2Y5a/vy3dX4taMm2k917AABtR1huyvkAAAzZSURBVA9i2lp9ynPZ+983ArBr1y5aWzuX9+zZk8GDBwOwY8cOvxb4aaedxqBBgwDYtm0bvue1T58+DBgwAGMM27Zt8/udCgoK6N+/Py6Xix07dviVFxYW8kHdCeY4PuLQ/i+m7Ym452cPOWMgsyoruOLc/uzevdvv9QMGDKBPnz6cOHGCPXs6T/t7a9MB/mvFAcjNp731BK5jn/u93tbrdLJyetDe0ozreOfBSwEe/P44/vP1LZxoasTVeNjv9dm9+yPZObSfOI6r6Yh/eZ8BiM2O68Rx+mefZNEP3eMaD7+1mb+u30d2wSAky4ar+Sjtze5vT7YsYeblZ/Ox8whv7DYgWbiajtB+4jgTLxjMzy4Z0bH/srIyAD777DOOHj3a6dhZWVkMHz7cLyalUlGoC4siSehjgDnGmPHW81kAxpgHgtS3AfXGmJBXaMTrStHhs96IelrbrnkTYhxN6ovFTJuCPDuNLW20uro/nTDPbuOBa0Z26peP5CrSQMoG9OSt6Rd3OyalUkmohB7JhUVrgTIRKRWRHGAysMTnAGd4Pb0K+CTaYLsr2mReNqBnjCNJHEeNk7Hz3qF05lLGznunSxcZRTodM5g8uw0RupXMRdz/FhXk+SVzcM9kkSj2u/VgI7MdtVHHpVS6CTsoaoxpE5HbgeW4py0+Y4zZKCL34e7LWQLcISJXAW1APe4+9aQId4HKwNNy6J1n77TkbTq35Bw1Tma9Uktzq7sLyNnQzKxX3Eksktkn3V2P5qvFfbo0f9xXpN+KIpllFMjza+pS/mIrpWIlonnoxpg3gDd8tv3c6/EsYFZsQ4vOmf3zA65PHul883Qzf/nmjmTu0dzqimg+uqPGGfFVq/YsuLDU/+Kf7iTzG0YXR1x3xvgRnT64ItXe/V4gpdJGxlwpGu4WcLs+T95iWvEUrNUayXz0+cs3h03m9ixh/ne+QmV5EcNnvRGmdmSiWRo31hdhKZWJ0j6hO2qc3PXShk5zzgPxTXCOGidzlmykoTnwhUfp0KIP1T9sgJKZSykb0JO9h5s7LtTJEviutQRtuKRfkGdnzlXndiTTro5PxHqZ38ryIqp313dpgDTPruvPqVNHWid0R42TGYvDJ3PofBWnexmADbSGeN3K7fVMfXJ1Sif1SFYd9O1+avdagjZYv3RRQR4rZ/rPOu3KAloQn+V9Pa16z40wwsnNtuGocYaNJdk311AqFtI6oc9fvjni2RWe1Q87XhfBh0B3+ocToTurDi5cs4cHr/+KX790nt3W6Vx5mzJqaMSt43g2jEMtdua73ENDcyszXgq9frxvd53nfq+eYymVLtL6+2ik65b0zLF1+s/clfVOSmYu5dyfv5mS643bJJrJfG4uY6gsL+KBa0ZSVJCHEHzaoMfcypERD2TO/84FUccWrcryIgJ9xrW2u9ePL7/Pf914R40z6Ad3vNZdVype0rqFnmfPoinMIk52m3D/1Z1bWV2dAtfY4uLOMK28ZOhKi9mXTQRHjZP5yzezr6E54oW2PK3jSx5aEXA2kT3LncyTdZ6CjYkAHG5q7bgZN7jHCEJ9JsZr3XWl4iXtErqjxsk9r9bS2BJ++lpRkCQ1Y/yIsH3ovlzthjlLNqZUQvd0ByyoquvyDTNGn9m3W/PX03XevrdQyR+69w1IqWRIq4TuqHFy50sbcEWQiENdsOJJWKFmuQTS0Nwa0QBbIvn2J/v2Bweb5fLup4einr+eyvrm28MumRypKaP0ZtEqvaRVQp+/fHNEybwognXJK8uLAiaucGvBpHrCi3RWTunMpQG3J/Pm17Fw75XnxmSuetmAnjogqtJOWiX0SJNNsFkakQjXL52uCS/SG35Ec5OOVBLNXPVAmlpicbsRpRIrrWa5RJpsutOCDjeTIx0Tnme+fiRdEd35MEwVcytH8ttJF9A33x71PtL1g1ud2tIqoc8YPwJbVuiBqq6sDxKMJyHk2W2dtoeao53KujJfP5W7k7qisryImp9fyq55E9g1b0KX/y7S8YNbqbTqcvEkm0CzXLwvaY/lseYv34yzoRmbSMegoXd5OtDW5hczgrzvjxqM3SZp+cGtVFoldAg+mBmvYwHdmt6XCqJdejbTeGYEOWqcQW/K3Tffzr1Xnps2761S3tIuoSdad5anTRUzxo9gxuINYbtdxg4vTFBEyZXIRoFSiaQJPYxg3RXp1I3hSV6hZrmUDeiZ0guRKaXC04QeRrDuinQbNPNulUZzyb9SKvVpQg8j0J1y0nW2i4d2OSiVmTShh+E920VbtEqpVKYJPQLaolVKpYO0urBIKaVUcJrQlVIqQ2hCV0qpDKEJXSmlMoQmdKWUyhBiknTfRBE5BOyO8uX9gH/GMJxY0bi6JlXjgtSNTePqmkyMa5gxpn+ggqQl9O4QkWpjTEWy4/ClcXVNqsYFqRubxtU1p1pc2uWilFIZQhO6UkpliHRN6E8kO4AgNK6uSdW4IHVj07i65pSKKy370JVSSvlL1xa6UkopH5rQlVIqQ6RsQheR74jIRhFpF5Gg03tE5DIR2Swi20Rkptf2UhFZY21fJCI5MYqrUETeEpGt1r99A9T5FxFZ7/VzQkQqrbI/i8hOr7ILEhWXVc/ldewlXtuTeb4uEJHV1vv9kYhM8iqL6fkK9vfiVZ5r/f7brPNR4lU2y9q+WUTGdyeOKOKaLiKbrPPztogM8yoL+J4mKK6bROSQ1/Fv9iq70Xrft4rIjQmO62GvmLaISINXWTzP1zMiclBEPg5SLiLyqBX3RyLyVa+y7p8vY0xK/gBfBkYAK4CKIHVswHbgTCAH2ACcY5W9CEy2Hj8O/DhGcf0XMNN6PBP4dZj6hUA9kG89/zNwXRzOV0RxAceDbE/a+QK+BJRZjwcD+4GCWJ+vUH8vXnVuAx63Hk8GFlmPz7Hq5wKl1n5sCYzrX7z+hn7siSvUe5qguG4Cfh/gtYXADuvfvtbjvomKy6f+vwPPxPt8Wfv+BvBV4OMg5VcAywABRgNrYnm+UraFboz5xBizOUy1C4FtxpgdxpgW4AVgoogIMA5YbNV7FqiMUWgTrf1Fut/rgGXGmKYYHT+YrsbVIdnnyxizxRiz1Xq8DzgIBLwSrpsC/r2EiHcx8C3r/EwEXjDGnDTG7AS2WftLSFzGmHe9/oaqgCExOna34gphPPCWMabeGHMYeAu4LElxTQEWxujYIRlj3sfdgAtmIvAX41YFFIjIGcTofKVsQo9QEbDH6/lea9vpQIMxps1neywMNMbstx4fAAaGqT8Z/z+m+62vWw+LSG6C4+ohItUiUuXpBiKFzpeIXIi71bXda3Oszlewv5eAdazzcQT3+YnktfGMy9sPcLfyPAK9p4mM61rr/VksIkO7+Np4xoXVNVUKvOO1OV7nKxLBYo/J+UrqHYtE5O/AoABF9xhj/proeDxCxeX9xBhjRCTovE/rk3cksNxr8yzciS0H91zUu4H7EhjXMGOMU0TOBN4RkVrcSStqMT5f/wPcaIxptzZHfb4ykYjcAFQA3/Ta7PeeGmO2B95DzL0GLDTGnBSRH+L+djMuQceOxGRgsTHG5bUtmecrrpKa0I0x/9rNXTiBoV7Ph1jbPsf9VSbbamV5tnc7LhH5TETOMMbstxLQwRC7uh541RjT6rVvT2v1pIj8CbgrkXEZY5zWvztEZAVQDrxMks+XiPQGluL+MK/y2nfU5yuAYH8vgersFZFsoA/uv6dIXhvPuBCRf8X9IflNY8xJz/Yg72ksElTYuIwxn3s9fQr3mInntRf7vHZFDGKKKC4vk4GfeG+I4/mKRLDYY3K+0r3LZS1QJu4ZGjm437wlxj3K8C7u/muAG4FYtfiXWPuLZL9+fXdWUvP0W1cCAUfD4xGXiPT1dFmISD9gLLAp2efLeu9exd23uNinLJbnK+DfS4h4rwPesc7PEmCyuGfBlAJlwIfdiKVLcYlIOfBH4CpjzEGv7QHf0wTGdYbX06uAT6zHy4FLrfj6ApfS+ZtqXOOyYjsb9wDjaq9t8TxfkVgC/Js122U0cMRqtMTmfMVrtLe7P8DVuPuRTgKfAcut7YOBN7zqXQFswf0Je4/X9jNx/4fbBrwE5MYortOBt4GtwN+BQmt7BfCUV70S3J+6WT6vfweoxZ2YngN6JSou4CLr2Busf3+QCucLuAFoBdZ7/VwQj/MV6O8FdxfOVdbjHtbvv806H2d6vfYe63Wbgctj/PceLq6/W/8PPOdnSbj3NEFxPQBstI7/LnC212v/r3UetwHfT2Rc1vM5wDyf18X7fC3EPUurFXf++gHwI+BHVrkAj1lx1+I1gy8W50sv/VdKqQyR7l0uSimlLJrQlVIqQ2hCV0qpDKEJXSmlMoQmdKWUyhCa0JVSKkNoQldKqQzxvwbvDlUCVqolAAAAAElFTkSuQmCC\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "tmYzmE5NTf_p",
        "colab_type": "text"
      },
      "source": [
        "## S-Learner"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "oLPw7a27TdTA",
        "colab_type": "code",
        "outputId": "55ed183f-853a-4efd-e217-2bd8ec8ebffe",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 194
        }
      },
      "source": [
        "# ランダムフォレストモデルを作成\n",
        "from sklearn.ensemble import RandomForestRegressor\n",
        "\n",
        "# モデルを学習\n",
        "reg = RandomForestRegressor(max_depth=4)\n",
        "X = df.loc[:, [\"x\", \"Z\"]]\n",
        "reg.fit(X, df[[\"Y\"]])\n"
      ],
      "execution_count": 11,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "/usr/local/lib/python3.6/dist-packages/ipykernel_launcher.py:7: DataConversionWarning: A column-vector y was passed when a 1d array was expected. Please change the shape of y to (n_samples,), for example using ravel().\n",
            "  import sys\n"
          ],
          "name": "stderr"
        },
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "RandomForestRegressor(bootstrap=True, ccp_alpha=0.0, criterion='mse',\n",
              "                      max_depth=4, max_features='auto', max_leaf_nodes=None,\n",
              "                      max_samples=None, min_impurity_decrease=0.0,\n",
              "                      min_impurity_split=None, min_samples_leaf=1,\n",
              "                      min_samples_split=2, min_weight_fraction_leaf=0.0,\n",
              "                      n_estimators=100, n_jobs=None, oob_score=False,\n",
              "                      random_state=None, verbose=0, warm_start=False)"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 11
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "IJmgnebQTfkB",
        "colab_type": "code",
        "outputId": "0463d635-75d8-4900-d3dc-e7c4211f972c",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 34
        }
      },
      "source": [
        "# 処置が0と1の状態を作成する\n",
        "X_0 = X.copy()\n",
        "X_0[\"Z\"] = 0.0\n",
        "\n",
        "X_1 = X.copy()\n",
        "X_1[\"Z\"] = 1.0\n",
        "\n",
        "# ATEの計算\n",
        "ATE = (reg.predict(X_1)-reg.predict(X_0)).mean()\n",
        "print(\"ATE：\", ATE)\n"
      ],
      "execution_count": 12,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "ATE： 0.6723197979076997\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "sv1NNZi8TjUJ",
        "colab_type": "code",
        "outputId": "47f5a676-a90d-4663-c812-f1402f97cc86",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 282
        }
      },
      "source": [
        "# 推定された治療効果を各人ごとに求めます\n",
        "t_estimated = reg.predict(X_1)-reg.predict(X_0)\n",
        "plt.scatter(df[[\"x\"]], t_estimated,\n",
        "            label=\"estimated_treatment-effect\")\n",
        "\n",
        "# 正解を描画\n",
        "plt.plot(x_index, t_ans, color='black', ls='--', label='Baseline')\n"
      ],
      "execution_count": 13,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "[<matplotlib.lines.Line2D at 0x7f6b85aff5f8>]"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 13
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3dfZzVZZ3/8ddnzhyGAZEbEYXhVgQKFwVFwajFNQ2TwslMIbRsU2vL9lcaCauPRH8UbFiZa7/adMs1WLxbm0WhJTdyKQKSGnREHblRBw4ixjCAzMDcXb8/zveMZ87NnJs5t8P7+XjMg3O+1/fmc77n8DnXub7XdX3NOYeIiBS/knwHICIimaGELiLSQyihi4j0EEroIiI9hBK6iEgPUZqvAw8ePNiNHj06X4cXESlKf/7zn//qnDs9VlneEvro0aPZunVrvg4vIlKUzOyteGVqchER6SGU0EVEeggldBGRHkIJXUSkh1BCFxHpIfLWy0VE5GRTVR1g+bpa9jU0MWxAOQtmTaBySkXG9q+ELiKSA1XVARY9XUNTSxsAgYYmFj1dA5CxpK4mFxGRHFi+rrYjmYc0tbSxfF1txo6hhC4ikgP7GppSWp4OJXQRkRwYNqA8peXpUEIXEcmBBbMmUO73dVpW7vexYNaEjB1DF0VFRHIgdOFTvVxERHqAyikVGU3gkdTkIiLSQyihi4j0EGpyERHJgWyPEgUldBGRrMvFKFFQk4uISNblYpQoKKGLiGRdLkaJghK6iEjW5WKUKCihi4hkXS5GiYIuioqIZF0uRomCErqISE5ke5QoKKGLiKRk/kOb2LirvtOy66ePZEnlJO6qqmHF5rqobULl2WbOuawfJJapU6e6rVu35uXYIiLpiJXM0zFj7CBW3nxxWtua2Z+dc1NjlSWsoZvZz4FPAAecc38To9yAHwFXAo3Ajc65v6QVqYhIAamqDvCtp16kuS2zFd+Nu+qZ/9CmtJN6PMn0cnkEuKKL8o8D47y/W4CfdD8sEZH8qqoO8I3Ht2U8mYdkoqYfKWEN3Tm3wcxGd7HKVcCjLth2s9nMBpjZUOfc2xmKUUQK2IEDBxgyZAgA+/bt4/Dhw53KfT4f48ePB2Dv3r0cPXq0U7nf7+fss88GoK6ujmPHjnUqLysr46yzzgLgjTfe4Pjx453Ky8vLGT16NAC7du2iubm5U3nfvn0ZOXIkADt27KC1tbVTeb9+/Rg+fDgAtbW1tLe3A3DPLzfTfOQEVtaH0n6nAdDy1z1Rr7+k9yn4ThmIc+20HgxEl5f3w9d3AK69DdfWQom/d9Q6GeOcS/gHjAZejlP2LPDhsOe/BaYm2ucFF1zgRKS4HTp0yPXr188dOHDAOefc5z73OQd0+hs8eHDH+p/61KeiykePHt1Rftlll0WVT5o0qaN82rRpUeUf+tCHOsonTpwYVT5r1qyO8pEjR0aVX3PNNR3lAwcOjCrv+zeXuVF3POtG3fGso8QXVd7vgk+6UXc860be/nRUGeBOvfhaN+qOZ93wr610Z3x2Wce+Rt3xbFrnHNjq4uTVnPZyMbNbCDbLdHxjikjxOnz4MEePHuXZZ5/lC1/4Al/5yle48sorO61TVlbW8fi2227juuuu61Tet2/fjsd33nknN910U6fy/v37dzxesmQJBw8e7FQ+ePDgjsf33XcfR44c6VR+5plndjx+8MEHaWxs7FQeqp0DXHrLt9nw2tuE9xUp7X/G+8f6xDcJ5umw8oHDgg9KfAye8y0i+U8bESzuVY5/0PvHmjF2UNS63ZVULxevyeVZF/ui6L8CzzvnVnnPa4FLXIImF/VyESl+u3bt4uyzz+bRRx/lhhtuyHc43ZKpHiwhRmTqf1/eerkkYTVwq5k9BkwDDidK5iLSM4Tao0tLi3tIy11VNSklc58Zbc5RkaURn+lKptviKuASYLCZ7QXuBvwAzrmfAmsJdlncSbDb4heyFayIFJZiT+hV1QFue3wb7Sls8+ay2VmLp7uS6eUyL0G5A76asYhEpGj4/X7OO+88Bg4cmO9QgGBNe+Xmuk5NHfFGaVZVB/j649tyF1wOFOfXqogUhPHjx7NtW/6TYlV1gDt/VcOx5raostBQ/Mikns7NJbJxITOTNH2uiBS10O3dYiXzkJUx5ldJ9eYS3bmQmSuqoYtI2rZt28Ytt9zCgw8+yEUXXZSXGGLd3i2SI5j4wy9eDhtQTqCLpN7HX8Ir//fjmQozJ1RDF5G0HT58mBdeeCFqdGcuJVvTvvNXNZ2eJ7q5xHevPjftmPJFCV1E0lYIvVySvY3bseY25j+0qeN55ZQK7r9uMv6ILOgvgfuvm1wwXRFToSYXEUlbIST0BbMmsOjpmoTNLhCcECu86SUXN53IJSV0EUlbPhN6vJtJJLJ8XW2PSuLhlNBFJG39+/dnxowZnHrqqTk9brrJHIJt7lXVgazf3zMfdMciESk6YxetpS3N3OUziDXFea5uE9ddXc3loouiIlJ0EiXzcr+vi21jL1+xuY6q6uj5zIuJErqIpO2ZZ55h4sSJvPXWWzk9rs+sy/KlV0+iIsneL+HSGT1aSJTQRSRt9fX1vPrqqx13+cmVedNGxC2bMXYQlVMq2LjwUrpO+9FSHT1aaJTQRSRt+erlsqRyEtdPj75JTuTw/GT7qKe7fqFRLxcRSVs+uy0uqZyU8CLmglkTUppRMdHo0UKnGrqIpK2lpQUo3vnQQ0qseEeHhivud0FE8mrEiBFcccUV9O6dxTvZd8Pi1dtjLi8rLaF2SXFNvJUMJXQRSdtVV13FVVddle8wYqqqDtDQ1BKzrLk1txdxc0UJXUSKTjIjPbvqgljsFz/jURu6iKTtRz/6ESNGjOhoS8+F0A0tAg1NOCDQ0MSip2uiBgV1Ndd5sV/8jEcJXUTS1tDQwN69e/H54o/MzLRYN7RoamlLaVBQsV/8jEcJXUTS1tLSQklJCSUluUsl8WrexT4oKBOU0EUkba2trTntslhVHYg7+jOyXTze9ACJpg0oZkroIpK21tZW/H5/zo63fF0tsebWMqLbxeNND9DVtAHFTr1cRCRt5557Lp/+9Kdzdrx4zSqO6Hbx0CjSVVv20OYcPjPmTRtRFFPkpiup+dDN7ArgR4APeNg5tyyifBTwc+B0oB643jm3t6t9aj50kZPXXVU1aSXaGcvWx2xDrxhQzsaFl2Yj1ILT1XzoCWvoZuYDfgxcDuwFXjCz1c65V8JWuw941Dn372Z2KbAUuKH7oYtITxDsavgSTS3RA3ranGPF5rqk7kAUeUNnCM593lO7IaYqmTb0i4Cdzrndzrlm4DEgcmjYRGC99/h3McpFpAf62te+xsSJE7tcp6o6wG2Pb4uZzFMVuYuBffwsvXpSj+2GmKpk2tArgD1hz/cC0yLWeRG4mmCzzKeAfmZ2mnPuYPhKZnYLcAvAyJHRU1+KSHE5evQox44di1lWVR3gnme2c6gxe4OO+vQqVTIPk6leLt8EZppZNTATCABtkSs5537mnJvqnJt6+umnZ+jQIpIv8botVlUHuP3JF7OazEF9zyMlU0MPAOH9fIZ7yzo45/YRrKFjZqcAn3bONWQqSBEpTPG6Ld7zzHba2rN/A/qeOidLupKpob8AjDOzMWbWC5gLrA5fwcwGm1loX4sI9ngRkR4uXg092zVz0MXQWBLW0J1zrWZ2K7COYLfFnzvntpvZvcBW59xq4BJgqZk5YAPw1SzGLCIF4tJLL+Wcc87J+XEHlPtZPOcctZ9HSKofejaoH7pIzzX5nt/EnYs8pJfPeP07VybcVzJT5Z5MutUPXUQkVYvnnMOCJ1+kJU47ugHfu+a8pPZVOaXipE7gqVBCF5G0XX755Zw4cYINGzZ0Wh5KwKGadf9yP2bQ0NiiWnYWKaGLSNpC0+fGopp17mm2RRFJW66nz5WuKaGLSNqU0AuLErqIpE0JvbDonRCRtM2bN4/+/fvnOwzxKKGLSNpuv/32qGXzH9rExl31Hc9njB3EypsvzmVYJy0ldBFJW2NjI2tefof717/BvoYmSkuip7jduKue+Q9tUlLPASV0EUnb2AkTea//GE77xDeB6GQeEl5jl+zRRVERSdvB947jzJfvMMSjhC4iaWtvbcVKlNALhRK6iKTNuTaIM1I03Iyxg3IQjSihi0jarL09YQ29l890QTRHlNBFJG3zbvoKfc6KOZNrh2RnVZTuU0IXkbT98oGl/PSfbmZAefRt6HwG9183WRN05ZC6LYpI2vbv389l4/pTeffH8h2KoBq6iHTD8OHDWbZsWb7DEI8SuoikxTlHW1ubJucqIEroIpKWtrY2ACX0AqKELiJpaW1tBZTQC4kSuoikRQm98Cihi0ha/H4/y5cvZ+bMmfkORTxJJXQzu8LMas1sp5ktjFE+0sx+Z2bVZvaSmV2Z+VBFpJCUlZXxzW9+kwsvvDDfoYgnYUI3Mx/wY+DjwERgnplNjFjtLuAJ59wUYC7w/zIdqIgUltbWVnbs2MGRI0fyHYp4kqmhXwTsdM7tds41A48BV0Ws44BTvcf9gX2ZC1FECtH+/fsZP348TzzxRL5DEU8yCb0C2BP2fK+3LNxi4Hoz2wusBb4Wa0dmdouZbTWzre+++24a4YpIodBF0cKTqYui84BHnHPDgSuBX5pZ1L6dcz9zzk11zk09/fTTM3RoEckHJfTCk0xCDwAjwp4P95aF+yLwBIBzbhPQGxiciQBFpDApoReeZBL6C8A4MxtjZr0IXvRcHbFOHfBRADP7IMGErjYVkR5MCb3wJEzozrlW4FZgHfAqwd4s283sXjOb4612O3Czmb0IrAJudM65bAUtIvk3dOhQfvrTn3L++efnOxTxWL7y7tSpU93WrVvzcmwRkWJlZn92zsW8q4hGiopIWt577z3+8pe/cPTo0XyHIh4ldBFJy7Zt27jgggvYsmVLvkMRjxK6iKRFF0ULjxK6iKRFCb3wKKGLSFqU0AuPErqIpEUJvfAooYtIWiZPnsyKFSsYM2ZMvkMRj75aRSQtw4cPZ/78+fkOQ8Kohi4iaXnnnXfYsGEDTU1N+Q5FPEroIpKW3/zmN8ycOZNAIHKuPskXNbmIZMldVTWs2rKHNufwmTFv2giWVE7Kd1gZ09bWBuiiaCHROyGSBXdV1bBic13H8zbnOp73lKSuXi6FR00uIlmwasuemMv/Y0tdzOXFZv5Dm7jjyWoALl72PPMf2pTniASU0EWyoi3OLKbtDqqqi7vNef5Dm9i4qx7XHmxyoaSEjbvqldQLgH4riWSBzyxuUl++rpbKKZG35S0eG3fVA1A+9kJKTx1CSVmfTsslf1RDF8mCedNGxC3b19Azuvn5B5xJn3HTMJ8/36GIRwldJAuWVE6i3B/7v9ewAeU5jiY7WuoDNO56Aefa8x2KeNTkIhImVlfDqaMGsejpl2hq6Zy4YnVFrKoO8K2nXqS5Lf6dwAINTYxeuIaBffzc/clziq75ZcbYQWzcVc+xVzdw+A8rGbngv8CCyyW/dAs6EU9kV8NkjRvSl8bmdgLdaEoptuQ+/6FNrHnkXzj8x1WM/NYzfPjs01h588X5Duuk0NUt6FRDF/HE62qYyI4Dx7p97EONLSx46kWAokjqK2++mDvffJbv/amUt/75E/kORzxqQxfxxOuVkistbY7Fq7fnNYZUtLa24vP58h2GhFFCFykgDU0tRdNPvbW1VaNEC4za0EU8oxeuyXcInRT6/C87d+4kEAgwc+bMfIdyUumqDT2pGrqZXWFmtWa208wWxij/oZlt8/5eN7OG7gYtkmsVBdadMDT/y11VNfkOJaazzz5bybzAJEzoZuYDfgx8HJgIzDOzieHrOOe+4Zyb7JybDPwL8HQ2ghXJpgWzJlDuz1yb8IByP5aB/axMo+dNLmzZsoWqqqp8hyFhkmkAuwjY6ZzbDWBmjwFXAa/EWX8ecHdmwhPJncopFTjn+MaDT9HU2Lnnin/QcEpPHUz7iUYGN9URaDjeqbzX4FH4ThlIacsxbhjXzoyzB3slZWzc+VfW7O3Fu61lnF56gtnDm/nlpjc5eqLt/e3PGIuvvB9t7x2i+a9vRcU2/lvH+d686Uw93fHaa69FlV988cX07duXt956ix07dkSVf+QjH6GsrIzdu3eze/fuqPJLLrmE0tJSXn/9derqor9APvrRj2JmvPrqqx3zn3//+9/n5ZdfprKyMmp9yRPnXJd/wDXAw2HPbwAejLPuKOBtwBen/BZgK7B15MiRTqTQ3PSDJx0Q9Tdo1lfdWQvXuL+/77GY5YM/+U332Z/90T3//PMxy1evXu2cc+6ZZ56JWX7G3O+6UXc86wZ/ckHM8jM/90M36o5n3T/cuSxm+Wuvveacc+6+++6LWR4IBJxzzi1evDhm+ZEjR5xzzt1+++0xy0d+6xk36o5n3SnnXdFp+ZBR4/LzRp3EgK0uTr5OeFHUzK4BrnDO3eQ9vwGY5py7Nca6dwDDnXNf63Kn6KKoFJ67qmp49Pc7aN6/A9dyAuvVG4BZE89g2Rc+xtChQ3nvvffYtm1b1Lbjx49nyJAhHD58mJqa6DbvD37wg5x22mnU19fzyivBH7cbXn+XB9YHa9O9Bo+ipPcptDUepqV+b9T2vU4fQ0lZH/wnGvjFp0dFlZ9//vn06dOHQCDAG2+8EVV+4YUXUlZWRl1dXaca+IbX32XVn+o4euoYKgadwuf+ppxzBgR/Odyz+mVq9h0BoKxiImZGy6F9tB071LF96YCh3HjZ5IK9cNsTdXVRNJmEfjGw2Dk3y3u+CMA5tzTGutXAV51zf0wUlBK6FJqxi9bG7IvuM2PX0iuzcsx0Rqe+uWx2Ro5dVR1g0dM1NLW83/RT7vex9OpJVE6pSLrXTzbPj0Trbi+XF4BxZjbGzHoBc4HVMQ7yAWAgoEmRpSi1OUfLoX0c+dOvOtVCszngaEnlJK6fPpKSTFw9TdHydbWdkjlAU0sby9fVprSffA/IkvclTOjOuVbgVmAd8CrwhHNuu5nda2ZzwladCzzmElX5RQqUz4yWd9/k0O/+rVNC91l2s+2SyknsXjqb+6+bnLCXTSYTf7xpfFOd3jfb50eSl9QwL+fcWmBtxLJvRzxfnLmwRHJv3rQR/Ov2FgCsxN9peS6E5nBZvq427kRfn502MmpZVXWAe57ZzqHGlo5lffwllPl9NDS2MGxAOQtmTYiaI2bYgPKYx3GkNsgqV+dHEtPQfxHPkspJXDymf/CJrxSfGddPH5nTC36VUyrYuPBS3lw2m+unj+yo/caLpao6wIKnXuyUzAEaW9o51NiCIzhd79cf38bohWsYs2hNx0Cl7va7NyPn50e6pqH/ImEefvhhbr75Zurq6hgxIv81z67mV/cZdDHtekJ9/CVcfcHwlC7KGvBGhi7KSno0fa5IkpqbmwEoKyvLcyTBZP6Nx7cRL2d3J5lDsBafag+bnnK3pZ5KTS4iYT7/+c+zZ88eTjvttHyHwvJ1tXGTeT4YwWYaKVyqoYuE6du3L3379s13GEDh3Ux6/vSRRXHzjZOZaugiYZ5//nnuvfde2tvzf+PjQmne8JfA/ddpNGgxUEIXCbN+/XruvvturAD6Vi+YNSEjszV2x5vLZrPju7NVMy8SanIRCdPc3ExZWVlBJPRQEo3XyyWWAeV+Fs+JfbPpu6pqWLmljmQ7tmnAUPFRQhcJ09zcTK9evfIdRofKKRUdyXnGsvVxBxzdf93khLXoJZWTOppNxixck/CCqwYMFR81uYiEKbSEHi7WQKByvy+pZB4pmUq62syLj2roImEKOaGHTw2wr6Ep7pD+ZPjMupxUS80txUkjRUXCnDhxghMnTnDqqafmO5SsSjRtr4b0Fy6NFBVJUllZWUGMEs22ULJetWVPp5q6z4x500YomRcp1dBFwjzyyCPU19dz22235TsUkZi6e4MLkZPG008/zYoVK/IdhkhalNBFwhTyRVGRRJTQRcIooUsxU0IXCaOELsVMCV0kTEtLixK6FC11WxQJs2XLloKYaVEkHaqhi0QoKdF/CylO+uSKhFm0aBG/+MUv8h2GSFqU0EXCrFy5kt///vf5DkMkLUroImFOnDihi6JStJJK6GZ2hZnVmtlOM1sYZ51rzewVM9tuZv+R2TBFckPdFqWYJezlYmY+4MfA5cBe4AUzW+2ceyVsnXHAImCGc+6QmQ3JVsAi2aSELsUsmRr6RcBO59xu51wz8BhwVcQ6NwM/ds4dAnDOHchsmCK50bt3b/r27ZvvMETSkkw/9ApgT9jzvcC0iHXGA5jZRsAHLHbO/XfkjszsFuAWgJEjR6YTr0hWHTx4MN8hiKQtUxdFS4FxwCXAPOAhMxsQuZJz7mfOuanOuamnn356hg4tIiKQXEIPAOF3ix3uLQu3F1jtnGtxzr0BvE4wwYsUjcbGRubOncvatWvzHYpIWpJJ6C8A48xsjJn1AuYCqyPWqSJYO8fMBhNsgtmdwThFsq6pqYnHH3+cnTt35jsUkbQkTOjOuVbgVmAd8CrwhHNuu5nda2ZzvNXWAQfN7BXgd8AC55waI6WoNDc3A5wUt6CTnimpybmcc2uBtRHLvh322AG3eX8iRSmU0NVtUYqVRoqKeJTQpdgpoYt4nHMMGzaMfv365TsUkbRoPnQRz/jx4wkEIjtwiRQP1dBFRHoIJXQRT01NDbNnz+all17KdygiaVFCF/EcOHCAtWvXcvjw4XyHIpIWJXQRj3q5SLFTQhfxKKFLsVNCF/EooUuxU0IX8fTp04fx48fTp0+ffIcikhYLjtrPvalTp7qtW7fm5dgiIsXKzP7snJsaq0w1dBGRHkIJPQuqqgPMWLaeMQvXMGPZeqqqNfqwGDz55JPMnDmTI0eO5DsUkbRo6H+GzX9oExt31Xc8DzQ08fXHt/H1x7cBcP30kSypnJSv8KQLdXV1bNiwATPLdygiaVFCz5DIRB7Pis11AErqBUi9XKTYqcklA5JN5iErNtdx1qI13FVVk8WoJFWhhO73+/MciUh6VEPPgFSSeUi7Cyb2UI09ZMbYQay8+WKqqgP809Mv0djSDoAZzJ+m5ppsam5uprS0lJIS1XOkOCmhd0NVdYDbn9iW0X1u3FXP6IVropY7p+aabDvzzDO56KKL8h2GSNrUDz1Fl//geXYcOJa34/vM2LX0yrwdP5vuqqrhP7bU0e59JMv9JSy9+lwqp1TkNzCRAtJVP/QeX0Of9p3neOdoc5frjBvSl+duuyQj+8q2tjx9AWfbXVU1Uc1PTS3tnXoIqdlJpGs9trGwqjrA6IVrkkrAOw4c4/IfPJ9wf6km8zeXzeb66SNT2iYZPfFi6qotexKuE2p2ytbrv/vuu7nyyp7560dODj0yoVdVB1jw1IspbZOoGWX5utqU9hdK5EsqJ3H99JFksmvzis11nPPt/+5RA5ZS+eWRraS+e/duamtTe59FCkmPbHJZvq6WlrbMNk3sa2hKet1xQ/p2ahZYUjmp43lVdYDl62oJNDRRYnS0F6fqWHNbR3NET2hr9pmlnNQhsxeIm5ub1QddilpSNXQzu8LMas1sp5ktjFF+o5m9a2bbvL+bMh9q8gIpJN9kDRtQnnAdI1gz76o9vnJKBRsXXsqby2aze+ls3lw2m/uvm9ytn0qhtuZiboqZN21Eytsk00yTCiV0KXYJ84iZ+YAfAx8HJgLzzGxijFUfd85N9v4eznCcSUs3qY0b0rfL8gWzJlDu98UsG9jHz/3XTeaNZbPTqjFWTqlg97LZVCTxpdGVlZvrirYZJtQ0lYpMXyBWQpdil0yTy0XATufcbgAzewy4Cnglm4F15bvf/S5/+ctfOi0bNmwYDzzwAKu27KFhwy9pOdi59lY6cCgDL/kCAIfW/xuth9/pKDuldykzPvV3wCUA3Hrrrezfv7/T9hdeeCFLr76e5etqeWnVP9OrvYkJZ/TrqLm/1f8jMOX/ADB//nxOnDjRafvLL7+cL33pS7S3t3PttddGvaY5c+awYNZHueOxF9hbtfz9AgMc9Png39L3Ax+mrekI9f/9YNT2fSddRp+zL2LJE39gxXeejCr/8pe/zGWXXcbOnTtZuDDqRxZf//rX+fCHP0xNTQ333HNPVPmiRYu44IIL+NOf/sT3vve9qPJ7772XiRMnsmHDBh544IGo8uXLlzNmzBg+cONS9mx6Nqp84OVfZtTwCuad8Q7/u+Y/qX3nKMeb2zrKB338H/H1PoVjr/wvjbUbO5Z/cMP97Klv4rRP3E6pv4yJ7/2Zsrc7Xz8xM558MnhOfvKTn/Db3/62U3mfPn149NFHmTRpEkePHo2KTaRYJJPQK4Dw7LgXmBZjvU+b2d8CrwPfcM5F/R42s1uAWwBGjky/98fevXt57bXXOi0LJdA252g9+i4t9Xs7b1Ti4/7rJlM5pYJrrvkFr70WdiPgtuDETCG7d+/u9Bxg6NCh3DGlgsopFcz89T0cPFjPkf31HPHy/tixYzvWff3112lq6tzsc+6553Y8jowdYPr06XxuSgVNx97ji/++j9a2dkp9JQw+pYzjza20NXkzALa3Rb82oP14MBHtP/Qex2Psv6GhAYDjx4/HPH4okTU2NnYqP3DkOA1NLfyh13P0HrmflroXqd+wlTbnOuI7tXcpjY2NABw5ciTm/o8fP8607zzHwUMNMeOnrZVAQxOPbH8Z/5u7GNq7FPzvH5/2YHJvazrSafvQIF3n2mlzjj/U7Ma/extDTu3dsU74ZFv79++Piu+UU04BYNmyZdFxiRSRhAOLzOwa4Arn3E3e8xuAac65W8PWOQ14zzl3wsy+BFznnLu0q/1mY2BRVXWgo89ypBKD3UtnZ/R4uVRVHeC2J7YlvIhaMaCcjQu7PPUx97149fZg4iTYhDT4lF5JDaAq9/tYevWkpC7IxhoBG0vka4jVRz0RzWopPVV3b3ARAMKvWA33lnVwzh10zoXaGB4GLkgn0O6Y/9CmuMkcoKy0uHtoVk6p4AfXTmZgn/gTR5X7fSyYNSGl/VZVB1jw5IsdyRzgUGNL0qNhm1rakurSmUrbfmSPoiWVk0i11+eKzXVM+85zKW4lUtySaXJ5ARhnZmMIJvgAM/EAAAsNSURBVPK5wGfDVzCzoc65t72nc4BXMxplAsnMdnjcm+Qqllg1wNAkWYWk0mvyCQnvAukz65Rck+3CuHj1dlrS7TvpSdSls6o6wG1dfNlGitWjaNiA8pR7L71ztJn5D20quPdRJFsSJnTnXKuZ3QqsA3zAz51z283sXmCrc2418I9mNgdoBeqBG7MYc5RkZjuMTBKRzQyx9lnoySCUtBc9XUNTS7CNOdDQxKKnazqVx1NVHYj7+lMRfm4j52NJVbxfGQtmTej0OpOVzkyYIsUqqYFFzrm1wNqIZd8Oe7wIWJTZ0DIrPEmEmhkS1UyLIRksX1cbleSaWt4fdDRuSF/2HmqiyfuFUmLwWW8+lFRHv8YSnoDTaesOVzGgnAWzJsT8Inr/y+uljteSKZqqWHqKoh8pmky/84F9/J2SxPJ1tUk3M4xdtJZ500YU7H/uRM0dkW3h7WHT8KYy+jXEZ9Cvt5/DTS0MC0vAVdWBtJJ5KrNHhpqc7qqqYdWWPRnphx7rYrOmKpZiVdQJvao6wMoESaTc7+PuT57TaVkqiazNuYL+z92/3J9Ws8mqLXtSapf2l8Dyz0yOWXtOZ+6ckHRGiIZPpQCJr6HcVVUT971bvq42bvPQqi17CvI9F4mn6BJ6VXWAO39Vw7HmxG2p8X7Cp3OBbcXmOqaOGlRw86Wk2qYc0uZcwnbpAeV+Fs85p9NrDq8d+8yYftbAtJumjMx8SYauc8RL7F19IXf15d5TpyqWnquoEnpVdYDbn3yRtiSaS7rqj71g1oSk2tCjtvNqoYWS1O+qquFEa3rtyT5vsE1ZaUlHQh/Yx8/dnzwn7uuLbCNvc65b1xnmZ3hq4ZU3X8zYRWtjJuJ4te2uvtx9mZwiUyQHiiqhL19Xm1QyB7rsjx1KWJGDae7+5Dnc/sSLcWtmLW2O5etqCyahd2dyqulnDYyqnXfVtbO7xwvnM8vadYl4712bc0y469edvgD7+Eu4+oLhcXvlpNMcJJJPRZXQU2n7TpR0I/t0h2x9q77Li3vpXEjMlmSaBCJ7uYT8cVc9kVuH+rHHO3epNkH4S4zlnzkvp1+AXU3DG/lrprGlnRWb6+jlM5ojplvWSFMpRkU1fDKZKWyBbs1auKRyEn17xZ5VEaDErGBmNIzXJGAE75b05rLZPHfbJSy9+lz8vs7rxkvNXX1hpdoEketkDunVqiOTud9nTB01KFMhieRMUSX0BbMm4CvpOqmUWNfNLcn4zqcm4Y9znDbnWPR0TUEk9XjJK7JtOpUbfnT1pZlKshw3pG9emqYyUasONa2JFJuiSuiVUyr4/mfOi1uDList4QfXxu5al+pxln/mPAaUx543Jdn5S7ItNId4qObsM4vZVJBKM1FXX4bJzll+apkvqZtuZ0t355WHwmpaE0lWwtkWsyUbsy1mw5iFa2I2TxjwxrLimL1xxrL1SXfTfDPJ13T5D56POYFXIbQ9h/rFd+c2hOnMWimSC13NtlhUF0XzIV63tmTb8wvBglkTup3gIuWzBp5I6BfaPc9s51Bj6oOu/D7rdrOdSD4ooScQa/BNOtPU5lOyCa4n9bqO1YupqjqQ8Bwk6osvUsiU0BMI/cdevq6WfQ1NneYvKSbhCS5ec0mmB/oUmtA5iDeitBCai0S6Q23oJ6nIIfyFPAFZNpzsr1+KV1dt6EroIiJFpLu3oBMRkSKghC4i0kMooYuI9BBK6CIiPYQSuohID5G3Xi5m9i7wVpqbDwb+msFwMkVxpaZQ44LCjU1xpaYnxjXKOXd6rIK8JfTuMLOt8brt5JPiSk2hxgWFG5viSs3JFpeaXEREeggldBGRHqJYE/rP8h1AHIorNYUaFxRubIorNSdVXEXZhi4iItGKtYYuIiIRlNBFRHqIgk3oZvYZM9tuZu1mFrd7j5ldYWa1ZrbTzBaGLR9jZlu85Y+bWa8MxTXIzJ4zsx3evwNjrPN3ZrYt7O+4mVV6ZY+Y2RthZZNzFZe3XlvYsVeHLc/n+ZpsZpu89/slM7surCyj5yve5yWsvMx7/Tu98zE6rGyRt7zWzGZ1J4404rrNzF7xzs9vzWxUWFnM9zRHcd1oZu+GHf+msLLPe+/7DjP7fI7j+mFYTK+bWUNYWTbP18/N7ICZvRyn3MzsAS/ul8zs/LCy7p8v51xB/gEfBCYAzwNT46zjA3YBZwG9gBeBiV7ZE8Bc7/FPgX/IUFzfAxZ6jxcC/5xg/UFAPdDHe/4IcE0WzldScQHvxVmet/MFjAfGeY+HAW8DAzJ9vrr6vISt8xXgp97jucDj3uOJ3vplwBhvP74cxvV3YZ+hfwjF1dV7mqO4bgQejLHtIGC39+9A7/HAXMUVsf7XgJ9n+3x5+/5b4Hzg5TjlVwK/JniDsOnAlkyer4KtoTvnXnXO1SZY7SJgp3Nut3OuGXgMuMrMDLgUeMpb79+BygyFdpW3v2T3ew3wa+dcY4aOH0+qcXXI9/lyzr3unNvhPd4HHABijoTrppifly7ifQr4qHd+rgIec86dcM69Aez09peTuJxzvwv7DG0Ghmfo2N2KqwuzgOecc/XOuUPAc8AVeYprHrAqQ8fuknNuA8EKXDxXAY+6oM3AADMbSobOV8Em9CRVAHvCnu/1lp0GNDjnWiOWZ8IZzrm3vcf7gTMSrD+X6A/Td7yfWz80s7Icx9XbzLaa2eZQMxAFdL7M7CKCta5dYYszdb7ifV5iruOdj8MEz08y22YzrnBfJFjLC4n1nuYyrk97789TZjYixW2zGRde09QYYH3Y4mydr2TEiz0j5yuv9xQ1s/8BzoxRdKdz7r9yHU9IV3GFP3HOOTOL2+/T++adBKwLW7yIYGLrRbAv6h3AvTmMa5RzLmBmZwHrzayGYNJKW4bP1y+Bzzvn2r3FaZ+vnsjMrgemAjPDFke9p865XbH3kHHPAKuccyfM7EsEf91cmqNjJ2Mu8JRzri1sWT7PV1blNaE75y7r5i4CwIiw58O9ZQcJ/pQp9WpZoeXdjsvM3jGzoc65t70EdKCLXV0L/Mo513Gb+bDa6gkz+wXwzVzG5ZwLeP/uNrPngSnAf5Ln82VmpwJrCH6Zbw7bd9rnK4Z4n5dY6+w1s1KgP8HPUzLbZjMuzOwygl+SM51zJ0LL47ynmUhQCeNyzh0Me/owwWsmoW0vidj2+QzElFRcYeYCXw1fkMXzlYx4sWfkfBV7k8sLwDgL9tDoRfDNW+2CVxl+R7D9GuDzQKZq/Ku9/SWz36i2Oy+phdqtK4GYV8OzEZeZDQw1WZjZYGAG8Eq+z5f33v2KYNviUxFlmTxfMT8vXcR7DbDeOz+rgbkW7AUzBhgH/KkbsaQUl5lNAf4VmOOcOxC2POZ7msO4hoY9nQO86j1eB3zMi28g8DE6/1LNalxebB8geIFxU9iybJ6vZKwGPuf1dpkOHPYqLZk5X9m62tvdP+BTBNuRTgDvAOu85cOAtWHrXQm8TvAb9s6w5WcR/A+3E3gSKMtQXKcBvwV2AP8DDPKWTwUeDltvNMFv3ZKI7dcDNQQT0wrglFzFBXzIO/aL3r9fLITzBVwPtADbwv4mZ+N8xfq8EGzCmeM97u29/p3e+TgrbNs7ve1qgY9n+POeKK7/8f4fhM7P6kTvaY7iWgps947/O+ADYdv+vXcedwJfyGVc3vPFwLKI7bJ9vlYR7KXVQjB/fRH4MvBlr9yAH3tx1xDWgy8T50tD/0VEeohib3IRERGPErqISA+hhC4i0kMooYuI9BBK6CIiPYQSuohID6GELiLSQ/x/G4aHXyA4sfQAAAAASUVORK5CYII=\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "gCPyyuXWb_ym"
      },
      "source": [
        "## X-Learner"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Aq7lb2fkcDDU",
        "colab_type": "code",
        "outputId": "0b0b0fdf-6d2d-4d17-8be5-961d3e8ef00e",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 229
        }
      },
      "source": [
        "# T-LearnerでM0とM1を求める\n",
        "from sklearn.ensemble import RandomForestRegressor\n",
        "\n",
        "# 集団を2つに分ける\n",
        "df_0 = df[df.Z == 0.0]  # 介入を受けていない集団\n",
        "df_1 = df[df.Z == 1.0]  # 介入を受けた集団\n",
        "\n",
        "# 介入を受けていないモデル\n",
        "M0 = RandomForestRegressor(max_depth=3)\n",
        "M0.fit(df_0[[\"x\"]], df_0[[\"Y\"]])\n",
        "\n",
        "# 介入を受けたモデル\n",
        "M1 = RandomForestRegressor(max_depth=3)\n",
        "M1.fit(df_1[[\"x\"]], df_1[[\"Y\"]])\n"
      ],
      "execution_count": 14,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "/usr/local/lib/python3.6/dist-packages/ipykernel_launcher.py:10: DataConversionWarning: A column-vector y was passed when a 1d array was expected. Please change the shape of y to (n_samples,), for example using ravel().\n",
            "  # Remove the CWD from sys.path while we load stuff.\n",
            "/usr/local/lib/python3.6/dist-packages/ipykernel_launcher.py:14: DataConversionWarning: A column-vector y was passed when a 1d array was expected. Please change the shape of y to (n_samples,), for example using ravel().\n",
            "  \n"
          ],
          "name": "stderr"
        },
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "RandomForestRegressor(bootstrap=True, ccp_alpha=0.0, criterion='mse',\n",
              "                      max_depth=3, max_features='auto', max_leaf_nodes=None,\n",
              "                      max_samples=None, min_impurity_decrease=0.0,\n",
              "                      min_impurity_split=None, min_samples_leaf=1,\n",
              "                      min_samples_split=2, min_weight_fraction_leaf=0.0,\n",
              "                      n_estimators=100, n_jobs=None, oob_score=False,\n",
              "                      random_state=None, verbose=0, warm_start=False)"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 14
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "KPpIUdNtQK5-",
        "colab_type": "code",
        "outputId": "eae9b81e-a3c6-46a1-eb59-4a480fee3987",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 139
        }
      },
      "source": [
        "# 推定された治療効果を各人ごとに求めます\n",
        "tau_0 = reg_1.predict(df_0[[\"x\"]]) - df_0[\"Y\"]\n",
        "tau_1 = df_1[\"Y\"] - reg_0.predict(df_1[[\"x\"]])\n",
        "\n",
        "# ATTとATUを求めるモデルを作成します\n",
        "M2 = RandomForestRegressor(max_depth=3)\n",
        "M2.fit(df_0[[\"x\"]], tau_0)\n",
        "\n",
        "M3 = RandomForestRegressor(max_depth=3)\n",
        "M3.fit(df_1[[\"x\"]], tau_1)\n"
      ],
      "execution_count": 15,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "RandomForestRegressor(bootstrap=True, ccp_alpha=0.0, criterion='mse',\n",
              "                      max_depth=3, max_features='auto', max_leaf_nodes=None,\n",
              "                      max_samples=None, min_impurity_decrease=0.0,\n",
              "                      min_impurity_split=None, min_samples_leaf=1,\n",
              "                      min_samples_split=2, min_weight_fraction_leaf=0.0,\n",
              "                      n_estimators=100, n_jobs=None, oob_score=False,\n",
              "                      random_state=None, verbose=0, warm_start=False)"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 15
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Ou5_JAZQQ0QX",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# 傾向スコアを求めます\n",
        "from sklearn.linear_model import LogisticRegression\n",
        "\n",
        "# 説明変数\n",
        "X = df[[\"x\"]]\n",
        "\n",
        "# 被説明変数（目的変数）\n",
        "Z = df[\"Z\"]\n",
        "\n",
        "# 回帰の実施\n",
        "g_x = LogisticRegression().fit(X, Z)\n",
        "g_x_val = g_x.predict_proba(X)\n",
        "\n",
        "# それぞれのモデルで全データの効果を予測し、傾向スコアで調整\n",
        "tau = g_x_val[:, 1]*M2.predict(df[[\"x\"]]) + g_x_val[:, 0]*M3.predict(df[[\"x\"]])\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "MVuQlQCrfvjB",
        "colab_type": "code",
        "outputId": "66545a99-f84d-49b6-8979-7886ce5d1c93",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 282
        }
      },
      "source": [
        "# 推定された治療効果を各人ごとに求めます\n",
        "plt.scatter(df[[\"x\"]], tau, label=\"estimated_treatment-effect\")\n",
        "\n",
        "# 正解を描画\n",
        "plt.plot(x_index, t_ans, color='black', ls='--', label='Baseline')\n"
      ],
      "execution_count": 17,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "[<matplotlib.lines.Line2D at 0x7f6b85a81ba8>]"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 17
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3de3xU9Z3/8dcnFyCAEC7icgugUlFAoRsFde2yuC3WVqHWrqj09qu67VZ3lZYira3Uh7ta6a+1rfaiPtx2iwtYLymia0rx9tMKGkgAQUFAuQwqIBeLieT2+f0xJ9lJZiaZTCZzCe/n45GHmfP5nnM+nBk/+c73fM855u6IiEjuy8t0AiIikhoq6CIi3YQKuohIN6GCLiLSTaigi4h0EwWZ2vHgwYN99OjRmdq9iEhOWrt27QF3PzFWLGMFffTo0VRUVGRq9yIiOcnMdsaLachFRKSbUEEXEekmVNBFRLoJFXQRkW5CBV1EpJvI2CwXEZHjTVlliEXlW9h7uIZhxUXMm3EasyYPT9n2VdBFRNKgrDLEgsc2UlPXAEDocA0LHtsIkLKiriEXEZE0WFS+pbmYN6mpa2BR+ZaU7UMFXUQkDUKHa2Iu3xtneTJU0EVEulhZZQiLExtWXJSy/WgMXUSkA66+/2Ve2n4wavmcqSUALF69K+FtGTBvxmmpSk0FXUQkUfGKOXSskDdxUndCFBIYcjGzB81sn5m9FiduZvZzM9tmZhvM7OMpy05EJIvEK+bJGp7C4RZIbAz9t8BFbcQ/DYwNfq4DftX5tEREur9UDrdAAkMu7v6CmY1uo8lM4L/c3YHVZlZsZkPd/Z0U5SgiWWz9+vXs27ePcePGMXLkSD744APWrFkT1W7ChAkMHTqUQ4cOxbx19llnncWQIUM4cOAAlZWVUfG//du/ZeDAgbz77rts3LgxKn7OOefQv39/QqEQmzdvjoqfe+659O3bl507d7J169ao+AUXXECvXr3YsWMH27dvj4pPmzYNgLqDIeqP7IuK9xo9CTOj7sBu6v96oGXQjKLRk6LWSeVwCwDu3u4PMBp4LU5sBfB3Ea9XAaVx2l4HVAAVJSUlLiK57cCBA25mDvgvf/lLd3dft26dEx4ebvHz+9//3t3dn3/++ZjxsrIyd3d/8sknY8ZXrVrl7u5Lly6NGV+zZo27u99///0x45s3b3Z395/85Ccx47t373Z39x/+8Icx4yNvXOaj5q/wfudcFjNe8p3lPmr+Cu876dNRMSvo6aPmr2jxc9V9f0nqmAMVHqdWWzjetqCHvsLdJ8SIrQDudPcXg9ergPnu3ubTK0pLS10PuBDJbTt37mT06NEsWLCAG264gaFDh3L06FHWr18f1Xbs2LEMGTKEI0eO8Npr0afkxo0bx6BBgzh48CCvv/56VHz8+PEUFxezf//+mD3sM888kxNOOIF33303Zg978uTJ9O7dm1AoxNtvvx0VLy0tpWfPnuzatYvdu3fzwtb9/HzVmzRVyJ7DTsPy8qk/8h71f30/av2ew08P99APvUPDh4daxMyMnsNPb359/ikDeejac6O2kQgzW+vupTFjKSjovwGec/clwestwDRvZ8hFBV0k923fvp1TTz2V3/3ud3zpS1/KdDopdf6dz8S9GKgj5kwt4fZZE1OQUVhbBT0VFxYtB74UzHaZChxpr5iLSPdQX18PQEFB95sB3dkrOPPNUl7M29Puu2BmS4BpwGAz2wPcChQCuPuvgaeAi4FtQDXw1a5KVkSyy4gRI/jzn//M+PHjM51Kyg0rLkq4hz68uIiXbp7exRm1L5FZLle2E3fgmynLSERyRp8+fbjwwgsznUbKlVWGEi7mRYX5KZ9+mCzdy0VEknbw4EGWLVvG3r17M51KypRVhrhxWVWbbfItfGeW4cVF3HHZxNRPP0xS9xv4EpG02b59O7Nnz2bFihUMGzYs0+kkrawyxA+f2MSh6rp2275952fSkFFyVNBFJGlNJ0Xz8/MznElyyipDzF1WRWOmE0kRFXQRSVpDQ/iBDbk4yyWRoZVcozF0EUlark5bLKsM8a2Hoy9+ak+/ntn9TSS33gURySq5WNDbugVuW/r1zGfDD9u6T2Hm5c67ICJZ5+yzz2b16tWcfvrp7TfOAreUbexwMc8DfnLFpKyZydIWFXQRSVr//v2ZMmVKptNoU1lliPmPbuBYfcdPfRYV5nHHZWfmRDEHFXQR6YS3336bVatWMWvWLAYNGpTpdKKUVYaY+3AVje3fsqpZYR4s+kJu9MhbU0EXkaStXbuWa665hrPPPjsrC/qi8i0dKubZcgl/sjTLRUSSlu0nRTt6g61suYQ/WSroIpK0bJ+HPizBZ3bmGdydIyc+25Kd74KI5IRsuVK0rDLEovItHbp/eWceMpGtVNBFJGmZ7KHfUraRxat3JbVuz4K8blfMQUMuItIJn/vc59i0aRNDhw5N6347U8wBapOYwpgL1EMXkaQVFxdTXFyc9v0uWbO7U+snOraea9RDF5GkVVZWcvfdd1NT0/lnb3ZEQwLPQm5Lrs9miUcFXUSS9vzzz3PTTTdx7NixtO636QETyTjphB45P5slHhV0EUlapuahXzllZFLrjR3ShzXf+2SKs8keGkMXkaRlatri7bMmAkSdGO2OUxE7QgVdRJKWqR56WWWIZ9/YjxE+wTlvxmnddhilI1TQRSRpTfPQ09lDL6sMseCxjdTUhfcdOlzDgsc2Ahz3RV1j6CKStJtuuomdO3eSl5e+UrKofEtzMW9SU9fAovItacshW6mHLiJJ69evH/369UvrPuNd3t/RG3F1RyroIpK0lStXsm7dOubPn5/S7ZZVhvjuYxuorkv8is7uerFQR6igi0jSnn76ae6///4OF/Smm2ntPVwTdVIzmWd+Gt33YqGOUEEXkaTV19d3+IRoWWWIeX9YT13w5InQ4RpuXFbFjcuqks7D0QlRSPCkqJldZGZbzGybmd0cIz7KzFaZ2QYze87MRqQ+VRHJNvX19R2esrhw+abmYp4qwzXcAiTQQzezfOBe4JPAHuBVM1vu7psjmv0Y+C93/52ZTQfuAL7YFQmLSPZItKCHpxpuoKYDY+KJKirM13BLIJEe+jnANnff4e61wFJgZqs2ZwDPBL8/GyMuIt1QIkMuZZUh5i6r6pJiPqB3IXdcNlHDLYFEvisNByLvVbkHmNKqzXrgMuBnwOeAE8xskLu/H9nIzK4DrgMoKSlJNmcRyRL33HNPuzfmWlS+hVSVciM8Xj5cV4fGlKqTot8G7jGzrwAvACGgoXUjd78PuA+gtLQ0tYNoIpJ2RUVFFBW1PX7d2fnhY4f0YeXcaZ3axvEikSGXEBB5a7MRwbJm7r7X3S9z98nA94Jlh1OWpYhkpcWLF/PjH/+4zTbJzg/PN2PO1BIV8w5IpIf+KjDWzMYQLuSzgasiG5jZYOCguzcCC4AHU52oiGSfP/7xj2zevJlvf/vbcdvMm3Eac5dVxR12mTO1pPnuidI57fbQ3b0euB4oB14HHnb3TWZ2m5ldGjSbBmwxs63AScC/d1G+IpJFEpnlUrHzoIp5miQ0hu7uTwFPtVr2g4jfHwEeSW1qIpLt2pvl0tbDnMcO6aNinmK626KIJK2hoSFuD72sMhS3mAO8ue/DrkrruKWCLiJJi9dDL6sM8a0/rM9ARsc33ctFRJJWXl7e/JCLSD98YhMNKb68X9qnHrqIJM3MYg65HKqua3fdsUP6dEVKxzUVdBFJ2l133cUvfvGLDq930gk9NL+8C6igi0jSHn30UZ588smo5cVFhXHXmTO1hDXf+2RXpnXcUkEXkaTFm4e+8NLxFOZZ1HLNO+9aOikqIkmLN22x6aZZ8Z5KJF1DBV1EktbWlaKzJg9XAU8zDbmISNJ69OjR7t0WJX3UQxeRpK1bty7TKUgE9dBFRLoJFXQRSdr111/PAw88kOk0JKCCLiJJe/TRR3nllVcynYYEVNBFJGlt3W1R0k8FXUSSlsgDLiR9VNBFJGntPeBC0ksFXUSSNnjwYPr375/pNCSg70oikrQdO3ZkOgWJoB66iEg3oYIuIklxdy655BKWLl2a6VQkoIIuIklpaGhgxYoVbNu2LdOpSEBj6CKSlPr6eoCoWS5X3/8yL20/2Pz6/FMG8tC156Y1t+OVeugikpSmh0NHzkNvXcwBXtp+kKvvfzmtuR2vVNBFJCmxeuiti3l7yyW1VNBFJGmnnnoqAwcOBFAvPAskNIZuZhcBPwPygQfc/c5W8RLgd0Bx0OZmd38qxbmKSBbp378/b775JgC3lG1ULzwLtNtDN7N84F7g08AZwJVmdkarZrcAD7v7ZGA28MtUJyoi2Wvx6l1txs8/ZWCaMjm+JdJDPwfY5u47AMxsKTAT2BzRxoF+we/9gb2pTFJEss97773H7NmzYeIl0Pu0uO165JtmuaRJImPow4HdEa/3BMsiLQTmmNke4CnghlgbMrPrzKzCzCr279+fRLoiki2qq6t57rnn2Lh9T5vt7rr8rDRlJKk6KXol8Ft3HwFcDPzezKK27e73uXupu5eeeOKJKdq1iGRC0ywX8touI7Mmt+7/SVdJpKCHgJERr0cEyyJ9DXgYwN1fBnoBg1ORoIhkp6aCbnnxb587Z2pJutIREivorwJjzWyMmfUgfNJzeas2u4ALAczsdMIFXWMqIt1Ycw89+ss4ED4RevusiWnMSNot6O5eD1wPlAOvE57NssnMbjOzS4Nm3wKuNbP1wBLgK+7uXZW0iGRez549OenkceQVnRAV0+X+mWGZqrulpaVeUVGRkX2LSOrcUraRJWt20+BOvhlXThmpnnkXMrO17l4aM6aCLiKSO9oq6Lr0X0SSUlVVRWlpKa+88kqmU5GACrqIJOXIkSOsXbuWDz/8MNOpSEAFXUSS0jTLJfL2uZJZKugikpR4D7iQzFFBF5GkxHrAhWSWCrqIJKW4uJgLLriAfv36td9Y0kJ/WkUilFWGWFS+hb2HaxhWXMS8GafpXiRxnHfeebzwwguZTkMiqKCLBMoqQ8z7w3rqGsPXZoQO1zDvD+sB3WBKcoOGXEQCC5dvai7mTeoanYXLN2Uoo+y2YsUKxo0bx1tvvZXpVCSggi4SOFxT16Hlx7tDhw6xZcuW5pOjknkq6CKSFM1yyT4q6CKBAb0LO7T8eFZWGeL2JzYCcPmvV1NW2foRCZIJKugigVsvGU9hvrVYVphv3HrJ+AxllJ3KKkPMXVbFoQ+PAfDu0TrmLqtSUc8CKugigVmTh7Po8rMYXlyEAcOLi1h0+Vma4dLKgsc20AgUnDCYopNLySvoQWOwXDJLg18iEWZNHq4C3o6aukYAik45m6JTzo5aLpmjgi4SQRcWSS7TkItIoKwyxILHNhI6XIMTvrBowWMbNTbcSl5wmuGv61aw594v0Vh3rMVyyRwVdJHAovIt1NS1nFNdU9fAovItGcooO101pQSAxmPVNBw9iOXltVgumaMhF5HA3sM1HVp+vGp6XugvXwr/8cvPy+eqqSV6jmgWUA9dJDCsuCjmcif8IGT5X7fPmsgN008BYPudn1UxzxIq6CKBeTNOo6gw9sMaFq/epaLeSn19PQUFBZhp8DxbaMhFJNA0m+XGZVUx40vW7O5QT/SWso3895pdNN3vq6gwjzsuO7PbzJqZMGECV1xxRabTkAjm7u236gKlpaVeUVGRkX2LtGX0zU926fbvvmJStynqkn5mttbdS2PFNOQi0kp+Fw8h3LisitO//z+aDikppx66SCu3lG1k8epdGdn3+acM5KFrz83Ivjtq7ty5PP7447ofepq11UNPaAzdzC4CfgbkAw+4+52t4j8F/iF42RsY4u7FyacskjlN4+QPrdlFuvs7L20/2DzkM6B3IbdeMj5rh2eqq6uprq7OdBoSod0hFzPLB+4FPg2cAVxpZmdEtnH3m9x9krtPAn4BPNYVyYqky+2zJtKrILMjkoeq65j7cPbexbBplotkj0TejXOAbe6+A8DMlgIzgc1x2l8J3Jqa9EQyZ+ey26Cx5ZWjRaeewwmTPo031LH/8f+IWqf3uAvoO2E6jcc+5MATP46K95lwIX3G/R0NRw/x/tM/j4r3nXwxvU85m7rD73Loz78B4Ko/hGMD+vTg9z+9jenTp7N582bmz58ftf6CBQs477zzWLt2LQsXLoyK33bbbUyePJkXX3yRH/3oR1HxRYsWMW7cOFauXMnPfx6d3z333MOoUaN44oknKC8vJz8/9jRPyYxECvpwYHfE6z3AlFgNzWwUMAZ4Jk78OuA6gJISXSYs2ae+vp5jx45RVFREw9GDUQW98VgwxOCE4614bTju7nHi4atOvbEhTvyjYEfR8QNHYe5/r6Fq+nRqa2vZu3dv1PoffRRe/9ixYzHjx46F77tSU1MTM15XF37cXnV1dcx4fX09AEePHmXIkCFceOGFUW0kc9o9KWpmlwMXufs1wesvAlPc/foYbecDI9z9hvZ2rJOiko2eeeYZLrzwQp577jm++vRRGmP872HAW3d+JuFtllWG+O5jG6hO4e1lc+nkqaRWZ0+KhoCREa9HBMtimQ18s2PpiWSPph5qjx49uGpKSczZLldP7di3y8h7rDfdnjfUyfvDRJ48jaWrTqjq9sLZLZGC/iow1szGEC7ks4GrWjcys3HAAODllGYokka1tbUAFBYWNs92WbJmNw3u5Jtx5ZSRnbpvSawHaLS+ojQVDlXXceOyKr7zyHruStFTl5puL9x0R8qm2wsDKupZot2C7u71ZnY9UE542uKD7r7JzG4DKtx9edB0NrDUMzWxXSQFInvoEJ7t0tU3nmraR9OzOlP53J/aBmfuw+FbGXS26LZ1e2EV9OyQ0Jwjd38KeKrVsh+0er0wdWmJZEbrgp5OTUUxFUMykRodFi7flFTRvfr+l3lpe/TJ20ipzFU6R5NIRSKcccYZfP/73+fEE0/MyP5bD8l88ifP8ea+Dzu93cM1dV12j5quvlWCJE4FXSTCxIkTmTgxe+7tvXLutJQV9a7SoFHWrKGbc4lEOHr0KO+99x6NjdnzBPuVc6dx9xWTyM/SjvDwOA8GkfRTD10kwoMPPsi//du/8f777zNw4MBMp9Ms1uyYWNJ9YzEj/GAQyQ7qoYtEaDopWlhYmOFMknP7rIkM6J2+3K+eWqIZLllEBV0kQtM89EzMckmVWy8Z3+X7KC4q5O4rJulZollGQy4iEZ7esAeAj/3gTxTk5Xf6QqJMaOoxf+eR9dQ2pOaEZbbfylfCVNBFAreUbaRq5/tgeZjl0eDePB6di0Vdxff4oyEXkcCSNbspGjOZAdO+ErVcJBeohy4SaHCn18gJ9Bo5IWq5SC5QD10kkG9Gw9FD1B1+N2q5SC5QQRcJXDllJIdfXMx7D30narlILlBBFwncPmsiJw/sieWFH6uWb8acqSU5d0JUjl8aQxeJMGFoX6qH9OfNDjyRSCRbqIcuEqG2tjanLyqS45sKukiEurq6nL3sX0RDLiIRvvGNb/DBBx9kOg2RpKigi0SYMWNGplMQSZqGXEQivPHGG2zbti3TaYgkRT10kQhf/epXOeGEE/jTn/6U6VREOkw9dJEIdXV1muUiOUsFXSSCZrlILlNBF4mgeeiSy1TQRSKohy65TCdFRSL89Kc/ZdCgQZlOQyQpKugiES655JJMpyCStISGXMzsIjPbYmbbzOzmOG3+ycw2m9kmM/vv1KYpkh4vvPAC27dvz3QaIklpt6CbWT5wL/Bp4AzgSjM7o1WbscAC4Hx3Hw/c2AW5inS5z3zmM9x7772ZTkMkKYn00M8Btrn7DnevBZYCM1u1uRa4190PAbj7vtSmKZIemuUiuSyRgj4ciHxK7p5gWaSPAR8zs5fMbLWZXRRrQ2Z2nZlVmFnF/v37k8tYpIu4uy4skpyWqmmLBcBYYBpwJXC/mRW3buTu97l7qbuXnnjiiSnatUhqNDQ04O6atig5K5GCHgIiH6o4IlgWaQ+w3N3r3P0tYCvhAi+SM+rq6gDUQ5eclUhBfxUYa2ZjzKwHMBtY3qpNGeHeOWY2mPAQzI4U5inS5QoLC3niiSf4/Oc/n+lURJLS7jx0d683s+uBciAfeNDdN5nZbUCFuy8PYp8ys81AAzDP3d/vysRFUq2goIDPfvazmU5DJGnm7hnZcWlpqVdUVGRk3yKx1NTU8Kc//YnJkydTUlKS6XREYjKzte5eGiume7mIBPbt28esWbNYtWpVplMRSYoKukigtrYWQLNcJGepoIsENMtFcp0KukhAPXTJdSroIoGmgq4euuQq3T5XJDBu3DieffZZJk6cmOlURJKigi4S6NevH9OmTct0GiJJ05CLSCAUCvHQQw9x4MCBTKcikhQVdJFAZWUlc+bM4a233sp0KiJJUUEXCWjaouQ6FXSRgKYtSq5TQRcJqIcuuU4FXSSgHrrkOk1bFAnMnDmTs846i6FDh2Y6FZGkqKCLBAYNGsSgQYMynYZI0jTkIhKoqqriN7/5TfNYukiuUUEXCTz99NN8/etfp76+PtOpiCRFBV0koFkukutU0EUCtbW1mBn5+fmZTkUkKSroIoHa2lr1ziWnqaCLBOrq6lTQJaepoIsEvvvd77Ju3bpMpyGSNM1DFwkMHjyYwYMHZzoNkaR124JeVhnih09s4lB123OKzz9lIA9de27K9ntL2UYeWrML9/Dr3oV5/MdlZzJr8vCU7UO6xpNPPsnevXu59tprM52KSFLMmypPmpWWlnpFRUWXbPvq+1/mpe0HE27fK994498vbrddWWWIbz1cRUMKDlmewVVTSrh9lh53li3mzJnD6tWr2bZtW6ZTEYnLzNa6e2msWLfroXe0mAN81OCMvvnJFsvmTG1ZbMsqQ9y4rColOQI0OixevYvFq3e1WJ7qbwySuNraWt2YS3JaQgXdzC4CfgbkAw+4+52t4l8BFgGhYNE97v5ACvNsVzKFvC2ti61Zyjbdppe2H4z64wJQmAeLvjBJQzddSLNcJNe1W9DNLB+4F/gksAd41cyWu/vmVk2Xufv1XZBju1JdzGPJ0MhUs7pGuHFZVcxvCUWFedyhcfpOKasM8dzr71B9pJrz73yGeTNO0/GUnJNID/0cYJu77wAws6XATKB1QU+bv/zlL7z77rvNr1c+tRbr0YuiMR8H4KNdG2ms+WuLdaxXH4pGnQVAzdtV+LHqFvG83v3oNXJCOL5jLV53rEU8v+8Aeg4/HYDq7a9CfcuTrfn9BtNz6MfC8TdXQ2Nji3hB8Un0OOkU3J2arS9H/ZsKBg6jx4mj8YZ6ara9EhUvHDSSwsEjaaw7xkc71raIVQPffG83FRefy/wLR7Ny5cqo9SdPnsyYMWM4dOgQzz77bFT87LPPZuTIkezbt48XX3wxKn7uuecydOhQ9u7dy+rVq6Pin/jEJxg8eDC7du0i1rmR6dOnU1xczI4dO6iqiv6j9KlPfYq+ffuydetWXnvttaj4xRdfTK9evdi8eTNvvPFGVPzSSy+loKCADRs2xBwDv+yyywBYt24db7/9dotYxa4jPP7+UD46VovlFRA6XMOCxzYCqKhLbnH3Nn+AywkPszS9/iLhIZXINl8B3gE2AI8AI+Ns6zqgAqgoKSnxZF188cUOtPgpGDjCR81f4aPmr/CeI8ZHxXv8zdjmeOGQk6PiPUvObI4XDBgaFS86dUpzPK9PcVS89xl/3xy3wp5R8b6TLvJR81d4yXeWR8UA73fOZT5q/gofeeOymPH+f3e1j5q/wof/y29jxgdMv8ZHzV/hp33z/pjx++67z93dX3nllZjxJUuWuLv7qlWrYsZXrFjh7u5lZWUx488//7y7uy9evDhmfN26de7u/qtf/SpmfOvWre7uftddd8WMv/POO+7u/v3vfz9m/OjRo+7ufuONN0bFzKz5s3PNNddExfN79gkf+39b6iNueKj5fTzvjlVJf0ZFugpQ4XHqdbuzXMzscuAid78meP1FYIpHDK+Y2SDgqLsfM7N/Bq5w9+ltbbczs1zeeust/vrXcA/8mdff467yLVh+IYWDRgBQd+gdvO6jFuvcMvMs/mXmBQCcedN/cuDIhy3/nYW9KBwQfrBB3fu78YaWd9zL69mbgv4nAVB7YBc0NrSM9+pDQb8h4fj+t+mRB7UR02Hyep1AQb/BuDt1+9+O+jfl9e5PQd+BeGMDdQd2RcXz+xST32cA3lBP3fu7o+N9B5Lfuz9eX0vdwVDz8pKBRXxU18hB+jJy6En86ydKOL1vddT6JSUlFBcXc/ToUXbs2BEVHz16NP369ePIkSPs3LkzKn7yySfTt29fDh06xO7d0fm9/mFvFj71Jkc/OEzDX9+PihcOHEHvop58ePgg9R8eio4PGonlF9Bw9BAN1Yf5zJlDuWH62Ob4hAkTyMvLY+/evRw4cCBq/TPPPBOAPXv2cPBgy+G5i372/+hx4uiodQx4687PRC0XyaS2ZrkkUtDPBRa6+4zg9QIAd78jTvt84KC7929ru6mYtlhWGWLeI+upa2ceYesZK03rLly+icM1Hb/3dVFhPndcNjGpr+NllSEWlW8hdLiGPAvPdkm3ngV5/Ojz4TH3ROfrNzHC3dpskW/GlVNGJj39s63ZSwN6F1L5g091Jj2RlOtsQS8AtgIXEp7F8ipwlbtvimgz1N3fCX7/HDDf3ae2td1UFPTz73yG0OGauPHhxUUJn9xqKrR7D9fQu0c+H9Y2tNl+eHERL93c5peQhJVVhpi7rIrG9ptKByQyBbStz1BxUSFVt6qgS3bp1Dx0d683s+uBcsLTFh90901mdhvhsZzlwL+a2aVAPXCQ8Jh6l9vbRjE36FDBnTV5eMzCP+bmJ2P2SNvad0fF23dkb1467qXtB/nkT56juraRvYdr6FGQx7H6xP9sHkni25tIJiU0D93dnwKearXsBxG/LwAWpDa19g0rLopb7IYVF3XpPlK1/bbEK/RNTl7wZEaGbHLJm/v+91xJR4o5pOc9FkmlnL7b4rwZp1GYH33FT2GeMW/GaSnbR1FhywceFBXmp2z7nXHVlJJMp9CtZcN7LNIROX3pf1PvNfKkXnFRIQsvHT+2BDUAAAeASURBVJ+y+cNN22kaXx/WgXH5rtZ0IvCh1buah4Van2jt1zOfD461fT4gWxQV5vFRXWPWnHTNhvdYpCO65c25pKVbyjZG3TMmWamY5ZLszJR0nlNI5UlvkVTq1CyXrqKCnn5tFcRcvFfMLWUbW3w7SaW7r8itYyHHDxV0OW6k4tuIpitKNutWt88tqwzxvcc3xpwnPqB3Ibdekrrxc8k9sc4rQMtzK23NPS8qzGfhpePTkKlI6uVUD72sMsS3/rCehjbm6hXmG4suP0tFXeIqqwyx4LGN1NS17BSk+oS6SFfoNj30ReVb2izmAHUNzqLyLfqfUuLK5plLIp2RUwU90aszU3kVp3RP7V20JZKLcurCokSv3NMVfiJyPMqpgj5vxmnk57X/LDhd4Scix6OcKuizJg/n/37hLPr0yI/b5vxTBuqrtIgcl3JqDB1ajn3eUraRJWt20+De6ftii4jkupyatigicrxra9piTg25iIhIfCroIiLdhAq6iEg3oYIuItJNqKCLiHQTGZvlYmb7gZ1Jrj4YOJDCdFJFeXVMtuYF2Zub8uqY7pjXKHc/MVYgYwW9M8ysIt60nUxSXh2TrXlB9uamvDrmeMtLQy4iIt2ECrqISDeRqwX9vkwnEIfy6phszQuyNzfl1THHVV45OYYuIiLRcrWHLiIiraigi4h0E1lb0M3sC2a2ycwazSzu9B4zu8jMtpjZNjO7OWL5GDNbEyxfZmY9UpTXQDNbaWZvBv8dEKPNP5hZVcTPR2Y2K4j91szeiohNSldeQbuGiH0vj1ieyeM1ycxeDt7vDWZ2RUQspccr3uclIt4z+PdvC47H6IjYgmD5FjOb0Zk8kshrrpltDo7PKjMbFRGL+Z6mKa+vmNn+iP1fExH7cvC+v2lmX05zXj+NyGmrmR2OiHXl8XrQzPaZ2Wtx4mZmPw/y3mBmH4+Idf54uXtW/gCnA6cBzwGlcdrkA9uBk4EewHrgjCD2MDA7+P3XwDdSlNddwM3B7zcDP2qn/UDgINA7eP1b4PIuOF4J5QUcjbM8Y8cL+BgwNvh9GPAOUJzq49XW5yWizb8Avw5+nw0sC34/I2jfExgTbCc/jXn9Q8Rn6BtNebX1nqYpr68A98RYdyCwI/jvgOD3AenKq1X7G4AHu/p4Bdv+BPBx4LU48YuB/wEMmAqsSeXxytoeuru/7u5b2ml2DrDN3Xe4ey2wFJhpZgZMBx4J2v0OmJWi1GYG20t0u5cD/+Pu1SnafzwdzatZpo+Xu2919zeD3/cC+4CYV8J1UszPSxv5PgJcGByfmcBSdz/m7m8B24LtpSUvd3824jO0GhiRon13Kq82zABWuvtBdz8ErAQuylBeVwJLUrTvNrn7C4Q7cPHMBP7Lw1YDxWY2lBQdr6wt6AkaDuyOeL0nWDYIOOzu9a2Wp8JJ7v5O8Pu7wEnttJ9N9Ifp34OvWz81s55pzquXmVWY2eqmYSCy6HiZ2TmEe13bIxan6njF+7zEbBMcjyOEj08i63ZlXpG+RriX1yTWe5rOvD4fvD+PmNnIDq7blXkRDE2NAZ6JWNxVxysR8XJPyfHK6CPozOzPwN/ECH3P3f+Y7nyatJVX5At3dzOLO+8z+Ms7ESiPWLyAcGHrQXgu6nzgtjTmNcrdQ2Z2MvCMmW0kXLSSluLj9Xvgy+7eGCxO+nh1R2Y2BygF/j5icdR76u7bY28h5Z4Alrj7MTP7Z8Lfbqanad+JmA084u4NEcsyeby6VEYLurv/Yyc3EQJGRrweESx7n/BXmYKgl9W0vNN5mdl7ZjbU3d8JCtC+Njb1T8Dj7l4Xse2m3uoxM/tP4NvpzMvdQ8F/d5jZc8Bk4FEyfLzMrB/wJOE/5qsjtp308Yoh3uclVps9ZlYA9Cf8eUpk3a7MCzP7R8J/JP/e3Y81LY/znqaiQLWbl7u/H/HyAcLnTJrWndZq3edSkFNCeUWYDXwzckEXHq9ExMs9Jccr14dcXgXGWniGRg/Cb95yD59leJbw+DXAl4FU9fiXB9tLZLtRY3dBUWsat54FxDwb3hV5mdmApiELMxsMnA9szvTxCt67xwmPLT7SKpbK4xXz89JGvpcDzwTHZzkw28KzYMYAY4FXOpFLh/Iys8nAb4BL3X1fxPKY72ka8xoa8fJS4PXg93LgU0F+A4BP0fKbapfmFeQ2jvAJxpcjlnXl8UrEcuBLwWyXqcCRoNOSmuPVVWd7O/sDfI7wONIx4D2gPFg+DHgqot3FwFbCf2G/F7H8ZML/w20D/gD0TFFeg4BVwJvAn4GBwfJS4IGIdqMJ/9XNa7X+M8BGwoVpMdA3XXkB5wX7Xh/892vZcLyAOUAdUBXxM6krjleszwvhIZxLg997Bf/+bcHxODli3e8F620BPp3iz3t7ef05+P+g6fgsb+89TVNedwCbgv0/C4yLWPf/BMdxG/DVdOYVvF4I3Nlqva4+XksIz9KqI1y/vgZ8Hfh6EDfg3iDvjUTM4EvF8dKl/yIi3USuD7mIiEhABV1EpJtQQRcR6SZU0EVEugkVdBGRbkIFXUSkm1BBFxHpJv4/Nl4hd0uBqSgAAAAASUVORK5CYII=\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "3oYzELNxl04e",
        "colab_type": "text"
      },
      "source": [
        "以上"
      ]
    }
  ]
}