6_multiple_interventions_in_2D.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "76bc83dc-72d6-471e-be74-252a257c6a65",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "90e2880e-3dca-4a46-81e6-8384346895c5",
   "metadata": {},
   "source": [
    "# Original Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "4606cc95-4fce-4e95-8a96-ed47d9ba6a20",
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 30\n",
    "C = np.linspace(-3, 3, N) ** 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "f6b050ef-eaef-4993-ae3d-206f5a93151d",
   "metadata": {},
   "outputs": [],
   "source": [
    "x01 = np.linspace(0, 10, N) + 2*C\n",
    "x02 = x01 + 2 * x01 - C\n",
    "y0 = x01**2 + x02 + C"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "292c80fb-0ba6-48d6-bee6-41370cdb8a5a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEGCAYAAACUzrmNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAAAT8ElEQVR4nO3dfYxc133e8e/DlVcW4gp6oyRClEq1JZpIbvPiraqFgWJhJhXrGJZQQABbJGRbAUQNObEAtY7o1G1TF5XaAKls1DIgSK4p2K3AOilEGFFcdZNFAXf1soziKpSiirBimRH1YhlGVBTV8uXXP+ZSGa+WPEtmZ2dn9vsBFnPnzL2X54BcPnPOuffcVBWSJJ3NpmFXQJK0/hkWkqQmw0KS1GRYSJKaDAtJUtMFw67AoFxxxRW1bdu2YVdDkkbKoUOHvl9Vm5eWj21YbNu2jYWFhWFXQ5JGSpLvLlfuMJQkqcmwkCQ1GRaSpCbDQpLUZFhIkpoMC0lSk2EhSWNifn6ee++9l/n5+VU/99jeZyFJG8n8/Dw7duxgcXGRyclJZmdnmZ6eXrXz27OQpDEwNzfH4uIiJ0+eZHFxkbm5uVU9v2EhSWNgZmaGyclJJiYmmJycZGZmZlXP7zCUJI2B6elpZmdnmZubY2ZmZlWHoMCwkKSxMT09veohcZrDUJKkJsNCktRkWEiSmgwLSVKTYSFJajIsJElNAw+LJBNJnk3yje79ZUmeSPJS93pp3777khxJ8mKSW/rKP5Tkue6zLyTJoOstSfoza9Gz+BTwQt/7e4DZqtoOzHbvSXIDsAu4EdgJPJBkojvmS8BeYHv3s3MN6i1J6gw0LJJsBX4eeKiv+FZgf7e9H7itr/zRqnqnql4GjgA3JdkCXFxV81VVwCN9x0iS1sCgexb3A58GTvWVXVVVxwC61yu78muA7/Xtd7Qru6bbXlouSVojAwuLJB8D3qiqQys9ZJmyOkv5cn/m3iQLSRbefPPNFf6xkqSWQfYsPgx8PMkfA48CH0nyVeD1bmiJ7vWNbv+jwLV9x28FXu3Kty5T/h5V9WBVTVXV1ObNm1ezLZK0oQ0sLKpqX1Vtrapt9Cauf7eqfgE4COzpdtsDPNZtHwR2JbkwyfX0JrKf7oaq3k5yc3cV1O6+YyRJa2AYq87eBxxIcgfwCnA7QFUdTnIAeB44AdxZVSe7Yz4BfAW4CHi8+5EkrZH0LjAaP1NTU7WwsDDsakjSSElyqKqmlpZ7B7ckqcmwkCQ1GRaSpCbDQpLUZFhIkpoMC0lSk2EhSWoyLCRJTYaFJKnJsJAkNRkWkqQmw0KS1GRYSJKaDAtJUpNhIUlqMiwkSU2GhSSpybCQJDUZFpKkJsNCktRkWEiSmgwLSVKTYSFJajIsJElNhoUkqcmwkCQ1GRaSpCbDQpLUZFhIkpoMC0lSk2EhSWoyLCRJTYaFJKnJsJAkNRkWkqQmw0KS1GRYSJKaBhYWSd6f5Okk305yOMmvdeWXJXkiyUvd66V9x+xLciTJi0lu6Sv/UJLnus++kCSDqrck6b0G2bN4B/hIVf0k8FPAziQ3A/cAs1W1HZjt3pPkBmAXcCOwE3ggyUR3ri8Be4Ht3c/OAdZbkrTEwMKiev5P9/Z93U8BtwL7u/L9wG3d9q3Ao1X1TlW9DBwBbkqyBbi4quarqoBH+o6RpHVtfn6ee++9l/n5+WFX5c/lgkGevOsZHAL+CvDFqnoqyVVVdQygqo4lubLb/Rrgyb7Dj3Zlx7vtpeWStK7Nz8+zY8cOFhcXmZycZHZ2lunp6WFX67wMdIK7qk5W1U8BW+n1Ej54lt2Xm4eos5S/9wTJ3iQLSRbefPPNc66vJK2mubk5FhcXOXnyJIuLi8zNzQ27SudtTa6GqqofAnP05hpe74aW6F7f6HY7Clzbd9hW4NWufOsy5cv9OQ9W1VRVTW3evHk1myBJ52xmZobJyUkmJiaYnJxkZmZm2FU6b4O8Gmpzkku67YuAnwX+CDgI7Ol22wM81m0fBHYluTDJ9fQmsp/uhqzeTnJzdxXU7r5jJGndmp6eZnZ2ls997nMjPQQFg52z2ALs7+YtNgEHquobSeaBA0nuAF4BbgeoqsNJDgDPAyeAO6vqZHeuTwBfAS4CHu9+JGndm56eHumQOC29C4zGz9TUVC0sLAy7GpI0UpIcqqqppeXewS1JajIsJElNhoUkqcmwkCQ1GRZLjMut+ZK0mga63MeoGadb8yVpNdmz6DNOt+ZL0moyLPqM0635krSaHIbqc/rW/Lm5OWZmZhyCkqSOYbHEuNyaL0mryWEoSVKTYSFJajIsJElNhoUkqcmwkCQ1GRaSpCbDQpLUZFhIkpoMC0lSk2EhSWoyLCRJTYaFJKnJsJAkNRkWkqQmw0KS1GRYSJKaDAtJUpNhIUlqMiwkSU2GhSSpybCQJDUZFpKkJsNCktTUDIskn0xy6VpURpK0Pq2kZ3E18EySA0l2JsmgKyVJWl+aYVFV/wzYDjwM/APgpST/JslfHnDdJEnrxIrmLKqqgNe6nxPApcDXk/y7AdZNkgZqfn6ee++9l/n5+WFXZd27oLVDkl8G9gDfBx4C/mlVHU+yCXgJ+PRgqyhJq29+fp4dO3awuLjI5OQks7OzTE9PD7ta69ZKehZXAH+3qm6pqv9SVccBquoU8LEzHZTk2iS/l+SFJIeTfKorvyzJE0le6l4v7TtmX5IjSV5Mcktf+YeSPNd99gXnTST9ec3NzbG4uMjJkydZXFxkbm5u2FVa11YyZ/HPq+q7Z/jshbMcegK4u6p+ArgZuDPJDcA9wGxVbQdmu/d0n+0CbgR2Ag8kmejO9SVgL725k+3d55J03mZmZpicnGRiYoLJyUlmZmaGXaV1rTkMdb6q6hhwrNt+O8kLwDXArcBMt9t+YA74la780ap6B3g5yRHgpiR/DFxcVfMASR4BbgMeH1TdJY2/6elpZmdnmZubY2ZmxiGohoGFRb8k24CfBp4CruqChKo6luTKbrdrgCf7DjvalR3vtpeWL/fn7KXXA+G6665bxRZIGkfT09OGxAoN/A7uJB8AfhO4q6r+9Gy7LlNWZyl/b2HVg1U1VVVTmzdvPvfKSpKWNdCwSPI+ekHxtar6ra749SRbus+3AG905UeBa/sO3wq82pVvXaZckrRGBhYW3RVLDwMvVNVv9H10kN6luHSvj/WV70pyYZLr6U1kP90NWb2d5ObunLv7jpEkrYFBzll8GPhF4Lkkf9CVfQa4DziQ5A7gFeB2gKo6nOQA8Dy9K6nurKqT3XGfAL4CXERvYtvJbUlaQ+ndnD1+pqamamFhYdjVkKSRkuRQVU0tLXeJcklSk2EhSWoyLCRJTYaFJKnJsBgyl0iWNArWZLkPLc8lkiWNCnsWQ+QSyZJGhWExRC6RLGlUOAw1RC6RLGlUGBZD5hLJkkaBw1CSpCbDQpLUZFhIkpoMC0lSk2EhSWoyLCRJTYaFJKnJsJAkNRkWkqQmw0LSyHBJ/+FxuQ9JI8El/YfLnoWkkeCS/sNlWEgaCS7pP1wOQ0kaCS7pP1yGhaSR4ZL+w+MwlCSpybCQJDUZFpKkJsNCktRkWEiSmgwLSVKTYSFJajIsJElNhoUkqcmwkCQ1GRaSpCbDYsz5sBhJq2FgYZHky0neSPKHfWWXJXkiyUvd66V9n+1LciTJi0lu6Sv/UJLnus++kCSDqvO4Of2wmM9+9rPs2LHDwJB03gbZs/gKsHNJ2T3AbFVtB2a79yS5AdgF3Ngd80CSie6YLwF7ge3dz9Jz6gx8WIyk1TKwsKiq/wH8YEnxrcD+bns/cFtf+aNV9U5VvQwcAW5KsgW4uKrmq6qAR/qOUYMPi5G0Wtb6eRZXVdUxgKo6luTKrvwa4Mm+/Y52Zce77aXly0qyl14vhOuuu24Vqz2afFiMpNWyXh5+tNw8RJ2lfFlV9SDwIMDU1NQZ99tIfFiM1oP5+Xm/tIy4tQ6L15Ns6XoVW4A3uvKjwLV9+20FXu3Kty5TLmlEnL7QYnFxkcnJSWZnZw2MEbTWl84eBPZ023uAx/rKdyW5MMn19Cayn+6GrN5OcnN3FdTuvmMkjQAvtBgPA+tZJPnPwAxwRZKjwL8A7gMOJLkDeAW4HaCqDic5ADwPnADurKqT3ak+Qe/KqouAx7sfSSPi9IUWp3sWXmgxmtK7yGj8TE1N1cLCwrCrIQnnLEZJkkNVNbW0fL1McEsaY15oMfpc7kOS1GRYSJKaDAtJUpNhIUlqMiwkSU2GhSSpybCQJDUZFpKkJsNC0or5mN6Nyzu4dU5ctmHjcvXYjc2w0Ir5n8XGttzqsf79bxwOQ2nFXGp6Y/MxvRubPQutmEtNb2w+pndjc4lynRPnLKTx5hLlWhUuNS1tTM5ZSJKaDAtJUpNhIUlqMiykDcw7srVSTnBLG5Q3Wepc2LOQNihvstS5MCykDco7snUuHIaSNijvyNa5MCw0FN4Jvj54k6VWyrDQmnNidfAMY602w0JrzqWuB8sw1iA4wa0158TqYHmVkwbBnoXWnBOrg+VS8hoElyiXRlBrTsI5C50vlyiXxsRK5iS8ykmrzTkLjRTXMnJOQsNhz0IjY6Ne5bN0SMk5CQ2DYaGRsREvuT1TQHqBgNaaYaGRsRG+US/tRZwpIJ2T0FozLDQyzucb9Xq+Kuh03S6//HLeeustLr/8cu66664f6UVshIDUaDAsNFLO5Rv1epvjmJ+f55FHHuG1117jBz/4Ad/61rc4deoUVcWmTZvYtGkTp06d4tSpU+/2Ivbt2+eQk9aFkQmLJDuBzwMTwENVdd+Qq6R17lzmOFarB9J/ntN1uPzyy3n22Wd5+OGHOX78+LLHnTp1CoCJiQmS/EgvwiEnrQcjERZJJoAvAj8HHAWeSXKwqp4fbs20nq10CGclPZAzhcnScDh9ntP/6R8/fpxTp06RhLPdALtp0yYuvPBC7r//ft566y17EVp3RiIsgJuAI1X1HYAkjwK3AoaFzmilcxytHsiZwmRp+Z49e949z+mewumAOFNQTExMcPfdd3PJJZcYEFrXRiUsrgG+1/f+KPA3l+6UZC+wF+C6665bm5ppXVvJEE6rB3KmMFlaDrx7nqU9i02bNnHBBRfw0Y9+9N3zXn311ezevduA0EgYlbDIMmXv+apWVQ8CD0JvbahBV0rjodUDOVOYLC3fvXs3u3fvXnbOwqEljbqRWEgwyTTwL6vqlu79PoCquvdMx7iQoFbTSuYsDAKNgzMtJDgqYXEB8L+BHcCfAM8Af7+qDp/pGMNCks7dSK86W1UnknwS+Ca9S2e/fLagkCStrpEIC4Cq+m3gt4ddD0naiFyiXJLUZFhIkpoMC0lSk2EhSWoaiUtnz0eSN4HvDrseq+QK4PvDrsQqGqf2jFNbYLzaM05tgbVrz1+sqs1LC8c2LMZJkoXlrnseVePUnnFqC4xXe8apLTD89jgMJUlqMiwkSU2GxWh4cNgVWGXj1J5xaguMV3vGqS0w5PY4ZyFJarJnIUlqMiwkSU2GxQhI8k+SVJIr+sr2JTmS5MUktwyzfiuR5NeT/FGS/5Xkvya5pO+zkWrLaUl2dnU+kuSeYdfnXCS5NsnvJXkhyeEkn+rKL0vyRJKXutdLh13XlUoykeTZJN/o3o9yWy5J8vXud+aFJNPDbo9hsc4luRb4OeCVvrIbgF3AjcBO4IEkE8Op4Yo9AXywqv46vWeT7IORbQtdHb8I/B3gBuDvdW0ZFSeAu6vqJ4CbgTu7+t8DzFbVdmC2ez8qPgW80Pd+lNvyeeB3qurHgZ+k166htsewWP/+PfBpfvQxsrcCj1bVO1X1MnAEuGkYlVupqvpvVXWie/sksLXbHrm2dG4CjlTVd6pqEXiUXltGQlUdq6rf77bfpvef0TX02rC/220/cNtQKniOkmwFfh54qK94VNtyMfC3gIcBqmqxqn7IkNtjWKxjST4O/ElVfXvJR9cA3+t7f7QrGxX/CHi82x7Vtoxqvd8jyTbgp4GngKuq6hj0AgW4cohVOxf30/tSdaqvbFTb8peAN4H/2A2rPZTkxxhye0bm4UfjKsl/B65e5qNfBT4D/O3lDlumbOjXQJ+tLVX1WLfPr9IbAvna6cOW2X/obVmBUa33j0jyAeA3gbuq6k+T5Zq1viX5GPBGVR1KMjPk6qyGC4CfAX6pqp5K8nnWwRCaYTFkVfWzy5Un+WvA9cC3u1/grcDvJ7mJ3rfYa/t23wq8OuCqNp2pLacl2QN8DNhRf3aDz7psywqMar3fleR99ILia1X1W13x60m2VNWxJFuAN4ZXwxX7MPDxJB8F3g9cnOSrjGZboPdv62hVPdW9/zq9sBhqexyGWqeq6rmqurKqtlXVNnr/gH6mql4DDgK7klyY5HpgO/D0EKvblGQn8CvAx6vq//Z9NHJt6TwDbE9yfZJJepP0B4dcpxVL7xvIw8ALVfUbfR8dBPZ023uAx9a6bueqqvZV1dbu92QX8LtV9QuMYFsAut/x7yX5q13RDuB5htweexYjqKoOJzlA7x/QCeDOqjo55Gq1/AfgQuCJrqf0ZFX94xFtC1V1IskngW8CE8CXq+rwkKt1Lj4M/CLwXJI/6Mo+A9wHHEhyB70r8G4fTvVWxSi35ZeAr3VfRL4D/EN6X+6H1h6X+5AkNTkMJUlqMiwkSU2GhSSpybCQJDUZFpKkJsNCWmNJfifJD0+vjiqNAsNCWnu/Tu8eB2lkGBbSgCT5G93zO96f5Me650Z8sKpmgbeHXT/pXHgHtzQgVfVMkoPAvwYuAr5aVX845GpJ58WwkAbrX9FbR+r/Ab885LpI581hKGmwLgM+APwFeiuiSiPJsJAG60Hgs/Se3/Fvh1wX6bw5DCUNSJLdwImq+k/dM7v/Z5KPAL8G/DjwgSRHgTuq6pvDrKvU4qqzkqQmh6EkSU2GhSSpybCQJDUZFpKkJsNCktRkWEiSmgwLSVLT/wfVUR7kWheXbAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(x01, y0, 'k.')\n",
    "plt.xlabel(\"x1\")\n",
    "plt.ylabel(\"y\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "a7075054-b563-468a-9876-81eab14f9354",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEGCAYAAACUzrmNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAAAUdklEQVR4nO3df4zc953X8ecr69s04i4icZzE2DF2wYeaUFHIErJUnFbng5jjIDmqCAuBXRHJUpWKiwocCb1CURAuJ0FDxCUo9Ko6/LjIUCCmXHSEheWqsr1kTdqmTi5nt+m1vrhxmlPUIIHHP978MV+H6WZ3v2t3Zmdm9/mQRjPzme93/P5od/zaz/fz/X4mVYUkSSu5atgFSJJGn2EhSWplWEiSWhkWkqRWhoUkqdWmYRcwKDfccEPt3Llz2GVI0lg5duzY96pqy+L2dRsWO3fuZGFhYdhlSNJYSfLbS7V7GEqS1MqwkCS1MiwkSa0MC0lSK8NCktTKsJAktTIsJGmdmJ+f59ChQ8zPz/f9vdftdRaStJHMz8+zZ88eOp0Ok5OTzM7OMj093bf3d2QhSevA3NwcnU6HCxcu0Ol0mJub6+v7GxaStA7MzMwwOTnJxMQEk5OTzMzM9PX9PQwlSevA9PQ0s7OzzM3NMTMz09dDUGBYSNK6MT093feQuMTDUJKkVoaFJKmVYSFJamVYSJJaGRaSpFaGhSSp1cDDIslEkheSfKF5fn2SZ5OcaO6v69n2oSQnk7yS5K6e9tuTvNi89miSDLpuSdL/txYji58DXu55/iAwW1W7gdnmOUluBfYBtwF7gceSTDT7PA4cBHY3t71rULckqTHQsEiyHfhzwGd6mu8GDjePDwP39LQ/VVVnq+pV4CRwR5KtwLVVNV9VBTzZs48kaQ0MemTxCPDzwMWetpuq6jRAc39j074N+E7Pdqeatm3N48XtkqQ1MrCwSPIzwJmqOrbaXZZoqxXal/o3DyZZSLLwxhtvrPKflSS1GeTI4oPAX0jyLeAp4CeT/Cvg9ebQEs39mWb7U8AtPftvB15r2rcv0f4uVfVEVU1V1dSWLVv62RdJ2tAGFhZV9VBVba+qnXQnrv9bVf0V4ChwoNnsAPB08/gosC/J1Ul20Z3Ifq45VPV2kjubs6D29+wjSVoDw1h19lPAkST3Ad8G7gWoquNJjgAvAeeB+6vqQrPPR4DPAdcAzzQ3SdIaSfcEo/VnamqqFhYWhl2GJI2VJMeqampxu1dwS5JaGRaSpFaGhSSplWEhSWplWEiSWhkWkqRWhoUkqZVhIUlqZVhIkloZFpKkVoaFJKmVYSFJamVYSJJaGRaSpFaGhSSplWEhSWplWEiSWhkWkqRWhoUkqZVhIUlqZVhIkloZFpKkVoaFJKmVYSFJamVYSJJaGRaSpFaGhSSplWEhSWplWEiSWhkWkqRWhoUkqZVhIUlqZVhIkloZFpKkVoaFJKmVYSFJamVYSJJaDSwskrwnyXNJvprkeJK/37Rfn+TZJCea++t69nkoyckkryS5q6f99iQvNq89miSDqluS9G6DHFmcBX6yqv4I8AFgb5I7gQeB2araDcw2z0lyK7APuA3YCzyWZKJ5r8eBg8Du5rZ3gHVLkhYZWFhU1/9unv5IcyvgbuBw034YuKd5fDfwVFWdrapXgZPAHUm2AtdW1XxVFfBkzz6SNNLm5+c5dOgQ8/Pzwy7lh7JpkG/ejAyOAX8Q+KWq+o0kN1XVaYCqOp3kxmbzbcCXe3Y/1bSdax4vbpekkTY/P8+ePXvodDpMTk4yOzvL9PT0sMu6IgOd4K6qC1X1AWA73VHCH15h86XmIWqF9ne/QXIwyUKShTfeeOOy65Wkfpqbm6PT6XDhwgU6nQ5zc3PDLumKrcnZUFX1FjBHd67h9ebQEs39mWazU8AtPbttB15r2rcv0b7Uv/NEVU1V1dSWLVv62QVJumwzMzNMTk4yMTHB5OQkMzMzwy7pig3ybKgtSX5v8/ga4KeA3wSOAgeazQ4ATzePjwL7klydZBfdieznmkNWbye5szkLan/PPpI0sqanp5mdneXhhx8e60NQMNg5i63A4Wbe4irgSFV9Ick8cCTJfcC3gXsBqup4kiPAS8B54P6qutC810eAzwHXAM80N0kaedPT02MdEpeke4LR+jM1NVULCwvDLkOSxkqSY1U1tbjdK7glSa0MC0lSK8NCktTKsJAktTIsFlkvl+ZLUj8NdLmPcbOeLs2XpH5yZNFjPV2aL0n9ZFj0WE+X5ktSP3kYqselS/Pn5uaYmZnxEJQkNQyLRdbLpfmS1E8ehpIktTIsJEmtDAtJUivDQpLUyrCQJLUyLCRJrQwLSVIrw0KS1MqwkCS1MiwkSa0MC0lSK8NCktTKsJAktTIsJEmtDAtJUivDQpLUyrCQJLUyLCRJrQwLSVIrw0KS1MqwkCS1MiwkSa0MC0lSq9awSPLRJNetRTGSpNG0mpHFzcDzSY4k2Zskgy5KkjRaWsOiqn4B2A38MvBh4ESSf5jkDwy4NknSiFjVnEVVFfDd5nYeuA74d0l+cYC1SdLAzc/Pc+jQIebn54ddykjb1LZBkr8OHAC+B3wG+FtVdS7JVcAJ4OcHW6IkDcb8/Dx79uyh0+kwOTnJ7Ows09PTwy5rJK1mZHED8Ber6q6q+rdVdQ6gqi4CP7PcTkluSfLfk7yc5HiSn2var0/ybJITzf11Pfs8lORkkleS3NXTfnuSF5vXHnXeRFI/zM3N0el0uHDhAp1Oh7m5uWGXNLJWM2fxd6vqt5d57eUVdj0P/I2qeh9wJ3B/kluBB4HZqtoNzDbPaV7bB9wG7AUeSzLRvNfjwEG6cye7m9cl6YcyMzPD5OQkExMTTE5OMjMzM+ySRlbrYagrVVWngdPN47eTvAxsA+4GZprNDgNzwN9u2p+qqrPAq0lOAnck+RZwbVXNAyR5ErgHeGZQtUvaGKanp5mdnWVubo6ZmRkPQa1gYGHRK8lO4I8CvwHc1AQJVXU6yY3NZtuAL/fsdqppO9c8Xty+1L9zkO4IhB07dvSxB5LWq+npaUNiFQZ+BXeSHwU+DzxQVd9fadMl2mqF9nc3Vj1RVVNVNbVly5bLL1aStKSBhkWSH6EbFP+6qv590/x6kq3N61uBM037KeCWnt23A6817duXaJckrZGBhUVzxtIvAy9X1T/peeko3VNxae6f7mnfl+TqJLvoTmQ/1xyyejvJnc177u/ZR5K0BgY5Z/FB4K8CLyb5StP2d4BPAUeS3Ad8G7gXoKqOJzkCvET3TKr7q+pCs99HgM8B19Cd2HZyW5LWULoXZ68/U1NTtbCwMOwyJGmsJDlWVVOL212iXJLUyrCQJLUyLCRJrQwLSVIrw2KIXBpZ0rhYk+U+9G4ujSxpnDiyGBKXRpY0TgyLIXFpZEnjxMNQQ+LSyJLGiWExRC6NLGlceBhKktTKsJAktTIsJEmtDAtJUivDQpLUyrCQJLUyLCRJrQwLSVIrw0KS1MqwkDQ2XNZ/eFzuQ9JYcFn/4XJkIWksuKz/cBkWksaCy/oPl4ehJI0Fl/UfLsNC0thwWf/h8TCUJKmVYSFJamVYSJJaGRaSpFaGhSSplWEhSWplWEiSWhkWkqRWhoUkqZVhIUlqZVhIkloZFuuYXxQjqV8GFhZJPpvkTJKv97Rdn+TZJCea++t6XnsoyckkryS5q6f99iQvNq89miSDqnk9ufRFMZ/4xCfYs2ePgSHphzLIkcXngL2L2h4EZqtqNzDbPCfJrcA+4LZmn8eSTDT7PA4cBHY3t8XvqSX4RTGS+mlgYVFVvw787qLmu4HDzePDwD097U9V1dmqehU4CdyRZCtwbVXNV1UBT/bsoxX4RTGS+mmtv8/ipqo6DVBVp5Pc2LRvA77cs92ppu1c83hx+5KSHKQ7CmHHjh19LHv8+EUxkvppVL78aKl5iFqhfUlV9QTwBMDU1NSy220UflGMRsn8/Lx/vIyxtQ6L15NsbUYVW4EzTfsp4Jae7bYDrzXt25dolzRGLp1w0el0mJycZHZ21sAYM2t96uxR4EDz+ADwdE/7viRXJ9lFdyL7ueaQ1dtJ7mzOgtrfs4+kMeEJF+NvYCOLJL8CzAA3JDkF/D3gU8CRJPcB3wbuBaiq40mOAC8B54H7q+pC81YfoXtm1TXAM81N0hi5dMLFpZGFJ1yMn3RPMlp/pqamamFhYdhlSGo4ZzEekhyrqqnF7aMywS1pnfOEi/Hmch+SpFaGhSSplWEhSWplWEiSWhkWkqRWhoUkqZVhIUlqZVhIkloZFpIui1/XuzF5Bbcui0s2bGyuHrtxGRZaNf+j0FKrx/o7sDF4GEqr5jLT8ut6Ny5HFlo1l5mWX9e7cblEuS6LcxbS+uYS5eoLl5mWNibnLCRJrQwLSVIrw0KS1MqwkDY4r8jWajjBLW1gXmip1XJkIW1gXmip1TIspA3MK7K1Wh6GkjYwr8jWahkWWnNeBT5avNBSq2FYaE05obp2DGX1k2GhNeUS12vDUFa/OcGtNeWE6trwLCf1myMLrSknVNeGy8mr31yiXBpTbXMSzlnoSrhEubSOrGZOwrOc1E/OWWisuI5Rl3MSWmuOLDQ2NuoZPksdTnJOQmvNsNDY2Iin3S4XkJ4ooLVmWGhsbJS/pntHEisFpHMSWkuGhcbGlfw1PcpnBPXWBt2R0+bNm3nggQfeCcRHHnlkQwSkRp9hobFyOX9Nj9Icx6Vg2Lx5My+88AIvvfQSX/rSl6gqNm3aRBLOnz9PEi5evMjFixfpdDq8+eabHm7SSBibsEiyF/inwATwmar61JBL0ohb7RxHv0Yfi9/n0vO33nqLT3/605w/f56lrms6d+4cAFXFVVddxcTEBEneGUl4uEmjYCzCIskE8EvAnwZOAc8nOVpVLw23Mo2y1cxxrGb0sVyYLD6M1Ps+jzzyCA888ABnz57l4sWLK9Z51VVXsWnTJs6fP//Ovm+++aYjCY2UsQgL4A7gZFV9EyDJU8DdgGGhZa1mjqNt9LFcmCxuP3DgwA+8z+c//3k6nU5rUExMTPDYY4/x/ve/30NNGmnjEhbbgO/0PD8F/InFGyU5CBwE2LFjx9pUppHWdginbfSxXJgsbgd+4H0+9KEP8cUvfvGdkUUSJiYm+NjHPsb3v/99vvvd73LzzTezf//+Hzi7SRpV4xIWWaLtXQd/q+oJ4Anorg016KI0/tpGH8uFyeL2/fv3s3///h94n0ujhc2bN3tYSWNvLBYSTDINfLKq7mqePwRQVYeW28eFBNUvq5mzMAS0Xiy3kOC4hMUm4LeAPcDvAM8Df7mqji+3j2EhSZdvrFedrarzST4K/BrdU2c/u1JQSJL6ayzCAqCqfhX41WHXIUkbkUuUS5JaGRaSpFaGhSSplWEhSWo1FqfOXokkbwOvDLuOPrgB+N6wi+gT+zKa1ktf1ks/YLh9+f1VtWVx49icDXUFXlnqXOFxk2RhPfQD7MuoWi99WS/9gNHsi4ehJEmtDAtJUqv1HBZPDLuAPlkv/QD7MqrWS1/WSz9gBPuybie4JUn9s55HFpKkPjEsJEmtxj4sktyb5HiSi0mmetp3Jvk/Sb7S3P55z2u3J3kxyckkjyZZ6suV1txyfWlee6ip95Ukd/W0j2RfeiX5ZJLf6flZ/HTPa0v2a1Ql2dvUejLJg8Ou53Il+Vbz+/KVJAtN2/VJnk1yorm/bth1LiXJZ5OcSfL1nrZlax/l361l+jLan5OqGusb8D7gDwFzwFRP+07g68vs8xwwTfcb+J4B/uyw+9HSl1uBrwJXA7uAbwATo9yXRf36JPA3l2hftl+jeKO7PP43gPcCk03ttw67rsvsw7eAGxa1/SLwYPP4QeAfDbvOZWr/CeCP9X6ul6t91H+3lunLSH9Oxn5kUVUvV9Wqr9ROshW4tqrmq/uTeBK4Z1D1XY4V+nI38FRVna2qV4GTwB2j3JdVWrJfQ65pJXcAJ6vqm1XVAZ6i24dxdzdwuHl8mBH9HaqqXwd+d1HzcrWP9O/WMn1Zzkj0ZezDosWuJC8k+R9J/lTTtg041bPNqaZtlG0DvtPz/FLN49SXjyb5WjP8vnSoYLl+japxq3cpBfyXJMeSHGzabqqq0wDN/Y1Dq+7yLVf7uP6sRvZzMhbLfST5r8DNS7z08ap6epndTgM7qurNJLcD/zHJbXQP1yy2ZucPX2Fflqt5qH3ptVK/gMeBh+nW9jDwj4G/xgjVv0rjVu9SPlhVryW5EXg2yW8Ou6ABGcef1Uh/TsYiLKrqp65gn7PA2ebxsSTfAH6cbipv79l0O/BaP+pcZV2X3Re6Nd/S8/xSzUPtS6/V9ivJvwC+0Dxdrl+jatzqfZeqeq25P5PkP9A9nPF6kq1Vdbo5tHlmqEVenuVqH7ufVVW9funxKH5O1u1hqCRbkkw0j98L7Aa+2QxV305yZ3Pm0H5gub/oR8VRYF+Sq5PsotuX58alL82H+JKfBS6dAbJkv9a6vsvwPLA7ya4kk8A+un0YC0l+T5Ifu/QY+DN0fxZHgQPNZgcYwd+hFSxX+7j9bo3+52TYZwX04ayCn6WbvGeB14Ffa9o/BBynexbB/wL+fM8+U3R/EN8A/hnNlezDvi3Xl+a1jzf1vkLPGU+j2pdF/fqXwIvA1+j+4m9t69eo3oCfBn6rqfnjw67nMmt/b/N5+Grz2fh4074ZmAVONPfXD7vWZer/FbqHl881n5P7Vqp9lH+3lunLSH9OXO5DktRq3R6GkiT1j2EhSWplWEiSWhkWkqRWhoUkqZVhIa2hJB9IMt+sLvy1JH9p2DVJq+Gps9IaSvLjQFXViSS/DzgGvK+q3hpuZdLKHFlIA5Lkjzejh/c0V08fByar6gS8s/TGGWDLUAuVVsGRhTRASf4B8B7gGuBUVR3qee0Oustq31ZVF4dUorQqhoU0QM0aUs8D/xf4k1V1oWnfSvdLrg5U1ZeHV6G0Oh6GkgbreuBHgR+jO8IgybXAfwZ+waDQuHBkIQ1QkqN0v1FvF7AV+Bjdr7/9T1X1yBBLky7LWHyfhTSOkuwHzlfVv2mWy/+fdJc1/wlgc5IPN5t+uKq+MpwqpdVxZCFJauWchSSplWEhSWplWEiSWhkWkqRWhoUkqZVhIUlqZVhIklr9P45O8nOo+wAvAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(x02, y0, 'k.')\n",
    "plt.xlabel(\"x2\")\n",
    "plt.ylabel(\"y\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "d6e5e822-6c0a-4b6e-b7fe-ce59559da296",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEGCAYAAACZ0MnKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAAAUZElEQVR4nO3df6zd9X3f8efL1zOJkqAmAxKGYaYVawdZy5obL1akyR1O8boqQNVUrrSABotLAm0ndVox2TT2A5G2WrN0CdlMkgXyowwlpVgpDSXeUCL1JnBp8wNDWBzIwMUNTtOqTGvs2Lz3x/ne5mAf21+Te873fM95PqSre87n+z3X74/E5XU/n+/n+/mmqpAkqY01XRcgSeoPQ0OS1JqhIUlqzdCQJLVmaEiSWlvbdQHjdsYZZ9SGDRu6LkOSeuXhhx/+VlWdeXT7zIfGhg0bWF5e7roMSeqVJP9nVLvTU5Kk1gwNSVJrhoYkqTVDQ5LUmqEhSWrN0JAktWZoSNKMWVpa4pZbbmFpaWnVf/bM36chSfNkaWmJSy65hEOHDrFu3Tp2797Npk2bVu3nO9KQpBnywAMPcOjQIY4cOcKhQ4d44IEHVvXnGxqSNEM2b97MunXrWFhYYN26dWzevHlVf77TU5I0QzZt2sTu3bt54IEH2Lx586pOTYGhIUkzZ9OmTaseFiucnpIktWZoSJJaMzQkSa11GhpJPpTk2SSPDLXdlORPknyx+fqpoWM7kuxN8niSS7upWpLmV9cjjQ8DW0e0v7uqLm6+7gVIciGwDbio+cytSRYmVqkkqdvQqKrPAt9uefplwJ1VdbCqngT2AhvHVpwk6RhdjzSO5/okX26mr17ZtJ0DPD10zr6m7RhJtidZTrJ84MCBcdcqSXNjGkPj/cAPARcD+4H/1LRnxLk16gdU1c6qWqyqxTPPPOa56JKkF2nqQqOqvllVR6rqeeA2vjcFtQ84d+jU9cAzk65PkubZ1IVGkrOH3l4BrKys2gVsS3JakvOBC4AHJ12fJM2zTrcRSfLbwGbgjCT7gH8LbE5yMYOpp28AvwBQVXuS3AU8ChwGrquqIx2ULUlzK1UjLwvMjMXFxVpeXu66DEnqlSQPV9Xi0e1TNz0lSbNonE/TmyR3uZWkMRv30/QmyZGGJI3ZuJ+mN0mGhiSN2bifpjdJTk9J0piN+2l6k2RoSNIEjPNpepPk9JQkqTVDQ5LUmqEhSWrN0JAktWZoSJJaMzQkSa0ZGpKk1gwNSVJrhoYkqTVDQ5LUmqEhae7NyrMuJsG9pyTNtVl61sUkONKQNNdm6VkXk2BoSJprs/Ssi0lwekrSXJulZ11MgqEhae7NyrMuJsHpKUlSa4aGJKm1TkMjyYeSPJvkkaG2VyW5P8nXmu+vHDq2I8neJI8nubSbqiVpfnU90vgwsPWothuA3VV1AbC7eU+SC4FtwEXNZ25NsjC5UiVJnYZGVX0W+PZRzZcBtzevbwcuH2q/s6oOVtWTwF5g4yTqlCQNdD3SGOXVVbUfoPl+VtN+DvD00Hn7mrZjJNmeZDnJ8oEDB8ZarCTNk2kMjePJiLYadWJV7ayqxapaPPPMM8dcliTNj2kMjW8mORug+f5s074POHfovPXAMxOuTdIUcIPB7kzjzX27gKuAdzXf7xlq/3iS3wT+FnAB8GAnFUrqjBsMdqvrJbe/DSwBP5xkX5JrGITFm5J8DXhT856q2gPcBTwKfBq4rqqOdFO5pK64wWC3Oh1pVNXPH+fQJcc5/2bg5vFVJGnarWwwuDLScIPByZrG6SlJOi43GOyWoSGpd9xgsDvTuHpKkjSlDA1JUmuGhiSpNUNDktSaoSFJas3QkDQxbv/Rfy65lTQRbv8xGxxpSJoIt/+YDYaGpIlY2f5jYWHB7T96zOkpSRPh9h+zwdCQNDFu/9F/Tk9JklozNCRJrRkakqTWDA1JUmuGhqRT5p3d88vVU5JOiXd2zzdHGpJOiXd2zzdDQ9Ip8c7u+eb0lKRT4p3d883QkHTKvLN7fjk9JUlqbWpHGkm+ATwHHAEOV9ViklcB/wPYAHwD+Lmq+vOuapRmxdLSktNNamVqQ6PxE1X1raH3NwC7q+pdSW5o3v9qN6VJs8EltDoVfZueugy4vXl9O3B5d6VIs8EltDoV0xwaBfxBkoeTbG/aXl1V+wGa72eN+mCS7UmWkywfOHBgQuVK/eQSWp2KaZ6eemNVPZPkLOD+JF9t+8Gq2gnsBFhcXKxxFSjNApfQ6lRMbWhU1TPN92eT3A1sBL6Z5Oyq2p/kbODZTouUZoRLaNXWVE5PJXlZklesvAZ+EngE2AVc1Zx2FXBPNxVK0nya1pHGq4G7k8Cgxo9X1aeTPATcleQa4CngLR3WKE09l9JqtU1laFTVE8CPjWj/M+CSyVck9Y9LaTUOUzk9Jen751JajYOhIc0ol9JqHKZyekrS98+ltBoHQ0PqsZNd6HYprVaboSH1lBe61QWvaUg95YVudcHQkHpiaWmJW265haWlJcAL3eqG01NSDxxvKsoL3Zo0Q0PqgVFTUSsXuQ0LTZLTU9IUcipK08qRhjRlnIrSNDM0pCmxcs/FU0895VSUppahIU2B4dHFwsICa9cOfjWditK0MTSkji0tLXHTTTdx8OBBnn/+eQDe9ra3cd555zkVpalzwtBIcjpwZlV9/aj2H62qL4+1MmkO7Ny5k+uvv57Dhw9TVaxZs4Z169Zx5ZVXGhaaSsddPZXk54CvAp9MsifJ64cOf3jchUmzbGlpiSuuuIJrr72W7373u1QVSdiyZYvbgWiqnWikcSPwuuZ53BuBjyS5sap+B8hkypNmz86dO3nHO97BkSNHXtC+du1abrrpJgNDU+1EobFQVfsBqurBJD8BfCrJeqAmUp00Q5aWlrjjjju47bbbjgmMNWvW8N73vtfA0NQ7UWg8l+SHVq5nNCOOzcDvAheNvzRpdqysjvrOd75D1Qv/5lpYWODWW29l+/btHVUntXeiO8LfDqxJcuFKQ1U9B2wF/vm4C5NmxfDqqOHAWFhY4PLLL+dzn/ucgaHeOO5Io6q+BJDkkSQfAX4deEnzfRH4yEQqlHpsZYSxspx2zZo1rF27lquvvtoVUuqlNvdp/APg14A/BF4BfAx44ziLkvru6Lu7VwJjy5YtXuxWr7UJje8CfwW8lMFI48mqen6sVUk9dqK7uw0M9V2b0HgIuAd4PfA3gf+W5Ger6mfHWpnUU8PbmIN3d2u2tAmNa6pquXn9p8BlSd46xppOKMlW4D3AAvCBqnpXV7VIK1amozZv3vzX25iv7FLrtQvNkpOGxlBgDLd1chE8yQLwPuBNwD7goSS7qurRLuqRYPRW5m5jrlnVtw0LNwJ7q+oJgCR3ApcBhoY6M+qpejt27DAsNJP69uS+c4Cnh97va9peIMn2JMtJlg8cODCx4jSffKqe5knfRhqj9rw6ZkuTqtoJ7ARYXFx0yxOtmuFrFysjCZ+qp3nSt9DYB5w79H498ExHtWjOHO8xrIBP1dPc6Nv01EPABUnOT7IO2Abs6rgmzYlR1y6kedOrkUZVHU5yPXAfgyW3H6qqPR2XpTlx9FJar11oHvUqNACq6l7g3q7r0Pzx2oXUw9CQxmnUhe5hXrvQvDM0pMaJLnRLGujbhXBpbLzQLZ2coSE1vElPOjmnp6SGF7qlkzM0pCFe6JZOzOkpzYWlpSVuueUWlpaWui5F6jVHGpp5roqSVo8jDc08V0VJq8fQ0MxzVZS0epye0sxzVZS0egwNzQVXRUmrw+kpSVJrhoYkqTVDQ73kfRdSN7ymod7xvgupO4401DvedyF1x9BQ73jfhdQdp6fUO953IXXH0FAved+F1A2npyRJrRkakqTWDA1JUmuGhjrlTXpSv0zdhfAkNwFvAw40TTdW1b3NsR3ANcAR4Jeq6r5OitSq8CY9qX+mdaTx7qq6uPlaCYwLgW3ARcBW4NYkC10Wqe+PN+lJ/TOtoTHKZcCdVXWwqp4E9gIbO65J3wdv0pP6Z+qmpxrXJ7kSWAZ+par+HDgH+PzQOfuatmMk2Q5sBzjvvPPGXKpeLG/Sk/onVTX5fzT5DPCaEYfeySAYvgUU8B+As6vq6iTvA5aq6qPNz/ggcG9VffJE/9bi4mItLy+vav2SNOuSPFxVi0e3dzLSqKotbc5LchvwqebtPuDcocPrgWdWuTRJ0glM3TWNJGcPvb0CeKR5vQvYluS0JOcDFwAPTro+SZpn03hN49eTXMxgeuobwC8AVNWeJHcBjwKHgeuq6khXRUrSPJq60Kiqt57g2M3AzRMsR5I0ZOqmpyRJ08vQ0Ivi9h/SfJq66SlNP7f/kOaXIw2dMrf/kOaXoaFT5vYf0vxyekqnzO0/pPllaOhF8Rnd0nxyekqS1JqhIUlqzdCQJLVmaEiSWjM0JEmtGRqSpNYMDUlSa4aGJKk1Q2NOuCutpNXgHeFzwF1pJa0WRxpzwF1pJa0WQ2MOuCutpNXi9NQccFdaSavF0JgT7koraTU4PSVJas3QkCS1ZmhIklrrJDSSvCXJniTPJ1k86tiOJHuTPJ7k0qH21yX5SnPst5Jk8pVL0nzraqTxCPAzwGeHG5NcCGwDLgK2ArcmWWgOvx/YDlzQfG2dWLWSJKCj0Kiqx6rq8RGHLgPurKqDVfUksBfYmORs4PSqWqqqAu4ALp9cxZIkmL5rGucATw+939e0ndO8Prp9pCTbkywnWT5w4MBYCpWkeTS2+zSSfAZ4zYhD76yqe473sRFtdYL2kapqJ7ATYHFx8bjnSZJOzdhCo6q2vIiP7QPOHXq/HnimaV8/ol2SNEHTNj21C9iW5LQk5zO44P1gVe0HnkvyhmbV1JXA8UYrveTW5ZL6oJNtRJJcAfwX4Ezg95J8saourao9Se4CHgUOA9dV1ZHmY28HPgy8FPj95msmuHW5pL7oJDSq6m7g7uMcuxm4eUT7MvDaMZfWiVFblxsakqbRtE1PzSW3LpfUF+5yOwXculxSXxgaU8KtyyX1gdNTkqTWDA1JUmuGhiSpNUNDktSaoSFJas3QkCS1ZmhIklozNCRJrRkakqTWDA1JUmuGhiSpNUNDktSaoXEcPklPko7lLrcj+CQ9SRrNkcYIo56kJ0kyNEbySXqSNJrTUyP4JD1JGs3QOA6fpCdJx3J6SpLUmqEhSWrN0JAktdZJaCR5S5I9SZ5PsjjUviHJXyX5YvP1X4eOvS7JV5LsTfJbSdJF7ZI0z7oaaTwC/Azw2RHHvl5VFzdf1w61vx/YDlzQfG0df5mSpGGdhEZVPVZVj7c9P8nZwOlVtVRVBdwBXD6u+iRJo03jktvzk/wx8JfAv66qzwHnAPuGztnXtI2UZDuDUQnA/03SOqCm3BnAt7ouYhXNUn9mqS8wW/2Zpb7A5Przt0c1ji00knwGeM2IQ++sqnuO87H9wHlV9WdJXgf8bpKLgFHXL+p4/3ZV7QR2nmrN0y7JclUtnvzMfpil/sxSX2C2+jNLfYHu+zO20KiqLS/iMweBg83rh5N8Hfg7DEYW64dOXQ88sxp1SpLam6olt0nOTLLQvP5BBhe8n6iq/cBzSd7QrJq6EjjeaEWSNCZdLbm9Isk+YBPwe0nuaw79Q+DLSb4EfAK4tqq+3Rx7O/ABYC/wdeD3J1z2NJi1KbdZ6s8s9QVmqz+z1BfouD8ZLEaSJOnkpmp6SpI03QwNSVJrhkaPJPmXSSrJGUNtO5qtVR5PcmmX9bWR5DeSfDXJl5PcneQHho71qi8rkmxtat6b5Iau6zkVSc5N8r+SPNZs7fPLTfurktyf5GvN91d2XWtbSRaS/HGSTzXv+9yXH0jyieZ35rEkm7ruj6HRE0nOBd4EPDXUdiGwDbiIwbYqt66sPpti9wOvraofBf43sAN62xeaGt8H/GPgQuDnm770xWHgV6rq7wJvAK5r6r8B2F1VFwC7m/d98cvAY0Pv+9yX9wCfrqofAX6MQb867Y+h0R/vBv4VL7yp8TLgzqo6WFVPMlhZtrGL4tqqqj+oqsPN28/zvftveteXxkZgb1U9UVWHgDsZ9KUXqmp/Vf1R8/o5Bv9TOodBH25vTrudnmzbk2Q98E8YrLRc0de+nM5gRekHAarqUFX9BR33x9DogSRvBv6kqr501KFzgKeH3p9we5UpdDXfWzrd1770te5jJNkA/H3gC8Crm/ujaL6f1WFpp+I/M/jj6vmhtr725QeBA8B/b6bbPpDkZXTcn2nce2ounWjbFeBG4CdHfWxEW+drqNtsIZPknQymRj628rER53felxb6WvcLJHk58EngX1TVX/bxyQNJfhp4ttlNYnPH5ayGtcCPA79YVV9I8h6mYGrN0JgSx9t2JcnfA84HvtT8Iq8H/ijJRgZ/1Z47dPpUbK9ysi1kklwF/DRwSX3vRqGp7EsLfa37ryX5GwwC42NV9TtN8zeTnF1V+5tdpp/trsLW3gi8OclPAS8BTk/yUfrZFxj8t7Wvqr7QvP8Eg9DotD9OT025qvpKVZ1VVRuqagOD/5B+vKr+FNgFbEtyWpLzGWy78mCH5Z5Ukq3ArwJvrqr/N3Sod31pPARckOT8JOsYXMzf1XFNrTXb8nwQeKyqfnPo0C7gqub1VfRg256q2lFV65vfk23A/6yqf0oP+wLQ/I4/neSHm6ZLgEfpuD+ONHqsqvYkuYvBf0iHgeuq6kjHZZ3Me4HTgPubkdPnq+ranvaFqjqc5HrgPmAB+FBV7em4rFPxRuCtwFeSfLFpuxF4F3BXkmsYrNh7SzflrYo+9+UXgY81f5A8AfwzBn/sd9YftxGRJLXm9JQkqTVDQ5LUmqEhSWrN0JAktWZoSJJaMzSkjiT5dJK/WNmNVeoDQ0Pqzm8wuEdC6g1DQxqzJK9vnh/ykiQva55b8dqq2g0813V90qnwjnBpzKrqoSS7gP8IvBT4aFU90nFZ0otiaEiT8e8Z7FP1HeCXOq5FetGcnpIm41XAy4FXMNiBVeolQ0OajJ3Av2Hw/JBf67gW6UVzekoasyRXAoer6uPNM8X/MMk/Av4d8CPAy5PsA66pqvu6rFU6GXe5lSS15vSUJKk1Q0OS1JqhIUlqzdCQJLVmaEiSWjM0JEmtGRqSpNb+P6rTW4WIuvLHAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(x01, x02, 'k.')\n",
    "plt.xlabel(\"x1\")\n",
    "plt.ylabel(\"x2\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e3a55090-b99b-4980-92ed-a2fe297db360",
   "metadata": {},
   "source": [
    "So, there is likely a path x1 -- x2 or a counfounder of the two or a confounder of y or both."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "717a547c-46c1-40b0-8496-53cd929db54a",
   "metadata": {},
   "source": [
    "# do(x1 = a, x2 = b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "35d3b0aa-41ff-4159-828e-1e13e87a8c17",
   "metadata": {},
   "outputs": [],
   "source": [
    "a = 35\n",
    "b = 10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "05a3e630-329f-4705-8ab2-4342094d2096",
   "metadata": {},
   "outputs": [],
   "source": [
    "x11 = a * np.ones(N)\n",
    "x12 = b * np.ones(N)\n",
    "y1 = x11**2 + x12 + C"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "33c2adb9-d7e0-43b0-a4ef-24a15743f6b6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEGCAYAAACUzrmNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAAAa2UlEQVR4nO3df5RU9Znn8fcnYgSVXyJGYuvguIgKgj86EXQE1GDYSIhukpN43IEsRoKrMya7UdPDEp0k6rrImHHdNRJhFUEmxlmOhBkU1kgwG1DbyE81BOPE6QUiggiRYOjw7B91G4ui6NtIVd2uqs/rnDrd96lbt5+vAk9/f9zvVURgZmbWno9knYCZmXV+LhZmZpbKxcLMzFK5WJiZWSoXCzMzS9Ul6wTK5fjjj4/+/ftnnYaZWVV56aWX3o6IvoXxmi0W/fv3p7m5Oes0zMyqiqTfFot7GMrMzFK5WJiZWSoXCzMzS1WzcxbF7Nmzh5aWFnbv3p11KlZhXbt2paGhgSOPPDLrVMyqUl0Vi5aWFrp3707//v2RlHU6ViERwdatW2lpaeHUU0/NOh2zqlRXw1C7d++mT58+LhR1RhJ9+vRxj9LsMNRVsQBcKOqU/79bXVi+HO66K/e1xOpqGMrMrGYtXw6XXQZ//CN89KPwzDMwfHjJLl93PQszs5q0dGmuUPzpT7mvS5eW9PIuFhm6/fbbueeee1LP+/rXv86yZcvaPWfKlCmcfPLJHHvssaVKr2RuvvlmzjjjDIYMGcJVV13F9u3bD/uao0aNKnqH/sKFC7ntttsO+/pmVWfUqFyP4ogjcl9HjSrp5V0s0pRxDLAjtm3bxooVKxgxYkS75332s5/lhRdeqFBWh2b06NGsXbuW1atXc/rpp3PXXXeV7WddccUVLFiwgF27dpXtZ5h1SsOH54aevvvdkg9BgYtF+9rGAKdOzX0tQcG44447GDhwIJ/61Kf41a9+tS++cuVKhg0btu+373feeQeAJ554gjFjxgDw7rvvMnDgwH2fu/rqq/nhD38IwLBhw+jXr99h55fv+uuvp7GxkUGDBnXot/XPfe5zzJ49G4AHH3yQa665BoDLL7+cLl267MuzpaWlJPnNmTOHCy+8kMGDB+8rlJIYNWoUCxcuLMnPMKsqw4dDU1PJCwWQW4NejhcwC3gLWJsXmwa8BqwG5gO98t4bAiwH1gFrgK5J/PzkeANwH6CO/Pzzzz8/Cr3yyisHxNp1550RRxwRAbmvd955aJ8v0NzcHIMHD4733nsv3n333TjttNNi2rRpERFx9tlnx9KlSyMiYurUqXHTTTdFRMT48eNjwYIF+66xePHiGDZsWMybNy8+/elPH/AzjjnmmIP+/J/+9KcxdOjQA17Dhw8vev7WrVsjIqK1tTVGjhwZq1atard9mzdvjtNOOy2WLVsWAwYM2Pf5fGPHjo1HH330gPiOHTuK5jZ06NBYt27dAeePHDkyvvrVr0ZExM9+9rMYNGjQvvfmzJkTN9544wGfOeT//2Z1CGiOIv+mlnM11MPA/cDsvNgSoCkiWiXdDTQBt0rqAswB/jIiVknqA+xJPvMAMAlYAfwzMAZYVMa8P9A2Bti2uuAwxwCfe+45rrrqKo4++mgAxo0bB+R6DNu3b2fkyJEATJgwgS9+8YsAbNq0ib59P9gtePTo0fz4xz/mhhtuYNWqVYf08y+55BJWrlzZ4fMff/xxZsyYQWtrK5s2beKVV15hyJAhBz3/Yx/7GN/5zne45JJLmD9/Pscdd9x+799xxx106dJlX48jX/fu3Q8pN8j1rABGjBjBjh072L59O7169eKEE05g48aNh3QtM2tf2YpFRCyT1L8gtjjvcAXwheT7y4HVEbEqOW8rgKR+QI+IWJ4czwaupFLFom0McOnSXKEoQdfuUNf7d+vWbb+byfbu3curr75Kt27d2LZtGw0NDR2+1rPPPss3vvGNA+JHH300v/jFL/aLvfHGG9xzzz28+OKL9O7dm6985SsduqltzZo19OnT54B/rB955BEWLlzIM888U/S/wc6dO7n44ouLXvOxxx7jrLPOOiBeeJ224927d9OtW7fUXM2s47Kcs5jIB//onw6EpKcl/VLSLUn8JCB/gLsliVVOCccAR4wYwfz58/nDH/7Azp07+clPfgJAz5496d27N8899xwAjz766L5explnnsmGDRv2XePee+/lzDPPZN68eUycOJE9e/Yc+IMOoq1nUfgqLBQAO3bs4JhjjqFnz5787ne/Y9GiD+pzU1MT8+fPP+AzL7zwAosWLeLll1/mnnvu4Y033gDgqaee4u6772bBggX7elWF2noWxV7FCgXAj370IwB+/vOf07NnT3r27AnA+vXrGTx4cIf/u5hZukyKhaQpQCswNwl1Af4CuCb5epWky4Biv4ZHO9edJKlZUvOWLVtKnPXhO++88/jSl77EOeecw+c///n9fpN+5JFHuPnmmxkyZAgrV67k29/+NpBb3bM0WS+9fv16HnroIaZPn87FF1/MiBEj+N73vgfALbfcQkNDA7t27aKhoYHbb7/9sHIdOnQo5557LoMGDWLixIlcdNFF+95bs2YNJ5544n7nv//++1x33XXMmjWLj3/840yfPp2JEycSEdx4443s3LmT0aNHc8455zB58uTDyq1N7969ufDCC5k8eTIzZ87cF3/22We54oorSvIzzCxRbCKjVC+gP3kT3ElsArmJ7KPzYl8GHs47ngrcDPQDXsuLXw082JGfXZIJ7k7ioosuinfeeSfrNPa5/PLLs07hoDZv3hyXXnpp0feq9f+/WSVxkAnuivYsJI0BbgXGRUT+QvingSGSjk4mu0cCr0TEJmCnpGHKDUiPB56sZM6dwfTp03nzzTezTmOfp59+OusUDurNN99k+vTpWadhVnPKNsEtaR4wCjheUgtwG7nVT0cBS5LJyBURMTki3pH0d8CL5IaZ/jki/im51PXkVlZ1IzfHcViT2xFRdZvKXXDBBVmnUDU+8YlPFI3nfmEysw+rnKuhri4Snlkk1nb+HHLLZwvjzUBJZiu7du3K1q1bvU15nYnkeRZdu3bNOhWzqlVXu842NDTQ0tJCZ5z8tvJqe1KemX04dVUsjjzySD8pzczsQ/DeUGZmlsrFwszMUrlYmJmVU8aPOSiVupqzMDOrqDI/6rSS3LMwMyuXMj/qtJJcLMzMyqXMjzqtJA9DmZmVSxkec5AVFwszs3IaPryqi0QbD0OZmVkqFwszM0vlYmFmZqlcLMzMLJWLhZmZpXKxMDOzVC4WZmaWysXCzMxSuViYmVkqFwszM0vlYmFmZqlcLMysftXIg4kqwRsJmll9qqEHE1WCexZmVp9q6MFEleBiYWb1qYYeTFQJHoYys/pUQw8mqgQXCzOrXzXyYKJK8DCUmZmlKluxkDRL0luS1ubFpkl6TdJqSfMl9Uri/SX9QdLK5PWDvM+cL2mNpA2S7pOkcuVsZmbFlbNn8TAwpiC2BBgcEUOA9UBT3nuvR8Q5yWtyXvwBYBIwIHkVXtPMzMqsbMUiIpYB2wpiiyOiNTlcATS0dw1J/YAeEbE8IgKYDVxZhnTNzKwdWc5ZTAQW5R2fKullST+TdHESOwloyTunJYkVJWmSpGZJzVu2bCl9xmZmdSqTYiFpCtAKzE1Cm4BTIuJc4D8Bj0nqARSbn4iDXTciZkREY0Q09u3bt9Rpm5nVrYovnZU0ARgLXJYMLRER7wPvJ9+/JOl14HRyPYn8oaoGYGNlMzYzs4r2LCSNAW4FxkXErrx4X0lHJN//ObmJ7N9ExCZgp6RhySqo8cCTlczZzMzKu3R2HrAcGCipRdK1wP1Ad2BJwRLZEcBqSauAJ4DJEdE2OX498BCwAXid/ec5zKyeeJfYzCgZCao5jY2N0dzcnHUaZlYq3iW2IiS9FBGNhXHfwW1m1cG7xGbKxcLMqoN3ic2UNxI0s+rgXWIz5WJhZtXDu8RmxsNQZmaWysXCzMxSuViYmVkqFwszM0vlYmFmZqlcLMzMLJWLhZmVn/d0qnq+z8LMyst7OtUE9yzMrLy8p1NNcLEws/Lynk41wcNQZlZe3tOpJrhYmFn5eU+nqudhKDMzS+ViYWZmqVwszMwslYuFmXWcb66rW57gNrOO8c11dc09CzPrGN9cV9dcLMysY3xzXV3zMJSZdYxvrqtrLhZm1nG+ua5ueRjKzMxSuViYmVkqFwuzeub7JqyDylYsJM2S9JaktXmxaZJek7Ra0nxJvQo+c4qk30v6Zl7sfElrJG2QdJ8klStns7rSdt/E1Km5ry4Y1o5y9iweBsYUxJYAgyNiCLAeaCp4/15gUUHsAWASMCB5FV7TzD4M3zdhh6BsxSIilgHbCmKLI6I1OVwBNLS9J+lK4DfAurxYP6BHRCyPiABmA1eWK2ezuuL7JuwQZLl0diLwIwBJxwC3AqOBb+adcxLQknfcksSKkjSJXC+EU045pcTpmtUY3zdhhyCTYiFpCtAKzE1CfwvcGxG/L5iSKDY/EQe7bkTMAGYANDY2HvQ8M0v4vgnroIoXC0kTgLHAZcnQEsAFwBck/TegF7BX0m7gH8kbqkq+31jBdM2q0/Ll7jFYSVW0WEgaQ264aWRE7GqLR8TFeefcDvw+Iu5PjndKGgY8D4wH/nslczarOt4d1sqgnEtn5wHLgYGSWiRdC9wPdAeWSFop6QcduNT1wEPABuB1DlwtZWb5vMrJyqBsPYuIuLpIeGYHPnd7wXEzMLhEaZnVvrZVTm09C69yshLwRoJmtcarnKwMXCzMqlHaBLZXOVmJuViYVRtPYFsGvJGgWbXxBLZlwMXCrNp4mw7LgIehzDq7wvkJT2BbBlwszDqzg81PeALbKszDUGadSeHDiDw/YZ2EexZmnUWxXoRvsLNOwsXCLGttcxJvvnlgL6KpyfMT1im4WJhVyvLlMHs2bN68f3zRImhtza1u6pL8lczvRXh+wjqB1GIh6UZgbkS8U4F8zGpHW4+hTx94+WWYORP27Gn/M9ddB6ec4l6EdTod6VmcCLwo6ZfALODpvOdQmFkxbfMP778Pe/eCBO39tZFyvYnx410krFNKXQ0VEf8FGEBux9ivAL+WdKek08qcm1n1aVvNNHt2bt5h795cvL1CceSR8LWvedsO69Q6NGcRESFpM7CZ3ONQewNPSFoSEbeUM0GzqpG/mqlt/iEiVzA+8pHc8Wc+s/9nTjzRvQmrCh2Zs/hrYALwNrmHEN0cEXskfQT4NeBiYfUr/+7q/Hsi4IP5hz59YOtWz0NYVetIz+J44N9FxG/zgxGxV9LY8qRlVgUK74v4/vf3vyfCPQarIanFIiK+3c57r5Y2HbMqUnh39datvifCapbvszD7sIrdXe17IqxGuViYdUSxJ9N591erIy4WZmnaezKdexJWJ7zrrFka7/xq5mJhlspPpjPzMJRZKs9NmLlYmAHFJ7DzeW7C6pyLhVl7E9hmBnjOwswT2GYdULZiIWmWpLckrc2LTZP0mqTVkuZL6pXEPylpZfJaJemqvM+cL2mNpA2S7pOkcuVsdcoT2GapytmzeBgYUxBbAgyOiCHAeqApia8FGiPinOQzD0pqGyJ7AJhEbpv0AUWuaXZ42iawv/tdD0GZHUTZ5iwiYpmk/gWxxXmHK4AvJPFdefGuQABI6gf0iIjlyfFs4EpgUbnytjrlCWyzdmU5ZzGRvH/0JV0gaR2wBpgcEa3ASUBL3mdaklhRkiZJapbUvGXLljKlbVWl7WFEy5dnnYlZVctkNZSkKeQeojS3LRYRzwODJJ0JPCJpEVBsfuKgjxyLiBnADIDGxkY/+rXeeZWTWclUvGchaQIwFrim2LO8k23P3wMGk+tJNOS93QBsrESeVgO8ysmsZCpaLCSNAW4FxuXPU0g6tW1CW9KfAQOBf4mITcBOScOSVVDjgScrmbNVMa9yMiuZsg1DSZoHjAKOl9QC3EZu9dNRwJJkBeyKiJgM/AXwLUl7gL3Af4yIt5NLXU9uZVU3cnMcnty2jvE2HWYloyIjQTWhsbExmpubs07DzKyqSHopIhoL476D28zMUrlYmJlZKhcLqy6+b8IsE9511qqH75swy4x7FlY9fN+EWWZcLKx6+L4Js8x4GMqqh++bMMuMi4VVF+8Oa5YJD0OZmVkqFwszM0vlYmFmZqlcLMzMLJWLhWXDd2KbVRWvhrLK853YZlXHPQurPN+JbVZ1XCys8nwntlnV8TCUVZ7vxDarOi4Wlg3fiW1WVTwMZWZmqVwszMwslYuFmZmlcrEwM7NULhZmZpbKxcIOjbfpMKtLXjprHedtOszqlnsW1nHepsOsbrlYWMd5mw6zuuVhKOs4b9NhVrfK1rOQNEvSW5LW5sWmSXpN0mpJ8yX1SuKjJb0kaU3y9dK8z5yfxDdIuk+SypWzdcDw4dDU5EJhVmfKOQz1MDCmILYEGBwRQ4D1QFMSfxv4bEScDUwAHs37zAPAJGBA8iq8ppmZlVnZikVELAO2FcQWR0RrcrgCaEjiL0fExiS+Dugq6ShJ/YAeEbE8IgKYDVxZrpzNzKy4LCe4JwKLisQ/D7wcEe8DJwEtee+1JLGiJE2S1CypecuWLSVN1sysnmVSLCRNAVqBuQXxQcDdwNfaQkU+Hge7bkTMiIjGiGjs27dvqdI1M6t7FV8NJWkCMBa4LBlaaos3APOB8RHxehJuIRmqSjQAGzEzs4qqaM9C0hjgVmBcROzKi/cC/gloioj/2xaPiE3ATknDklVQ44EnK5mzmZmVd+nsPGA5MFBSi6RrgfuB7sASSSsl/SA5/Ubg3wBTk/hKSSck710PPARsAF6n+DyHmZmVkfJGgmpKY2NjNDc3Z52GmVlVkfRSRDQWxr3dR63zLrFmVgLe7qOWeZdYMysR9yxqmXeJNbMScbGoZd4l1sxKxMNQtcy7xJpZibhY1Lrhw10kzOyweRjKzMxSuViYmVkqFwszM0vlYmFmZqlcLMzMLJWLhZmZpXKxMDOzVC4WZmaWysXCzMxSuViYmVkqFwszM0vlYpE1P5zIzKqANxLMkh9OZGZVwj2LLPnhRGZWJVwssuSHE5lZlfAwVJb8cCIzqxIuFlnzw4nMrAp4GMrMzFK5WJiZWSoXCzMzS+ViYWZmqcpWLCTNkvSWpLV5sWmSXpO0WtJ8Sb2SeB9Jz0r6vaT7C65zvqQ1kjZIuk+SypWzmZkVV86excPAmILYEmBwRAwB1gNNSXw3MBX4ZpHrPABMAgYkr8JrmplZmZWtWETEMmBbQWxxRLQmhyuAhiT+XkT8nFzR2EdSP6BHRCyPiABmA1eWK2czMysuyzmLicCilHNOAlryjluSmJmZVVAmxULSFKAVmJt2apFYtHPdSZKaJTVv2bLlcFI0M7M8FS8WkiYAY4FrkqGl9rSQDFUlGoCNBzs5ImZERGNENPbt2/fwkzUzM6DCxULSGOBWYFxE7Eo7PyI2ATslDUtWQY0HnixzmmZmVqBse0NJmgeMAo6X1ALcRm7101HAkmQF7IqImJyc/y9AD+Cjkq4ELo+IV4Drya2s6kZujiNtnsPMzEqsbMUiIq4uEp7Zzvn9DxJvBgaXKC0zM/sQfAd3IT/m1MzsAN6iPJ8fc2pmVpR7Fvn8mFMzs6JcLPL5MadmZkV5GCqfH3NqZlaUi0UhP+bUzOwAHoYyM7NULhZmZpbKxcLMzFK5WJiZWSoXCzMzS+ViYWZmqZT+SInqJGkL8Nus8yiR44G3s06ihGqpPbXUFqit9tRSW6By7fmziDjggUA1WyxqiaTmiGjMOo9SqaX21FJboLbaU0ttgezb42EoMzNL5WJhZmapXCyqw4ysEyixWmpPLbUFaqs9tdQWyLg9nrMwM7NU7lmYmVkqFwszM0vlYlEFJH1TUkg6Pi/WJGmDpF9J+nSW+XWEpGmSXpO0WtJ8Sb3y3quqtrSRNCbJeYOkb2Wdz6GQdLKkZyW9KmmdpJuS+HGSlkj6dfK1d9a5dpSkIyS9LGlhclzNbekl6Ynk78yrkoZn3R4Xi05O0snAaODNvNhZwJeBQcAY4H9KOiKbDDtsCTA4IoYA64EmqNq2kOT4P4B/C5wFXJ20pVq0Av85Is4EhgE3JPl/C3gmIgYAzyTH1eIm4NW842puy98DT0XEGcBQcu3KtD0uFp3fvcAtQP5KhM8B/xAR70fEG8AG4JNZJNdREbE4IlqTwxVAQ/J91bUl8UlgQ0T8JiL+CPwDubZUhYjYFBG/TL7fSe4fo5PIteGR5LRHgCszSfAQSWoArgAeygtXa1t6ACOAmQAR8ceI2E7G7XGx6MQkjQP+X0SsKnjrJOBf845bkli1mAgsSr6v1rZUa94HkNQfOBd4HvhYRGyCXEEBTsgwtUPxfXK/VO3Ni1VrW/4c2AL8r2RY7SFJx5Bxe/xY1YxJ+j/AiUXemgL8DXB5sY8ViWW+Brq9tkTEk8k5U8gNgcxt+1iR8zNvSwdUa977kXQs8I/A1yNih1SsWZ2bpLHAWxHxkqRRGadTCl2A84C/iojnJf09nWAIzcUiYxHxqWJxSWcDpwKrkr/ADcAvJX2S3G+xJ+ed3gBsLHOqqQ7WljaSJgBjgcvigxt8OmVbOqBa895H0pHkCsXciPjfSfh3kvpFxCZJ/YC3ssuwwy4Cxkn6DNAV6CFpDtXZFsj92WqJiOeT4yfIFYtM2+NhqE4qItZExAkR0T8i+pP7A3ReRGwGFgBflnSUpFOBAcALGaabStIY4FZgXETsynur6tqSeBEYIOlUSR8lN0m/IOOcOky530BmAq9GxN/lvbUAmJB8PwF4stK5HaqIaIqIhuTvyZeBn0bEv6cK2wKQ/B3/V0kDk9BlwCtk3B73LKpQRKyT9Di5P0CtwA0R8aeM00pzP3AUsCTpKa2IiMlV2hYiolXSjcDTwBHArIhYl3Fah+Ii4C+BNZJWJrG/Af4r8Lika8mtwPtiNumVRDW35a+AuckvIr8B/gO5X+4za4+3+zAzs1QehjIzs1QuFmZmlsrFwszMUrlYmJlZKhcLMzNL5WJhVmGSnpK0vW13VLNq4GJhVnnTyN3jYFY1XCzMykTSJ5Lnd3SVdEzy3IjBEfEMsDPr/MwOhe/gNiuTiHhR0gLge0A3YE5ErM04LbMPxcXCrLy+Q24fqd3AX2eci9mH5mEos/I6DjgW6E5uR1SzquRiYVZeM4Cp5J7fcXfGuZh9aB6GMisTSeOB1oh4LHlm9y8kXQr8LXAGcKykFuDaiHg6y1zN0njXWTMzS+VhKDMzS+ViYWZmqVwszMwslYuFmZmlcrEwM7NULhZmZpbKxcLMzFL9f7GxEl2spe4iAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(x01, y1, 'r.', label = 'do(x1 = a, x2 = b)')\n",
    "\n",
    "plt.legend()\n",
    "plt.xlabel(\"x1\")\n",
    "plt.ylabel(\"y\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "f34a2b83-5b53-4e3e-87df-86e44dab93bc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEGCAYAAACUzrmNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAAAbNUlEQVR4nO3de5CU9b3n8ffngAoabiJG4sQDcZEoBDVMIugRMAbCHq+sSSWWteCSyOLG2phU1FAcIhtjLI+yplyrjEQpxQsniac4EnKIsEbEbIajQ+TqBSHmuLOAIogQEQLy3T/6GWybnnka6Ht/XlVdPf3t7me+P2aa7/wuz+9RRGBmZtaZv6l0AmZmVv1cLMzMLJWLhZmZpXKxMDOzVC4WZmaWqmulEyiVk046KQYMGFDpNMzMasqKFSveiYh+ufG6LRYDBgygtbW10mmYmdUUSf+eL+5hKDMzS+ViYWZmqVwszMwsVd3OWeSzb98+2tra2LNnT6VTsTLr1q0bTU1NHHPMMZVOxawmNVSxaGtro0ePHgwYMABJlU7HyiQi2LZtG21tbQwcOLDS6ZjVpIYahtqzZw99+/Z1oWgwkujbt697lGZHoaGKBeBC0aD8c7eG0NICd9yRuS+yhhqGMjOrWy0tcPHF8Ne/wrHHwjPPwMiRRTt8w/UszMzq0tKlmULx4YeZ+6VLi3p4F4sKmjlzJnfffXfq62688UaWLVvW6WumT5/Opz/9aT7xiU8UK72iuemmm/jsZz/LsGHDmDBhAjt27DjqY44ZMybvGfoLFy7k1ltvPerjm9WcMWMyPYouXTL3Y8YU9fAuFmlKOAZYiO3bt7N8+XJGjRrV6esuu+wyXnjhhTJldXjGjh3L2rVrWb16NWeccQZ33HFHyb7XJZdcwoIFC9i9e3fJvodZVRo5MjP0dNttRR+CAheLzrWPAc6YkbkvQsG4/fbbGTx4MF/+8pd57bXXDsZXrlzJiBEjDv71/e677wLw5JNPMn78eADee+89Bg8efPB9V199NT//+c8BGDFiBP379z/q/LJdf/31NDc3M2TIkIL+Wr/iiiuYO3cuAA888ADXXHMNAOPGjaNr164H82xraytKfo899hjnn38+Q4cOPVgoJTFmzBgWLlxYlO9hVlNGjoRp04peKIDMGvRS3IA5wNvA2qzYXcCrwGpgPtA767lhQAuwDlgDdEviw5PHG4B7ARXy/YcPHx65Xn755UNinfrJTyK6dImAzP1PfnJ478/R2toaQ4cOjffffz/ee++9OP300+Ouu+6KiIjPfe5zsXTp0oiImDFjRnznO9+JiIiJEyfGggULDh5j8eLFMWLEiJg3b1585StfOeR7nHDCCR1+/9/97ndx9tlnH3IbOXJk3tdv27YtIiL2798fo0ePjlWrVnXavi1btsTpp58ey5Yti0GDBh18f7ZLL700Hn300UPiO3fuzJvb2WefHevWrTvk9aNHj45vfetbERHx3HPPxZAhQw4+99hjj8UNN9xwyHsO++dv1oCA1sjzf2opV0M9DNwHzM2KLQGmRcR+SXcC04BbJHUFHgP+c0SsktQX2Je8535gCrAc+FdgPLCohHl/pH0MsH11wVGOAT7//PNMmDCB448/HoDLL78cyPQYduzYwejRowGYNGkSX/va1wDYvHkz/fp9tFvw2LFj+dWvfsW3v/1tVq1adVjf/6KLLmLlypUFv/6Xv/wls2fPZv/+/WzevJmXX36ZYcOGdfj6T37yk/zoRz/ioosuYv78+Zx44okfe/7222+na9euB3sc2Xr06HFYuUGmZwUwatQodu7cyY4dO+jduzcnn3wymzZtOqxjmVnnSlYsImKZpAE5scVZD5cDX02+HgesjohVyeu2AUjqD/SMiJbk8VzgSspVLNrHAJcuzRSKInTtDne9f/fu3T92MtmBAwd45ZVX6N69O9u3b6epqangYz377LN897vfPSR+/PHH84c//OFjsTfeeIO7776bF198kT59+nDttdcWdFLbmjVr6Nu37yH/WT/yyCMsXLiQZ555Ju+/wa5du7jwwgvzHvOJJ57grLPOOiSee5z2x3v27KF79+6puZpZ4So5ZzGZj/7TPwMISU9L+qOkm5P4qUD2AHdbEiufIo4Bjho1ivnz5/PBBx+wa9cufv3rXwPQq1cv+vTpw/PPPw/Ao48+erCXceaZZ7Jhw4aDx7jnnns488wzmTdvHpMnT2bfvn2HfqMOtPcscm+5hQJg586dnHDCCfTq1Yu33nqLRYs+qs/Tpk1j/vz5h7znhRdeYNGiRbz00kvcfffdvPHGGwD89re/5c4772TBggUHe1W52nsW+W75CgXAL37xCwB+//vf06tXL3r16gXA+vXrGTp0aMH/LmaWriLFQtJ0YD/weBLqCvwdcE1yP0HSxUC+P8Ojk+NOkdQqqXXr1q1Fzvroff7zn+frX/8655xzDlddddXH/pJ+5JFHuOmmmxg2bBgrV67khz/8IZBZ3bM0WS+9fv16HnzwQWbNmsWFF17IqFGj+PGPfwzAzTffTFNTE7t376apqYmZM2ceVa5nn3025557LkOGDGHy5MlccMEFB59bs2YNp5xyysdev3fvXq677jrmzJnDpz71KWbNmsXkyZOJCG644QZ27drF2LFjOeecc5g6depR5dauT58+nH/++UydOpWHHnroYPzZZ5/lkksuKcr3MLNEvomMYt2AAWRNcCexSWQmso/Pin0DeDjr8QzgJqA/8GpW/GrggUK+d1EmuKvEBRdcEO+++26l0zho3LhxlU6hQ1u2bIkvfelLeZ+r1Z+/WTnRwQR3WXsWksYDtwCXR0T2QvingWGSjk8mu0cDL0fEZmCXpBHKDEhPBJ4qZ87VYNasWbz55puVTuOgp59+utIpdOjNN99k1qxZlU7DrO6UbIJb0jxgDHCSpDbgVjKrn44DliSTkcsjYmpEvCvpfwIvkhlm+teI+E1yqOvJrKzqTmaO46gmtyOi5jaVO++88yqdQs34whe+kDee+YPJzI5UKVdDXZ0n/FCeWPvrHyOzfDY33goUZbayW7dubNu2zduUN5hIrmfRrVu3SqdiVrMaatfZpqYm2traqMbJbyut9ivlmdmRaahiccwxx/hKaWZmR8B7Q5mZWSoXCzMzS+ViYWZWShW+zEGxNNSchZlZWZX4Uqfl5J6FmVmplPhSp+XkYmFmViolvtRpOXkYysysVEpwmYNKcbEwMyulkSNruki08zCUmZmlcrEwM7NULhZmZpbKxcLMzFK5WJiZWSoXCzMzS+ViYWZmqVwszMwslYuFmZmlcrEwM7NULhZmZpbKxcLMGludXJyo1LyRoJk1rjq6OFGpuWdhZo2rji5OVGouFmbWuOro4kSl5mEoM2tcdXRxolJzsTCzxlYnFycqNQ9DmZlZqpIVC0lzJL0taW1W7C5Jr0paLWm+pN5JfICkDyStTG4/y3rPcElrJG2QdK8klSpnMzPLr5Q9i4eB8TmxJcDQiBgGrAemZT23MSLOSW5Ts+L3A1OAQckt95hmZlZiJSsWEbEM2J4TWxwR+5OHy4Gmzo4hqT/QMyJaIiKAucCVJUjXzMw6Uck5i8nAoqzHAyW9JOk5SRcmsVOBtqzXtCWxvCRNkdQqqXXr1q3Fz9jMrEFVpFhImg7sBx5PQpuB0yLiXOB7wBOSegL55ieio+NGxOyIaI6I5n79+hU7bTOzhlX2pbOSJgGXAhcnQ0tExF5gb/L1CkkbgTPI9CSyh6qagE3lzdjMzMras5A0HrgFuDwidmfF+0nqknz9GTIT2X+KiM3ALkkjklVQE4GnypmzmZmVdunsPKAFGCypTdI3gfuAHsCSnCWyo4DVklYBTwJTI6J9cvx64EFgA7CRj89zmFkj8Q6xFaNkJKjuNDc3R2tra6XTMLNi8Q6xZSFpRUQ058Z9BreZ1QbvEFtRLhZmVhu8Q2xFeSNBM6sN3iG2olwszKx2eIfYivEwlJmZpXKxMDOzVC4WZmaWysXCzMxSuViYmVkqFwszM0vlYmFm5eF9nWqaz7Mws9Lzvk41zz0LMys97+tU81wszKz0vK9TzfMwlJmVnvd1qnkuFmZWHt7XqaZ5GMrMzFK5WJiZWSoXCzMzS+ViYWaHxyfXNSRPcJtZ4XxyXcNyz8LMCueT6xqWi4WZFc4n1zUsD0OZWeF8cl3DcrEws8Pjk+sakoehzMwslYuFmZmlcrEwa3Q+b8IKULJiIWmOpLclrc2K3SXpVUmrJc2X1DvnPadJ+ouk72fFhktaI2mDpHslqVQ5mzWc9vMmZszI3LtgWAdK2bN4GBifE1sCDI2IYcB6YFrO8/cAi3Ji9wNTgEHJLfeYZnakfN6EFahkxSIilgHbc2KLI2J/8nA50NT+nKQrgT8B67Ji/YGeEdESEQHMBa4sVc5mDcfnTViBKrl0djLwCwBJJwC3AGOB72e95lSgLetxWxLLS9IUMr0QTjvttCKna1aHfN6EFagixULSdGA/8HgS+h/APRHxl5wpiXzzE9HRcSNiNjAboLm5ucPXmVkWnzdhBSh7sZA0CbgUuDgZWgI4D/iqpH8EegMHJO0B/pmsoark601lTNesdrW0uMdgRVPWYiFpPJnhptERsbs9HhEXZr1mJvCXiLgvebxL0gjg34CJwP8qZ85mNcm7w1qRlXLp7DygBRgsqU3SN4H7gB7AEkkrJf2sgENdDzwIbAA2cuhqKTPL5VVOVmQl61lExNV5wg8V8L6ZOY9bgaFFSsusMbSvcmrvWXiVkx0lbyRoVo+8ysmKzMXCrFalTWB7lZMVkYuFWS3yBLaVmTcSNKtFnsC2MnOxMKtF3qbDyszDUGbVLt/chCewrcxcLMyqWWdzE57AtjLyMJRZtcm+GJHnJqxKuGdhVk1yexI//alPrrOq4GJhVg3aexFvvvnxnsS2bZ6bsKrgYmFWDu3FoG9feOkl2LIlEz/lFDj3XLjxxkxx6NIFuiYfy/aehOcmrAqkFgtJNwCPR8S7ZcjHrH60F4gdO+Cee2D/fog8l1np0iUTP3Ag8/i66+C009yTsKpSSM/iFOBFSX8E5gBPZ12HwszyaZ972Lv3oyLQkQ8/hGOOASnTm5g40UXCqk5qsYiIf5A0AxgH/BfgPkm/BB6KiI2lTtCsZmSfD9G+iimtUAAcdxzce29mfsK9CatSBc1ZRERI2gJsIXM51D7Ak5KWRMTNpUzQrCZ0tIqpvWchZYabvvc92Lnz43MW7klYDShkzuK/A5OAd8hchOimiNgn6W+A1wEXC2tM+XoS+VYx9e3rXoPVvEJ6FicB/yki/j07GBEHJF1amrTMqlwh50N4FZPVkULmLH7YyXOvFDcdsxrRWU/CPQirQz7PwuxI5LtsqXsSVsdcLMzSeNdXMxcLs05511czwLvOmnXOu76aAS4WZp3zFenMAA9DmXXOcxNmgIuFWf4J7GyemzBzsbAG19kEtpkd5DkLa2yewDYrSMmKhaQ5kt6WtDYrdpekVyWtljRfUu8k/kVJK5PbKkkTst4zXNIaSRsk3StJpcrZGpAnsM0KUsqexcPA+JzYEmBoRAwD1gPTkvhaoDkizkne84Ck9iGy+4EpwKDklntMsyPXPoF9220egjLrRMnmLCJimaQBObHFWQ+XA19N4ruz4t2AAJDUH+gZES3J47nAlcCiUuVtDcgT2GapKjlnMZms//QlnSdpHbAGmBoR+4FTgbas97QlsbwkTZHUKql169atJUrbakpLC9xxR+bezI5YRVZDSZpO5iJKj7fHIuLfgCGSzgQekbQIyDc/0eElXSNiNjAboLm52Zd+bXRe6WRWNGXvWUiaBFwKXJPvWt7JtufvA0PJ9CSasp5uAjaVI0+rA17pZFY0ZS0WksYDtwCXZ89TSBrYPqEt6W+BwcCfI2IzsEvSiGQV1ETgqXLmbDXMK53MiqZkw1CS5gFjgJMktQG3kln9dBywJFkBuzwipgJ/B/xA0j7gAPDfIuKd5FDXk1lZ1Z3MHIcnt60w3qrDrGiUZySoLjQ3N0dra2ul0zAzqymSVkREc27cZ3CbmVkqFwszM0vlYmG1xedNmFWEd5212uHzJswqxj0Lqx0+b8KsYlwsrHb4vAmzivEwlNUOnzdhVjEuFlZbvEOsWUV4GMrMzFK5WJiZWSoXCzMzS+ViYWZmqVwsrPx8FrZZzfFqKCsvn4VtVpPcs7Dy8lnYZjXJxcLKy2dhm9UkD0NZefksbLOa5GJh5eezsM1qjoehzMwslYuFmZmlcrEwM7NULhZmZpbKxcLMzFK5WNjh8VYdZg3JS2etcN6qw6xhuWdhhfNWHWYNy8XCCuetOswaloehrHDeqsOsYZWsZyFpjqS3Ja3Nit0l6VVJqyXNl9Q7iY+VtELSmuT+S1nvGZ7EN0i6V5JKlbMVYORImDbNhcKswZRyGOphYHxObAkwNCKGAeuBaUn8HeCyiPgcMAl4NOs99wNTgEHJLfeYZmZWYiUrFhGxDNieE1scEfuTh8uBpiT+UkRsSuLrgG6SjpPUH+gZES0REcBc4MpS5WxmZvlVcoJ7MrAoT/wq4KWI2AucCrRlPdeWxPKSNEVSq6TWrVu3FjVZM7NGVpFiIWk6sB94PCc+BLgT+K/toTxvj46OGxGzI6I5Ipr79etXrHTNzBpe2VdDSZoEXApcnAwttcebgPnAxIjYmITbSIaqEk3AJszMrKzK2rOQNB64Bbg8InZnxXsDvwGmRcT/aY9HxGZgl6QRySqoicBT5czZzMxKu3R2HtACDJbUJumbwH1AD2CJpJWSfpa8/AbgPwAzkvhKSScnz10PPAhsADaSf57DzMxKSFkjQXWlubk5WltbK52GmVlNkbQiIppz497uo555h1gzKxJv91GvvEOsmRWRexb1yjvEmlkRuVjUK+8Qa2ZF5GGoeuUdYs2siFws6tnIkS4SZlYUHoYyM7NULhZmZpbKxcLMzFK5WJiZWSoXCzMzS+ViYWZmqVwszMwslYuFmZmlcrEwM7NULhZmZpbKxcLMzFK5WFSSL05kZjXCGwlWii9OZGY1xD2LSvHFicyshrhYVIovTmRmNcTDUJXiixOZWQ1xsagkX5zIzGqEh6HMzCyVi4WZmaVysTAzs1QuFmZmlqpkxULSHElvS1qbFbtL0quSVkuaL6l3Eu8r6VlJf5F0X85xhktaI2mDpHslqVQ5m5lZfqXsWTwMjM+JLQGGRsQwYD0wLYnvAWYA389znPuBKcCg5JZ7TDMzK7GSFYuIWAZsz4ktjoj9ycPlQFMSfz8ifk+maBwkqT/QMyJaIiKAucCVpcrZzMzyq+ScxWRgUcprTgXash63JTEzMyujihQLSdOB/cDjaS/NE4tOjjtFUquk1q1btx5NimZmlqXsxULSJOBS4JpkaKkzbSRDVYkmYFNHL46I2RHRHBHN/fr1O/pkzcwMKHOxkDQeuAW4PCJ2p70+IjYDuySNSFZBTQSeKnGaZmaWo2R7Q0maB4wBTpLUBtxKZvXTccCSZAXs8oiYmrz+z0BP4FhJVwLjIuJl4HoyK6u6k5njSJvnMDOzIitZsYiIq/OEH+rk9QM6iLcCQ4uUlpmZHQGfwZ3Llzo1MzuEtyjP5kudmpnl5Z5FNl/q1MwsLxeLbL7UqZlZXh6GyuZLnZqZ5eVikcuXOjUzO4SHoczMLJWLhZmZpXKxMDOzVC4WZmaWysXCzMxSuViYmVkqpV9SojZJ2gW8Vuk8iuAk4J1KJ1Ekbkt1qpe21Es7oLJt+duIOOSCQPV8nsVrEdFc6SSOlqTWemgHuC3Vql7aUi/tgOpsi4ehzMwslYuFmZmlqudiMbvSCRRJvbQD3JZqVS9tqZd2QBW2pW4nuM3MrHjquWdhZmZF4mJhZmapar5YSPqapHWSDkhqzooPkPSBpJXJ7WdZzw2XtEbSBkn3SlJlsv+4jtqSPDctyfc1SV/JildlW7JJminp/2X9LP4+67m87apWksYnuW6Q9INK53O4JP05+X1ZKak1iZ0oaYmk15P7PpXOMx9JcyS9LWltVqzD3Kv5d6uDtlT35yQiavoGnAkMBpYCzVnxAcDaDt7zAjASELAI+I+VbkdKW84CVgHHAQOBjUCXam5LTrtmAt/PE++wXdV4A7okOX4GODbJ/axK53WYbfgzcFJO7B+BHyRf/wC4s9J5dpD7KODz2Z/rjnKv9t+tDtpS1Z+Tmu9ZRMQrEVHwmdqS+gM9I6IlMj+JucCVpcrvcHTSliuAf4qIvRHxBrAB+GI1t6VAedtV4Zw680VgQ0T8KSL+CvwTmTbUuiuAR5KvH6FKf4ciYhmwPSfcUe5V/bvVQVs6UhVtqflikWKgpJckPSfpwiR2KtCW9Zq2JFbNTgX+b9bj9pxrqS03SFqddL/bhwo6ale1qrV88wlgsaQVkqYksU9GxGaA5P7kimV3+DrKvVZ/VlX7OamJ7T4k/W/glDxPTY+Ipzp422bgtIjYJmk48C+ShpAZrslVtvXDR9iWjnKuaFuyddYu4H7gNjK53QbMAiZTRfkXqNbyzeeCiNgk6WRgiaRXK51QidTiz6qqPyc1USwi4stH8J69wN7k6xWSNgJnkKnKTVkvbQI2FSPPAvM67LaQyfnTWY/bc65oW7IV2i5JPwcWJg87ale1qrV8DxERm5L7tyXNJzOc8Zak/hGxORnafLuiSR6ejnKvuZ9VRLzV/nU1fk7qdhhKUj9JXZKvPwMMAv6UdFV3SRqRrByaCHT0F321WAB8Q9JxkgaSacsLtdKW5EPcbgLQvgIkb7vKnd9heBEYJGmgpGOBb5BpQ02QdIKkHu1fA+PI/CwWAJOSl02iCn+HOtFR7rX2u1X9n5NKrwoowqqCCWQq717gLeDpJH4VsI7MKoI/ApdlvaeZzA9iI3AfyZnslb511JbkuelJvq+RteKpWtuS065HgTXAajK/+P3T2lWtN+DvgfVJztMrnc9h5v6Z5POwKvlsTE/ifYFngNeT+xMrnWsH+c8jM7y8L/mcfLOz3Kv5d6uDtlT158TbfZiZWaq6HYYyM7PicbEwM7NULhZmZpbKxcLMzFK5WJiZWSoXC7MyknSOpJZkd+HVkr5e6ZzMCuGls2ZlJOkMICLidUmfAlYAZ0bEjspmZtY59yzMSkTSF5LeQ7fk7Ol1wLER8Toc3HrjbaBfRRM1K4B7FmYlJOnHQDegO9AWEXdkPfdFMttqD4mIAxVK0awgLhZmJZTsIfUisAc4PyI+TOL9yVzkalJELK9chmaF8TCUWWmdCHwC6EGmh4GknsBvgH9wobBa4Z6FWQlJWkDminoDgf7A98hc/vbXEfHTCqZmdlhq4noWZrVI0kRgf0Q8kWyX/wcy25qPAvpKujZ56bURsbIyWZoVxj0LMzNL5TkLMzNL5WJhZmapXCzMzCyVi4WZmaVysTAzs1QuFmZmlsrFwszMUv1/1z/iX3p/sC0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(x02, y1, 'r.', label = 'do(x1 = a, x2 = b)')\n",
    "\n",
    "plt.legend()\n",
    "plt.xlabel(\"x2\")\n",
    "plt.ylabel(\"y\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "87effa40-0e7f-4620-a14d-78f3ceb57bb0",
   "metadata": {},
   "source": [
    "Since y ~ x1 and y ~ x2, we confirm that there is a confounder of y."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6a3f20d5-7549-47eb-b8ac-01ffa555485e",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}