5_3_doubly_robust_learning.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "colab": {
      "name": "5_3_doubly_robust_learning.ipynb",
      "provenance": [],
      "collapsed_sections": []
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "aoxI3DOK9vm2",
        "colab_type": "text"
      },
      "source": [
        "# 5.3 Doubly Robust Learningの実装\n",
        "\n",
        "本ファイルは、5.3節の実装です。\n",
        "\n",
        "5.2節と同じく、人事研修の効果について因果推論を実施します。"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "2XdIDbdlejUk",
        "colab_type": "text"
      },
      "source": [
        "## プログラム実行前の設定など"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "XZFKJwcu-_Oj",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# 乱数のシードを設定\n",
        "import random\n",
        "import numpy as np\n",
        "\n",
        "np.random.seed(1234)\n",
        "random.seed(1234)\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "hx1idArc_F15",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# 使用するパッケージ（ライブラリと関数）を定義\n",
        "# 標準正規分布の生成用\n",
        "from numpy.random import *\n",
        "\n",
        "# グラフの描画用\n",
        "import matplotlib.pyplot as plt\n",
        "\n",
        "# SciPy 平均0、分散1に正規化（標準化）関数\n",
        "import scipy.stats\n",
        "\n",
        "# シグモイド関数をimport\n",
        "from scipy.special import expit\n",
        "\n",
        "# その他\n",
        "import pandas as pd\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "AWqP6yeQlI_t",
        "colab_type": "text"
      },
      "source": [
        "## データの作成"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "DpnGB2KZ_L8x",
        "colab_type": "code",
        "outputId": "55ffed85-5df1-4963-9799-2476ad98ebc3",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 282
        }
      },
      "source": [
        "# データ数\n",
        "num_data = 500\n",
        "\n",
        "# 部下育成への熱心さ\n",
        "x = np.random.uniform(low=-1, high=1, size=num_data)  # -1から1の一様乱数\n",
        "\n",
        "# 上司が「上司向け：部下とのキャリア面談のポイント研修」に参加したかどうか\n",
        "e_z = randn(num_data)  # ノイズの生成\n",
        "z_prob = expit(-5.0*x+5*e_z)\n",
        "Z = np.array([])\n",
        "\n",
        "# 上司が「上司向け：部下とのキャリア面談のポイント研修」に参加したかどうか\n",
        "for i in range(num_data):\n",
        "    Z_i = np.random.choice(2, size=1, p=[1-z_prob[i], z_prob[i]])[0]\n",
        "    Z = np.append(Z, Z_i)\n",
        "\n",
        "# 介入効果の非線形性：部下育成の熱心さxの値に応じて段階的に変化\n",
        "t = np.zeros(num_data)\n",
        "for i in range(num_data):\n",
        "    if x[i] < 0:\n",
        "        t[i] = 0.5\n",
        "    elif x[i] >= 0 and x[i] < 0.5:\n",
        "        t[i] = 0.7\n",
        "    elif x[i] >= 0.5:\n",
        "        t[i] = 1.0\n",
        "\n",
        "e_y = randn(num_data)\n",
        "Y = 2.0 + t*Z + 0.3*x + 0.1*e_y \n",
        "\n",
        "# 介入効果を図で確認\n",
        "plt.scatter(x, t, label=\"treatment-effect\")\n"
      ],
      "execution_count": 3,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "<matplotlib.collections.PathCollection at 0x7f22891d3160>"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 3
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAASc0lEQVR4nO3df5BdZX3H8fc3mwQjVkPIaiUJJDgRTRsFvYO0zBT8SWA6SfzZZMoIlppqxc60yhRGxjIMjra049BKa5GhqNggUqXrFCejAuOMQ5ClCBiYQIgjJFBZQZxpoZAf3/5xT+R0c3fv3ey5e7MP79fMzt7zPM8955vn3nz27Dnn7onMRJI0+80ZdAGSpGYY6JJUCANdkgphoEtSIQx0SSrE3EFtePHixbl8+fJBbV6SZqW77rrrF5k53KlvYIG+fPlyRkdHB7V5SZqVIuJnE/V5yEWSCmGgS1IhDHRJKoSBLkmFMNAlqRBdr3KJiGuA3weeyMzf7tAfwBXAWcAzwLmZ+Z9NFyrp8HTxTffxtTseYbp/52/unGDv/slXEkACEUx7e4O2cME8Lln7W6w/aUlj6+xlD/1aYM0k/WcCK6uvTcA/Tb8sSbPBxTfdx3Vbpx/mQNcwh3aYw+wPc4Cnn93DBd+4h5vu3t3YOrsGemb+AHhqkiHrgK9k21ZgYUS8uqkCJR2+Nt/x6KBLmNX27E8u37K9sfU1cQx9CVB/VXdVbQeJiE0RMRoRo2NjYw1sWtIg7SthV3nAHnv62cbWNaMnRTPzqsxsZWZreLjjJ1clzSJDEYMuYdY7ZuGCxtbVRKDvBpbVlpdWbZIKt/Ety7oP0oTmzQkuOOOExtbXRKCPAB+MtlOAX2Xm4w2sV9Jh7rL1qzn7lGNpYkd97pzuKzkwooRfDBYumMfl739jo1e5RLd7ikbEZuB0YDHwc+CvgHkAmfnF6rLFL9C+EuYZ4EOZ2fWvbrVarfSPc0nS1ETEXZnZ6tTX9Tr0zNzYpT+Bjx1ibZKkhvhJUUkqhIEuSYUw0CWpEAa6JBXCQJekQhjoklQIA12SCmGgS1IhDHRJKoSBLkmFMNAlqRAGuiQVwkCXpEIY6JJUCANdkgphoEtSIQx0SSqEgS5JhTDQJakQBrokFcJAl6RCGOiSVAgDXZIKYaBLUiEMdEkqhIEuSYUw0CWpEAa6JBXCQJekQhjoklQIA12SCtFToEfEmojYHhE7IuLCDv3HRcT3I+LeiLgtIpY2X6okaTJdAz0ihoArgTOBVcDGiFg1btjfAl/JzDcAlwKfbbpQSdLketlDPxnYkZk7M/N54Hpg3bgxq4Bbqse3duiXJPVZL4G+BHi0tryraqu7B3hP9fjdwG9ExNHjVxQRmyJiNCJGx8bGDqVeSdIEmjop+kngtIi4GzgN2A3sGz8oM6/KzFZmtoaHhxvatCQJYG4PY3YDy2rLS6u2X8vMx6j20CPiZcB7M/PppoqUJHXXyx76ncDKiFgREfOBDcBIfUBELI6IA+u6CLim2TIlSd10DfTM3AucD2wBHgBuyMxtEXFpRKythp0ObI+IB4FXAZ/pU72SpAlEZg5kw61WK0dHRweybUmarSLirsxsderzk6KSVAgDXZIKYaBLUiEMdEkqhIEuSYUw0CWpEAa6JBXCQJekQhjoklQIA12SCmGgS1IhDHRJKoSBLkmFMNAlqRAGuiQVwkCXpEIY6JJUCANdkgphoEtSIQx0SSqEgS5JhTDQJakQBrokFcJAl6RCGOiSVAgDXZIKYaBLUiEMdEkqhIEuSYUw0CWpED0FekSsiYjtEbEjIi7s0H9sRNwaEXdHxL0RcVbzpUqSJtM10CNiCLgSOBNYBWyMiFXjhl0M3JCZJwEbgH9sulBJ0uR62UM/GdiRmTsz83ngemDduDEJvLx6/ArgseZKlCT1opdAXwI8WlveVbXVXQKcHRG7gJuBj3daUURsiojRiBgdGxs7hHIlSRNp6qToRuDazFwKnAV8NSIOWndmXpWZrcxsDQ8PN7RpSRL0Fui7gWW15aVVW915wA0AmXk78BJgcRMFSpJ600ug3wmsjIgVETGf9knPkXFjHgHeDhARr6cd6B5TkaQZ1DXQM3MvcD6wBXiA9tUs2yLi0ohYWw37BPDhiLgH2Aycm5nZr6IlSQeb28ugzLyZ9snOetuna4/vB05ttjRJ0lT4SVFJKoSBLkmFMNAlqRAGuiQVwkCXpEIY6JJUCANdkgphoEtSIQx0SSqEgS5JhTDQJakQBrokFcJAl6RCGOiSVAgDXZIKYaBLUiEMdEkqhIEuSYUw0CWpEAa6JBXCQJekQhjoklQIA12SCmGgS1IhDHRJKoSBLkmFMNAlqRAGuiQVwkCXpEIY6JJUCANdkgoxt5dBEbEGuAIYAq7OzM+N6/888NZq8aXAKzNzYZOFSv128U33cd3WRwZdxqx06msW8bUP/86gy3jR67qHHhFDwJXAmcAqYGNErKqPycw/z8wTM/NE4B+Ab/ajWKlfDPPp+eHDT/GHX7p90GW86PVyyOVkYEdm7szM54HrgXWTjN8IbG6iOGmmbL7j0UGXMOv98OGnBl3Ci14vgb4EqL/bd1VtB4mI44AVwC0T9G+KiNGIGB0bG5tqrVLf7MscdAnStDV9UnQDcGNm7uvUmZlXZWYrM1vDw8MNb1o6dEMRgy5BmrZeAn03sKy2vLRq62QDHm7RLLTxLcu6D9KkTn3NokGX8KLXS6DfCayMiBURMZ92aI+MHxQRrwOOAjwzolnnsvWrOfuUYwddxqzlVS6Hh66XLWbm3og4H9hC+7LFazJzW0RcCoxm5oFw3wBcn+nBSM1Ol61fzWXrVw+6DOmQ9XQdembeDNw8ru3T45Yvaa4sSdJU+UlRSSqEgS5JhTDQJakQBrokFcJAl6RCGOiSVAgDXZIKYaBLUiEMdEkqhIEuSYUw0CWpEAa6JBXCQJekQhjoklQIA12SCmGgS1IhDHRJKoSBLkmFMNAlqRAGuiQVwkCXpEIY6JJUCANdkgphoEtSIQx0SSqEgS5JhTDQJakQBrokFcJAl6RCGOiSVIieAj0i1kTE9ojYEREXTjDmAxFxf0Rsi4h/bbZMSVI3c7sNiIgh4ErgncAu4M6IGMnM+2tjVgIXAadm5i8j4pX9KliS1Fkve+gnAzsyc2dmPg9cD6wbN+bDwJWZ+UuAzHyi2TIlSd30EuhLgEdry7uqtrrXAq+NiB9GxNaIWNNpRRGxKSJGI2J0bGzs0CqWJHXU1EnRucBK4HRgI/CliFg4flBmXpWZrcxsDQ8PN7RpSRL0Fui7gWW15aVVW90uYCQz92TmT4EHaQe8JGmG9BLodwIrI2JFRMwHNgAj48bcRHvvnIhYTPsQzM4G65QkddE10DNzL3A+sAV4ALghM7dFxKURsbYatgV4MiLuB24FLsjMJ/tVtCTpYJGZA9lwq9XK0dHRgWxbkmariLgrM1ud+vykqCQVwkCXpEIY6JJUCANdkgphoEtSIQx0SSqEgS5JhTDQJakQBrokFcJAl6RCGOiSVAgDXZIKYaBLUiEMdEkqhIEuSYUw0CWpEAa6JBXCQJekQhjoklQIA12SCmGgS1IhDHRJKoSBLkmFMNAlqRAGuiQVwkCXpEIY6JJUCANdkgphoEtSIQx0SSqEgS5JhZjby6CIWANcAQwBV2fm58b1nwtcDuyumr6QmVc3WCcAN929m8u3bGf30882vepJHTl/iMzkmT37Z3S707Vg3hze++alfP1Hj9BE6XOArL6atmThApYfvYDbdz7F/n5sADhi7hye29t9Io6cP8Rn3r2a9Sct6U8hUp90DfSIGAKuBN4J7ALujIiRzLx/3NCvZ+b5fagRaIf5Rd+8j2f37OvXJib0P8/P/Dab8Oye/Vy39ZHG1tfPH2e7n3627z+oewlzaL/en/jGPQCGumaVXg65nAzsyMydmfk8cD2wrr9lHezyLdsHEuZ6cdq3P7l8y/ZBlyFNSS+BvgR4tLa8q2ob770RcW9E3BgRyzqtKCI2RcRoRIyOjY1NqdDHZvgwi+R7TrNNUydFvw0sz8w3AN8FvtxpUGZelZmtzGwNDw9PaQPHLFww/SqlKfA9p9mml0DfDdT3uJfywslPADLzycx8rlq8GnhzM+W94IIzTmDBvKGmVyt1NDQnuOCMEwZdhjQlvQT6ncDKiFgREfOBDcBIfUBEvLq2uBZ4oLkS29aftITPvmc1Swaw13Tk/CFeOm/2XeG5YN4czj7lWJoqfQ4QzazqIEsWLuDU1yxiTr82QPsql14cOX+Iv3v/Gz0hqlmn61Uumbk3Is4HttC+bPGazNwWEZcCo5k5AvxZRKwF9gJPAef2o9j1Jy3xP9khuGz96kGXIGkGRGafLvrtotVq5ejo6EC2LUmzVUTclZmtTn2z7ziCJKkjA12SCmGgS1IhDHRJKsTATopGxBjws0N8+mLgFw2W0xTrmhrrmhrrmrrDtbbp1HVcZnb8ZObAAn06ImJ0orO8g2RdU2NdU2NdU3e41tavujzkIkmFMNAlqRCzNdCvGnQBE7CuqbGuqbGuqTtca+tLXbPyGLok6WCzdQ9dkjSOgS5JhTgsAz0i3h8R2yJif0RMeGlPRKyJiO0RsSMiLqy1r4iIO6r2r1d/9rep2hZFxHcj4qHq+1Edxrw1In5c+/rfiFhf9V0bET+t9Z04U3VV4/bVtj1Sa+/LnPU4XydGxO3Va35vRPxBra/R+ZroPVPrP6L69++o5mN5re+iqn17RJwxnToOoa6/iIj7q/n5fkQcV+vr+JrOUF3nRsRYbft/XOs7p3rdH4qIc2a4rs/XanowIp6u9fVzvq6JiCci4icT9EdE/H1V970R8aZa3/TnKzMPuy/g9cAJwG1Aa4IxQ8DDwPHAfOAeYFXVdwOwoXr8ReCjDdb2N8CF1eMLgb/uMn4R7T8p/NJq+VrgfX2Ys57qAv57gva+zFkvdQGvBVZWj48BHgcWNj1fk71namP+FPhi9XgD7ZufA6yqxh8BrKjWMzSDdb219h766IG6JntNZ6iuc4EvdHjuImBn9f2o6vFRM1XXuPEfp/1nv/s6X9W6fw94E/CTCfrPAr5D+9YCpwB3NDlfh+UeemY+kJnd7tDb8ebVERHA24Abq3FfBtY3WN46XrjFXi/rfh/wncx8psEaOplqXb/W5znrWldmPpiZD1WPHwOeAKZ2j8Le9HLD83q9NwJvr+ZnHXB9Zj6XmT8FdlTrm5G6MvPW2ntoK+07h/XbdG4Qfwbw3cx8KjN/SfvWlGsGVNdGYHND255UZv6A9g7cRNYBX8m2rcDCaN8gqJH5OiwDvUcT3bz6aODpzNw7rr0pr8rMx6vH/wW8qsv4DRz8ZvpM9evW5yPiiBmu6yXRvlH31gOHgejvnE1pviLiZNp7XQ/Xmpuar15ueP7rMdV8/Ir2/PR6s/R+1VV3Hu29vAM6vaYzWVenG8QfFvNVHZpaAdxSa+7XfPViotobma+udyzql4j4HvCbHbo+lZn/PtP11E1WW30hMzMiJrzus/rJu5r23Z4OuIh2sM2nfS3qXwKXzmBdx2Xm7og4HrglIu6jHVqHrOH5+ipwTmbur5oPeb5KFBFnAy3gtFrzQa9pZj7ceQ2N+zawOTOfi4g/of3bzdtmaNu92ADcmJn7am2DnK++GligZ+Y7prmKiW5e/STtX2PmVntYB93Uejq1RcTPI+LVmfl4FUBPTLKqDwDfysw9tXUf2Ft9LiL+BfjkTNaVmbur7zsj4jbgJODfmMacNVFXRLwc+A/aP9C31tZ9yPPVQdcbntfG7IqIucAraL+nenluP+siIt5B+4fkafnCTdknek2bCKiebhBfW7ya9jmTA889fdxzb2ugpp7qqtkAfKze0Mf56sVEtTcyX7P5kEvHm1dn+wzDrbSPXQOcAzS5xz9SrbOXdR907K4KtQPHrdcDHc+G96OuiDjqwCGLiFgMnArc3+c566Wu+cC3aB9bvHFcX5Pz1fWG5+PqfR9wSzU/I8CGaF8FswJYCfxoGrVMqa6IOAn4Z2BtZj5Ra+/4ms5gXRPdIH4L8K6qvqOAd/H/f1Pta11Vba+jfYLx9lpbP+erFyPAB6urXU4BflXttDQzX/062zudL+DdtI8hPQf8HNhStR8D3FwbdxbwIO2frp+qtR9P+z/bDuAbwBEN1nY08H3gIeB7wKKqvQVcXRu3nPZP3Tnjnn8LcB/tYLoOeNlM1QX8brXte6rv5/V7znqs62xgD/Dj2teJ/ZivTu8Z2odw1laPX1L9+3dU83F87bmfqp63HTiz4fd8t7q+V/1fODA/I91e0xmq67PAtmr7twKvqz33j6p53AF8aCbrqpYvAT437nn9nq/NtK/S2kM7w84DPgJ8pOoP4Mqq7vuoXcXXxHz50X9JKsRsPuQiSaox0CWpEAa6JBXCQJekQhjoklQIA12SCmGgS1Ih/g9Sp439On0bGgAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "BHcdUlW9koTa",
        "colab_type": "text"
      },
      "source": [
        "## データをまとめた表を作成し、可視化する"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "1EMwdGIIIPrK",
        "colab_type": "code",
        "outputId": "f142dfd1-6479-4131-bd79-d3aee325e835",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 196
        }
      },
      "source": [
        "df = pd.DataFrame({'x': x,\n",
        "                   'Z': Z,\n",
        "                   't': t,\n",
        "                   'Y': Y,\n",
        "                   })\n",
        "\n",
        "df.head()  # 先頭を表示\n"
      ],
      "execution_count": 4,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/html": [
              "<div>\n",
              "<style scoped>\n",
              "    .dataframe tbody tr th:only-of-type {\n",
              "        vertical-align: middle;\n",
              "    }\n",
              "\n",
              "    .dataframe tbody tr th {\n",
              "        vertical-align: top;\n",
              "    }\n",
              "\n",
              "    .dataframe thead th {\n",
              "        text-align: right;\n",
              "    }\n",
              "</style>\n",
              "<table border=\"1\" class=\"dataframe\">\n",
              "  <thead>\n",
              "    <tr style=\"text-align: right;\">\n",
              "      <th></th>\n",
              "      <th>x</th>\n",
              "      <th>Z</th>\n",
              "      <th>t</th>\n",
              "      <th>Y</th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>0</th>\n",
              "      <td>-0.616961</td>\n",
              "      <td>1.0</td>\n",
              "      <td>0.5</td>\n",
              "      <td>2.303183</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>1</th>\n",
              "      <td>0.244218</td>\n",
              "      <td>1.0</td>\n",
              "      <td>0.7</td>\n",
              "      <td>2.668873</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>2</th>\n",
              "      <td>-0.124545</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.5</td>\n",
              "      <td>2.193123</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>3</th>\n",
              "      <td>0.570717</td>\n",
              "      <td>0.0</td>\n",
              "      <td>1.0</td>\n",
              "      <td>2.245229</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>4</th>\n",
              "      <td>0.559952</td>\n",
              "      <td>0.0</td>\n",
              "      <td>1.0</td>\n",
              "      <td>2.139868</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "</div>"
            ],
            "text/plain": [
              "          x    Z    t         Y\n",
              "0 -0.616961  1.0  0.5  2.303183\n",
              "1  0.244218  1.0  0.7  2.668873\n",
              "2 -0.124545  0.0  0.5  2.193123\n",
              "3  0.570717  0.0  1.0  2.245229\n",
              "4  0.559952  0.0  1.0  2.139868"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 4
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "L6Tb2Hjk9vno",
        "colab_type": "code",
        "outputId": "7650ce67-a6b7-42e1-9388-5335c313ae7f",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 282
        }
      },
      "source": [
        "plt.scatter(x, Y)\n"
      ],
      "execution_count": 5,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "<matplotlib.collections.PathCollection at 0x7f22864681d0>"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 5
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO2df5Qc1XXnv3d6WlKPyGqEURwYEJK9jghYQRO0QFY5a0ubSNgEmAAOkODgxD6svU42EKITYbNGUsiiRBuz3nV2bR3bazsmWBjBRBiIkCN5HbMWtsYjWR6QbGzMj4Y1Y9DISBqknpm7f1TV6HX1e69e/equ7rmfc+bMTNXrqlevq+57dX8SM0MQBEHoXLpa3QFBEAQhX0TQC4IgdDgi6AVBEDocEfSCIAgdjgh6QRCEDqe71R3QccYZZ/CiRYta3Q1BEIS2YWho6GfMvEC3r5CCftGiRdi7d2+ruyEIgtA2ENFzpn2iuhEEQehwRNALgiB0OCLoBUEQOhwR9IIgCB2OCHpBEIQOp5BeN4IgCO3G4HAVm3ccwktj4zirt4K1a5ZgoL+v1d0CIIJeEAQhNYPDVdz+4AGM1yYBANWxcdz+4AEAKISwj1TdENEcIvo2Ee0nohEi2qBp8z4iGiWiff7PB5R9NxHRD/2fm7K+AEEQhFazecehaSEfMF6bxOYdh1rUo3pcVvQnAKxi5qNEVAbwTSJ6jJn3hNptZeY/VjcQ0ekA7gSwHAADGCKi7cx8OIvOC4IgFIGXxsZjbW82kSt69jjq/1v2f1yrlawBsJOZX/OF+04AlyXqqSAIQkE5q7cSa3uzcfK6IaISEe0D8Ao8wf2kptk1RPQ9InqAiM7xt/UBeEFp86K/TXeOm4loLxHtHR0djXEJgiAIrWXtmiWolEt12yrlEtauWdKiHtXjJOiZeZKZlwE4G8DFRPT2UJOHASxi5l+Ft2r/QtyOMPMWZl7OzMsXLNDm5REEQSgkA/19uPvqpejrrYAA9PVWcPfVSwthiAViet0w8xgR7Yanfvm+sv1VpdlnAPyN/3cVwDuVfWcD+HqSjgqCIBSZgf6+wgj2MC5eNwuIqNf/uwLgtwAcDLU5U/n3SgBP+3/vALCaiOYT0XwAq/1tgiAIQpNwWdGfCeALRFSCNzHcz8xfJaKNAPYy83YA/4mIrgQwAeA1AO8DAGZ+jYj+EsB3/GNtZObXsr4IQRAEwQwxuzrQNI/ly5ez5KMXBEFwh4iGmHm5bp/kuhEEQehwRNALgiB0OJLrRhCEGU2Rk5FlhQh6QRBmLEVPRpYVoroRBGHGUvRkZFkhgl4QhBlL0ZORZYUIekEQZixFT0aWFSLoBUGYsRQ9GVlWiDFWEIQZS2BwFa8bQRCEDsaUjCxPt8tmu3SKoBcEQQiRp9tlK1w6RUcvCIIQIk+3y1a4dIqgFwRBCJGn22UrXDpF0AuCIITI0+2yFS6dIugFQRBCZO12OThcxYpNu7B43SM4dmIC5RJldmwXxBgrCIIQIku3y7DxdWy8hnIXYX5PGWPHa+J1IwiCEJesXBezqgGrM77Wphg9s7ox/LHVqY/vQqSgJ6I5AL4BYLbf/gFmvjPU5s8AfABeKcFRAH/EzM/5+yYBHPCbPs/MV2bXfUEQhFOkdV3Mw7+9CPl0XFb0JwCsYuajRFQG8E0ieoyZ9yhthgEsZ+bjRPQhAH8D4Dp/3zgzL8u224IgCI3YXBejBHZek8RZvRVUNUI9ML42I3gqUtCzV1T2qP9v2f/hUJvdyr97ANyYVQcFQeg88hJuplVydWwcKzbtsp4vr0li7ZoldfsAz/i68rwF6N/4OA4fr9X1M4/gKSevGyIqEdE+AK8A2MnMT1qavx/AY8r/c4hoLxHtIaIByzlu9tvtHR0ddeq8IAjtRyAUq2PjYJwSboPD1dTHNrkokn8e2/nSqFiiJom7r16Kvt4KCEBfbwXXXNSHbUPVOiEf/lyWOAl6Zp701S9nA7iYiN6ua0dENwJYDmCzsvlcvzL57wH4b0T0VsM5tjDzcmZevmDBglgXIQhC+5BnZKjOLZIQUkEYzpfGvz1qkhjo78MT61bh2U2X44l1q7D74GjDGLgcLymx/OiZeQzAbgCXhfcR0W8C+CiAK5n5hPKZqv/7xwC+DqA/RX8FQWhz8jRO6lbPYSFvOl8a3/m4k0TUtWYdPBUp6IloARH1+n9XAPwWgIOhNv0APg1PyL+ibJ9PRLP9v88AsALAU9l1XxCEdiPvyNDw6rnP8Xy6SeLuq5c66crjThK2a80jeMrF6+ZMAF8gohK8ieF+Zv4qEW0EsJeZt8NT1ZwG4CtEBJxyo/wVAJ8moin/s5uYWQS9IMxgTMbJvCJD45wvqe983AArXZ8AoLdSxvorL8jc64Y8p5pisXz5ct67d2+ruyEIQk40Ox972vPl0d+sj0lEQ749tHGfCHpBEAQzYddJwHsjcFXrNAuboJcUCIIgzDjirKbT+NcXBRH0giDMKOJGwBYhhUFaJE2xIAgzirh+/K3IH581IugFQWgL1JzuKzbtShxJq8s7Y9uedW76ViCqG0EQCk+WBbVLRJjUOKGUiDSts81N3ypE0AuCUDjCxtLjJycyM4jqhLxtO5A+N32z3UnDiKAXBKFQ6FbvJpIYRPsMaYNNEbRpyfJtJCmioxcEoVDojKUmkhhEm61zzzOJmyuyohcEoVC4rtKTCudm69yL4J4pgl4QhEJhqsjUWylj7uzuTIRzVvVgXYiqMNUMRNALgpA7cYyRpiRkrsm+Wm34DNPsJG46RNALgpArcY2RaVQrRTB8himCe6YkNRMEIVdWbNpl9HJ5Yt2qQp2raG8DcZCkZoIgZEocgdhMY2SacxXxbSArxL1SEBzJKgS/3Ylb3LuZuWLSnKsIbpB54VJKcA4RfZuI9hPRCBFt0LSZTURbiegZInqSiBYp+273tx8iojXZdl8QmkNc4VYU8pic4grEZvqtpzlXEdwg88JFdXMCwCpmPkpEZQDfJKLHmHmP0ub9AA4z878mousB/DWA64jofADXA7gAwFkAvkZEv8zMbtEQglAQ2jEneV6qiLgCMQ9jpEl1lOZcRXCDzItIQc+etfao/2/Z/wlbcK8CsN7/+wEAnySveOxVAL7MzCcAPEtEzwC4GMC30nddEJpHO6728pqckgjELP3WoyawpOcqghtkXjjp6ImoRET7ALwCYCczPxlq0gfgBQBg5gkARwC8Sd3u86K/TXeOm4loLxHtHR0djXcVgpAz7ZiT3DQJVcfGU6lw0qhHslAl5aVLH+jvw91XL0VfbwUEz1OnaOUCk+LkdeOrWpYRUS+Ah4jo7cz8/Sw7wsxbAGwBPPfKLI8tCGlpx9WeaeUNILYKJ6wqueaiPuw+OBpLPZKVKinPt6tmRsw2k1julcw8RkS7AVwGQBX0VQDnAHiRiLoBzAPwqrI94Gx/myC0FUUIeomLbnIKcFHhBMK9OjYOwil9bXVsHNuGqrFXu1mpkjpZl54XkYKeiBYAqPlCvgLgt+AZW1W2A7gJnu79WgC7mJmJaDuAfyCij8Mzxr4NwLezvABBaBbtttoL+nrL1n3a/bYVcHj1HX7FTiKgs1qJt+PbVatxWdGfCeALRFSCp9O/n5m/SkQbAexl5u0APgvg731j62vwPG3AzCNEdD+ApwBMAPiweNwIQja4BC0N9PdNr8rD2FbALqmCX/J1/a5vOVmtxNvx7arVSAoEQWhDwituwFvVhtUpg8NVbHh4BIeP1+o+r2ursnjdIw2r+DC9lTJOTExF9iFun4VkSAoEQSgYppWw6wrZRd+tE6yAJ6CjMkHaDLmAJ6CJEEvnLivx1iGCXhCazOBwFWu/sh+1KW/NXB0bx9qv7Mfe517DtqGqk1eKi77bpH6ZO7s7Urjq9OCBQbbPF9C3JtD9t5udo1MQQS8ITWb99pFpIR9Qm2Lcu+d5rdHztvv349at++pWwKYV97xKGSs27cJLfqoGHS7GT5fVdxLdv9AaRNALQpMZG69pt5sE8ySfWvkHK3zdirvcRTh2csJ4/ABXQRy1+hbvl/ZBslcKQhuh6sDDUZynzelGbdJuQs1SEHdyJGmnISt6QWgy83vKDV4wADC7uwtdRE5ujUCjesUm4gnIxfgpOvf2QAS9IDSZO6+4ALd9ZT8mQ3r6qSnGey4+ezq1QBfRtNpGJVC9mLxqwhCAe65blmlSMfGcaS9E0AtCkxno78P67SMNuvTaFGP3wdHpkncmv/NA9eIS1AR4uv8g4VdaAd3JVZg6GdHRC0ILOGIwmKoeMQP9fbjmoj6UiAAAJSJcc9EpVUmc1AGBQLYVTnHJLNnJVZg6GRH0gtACXNIeDw5XsW2oOq2+mWTGtqHqtACO48ZY0uj+VQHtWkGrHfPyCyLoBSE3bCtkXU53ArDyvAXT/9tWz4PDVRw/OeHUj0q5pNX1A6cEtOlc67eP1G1rx7z8ggh6QXAmTtGMqBVyoJYh5TMM1K3YbYVDbn/wgNZzR0V1eeyLENCmc42N1yInKAA4dmKi8PVzZzIi6AXBgbjFwV102bsPjhrT/wLmVbJODROmr7eCZzddjifWrcJAf19kVaioTJYBge/8/J5yXZux8VpbFEufqYigFwQH4hohXXTZUW1MwtmkhlHbhIOiooKbbEFU4X4O9PdB1wUxyhYXca8UBAfiGiFdcq9HtTHlmzHlmAFOJRwzZY80uUAO9Pdp0xkDXv4clcHhqjHNghhli4ms6AXBgbhGSJcC2i5tBvr78MS6VbjnumUAgFu37sOxExMol6jucwTgxksXTqtqknDnFReg3EUN24+drNe/21btYpQtJpGCnojOIaLdRPQUEY0Q0Z9q2qwlon3+z/eJaJKITvf3/YSIDvj7pJqI0Ja4CGUVlzwwrrliwvaBsfFaQ06bsCE3CQP9fThtTuNLfm2S64S7bdUuCc2KiYvqZgLAbcz8XSL6BQBDRLSTmZ8KGjDzZgCbAYCIrgBwKzO/phxjJTP/LMuOC0IzSVI0wyUPjEsb1wjYJHVcw4wZPHlU4W5SOc3vKUt0bEGJFPTM/DKAl/2/XyeipwH0wasDq+MGAPdl1kNhRtAO+VNcE3hlfS1x9N4ubW39c7EtmNIT33nFBc79FJpLLB09ES0C0A/gScP+HgCXAdimbGYAjxPREBHdbDn2zUS0l4j2jo6OxumW0ObEdV0sMnlcSxy9d1RbXf9u2boPyzY87lW+crQbSHri9sK5ODgRnQbg/wD4K2Z+0NDmOgA3MvMVyrY+Zq4S0S8C2AngT5j5G7ZzSXHwmcWKTbu0q8i+3sp0gq+iEl4dHzuhL/yR5lpcs1QC0fVgTWMNnCrUDUhd13YkdXFwIirDW6XfaxLyPtcjpLZh5qr/+xUiegjAxQCsgl6YWbRD/hSdugNAQyZHE2muJRCyJvdHlSBwSf2ciq2PgY4/jeeOUEwiBT0REYDPAniamT9uaTcPwDsA3Khsmwugy9ftzwWwGsDG1L0WOgoXvXAr0aXmvXXrPmuhjzBhX/S4TAc2KUXFTZjqzA4OV6cLfJsIJqR2sJkI7ris6FcAeC+AA0QUlH3/CICFAMDMn/K3/Q6Ax5n5mPLZNwN4yJsr0A3gH5j5n7LouNA5FL32qM7rJY6QBwBqdE/XYhOwm3ccihTyAbo6s5t3HIrsNwNYtuFxHDs5Me3CKTnn2x8Xr5tvAoi8TZn58wA+H9r2YwAXJuybMENI4rrYTLJQIY0dr0WukqOKeiTtR6CScf28zsaQheum0DokBYJQCIpce9SkWorDvEo5sjKTLZ/OQH9fqn4Ek0ua6yiSzUSIh6RAEIQITKl5XamUSyBCZFK0JEnOXAneIHQ58OMcQ2hPRNALQgSq3zgQTzgG5f9cI051qEnObLnlbQRFSnQ58F2up0g2EyE+zn70zUT86IUiE9a1Hz85Een22EWAzo6q+tebioHrgpEWrXvEej6dd02lXMLs7i6tDj7cvtxFOG1ON8aO1wpnMxH0pPajF4SZjs2Q6hLQpBPywSpZPXZvTxmzu7twZNwuYEtE1rz0uj3jtUljHxnepFNEY7iQHhH0ghBBlDeMakyNMnaWiDDFbAy6Ony8hkq5hHuuW2YVtDdccg6+tOf51NcW0A5RyEJyREcvCBG4VJcK8sZH6c+nmOtK/MWtXBVw18BS3HjpQq1+vVIuNZT6C5jfU46VblnoDGRFLwgKqhplXqUMIhj17zovGV3wl0rY4Jom/cNdA0tx18BSp/QMQH2GyaLGLAj5IIJeEHzCKhpTubwAnZdMIDDXbx9p+Lxu5ZxF+gdbDIJJoItgn1mIoBcEH9cCH0B0dakgt4xupb1i067pbSvPW4BtQ9Vc0j8UOQhNaC4i6AXBxzXy01aAWyUsaHVG3W1DVVxzUR92HxwVVYqQGyLohY4hbcZFlxQBabxTTIbX3QdHxeNFyBXxuhE6giwqO0WlGAirVAaHq1ixaRcWr3sEKzbtijxXO+TdFzoTWdELuZJVXvOo40QlBHMhnEUz8LrRRYdG+dbrKHrefaFzkRQIbUgS4ZlU4KYR1LqI0SDU3lXPbTpOODXA4nWPaKNBCcCzmy536m8ckpQ/jJPiQBDiYkuBIKqbNiOJiiKpWiOtOsRWsCPOsVyCiqISgmVNEjWMFNUWWoWobtqMJCqKpGoN0+d0Zep0ROmeXVUrLkJVF6hULhGOnZjA4nWPRPY17ptLUjWMuDwKrSByRU9E5xDRbiJ6iohGiOhPNW3eSURHiGif//MxZd9lRHSIiJ4honVZX8BMI8lK0rSvOjZuNSaaPjfJ7LTCd1lNuxgiXVbr4dXy/J4ywF7QU9DXW7buQ//Gxxv6q3tzMbUN0BluJZWAUFRcVDcTAG5j5vMBXArgw0R0vqbdvzDzMv9nIwAQUQnA3wF4F4DzAdxg+KzgSBIVhWkfAVa1TK8hX4qKLS+LS6GMs3orkd4rrkI1yDfz7KbL0TOrW1tf9fDxWsN1bnh4RBsopWurnkvUMEK74FIz9mUAL/t/v05ETwPoA/CUw/EvBvCMXzsWRPRlAFc5flbQkKSQtu4zunzlqiplcLiKo29MOPXJtCoPZ3UMn7NSLmHleQsivVdMNWWB+ihTVd1i84cPX6ctl7xNVSVqGKFdiKWjJ6JFAPoBPKnZ/etEtB/ASwD+nJlH4E0ILyhtXgRwieHYNwO4GQAWLlwYp1uFI5xfnBmR+cVdSVJIW/cZkyAMhPbmHYe0K2IdtrcJVRjq9OCu9gOXKNNgggD0E5npOqMI8r67uFAKQhFxFvREdBqAbQBuYeafh3Z/F8C5zHyUiN4NYBDA2+J0hJm3ANgCeO6VcT5bJMICSF0tZiEokro7hgWlyT0wENquQTy6ICJbIq1wX2/duk97XNP5g+Pr+q6qkaJuoLjXGT6HCHqhnXByrySiMjwhfy8zPxjez8w/Z+aj/t+PAigT0RkAqgDOUZqe7W/rWKISY7nkGjeRRfRnQJTe27RKr5S7UCIvC3pQDzUcRKT279at+3DH4AHtsWzn0W1Xj2/ipbHxyDQGLtdpQyJZhXbDxeuGAHwWwNPM/HFDm1/y24GILvaP+yqA7wB4GxEtJqJZAK4HsD2rzmdF3FB2Gy5CIKmgSFqkQkeUMVE3EZS7CBNTPK3KmGTGtqHq9HiZ/Obv3fO8cUxXnregoXgG+dvDuGSX7O0pW4tdu1xn8PlgQgsjkaxCu+GiulkB4L0ADhBR8J79EQALAYCZPwXgWgAfIqIJAOMArmcv5HaCiP4YwA4AJQCf83X3LcGlQENa9YpLYqywoHBVx2SdK8VmTAwbUktEWp29OtGYrpsB3Hb//rrjAt51bxuqNqhZGMC2oSqWn3t6Xfuo66yUS2DWq20I0Jbns9k8TJGs4kLZGWSVnqMdmDEpEEwP7Zxyl9brImmWwqhC0eGQ98HhKtY+sB+1yVPfQ7lE2HzthQ03nUmv3lspY+7s7lxuWJfC14B3XS5t1Gs3XU9A+DuwtQ9SKty6dZ9RP/+TBKkQZpIwaHfifFe2dBRAe1bgkhQIMKs94pSJc0EXuNNb8dQJvZUy5pS7cOvWfdMqog0Pj9QJeQCoTTI2PNz44mNSpxw7OZGJ3l6HazEO1zaqmilqjMP7TX75vZXy9MNoUqtE1XIVmkeWqlL1mHHsVyZ5sH77SGZ2sCIxY1IgxBXccfWwttXE4HAVH33oQF1pueAGMgnIw8drDaH7OjXD8ZMTDZOVKlDTrkyi1FBxCb6HweEquoim9f06wt9B0PcND4/UXfPYeG1a3ZYkzsBEkgyVgp28xjRumg+TPNCVjwxiKdL2sZXMGEFv0p1Xyl0AKJVgiPLpDqtmAqJWweqKAjilU1dvtsXrHtF+NvB4CScRC47jel1R/uhxCSJhb3/wgFXIA8CxExMYHK42+NNv3nHIOLkFqp4sXr2zSH0s1JPXmMa1X7nY0lQmmTOb5O8YPID7nnwBk8woEeGGS87BXQNLUx0zihkj6NeuWYK1X9nfYFCcmGJc92/OTlXKbf32xhB6dVWtE/JxCEdyqkKst6dsVD/ZIl9d2LzjUKZCvlwirDxvAW67f3+kkAfqV+ouRtmk6jbT25gUCsmevMY0bpI505ufyWYH6N+UbTULdNwxeABf2vP89P+TzNP/5ynsZ4yOfqC/D6fNaZzXapM8Xcrt2U2X44l1q2IJ+cHhqvZ1D/Bu3qgbeH5PGeUum0PgqWPp9JBH35hAuRT9efU4ebR1YXKSsfU7LzgJ+QCd+6jN9z6urtbWvtmpj2cCeY1p3CRzJvfiO6+4wJqfKbg/gvtlbLyGw8drzvr8+558Idb2rJgxgh7wZl0dgRBNYiCy+bCf1VuJvIGZgesuPmf6hrP5butee2tTjNokGz/XcD7A+friPnw9ZfvtNIVkbzcuRtngoY4ba2BrLxkqsyevMU2SZE5Nghcs8ILjmJ4mgl3lGhXXYlrkxFn8JKHjBb0qwLsMwnBepZzY0m5b9a5dswRr1yyxrrjHxmvYNlTF2jVL8Oymy/G3v3uh8UGwnWuS2RoopOJ6fSYvH51AL3cR/svVv4reSnTGy7jojLKmh9qWklk3idtUCZKhMnvyHFOd4E56HFPmVhdxrLunAjlkI0svpDAdraMPG0l1s2alXAJR4yw9XpvEhodH6nS3K89b0KDLN+kG5/eU6260sKdI+FyB7lznWTK72xOsUQYkRmMyL5Mx1UVfHxVMZDJ4uvjdu2Ja7ZmCvWxjpDNuR+l2JUPlKbKKKSjSmJquyfT270JvT7kuq+rK8xZg21A1kfNFVnR0wJQpwKZEhCnm6S/WFmRjo1Iu4ZqL+hq+RFMd0MHhKm4xJPFSa5uagjl059LR11upu3FN15dXPVXbdYYplwjlLsLx2hQAzwtqTrnkbNzSndtlogmCsaSOqx01iZwuzXQ7j5Puuy+XCHNndRvtblGUSwQw6pw+kniuJQnYtAVMdfSK3vRaPsVcJ+BM2RCjGK9NYvfBUdx99VKnlU6UPl9tp3vDuO/JF3DDJedg98FRY39LRA39MF1fUgOY+vCXfF/4OMW+A5J8JorwW4jpAQvujSRpn4tIHhG8YUGY1ouraGhtXpOcWMjP91OShz9vugfJsj/r+JWOFvSuLldr1ywx+rpHEehyXW72KH1+VLsgiVgQpq1buepyp+cZRKQ73/rt0emMCEiUYsIF9fuISsccbt+ONDMIKUyUZ1aRU0hk7VXWM6s71jGDe1B3fxLQEEOSho41xg4OV3H8ZGOFJKOAS6jBCr6sOwYP4K23P4pF6x7BW29/VJua17SCDuvzbSttteLR7O4uzPezNeq8btQVV1YGMNvDH4SQu6yImuWiOBM8Z7LMaqqSpp4vkG1a7TzI+h4MJjMXyiWadtbQOVEw3IriuNKRK3qTnra3Usb6Ky9oEHBxqimpBALDJQhicLiKYyf0E8+dV1xQt80U3KUeH/BeESvlEu65bplRJ66qKJIIdp2axoaLkG+moE2iminyKlRHs4OQAqK+xzhRsHETkiVVHarnmVcpo1yi1AGNAUG/XWxEE5M8XZ4ySr2YBR0p6E3FnufO7o6V9yKKOb6boS0I4q6BpcaJp4vqV151fXP0lQxW0SaDj2mFEb7hddF9Ll5LcTE9kHkK1ziTXDvmt4kbFeqKrdawi2B1nYB0Y37L1n1Yv32kYWHmojrU9SlIIBjOkVTuIszvKWPMD3py4c2/MAs/ff1k3bZg0nO1EampSeI+u0noOEFvK/as3mCqYDGNdJS1/PDxmjVnS7DdpO4IFuzhm3TzjkOxVhmmVTQB2hVX+GHRJVuz9duELYR8fk8Zwx9brf1ckYRrO+a3ydIGo5LWUG2agIKgPdVZQHef6VJgRKkOdalCenvKOPrGhPYNuTbFOHy8hr4YuW9+drSGGy9dOO0UUSJqWLAF/e3f+Li1+HwwHjqPpizfejtOR+/i2RLWHZoWqj2zSpHpCWyCMNCbu7wxqDdKVq9sDL2gdC136NKP4BptIeQ69VRUf7LQMSehHfPbFDWwy5RWGjg1md8xaC8Nqd4Hg8PVSGGsSxVy+HgtUjUbN8HZtqEqVp63AJVyqeGtIrBBDA5XccTRHz94S8rr+4tc0RPROQC+CODNfn+2MPMnQm1+H8BfwJuYXgfwIWbe7+/7ib9tEsCEyc8zK1w8W1xXqsdOTsbKIxPmhku8crmumfKCvsfJrBdVPMV2HhsuOnmbr2+cVWCRhGteahAdWaqrotRTSc4VlZU1KphQXYWbirnfq9i2TATCe+0D+yPbmlKFZE3g7hx+PlRnCZC7j0fSQkeuuKhuJgDcxszfJaJfADBERDuZ+SmlzbMA3sHMh4noXQC2ALhE2b+SmX+WXbfNmB7U3sopzxZXAVIiN0NNn3+jh1OPLj/39Gn3PpegiSDsOsqgUyl34Y3a1PQDtfe513DvnucbXv1WnregLkJPLc7hMpHYhHxUwqg4AquZwjWKvNQgYZqprkp6LlNW1vXbR3BiYqrueKozgi619uJ1j2jvfxdB2EWEjzz4vchnMfiebnUM1ktLZN6aGCat6th4nTorayIFPTO/DOBl/+/XiehpAH0AnlLa/L0d/cIAAB3lSURBVF/lI3sAnJ1xP50xPajrrzylOrCl9g0guBkfVSOMmmbUFmxiEvpv+G2jVkIATdc/1dVdJQC/tnBeXRRtlF+9jWBlnyY4ysbK8xbUCQp1e7NpVgBVlC0gvAK3rZjTnkuHLSuri2dV+Phx87+rTDLjeM3+LKr3ZNIAyFaT52QfyxhLRIsA9AN40tLs/QAeU/5nAI8TEQP4NDNvMRz7ZgA3A8DChQvjdKsOlwf1DQcBpzOQBIRTKMTRg9tW9uO1qekgieBn2YbHGx4sVW+py+3OAPb8+LD2tVJXnCPwujFNflPMudZb3X1wVPt50/a8aUYAlU1dpVuBm1bMQPSklEQ1loV9JDj+HYMH8PKR/ARv+N40LWSIPHtcb6WMkxOT02k3dKh1mG3G3KzJy/DvLOiJ6DQA2wDcwsw/N7RZCU/Q/4ay+TeYuUpEvwhgJxEdZOZvhD/rTwBbAC/XTYxraCDqQR23fMF1fdJsi8rvofr4mo5p032rX7JtVRU87HHTntr86l2iSF2Joy4oko6+WcyrlLXfrauOOUi690ZtKnKMk6jGshj7s3orDTEmWaPLlhpcezh4L3gkot5Iyl2kdeuMcpnMijzeRpy8boioDE/I38vMDxra/CqAzwC4iplfDbYzc9X//QqAhwBcnLbTaUgTlUcArrnIPIncMXgAt27dF/lF2VRC6gNmSyUQuHTZ9uuwPdxZRpHG8aSZaQU+BoerOKaJ2i53UWQ6apXDx2tOYxz3ew3q+aYhOH7eBTV++8IztdsH+vswd3Z87/HeShmb33NhwzM+0H8qBbJr7Yek5HH8SEFPRATgswCeZuaPG9osBPAggPcy8w+U7XN9Ay6IaC6A1QC+n0XHw7gWDknzSsrwgqB05xgcrjYYRE309VYw35DvWhVutpVHlKH0hkvOiS20s3TTi7NKzytNQdJiMnljipM4bU73tLE8DeExjvO9utbzVenrreDGSxdqj593QQ2bei/OWwnBUwHtu3N15P0eeNPlRR5j5jLlrQDwXgAHiCgwZ38EwEIAYOZPAfgYgDcB+J/evDDtRvlmAA/527oB/AMz/1OmVwC9mmDtV/Zjw8MjDdGeaV9JTZF4rvVV1SCmNN4dJvVPiWj6IVt+7umxjYpZ6afjqAvyMIAWKQgrjOmNL7CRuBjLK+USZnd3GdU/YVy/1ySuidWxcew+ONqSVBG25zmOATgqZ0/43tz67efhqAGOjcktOg0uXjffRERAPjN/AMAHNNt/DODCxL1zxFRiL3hw1Ifc9uXPnVVCb88svDQ2ji6HvC6q4cR1AgkHMdmE23yLd5CpiEog5FudryWum2LWBtAiR7jaJmlAP/HpvG4A+2IhyT2QdCHUqonUluJDl1tKhymCPDiObsFQLnWhNhVP0pe6CJMRBt288kB1RAqEOJGna9csMSYAO35yEiMbPY8U1wIWcYOc1NlaFW7BQxkkOlq7ZgnuvOIC5/TJ6kq+CKvZZrkpmmiGgTfpZOpiQA9PfIPD1QY1hW2Mk94DadwgdROpKbVAX28Fx09ORLo52zAJRddnN+D3L11oHBPTgiEJJiHv4sWXlo4Q9K43Z3Vs3Kqjn1epLwF2zUV906so0wo/WFG4vm673JjBQ3n31Uux+doL6x5k03VOMVvzgbRiNdsMN0UTprHqIsLidY+kfqjSTKY24ace35SrRReQFCaun34wFqZEZl0EuKRf0hVyN711RAU2ke+LrL7BuBb4sQUbBl5383vKuPOKxmy2tuvJg3AhpDzoiFKCcUrXmSh3EUCoWz2HVSFRJeeSBrmY3Bp1YdEubU1RiHmVDiwiLqu6NKXw4nxnLn2Lutd0qOUQwwLQVj7ynuuWGUvoHRmvYY4iDOOi+p9HCeiohF8ugliH6f4PjhkUCDFlbVUxfc9EQHdXNimOs0p/YCsl2BGCHoA2sMgV22uk+iXkpfeOI5hdJpw0QqiTUL8v0xtZ0jFJMpmGV+nMwJHxRiFj+v5059IJbZf8R3lFjpoCAsNBSnNnlXDsZLQKxFZ/2bSogqEPSc5lm3TLXYTT5nRPTxSHj52wBmHp6slmWXfXJug7Jnvl+isbsya68sS6VcaVRZCDYvG6R6Z1/PdctwwAcOvWfZm47cXxI3dxlZsJVZVcUH2fpyICyOIS1/dfl1HxxMQU7rluGZ5Yt6ru+3Ptkym4arw2CWYY74E81REmATs2XqsTgi5CHtDHBegqV31pz/OnstEm67r2XMHzpvNtr00xemZ149lNl+OJdausQr6vt+KpYd9zYUuyjHaEjh5wyQ9jx5aaIDhe4LapqniyMHRm7aHSakNoEck6cVrc7yyO3cTF5hSl5z4yXsM91y3T3gPtlgsmPDHlmZ1SNwkO9PcZx9l10lTfGlvxHHaMoFchmEPMwwQh1K6rAF2+i7SGzjwEcxJDaBrVVKvdOaPIOitl3O8sbgCZLVeLSwKvs3orxnsgblI7wBsrF+eEPAhPxnm+kdje1KIWCr0GmaNL09BsOkbQ26ommQhyWmRB2puvlR4qQDovkiK4c0bRzMlUN+klCSDT5WoJl6xLMoGFx0KXtCusfw6PVVwXxqSEryVIz5DHJGNL7e0yzuuvvKCh1nOWMiYNHSPo477OhQuFu+SLt9HueVnSuGQWxZ0zimZMpqZJ75qL+urSRgPR6rnNOw4ZM5cG15F0AtP56asFtwP9c1AZLBzjEXxWlz3VlahnTo0NCfromp5BnahsrXt8D6PAoGtK7e06znNnd09/Z0m9hvKgYwS964qa4AVIqLnj4W/TZdkLR7OZ3DBdVQBFVXGkCTAqUvbJVo+vadLbfXDUi4vIofJWWAiptUtdCdrqinSr6ISf0ZWTgHlzyjgyXsM8jdfNyYkpa+pfNTYkuD6XxVy4XoLJi6m3Usa+O0/VMV6xaZd1wRL8qMGNgYMG0Bil/EZeORIS0DGC3lZZKuzXq3sAAsGvqxIVfjiBZCqAIqs40hgri1IhKsvxTTphmIRzdWx8ekUcFI2Jwhb0FdQtCPqaRu0W5YIaJiz8TBXObB4lLi6kcXXzpnO6FCOyHV/tp2msZ3d3FfqttmP86F38y1tNkf3b04xfUcY+q/HVXY8aUJQk+E2FAPzbt56On7w6bp1IbHpwdXyTXncaPXs4XiDuxGgLagL0949tbKOqnrn0zxgcBUxPzq4xDupnXYIUs3gTnREBU0DrX9ujKHrEapG9bnTHB+rfrEwPYByhNDhcddY763SwdwwecE5XrWILDDL1JxDkSe+ruEJLZX5PGcMfW63dl0aoAqdqMLsmcUuyoDDdTyY1VNRYm8giUtoVm6DvGNUN0HrPlSiKouJQyUpAZzH2pr7oXpd1umOTcU8d38Hhap1nxHRshE+cXOyHj9fqVCS6+r2uBEW3w2Po4sNtqoEcdV+lsaGYhshVjWRSp9x9tadCNeV+imvniNPHu69eavzuopIXzu8p11X7Cq7HxXbXDGeGjhL0AUVd2Wfty52WItkMbH1xNcLp6vw2uMBtH2kwANamGOu3j2Du7O7Yagz1gUwbyDM2XkP/xscb3hJsC4TB4SqOvqGpVlWiBrdEV5fPLgKiyqMeMbgvm4TW+u0jDec3CW2bUTQcRZwEm2A1JZyzJS+slEt13klZ2XWydGboOEFfJOEVpmgRq0Vyi7T1Jc4Nz/Bel03ja4qvGBuvJc6VFKTJyCLaNPyWANgXCJt3HNJ6rsyd1W011gYqCt3EGKyqbRG0prcF03eljq+6gtapNfIWfLbjm3IHBZNm1DOc5Llpxpt+pKAnonMAfBFetSgGsIWZPxFqQwA+AeDdAI4DeB8zf9ffdxOAO/ymdzHzFzLrvYYiCS8dRVIvFckt0taXODnS0xi208RS2FRHcVFzrgQCZV6ljDnlroYAJpNaZ2y8hhWbdk1PBuFngpXfQb/DBk2d2gywv4W6fle2ZzJvwWc7vstiLOtnuBlv+i4r+gkAtzHzd/36r0NEtJOZn1LavAvA2/yfSwD8LwCXENHpAO4EsBzevTRERNuZ+XBmVxCiSMKr6DTTZhClTrP1JU7I/rETE3Wuh2FMVbuC9AJpsH08cPN1nbCCVa8a6V0plxpcM22CNXwMW79NE2Tct9A435Xpmcxb8EUdv9mLsWa86buUEnwZwMv+368T0dMA+gCogv4qAF9kz4VnDxH1EtGZAN4JYCczvwYARLQTwGUA7svsCkIU0eBZVJplM3BRp9n6oj4IUSvnsfFG1YeKrmpXuRQvr7ipFKAJ1Wc7nNbABFFjJSPdKjhKsI7XJp36a5uA4gg+ndAypQA3PZN5C76iqVCDPuV5/ljulUS0CMA3ALydmX+ubP8qgE1+fVkQ0T8D+At4gn4OM9/lb//PAMaZ+b9qjn0zgJsBYOHChRc999xziS6oKD7d7ULSYilxcPXzdjWip80zrzuPSR9t0mG7ZoAM1CFAo2ugCVttUVONgqj+VMol67lLRPjR3e+O7FsS5JlsDpm4VxLRaQC2AbhFFfJZwcxbAGwBPD/6pMcp4mxdZNSVRF6GbBd1WhxPKbXPi9c9kujYuklAJ4zUjI3hfkUJbsKp9LQ6TxIdUQU5TEnQogKnbJMZYK5jmwXyTLYeJ0FPRGV4Qv5eZn5Q06QK4Bzl/7P9bVV4q3p1+9eTdDQORTJ4thM217g04xmlTkszwWR5bK+EntcunPROh9re1LcAVxtRVH4Um1otSv0VNRnkiTyTrSWywpTvUfNZAE8z88cNzbYD+APyuBTAEV+3vwPAaiKaT0TzAaz2twkFxOYal6aKVlTFK5unVNpjb3h4JPLYwWSg6pFPTJgFrq59mLCtw9VGZFtZz+8pT3vCBFXP1ApnA/1SfUzQ47KiXwHgvQAOEFHgy/URAAsBgJk/BeBReK6Vz8Bzr/xDf99rRPSXAL7jf25jYJgViofNgyPPwippPKVMxwbsdYTVY8d1yTUFRgVGzxJR3WQy0K/PZx4HgmdIjnpDiVo5ixplZuLidfNNePeZrQ0D+LBh3+cAfC5R74TMcNGBr12zpCG1QEDaYCCbAErjKWXKWRIlVF3UKnG3TzLXGT1t+cx1xT5MBKm1o6JG43jGiGCfWXRcZKzQiKueeqC/DxseHtH7mfvHyUNARLl5xsmBY0oZqztnQNyJxtQ+WMmr6PKZB7h4y4SDmCROREhCpI5eaE9UPe5t9+930oHb9PAMOOnMTfpjG4FuWa2tOafcNX282x88gOrYOBinhHkgJHXX5eKnvnnHoem+xdVbm9qb9OsmITzQ34cn1q0yGkIDF9FwgJQOiRMRbMiKPidamVgtvNJ1EUAuucmjVo1p3TNVA2iQ80Xn2RJMUmnUSUnKxAWE28+rlLUBTgFRQjhO4FrREuMJ7YEI+hxodWI11yyKqgBy+UyUwMqj7qypTy+NjRsjPrsImN1tDxAK9y2u3jpoHzVBugjhOBONGFOFJIigz4FWJ1Zz0deGBZBLmbYogZVH3VkTNg+hKUZDClxT27S6bdsEqat6ZHrTi5tmQAS7EAcR9CkwPbStMpgF/TH5cZSIMMUcK6kYEF2mLeoYujeB8NiZimcA5hzzJvVNn5+FUI34NVVpSqvbNn2nalRsQKvf9ISZS8cZY5MYA5Oex2QkbIXBTO2Pjkq5hL/93Qvx7KbLtcUbdAZGArDiracD8EqsRY2n7hjAqYySur4GY3f0jQmUS3ov3iCVLlAfBORiRA3OpRPyWei243zXaQLDBCENHbWib+aKyfTQ3rJ1H3or5YaMiEmEShyDblwVQpiB/j7sfe61unqnDOCJH52Kb4saz2Bb2EUzyCi597nXsPvgqHYyqk2xNZWvLpWui77aFtyURVKtOMbRNOqjolZNE9qDjhL0zdSN2x7OsfEayl2E+T3lhkIRrpjqpG54eKSh1JytPzoVgondB0cjC2dEjedAv1dSL6yGGa9NRhbNPjJew747VxsLMOuuMUpfbRqXKX+Fv2LTrlTC09U4Ojhcdappq0NUPkJaOkrQN1M3HlVJpzbF6JnVjeGPrTa2sa3STCtRXak5W3/iqItcxymqnWl/1CQS9DXLmgKmY82rlDMTni7GUZPthGBPVBZ8tshV04Ti01E6+qx04y56fpM+WsUmEG06/qjP6vS6SZJVha+zt6dsbKsSNZ5JBLLa1ywTb5mOZSvskQe2yS9KWEs0rJCWjhL0WQiIKAEcoGYKNNFFZJwsogxzUcIyvEp1yVwYdZ02g2iAy3iaDLsmwn2Ney02TMcaM3j45CU8Td+nS3pgiYYV0hKrwlSzWL58Oe/duzfRZ9MarVyrIYXPGRVVGq6oY9JDB+eKKpmXtiKQ6ToDg2jailO6ylXbhqqFqTKU5HtOQ5oqS1KhSXAhkwpT7ULaYJIkr8lhg5yuvF1Yp2rSHRNOrdZtU3DaikCm6wkMomnRfQ/Lzz29MJ4jzU4lkCaiVaJhhbR0nKBPS1JDoCrYXMrb6QSNbQUfJm1FoFYUUW92RKft7a4VwjPN9Us0rJCGjhf0cVU5Waz0XISoTtC4JunKYuXZ6cmxXFwSRXgKM4WOFvRJ/I+zWOm5CtGwoDHpjefOKuGN2tR09aJrLkovoDpdHSAuiYJwikhBT0SfA/DbAF5h5rdr9q8F8PvK8X4FwAK/jOBPALwOYBLAhMlQkBdJH/a0K72kQlQ3QZRLhJMTU9M6+UlmbBuqYvm5p2ci7DtV6IlLoiCcwmVF/3kAnwTwRd1OZt4MYDMAENEVAG4N1YVdycw/S9nPRLTyYbcJ0TsGD+C+J1+YXqHfcMk5uGtgad0EUfXT8KppFAJkZRqtkmuFDUIQikqkHz0zfwOAa0HvGwDcl6pHGeLif9ysJGgBdwwewJf2PF+3Qv/SnuexyD8/cMoP3eZZExWM1cxrajYusQ5ZBl0JQruTWcAUEfUAuAzANmUzA3iciIaI6OaIz99MRHuJaO/o6GgmfYp62F2Do7LkvidfMO4Lzr/h4ZHERUBacU3NxiULZJZBV4LQ7mRpjL0CwBMhtc1vMHOViH4RwE4iOui/ITTAzFsAbAG8gKksOhSlK2+FwS7K/91WVSnAtjKdCUZIV5VcJ9sgBCEOWQr66xFS2zBz1f/9ChE9BOBiAFpBnxe2h73ZOvwsVtVRKYdnghFS9O+CEI9MBD0RzQPwDgA3KtvmAuhi5tf9v1cD2JjF+ZKgM941W2C4JszqrZRxYmIqUch7q4RgM/Old3oMgCBkTaSOnojuA/AtAEuI6EUiej8RfZCIPqg0+x0AjzPzMWXbmwF8k4j2A/g2gEeY+Z+y7LwrJr31yvMWNNVg51rLdf2VFyTWL7fCCNlsu4Do3wUhHh2X1EyHLYFVUHu0GStRUz9stVyT0OxqRM1OECYIQiMzKqmZDpveupkGO5PKIevVaLONkDPBLiAI7UxH5aM3UZR83p2qcijK+AqCoGdGrOiLZLzrRJe/Io2vIAiNzAhB3+kJvFqNjK8gFJsZYYyNotnGS0EQhKyZ8cZYG0lSGQuCILQTM8IYa8Mlb4ogCEI7M+MFvbgGCoLQ6cx41U3eKQNE/y8IQquZ8Sv6tWuWoFyium3lEmXiGjgTUgYLglB8ZrygB+Blzbf9nxDR/wuCUARmvKDfvOMQalP1kr02xZkIY9H/C4JQBGa8oM9TGEtqAEEQisCMF/R5CmOpWyoIQhGY8YI+T2HcqUnMBEFoL2a8e2XeeVo6MYmZIAjtRaSgJ6LPAfhtAK8w89s1+98J4B8BPOtvepCZN/r7LgPwCQAlAJ9h5k0Z9TtTRBgLgtDJuKhuPg/gsog2/8LMy/yfQMiXAPwdgHcBOB/ADUR0fprOCoIgCPGJFPTM/A0AryU49sUAnmHmHzPzSQBfBnBVguMIgiAIKcjKGPvrRLSfiB4jogv8bX0AXlDavOhv00JENxPRXiLaOzo6mlG3BEEQhCwE/XcBnMvMFwL4HwAGkxyEmbcw83JmXr5gwYIMuiUIgiAAGQh6Zv45Mx/1/34UQJmIzgBQBXCO0vRsf5sgCILQRFK7VxLRLwH4KTMzEV0Mb/J4FcAYgLcR0WJ4Av56AL/ncsyhoaGfEdFzCbpzBoCfJfhc3ki/4lPUvkm/4iH9ikeafp1r2uHiXnkfgHcCOIOIXgRwJ4AyADDzpwBcC+BDRDQBYBzA9ezVJ5wgoj8GsAOee+XnmHnEpbfMnEh3Q0R7TaW0Won0Kz5F7Zv0Kx7Sr3jk1a9IQc/MN0Ts/ySATxr2PQrg0WRdEwRBELJgxqdAEARB6HQ6TdBvaXUHDEi/4lPUvkm/4iH9ikcu/SJPnS4IgiB0Kp22ohcEQRBCiKAXBEHocNpO0BPRe4hohIimiMjohkRElxHRISJ6hojWKdsXE9GT/vatRDQro36dTkQ7ieiH/u/5mjYriWif8vMGEQ34+z5PRM8q+5Y1q19+u0nl3NuV7a0cr2VE9C3/+/4eEV2n7Mt0vEz3i7J/tn/9z/jjsUjZd7u//RARrUnTjwT9+jMiesofn38monOVfdrvtIl9ex8RjSp9+ICy7yb/u/8hEd3U5H7do/TpB0Q0puzLZcyI6HNE9AoRfd+wn4jov/t9/h4R/ZqyL/1YMXNb/QD4FQBLAHwdwHJDmxKAHwF4C4BZAPYDON/fdz88X38A+BSAD2XUr78BsM7/ex2Av45ofzq8ZHE9/v+fB3BtDuPl1C8ARw3bWzZeAH4ZwNv8v88C8DKA3qzHy3a/KG3+I4BP+X9fD2Cr//f5fvvZABb7xyk1sV8rlXvoQ0G/bN9pE/v2PgCf1Hz2dAA/9n/P9/+e36x+hdr/CbwYn1zHDMC/A/BrAL5v2P9uAI8BIACXAngyy7FquxU9Mz/NzFGVu7WZM4mIAKwC8IDf7gsABjLq2lX+8VyPey2Ax5j5eEbnNxG3X9O0eryY+QfM/EP/75cAvAIgj0RILplW1f4+AODf++NzFYAvM/MJZn4WwDP+8ZrSL2berdxDe+ClGmkGabLTrgGwk5lfY+bDAHYiOhV6Xv26AcB9GZ3bCEdnAb4KwBfZYw+AXiI6ExmNVdsJekdMmTPfBGCMmSdC27Pgzcz8sv/3/wPw5oj216PxBvsr/7XtHiKa3eR+zSEve+ieQJ2EAo0Xeek1ZsFbrQVkNV4umVan2/jjcQTe+MTK0ppDv1TeD29VGKD7TrPCtW/X+N/RA0QU5L4qxJj5aq7FAHYpm/McMxumfmcyVoUsJUhEXwPwS5pdH2Xmf2x2fwJs/VL/YWYmIqPfqj9TL4WXHiLgdngCbxY8X9q/ALCxif06l5mrRPQWALuI6AA8YZaYjMfr7wHcxMxT/ubE49WJENGNAJYDeIeyueE7ZeYf6Y+QCw8DuI+ZTxDRf4D3RrSqieeP4noADzDzpLKt1WOWC4UU9Mz8mykPYcqc+Sq8V6Juf1UWK6OmrV9E9FMiOpOZX/YF0yuWQ/0ugIeYuaYcO1jdniCi/w3gz5vZL2au+r9/TERfB9APYBtaPF5E9K8APAJvkt+jHDvxeGlwybQatHmRiLoBzIN3P+WZpdXp2ET0m/Amz3cw84lgu+E7zUpoRfaNmV9V/v0MPLtM8Nl3hj779Wb1S+F6AB9WN+Q8ZjZM/c5krDpVdfMd+JkzyfMSuR7AdvasG7vh6ccB4CZ49W6zYLt/PJfjNugFfWEX6MUHAGit83n0i4jmB6oP8lJMrwDwVKvHy//uHoKnu3wgtC/L8dLeL5b+Xgtglz8+2wFcT55XzmIAbwPw7RR9idUvIuoH8GkAVzLzK8p27XeaUb9c+3am8u+VAJ72/94BYLXfx/kAVqP+7TbXfvl9Ow+ecfNbyra8x8zGdgB/4HvfXArgiL+YyWas8rAw5/kD4Hfg6alOAPgpgB3+9rMAPKq0ezeAH8CbjT+qbH8LvAfxGQBfATA7o369CcA/A/ghgK8BON3fvhxeYfSg3SJ4s3RX6PO7AByAJ7C+BOC0ZvULwL/1z73f//3+IowXgBsB1ADsU36W5TFeuvsFniroSv/vOf71P+OPx1uUz37U/9whAO/K+H6P6tfX/OcgGJ/tUd9pE/t2N4ARvw+7AZynfPaP/LF8BsAfNrNf/v/rAWwKfS63MYO3sHvZv59fhGdP+SCAD/r7CV6N7R/5516ufDb1WEkKBEEQhA6nU1U3giAIgo8IekEQhA5HBL0gCEKHI4JeEAShwxFBLwiC0OGIoBcEQehwRNALgiB0OP8fZCI1W25tGuYAAAAASUVORK5CYII=\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "AeC7Uv29KsXC"
      },
      "source": [
        "## DR-Learnerの開始、まずはT-Learner"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "xp2P-IDT9vql",
        "colab_type": "code",
        "outputId": "e918b1c4-f126-4d79-8008-e99dff428e5a",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 229
        }
      },
      "source": [
        "# ランダムフォレストモデルを作成\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",
        "M_0 = RandomForestRegressor(max_depth=3)\n",
        "M_0.fit(df_0[[\"x\"]], df_0[[\"Y\"]])\n",
        "\n",
        "# 介入を受けたモデル\n",
        "M_1 = RandomForestRegressor(max_depth=3)\n",
        "M_1.fit(df_1[[\"x\"]], df_1[[\"Y\"]])\n"
      ],
      "execution_count": 6,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "/usr/local/lib/python3.6/dist-packages/ipykernel_launcher.py:10: DataConversionWarning: A column-vector y was passed when a 1d array was expected. Please change the shape of y to (n_samples,), for example using ravel().\n",
            "  # Remove the CWD from sys.path while we load stuff.\n",
            "/usr/local/lib/python3.6/dist-packages/ipykernel_launcher.py:14: DataConversionWarning: A column-vector y was passed when a 1d array was expected. Please change the shape of y to (n_samples,), for example using ravel().\n",
            "  \n"
          ],
          "name": "stderr"
        },
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "RandomForestRegressor(bootstrap=True, ccp_alpha=0.0, criterion='mse',\n",
              "                      max_depth=3, max_features='auto', max_leaf_nodes=None,\n",
              "                      max_samples=None, min_impurity_decrease=0.0,\n",
              "                      min_impurity_split=None, min_samples_leaf=1,\n",
              "                      min_samples_split=2, min_weight_fraction_leaf=0.0,\n",
              "                      n_estimators=100, n_jobs=None, oob_score=False,\n",
              "                      random_state=None, verbose=0, warm_start=False)"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 6
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "wAeHIJiqOF-h",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# 傾向スコアを求めます\n",
        "from sklearn.linear_model import LogisticRegression\n",
        "\n",
        "# 説明変数\n",
        "X = df[[\"x\"]]\n",
        "\n",
        "# 被説明変数（目的変数）\n",
        "Z = df[\"Z\"]\n",
        "\n",
        "# 回帰の実施\n",
        "g_x = LogisticRegression().fit(X, Z)\n",
        "g_x_val = g_x.predict_proba(X)\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "xTjMfuZTNrLO",
        "colab_type": "text"
      },
      "source": [
        "## DR法に基づく推定"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "jwEhxJQeNvhw",
        "colab_type": "code",
        "outputId": "7c24f63d-2e4e-4cc9-fbc1-96f8ea21daf7",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 406
        }
      },
      "source": [
        "# 処置群\n",
        "Y_1 = M_1.predict(df_1[[\"x\"]]) + (df_1[\"Y\"] - M_1.predict(df_1[[\"x\"]])) / \\\n",
        "    g_x.predict_proba(df_1[[\"x\"]])[:, 1]  # [:,1]はZ=1側の確率\n",
        "df_1[\"ITE\"] = Y_1 - M_0.predict(df_1[[\"x\"]])\n",
        "\n",
        "# 非処置群\n",
        "Y_0 = M_0.predict(df_0[[\"x\"]]) + (df_0[\"Y\"] - M_0.predict(df_0[[\"x\"]])) / \\\n",
        "    g_x.predict_proba(df_0[[\"x\"]])[:, 0]  # [:,0]はZ=0側の確率\n",
        "df_0[\"ITE\"] = M_1.predict(df_0[[\"x\"]]) - Y_0\n",
        "\n",
        "# 表を結合する\n",
        "df_DR = pd.concat([df_0, df_1])\n",
        "df_DR.head()\n"
      ],
      "execution_count": 8,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "/usr/local/lib/python3.6/dist-packages/ipykernel_launcher.py:3: SettingWithCopyWarning: \n",
            "A value is trying to be set on a copy of a slice from a DataFrame.\n",
            "Try using .loc[row_indexer,col_indexer] = value instead\n",
            "\n",
            "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
            "  This is separate from the ipykernel package so we can avoid doing imports until\n",
            "/usr/local/lib/python3.6/dist-packages/ipykernel_launcher.py:7: SettingWithCopyWarning: \n",
            "A value is trying to be set on a copy of a slice from a DataFrame.\n",
            "Try using .loc[row_indexer,col_indexer] = value instead\n",
            "\n",
            "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
            "  import sys\n"
          ],
          "name": "stderr"
        },
        {
          "output_type": "execute_result",
          "data": {
            "text/html": [
              "<div>\n",
              "<style scoped>\n",
              "    .dataframe tbody tr th:only-of-type {\n",
              "        vertical-align: middle;\n",
              "    }\n",
              "\n",
              "    .dataframe tbody tr th {\n",
              "        vertical-align: top;\n",
              "    }\n",
              "\n",
              "    .dataframe thead th {\n",
              "        text-align: right;\n",
              "    }\n",
              "</style>\n",
              "<table border=\"1\" class=\"dataframe\">\n",
              "  <thead>\n",
              "    <tr style=\"text-align: right;\">\n",
              "      <th></th>\n",
              "      <th>x</th>\n",
              "      <th>Z</th>\n",
              "      <th>t</th>\n",
              "      <th>Y</th>\n",
              "      <th>ITE</th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\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",
              "      <td>-0.001272</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>3</th>\n",
              "      <td>0.570717</td>\n",
              "      <td>0.0</td>\n",
              "      <td>1.0</td>\n",
              "      <td>2.245229</td>\n",
              "      <td>0.827965</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>4</th>\n",
              "      <td>0.559952</td>\n",
              "      <td>0.0</td>\n",
              "      <td>1.0</td>\n",
              "      <td>2.139868</td>\n",
              "      <td>0.962872</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>5</th>\n",
              "      <td>-0.454815</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.5</td>\n",
              "      <td>1.973293</td>\n",
              "      <td>0.218243</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>6</th>\n",
              "      <td>-0.447071</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.5</td>\n",
              "      <td>1.953387</td>\n",
              "      <td>0.292839</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "</div>"
            ],
            "text/plain": [
              "          x    Z    t         Y       ITE\n",
              "2 -0.124545  0.0  0.5  2.193123 -0.001272\n",
              "3  0.570717  0.0  1.0  2.245229  0.827965\n",
              "4  0.559952  0.0  1.0  2.139868  0.962872\n",
              "5 -0.454815  0.0  0.5  1.973293  0.218243\n",
              "6 -0.447071  0.0  0.5  1.953387  0.292839"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 8
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "XvOWVBt99vq7",
        "colab_type": "code",
        "outputId": "bcbc6cf8-e22b-4337-ec11-a157ab6ede25",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 337
        }
      },
      "source": [
        "# モデルM_DRを構築し、各人の治療効果をモデルから求める\n",
        "\n",
        "# モデルM_DR\n",
        "M_DR = RandomForestRegressor(max_depth=3)\n",
        "M_DR.fit(df_DR[[\"x\"]], df_DR[[\"ITE\"]])\n",
        "\n",
        "\n",
        "# 推定された治療効果を各人ごとに求めます\n",
        "t_estimated = M_DR.predict(df_DR[[\"x\"]])\n",
        "plt.scatter(df_DR[[\"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": 9,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "/usr/local/lib/python3.6/dist-packages/ipykernel_launcher.py:5: DataConversionWarning: A column-vector y was passed when a 1d array was expected. Please change the shape of y to (n_samples,), for example using ravel().\n",
            "  \"\"\"\n"
          ],
          "name": "stderr"
        },
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "[<matplotlib.lines.Line2D at 0x7f22813d1630>]"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 9
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3dfXyU5Z3v8c8vkwQSUAKCCgHkoRQfmlY0KyitVXcV2p4KPqzCwq7uVrFd6SnVQ4tHrNZDFcs5rd1T2x5xPW4ri6jr5lChm/qA7SkLHOIGRbCUBxUSVBCIqAQyJL/zx9wJk2QmmSQzmQe+79crL2au65r7/uWe4TdXrvu6r9vcHRERyX556Q5ARESSQwldRCRHKKGLiOQIJXQRkRyhhC4ikiPy07XjwYMH+6hRo9K1exGRrPTqq69+4O5DYtWlLaGPGjWKqqqqdO1eRCQrmdk78eo05CIikiOU0EVEcoQSuohIjlBCFxHJEUroIiI5Im2zXEREstHCis0s37CHRndCZsycOIJF08vSHRaghC4ikrCFFZt5cv3ulueN7i3PMyGpa8hFRCRByzfs6VJ5b1NCFxFJUGOc+0fEK+9tSugiIgkyi1Peu2HEpYQuIpKAiupa4nXEPahPN50UFRGJo6K6lu88+xoNjZ0PqSyp3Mb0CaW9EFV8SugiIjFUVNfy7RWbSHR0fG9dfUrjSYQSuohIDPNWbOpS+2ElRSmKJHEaQxcRaeOz9/5bl19z+dkxlyjvVUroIiJtHD7W2OXXPLl+NxPu/21aT452mtDN7HEz22dmb8SpNzP7BzPbYWavm9kFyQ9TRCT1KqprOfee33T79YeOhJm3YhMLKzYnMarEJTKG/gTwU+CXceq/BIwLfiYCPw/+FRHJeAsrNrNsw+64UxLbKgxZp7Nenly/m/KzBvX6rJdOE7q7/97MRnXQZBrwS3d3YL2ZlZjZUHd/N0kxikiGO3z4MF/4whfald9+++3MmTOHd999l6lTp7ar/853vsOsWbPYsWMH1113Xbv6++67j2uuuYbNmzcze/bs1vusD1P8+Zs5cvp5nHL4LT568WecWlTQqs0jjzzC5z//eV566SXuuOOOVq997/BRBn5pHoVnjOXIjg3U/f5X7fY/ZNoCCk4bzidbf8eH658hlGeceWrflhktp19/H/mnDubj13/L4aqVrV478wmjdmsVgwYN4qc//SlLly7l8ccf58ILL+zgSPZMMma5lALRCxnUBGXtErqZzQHmAIwcOTIJuxaRTJCXl8eYMWPalQ8cOBCA/Pz8mPUDBgwAoLCwMGb9qaeeCkCfPn1a1b/7YT279h6mX4PRF/jgqPNRXglnnHYqQwecmG1SXFwMQFXNJ9QcP4Wj4UYKQnmEacJOHQChyBdAXkER+SVntv/FQpEUmdenmPySM5lyXqTNvi3v4YDlhYL6/jFfHwqFWo7DmDFj6Nu3b/t9JJF5An9nBD305939MzHqngcWu/sfgucvAd919w7vAF1eXu66SbRIdtu1axdlZWX84z/+IzNmzOi1/U5e/DK1MeZ9l5YUsXbBFa3KKqprueu5zdSHu36is5kBP77x/JYhlLarLsbz9uKvdHufcWMxe9Xdy2PVJaOHXguMiHo+PCgTkRx37Ngxjhw5gsVb5CRF4l3EE6t8SeW2HiVzgFmTRrYaD29eKrejpD6wuKBdWUV1LUsqt7G3rp5hJUXMnzI+qePsyZi2uBL4m2C2yyTgQ42fi5wcwuEwAAUF7ZNXKsW7iCdWeU+v4Jw9aWTMtc4XTS/j7cVfYfLYQe3qCkLGvV89r1VZ818KtXX1OFBbV89dz21O6jTHRKYtLgfWAePNrMbMvmZmXzezrwdNVgO7gB3AUuDvkxadiGS0hoYGoPcTeryLeGKVl8ToKSeiX2GIh288v9MbVyy79WIevvF8SqO+TMKNzrwVmxi1YBWzlq4D4L6VW9r9pVAfbmRJ5bZuxRdLIrNcZnZS78DtSYtIRLJGunroa/64P+Hyj4+Gu7TtgcUF3PvV87o0FDJ9QilV7xyMOQSzdudBRi1YFfe1yVwDRmu5iEi3DRkyhNtuu42zzjqrV/eb6Bj6worNhJvibyeZJy27e9eiZK4Bo4QuIt32qU99il/84he9vt9hJUUxZ7k4dNgbbmvUglUUF+TxwLWf7fHJye7etWj+lPE92m80reUiIt3W1NREU1MHXeAUSWYSPBJu4o6nN/X45GSoGzN9BhYXZNwsFxE5ST3//POEQiFeffXVXt1vsi+pb3J6fHJy5sQRnTdqo+1MmJ5SQheRbkvXSdFU6O0bVIw7vV/Sv5iU0EWk25oTemFhYZoj6bmenpzsyknRyWMH8cIdl/Vof7HopKiIdFu65qFDJCmu3Xkwadvr6bh8RydFQ2bMnDii0zntPaUeuoh0WzqHXJbdenHMqzQ7c2qfULuyglDPly6Id1I0ZMbOB7+c8mQO6qGLSA985jOf4c4772xZNbGtVK9dsuzWizusn7V0Xate/OSxg3j7QD2Hj7UeLw83Oksqt/UotpkTR8S8sKg7J0u7K6HVFlNBqy2K5LZYqxwWFYR48NqyXr/xQ7SO5qn39EKjhRWbWb5hD43uKRtmSfVqiyJykjp69CgNDQ2ccsop7VZcjLXKYfPaJelM6Km0aHpZrwytxKMxdBHptocffpgBAwZw9OjRdnVdWeJWkkMJXUS6raNpi11Z4rY3dXTyMtspoYtItzU0NGBmLbdaizZ/yniKClqXFxWEknrZfnfEO0nZmycvU0Vj6CLSbeFwOO6UxeZx8u//eguHjkR68n3y09+HbB7jTvXJy3RIKKGb2VTgJ0AIeMzdF7epPwt4HBgCHARmu3tNkmMVkQzT0NDQ6VWih48eb3lcVx9m3opNzFuxKeF99CsM8YNrkjszJt0nL1Ol04RuZiHgEeBKoAbYaGYr3X1rVLP/DvzS3f/JzK4AHgT+OhUBi0jmmDp1Kqeffnrc+rv/dTONTT2bGv1JQyPzVmyi6p2DOZmEkymRHvpFwA533wVgZk8B04DohH4ucEfweA1QkcwgRSQzXXXVVVx11VVx6z9p6NnNmaMtW7+b8rMGMX1CacovWMpWiST0UiB61ZkaYGKbNq8B1xIZlrkGOMXMTnP3A0mJUkQy0oEDBwiHw5x55pkp35dzYonb6AuWmm+2DMlfVjfbJOsMxX8Bvmhm1cAXgVqg3Vezmc0xsyozq9q/P/Y9AUUke3z729/m4ovjX36f7ImAtXX1HV6wdLJLJKHXAtHzeYYHZS3cfa+7X+vuE4C7g7K6thty90fdvdzdy4cMiX3XbhHJHuFwuMOTorMmjUzq/kJmumCpA4kk9I3AODMbbWaFwAxgZXQDMxtsZs3buovIjBcRyXENDQ0drrS4aHoZs5OY1BvdM/aCpUzQ6Ri6ux83s7lAJZFpi4+7+xYzux+ocveVwGXAg2bmwO+B21MYs4hkiM566NC9KYKTF78c8ybQpcEJ0FiLfnXlgqVcPama0Bi6u69290+7+1h3/0FQ9r0gmePuz7r7uKDNLe5+LJVBi0hm6OjCop6IdZUpQN2RyA01Hry2jNKSIoxIku/KCo7Nq0DW1tXjnDip2tObRGcCXSkqIt02Z86cmAtz9dT0CaVUvXOw3frizXPSjcisl9Ju9K5zeRVIJXQR6bZrrrkmZdt+/rV349Y1X6rUnSmLuXxSNf0LK4hI1tq1axc1NalZ5aOuPpxQu65OWczlk6pK6CLSbTfccAO33XZbusPoUu86U1eBTAYNuYhItyWyOFd3DSwuaFmlsTNd6V03D83k4iwXJXSRKLk6nS1VOprl0tNjee9Xz2P+s68Rbux4ca/u9K6nTyjNyfdVCV0k0PamxlojpHPxEnoyjmXbnnRJcQHukbH1kBmN7t2a5ZLLlNBFArk8nS1VohP6worNLNuwG4/Toe7OsczVnnSqKKGLBGJdmdhRucBDDz3E0KFDmbV0HWt3Huy0fS5MDcxkmuUiEsjlmwenyowZMzh06qcSSuaQG1MDM5l66CKBxjhjBfHKBTZu3Miif9kJnJJQ+8vP1iqrqaQeukigNE7vMV65wKWXXsr2Nc8k3H7NH3UfhFRSQhcJ5PIFJ8lWUV3L5MUvc/RYA4QS/0NfY+ippSEXkUAuX3CSTM1TEo80HAdvwvLar4oYj8bQU0sJXSSKpsl1rmV6Z9NxACwvsTSiv3ZSTwldRLqkeRqnN0YSekdDLroAqHcllNDNbCrwEyJ3LHrM3Re3qR8J/BNQErRZ4O6rkxyriGQQCxUweNoCCoec1a5u9qSRXb5LkfRcpwndzELAI8CVQA2w0cxWuvvWqGYLgafd/edmdi6wGhiVgnhFJENYKJ9+Z3++VVnIjJkTRyiZp0kiPfSLgB3uvgvAzJ4CpgHRCd2BU4PHA4C9yQxSRDJH8zBKU/gYx2q2UjDkLPL7DyJkxs4Hv5zu8E5qiUxbLAX2RD2vCcqi3QfMNrMaIr3zb8bakJnNMbMqM6vav1/zUUWy0cyJIwBo/OQQ+56+h6Nv/UerckmfZM1Dnwk84e7DgS8DvzKzdtt290fdvdzdy4cM0RVjItlo0fQyZk8aSZ43AZAXyteYeYZIZMilFoj+6h0elEX7GjAVwN3XmVlfYDCwLxlBivQWrYeemEXTy5g5Lo/PPAo/nf1n3KBknhESSegbgXFmNppIIp8B/FWbNruBPweeMLNzgL6AxlQkq2g99M5Ff+GdeiRyL9FU3bFIuq7TIRd3Pw7MBSqBN4nMZtliZveb2dVBszuBW83sNWA5cLO7VjSS7NLReuhy4guvtq4eB/bVfQJA1e7D6Q1MWiQ0Dz2YU766Tdn3oh5vBSYnNzSR3hVvnRGtPxLR9guvYFApp19/H5X7+rEojXHJCVqcSyQQb50RrT8S0faLLa9vf4rGlvPB8b5pikjaUkIXCWi1xY61/WJr/PgQR7b9O0MKw2mKSNpSQhcJTJ9QyoPXllFaUoQRWQf9wWvLdEI00PYL79j7O9hf8QDXj9OSUJlC74RIFK22GF/b5YUH9c1jP3BV2fD0BiYtlNBFJGHRX3hPP/0JNz4JBQUFaY5KmmnIRUS6JRyOjJ0roWcO9dBFJCEV1bXct3ILdfWRRN60rRpQQs8k6qGLSKcqqmuZ/8xrLckcwIdPoHT2Q2x8vymNkUk0JXQR6dSSym2Em1pf/B3qV0J+6Xn85JV30hSVtKWELiKdqo1xtWzD/nf4+I2XqD2gS/8zhRK6iHSooroWi1Fev6uKA6t+zJmnaAw9Uyihi0iHllRuI+ZKe8FNou+cem6vxiPxaZaLSApETiJuIhx1vnDy2EEsu/Xi9AXVTfEWJ/OmSEK//s9G9WI00hEldJEku/JHr7B93yftytfuPMioBativqaoII8Hr/1sRl6lOqykKOYYev8C4+NQCLNYAzKSDkroIglovrFDrMSWDPXhJuat2MS8FZtalae7V19RXcuRhuPtyosKQpw98lT+TXPQM4oSukgnFlZs5sn1u9Oy77U7DzJr6bq0JPW2d3BqVlJUwH1Xn8fEMz/H/ju/0etxSXwJJXQzmwr8BAgBj7n74jb1PwYuD54WA6e7e0kyAxVJh4rq2rQl82Zrdx5My35j3cEJoF+f/JahoaFDh/Z2WNKBThO6mYWAR4ArgRpgo5mtDO5SBIC7fzuq/TeBCSmIVaTXLKzYzPINe2jMkDspjlqwisKQ8cPrP9dunL2iupbv/3oLh45EruJs7kH3dDw+3snQ5mGnl156iZqaGm666aYe7UeSJ5Ee+kXADnffBWBmTwHTgK1x2s8E7k1OeCK9a+/evZx3wUQOHz5xsUzBaSM4c9ZDALz/zL00vLu91WsKzxzHGTd8H4D3ln2X8IE9rer7jixjyPS7AHj3iW9x/HDr+6cXjb2IwV+ZB0Dto3NoOvpxq/p+51zKoCu/TkOjc/2ln6VfPhTmR2YcNxxvIu+cP2fApTfhjWFqHrmJPcB1D0H/PvkU5ufxrW99i4ULF3Lw4EHGj29/s467776befPmsXv3bi688MKW8kNHGnCHkkv/hlPOn0rDB7t5/58XRGL6eYjGY0c444wzlNAzSCIJvRSI/oTWABNjNTSzs4DRwMtx6ucAcwBGjhzZpUBFesP27dupe7+GonGTCPU/DYD8/oNa6vuO/Bz5A85s9Zr8AWe0PC4afQEFQ0a1qi847cR64UVjymlsk7ALTx/d8rh43CSawsda1w8dd6J+/GSKC4xrL4hs87n/qOH4kDGRSsuj+OwvtLR1g/KxgykrK4tsp7CQG264od3vfPbZZwPQr1+/VvW/WvcODhQMGgZAXp/ilu3nGdw86SwuueSSdtuT9DHv5E9KM7semOrutwTP/xqY6O5zY7T9LjDc3b/Z2Y7Ly8u9qqqqe1GLpEg4HGbMHU+T16cflt+zGRwGzJo0kkXTy2LWV1TX8t1/eZ1jx1O3uFVByFgSY5gmEfGmWDZ7e/FXuhuW9ICZveru5bHqEumh1wIjop4PD8pimQHc3rXwRDLHqjf2EerX+fn8h288v8dj1G3vjjR58ctJnxYZbnSWVG7rVqwhs4w5hyCJSeTS/43AODMbbWaFRJL2yraNzOxsYCCwLrkhivSOiupavv2zCg79/lc0HvkwbrvJYwel5AKgWDepToaufEnMWrqOUQtWMWrBKiXzLNRpQnf348BcoBJ4E3ja3beY2f1mdnVU0xnAU97ZGI5IhlpSuY2P9u7k8LoVNDXEToKpvNAn+ibVyVZRHe+P6oiFFZsZtWBVwlMkUxGj9FxC89DdfTWwuk3Z99o8vy95YYn0vr119XhjZOqfhU6MnxvwVi+NF7cdhulsHDtR81Zsouqdg+3G82ctXdflee5FBSHmT2k/W0bST1eKigSGlRTx4fEGACy/sFV5ukweOyhpFxY9uX53Ui6SevDasoxcc0a0fK5Ii/lTxpPvkXVLmme4pLs3uuzWi5k8dlDnDXtJaUmRknkGUw9dJDB9QikrPj2Qp16GvFABpSVFzJ8yPu0JrCtj9qlcd6YgZBpqyXCdzkNPFc1Dl0zU2NjIsWPHKC4uTnco3TbmrlU0Jfm/9cDiAu79as+XE5Ce6+k8dJGTRigUyupkDvBXE0cmpZee7qV7pes0hi4SZfny5dxzzz3pDqNHFk0v6/G4++xJI5XMs5ASukiUF198kSeeeCLdYfTYslsv5uEbz2dgceLLFxQX5PHwjefz9uKvxF2uQDKbhlxEohw9epS+ffumO4ykaDunXXKfeugiUXIpocvJRwldJIoSumQzJXSRKO5OUZHWKZHspDF0kSirV6/uvJFIhlIPXUQkRyihi0RZsGABP//5z9Mdhki3KKGLRHnmmWdYu3ZtusMQ6ZaEErqZTTWzbWa2w8wWxGlzg5ltNbMtZvbPyQ1TpHfUffQJL2w7yOgFq5i8+OVObwwhkkk6PSlqZiHgEeBKoAbYaGYr3X1rVJtxwF3AZHc/ZGanpypgkVSpqK7lw4/rKW7Ko4jIrdvuem4zgC7QkayQSA/9ImCHu+9y9wbgKWBamza3Ao+4+yEAd9+X3DBFUm9J5TaajjdgoRM3t6gPN7KkclsaoxJJXCIJvRTYE/W8JiiL9mng02a21szWm9nUWBsyszlmVmVmVfv37+9exCIpsreunrzCIqxPcbtykWyQrHno+cA44DJgOPB7Mytz97roRu7+KPAoRNZDT9K+RZJiWEkRPvdXMctFskEiPfRaYETU8+FBWbQaYKW7h939LeBPRBK8SNaYP2U8RQWhVmXpvgWdSFckktA3AuPMbLSZFQIzgJVt2lQQ6Z1jZoOJDMHsSmKcIil3+Zj+DFj7MMX73sCI3D9TN0SWbNLpkIu7HzezuUAlEAIed/ctZnY/UOXuK4O6q8xsK9AIzHf3A6kMXCTZPvroIzasqWTp0uu55ZavpDsckS5LaAzd3VcDq9uUfS/qsQN3BD8iWeno0aMAWm1RspauFBUJNCf0Pn36pDkSke5RQhcJqIcu2U4JXSRgZowePZoBAwakOxSRbtF66CKBCy+8kF27NDlLspd66CIiOUIJXSTwyiuvcOWVV/L222+nOxSRblFCFwnU1tby4osvEg6H0x2KSLcooYsENMtFsp0SukhACV2y3Uk1y2XW0nWs3XmwVdnsSSNZNL0sTRFJJlFCl2x30vTQYyVzgCfX72ZhxeY0RCSZZuDAgZSVlelKUclaFlmGpfeVl5d7VVVVSra9sGIzy9bvJtHfzAzeelCLMYlI5jOzV929PFZdzg25xOuJdyRN32kiIkmVU0MuFdW1XU7mIs0eeughpk6NefdEkayQ9T30iupavv/rLRw60v25w8UFOfW9Jt20Y8cOXn/99XSHIdJtWZ3QK6pruePpTTT1YMgkz+CBaz/b41i6M9QTy+Sxg1h268U93k6uqKiuZUnlNvbW1TOspIj5U8an7A5CR48e1QwXyWoJJXQzmwr8hMgdix5z98Vt6m8GlnDiXqM/dffHkhhnTPet3NKjZA7Q5DBvxSbmrdgUt02/whA/uCYytfHOpzfRmMIx97U7DzJqwap25QV5sOQvzz+pbodWUV3L/GdeIxy8ybV19cx/5jWAlBwHJXTJdp3OcjGzEJGbPl9J5GbQG4GZ7r41qs3NQLm7z010x8mY5RIr8Z2sSooKuO/q83Iq4Z///d9SV99+KK2kqIBN916VtP00z4p6/9nv0/jxQYbe/BNdnyAZq6ezXC4Cdrj7rmBjTwHTgK0dviqFfvjDH/L8mn9n/zt1LWWh/gMZdMUtANStXU74QE1LXWG+cfvVk1m0aBGjFqzi0O+e4PiH+1tts2DwCEoumQHAwZeW0vhJXav6PmeO5dSLrgXgQOUjNB070qq+7/BzOOWC/wTAB6sexhtbJ6K+Z32OUz53Fe5NfPDr/9HudyoaW07/8y6nKXyUA7/5n+3q+42fTPH4S2is/4iDL/yiXX39eZcxrz7M3Mde4tCax1vKh5xSyCVjB3PLLbdwxRVXsHPnTu655552r587dy6XXHIJW7Zs4Qc/+EG7+vnz5zNhwgSqqqr40Y9+1K7+nnvu4ZxzzuEPf/gDP/vZz9rVP/DAA4waNYoXX3yRxx8/EV/NoSO8+e5H9P3839J/4GkcfHM9H2/9XbvXn3bVN8jr259P/vgH9v9pHf1XPIADxYUh+vUJ4V/4BpZfyMdvvEz9rlcZPbiYzw4vASLrnC9btgyApUuXsmbNGl6vqePtA0dwB8svZPCXv0XBaSMI9R8ERK5PAJTUJaskktBLgT1Rz2uAiTHaXWdmlxLpzX/b3fe0bWBmc4A5ACNHjux6tIHt27ezcWMVx6PGW/JLzmx5HD5QQ8N725v3yYD+hbz55qAT9R/sIXygdXiWd+LEaHj/Oxw/vK9VfV6f4hP1+96isf5wq/pQ/4Etjxve34Efb2hVXzBw6In6ILZohaePjjxoaopZ36f07KD+eMz6xlHnA+DHG1rV174Hz78dYtq0aVRU13LfEy+z5d/aJ8wXw5+meOUhjtW+yQcx6l+mjL4j93L07U0c+G37+jUFf0afobs4sn0Dh9a0r3+l+DcUDh7JJ1t/R90f2tef0dhAfbiJ458civn7eVNj5Pf8+ECr+gagDhj6+SYMaPzoAxre286292DbG0EjM9YuWEUecOD3L3Hkj/+31bbzCosAGHj537UqX75hjxK6ZJVEhlyuB6a6+y3B878GJkYPr5jZacDH7n7MzG4DbnT3Kzrabk+HXEYvWNXphUPNY9/RwxAappGueHuxLjiTzNLRkEsi8/VqgRFRz4dz4uQnAO5+wN2PBU8fAy7sTqBdMaykKG5dyIzZk0ay5f6p7caUx53eL9WhSY4ImaU7BJEuSWTIZSMwzsxGE0nkM4C/im5gZkPd/d3g6dXAm0mNMob5U8Zz13ObqQ83tpQVFYR48NqyDk8MvnDHZVz5o1fYvu+TlMRlwKwunlCrqK7ljhWbaEpJRNJdMyeO6LyRSAbpNKG7+3EzmwtUEpm2+Li7bzGz+4Eqd18J/Gczuxo4DhwEbk5hzEBk2lrVOwdZvmEPje6EzLjuwtKEZnm8cMdlCe2jeQ50bV09ITMa3SlNwVzo6RNixx29f+ldmuUi2ShrF+dqO0cZoCDPWPKXn8upqXudWVixuWVGRrYryIPjTbQ6N5Jn9Phag64qLSli7YIOTwGJpE1HY+hZm9B7a45yNkvGlbRtGSS8imWiYp287kh3/3LJg06HtRIZthNJp5xcbTFWMu+o/GTUnJT+63OvcyR8IpUVF+TxwLWfzdqkFW+Iqit6azhNpDdlbUKXxCQj+eUiHRfJRVm7zODA4oIulYuI5LqsTej3fvU8CkKt5wkXhIx7v3pemiISEUmvrB1yaf5zubeWVhURyXRZm9BB46AiItGydshFRERaU0IXEckRSugiIjlCCV1EJEcooYuI5AgldBGRHKGELiKSI5TQRURyhBK6iEiOSCihm9lUM9tmZjvMbEEH7a4zMzezmGv1iohI6nSa0M0sBDwCfAk4F5hpZufGaHcK8C1gQ7KDFBGRziXSQ78I2OHuu9y9AXgKmBaj3X8DHgKOJjE+ERFJUCIJvRTYE/W8JihrYWYXACPcfVVHGzKzOWZWZWZV+/fv73KwIiISX49PippZHvAj4M7O2rr7o+5e7u7lQ4YM6emuRUQkSiIJvRYYEfV8eFDW7BTgM8ArZvY2MAlYqROjIiK9K5GEvhEYZ2ajzawQmAGsbK509w/dfbC7j3L3UcB64Gp3r0pJxCIiElOnCd3djwNzgUrgTeBpd99iZveb2dWpDlBERBKT0B2L3H01sLpN2ffitL2s52GJiEhX6UpREZEcoYQuIpIjlNBFRHKEErqISI5QQhcRyREJzXLJJBXVtSyp3MbeunqGlRQxf8p4pk8o7fyFIiI5LqsSekV1LXc9t5n6cCMAtXX13PXcZgAldRE56WXVkMuSym0tybxZfbiRJZXb0hSRiEjmyKqEvreuPmZ5bZxyEZGTSVYl9GElRXHrFlZs7sVIREQyT1Yl9PlTxsetW75hT9w6EZGTQVYl9I5OfDa692IkIuy2LmkAAAg6SURBVCKZJ6sSOkDIrEvlIiIni6xL6DMnjuhSuYjIySKr5qEDLJpeBkTGzBvdCZkxc+KIlnIRkZOVeQJjz2Y2FfgJEAIec/fFbeq/DtwONAIfA3PcfWtH2ywvL/eqKt3USESkK8zsVXePeYvPTodczCwEPAJ8CTgXmGlm57Zp9s/uXubu5wM/JHLTaBER6UWJjKFfBOxw913u3gA8BUyLbuDuh6Oe9gM05UREpJclMoZeCkRP8q4BJrZtZGa3A3cAhcAVsTZkZnOAOQAjR47saqwiItKBpM1ycfdH3H0s8F1gYZw2j7p7ubuXDxkyJFm7FhEREkvotUD0nMDhQVk8TwHTexKUiIh0XSIJfSMwzsxGm1khMANYGd3AzMZFPf0KsD15IYqISCI6HUN39+NmNheoJDJt8XF332Jm9wNV7r4SmGtmfwGEgUPATakMWkRE2kvowiJ3Xw2sblP2vajH30pyXCIi0kVZd+m/iIjEpoQuIpIjlNBFRHKEErqISI5QQhcRyRFK6CIiOSLr1kNvVlFdy5LKbeytq2dYSRHzp4zv8BZ1IiK5LisTekV1LXc9t5n6cCMAtXX13PXcZqDj+46KiOSyrBxyWVK5rSWZN6sPN7KkcluaIhIRSb+sTOh76+q7VC4icjLIyoQ+rKSoS+UiIieDrEzo86eMp6gg1KqsqCDE/Cnj0xSRiEj6ZeVJ0eYTn5rlIiJyQlYmdIgkdSVwEZETsnLIRURE2lNCFxHJEUroIiI5QgldRCRHKKGLiOQIc/f07NhsP/BON18+GPggieEki+LqmkyNCzI3NsXVNbkY11nuPiRWRdoSek+YWZW7l6c7jrYUV9dkalyQubEprq452eLSkIuISI5QQhcRyRHZmtAfTXcAcSiursnUuCBzY1NcXXNSxZWVY+giItJetvbQRUSkDSV0EZEckbEJ3cz+0sy2mFmTmcWd3mNmU81sm5ntMLMFUeWjzWxDUL7CzAqTFNcgM3vBzLYH/w6M0eZyM9sU9XPUzKYHdU+Y2VtRdef3VlxBu8aofa+MKk/n8TrfzNYF7/frZnZjVF1Sj1e8z0tUfZ/g998RHI9RUXV3BeXbzGxKT+LoRlx3mNnW4Pi8ZGZnRdXFfE97Ka6bzWx/1P5viaq7KXjft5vZTb0c14+jYvqTmdVF1aXyeD1uZvvM7I049WZm/xDE/bqZXRBV1/Pj5e4Z+QOcA4wHXgHK47QJATuBMUAh8BpwblD3NDAjePwL4BtJiuuHwILg8QLgoU7aDwIOAsXB8yeA61NwvBKKC/g4TnnajhfwaWBc8HgY8C5Qkuzj1dHnJarN3wO/CB7PAFYEj88N2vcBRgfbCfViXJdHfYa+0RxXR+9pL8V1M/DTGK8dBOwK/h0YPB7YW3G1af9N4PFUH69g25cCFwBvxKn/MvAbwIBJwIZkHq+M7aG7+5vu3tldny8Cdrj7LndvAJ4CppmZAVcAzwbt/gmYnqTQpgXbS3S71wO/cfcjSdp/PF2Nq0W6j5e7/8ndtweP9wL7gJhXwvVQzM9LB/E+C/x5cHymAU+5+zF3fwvYEWyvV+Jy9zVRn6H1wPAk7btHcXVgCvCCux9090PAC8DUNMU1E1iepH13yN1/T6QDF8804JcesR4oMbOhJOl4ZWxCT1ApsCfqeU1QdhpQ5+7H25Qnwxnu/m7w+D3gjE7az6D9h+kHwZ9bPzazPr0cV18zqzKz9c3DQGTQ8TKzi4j0unZGFSfreMX7vMRsExyPD4kcn0Rem8q4on2NSC+vWaz3tDfjui54f541sxFdfG0q4yIYmhoNvBxVnKrjlYh4sSfleKX1jkVm9iJwZoyqu939//R2PM06iiv6ibu7mcWd9xl885YBlVHFdxFJbIVE5qJ+F7i/F+M6y91rzWwM8LKZbSaStLotycfrV8BN7t4UFHf7eOUiM5sNlANfjCpu9566+87YW0i6XwPL3f2Ymd1G5K+bK3pp34mYATzr7o1RZek8XimV1oTu7n/Rw03UAiOing8Pyg4Q+VMmP+hlNZf3OC4ze9/Mhrr7u0EC2tfBpm4A/tXdw1Hbbu6tHjOz/w38l96My91rg393mdkrwATgX0jz8TKzU4FVRL7M10dtu9vHK4Z4n5dYbWrMLB8YQOTzlMhrUxkXZvYXRL4kv+jux5rL47ynyUhQncbl7geinj5G5JxJ82sva/PaV5IQU0JxRZkB3B5dkMLjlYh4sSfleGX7kMtGYJxFZmgUEnnzVnrkLMMaIuPXADcByerxrwy2l8h2243dBUmtedx6OhDzbHgq4jKzgc1DFmY2GJgMbE338Qreu38lMrb4bJu6ZB6vmJ+XDuK9Hng5OD4rgRkWmQUzGhgH/L8exNKluMxsAvC/gKvdfV9Uecz3tBfjGhr19GrgzeBxJXBVEN9A4Cpa/6Wa0riC2M4mcoJxXVRZKo9XIlYCfxPMdpkEfBh0WpJzvFJ1trenP8A1RMaRjgHvA5VB+TBgdVS7LwN/IvINe3dU+Rgi/+F2AM8AfZIU12nAS8B24EVgUFBeDjwW1W4UkW/dvDavfxnYTCQxPQn07624gEuCfb8W/Pu1TDhewGwgDGyK+jk/Fccr1ueFyBDO1cHjvsHvvyM4HmOiXnt38LptwJeS/HnvLK4Xg/8HzcdnZWfvaS/F9SCwJdj/GuDsqNf+XXAcdwB/25txBc/vAxa3eV2qj9dyIrO0wkTy19eArwNfD+oNeCSIezNRM/iScbx06b+ISI7I9iEXEREJKKGLiOQIJXQRkRyhhC4ikiOU0EVEcoQSuohIjlBCFxHJEf8fgo8bWEYeUKAAAAAASUVORK5CYII=\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "4riPBjmmWX__"
      },
      "source": [
        "以上"
      ]
    }
  ]
}