5_2_meta_learners_issue18_issue36.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "colab": {
      "name": "5_2_meta_learners_issue18_issue36.ipynb",
      "provenance": [],
      "collapsed_sections": [],
      "toc_visible": true
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "aoxI3DOK9vm2"
      },
      "source": [
        "# 5.2 Meta-learners（T-Learner、S-Learner、X-Learner）の実装\n",
        "\n",
        "本ファイルは、5.2節の実装です。\n",
        "\n",
        "第1章から取り上げている、人事研修の効果について因果推論を実施します。"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "2XdIDbdlejUk"
      },
      "source": [
        "## プログラム実行前の設定など"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "XZFKJwcu-_Oj"
      },
      "source": [
        "# 乱数のシードを設定\n",
        "import random\n",
        "import numpy as np\n",
        "\n",
        "np.random.seed(1234)\n",
        "random.seed(1234)\n"
      ],
      "execution_count": 1,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "hx1idArc_F15"
      },
      "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": 2,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "AWqP6yeQlI_t"
      },
      "source": [
        "## データの作成"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "DpnGB2KZ_L8x",
        "outputId": "e35f7dd5-7d94-44e2-df82-a5a43422fa1a",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 283
        }
      },
      "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(-1*-5.0*x+5*e_z)  # xの効果が反対になっていたのを修正Issue：#18\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 0x7f04dd963ad0>"
            ]
          },
          "metadata": {},
          "execution_count": 3
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAASc0lEQVR4nO3df5BdZX3H8fc3mwQjVkPIaiUJJDgRTRsFvYO0zBT8SWA6SfzZZMoIlppqxc60yhRGxjIMjra049BKa5GhqNggUqXrFCejAuOMQ5ClCBiYQIgjJFBZQZxpoZAf3/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"
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "BHcdUlW9koTa"
      },
      "source": [
        "## データをまとめた表を作成し、可視化する"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "1EMwdGIIIPrK",
        "outputId": "cdb92fab-894a-42a8-d031-960bf40bd47c",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 206
        }
      },
      "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/plain": [
              "          x    Z    t         Y\n",
              "0 -0.616961  0.0  0.5  1.803183\n",
              "1  0.244218  1.0  0.7  2.668873\n",
              "2 -0.124545  0.0  0.5  2.193123\n",
              "3  0.570717  1.0  1.0  3.245229\n",
              "4  0.559952  1.0  1.0  3.139868"
            ],
            "text/html": [
              "\n",
              "  <div id=\"df-6a7b716a-1935-4790-b88a-6f0cd44ab2b8\">\n",
              "    <div class=\"colab-df-container\">\n",
              "      <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>0.0</td>\n",
              "      <td>0.5</td>\n",
              "      <td>1.803183</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>1.0</td>\n",
              "      <td>1.0</td>\n",
              "      <td>3.245229</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>4</th>\n",
              "      <td>0.559952</td>\n",
              "      <td>1.0</td>\n",
              "      <td>1.0</td>\n",
              "      <td>3.139868</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "</div>\n",
              "      <button class=\"colab-df-convert\" onclick=\"convertToInteractive('df-6a7b716a-1935-4790-b88a-6f0cd44ab2b8')\"\n",
              "              title=\"Convert this dataframe to an interactive table.\"\n",
              "              style=\"display:none;\">\n",
              "        \n",
              "  <svg xmlns=\"http:https://www.w3.org/2000/svg\" height=\"24px\"viewBox=\"0 0 24 24\"\n",
              "       width=\"24px\">\n",
              "    <path d=\"M0 0h24v24H0V0z\" fill=\"none\"/>\n",
              "    <path d=\"M18.56 5.44l.94 2.06.94-2.06 2.06-.94-2.06-.94-.94-2.06-.94 2.06-2.06.94zm-11 1L8.5 8.5l.94-2.06 2.06-.94-2.06-.94L8.5 2.5l-.94 2.06-2.06.94zm10 10l.94 2.06.94-2.06 2.06-.94-2.06-.94-.94-2.06-.94 2.06-2.06.94z\"/><path d=\"M17.41 7.96l-1.37-1.37c-.4-.4-.92-.59-1.43-.59-.52 0-1.04.2-1.43.59L10.3 9.45l-7.72 7.72c-.78.78-.78 2.05 0 2.83L4 21.41c.39.39.9.59 1.41.59.51 0 1.02-.2 1.41-.59l7.78-7.78 2.81-2.81c.8-.78.8-2.07 0-2.86zM5.41 20L4 18.59l7.72-7.72 1.47 1.35L5.41 20z\"/>\n",
              "  </svg>\n",
              "      </button>\n",
              "      \n",
              "  <style>\n",
              "    .colab-df-container {\n",
              "      display:flex;\n",
              "      flex-wrap:wrap;\n",
              "      gap: 12px;\n",
              "    }\n",
              "\n",
              "    .colab-df-convert {\n",
              "      background-color: #E8F0FE;\n",
              "      border: none;\n",
              "      border-radius: 50%;\n",
              "      cursor: pointer;\n",
              "      display: none;\n",
              "      fill: #1967D2;\n",
              "      height: 32px;\n",
              "      padding: 0 0 0 0;\n",
              "      width: 32px;\n",
              "    }\n",
              "\n",
              "    .colab-df-convert:hover {\n",
              "      background-color: #E2EBFA;\n",
              "      box-shadow: 0px 1px 2px rgba(60, 64, 67, 0.3), 0px 1px 3px 1px rgba(60, 64, 67, 0.15);\n",
              "      fill: #174EA6;\n",
              "    }\n",
              "\n",
              "    [theme=dark] .colab-df-convert {\n",
              "      background-color: #3B4455;\n",
              "      fill: #D2E3FC;\n",
              "    }\n",
              "\n",
              "    [theme=dark] .colab-df-convert:hover {\n",
              "      background-color: #434B5C;\n",
              "      box-shadow: 0px 1px 3px 1px rgba(0, 0, 0, 0.15);\n",
              "      filter: drop-shadow(0px 1px 2px rgba(0, 0, 0, 0.3));\n",
              "      fill: #FFFFFF;\n",
              "    }\n",
              "  </style>\n",
              "\n",
              "      <script>\n",
              "        const buttonEl =\n",
              "          document.querySelector('#df-6a7b716a-1935-4790-b88a-6f0cd44ab2b8 button.colab-df-convert');\n",
              "        buttonEl.style.display =\n",
              "          google.colab.kernel.accessAllowed ? 'block' : 'none';\n",
              "\n",
              "        async function convertToInteractive(key) {\n",
              "          const element = document.querySelector('#df-6a7b716a-1935-4790-b88a-6f0cd44ab2b8');\n",
              "          const dataTable =\n",
              "            await google.colab.kernel.invokeFunction('convertToInteractive',\n",
              "                                                     [key], {});\n",
              "          if (!dataTable) return;\n",
              "\n",
              "          const docLinkHtml = 'Like what you see? Visit the ' +\n",
              "            '<a target=\"_blank\" href=https://colab.research.google.com/notebooks/data_table.ipynb>data table notebook</a>'\n",
              "            + ' to learn more about interactive tables.';\n",
              "          element.innerHTML = '';\n",
              "          dataTable['output_type'] = 'display_data';\n",
              "          await google.colab.output.renderOutput(dataTable, element);\n",
              "          const docLink = document.createElement('div');\n",
              "          docLink.innerHTML = docLinkHtml;\n",
              "          element.appendChild(docLink);\n",
              "        }\n",
              "      </script>\n",
              "    </div>\n",
              "  </div>\n",
              "  "
            ]
          },
          "metadata": {},
          "execution_count": 4
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "L6Tb2Hjk9vno",
        "outputId": "3e0ec3d4-976a-4813-9284-c01aa9015fd2",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 283
        }
      },
      "source": [
        "plt.scatter(x, Y)\n"
      ],
      "execution_count": 5,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "<matplotlib.collections.PathCollection at 0x7f04dd4892d0>"
            ]
          },
          "metadata": {},
          "execution_count": 5
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO2de5RV133fv7+5c4A7UsIM1iSWrnnIrgs1wcxEU5mYtRKjJsKxLHkqZCMvObFTZ7lx0wdYmVWUqAKpyoKUOiip27pu4saOtZRBoIyRkRdWAmkaxSiBzIwJFsRSsCRfKdHEcInFXODOzO4f9+zh3HP367zua36ftVjMPc999jnnt3/nt38PEkKAYRiGaX+6mt0AhmEYJh1YoDMMw3QILNAZhmE6BBboDMMwHQILdIZhmA6hu1knvuGGG8SqVauadXqGYZi25OTJk/8ghOhXrWuaQF+1ahVOnDjRrNMzDMO0JUT0sm4dm1wYhmE6BBboDMMwHQILdIZhmA6BBTrDMEyHwAKdYRimQ2ialwvDMEw7MjZexN4jZ/FaqYybevMY2bwaw4OFZjcLAAt0hmEYZ8bGi3jgqVMoV2YBAMVSGQ88dQoAWkKoW00uRLSEiP6CiCaJ6DQRPazY5hNENEVEE/6/X8ymuQzDMM1j75Gz88JcUq7MYu+Rs01qUS0uGvoVALcJId4kIg/AnxHR14UQx0PbjQoh/m36TWQYhmkNXiuVIy1vNFYNXVR50//p+f+4KgbDMAuOm3rzkZY3GicvFyLKEdEEgDcAPCuEeF6x2RYi+hYRHSCi5ZrjfIqIThDRiampqQTNZhiGaTwjm1cj7+VqluW9HEY2r25Si2pxEuhCiFkhxACAtwG4lYh+LLTJ0wBWCSHeDeBZAF/SHOcLQoghIcRQf78ytwzDMEzLMjxYwO6716HQmwcBKPTmsfvudS0xIQoAFLWmKBE9BGBaCPFfNetzAM4LIZaajjM0NCQ4ORfDMEw0iOikEGJItc7Fy6WfiHr9v/MAfgbAmdA2NwZ+3gXghfjNZRiGYeLg4uVyI4Av+Zp3F4D9QoivEdEjAE4IIQ4B+PdEdBeAGQDnAXwiqwYzDMMkpZWDg5IQ2eSSFmxyYRimGYSDg4DqxGYr2cJNJDK5MAzDdBKtHhyUBA79Zxim4wmaWHQ2iWKpjI17jra1+YUFOsMwHY3KxKKj1XKzRIUFOsMwHY3KxGJCml+kQG+nCVQW6AzDdDRx8qzIfVo9u2IYnhRlGKaj0eVZKfTmUbDkZok6gTo2XsTGPUdx847D2LjnKMbGiwlaHh3W0BmG6WhGNq9WuinK/Csq+/r01Rnc97+/iaIhu2LYFLNpTT8Oniw2VZtnP3SGYToelfA9dmYKr5XKyHtdKFfmIqWQ7evxcLkyVzMQENRpaAu9eTy347akl3DtPAY/dBboDMMsKKJ4vajIezks7u5CqVxx3oeA1CZUObCIYRjGJ6rXS5jdd6/DxQjCHKhq7tIEk6VdnW3oDMO0JNJMUiyVkSPCrBAopKDl6uzirgwPFubbFZWwS2TasIbOMEzLIc0iUmjO+qbhpFru2HgRlKBdea8qMlWFLlzJslwda+gMw6RKGoE4JrOISctVnVse77VSGV1ETpOfPV4XpitzNcu6AOy++90ArnmtxNHUsyxXxwKdYZjU0AXinHj5/LxXiYuQt2mx0m3w4adP48J01Z7d43WhMitQmbumzY88OQkQUJmtLpt1dAL59n/+WaXJR/qfDw8W5tv/mdEJzJkOFiDrcnUs0BmGSQ1dIM7jx1+Z14xt/tlj40V0+QJUR2+Ph5EDk/OCGkCdRg1gXrhHQQYbybaZIkV3HTptFeY5IswJ0ZC0ASzQGYZJDZ1mHRarOrOJ1PBNwjzv5SAEaoR5WhAwr0E/OHYKXzn+St025cos7t8/ie2jE07mmzkhcG7PHek2VAMLdIZhUsFFsw6isj3bXAqll8v20YnY7dRBAO7bsALDgwWtMJe4XiNQHcze9Z++Pv8F0Zv3sOuutZlo6izQGYZJzNh4ESMHJiMJuhxR3SSmboKRgBotN67boIl9WwfmhewTz7+a6rGD5qBSuVK17SP9lADstsgwTGIefvp0ZBPIrBAYOTCJol90olgqa10Kw54hI5tXw8slcUCspdCbrxGuUQamOFTmRCYVkligMwyTGOlpEpXwICCAOqGu8gwZHixg7z3r0dfjRT6ny/FzlN5goSMLf3QW6AzDtBQCVY2Z/P91xZuHBwsYf+j2SIFChKppxXb8j75neYIrcCMLf3S2oTMMk5jevBcpWZWJqNkJTbZ31bZBH3KJLiAp6G6pglCNHlW5TNrIwh+dNXSGYRKz66618LpqdeUuANctihYeHyfwxjUMX3fsYJoBacvfNjqBr02+jl6LSUcAEDGSCXzM96ZJG9bQGYZJTDAUPhwNevOOw1ot1+siXL+kG6XpSuzAm3AYvozq7OvxIARwsWw+ts5V0vWLwzVzY5opdLXn4HzoDMNkyeAj3zBOmibJoGjL3dLb4+FyZRZl3yTS1+Nh5521PuCmASct0ixyYcqHzho6wzCJUOU8CWrJNg+YuKXaVHljPrN/AsFo//C5L0xXMHKg1gc8ig0+DlnnbwnCNnSGYWKjS3Mr/3d1Z5SpAKIUWVaZSlxSt1Rma33AXWzwOSIQ3NwZvRyhN+9ZvXSygDV0hmFik7T6TxCpqQc17pEDk9h16LTSDp7Ejzu4rzxeMHNjGJmP5eYdh43HTaMARxJYQ2cYJjZpBsd0Uf0EY2VWoFSuKEu4JfHjDu8rfdp782qvFrm97pwE4LGtAxjZvBp7j5x1+sLIAhboDMPEJs3gGBdziTTNAPGrBnk5wqY1/UrTzq671tYdM2gDV51TJvUCUOf+mHUN0TAs0BmGiY2T/bkr3TB6+VUwPFjA7rvXzUd99uY9a36XLgK2/vPlOHiyqBS8wWMCVZt50L4/PFjAllsK87b0HBHu27ACjw6v0+aCzyJniw62oTMMExudD3i4qLPcxpZV0YXgV0E46jPocUOozcOe93LYfbdZ8AaPp6u8dPBksWby9+DJIoZWLtOan7KsIRqG/dAZhkmFKLVEN+45qhTqvXkPP7g8o812KIUyoA5icmnPKsPE5nf9FL269uU0+d6lRq/aJ+iDnka9VfZDZxgmU3S1RAG1b/nI5tU12wNVW/QH19+IoZXL6tYB1wpDAOaycBJVzhZAL5SDLok6rVo30LxWKmPf1oG6dgft71H7KA5sQ2cYJpL/t4qo9mNpiw5avAWAgyer5w3axgu9eTy2dQATO2+vps1NaKvWCeXgct1kr84PXSb9Crc76IPeCBs7a+gMs8BJQ3OMYz8+dmZKW2v0uR23ac+d1FZd0NjwCwEhrvqCyHs5bLmlgIMni1otXPdVkEa7XWCBzjAB0rBxths2zdGlP3QTnSa3xrgCLs65gmxa06+sF7ppTf/836ZkY0Mrl8V6RpK22wWrQCeiJQD+FMBif/sDQoidoW0WA/gygFsAfB/AViHEd1NrJcM0gEbYOLMiyUCkE6CqyE1df+g0WlMOk7gCLs65ghw7M+W0XKdtm7RwE0nb7YKLDf0KgNuEEOsBDAB4PxFtCG3zSQAXhBD/BMA+AL+RWgsZpkG0gh9xHFT5vKMEtJjsxa79YbMfq1D5sIcFnMq2H+dcQZrlXpi03S5YNXRR9Wt80//p+f/Cpq8PAdjl/30AwOeIiESzfCIZJgat4EccB5tftQ2d5qjL0aLrj6iaq8msAdi/mOIKwkaYPnQkabcLTl4uRJQjogkAbwB4VgjxfGiTAoBXAUAIMQPgIoC3KI7zKSI6QUQnpqbUnz0M0yx0L3QjXvQkmEwmLlq6TnMsxOyPKB4zw4MFPLfjNpzbc0fdRGhWX0wuXwbtipNAF0LMCiEGALwNwK1E9GNxTiaE+IIQYkgIMdTf32/fgWEaSLu+6CYBGzWXiADwdxcvY9voBC5dmakLpbf1R1LzT5CsvpgaYfpoFpG8XIQQJSI6BuD9AP46sKoIYDmA7xFRN4ClqE6OMkzbYDMBtCoqk4nEZnoZGy/WpY2V/tilcgVeF6Gvx5svEbdpTT/2HjmL7aMTyv5Jav4JkqVpJGvTR7Nw8XLpB1DxhXkewM+gftLzEICPA/gmgHsAHGX7OdOOtOOLLtu7bXRCuV6n0YZt1CoqcwI9i7ox/tDtTl5AaWrVjfAK6TRcNPQbAXyJiHKommj2CyG+RkSPADghhDgE4HcB/D4RvQjgPIB7M2sxwywgXNwR5TY6dBqta3GK13xb/P37J+uiLMPad5padbt+MTUTFy+XbwEYVCx/KPD3ZQAfTrdpDNPemISxq6C2acQ2Lduk0bpqzUvzHh546pQxj4kkba26Hb+YmglHijJMBqiE8fbRCZx4+Xxd8ildwI6LPdqkZdvKobmksSUAV2dmUa7MabcJp7OV7WKtuvFw+lyGyQBd+lVCVeMtldW1K4NC+OYdh+sCPoLbveZ7kqggAOf8VLA6xsaL2D46oT2GCzKdLQvsxsHpcxmmwejMGQLQCnOgVlvXadAEdd7tIC426+HBAk68fB6PH38lllDPEbEwbzE4fS7DZEAS1zppVtHVr7QJ3yg260eH12Hf1oF5n2xXvC7CZz+ynoV5i8ECnWEyYGTzaq2A7OvxrHU4i6Uy9h45iy23FGrqW9qEeZwgmWC0pi46NMz1S7pZmLcgLNAZJgOGBwt47zuWKdfd8e4bjaH1kmKpjIMni9i0ph95L6f1MpHIUmfhGptRCle4FH0GgNK03mzENA+2oTNMRnz3+2o797EzU3h0uKpF29wOy5VZPPH8q1ZhDlzzF5ceJkvzHi5dnUFltrqvSzpguVzlcx5kad7Dxj1H2ZOlxWANnWEywiVqUlWKLYyLMAeA3h6vJo9KqVyZF+aScmUW9++fNGrsw4MFzFnO+YMrM6nka2HShQU6w2SEbmK0i6hG+KlKsUUl7+UgBJwiP2eFmBfE20YnMPjIN+qEsW1Sd3ZOHTHKNBcW6AwTIE6xZN0+Onv0rBA1Gq0pYtPmeRLMFnjR4A5p4sJ0pU7DdrWlB2n1vPELAbahM4xPnBJ0LvvYcqCYAo1MmnuOCC/t/sD8771Hzlr903WEI1Dl/5/ZP4E5x8+HVs8bvxBgDZ1hfOIUVLDtY7JHy0nMS1dnYrX3o+9ZXvNbpVXL9LeE6gBgQqlhK5reRYicJ51pDCzQGcYnTupXl31MlZD2HjlbN3FpI0eEj21YgUeH19UsVxVu2Pvh9Rh/6Hac23MHPvuR9VYzStDssuvQaagyuCzu7sLee9Z3ZIGIdodNLgzjEyf1q8s+pgyE2zU5zHXkyByhacpOKJfvOnRaaeIRAEaenJzfVmcGKlfmOAtii8IaOsP4xClB57KPqeRZVLtzeEI1KsODBUzsvB19PZ5yfWVO4P79k9bjx5k8ZrKHsy0yTACXPOVp7BPc11Y1SIWMCo2LKZMjUB2Uugi4dLW+XdctymEu5CLJWRcbB2dbZBhHopgSwoJ839aByAItmD88ioeKynYfjhIlqobo9/Z4EAK4WK7MDzi2XOjlyix68x683FyNjd/LEbxcV505Jm7dUCZd2OTCMDFIs7q9TI712NYBZ9/vcHBSuD2lcgUXpisQqPqZl8qVmnZuWtNf56kS5mK5Ujf5ufee9Vp/d/ZDbz6soTNMAJX5BKivwLPr0OnUqttLXPOoANds6XI/1/qgsp3Hzkxh7z3r8fDTp3FBk2jrpt688otF9zXBfujNhwU6w/iogoS2hbxQiqUyRp6cREUTbZNUS5XCM3xeFeXKLLaNTsQKKHqtVJ4/l0qomyaD064byqQHC3SG8XHVcnXCHHDXUpNMpIYplspOhS+CdBHhwbFTOHiyWHfNvXkPu+5aa3V/5LqhrQd7uTCMj83zw4XHtg4AMAs7nWeLFKRxQ/ijCnXd9kk9aJhsMXm58KQow/gktQFL327bZKnuS6BUrszvGwcBOFccktur4MnN9oUFOsP4xMkwKMl7Oey8c601t8vYeNHqLmjLuaJDatZRhLoKntxsX1igM4yPjOjszddHUdpE7OLu6qtkyu0iTS02ZoWIPLAQql8DG/ccdXJJDO4XhCc32xsW6AwTQIbGP7Z1oMb/ep//W0epXMH20QnkPfUrJRNxuUy6ytQArpp20BYu65B2d9kFet7L4b4NKzjJVgfBXi4ME8AW/WkK0xcApitz8LqoxhMm7+WwaU0/vnL8Fev5pYYs/b/f8cAzVp/08FrboEEAe6Z0KCzQGcbHVqzCNUz/+iXd6FnUPT8obFrTj4Mn9RGkOSLMCaEUsh99z3KngcAV9mDpbFigMwsOXc6TLiJjZSHgmmDfuOeoVqiXpisYf+j2+d8b9xzVas22pFYy5/njx1+p08TzXg5LvC5lpGdfj4fLlTkO/llgsA2dWVCYcp7oTBuqic6Rzau1E6VhLxGTG6CLzfrR4XU4t+eOOrv+7rvXYeeda5Xpe3feuVabspfpXFhDZxYUUXKeSFRufMODBZx4+Xyd5iy14OBXgErzB6pCNoqANWWC1AUysQBfWLBAZxYUUYNmTGaKR4fXYWjlMmUyr6AtXiXM0zR/cPUgRsICnWkrkuZAseUBB2onKTet6cfeI2exfXRCeT6VMNXZzE2TnwyTBizQmbbB5oXigipTYJDgJGXc8+m+AuaEwLk9dzi1k2HiwAKdSYU0sgfajmEKq3c9VzhTYNDLJXzOuOeLU2yaYdKABXqHkWZa1ijnTKo5uxzDFFYfBVebc9zzcb5wplmw22IHkWZZtCjYElKldQydhhslB3mUSvVxzydzwrDLINNorBo6ES0H8GUAP4pqlPEXhBC/FdrmfQC+CuCcv+gpIcQj6TaVsZGGSSIOaWjOLsfQab6b1vTPB/rkfBfBQujrRPUFsH10AttGJ+q2tZ3PRdNmzxOmGbiYXGYA3C+E+Csi+iEAJ4noWSHEt0Pb/T8hxAfTbyLjSlomiaikYTN2OYaqUo4Mqw+7CIZNNqrBLpjQSmUi4so8TLthFehCiNcBvO7//QMiegFAAUBYoDNNplmTcS6arM2276oNhzVfU1h98OvE5qoYrM8ZDsxhAc60C5Fs6ES0CsAggOcVq3+CiCaJ6OtEtFaz/6eI6AQRnZiamorcWMaMqkBDIybjbDZjF9u+6Rgm27dNUMuvE9eiEY2ad2CYLHCuKUpE1wP4vwB+XQjxVGjdDwOYE0K8SUQfAPBbQoh3mo7HNUWzoRFeLlHPoUtkpcv8F06edenqDCqztelod99dTVq1fXTCWEdTnmPVjsPuF2hoG8M0G1NNUSe3RSLyABwE8HhYmAOAEOIfA38/Q0T/g4huEEL8Q9xGM/FIw0RgEthxXBSj2PbDxy+V6zMJSlPKpSszRmEe/DopOESIurSZYVoZq8mFiAjA7wJ4QQjxm5pt3upvByK61T/u99NsaLsT1WWuWdjMI3FcFKO4/7kmzyqWykphLwmbfVTmKJMRhoOAmHbERUPfCODnAJwiogl/2a8CWAEAQojPA7gHwKeJaAZAGcC9wtWWswBII/CmUdhcH+N40kRx/3PVjHOaDIaA2lyi81gB6qsQcRBQaxPF5NeMQLtm4uLl8mew1MgVQnwOwOfSalS7E36Ipq/ONMU/PA42gR3HkyaK+59L8qy8lzNq8TphHCf9LBOfLIRpFOWonRSptOBI0Qi4mE1UJgtVRRmgNe20NvNIVE8a2WfbR6sfd/u2DuC5HbdpXyjV8YMQgC23FLQFlPt6vMgv6/BgASObV+Om3jxeK5Wx98jZljWJtQtZRS1HMfmlEcHcbrBAd8T1AY1SQKEV7bQ2gR0lrF3VZ9tGJzD4yDe0L3bw+CoEgGNnprTt3Hmn0mPWSLNSJnQyWQnTKCa/ZgXaSR4cO4V3PPAMVu04jHc88AweHDuV+Tk5OZcjrmH1rg9Lq9ppXcwjLp40Y+NF3L9/UmnnvjBdwciTkzXnC7dheLCAm3ccVnqxvFYqx4ri1JkAmpUyoZPJSphGMfklDbRLYjJ6cOxUTXHvWSHmf8s6sVnAAt0RnV03/IDqHqLevIfrFne3hZ02qeuj1Hh1k5YAUJkT2HXodKI0tFHaabKnNluT60SyilqOMsHuuq1KcANIZH9/4vlXtctZoDeZsfEiCFBqi+EHVPcQ7bprbWwh2eoz9S6TwCpMbodAumloTVo45y9Pn6xSCEf5MnPZVjfQL+7uSvTVplNmTEpOGrBAd2DvkbNKYU6o96hIO6FTq8/Uq9qXFlH60jbombTwfVsH2HUxZbJMbBbly8y2rW6g1ykkabjVbtxzNDOlbMELdBftV3cTBcw24DRoZfuuyU7uQl+PZ93G1V5vG/RMWjhnVbxGml+DrZTYTHddUc1qrl9tH33P8hobepAslbIFLdBdtV+dMNB5YqSJ7oErlsq4ecfhzIWP7kVwsZOb8HIUyyNFhcugZzMBtJLwaRat/jUYF9N16d7tvh4Plytzsb/apJ38iedfVb4jWSllC9pt0dW1qllZDAGzRpC1i53Jnc/mntmb92pcGz+2YUXN7733rE/tYXaZ1OyUKkJZppCI62rY6mktdNd1//5JFEvluqhJ6f6a9Hl5dHgdXtr9AW1UZprmScmC0tDD2qar50ozP8ltVeqB7EZ70wtu+lRNOgkcFddJzXbXwrPWoON4+7SDVq9rv9ScBTDv9BCuXpXGNeieT0K1/9Lsp7YX6K42P9WD5+q5AjRPGIQHE52BIwsXO9MLrntIc0QN13zjeFS0uueQiqznU+J4+0RpU7P63CWdhED12ZWRwkAyYR5OAa07Z9qKWNsK9LHxInYdOl3j+mbSDnQlyMJCPakpJXwjiYDSdCXRAywHk7Hxojb/d1IXO9XLZnrBdUK0GWaMqF9Q7aBVqsjaXz7OwOjapkb3efg99HJUk1Nfha58YZxz21JAS9JWxNpSoIc7LEhYO5A3VjdCy8+sLFwMXQcbV6K4TwbbJK8/SgHlB546hS23FGrqdQLXXvBW8wyJ8gXVyp5DJrL2l49zT13bpOtzVdm/pKjeQ6+L0NfjoTRdQZfBpTDYtvv366OZTew6dLpp6T/aUqDbJuTkqGcS/JI0K9PY2pVUaER1nwxff5QCyuXKLI6dmcLuu9dpX/B2tUm3a2RoVsE6QaLeU9c2mfo2bW1d9TxX5gQuTFdQUBQW1zErROR2jY0XrQFzkiwcK9pSoNtePDnq2QRs2h3qIhCSCI2o7pOm65caiKmEm8yZ0o5C20QjI0PT9usG7JGPaX412Y7nqtXb7NhBZSfpNdgGj4Mni9hySwHHzkzhtVLZqLFHnQ8weQQ1Iv1HWwp008MRFNKmByhsdsi6XcFt4hJVQ7MNHrbPzlYXcHFphKYLZGM3Ng2waZ/P9Xgug76Lt1axVE7lGlwGj2Nnpua/zG1f8i7zAdtHJ3Di5fPGd44ImT/3bSnQdQ9HX4+HnXeunR/ldV4sMkJxe8r2O5eH9tKVmUiuSmEBGNQsbA+GywCjIwsB9+DYKTx+/JX5e9KsychG2f9ttvrwvd20pt/53sY5nw3VfEuYuGbDYJ+bnsltoxN1y6Ke0+U9DMcoANBGPbvMBwgAjx9/Bb09nrb+gVzOkaIhXF5I3QQiAFycrmTSueF2Lc17qMzO4dLVaze/VK44n0+lCRw8WXT2JHF5sMMQEEmYRHEbDQpzSbMmIxthSjLZ6lX3NhgqHnwuAbfBJ+7cwNh4EQ8/fbpGEJm+3qKaDcPPyMc2rFA+CyainnOJV59cK4gqRgFwK0domssSor6ilkqxzOq5b0uBDthfSNMDMBf6HbQnp52/YuOeo7h0tbYtrjczqcYV1oqk1qXTvqJOEEf5PDYNsK0+GRmXpXlPOUG2NO85FUIpV2bx8NOna0LQTX0cZ27AxXEgyvFsx5dKSdSEEVFymNuuR/f1mcZ8wMVyBfu2DsQKYEyDthXoNqKaG5L4oJq01CQeFWl4Y6gGPtVDH8fEEmXAMbU5ia2+FWzyOkgT8y1tqS6oPt91fRxlbsDmzqvDdjzX4iFR8LrI+dm0DZS2uTPX+QBTPIhKqWvUJHzH5nIZ2bzaXNnagEv+ComtfJmtRqeJJPuaGB5MJ69JlAFH12aTD72NVi8dZ7KlJr2Hqj52va/BfouCy/HC9yENLfT6Jd3Oz6bpfPILNA3T6n0bVihzwOgKbTQqF1RbauguWtnwYAEnXj5fZ6vzcgSIql+qCdcH0aalJvGoyNIbIw0bcpRPfNW1EID7NqxIFGjVygFCOtNWjshpfiPv5bC4u0tpttENCC731cXcE6Y372nNcXGKh4TROTAA1UhrFVGim4F0TRyPDq/D0MplqRXaSIu2E+hR7La6Tgeuda7OB9VVg7JpqUluZqtFY4aJMuBkcS2tHiBkqlqj6g+VlwtgnqiLY3KK0z8XY4Sv64qHhCEA+7YOOHuZAHrXwfe+Y5k251HaJo4oSlGj4jnaTqBHDSHWdaRcltSe7KKlJrmZrRzYE1VIp30tur7vIkotV3wSG33BEggW7r9jZ6a0xzflpI/qsx3HnVWXYMp0vHDxEN05g9u5vos618E/f+k83vuOZfjzl86nmqOpXSCRcY07HUNDQ+LEiROR99NVgpd4OcJ1i7pxseyeECvJS6sbENox13a74erREPdeJL23tv2THl832abzVgpOhJpMHCr6ejyMP3S7cp3rdbhs5/oumuSAnPhs1S/bpBDRSSHEkGpd22noNu2iMivmbY46jUX10MTN59LqZpFOJtz3KvNZEpt6HBt9lECwpHMAUUxOYWEaVY0L27HjBLy5vCvBrzh5DulOHDRJmcL1OzVlhQttp6HH8ZsNaizhaEWANepOQae1EYBze+7I/HhRNe6k7Y2ioeu2dSWo9ao0/LTfoTjvuSRHhDkhOla5Mmnobee2GHTNciWYfdEUrcjYaeVyY2m7eUY9nknjTuP4YaK4wyWZKM57OWxa01/j6hjnHYry7MTxxJHMCtGSbqyNoO0EOlAV6s/tuA2PbR2oe6BVBLMvLrRoxTRpdb/vtP19bccLC6io7nJJ2xslnsAUBxA+f7j+68g7z3UAABhVSURBVO671+HYmSmrgDW9Q1GfnSjvo2xrThHJtdCUtbazoQdxmUEPBq5kFa24UGh1v++05zNMx0ta0jCt9rrainVxALL0mkwJITMRhtuxXZE0K4zuOsfGi0qXRNOz4+qJEzQv3bzjsHKbhaSstbVAB2rLs9kCV3QPSZJoxYVEq/h9mzwh0p4M0x1P5zYXtaRhlPYm8cYKKz/BdkpBa0p/YROwppQADzx1yjiBqcI18Cp4zkbmuW9V2tLkokL1+blv6wAeHV43v43qEzdptOJCIqtUBFFI0+yTZD7AlHFP0pv3UpsoTHLd8jqllt2b96xeLmFThSmVhqkwuM0Wbvp6Cb/PKlNQ8JyNDLFvVdrOyyUprZzMqdVpBZ/7qL7XOnTX4ppv3sVrJNdF+KHF7jERpmcz7nUn8RYJe9vE8RAz+Ytn8ewshPe7o/zQk7JQ/VPToBE+96oXMnhO10lt04ttsumGC3BsG53Aw0+fni+cIve/dGXGei2zc/aYiGB7TVGfcc1dSbxFwtpzlPwlwWOoBiKTVp+E8Pstv046WcAHWXACfaFiE3CNDN/XnU8l1FQVbFQEhc/YeBEjT07OJ2ArlsoYeXJyfr3JpqtaemH6WlESuX8cIVmuzGLXodPOpongpGFc+3Dc+Q1TTp4o91+X76cRX3VZlABsddpaoHfy51Wa12Z6sAGk+tDb2m1qS1xtMix8dh06XZdNszInsOvQaVy3uDu2MJY25bgaL1CtWLVqx+G6vNw2DXzTmv6aikaSTWv6jefTDQR9PR56Fl0rWJy0/J1Edf93372uKe9pq3tlZUHbCvROHn3TvjZbwEtaD71Lu01tiapN6srlqdLNyuWmrIE2kkRaqo71wFOncOLl8zh2ZkprSpIa+LEzU8r1uuWA3jSU93I1JqS00N3/3Xevi51aIwmt4pXVSKxeLkS0nIiOEdG3ieg0Ef0HxTZERL9NRC8S0beI6Mezae41okbltRNpX5vpwdatK5bKkT1AXNptaksUb5lCbx7n9twRuWBBb48+a2CjkTZ73UAR/PIw3SfVvZHCNTy49fWk53kTptXeyVbwymo0Lm6LMwDuF0K8C8AGAL9MRO8KbfOzAN7p//sUgP+ZaisVdPLom/a16YTYTb15YwRhVBc5nWAKLje9ZCq3Mx2Xrsxo29Onud6+Hg9ZOnX15r15tzrXalmmjIFBwWsSQqp7ozNf9Sxyr/4TlVZ7JxeiG6NVoAshXhdC/JX/9w8AvAAg/ER8CMCXRZXjAHqJ6MbUWxugk0ffNK9tbLyINy/Xf3Z7uWrVHJ0NNk6uDlXodXi56SWTvse9htzbklK5oh1kdt65tlqZKoCXI+y8c62zyUUKZ1fyXg677lqL53bchn1bB9wlugIC6r48TIOd6t64DK5p02rvZJTUCJ1CJBs6Ea0CMAjg+dCqAoBXA7+/5y97PbT/p1DV4LFixYpoLQ2RZXk2V4L5pWX4tK0IrQtpXtveI2eV5fa8LopcJNimaZkq9Ehsro/ByN/gNpeuzNSZD8K2/uA+S/MeiKppX8NFi23XLIXz8GDBOUthOMd5ki8BlQCU16jz+gnfG1P5u6xohXcyzEJzU3aOFCWi6wEcBLBNCPGPcU4mhPiCEGJICDHU32+enbfR7NE3XGg3HDadJGHV8GABW24pzL98OSJsuSXeg6kTwtOVucjamk3T0mm0weWu3jsyAZu0k+s062AmzWAkZalcweXKHPZtHajRdjet6a9Tnr0uQl+Pp3yOVNurrs806atCJ1hNaSiGBwvaPg7fG5fBNW2a/U4yjho6EXmoCvPHhRBPKTYpAlge+P02f1mmNHP0Nb24SV2jxsaLOHiyWJNj4+DJIoZWLot8zDjlxgC3nCRh4bxpTT8Oniwa61/q/MNt12Xzw3746dNWbx3Zr8HrIgBbb11ekyIieH3h7cOE+8XFXlzozRtTB5jiA1y1YFv5u6xYaBpxq+Hi5UIAfhfAC0KI39RsdgjAz/veLhsAXBRCvK7Ztmmkmcvb9uImmQhK01sgykSjJO/lcJ8mb4bsw1U7DmP76ETNxOnBk0VsuaWg3W/76ITWPzzOdcg83QMPfwMXps0aPKBPqKVz/dMN2kHdeolX+wq52ItNHj1S4OpytwBw0oIX4oQg46ahbwTwcwBOEZE04P0qgBUAIIT4PIBnAHwAwIsApgH8QvpNdSNKFGIS326b5ptkIsjkomYrfuwS2KGyRwdZ3N2FoZXL6rRWWxmzcmUWh7/1el3tyXm7suZ8prYEr6lcma2Zq1B9EYQJ3oeoXhguCbiCUaTDgwWnLIHyvpg0bdOg7uKq2Yg0DUzrYRXoQog/g2XOXlQzfP1yWo2KS9QoxCSmEdOLa0ol6vKCmQaLsLbmEoUZDuwYGy8aQ+qlB4kMegkOBDb78IXpCsbGizXtcrUr6/K4BK9pVoj5/nU5bpL0qq7mquBzZEpTC9R69MjtVM9DGi6AbP5YeHRUtkVTRjpdYqe49SYBtZdLr8a7IkqmQtcMeTJ8O3h+FcGMfME2Z0U4A6Ap4x5QvY6dd65V9s/i7i6lBt/X46E0XbGmgS3EvAdAtEyFpjqjcbTktLJKMp3Hgsm2aDJV9OY9pWBIw0eWALx16RKseksef/7S+ZpsfbYvhPv3T85XNZcve1h70wmtC9OVeduxyXsh7AmSJBeJC+H7sFTT95LLlVntpKaurRemK9p7GkSXfsA1EZnc3jZomvJ6p/UFyDZwxkZHCXTTJ/KlqzPwuqhmUi7JC6Iyb6jObctTYqoSIzEJEheCNVWzFuZAfWSqzfW5XJlDuTIX+TxXZ2a1Zd9qj19rEokiYOW2poHQxQMoTmk5gG3gTDQ6SqCb7NqVWVGXYS7JCxJFOMrz2cwcQU+WsN04Li75QNLmzcszNXb0ksYDxYW+Hk/rwTIdYRBI2+tIogokS2sCnm3gTFTatgSdygVRBjbo0AmGOEQREALV3CPhcHTdcXUCJGqUX7gEWqNCsCtzYn5gGhsvoitmdKLMCuiSCsBGFl5HqhB9oPWSVDELh7bU0FUa0PbRCWwbnUChN6+1rcqEU3IfnaeILZ/33iNnrZ/5YUrlynxEYmm6gi6NGWVp3tMKkDnfXc91QvPKTK0GO7J5NUYOTKIyW3veLt9uEd3ooadYKmPVjsNOJpEgsn/CfZ/E9p/U9hzVO0Z3fzohcRzT2rSlhq4LEAGqL5O0lwdRCZaw1mQrxBsO949KZU6gZ1E3zu25A5/9yPq6NgLAD67MGLMjRhEK4esbHizgukX1Y/icQZiHW5j3cnjnj1zn3IaoA5/sn6Dmqwop12VUDJNG+HmUIJ2x8aLWx9flKyHN4Ddm4dGWAt0m1CqzAtcv6a4RAC61KG2fyjZbqqxKbkIOBsODBSzqru/+2TmBC9MVpSAd2bw6sukg3FdRCjzkiCBwzdQjheOzn3lfpDaEMRlgdPc2nNtl551rrRGw0r6998jZRAIySo4S3debKUeLxKZQMIyNtjS5uEwwlqYrNRGLOr/eLqL56Evbp7LNlmo7l9xW2vsvXdWbEIJCITzxFsVf+6befI0ZSWfqURHMJRMOiIli+lG1CVCbJlwHLJcAnk1r+lOLDnadoHTJ0aJjIZZMY9KlLTV0l/wkYcGg22dWiHltyHaspZrJOQHUaH+mLwgBRJocI0JdelmVtrjrrnqNNSjUpNanEuYuU5Zh842tlqUOOTCkkWtEau3f3XMH9m0dqOuTY2emGj45acvRYqLVCkQw7Udbauiu4dW6faJoqvJYY+NFXLpaXyhCEtT+bF8QsrybC0KgTqs0aYvhCV2Tx8ycENosiSqCgsVUyzKMvD8qF7+0/KxVfbLdMXd4miQJCIo6+cowYdpSoAO1L3CUHNvAtYHABZlNb++Rs3XeIWGk9mdL0BT0tnHB9bM7ilCbE6ImVH1o5TI8/PRpo2unS6Ir4Fqpt4vlem8VOenXiGCZZgjIJAFBHB3KJKVtBXqQ8EskP6nj5kgJIrPpRQkikudVCciobnzB48bBVajJaj46gR4WLLrjmnKNpJ3xUh5TJzybJSDjBgRxdCiTlI4Q6K6CIm7oezBtqw0pKOVLHRY4SScSoxJFqJkGjbBXRxxhmfakn+2+t6OA5OhQJgkdIdBdBUUS2+msEFbt2kX7M4WyJzmujihCzaR1h7ePIyzTnvRzue8sIJmFREcIdFdBkURDlgUVHj/+Sl0JM+H/X67MYtvoBB5++jTuePeNeOrk92ryjRRLZXhdBC9HVnu8JEc0n85AZ3u2zSG4CrWoWndUYZm2TZu9Qhimlo4Q6CZBMTZetE722Qj6YA+tXFZXR3P0L16tyeJ4YbqCrxx/RXmsypxAj9eFH/mhJU6Dy5xv5lGZFk68fB5fm3y9xv88qa81kJ2JIo6ZxjRYsVcIw9TSEQUudIULttxSwOhfvqrVhguGcmxBtz6VUEtaKOKxrQNO+0v/ZdV2JhNQqxZCiJJW1laQImrBCobpBDq+wIVKs9y0ph9PPP+qtYpPHKHw4NipOtNLVEw50oOMbF6tdT00nb9VzQ5RzDQ2G3k7TnoyTJZ0hEAH6v3SH3jqlFMVH91gsPfI2bpKQvLYSYW5PL+LTV+6E0b9EsjC7JC0aENUXGzkPOnJMNdoK4EeFiib1vTXFDGWttj7909aXQyDAk81GOhc4eKkztWd36VCPKC2PZvMLVn4WmfhQ26DbeQME422yeWiykT3leOv1Py+/8lJjDxpF+Zd0Ge+s2VcTMOUEZxkNRXkAKqeLdtHJ7C4uwt9Pd58rpL7NqxQ5qa5blE2NuRmFG1II98Lwywk2kaguwQFzc6JGm8THUt7PK3As33m67RDQnWiU5eEKUcEQrWK0BKvC9tHJ+bzuZhye8sBq1Su4HJlDvu2DuC5Hbfh0eF12HJLoS6xlsPlx6IZLoJR0tYyDNNGJpc0BceF6YrWp9v2ma8zf9y3YYUxva3UxFVmC5s3jiQcNHPszJS2aEfaQq9Z5g+2kTOMO22joacpOGRyLFURAdtnvkpr3Ld1AI8Or9Oul1qlzmxx7MwU9t6zvmYfHcGBrZFaM5s/GKb1aRsN3XUC0YapFJ2rK5xNa9StNwng8D66IhnBga2RWjO7CDJM69M2Aj0sUHp7PFycrjgXNiZUC1SogoiAaznKwx4z0n1RpsVNIsCiCGCXqMpGZxNk8wfDtDZtY3IBrlWo2bd1AD2LujEHt2o70ixyZcYs/oNmmJEnJzFyYDLV+o5RzBYuE4I8acgwTJC2C/2PmtNcTkgmCdMPkjSkvtHBOa12foZhktFRof9RcpoHS57pwuejknTCsZlmi2YEBzEM0zjaTqC7CtSwJp0kdW4Q3YRjO2i+XFWeYTqbtrKhA24eHCq79Mjm1fByaou7zFFuW6azd6uiWEeenMTgI9/AzTsOY+Oeo4ls72nB+cMZprNpO4Gumlj0uqgmLF47MaiYLujNe9j74fV1fuCqZbrjqjTfypzAhelKahOqaaAbDDk3CsN0Bm1nconrD733yFllWoDrFnfXlCvTnc+Ei4bbCqYNrirPMJ1N2wl0IN7EYpbmBlf7fLNNGxwcxDCdTVsK9DhkGVXpGsXaCqYNDg5imM6l7WzocckyF0k4wKc37zlPqDIMw6SFVUMnoi8C+CCAN4QQP6ZY/z4AXwVwzl/0lBDikTQbmQZZmxvCmm87uDEyDNNZWCNFiegnAbwJ4MsGgf4rQogPRjlxmkWiGYZhFgqmSFGryUUI8acAzqfeKoZhGCZV0rKh/wQRTRLR14lorW4jIvoUEZ0gohNTU1MpnZphGIYB0hHofwVgpRBiPYD/BmBMt6EQ4gtCiCEhxFB/f38Kp2YYhmEkiQW6EOIfhRBv+n8/A8AjohsSt4xhGIaJRGKBTkRvJSLy/77VP+b3kx6XYRiGiYaL2+ITAN4H4AYi+h6AnQA8ABBCfB7APQA+TUQzAMoA7hXNSrLOMAyzgLEKdCHERy3rPwfgc6m1iGEYhonFgokUZRiG6XRYoDMMw3QILNAZhmE6hAWTbTEJnJeFYZh2gAW6BS6szDBMu8AmFwumwsoMwzCtBAt0C1xYmWGYdoEFugUurMwwTLvAAt1ClpWOGIZh0oQnRS1wYWWGYdoFFugOcGFlhmHaATa5MAzDdAgs0BmGYToEFugMwzAdAgt0hmGYDoEFOsMwTIdAzSouRERTAF6OufsNAP4hxeakBbcrGtyuaLRqu4DWbVsntmulEKJftaJpAj0JRHRCCDHU7HaE4XZFg9sVjVZtF9C6bVto7WKTC8MwTIfAAp1hGKZDaFeB/oVmN0ADtysa3K5otGq7gNZt24JqV1va0BmGYZh62lVDZxiGYUKwQGcYhukQWlagE9GHieg0Ec0Rkda9h4jeT0RniehFItoRWH4zET3vLx8lokUptWsZET1LRN/x/+9TbLOJiCYC/y4T0bC/7veI6Fxg3UCj2uVvNxs496HA8mb21wARfdO/398ioq2Bdan2l+55Caxf7F//i35/rAqse8BffpaINidpR4x2fYaIvu33zx8T0crAOuU9bVC7PkFEU4Hz/2Jg3cf9+/4dIvp4g9u1L9CmvyGiUmBdlv31RSJ6g4j+WrOeiOi3/XZ/i4h+PLAueX8JIVryH4B/BmA1gD8BMKTZJgfgJQBvB7AIwCSAd/nr9gO41//78wA+nVK7/guAHf7fOwD8hmX7ZQDOA+jxf/8egHsy6C+ndgF4U7O8af0F4J8CeKf/900AXgfQm3Z/mZ6XwDb/BsDn/b/vBTDq//0uf/vFAG72j5NrYLs2BZ6hT8t2me5pg9r1CQCfU+y7DMDf+v/3+X/3Napdoe3/HYAvZt1f/rF/EsCPA/hrzfoPAPg6AAKwAcDzafZXy2roQogXhBC2Ssy3AnhRCPG3QoirAP4AwIeIiADcBuCAv92XAAyn1LQP+cdzPe49AL4uhJhO6fw6orZrnmb3lxDib4QQ3/H/fg3AGwCUkXAJUT4vhvYeAPAv/P75EIA/EEJcEUKcA/Cif7yGtEsIcSzwDB0H8LaUzp2oXQY2A3hWCHFeCHEBwLMA3t+kdn0UwBMpnduIEOJPUVXgdHwIwJdFleMAeonoRqTUXy0r0B0pAHg18Pt7/rK3ACgJIWZCy9PgR4UQr/t//x2AH7Vsfy/qH6Zf9z+39hHR4ga3awkRnSCi49IMhBbqLyK6FVWt66XA4rT6S/e8KLfx++Miqv3jsm+W7QrySVS1PInqnjayXVv8+3OAiJZH3DfLdsE3Td0M4GhgcVb95YKu7an0V1MrFhHRHwF4q2LVrwkhvtro9khM7Qr+EEIIItL6ffoj7zoARwKLH0BVsC1C1Rf1PwJ4pIHtWimEKBLR2wEcJaJTqAqt2KTcX78P4ONCiDl/cez+6kSI6GMAhgD8VGBx3T0VQrykPkLqPA3gCSHEFSL616h+3dzWoHO7cC+AA0KI2cCyZvZXpjRVoAshfjrhIYoAlgd+v81f9n1UP2W6fS1LLk/cLiL6eyK6UQjxui+A3jAc6iMA/lAIUQkcW2qrV4jo/wD4lUa2SwhR9P//WyL6EwCDAA6iyf1FRD8M4DCqg/nxwLFj95cC3fOi2uZ7RNQNYCmqz5PLvlm2C0T006gOkj8lhLgil2vuaRoCytouIcT3Az9/B9U5E7nv+0L7/kkKbXJqV4B7AfxycEGG/eWCru2p9Fe7m1z+EsA7qeqhsQjVm3dIVGcZjqFqvwaAjwNIS+M/5B/P5bh1tjtfqEm79TAA5Wx4Fu0ioj5psiCiGwBsBPDtZveXf+/+EFXb4oHQujT7S/m8GNp7D4Cjfv8cAnAvVb1gbgbwTgB/kaAtkdpFRIMA/heAu4QQbwSWK+9pA9t1Y+DnXQBe8P8+AuB2v319AG5H7Zdqpu3y27YG1QnGbwaWZdlfLhwC8PO+t8sGABd9pSWd/spqtjfpPwD/ElU70hUAfw/giL/8JgDPBLb7AIC/QXWE/bXA8rej+sK9COBJAItTatdbAPwxgO8A+CMAy/zlQwB+J7DdKlRH3a7Q/kcBnEJVMH0FwPWNaheA9/rnnvT//2Qr9BeAjwGoAJgI/BvIor9UzwuqJpy7/L+X+Nf/ot8fbw/s+2v+fmcB/GzKz7utXX/kvweyfw7Z7mmD2rUbwGn//McArAns+6/8fnwRwC80sl3+710A9oT2y7q/nkDVS6uCqvz6JIBfAvBL/noC8N/9dp9CwIMvjf7i0H+GYZgOod1NLgzDMIwPC3SGYZgOgQU6wzBMh8ACnWEYpkNggc4wDNMhsEBnGIbpEFigMwzDdAj/H3oJBsohB0bfAAAAAElFTkSuQmCC\n"
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "AeC7Uv29KsXC"
      },
      "source": [
        "## T-Learner"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "EZZFqyFB9vn1"
      },
      "source": [
        "# 集団を2つに分ける\n",
        "df_0 = df[df.Z == 0.0]  # 介入を受けていない集団\n",
        "df_1 = df[df.Z == 1.0]  # 介入を受けた集団\n"
      ],
      "execution_count": 6,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "xp2P-IDT9vql",
        "outputId": "188bf4e2-dc89-4e83-ba99-b51f625f3fb8",
        "colab": {
          "base_uri": "https://localhost:8080/"
        }
      },
      "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",
          "name": "stderr",
          "text": [
            "/usr/local/lib/python3.7/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.7/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"
          ]
        },
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "RandomForestRegressor(max_depth=3)"
            ]
          },
          "metadata": {},
          "execution_count": 7
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Rx2_b0CD9vqw",
        "outputId": "788c7814-6ee3-40b2-9763-1c9529401781",
        "colab": {
          "base_uri": "https://localhost:8080/"
        }
      },
      "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",
          "name": "stdout",
          "text": [
            "ATE： 0.6678485877933638\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "XvOWVBt99vq7",
        "outputId": "ab3769f1-30c3-4f56-88b0-600bab6f092a",
        "colab": {
          "base_uri": "https://localhost:8080/"
        }
      },
      "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",
          "name": "stdout",
          "text": [
            "ATT： 0.7338369884256872\n",
            "ATU： 0.6090862934128775\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "8WXoUe6E9vrE",
        "outputId": "a7095859-00ee-4b6c-e039-75ad6e02dbb2",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 283
        }
      },
      "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 0x7f04dbb3ea10>]"
            ]
          },
          "metadata": {},
          "execution_count": 10
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXxUZZov8N+TohISFsOqEogsTeOgQaJRUHpxD+oIcQWUuaLdYjujdxQnMzAybgMtd9Jj2/Z424WLtg2NQNQitGhc0O5Po9DECZBGJ8hmoEBZgyyBVCrP/aNOxZPaUzlVp6ry+34++VB13lPnPDkVnnrrPe8iqgoiIkp/WXYHQERE1mBCJyLKEEzoREQZggmdiChDMKETEWWIbnaduH///jp06FC7Tk9ElJY+++yzg6o6IFSZbQl96NChqKmpsev0RERpSUS+ClfGJhciogzBhE5ElCGY0ImIMgQTOhFRhmBCJyLKEFETuogsEpH9IvLXMOUiIs+JyDYR2SwiF1ofJhERRRNLDf1VABMjlF8HYKTxMxPAbzofFhERdVTUfuiq+icRGRphl8kAXlPfPLzrRCRfRM5W1X0WxUhEKUpV8cILL2Dfvu/+u1966aW47rrrcPr0acyfPz/oNZdffjmuvPJKHDt2DBUVFUHlpaWlmDBhAg4dOoRf/epXQeU33ngjLr74YuzduxcvvPBCUPktt9yCCy64ALt27cKiRYuCyu+44w6ce+65qK+vx5IlS4LKZ8yYgeHDh6Ourg4rVqwIKr/vvvtQUFCAmpoaVFVVBZU/+OCDGDBgAD755BO8++677cruueceJHRApapG/QEwFMBfw5T9AcAPTM8/BFASZt+ZAGoA1BQWFioRpTe3260AFICKiIqIzpo1S1VVjx071rbN/PP444+rquq+fftClldUVKiq6tatW0OW/+Y3v1FV1c8++yxk+eLFi1VV9eOPPw5ZvnLlSlVVXbVqVcjyNWvWqKrq0qVLQ5Zv2LBBVVVfeumlkOVffPGFqqr+53/+Z1DZH//4x05fcwA1GiZXi8awwIVRQ/+Dqp4fouwPABao6p+N5x8C+BdVjTgMtKSkRDlSlCi97dmzB5dddhl+8Ytf4Pbbb7c7nC5BRD5T1ZJQZVYM/XcDGGJ6PtjYRkQZbvDgwWhoaLA7DLhq3aioroe7saltm0ME08YNwbyyIhsjSy4rEnoVgAdE5HUA4wAcVbafE1ESzHXVYcm6BoRqZ/CqYvG6BqzfcQjvz7o82aHZIpZui0sBfApglIjsEZGfiMjPRORnxi6rAewAsA3AywD+PmHRElFK2bVrFyZOnIg///nPST/3XFcdFodJ5mZf7j+BO1/+NCkx2S2WXi7TopQrgH+wLCIiShtHjx5FdXU17rvvvqSfe/G62Jt61m4/DFetG2XFBQmMyH4cKUpEcWtpaQEAdOuW3Jm457rqOvyaiur6BESSWpjQiShuXq8XQPIT+tL1uzv8GvMN00zFhE5EcfPX0B0OR1LP642hu3Ugga83TCazbcUiIkp/OTk5GD16NHr37g3AlzCfXLUFR056AAD5uU48Mem8TrVd+7sk7m1swqD8XFxxbsjV16JS+JpdMrkdPaaBRYnAgUVEmcFV68YTVVvQ2OQJu8+zU8Z2OJEGfjjEIkuAO8YVRrxhumvBDR2KI9UkemAREXVRrlo3Zi3biNYo+z20bCMeWrYxaPv08YXtBv5E6lcei2duHxv15mcm93ZhDZ2I4jbyp89i97svo+/EB5Ddv9DWWATAzgU3YNjstyN+IBTk52Lt7CuTFZblItXQeVOUiOJ2tPEITrs/h3pO2x0KBuXntvs3nL0Z3NuFCZ2I4tfq67YoWcnt5RLImSUoLx0FACgvHQVnloTdN1rCT2dM6EQUtx5OI3Fm2ZdK8nOdqLjtgrZ28bLiAlTcdgFyncEx5TodbYk/E/GmKBHF7ebis/HMckAkuTX0gvxclJeOCntzs6y4AGXFBUFdHiO9JhMwoRNR3K4tHoFV549F7/5n4GCSzhnYMyYSf2LvKpjQiShupaWl2FpaGrJsrqsOS9Y3IFRHupEDe2Db/hMhe6OYa99zXXVYun43vKpdcn7zjmK3RSLKWJnY5MKBRUSUEMuXL8f8+fPxwQcfYMCA+IbkJ4qr1o05b9ahyePrieNubMKcN32zNMaS1MfNfx/fHGuOut+EEX2x5N5LOxesRZjQiShuBw4cwObNm2HXN/1IKqrr25K5X5PH2zafSzxTC4SydvthDJ39dsR9OtLu3xlM6EQUN7umz41FuAFE7sammGvfVlm8rqHd/DKJqtWzHzoRxc2uBS5iEWkAUTKTeShrtx9OyLJ4TOhEFLdUTuhD+6X2iNC12w9bfkwmdCKKW2FhIa666qqUTOjrdhyxO4SkY0InorhNnToVH3zwAbKzs+0OJUg8qxqlOyZ0IspIDgk/QVdH5ec68eyUsdi14AY8O2VsyHliOmrCiL4WRNZe6n1PIqK0UVFRgYULF6K+PvKiEnaYNm5IxJWLQollNaNYpxOItFhHonq5MKETUdwOHjyIhoaOJc1kmVdW1KGEbl19/rvzJ3uaAja5EFHcWlpa4HDYOxd6JB1pdrlzvL0rLlmBNXQiipvX603JHi5+4ZpdRg7sgR0HTiZ10q/Akan5uU48Mek8S+eWiemdEJGJAH4FwAFgoaouCCg/B8AiAAMAHAYwXVX3WBYlEaWklpaWlE7o/iRt94yNrlo3yis3weP9rkW9scmD8hWbAMQ2t0wsor4T4pu5/nkA1wDYA2CDiFSp6uem3X4B4DVV/a2IXAngaQB/Z0mERJSyxowZg+PHj9sdRkR2tGUHqqiub5fM/Tyt2ja3jBViaUO/BMA2Vd2hqs0AXgcwOWCf0QDWGI8/ClFORBlo5syZePXVV+0OI+VFWpjaykWrY0noBQB2m57vMbaZbQJws/H4JgC9RKRf4IFEZKaI1IhIzYEDB+KJl4go7USaV8bKRaut6uXyTwB+LCK1AH4MwA3AG7iTqr6kqiWqWpJqcycTUcfdf//9uPjii+0OI+WVl46C0xHc48aZJZYuWh3L3Qw3gCGm54ONbW1UdS+MGrqI9ARwi6o2WhUkEaWmxsZGHD161O4wUp6/jTwVerlsADBSRIbBl8inArjDvIOI9AdwWFVbAcyBr8cLEWW4VO+2mEqSsWB11CYXVW0B8ACAagBfAFiuqltE5CkRmWTsdjmAehHZCuBMAPMTFC8RpZBU77bY1cT0TqjqagCrA7Y9ZnpcCaDS2tCIKNWl+kjRroYfrUQUtyuuuAJHjnS9ecdTFRM6EcXt4YcftjsEMuHkXEREGYIJnYjidvXVV+Paa6+1OwwyMKETUdxOnTqF1tZWu8MgAxM6EcWN3RZTCxM6EcWN3RZTCxM6EcWNNfTUwneCiOI2ZcoU9OsXNLEq2YQJnYjiNmfOHLtDIBM2uRBR3Jqbm+H1Bs2UTTZhQieiuI0aNQp333233WGQgQmdiOLG6XNTCxM6EcWNvVxSCxM6EcWNCT21MKETUdw4sCi18KOViOL2j//4jyguLrY7DDIwoRNR3B5//HG7QyATNrkQUdwOHjyIpqYmu8MgAxM6EcXtzDPPxPz5XBM+VTChE1FcVBWtra3s5ZJCmNCJKC7+If9M6KmDCZ2I4tLS0gIA7LaYQpjQiSgurKGnHiZ0IopLt27dMH/+fPzwhz+0OxQy8KOViOKSk5ODf/3Xf7U7DDKJqYYuIhNFpF5EtonI7BDlhSLykYjUishmEbne+lCJKJW0tLRg586dOHbsmN2hkCFqQhcRB4DnAVwHYDSAaSIyOmC3uQCWq2oxgKkA/q/VgRJRatm/fz+GDx+OpUuX2h0KGWKpoV8CYJuq7lDVZgCvA5gcsI8C6G08PgPAXutCJKJUxF4uqSeWhF4AYLfp+R5jm9kTAKaLyB4AqwE8GOpAIjJTRGpEpObAgQNxhEtEqcKf0NnLJXVY1ctlGoBXVXUwgOsB/E5Ego6tqi+paomqlgwYMMCiUxORHdhtMfXEktDdAIaYng82tpn9BMByAFDVTwF0B9DfigCJKDWxhp56YknoGwCMFJFhIpIN303PqoB9GgBcBQAi8jfwJXS2qRBlsIEDB+K5557DhRdeaHcoZIj60aqqLSLyAIBqAA4Ai1R1i4g8BaBGVasAPALgZRF5GL4bpDNUVRMZOBHZq1+/fnjwwZC3y8gmMX1XUtXV8N3sNG97zPT4cwATrA2NiFLZyZMnsX37dgwdOhS9evWyOxwCh/4TUZw2b96MMWPGYO3atXaHQgYmdCKKC2+Kph4mdCKKCwcWpR4mdCKKC2voqYcJnYjiwoFFqYcJnYjicv755+OVV17B9773PbtDIQM/WokoLgUFBZgxY4bdYZAJa+hEFJeDBw/ik08+wcmTJ+0OhQxM6EQUlzVr1mDChAnYtWuX3aGQgQmdiOLCXi6phwmdiOLChJ56mNCJKC4cWJR6+NFKlCBzXXVYun43vKpwiGDauCGYV1Zkd1iWYQ099fCdIEqAua46LF7X0Pbcq9r2PFOS+jXXXIPKykr07dvX7lDIwCYXogRYun53h7ano2HDhuGWW25Bbm6u3aGQgQmdKAG8YdZ3Cbc9He3cuRPV1dVtTS9kPyZ0ogRwiIQtc9UGLsmbniorKzFx4kScOnXK7lDIwDZ0ogSYNm5IuzZ0s4rqepQVFyQ5ovBctW48+lYdTjT7JtsSAHeOL4za1s/JuVIP3wkiE6t6pswrKwqb0Pc2NnU2TMu4at14ZMUmeFu/awpSIKYbuOy2mHrY5EJk8PdM8bdz+3umzHXVxXW8gvzQNwsHhdluh4rq+nbJ3CzaDVx/DZ0JPXUwoRMZrO6ZUl46qkPb7eCO8G0h2g3clpYWiAiysphGUgWbXIgMVvRMcdW68c+Vm9DsDf+aFTUNYdvQ73z5U6zdfrjdtkQOSnKIhP39It3YBYC77roLP/zhDy2PieLHhE5kyBIgVOtDVuS8BletGxXV9RFru2aBCdsvVDIHvmv62XngOJbce2lM54hVpA+raeOGRHzt9773PS5ukWL4XYnIkNMt9H+HcNsBXzKf82ZdzMk8knCJ3lxudZfHcO38PbIdUb8RbNy4EStXrrQ0HuocJnQiwylPa8jtTZ7WsDdGK6rr0eTxdvhcrlo3xj75HobOfrvtJxYV1fUdPlck5aWjkOtsf1Mz1+nA/JuiN+8sWrSIKxalmJiaXERkIoBfAXAAWKiqCwLKfwngCuNpHoCBqppvZaBEiTYoPzdsTXvxuoaw3RDj8dCyjXG9zopvAoG6O7PaPpTyc514YtJ5MfWT93q97IOeYqLW0EXEAeB5ANcBGA1gmoiMNu+jqg+r6lhVHQvg1wDeTESwRIkUqraayfzNRUdOetq2NTZ58PCyjVG7arpq3Xij5iscafJiwoI1GTP6Nd3F8vF6CYBtqroDAETkdQCTAXweZv9pAB63Jjyi5PHXSv21Z21pRvM324P2c5xxJrr17ItWzyl49u8MKu+WfzYcPfLRevokPAe/Ci7vWwBHbm+0njoOz6HgLpHOfkOQ1b0nvE3fouVwcKI8duxH6NWrF/bv34/t24PjGzNmDHr06IF9+/aFXB6uuLgY3bt3x7zlf0bjnuDzZ581EovXNeDogX244/xeQeX7sgfjUdcWHP1mNyACd2MT5rzp+wBIpRGwXVEsCb0AgPld3wNgXKgdReQcAMMArOl8aETJd8mZguOfrUL3keMBbcXXi8uD9ulz9X3ofdGNaGn8OmR5v+sfQs+iq5FzzI3dIcr7T56NHuf+AKf3bcX+5Y8FlQ+87UnkDr8IpxrqcND1dFB53f2X4bLLLsPq1atx9913B5Vv3rwZRUVFeOONN/Dggw8Gle/cuRNDhw7F9k9Wo/FPrwWVD35wCRx5Z2DxK/8P//Xp8qDy8f/+jq+JxuFEVk5PAECTx5tyUxp0RVY3gE0FUKmqIe8SichMADMBoLCw0OJTE3Xel19+iUMfvIiB/QuRM+hcDLztyaB9nP183fm69R4YunzAOSjIz8Xbs2fgstbTOHT8dLvy7IHDff+eOSLk67PP8nUF7F4wOmT5lOW74ah6G5NGDsa7774bVD506FAAwKRJkzBy5Mig8jPPPBMAMHzcNdhz5oig8qycPABA3vlXIXvweW3be3bvhunjzsFvGzyAZKHvFT+FOLPbylNpSoOuKpaE7gZg7pA62NgWylQA/xDuQKr6EoCXAKCkpCRz5hGljPHiR75eJOLohixnDnKHXxR236ycvLDlexubkJ+fj58/OB2zlm8M2b/dkXdGxOM7evZBbs/w5VVfnkbvfuEn0SosLIxYcfq3aVdgzpv9w/bScfYtgLPvdzVuL4Df7xXk53VDY5MHzv7t+6mn0pQGXVUs3RY3ABgpIsNEJBu+pF0VuJOInAugD4BPrQ2RKDnmuurw4Za9AADJ6tyXV39yKysuwDO3j0WeM7Yewg7xzXZYkJ+L6eMLEWVMU6d63pQVF+Dpm4vQIzv2G8Eer0IEIbs6ptKUBl1V1L9aVW0RkQcAVMPXbXGRqm4RkacA1KiqP7lPBfC6agbN4E9dytL1u6Fe3wyC4gj+r+Efgg9ETqSBya2suCBk27J/hOnexiYMys9FeemodvtNWLAGif7P5I8t1n7wANB40oM7xxe2m5XylotC/46UXDFVQ1R1NYDVAdseC3j+hHVhESWfVxXaajQ/BNTQdy24IWj/UEm9IERiDidcovdLZpt0QYQ++IHM0+sCvuu2bMNulJzTl0ndZhwVQGRwiCB3RAkG3bcQ3Xr1a7c90LyyooQv9hxpoJNftCaZWJWXjop7sBPga4p5ctUWJnSbceg/kWHauCHIcnaHM/8siMPZbrsdyktHRU3YClgysKesuADPThmLGJv6QzIPUCJ7sIZOZJhXVoR92z/HW6veQc8L/xbOnLyETVsbC39tN9p0vO7GJpRXbmr3mnjPZ379iDmrM2pR666ACZ3IpDjnABb98bfYWvkL9O/f3+5w2iVZV60bs5ZtRKgpxDxexazlG9teY4VI66KGkp/rjL4TJRQTOpFJc3MzAMDpTL3kVFFdHzKZ+7Wqb9qCSG3hHZl8K9K6qIGcWYInJp0XfUdKKLahE5l4PL524Ozs7Ch7Jp8VvV4amzyYtWxjzG3uE0b0jbpPQX4uKm67gDdEUwBr6EQmqVxDj6XXSyxa4avJ13x1OOr9gSX3Xhq0klK2Q/AftzKBpyImdCKT5uZmiEhKrmRfXjoqbBt6PPzL2u061BR2cBMAy5e9o8QRuwZ2lpSUaE1NjS3nJgqnqakJJ06cSIkboqH45jDfjKYwqyt1Vq7TgadvLmLtO4WJyGeqWhKyjAmdKH2N/rd3cNLi5F6Qn4u1s6+09JhknUgJnTdFiUxcLhfmz59vdxgx+/nNY+B0WDVe1IfT4KYvtqETmbzzzjtYuXIlHn30UbtDiYm/aSTSJF+uWjcefasOJ5pjW8ya0+CmLyZ0IhOPx5OSXRYjiTbJl788luYZToOb3tjkQmTS3Nyckl0WrfDzm8dEXQSbN0TTG2voRCbpWEOPVeAi2IEcIkzmaY41dCKTTK6hA76kPn186GXpsrtJp2dtJHuxhk5kUllZiZaWFrvDSCj/6NDfr29ot9Zpk6cV5Ss6P2sj2Yc1dCITh8OBnJwcu8NIuHllRQjV29HTqniiakvyAyJLMKETmTzzzDN47rnn7A4j4ea66hCuw0tjExeqSFdM6EQmlZWVWLVqld1hJNzS9bvtDoESgAmdyCSTe7mYRVqJqE9e5t4UznRM6EQmmd7LxS/Uwtd+j9/IhSrSFRM6kUlzc3OXqKGHW/h6woi+7OGSxthtkcjE4XCge/fudoeRcP6ui0vX74ZXFQ4RWxfEJmtw+lwiojTC6XOJiLqAmBK6iEwUkXoR2SYis8Psc7uIfC4iW0Tk99aGSZQc999/PxYvXmx3GERxiZrQRcQB4HkA1wEYDWCaiIwO2GckgDkAJqjqeQAeSkCsRAm3ZMkS/Pd//7fdYRDFJZYa+iUAtqnqDlVtBvA6gMkB+9wL4HlVPQIAqrrf2jCJkqOrdFukzBRLQi8AYB5WtsfYZvZ9AN8XkbUisk5EJoY6kIjMFJEaEak5cOBAfBETJVBXGVhEmcmqm6LdAIwEcDmAaQBeFpH8wJ1U9SVVLVHVkgEDBlh0aiJreL1etLa2soZOaSuWhO4GYB6FMNjYZrYHQJWqelR1J4Ct8CV4orTR0tKCAQMGoFevXnaHQhSXWAYWbQAwUkSGwZfIpwK4I2AfF3w181dEpD98TTA7rAyUKNFycnKwfz9v/1D6ilpDV9UWAA8AqAbwBYDlqrpFRJ4SkUnGbtUADonI5wA+AlCuqocSFTQREQWLqQ1dVVer6vdVdYSqzje2PaaqVcZjVdVZqjpaVYtU9fVEBk2UCIcOHUJZWRnef/99u0MhigtHihIZjh8/jpUrV2LPnj12h0IUFyZ0IkNzczMAsJcLpS0mdCKDP6GzHzqlKyZ0IoPH41tLkwmd0hUTOpEhKysLw4YNQ+/eve0OhSguXOCCyDBmzBjs2MHhE5S+WEMnIsoQTOhEhg0bNuDqq6/G559/bncoRHFhQicyfPPNN/jwww9x4sQJu0Mhigvb0GPgqnWjoroeexubMCg/F+Wlo7gyegZiLxdKd0zoEbhq3Xhy1RYcOelp2+ZubMKcN+sAgEk9w7AfOqU7NrmE4ap1Y86bde2SuV+Tx4uK6noboqJEYg2d0h0TehgV1fVo8njDlu9tbEpiNJQMvXr1QlFREfLy8uwOhSguTOhhREvYg/JzkxQJJcvkyZOxefNmnH322XaHQhQXJvQwIiXsXKcD5aWjkhgNEVF0TOhhlJeOQq7TEbQ9P9eJp28u4g3RDPS73/0O48aNw8mTJ+0OhSgu7OUShj9hs7ti17Fnzx785S9/QVYW6zmUnpjQIygrLmAC70I4HzqlOyb0EOa66rBkXQPUeN4j24H5N7GZJdN5PB44HA44HMFNbUTpoEsndP8IUHeUHi0nmr14ZMUmABxMlMmam5tZO6e01mUTun/gUKS+5mbeVsWTq7Z0qk19rqsOS9Y3QDX6vgAwfXwh5pUVxXx86pxzzjkHP/jBD+wOgyhuorFmF4uVlJRoTU2NLecGgAkL1kStmUeT63TE3ONlrqsOi9c1dPgcE0b0xZJ7L40nPCLKQCLymaqWhCrrsrfzrRjp2eTx4qFlG1H81Htw1boj7rt0/e64zrF2++GoxyYiArpwk8ug/NxO19D9jpz0oLyyfRv7nS9/irXbD1ty/IrqerbdJ8GcOXPw2Wef4b333rM7FKK4pF1Cd9W68ehbdTjR7Gv7FgB3xtHWXF46KuY29B7ZDjgdWWhsCp6oy8/j1bbEa2UyBzhvTLI0NDRg165ddodBFLe0SuiuWjceWbEJ3tbv2v0VwOJ1De3ap2NJ8uaBQ6Fq6oFt165aNx5atjFifO7GJlzzzMf4cr+1CyRkicBV62YtPcHYy4XSXUwJXUQmAvgVAAeAhaq6IKB8BoAKAP7G3v9S1YUWxgnAl3zNyTwcf5Jfsq4BOQ7BKW/wawRAd2cWTnlaURBDj5Wy4oKoCR2A5ckcALyqaTMH+1xXHX6/vgGh3iaHCKaNG5KyPXeam5s5dS6ltagJXUQcAJ4HcA2APQA2iEiVqgYuvLhMVR9IQIxtOtr0oEDIZO4va/K0Aoh90YqRA3tYmrAFwM4FN7Q9d9W68UTVlpBNO/452EPFlworKrlq3fjnyk1oDnO9Ad8Hk/+DNp5mskRjDZ3SXSw19EsAbFPVHQAgIq8DmAwg6Svp+m9kNv7592j++st2ZY5e/dCv1Pd5cuTjV+E5+FW78m59BqHvVfcCAA5/8CJaGr9uV+4cOAwVPe5DWXEBZs6ciX379rUrHzduHN6fOxfXPPMxPl30FFpPHW9XnjOkCGeMuxkAcOCtn0O97ZNy7vCL0OvCv4VqKw688e8AgO5OB27c8gIA4KabbsI999yDa76fjwEXXBH0u/cYfTn2jv4xDh06hBkzZrRt//roKXzx9bfIO/8a5H3/UnzVsBvTp8zB35zVG2ed0b1tvwceeAClpaXYunUrHnnkkaDjl5eX40c/+hE2bdqEuXPnBpU/9thjuPjii7Fu3TrMnz+/XdnXR0/hm1E3IavfOTjVsBnf/uWtoNf3uWomnH3ORtP2GhyrfRvPVALL5ufi3LN6AwBefPFFDBo0CG+++SZeeeWVoNe/9tpr6NOnD5YsWYLXX389qHzFihXo3r07Fi5ciJUrV7YrExFUVVUBAH79618H3fTMy8vDsmXLcNFFF3FiLkprsST0AgDmPnd7AIwLsd8tIvIjAFsBPKyqQf30RGQmgJkAUFhY2OFgy0tH4ZEVm9B66hi8x9vfdBTHdzWr1qZvg8qzuvdse+w9GVzu6Nmv7RvA/v37sXfv3nblhw/79n9/1uW49t2f45MtjfB4W7875+nvErz3xBFoS3O717ee/i5ReI8fhgjQJz8Pe/f69vv22299+7W2Iqup/bEBoLXZV/v2er3tYvti37fwtipaPacAANraAs+3h/DFiSNoPatX237+ROXxeIJ+NwBoavL97s3NzSHLT5061fZvYHn918fQa7gHOQDU0xx0bX2/QIvvH8+ptvKG40DvVt91a2nxlR8/fjzk+b1eb9t1ClXuH09x9OjRoHIRaXvc2NgYVN6zp+9vY968ecFxE6WRqAOLRORWABNV9afG878DMM7cvCIi/QAcV9XTInIfgCmqemWk48Y7sMhV68a/vLEZp1tao+/cQQX5uVg7O2LY7eKYtXxjyLbiaPrkOfH4jeeFbRYJNYpVAHTkVIHNOYk0bPbbHYrN7NkpY9uuQ6T2dwAQAe4cl3pNNUTJFGlgUSw1dDeAIabng/HdzU8AgKoeMj1dCOA/OhpkrPwzILpq3ShfsREeC/N6Rxat8CehR5ZvRIRmY+Q5s/Dzm8d0qE07sAdOR5M5kNwVlTrTp//JVVtQVlwQ00haVbTtw6ROFCyWhL4BwEgRGQZfIp8K4A7zDiJytqr6G0Aifu8AAAsoSURBVJ0nAfjC0ihDCJzaNtaJtqIdM54YQiWjzg7Z9x873ikKkrmiUnnpKJRXboIn0idbGP5FuDsykpZJnSi0qAldVVtE5AEA1fB1W1ykqltE5CkANapaBeB/i8gkAC0ADgOYkcCYQzIn+FADe3KdWW29WkIZObBH3OeeV1aUkOTiqnVbNpo1kfzX/clVW9oSdH6uE09MOi/mDyVvB+cUWryuATsPHOc8N0QmXWpyrnCDfro7BP8z//qkxhJNZ0ebduR+QKJFG5T17JSxeGT5pg4n9VA4QyVlOk7OZXh/1uWYMKJvu20TRvRNuWQ+11XX6akDUmm6gLLiAkwfH75X08PLNmL88D6WnGvxugaM/rd3OKEZdUldqoaeLkbMWW1JbRXwNSUdONbcNlipT54TN4w5Gx/9z4G4BiKFGsRU89VhLF2/G17ViKNBh85+O+Kxc42Ru1b85s4sQcVtF6T8yFqijupsLxdKsliSub9JJdpCHYFNTEdOetrdwI11lCxgdNVcthH+OxHuxqagphT/aFCg4zctmzytQXPMxzuPvKdVMWvZxpgT+lxXXUwfSkSprEs1uaQLh2kgTCi5TkdbL5ay4gI8fXMRCjrRTdE/rUA0c97cjFh7icY7/3tgLPPKiiI210TSCt+9iGjufPlTLF7X0PZB6v9Qmuuqi+u8RHZhQk9B08YNCVvWJ88ZtEpSWXFBp2+AxtLmHqmXUKDONBkFxjKvrAi7FtyAZ6eMRa6zY3+y0e5FuGrdYfeJ90OJyC5scklB/q/65lGTuc4sPN3BAUodETgQyeo53TsTi5+537+/eaSzIn0zseo+BlGyMKGnqHj6tk8Y0TfuJHzFuQPaHicqmRfEMKLU3JwUjvnadDbWSN9MojV9EaUaNrlkkCX3XhrULdMRY0766H8OtD22IpmHSoblpaOQ63QEbffvWZCfG/Oi236hfmczZxYidmGMNEVCpKYvolTEGnqGCRw56ap14+FlG6N2BbS633qoZGieo8bKudv9v3Oo2rqnFShf0X69V7Py0lEoX7EJnoAZwSaM6MteLpR2mNAzXFlxAWq+Oowl6xoiJvX8vI4t7GAeidqRLn+Bc/BYacm9l2Lsk+8FLRDiaVU8UbUl5HlrvjoclMydDsFtJfH1rCGyE5tcuoB5ZUX45ZSxEbs2Hj/V0tY0EakJA/D90ZjbueeVFWH709dj14IbsP3p622t2YZbyLuxyYPip95r1/ziqnVjSYg+7v4Fv4nSDUeKdkGharGxSnRvm86KNhrVLNshYZfMS+Z88kQdwZGi1M7RDiTzdBtC3yfP2TbjYzSR1j9N5nzyRFZhk0sXlJcd3NMkHE9rejU/PH7jeZ0+hiC588kTWYUJvQs62Rx63pdwUmnmxmjKigvQp4M3eM0EwJ3jC9PmGwmRGRN6F5TKy9lZ4fEbz4Mz1g74Jj2yHfjllLHsrkhpiwm9C+rICEhnlqRd80NZcQEqbr2gwzX1E81e1swprfGmaBc0bdyQmKakFUFa3RA1C9XfPdyKVUSZgjX0LmheWVHUvuZOh+CXt49Ny2QezvuzLkd+bvzt60Spjgm9i1py76V4dsrYkAmuT54TFbemZ808micmhe8Fw6m4KN2xyaULS+Qw/FRVVlyAFTUNIScguzPOhTSIUgVr6NTlLLn3UkwfX9h2c9ghgunjC9m7hdIeh/4TEaWRSEP/WUMnIsoQTOhERBmCCZ2IKEMwoRMRZQgmdCKiDGFbLxcROQDgqzhf3h/AQQvDsQrj6phUjQtI3dgYV8dkYlznqOqAUAW2JfTOEJGacN127MS4OiZV4wJSNzbG1TFdLS42uRARZQgmdCKiDJGuCf0luwMIg3F1TKrGBaRubIyrY7pUXGnZhk5ERMHStYZOREQBmNCJiDJEyiZ0EblNRLaISKuIhO3eIyITRaReRLaJyGzT9mEist7YvkxEsi2Kq6+IvC8iXxr/9gmxzxUistH0c0pEyoyyV0Vkp6lsbLLiMvbzms5dZdpu5/UaKyKfGu/3ZhGZYiqz9HqF+3sxlecYv/8243oMNZXNMbbXi0hpZ+KII65ZIvK5cX0+FJFzTGUh39MkxTVDRA6Yzv9TU9ldxvv+pYjcleS4fmmKaauINJrKEnm9FonIfhH5a5hyEZHnjLg3i8iFprLOXy9VTckfAH8DYBSAjwGUhNnHAWA7gOEAsgFsAjDaKFsOYKrx+AUA91sU138AmG08ng3g/0TZvy+AwwDyjOevArg1AdcrprgAHA+z3bbrBeD7AEYajwcB2Acg3+rrFenvxbTP3wN4wXg8FcAy4/FoY/8cAMOM4ziSGNcVpr+h+/1xRXpPkxTXDAD/FeK1fQHsMP7tYzzuk6y4AvZ/EMCiRF8v49g/AnAhgL+GKb8ewDvwLZA1HsB6K69XytbQVfULVa2PstslALap6g5VbQbwOoDJIiIArgRQaez3WwBlFoU22TherMe9FcA7qnrSovOH09G42th9vVR1q6p+aTzeC2A/gJAj4Top5N9LhHgrAVxlXJ/JAF5X1dOquhPANuN4SYlLVT8y/Q2tAzDYonN3Kq4ISgG8r6qHVfUIgPcBTLQprmkAllp07ohU9U/wVeDCmQzgNfVZByBfRM6GRdcrZRN6jAoA7DY932Ns6wegUVVbArZb4UxV3Wc8/hrAmVH2n4rgP6b5xtetX4pITpLj6i4iNSKyzt8MhBS6XiJyCXy1ru2mzVZdr3B/LyH3Ma7HUfiuTyyvTWRcZj+Br5bnF+o9TWZctxjvT6WIDOngaxMZF4ymqWEA1pg2J+p6xSJc7JZcL1vXFBWRDwCcFaLoUVVdmex4/CLFZX6iqioiYft9Gp+8RQCqTZvnwJfYsuHri/ovAJ5KYlznqKpbRIYDWCMidfAlrbhZfL1+B+AuVW01Nsd9vTKRiEwHUALgx6bNQe+pqm4PfQTLrQKwVFVPi8h98H27uTJJ547FVACVquo1bbPzeiWUrQldVa/u5CHcAIaYng82th2C76tMN6OW5d/e6bhE5BsROVtV9xkJaH+EQ90O4C1V9ZiO7a+tnhaRVwD8UzLjUlW38e8OEfkYQDGAN2Dz9RKR3gDehu/DfJ3p2HFfrxDC/b2E2mePiHQDcAZ8f0+xvDaRcUFErobvQ/LHqnravz3Me2pFgooal6oeMj1dCN89E/9rLw947ccWxBRTXCZTAfyDeUMCr1cswsVuyfVK9yaXDQBGiq+HRjZ8b16V+u4yfARf+zUA3AXAqhp/lXG8WI4b1HZnJDV/u3UZgJB3wxMRl4j08TdZiEh/ABMAfG739TLeu7fga1usDCiz8nqF/HuJEO+tANYY16cKwFTx9YIZBmAkgL90IpYOxSUixQBeBDBJVfebtod8T5MY19mmp5MAfGE8rgZwrRFfHwDXov031YTGZcR2Lnw3GD81bUvk9YpFFYD/ZfR2GQ/gqFFpseZ6Jepub2d/ANwEXzvSaQDfAKg2tg8CsNq03/UAtsL3Cfuoaftw+P7DbQOwAkCORXH1A/AhgC8BfACgr7G9BMBC035D4fvUzQp4/RoAdfAlpsUAeiYrLgCXGefeZPz7k1S4XgCmA/AA2Gj6GZuI6xXq7wW+JpxJxuPuxu+/zbgew02vfdR4XT2A6yz+e48W1wfG/wP/9amK9p4mKa6nAWwxzv8RgHNNr73HuI7bANydzLiM508AWBDwukRfr6Xw9dLywJe/fgLgZwB+ZpQLgOeNuOtg6sFnxfXi0H8iogyR7k0uRERkYEInIsoQTOhERBmCCZ2IKEMwoRMRZQgmdCKiDMGETkSUIf4/yDSxu6xot4IAAAAASUVORK5CYII=\n"
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "tmYzmE5NTf_p"
      },
      "source": [
        "## S-Learner"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "oLPw7a27TdTA",
        "outputId": "19e713a5-bac7-4c09-f6e0-3e863fe4dc3a",
        "colab": {
          "base_uri": "https://localhost:8080/"
        }
      },
      "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": 14,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stderr",
          "text": [
            "/usr/local/lib/python3.7/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"
          ]
        },
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "RandomForestRegressor(max_depth=4)"
            ]
          },
          "metadata": {},
          "execution_count": 14
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "IJmgnebQTfkB",
        "outputId": "4d59898f-a47e-4f31-8ee7-e15fb138374c",
        "colab": {
          "base_uri": "https://localhost:8080/"
        }
      },
      "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": 15,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "ATE： 0.6686144361959908\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "sv1NNZi8TjUJ",
        "outputId": "1b660455-1f1e-43e2-e548-1aec5d873a02",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 283
        }
      },
      "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": 16,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "[<matplotlib.lines.Line2D at 0x7f04db8beed0>]"
            ]
          },
          "metadata": {},
          "execution_count": 16
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dfXxU5Znw8d+VSQIDKiHyngBBBAQLJTUfEek+Wm0FX4opqwLF3Wqr2K22D9WmQqUttSh02aU+2/q0ik/rurKAujBS0aWIoiuVlEjAFDW8KjAgRGJQSSCT5H7+mDPxJPM+mfdc388nH2bOfeacK2fCNffcb0eMMSillMp8OakOQCmlVHxoQldKqSyhCV0ppbKEJnSllMoSmtCVUipL5KbqxP369TMlJSWpOr1SSmWkt9566yNjTP9AZSlL6CUlJVRVVaXq9EoplZFE5INgZdrkopRSWUITulJKZQlN6EoplSU0oSulVJbQhK6UUlkibEIXkT+IyAkR+VuQchGRfxORfSLytoh8Kf5hKqWUCieSGvqTwLQQ5dcCo6yfucDvuh6WUkqpaIUdh26MeV1ESkLsciPwlPGuw7tNRApEZLAx5licYlRKpSljDL///e85dsz/v/v9999P7969efnll3n99df9yn/2s5+Rm5vLhg0bqKys7FCWk5PDokWLAFi3bh3V1dUdyp1OJwsWLABg9erVvPPOOx3K+/Tpw3333QfAU089xb59+zqUDxgwgHvuuQeAFStWcPjw4Q7lxcXFzJ07F4BHH32U48ePdygfOXIk3/rWtwBYvnw5DQ0NHcrHjh3L7NmzAVi6dCmNjY0AfPvb3yahEyqNMWF/gBLgb0HKXgC+bHu+GSgLsu9coAqoGjZsmFFKZbYzZ86Y3/zmN8bpdBoR6fBTV1dnjDHmgQce8CsTEXPmzBljjDE/+MEP/Mry8vLaz3H77bf7lRcWFraX33TTTX7lw4cPby+fOnWqX/kXvvCF9vIpU6b4lU+ePLm9fMKECX7lX/va19rLR4wY4Vc+Y8aM9vJ+/fq1b3/ttde6fM2BKhMkV4uJ4AYXVg39BWPMFwKUvQAsNca8YT3fDNxvjAk5DbSsrMzoTFGllIqOiLxljCkLVBaPqf9uYKjtebG1TSmV5RobG6mqquKiiy5iwIABST//QlcNqyoP0xqkYpoDLJ85kfLSouQGliLxGLa4HvhHa7TLZcApo+3nSnULhw8f5oorrmDz5s1JP/dCVw1PbzsUNJkDtAHz1uzEVd096piRDFtcBbwJjBGRIyLyHRH5roh819rlReAAsA9YAXwvYdEqpdJKS0sLALm5yV/n7z8rD0W874K1bycwkvQRySiX2WHKDXB33CJSSmWMVCV0V7Wbtijub9/kacNV7c76phedKaqUilmqEvqi9bujfs2yjbUJiCS9pGw9dKVU5vMldIfDkbBzuKrdLNtYy9GGJoYUOKmYOoaGJk/Ux3E3NCUguvSiCV0pFbPRo0fjcrm45JJLgMDJN9ZmDl+np527oYl5a3bGdDyx4svmZpeIxqEngo5DVyo7uKrdLFq/O2itedSA3nzS5OH4p80dtvfKy+HhGRP8EmygRB6NUQN6s/fE6YBlRQVOts6/KuZjp4NEj0NXSnVT//7K2yz4/Tocg8fgcJ4bcJ9gybXR08a8NTtjrnEHcutlw3j1vbqg5UezvNlFO0WVUjH756f/m2PPLsJz8nD4nZNgcfn4kEm7jzMvidEknyZ0pVTMPvrUmzwlJ3GdopEqKnACMMT6N5BPzniyepKRJnSlVMzEtHkfpDih5+UIFVPHAFAxdQwSZL82k93DFzWhK6Vi1tbmHbYoOclJJb3yvOcRW8YucOax7OYvtneulpcWEWqoRza3o2unqFIqZn2dudQBSOISukPgX2+JboGtogJn0HHnoZpkMp3W0JVSMVvw7RkMnfMwuX0Gdthe4MxjysjCDjVpn7wIs45DhFsvG8b+JddHPXa8YuoY8hz+J7c3zWQjraErpWJ229VfpKCwX9wmE8WL7/y/+NNuPm70jo8vcOaxaPrFKY8tkTShK6ViduDAAc7u2c6fv38DvXv3TnU4HZSXFmV18g5Em1yUUjF75ZVXmDVrFh9//HGqQ1FoDV0p1QXJWJwrVSY9tMlvuYJApowsZOWdk5MQUXia0JVSMWttbQVSc4OLeHFVu1mw9m2aPG0xvX7r/npK5m8IWp7MhJ+574JSKuVSeceieOjqQmCRsCd8hwizJw1lcfn4hJxL29CVUjHL5ISejGTeWasxPL3tEJMe2pSQ42tCV0rFbM6cOVRWVtKrV69UhxIVV7WblUlO5nbHP21mzoo3437czPtYVUqljUGDBjFo0KBUhxFUsBtuPLCuJuTyAMmwdX993I+pNXSlVMzeeustnnjiCVJ1o5xQvJ2dNbgbmjB473a0YG0Nc1a8yenm1lSHlxCa0JVSMXv++eeZO3cuEmiOf4ot21hLk6dj4m7ytMZUM+6Vl8MjMyfy/tLreWTmRJyRrl+QZNrkopSKWUtLS9p2iMZjVUVnXg5LOt0mL5IZqF9bviXonZp8pows7HJ8naXnO6GUyggtLS1pO6moV74jqqYVAQ4uvT4u595075UsdNWwctuhgG31iRqbrgldKRWz1tbWtK2hN0bZTn5uj/h+MC0uH5+w8ebBpGdDkFIqI6Rzk0u03bSfnM38jtKIErqITBORWhHZJyLzA5QPF5HNIvK2iGwRkeL4h6qUSjc/+clPqKysTHUYATnSsKM20cImdBFxAI8C1wLjgNkiMq7Tbv8CPGWMmQA8CCyJd6BKqfQzcOBARo8eneowApo9aWiqQ0i6SL4rXQrsM8YcABCR1cCNwDu2fcYB91qPXwVc8QxSKZWeNm7cyOHDh7njjjtSHYqfxeXjg3ZKBjJqQGLXc3dVuxN+w41ImlyKgMO250esbXa7gBnW428A54rI+Z0PJCJzRaRKRKrq6upiiVcplUZWrVrF4sWLUx1GUNEk8033XpmwOFzVbiqe29WezAEamjxUPLsLV7U7bueJV2/Gj4DfishtwOuAG/DrYTDGPA48DlBWVpZ+U8uUUlFJ505RCH6z6KICJ1vnX5W0OJZtrMXT6p/yPG2GZRtr41ZLj6SG7gbsjVHF1rZ2xpijxpgZxphS4AFrW0NcIlRKpa10HocO3ptFO/M6xufMcyT9RtGhJjnFYwKUTyQJfTswSkRGiEg+MAtYb99BRPqJiO9YC4A/xC1CpVTaSvcaenlpEUtmjKeowIngrZkvmTE+6fcaHVLgjKksWmHfCWNMi4jcA2wEHMAfjDG7ReRBoMoYsx64ElgiIgZvk8vdcYtQKZW20nlikU863Cy6YuoYKp7b5dfskpcjcf22IKlaJa2srMxUVVWl5NxKqfior6/H4/EwcODAVIeS9uI1ykVE3jLGlAUqS++PVqVUWissjP8CU9kqGd8UdOq/UipmTz31FH/84x9THYayaEJXSsXsySef1ISeRjShK6Vilu6jXLobTehKqZhpQk8vmtCVUjFL94lF3Y0mdKVUzLSGnl70nVBKxex//ud/Uh2CstGErpSKWe/eiV1yVkVHm1yUUjFbunQpK1euTHUYyqIJXSkVsxUrVvDSSy+lOgxl0YSulIpZJizO1Z1oQldKxUxHuaQXTehKqZhpQk8vmtCVUjHTiUXpRT9alVIxq6urI1X3VFD+NKErpWImIohIqsNQFm1yUUrF7O6772bdunWpDkNZNKErpWL2+OOPs3379lSHoSya0JVSMdNRLulFE7pSKiZtbW0AmtDTiCZ0pVRMWlpaAE3o6UQTulIqJq2trZxzzjn07Nkz1aEoi360KqVi4nQ6+fTTT1MdhrLRGrpSSmWJiBK6iEwTkVoR2Sci8wOUDxORV0WkWkTeFpHr4h+qUiqdnDp1itmzZ7Np06ZUh6IsYRO6iDiAR4FrgXHAbBEZ12m3hcAzxphSYBbwf+MdqFIqvTQ2NrJ69WoOHDiQ6lCUJZIa+qXAPmPMAWNMM7AauLHTPgY4z3rcBzgavxCVUunIN8pFF+dKH5Ek9CLgsO35EWub3SLgVhE5ArwIfD/QgURkrohUiUhVXV1dDOEqpdKFDltMP/HqFJ0NPGmMKQauA/5DRPyObYx53BhTZowp69+/f5xOrZRKhdbWVkATejqJJKG7gaG258XWNrvvAM8AGGPeBHoC/eIRoFIqfQ0ePJjevXunOgxlieSjdTswSkRG4E3ks4BvdtrnEHA18KSIjMWb0LVNRaksduGFF3L0qHaXpZOwNXRjTAtwD7AReBfvaJbdIvKgiEy3drsPuFNEdgGrgNuMrnqvlFJJFVEbujHmRWPMaGPMSGPMQ9a2nxlj1luP3zHGTDHGfNEYM9EY8+dEBq2USr29e/dyww036PK5aURniiqlYlJfX8+GDRv46KOPUh2KsmhCV0rFRIctph9N6EqpmOjEovSjCV0pFROtoacfTehKqZj07NmTUaNG6Tj0NKIfrUqpmEyZMoU9e/akOgxlozV0pZTKEprQlVIxef3117nyyivZv39/qkNRFk3oSqmYHD9+nNdee40zZ86kOhRl0YSulIqJjnJJP5rQlVIx0YSefjShK6ViohOL0o8mdKVUTPr27cuXvvQlevbsmepQlEVStcptWVmZqaqqSsm5lVIqU4nIW8aYskBlWkNXSqksoQldqQRZ6Kph5IIXKZm/gZELXmShqybVIcXVypUrueSSS/jkk09SHYqyaPe0Ugmw0FXD09sOtT9vNab9+eLy8akKK66OHTvGjh07EJFUh6IsWkNXKgFWVR6Oansmam1tBXTYYjrRd0Ipm4WuGlZVHqbVGBwizJ40NKYadWuQwQbBtqeSq9rNA+tqON3sTdACzLlsWNjfW8ehpx+toStl8TWT+JKur5kklrZvR4hmCFe1O+YY481V7ea+Z3e1J3MAAxH93joOPf1oQlfKEo9mEle1m9EPvBiyJr5sY23UsSXKso21tLYFjjXc7z1s2DCuuOIKcnI0jaQLfSeUsnS1mcRV7eaHa3bS3Bp6/6MNTVHHliihYgn3e99+++1s2bIlzhGprtDGL6UsOQJBKquUzN9AUYGTiqljKC8t6lDmqnazbGMt7ggT9ZACZ1dDjZuCXnl83OhJdRgqTjShK2XpkZtDk6ctaLm7oYl5a3bywLoaHvrGeMpLi3BVu1mwtoYmT2vQ19kJUDF1TNDyOSveZOv+er/tzrwclsyY4Pdh0lVd6aN9+OGHee6559ixY0f8AlJdogldKcuZEMnc7nRzK/PW7GTemp0I3k7ESBkImpSDJXOAJk8b967ZCSFeH4tTTaFr565qd9Dzffjhhxw8eDBusaiui6gNXUSmiUitiOwTkfkByn8tIjutnz0i0hD/UJVKrFiaQmKp4HYePbLQVcMFCzYETeY+bcS/Q7WPMy9keajztba26pDFNBM2oYuIA3gUuBYYB8wWkXH2fYwxPzTGTDTGTAR+A6xNRLBKJVLF1DE48xI/BG9l5eczSH1DJYO13XcWzw5VV7Wb080tIfcJ1S/Q0tKiCT3NRPJuXArsM8YcABCR1cCNwDtB9p8N/Dw+4SmVPL6mhYqnXqM57xwkx0Fr0ye0nTntt29uwUBEcmhtPEXb2Ua/8ry+gwFoPd1AW3OnpCiCq9rNovW7OVl3gjZPx1u4SU4OuX0GAtDy6UlMS3N7Wc/cHA4dOsSwYcMAcLvdfreA69GjB8XFxQAcPnyY5ubmDuVOp5MhQ4awbGMtTfXHMW0d2/8lrwe55xR6z3/qOEV3PQFAH2cu379qFDMmXcjAgQNpaWnRIYtpJpKEXgTYB6QeASYF2lFEhgMjgFe6HppSyTe212fse2QO3/3t87x0GD7Z/jyfvLnGb7+h9/4XkteDU39Zw6dvre9YKDkM/7F328ev/TunazZ1LO7Rm4rCZ/C0GepffozG2jc6lDvO7Ufx954E4ORL/8aZg291KJ/26ljeecdbn5o5cyZbt27tUH7ppZdSWVkJwNe//nV27drVofzqq6/m5Zdf5mhDE8dX/4SWhg87lDtHXcaAGQsBOPbUvbQ1ngLgKPC9/wNrrptBUfmPeNVM5KOWVymZv4EpIwtZeedkv+ukkive35dmAc8ZYwJ2+YvIXGAu0F7DUCqdLH72LwCs3bKDcy68lLGTv8r+wiK/tnJxeP/r9L74SvIHXdipUCgqcOJuaOKcCdfQc1jHKfSS48BjtbGce8kNOC+8tGN5Xo/2x+dd+g16j7uiQ/nSuZ8/X7hwIXV1dR3Kzz///M9/n8WL+fjjjzuUDxo0CPD2GZz+yrdpa+5Yw889r1/748KvfhfT2rHjdH/hIN7fX0/+gAvod8N9AGzdX8+cFW9qUk+xSBK6Gxhqe15sbQtkFnB3sAMZYx4HHgfvDS4ijFGppFjoquHPfzsKgPToRasxfJg3mLvumETZ8MKAY817DB5Nj8GjO2xz5jmomDqGeWt20rN4LBSPDXrOnkO/0PF/VyfOkokdnjtEmD79uvbn06ZNC/k73XDDDUHLKqaOYcHp5pBDLnuP/buQx8/vX9L+OFynrkq8SBrAtgOjRGSEiOTjTdrrO+8kIhcBfYE34xuiUsmxqvIwWLVRycntsL28tIit86/i/aXXc+tlwb9d9sjNYcmM8XEfL+4ze1KI7B+l8tIilswYT1GBEwHyHboMbqYLW0M3xrSIyD3ARsAB/MEYs1tEHgSqjDG+5D4LWG1SdU87pbqo1RhMq3fUh69JxbfdzrcKoX29c8CvHXnUgN7sPeHfoRqJogInJec7+cuB+g6Tf57edqj9vKMG9GbTvVfiqnbziz/tDjjjs1deDg+HmJBUXlrUoWzCz/+bT85GNklKpR+9p6hSlpELXqTx2F4+2/kSfS6f3d6W7BBh/5Lrwrw6sJL5G6J+je+Dwbc2TKj/oQPPzae+0YMnzPoxPgXOPBZNvzjkN4hQE5xC0Y7R5NB7iioVgdmThtJj0IWcP+37HToGu9LMURTlZKVRA3q3J8VlG2vDTlw6/mlzxMkcoKHJQ8Wzu0Iu4bvyzslhm5Y602SeHjShK2VZXD6eOZcWt/+ncIhwawQ3egilYuoYciJomu7bK49HZk5k071Xtm9L1KqMnjYT0YzTSJcNdohoMk8TOs1LKZvhH1Vy8Fd34Xa7GTJkSJeP52va+Mnat2m0rRUTyWJbQ6yhj4kQyXEjXTY4nh21qms0oStl4/F4Oxbz8/PjdszOHY+Rqpg6JmwbelcsdNWE/PbhEAmZ1EVgzqSufYNR8aVNLkrZ+KbJ5+WFXrQqGcpLi/j1zIkJG0749LZDzFkRfJRxqJq3My+Hg0uu12SeZjShK2XjS+jxrKF3RXlpEXseuo73l17PIzMnRvQfNi9HmDKyMKLjb91fH7SDdHH5+IDHyQGWzJgQ0fFVcmlCV8omnWronZWXFrF85kScecH/2xY481h28xfbR6o8MnNi2BUk563Z6V2PZekrfsl95Z2TeWTmxPbJR0UFTpbPnJiwiVOqa7QNXSmbSZMm8cMf/jBt72QfbXu8b9951s0xQnE3NLFgbU2H18VyTpU6WkNXyuaaa65h+fLliGTPNPjy0qKIx5Q3eVrjfhMNlTya0JWyaWpq4vTp2Kbrp7Ng7eGBJGr8u0o8TehK2fz4xz/O2qWdfe3h4cRyKz6VHjShK2XT3NycNiNcEqG8tIiCEPcR9S39qzKTJnSlbDweT1qOcImnRdMvJi/IegRfGtZHO0AzmCZ0pWyyvYYO3lr6spu/SK8Awx+37q9noasmBVGpeNCErpRNd0jo4E3qTba1ZewiXZRLpR8dh66UzU033cQVV1wRfscMt9BVE3SNmEgX5VLpRxO6Uja33HJLqkNIilC1cEcWjcHvbrTJRSmburo6Tp48meowEi5ULVyXw81cWkNXyubmm28GYMuWLakNJMFCLY2rKyhmLq2hK2XTXTpFg9XCo7ntnEo/WkNXyqa5uTnrx6HD57XwVZWHaTUGhwizJw3V2nmG04SulI3H4+kWNXTwJnVN4NlFm1yUsukuTS4qO2kNXSmbH//4x/Tv3z/VYSgVE03oStncfvvtqQ5BqZhpk4tSNnv27OHEiROpDkOpmESU0EVkmojUisg+EZkfZJ9bROQdEdktIv8Z3zCVSo7Jkyfzy1/+MtVhKBWTsE0uIuIAHgW+BhwBtovIemPMO7Z9RgELgCnGmI9FZECiAlYqkbRTVGWySGrolwL7jDEHjDHNwGrgxk773Ak8aoz5GMAYo99ZVUbShK4yWSQJvQiwr+RzxNpmNxoYLSJbRWSbiEwLdCARmSsiVSJSVVdXF1vESiWIMabbTCxS2SlenaK5wCjgSmA2sEJECjrvZIx53BhTZowp06FhKt20tLQAaA1dZaxIhi26AfvCD8XWNrsjQKUxxgMcFJE9eBP89rhEqVQSiAgrVqzgkksuSXUoSsUkkhr6dmCUiIwQkXxgFrC+0z4uvLVzRKQf3iaYA3GMU6mEy83N5Y477qC0tDTVoSgVk7AJ3RjTAtwDbATeBZ4xxuwWkQdFZLq120bgpIi8A7wKVBhjsn9RaZVVmpub+etf/4r276hMJSZFt5sqKyszVVVVKTm3UoG43W6Ki4t57LHHmDt3bqrDUSogEXnLGFMWqExniiplaW5uBrRTVGUuTehKWTweD6AJXWUuTehKWbSGrjKdJnSlLL6ErhOLVKbShK6UpaSkhFWrVlFWFrC/Sam0p+uhK2UpLCxk1qxZqQ5DqZhpQg/BVe1m2cZa3A1NOERoNYaiAicVU8dQXtp5ORuV6U6ePMnOnTspKyujT58+qQ5Hqahpk0sQrmo3C9bW4G5oAqDVGq/vbmhiwdoaXNWdVz9QmW779u189atf5d133011KErFRBN6EMs21tLkaQ1Y1uRpZdnG2iRHpBJNO0VVptOEHsRRq2Yea7nKPDpsUWU6TehBDClwdqlcZR6dWKQynSb0ICqmjsGZ5whafvpsi7ajZxmtoatMp6NcAljoqmHltkPYly0T6PC8ocnDgrU1ADriJUtcffXVvPDCCwwcODDVoSgVk26d0O3DEsMJtCalr3NUE3p2KC4upri4ONVhKBWzbtvk0nlYYqzcDU2UPvhnbX7JAnv37uX5559vb0tXKtN024QealhitD5u9FDx3C5N6hnO5XJRXl7e3pauVKbptgk93sMOPa1Gx6ZnOO0UVZmu2yb0aIcdOnIk7D46Nj2z+RJ6bm637lpSGazbJvRwwxLteuc7+Nebvxh2v3AfEgtdNYxYsIGS+d6fcT99SZtp0khzczN5eXmIhP/wVioddduqiG9kSrBRLlNGFrLyzskdtj366l72njgd8Hh5DqFi6pig51voquHpbYc6bGv0tDFvzU7mrdnpt38OsHzmRB1Bk0Qej0ebW1RG05tER+lry7cETep25/Vw8MnZrne6PqJJPWkOHjyI2+3my1/+cqpDUSqoUDeJzriE7qp288C6Gk43e5OlAHMuG8bi8vEB9w+WgHMEvjkp+OsiiePeZ3bSluDLl5cDex++PrEnUUpljFAJPaPa0F3Vbu57dld7MgfvhJ+ntx2iZP4GRi54kYWumvaySQ9tClqbbjPe19n3j8ayjbUJT+YAnja0nT1J3njjDdatW5fqMJSKWUYl9GUba2kNkUVbjWlP7iXzN3D80/DjiVdVHo4plmSOaNHhkMnx2GOPcd9996U6DKVillEJPRFJtNUYLpi/IepacDJXW9ThkMmhnaIq00WU0EVkmojUisg+EZkfoPw2EakTkZ3Wzx3xDzVxSbQNmLdmZ1RJvWLqGCIYmh5SpKPjdKne5GhubtaErjJa2IQuIg7gUeBaYBwwW0TGBdh1jTFmovXzRJzjBLxJNJIJPrGKpmmjvLSI5bdMpFdebF9ycgR+fctE3l96ffvPlJGFAffVpXqTQxO6ynSRjEO/FNhnjDkAICKrgRuBdxIZWCC+4Xt3P7qe1jMdOztz8vLJHzgSgOa692k727GZIie/J/kDRnjLTxygrflsx/IeTo5SAsCOHTs4c+ZMh/K+ffsyduxYwHvvSY/HwwDgiWu9NxPu168fo0ePZqGrhv+39s+YTm39jt4F5PUdjDEGOV7L3CsvZEDTB/zlLx8AMGTIEFbeOZn/2v4BCx5byydnPu/4/RD40R+PA1dx7bh+vPXWW37XpqSkhCFDhtDY2MjOnf7j2keOHMnAgQP59NNPqanx7wgePXo0/fr149SpU+zevduvfOzYsfTt25f6+nree+89v/KDLX35ly2H+eijj/DUH+XcHg6+83cX8L9G9wdgwoQJnHPOOXz44YccOHDA7/WlpaU4nU7cbjcffPCBX3lZWRn5+fkcOnSII0eO+JVPmjQJh8PBwYMHOXbsmF/55ZdfDsC+ffs4ceJEhzKHw8GkSZPaJxYplbGMMSF/gJuAJ2zP/wH4bad9bgOOAW8DzwFDwx33kksuMbG65MtXGbwDXNp/cguLzfD7XzDD73/B9Ci+2K88f9Co9vK8ARf4lfcYNsFcvmSzMcaYCy+80K98+vTp7ecfOHCgX/k3v/nN9vJevXr5ld91113GGGNaW1v9ygDzox/9yBhjTENDQ8DyPl+eYy5fstkcOXIkYPny5cuNMca8++67Acsff/xxY4wxf/3rXwOWr1q1yhhjzObNmwOWv/DCC8YYY1wuV8Dyolt/ZYbf/4I5/4b7Apbv2LHDGGPM7373u4Dle/bsMQ+se9sUfuXbAcuPHTtmjDHmpz/9acDyzz77zBhjzLx58/zKRKT9vbnjjjv8ys877zxjjDH79+83tbW1Mf9dKpUMQJUJklfDjkMXkZuAacaYO6zn/wBMMsbcY9vnfOAzY8xZEbkLmGmMuSrAseYCcwGGDRt2SaCaWCR27NjB+m3vseL1/bRY4UteD3oWXwzA2aO1tJ3tWIPv17cPpwu8Nfiz7ndpa+5Ug+95Lr+bdzPlpUW88cYbNDY2dijv378/paWlAGzZssVvRb5BgwYxYcIEADZv3kxra8dJRcXFxYwbNw5jDJs2bfL7nYYPH86YMWPweDwMvXWp3/rruQWDyW8yOqsAAA14SURBVO87mHcXXc3rr7/u9/oxY8YwfPhwTp8+zdatW/3Kx40bR3FxMadOnaKystKvfPz48QwePJiTJ08G/AZQWlpK//79OXHihN83gB+u3skn5w3H0fMcWj6rx1P3vt/rewy5iJwevWj5pI7W+iN85aL+3Hb5iPbyl+v78MzOj2g5dRxP/efNS0UFPTl26iz5xV8gNy+fqcMM5SP9a9FXXXUVubm5vPfeexw6dMiv/JprrgFg9+7duN0dm69yc3O56iq/P1el0lKXJhaJyGRgkTFmqvV8AYAxZkmQ/R1AvTGmT6jjJnKm6JwVb7J1f337c980/lCzPNNpRuaUpa8EXI6gqMDJ1vnpl3hGzN8Q8AYg4dxqmxA2csGLtEY4ye3WEBPJlMp2XZ1YtB0YJSIjRCQfmAWs73SCwban04F3Yw02HlbeOblDZ6NvTZZN917p1/E4ZWQh7y+9Pm2SOQReOMyZ5wi5VkwqxToK5z8rP69JR5rMIfa5A0plu7CdosaYFhG5B9gIOIA/GGN2i8iDeNty1gM/EJHpQAtQj7dNPS11XnArHdkXDjva0MSQAicVU8cE/dDx3Uovkn0ToWLqGCqe24WnNbp6epvxxl5eWoRDJOKkHk3yV6o7ybi1XLoLV7WbX/xpNx83em+H5rtJtS/xFVmJu+qDer8bWjvzHCyZMT6pSd1V7Q64amQ4vmakQKtRhhLPZhdXtZufrH2bRk8b4J0fMKcL6/wolUhZs5ZLd+Fbs8aXzOHzm1T7aqfuhibmrdnJ052SOXx+8+pkKi8t4tbLhkX9Ot8s2MXl43FGMab/6W2HmLPizajP19lCVw3z1uxsT+YApovr/CiVKprQ09Av/rQ75Jo1kUjFcgGLy8dz62XDoppBW9Dr8xErS2ZMIC+KF2/dX9++bk/nhdki4ap2szLEtwJtq1eZptve4CKd2WvmsUrUcgGB2uurPqhnVeVhWo3BIeK3LHGo5hh7i5+viWjR+t00NEV3DXwLs63b4aaxuTWivoRlG2tDjs7RtnqVaTShZ6mvXNQ/7sd0VbupeHYXnraOzT52vsQKtCf18tKioAm9c+IuLy3qkISjGc4ItC+t7G5oouK5Xe3HDCTctxiH3opOZRhN6GmowJkXdQ21s1ffq+vS6zt3ykbr6W2HOtTSiwqcAcfWA5TM39C+T+da9exJQ6PqLLXztBp+8afdQRP6kBAx+c6tVCbRNvQ0tGj6xV1eydHd0ETJ/A0xdRz67sYUj6Yfn0jG0LsbmliwtqbDQmS+dvlYhfodKqaOIdhl7p3v0FEuKuNoQk9DvpUc+9o6DH3f/qNN9Fv313OB1XFYMn8DpQ/+OezKjQ+sq4n73ZgiHUIZaITO4vLxvL/0eh6ZOTGqkTCRxDTnsmF+Sd2Z5+Chb2gyV5lHm1zSVOe25GBc1W4WrK2hyRP8htRttscfN3q495md/OJPu2lo9ATsPLTf4i9WXfmCEaxt23dNFrpq2jthu2px+XjKhhemdGKWUvGiCT3D2WeVhmoPtmsznzdF+Jo57MeKhzldaCYJN0Jncfl4FpePj2gyUoEz/HK4kX54KpXutMklC5SXFnVp0a7OzRxdHdsx8Nz8mNufo1mzZnH5+KA3BfE51eTRCUKq29CEnkXCJbdQ7M0coWrXeTmha71TRhZS+cDXApYVhal5FxU4o16yYOWdk3lk5sT2Y3f+MDLQfuPwSPoPlMpkupZLlum8dHCkOi/Na2+ndogwe9LQLo/6CNTeH+91ZyIdtx5oiKRSmaBL66Eniib05Ch98M9hhx/m5QjLbv5iUpJboleG9I1pj0aO4De7Val0FSqha6dolvv51y/uMLszoCROiEx0B2Q0y/D6tFmLcQGa1FVG0xp6N+CqdoddH6VHbg5nW9o6bCtw5rFo+sUZ1SwR7TK8nb2/9Po4RqNU/GkNvZuz14qD3S6uczIH7zorFc+GXg8l3SwuH8/Bus9i6kdQKtPpKJduJtpVGD1tJulrq3eVb+SLfaatUt2B1tC7mYqpY6K+s1Aq1lbvqs5t9ZE0xejaiirTaQ29m4ml6SRRa6snUyQ33+jK7Fal0oEmdBVSXo5EPHMz3S0uH8+BJdfz/tLrufWyYe3rnTtE4nqPUqVSRZtcuiHfDafD7ickbXx6svnWg1Eqm2gNvRuKpGkhzyH8+paJWZnMlcpWmtC7oXCLWvXOd7DspuysmSuVzbTJpZtaeedkvwlHfXvl8fOvZ9ZEIqXU5zShd2O6DrhS2UWbXJRSKktElNBFZJqI1IrIPhGZH2K/vxcRIyIB1xlQSimVOGETuog4gEeBa4FxwGwRGRdgv3OB/w1UxjtIpZRS4UVSQ78U2GeMOWCMaQZWAzcG2O+XwK+AM3GMTymlVIQiSehFwGHb8yPWtnYi8iVgqDEm5N0FRGSuiFSJSFVdXV3UwSqllAquy6NcRCQHWA7cFm5fY8zjwOPW6+pE5IMYT9sP+CjG1yaSxhWddI0L0jc2jSs62RjX8GAFkSR0NzDU9rzY2uZzLvAFYIt418YYBKwXkenGmKB3sDDG9I/g3AGJSFWwBd5TSeOKTrrGBekbm8YVne4WVyRNLtuBUSIyQkTygVnAel+hMeaUMaafMabEGFMCbANCJnOllFLxFzahG2NagHuAjcC7wDPGmN0i8qCITE90gEoppSITURu6MeZF4MVO234WZN8rux5WWI8n4Ryx0Liik65xQfrGpnFFp1vFlbKbRCullIovnfqvlFJZQhO6UkplibRN6CJys4jsFpG2UGvDBFtnxhqVU2ltX2ON0IlHXIUisklE9lr/9g2wz1dEZKft54yIlFtlT4rIQVvZxGTFZe3Xajv3etv2VF6viSLypvV+vy0iM21lcb1e4dYlEpEe1u+/z7oeJbayBdb2WhGZ2pU4YojrXhF5x7o+m0VkuK0s4HuapLhus+aU+M5/h63sW9b7vldEvpXkuH5ti2mPiDTYyhJ5vf4gIidE5G9BykVE/s2K+23xTsr0lXX9ehlj0vIHGAuMAbYAZUH2cQD7gQuAfGAXMM4qewaYZT3+PfBPcYrrn4H51uP5wK/C7F8I1AO9rOdPAjcl4HpFFBfwWZDtKbtewGhglPV4CHAMKIj39Qr192Lb53vA763Hs4A11uNx1v49gBHWcRxJjOsrtr+hf/LFFeo9TVJctwG/DfDaQuCA9W9f63HfZMXVaf/vA39I9PWyjv2/gC8BfwtSfh3wEt47QV4GVMbzeqVtDd0Y864xpjbMbgHXmRERAa4CnrP2+3egPE6h3WgdL9Lj3gS8ZIxpjNP5g4k2rnapvl7GmD3GmL3W46PACSDmiWchRLIukT3e54CrretzI7DaGHPWGHMQ2GcdLylxGWNetf0NbcM7wS/RIl3HKZCpwCZjTL0x5mNgEzAtRXHNBlbF6dwhGWNex1uBC+ZG4CnjtQ0oEJHBxOl6pW1Cj1CwdWbOBxqMdwy9fXs8DDTGHLMefwgMDLP/LPz/mB6yvm79WkR6JDmunuJdT2ebrxmINLpeInIp3lrXftvmeF2vsOsS2fexrscpvNcnktcmMi677+Ct5fkEek+TGdffW+/PcyLim1WeFtfLapoaAbxi25yo6xWJYLHH5Xql9I5FIvIy3qUCOnvAGPN8suPxCRWX/YkxxohI0HGf1ifveLyTsnwW4E1s+XjHot4PPJjEuIYbY9wicgHwiojU4E1aMYvz9foP4FvGmDZrc8zXKxuJyK1AGXCFbbPfe2qM2R/4CHH3J2CVMeasiNyF99vNVUk6dyRmAc8ZY1pt21J5vRIqpQndGPPVLh4i2DozJ/F+lcm1almd15+JOS4ROS4ig40xx6wEdCLEoW4B1hljPLZj+2qrZ0Xkj8CPkhmXMcZt/XtARLYApcB/keLrJSLnARvwfphvsx075usVQLh1iez7HBGRXKAP3r+nSF6byLgQka/i/ZC8whhz1rc9yHsajwQVNi5jzEnb0yfw9pn4Xntlp9duiUNMEcVlMwu4274hgdcrEsFij8v1yvQml4DrzBhvL8OreNuvAb4FxKvGv946XiTH9Wu7s5Kar926HAjYG56IuESkr6/JQkT6AVOAd1J9vaz3bh3etsXnOpXF83qFXJcoQLw3Aa9Y12c9MEu8o2BGAKOAv3YhlqjiEpFS4DG86ySdsG0P+J4mMa7BtqfT8S4PAt5vpddY8fUFrqHjN9WExmXFdhHeDsY3bdsSeb0isR74R2u0y2XAKavSEp/rlaje3q7+AN/A2450FjgObLS2DwFetO13HbAH7yfsA7btF+D9D7cPeBboEae4zgc2A3uBl4FCa3sZ8IRtvxK8n7o5nV7/ClCDNzE9DZyTrLiAy61z77L+/U46XC/gVsAD7LT9TEzE9Qr094K3CWe69bin9fvvs67HBbbXPmC9rha4Ns5/7+Hietn6f+C7PuvDvadJimsJsNs6/6vARbbXftu6jvuA25MZl/V8EbC00+sSfb1W4R2l5cGbv74DfBf4rlUueO8At986f5nttV2+Xjr1XymlskSmN7kopZSyaEJXSqksoQldKaWyhCZ0pZTKEprQlVIqS2hCV0qpLKEJXSmlssT/B/k40HoLpdVsAAAAAElFTkSuQmCC\n"
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "gCPyyuXWb_ym"
      },
      "source": [
        "## X-Learner"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Aq7lb2fkcDDU",
        "outputId": "28185497-7706-4bf7-edac-495a63df06ce",
        "colab": {
          "base_uri": "https://localhost:8080/"
        }
      },
      "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": 17,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stderr",
          "text": [
            "/usr/local/lib/python3.7/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.7/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"
          ]
        },
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "RandomForestRegressor(max_depth=3)"
            ]
          },
          "metadata": {},
          "execution_count": 17
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "KPpIUdNtQK5-",
        "outputId": "ff5958ba-1f59-4ac6-f12f-639af219cf17",
        "colab": {
          "base_uri": "https://localhost:8080/"
        }
      },
      "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",
        "# Issue 38ご指摘で判明：修正（2022/04/11）\n",
        "tau_0 = M1.predict(df_0[[\"x\"]]) - df_0[\"Y\"]\n",
        "tau_1 = df_1[\"Y\"] - M0.predict(df_1[[\"x\"]])\n",
        "\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": 18,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "RandomForestRegressor(max_depth=3)"
            ]
          },
          "metadata": {},
          "execution_count": 18
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Ou5_JAZQQ0QX"
      },
      "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": 19,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "MVuQlQCrfvjB",
        "outputId": "7ebb7b38-8b32-447d-8aeb-bc6d633689fe",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 283
        }
      },
      "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": 20,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "[<matplotlib.lines.Line2D at 0x7f04db7f8650>]"
            ]
          },
          "metadata": {},
          "execution_count": 20
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3de3xU9Z3/8dcnIeEiKCBaIdxZoNji4poVLdrarlaqqFTdCsXd0tZbd9VVViwu1Au1irrqrlZ/tfhQWy+A9ZIiXtBFxcqPqNFw8cZFVCBaoWKsYCC3z/4xJ9mTzCSZHGYyM+H9fDzyIHO+Z858cmZ4z3e+53vOmLsjIiK5Ly/TBYiISGoo0EVEOgkFuohIJ6FAFxHpJBToIiKdRJdMPXC/fv186NChmXp4EZGc9Prrr//F3Q9K1JaxQB86dChlZWWZengRkZxkZh+21KYhFxGRTkKBLiLSSSjQRUQ6CQW6iEgnoUAXEekk2gx0M7vHzLaZ2ZsttJuZ3WZmG81sjZn9XerLFBGRtiTTQ78PmNhK+/eAkcHPecD/2/uyRESkvdqch+7uL5nZ0FZWOQ34vceuw1tqZr3NrL+7f5yiGkUkS7k7v/nNbzjmmGMYO3YsmzZt4r777otbb9q0aYwePZp33nmHBQsWxLX/5Cc/YejQoaxevZpHH300rv1nP/sZ/fv359VXX2XJkiVx7f/2b//GgQceyMsvv8yzzz4b13755ZfTs2dPli1bxvLly+Pa58yZQ2FhIU8//TQrV65s0mZmXHPNNQCUlJTwxhtvNGnv2rUrs2fPBuDhhx/mzTebDmbsv//+XHbZZXGPmRbu3uYPMBR4s4W2JcAxodvLgOIW1j0PKAPKBg8e7CKS2yoqKhzwSy+91N3dly1b5mYW97NkyRJ3dy8pKUnYvnz5cnd3v//++xO2v/HGG+7ufueddyZs37Bhg7u733jjjQnb//znP7u7+5VXXpmwfefOne7ufumll8a15eXlNf6955xzTlz7AQcc0Nh+1llnxbUPHDgwpfscKPMWsto8iS+4CHroS9z96wnalgDz3P3l4PYy4Ofu3uppoMXFxa4zRUVy2wcffMCwYcO49957mT59eqbL2SeY2evuXpyoLRWzXCqAQaHbA4NlItLJ1dbWApCfn5/hSgRSE+iLgX8OZrscBXzuGj8X2SfU1dUB0KVLxi4LJSFtPgtmtgA4DuhnZluBq4ACAHf/DfAUcBKwEfgS+HG6ihWR7KIeenZJZpbL1DbaHfjXlFUkIjlj9OjRfPzxxxxwwAGZLkXI4OVzRST3denShUMOOSTTZUhAp/6LSGRbt27lqquuYv369ZkuRVAPXUT2wtatW5k7dy5HH300o0aNanG9E255kQ3bdjXeHnnwfjw347gOqHDfoh66iETW1kHRkvIKhs56skmYA2zYtosTbnkx3eXtcxToIhJZQ6AnmrZYUl7BpYtWtXjfDdt2Madkbdpq2xdpyEVEImttHvpNS9fR1nnoD5Ru5sHSzTiQb0ZdG2euN6xT1Ls7M08czeTDiyJW3jkp0EUkstaGXD6qrEpqGw0R3laYh9epqKzikkWruKSVTwAAE0b05cFzj06qjs5AQy4iEtnxxx/PF198wZFHHhnXNqB39wxU1NSK93Ywbf7KtlfsJNRDF5HI8vPz6dmzZ+PtkvIK/uOxNXxZU5/Bqppa8d4Oxl3zLFef+rVOP0SjHrqIRLZmzRpmzJjBRx99xLT5K7lk0aqsCvMGlVU1XLJoVac/CKseuohEtn79em699VZ2DjmGFR93bXE9A7oV5FGV4bB/oHQzxUP6dtqeunroIhJZw0HRJ9/c1up6Drzzy+/xwbyT+a+zxpFnHVBcC25aui5zD55m6qGLSGQN0xax5PuGDb3j2Y+vZVd1XePyHgV5XHf6YU16z3NK1rLglS3UuZNvxlHD+/DG5s+pqqmL226ykp19k4sU6CISWUMPnby2L587p2Qt104eC8RCPZlhj2snj228T4OS8gpuWrqOjyqrGNBsPvq0+StZ8d6OVreZDbNv0kWBLiKR1dfHxsQtr+0e+gOlmwHiArq9WnszePDco9sM9Zknjt6rx89mSX2naDroO0VFOochP18CgFnbA+P5Zrx3/UnpLomS8gpmLFpFlEOwDScjtfZJAOKHg6aOH7TXb1bJaO07RRXoIhJZw/Va2pMiH8w7OW31hA2d9WSHPE4yzj5qcMrCPt1fEi0i+6jZdyzgL0/fRn11cgca85PoxXdGD5Ru7pA58Ap0EYns4/fXsXPNs1Cf3KyTqeMHpbmi7PVA6WaGznqSobOeTNvlCBToIhLZ/l2D2S3Npi0W9e7O2UcNbuyR55uldNghGRNG9O2wx2qvdF1jRrNcRCSyb47sy/1P02TaYveC/MYDiB0Z4M09eO7RjP/Vc3zyRXXGamhNW9Mro1APXUQiG33wfgAU9dkPI9Yzv/70sVlzav0rs0/g7KMGZ7qMDqMeuohEVlBQQK9evXj5iuNb/Bq6TGt+clI2zX5JNQW6iER2+eWXM+qEaXzzpuUtztfONmaQaLa2GUwbP5iHXtlMfbP23t0L2nX53TklaxtPpGpJvsWmfaZyXynQRSSykvIKrnhsbeO1VSoqq7jisdj0vGwN9WnjBycM22njBye81EAUDdto+Hq9ROocZv5hNZC6faUxdBGJ7PJ5v2ZLyU1NllXV1GX1FQ2vnTy2Q2bgXDt5LO/PO5kP5p1MUQvXj6mp95TuK/XQRSSyT95/ly/X/X84eUaT5dl+RcNU9cST1dr+SOW+Ug9dRCLbryAv4ZUWO/MVDaNobX+kcl8lFehmNtHM1pnZRjOblaB9iJktM7M1ZvaimQ1MWYUikrXGDewVd6XFhnno8n9mnjiagvz4yx4U5FlK91WbgW5m+cAdwPeAQ4GpZnZos9X+E/i9ux8GzAWuT1mFIpK1hvTtTs9uhRT17p6V89CzxeTDi7jpzL+lT4+CxmW9uxdw0z/+bYfPcjkS2OjumwDMbCFwGvB2aJ1DgYZBtBeAkpRVKCJZa//992fY4IGsmPWdTJeS9ZL9Uo+9kcyQSxGwJXR7a7AsbDVwevD794FeZnZg8w2Z2XlmVmZmZdu3b49Sr4hkkRtuuIFVq1ZlugwJpOqg6GXAt8ysHPgWUAHEXX7N3X/r7sXuXnzQQQel6KFFRASSC/QKIHzNy4HBskbu/pG7n+7uhwOzg2WVKatSRLLSjTfeyHnnnZfpMiSQTKC/Bow0s2FmVghMARaHVzCzfmaN18+8ArgntWWKSDZ6/fXXeemllzJdhgTaDHR3rwUuBJYC7wAPu/tbZjbXzE4NVjsOWGdm64GvAL9KU70ikkXq6uro0kXnJ2aLpJ4Jd38KeKrZsitDvz8CPJLa0kQk29XW1mbtVRb3RTpTVEQiUw89u+iZEJHIBgwYQNeuXTNdhgQU6CIS2V133ZXpEiREQy4iIp2EAl1EIrvwwgu55JJLMl2GBDTkIiKRlZeX0727LpWbLdRDF5HIamtrNcsliyjQRSSyuro6zUPPIgp0EYlM89Czi54JEYlszJgxDBgwINNlSECBLiKRPfTQQ5kuQUI05CIi0kko0EUkskmTJvGLX/wi02VIQEMuIhLZ2rVr0bePZQ/10EUkMl0+N7so0EUkMs1Dzy4KdBGJTGeKZhcFuohEdswxxzB69OhMlyEBvbWKSGQlJSWZLkFC1EMXEekkFOgiEtmYMWO4+eabM12GBBToIhLZ+vXrqayszHQZElCgi0gk7k59fb2mLWYRBbqIRFJXVwegaYtZRIEuIpE0BLp66NlDgS4ikZ122mmMGjUq02VIQJ+VRCSSrl27ah56llEPXUSkk0gq0M1sopmtM7ONZjYrQftgM3vBzMrNbI2ZnZT6UkUkm/zlL3/h4IMP5t577810KRJoM9DNLB+4A/gecCgw1cwObbbaHOBhdz8cmALcmepCRSS71NTUsH37dvbs2ZPpUiSQTA/9SGCju29y92pgIXBas3Uc2D/4/QDgo9SVKCLZSNMWs08ygV4EbAnd3hosC7saONvMtgJPARcl2pCZnWdmZWZWtn379gjliki2qK2tBRTo2SRVB0WnAve5+0DgJOB+M4vbtrv/1t2L3b1YX1slktsaAl3z0LNHMoFeAQwK3R4YLAv7KfAwgLuvBLoB/VJRoIhkp549e3L22WczfPjwTJcigWQ+K70GjDSzYcSCfArww2brbAb+AbjPzMYQC3SNqYh0Yocccgj3339/psuQkDZ76O5eC1wILAXeITab5S0zm2tmpwar/TtwrpmtBhYA093d01W0iIjES2oM3d2fcvdR7j7C3X8VLLvS3RcHv7/t7hPc/W/dfZy7P5vOokUk81atWkX37t1ZsmRJpkuRgM4UFZFIampq2L17N2aW6VIkoEAXkUg0yyX7KNBFJBKdWJR9FOgiEol66NlHgS4ikfTv358LLriAoqLmJ45LplimZhcWFxd7WVlZRh5bRCRXmdnr7l6cqE09dBGJpL6+nrq6OnTKSfZQoItIJEuWLKFLly688cYbmS5FAgp0EYlEs1yyjwJdRCLRLJfso0AXkUjUQ88+CnQRiURfcJF9FOgiEsmYMWOYOXMmffr0yXQpEtBbq4hEcsQRR3DEEUdkugwJUQ9dRCLZvXs3lZWV1NfXZ7oUCSjQRSSSe+65hz59+qAvfM8eCnQRiUSzXLKPAl1EItE89OyjQBeRSDRtMfso0EUkEg25ZB8FuohEcuyxxzJ37lwKCgoyXYoE9NYqIpFMmDCBCRMmZLoMCVEPXUQiqaysZOvWrZkuQ0IU6CLSbiXlFXz9+//CoCFDmTDveUrKKzJdkqBAF5F2Kimv4IrH1vLFl3uwvHwqKqu44rG1CvUsoEAXkXa5aek6qmrqcK+HvNgc9KqaOm5aui7DlYkOiookoaS8gpuWrqOisqrJ8v0K8/nV98cy+fCixmVzStay4JUt1LmTb8bU8YO4dvLYji45bT5q2Af1dWB58cslYxToIm0oKa9g5iOrqamL/zLkXdV1XLJoFZcsWpXwvnXuPFC6mQdKNzNhRF8ePPfodJebdgN6d4+9sXk9lpffZLlkVlKBbmYTgf8G8oG73X1es/ZbgW8HN3sAB7t771QWKtIRTrjlRTZs25WWba94bwcn3PIiz804Li3b7wgl5RVs+2usJ95j1AQK+g1pbBt6oAI909ocQzezfOAO4HvAocBUMzs0vI67X+ru49x9HHA78Fg6ihVJp3SGeYMN23bl7MHDkvIKZv5hNTXB1XK7DTmMXoef1NheuumzDFUmDZLpoR8JbHT3TQBmthA4DXi7hfWnAlelpjyRjrNh2y68roa/vlZC/Z4vAbC8LvQ+dhoAu975E9XbNjW5T15hdw44+gcA7HzzeWo+3dK0vVsvDhh/OgBfrH6W2sqPmb78d/zLt/8GgP79+3PRRRcBcOedd8bN6x4yZAjnn38+ALfeemvcpWpHjRrF9OnTAbjhhhv4/PPPm7SPHTuWqVOnAjB37lx2797dpP2II47gjDPOAGD27Nm4Nx1W+sY3vsGkSZOorq7mnAtnUFMfaq+rpee4iRT0jR0/qPP4ISnpWMkEehEQfpVuBcYnWtHMhgDDgOdbaD8POA9g8ODB7SpUpCNU/3kjlct/FzvYZ3lYl8LGQK/aVMaut5c3WT+/xwGNgV61oZQvN77apL1L70MaA/3Ld//E7s1rAbg+WK3w4KHcXDGcHgV51C++hw1vrSacmV0HjuH69wfSvSCPqoW/5aMP32uy/YkTJzYG+l133RX3hnDmmWc2Bvrtt98eF/jTp09vDPSbb7457ssqLrroIiZNmsTjZZv5tLTZB+/6Wmo++4iDz/hFbF+YIZllzd+R41YwOxOY6O7nBLf/CRjv7hcmWPfnwEB3v6itBy4uLvaysrJoVYukwdBZT7L7wzV8svA/+MrU6+g2+LBMlxTnv84a12RGTUcZd82zVFbVtLrO2UcN7lSzebKVmb3u7sWJ2pKZh14BDArdHhgsS2QKsKB95Ylkj4J+g+h3ykwKDmz5E6QBX+lV2HFFhVzx2JoOf8yS8gqFeY5IZsjlNWCkmQ0jFuRTgB82X8nMvgr0AVamtEKRDlLUuzsVwH6HfqvV9f7m4P2SmqkydNaTqSkspKqmnqGznmRkkjWkQlsnDGXqU4PEazPQ3b3WzC4ElhKbtniPu79lZnOBMndfHKw6BVjobY3hiGSpmSeOZub9f+KvH2+ia/9R5BUmnoaX7EyYkQfvl7ZZMxu27UrqDaNHQR7XnX7YXgVu85Opws4+arDCPIu0OYaeLhpDl2w088a7+M+fX8CAn9xBwUFDWlzvg3knJ7W9jpgK2R69uxdw9alfSzqEx//qOT75ojphW58eBZRf+d1UlidJaG0MXWeKioQcPrAXAC/8/Hi+e8+Gvd5ew7DInJK1PPTK5iYzWDKhsqqm1TNb2+OqU76WgooklRToIiHV1bHeaEFBARNG9GXFezvi1pkwom+7t3vt5LFxBw1Lyiv4+aNr2FNb38K9spuGWrKPAl0kpKYmNpujsLCQB889mmnzVzYJ9VRej2Xy4UWthmJJeQWXLlqFDkpJshToIiHhHjqQ0YtpNYT95Y+spjrBhcEyKcqnFEk/BbpIyKRJkxg6dCi9e2fHteXCvfhsOcCaZ5l9o5OWKdBFQoYMGcKQIS3Pbsmk9sw7n1OylgdKN6e8hjyDW34wLuXbldTQtEWRkHfffZd169ZxyimnkJfXOb7Qq6Uv52ivPj0KuOqU5Kc8Snpo2qJIkhYtWsTVV19NXV1dpktJmbYOvkrn0Tm6ICIpUl1dTX5+fqfpncu+Ra9akZCamprGGS4iuUaBLhJSU1NDYWFmrqQosrcU6CIh1dXV6qFLzlKgi4TMmDGDJ554ItNliESiWS4iISNGjGDEiBGZLkMkEvXQRUJeeukl/vjHP2a6DJFIFOgiIXfccQezZs3KdBkikSjQRUJ0UFRymQJdJETTFiWXKdBFQnRikeQyBbpIiAJdcpmmLYqEzJ8/v1NdmEv2LQp0kRDNQZdcpiEXkZCFCxfyzDPPZLoMkUgU6CIh1157LfPnz890GSKRKNBFQnRQVHKZAl0kRIEuuUyBLhKiQJdcpkAXCamurtaZopKzkpq2aGYTgf8G8oG73X1egnV+AFwNOLDa3X+YwjpFOsSrr75K9+7dM12GSCRtBrqZ5QN3ACcAW4HXzGyxu78dWmckcAUwwd0/M7OD01WwSDoNGTIk0yWIRJbMkMuRwEZ33+Tu1cBC4LRm65wL3OHunwG4+7bUlinSMW688UZeeumlTJchEkkygV4EbAnd3hosCxsFjDKzFWZWGgzRxDGz88yszMzKtm/fHq1ikTSaNWsWzz33XKbLEIkkVQdFuwAjgeOAqcB8M+vdfCV3/627F7t78UEHHZSihxZJjbq6Otxds1wkZyUT6BXAoNDtgcGysK3AYnevcff3gfXEAl4kZ9TU1ABolovkrGQC/TVgpJkNM7NCYAqwuNk6JcR655hZP2JDMJtSWKdI2lVXVwOohy45q81Ad/da4EJgKfAO8LC7v2Vmc83s1GC1pcCnZvY28AIw090/TVfRIunQ0ENXoEuuMnfPyAMXFxd7WVlZRh5bJJH6+np27NhBjx496NGjR6bLEUnIzF539+JEbboeukggLy+Pfv36ZboMkch06r9I4NNPP2X27NmsXr0606WIRKJAFwls376d6667jrfffrvtlUWykAJdJNAwy0XTFiVXKdBFAprlIrlOgS4SUKBLrlOgiwQU6JLrNG1RJDBhwgSqqqoU6JKzFOgigby8PLp165bpMkQi05CLSGDNmjVcfPHFbNmype2VRbKQAl0ksGHDBm6//XY+++yzTJciEokCXSSgg6KS6xToIgEFuuQ6BbpIQIEuuU6BLhKor68nPz9fgS45S4EuEjjnnHOora1lwIABmS5FJBIFuohIJ6FAFwksWbKEH//4x41j6SK5RoEuEigvL+e+++7DzDJdikgkCnSRQEPPPD8/P8OViESjQBcJ1NTUUFBQoB665CwFukigpqZG31YkOU2BLhIoLCykX79+mS5DJDJz94w8cHFxsZeVlWXksUVEcpWZve7uxYna1EMXEekkFOgigdtuu40LL7ww02WIRKZvLBIJlJaWomFAyWUK9EBJeQU3LV3HR5VV9CjMZ1d1XZP2PIP6Djjc0LVLHjeccRiTDy9K/4NJE9XV1bowl+S0pALdzCYC/w3kA3e7+7xm7dOBm4CKYNGv3f3uFNaZciXlFVy9+C0qq+JP824e5tAxYQ6wp7aeSxat4pJFq/ivs8Yp2DuQpi1Krmsz0M0sH7gDOAHYCrxmZovd/e1mqy5y96wZgCwpr+CaJ97isy9z97oclyxaBaBQ7yANJxaJ5KpkeuhHAhvdfROAmS0ETgOaB3qHuffee3n33XebLDvooIO47LLLALhg9jwWPV/WpFed36sf+x9xCgB/LfsjdTt3NLl/l9796TVuIgCfv/II9VVfNGkvOHAwPcf+AwCVKxbgNbubtBcePJz9Dv0WAJ+9dD/U1zZp79p/ND1GfwP3eiqX/y7ub+o66Ov0GPH31Nfs4fMVDzUu//FyWPnNEZxwwgkcf/zxfP7551x33XVx9580aRLHHnss27Zt4+abb45rP/300xk/fjxbtmzh17/+dVz71KlTGTduHBs3bmT+/Plx7dOnT2fMmDG89dZb/P73v49rP//88xk+fDjl5eUsXLgwrv3iiy+mqKiI0tJSHn/88bj2mTNn0q9fP5YvX85TTz0V1z5nzhx69erFs88+y7Jly+Laf/nLX1JYWMgTTzzByy+/3KTNzJg3L/ah8tFHH+XVV19t0t6tWzcOn3weZR/XsOvLLkyY9zwzTxytN1LJOckEehEQ/hr0rcD4BOudYWbfBNYDl7p73Fenm9l5wHkAgwcPbn+1gcWLF/PMM880WTZq1KjGQF+46A/89cM3m7QXfuVvGgN919svUbP9/SbtXQd9vTHQd735PLWVf27S3n3E3zcG+s41z1H/ZWWT9h5fPbYx0HeWP4XX7mnS7oftocfobwDwxetPxP9RZvQY8fdQXxvXflMZVOwyjj/+eHbu3Mltt90Wd/eioiKOPfZYPvvss4Tto0ePZvz48XzyyScJ24uLixk3bhxbt25N2H7ccccxZswYNm3alLD9lFNOYfjw4axbty5h+9SpUykqKmLt2rUJ288991z69evHG2+8kbD90ksvpVevXpSWliZsv/LKKyksLORPf/oTt99+e5O2cKAvW7aMe++9t0l71x496V9/ND1OvIQeQEVlFVc8thbQpyPJLW2eWGRmZwIT3f2c4PY/AePDwytmdiCw0933mNn5wFnu/p3WtpuqE4vmlKxlwStbqMvQCVIdpXtBPtefPlYBkwYT5j1PRWVV3PKi3t1ZMavVl7FIh2vtxKJkeugVwKDQ7YH838FPANz909DNu4Eb21tkFHNK1vJA6eaOeKgOm+XSkqqausaDpRALm5knjgZonJ0zIFim0G+fjxKEeWvLRbJVMoH+GjDSzIYRC/IpwA/DK5hZf3f/OLh5KvBOSqtswYJX4kZ1UubsowZz7eSxadt+c4dd9Qx/3RM/u6YlFZVVjeHefFnZhzs6tPZcN6B394Q99AG9u2egGpHo2gx0d681swuBpcSmLd7j7m+Z2VygzN0XAxeb2alALbADmJ6ugkvKK5j9+NqEUwujGnnwfjw347iUbS+KNddMZNr8lax4b0fbK7fhgdLNFA/pq556kr791YPiPukV5FvjJyCRXJFTF+cqKa/g3/+wmrqIYx8fzDs50v06Ukl5RVzPO4psH/8Nn8jVMFRU9uEOHnplc6tDWwZM28tPTyXlFVzx2BqqaupbXCcPuEXnAUgWam0MPacCvaWDV0ndd0RfHjz36Ej37WgjrngqZQd50/V3n3DLi2zYtivl202X3t0LuPrUrwEwY9EqWo7y/5Ptb4qyb9rbg6JZI+pBqlwKc4Cp4wel7GDvivd2MG3+yqT//oaec0tvnAbk4nyiyqqadn/y0UFRyTU5FegtHbxqrqE3lqsflxuGEx4s3ZyS8Fzx3g6GznoyYVu+GVPHD+L97TuTGr/PxTCPSgdFJdfk1JBLMmPoudYbb4/OcDmDXKExdMlWnWbIpeE/V0uzXDpzmEPs728eMLEDfGupqkndrB+BA3oUKMwl5+RUoEPiUNuXNeyL1sa9pf0q9SlIclDOBbrEC7/JpWoue0fqXpDH7pr6Vs90nVOyNmXHFJKh8XPJRTk1hi7tU1JekfQUvWQ1n+WS7cNc4fnuhV3y2FPb9t4w4FaNn0uW6jRj6NI+DYHU1kk0Lenoyx+kQ6IhutbOZ2g4cUlhLrlIgd7JtXTMobWhmT49CrjqlNyd9tmWmSeOTnggOdenu4oo0PdR2TxMkm7hA8m6SqV0Jgp02SdptpR0RnmZLkBERFJDgS4i0kko0EVEOgkFuohIJ6FAFxHpJDJ2pqiZbQc+jHj3fsBfUlhOqqiu9snWuiB7a1Nd7dMZ6xri7gclashYoO8NMytr6dTXTFJd7ZOtdUH21qa62mdfq0tDLiIinYQCXUSkk8jVQP9tpgtogepqn2ytC7K3NtXVPvtUXTk5hi4iIvFytYcuIiLNKNBFRDqJrA10M/tHM3vLzOrNrMXpPWY20czWmdlGM5sVWj7MzF4Jli8ys8IU1dXXzJ4zsw3Bv30SrPNtM1sV+tltZpODtvvM7P1Q27iOqitYry702ItDyzO5v8aZ2crg+V5jZmeF2lK6v1p6vYTauwZ//8ZgfwwNtV0RLF9nZifuTR0R6pphZm8H+2eZmQ0JtSV8Tjuorulmtj30+OeE2n4UPO8bzOxHHVzXraGa1ptZZagtnfvrHjPbZmZvttBuZnZbUPcaM/u7UNve7y93z8ofYAwwGngRKG5hnXzgPWA4UAisBg4N2h4GpgS//wb4WYrquhGYFfw+C7ihjfX7AjuAHsHt+4Az07C/kqoL2NnC8oztL2AUMDL4fQDwMdA71furtddLaJ1/AX4T/D4FWBT8fmiwfvFU4qUAAARUSURBVFdgWLCd/A6s69uh19DPGupq7TntoLqmA79OcN++wKbg3z7B7306qq5m618E3JPu/RVs+5vA3wFvttB+EvA0sS/HOgp4JZX7K2t76O7+jruva2O1I4GN7r7J3auBhcBpZmbAd4BHgvV+B0xOUWmnBdtLdrtnAk+7+5cpevyWtLeuRpneX+6+3t03BL9/BGwDEp4Jt5cSvl5aqfcR4B+C/XMasNDd97j7+8DGYHsdUpe7vxB6DZUCA1P02HtVVytOBJ5z9x3u/hnwHDAxQ3VNBRak6LFb5e4vEevAteQ04PceUwr0NrP+pGh/ZW2gJ6kI2BK6vTVYdiBQ6e61zZanwlfc/ePg9z8DX2lj/SnEv5h+FXzcutXMunZwXd3MrMzMShuGgcii/WVmRxLrdb0XWpyq/dXS6yXhOsH++JzY/knmvumsK+ynxHp5DRI9px1Z1xnB8/OImQ1q533TWRfB0NQw4PnQ4nTtr2S0VHtK9ldGv7HIzP4HOCRB02x3/2NH19OgtbrCN9zdzazFeZ/BO+9YYGlo8RXEgq2Q2FzUnwNzO7CuIe5eYWbDgefNbC2x0IosxfvrfuBH7t7wrdaR91dnZGZnA8XAt0KL455Td38v8RZS7glggbvvMbPziX26+U4HPXYypgCPuHv4C2Qzub/SKqOB7u7H7+UmKoBBodsDg2WfEvso0yXoZTUs3+u6zOwTM+vv7h8HAbStlU39AHjc3WtC227ore4xs3uByzqyLnevCP7dZGYvAocDj5Lh/WVm+wNPEnszLw1tO/L+SqCl10uidbaaWRfgAGKvp2Tum866MLPjib1Jfsvd9zQsb+E5TUVAtVmXu38aunk3sWMmDfc9rtl9X0xBTUnVFTIF+NfwgjTur2S0VHtK9leuD7m8Boy02AyNQmJP3mKPHWV4gdj4NcCPgFT1+BcH20tmu3Fjd0GoNYxbTwYSHg1PR11m1qdhyMLM+gETgLczvb+C5+5xYmOLjzRrS+X+Svh6aaXeM4Hng/2zGJhisVkww4CRwKt7UUu76jKzw4G7gFPdfVtoecLntAPr6h+6eSrwTvD7UuC7QX19gO/S9JNqWusKavsqsQOMK0PL0rm/krEY+OdgtstRwOdBpyU1+ytdR3v39gf4PrFxpD3AJ8DSYPkA4KnQeicB64m9w84OLR9O7D/cRuAPQNcU1XUgsAzYAPwP0DdYXgzcHVpvKLF33bxm938eWEssmB4AenZUXcA3gsdeHfz702zYX8DZQA2wKvQzLh37K9HrhdgQzqnB792Cv39jsD+Gh+47O7jfOuB7KX69t1XX/wT/Dxr2z+K2ntMOqut64K3g8V8Avhq670+C/bgR+HFH1hXcvhqY1+x+6d5fC4jN0qohll8/BS4ALgjaDbgjqHstoRl8qdhfOvVfRKSTyPUhFxERCSjQRUQ6CQW6iEgnoUAXEekkFOgiIp2EAl1EpJNQoIuIdBL/C82D2bSx5l4HAAAAAElFTkSuQmCC\n"
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "3oYzELNxl04e"
      },
      "source": [
        "以上"
      ]
    }
  ]
}