5_3_doubly_robust_learning_issue18.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "colab": {
      "name": "5_3_doubly_robust_learning_issue18.ipynb",
      "provenance": [],
      "collapsed_sections": []
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "aoxI3DOK9vm2"
      },
      "source": [
        "# 5.3 Doubly Robust Learningの実装\n",
        "\n",
        "本ファイルは、5.3節の実装です。\n",
        "\n",
        "5.2節と同じく、人事研修の効果について因果推論を実施します。"
      ]
    },
    {
      "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": "b67517df-c4d8-4c40-851b-669d4aefabdc",
        "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(-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 0x7fb171d483c8>"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 3
        },
        {
          "output_type": "display_data",
          "data": {
            "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",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "BHcdUlW9koTa"
      },
      "source": [
        "## データをまとめた表を作成し、可視化する"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "1EMwdGIIIPrK",
        "outputId": "c3575b0b-6d6d-4bbf-b59c-321d7e5978e2",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 195
        }
      },
      "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>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>"
            ],
            "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"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 4
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "L6Tb2Hjk9vno",
        "outputId": "dd55119a-872c-4dc1-ebc7-0f66a3bab881",
        "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 0x7fb171858438>"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 5
        },
        {
          "output_type": "display_data",
          "data": {
            "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",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "AeC7Uv29KsXC"
      },
      "source": [
        "## DR-Learnerの開始、まずはT-Learner"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "xp2P-IDT9vql",
        "outputId": "122ee358-9b12-4166-94aa-38ad5ef19208",
        "colab": {
          "base_uri": "https://localhost:8080/"
        }
      },
      "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"
      },
      "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": 7,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "xTjMfuZTNrLO"
      },
      "source": [
        "## DR法に基づく推定"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "jwEhxJQeNvhw",
        "outputId": "a56db28c-a884-489c-b597-0d7190281fbf",
        "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>0</th>\n",
              "      <td>-0.616961</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.5</td>\n",
              "      <td>1.803183</td>\n",
              "      <td>0.514190</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",
              "      <td>0.081865</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.333970</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.364906</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>9</th>\n",
              "      <td>0.751865</td>\n",
              "      <td>0.0</td>\n",
              "      <td>1.0</td>\n",
              "      <td>2.289369</td>\n",
              "      <td>0.776072</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "</div>"
            ],
            "text/plain": [
              "          x    Z    t         Y       ITE\n",
              "0 -0.616961  0.0  0.5  1.803183  0.514190\n",
              "2 -0.124545  0.0  0.5  2.193123  0.081865\n",
              "5 -0.454815  0.0  0.5  1.973293  0.333970\n",
              "6 -0.447071  0.0  0.5  1.953387  0.364906\n",
              "9  0.751865  0.0  1.0  2.289369  0.776072"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 8
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "XvOWVBt99vq7",
        "outputId": "97da2ae4-09f7-4411-a4ef-1f177442fa60",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 338
        }
      },
      "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 0x7fb16b4c6d68>]"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 9
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3de3iU5Z3/8fc3kwQip3AIp0AAEUEUC25+otW2irVQvUSqbsFqt+5a7clu1ZafWA+1alda9rJ2W7uuutrt1oraakrVbtZj+S0VJDQogpuKhyJB0XKySAg5fH9/zDNhJjOTTMJM5sDndV25mLnve+b55pnwzZ37ue/7MXdHRETyX1G2AxARkfRQQhcRKRBK6CIiBUIJXUSkQCihi4gUiOJsHXjEiBE+ceLEbB1eRCQvrVu37i/uXpGoLmsJfeLEidTV1WXr8CIiecnM/pysTkMuIiIFQgldRKRAKKGLiBQIJXQRkQKhhC4iUiCyNstFRCTTauobWVbbwLbdTYwtL2Px3KksmFWZ7bAyRgldRApSTX0j1z66gaaWNgAadzdx7aMbAAo2qWvIRUQK0rLaho5kHtHU0say2oYsRZR5SugiUpC27W7qUXkhUEIXkYI0trysR+WFQAldRArS4rlTKSsJxZSVlYRYPHdqliLKPF0UFZGCFLnwqVkuIiIFYMGsyoJO4J11O+RiZveZ2Xtm9kqSejOzfzGzzWb2spmdkP4wRUR6pqa+kVOWPsukJU9wytJnqalvzHZIGZfKGPrPgHld1H8amBJ8XQ7866GHJSLSe5E56I27m3AOzkEv9KTe7ZCLu680s4ldNDkX+Lm7O7DazMrNbIy7v5OmGEUkh/3Xf/0X7777Lueddx6DBw/mlVdeSXivg4ULF1JWVkZ9fT0vvfRSXP1FF11ESUkJL774Ips2bYqrv+SSSwBYtWoVr732WkxdSUkJF110EQC///3vWfyTJ9i1r6Wj3or7wTEfY1ltAwN3vMrWrVtjXj9kyBA+85nPAPC73/2O7du3x9QPHz6cc845B4AVK1awc+fOmPrRo0czb1643/vrX/+av/71rzH148eP54wzzoj7ntLO3bv9AiYCrySpexw4Ner5M0B1kraXA3VAXVVVlYtIftu+fbsDDvhrr73m7u4/+MEPOsqiv9555x13d7/hhhsS1u/du9fd3a+88sq4OjPrOOall14aVz948OCO+oULF8bVhwYO9wnXPO4Tr3nczzrrrLj6qVOndrz+Yx/7WFx9dXV1R/3MmTPj6k8//fSO+smTJ8fVz58/P23nHKjzJLnawvVdC3roj7v7cQnqHgeWuvv/BM+fAa5x9y5vR1RdXe26Y5FIfnvjjTeYPHkyy5Yt4xvf+AYlJSXs2bOHXbt2xbUdP348oVCIXbt2sWfPnrj6qqoqioqK2LlzJx988EFcfeSWlX/5y1/Yu3dvTJ2ZMWHCBADef/995v3zU2z/YH9UgyKKB1dQWV7Go39/LE1NsYuLSkpKqKwMXzx95513aG5ujqkvLS1l7NixAGzbto0DBw7E1Pfv35/Ro0cDsHXrVlpbW2Pqy8rKGDVqVNz31Btmts7dqxPVpWOWSyMwPur5uKBMRApcS0t4WGPMmDGUlJQA4eGLIUOGAOGx7Ose28CHB9qAjRhw0UlV3LpgRtL3HDZsGMOGDUtaP2LECEaMGJG0vqKighsWfSJmHxc4OAe9u8Q6ZsyYLusjiT2ZcePGdVmfSelI6CuAK8xsOTAb2OMaPxc5LER6opFkHr274ZCyEj7Y30J71CCAA79YvQWgy6R+qA7HOehA90MuZvYgcBowAtgOfAcoAXD3u8zMgJ8QngmzD/j77oZbQEMuIoVgx44dPPjgg8ybN49X/loW1ytOJmTG67ed1QcRFp6uhlxSGkPPBCV0kcJyytJnaezFxleVveg919Q38u1HX2ZfS3vXMU0exgOXndzjmHKZErqIZMTevXt58803mTRpEjNu/T2ZziYGPT6GAT9cOLNghlu6SujanEtEem3dunUcf/zxvPjiiwwpK8n48XrzC8OBKx9az/U1G9IdTs5RQheRXovMcln95m52N7V00zq7Hli9peBXiiqhi0ivReZj/+v/eythvfVhLN1xKOi7FYESuogcgkgPvZVQwnoH+odyJ60X8t2KQAldRA5BJKFbKPmSlv1t2Zl4kUgh360IlNBF5BBUV1dz1AWLCQ1KvnIzoqwkxB0LZ/LW0rO5Y+FMjihJT/rp/L7J/h4wKOi7FYFucCEih2DixIm0Tv5EkgGXWE0tbSyrbei46URvpxFGr0btvAI08m+iBU7RY+g9OfZF97zAqtd3dt8wgSkjB/DU1af16rW9oXnoItJr7777Lmd+9xE+GDgeC3U/bdGAN5eenfG4Ikk/2UKnIqAdKDJitibIhMH9Qrz83a5uKdEzmocuIhlRU1PDK3f9I6Ut+1Jq31dj2AtmVXY5vBJZX5rpZA7wQXNbn82BV0IXkV6LXBS9Yf4MKrtJ1pHdDvtKLk1RfHDN231yHCV0Eem1SEI/92+qWLVkTtKkHjLjtvNm9Ony+1yaotjm3if3N1VCF5FeiyT0yPa5ycas29z7fC+VXJui2Bf3N1VCF5FeiyT00tJSINwTTyRZeSYtnjuVotxZ0xQjMuMn3TRtUUR67YILLmDatGmEQuGJi21JZs0lK8+kuj/vTPtFzzu62bWxJ1sIZ2JISAldRHpt2rRpTJs2reN5ZXlZwoTW3QXTTEh2ITLRFrzpmlq4eO5UrnxofUptizMwPqKELiK99uqrr7Jt2zbOOOMMIJzQkt3Ls68l+6vAgbcyNBd+waxKHqnbktJCpJb28Hz5dF5bSOl3hJnNM7MGM9tsZksS1E8ws2fM7GUze97MsneXVBHpM3fddRfnn39+x/MFsyq57bzwFEYj3DPv69ktEdkaz3/gspO5+KQqUjlMusfRu+2hm1kIuBM4E9gKrDWzFe6+KarZPwM/d/f/MLM5wG3A59MaqYjknJaWlo4ZLhGHsqw/nS6cPb7jhtSdyzPt1gUzeOyPjXx4oOv7q/bmln1dSaWHfiKw2d3fcPcDwHLg3E5tpgPPBo+fS1AvIgUoUULPFbcumMHFJ1V19MhDZlx8UhW3LpiR8WNfX7Oh22QO4fH8dE5fTGUMvRKIvrqwFZjdqc1LwHnAj4DPAIPMbLi774huZGaXA5cDVFVV9TZmEckRuZzQIZzU+yKBd5bqytDIhmHp+osmXddZvwV8wszqgU8AjUDcryd3v9vdq929uqKiIk2HFpFsyfWEni09maaZzumLqfTQG4HoQadxQVkHd99GuIeOmQ0Eznf33ekKUkRy05IlS9i9W//VOwuZpZzU07miNZWEvhaYYmaTCCfyRcDnohuY2Qhgp7u3A9cC96UtQhHJWTNmJB/O6Grf8kKX7IJsqMhoi1rtlO4pnd0mdHdvNbMrgFogBNzn7hvN7Gagzt1XAKcBt5mZAyuBr6UtQpE+dDgnod5YuXIlra2tzJkzJ6a8pr4xZj56ZP8S6NnNJfJVZNz+wTVv0+ZOyIwLZ4+nesKwjP586QYXIoHOSQjCPahszaPOB3PnzmXPnj2sXr06pjzZEvjK8jJWLZkTVy6p0w0uRFKwrLYh7rZlmdpEqVAkuyia7EJfLm1pW4i09F8koCTUcy0tLeze384pS5+NGUYYm2RPl1zb0rbQqIcuEkiWbJSEkntv94e8uXN/3F7fp0+roKwk9tbR2drT5XCihC4SWDx3qpJQDzXu3Eu7xZ6zppY2nvvf93NmT5fDiYZcRAKRZKNZLqkbPO8bWFEorrxxd1PO7OlyOFFCF4miJNQzZSMnJVxAk407FImGXETkEHywaSX7334lrjwbdygS9dBFYmhhUc98sPJ+SsbNoP/442LKs3GHIlFCF+lwuK9u7I0BxdBeHDsPXReSs0dDLiIBLSzquSJv4+QpIzWbJUeohy4S0MKinmtpaWHq2KHcoeX8OUE9dJGAFhb1nPZDzy3qoYsEcumO9fli3bp1lJeXZzsMCSihiwS0sKjnjjnmmGyHIFGU0EWiaGFR6tra2vjxj3/MqaeeSnV1wt1cpY9pDF1EeqW5uZmrrrqKZ555JtuhSEA9dJEoWliUupaWFgBdFM0hKfXQzWyemTWY2WYzW5KgvsrMnjOzejN72czOSn+oIpkVWVjUeSvYmvrGbl97OFJCzz3dJnQzCwF3Ap8GpgMXmtn0Ts2uBx5291mEbyL903QHKpJpWljUM0rouSeVHvqJwGZ3f8PdDwDLgXM7tXFgcPB4CLAtfSGK9A0tLOoZJfTck8oYeiXwdtTzrcDsTm1uAv7bzL4ODAA+meiNzOxy4HKAqqqqnsYqklG6bVrPjB07li1btmgeeg5J1yyXC4Gfufs44CzgP80s7r3d/W53r3b36oqKijQdWiQ9dMeinikuLmb8+PEMGjQo26FIIJWE3giMj3o+LiiLdinwMIC7vwD0B0akI0CRvrJgVqVum9YD27dv55ZbbqGhQdcYckUqQy5rgSlmNolwIl8EfK5Tmy3AGcDPzOwYwgn9/XQGKtIXtLAodY2Njdx4440cf/zxTJ2qv2JyQbc9dHdvBa4AaoFXCc9m2WhmN5vZ/KDZN4HLzOwl4EHgEnfdskSkkOmiaO5JaWGRuz8JPNmp7Maox5uAU9IbmojkMiX03KOl/yLSK0rouUcJXUR6RQk992gvFxHplTlz5rBz504GDhyY7VAkoIQuIr1SXFzM0KFDsx2GRFFCF5GU1NQ38u1HX2ZfSzsALe+9wYTd9dTctZSRI0dmOToBjaGLSApq6hu5+uH1Hckc4MBftvDCY/fznUfWZDEyiaaELiLdWlbbQHunlSXeFt6ZcsWG97TFcI5QQheRbiXacdLbWwGwomIWP/KSknoOUEIXkW4l3HGyPdxDt1CIlnbnphUb+zgq6UwJXUS6tXjuVIostszbwj10isJzK3Y3tfRxVNKZZrmISLciG5ZFz3IZ9DfnMOiEsyF+p2zJEn0SIpKSBbMq2XTLpxlQGt4z3sywohBm4a57pFyyRwldRHqkJBROG/s2v8iO//4pkY1VI+WSPRpyEZFu1dQ3sqy2gW27m4jMXmze1sDe+t8x/FNfBWCPxtCzTgldJIHoBDa2vIzFc6f26MYXNfWNLH5kPVHrcAC4+KQqbl0wI83RZlZNfSPXPrqBppa22Ir2VggdHGbRvVezTwldpJNwMn6JlmAlTePuJhY/8hJASkn9zNuf57X3PkxY94vVWwDyKqkvq22IT+aEZ7lYKLzTou69mhtSSuhmNg/4ERAC7nX3pZ3qfwicHjw9Ahjp7roVuOSlm1Zs7EjmES3tzpUPrefKh9Yf8vs/sHpLXiX0RIuKAGhvw4pCVPbiLxjJjG4TupmFgDuBM4GtwFozWxHcpQgAd78qqv3XgVkZiFWkT2R6PnW+3ZtxbHkZjQmS+sD+pZRXDGXVkjlZiEoSSaWHfiKw2d3fADCz5cC5wKYk7S8EvpOe8ET61mN/3Epz4//SfmAfoSOGUDpqMgD7t7x8cCFNIDRwGKUVEwFoems9eOyAeWjQCEpHVOHu7H+rPqbupp/u4XOfrOa7v9/B//zpPfZveTkuluLyMUyfehS//epJrFy5Mq5+6tSpTJgwgb179/KHP/whrn769OmMGzeOPXv2sGZN/AZaxx9/PKNHj2bHjh2sW7curn7WrFlUVFRw2f8Zzo3//jjNbQe/v36hIn7wTzdz8SeOjXudZJG7d/kFXEB4mCXy/PPAT5K0nQC8A4SS1F8O1AF1VVVVLpJLHvvjVp942Z1OuBPtZUfN9gnXPO4Trnnci44o7yiPfA2YflpHvRX3i6sfOPPTPuGax71q8W/i6gAffOJ5PuGax338lQ8lrB9y6kU+4ZrHvfIrP0tYf/vtt7u7+6ZNmxLW33PPPe7uvmbNmoT1y5cvd3f3p59+OmH9E0884Y/9catP+/zNCetXrlyZtc/qcAbUeZJ8ne6LoouAX7l7/BUUwN3vBu4GqK6uzre/PKXALattYN9fdwMw7Mwv07/qIx11Iy+4sWN3wYjQEYM7Ho9a9L2O+dgd9QOCy0hWxKiLlsUdr3jQsHB1Sf/E9YNHBMcZ0lE/77hRXP7x8F8NEydOBGDChAmsWrUq7vWTJ4fbHXPMMQnrjz76aGrqG1lSuyvu+POOG8XOARP43qMb2Dt0Skd9v+IivnLaZD5+dAXHHXdc3HtKdqWS0BuB8VHPxwVliSwCvnaoQYlkw7bdTdAaHj8vHXUUJSMO/tj3G3N0l6/tVzktaZ2Z0X/cMcnri0Jd1xeXdNQ/vxueX7ErqNkF1Cd93dAjXuI75xzLglmVfPSjH42rr6lv5JuPvERbUVnc8Z/fDc//7i0AQmWDCEXV37+5iCWXxL+fZF8qCX0tMMXMJhFO5IuAz3VuZGbTgKHAC2mNUKSPjC0vY3OomOIho7DS/J9TvWtfS9pm5kRrbm3nonte4IHLTk7r+8qh6zahu3urmV0B1BKetnifu280s5sJj+WsCJouApZ75787RfLE4rlTufbDA5RN/PeOsiKg09ogykpC3HbejITT9DovSDp9WkXH3PNCsur1ndkOQRJIaQzd3Z8EnuxUdmOn5zelLyyRvhdJ0JGEPKSshA8PtNLedrCPYsD5f1OZdM71glmJ6woxqUvu0W46IlFsSx1lT/8T6/7vyQzoV0xLW+wfnA4897/v9+g9b10wgzsWzmToESVdtrv4pCreWno2by09mykjB/Q0dBEt/ReJ9tZbb/HMM88AJFxM01V5V5L13JN56urTuL5mQ8727E+ZPCzbIUgC6qGLRNm/fz8A/fv3z3Ik4Z79W0vPzrnkecrkYbogmqPUQxeJ0tzcDEC/fv2yHMlBkeR50T0vxFyMLCmCZX87M6bnX1PfyE0rNvZ6+4JRg0rZ/tcDceVDjyjpmAIpuUsJXSTK/v37KS4uJhQKETKjLcGkrZBZgldmXiq94p4O7Uhh0ZCLSJSRI0dywgknAHDh7PEJ2yQrF8k2y9a08erqaq+rq8vKsUVSdX3NBh5c8zZt7oTMuHD2+Lza+lYKj5mtc/fqhHVK6CIi+aOrhK4hF5Eo1113HRdccEG2wxDpFSV0kSgNDQ28+uqr2Q5DpFeU0EWiNDc359SURZGeUEIXibJ///6cWFQk0htK6CJRmpubldAlb2lhkUiU4447Tgld8pYSukiUn/70p9kOQaTXlNBFonS+QcXiuVO1lF7yhsbQRQI19Y18bsGneWXFv+GEt8m99tEN1NQnu4WuSG5JKaGb2TwzazCzzWa2JEmbz5rZJjPbaGa/TG+YIpm3rLaB5h1v077vg46yppY2ltU2ZDEqkdR1O+RiZiHgTuBMYCuw1sxWuPumqDZTgGuBU9x9l5mNzFTAIpmybXcT3tqCFZfElYvkg1R66CcCm939DXc/ACwHzu3U5jLgTnffBeDu76U3TJHMG1teBm0tWKgkvlwkD6SS0CuBt6Oebw3Koh0NHG1mq8xstZnNS/RGZna5mdWZWd377/fsvowimfatTx2Ntx7Aiks7yspKQiyeOzWLUYmkLl0XRYuBKcBpwIXAPWZW3rmRu9/t7tXuXl1RUZGmQ4ukxznHj6b6Y59k5PjJGFBZXsZt583QLBfJG6lMW2wEonf0HxeURdsKrHH3FuBNM/sT4QS/Ni1RivSB4uJi1q58KtthiPRaKj30tcAUM5tkZqXAImBFpzY1hHvnmNkIwkMwb6QxThER6Ua3Cd3dW4ErgFrgVeBhd99oZjeb2fygWS2ww8w2Ac8Bi919R6aCFsmExsZGxowZw/Lly7MdikivpLRS1N2fBJ7sVHZj1GMHrg6+RPJSU1MT7777Lq2trdkORaRXtFJUJLB//34A7YcueUsJXSQQSejabVHylRK6SKC5uRlQD13ylxK6SGDo0KFceOGFVFZq3rnkJ22fKxKYPn06v/yl9pWT/KUeuohIgVBCFwksX76cQYMGsXnz5myHItIrSugigX379rF3715KSkq6byySg5TQRQKRWS6atij5SgldJKCFRZLvlNBFAlpYJPlOCV0k8JGPfIQvfvGLlJaWdt9YJAdZeF+tvlddXe11dXVZObaISL4ys3XuXp2oTj10kUBbWxvZ6uCIpINWiiZQU9/ITSs2sruppaOsyKDdwYDIf3kDLjqpilsXzMhGmJJmV1xxBY8++ijbt2/PdigivaKEHrjonhdY9frOpPXtQRaP7r858IvVWwCU1PPY9TUbeGD1Ft5ftZn9Te1MXPIEF+sXteQhJXS6T+bd+cXqLR2JPRVHlBTxT+cdr5sP54DrazZ0fHbe2oIVhxcV6Re15KOUErqZzQN+BISAe919aaf6S4BlHLx59E/c/d40xplRh5LMe2NfSztXPrSeKx9a32W7UyYP44HLTu6jqA5PD655++CTthYsVBpTp4Qu+aTbhG5mIeBO4ExgK7DWzFa4+6ZOTR9y9ysyEOMhmf29p9j+1wNx5UUGn5tdRfWEYVmIKjWrXt/JxCVPUFlexuK5U1Pq0dfUN7KstoFtu5sY24PXHa7aoi6CeuuBjh565zqRfJBKD/1EYLO7vwFgZsuBc4HOCb3PXHrppTz33HMxZUceeSRPP/00AAsXLmTt2rVs291EWzD4XVIxkZHn3wDAe7/6Li1/2cL3//Xg60srp1FxzmIA3n3w27Ttib0w1n/iTIbP+zoA7/z8Ktr3fRBTX3bUbIZ98nIAtt37Vby1Oab+iOmnMfTjn8e9nW3/dlnc9zTwI3MZcvJnaW/exzv3fz2mrhH4h6fP5crq+bTt3cW7v/hW3OuHfHQhA4//FC273uG9h64Hwr99P/t9Y9iAUu745++zaNEiXn75ZRYsWBD3+h/96Eecc845/OEPf+Diiy+Oq7/33nuZM2cOTz31FF/60pfi6pcvX86JJ57Ib37zG6666qq4+t/+9rcce+yxPPDAA9xwww1x9c888wyTJk3i7rvvZunSpXH1a9asoaKigh/+8If8+Mc/jqt/5ZVXOOKII7jlllu4//77Y+rMjNdffx2Aa6+9loceeoh9B9rYva+FtvZ2rLSMsf/wE8qOPhlvOfi5hczijiOSy1JJ6JVA1N+lbAVmJ2h3vpl9HPgTcJW7v925gZldDlwOUFVV1fNoA9OnT6elpSWmbMyYMR2Pjz/+ePr168ejf9za8Q0WDxndUV868kiK+g+MeX3J8PEdj/uNPoq2QcOT14+ZSnvzh53qxx2sr5yGt8XGV1J+8Pj9xk2P+56Kh4wKPygqSlw/uCL8IFScsD40MByvFZfG1beXFjN6dPj4AwcO5NRTT417fUVF+P3Ly8sT1g8bFv5LZvjw4QnrBw8e3PE+ieoHDBgAhD+nRPVlZWUAVFZWJqyPbJhVVVWVsD4UCgHhX+yd6y0qMR911FFMmH4C6/68i+Jypxiw4vBS/0EfmRvzugtnj0ckn3S7sMjMLgDmufsXg+efB2ZHD6+Y2XBgr7s3m9mXgIXuPqer9830wqKa+sZux6i7YxAzbBF9AS3fvLX07GyH0OH6mg08sGYLnX/0ykqK2N/SzqEMdFx8UngYrfO0Uwgvumjvwfto/FxyUVcLi1LpoTcC0V2VcRy8+AmAu++Ienov8IOeBpluy2obDun15WUlrP/Op2LKIv/Bf7lmS8c0xnwxcckTaXmfyHz8TGhqSTXdJtfVjKNU3r2sJMRt583QdQfJS6kk9LXAFDObRDiRLwI+F93AzMa4+zvB0/nAq2mNshcadzf1+rUlRcZN849NWHfrghm97rkl65nmk3z7RdZTTS1tLKttUEKXvNRtQnf3VjO7AqglPG3xPnffaGY3A3XuvgL4RzObD7QCO4FLMhhzt2rqG2NWdPbE0CNK+M45x2bkP3Qqvwxq6hv59qMvsy8NvVXpnW2H0BkQyaaC3Jxr1s3/za59Ld037OSOhTNzpmcWWb1Y4B3inFRZXsaqJV1eAhLJmkMdQ89p19ds4ME1b9PmTsiMk44c2mUyv2PhTK5+eH3M0EGRwe2fzZ1kDj0b2omee96/pIjm1vaCHxrJlLKSEIvnTs12GCK9ktcJvfOskzb3Lld9VpaXdSTtQlp8s2BWZcL4a+ob+eYjL3XMxc916Zjl0pXuZrn0ZAGXSC7K64Qes2w7BZGeV7IEWGgi3+N1j23gwwNtaXnPdM9yGTWolDXXnXlI79HVXjxK0nI4yeuE3pOl2eVlJYflf+rD4ZeX9rsRCcvrG1ykujS7yEg6DVFEpFDkdUJPdWl2u0Pdn/t2R0URkb6W1wm9Jwt8ejreLiKSb/I6oUP4olcqtBWqiBS6vE/oqc4Z1laoIlLo8j6hL5hVyZSRA7ptp61QRaTQ5X1CB3jq6tM4ZXLyOw8NKA1pK1QRKXgFkdAhPBf5joUzKSsJxZSXlYT43meUzEWk8OX1wqLOCnFZv4hIqgoqocPhsTJSRCSRghlyERE53Cmhi4gUiJQSupnNM7MGM9tsZku6aHe+mbmZJdx8XUREMqfbMXQzCwF3AmcCW4G1ZrbC3Td1ajcI+AawJhOBRnS+RZsZXDRbd2gXEUmlh34isNnd33D3A8By4NwE7W4Bvg/sT2N8MWrqG7n64fUx99t0D9/p/fqaDZk6rIhIXkgloVcC0TtbbQ3KOpjZCcB4d38ijbHFWVbbkPTmCtp8S0QOd4d8UdTMioDbgW+m0PZyM6szs7r333+/x8fq6m7s2nxLRA53qST0RiB6I5RxQVnEIOA44Hkzews4CViR6MKou9/t7tXuXl1RUdHjYMd2sbOiNt8SkcNdKgl9LTDFzCaZWSmwCFgRqXT3Pe4+wt0nuvtEYDUw393r0h3s4rlTKUqSt4vMqalvTFwpInIY6Dahu3srcAVQC7wKPOzuG83sZjObn+kAoy2YVcntn51JaSg+q7e0w+JHXlJSF5HDlnmWxp6rq6u9rq53nfhTlj5LY5Lx9MryMlYtmXMooYmI5CwzW+fuCdf65OVK0a4ujnZVJyJSyPIyoXd1cbSrOhGRQpaXCX3x3KmUJBhHLymylG9JJyJSaPJy+9zI9rjf/e1Gdu1rASFeU4UAAAhLSURBVKC8rISb5h+rrXNF5LCVlwkdtO+5iEhneTnkIiIi8ZTQRUQKhBK6iEiBUEIXESkQSugiIgVCCV1EpEAooYuIFAgldBGRApG3C4sgfI/RZbUNbNvdxNjyMhbPnarFRiJy2MrbhF5T38i1j26gqaUNgMbdTVz7aPhG0UrqInI4ytshl2W1DR3JPKKppY1ltQ1ZikhEJLvyNqEn2/dc+6GLyOEqpYRuZvPMrMHMNpvZkgT1XzazDWa23sz+x8ympz/UWMn2Pdd+6CJyuOo2oZtZCLgT+DQwHbgwQcL+pbvPcPeZwA+A29MeaSeL506lrCQUU1ZWEtJ+6CJy2ErlouiJwGZ3fwPAzJYD5wKbIg3c/YOo9gOAjN+oNHLhU7NcRETCUknolcDbUc+3ArM7NzKzrwFXA6VAwrs0m9nlwOUAVVVVPY01jvZEFxE5KG0XRd39TnefDFwDXJ+kzd3uXu3u1RUVFek6tIiIkFpCbwTGRz0fF5QlsxxYcChBiYhIz6WS0NcCU8xskpmVAouAFdENzGxK1NOzgdfSF6KIiKSi2zF0d281syuAWiAE3OfuG83sZqDO3VcAV5jZJ4EWYBfwhUwGLSIi8VJa+u/uTwJPdiq7MerxN9Icl4iI9FDerhQVEZFYebs5F2i3RRGRaHmb0LXboohIrLwdctFuiyIisfI2oWu3RRGRWHmb0LXboohIrLxN6NptUUQkVt5eFNVuiyIisfI2oYN2WxQRiZa3Qy4iIhJLCV1EpEAooYuIFAgldBGRAqGELiJSIMw94/dzTnxgs/eBP/fy5SOAv6QxnHRRXD2Tq3FB7samuHqmEOOa4O4J7+GZtYR+KMyszt2rsx1HZ4qrZ3I1Lsjd2BRXzxxucWnIRUSkQCihi4gUiHxN6HdnO4AkFFfP5GpckLuxKa6eOaziyssxdBERiZevPXQREelECV1EpEDkbEI3s781s41m1m5mSaf3mNk8M2sws81mtiSqfJKZrQnKHzKz0jTFNczMnjKz14J/hyZoc7qZrY/62m9mC4K6n5nZm1F1M/sqrqBdW9SxV0SVZ/N8zTSzF4LP+2UzWxhVl9bzleznJaq+X/D9bw7Ox8SoumuD8gYzm3socfQirqvNbFNwfp4xswlRdQk/0z6K6xIzez/q+F+MqvtC8Lm/ZmZf6OO4fhgV05/MbHdUXSbP131m9p6ZvZKk3szsX4K4XzazE6LqDv18uXtOfgHHAFOB54HqJG1CwOvAkUAp8BIwPah7GFgUPL4L+Eqa4voBsCR4vAT4fjfthwE7gSOC5z8DLsjA+UopLmBvkvKsnS/gaGBK8Hgs8A5Qnu7z1dXPS1SbrwJ3BY8XAQ8Fj6cH7fsBk4L3CfVhXKdH/Qx9JRJXV59pH8V1CfCTBK8dBrwR/Ds0eDy0r+Lq1P7rwH2ZPl/Be38cOAF4JUn9WcDvAANOAtak83zlbA/d3V919+7u+HwisNnd33D3A8By4FwzM2AO8Kug3X8AC9IU2rnB+6X6vhcAv3P3fWk6fjI9jatDts+Xu//J3V8LHm8D3gMSroQ7RAl/XrqI91fAGcH5ORdY7u7N7v4msDl4vz6Jy92fi/oZWg2MS9OxDymuLswFnnL3ne6+C3gKmJeluC4EHkzTsbvk7isJd+CSORf4uYetBsrNbAxpOl85m9BTVAm8HfV8a1A2HNjt7q2dytNhlLu/Ezx+FxjVTftFxP8wfS/4c+uHZtavj+Pqb2Z1ZrY6MgxEDp0vMzuRcK/r9ajidJ2vZD8vCdsE52MP4fOTymszGVe0Swn38iISfaZ9Gdf5wefzKzMb38PXZjIugqGpScCzUcWZOl+pSBZ7Ws5XVu9YZGZPA6MTVF3n7r/p63giuoor+om7u5klnfcZ/OadAdRGFV9LOLGVEp6Leg1wcx/GNcHdG83sSOBZM9tAOGn1WprP138CX3D39qC41+erEJnZxUA18Imo4rjP1N1fT/wOafdb4EF3bzazLxH+62ZOHx07FYuAX7l7W1RZNs9XRmU1obv7Jw/xLRqB8VHPxwVlOwj/KVMc9LIi5Yccl5ltN7Mx7v5OkIDe6+KtPgs85u4tUe8d6a02m9n9wLf6Mi53bwz+fcPMngdmAb8my+fLzAYDTxD+Zb466r17fb4SSPbzkqjNVjMrBoYQ/nlK5bWZjAsz+yThX5KfcPfmSHmSzzQdCarbuNx9R9TTewlfM4m89rROr30+DTGlFFeURcDXogsyeL5SkSz2tJyvfB9yWQtMsfAMjVLCH94KD19leI7w+DXAF4B09fhXBO+XyvvGjd0FSS0ybr0ASHg1PBNxmdnQyJCFmY0ATgE2Zft8BZ/dY4THFn/VqS6d5yvhz0sX8V4APBucnxXAIgvPgpkETAFePIRYehSXmc0C/g2Y7+7vRZUn/Ez7MK4xUU/nA68Gj2uBTwXxDQU+RexfqhmNK4htGuELjC9ElWXyfKViBfB3wWyXk4A9QaclPecrU1d7D/UL+AzhcaRmYDtQG5SPBZ6MancW8CfCv2Gviyo/kvB/uM3AI0C/NMU1HHgGeA14GhgWlFcD90a1m0j4t25Rp9c/C2wgnJh+AQzsq7iAjwbHfin499JcOF/AxUALsD7qa2YmzleinxfCQzjzg8f9g+9/c3A+jox67XXB6xqAT6f55727uJ4O/h9Ezs+K7j7TPorrNmBjcPzngGlRr/2H4DxuBv6+L+MKnt8ELO30ukyfrwcJz9JqIZy/LgW+DHw5qDfgziDuDUTN4EvH+dLSfxGRApHvQy4iIhJQQhcRKRBK6CIiBUIJXUSkQCihi4gUCCV0EZECoYQuIlIg/j8qlDYa6ALtuAAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "4riPBjmmWX__"
      },
      "source": [
        "以上"
      ]
    }
  ]
}