logTransformedFit_errors.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "nteract": {
     "transient": {
      "deleting": false
     }
    }
   },
   "source": [
    "# Protein stability as a log-linear problem\n",
    "## Assessing model error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "jupyter": {
     "outputs_hidden": false,
     "source_hidden": false
    },
    "nteract": {
     "transient": {
      "deleting": false
     }
    }
   },
   "outputs": [],
   "source": [
    "# Modules\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "from scipy.stats import linregress\n",
    "import matplotlib.pyplot as plt\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nteract": {
     "transient": {
      "deleting": false
     }
    }
   },
   "source": [
    "## Same mock dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "jupyter": {
     "outputs_hidden": false,
     "source_hidden": false
    },
    "nteract": {
     "transient": {
      "deleting": false
     }
    }
   },
   "outputs": [
    {
     "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>GdmCl</th>\n",
       "      <th>Intensity</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.25</td>\n",
       "      <td>101.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.50</td>\n",
       "      <td>97.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.75</td>\n",
       "      <td>102.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1.00</td>\n",
       "      <td>96.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1.50</td>\n",
       "      <td>88.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>1.75</td>\n",
       "      <td>62.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>2.00</td>\n",
       "      <td>52.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>2.50</td>\n",
       "      <td>48.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>2.75</td>\n",
       "      <td>28.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>3.00</td>\n",
       "      <td>20.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>3.50</td>\n",
       "      <td>26.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>4.00</td>\n",
       "      <td>23.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    GdmCl  Intensity\n",
       "0    0.25      101.0\n",
       "1    0.50       97.0\n",
       "2    0.75      102.0\n",
       "3    1.00       96.0\n",
       "4    1.50       88.0\n",
       "5    1.75       62.0\n",
       "6    2.00       52.0\n",
       "7    2.50       48.0\n",
       "8    2.75       28.0\n",
       "9    3.00       20.0\n",
       "10   3.50       26.0\n",
       "11   4.00       23.0"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df = pd.DataFrame({\n",
    "    \"GdmCl\": np.array([\n",
    "        0.25,\n",
    "        0.50,\n",
    "        0.75,\n",
    "        1,\n",
    "        1.5,\n",
    "        1.75,\n",
    "        2,\n",
    "        2.5,\n",
    "        2.75,\n",
    "        3,\n",
    "        3.5,\n",
    "        4\n",
    "    ], dtype=float),\n",
    "    \"Intensity\": np.array([\n",
    "        101,\n",
    "        97,\n",
    "        102,\n",
    "        96,\n",
    "        88,\n",
    "        62,\n",
    "        52,\n",
    "        48,\n",
    "        28,\n",
    "        20,  # 20 or 19\n",
    "        26,\n",
    "        23\n",
    "    ], dtype=float)\n",
    "})\n",
    "\n",
    "df.head(12)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nteract": {
     "transient": {
      "deleting": false
     }
    }
   },
   "source": [
    "## Same data manipulations and same transition region."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "jupyter": {
     "outputs_hidden": false,
     "source_hidden": false
    },
    "nteract": {
     "transient": {
      "deleting": false
     }
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/lib/python3.9/site-packages/pandas/core/base.py:937: RuntimeWarning: invalid value encountered in log10\n",
      "  new_values = map_f(values, mapper)\n"
     ]
    }
   ],
   "source": [
    "def rescale(y, ref_min, ref_max):\n",
    "    return (y - ref_min) / (ref_max - ref_min)\n",
    "\n",
    "\n",
    "df['N_ratio'] = df.Intensity.apply(rescale, args=(20, 98))\n",
    "df['K'] = (1 - df.N_ratio) / df.N_ratio\n",
    "df['logK'] = df.K.map(np.log10)\n",
    "transition = df[df.GdmCl.between(1.25, 2.80)].copy()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nteract": {
     "transient": {
      "deleting": false
     }
    }
   },
   "source": [
    "Recover paramaters estimated from fitted model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "jupyter": {
     "outputs_hidden": false,
     "source_hidden": false
    },
    "nteract": {
     "transient": {
      "deleting": false
     }
    }
   },
   "outputs": [],
   "source": [
    "INTERCEPT = np.array(-2.293263723543123)\n",
    "SLOPE = np.array(1.1350753137863376)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nteract": {
     "transient": {
      "deleting": false
     }
    }
   },
   "source": [
    "For any model fitted to empirical data, there will be deviations between the\n",
    "model predictions and the actual measurements. These are due to any factor\n",
    "outside of the control of the experimental design, think anything from\n",
    "thermal fluctuations, sample heterogeinity, instrument inaccuracy, pipetting...\n",
    "The distances between observations and model $y_{obs} - y_{pred}$ are called\n",
    "**residuals** and allow us to quantify the model errors. \n",
    "\n",
    "Below, we calculate the residuals and squared residuals for the transition\n",
    "region."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "jupyter": {
     "outputs_hidden": false,
     "source_hidden": false
    },
    "nteract": {
     "transient": {
      "deleting": false
     }
    }
   },
   "outputs": [],
   "source": [
    "transition['logK_pred'] = INTERCEPT + SLOPE * transition.GdmCl\n",
    "transition['logK_res'] = transition.logK - transition.logK_pred\n",
    "transition['sqres'] = transition.logK_res ** 2\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nteract": {
     "transient": {
      "deleting": false
     }
    }
   },
   "source": [
    "Here we represent residuals as vertical arrows between observations and\n",
    "predictions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "jupyter": {
     "outputs_hidden": false,
     "source_hidden": false
    },
    "nteract": {
     "transient": {
      "deleting": false
     }
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1288efe50>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAFxCAYAAABgNZ+NAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAA0E0lEQVR4nO3deXhU5fnG8e8TCATEskRcECEqKigqShBZhIC4L1VwAVFQUFyq1p8bttFqtWkpglJtFRFaRQMUxY0qKAgIiKhhERVpq1UwUhQjskgICXl+f5yBxhiSCSRzJpn7c11zceacMzPPJEzuec97zvuauyMiIlKRpLALEBGRmkGBISIiUVFgiIhIVBQYIiISFQWGiIhERYEhIiJRqRt2AdVhv/3287S0tLDLEBGpMZYsWfKtuzcvb59aGRhpaWnk5OSEXYaISI1hZqsr2ieUQ1JmdqmZjShj/W1mtszMlphZ38i6lmY238zeNbMXzWyf2FcsIiIxDQwzSzKzWcDTZWw7AugPdAJOBUaaWT3gD8DD7t4ZWAbcEMOSRUQkIqaB4e7FwFmU/Ue/N/APdy9y9++BVUB7oCcwPbLPa8ApMShVRERKiXkfhrsXmVlxGZtSgbwS9zcCTYC67l5Uat1PmNkwYBhAq1atfrK9sLCQ3Nxctm3btse1S3RSUlJo2bIlycnJYZciIlUonjq9vyMIjZ2aAeuB7WZWx913lFj3E+4+DhgHkJ6e/pMRFXNzc9l3331JS0vDzKq8eAm4O3l5eeTm5nLooYeGXY6IVKF4ug5jDnBBpJ+jOZAGfAy8DZwT2edC4PU9efJt27aRmpqqsKhmZkZqaqpaciK1UOiBYWaZZnamu/8LmEbQsT0DuCnS55EJDDezd4FWwN/24rWqomSpgH7OIrVTKIHh7k+5+12R5Sx3nxlZHuHux7t7urvPjqz7wt27uXtnd7/M3QvDqLk6pKWlhfJN/KWXXuLTTz8FoG/fvjF/fRGpmUJvYcSzvC0FfPDl9+RtKQi7lHIVF5d1DsHulQyMF154oTpKEpFaKJ46vePKy8u/Yvi0FSQnJVFYXMzIfsdxfoeD9/j5Nm/ezLBhw1i7di2FhYX89re/BeDGG2/ks88+A+DZZ5+lTp06XH755RQWBg2pxx9/nKOPPpoHHniA2bNns3XrVm655RYGDhxIRkYGhx12GMXFxaxevZr77ruPnj17Mn/+fEaPHk12djaXXXYZmzZtYvPmzdxxxx0cdNBBzJw5k6VLl9KhQwc6dOjAunXrWLRoEXfddRcAqampTJgwgRUrVnDvvffSuHFjvvrqK0477TRGjBjBc889x4MPPkiDBg047LDDmDBhAklJ+u4hEoa8LQXkbsinZdMGpDaqX70v5u617taxY0cvbeXKlT9Ztzvfbt7mR939mrce/o9dt6Pufs2/3bwt6ucoLTMz00ePHu3u7l9//bW3aNHCDznkEF+4cKG7u48dO9avv/56f+qpp/yqq67yHTt2+JIlS3zGjBk+e/Zs79+/v7u75+fne9u2bT0vL8979uzpzz77rLu7T5w40a+//np3d7/66qt9+vTpvmLFCp80aZK7u7/zzjveu3dvd3cfPHiwz5gxw93dDzjgAHd3P/LIIz03N9fd3R955BG/+eabfe7cuX7UUUd5fn6+FxQUeJMmTXzHjh2ekZHhc+bMcXf38ePH+7p1637yfivz8xaRPfPSslw/6u7XvP1vZvpRd7/mLy/L3ePnAnK8gr+t+lpYhtwN+SSX+sacnJRE7ob8PX7OZcuW0bt3bwD2339/UlNTqVOnDp06dQKgc+fOfPHFFwwYMICjjjqK/v378+CDD/Kzn/2MpUuXkpOTQ0ZGBmeeeSZmxpo1awDo2LEjABdddNGuFsjixYs566yzaNy4MbNmzeK6667jiSeeYMeOHWXWtn79eho0aMDBBwctqG7duu06ZHX88ceTkpJCvXr1aNy4Mdu3b2fs2LE8//zzXHLJJSxfvpyGDRvu8c9FRPZM3pYChk9bwbbCYjYXFLGtsJg7p62o1kPoCowytGzagMJS/QKFxcW0bNpgj5+zTZs2LF68GIC1a9dSUFCAu7NgwQIAFixYQIcOHXj++ec59thjmTp1KjfeeCP3338/Rx11FGeccQbz5s1j1qxZXHDBBbRt2/ZHz9+gQQN69+7NrbfeyoUXXkidOnUYPXo0GRkZjB07ltNPP33XvqXPYmrWrBnff/8969cHl7i89dZbnHjiiWXuC/DQQw/x4IMPMnXqVLZt28arr766xz8XEdkzJb/YFm5YC+z9F9uKKDDKkNqoPiP7HUdKchL71q9LSnISI/sdt1fHB++++26mT59Oz549ufTSSxk/fjx169YlOzub7t2788Ybb3DHHXfQrVs3RowYQY8ePbj55pu59dZbOf/886lXrx4ZGRn07NmTww8/nJSUlJ+8xtChQ3nyyScZMmQIABdccAG///3v6dWrF19++SVr167l7bff5qSTTuLOO+/cFRB16tThscce49xzz6VXr14sXryYO+64Y7fvpVu3bvTo0YOuXbuydu1azjrrrD3+uYjIntn5xbbw2y9ZP+13uBfv9Rfbilhw6Kp2SU9P99LDm3/yySe0a9euUs8T086kWmZPft4iUjmvLP+KKwYNZtOHb9Ki3695/O4b9vjkHDNb4u7p5e2jFkY5UhvV5/hDmigsRCQuHdt4Oz+snAdA6qevcd7xLar19RQYIiI1VIMGDRgzZgwpKSmMfewvu07Hry66DkNEpIY68MADad++PXXq1KFr167V/npqYYiISFQUGCIiEhUFhoiIREWBsRuTs7Npn5ZGnaQk2qelMTk7O+ySuO+++xg7dizLly/n/vvv3+1+Y8aMoaioqML9REQqQ53eZZicnU3msGFM2LqV7sDC1asZOmwYAAMGDqzy1ysuLq7U4H07Bw3cnTFjxnDddddVuJ+ISGUkbgsjI+Ont8ceAyDrV79iwtat9AKSgV7AhK1byfrlL4PHfvvtTx9bgaeeeoqzzjqLPn36cOyxx/L444+TkZHBkCFDGDJkCDt27OAXv/gFvXr1onPnzsyePRuA7OxsTjzxRDIyMpg/fz4A8+bNo3///gA8+uijdO7cmU6dOjFmzBj+9Kc/sW7dOgYNGvSj/caNG0fXrl3p2rUrt99+OxC0WC6++GLOOussjjnmGA11LiLlUgujDJ/k5tK91LruwCd5eXv1vJs2bWLBggXk5+dz7LHHUq9ePa655hoGDhzIhAkT+NnPfsbcuXPZsGEDnTp14t133+WBBx5gyZIlNGzYkHPOOedHz7dq1SqeeeYZFi1aRHFxMYMGDWLSpEk8/PDDTJw4cdfYVf/6178YP348b7/9NsnJyfTr148XX3wRgH322YfnnnuOd955h8zMTE2oJCK7lbiBMW/ebje1a9WKhatX06vEuoVAu9atgzv77Vfu43ena9euJCUlsc8++9C+fXvmzJmza7TZpUuXMn/+fN555x0AkpKS+Oyzz2jbti377LMPACeddNKPnm/FihV0796dunWDX+OUKVPKfN0PPviArl27kpycDPx4NNrOnTsDcNBBB7F169ZKvycRSRyJe0iqHJlZWQxt2JC5QCEwFxjasCGZWVl79bzLli3D3fnhhx9YuXIlRxxxxK5tRx11FIMHD2bevHm88sorXHzxxRx66KGsXLmS/Px83J1Fixb96PmOOOKIXc+5ceNGevfuTXFx8U9GmG3Tpg05OTm7hjefP39+uaPRioiURYFRhgEDB5I1bhw3tW5Nihk3tW5N1rhxe93hXVRURK9evejSpQt33303jRs33rXt2muv5cMPP6RXr16cc845dO7cmebNm/OrX/2KLl26cN5555Gamvqj5zvhhBPo3bs3Xbt25bTTTuOaa64hKSmJLl26MGjQoB/td/bZZ9OlSxd69OhB+/btOfXUU/fqvYhI4tFotTHy1FNPsWrVKkaMGBFaDbEU9s9bJBFMzs7mnttu4z9ff83RrVuTmZW1x19soxmtNnH7MEREarBYn/4PamFINdHPW6R6tU9L49FSJ+fMBW5q3ZqPvvii0s+n+TBERGqpT9asKfv0/zVrqu01EyowamNrKh7p5yxS/dq1asXCUusWRtZXl4QJjJSUFPLy8vTHrJq5O3l5eWXOOS4iVcAdxowh85prquX0//IkTKd3y5Ytyc3NZf369WGXUuulpKTQsmXLsMsQqX1++AGuvhqmTGHAjTfCuHHclJnJJ2vW0K5VK7L24iypaCRMp7eISI326adw4YXw8ceQlQXDh0MlBi2tiE6rFRGpDT78EE45BerUgZkz4fTTQykjYfowRERqrLZt4fLLYcmS0MICFBgiIvHp+++D/opvvoHkZPjznyEtLdSSFBgiIvHmww8hPR2efhoiI1jHAwWGiEg8mTwZTj4Ztm6Ft96Cn/887Ip2UWCIiMSL8ePhssvgxBNh6VLo2jXsin4kZmdJmVky8DfgCGAHMMTdV0W2nQncVWL3ZsBfIvsvBb6NrJ/j7vfHqmYRkZi68ELIzYXMzKDfIs7EsoUxCMhz984E4TB65wZ3n+nuGe6eAVwCbAAmAW2ABTu3KSxEpNZZvBguuQS2b4fUVLjvvrgMC4htYPQBXogsLwA67Ga/R4Asd98MtAU6mdk8M3vFzI7YzWNERGoWd3jiCejRIzhddu3asCuqUCwDIxXIA/Dg8nI3sx+9vpl1AJq6+xuRVRuAUZGWx8MErY4ymdkwM8sxsxwN/yEicW3btuCU2euugz59ICcn9FNmoxHLwPgOaAxgwUTS7u7Fpfa5HfhziftvufsUgp3nAi1sN5NQu/s4d0939/TmzZtXffUiIlXlyivhr3+Fe+6B6dOhadOwK4pKLAPjTaBfZPkMgsNSu5hZfSADmFFi9RNmdlVk+3HAaq+Ng1+JSGLY+efr7rvh5Zfh/vuD4T5qiFiOJfU0MNHMcoAfgMvNLBNY4u4zgW5AjrsXlXjM/cDTkdAoAq6JYb0iIlXDHR58EP7zHxg7Ftq3D241TMwCw923A/1Lrc4qsX0OMKfUY74Eeld/dSIi1WTzZrjqKpg2LTgbqrAwbs+Cqogu3BMRqS6rVkHnzvDSSzBqFEyZUmPDAjS8uYhI9SgoCM6A2r4dZs2CXr3CrmivKTBERKrSjh3BxEb16weDBx55JBxySNhVVQkdkhIRqSp5eXD22cFQ5ACnnlprwgIUGCIiVWPp0mBI8nnzYJ99wq6mWigwRET21tNPQ7duUFQECxbAkCFhV1QtFBgiIntj5crgtNkuXYIxoU46KeyKqo06vUVE9sS2bZCSAkcfDa+/HpwFVbd2/0lVC0NEpLIWLIA2bWBO5Frj006r9WEBCgwRkei5w6OPQu/e0LAhHHBA2BXFlAJDRCQaW7fCoEFw883BqbPvvw/HHBN2VTGlwBARicbkyZCdDQ88AC++CI0bh11RzNX+g24iInvj+++hSZPgVNnjjw+utUhQamGIiJSluBh+97ugc/vzz8EsocMC1MIQEfmpjRvhiiuC2fAGDky4zu3dUWCIiJT00UfQt2/QqnjkEbjxxqB1IQoMEZEfGTMGNm0KrrE45ZSwq4kr6sMQESkqgv/+N1h+5JFgIEGFxU8oMEQksa1fD6efHgxFXlAQXJDXokXYVcUlHZISkcT1/vvQr18QGmPHBpMeyW6phSEiiWn8eOjePZgd7+23YfDgsCuKewoMEUk8hYXwxBPQs2cwJPmJJ4ZdUY2gQ1Iikji+/BJ+9rNgWI8ZM6BpU6hTJ+yqagy1MEQkMcydCx07wg03BPf3209hUUkKDBGp3dxh9Gjo0wdSU+Gee8KuqMbSISkRqb22bIGhQ2Hq1OBsqL/9DfbdN+yqaiy1MESk1picnU37tDTqJCXRPi2NyU89BQsXwogR8NxzCou9pBaGiNQKk7OzyRw2jAlbt9IdWLh6NUOHD4dHH2XAkCFhl1crmLuHXUOVS09P95ycnLDLEJEYap+WxqOrV9OrxLq5wE2tW/PRF1+EVFXNYWZL3L3c8dt1SEpEaoVP1qyhe6l13SPrpWooMESk5luxgnZ16rCw1OqFQLtWrcKoqFZSYIhIzZefT2bTpgxNSWEuUEhwOGpow4ZkZmWFXFztocAQkZqpsBBefjlY7tyZAV99Rdb48VzVrBn1CfoussaNY8DAgaGWWZvE7CwpM0sG/gYcAewAhrj7qhLbxwAnAsWRVacDDYDJQDNgC3C5u6+LVc0iEqfWrYNLLoEFC4K5K044AZKTGTBwIJt/+IGxY8eydOnSsKusdWLZwhgE5Ll7Z+AuYHSp7ccDGe6+87YduB2Y5e4nAxMBXaIpkujeeScYLDAnByZNCsJCYiKWgdEHeCGyvADoUGr7IcCrZrbQzAaV8ZjXAE2BJZLInnwyGGG2YUNYvBgGDAi7ooQSywv3UoE8AHd3M3MzS3L3YjNrAPwdyAKSgXlm9kHJxwAbgSa7e3IzGwYMA2ilsyJEaqd69YLZ8Z55JhhpVmIqli2M74DGAGZmBLmxs79iG3CPu291943ALODYko8h6MdYv7snd/dx7p7u7unNmzevrvcgIrH2xRcwfXqwPHhwsKywCEUsA+NNoF9k+QyCw1I7dQXetEBdoBuwtNRjLgRej1GtIhIPZs0KhiQfNgzy84N1ZuHWlMBiGRhPAy3MLAf4FTDczDLN7Ex3fxt4P3KbB0x295XAKOAcM3sPOB8YGcN6RSQs7sGAgWeeCS1aBGdDNWgQdlUJL2Z9GJGznvqXWp1VYvudZTxmA0FrREQSxY4dwSmzL7wA/fsHc2/vs0/YVQm6cE9E4k2dOtCmDTz8cHDarMIibmh4cxGJD9OmQcuW0Lkz/PGPYVcjZVALQ0TCtWMH3HUXXHQRjFQ3ZTxTC0NEwvPtt8HFd7Nnw7XXwp/+FHZFUg4FhoiEY80aOOUU+PprmDABNCte3FNgiEg4Dj4Y+vSB66+H9HInepM4oT4MEYmd7duD/oq1a4OzoSZMUFjUIAoMEYmNr74KBg784x//N9SH1Cg6JCUi1W/+fLj4YvjhB5g6NViWGkctDBGpXtOnQ+/e0KQJvPeewqIGU2CISPU65ZSgY/u99+Doo8OuRvaCAkNEqt5nn8GVV8K2bUHL4tFHoXHjih4lcU6BISJV67XXgjOfXnkF/vnPsKuRKqTAEJGqUVwM998P554LaWmwZAkcf3zYVUkVUmCISNW47Ta491644gpYtAgOPTTsiqSK6bRaEaka110HRx4Z/KtZ8WolBYaI7LkpU2DuXBg7Fo46KrhJraVDUiJSeYWFcOutwUizK1cGF+RJrafAEJHK+fprOO20YEa8m2+GOXOgUaOwq5IY0CEpEYlecXEQFp9+Cs88A5dfHnZFEkMKDBGpmHvwb1ISPPQQ7LcfdOgQakkSezokJSLl27YNrr4aRo0K7vfpo7BIUAoMEdm9nbPi/fWvsGVL2NVIyHRISkTK9uab0L9/MOnRSy/Bz38edkUSMrUwROSncnPh7LNh//3h/fcVFgKohSEiJRUVQd260LIlPP889OqlU2ZlF7UwRCTwz3/CCSfAjBnB/fPOU1jIjygwRCToo+jUCdatg5SUsKuROKXAEElkO3bA3XfDhRdC27bBkOS9eoVdlcQpBYZIInv5ZcjKgqFDYf58aNUq7IokjqnTWyQRbd0KDRsGLYvZs+HUU8OuSGoAtTBEEs0zz8BhhwWd3GYKC4maAkMkUWzfDjfdBIMGQbt20LRp2BVJDaPAEEkEa9cGndl//nMwleqsWcFFebXMM888ww033MCyZcvo2bMnvnPQRKkSMQsMM0s2s2fN7F0zW2RmbUtt/7WZfWBm75jZ8Mi6NDNbYWbzIrdhsapXpFZ56CFYvjyYIW/UqODivFro0ksv5cCDDgLgmhtuxjRVbJWKZQtjEJDn7p2Bu4DROzeY2aHApUAnoBvQ18yOAI4Gst09I3IbF8N6RWo2d8jLC5Z/9zvIyYFLLw23pmo2Y+V6Co85j/r7p/HAhym8svyrsEuqVWIZGH2AFyLLC4AOpbb/xt23u3sxkA80AdoC55nZW2Y2xcwOiFWxIjXa1q0weDCcfDJs2hRcjNeuXdhVVau8LQUMn7aC+sf0oXHGUAqK4M5pK8jbUhB2abVGLAMjFcgD8ODAoptZUuT+5+7+spmlmtlTwCYgB1gN3O/uPYGZwF929+RmNszMcswsZ/369dX8VkTi2H/+A127wrPPBh3cCTK8R+6GfJKTkkhKrk+DQ08AIDkpidwN+SFXVnvEMjC+AxoDWHBg0SOtCSLrzgUWA+8BF0RC5SV3fyOyy1R+2irZxd3HuXu6u6c3b968mt6CSJybORPS02H1avjHP+Cee4JZ8hJAy6YNKCwu/tG6wuJiWjZtEFJFtU8s/ye9CfSLLJ9BcFgKCDq3gVFAT3d/rESQzDSzneMU9CZodYhIWdyDq7YPOSTorzj77LAriqnURvUZ2e84UpKT2Ld+XVKSkxjZ7zhSG9UPu7Raw2J12pmZ1QMmAm2AH4DLCTrClwAtgPuBT0s85BagGHgcKAK2ANe6e25Fr5Wenu45OcoWSRCbNgXDkjdrBt98A/vsE9wSVN6WAnI35NOyaQOFRSWY2RJ3Ty93n9p4nrICQxLGypXB8B5HHBEcghLZQ9EERmIc3BSpjZ5/Hk46CTZuhDvvDLsaSQAKDJGapqgoCIiLL4Zjjw2GJO/RI+yqJAEoMERqmu+/h8mT4frr4a234OCDw65IEkTtHB9ApDb66KNgkqP99oNly4J/RWJILQyRmmDCBOjYEf7wh+C+wkJCoMAQiWcFBXDttXD11UE/xfXXh12RJDAFhki8ys0NQmLcOPjVr4KruNWykBCpD0MkXq1fD198AdOmQd++YVcjohaGJJ68LQV88OX38TmKqTu8+WawfMIJ8PnnCguJG2phSEJ5eflXDJ+2guSkJAqLixnZ7zjO7xAnp6Vu2RL0Vfz97zB7djDXdsOGYVclsotaGJIwds6X8MOGPDblF7CtsDh+5kv497+hSxd47jkYMQJ69w67IpGfUGBIwtg5X8LGxVP5YeVbQJzMl/Dqq9CpE/z3v0HH9vDhoKlFJQ4pMCRhtGzagPxNeWz54HU2LX4O9+L4mC9hyxZo0yYY4uO008KtRaQcCgxJGKmN6nP892/jRdspzPuSos8WhzdfwoYNMGNGsHzppbB4MbRuHfs6RCpBnd6SUG6+8hL+s+JdDjwkjeG3/pyMMDq8V6wIhiT/5ptgZrxmzaCuPooS/9TCkITSs2dPDj7oQDp3aE9G15NiX8CkSXDyybBtG7zxRhAWIjWEAkMkFtzh//4PBg4MOriXLAnOihKpQRQYIrFgFkybesstwTUWBx4YdkUilaYDpyLV6Z13oLgYunWDBx7Q6bJSo6mFIVId3OHxx6Fnz+C6CneFhdR45QaGmR27m/UXVU85IrVAfj4MGQI33BBcVzF9usJCaoWKWhgvm1mbnXfMrL6ZPQH8oXrLEqmhvvsOTjkFnnoK7r03CIumTcOuSqRKVBQY1wKvmllLM2sH5AD1gROqvTKRmqhJEzj66CAo7rsPknTUV2qPcju93X2Wmd0KzCMIitvd/e+xKEykxnCHhx+Giy6CVq1g4sSwKxKpFhV+/XH3V4HbgK+BGdVekUhNsmlTEBS33RbMuy1Si5XbwjCzdwAHDDgM+MDM1gHu7l1jUJ9I/Fq1Khji49//htGjgwvzRGqxiq7D6B+TKkRqmnfegTPOgJSU4EK8jIywKxKpdhX1YawGMLPflNq0w8y+AKa4+45qqk0kfrVvDxdcAL//PbRsGXY1IjER7SkcxxEclsoBmgLHA4cDj1VTXSLx59tv4cYbYetW2HffoHNbYSEJJNrASHX337r7a+7+f0ALd7+fIDREar+lSyE9HZ58Et59N+xqREIRbWDUj1yHgZm1itxvDjSutspE4sXTTwdjQRUXw8KF0KtX2BWJhCLawQdvAP5qZgcDm4GbgbPQFd9S2z34INx5J/TuDVOmQPPmYVckEpqoAsPdl5vZ9UAbYLm7f1q9ZYnEiX79gmst7r1Xs+JJwovqkFTkau/HgVOAZ83s2mqtSiRM8+cHAwe6w2GHBcOSKyxEou7DuBjo5u6/BLoDl+/Ji5lZspk9a2bvmtkiM2tbavttZrbMzJaYWd/IupZmNj/ymBfNbJ89eW2RCrnDI4/AqacG11asXx92RSJxJdrAKHL3YgB3LwLq7OHrDQLy3L0zcBcweucGMzuC4ELBTsCpwEgzq0fQT/Jw5DHLCPpTRKrW1q1wxRXwy1/COefA++/D/vuHXZVIXIk2MGZHvt3fbGbZBNdj7Ik+wAuR5QVAhxLbegP/cPcid/8eWAW0B3oC0yP7vEZwWEykav385zBpEvzud/DCC9BYJwCKlFbRWFJ/4H9jSbUEHgI+BV7cw9dLBfIgGIzKzNzMkiKtl13bIjYCTYC6kVZNyXVl1ToMGAbQqlWrPSxPEtavfw233x4M9yEiZaqoJ2/Vbpb31HdErt0wMyPIjeIS21JL7NsMWA9sN7M6kSFIdq77CXcfB4wDSE9P9yqoVWqhydnZvDdnDjNefZXnx4wh87TTGPDcc7q2QiQKFY0l9XQVv96bQD/gbeAMgsNSO80BsiOtmlQgDfg4su85wCvAhcDrVVyTJIjJ2dlkDhvG3/Pz6Q4s3LiRoS++CNnZDBg4MOzyROKeucfuy3ikE3siwfUcPxCcbTUIWOLuM83sLmAAUAjc5e6zzSwNyCYIt8+Awe5eWN7rpKene07OnnazSG3VPi2NR1evpmRbYi5wU+vWfPTFFyFVJRIfzGyJu6eXu08sAyNWFBhSljpJSWxzJ7nEukIgxYwdxcW7e5hIQogmMDThsNR+kS9F7Vq1YmGpTQsj60WkYgoMqd2++Qb69IGXXiIzK4uhDRsyl6BlMRcY2rAhmVlZIRcpUjMoMKT2evdd6NgRFi2C/HwGDBxI1rhx3NS6NSlm3NS6NVnjxqnDWyRK6sMQ8rYUkLshn5ZNG5DaqH7Y5VSNJ58MJjtq0SK4EO+EE8KuSCSuRdOHoRHVEtzLy79i+LQVJCclUVhczMh+x3F+h4PDLmvvzJ8Pw4YFF+FNmgTNmoVdkUitoENSCSxvSwHDp61gW2ExmwuK2FZYzJ3TVpC3pSDs0vZMYeRs6x49glbFq68qLESqkAIjgeVuyCc5KYnNH7xO8fZ8AJKTksjdkB9yZXtgzhw48kj46KPg/oUXQp09HSNTRMqiwEhgLZs2oKBgK9/Pn8iW5TMBKCwupmXTBiFXVgnuMGoUnHYaNGwI9eqFXZFIraXASGCpjerTs/hDirduZPP7L1LPihjZ77ia0/G9eTNceinccQf07QuLFwetDBGpFgqMBOburPvkferVq0fzxg35bUevWR3ejzwC06bByJEwdSrsu2/YFYnUajqtVjjyyCO57777uOyyy8IuJTpbtkCjRrB9OyxZAl26hF2RSI2noUGkdtmxA+65B447DvLygv4KhYVIzOg6DKkZvvsOLrsMXn8dhgyBfTS1u0isKTAk/i1fHnRq5+bCE0/ANdeAWdhViSQcBYbEv8zMoL9iwQLo3DnsakQSlgJD4lNhYdC53bQpPPUUFBfDAQeEXZVIQlNgSPz573/hkkugbl14801o3jzsikQEBYbEm0WL4KKLYONGmDABknQin0i80KdR4oM7/OUv0LNnMMTH4sXQv3/YVYlICQoMiQ9btgRjQp1xBuTkwLHHhl2RiJSiQ1ISrjVr4MADg2E9Fi6Egw7SYSiROKVPpoTnjTeCmfB+/evg/sEHKyxE4pg+nRJ7xcXw+9/DmWcGIXH99WFXJCJR0CEpia1Nm2DwYHjppWCoj3HjNMyHSA2hFobEVm4uzJ0LY8bAs88qLERqELUwJDbefx/S0+Hoo+Hzz4MruEWkRlELQ6pXURHceSecdFIw2REoLERqKLUwpPqsXx9cfDdnDlx3HZx3XtgVicheUGBI9cjJgX794Ouv4a9/hauuCrsiEdlLCgypHrm5wZwVb78NHTuGXY2IVAH1YUjVKSgIDj8BXHABrFqlsBCpRRQYUjVyc4OBA884A774IliXkhJqSSJStRQYsvfmzQtaEh9/DFOmQFpa2BWJSDVQYMje+dOfoE+f4FTZ994LOrpFpFaKWWCYWWMze83MFpvZbDM7sNR2M7NHzWypmb1rZldE1nePrJsXuZ0fq5olCps3B6fLvvcetGsXdjUiUo1ieZbU7cAsd3/YzAYB9wC/KLG9J9AG6AjUBz40sxeB44A/uvvfY1irlOff/w5Ol+3e/X8jzWqUWZFaL5af8j7AC5Hl14BTSm3fCIxwdwcKAQMaAG2Ba8xsvpmNM7NGsSpYyvCPf0CnTnD11bBjRxAUCguRhBDLT3oqkBdZ3gg0KbnR3Ze5+1tm1gp4BZjh7uuBj4Dh7t4DWAvcV9aTm9kwM8sxs5z169dX01tIYMXFcN99weGnww+HmTOhTp2wqxKRGKqWwDCze8zso5I3goBoHNmlGfCTv+pmdjUwC5jo7jdFVk9w9yWR5alAh7Je093HuXu6u6c3b968Ct+NkJ8P558Pv/0tXHllMDOezoQSSTjVEhju/oC7ty95A54Edp5CcyHwesnHmFlXYBhwUqn+ipVmdnhkuTeQUx01SzlSUqBJE3j88WCYjwYNwq5IREIQy07vUcAUM7sc+Ba4DMDMHgUmAGcStDxeNrOdj+kPXA9MNrMfCFol18Sw5sQ2ZUrQX3H44fDMM8FQHyKSsGIWGO6+ATijjPU7Dz0tB35TxkPXASdVX2XyE4WFwZDkY8bAsGHwxBMKCxHR4INSytdfwyWXwPz5cPPNMGpU2BWJSJzQ+ZAJbnJ2Nt9+/jlXDBxI+xYtmNy2bTA73rPPBldxJyeHXaKIxAm1MBLY5OxsMocNY1pREd2Bhf/9L0Pr1IEHHmDAwIFhlycicUYtjASWlZnJhK1b6QUkA72ACTt2kPXEEyFXJiLxSIGRwD5Zs4bupdZ1j6wXESlNgZGoZs+mnRkLS61eCLRr1SqMikQkzikwEtGzz8IZZ5B50EEMTUlhLsHgXXOBoQ0bkpmVFXKBIhKPFBiJqEcPuOYaBqxaRdb48dzUujUpZtzUujVZ48apw1tEymTB4LC1S3p6uufkaASRH1m1CsaOhYce0uiyIvITZrbE3dPL20d/ORLBiy/CSSfBpEn/m29bRKSSFBi12Y4dkJkJffsGs+EtWQKHHRZ2VSJSQ+nCvdpsyBCYODEYD+qRR6B+/bArEpEaTIFRm11zDZxySjA7nojIXlJg1DZPPw2rV8NvfhPMud299KV5IiJ7Rn0YtcX27fCLXwQz4r31VjBEuYhIFVJg1AZr10JGBjz2GNx+O7z+ukaZFZEqp0NSNd22bdC1K3z7Lfz978FcFiIi1UCBUdOlpMDIkXDMMcFNRKSa6JBUTbR1K1xxRTDnNgStCoWFiFQztTBqmv/8J7gQb8UKOO64sKsRkQSiwKhJZsyAyy4DM3jtNTjzzLArEpEEokNSNcVHH8E550Dr1pCTo7AQkZhTYMS74uLg3/btg2E+Fi3SeFAiEgoFRjz7+GM44QRYujS4f/nl0LBhuDWJSMJSYMSrqVOhc2f45hsoKAi7GhERBUbcKSqCO+6ASy+FDh2C1kWXLmFXJSKiwIg7Tz4Jo0bBjTfCnDlw0EFhVyQiAui02vixfTvUqxcMSd66NZx9dtgViYj8iFoY8eDJJ4MZ8b7+GurWVViISFxSYIRp27agRTFsGLRpE4SFiEicUmCE5csvoUcPGD8efv3r4Mrt1NSwqxIR2S19pQ3LXXfBqlXw4otwwQVhVyMiUiG1MGLJHTZtCpYffRTef19hISI1RsxaGGbWGJgMNAO2AJe7+7pS+0wF9o/c/c7d+5pZS2ASUB9YG3ncD7Gqu8ps2QJDhwaHoubNg2bNgpuISA0RyxbG7cAsdz8ZmAjcU8Y+B7h7RuTWN7LuD8DD7t4ZWAbcEJtyq9C//gUnnwzPPx+0KDR9qojUQLEMjD7AC5Hl14BTSm40s2ZAKzObbWZvmdnO4Vh7AtN397i498or0KkTrFsXzLV9553B8OQiIjVMLDu9U4G8yPJGoEmp7Q2BccCDBIel5ptZV6CuuxeV8zgAzGwYMAygVatWVVn3nissDALiiCNg2rTggjwRkRqqWgLDzO4BLi21ugnQmKD/ohmwvtT2r4AR7u7AWjNbChwJbDezOu6+YzePA8DdxxEEDunp6V5Fb2XPbNgADRoE823PnAkHHBDcFxGpwarlkJS7P+Du7UvegCeBfpFdLgReL/Wwy4CnAMysEXAM8AnwNnBOOY+LLx98AB07wi9/GdxPS1NYiEitEMs+jFHAOWb2HnA+MBLAzB41sw7A1Mj9xQSh8Bt3zwMygeFm9i7QCvhbDGuunGefDUaW3b4drroq7GpERKpUzPow3H0DcEYZ628qcXdwGdu/ALpVX2VVoLAQbrstuLaiR49gLosDDgi7KhGRKqUL96rCl1/C00/DLbfA7NkKCxGplTQ0yN745z/hyCODObY/+QRatAi7IhGRaqMWxp5wh8ceg2OPhQkTgnUKCxGp5RQYlZWfD1deCb/4BZx+Olx0UdgViYjEhAKjMj7/HLp1g4kT4b77gqu4mzQJuyoRkZhQH0Zl/POfsHo1TJ8O554bdjUiIjGlwKiIO7z7bjB44JlnBq2Mn/0s7KpERGJOh6TKs2kT9O0bHIb64INgncJCRBKUAqOUydnZtE9Lo05SEu1TU5n88sswejQcd1zYpYmIhEqHpEqYnJ1N5rBhTNi6le7AwqIihtavD82bM0BDkotIglMLo4SszEwmbN1KLyAZ6AVMKCggKzMz5MpERMKnwCjhkzVr6F5qXffIehGRRKfAKKFdq1YsLLVuYWS9iEiiU2CUkJmVxdCGDZkLFAJzgaENG5KZlRVyZSIi4VOndwkDBg7EgfMvv5ytZrRr1YqsrCwGDBwYdmkiIqGzYEbU2iU9Pd1zcnL26LHuTlJSEjk5OXTs2LGKKxMRiU9mtsTd08vbR4ekREQkKgoMERGJigJDRESiosAQEZGoKDBERCQqCgwREYmKAkNERKKiwBARkagoMEREJCoKDBERiYoCQ0REoqLAEBGRqCgwREQkKgoMERGJigJDRESiosAQEZGoKDBERCQqCgwREYlKzOb0NrPGwGSgGbAFuNzd15XYfiVwZYmH7A/cAXwMvAJ8F1k/yd3HxaBkEREpIZYtjNuBWe5+MjARuKfkRnd/yt0z3D0DuA5YA7wOHA1k79ymsBARCUcsA6MP8EJk+TXglLJ2MjMDxgH/5+5FQFvgPDN7y8ymmNkBMalWRER+JJaBkQrkRZY3Ak12s9/PgVXu/knk/mrgfnfvCcwE/lLWg8xsmJnlmFnO+vXrq65qEREBqikwzOweM/uo5I0gIBpHdmkG7O6v+u3AoyXuv+Tub0SWpwIdynqQu49z93R3T2/evPlevwcREfmxagkMd3/A3duXvAFPAv0iu1xI0D/xI2bWCtjX3T8ssXqmmfWKLPcGcqqjZhERKV/MzpICRgFTzOxy4FvgMgAzexSY4O7LCfo55pZ63G3A42ZWRHB21bUxq1hERHaJWWC4+wbgjDLW31Ri+a9lbF8BdKve6kREpCK6cE9ERKKiwBARkagoMEREJCoKDBERiYoCQ0REoqLAKGHLli307NkTgMGDB/Pxxx+HXJGISPxQYJTQqFEjjj76aAAsOYUDWh0eckUiIvFDgVHKiecNBkti01Hn0n3kXF5Z/lXYJYmIxAUFRgl5Wwp4aPH3NOl5JbTuyLbCYu6ctoK8LQVhlyYiEjoFRgm5G/JJTkqicee+mAU/muSkJHI35IdcmYhI+BQYJbRs2oDC4uIfrSssLqZl0wYhVSQiEj8UGCWkNqrPyH7HkZKcxL7165KSnMTIfseR2qh+2KWJiIQulqPV1gjndziYbm32I3dDPi2bNlBYiIhEKDDKkNqovoJCRKQUHZISEZGoKDBERCQqCgwREYmKAkNERKKiwBARkagoMEREJCoKDBERiYoCQ0REoqLAEBGRqCgwREQkKubuYddQ5cxsPbA67DoqsB/wbdhF7AXVHy7VH67aWH9rd29e3oNqZWDUBGaW4+7pYdexp1R/uFR/uBK1fh2SEhGRqCgwREQkKgqM8IwLu4C9pPrDpfrDlZD1qw9DRESiohaGiIhERYERI2Z2qZmNKGP9mWa2zMzeM7PfhlFbNMqpf4iZLTKzHDP7ZRi17Y6Z1Tezv0d+tovN7PRS22+L/OyXmFnfsOrcnSjq/7WZfWBm75jZ8LDqLE9F76HEfiPL+v8Vtih+B3H9+Y2i/sp9ft1dt2q8EYTyLGAbMKLUtgbAF0BroA6QAxwads2VqL8h8B8gOXL7FNgv7JpL1Hcl8HhkuTnw7xLbjgDeJ5imuEmk9nph11yJ+g8FPgDqRX5H7wJHhF1zZd5DiX1OBNaX/v8VD7cKfgc14fNbXv2V/vxqTu9q5u7FZnYWMAg4stTmk4Bl7r4awMz6Ad/FuMRyVVC/RW77AMUEf7i2x7bCcq0GlkWW84FGZmYefFp6A/9w9yLgezNbBbQHloZTapnKqx/gN+6+HcDM8gmCL96U+x7MrC7wIDAKaBpOieUqr/64//xSfv2V/vwqMGLA3YvMrLiMTS0AM7PngAOBucC9MS0uCrur391/MLMZwCeRVW+6+6bYVrd77j4XwMzaA08Co0v8sU0F8krsvpE4+4NbXv3u/jnwuZmlAqOBTQTfcONKBb8DgNuAKUAhcRgYFdQf95/fCv4PVfrzqz6McG0maM4OBjIImubnhFlQZZhZT6AtkEbwPvaLtEbihpn9BpgE3O3uo0ps+g5oXOJ+M4LDInGlnPoxs3OBxcB7wAWl/hDHjd29BzNrA5wKjA+rtmiU8zuoEZ/fcn7+lf78qoURrqUE32y3RQ79bCL4plVT1Ac2unsB7BrDy8It6X/MbADQCei0s8YS5gDZZvYHgtZGGvBxbCssX3n1m1kawWGcnu6+NoTyolLB7+AUgp/9XIJv6A3NLNfd/xzjMnergvrj/vNbQf2V/vwqMEJgZpnAEnefaWYTgPlmtgN4z91fD7m8Cu2sH3gdON3MFgE7CD5AM8KsrZSzCDqHXzfb9TmYxf9+9tMIju8WAje5e1mHDcO02/oJDoc0AiaV2HaLuy+PcY0VKe938DfgbwBmdiXQNp7CIqKi/0Px/vkt7/9QpT+/unBPRESioj4MERGJigJDRESiosAQEZGoKDBERCQqCgwREYmKAkNkD5hZmplNM7N5kcH/Jkeuui65T4aZTanEc15gZu+a2cLIgHA3RtZfGY8D80ni0XUYIpVkZinAa8AN7j4vsm4o8ARw0R4+5wnACKC3u681s/rADDNbXTVVi+w9BYZI5Z0NfLQzLCL+CrxoZscDYwkGccuL/IuZ5QEvEgwfsQxIAQ6PPM/VwHXAn3Zete3uBZHB7OoA58biTYlURIEhUnmHEwwLjZk1A16IrG9JMB7Vte6+wszuAjpEtjUFHgC+BDYAJ7v7J2b2LzNrSjA0yYslX8TdN0Reo1rfjEi0FBgilZcLnADg7t8BGRb8Vf8BaObuKyL7LeZ/gfFtiWGw89195wihWwnG9PkvQeDsYmZ9gH2r8X2IVIo6vUUq73Wgq5l1LLHuFoIJdT43s2Mi67qV2F7ROFVPA7eYWXMAM2sMjAS+qZKKRaqAWhgileTu35nZxcBIM2tAMHDbW8BMIAv4m5l9B3xWieeca2b3A9PNbOckNr9z97fN7Igqfgsie0SDD4qISFR0SEpERKKiwBARkagoMEREJCoKDBERiYoCQ0REoqLAEBGRqCgwREQkKgoMERGJyv8DAJXCcShaJVYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize=(6, 6))\n",
    "ax = fig.subplots(1, 1)\n",
    "# observations\n",
    "transition.plot.scatter(\n",
    "    \"GdmCl\",\n",
    "    \"logK\",\n",
    "    label=\"observations\",\n",
    "    zorder=1,\n",
    "    ax=ax,\n",
    ")\n",
    "# predictions\n",
    "transition.plot(\n",
    "    \"GdmCl\",\n",
    "    \"logK_pred\",\n",
    "    color=\"red\",\n",
    "    label=\"prediction\",\n",
    "    markeredgecolor=\"k\",\n",
    "    linestyle=\"--\",\n",
    "    marker=\"o\",\n",
    "    zorder=3,\n",
    "    ax=ax\n",
    ")\n",
    "# residual arrows\n",
    "for i in range(transition.shape[0]):\n",
    "    ax.arrow(\n",
    "        x=transition.iloc[i].GdmCl,\n",
    "        y=transition.iloc[i].logK_pred,\n",
    "        dx=0,\n",
    "        dy=transition.iloc[i].logK_res,\n",
    "        width=0.005,\n",
    "        length_includes_head=True,\n",
    "        facecolor=\"k\",\n",
    "        edgecolor=\"w\",\n",
    "        lw=4e-5,\n",
    "        zorder=2,\n",
    "    )\n",
    "ax.legend()\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nteract": {
     "transient": {
      "deleting": false
     }
    }
   },
   "source": [
    "Residuals are often associated with the **unexplained** variation in the data.\n",
    "For this particular case, this means that while the model explains\n",
    "something about how $\\log{K}$ changes as a function of [GdmCl], we can't predict\n",
    "with absolute certainty what the exact value of $\\log{K}$ will be for any one\n",
    "particular concentration of denaturant, we can only provide an approximation.\n",
    "\n",
    "A probabilistic interpretation of this situation is that, in repeating this\n",
    "experiment, we will be willing to accept similar deviations between model and\n",
    "observations, residuals of similar size. Consequently, different experiments\n",
    "could have given us sligthly different datasets, yielding different parameter\n",
    "estimates. It follows that our current model is only one of many possible\n",
    "lines, so our parameter estimates carry some **uncertainty**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fitted equation:\n",
      "log K = -2.29 + [GdmCl] • 1.14\n"
     ]
    }
   ],
   "source": [
    "print(\n",
    "    f\"Fitted equation:\\n\"\n",
    "    f\"log K = {INTERCEPT.round(2)} + [GdmCl] • {SLOPE.round(2)}\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What residuals can be accepted as typical we need to asses depending on the\n",
    "residuals we've observed so far. By definition, residuals have mean 0 and their\n",
    "standard deviation can be calculated as follows:\n",
    "\n",
    "$$\\sigma_{res} = \\sqrt{\\frac{\\sum{(y_{res} - \\bar{y}_{res})^{2}}}{\\mathrm{df}_{res}}},$$\n",
    "\n",
    "where $\\mathrm{df}_{res}$ is the number of degrees of freedom of the residuals.\n",
    "For $n_{obs}$ observations and a linear model with one variable,\n",
    "$\\mathrm{df}_{res} = n_{obs} - 2$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "jupyter": {
     "outputs_hidden": false,
     "source_hidden": false
    },
    "nteract": {
     "transient": {
      "deleting": false
     }
    }
   },
   "outputs": [],
   "source": [
    "dfres = transition.shape[0] - 2\n",
    "sigma = np.sqrt(transition.sqres.sum() / dfres)  # sigma = sqrt(MSE) = sqrt(SSE/df)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nteract": {
     "transient": {
      "deleting": false
     }
    }
   },
   "source": [
    "Rather than actually repeating the experiment, we can simulate what measurements\n",
    "we are likely to obtain by adding appropriately *typical* residuals to our\n",
    "current model. A common assumption is that the residuals follow a Gaussian, so with\n",
    "mean 0 and $\\sigma_{res}$ we can generate datasets that would be indistinguishable\n",
    "from the one we measured."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "jupyter": {
     "outputs_hidden": false,
     "source_hidden": false
    },
    "nteract": {
     "transient": {
      "deleting": false
     }
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAFUCAYAAAAgdpvBAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAA/8ElEQVR4nO3deXhU5dn48e8zk41FAoSgIEtAAUERkDWRkomogFurVAVpg0gL4vKqrS+1RavVFxfenxW1b8UoVaNIqwXR2qJgYFjjwioaEAUDiZQawiaQzCQzz++Pk4RMSA5JmPWc+3Nduch9zsycOyGZO89ynkdprRFCCCGqOSKdgBBCiOgihUEIIUQAKQxCCCECSGEQQggRQAqDEEKIAFIYhBBCBIiLdAJnYuzYsfqDDz6IdBpCCBFrlNnJmG4xHDhwINIpCCGE5cR0YRBCCBF8UhiEEEIEkMIghBAiQEwPPtenoqKC4uJiysvLI52KLSQlJdGlSxfi4+MjnYoQIkgiUhiUUjcDg7TWD9Q5Phe4BPBXHbpSa+1tymsXFxdz1llnkZaWhlKmA+/iDGmtKS0tpbi4mB49ekQ6HSFEkIS1K0kp5VBKLQdea+AhAwCX1rr6o0lFAaC8vJyUlBQpCmGglCIlJUVaZ0JYTFgLg9baD4wD7mjgIV2Bfyql1iqlspt7HSkK4SPfayGsJ+yDz1rrSk52FdVQSrUA/gaMB64G7lNKDajncdOUUhuUUhtKSkpCnm8wpKWlReSv6iVLlvDNN98AcMMNN4T9+kKI2BRNs5LKgYe01ie01keA5UD/ug/SWudorYdorYekpqYG5cKlxzxsLTpM6TFPUF4vVPz+U+qpqdqFYfHixaFISQhhQdFUGDKAPGWIAy4FNoX6ou9u+Y5Ln1rBz17+hEufWsF7W747o9f74YcfmDhxIpmZmWRkZLB8+XIA7rrrLrKyssjKyuK7775j//79XH755WRmZpKZmUlBQQEAjz32GJmZmQwdOpQFCxYA4HK5uO2227jtttvIyspi1apVAKxevZof//jHHDt2jOuuuw6Xy8XgwYP561//yqpVq/jggw+YOXMm+/fv55xzzgFg/fr1jBo1ilGjRnH99ddz8OBB3G43mZmZXHfddQwePJgHHjDmBLz99tsMGzaMzMxMpkyZ0uTCJISIUVrrsH8AtwJPVn0+Cxhb9fkcYAOwFrjrdK8zePBgXVdBQcEpxxpy4Idy3efBf+nuv3m/5qPPg//SB34ob/Rr1DVr1iz99NNPa621/s9//qM7d+6su3btqteuXau11nrevHl6xowZ+tVXX9VTpkzRPp9Pb9y4US9dulR/9NFHesKECVprrcvKyvQFF1ygS0tLdWZmpn7jjTe01lrn5ubqGTNmaK21/sUvfqH/8Y9/6M8//1y/+eabWmut8/Pz9WWXXaa11nry5Ml66dKlWmutzz77bK211r1799bFxcVaa62fe+45/V//9V965cqVuk+fPrqsrEx7PB7dtm1b7fP5tMvl0itWrNBaa/3yyy/r/fv31/s1N+V7LoRovvV71+vHVz+u1+9df6YvZfreGpHpqlrrV2t9PrvW5zPDmUfxoTLiHQ7Kaw15xDscFB8qI6V1YrNec/PmzcyebXxJHTt2JCUlhR9++IGhQ4cCMHz4cN59913mzp3L/v37mTBhAk6nk7vvvpt169axYcMGXC4XYAzs7t27F4DBgwcD8NOf/pTHHnuMEydO8PHHHzNv3jy+++47li9fzqpVq/B4PPh8vnpzKykpoUWLFpx77rkAXHrppVQvQjhgwACSkpIASE5Oxuv1Mm/ePJ577jleeOEFzj77bG666aZmfU+EEGcuvyif0bmj8fq8JDgTyMvOI71rekiuFU1dSWHXpV0LKup0j1T4/XRp16LZr3n++efz8ccfA7Bv3z48Hg9aa9asWQPAmjVrGDhwIH//+9/p378/b731FnfddRePPvooffr0YcyYMbjdbpYvX85PfvITLrjggoDXb9GiBZdddhm/+tWvuP7663E6nTz99NO4XC7mzZvHlVdeWfPYujOG2rdvz+HDh6ketF+1ahWXXHJJvY8F+OMf/8j//u//8tZbb1FeXs4///nPZn9fhBBnxl3oxuvz4tM+vD4v7kJ3yK5luTufmyKldSJzxl/MzEWfE+9wUOH3M2f8xc1uLQA8+OCD3HrrrSxcuBC/38/LL7/MlClTWLBgAQ8//DDt2rUjNzeXo0eP8vOf/5wnn3yS48eP88QTT3DllVfidrtxuVx4vV6mTp1a81d8bVOnTmXEiBHs2rULgJ/85CfMmDGDV155hXHjxrFv3z7WrVvHsGHDmDlzZk1rw+l08uc//5lrrrmGli1b0rFjR1566SU2bap/KOfSSy9l1KhRJCQk0K5dO8aNG9fs74sQ4sy40lwkOBNqWgyuNFfIrqW00bcfk4YMGaI3bNgQcGz79u307du3Sa9TesxD8aEyurRrcUZFwa6a8z0XQjRdflE+7kI3rjTXmXYjmd6AZOsWQ7WU1olSEIQQUS+9a3rIxhVqs/UYgxBCxJL8onyeWPME+UX5Ib2OtBiEECIGyKwkIYQQAdyFbjw+Dz7tw+PzhHRWkhQGIYSIASktU/BrY3q9X/tJaZkSsmtJYRBCiBhQeqIUR9VbtgMHpSdKQ3Yt2xeGhQsWcFFaGk6Hg4vS0lhYtT5RpD3yyCPMmzePLVu28Oijjzb4uLlz51JZWXnaxwkhYpsrzUViXCJO5SQxLjGk9zHYevB54YIFzJo2jfknTjASWLtnD1OnTQNg4qRJIbmm3+/H4Wh8PR44cCADBw5s8PzcuXO5/fbbT/s4IURsS++aTl52XrDuYzBl/RaDy3Xqx5//DMDs3/6W+SdOkAXEA1nA/BMnmH3PPcZzDxw49bmN8OqrrzJu3Dguv/xy+vfvzwsvvBCwQqrP5+POO+8kKyuL4cOH89FHHwGwYMECLrnkElwuF6tXrwbA7XYzYcIEAJ5//nmGDx/O0KFDmTt3Ls8++yz79+8nOzs74HE5OTlkZGSQkZHB/fffDxgtkBtvvJFx48Zx4YUXyjLcQogG2brFsL24mJF1jo0Etpeeed/d0aNHWbNmDWVlZfTv35+EhAR++ctfMmnSJObPn0+bNm1YuXIlhw4dYujQoXzyySc89thjbNy4kZYtW3L11VcHvN6OHTt4/fXXWb9+PX6/n+zsbN58802eeeYZcnNza9Zn2rlzJy+//DLr1q0jPj6e8ePH88477wDQqlUr3n77bfLz85k1a5Zs3iNEDMkvyifrtaya6aorJ68MWavB+oXB7W7wVN9u3Vi7Zw9ZtY6tBfp2724EHTqYPt9MRkYGDoeDVq1acdFFF7FixYqaNYs2bdrE6tWryc83blJxOBzs2rWLCy64gFatWgEwbNiwgNf7/PPPGTlyJHFxxn/ZX//613qvu3XrVjIyMoiPjweM9Y6qN+sZPnw4AJ06deLEiRPN+rqEEJGRuzUXj8/YTMzj85C7NVfuYwiFWbNnM7VlS1YCFcBKYGrLlsyaPfs0zzy9zZs3o7Xm+PHjFBQU0KtXr5pzffr0YfLkybjdbt577z1uvPFGevToQUFBAWVlZWitWb9+fcDr9erVq+Y1jxw5wmWXXYbf7z9lVdTzzz+fDRs21Cy9vXr1atMVVIUQsaGgpMA0DiZbF4aJkyYxOyeHu7t3J0kp7u7endk5OUEZeK6srCQrK4v09HQefPBBkpOTa85Nnz6dbdu2kZWVxdVXX83w4cNJTU3lt7/9Lenp6Vx77bWkpATOUR40aBCXXXYZGRkZXHHFFfzyl7/E4XCQnp5OdnZ2wOOuuuoq0tPTGTVqFBdddBGjR48+469HCBFZ5ZXlpnEwyeqqIfDqq6+yY8cOnnzyyYjmES7R8D0XwupyNuYw/f3pNfGL17zItMHTmvtypt0Htm4xCCFErOjfsT/xDmPsMN4RT/+O/UN2LesPPkfArbfeGukUhBAW4y50ByyJ4S50y+CzEELYWfUObk7lDPkObpZsMWitZQZOmMTyGJUQsUTufD4DSUlJlJaWyhtWGGitKS0trXdfaiFE7LJci6FLly4UFxdTUlIS6VRsISkpiS5dukQ6DSEsL5wb9ViuMMTHx9OjR49IpyGEEEHlLnTj9XnxaR9en1cGn4UQwu5k8FkIIUSAcA4+S2EQQsSs/KL8sLxRRov0rulh+TqlMAghYlI4B2PtRsYYhBAxqb7BWKvL2ZjDmNfHkLMxJ6TXkRaDECImVQ/GVrcYQjkYGw1qL6K3bPcygDNZRM+UtBiEEDGpejD2sazHbNGNNH/TfNM4mKTFIISIWeEajI0GSXFJpnEwSYtBCCFiQd3l30K4HJwUBiGEiAEHjh8wjYNJCoMQQsSA3im9TeNgksIghBAxYFyvcaZxMElhEEKIGLD535tN42CSwiCEEDGg4ECBaRxMUhiEECIGlBwvMY2DSQqDEELEgD4pfUzjYJLCIIQQMUAGn4UQQgQoPVGKqrqrzYGD0hOlIbuWFAYhhIgBKS1T0GgA/PhJaZkSsmtJYRBCxKz8onyeWPME+UX5kU4l5Gq3GBQqpC0GWURPCBGT8ovycb3mosJXQbwzHvdkt6UX1KvdYtBo67UYlFI3K6WerOf4r5VSm5VSG5VSN0QiNyFEbMjdmovX50Wj8fq85G7NjXRKIbX066WmcTCFtTAopRxKqeXAa/Wc6wVMAIYCo4E5SqmEcOYnhBDR6pPvPjGNgymshUFr7QfGAXfUc/oy4H2tdaXW+jCwA7gojOkJIWJI9oBsEp2JKBSJzkSyB2RHOqWQqvRXmsbBFPauJK11JeCv51QKUHs05QjQtu6DlFLTlFIblFIbSkpCd+efECK6pXdNZ+Xklcy+bDYrJ6+09PgCwJXnXWkaB1M0zUo6CCTXitsDp7zza61ztNZDtNZDUlNTw5acEEJEUtHRItM4mKJpVtIKYIFS6gmM1kMa8GVEMxJCRK38onxG547G6/OS4Eyw/L7P20u2m8bBFPEWg1JqllJqrNZ6J7AI2AwsBe6uGpMQQohTuAvdlFeW49M+yivLcRe6I51SSCUnJpvGwRSRFoPW+tVan8+u9fmTwCnTWIUQoq7DnsMB8/oPew5HNqEQK69QoDH2etZVcYhEvMUghBDNseXfW0xjKyk95uHwkQ4Bxw4f6UDpMU9IrieFQQgRkwZ2GmgaW0nxoTKcqgXsBLzGMadqQfGhspBcTwqDECIm7Tyw0zS2ki7tWnB8/zZ4C/jAOHaCHXRp1yIk15PCIISISft+2GcaW0m7lvF03NnaCIqB/ZDRdQQprRNDcj0pDEKImNQrpZdpbCUOh4Nhtw+A1kAm0AnanRW6SZtSGIQQMenr0q9NY6v56sBXpnEwSWEQQsSkpPgk09hqZM9nIYQ4jX4d+pnGVtO7Q2/TOJikMAghYlL2gGwcVW9hDhyWX101nPdtSGEQQsSkJV8twV+1ULMfP0u+WhLZhEIstVWqaRxMUhiEEDFpccFi09hqti3fSuujoN6G1s8YcahIYRBCxKThXYabxlaycMECSl/ZznsaPMB7R6D0le0sXLAgJNeTwiCEiEkXpl5oGlvJ7FmzeN3jIwuIB7KA1z0+Zs+aFZLrSWEQQsQkV5qLBGcCCkWCMwFXmivSKYXM9r17GVnn2Miq46EghUEIEbN8fh8ajc/vi3Qqwef3g9dYMa9v+/asrXN6LdC3W7eQXFoKgxAiJs1ZNwefNgqCT/uYs25OhDMKEr8fFi+GQYNg7lwAZj3zDFNbtmQlUAGsBKa2bMms2bNNXqj5pDAIIWLSzoM7TeOYU7sgjB8P5eVw/vkATPz5z5mdk8O0c88mScG0c89mdk4OEydNCkkq0bTnsxBCNFqHFh1M45gzYwbk5EDv3vD66zBhAsSdfItOG9WTvbcfQvsUe52HSBvVM2SpSItBCBGT+qX2M42jXnUL4bvvjPi224yC8OWX8LOfBRQFgNytuXh9XjQar89L7tbckKUmhUEIEZPaJLUxjaNW3S6jl14yjg8fXm9BqFZwoMA0DiYpDEKImPTWF2+ZxlFpyZLAMYQ33oCHHmrUU8O5Y52MMQghYtLxiuOmcdTQGpQyPl+8+GRBmDABnM5Gv4xf+03jYJIWgxAiJvVo28M0jrjaXUabNxvHnnsOCgpg0qQmFQWAtLZppnEwSWEQQsSkqN3zue4YQlkZHDlinGvbtskFodq3h781jYNJCoMQIiZ1aNXBNI4IrSEz8+QYQvUsI5frjF/a6XCaxsEkhUEIEZPaJLQxjcPG74cPPzw5lnDzzabTTptrxLkjTONgksIghIhJe4/uNY1DrnaX0dixsGyZcfyuu4JaEKqFs+tMCoMQIia1iGthGodMfUtXvPEGXH55SC+77fttpnEwSWEQQsSka/tcaxqHTEUF3HPPyYLQzFlGTZWcmGwaB5PcxyCEiElHy4+axkHj9xs3ps2fb7QUEhNhxQro2TPkxaC2ExUnTONgkhaDECImhXyJiLpdRt98A3v2GOd69QprUQDo2a6naRxMUhiEEDGp5HiJaXxG/vOfU8cQCgqMlU8jZMbQGaZxMElhEELEpD4d+pjGTeb3G2/+AB07Qr9+YR1DOJ2lXy81jYNJCoMQIiaNO3+cadxotbuMhg+H0lLjfoSFC6OiIFT7+LuPTeNgksIghIhJpSdKURiL0ykUpSdKm/YC9U07feEFSA7dbJ8zEbbBdmRWkhAiRqW0TEGjAdBoUlqmNO0FCgqMgtDAjmnRRlWv0NpAHEzSYhBCxKSl3yw1jU9R3UJ49FEjvugiyMsL+tIVodKpdSfTOJikMAghYtLHRR+bxjXqdhn97W/g8RjnLrss6gtCtYNlB03jYJLCIISISeW+ctMYMPZBqDvt9PPPjZvUYsyximOmcTBJYRBCxKS05LT6Y78fDhwwPu/Y0fg3iqadNlf7pPamcTBJYRBCxKS4L520fgbUI9D6GYjf5jjZZXTTTcaDzj0XtmyJ6YJQLXtgtmkcTFIYhBAxZ+GCBRS/uIn3joAHeO8IlLywiYXVXUZTpxr7I8DJ/ZZjnExXFUIIE7NnzeLNCk1WVZwF/EXD3SkpTCwoiPnWQaRJi0EIEVv8frbv3cvIOodHAtsPHrRsUfjB+4NpHExSGIQQsaHWtNO+WrO2zum1QN9u3SKRWVh8UvyJaRxMYSsMSql4pdQbSqlPlFLrlVIX1Dk/Vym1WinlrvpICFduQogot3x5wLTTWTNmMClesRKoAFYCkxIUs2bPjnCioXNDvxtM42AK5xhDNlCqtf6ZUmoU8DRwda3zAwCX1tofxpyEENHK74fKSkhIgJKSk/chTJjARKeTW//zMtctruA40CoZKq6IY+KkSZHOOmSeuvwpvjv6HUu/Xsq4XuN46vKnQnatcHYlXQ4srvp8DTCwzvmuwD+VUmuVUqGbhyWEiG6171R+5hnj2M03n3IfQurwjhwD9O1w7D5ITT87cjmHQX5RPou3L+aI5wiLty8mvyg/ZNcKZ2FIAUoBtNYa0EopB4BSqgXwN2A8RiviPqXUgPpeRCk1TSm1QSm1oaQkiBtz2EB+UT5PrHkipD9QQjRbfaud9uplnHM6TxlU7nRW4FpBnc/qHK5MI8Jd6Mbr8+LTPrw+L+5Cd8iuFc7CcBBIBlDGsoC6VrdROfCQ1vqE1voIsBzoX9+LaK1ztNZDtNZDUlNTw5G3JeQX5TM6dzQPrXyI0bmjbVEccjbmMOb1MeRszIl0KqIx7rjj1B3Tbmi4H/2YN3BJiFDO0okGrjQXDuNvaRzKgSvNFbJrhbMw5GG0CADGYHQnVcsA8pQhDrgU2BTG3CwvnH9tRIOcjTlMf386y3YvY/r706U4RKPqFkJxsRFPmdKkpSt2H9ptGlvNtu+3UeGvAKDCX8G277eF7FrhLAyvAZ2VUhuA3wK/UUrNUkqN1VqvAz6r+nADC7XWQd7Z295caS4SnAk4lZMEZ0JI/9qIBvM3zTeNRQTV7TJ66SXj+PDhTVq6oltyN9PYahYVLDKNgylss5K01l5gQp3Ds2udnxmuXOwovWs6edl5uAvduNJcpHdNj3RKIZUUn2QaiwhZsgQefthY4bR375pZRs1xX/p93MmdNfH9GfcHKcnoNL7feJbtXhYQh4osiWEj6V3TLV8QqvXr0I/Ve1YHxCJCtD65XtE77wRMOz2Tu5Q/3PVhQLz0m6VMGzztTDKNatVf26KCRYzvNz6kX6vc+SwsKXtANvGOeADiHfFkD5AZ0GFXu8to82bj2LPPBm35631H95nGVjRt8DQ+/PmHIS+AUhiEZdXeD1iEUd0xhLIyOHLEONe2bdDWMspMywyIXT1cQXndaBauKefSlSQsKXdrLpX+SgAq/ZXkbs21TTdaRGkNmZmwdm3gGEIIFrb7+uDXAfHO0p1Bv0Y0yS/Kx/WaiwpfBfHOeNyT3SH7mZYWgxDizPj98OGHJ8cSJkwIy45pOw8EFoKvDnwVkutEi9ytuXh9XjQar89L7tbckF1LCoOwpEGdBpnGIghqdxmNHWsUB4A77wzLjmkdWnUIiFNbyQ2vwSKFQVhS6YlSFMZMGIWi9ERphDOykPqWrnjjDbjiirCm0b5Fe9PYasI5oULGGIQlpbRMCRh8TmmZEuGMLKSyEu69F1q0COkYwukcPHHQNLYipRQKhQrxdqXSYhCWVLvF4MAhLYYzUd1CuOoq8HiMZbBXrAj5GMLplFeWm8ZW4y504/P70Gh8fp9lFtETImxqtxj8+KXF0Bx1u4x27YI9e4xz558f8S002ya1NY2txpXmwulwolA4HU7LLKInRNjUbSFIi6GJvv/+1DGEggJjCmqU+PbQtwGx1RfRA/BrPxqNP8T7mUlhEJZ02HPYNBb18PuNN3+A1FS48MKwTDttrh7tegTEPdv1jFAm4TFn/ZyAe3PmrJ8TsmtJYbARO23Us+XfW0xjUUvtLqPhw6G01Lgf4c03o7Ig2FU4lwCRwmATdtuop+7Kk6FciTJm1TftdN48Y9mKGDDgnMBNHgd2GhiZRMJk6iVTTeNgkumqNlHfRj2yRITNbd9uFIQQL10RKsmJyTWfO3DQNrFt5JIJg3CuriqFwSaqN+rx+ry22Kinvk1NrLwkc6P4/cay19u2wSOPGGMIK1fCj34UUwWhWmZ3YxE9h3KQGJdo+Z9pMIpDOH6OpSvJJqo36nks6zHysvMs31qo261g9W4GU34/LFoEAwfCT38Kb79tdBsBuFwxWRQARnQdAcC5Z53L3cPvtvzPdDhJYRCWdLT8qGlsG5s2nSwIHo/RZfT555AU+zvavbzpZQCKjhYxZ90c2dc7iKQrySaqB5+ru5Ks3mrYf2y/aWxpfj8cPAgdOsDZZ4PDEZNjCKfzzvZ3AmLpLgweaTHYRH2Dz5ZWdymZ0C4tEx1qzzK68Ubj2LnnGrunWWzaaVlZGZ8+XrUr3FvAUfB4Wkc0JyuRwmAT1YPPTuW0xeDzOa3OMY0tpb5pp7/4hbE/Apzcb9lCTvgceFwdwAmMBs6Cz4oKKT3miXRqliCFwSbSu6Yzd+xcRvcYzdyxcy3djQTQJqmNaWwp8+efunTFpEmWLAjVig+VkdT5XPgF0M84Fk9big+VRTQvq5AxBpvIL8rn3g/uxevzsmbvGvp37G/p4uD+1m0axzS/H5YsMQaQr7oKbrkFWrWCm2+2VHeRmS7tWlBOMZyD0U2ooZxiurRrEenULEFaDDZhtzGGzm06m8YxqW6X0QsvGMdbtTKKg02KAkBK60TO69Au4Nh5HdqR0joxQhlZixQGm3CluYhzxKFQxDniLD/GMO78caZxzFm+/NTVTpcsiXRWETWqx6CTkwpUVWxxORtzGPP6mJBPzbV1YbDTonIAPl21yYf2RTqVkFvw+QLTOCb4/eD1Gp8fOHDqGIKNWgj1yR6QTaIzEYUi0ZkY0q0uo0HOxhymvz+dZbuXMf396SEtDrYtDHZbVC53a27Akr25W3MjnFFobT+w3TSOarXvVP7jH41jN90kBaGO9K7pPDfuOa7oeQXPjXvO0mNmUP8yL6Fi28Jgtz53u+mb2tc0jkp1l67weE5ujON0SkGoI78on7uX3s3y3cu5e+ndlv/jLpwrBtu2MNhtXr/dmt1Pjn4SJ8YbqRMnT45+MsIZNcIddwQuXVFQADfcEOmsolbu1ly8Pi8ajdfntXwreNrgabx4zYtc2fNKXrzmRVldNRSqF5VzF7pxpbks3wxN75rOyskrbfP1rv66BF/V/V0+YM3XJaR3jWhKp6pe7XT4cOjSBW67zVjp1GJLV4jgCdfqqrYtDGC8WVr9DbI2u3y9pcc8POGeBw5f1Rx3H4+75zF12JjomM5YXRD+8AdjCeyHHoJHH4Vhw4wP0SiDOg0yjUXz2bYryY7sMgur+FAZfnUY/EDV7od+dTg67op9551TVzt9+OFIZxWTSk+Umsai+WzdYrCT/KJ8Ml/NpMJfQbwjnlW3rrJs66FLuxY4dDtYBXwKTAdHm3aRuytW65PLU7z33smCIF1GZ+TLki9NY9F80mKwiTnr5lDhrwCgwl/BnHVzIpxR6KS0TmRwYivYBHiAFXBZr+7h70aqPcto0ybj2Ny5Mu00SFYVrjKNRfNJYbCJfT/sM42tJq7bERgHtAHGw3F2he/idaedlpfD0aqNgpKTpSAESc/2PU1j0XymhUEp1b+B4z8NTToiVFw9XKax1YRzzncArSEz89Rppy5XeK5vIyO6jDCNRfOdrsXwrlLq/OpAKZWolHoReCK0aYlga5vYFlW1sIxC0TaxbWQTCrH+HfvXfL0OHPTvWO/fOMHh98MHH5wcS5g4UZauCANLr6AbYacrDNOBfyqluiil+gIbgERA5oXFGFeai6S4JJzKSVJckuVv6Hsg7wE0xo0Mfvw8kPdA8C9Su8to3DijOIBxo5oUhJDrfFZn01g0n2lh0FovB34FuIFlwP9orW/VWh8LQ24iiOy2Uc/ug7tN4zNS39IVb7wBV14ZvGuI05p56UziHMbEyjhHHDMvnRnhjKzjtNNVtdb/VErFAQ8BS0OfkgiF/KJ8/mvpf+H1eVm1Z5XlN+rJTMtkwccLAuKgqayEX/3K2ChHpp1GTHrXdFbfuto2d/OHk2lhUErlAxrj/tGewFal1H5Aa60zwpBfSOUX5dvmhyp3ay4en7EfrsfnIXdrrqW/5uMVx03jJqm+U/mll07unJaXBz16SEGIMLvczR9up2sxTAhLFhFQvey21+clwZlAXnaepX/A9h/bbxpbTd7uPNO4UeouXdG7N+zda/x7/vmnf74QMcq0MGit9wAopX5f55RPKVUI/FXr2Nz1pb5lt61cGM5pfY5pbDVn3GL4/nu4/PKTBUG6jKKSnVr94dTYG9wuxuhO2gC0AwYA5wF/DlFeIWfHZbdrD9RZfdltR50f7bpxvfx+Y4opQGoq9O8v006jWH5RPlmvZTFrxSyyXsuy/Bpg4dTYwpCitf6D1vpfWuv7gM5a60cxikNMql52+7GsxyzfjQSw7fttATu4bft+W4QzCq0uyV1M4wB+PyxebOypPGyYsY2mUrBggRSEKFY9bqbRNeNmIjgaWxgSq+5jQCnVrSpOBZIbeyGlVLxS6g2l1CdKqfVKqQvqnP+1UmqzUmqjUiosu5Okd03ntz/6reWLAsD8TfNNY6v57cjfmsZAYEEYP96YdpqTA+3ahSlLIaJTYwvDHcBflFJ7MaasPoCxEk1T7oDOBkq11sOrnv909QmlVC+Mge6hwGhgjlIqoQmvLU4jKT7JNLaa/h3741DGj7dTOeu/83nHjpMFYcEC+PJLuOUWaSHEiOwB2SQ4E1AoEpwJlu8eDadGLbuttd6ilJoBnA9s0Vp/04xrXQ7Mq/p8DbCw1rnLgPe11pXAYaXUDuAijPUxRRAc9Rw1ja0md2sufu0HwKd9xvTcc4cbs4w+/9yYadSvH7jdMHKkFIMYlN41nefHPc+igkWM7zfeFi3/cFFa69M/SKlfATdirG4/HHhFa/1iky6k1DLgV1rrL6riYqCb1tqvlPodcFRr/aeqcwuA+VrrFfW8zjRgGkBSUtLgCy+8sClphF1JSQmpqamRToPP//M5Fb6KmjjeGc/FZ1982udFS/5NtWNvAaq0jON+aOWAQa3imV8ZR6+yMnYlJfHzvn3xOqJ/ceFY/f7XFqqv4bj3OF+VfoVGo1D0SelDq4RWQb9OrP8f1Jf/xo0bP9Raj23wSVrr034A+YCj6vM4YE1jnlfnNf4KXFr1uQKKap27HZhVK14K9D/daw4ePFhHu2jJcdKiSZpHqPmYtGhSo54XLfk3xZtvvKE7JSi9ArQX9ArQPUC/ec45Wi9YoHVlZaRTbLRY/P7XFaqv4fHVj2vHHxyaR9COPzj046sfD8l1Yv3/oIH8Td9bG/snU6XWRrtcG909zWl35wHVax+PwehOqrYC+IlSylE1qJ0GyHZMQXRh6oUBq6temBrdLa0zMft3v2OBV5MFxANZwHxgdkKCjCFYSErLlJruQr/2k9IyJcIZWUdjt/b8SCn1DrASoytpQzOu9RqQq5TaABwHfqaUmgVs1Fp/oJRaBGwGKoC7qwuRCA5Xmos4RxwV/griHHHWvG+j6k7l7Xv3MrLOqZHA9qKiSGQlQqT0RCkOHPjx48Ahez4H0enWSnqCk2sldQH+CHwDvNPUC2mtvZy6xMbsWuefBJ5s6utGu2nTpkU6BcC4j6H21p7bvt/WqMG6aMnfVJ2lK/rGxbG2spKsWg9ZC/Tt1i1SGTZbTHz/TyNUX4MrzUViXGLNsjah+mMn1v8PmpO/6eCzUmpyQ+e01q81+WpBNmTIEL1hQ3MaLwY73U4/6M8utny/yijxGgZ2zGTzHe5IpxUcL78Mv/wl9OkDv/89C30+7puazcIKo6WwFpgYD8+88gYTJ02KdLYiiOz0OxxkyvSkWWGIdmdSGKpvp6/+a2Pl5JWW/cEqPeah11M/45Dz7zXH2vl+yte/eYOU1okRzKyZqlsISUlw9dVw/Di8+y7cfHPN+IFjvKJVHhw/Aq2S4fho8C+K3Z91IYLMtDBE/3y9ELHT7fTFh8qodHwHhzE6AoFKx3cUHyqLZFpNV3eDnBdeMI63anXKoHLSJS04dh/oR+DYfUYshGgc2xaGUPnb3/7GAw+cuo3kBx98wKBBgxg2bBgPP/xwWHPq0q4FHt8B+DuwCDgMXkrp0u7UN8uG8v/LX/5CRkYGQ4YM4dlnnw15zqdYvvxkQfB6jTuV33235rTH4+Hmm29m2LBhjBgxgr5H+gY8veuXXRk0aBCDBw9m8eLFYU7+9Ormv2zZsoDzjz/+OAMGDCA9PZ2nnnoqQlk27HT5V5s5c2a9P1+Rdrr8I/n72xiny7/Jv7+nm88azR9nMr94/d71Ou7ROM0j6LhH4/T6veub/Vpaa+3z+fTll1+uExMT9W9+85uAcydOnNDdu3fXhYWFurKyUg8ePFjv3r37jK7XVFc/OEXTEU0rNFeg73jnqYDzZvkfP35c9+jRQ3u9Xu31evV5552nS0pKQp+0z6e1x2N8vnCh1n36NHgfwiuvvKJvv/12rbXW33//vU5ITTh538bd6KRuSbqiokIfOnRIn3feedpT/bpRom7+559/fs253bt364svvlh7PB7t8/n0sGHD9M6dOyOVar3M8q+2ceNG3aFDh1N+vqKBWf7R8Pt7Omb51/f7C3TQJu+tjZ2uaklO5cSHD6c683ntDoeDpUuXkpuby86dOwPOffrppwwaNIju3bsDsGjRItq3b3/G12yK9x/7C1d/f4B1a9by8CP/w30ZdwScN8u/+ofl+PHjOBwO/H4/CQkhXMqq9iyjCRPgd7+Dm26CG29s8B6E7t27M2jQIABatGhBZXnlyfl030Jcnzji4uJo27YtF1xwAV988QWXXHJJ6L6GJqqb/7Fjx9Bao5TRFfzoo4/WfM9btGjB4cOHI5VqvU6Xf2VlJf/93//N/fffz6FDhyKZar3M8o+G39/TMcu/vt9fwGv2erbtSnIXuqn0V6LRVPorcRe6z/g14+LicNSzzMK+ffvQWnPjjTfyox/9iPnz59O6deszvl5Txcc5iHM6aJVY/98DDeXfqlUrxo0bR9++fenTpw8ZGRm0adMm+AnWHUPweuGCqkV4HQ7TG9OysrIYMGAAX3zxBVdccQXJWcknh9dOQEKbk4UsOTk56t5Y6+b/61//uuZNtUePHvz4xz+mtLSUW2+9lTZt2jBkyJAIZxzILH+Ap59+mgkTJnD22WdHMMuGmeUfLb+/Zszyr+/3V2ttuliabQuDK82F0+FEoXA6nCG94euss85iz549vPbaa7jdbjZt2sQ///nPkF2vPjkbc3h3x7uUnihl+vvTydmY0+jnrlq1ih07dlBYWMiePXs4cOAAS5cuDX6Sd9wROIbw5ZdwQ+NXYH/00Ue55ZZb+J//+R8yJ2aePNECznaefEM6ePBgVK59Uzv/+++/P+Dc+++/z4gRIxg2bBhLliwJeNONFg3l/80335CXl8cvfvGLCGZ3eg3lHw2/v43RUP71/f4qpcaZvZatu5K01mhq1mcKmUsuuYTk5GSSkpJwOBy0adOG+Pj4kF6zrmc/efaUeNrgxt344vF4SE5OJjHRmNqampoanO9ZdZfRsGHQtStMnQqjRgVMO22shQsX8tlnn/HZZ5+RmJhIy6KW/GPnP/BpH46eDvQqjd/vp7S0lMLCQqJt8cW6+ddWWFjI/fffz6pVq+jcuXOEMjRnlv+aNWsoLS0lKyuL/fv3c+LECbp06cJdd90VoWxPZZZ/NPz+no5Z/vX9/nKa6aq2LQy5W3MD7gTO3Zob9PsYZs+ezeDBgxk7dixTp05l1KhROJ1Ohg0bxpgxY4J6rdM5eOKgaVyf6vzHjBnDsmXLyMjIwOl0cskllzBunOkfHObq3KnMQw/Bo4/C0KHGRzMsXbqUb7/9tub7etRzFH9LP3QC1VvhGudi0KBBxMfH8/zzz9fbZRZJdfMHuOKKKxg8eDD79u3j2LFj3HLLLTXn5s6dy8CBAyOQaf3M8p8yZQpTpkwB4NVXX2XHjh1RVRTAPP9o+P09HbP86/v9xViotEG2vcHt+r9dz5IdS2rin1zwE965uckrfcSMTv+vE/uX7odCIBvOaXUO/77/3+FPZMkS+P3vjYLQu7fx+YQJQV/YLvPVTFbvWV0Tj+o+ilW3rgrqNYSIYdJiqM85rc4xja2mfcv27Gd/QBw2Wht7KAO8956xY9obb4SkIFTbXrLdNBZCNCy62tNhNKjTINPYarqc1cU0Donas4yqW3bPPAMFBTBpUkiXv26b1NY0FkI0zLaFoe4SvVZfsnf3od2mcVDVnXbq8cCxY8a55OSw7IfQLqmdaSyEaJhtC8Nhz2HT2GqGdxluGgeN1pCZeeq0U5crNNdrgKuHyzQW1pBflM8Ta54gvyg/0qlYim3HGLb8e4tpbDVFR4tM4zPi98OHH8LYscZYwi23wIwZzZp2GixHy4+axiL25RflMzp3dM0KyXnZeZZdITncbNtiGN9vvGlsNSEZjK3dZXTVVVB909uMGZbZQrOwsJDx48fjcrlIT09n4sSJlJYGdju63W4mTKi7B1XDlixZwvDhwxk5ciQZGRn86U9/AoypnNG4wFy0che68fq8+LQPr88blNULhMG2LYb+HfsbayVpY62k/h37RzqlkDq3zbmUUBIQN1vd+xD69DG6jKJobnebpDamcWOUl5dz1VVX8ec//xlXVVfY/PnzmT59On//+9/Nn9yAzZs388ADD7BixQo6d+6Mx+Nh3LhxNevwiMZzpblIcCaEfAc3O7JtYcjdmotP+wDwaV9IbnCLJmnJaWxhS0DcbJWV8OtfGxvlLFgQ0S6jhri/dZvGjfGvf/2Liy66qKYoANx2221cf/31bN26ldtvv52EhARSUlJqFrhLSUnh+uuvZ9OmTQwaNIjy8nJ27drFRRddxMsvv8y8efO45557au5gTkxMZNGiRfh8Pt5///3mfrm2lN41nbzsPNnBLQRsWxjspvBIoWlsqrqF8NJLxg1qSUmQlwdpaVFXEKp1btMZ9tWJm2jXrl307NkTMNZXuqFq3abi4mJSU1N58cUXufjii3nyySfZsmULAIcOHeKhhx6ia9eutGvXjo8//pi+ffvSu3dvDh06RGFhIddff33Addq1kxlTzZXeNV0KQgjYdowhe0A2ic5EFIpEZyLZA7IjnVJIfVP6jWlcr7rTTgsLYc8e49x550VtUQCYmTGTOIfxd0+cI46ZGTOb/BpdunRh7969ALRv3x63283KlSvZt28fX331FRdffDEAI0aMqHlOhw4d6N69Ow6HgxYtWtC3r7FhUMuWLfF4PHTq1Ini4uKA63z00Ue8845177oXsce2hSG9azorJ69k9mWzLb3fc7VyX7lpfIqSklN3TPvyS2M8wSbGjBnD+vXr2bhxY82xuXPnUlZWRo8ePfjyyy8BWLduXc35063BNHnyZObOnUtJiTHec+TIEWbOnEnHjh1D8BUI0TzSlWQTPr/PNAaMFsL27XDhhdChA1x8MTzwQFSOIZxO7tZcKv2VAFT6K5s1htS+fXvefvttZs6cSVlZGU6nk8zMTMaOHcusWbOYMmUK7du357zzzmv0a2ZlZfH73/+ea6+9tmZc4sEHH+TSSy/l66+/blJ+QoSKbQuDzIGupfYso9274dtvITXVWM/I5oYOHUpeXl695z799NNTju3fv7/ez6vHIABuuukmbrrpplOee+uttzY/USGCyLZdSe5CNx6fB5/24fF5LD8HukVci1Pj+nZMy8mBKNu2sDnsNoYkRDDZtjCktEzBr/0A+LWflJYpEc4otNp93ZbW60DthtbPGDE7dpw6hmCRG9PSu6Zzz4h7OK/dedwz4h77tgaFaAbbFobSE6WoqiXJFcrSi+gtXLAAx98O8J4HPMB7R8DxtwMs3LwZ3G5LFYRqORtzmLNuDt8c+oY56+Y0aStTIezOtoUhpWUKGmOTIo22dIth9u9+x2vlXrKAeCALeK3cy+xZs4wF7yxUEKo9+/GzprEQomG2LQy2WXZ782a2793LyDqHRwLbq+boW1F10W8oFtYgq6uGhm0LQ93BZksNPvv9xn0IAOecQ9/4eNbWechaoG+3buHOLGzuHXGvaSxiX/XMwodWPsTo3NFSHILItoVhw74NpnFMqjvLSGvo1IlZr7zCpHjFSqACWAlMilfMmj07svmG0LTB03jxmhe5sueVvHjNi0wbPC3SKYkgk9VVQ8e2hWFcr3GmcUypb9rp9Ok1pydOmsSQX1/NdcmQCFyXDEN+fTUTJ02KWMrhMG3wND78+YdSFCyqenVVp3LK6qpBZtsb3O4ceicLti0IiGPWX/4Cv/zlyeWv67lTOf2qDP6RZKzeeQxIvywjAomGV35Rvqy8aWGyumro2LYwPJD3wCnxqltXRSibJvL7YfFiY5XTa66BiROhVSu46aYGZxjZbStTubPdHmR11dCwbVfSroO7TOOo5PfD3/8OAwbAjTfCvHnG8VatjOJgMu3UbluZ2u3OdiGCybaFoW5/ZNT3Ty5ffrIgVFQYXUbvvhvprKKW3e5sFyKYbFsYjnuPm8ZRwe83BpIBDh06WRCacadyQUmBaWw1S79eahoLIRpm28KwZf8W0ziias8y+t//NY799KdntHRFz/Y9TWOr2ffDPtNYCNEw2xaGlvEtTeOIqG/aadUOYDgcZ7R0xZOjnzSNrWbqJVNNYyFEw2xbGK7pc41pHBF33nnqaqdV+wyfqf/77P9MY6vp37E/8Y54AOId8fTv2D/CGQkRO2xbGJbtWmYah0X1LKPqNYumTg3Z8tf/+OofprHVuAvdAYPPMitJiMazbWHYfXC3aRxSdaedvvSScXzIkJAtf90msY1pbDWuNBdxjjgUijhHXPTPOhMiiti2MFzb51rTOGSWLDHGEGpPO33kkZBftvNZnU1jK6q9rLoQovFsWxjObXOuaRxUutYb0z/+EZEd03ql9DKNrcZd6KbSV4lGU+mrlK4kIZrAtoVhccFi0zgoancZffaZceyZZyKyY1rJ8RLT2GpSWqbgp2qMAbnBTYimCFthUEolK6X+pZT6WCn1kVLqnHoe85ZSyl31EYJ36pOGdxluGp+RumMIXi8cO2aca9MmIjumje833jS2Gjtt3SpEsIVzEb37geVa62eUUtnAQ0DdJU3P1lpnhiOZC1MvNI2bTWtju8y1a01XOw236qWnFxUsYny/8ZZfitpOW7cKEWzhLAyXAxOqPv8XRqGooZRqD3RTSn2EsTXxE1rrD0KVTFBXG/X74cMPYcwY40a0SZNgxoyoKAi1TRs8zfIFodrmf282jYUQDQvnGEMKUN2ePwK0rXO+JZADjAUmAn9SSnWs+yJKqWlKqQ1KqQ0lJc3vJ3d/6zaNG6X2ncpXXQUfVNWx228P+xhCY8j+uEKIxghJYVBKPaSU+qL2B0YhSK56SHug7rv6d8CTWutKrfU+YBPQu+5ra61ztNZDtNZDUlNTm53jobJDprGp+pauWLDAaDFEqfyifLJey2LWillkvZZl+eKQPSCbBGcCCkWCM4HsAdmRTkmImBGSriSt9WPAY7WPKaVmA+OB54DrgQ/rPO0W4EpgslKqNXAhsD0U+QGUlJWYxqYqK+H++yExMWrGEE4nd2suHp8HAI/PQ+7WXEtvcJLeNR33ZLfs7iVEM4SzK+n/AVcrpT4FrgPmACilnldKDQTeqoo/xigav9dah2wqSZIzyTQOUN1CGDMGysshIQHy8iIy7bS59h/bbxoLIUS1sA0+a60PAaf0tWit764VTg5XPq0TWsPxOnFdfj+88w784Q+wbZsxy2jPHuPfnrG1bPU5rc8xja2muuusemvPlZNXSqtBiEay7Q1uRzxHTGMOHDh1DOHLL42iEIMGdRpkGltNddeZRtd0nQkhGse2heHcs849Nfb74YsvjAMpKcYNamFeuiJUZPqmEKKxwnkfQ1TxbCqj9Vtw/Ai0SoaeQ7+H/xsAu3fDt99Cx47w+uuRTjNo9h/fbxpbTfaAbP6y5S9U+CqId8bLrCQhmsCWhWHhggUcfu1r3vPCSGDtEZj60T4WnuNn4ksvGa0Fizl44qBpbDUyK0mI5lNax+6SxEOGDNEbNmxo8vMuSkvj+T17yKp1bCVwd7dufLFnT9DyiyY95vag8EhhTZyWnMa3934buYSEEJGkzE7acoxh+969jKxzbCSwvagoEumERd09CeywR4Hc6S1E89iyK6lvt26srdNiWFt13KoGdRrEniN7AmIryy/Kx/Waq2aMwT3ZLd1JQjSSLVsMs2bPJjsxjpVABUY3UnZiHLNmz45wZqEzM2Mm8Y54AOId8czMmBnhjEIrd2suXp8Xjcbr88p0VSGawJaFYeKkSbSfcgHXxUMicF0ytJ9yARMnTYp0aiGT3jWdP131J67seSV/uupP8tezEKJBtiwMAP2vGMCxfqCHwbH7jNjK8ovyufeDe8n7No97P7jX8v3u2QOySXQmolAkOhNluqoQTWDLMQaw31aX7kI3Xp8Xn/bh9XlxF1q7zz29azorJ6+U6apCNINtC8P4fuNZxrKA2MpcaS4SnAk1awe50lyRTink0rumS0EQohns25XUsX/N5w4cAbEVpXdN5+7hd9OjbQ/uHn63vGEKIRpk28IwZ92cms/9+ANiK8rZmMOcdXP45tA3zFk3h5yNOZFOSQgRpWxbGPb9sM80tppFBYtMYyGEqGbbwtA2qa1pbDV1x1CsPqYCcuezEM1l28HnT7/71DS2mmmDp7Hr0C4WFyzmhn43MG3wtEinFFL5RfmMzh1dM9iel50n4ypCNJJtWwxllWWmsdXkF+Xz/CfP8+3hb3n+k+ct/1d0fdNzhRCNY9vC0CKuhWlsNe5CN+WV5fi0D0+lx/JvlK40F3GOOBSKOEecLabnChEsti0MvVN6m8ZWc9hzuGZFVT9+DnsORzahMKj+eu2wkqwQwWTbwmC3WUlb/r3FNLYad6Ebn9+HRuPz+yzfQhIimGxbGDq06mAaW83ATgNNY6upvtPbqZy2udNbiGCxbWEYce4I09hqdh7YaRpbTXrXdOaOncvoHqOZO3auzEgSoglsO131B+8PprHV2K3rrHo1Wa/Py5q9a+jfsb8UByEaybYthk+KPzGNrWbqJVNNY6uR6apCNJ9tWwzxznjT2Gqqb2hbVLCI8f3GW/4GNzuuJitEsCitY3cq35AhQ/SGDRua9dzWj7fm+FvHjS3croJW8a049rtjwU1QRFR+Ub7sxyBE/ZTZSdu2GBzKYRpbkd3eKGU/BiGax7aFoWKNhq1VQUdwZlj7zuf8onxcr7mo8FUQ74zHPdnaO7gJIZrP+n8m16P0mAfPkDg4H7gQuASOlldQeswT6dRCJndrLl6fF43G6/OSuzU30imFnKyuKkTz2LLFUHyoDIdS+K4HnIACpY3jKa0TI52eCAJZXVWI5rNli6FLuxZAPLQEkqqPxlcdt6bsAdkkOhNRKBKdiWQPyI50SiEl01WFaD5bthhSWifSp0MPCkq/rznWp0MPS7cW0rums3LyStsMPst0VSGaz5aFAaDDWYlwsCpQVbGwjPSu6eRl59mmEAoRTLYtDHsO7zGNrcaOfe4yXVWI5rHlGAOA1+c1ja3GXejG4/MYG/X4rL9RjxCi+WxbGBKcCaax1aS0TMGv/QD4tZ+UlikRzkgIEa1sWxgGdRpkGlvN5n9vNo2FEKKabQuD3bb2FEKIxrJtYbDbVpd2ayEJIZrPtoVhfL/xprHVlJ4oxVH13+3AQemJ0ghnJISIVradrmrH/QkS4xLlhi8hxGnZdj8GO7LbsttCiAaZ7sdg264kIYQQ9bNtV5Ld2PHOZyFE84S9xaCU+o1S6vYGzv1RKbVBKfWZUmpkqHOx03r9stqoEKKxwtZiUEqdDbwNZAB31XN+NHCe1nqIUqoH8C5wcajysdtf0LLaqBCiscJWGLTW/1FKXQb8voGHXA4srnrst0oph1Kqvdb6YAOPPyPVawf5tb9m7SArFwZZbVQI0VhhHWPQWlcqpfwNnE4Bak+uPwK05eTi2AAopaYB0wC6devW7FzsuHaQrDYqhGiMkIwxKKUeUkp9UefjjtM87SCQXCtuCxyo+yCtdY7WeojWekhqamqzc5QbvoQQon4haTForR8DHmvi0/KAGcACpVRf4LDW+mjQk6siN3wJIUT9Ij5dVSn1PDBfa71cKXWtUmoz4KWquyhUpM9dCCHqJ3c+CyGE/cidz0IIIRpPCoMQQogAUhiEEEIEkMIghBAigBQGIYQQAaQwCCGECCCFQQghRAApDEIIIQJIYRBCCBFACoMQQogAUhiEEEIEsHVhsNPWnkII0VgRX101Uuy2tacQQjSWbVsM7kI3Xp8Xn/bh9XlxF7ojnZIQQkQF2xYGV5qLBGcCTuWUjXqEEKIW23YlpXdNZ+7YuSwqWMT4fuOlG0kIIarYtjDkF+Vz7wf34vV5WbN3Df079pfiIIQQ2LgrScYYhBCifrYtDDLGIIQQ9bNtV1J613TysvNwF7pxpbmkG0kIIarYtjCAURykIAghRCDbdiUJIYSonxQGIYQQAaQwCCGECCCFQQghRAApDEIIIQJIYRBCCBFACoMQQogAUhiEEEIEkMIghBAigK0Lg2ztKYQQp7LtkhiytacQQtTPti0GWXZbCCHqZ9vCIMtuCyFE/WzblSTLbgshRP1sWxhAlt0WQoj62LYrSQghRP2kMAghhAgghUEIIUQAKQxCCCECSGEQQggRQAqDEEKIAFIYhBBCBJDCIIQQIoAUBiGEEAGkMAghhAgghUEIIUQApbWOdA7NppT6QGs9NtJ5CCGElcR0YRBCCBF80pUkhBAigBQGIYQQAaQwCCGECCCFQQghRAApDEIIIQL8fx5MV0zUaHHgAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "aa = []  # intercepts\n",
    "bb = []  # slopes\n",
    "for _ in range(100): # repeat X times\n",
    "    sim_error = np.random.normal(0, sigma, transition.shape[0]) # Guassian residuals\n",
    "    simulated = transition.logK_pred + sim_error\n",
    "    ax.scatter(\n",
    "        transition.GdmCl,\n",
    "        simulated,\n",
    "        color=\"g\",\n",
    "        marker=\".\",\n",
    "    )\n",
    "    # fit for each dataset\n",
    "    sim_res = linregress(transition.GdmCl.values, simulated) \n",
    "    # track intercepts and slopes\n",
    "    bb.append(sim_res.slope)\n",
    "    aa.append(sim_res.intercept)\n",
    "\n",
    "ax.spines.bottom.set_position((\"data\", 0))\n",
    "ax.spines.right.set_visible(False)\n",
    "ax.spines.top.set_visible(False)\n",
    "# show updated figure again\n",
    "fig\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For each simulation we refitted the dataset to a line and kept\n",
    "track of the intercepts and slopes we were getting."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0, 'simulated slopes')"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA2UAAAE7CAYAAABdbSn9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAZp0lEQVR4nO3dfZQleVkf8O+zDAsBPMDCoGfVcdajskEOymZAjiIioLAOSdQEdzWEIJhxJSAYQpiIrxHNqIAYCOoadSFKEBN8CRNAgrALoiz7AhpP8OXI4DEo4m4AF4F12Sd/VA3bM8xMd9++M7++3Z/POX3m9q3quk/9bnU9/a2quVXdHQAAAMY4b3QBAAAAu5lQBgAAMJBQBgAAMJBQBgAAMJBQBgAAMJBQBgAAMNCec/Ei973vfXv//v3n4qUAGOz666//6+7eO7qOVaFHAuwOZ+qP5ySU7d+/P9ddd925eCkABquq942uYZXokQC7w5n6o8sXAQAABhLKAAAABhLKAAAABhLKAAAABhLKAAAABhLKAAAABhLKAAAABhLKAAAABhLKAAAABhLKAAAABhLKAAAABtozugDYjvYfPrrlZRw7cnAJlQDAciyjtyXL6W/6LJzImTIAAICBhDIAAICBhDIAAICBhDIAAICBhDIAAICBhDIAAICBhDIAAICBhDIAAICBhDIAAICBhDIAAICBhDIAAICBhDIAAICBhDIAAICBhDIAAICBhDIAAICBhDIAAICBNhzKquqyqjoyP760qt5VVddU1Sur6vyzVyIAbG96JABbsW4oq6rzquqNSV6+5umfTPKPuvsRSd6f5ElnqT4A2Lb0SACWYd1Q1t23J7k0ydOSpKoqyUu7+8/mWW5Jcu+zViEAbFN6JADLsGcjM3X3bVV1+/y4k/zHqrpzkmcluSzJI07+mao6lORQkuzbt29Z9QLAtqJHArBVC33QR1XdP8k7knx2ki/r7g+ePE93X9ndB7r7wN69e7dYJgCsBj0SgM3a0JmyteZLM/5bkqd399XLLwkAVpMeCcAiNh3Kklw0f/3g1HuSJFd191XLKgoAVpQeCcCmbTiUndRQ7rH8UgBgNemRAGyFm0cDAAAMJJQBAAAMJJQBAAAMJJQBAAAMJJQBAAAMJJQBAAAMJJQBAAAMJJQBAAAMJJQBAAAMJJQBAAAMJJQBAAAMJJQBAAAMJJQBAAAMJJQBAAAMJJQBAAAMtGd0AQAAsFn7Dx9dynKOHTm4lOXAVjhTBgAAMJBQBgAAMJBQBgAAMJBQBgAAMJBQBgAAMJBQBgAAMJBQBgAAMJBQBgAAMJBQBgAAMJBQBgAAMJBQBgAAMJBQBgAAMJBQBgAAMJBQBgAAMJBQBgAAMNCGQ1lVXVZVR+bHj66qG6rqnVX1w2evPADY/vRIALZi3VBWVedV1RuTvHz+vpK8LMnju/shSR5WVQ85u2UCwPajRwKwDOuGsu6+PcmlSZ42P/VFSd7f3e+fv39dkq88O+UBwPalRwKwDBu6fLG7b0ty+/ztfZLctGbyh5Pca7llAcBq0CMB2Ko9C/zMzUnuueb7C5J88OSZqupQkkNJsm/fvoWKA5L9h49ueRnHjhxcQiXABuiRfBr7cWA9i3z64h8l+ZyqurCq7pTk8UneePJM3X1ldx/o7gN79+7dap0AsAr0SAA2bdNnyrr79qp6Zqbr5G9L8qrufs/SKwOAFaNHArCIDYey7r5qzePfTPKbZ6MgAFg1eiQAW+Hm0QAAAAMJZQAAAAMJZQAAAAMJZQAAAAMJZQAAAAMJZQAAAAMJZQAAAAMJZQAAAAMJZQAAAAMJZQAAAAMJZQAAAAMJZQAAAAMJZQAAAAMJZQAAAAMJZQAAAAMJZQAAAAPtGV0ArLX/8NEtL+PYkYNLqGRnMa4AANuXM2UAAAADCWUAAAADCWUAAAADCWUAAAADCWUAAAADCWUAAAADCWUAAAADCWUAAAADCWUAAAADCWUAAAADCWUAAAADCWUAAAADCWUAAAADCWUAAAADCWUAAAADLRTKavJTVXVNVV1bVV+97MIAYNXojwAsYs+CP/eYJBd09yOq6guSvCbJg5ZXFgCsJP0RgE1b9PLFTyb5jKo6L8m9k9yyvJIAYGXpjwBs2qKh7LeTfFaS9yR5c5JfWFpFALC69EcANm3RyxcPJ3lddz+vqu6X5Nqq+pXu/tDxGarqUJJDSbJv374tFwqMtf/w0aUs59iRg0tZDmxT6/bHRI8E4ESLnim7S5K/mh9/KMnHk/TaGbr7yu4+0N0H9u7du3iFALA61u2PiR4JwIkWPVP2giQ/X1XfkKkB/Uh3f3h5ZQHAStIfAdi0hUJZd9+c5OuXWwoArDb9EYBFuHk0AADAQEIZAADAQEIZAADAQEIZAADAQEIZAADAQEIZAADAQEIZAADAQEIZAADAQEIZAADAQEIZAADAQEIZAADAQEIZAADAQEIZAADAQEIZAADAQEIZAADAQHtGFwAAAKPsP3x0y8s4duTgEiphN3OmDAAAYCChDAAAYCChDAAAYCChDAAAYCChDAAAYCChDAAAYCChDAAAYCChDAAAYCChDAAAYCChDAAAYCChDAAAYCChDAAAYCChDAAAYCChDAAAYCChDAAAYCChDAAAYKCFQ1lVPbeq3lVV11fV45dZFACsMj0SgM3Ys8gPVdVDknxTkocmuXeSq5O8dol1AcBK0iMB2KxFz5R9XZJXdPet3f2BJP90iTUBwCrTIwHYlIXOlCW5MMn9qup/JrlHkpck+d9rZ6iqQ0kOJcm+ffu2UiNn2f7DR7e8jGNHDi6hkuVYxvosw3apY7vZadvbMhiTHUePZOn0FDZiWduJnnLuLXqm7G+S3CnJwST/OMkLq+pea2fo7iu7+0B3H9i7d+/WqgSA1aFHArApi4ay30ny4e7uJB9N8rdJbl9aVQCwuvRIADZl0csXX5Pk4VX1lnkZP9jdH1laVQCwuvRIADZloVA2H/37riXXAgArT48EYLPcPBoAAGAgoQwAAGAgoQwAAGAgoQwAAGAgoQwAAGAgoQwAAGAgoQwAAGAgoQwAAGAgoQwAAGAgoQwAAGAgoQwAAGAgoQwAAGAgoQwAAGAgoQwAAGAgoQwAAGCgPaMLAABYtv2Hjy5lOceOHFzKcnaSZY0tcAdnygAAAAYSygAAAAYSygAAAAYSygAAAAYSygAAAAYSygAAAAYSygAAAAYSygAAAAYSygAAAAYSygAAAAYSygAAAAYSygAAAAYSygAAAAYSygAAAAYSygAAAAYSygAAAAbaUiirqvOq6neq6nHLKggAVp3+CMBmbPVM2TOSXLyMQgBgB9EfAdiwhUNZVe1LcmmS31heOQCw2vRHADZrzxZ+9iVJnp3kOaeaWFWHkhxKkn379m3hZXam/YePLmU5x44cXMpyYJUs4/fH7w5n0Rn7Y6JHwk6zrL/r2L0WOlNWVU9M8vvd/Qenm6e7r+zuA919YO/evQsXCACrYiP9MdEjATjRomfKvjLJA6vqLZmumb+kqm7p7rctrTIAWD36IwCbtlAo6+5vP/64qq5K8ioNB4DdTn8EYBFb+T9lSZLufvIS6gCAHUV/BGCj3DwaAABgIKEMAABgIKEMAABgIKEMAABgIKEMAABgIKEMAABgIKEMAABgIKEMAABgIKEMAABgIKEMAABgIKEMAABgIKEMAABgIKEMAABgIKEMAABgIKEMAABgoD2jC1hF+w8fHV0CAGxLeiTA5jlTBgAAMJBQBgAAMJBQBgAAMJBQBgAAMJBQBgAAMJBQBgAAMJBQBgAAMJBQBgAAMJBQBgAAMJBQBgAAMJBQBgAAMJBQBgAAMJBQBgAAMJBQBgAAMJBQBgAAMJBQBgAAMNBCoayq7lJVv1xV11bV71bV1y67MABYRXokAJu1Z8Gf++YkN3f3ZVW1N8nbk3zh8soCgJWlRwKwKYuGsvcluXF+/LEk96iq6u5eTlkAsLL0SAA2ZaFQ1t1vTpKqemCSn03ywpObTVUdSnIoSfbt27fFMpP9h49ueRlJcuzIwaUshxMt6/2Bc8U2y9kyokdy9thXsBttl+1+WX+3L2N9znaGWPiDPqrq+5K8Msn3dPcLTp7e3Vd294HuPrB3796t1AgAK0WPBGAzFjpTVlXfnOQhSR7S3Z9YbkkAsLr0SAA2a9H/U3ZpkouSvKGqkiTd/cgl1QQAq0yPBGBTFv0/ZU9adiEAsBPokQBslptHAwAADCSUAQAADCSUAQAADCSUAQAADCSUAQAADCSUAQAADCSUAQAADCSUAQAADCSUAQAADCSUAQAADCSUAQAADCSUAQAADCSUAQAADCSUAQAADCSUAQAADLRndAFszf7DR0eXAGwDy9gXHDtycAmVMIJeALDanCkDAAAYSCgDAAAYSCgDAAAYSCgDAAAYSCgDAAAYSCgDAAAYSCgDAAAYSCgDAAAYSCgDAAAYSCgDAAAYSCgDAAAYSCgDAAAYSCgDAAAYSCgDAAAYSCgDAAAYSCgDAAAYaKFQVlV3rqpfrKp3VNXbq+riZRcGAKtGfwRgEYueKXtSkpu6+8uSHE7ywuWVBAArS38EYNMWDWWPSfKa+fFbk3zpUqoBgNWmPwKwaXsW/Ln7JLkpSbq7q6qr6rzuvv34DFV1KMmh+dtbquoPt1bqctSPfurhfZP89bhKthVjcSLjcQdjcaIdPR5r9o8bcaax+LwtF7O61u2PyfbtkVu0034/dtr6JDtvnXba+iQ7b522tD6b7Etn1VzLVt+f0/bHRUPZzUnumSRVVZl6zwkNp7uvTHLlgss/66rquu4+MLqO7cBYnMh43MFYnMh43MFYnNa6/THZ/j1yETttm9hp65PsvHXaaeuT7Lx1sj4bt+jli29K8k/mx4/NdIkGAOx2+iMAm7bombKXJ3lFVV2X5KNJnri8kgBgZemPAGzaQqGsu29NcvmSaznXdtRlI1tkLE5kPO5gLE5kPO5gLE5hh/THRe20bWKnrU+y89Zpp61PsvPWyfpsUHX32Vo2AAAA61j0/5QBAACwBLsilFXVhVX1v6rq6vnrolPM8+qqesv89ZpTLWen2OB4vKiqrquqd1bVw0fUeS5U1QVV9Yaq+p2qeltVPfgU8+ymbWMj47Erto3jquq5VXXFaabtmm3juHXGY1dtG7tRVd25qn6xqt5RVW+vqotPmv64qrqxqq6tqh8cVedmbGCdnjI/f11VPXNUnYuoqsuq6sgpnn/2/D5dX1XfOKK2RZxhfb67qt49967njqhtUadbpzXTf+xM07ebM7xHK7dvSM64PsvfL3T3jv9K8gtJvmF+/JQkLz3FPFePrnO7jEeSRyf59fnxRUl+b3TNZ3EsfiDJc+fHD0zytl2+bZxxPHbZtvGZSa5JcluSK04zz27aNs44Hrtp29jNX0memuQn58ePSHJ0zbS/l+RYpvvw3CnJdUkuGl3zFtfpbkn+NMmd568/SXLf0TVvYJ3OS/LGJB9PcuSkaV+Y5J2ZPlfgXvM6nT+65i2sz0VJ3p3k/Hm+dyT5wtE1b2Wd1sxzSZIPnm76dvpa5z1auX3DOutzVvYLu+JMWZJXZRrYJLkgyd+snVhVFyTZt+bs0ePOdYHn2BnHI8ljkrwmSbr7vUnOm8doJ3p3pvFIkluS3HvtxF24bZxxPLKLto3u/kCSRyX5kVNN323bxnrjkV20bexyn3qfM33c/5eumfbQJDd29/u6+5OZbg2wCjfBPdM61fx190x/WJ6X5NZzWdwiero33qVJnnaKyY9K8truvq27P5TkPZkOwm1b66xPknxfd986z/exTGFzW1tvnapqT5IfT/KCc1nXotZZn5XbN6yzPmdlv7DoR+KvlO5+Q1WdX1U3Jrl/ph3wWnfL9GkqP57kfkmuqaov7+6/OselnhMbGI/7JLlpzfcfzrSDu/ncVHjudPevJklVfUWSn0ny70+aZbdtG+uNx67ZNpKku2+rqk+78e9sV20bybrjsau2jV3sU+9zd3dVdVWdN/8Bc2Gme2b/SpLPSvLmJN8/rtQNO+06dfdHq+p1Sf7PPO+buvsjwyrdhDP8vp7ud3VbO936zAeB3ltV90nywiQfyXQmZttbZ5/67EwHSf8un36AdFs6w/qs5L7hDNvcWdkv7LhQVlXfm+Syk57+qe7+T0keXFUPy/TH5pesmf5/M52a7CTvr6obknxRkpX/42rB8bg5yT3XfH+vbPMjGhtxmrG4MlMwfXCSJ3X3DSdN323bxnrjsZu2jZd198vO8GO7bdtYbzx25LbBp/nU+1xVlSnHHP+j5W8yXZ70FUk+keTXkxxM8toBdW7Gadepqr4qycVJ9ifpJL9RVZd29+sG1boMN2cKZsddkOkSuZVVVY9P8hPz10/P++WVVVVfkOmS8Mcm+ReDy1mGVd03nNLZ2i/suMsXu/uHuvuBa7+SfFtVffE8yy2Zbui51rckuSpJquoeSb44d6TflbbgeLwp06nlVNXfT/KhVTkyeCanGYu7Z7q++eGnCCDJ7ts21huPXbNtrBNAkl22bWxgPHbktsGn+dT7nOkPxreumXZDprMuH58vUfpIpqP8292Z1ukuST7c3Z/o6R50H8x02dIq+60kX19V51XV3kx/WP7B2JIWV1X7M13i91Xd/bI1BwlW2VdmCs5vTnI4ybdU1dPHlrQlq7pvOJ2zsl/YcWfKTuMZSX6+qj42f/+0JKmqlyT5uSSvTvK1VfW7ST6Z6drkm065pJ3hjOPR3W+sqn84X954a5JDg+o8Fy7NdBnab00HSPOX3X35Lt42zjgeu2zbOKVdvG2ckm1j13l5kldU1XWZDug9saqel+T67n59Vf1cpkt5P5nk2u5+w8hiN+i065TkDZl+z9+e6ff8hiQreZbspPfpvye5MdMfxs9YxSCz5j26MMk9krxy7ltJ8qzufteg0ha25j36hUwfypaqenKSi7v7pSNrW8QO2Dec4GzvF9w8GgAAYKAdd/kiAADAKhHKAAAABhLKAAAABhLKAAAABhLKAAAABhLKWGlV9biqWvijt6tq//yR5qebfq+q+rZNLO+qqnrcSc8drqqHnuFnHllVBzb6GltVVU+uqvueq9cDYPtZhf65lflg1QhlrLTufn13X3kWX+JeSTbcVE6lu49097VnmOWRSTYUyqpqGb+zT04ilAHsYqvQP2E32S03j2YHqKoHJPmZTDfquzVTuPjaJBcn+ekkr0/yriQPzHRjv4uSfF6S/9HdP1BVb0lyRXe/p6quSPJZSa5as/zLkzxzXv5HkjwhyX9I8oCqOpzkpfPrX5jkzkme3t3vqqrnJLks093q736Kuq9K8qr59Z4wP/25SX42yW/O63HrfMTxrvNrnp/kjzLd2PuyJJdnuoP8d1TV1yR5YqaDKr/U3S+uqqfOy7lrkl/u7hfMr7snyb5MN9Z8RpLPSfKlmW6U+tVJfinJfZLcKclzu/ut670PAKyWFe6fT0jynCQfS/KnSZ66ZtpnJLlyzTK/f76J/bEkr03yD5L8bZJ/1t1/WVXfm+QxSe6W5MXd/UtV9YxM/fS2JFd393dvfnRhOYQyVskTkrxpbhCPyhQ21vrMJN+eabv+8ySfnWlH/+dJfmADy//8JF/T3bdU1euTXJLk3yW5qLuPVNUPJbmmu3+mqu6f5L9U1ZOSfGOShyWpTHd6P5O7ZmqE90vy2939kjk8/eXcoH4/yaO7+6+q6keSPCXJRzMFsEcluX+Sf57ky+fnXlFVX5TkSZnOuHWSq6vq+J3l/7y7n1hVD0ryyu5+YFV9R5IrMp2du1uSR89j9bANjBEAq2dV++fTMh0wfPN88HHvmmnPTfLO7n5RVd0vyY1V9bnztNd299Or6juTfE9V/WqSB3T3V1XVXed5X5fkuzL13PdW1bOr6i7d/YkNrC8snVDGKnlxkn9TVa9O8vFMR+HW+uPu/vC8w/1Ad9+cJFVVp1jWqS4D/FCSl1TV32Y6k3Wnk6ZfkuRxVfXN8/efkeRBSd7R3bfNr3XDOutwXXd/MslfVNXd1k6oqr1J9id59Vzy3TKdSfuTJDd0d8/h6m3HXy/J5fORxM9P8qb5uXtmOsqZJG9Nku7+vaq6z0m1XJPkAUlekTngrVM7AKvpxVnN/nlFku+cDyZ+IMmr10x7cJLnJcl8IPOmTFd+JHPvS/K2TAdCL0lyYD7jl0wHMPdlugrleVV1z0xXp9x+ihrgnPB/ylgl35rk17r7mzJdXvHMk6avtzP9RO64POJL1k6oqguSPLu7vzXJszNdypBMO+7j/jDJj3b3IzNd7nDV/NxDq+pOVXV+koesU0Of4bmbkvxZkq+bX+PHkrzxpHn/OMmDa3LPqvqtTI3kXfPPfHWS/5o7jjheMq/fg5K876TXe1SSj3X35ZmOFr5gndoBWE2r2j//dZLnzHXfNcnBNdP+JPMVHlV1YZK7dPcH52mXzP8+PMkfzK/1hvn1vybJr2Xqp1ckOdTdT0jy0EyXPMIQzpSxSq5J8rKq+rtM149/R05qDut4WZIfrar3Zbo04y/WTPt/SW6oquszhZc3Zfo/WJcnuaCqnpXk+Un+c1X9q0zN5nu6+91VdTTJtZkC1XsXWK/rkrwo05G9ZyV5/Xx08r3zOh4/65XuvnEOYm/PdCTyJ+Yarq6qt87PHe3uv5gPcF4yP3/3JP9yXszvZmqIj03yb6vqKZmOWj5/gdoB2P5WtX/+dpJrqurW+XVel+T4Jy8+P8lV89m383Lih4o8taqOZAqT3zCfBXzkfKbs/CQ/190frao/TPKOqro9U8hb778gwFlT3ac6cA+suuMfMNLdrx9dCwCcC/MHfVzc3R8fXQtshssXAQAABnKmDAAAYCBnygAAAAYSygAAAAYSygAAAAYSygAAAAYSygAAAAYSygAAAAb6/xZYjxKXOkC3AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1080x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "histograms = plt.figure(figsize=(15, 5))\n",
    "ax0, ax1 = histograms.subplots(1, 2)\n",
    "ax0.hist(aa, bins=20)\n",
    "ax0.set_xlabel(\"simulated intercepts\")\n",
    "ax1.hist(bb, bins=20)\n",
    "ax1.set_xlabel(\"simulated slopes\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To express the uncertainty around our parameter, we can report them\n",
    "together with the standard deviation of their simulated distributions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "jupyter": {
     "outputs_hidden": false,
     "source_hidden": false
    },
    "nteract": {
     "transient": {
      "deleting": false
     }
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fitted equation with erros:\n",
      "log K = -2.3 (± 0.6) + [GdmCl] • 1.1 (± 0.2)\n"
     ]
    }
   ],
   "source": [
    "std_intercept = np.array(aa).std()\n",
    "std_slope = np.array(bb).std()\n",
    "print(\n",
    "    \"Fitted equation with erros:\\n\"\n",
    "    f\"log K = {INTERCEPT.round(1)} (± {std_intercept.round(1)}) + \"\n",
    "    f\"[GdmCl] • {SLOPE.round(1)} (± {std_slope.round(1)})\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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.9.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}