Sample Code- Hourly Analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# An Application of Gaussian Mixture Regression for the Correction of Low Cost $PM_{2.5}$ Monitoring Data in Accra, Ghana\n",
    "## Hourly Analysis\n",
    "*Celeste McFarlane, Garima Raheja, Carl Malings, Emmanuel K. Appoh, Alison Felix Hughes, Daniel Westervelt*\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import required packages\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "from matplotlib import colors as mpc\n",
    "from sklearn.mixture import GaussianMixture\n",
    "import pandas as pd\n",
    "import math\n",
    "import seaborn as sns; sns.set()\n",
    "from sklearn.metrics import r2_score \n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.metrics import mean_absolute_error\n",
    "from sklearn.model_selection import KFold\n",
    "from sklearn import linear_model\n",
    "from sklearn.ensemble import RandomForestRegressor\n",
    "from gmr import GMM\n",
    "\n",
    "from IPython.display import clear_output\n",
    "\n",
    "%load_ext autoreload\n",
    "%autoreload \n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import data\n",
    "# INSERT PATH TO FILE HERE\n",
    "path_h = \"C:/Users/cmcfa/OneDrive/Desktop/2020 Research/Accra/Full Time Series Scripts/hourlyPA_embassy.csv\"\n",
    "data_h = pd.read_csv(path_h, header = 0, index_col = 0)\n",
    "\n",
    "# subset just to get columns for model\n",
    "model_data_h = data_h[[\"purpleair_pm2.5\", \"avg_temp\", \"RH\", \"embassy_pm2.5\"]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define functions for model evaluation\n",
    "def n_bias(estimated, true):\n",
    "    estimated_ = estimated.reshape(estimated.shape[0])\n",
    "    \n",
    "    true_ = true.reshape(true.shape[0])\n",
    "    return np.sum(np.subtract(estimated_, true_))/estimated_.shape[0]\n",
    "\n",
    "def cvMAE(estimated, true):\n",
    "    estimated_ = estimated.reshape(estimated.shape[0])\n",
    "    true_ = true.reshape(true.shape[0])\n",
    "    return np.sum(abs(np.subtract(estimated_, true_) - n_bias(estimated_, true)))/np.sum(true_)\n",
    "\n",
    "def mnbias(estimated, true):\n",
    "    estimated_ = estimated.reshape(estimated.shape[0])\n",
    "    true_ = true.reshape(true.shape[0])\n",
    "    return np.sum(np.subtract(estimated_, true_))/np.sum(true_)\n",
    "\n",
    "def nmae(estimated, true):\n",
    "    estimated_ = estimated.reshape(estimated.shape[0])\n",
    "    true_ = true.reshape(true.shape[0])\n",
    "    return np.sum(np.abs(np.subtract(estimated_, true_)))/np.sum(true_)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Build Model With 80:20 Training:Testing Data Split"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Model Selection via Bayesian Information Criterion"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZkAAAEJCAYAAABGw1qNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABPVElEQVR4nO3deVxUVR/H8c+wgwPiAgKiuCMuITVaLklqqYi7aSoulWXulrmMe/pIro/yuJWWlmupleCCu6Fp5pJb7omCICIgKoKsw33+QCdJcHCBAfy9X695JWfu8p2p+HnuufcclaIoCkIIIUQ+MDF2ACGEEMWXFBkhhBD5RoqMEEKIfCNFRgghRL6RIiOEECLfSJF5REZGBpGRkWRkZBg7ihBCFAtSZB4RHR1NixYtiI6ONnYUIYQoFqTICCGEyDdSZIQQQuQbKTJCCCHyjRQZIYQQ+cbM2AGEEMVLeno6kZGRpKSkGDuKeMGsrKxwdXXF3Nw8z/tIkRFCvFCRkZHY2tpSqVIlVCqVseOIF0RRFG7dukVkZCSVK1fO835yuUwI8UKlpKRQpkwZKTDFjEqlokyZMk/dQ5Ui84IoioKsmiBEFikwxdOz/HuVIvOCnP10ARHfbTN2DCGEKFSkyLwgik5H/MG/jB1DCPGIKVOm0KFDB9q0aUOdOnXo0KEDHTp04Oeff87zMTp06PDE9/fs2cP//ve/543KggULaNy4MR06dKB9+/b4+PjwzTffGNzv9OnTzJ49+7nP/6h169axZcuWF3IsGfh/QdQebkStDyEzPQMTc/lahSgMJk+eDGTdjNCnTx+CgoKe+hiG9mnRogUtWrR4pnz/1r17d4YOHQpAfHw8ffv2xd7enq5du+a6z+XLl7l169YLOf9DJ06coEGDBi/kWPLb8AWx9aiEkp5B0uXr2Hq4GTuOEMKA5s2b88orr3D+/HnWrl3LypUrOXToEHfv3qVUqVIsWLAABwcH3N3duXjxIgsWLODmzZuEh4dz/fp1unbtysCBA/nll184cuQIM2bMoHnz5rRv354DBw6QnJzMzJkzqVOnDpcuXUKr1aLT6dBoNOzfv59du3Y9MV/p0qUZOHAg33zzDV27duXSpUv85z//4f79+8THx/PBBx/QsWNH5s+fz/379/nqq6/o3bs348aN4+bNm8TExKDRaJg1axY3b95k5MiR3L9/HxMTEyZMmEC9evU4ffo006dPJyUlhVKlSjFlyhQiIiLYu3cvf/zxBw4ODrz55pvP9T1LkXlB1A8KS+L5cCkyQjxwc9sf3NxyKF+OXa5tQ8r5vPFcx2jatCkBAQGEh4dz5coVfvzxR0xMTBg9ejSbN2/mww8/zLb9xYsXWbNmDffu3ePtt9/Gz8/vsWPa29vz008/sWrVKpYsWcKCBQvQarUMHz4cb29vvv/+e3Q6XZ7y1ahRgytXrgCwYcMGBg0aRMOGDYmIiKB9+/b06dOHYcOGceTIEQYOHMiWLVvw8PBg/vz5pKWl4evry9mzZwkJCeGtt97io48+4vDhw/z555/UqlWLCRMm8PXXX+Pi4sJvv/3GxIkT+f7772nevDkNGjR47gIDUmReGEun0pjZq7l3PhxnY4cRQuSJp6cnAG5ubowZM4YNGzZw9epVTp48ScWKFR/b/vXXX8fCwoIyZcpgb2/PvXv3Htvm4S/m6tWrs3PnTu7cucP169fx9vYGoEuXLqxcuTJP+VQqFVZWVgBotVp+++03lixZwsWLF7l///5j27dt25bTp0/z/fffc+XKFe7cucP9+/dp2LAhQ4cO5fz583h7e9OrVy/CwsKIiIhg4MCB+v0TExPzlOtpSJF5QVQqFbYebiSeDzd2FCEKjXI+bzx3byM/WVpaAnDmzBk+//xz3n//fVq1aoWJiUmOjyQ83B6y/p9/0jYPb/c1NTV95scbLl68SNWqVQH49NNPsbOzo1mzZrRp04atW7c+tv2qVavYsWMH3bp1o1GjRly6dAlFUXjttdfYunUrISEhBAcHs3HjRsaMGYOrq6t+zEmn0xEXF/dMOZ8kX+8uS0xMpG3btkRGRrJv3z79nR0dOnTgjTfe4JNPPsm2fUhICM2bN9f/nJCQQP/+/fHx8cHPz4/Y2FgA0tLSGDVqFD4+PnTq1InQ0FAg61mVmTNn0rp1a9q0acOff/6Znx/vMWoPN+6H3UCXnFqg5xVCPJ+jR4/SoEEDevToQbVq1Th48GCeL2kZYmtrS8WKFdm3bx8AmzdvztN+MTExfP311/pLcgcPHmTYsGG8/fbbHD16FMgqDKampvqFFg8ePMh7771H+/btUalUXLhwgczMTGbNmkVQUBCdOnVi0qRJnDt3jipVqnD37l2OHTsGwM8//8zIkSOBrML4oj5/vvVkTp06xYQJEwgLCwPA29tb312MjY2lR48ejB07Vr99XFwcM2fOzHaMgIAANBoNS5cuJTAwEH9/fwICAli1ahXW1tZs27aNo0ePotVq2bBhAzt27CA0NJTg4GDCw8Pp378/27Ztw8ysYDpsth5ukKmQeDGCkvWqFcg5hRDPr02bNgwZMoR27dphbm6Ou7s7kZGRL+z4M2fOZNy4cQQEBODu7q6/BPZvP/74I7t379b3kt577z18fX0BGDp0KD179sTOzo7KlStTvnx5IiMjeeWVV1i4cCFz5syhb9++fPHFFyxfvpwSJUrg5eVFZGQkvXv35vPPP2fjxo2YmpoyefJkLCws+N///oe/vz+pqamo1Wr97+BGjRoxd+5cbG1tad269fN9eCWfjBs3Tjl69KjSrFkzJSIiItt7I0eOVL755ptsbZ988omyZcsWpVmzZvq2Zs2aKVFRUYqiKEp6erri5eWlpKWlKb169VKOHj2q365FixbK9evXFa1Wq2zcuFHf3qdPH+XIkSM55rt7964SERGR7XX06FGlRo0aj+XNq9Rbd5XfGg5UItfufqb9hSgOzp07Z+wIhc6CBQuUmzdvKoqiKDt27FCGDBli5ETP7mn//ebbX/H9/f1zbA8LC+PIkSPZ3l+5ciW1atXSD8I9FBMTg4ODAwBmZmao1Wri4+OztQM4ODgQHR1NTEwMjo6Oj7XnZMWKFSxcuPCZP19OLErbYVmuFPdkXEYI8QgXFxc+/PBDzMzMsLOzy/X3Y3FU4AP/69ato2fPnlhYWABw6dIldu7cyffff59rQXiUiUnOw0i5DdTltn3fvn3p1KlTtrbo6Ogcb0l8GmoZ/BdC/Evnzp3p3LmzsWMYRYEXmT179rBs2TL9z9u3byc2NpYuXbqQnp5OTEwMPXv2ZO3atTg6OhIXF4eTkxMZGRkkJiZib2+Po6MjsbGxuLllPY8SGxuLo6Mj5cqV098c8Gh7Tuzs7LCzs3vhn09d041bISdJT0jC3K7ECz++EEIUJQU6d1l8fDwpKSlUqFBB3zZs2DB27NhBUFAQS5cuxdHRkbVr1wJZNwsEBgYCEBwcjEajwdzcHG9vb/1td8eOHcPS0hIXFxeaNm3K5s2b0el0hIeHExYWRt26dQvyI2Jb65+HMoV4WeV0VUEUfc/y79VgTyYtLY1du3YRHx+f7QR9+vR56pNFRkbi5OSU5+2HDx+OVqvF19cXW1tb5syZA0Dv3r2ZNGkSvr6+WFhYMGvWLABat27N6dOnad++PZA1LpTbXRz5Re2e9QBX4vlwSr1eq0DPLURhYGVlxa1bt2RNmWJGebBo2dP+TlUpBkrTsGHD2LVrV7YCo1KpOH/+/LMlLcQiIyNp0aIFe/bswdXV9ZmP82f3L7B2c6LWzAEvMJ0QRYMsv1x85cvyywcOHODtt9+me/fuT3Xgl5m6pht3j18ydgwhjMLc3PyplucVxZvBIuPq6krLli1p3LhxQeQpFtQebsTuPEpq7B0sHeyNHUcIIYzGYJFp1qwZs2bNIjIykhIlsu6WUqlU9O7dO9/DFVWPDv5LkRFCvMwMFpklS5YAZFv5TYrMk5WoXgFMTUg8H06Zpp6GdxBCiGLKYJHJ65TU4h+mVhaUqOIiT/4LIV56BotMgwYNSEhI4ODBg6hUKho3boytrW1BZCvS1DXduLXvJIqiyG2cQoiXlsEic/nyZT744APi4uJQFAUnJyeWL19OlSpVCiJfkWXr4cbNzQdJuR6HtauD4R2EEKIYMvjE//Tp07GxsWHu3LnMnTsXKysrpk+fXhDZirRHl2MWQoiXlcGezMmTJ/H399evKaBSqZg4cWK+ByvqbKq6YGJhzr3z4Ti8ozF2HCGEMAqDPRlra2suX76s/zk0NLTAp2opikzMTClR3VV6MkKIl5rBnoyvry+LFi3i4MGDQFbPplevXvkerDhQe7hxc+shFF0mKtMCnYtUCCEKBYNF5vPPPycjI4NNmzZhYmLCu+++y+eff14Q2Yo8Ww83bvwUwv2waEpUdTF2HCGEKHAGi4yFhQUTJ06UcZhnoPb4Z0ZmKTJCiJdRrkVmyJAhDBs2jPnz5z/2nkqlYsGCBfkarDiwrlgOUxsrEi+EU65tQ2PHEUKIApdrkdm9ezd+fn7s3r37sffk4cK8UZmYoK5ZQZ78F0K8tHItMhcuXACylkeuVKlSQeUpdtQebkStDyEzPQMT8wJf7VoIIYzK4C1PXbp0ITg4uCCyFEu2HpVQ0jNIunzd2FGEEKLAGSwyderU4cyZM6SnpxdEnmJHnvwXQrzMDF6/uXHjBkeOHGHFihVYW1tjYmKCSqXi8OHDeTpBYmIi3bt35+uvvyY0NJS5c+fq37t58yaenp4sWbKE3bt3s2DBAhRFwdXVlenTp1OyZEmioqIYNWoUt27donLlysyZM4cSJUqQkJDAyJEjiYiIoHTp0gQEBODg4EBaWhrjx4/nzJkzWFlZMWfOHKpWrfrs39BzsnQqjZm9mnvnw3E2WgohhDAOgz2ZcuXKodFoePXVV/Hw8MDd3Z0aNWrk6eCnTp2iR48ehIWFAeDt7U1QUBBBQUF8++23qNVqxo4dS2JiIl988QVLly5l06ZNuLu76+9emzJlCj179mT79u3UqVOHxYsXAxAQEIBGo2Hbtm107doVf39/AFatWoW1tTXbtm1j3LhxaLXaZ/leXhiVSoWth5v0ZIQQLyflKWRmZj7N5sq4ceOUo0ePKs2aNVMiIiKyvTdy5Ejlm2++URRFUeLj45Vdu3bp39u2bZsyaNAgJS0tTfHy8lLS09MVRVGUqKgopXnz5oqiKEqzZs2UqKgoRVEUJT09XfHy8lLS0tKUXr16KUePHtUfq0WLFsr169fzlDciIkKpUaPGY1mfV9g3m5XfGg9SMu6nvNDjCiFEYWewJxMfH8/w4cN59dVX+euvv/jkk0+4efNmngqYv78/Gs3jk0OGhYVx5MgR+vTpA0CpUqV4++23AUhJSWHp0qW8/fbb3L59G7VajZlZ1lU9BwcH/bljYmJwcMiaQt/MzAy1Wk18fHy29of7REdHP5YhISGByMjIbK+ctnsRbD3cIFMh8WJEvhxfCCEKK4NjMtOmTePcuXMkJyeTkJDAyZMnGTduHMuWLXvmk65bt46ePXtiYWGRrf3evXsMGjSImjVr0qlTpxyL2ZOe0TExyblm5tS+YsUKFi5c+JTJn82jg/8l61UrkHMKIURhYLAnc/DgQQYNGgSAnZ0d48eP5/jx48910j179tCmTZtsbTExMfTs2ZOaNWvqx1dKly5NYmIiOp0OgNjYWBwdHQFwdHQkLi4OgIyMDBITE7G3t8fR0ZHY2Fj9cR/d51F9+/Zlz5492V5r1qx5rs+VG4vSdliWKyUPZQohXjoGi4yVlZX+l7ZOp+PkyZOULVv2mU8YHx9PSkoKFSpU0LfpdDoGDBiAj48P48eP1/dWzM3N0Wg0+ud0AgMDadq0KZB1E0FgYCAAwcHBaDQazM3N9TcXABw7dgxLS0tcXB6fN8zOzg5XV9dsLycnp2f+XIaoZfBfCPESMni57L333tPfdty7d290Oh3Dhw9/5hNGRkY+9st87969nDt3Dp1Ox44dO4Cs53P8/f2ZPHkyWq2Wr776CmdnZ32W4cOHo9Vq8fX1xdbWljlz5ugzTpo0CV9fXywsLJg1a9YzZ32R1DXduBVykvSEJMztShg7jhBCFAiVoiiKoY0CAwMJCQkhIyODN998k/fee68gshW4yMhIWrRowZ49e3B1dX2hx75z7AJnhs2n9rwhlHq91gs9thBCFFYGL5eNHTuWunXrEhAQwMKFC3F3d2fSpEkFka1YUbv/M+2/EEK8LHK9XLZgwQISExMJDAwkKSkJZ+es59UvX77MiRMnmDp1aoGFLA7MbG2wrugog/9CiJdKrkXGxMSEFStWoFKp2LlzZ7b3WrRoke/BiiN1TTfuHr9k7BhCCFFgci0y/fv3p0OHDvTp0wetVkvt2rWBrDu+crolWBim9nAjdudRUmPvYOlgb+w4QgiR73ItMqGhobi5uennCrt3757+vfj4eGrWrJn/6YoZ21r/PJQpRUYI8TLItch06tSJ5cuX88EHH+T4lP358+fzNVhxVKJ6BTA1IfF8OGWaeho7jhBC5Ltci8zgwYNxdXVlyJAhBZmnWDO1sqBEFRcZ/BdCvDRyLTIPi0t0dDS9evWSy2MviLqmG7f2nURRlCfOwyaEEMWBwedkdu3axbVr1woiy0vB1sONjIQkUq7HGTuKEELkO4PTyjRs2JDZs2dz7tw57O3t9TMaP5ymXzydR2dktnZ1MLC1EEIUbQaLzPbt2wH4+uuv9W0qlUqKzDOyqeqCiYU5986H4/DO42vtCCFEcWKwyKxcubIgcrw0TMxMKVHdVaaXEUK8FJ5YZG7evEmVKlX0U/tfvHgRV1dXSpSQWYSfh9rDjZtbD6HoMlGZGhwWE0KIIivX33AXLlygbdu2+qn3AWbMmEGbNm2IiJBlhJ+HrYcbmcmp3A/Pn+WehRCisMi1yMydOxdLS0vc3d31bT4+Pty/f5/Zs2cXSLji6tHBfyGEKM5yLTKnTp3i008/RaP5Z3C6W7dufPrppxw5cqRAwhVX1hUdMbWxkiIjhCj2ci0yaWlpqNXqx9rLlClDSkpKvoYq7lQmJqhrVpAn/4UQxV6uRaZu3bp8//333L17V9+WkJDAypUrqVVLVnZ8XmoPN5IuXyczPcPYUYQQIt/kenfZsGHDeP/993nrrbeoUqUKiqJw9epV0tLS+PbbbwsyY7Fk61EJJT2DpMvXsX0wRiOEEMVNrj0ZjUbD8uXLqVOnDleuXCE0NJSaNWuydOlSGjZsmOcTJCYm0rZtWyIjI9m3bx8dOnTQv9544w0++eQTIGtW5y5dutCqVSvGjx9PRkbW3/CjoqLw8/OjdevWDBw4kKSkJCCrV9W/f398fHzw8/MjNjYWyLrMN2rUKHx8fOjUqROhoaHP/OXkJxn8F0K8FJR8dPLkSaVt27ZK7dq1lYiIiGzvxcTEKC1atFCuXr2qKIqi+Pr6KidOnFAURVHGjh2rrFmzRlEURenfv7+yZcsWRVEUZeHChcqsWbMURVGUKVOmKEuWLFEURVE2btyoDB8+XFEURfn222+ViRMnKoqiKEeOHFHefffdPOeNiIhQatSo8VjW/JCZmakc8hmlXJy2Mt/PJYQQxpKvTwKuX7+eyZMn57iS5qxZs+jevTuVKlXi+vXrpKSkUK9ePQA6d+7M9u3bSU9P5+jRo7Rq1SpbO0BISAjt2rUDoG3btuzfv5/09HRCQkJo3749APXr1+f27dtERUXl58d8JiqVClsPN+nJCCGKNYPTyjwPf3//HNvDwsI4cuSI/v2YmBgcHP6ZLNLBwYGbN29y+/Zt1Go1ZmZm2dr/vY+ZmRlqtZr4+PgcjxUdHY2Li0u2DAkJCSQkJGRri44u2Icj1R5u3D58Dl1yKqbWlgV6biGEKAj5WmRys27dOnr27ImFhQUAiqI8to1Kpcq1PTcPZ4jOS/uKFStYuHBhXiPnC1sPN8hUSLwYQcl61YyaRQgh8kOeikxMTAy3b9/O9kv/eRYx27NnD8uWLdP/XK5cOeLi/llfJTY2FkdHR0qXLk1iYiI6nQ5TU1N9O4CjoyNxcXE4OTmRkZFBYmIi9vb2ODo6Ehsbi5ubW7Zj/Vvfvn3p1KlTtrbo6Gj8/Pye+XM9rUcH/6XICCGKI4NFZuHChSxatOix9vPnzz/TCePj40lJSaFChQr6tvLly2Npacmff/7Ja6+9RmBgIE2bNsXc3ByNRkNwcDDt2rXTtwN4e3sTGBjIgAEDCA4ORqPRYG5ujre3N0FBQWg0Go4dO4alpeVjl8oA7OzssLOze6bP8KJYlLbDslwpeShTCFFsGSwyq1evpnLlyrzzzjuYm5s/9wkjIyNxcnJ6rH3OnDlMmDCBpKQkatWqpV+vZvLkyWi1Wr766iucnZ2ZO3cuAMOHD0er1eLr64utrS1z5swBoHfv3kyaNAlfX18sLCyYNWvWc2fOT2oZ/BdCFGMqJaeBj0c0bdqUcePG0bp164LKZDSRkZG0aNGCPXv24OrqWiDnjFi5g/Cvg3h9+2zM7WQJBSFE8WKwJ/Phhx+yatUqHBwcsq0j8zxjMuIftrX+GZcp9bpM1yOEKF4MFpkZM2YA0KtXr2ztzzomI7JTu1cEpMgIIYong0Vm8ODBT7xtWDwfM1sbrCs6yuC/EKJYMlhkhg4dCmTdFQZQunTp/E30ElLXdOPu8UvGjiGEEC+cwWll7t69S79+/WjcuDGNGzdmwIAB3Lt3ryCyvTTUHm6kxd0lNfaOsaMIIcQLZbDIzJo1i+PHj9O9e3e6devG4cOHZfnlF+zRwX8hhChODF4u27dvH2PGjKF79+4AeHh4sHDhQqZOnZrv4V4WJapXAFMTEs+HU6app7HjCCHEC2OwJ5OcnKyfYwzA3Nyc5OTkfA31sjG1sqBEFRcZ/BdCFDsGezINGzZk5syZXLt2DUVR+PHHH59q0TKRN+qabtzadxJFUeRuPiFEsWGwJzNhwgTKlSvH119/zZIlSyhdujRjx44tiGwvFVsPNzISkki5Hmd4YyGEKCIM9mScnJzYtGkTFy5cQKVSUaNGDfmbdj54dEZma1cHA1sLIUTRkGuRWblyJa1atWLHjh3Z2g8fPgygn8BSvBg2VV0wsTDn3vlwHN7RGDuOEEK8ELkWmS+//JLq1avz5ZdfPvaeSqWSIvOCmZiZUqK6q9zGLIQoVp7Yk6lZsyYrV64syDwvNbWHGze3HkLRZaIyNThcJoQQhV6uv8kaNGiAnZ0dGzdupEyZMjRo0IAGDRpgYWHBli1bCjLjS8PWw43M5FTuh0cbO4oQQrwQufZkFixYQGJiIoGBgSQlJeHs7AzA5cuXOXHihDyMmQ8eHfwvUeXx1TyFEKKoybXImJiYsGLFClQqFTt37sz2XosWLfI92MvIuqIjpjZWJJ4Pp5yvPIskhCj6ci0y/fv3p0OHDvTp0wetVkvt2rWBrCf+HR0dCyzgy0RlYoK6ZgV58l8IUWzkOiZjbm6Oq6sre/fupXLlyty4cYOoqCiuXLkiNwPkI7WHG0mXr5OZnmHsKEII8dwMPow5bdo01qxZ81h7Xm5hTkxMpHv37nz99de4urpy4sQJpk+fTlJSEu7u7syYMQMLCwvOnj3LpEmTSE9Px9nZmdmzZ2NnZ0dCQgIjR44kIiKC0qVLExAQgIODA2lpaYwfP54zZ85gZWXFnDlzqFq1KoqiMGvWLH799VdMTEz4z3/+w2uvvfZs34yR2HpUQknPIOnydWwfjNEIIURRZfA+2S1btuDr60uVKlUYM2YMdevWpV27dgYPfOrUKXr06EFYWBiQVXCGDh3K1KlT2bp1KwA//fQTAP7+/gwbNoxNmzZRuXJlli1bBkBAQAAajYZt27bRtWtX/P39AVi1ahXW1tZs27aNcePGodVqAdixYwehoaEEBwezaNEitFotGRlFq0fw6OC/EEIUdQaLTFJSEr6+vjRq1IgKFSrw/vvvc+DAAYMHXr9+PZMnT9aP3xw8eJB69epRs2ZNIGtOtHfeeQeAzMxMkpKSgKxZn62srAAICQnRF7S2bduyf/9+0tPTCQkJoX379gDUr1+f27dvExUVxb59+2jTpg0mJiZUrlwZFxcXTpw4kWO+hIQEIiMjs72io41/67ClU2nM7NUyLiOEKBYMXi6rUKECa9eupV27dnz77bfY2NiQmppq8MAPex0PhYeHY2Njw+DBg7l27RoajUbfA9FqtXzwwQd8+eWXWFtbs379egBiYmJwcMiax8vMzAy1Wk18fHy2dgAHBweio6OJiYnJdlPCw/acrFixgoULFxr8HAVNpVJh6+EmPRkhRLFgsCczcuRI0tPTad68OTExMRw8eBBfX9+nPpFOp+PAgQNotVoCAwNJTk5m6dKlpKSkMH78eFasWMGBAwfo2bMnY8aMyT2wSc6RTUxMUBQlz9v37duXPXv2ZHvlNPZkDGoPN+6H3UCXbLiYCyFEYWawJ9O8eXOaN28OQGBgIBEREdSqVeupT1S2bFk8PT2pUKECAD4+PqxevZpLly5haWnJK6+8AsB7773H//73PwAcHR2Ji4vDycmJjIwMEhMTsbe3x9HRkdjYWNzcssYvYmNjcXR0pFy5csTGxurP+bA9J3Z2dtjZ2T315ygIth5ukKmQeDGCkvWqGTuOEEI8M4M9mbi4OL766itmzJjBwoULCQoKYsaMGU99oiZNmnD27Flu3LgBwK+//krt2rVxc3MjOjqaK1euALBnzx7q1q0LgLe3N4GBgQAEBwej0WgwNzfH29uboKAgAI4dO4alpSUuLi40bdqUzZs3o9PpCA8PJywsTH+sokQG/4UQxYXBnsyQIUM4efJktjaVSqUfT8krZ2dnpk6dyoABA0hNTcXDw4MxY8ZgbW3N9OnT+fTTT1EUhTJlyuhnfh4+fDharRZfX19sbW2ZM2cOAL1792bSpEn4+vpiYWHBrFmzAGjdujWnT5/W3xTg7++vv4mgKLEobYdluVIy+C+EKPJUSk4DGY/w8vKiT58+vPvuu9nGN8qXL5/v4QpaZGQkLVq0YM+ePbi6uho1y/lxS0n6+zqaDVOMmkMIIZ6HwZ5M/fr1sbCw0I+liIKhrunGrZCTpCckYW5XwthxhBDimRgsMh07duTzzz8nKCiIEiWyftmpVCp++eWXfA/3MrOt9c+4TKnXn/5GCyGEKAwMFpn58+djYmJCWlpakXt6vihTu1cEpMgIIYo2g0Xm3r17+Pv707FjxwKIIx4ys7XBuqKjDP4LIYo0g7cw+/n5sW/fPmJiYkhMTNS/RP5T15Qn/4UQRZvBnszy5ctJSkpi+/bt+jaVSsW5c+fyNZjIel4mdudRUmPvYOlgb+w4Qgjx1AwWmYeTWIqC93DwP/630zh3bmrkNEII8fQMFhlTU1N69eqlnz1ZFBzbWpWxq1eNK/N/okR1V+zqVjF2JCGEeCoGx2R27drFtWvXCiKL+BeVqQkeX/bHsqw957VLSLlxy9iRhBDiqRjsyTRs2JDZs2dz7tw57O3t9U/952VlTPH8zO3V1Jo9kFP9Z3NuzNd4fv05pjZFb6ocIcTLyeC0MjldJlOpVJw/fz7fQhlLYZpW5t9u/3GWsyMXU7pRXTxm9EeVyxIGQghRmBjsyaxcubIgcggDSr1RmyrD3uVKwAbCv95EpUEdjR1JCCEMMlhkGjRowLVr19i5cycqlYqWLVvKPGZG4tz1Le6HRxO5eifWlZwo1+YNY0cSQognMlhkjh8/Tr9+/UhOTgZg8eLFfPfdd/pFxkTBUalUVPmsG8nXYrg8Yw1W5ctS0lMWNRNCFF4GL+zPnj2b6tWrs379etatW0e1atWeadEy8WKYmJlS0/8jLJ3LcH7sUrnjTAhRqBksMhcuXKB///688soreHp68vHHHxfLQf+ixNyuBLVmDUTR6Tg3cjEZScnGjiSEEDkyWGRKlSrFoUOH9D8fOnSIUqVK5WsoYZiNWzk8pn3M/Ws3uThpOYou09iRhBDiMQaLzLvvvsuaNWvw9vbG29ubtWvX0rlz54LIJgywr1+Tqp914/ahs1xdtNHYcYQQ4jEGi8zAgQMZPXo0pUqVokyZMgwfPpxBgwbl+QSJiYm0bduWyMhIAE6cOEG3bt3w9fVlxIgRpKWlAXDlyhV69+5N+/bt6devH3fv3gUgKioKPz8/WrduzcCBA0lKSgIgISGB/v374+Pjg5+fH7GxsQCkpaUxatQofHx86NSpE6GhoU/3jRQxzp2b4vyuN1E/7iF600FjxxFCiGxyLTL9+vXj9OnTLFq0iJYtWxIYGMgvv/zCgAED9E/9G3Lq1Cl69OhBWFgYkFVwhg4dytSpU9m6dSsAP/30E4qiMHDgQD7++GM2bdqEh4cHS5cuBWDKlCn07NmT7du3U6dOHRYvXgxAQEAAGo2Gbdu20bVrV/z9/QFYtWoV1tbWbNu2jXHjxqHVap/5yykqqgx7F/sGHoTO/oE7xy8ZO44QQujlWi2OHz/O1q1bWbRoEZs3b2bPnj3ZXnmxfv16Jk+ejKOjIwAHDx6kXr16+lkEJkyYwDvvvMPZs2exsbGhadOsmYYHDBiAn58f6enpHD16lFatWgHQuXNn/ZIDISEhtGvXDoC2bduyf/9+0tPTCQkJoX379gDUr1+f27dvExUV9SzfTZGhMjOl5n/6YeXqwIVx35AcGWvsSEIIATzhOZkaNWqwYsUKVCoV8+fP17cripLnaWUe9i4eCg8Px8bGhsGDB3Pt2jU0Gg1arZajR49StmxZxowZw7lz56hRowYTJ07k9u3bqNVqzMyyYjo4OHDz5k0AYmJicHBwyPoQZmao1Wri4+OztT/cJzo6GhcXl2xZEhISSEhIyNYWHR1t8DMVVma2NllznH00m3OjFuP5zWjM1NbGjiWEeMnlWmTmz5/PH3/8wZgxY+jXrx/Vqj3/Q386nY4DBw6wbt06XFxcGD9+PEuXLsXNzY0jR46wevVq6tatS0BAADNmzOCzzz577BgqlSrX4+d2GS+n9hUrVrBw4cJn/zCFkLWrIzW//Jizw+dzYeIyas8eiMrM1NixhBAvsVyLTLly5ejQoQPOzs54eHhga2v73CcrW7Ysnp6e+mlpfHx8WL16NRqNBjc3N+rWrQtkXf4aNmwYpUuXJjExEZ1Oh6mpKbGxsfpLb46OjsTFxeHk5ERGRgaJiYnY29vj6OhIbGwsbm5ZC349us+j+vbtS6dOnbK1RUdH4+fn99yf05jsX61B1VE9uDxjDVcW/EzVz7oZO5IQ4iVmcATf0tKSTz/9lM6dO9OpUyc6der0zLcwN2nShLNnz3Ljxg0Afv31V2rXro2Xlxfx8fFcuHABgL1791K7dm3Mzc3RaDQEBwcDEBgYqB+38fb2JjAwEIDg4GA0Gg3m5uZ4e3sTFBQEwLFjx7C0tHzsUhmAnZ0drq6u2V5OTk7P9LkKG6f2jXF5rzk3NoRwY+N+Y8cRQrzEDE7137ZtW0JDQylbtizm5ub69r179+b5JM2bN2flypW4uroSEhLCvHnzSE1NxcPDgy+//BJra2tOnTrFf/7zH5KTk3FycmLWrFmUKVOG69evo9VquXXrFs7OzsydO5eSJUty584dtFotERER2NraMmfOHFxdXUlNTWXSpEmcOXMGCwsLpk2bRu3atfOUszBP9f+0FF0m50Z/xe0j56kzbwj2GlnZVAhR8AwWGS8vL0aOHFnkLyPlRXEqMgAZScmc7j+H1Lg71PtmNNYVyxk7khDiJWPwclmLFi2IiYkpiCziBTMrYU2t2QNRmZpydtRXZCTcN3YkIcRLxuBU/7dv32br1q0EBwejVquBrDu8fvnll3wPJ56flUtZPKb358zQ/3F+wjfUnjsEE7njTAhRQAz2ZK5evYqzszM6nY67d+9y9+5d7ty5UwDRxItS0rMa1cb05O6xi1yZtx4DV0iFEOKFMdiTeZoBflF4lfNtyP2waK6v2YVNZWdc3n3L2JGEEC+BXIvMypUrc91JpVLRu3fvfAkk8k+lAR1IDr/JlYANWFdwpNTrtYwdSQhRzOV6d9nD+cVy3CmP08oUNcXt7rKcZCSlcHrgf0mNjsdr5XisnEobO5IQohh7pp6MKLrMSlhRa/onHO89jdDZP1BrzqAnTtUjhBDPI9ci06BBg4LMIQqQVfmyVBrQgSsBG4jdfgRHn9eNHUkIUUzlbWEYUew4d/HGtm4VrvzvJ9LiEwzvIIQQz0CKzEtKZWpC9bF+6JJTCf3vOmPHEUIUU1JkXmI2lZyp+GEbbv16griQE8aOI4Qohgw+JxMeHs6yZcu4ffu2/iE+lUrFggUL8j2cyH/l/d4hbu9xQv+7DvtX3TGzszF2JCFEMWKwyIwaNYrTp09na5O7kYoPEzNTqo/rzcmPZnJl/k/UmNDH2JGEEMWIwctlly5dYsiQIfz1119cuHCBCxcuFMtnZF5mavcKuPq9Q0zwH9z+45yx4wghihGDRaZBgwao1epsa8mI4qfiB22wdivH5ZlryUhKMXYcIUQxYfBy2cMFwUJCQrLNwixjMsWLiaU51cf14vSAuYQvCaLqiPeMHUkIUQwYLDJbt24F4I8//tC3yZhM8WRXtyrO73pzY0MIZVu8RknPasaOJIQo4gwWmQsXLhREDlFIVPqkPfEH/uLy9NXUWzEOU0sLY0cSQhRhBotMWloav//+u34NGZ1Ox99//41Wq83vbMIITG2sqDamJ2c/XUDEsmAqDepo7EhCiCLM4MD/2LFjGTBgAGPHjkWr1TJhwgTWrFmTp4MnJibStm1bIiMjAThx4gTdunXD19eXESNGkJaWlm37kJAQmjdvrv85ISGB/v374+Pjg5+fH7GxsUBW4Rs1ahQ+Pj506tSJ0NBQABRFYebMmbRu3Zo2bdrw559/5u1bENmUauBBubYNifxhN4kXrhk7jhCiCDNYZPbv30/fvn2pXbs2//3vf2nSpAk+Pj4GD3zq1Cl69OhBWFgYkFVwhg4dytSpU/XjPD/99JN++7i4OGbOnJntGAEBAWg0GrZt20bXrl3x9/cHYNWqVVhbW7Nt2zbGjRun71Xt2LGD0NBQgoODWbRoEVqtloyMjLx9EyKbykO7YFHKlr+nryYzQ2fsOEKIIspgkUlNTaV+/fo0aNAAlUpF586d2bdvn8EDr1+/nsmTJ+Po6AjAwYMHqVevnn6dmgkTJvDOO+/ot58wYQJDhgzJdoyQkBDatWsHQNu2bdm/fz/p6emEhITQvn17AOrXr8/t27eJiopi3759tGnTBhMTEypXroyLiwsnTsh0Kc/CzNaGqiO7k/R3JJGrdho7jhCiiDI4JlOtWjW+/fZbPvjgA+bPn4+lpWWeDvyw1/FQeHg4NjY2DB48mGvXrqHRaPQ9kJUrV1KrVi08PT2z7RMTE4ODg0NWUDMz1Go18fHx2doBHBwciI6OJiYmRl/UHm3PSUJCAgkJ2Wcfzm3bl1WZpp6UbfEaEd9vo+xb9bCp7GzsSEKIIsZgT2bSpEm4ubnRrFkzrK2tOXfuHL169XrqE+l0Og4cOIBWqyUwMJDk5GSWLl3KpUuX2LlzJ4MGDcpbYJOcI5uYmJDTIp+5bb9ixQpatGiR7eXn55f3D/SSqPJZN0xtLPn7y9UoukxjxxFCFDEGezL16tWjXr16QNYlsHv37lGyZMmnPlHZsmXx9PSkQoUKAPj4+LB69WoURSE2NpYuXbqQnp5OTEwMPXv2ZO3atTg6OhIXF4eTkxMZGRkkJiZib2+Po6MjsbGxuLm5ARAbG4ujoyPlypXT3xzwaHtO+vbtS6dOnbK1RUdHS6H5F4vStlT5tCuXpnxP1E8hlH+vueGdhBDigVyLzJAhQxg2bBjz589/7L1neeK/SZMmLFiwgBs3buDs7Myvv/5K7dq1GTZsGMOGDQMgMjKSPn36sHbtWgC8vb0JDAxkwIABBAcHo9FoMDc3x9vbm6CgIDQaDceOHcPS0hIXFxeaNm3Kzz//rL+jLSwsjLp16+aYx87ODjs7u6f6DC8rh5b1id11jPAlmyjT5BWsypc1diQhRBGRa5HZvXs3fn5+7N69+7H3nuWJf2dnZ6ZOncqAAQNITU3Fw8ODMWPGPHGf4cOHo9Vq8fX11U9vA9C7d28mTZqEr68vFhYWzJo1C4DWrVtz+vRp/U0B/v7+WFlZPXVWkZ1KpaLa6B4c7/kf/p6xhjrzh8msD0KIPFEpOQ1kvKQiIyNp0aIFe/bswdXV1dhxCp0bgb8ROusHqo3piVOHJsaOI4QoAvK0MmZERAQA27dvZ/ny5Y89RCleDk7tG1Py1epcXfgLqbF3jB1HCFEEGCwyCxYsYOTIkZw7d45PP/2U2bNnM23atILIJgoZlYkJ1bS9UDJ0hM7+Ice7+YQQ4lEGi8z69eupW7cuu3bton79+nz22Wfs2bOnILKJQsja1QG3/lmTaMbtOmbsOEKIQs5gkbl79y7e3t78+eefvP7661SvXp2kpKSCyCYKKZduzVDXqkRowAbSb98zdhwhRCFmsMi4uLiwdOlSjh8/jpeXFwEBAVStWrUgsolCSmVqQvVxvdAlJnMlYIOx4wghCjGDRWbYsGFcunSJli1b0qhRI6KjoxkxYkRBZBOFWIkqLlR434fYXce4deC0seMIIQqpp76F+d69e9ja2uZXHqOSW5ifTmZ6Bic/nEHG3SReXTMRM1sbY0cSQhQyBnsy8fHxDB8+nFdffZXTp08zcuRImUhSAGBibkb1cb1Ji0/g6sKNxo4jhCiEDBaZadOmce7cOZKTk0lISODkyZOMHz++ILKJIsDWw43yPd7m5uaD3DkmS3ULIbIzWGQOHjyonyHZzs6O8ePHc/z48XwPJoqOih/5YlXBkb9nrEGXnGrsOEKIQsRgkbGystLPbKzT6Th58iRly8oEieIfppYWVNf6kRp1iysBG1BkJU0hxAMGp/p/7733mDt3LpA1MaVOp2P48OH5HkwULSW9qlPe7x2ur9lFUmgUNSb2xcatnLFjCSGMzGCRGTRoEC4uLoSEhJCRkcGbb75Jt27dCiKbKGIqD+6E2r0CobN/5OT7X1JpcCecOzdFlcvCcUKI4s9gkZk7dy5du3alY8eOBRBHFHUOb2uw86zG5RlruDJ3PfG/nab6uF5Ylitt7GhCCCMw+FfMpUuX0qpVKz766CN2796NTifX28WTWTrYU2vOIKqN7knCmasc7+1PzPbDMqGmEC8hg0Vm+/bt9O/fn7CwMIYOHcpbb71FQEAAUVFRBZFPFFEqlQqnjk3wWjGOElVcuDR1BRfGf0v6nURjRxNCFKA8P/GfmJjIl19+yS+//AKAubk5kyZNomvXrvkasCDJE//5Q9Flcn3tbsK/2YyZXQmqj/WjdOOcl8UWQhQvBovMH3/8wcaNG9m5cyfJycnUrVuXjh07snHjRm7evMlvv/1WUFnznRSZ/JV0OZKLU1dw//J1yrVrROVhXTArYW3sWEKIfGRw4P/999/H0tKSNm3a0LNnT+rWzfobaOXKlRk4cGC+BxTFR4lqrtT7djTXlgcTuXond45dpMaEPpT0qm7saEKIfGJwTGbkyJHs27eP6dOn6wsMQO3atfWXzp4kMTGRtm3bEhkZCcCJEyfo1q0bvr6+jBgxQr+U8+7du+nQoQPt27dn0KBB3L17F4CoqCj8/Pxo3bo1AwcO1K9lk5CQQP/+/fHx8cHPz0//wGhaWhqjRo3Cx8eHTp06ERoa+pRfichPJhbmVBrQgVcWj0BlYsJfQwK4uuBnMlPTjR1NCJEPDBaZzp07s3XrVhYtWsTChQuZP38+Q4YMoWTJkgbXlTl16hQ9evQgLCwMyCo4Q4cOZerUqWzduhWAn376icTERL744guWLl3Kpk2bcHd3Z8GCBQBMmTKFnj17sn37durUqcPixYsBCAgIQKPRsG3bNrp27Yq/vz8Aq1atwtramm3btjFu3Di0Wu0zfzki/9i9UhWvFWNx6tiE6z/s4eSHM0i8GGHsWEKIF00x4KOPPlLc3d2VmjVrKu7u7oq7u7vSqFEjQ7spiqIo48aNU44ePao0a9ZMiYiIULZv364MHjxY//6tW7eUmJgYJT4+Xtm1a5e+fdu2bcqgQYOUtLQ0xcvLS0lPT1cURVGioqKU5s2bK4qiKM2aNVOioqIURVGU9PR0xcvLS0lLS1N69eqlHD16VH+sFi1aKNevX38s2927d5WIiIhsr6NHjyo1atRQIiIi8vT5xItx6/czyuF2WuVAk8FK+PJgJTM9w9iRhBAviMExmWPHjjFixAgOHDjA4MGD+emnn7h//36eCtjD3sVD4eHh2NjYMHjwYK5du4ZGo0Gr1WJpacnbb78NQEpKCkuXLqV3797cvn0btVqNmVlWTAcHB27evAlATEwMDg4OAJiZmaFWq4mPj8/W/nCf6OhoXFxcsmVZsWIFCxcuzNPnEPmrdMPaeK2eQOicdVz7ZjO3f/+LGhP7Yl1RpqURoqjL03wfFStWxMvLixs3btCsWTP++OOPZzqZTqfjwIEDaLVaAgMDSU5OZunSpfr37927x8cff0zNmjXp1KlTjg/vqVSqXI9vksv0JTm19+3blz179mR7rVmz5hk+lXgRzO1KUHPqh7hP+ZDkazGc6PslUT/vQ8nMNHY0IcRzMNiT8fDwYNGiRYwcOZLx48djYWGBlZXVM52sbNmyeHp6UqFCBQB8fHxYvXo1kNUz6devH2+88Qbjxo0DoHTp0iQmJqLT6TA1NSU2NhZHR0cAHB0diYuLw8nJiYyMDBITE7G3t8fR0ZHY2Fjc3NwAsu3zKDs7O+zs7J7pc4j84/COBrt61bg8fTVX/ruO+P2nqD6+N5aOpYwdTQjxDAz2ZKZPn06zZs148803qVevHnfv3n3mWZibNGnC2bNnuXHjBgC//vortWvXRqfTMWDAAHx8fBg/fry+t2Jubo5GoyE4OBiAwMBAmjZtCoC3tzeBgYEABAcHo9FoMDc3x9vbm6CgICDrUp+lpeVjl8pE4WbpYE+t/w6m6ugeJPx1heO9phG2OJD4g3+RcS9vl2qFEIVDnp/4fx7Nmzdn5cqVuLq6EhISwrx580hNTcXDw4Mvv/ySAwcOMHToUNzd3fX71KlTB39/f65fv45Wq+XWrVs4Ozszd+5cSpYsyZ07d9BqtURERGBra8ucOXNwdXUlNTWVSZMmcebMGSwsLJg2bRq1a9fOU055GLPwSY6MIXTOOu7+eRFFlwkqFTZVnCnpWQ27V6pi51lVJt8UohDLtcgkJSXx9ddf4+DgQPfu3Rk9ejT79u3Dw8OD2bNnU758+YLOmu+kyBReupQ07p29SsLpUBJOhXLvzBV097NW4bR0Kq0vOHae1bCp5CTLCwhRSORaZD7//HO2bt2KSqWiTp06nD9/nho1anDp0iXq1q3LDz/8UNBZ850UmaJDydCRFHqdhFOh3D11mYRToaTHJwBgZmuD7StVsno7nlVRu1fExMLcyImFeDnlOvC/b98+Pv74Yzw8PBgxYgQTJ07Ez8+P7777joCAgAKMKMTjVGamqN0ronaviEu3ZiiKQsr1WBJOhT7o7Vzm9sEzQNYsA+pabg96O9Wwq1sFM7XMmSZEQci1yCQlJdGoUSO8vLwAqFmzJgA1atTQTwUjRGGhUqmwdnXE2tWRcr4NAUiLv0fCX6EknLxMwulQItfsgpU7QKWiRLXy2NWtgk1lZ6wrlsPK1QHLcqXkMpsQL1iuRUZRFFQqlf4ZE3PzrMsNuT2LIkRhY1HalrLe9SjrXQ8AXXJq1rjOg95OzPbD+nEdyOrxWJUvi3VFR6xcHbGu8M/LvIzdE5/REkLk7InPyfzwww/s3r0blUrFihUrKFu2rP72YyGKGlNrS+w1NbHXZPXKFUUhLe4uyddukhIRS3LETZIjYrkfFk38wTMoGf+sAmtqY/lY4bGq4Ih1RUfM7UoY6yMJUeg9scjs2LFD/+eHE1rCk5+6F6KoUKlUWDrYY+lgD6+5Z3tPydCRevM2yRExJEfGkHwthpTIGBLPhxP363HI/Od+GTO7ElhXcMgqOhUcsaniQqn6NTG1ebaHloUoTnItMnv27CnIHEIUKiozU6zKl8WqfFlKUSvbe5npGaRcj8tWfJKvxXD3+CVitx/J2t/CjFINalGmWT3KNH4FMzsbY3wMIYwu1yJTHJ+DEeJFMDE3w6aSEzaVnB57T5eSxr1zYdzad5Jb+04Sf+A0l01NKPlqDco086LMm69gUaakEVILYRwF8sR/USHPyYgXScnMJPHCNW6FnCQu5AQpkbGgUmH3SlXKvFWPMt71sHKS2QpE8SZF5hFSZER+URSF+1eiiPv1BLf2neR+aBQA6poVKfNWPcq+5SVLG4hiSYrMI6TIiIKSHBGT1cPZd5LEc2EA2FR2zrqk5l2PEtXKyw02oliQIvMIKTLCGFJvxnNr3yniQk6QcCoUFAWr8g6UecuTsm95ofZwk4dERZElReYRUmSEsaXFJxD/22niQk5y99gFFF0mFg72lPH2xO6VqpjZlcDcrgRmtjaY2dpgqraSApSLzNR00u8mkn47kfS7iWTcSST9Ttaf0x/8OePB+4pOh4mlBSaW5o/8M+tlammBiVX2tofbmFqYP/Kexb/2Mce8pBqVmamxvwqDHpaB/Og9S5F5hBQZUZhkJNwn/uBfxO07yZ0/zpGZlv74RioVpmorfdHJ8WWXS1sJa1SmL6ZAKYoCigLKo39+/GflQds/72Xqf87tGIp+e4XMtIwcCkYS6XfukXE3ifTb90i/m0TG3cRsszn8+zszK1kC85IlMLdXY26fVQgyU9P1L11qWrafM1PSyExNy1pu4mmYmmQ9i+VcBivnMlg6lcbK6cE/nctgUa4UJgVQhJTMTNJuJZAaHa9/pUTf+ufnm/GUqObKK19//sLPbXBlTCGEcZjZ2eDo8zqOPq+jS04lJeoWGffuP/5KyP7z/bgb+j8raRlPPIep+kGhefiLHiXrQVMF4MEv+MysQqDwoD0z88HbjxQMIzKxtsxWMKzdnPR/Ni+pxqxU1j8ftpnZ2jxzcVUydOj0xefRQpT2eHFKTiUt7i4p0fGkRt/izrELpMXezf59maiwKFvyQQEqg5VT6WwFybJcqTzNIJ6ZnkFazJ1shSPlkYKSGnMbJT37fwtmtjZZxc7VgZIad+zr13ym78QQKTJCFAGm1paUqPr0K7zqUtMeK0L//pkHi8FlvR5cMnlw2URl8k87qP75mQfb5fQzqqymB8dRPdxf/96Dn7O9Z/LgFA8zmOi3wUSFigft5mbZCoqZvRpTS4vn/HbzTmVmipmZKZR4ttkc9MXgxq1svYmUG7dIOHWZ2N13sv59PMKibEl9z8fSqTTmpe1Iv5WQVURuZhWRtLh/FS+VCosydliWK62/g9HKqXRW4XrQkzJ7xs/wtKTICFGMmVpaYOpgkTV1jjA6E3Mz/UwSOVEydKTG3iE1+lZWEbkRry9I986GEbf3OIouE5WZKZaOpbB0Ko19/ZrZLsNZOpXG0tG+0KyhJEVGCCEKCZWZKVYPLpflNC+Eossk497957rkV9DyPWViYiJt27YlMjISgBMnTtCtWzd8fX0ZMWKEfm2a8+fP06VLF1q1asX48ePJyMi6fhgVFYWfnx+tW7dm4MCBJCUlAZCQkED//v3x8fHBz8+P2NhYANLS0hg1ahQ+Pj506tSJ0NDQ/P6IQghRIFSmJlk3KhSRAgP5XGROnTpFjx49CAsLA7IKztChQ5k6dap+VueffvoJgFGjRjFx4kR27NiBoiisX78egClTptCzZ0+2b99OnTp1WLx4MQABAQFoNBq2bdtG165d8ff3B2DVqlVYW1uzbds2xo0bh1arzc+PKIQQ4gnytcisX7+eyZMn4+joCMDBgwepV6+efpXNCRMm8M4773D9+nVSUlKoV68eAJ07d2b79u2kp6dz9OhRWrVqla0dICQkhHbt2gHQtm1b9u/fT3p6OiEhIbRv3x6A+vXrc/v2baKiovLzYwohhMhFvo7JPOxdPBQeHo6NjQ2DBw/m2rVraDQatFot586dw8HBQb+dg4MDN2/e5Pbt26jVaszMzLK1A8TExOj3MTMzQ61WEx8fn6394T7R0dG4uGS/MychIYGEhIRsbdHR0S/uwwshhCjYgX+dTseBAwdYt24dLi4ujB8/nqVLl9K4cePHtlWpVOT0nOiTnkjNbWnonNpXrFjBwoULnyK9EEKIp1WgRaZs2bJ4enpSoUIFAHx8fFi9ejWdO3cmLi5Ov11sbCyOjo6ULl2axMREdDodpqam+nYAR0dH4uLicHJyIiMjg8TEROzt7XF0dCQ2NhY3N7dsx/q3vn370qlTp2xt0dHR+Pn55dfHF0KIl06B3qLQpEkTzp49y40bNwD49ddfqV27NuXLl8fS0pI///wTgMDAQJo2bYq5uTkajYbg4OBs7QDe3t4EBgYCEBwcjEajwdzcHG9vb4KCggA4duwYlpaWj10qA7Czs8PV1TXby8np8UWohBBCPLsCmbusefPmrFy5EldXV0JCQpg3bx6pqal4eHjw5ZdfYm1tzYULF5gwYQJJSUnUqlWL6dOnY2FhwfXr19Fqtdy6dQtnZ2fmzp1LyZIluXPnDlqtloiICGxtbZkzZw6urq6kpqYyadIkzpw5g4WFBdOmTaN27dp5yhkeHk7Lli1Zs2aNFBwhhHgKTk5O+vHzR8kEmY84duyYXC4TQohnkNvEwlJkHpGSksKZM2dwcHDA1LRwTM/9cJyoKPSuilJWKFp5i1JWKFp5JeuLkVtPRqaVeYSVlRUajcbYMXLk5ORUZJYfKEpZoWjlLUpZoWjllaz5o+jMTSCEEKLIkSIjhBAi30iREUIIkW+kyBRydnZ2DBkyBDs7O2NHMagoZYWilbcoZYWilVey5i+5u0wIIUS+kZ6MEEKIfCNFRgghRL6RIlOILVy4EF9fX3x9fZk1a5ax4+TJzJkzi8RCcXv37qVz5860bt2aadOmGTvOEwUFBen/O5g5c6ax4+Tq36vg/v7777Rr146WLVsyb948I6fL7t9Z161bR9u2bWnXrh1jx47Vr9hbGPw760Nr1qyhd+/eRkqVd1JkCqnff/+dAwcOsHHjRgIDAzl79iy7du0ydqwnOnToEBs3bjR2DIMiIiKYPHkyixcvZvPmzZw7d459+/YZO1aOkpOT8ff3Z9WqVQQFBXHs2DF+//13Y8d6zL9XwU1JSWHcuHEsXryY4OBgzpw5U2i+439nvXr1KsuWLePHH39k06ZNZGZmsnbtWuOGfODfWR+6fPkyS5YsMU6opyRFppBycHBAq9ViYWGBubk5VatWLdQrfN65c4d58+YxYMAAY0cxaNeuXbRp0wYnJyfMzc2ZN28enp6exo6VI51OR2ZmJsnJyWRkZJCRkYGlpaWxYz3m36vgnj59Gjc3NypUqICZmRnt2rXTr2prbP/OamFhwRdffIFarUalUlGjRo1C8//av7MCpKWlMWnSJIYPH27EZHkn08oUUtWrV9f/OSwsjODgYH788UcjJnqySZMm8dlnn+mXcSjMwsPDMTc3p1+/fsTGxtKsWTM+/fRTY8fKkVqtZvjw4fj4+GBlZUWDBg149dVXjR3rMf9eBfffK9Q6OjrqV7U1tn9nLV++POXLlwcgPj6eNWvWMH36dGNEe8y/swL897//pUuXLjKtjHgx/v77bz788EPGjBlDpUqVjB0nRxs2bMDZ2ZmGDRsaO0qe6HQ6Dh06xOzZs1m/fj1//fVXob3Md+HCBX7++Wd+/fVXDhw4gImJCcuWLTN2LIOedlXbwuDmzZv07duXLl268Prrrxs7To4OHjzIjRs36NKli7Gj5JkUmULszz//5P333+fzzz9/bBXPwiQ4OJiDBw/SoUMH5s+fz969e/nyyy+NHStXZcuWpWHDhpQuXRorKytatGjB6dOnjR0rRwcOHKBhw4aUKVMGCwsLOnfuzJEjR4wdy6By5cplW+02JiYmxxVqC4vQ0FB69OhBp06dGDx4sLHj5GrLli38/fffdOjQgQkTJnDmzJlC2wt/SC6XFVI3btxg8ODBzJs3r9D3EL777jv9n3/55ReOHDnCuHHjjJjoyZo1a8aYMWNISEigRIkS/Pbbb7Ro0cLYsXJUs2ZNZs+ezf3797G2tmbv3r3UrVvX2LEM8vT05OrVq4SHh+Pq6sqWLVsK7d++ExMT6devH5999hkdOnQwdpwnevQy3uHDh1m4cCEBAQHGC5QHUmQKqWXLlpGamsqMGTP0bd27d6dHjx5GTFU8eHp68tFHH9GzZ0/S09Np3Lhxof0F2KRJE86dO0fnzp0xNzenbt269O/f39ixDLK0tGTGjBkMHTqU1NRUvL29ad26tbFj5einn34iLi6O5cuXs3z5ciBrNd+iMrBe2Mm0MkIIIfKNjMkIIYTIN1JkhBBC5BspMkIIIfKNFBkhhBD5RoqMEEKIfCNFRhRakZGRuLu78+abb3Lv3j19u0ajeeGzz/7yyy+4u7uze/fuF3pcQ7799luaNGlCgwYN2L9//2PvX7x4kQEDBvDaa6/x2muv0atXLw4fPlygGY3h999/p1evXsaOIV4AKTKi0IuJiWH27NnGjpEv1q1bh6IozJgxg1deeSXbe6GhoXTv3p2rV68yceJEvvjiC2JjY+nXrx9nzpwxUuKC8dVXX3HhwgVjxxAvgBQZUeiZmpqyfv16jh49+th7/+6BPPrzw57Q5MmT6d69O/Xr12f06NGsXbuWZs2a0aBBA1auXJnteAcOHKBly5a88cYbzJo1Sz8H1549e2jTpg1eXl706NFD/0t+wYIFuLu7o9Vq8fLy4pdffnks448//kirVq3w9PSkQ4cO7NmzB4DevXtz7do14uLiGDhwIPb29tn2++abb7h//z6LFy+mY8eOtGvXjvnz5zNixAhsbW0BuHTpEu+//z716tXD29ubmTNnkpKSoj++r68v06ZN4/XXX6djx47s2LGDnj174unpyccff8z9+/eBrIcP+/Tpw4gRI/D09KRr165cvHgRgMzMTBYvXkyzZs3w8vKie/fu/Pnnn8A/vc0pU6bQvXt36tWrx/vvv8+tW7eArCnpe/fujZeXF61bt2bz5s0G91uwYAFHjhzh3r17uLu7AzBv3jyaNGlCnTp1aN26NTt27Mjzfz/CuKTIiEKvXbt2lClThokTJz7TYlJbtmzhvffew8vLi6CgIH744QdGjx6NWq1m1qxZZGRk6Lc9evQoI0eOpHHjxixbtoygoCDCw8MZPnw4devWZerUqVhbWzN48GD9L3OA2NhY/P39eeutt7Kde8eOHUyePBl3d3f++9//Uq5cOYYMGcKxY8cYM2YMZcuWpWLFitmm5nnozJkz2NvbU7VqVX2bu7s7H374IW5ubiQlJfHhhx8SFhbGtGnT8PPz4/vvv882c+/ly5fR6XQMHDiQ8+fPM2LECNq0aUObNm3Yv38/O3fu1G97+PBhypcvz5QpU7h69SrDhw9HURS+//57/ve//+Ht7c2sWbPQ6XT069cv2xongYGBdOnShfbt23Po0CE2bNhARkYGgwcPJi0tjSlTptC4cWNGjx6drYeS034dO3bE3d0dGxsbvvvuOy5cuMDXX3/NO++8w6JFi3Bzc+Pnn39Gp9M99X8LouDJtDKi0LOzs2PixIkMHz6cRYsWPXHbnCawaN26NZ06dSIjI4N9+/bRr18/fHx82Lt3L5s2bSI5OVm/7ccff0zLli1p0KABW7Zs4fDhwyQlJZGenk5gYCCBgYH6bS9duqT/84cffkjjxo0fO3dQUBCmpqbMnDkTa2trvLy8aNSoEZs3b2bKlClYWlpiY2NDo0aNHtv34VoyuTl06BCxsbFMnTqVtm3bAllF8uGxAczNzZk4cSKKojB9+nQ0Gg29evXi0KFD/PLLL9y9e1d/vPLly/P5558D8Ndff7F69Wqio6MJCgrC0dGRyZMno1KpKFeuHF27dmXnzp20adNG/x137doVjUbDunXriI+PJywsTF+ITp48qT/Pb7/9ho+PT677VahQgZIlSxIVFUWjRo1ITEykYsWKbNq0iRs3bvD666/Tpk0bTE1Nc/1uROEhRUYUCa1bt+add95h2bJl2QrJw+njH/ZGHi0YD9nY2ADofymp1epsPz96vPT0dABMTLI6+WZmZvpjT548mVq1anH//n0sLS2pVKmSfrXHf1/qeujhcR7mfHiuvEx7X6tWLbZs2UJoaKi+N7N//34WL17M8OHDczyGoijZ2q2srPQZAEqUKJEtV06f/dF8JiYm2fbP7TM8/E7Nzc312zz83rp27cq7775LWloamZmZVKxYUV88c9rv39RqNRs3buS3337j1KlTrF+/nvnz57Nly5Yis6bKy0yKjCgyJk2axOHDh0lISNC3lS5dGsha7bJUqVJs2LDhuc6xZMkSSpQoQUhICACNGjWievXqmJiYsG3bNsqUKcOSJUu4ceOGfmwFeOwX8UMtW7Zk165daLVa2rZty4YNGzA1NaV9+/YGs/Tv3589e/YwbNgwPvnkE3Q6HfPmzSMxMRFXV1dKlSqFvb29PvONGzc4ePAg7733Xq55niQmJoaZM2dSq1YtNm/eTLVq1XB0dKRly5YEBATwn//8h0aNGvHNN99gY2NDq1atnni8KlWq4OLiwm+//UajRo3Yu3cvW7du5YcffqBs2bJP3Nfc3JyUlBSCg4MpVaoUH330ES1btqRjx47ExsZy9epV4uPjpcgUATImI4oMR0dHRo8ena2tUaNG+Pj4sHPnTr744gv9ZZhn5efnx5w5c9i/fz8DBw7Ex8eHatWqsWDBAu7cucPo0aNJT09n7ty5+h7Sk7Rv356JEyfqx0Nu3rzJokWL8rS6pbu7O6tXr8bJyYlJkyYxbdo0qlatyooVK6hQoQJqtZrvvvuOChUqMH78eFatWsUHH3zA2LFjn+mzV6pUiaioKCZPnoybmxsBAQGoVCr69+/PkCFD2Lt3LyNHjkSlUrF8+XIqVqz4xONZWFjw7bff4ubmxoQJE/RLQNSrV89glg4dOmBtbc2sWbPw8PBg1KhRnD59msGDB3PixAnGjh372N14onCSWZiFEDRv3hxbW1uCgoKMHUUUM9KTEUIIkW+kJyOEECLfSE9GCCFEvpEiI4QQIt9IkRFCCJFvpMgIIYTIN1JkhBBC5BspMkIIIfLN/wFBuYs+rc4uhwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# generate training and testing data splits\n",
    "train_h, test_h = train_test_split(model_data_h, test_size = 0.2, random_state = 13)\n",
    "\n",
    "# select models with lowest BIC\n",
    "n_components = np.arange(1, 16)\n",
    "models_train = [GaussianMixture(n, covariance_type = 'full', random_state = 2).fit(train_h) for n in n_components]\n",
    "\n",
    "sns.set_style(\"ticks\", {\"xtick.major.size\": 7, \"ytick.major.size\": 7})\n",
    "\n",
    "plt.plot(n_components, [m.bic(train_h) for m in models_train], color = \"#C33969\", label = 'Training Dataset')\n",
    "plt.legend(loc = 'upper right')\n",
    "plt.xlabel('Number of Components', fontweight = 'bold')\n",
    "plt.ylabel('Bayesian Information Criterion', fontweight = 'bold')\n",
    "\n",
    "sns.despine()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.7576008985854212\n",
      "4.136912735720333\n",
      "-0.24279327730376613\n",
      "0.1649208802738954\n",
      "-0.00961128670731657\n",
      "0.16376505489652082\n"
     ]
    }
   ],
   "source": [
    "gmm_h = GMM(n_components = 8, random_state = 11)\n",
    "gmm_h.from_samples(train_h.to_numpy())\n",
    "\n",
    "Y_h = gmm_h.predict(np.array([0,1,2]), test_h[[\"purpleair_pm2.5\", \"avg_temp\", \"RH\"]].to_numpy())\n",
    "\n",
    "print(r2_score(Y_h, test_h[[\"embassy_pm2.5\"]].to_numpy()))\n",
    "print(mean_absolute_error(Y_h, test_h[[\"embassy_pm2.5\"]].to_numpy()))\n",
    "print(n_bias(Y_h, test_h[[\"embassy_pm2.5\"]].to_numpy()))\n",
    "print(cvMAE(Y_h, test_h[[\"embassy_pm2.5\"]].to_numpy()))\n",
    "print(mnbias(Y_h, test_h[[\"embassy_pm2.5\"]].to_numpy()))\n",
    "print(nmae(Y_h, test_h[[\"embassy_pm2.5\"]].to_numpy()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cross Validation\n",
    "\n",
    "We will now perform a 10 fold cross validation within the training data. We will find the minimum of the BIC for each fold, build a GMM, and then perform a regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "kf = KFold(n_splits = 10, random_state = 8, shuffle = True)\n",
    "\n",
    "n_components = np.arange(1, 11)\n",
    "\n",
    "CVs = []\n",
    "\n",
    "training_data = []\n",
    "\n",
    "for train_index, test_index in kf.split(train_h):\n",
    "    train, test = train_h.to_numpy()[train_index], train_h.to_numpy()[test_index]\n",
    "    \n",
    "    n = 8\n",
    "    \n",
    "    #GMR\n",
    "    gmm_ = GMM(n_components = n, random_state = 11)\n",
    "    gmm_.from_samples(train)\n",
    "    Y_ = gmm_.predict(np.array([0,1,2]), test[:, [0,1,2]])\n",
    "                     \n",
    "    r2 = r2_score(Y_, test[:, [3]])\n",
    "    mae = mean_absolute_error(Y_, test[:, [3]])\n",
    "    bias_gmr = n_bias(Y_, test[:,3])\n",
    "    cvMAE_gmr = cvMAE(Y_, test[:,3])\n",
    "    mnbias_gmr = mnbias(Y_, test[:,3])\n",
    "    nmae_gmr = nmae(Y_, test[:,3])\n",
    "    \n",
    "    #MLR\n",
    "    regr = linear_model.LinearRegression()\n",
    "    regr.fit(train[:,[0,1,2]], train[:,[3]])\n",
    "    mlr_pred = regr.predict(test[:,[0,1,2]])\n",
    "    \n",
    "    mlr_r2 = r2_score(test[:,3], mlr_pred)\n",
    "    mlr_mae = mean_absolute_error(test[:,3], mlr_pred)\n",
    "    bias_mlr = n_bias(mlr_pred,test[:,3])\n",
    "    cvMAE_mlr = cvMAE(mlr_pred, test[:,3])\n",
    "    mnbias_mlr = mnbias(mlr_pred, test[:,3])\n",
    "    nmae_mlr = nmae(mlr_pred, test[:,3])\n",
    "        \n",
    "    #Random Forest\n",
    "    RF_regr = RandomForestRegressor(max_features = 1, random_state = 0)\n",
    "    RF_regr.fit(train[:,[0,1,2]], train[:,[3]])\n",
    "    RF_pred = RF_regr.predict(test[:,[0,1,2]])\n",
    "    \n",
    "    RF_r2 = r2_score(test[:,3], RF_pred)\n",
    "    RF_mae = mean_absolute_error(test[:,3], RF_pred)\n",
    "    bias_RF = n_bias(RF_pred, test[:,3])\n",
    "    cvMAE_RF = cvMAE(RF_pred, test[:,3])\n",
    "    mnbias_RF = mnbias(RF_pred, test[:,3])\n",
    "    nmae_RF = nmae(RF_pred, test[:,3])\n",
    "                   \n",
    "    row = np.array([r2, mlr_r2, RF_r2, \n",
    "                    mae, mlr_mae, RF_mae,\n",
    "                    nmae_gmr, nmae_mlr, nmae_RF,\n",
    "                    cvMAE_gmr, cvMAE_mlr, cvMAE_RF,\n",
    "                    bias_gmr, bias_mlr, bias_RF,\n",
    "                    mnbias_gmr, mnbias_mlr, mnbias_RF])\n",
    "    CVs.append(row)\n",
    "\n",
    "CVs = pd.DataFrame(CVs, columns = ['GMR R2', 'MLR R2', 'RF R2',\n",
    "                                   'GMR MAE', 'MLR MAE', 'RF MAE',\n",
    "                                  'GMR nMAE', 'MLR nMAE', 'RF nMAE',\n",
    "                                  'GMR cvMAE', 'MLR cvMAE', 'RF cvMAE',\n",
    "                                  'GMR Bias', 'MLR Bias', 'RF Bias',\n",
    "                                  'GMR mnBias', 'MLR mnBias', 'RF mnBias'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "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>GMR R2</th>\n",
       "      <th>MLR R2</th>\n",
       "      <th>RF R2</th>\n",
       "      <th>GMR MAE</th>\n",
       "      <th>MLR MAE</th>\n",
       "      <th>RF MAE</th>\n",
       "      <th>GMR nMAE</th>\n",
       "      <th>MLR nMAE</th>\n",
       "      <th>RF nMAE</th>\n",
       "      <th>GMR cvMAE</th>\n",
       "      <th>MLR cvMAE</th>\n",
       "      <th>RF cvMAE</th>\n",
       "      <th>GMR Bias</th>\n",
       "      <th>MLR Bias</th>\n",
       "      <th>RF Bias</th>\n",
       "      <th>GMR mnBias</th>\n",
       "      <th>MLR mnBias</th>\n",
       "      <th>RF mnBias</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.780633</td>\n",
       "      <td>0.742420</td>\n",
       "      <td>0.727646</td>\n",
       "      <td>3.946637</td>\n",
       "      <td>4.438719</td>\n",
       "      <td>4.128082</td>\n",
       "      <td>0.160276</td>\n",
       "      <td>0.180260</td>\n",
       "      <td>0.167644</td>\n",
       "      <td>0.159767</td>\n",
       "      <td>0.180383</td>\n",
       "      <td>0.169409</td>\n",
       "      <td>0.108338</td>\n",
       "      <td>-0.026027</td>\n",
       "      <td>-0.350244</td>\n",
       "      <td>0.004400</td>\n",
       "      <td>-0.001057</td>\n",
       "      <td>-0.014224</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.704466</td>\n",
       "      <td>0.619090</td>\n",
       "      <td>0.783646</td>\n",
       "      <td>4.288418</td>\n",
       "      <td>5.039847</td>\n",
       "      <td>4.468371</td>\n",
       "      <td>0.172419</td>\n",
       "      <td>0.202630</td>\n",
       "      <td>0.179654</td>\n",
       "      <td>0.172443</td>\n",
       "      <td>0.204447</td>\n",
       "      <td>0.179248</td>\n",
       "      <td>-0.003412</td>\n",
       "      <td>-0.269339</td>\n",
       "      <td>0.066210</td>\n",
       "      <td>-0.000137</td>\n",
       "      <td>-0.010829</td>\n",
       "      <td>0.002662</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.660764</td>\n",
       "      <td>0.744299</td>\n",
       "      <td>0.768464</td>\n",
       "      <td>3.740364</td>\n",
       "      <td>4.246531</td>\n",
       "      <td>3.948843</td>\n",
       "      <td>0.149770</td>\n",
       "      <td>0.170037</td>\n",
       "      <td>0.158117</td>\n",
       "      <td>0.150502</td>\n",
       "      <td>0.169882</td>\n",
       "      <td>0.158703</td>\n",
       "      <td>-0.246906</td>\n",
       "      <td>0.040511</td>\n",
       "      <td>-0.242390</td>\n",
       "      <td>-0.009886</td>\n",
       "      <td>0.001622</td>\n",
       "      <td>-0.009706</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.804988</td>\n",
       "      <td>0.788287</td>\n",
       "      <td>0.810798</td>\n",
       "      <td>3.841505</td>\n",
       "      <td>4.378033</td>\n",
       "      <td>3.906195</td>\n",
       "      <td>0.159375</td>\n",
       "      <td>0.181635</td>\n",
       "      <td>0.162059</td>\n",
       "      <td>0.158409</td>\n",
       "      <td>0.179757</td>\n",
       "      <td>0.160363</td>\n",
       "      <td>0.292031</td>\n",
       "      <td>0.569016</td>\n",
       "      <td>0.332253</td>\n",
       "      <td>0.012116</td>\n",
       "      <td>0.023607</td>\n",
       "      <td>0.013784</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.649269</td>\n",
       "      <td>0.629398</td>\n",
       "      <td>0.791738</td>\n",
       "      <td>4.284543</td>\n",
       "      <td>5.037121</td>\n",
       "      <td>4.304947</td>\n",
       "      <td>0.167317</td>\n",
       "      <td>0.196706</td>\n",
       "      <td>0.168114</td>\n",
       "      <td>0.168088</td>\n",
       "      <td>0.198650</td>\n",
       "      <td>0.168476</td>\n",
       "      <td>-0.125808</td>\n",
       "      <td>-0.288094</td>\n",
       "      <td>-0.069939</td>\n",
       "      <td>-0.004913</td>\n",
       "      <td>-0.011250</td>\n",
       "      <td>-0.002731</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>0.649158</td>\n",
       "      <td>0.682369</td>\n",
       "      <td>0.726254</td>\n",
       "      <td>4.413161</td>\n",
       "      <td>5.031743</td>\n",
       "      <td>4.534429</td>\n",
       "      <td>0.173153</td>\n",
       "      <td>0.197423</td>\n",
       "      <td>0.177911</td>\n",
       "      <td>0.173760</td>\n",
       "      <td>0.197259</td>\n",
       "      <td>0.178930</td>\n",
       "      <td>-0.090909</td>\n",
       "      <td>0.029544</td>\n",
       "      <td>-0.137900</td>\n",
       "      <td>-0.003567</td>\n",
       "      <td>0.001159</td>\n",
       "      <td>-0.005411</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>0.775698</td>\n",
       "      <td>0.715695</td>\n",
       "      <td>0.824508</td>\n",
       "      <td>4.369453</td>\n",
       "      <td>4.834938</td>\n",
       "      <td>4.284810</td>\n",
       "      <td>0.174661</td>\n",
       "      <td>0.193268</td>\n",
       "      <td>0.171278</td>\n",
       "      <td>0.172947</td>\n",
       "      <td>0.192462</td>\n",
       "      <td>0.170202</td>\n",
       "      <td>0.419499</td>\n",
       "      <td>0.179155</td>\n",
       "      <td>0.264566</td>\n",
       "      <td>0.016769</td>\n",
       "      <td>0.007161</td>\n",
       "      <td>0.010576</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>0.733333</td>\n",
       "      <td>0.734988</td>\n",
       "      <td>0.753833</td>\n",
       "      <td>4.086275</td>\n",
       "      <td>4.648985</td>\n",
       "      <td>4.159055</td>\n",
       "      <td>0.165684</td>\n",
       "      <td>0.188500</td>\n",
       "      <td>0.168635</td>\n",
       "      <td>0.164585</td>\n",
       "      <td>0.187034</td>\n",
       "      <td>0.168651</td>\n",
       "      <td>0.231299</td>\n",
       "      <td>0.386331</td>\n",
       "      <td>-0.005549</td>\n",
       "      <td>0.009378</td>\n",
       "      <td>0.015664</td>\n",
       "      <td>-0.000225</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>0.738012</td>\n",
       "      <td>0.755563</td>\n",
       "      <td>0.779416</td>\n",
       "      <td>3.979185</td>\n",
       "      <td>4.519452</td>\n",
       "      <td>4.205381</td>\n",
       "      <td>0.155341</td>\n",
       "      <td>0.176432</td>\n",
       "      <td>0.164171</td>\n",
       "      <td>0.155611</td>\n",
       "      <td>0.176872</td>\n",
       "      <td>0.165200</td>\n",
       "      <td>-0.373508</td>\n",
       "      <td>-0.158623</td>\n",
       "      <td>-0.468369</td>\n",
       "      <td>-0.014581</td>\n",
       "      <td>-0.006192</td>\n",
       "      <td>-0.018284</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>0.707145</td>\n",
       "      <td>0.670499</td>\n",
       "      <td>0.809074</td>\n",
       "      <td>4.532837</td>\n",
       "      <td>5.070728</td>\n",
       "      <td>4.241189</td>\n",
       "      <td>0.175472</td>\n",
       "      <td>0.196294</td>\n",
       "      <td>0.164182</td>\n",
       "      <td>0.175753</td>\n",
       "      <td>0.197969</td>\n",
       "      <td>0.165527</td>\n",
       "      <td>-0.076466</td>\n",
       "      <td>-0.449154</td>\n",
       "      <td>-0.327348</td>\n",
       "      <td>-0.002960</td>\n",
       "      <td>-0.017387</td>\n",
       "      <td>-0.012672</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     GMR R2    MLR R2     RF R2   GMR MAE   MLR MAE    RF MAE  GMR nMAE  \\\n",
       "0  0.780633  0.742420  0.727646  3.946637  4.438719  4.128082  0.160276   \n",
       "1  0.704466  0.619090  0.783646  4.288418  5.039847  4.468371  0.172419   \n",
       "2  0.660764  0.744299  0.768464  3.740364  4.246531  3.948843  0.149770   \n",
       "3  0.804988  0.788287  0.810798  3.841505  4.378033  3.906195  0.159375   \n",
       "4  0.649269  0.629398  0.791738  4.284543  5.037121  4.304947  0.167317   \n",
       "5  0.649158  0.682369  0.726254  4.413161  5.031743  4.534429  0.173153   \n",
       "6  0.775698  0.715695  0.824508  4.369453  4.834938  4.284810  0.174661   \n",
       "7  0.733333  0.734988  0.753833  4.086275  4.648985  4.159055  0.165684   \n",
       "8  0.738012  0.755563  0.779416  3.979185  4.519452  4.205381  0.155341   \n",
       "9  0.707145  0.670499  0.809074  4.532837  5.070728  4.241189  0.175472   \n",
       "\n",
       "   MLR nMAE   RF nMAE  GMR cvMAE  MLR cvMAE  RF cvMAE  GMR Bias  MLR Bias  \\\n",
       "0  0.180260  0.167644   0.159767   0.180383  0.169409  0.108338 -0.026027   \n",
       "1  0.202630  0.179654   0.172443   0.204447  0.179248 -0.003412 -0.269339   \n",
       "2  0.170037  0.158117   0.150502   0.169882  0.158703 -0.246906  0.040511   \n",
       "3  0.181635  0.162059   0.158409   0.179757  0.160363  0.292031  0.569016   \n",
       "4  0.196706  0.168114   0.168088   0.198650  0.168476 -0.125808 -0.288094   \n",
       "5  0.197423  0.177911   0.173760   0.197259  0.178930 -0.090909  0.029544   \n",
       "6  0.193268  0.171278   0.172947   0.192462  0.170202  0.419499  0.179155   \n",
       "7  0.188500  0.168635   0.164585   0.187034  0.168651  0.231299  0.386331   \n",
       "8  0.176432  0.164171   0.155611   0.176872  0.165200 -0.373508 -0.158623   \n",
       "9  0.196294  0.164182   0.175753   0.197969  0.165527 -0.076466 -0.449154   \n",
       "\n",
       "    RF Bias  GMR mnBias  MLR mnBias  RF mnBias  \n",
       "0 -0.350244    0.004400   -0.001057  -0.014224  \n",
       "1  0.066210   -0.000137   -0.010829   0.002662  \n",
       "2 -0.242390   -0.009886    0.001622  -0.009706  \n",
       "3  0.332253    0.012116    0.023607   0.013784  \n",
       "4 -0.069939   -0.004913   -0.011250  -0.002731  \n",
       "5 -0.137900   -0.003567    0.001159  -0.005411  \n",
       "6  0.264566    0.016769    0.007161   0.010576  \n",
       "7 -0.005549    0.009378    0.015664  -0.000225  \n",
       "8 -0.468369   -0.014581   -0.006192  -0.018284  \n",
       "9 -0.327348   -0.002960   -0.017387  -0.012672  "
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "CVs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model Comparison"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Multiple Linear Regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.7332640382598803\n",
      "4.699392865443719\n",
      "-0.16678494718813303\n",
      "0.18658662871665996\n",
      "-0.0066023984011889005\n",
      "0.18603156018849507\n"
     ]
    }
   ],
   "source": [
    "# Create linear regression object\n",
    "regr_h = linear_model.LinearRegression()\n",
    "\n",
    "# Train the model using the training sets\n",
    "regr_h.fit(train_h[[\"purpleair_pm2.5\", \"avg_temp\", \"RH\"]], train_h[[\"embassy_pm2.5\"]])\n",
    "\n",
    "mlr_pred_h = regr_h.predict(test_h[[\"purpleair_pm2.5\", \"avg_temp\", \"RH\"]]).reshape(test_h.shape[0])\n",
    "\n",
    "print(r2_score(test_h[[\"embassy_pm2.5\"]], mlr_pred_h))\n",
    "print(mean_absolute_error(test_h[[\"embassy_pm2.5\"]], mlr_pred_h))\n",
    "print(n_bias(mlr_pred_h, test_h[[\"embassy_pm2.5\"]].to_numpy()))\n",
    "print(cvMAE(mlr_pred_h, test_h[[\"embassy_pm2.5\"]].to_numpy()))\n",
    "print(mnbias(mlr_pred_h, test_h[[\"embassy_pm2.5\"]].to_numpy()))\n",
    "print(nmae(mlr_pred_h, test_h[[\"embassy_pm2.5\"]].to_numpy()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Random Forest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.8001335405790052\n",
      "4.2394457978075515\n",
      "-0.26437880633373934\n",
      "0.16902470701241804\n",
      "-0.010465777863497191\n",
      "0.16782395911183973\n"
     ]
    }
   ],
   "source": [
    "RF_regr_h = RandomForestRegressor(max_features = 1, random_state = 0)\n",
    "RF_regr_h.fit(train_h[[\"purpleair_pm2.5\", \"avg_temp\", \"RH\"]], train_h[[\"embassy_pm2.5\"]])\n",
    "\n",
    "RF_pred_h = RF_regr_h.predict(test_h[[\"purpleair_pm2.5\", \"avg_temp\", \"RH\"]])\n",
    "\n",
    "print(r2_score(test_h[[\"embassy_pm2.5\"]], RF_pred_h))\n",
    "print(mean_absolute_error(test_h[[\"embassy_pm2.5\"]], RF_pred_h))\n",
    "print(n_bias(RF_pred_h, test_h[[\"embassy_pm2.5\"]].to_numpy()))\n",
    "print(cvMAE(RF_pred_h, test_h[[\"embassy_pm2.5\"]].to_numpy()))\n",
    "print(mnbias(RF_pred_h, test_h[[\"embassy_pm2.5\"]].to_numpy()))\n",
    "print(nmae(RF_pred_h, test_h[[\"embassy_pm2.5\"]].to_numpy()))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}