Part III - Predicting Which Houses Own EVs.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:9a1e4d4c886f4830bc1b3c2ce72eb584cc810cabc194f58c23045fbea1063f3d"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Predicting Which Houses Own EVs"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The goal of this section is the following: Given the power consumption, is it possible to predict which houses own EVs?\n",
      "\n",
      "This section takes two main approaches:\n",
      "- In part a, predictions are made using aggregate statistics\n",
      "- In part b, predictions are made using the Markov Model from Part II"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%matplotlib inline\n",
      "\n",
      "import pandas as pd\n",
      "import matplotlib as mpl\n",
      "import matplotlib.pyplot as pl\n",
      "import matplotlib.colors as colors\n",
      "import numpy as np\n",
      "import scipy\n",
      "import scipy.stats as stats\n",
      "import scipy.stats.mstats as mstats\n",
      "import scipy.interpolate as interpolate\n",
      "import sklearn.neighbors as neighbors\n",
      "import sklearn\n",
      "import random\n",
      "\n",
      "import util\n",
      "\n",
      "# DataFrame where index = intervals, columns = house ids\n",
      "usages = pd.read_csv('./data/EV_train.csv', index_col='House ID').transpose().dropna()\n",
      "when_charging = pd.read_csv('./data/EV_train_labels.csv', index_col='House ID').transpose()[0:len(usages.index)]\n",
      "\n",
      "# Only use houses that own EVs\n",
      "usages.index = range(len(usages.index))\n",
      "when_charging.index = range(len(when_charging.index))\n",
      "houses = usages.columns\n",
      "\n",
      "# Series where index = house ids, values = Bool\n",
      "houses_with_ev = houses[when_charging.sum() > 0]\n",
      "houses_without_ev = houses[when_charging.sum() == 0]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Part III a - Using Aggregate Statistics"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "As noted in Part I, houses with EVs tend to use more power on average, and also tend to have higher variances in power consumption.\n",
      "This section elaborates by looking for similar aggregate statistics, but with a focus on separating houses with EVs from those without EVs as reliably as possible. Much of this section is experimental and guided by intuition."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "First split into train/test data"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "random.seed(0)\n",
      "train_houses, validate_houses, test_houses = util.split_data(houses, ratios=[0.6, 0.2, 0.2])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Next, get some statistics for each house.\n",
      "The main statistics here are averages of each derivative, weighted by some polynomial k.\n",
      "The statistics are also split by normalization, to see if it helps or not."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "derivatives = util.derivatives(usages)\n",
      "\n",
      "normed_usages = util.normalize_usages(usages)\n",
      "normed_derivatives = util.derivatives(normed_usages)\n",
      "\n",
      "house_stats = pd.DataFrame(index = houses)\n",
      "house_stats['owns_ev'] = when_charging.sum() > 0\n",
      "house_stats['avg'] = usages.apply(np.average)\n",
      "house_stats['std'] = usages.apply(np.std)\n",
      "for k in range(1, 6):\n",
      "    for i,d in enumerate(derivatives):\n",
      "        key = \"d%i_k%i\" % (i, k)\n",
      "        values = d.apply(abs).apply(lambda x: np.power(x, k)).apply(np.average).apply(lambda x: np.power(x, 1.0/k))\n",
      "        house_stats[key] = values\n",
      "for k in range(1, 6):\n",
      "    for i,d in enumerate(normed_derivatives):\n",
      "        key = \"normed_d%i_k%i\" % (i, k)\n",
      "        values = d.apply(abs).apply(lambda x: np.power(x, k)).apply(np.average).apply(lambda x: np.power(x, 1.0/k))\n",
      "        house_stats[key] = values"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "for (feature, values) in house_stats.iteritems():\n",
      "    if feature != 'owns_ev':\n",
      "        values = (values - np.mean(values)) / np.std(values)\n",
      "        house_stats[feature] = values"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Out of all the features generated above, see which statistic is the most accurate using 1D KNN."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from sklearn.neighbors import KNeighborsClassifier\n",
      "\n",
      "accs_1D = {}\n",
      "stats = house_stats.columns[1:] # take out the 'owns_ev' column\n",
      "for stat in stats:\n",
      "  \n",
      "    X = house_stats.T[train_houses].T[stat].values\n",
      "    y = house_stats.T[train_houses].T['owns_ev'].values.astype('int')\n",
      "\n",
      "    clf = KNeighborsClassifier(n_neighbors=20)\n",
      "    clf.fit(np.array([X]).T, y)\n",
      "    \n",
      "    X = house_stats.T[validate_houses].T[stat].values\n",
      "    y = house_stats.T[validate_houses].T['owns_ev'].values.astype('int')\n",
      "\n",
      "    pred = clf.predict(np.array([X]).T)\n",
      "\n",
      "    accs_1D[stat] = util.acc_info(pred, y)\n",
      "\n",
      "accs_1D = [(k,v) for k,v in accs_1D.iteritems()]\n",
      "accs_1D.sort(key = lambda (k,v): -v['acc'])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "for (k,v) in accs_1D[:10]:\n",
      "    print \"key: %s, acc: %.3f (%.3f), pos acc: %.3f (%.3f), neg acc: %.3f (%.3f)\" % (k, v['acc'], v['acc_unc'], v['pos_acc'], v['pos_acc_unc'], v['neg_acc'], v['neg_acc_unc'])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "key: d1_k5, acc: 0.840 (0.021), pos acc: 0.724 (0.045), neg acc: 0.891 (0.021)\n",
        "key: d1_k4, acc: 0.827 (0.021), pos acc: 0.684 (0.047), neg acc: 0.891 (0.021)\n",
        "key: d3_k5, acc: 0.818 (0.022), pos acc: 0.704 (0.046), neg acc: 0.868 (0.023)\n",
        "key: d2_k5, acc: 0.805 (0.022), pos acc: 0.684 (0.047), neg acc: 0.859 (0.023)\n",
        "key: d1_k3, acc: 0.802 (0.022), pos acc: 0.643 (0.048), neg acc: 0.873 (0.022)\n",
        "key: d3_k4, acc: 0.789 (0.023), pos acc: 0.582 (0.050), neg acc: 0.882 (0.022)\n",
        "key: d1_k2, acc: 0.767 (0.024), pos acc: 0.582 (0.050), neg acc: 0.850 (0.024)\n",
        "key: d0_k5, acc: 0.764 (0.024), pos acc: 0.592 (0.050), neg acc: 0.841 (0.025)\n",
        "key: d2_k3, acc: 0.761 (0.024), pos acc: 0.643 (0.048), neg acc: 0.814 (0.026)\n",
        "key: d2_k4, acc: 0.755 (0.024), pos acc: 0.571 (0.050), neg acc: 0.836 (0.025)\n"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "X = house_stats.T[train_houses].T['d1_k5'].values\n",
      "y = house_stats.T[train_houses].T['owns_ev'].values.astype('int')\n",
      "clf = KNeighborsClassifier(n_neighbors=20)\n",
      "clf.fit(np.array([X]).T, y)\n",
      "X = house_stats.T[test_houses].T['d1_k5'].values\n",
      "y = house_stats.T[test_houses].T['owns_ev'].values.astype('int')\n",
      "pred = clf.predict(np.array([X]).T)\n",
      "acc = util.acc_info(pred, y)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "acc: 0.8365 (0.0207), pos acc: 0.6429 (0.0523), neg acc: 0.9060 (0.0191)\n"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "bins = np.linspace(-3, 5, 41)\n",
      "xs = (bins[1:] + bins[:-1])/2\n",
      "h0 = np.histogram(house_stats[house_stats['owns_ev']]['d1_k5'], bins=bins)\n",
      "h1 = np.histogram(house_stats[~house_stats['owns_ev']]['d1_k5'], bins=bins)\n",
      "pl.plot(xs, h1[0], marker='o', color='blue', label='No EV')\n",
      "pl.plot(xs, h0[0], marker='o', color='red', label='Owns EV')\n",
      "pl.xlabel('Average normalized first derivative, weighted k=5')\n",
      "pl.ylabel('Counts')\n",
      "pl.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEPCAYAAABCyrPIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8XFX9//FXmrRJ952udEtblUUoILKmKdK0kLIoaxG+\n8GVXbAooW9N+E4VqEVFo/IqCCKhf9iI/SKBNEZoUUBC6sFmBNGHvhtQukNKk5/fHudOZSe7M3Mls\nd5L38/GYR+7c9eTOzP3cs9xzQEREREREREREREREREREREREREREJKMKgJeANcBbwM+c+YOA5cDb\nQC0wIGSbG4B3gHVASdpSKiIiKdfL+ZsH/B04Bvg5cK0z/zpgkTO9HzZ4dAfGAe8C3dKVUBERSY9e\nwD+A/bF3/8Oc+cOd92BzB9eFbLMUOCJdCRQR6epSfQfeDXvXvxF4DngTGww2Oss3EgwOI4EPQ7b9\nEBiV4vSJiIgjL8X73wMcDPQHlgHT2iw3ziuSaMtERCSJUh0QAv4D1ACHYnMFw4ENwAhgk7POR8C+\nIduMduaFKSwsNA0NDSlNrIhIJ9QATIy2QiqLjIYQbEHUE5gOrAaeAM535p8PPO5MPwGcDfQAxgOT\ngJfb7rShoQFjjO9fFRUVGU9DutJZXV1HYeE8ghk+Q2HhPKqr63yTxmw5l0qn0pmqF1AY66KdyoAw\nAngWW4fwEvAk8Fdsq6Lp2GanxxFsZfQW8LDz92ng+6jIKCssXlxLQ8PCsHkNDQupqlqeoRSJSEek\nssjodeAQl/n/Bo6PsM1PnZdkkV273L9Gzc25aU6JiCRC7fxTpLi4ONNJ8CQZ6czPb3GdX1DQmvC+\noWudy3RQOpMrW9LpRU6mE9ABxikPE5+oqaln7txlYcVGhYXzuP32mZSWFmUwZSISkJOTAzGu+QoI\nkhQ1NfVcffVy3n47lylTWrnxxukKBiI+ooAgaXX55XD33VBVZadFxD+8BATVIUjSrF8Phx4KH7V7\nekREsoECgiTN+vVw7LEKCCLZSgFBkqKlBT74AI4+WgFBJFspIEhSfPgh7LMPTJiggCCSrRQQJCnW\nr7fBYNQo+PjjTKdGRDpCAUGSIhAQBg2CL76wLxHJLgoIkhSBgJCTAyNHqthIJBspIEhSNDbC+PF2\netQoBQSRbKSAIEkRyCGAAoJItlJAkKRQQBDJfgoIkrBt22DnThjmjI6tOgSR7KSAIAlrbAxWKINy\nCCLZSgFBEhZaXAQKCCLZSgFBEqaAINI5KCBIwtoGhJEjYcMG2LMnc2kSkfgpIEjCAnUIAfn50K8f\nbN6cuTSJSPwUECRh69cHH0oLUJ9GItlHAUESsmcPNDW5BwTVI4hkFwUEScjHH8PAgdCrV/h8BQSR\n7KOAIAlpW6EcoIAgkn0UECQhCgginYcCgiREAUGk81BAkIRECgjqz0gk+yggSELaPoMQoByCSPZJ\nZUDYF3gOeBN4Ayhz5lcCHwKrndcJIdvcALwDrANKUpg2SZJIOYTBg+HzzzWUpkg2yUvhvncDVwFr\ngD7Aq8BywAC/dF6h9gPOcv6OAp4BJgPqAMGnPv8cPvsMRoxovyx0KM2JE9OfNhGJXypzCBuwwQBg\nB/BP7IUeIMdl/VOAB7CBpAl4Fzg8hemTBDU2wrhx0C3Ct0jFRiLZJV11COOAKcDfnfdzgLXA3cAA\nZ95IbFFSwIcEA4j4UKTiogAFBJHsko6A0Ad4FJiLzSncAYwHDgY+AW6Nsq1JeeqkwxQQRDqXVNYh\nAHQHlgB/Bh535m0KWf574Eln+iNsRXTAaGdeO5WVlXuni4uLKS4uTkpiJT5eAsIHH6QvPSIStGLF\nClasWBHXNm5l+cmSA9wHfIqtXA4Ygc0Z4Mz/BnAOtjL5fmy9QaBSeSLtcwnGGGUc/OCkk+Cii+DU\nU92XP/QQLFkCDz+c3nSJSHs5dozbqNf8VOYQjgbOBV7DNi8FmAfMxhYXGaARuMxZ9hbwsPO3Bfg+\nKjLytUjPIASoyEgku6Qyh5AqyiH4gDHQp48dGa1vX/d1Ghth2jTbPbaIZJaXHIKeVJYO2bgReveO\nHAzAPofwyScaSlMkWyggSIe4jZLWVn6+DRhbtqQnTSKSGAUE6ZBYLYwCVI8gkj0UEKRDFBBEOh8F\nBOkQBQSRzkcBQTpEAUGk81FAkA6J9QxCgAKCSPZQQJC4NTfD5s0wenTsdRUQRLKHAoLE7b33YN99\nITc39rqjRsHHH6c+TSKSOAUEiZvX+gNQDkEkmyggSNy8PJQWMHgw7NypoTRFsoECgsQtnhxCYChN\nFRuJ+J8CgsQtnoAAKjYSyRYKCBK3eAPCyJEKCCLZQAFB4mKMcgginZUCgsTl00+he3cYMMD7NgoI\nItlBAUHiEm/uABQQRLKFAoLERQFBpPNSQJC4KCCIdF4KCBKXeB5KC9BQmiLZQQFB4tKRHEJBgYbS\nFMkGCggSl44EBFAndyLZQAFBPNu92xb9jBkT/7aqRxDxPwUE8aSmpp7jjptPTk4ls2bNp6amPq7t\nFRBE/C8v0wkQ/6upqWfu3GU0NCwEoLYWGhrKASgtLfK0DwUEEf9TDkFiWry4dm8wCGhoWEhV1XLP\n+1B/RiL+p4AgMe3a5Z6RbG72MGSaQzkEEf9TQJCY8vNbXOcXFLR63ocCgoj/KSBITGVlJRQWlofN\nKyycx5w50z3vQwFBxP9yUrjvfYE/AvsABrgTWAwMAh4CxgJNwJnAVmebG4ALgVagDKh12a8xxqQw\n2eKmpqaeiy5azqBBuYwZ08qcOdM9VyiD7Ta7Z0/47DP7V0TSKycnB2Jc81MZEIY7rzVAH+BV4FTg\nv4EtwM+B64CBwPXAfsD9wDeAUcAzwGSgbYcHCggZMn06XHMNlJR0bPvx4+GZZ6CwMLnpEpHYvASE\nVBYZbcAGA4AdwD+xF/qTgfuc+fdhgwTAKcADwG5szuFd4PAUpk/i9NlnMHBgx7dXsZGIv6WrDmEc\nMAV4CRgGbHTmb3TeA4wEPgzZ5kNsABGfUEAQ6dzS8WBaH2AJMBfY3maZcV6RuC6rrKzcO11cXExx\ncXFCCRRvPvssvpHS2lJ/RiLps2LFClasWBHXNqmsQwDoDlQDTwO3OfPWAcXYIqURwHPAV7H1CACL\nnL9LgQpsriKU6hAyYM8e6NEDmpshr4O3EbfeanMIv/xlctMmIrFlug4hB7gbeItgMAB4AjjfmT4f\neDxk/tlAD2A8MAl4OYXpkzhs3w69enU8GICKjET8LpVFRkcD5wKvAaudeTdgcwAPAxcRbHYKNnA8\n7PxtAb5P9OIkSaNE6w9AAUHE71IZEJ4ncg7k+Ajzf+q8xGcSrT8A9Wck4nd6Ulk82bo18RxCYChN\nVQGJ+JO6vxZPklFk9Oyz9bS21nL00Xn07dtCWVlJXE87i0hqKSCIJ4kGhMCYCi0tC/nb3+y8eMdU\nEJHUUpGReJJoHUIyxlQQkdRSQBBPEq1DSMaYCiKSWgoI4kmiRUbJGFNBRFJLAUE8STQgJGNMBRFJ\nLVUqiyeJ1iEEKo6rqhbwwgu57LdfK//zPzNVoSziIwoI4kkynkMoLS2itLSIa6+FPn2gtDQ5aROR\n5FCRkXiSjOcQAoqLIc5OGEUkDRQQxJNkBoRjjoGXX7Y9p4qIfyggSEzGJKcvo4B+/WC//WxQEBH/\nUECQmL74AnJzoaAgeftUsZGI/yggSEzJLC4KUEAQ8R8FBIkpmcVFAapHEPEfBQSJKRU5BNUjiPiP\nAoLElIxnENyo2EjEX+INCIOAr6ciIeJfqcghgAKCiN94CQh1QD9sMHgV+D3wq1QmSvwlFXUIoHoE\nEb/xEhD6A9uA7wB/BA4n8pjI0gmlKoegegQRf/ESEHKBEcCZQI0zT6PidiGpqkMAFRuJ+ImXgPAT\nYBnQALwMFALvpDJR4i+pyiGAAoKIn3jp7fQTwiuSG1AdQpeSqjoEsPUIZ55p6xGS+SS0iMTPSw6h\nymXe4mQnRPwrlTkE1SOI+Ee0HMKRwFHAUOBqIMeZ3xdbryBdRCrrECBYbFSksXJEMipaDqEHwYt/\nX6CP89oGnJ76pIlfpDKHAKpHEPGLnNirMA5oSm0y4mKMUSOndOrdGzZsgL59U7P/bdtg5EjYskX1\nCCKpkpOTAzGu+V4qlfOBu7CBIbC+AY5LIG2SJb780r769EndMULrEVRsJJI5XiqVHwFWAfOBa0Je\nXvwB2Ai8HjKvEvgQWO28TghZdgO2Ses6oMTjMSSFtm61LYxyvOQlE6BiI5HM8xIQdgN3AC8Brziv\nVz3u/x5gZpt5BvglMMV5Pe3M3w84y/k7E/iNx/RJCqW6/iBAAUEk87xccJ8ErsA+rTwo5OXFSuAz\nl/lu95unAA9gA1AT8C62mwzJoFQ+gxBK/RqJZJ6XgHAB8CPgRWzOIPBKxBxgLXA3ELjcjMQWJQV8\nCIxK8DiSoFQ3OQ3Q8wgimeelUnlcko95B7Y7DIAbgVuBiyKs69qcqLKycu90cXExxcXFyUudhElX\nkRHoeQSRZFqxYgUr4iyH9VJVeD7uF+Y/ejzGOGyx04Exll3vzFvk/F0KVGDrLkKp2Wka/eY38Prr\ncMcdqT/WU0/BL34Bzz6b+mOJdDXJanb6DYIBoSe2uekqvAeEtkZg+0cC+DbBFkhPAPdjK5xHAZOw\nnelJBqWrDgFgx4566upqKSrKo2fPFsrKSigtVXZBJF28BIQftHk/AHjI4/4fAKYCQ4APsHf8xcDB\n2CDTCFzmrPsW8LDztwX4PupmO+O2boWhQ1N/nJqaeubNW8aePQtZudLOa2goB1BQEEmTjrQu7wG8\nAUxOclq8UpFRGl18MXzzm3DJJak9zowZ86mtvcll/gKWLr0xtQcX6QKSVWT0ZMh0N+xzAg93PFmS\nTdJVqbxrl/tXsblZ/SiKpIuXgHCr89dgi3Lexxb/SBeQrjqE/PwW1/kFBa2pP7iIAN6eQ1iB7Uqi\nHzAQ2JXKBIm/pOs5hLKyEgoLy8PmFRbOY86c6ak/uIgA3nIIZwK3AHXO+19j+zJ6JFWJEv9IV5FR\noOK4qmoB//hHLmPGtHLTTTNVoSySRl4qlV8Djgc2Oe+HAn8lfFjNdFKlchoNGACNjel7OA1gwQLb\nmd5PfhJ7XRHxxkulspcioxxgc8j7T2PtVDqH1lbYvt12K5FOhxwCq1al95gi4q3IaCmwDPvQWA62\nR9Kno24hncK2bXZQnNw0N/SZMgVWr07vMUUkekCYBAzD1hecBhztzH8RGxykk0tnP0ahxo6FL76w\no7QNH57+44t0VdGKjG7Djp8MsAS42nk9DvwqxekSH8hUQMjJscVGyiWIpFe0gDAMW6Hc1mvA+NQk\nR/wknf0YtaViI5H0ixYQol0KNBR6F5CuZxDcqGJZJP2iBYRXgEtd5l9C4gPkSBbIVJERKCCIZEK0\nSuUrgb8A3yUYAA4F8rHdVksnl8mAMGkSbN6c2TSIdDXRAsIG4ChgGnAAti+jakDDl3QRmaxD6NYN\nDjoI1qyBadMykwaRribWcwgGGwAUBLqgrVth330zd/xAsZECgkh6eHlSWbqoTBfXqKWRSHopIEhE\nmQ4IqlgWSS8FBIkok3UIAPvtB01NsHNn5tIg0pV46ctIuqhMPocA0L27DQqvvQZHHultm/qaGmoX\nLyZv1y5a8vMpKSujqLQ0tQkV6SQUECSiTBcZQbDYyEtAqK+pYdncuSxsaNg7r9yZVlAQiU1FRuLK\nGJtDyGSREcRXj1C7eHFYMABY2NDA8qqqFKRMpPNRDkFc7dgBPXrYVyZNmQJ33hl8H7FIaNMm8toE\ng4Dc5uY0pVYkuykgiKtM1x8EfP3rsG4d7NoFLz3jUiS0Zg2MHk1RQwMtvXq57qO1QF1viXihIiNx\n5Yf6A4CePaGwEN58M0KR0KZNLG9pgY8+ouSuuygvLAxbPq+wkOlz5qQzySJZSzkEceWXgADBB9Ty\ndu1yXZ47cCD07r234njBL35B7vPP0/qtbzFzzhxVKIt4pIAgrjL9DEKoQMXywPx81+WhRUJFpaUU\nnXgi9OkDDz0E/funK5kiWU9FRuLKL3UIYHMIq1ZBSVkZ5aNGhS1zLRLKyYGJE+Hdd9OYSpHspxyC\nuPJTkdHBB8Prr8PRM0th9mwW/OlP5H71q7QWFEQuEpo0Cd55Bw49NP0JFslSqQ4IfwBKgU3Agc68\nQcBDwFigCTgT2OosuwG4EGgFyoDaFKdPIvBTQOjfH0aMgH/9C4oGDKDoggtg0aLoGymHIBK3VBcZ\n3QPMbDPvemA5MBn4q/MeYD/gLOfvTOA3aUifROCnOgQI6fm0sRHGexjSe9IkBQSROKX6grsS+KzN\nvJOB+5zp+4BTnelTgAeA3dicw7vA4SlOn0TgpzoECHli2WtAmDjRFhmJiGeZuAMfBmx0pjc67wFG\nAh+GrPchEF6DKGnjpyIjCAkITU0wblzsDZRDEIlbpiuVjfOKtrydysrKvdPFxcUUFxcnNVHiv4Aw\nZQqsXdWKaf6QnLFjY28wYoTtf2PbNujXL/UJFPGZFStWsGLFiri2yURA2AgMx47ZPAJb4QzwERA6\nYONoZ147oQFBUsNvdQhDh8LkXh/S2msoeRGeRwgT2vT0kENSn0ARn2l7s/zjH/845jaZKDJ6Ajjf\nmT4feDxk/tlAD2A8MAl4Oe2pE8B/dQgAx01oYusAD/UHAapHEIlLqnMIDwBTgSHAB8D/AIuAh4GL\nCDY7BXjLmf8W0AJ8n+jFSZJCfisyAjh8aCMfbBvHEK8bqB5BJC6pDgizI8w/PsL8nzovyaDmZtiz\nx3Ys5ydfK2jk1S/GM8XrBhMnwvPPpzJJIp2K2vlLO4H6g5ycTKck3OiWRl7aNB7jNd+oHIJIXBQQ\npB0/1h8A9NrUxHs54/jkE48b6GllkbgoIEg7fqw/AMhpbKTfQeM9D6nJiBGwfbtteioiMSkgSDu+\nDAi7dsGmTex75GjvAaFbNzu6ToShNUUkXKYfTBMf2rrVX88gAPD++zBqFC28yF2/ruXZZ/PIz2+h\nrKyE0tKiyNsFmp5O8VwVLdJlKSBIO77MITQ1saXvAB58cBmbNy+krs7ObmgoB4gcFFSxLOKZioyk\nHV8GhMZGXvl0D++/vzBsdkPDQqqqlkfeTg+niXimgCDt+DUgfNR9kOui5ubcyNsphyDimQKCtOPL\nOoSmJrb06eu6qKCgNfJ2yiGIeKaAIO34NYdw7H+VUFhYHja7sHAec+ZMj7zdyJG22en27SlOoEj2\nU6WytOPXgHDUd7/N7fu9S1XVAp5/PpcDDmhlwYKZ0VsZBZqevvuuWhqJxKCAIO34LiDs3Gnv8ocP\np7R0JKWlRcybZ7vWKC31sH2gHkEBQSQqFRlJO76rQ3jvPRgzxt7tO048EZ56yuP26sJCxBMFBGnH\ndzkEl3GUjzjCxomPXIdQamPSJFUsi3iggCBhWlrg88+hr3uDnsxwCQh5eTB9Oixd6mF75RBEPFFA\nkDBbt0L//mGlM5nX1ATjxrWb7bnYSE1PRTzx089efMB39QfgmkMAmDkT/vpX2L07xvajRsF//gM7\ndqQmfSKdhAKChPFd/QFEDAjDhtnqgRdeiLF9t24wYYKKjURiUECQMNkUECCOYiN1YSESkwKChPFd\nQNi61ZYJDR7sulj1CCLJo4AgYXxXh9DUZHMHEQZ4Puww2LjRNkGNSjkEkZgUECSM73IIUYqLAHJz\nYcYMePrpGPtR01ORmBQQJIzvAkKEJqehPBUb6eE0kZgUECSM7wJCjBwC2BzCihV22OWIRo2y5WE7\ndyY1eSKdiQKChPFdHYKHgDB4MBxwANTXR1lJTU9FYlJAkDDZmEMAj8VGqkcQiUoBQcL4KiAY46kO\nAVSPIJIMGg9BwvgqIGzZAj162M6VYjj4YNs7xbvv2oyAq4kT4R//iLmv+poaahcvJm/XLlry8ykp\nK6PI08ALItktkwGhCdgGtAK7gcOBQcBDwFhn+ZnA1swkr2vyVR2Cx+IisFUEJ5xgm5/OmRNhpYkT\n4f77o+6nvqaGZXPnsrChYe+8cmdaQUE6u0wWGRmgGJiCDQYA1wPLgcnAX533kiZ79ti77FQFhPqa\nGubPmEFlcTHzZ8ygvqYm+gaNjZ6KiwJiFht5eDitdvHisGAAsLChgeVVVZ7TIZKtMl1k1Pbx05OB\nqc70fcAKFBTSZvt26NXLjjWQbB268w48pezR9Olw4YV2PIdevVxWGD0a/v1v2/S0d2/XfeRFaLua\n29zsOR0i2SrTOYRngFeAS5x5w4CNzvRG572kSSrrDzp05x1HkRHYnM0hh9hnElwFmp62SUeolj17\nXOe3quts6QIyGRCOxhYXnQBcARzbZrlxXpImqaw/6NCdd5wBATwUG0VremoMJV98QfmgQWGz540c\nyfT33oOyMlBOQTqxTBYZfeL83Qz8BVuPsBEYDmwARgCb3DasrKzcO11cXExxcXEKk9l1pDKH0JKf\n7zq/taAg8kYem5yG6t27nrvvruX11/MoKGihrKyE0tKi4ArRmp4uWULR55/D3Xez4Le/Jbe5mdaC\nAmbOmUPRUUfBZZfBN75B/eWXU/vEE2qFJL62YsUKVkTMLvtLLyAwam9v4AWgBPg5cJ0z/3pgkcu2\nRlJjyRJjTj01Nfuuq64283r1MsY+XWAMmBsKC01ddbX7Bq2txuTnG7Nzp+djVFfXmcLCeaGHMIWF\n80x1dV1wpTvuMObii9tvvHWrMSNHGrNyZeQD7Nlj6q680szr1i3s/5gX7f8Q8Qk8lLhkqshoGLAS\nWAO8BFQDtdgAMB14GzgO94AgKZLKHELRYYcxA1hw/PFUHnQQC3r3ZuZtt0W+s/7kE1t+5Vo77G7x\n4loaGhaGzWtoWEhV1fLgjEg5hPJyKC2FY46JfICcHGrfeouFbeoZ1ApJOotMFRk1Age7zP83cHya\n0yKOlD6DcP/9FJ1xBkX33mvvq/fbD/r1i7x+nE1OAXbtcv86NzfnBt+41SG89BIsWQJvvhnzGGqF\nJJ2Zuq6QvVKWQzAG7rkHLrjAvs/JgUsvhd/9LvI2cTY5BcjPb3GdX1DQGnyz777w6ae2bSrY0dgu\nvRRuvRXaVCa76VBdiEiWUECQvVIWENassQ85FIVU7p5/PtTU2Iuzmw60MCorK6GwsDxsXmHhPObM\nmR6c0a2b3W+g6entt8OwYTB7tqdjlJSVUV5YGDZvXmEh0yM+Hi2SPTL9YJr4RE1NPY89VsuAAXk8\n+aRL65xE3HuvDQDdQu4/Bg2Ck06C++6Dq69uv01jIxx5ZFyHCaS3qmoBn32Wy6pVrdxyy8z2/0eg\nHqFvX1i0yBYZRRiis61AnceCqipyN2ygdf16Zt5+u1oZiWRIpivrOx1PrXM6atcuY4YONaahof2y\nlSuN+cpXjNmzp/2yadOMqa1N6NDFxcY8+mj4vLrqalM+dqypGD/elA8ZYurOO6/jB9i925ghQ4xp\nakoonSLpgI9bGYmPeGqd01E1NbYCecKE9suOPtoOilxX135ZB4qM2rrgAps5CQh0n3HTe+9R2djI\nTVu2sOyFF2L3qRRJXp7N5fy//5dQOkX8QgFBvLXO6ah77w1WJreVk2Mf9mpbudzSAh9/DGPGJHTo\n006DlSthwwb73rX7jPXrE2syeuqp8PjjCaRSxD8UEMRb65yO2LjR3v2ffnrkdc47z/ZZvXlzcN4H\nH9iK3h49Ejp8nz7w7W/D//2ffZ+SJqPTp8Orr0auHBfJIgoIQllZCfn5MVrndMT999s76D59Iq8z\ncKBd5777gvM60GVFJIFiI2NS1GS0Z0/41rds0ZhIllNAEPbfv4gePWZQUrKAqVMrmTFjAbff7tI6\nJx5tnz2I5tJL4c477TaQlPqDgGOPtb1dr1qVwiajKjaSTkLNToUlS+Css4q4664kNTMF92cPIjny\nSMjPh+eeg+OOS2pA6NbNtni9916oqgppMhracV2iTUZLS+0wbREHYhDJDt4aX/uL04JKkuWII+An\nP4GSkiTudO5cWxwU0jNtVL/+ta0BfughOPdcOP54b7kLDxob4RvfgI8+snEnJY47Dq68Ek4+OUUH\nEElMjn3WJuo1X0VGXdz779uufaZNS+JOv/wSHngA/uu/vG9z7rmwbBls2tShbiuiGT8eDjwQqquT\ntsv2VGwknYACQhe3ZAmccgp0757EnUZ79iCSAQPgO9+x9Q5JLDIKaPtMQtKdcgo8+aRtMiuSpRQQ\nurhHHoEzzkjyTqM9exBF/f77M7+yksqPP2b+xRd3/IExF6edBs8/H3wmIenGjrUd5734YooOIJJ6\nqlTuwj74AP71L1v8naj6mhpqFy8mb/t2Wl5+mZLzziOeKur6mhqW3XEHCwPPBCxfTvn69QBJ6Sco\n9JmEH/4w4d25CxQbealIF5GkyHSXIJ3Gr35lzAUXJL6fuupqM6+wMKFRxMpLSsK2D7zmz5iReAID\n6awz5oAD3LtOSoq1a40ZPz6FBxDpONSXkUTz6KPJKS5y7RIizlHE0jHwzDHH2Jahq1YlbZfhDjzQ\n/n399RQdQCS1FBC6qI8+grfesq07E5WMi3k6Bp4JfSYhJXJy1NpIspoCQhe1ZIntqDPB7oKA5FzM\n0zXwzOjR9fz2t/MpKqpkxoz51NTUJ3X/CgiSzVSp3EU9+ihce21y9lVSVkb52rUs3Lhx77x5hYXM\njONiHjbwTDKfIg5RU1PPT3+6jJaWhaxcaec1NNg+nJI2GNBRR9na+vfesy2PRLKInlTugj75xD4m\nsGFDkp7c/fJL6seMYfmYMeT26kVrQQHTk3wxT4YZM+ZTW3uTy/wFLF16Y/IOdOGFcPDBUFaWvH2K\nJMjLk8rKIXRBS5bArFlJ7MbhzjspmjKFoqefTtIOUyPSuA///Gcub78Nkyfb9zU19SxeXMuuXXnk\n53dgONFTT4XbblNAkKyjgNAFPfJIEtvi79gBCxfaMQ18LtK4Dzk5rUydCqNHw8EH1/PMM8toagqO\nIBdvsdKtf1/L5udW8E6/fdnZozvH/+AyflR5XeL/gIi0k+nmvFntk0+M6d/fmC++SNIOf/xjY845\nJ0k7Sy3SW7wHAAASDUlEQVT3saNvMNXVdWb3bjuE88iR5W6PQ5gZM+Z7OsYtFYvMd/MGhG383bwB\n5paKRSn+7+JXXV1nSkrKzdSpFaakpDw5Y2iLb+HhOQTlELqYxx6zvTUnpTXn5s2weDG8/HISdpZ6\ngTv8qqoFNDfnUlDQypw5wXEfpk+HSZPy+Pjj9tt6HU70mV//jqUtW8Pm/bllKyf8+k5f5RJqauqZ\nO3dZ2FjaSa9gl6yjgNDFPPKI7Zk6KW66CWbPjq8TuwwrLS2KesGLVKy0bl3r3nqGaHUM3Xftcd2+\nx+7kd3qXSF3H4sW1YcEAoKFhIVVVCxQQujAFhC5k40ZYvRpmzEjCzhob4c9/tk+3dSJlZSU0NJSH\nXSzHj5/HscfO5Kij4LDD6lm3bhnvvRdc/vbb5TzxBLz2WhGtO9wf7Rm/7RNWXngPh992Dvn98vlF\n5c088+vfUdCyh+a8bu3qGfb2DbVrFy35+ZSUlYW12kr0Dj9SBbvXnJCIX2S6KC7rBMqKJ02qMMOH\nJ6ms+NxzjamoSHw/PlRdXWdmzJhvpk6tMDNmzN97vrZsMWbcOPc6hhEj5punnjLm5vnt6xDOyRtg\nKqZdaF4eMtNs6DbcLBx7vDkvt3/Eeoa66mozd/jIsOVzh48M6xuqpKTjdR3NzZH/j2nTvNWVJJPq\nMtIDD3UIfjQTWAe8A7gVumb6vGYV94rUeYn96NauNWbYMGO2bUteQrPE1KkVrhfSqVMr9q5zS8Ui\nM3PwBHNK/7Fm5uAJYRXK7zz+hrmwW1/XjvzO7TnYtNxWZcpGj3FdfvbEw83Pf25jce/esdPh5q23\njJkyxZjDD68z48aFfy8GDLjBDB1aZ159Nbh+XXW1KS8pMRVTp5rykpJ2HRbGuph7WZ7o99NLQEk0\nncmQ6cBHFgaEXOBdYBzQHVgDfK3NOmk9iR313HPPZToJxpjYd5IdSueJJxpz++3JTWgUfjmXxkQ/\nn17TeUr/sa4X/IvpY/435wpzMX1cl5/Wbai5qqzF3H23MUccUW76UG0Oo8RMZao5jBLTh2ozefJ8\ns3WrPU7bC9D3vldnBg825qqrnjN79rjnhB55xJghQ4z5059i92Ib62IebfmePcZ88IExhxwS+3xG\nC0peAkoi6QwV6YLu5XNP9Bhel0c7V2RhQDgSWBry/nrnFcoYY+/CZgweb07pP9bMGDy+XbO+WMtT\nvY8KpzglU+ncudOYBx80ZtCgCtcLR+BOMpDOWHeCe5cfdJApLygwdY891i6dqVLho6KpaE1XvaZz\nxuDxrhf8mYMnmJYWY47vP8F1+VXdehgzfLgxV15p/vfiK9oXTeUOMEfsv8gMHGjMSSfVmfFDzwj7\n3AfknmHuuKMu5mf++uvGFBYac9qYyF2S79hhL+Zu362RI+eba64xZvJk9+UDB843Awcas88+xgwY\n4P79nDy5wlx3XUXM4rOSEvdjHHnkfLN6tTGrVxtz5JHR14m0PLT4rbq6zkwaHn4+Jw0/I+xzj/Y7\njZROr8fwsjzWuSILA8LpwF0h788F2vahHLOtt5e24KneR0VFRcbSeXbOADO4xyIzfbox+48505xJ\n+F3emRSaoilnGWPsxTbWnWAyxjtIhJ8CgjGR6xi8ptPtMzsn5DMtmtL+Mzsj8JmtW2fMggWmvGdP\n94t1cbH5oHG3mTA08uce8TMfO9bUVVQYc/PNpvmM75prcgpcj3FJ3j7mvPyHzMH5J7oe4ysjZ5uf\n/cyYsYNmuy4/YOxss2GDifi/nkmhGd/3dDO0+1xzap+vuabh3H0PNDfP22pG9jnbdfsh+bPNQQcZ\nc9BBxgzJd09HYJ1Iy/vnzjZf/7ox06cbM2agezqP2P8sM3+++2/9nNwB5oLTFpnycmOG94lwjG6z\nzQEHGHPcccbsO8D9GF8ddZa59VZjvjIyQhq+dobZ9q+Pzfe+ur/rubr0kMONMdkZEE7DQ0CIdId1\nOb3MS/vMMpfh/mMJLPeyTqL7uLjXpJQfI9ryK/J6GTNrlvlRv/6uy6/p19+YWbNMxaRJpnzIEPeL\ny5AhxsyaFXl5EgevicZvASGSeNIZrZ4heCc4w7kTnGEmDT89rHigoqjI9TOpyMszJi/PXEu3iN+L\nqJ/5oEHGXHWVMX/4g/l2v1Gu65zfY6DZPetU86Nuea7Lf9Ddfveu6O7+3Qwsj/b9LCfHLOjbz/wo\np7vr8mvJNV/06Guud1nW9hix0hFp+RV5vczWY2eZjYfPMpd36+V+saWXmc0kc0mk32lOL7Nu0izz\nvdwoxyiaZTZ9M/IxLs/pZd4YP8tcluO+fB6YjTn7mB/i/nl8e+AYY4wxEDsg+K1zuyOASmzFMsAN\nwB7g5pB13gXC+0kWEZFYGoCJmU5EPPKwiR4H9MC9UllERLqIE4B/YXMCN2Q4LSIiIiIi4nc3Amux\nxUl/BfbNbHIiugX4JzatjwH9M5uciM4A3gRagUMynBY3sR5U9IM/ABuB1zOdkBj2BZ7Dft5vAH4d\nsKEAeAn7G38L+FlmkxNVLrAaeDLTCYmiCXgNm87s6IkyDn1DpucAv89UQmKYTnDM6kXOy4++CkzG\nXij8FhC8PKjoB8cCU/B/QBgOHOxM98EWzfrxfAL0cv7mAX8HjslgWqK5Gvg/4IlMJySKRmCQlxXd\ne+Lyt+0h032ALZlKSAzLsS2kwN7tjM5gWqJZB7yd6UREcDg2IDQBu4EHgVMymaAIVgKfZToRHmzA\nBlWAHdgc7MjMJSeqz52/PbA3Bv/OYFoiGQ2ciL0p9VuLzbY8pS8bAwLAQuB94Hz8e+cd6kLgqUwn\nIguNAj4Ief+hM08SNw6bq3kpw+mIpBs2eG3E5l792K3ur4BrCN74+ZUBngFeAS6JtqJfA8JybPa7\n7eskZ3k5MAa4F/uhZEqsdIJN65fA/WlPXZCXdPpRzAdppEP6AI8Cc7E5BT/agy3eGg0UAcUZTU17\ns4BN2HJ5v+cOjsYG/xOAK7BFnJ3SGGzlmF9dALyArSTzOz/WIRxBeN9WN+DfiuVx+L8OAWxdzDLg\nykwnJA4LgB9lOhFt/BSbe20EPgF2An/MaIq8qQCSNaK6L0wKmZ4D/ClTCYlhJrY1x5BMJ8Sj54BD\nM52INrLpQcVx+D8g5GAvWpnMVXsxBBjgTPcE6oFvZS45MU3Fv62MehFsiNMbe4NakrnkJN+j2B/e\nGmAJsE9mkxPRO8B72CzlauA3mU1ORN/G3ul8ga10fDqzyWknGx5UfAD4GNiFPZf/ndnkRHQMtihm\nDcHv5cyoW2TGgcAqbDpfw5bT+9lU/NvKaDz2PK7Blqb49TckIiIiIiIiIiIiIiIiIiIiIiIiIiIi\nQadi235/JdMJyQIrCD4tXQP0S3B/xUR+YOgBbNfkVwI/Jr4Hn6YCR3pctwmPPU2GuIuOP4A3r837\nFzq4n1Twcp4rcX+Stj/wvQ4cM9L+7sWO4R6vYuA/BJ/jmN+BffhWXqYT0AXMBqqdv5VJ2F83/NmZ\nVi52TIVEhPZdVJrgvqIZDhxG+FPvbiKd62nYXnf/5uFY8fbH1I0YHZDFcAO2W4WAoxPYV7JVeFgn\n0vkaCHwfuCPOY0baXyL9ZNUBJyewvW/5tXO7zqIP8E3gB8BZzryZwMMh6xQTvIstAV4EXnXW6e3M\nb8L26voqdkCbi7EDXazBPrnd01mvENt3/GvATYR3FX6Ns81aIgemHc52a7AXu8BT4OOAZ51tnyE4\nKNG9wG+dY/4cuAf7g/0btsuJYuA+bE+V94Qc5zfAP7BPTkZKSxMwGLic4N1Yo5MOiHyuZmK7dX4V\n+xS2m1psr6mrsU/v3kvwbrGJ8HNdhu2CZC22g8KxwGXAVSHbhxrs7P8N7J1+aMdn52J7F12NPW+B\n398O4BfY834kNqd0qHOcn4dsfwFQ5Uw/ju298g2CAWQR9ruwmmCXLoHO6x7EdtUccC/wHScNtxD8\nblxKdGcAtzrTc7GfM8AE4Hln+lDnf3gF2xfV8JBjBs7zidjP6RVgMeE5uf2wXak0YLunCfxvhc7/\ndrMzL9J3uhz7dPtKoufMA0HhRuz30+v10O+d2YlPfRf7wwfbH8sh2Dvp9whexO8AzsH231IXMv86\nbKdeYC+EoZ17hRZB3IgNOGBzIoHAcxnBgFAC/M6Z7ob98bn1eLiH4J35zdgfFs765znT/w38xZm+\nF/vIfuAHcg/BXl1PBrYB+zvLXwEOcpYNdP7mYn/4BzrvQzvYazuoRx72HJYS+VwVYLtFL3TmP4R7\nlwJjCe936B7sxTFw3NBz/RG2QzgIFmFVYAdGcbOYYDHCidhzOghbBPQE9n8GGxQD53QPcHrIPgLn\nYQi2C5SAp4CjnOnAOezp/C+B96E3AaHvT8V+XmD7hXofyMcGgMDnnI8N1OMi/G8AwwiOuvUoNsCN\nxHZFvxD7Ob2IDYxgv493O9OB8xz4nMY68+8n+DlVYou5ujv72II9Z20/s0jf6UOxN0QF2D583sH9\ns7oHG5xuIbxbmV8SvAEJfV3rLC8GPsUGoaewwavTUJFRas0m2JHYI877Vdi7ppOxfTGdiL0ATcN+\nuV501u8RMg324hZwIPZOvj82FxLoEfQIglnZB7B3nWB/PCXYLzbYu+mJ2DuoUF9iy+7B3iFPD9nv\nqc70nwnetRrn/wrNfgfu9N7A9o30pvP+TeyFZi32InEJ9vs3AnuxjNUx3GLskKk12K6H3c7VV7AX\n9MBd659xv+ONdYcXeq5fw16wHndesfZxLMGcyVPYgXNysGXnh2IDI9gL+QZnuhX7XWhrC7Aem8t8\nFzu6XeB/nkvwM9kXW/wVbXjEpcDt2HN1Ajag7sJ+Lw4kGJD6Yb8bTRH2sxH7neuD7Zr6fmz31Mc4\n/8NXsTcBzzjr52L7eQrIcdZZj70xAvtdDXxOBntjsxt74d2EDUJtz7fbd3oSNgg8BjQ7rydctg2k\nYwE2oF0WMj9SoA94FXu+P8eex8exIw52CgoIqTMIe5E/APslz3X+XoPNvv8AOwrUP7Bd54Idt+Cc\nCPvbGTJ9L/bC/zr2zmyqh/T8DLgzxjq7Q6b3EP79iHQB/LzN+y9Dtt/VZn+52M62fogtw/8P9k4t\nVvfgF2B/hN8Pmed2rg5q876jWfvQc12KveCdhL2TPtB1C2/HvY/2lb5gL1yRyrQfBM7Ejmz3mDOv\nGBtgjnC2fY7Y57AZW4wzw9nfAyHLfoA9n169iM0p/gtbTHQRtqjramzQf5NgTsZN2/+17fn6MmS6\nlcjXKbfv9Nw2+4v0WRjsb+9QbO4qMOLdr3Afe+FBbK45NAf2NDZ3MQh/jugWN9UhpM7p2K6Gx2Ev\ngmOwd6/HYu/ODsHeJT/orP8StgIwUNwRuONx0wd7d9kdWy4d8HeCd3pnh8xfhh21LVDOPgoYGsf/\n8mLI/r6LLbrpiBzsHdxObHHSMOxdVjSHYgPIeSHz/o77uVqHPd8TnPmzO5jO0PSOwV5IryeYI9tO\n+NjeoeoJBqoTsBcbg83dnE7wvA9y9h3LX7A5gdkEvyv9sBewZuzd9hEh6+8m8gX0Iez34FiCucpl\n2EAb2GYywfGM10XYz0rsjU0d9g59mpOW7dggMTQkTd0JL1YxzjoTCBYZnUUwSES6gLc955G+0/XY\n8xUoMppF5GC7FFs3UYP9XMHWDU1xeQXqLUJzK4c7050iGIACQiqdTbCsPWCJM38PNls80/kLsBl7\nJxxoDhkoAnETyOo+j62YC7gSe5e2Bnux/I8zfzk2a/83bBHIwwR/AKFMm+nA+znYO8K12IAwN8I2\nbvtou+w17EVkHXZw8udxZ7A/tiuwF9XnnO3uxBalXED7c7ULW/RQg83ab3RJQ6S0ucnFVs6+hi3q\nux17Tp/EFgutpn0rnh9jcxRvOOsEikX+ia1bqHXSXEuwsjVaWrZiK+XHECxuWoq9gL+FvUsObe10\np5PeQKVy6L5rnbQtB1qceb939rMKm+O8w/m/o43j8Tz2AlyP/S6/T/Bz/BIb+G4m2M122ya6zdgg\ntNT5n7YR/K6Gfu9CfYqtW3jd2Xek7/RqbOALlPFHK0Yz2HqQu7BFS/lR1g04nWD3+7cRfuMl4is9\nQ6bdApKIV6UEGyukQu+Q6f8l/CZDRJLgGOydy1psMceEqGuLZM6V2Lv5N7G5mWwYZlZERERERERE\nRERERERERERERERERETg/wPtrhmmEGLBzwAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0xae31fa0c>"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "It would appear that high k-values provide the most robust classifiers.\n",
      "It would also appear that first and second derivatives provide the most information.\n",
      "\n",
      "Next, see which combination of statistics is most accurate using 2D KNN."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from sklearn.neighbors import KNeighborsClassifier\n",
      "\n",
      "accs_2D = {}\n",
      "stats = house_stats.columns[1:] # take out the 'owns_ev' column\n",
      "for s1 in range(0, len(stats)-1):\n",
      "    for s2 in range(s1, len(stats)):\n",
      "        stat1 = stats[s1]\n",
      "        stat2 = stats[s2]\n",
      "\n",
      "        X = zip(house_stats.T[train_houses].T[stats[s1]].values, house_stats.T[train_houses].T[stats[s2]].values)\n",
      "        y = house_stats.T[train_houses].T['owns_ev'].values.astype('int')\n",
      "\n",
      "        clf = KNeighborsClassifier()\n",
      "        clf.fit(X, y)\n",
      "\n",
      "        X = zip(house_stats.T[validate_houses].T[stats[s1]].values, house_stats.T[validate_houses].T[stats[s2]].values)\n",
      "        y = house_stats.T[validate_houses].T['owns_ev'].values.astype('int')\n",
      "\n",
      "        pred = clf.predict(X)\n",
      "\n",
      "        accs_2D[(stat1,stat2)] = util.acc_info(pred, y)\n",
      "\n",
      "accs_2D = [(k,v) for k,v in accs_2D.iteritems()]\n",
      "accs_2D.sort(key = lambda (k,v): -v['acc'])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 9
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "for (k,v) in accs_2D[:10]:\n",
      "    print \"key: %s, acc: %.3f (%.3f), pos acc: %.3f (%.3f), neg acc: %.3f (%.3f)\" % (k, v['acc'], v['acc_unc'], v['pos_acc'], v['pos_acc_unc'], v['neg_acc'], v['neg_acc_unc'])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "key: ('d3_k4', 'd2_k5'), acc: 0.874 (0.019), pos acc: 0.714 (0.046), neg acc: 0.945 (0.015)\n",
        "key: ('d1_k3', 'd1_k5'), acc: 0.865 (0.019), pos acc: 0.765 (0.043), neg acc: 0.909 (0.019)\n",
        "key: ('d1_k4', 'd1_k5'), acc: 0.862 (0.019), pos acc: 0.796 (0.041), neg acc: 0.891 (0.021)\n",
        "key: ('d1_k3', 'd1_k4'), acc: 0.852 (0.020), pos acc: 0.745 (0.044), neg acc: 0.900 (0.020)\n",
        "key: ('d2_k1', 'd1_k5'), acc: 0.852 (0.020), pos acc: 0.724 (0.045), neg acc: 0.909 (0.019)\n",
        "key: ('d1_k1', 'd1_k3'), acc: 0.852 (0.020), pos acc: 0.724 (0.045), neg acc: 0.909 (0.019)\n",
        "key: ('d1_k5', 'normed_d2_k2'), acc: 0.852 (0.020), pos acc: 0.776 (0.042), neg acc: 0.886 (0.021)\n",
        "key: ('d3_k1', 'd1_k4'), acc: 0.852 (0.020), pos acc: 0.735 (0.045), neg acc: 0.905 (0.020)\n",
        "key: ('d2_k1', 'd1_k3'), acc: 0.852 (0.020), pos acc: 0.724 (0.045), neg acc: 0.909 (0.019)\n",
        "key: ('d2_k5', 'd3_k5'), acc: 0.852 (0.020), pos acc: 0.765 (0.043), neg acc: 0.891 (0.021)\n"
       ]
      }
     ],
     "prompt_number": 10
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "X = zip(house_stats.T[train_houses].T['d3_k4'].values, house_stats.T[train_houses].T['d2_k5'].values)\n",
      "y = house_stats.T[train_houses].T['owns_ev'].values.astype('int')\n",
      "clf = KNeighborsClassifier()\n",
      "clf.fit(X,y)\n",
      "X = zip(house_stats.T[test_houses].T['d3_k4'].values, house_stats.T[test_houses].T['d2_k5'].values)\n",
      "y = house_stats.T[test_houses].T['owns_ev'].values.astype('int')\n",
      "pred = clf.predict(X)\n",
      "acc = util.acc_info(pred, y)\n",
      "print util.acc_str(acc)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "acc: 0.8145 (0.0218), pos acc: 0.5476 (0.0543), neg acc: 0.9103 (0.0187)\n"
       ]
      }
     ],
     "prompt_number": 11
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "After doing cross-validation, it appears the accuracy of using 2D is only around 81%.\n",
      "It turns out 1D KNN is better than any pair of 2D KNN!"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Part b - Using the Markov Model"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The Markov Model from Part II was put into a file `markov.py`, which provides an easy interface for training and running."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from markov import EVMarkovModel\n",
      "\n",
      "random.seed(0)\n",
      "# want to only train the Markov Model with houses that own EVs, while also keeping the ratio (with/without ev) correct in the test data\n",
      "train_markov, train_classifier, test_classifier = util.split_data(houses, ratios=[0.4, 0.4, 0.2])\n",
      "run_markov = train_classifier.append(test_classifier)\n",
      "#train_markov, train_classifier, test = util.split_data(houses, ratios=[0.4, 0.4, 0.2])\n",
      "\n",
      "# train the markov model transition probabilities\n",
      "model = EVMarkovModel(usages, when_charging)\n",
      "model.train(train_markov, normed=True)\n",
      "\n",
      "# run the markov model on houses to be classified (this part takes a while)\n",
      "preds = pd.DataFrame()\n",
      "probs = pd.DataFrame()\n",
      "for i,house in enumerate(run_markov):\n",
      "    #print \"house %i out of %i\" % (i+1, len(run_markov))\n",
      "    (p0, p1) = model.run(house)\n",
      "    probs[house] = p1\n",
      "    preds[house] = (p1 > 0.5).astype('float')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 12
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now that the probabilities of charging are known at each timestep, it should be possible to perform aggregate statistics to see which houses own EVs.\n",
      "For example, houses with EVs should report more charging than houses without EVs.\n",
      "The plot below looks at a histogram of the average fraction of time the Markov Chain reported charging, separated into houses with or without EVs."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "avg_probs_with_ev = probs[train_classifier.intersection(houses_with_ev)].apply(np.average)\n",
      "avg_probs_without_ev = probs[train_classifier.intersection(houses_without_ev)].apply(np.average)\n",
      "\n",
      "bins = np.linspace(0, 1, 101)\n",
      "xs = (bins[1:] + bins[:-1])/2\n",
      "ys1 = np.histogram(avg_probs_with_ev, bins=bins)[0]\n",
      "ys2 = np.histogram(avg_probs_without_ev, bins=bins)[0]\n",
      "pl.plot(xs, ys2, linestyle='None', marker='.', color='blue')\n",
      "pl.plot(xs, ys1, linestyle='None', marker='.', color='red')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 13,
       "text": [
        "[<matplotlib.lines.Line2D at 0xae00eb4c>]"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXMAAAEACAYAAABBDJb9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAD+xJREFUeJzt3W2MXNdBh/FnE7stQ+tszK4cKy+15BJCo7Yu3bQBUvm2\nFJquII1AGwlQG2jF8IlGSJHj8IEsH/oSI6EKRSBtq1ZWxRYcSkMCGLwUj0OApjKKU8clLDGxVUps\nx443FJUXmy4fzl3vdHd298ydmTt3zjw/abR3Zu69c/bI+/edc+45ByRJkiRJkiRJkiRJkiRJkrpm\nJHK/U8B/AP8HXALeCWwF/gh4Y/7+PcBC10soSeqaFwnh3WwfsCfffgD4VKklkiS17UXgB1a89jyw\nLd++Ln8uSaqwfwWeAY4Cv5K/drHp/ZEVzyVJFbQ9/zkOHAPezerwfqXUEkmSrtgUud9L+c+XgS8T\nOkDPEppXzhDC/tzKg3bu3Ll48uTJLhRTkobKSeBN7RxwVcQ+NeAN+fb3Az8FHAceB+7NX78XeGxV\naU6eZHFx0cfiIg899FDfy1CVh3VhXVgX6z+Ane0EOcRdmW8jXI0v7f8HwCFC+/kB4KMs35ooSeqD\nmDB/EdjV4vVXgPd1tziSpCJimlnUBVmW9bsIlWFdLLMullkXnYkdAVrUYt7+I0mKNDIyAm3ms1fm\nkpQAw1ySEmCYS1ICDHNJSoBhLkkJMMwlKQGGuSQlwDCXpAQY5pKUAMNckhLQ8zDPMpichIU2lnqu\n14sdJ0nDqudzs0CYm2VqCg4ciDsoy+DIEdo+TpJSUNm5WSYmYGYmfv9ardhxkjSseh7mU1MwNwej\no/HHzM4WO06ShpVT4EpSxVS2mUWS1FuGuSQlwDCXpAQY5pKUAMNckhJgmEtSAgxzSUqAYS5JCTDM\nJSkBhrkkJcAwl6QEGOaSlADDXJISYJhLUgIMc0lKgGEuSQkwzCUpAYa5JCVg4MK8Xocsg8lJWFjo\nd2kkqRpiw/xq4Bngifz5VmAOmAcOAaUtuzw/D0eOwMGDIdglSfFhfh/wDWBpdea9hDC/GfhK/rwU\ntVr4OTEBMzNlfaokVVtMmN8ATAKfZXm16LuA/fn2fuDu7hettdlZmJqCuTkYLe37gCRV28jGu/Ao\n8AlgC3A/8DPAReDapnO80vS82eLi4mKLlyVJaxkZGYG4fL5i0wbv/zRwjtBenq2xzyLLzS+rTE9P\nX9nOsowsW+s0kjScGo0GjUajo3NslPyfAD4EXAZeR7g6/xPgNkK4nwG2A4eBW1oc75W5JLWpyJV5\nOzvvZrmZZR9wAXiY0Pk5SutOUMNcktpUJMzbvc98KZk/Bfwk4dbE9+bPJUl90lbyF+CVuSS1qYwr\nc0lSBRnmkpQAw1ySEmCYS1ICDHNJSoBhLkkJMMwlKQGGuSQloPJh7spCkrSxyoe5KwtJ0sYqH+au\nLCRJG6v83CwLC+GKfGbGlYUkDYdeT4FbRGkTbdXroUmmVgtLyxn8kgbVUE+0Zdu6pGGWTJjbti5p\nmCXTzGLbuqRUDHWbuSSlYqjbzCVpmBnmkpQAw1ySEmCYS1ICDHNJSoBhLkkJMMwlKQGGuSQloBJh\nXnQBCheukKSgEmFedJIsJ9eSpKASYV50kiwn15KkoBJzsxSdJMvJtSSlyIm2JCkBwzfRlj2gkgQM\nepjbAypJwKCHuT2gkgQMepu5PaCSEmQHqCQloBcdoK8DngaOAd8APpm/vhWYA+aBQ4CXxZLURzHJ\nXwO+A2wCngLuB+4CzgP7gAeAa4G9LY71ylyS2tSrWxO/k/98DXA1cJEQ5vvz1/cDd7fzoZKk7ooJ\n86sIzSxngcPACWBb/pz857aelE6SFGVTxD7fBXYB1wB/BbxnxfuL+UOS1CcxYb7kVeDPgXcQrsav\nA84A24Fzax00PT19ZTvLMrIsK1BMSUpXo9Gg0Wh0dI6NGtjHgMvAAvB9hCvz3wLeD1wAHiZ0fI5i\nB6gkdUUv7jN/C6GD86r88QXgtwm3Jh4AbgJOAfcQAn8lw1yS2uSgIUlKwPDNmihJAgxzSUqCYS5J\nCTDMJSkBpYa5CwNJUm+UGuYuDCRJvVFqmLswkCT1Rqn3mbswkCRtbKAGDdXrodmlVoPZ2e6H+5O3\n1NlyZp5Lm2vcfHSWa97o/x6SBsNADRrqdfv5ljPz7Hr1CLedP8iJO2ygl5S2voV5r9vPL20OH3Ci\nNsGtT9lALyltfWtm6XX7+aunFzhxR51bn5qxiUXSQBmoNnNJUmsD1WYuSeoew1ySEmCYS1ICDHNJ\nSoBhLkkJMMwlKQGGuSQlwDCXpAQY5pKUAMNckhJgmEtSAgxzSUqAYS5JCUgnzOt1yDKYnAzz68a9\nJUlJSCfM11m6qNerGklSv6UT5ussXdTrVY0kqd/SWZxinaWLer2qkSR1kysNSVIChmOlIXszJWmV\nwQtzezMlaZXBC3N7MyVplcFrM7c3U1Li7ACVpAT0qgP0RuAwcAJ4DvhY/vpWYA6YBw4BXiZLUp/E\nJP91+eMY8HrgH4G7gV8GzgP7gAeAa4G9K471ylyS2tSrK/MzhCAH+E/gn4DrgbuA/fnr+wkBL0nq\ng3bvZtkBvB14GtgGnM1fP5s/lyT1QTth/nrgS8B9wLdXvLeYPyRJfbApcr/NhCD/AvBY/tpZQlv6\nGWA7cK7VgdPT01e2sywjy7JiJZWkRDUaDRqNRkfniGlgHyG0iV8Afr3p9X35aw8TOj5HsQNUkjrW\nq/vM7wCeBL7OclPKg8DXgAPATcAp4B5g5WQphrkktclBQ5KUgOGYNVGStIphLkkJMMwlKQGGuSQl\nwDCXpAQY5pKUAMNckhJgmEtSAqof5vU6ZBlMToYl4yRJq1Q/zOfn4cgROHgwBLskaZXqh3mtFn5O\nTIRFnCVJq1R/bpaFhXBFPjMDoy4zKil9TrQlSQlwoi1JGlKGuSQlwDCXpAQY5pKUAMNckhJgmEtS\nAgxzSUqAYS5JCdjU7wKUol4Pc7zUatw3Psuzp0ep1WB21kGlktIwHCNAsyxM1gUcHp/ivS8fAGBq\nCg4c6GO5JKkFR4CupWmyrt9728zSpvN2SUrGcFyZN03WtcCo83ZJqjQn2pKkBNjMIklDyjCXpAQY\n5pKUAMNckhJgmDep18Mt6ZOT4QYYSRoUhnmT+fkwtujgwRDskjQoDPMmTWOLHFAkaaB4n3mTprFF\nDiiS1DcOGpKkBDhoKIa9nJISFBPmnwPOAsebXtsKzAHzwCFgcBol7OWUlKCYMP88cOeK1/YSwvxm\n4Cv588FgL6ekBMW2yewAngDekj9/HthNuGK/DmgAt7Q4rnpt5vZySqq4XnaA7uB7w/wicG3TOV5p\net6sWJg3rQzU6+WASvwoSYpSJMy7sWzcYv5oaXp6+sp2lmVkWbbxGZfatSGkbQ+XAyrxoySppUaj\nQaPR6OgcnTSzZMAZYDtwmG42s0xOhg7KiQmYm+vp5XKJHyVJUcq8NfFx4N58+17gsYLnaW12NizQ\nWUK6lvhRktQzMcn/RUJn5xihw/M3gT8FDgA3AaeAe4BWN21XrwNUkirOEaCSlABHgErSkDLMJSkB\nhrkkJcAwl6QEGOaSlADDXJISYJhLUgIMc0lKQDcm2iqmebrC8XE4fXpgpi50pkVJVdO/MG+ernBs\nDM6fD9sDMHWhMy1Kqpr+NbM0r/iza9fy9gCs/uNiRZKqpn9zszSv+AMDtfqPixVJ6iUn2uqi2HZx\n288ldZsTbXXRUrv4wYMhsDvdT5J6yTBfQ2y7uO3nkqrAZpY1xLaL234uqduq2Wa+e7cNypLUhmqG\n+dLW1JQ3ZEtShOqG+cSEKyZLUqRq3s0yNWWQS1KP2QEqSRVTzSvzqqvXIctgcjLcmtLZbh0fI0lF\nGOaRo36KDA5yQJGkshjmkaN+igwOckCRpLLYZh456qfI4CAHFEkqopq3JlY9zCWpYuwAlaQhZZhL\nUgIMc0lKgGEuSQkoN8wHaRRNwbJWcXDRIFW7pGLKDfNBGkVTsKxVHFw0SNUuqZhyw3yQRtEULGsV\nBxcNUrVLKqbc+8wHaRRNwbJWcXDRIFW7JAcNda5eD20SrVZGWu+9AueLPV3zfuPjcPp0OQs3Ffl1\nJXVHkTDv1J3A88C/AA+0eH9xoOzevbgYFtRYXJyain+vwPliT9e839hY+0UoqsivK6k7gLavgjtp\nM78aeCQP9DcDPw/8cAfn67/1Gpc7bAxvfPjDhU7XvN+uXe0XoahetrM3Go3unnCAWRfLrIvOdBLm\n7wReAE4Bl4A/BD7YhTL1z+zs2isjrfdexPkaR48WOl3zfo8+Wt7CTUV+3Vj+0S6zLpZZF53Z1MGx\n1wPfbHr+b8C7OitOn42Orr3o9HrvFThf7OlW7lfWmthFfl1J/dNJmA9Qz2aXrdcruWdPy57D0488\nwbFPN7i0ucbNR2e55uN7Wp7j718Yp3budNv7Pfv+PWw5M8+lzTX+e8s4b7jQ+hzNx63cr/kcb33f\nOK99qb3PWu98ze+dvvyt76mLTs/X7f3K/Kyq10WZdZZKXXSjzoropLf0dmCa0GYO8CDwXeDhpn1e\nAHZ28BmSNHS28drLZ/mfzWV93ibgJLADeA1wjEHvAJWkIfUB4J8JV+AP9rkskiRJ0nDbaOAQwO/m\n7z8LvL2kcvXLRvXxi4R6+Drwd8BbyytaqWL+XQDcBlwGfraMQvVJTF1kwDPAc0CjlFL1x0Z1MQb8\nJaHZ9jngl0orWbk+B5wFjq+zT6m5eTWhiWUHsJnW7eaTwF/k2+8CvtrrQvVRTH38KHBNvn0nadZH\nTD0s7fc3wJ8BP1dW4UoWUxejwAnghvz5WFmFK1lMXUwDn8y3x4ALdHbXXVW9mxDQa4V527nZ6ayJ\nMQOH7gL259tPE/7hbuvwc6sqpj7+AXg1336a5T/glMQOKPs14I+Bl0srWfli6uIXgC8RxmoAnC+r\ncCWLqYuXgC359hZCmF8uqXxl+lvg4jrvt52bnYZ5q4FD10fsk2KAQVx9NPsoy//7piT238UHgd/P\nn6c6biGmLn4Q2AocBo4CHyqnaKWLqYvPALcC/05oXrivnKJVTtu52enXl9g/wJX3s6f6h9vO7/Ue\n4CPAj/eoLP0UUw+fBvbm+45Q8gxxJYqpi83AjwA/AdQI396+SmgvTUlMXfwGofklI4xRmQPeBny7\nd8WqrLZys9Mw/xZwY9PzG1n+qrjWPjfkr6Uopj4gdHp+htBmvt5XrUEVUw/vIHzNhtA2+gHCV+/H\ne166csXUxTcJTSv/lT+eJARYamEeUxc/Bnw83z4JvAj8EOEbyzApPTdjBg41N+TfTpodfkti6uMm\nQrvh7aWWrFztDij7POnezRJTF7cAf03oIKwROsXeXF4RSxNTF78DPJRvbyOE/daSyle2HcR1gJaW\nm60GDv1q/ljySP7+s4SvkynbqD4+S+jUeSZ/fK3sApYk5t/FkpTDHOLq4n7CHS3HgY+VWrpybVQX\nY8AThKw4TugcTtEXCf0C/0v4ZvYRhjs3JUmSJEmSJEmSJEmSJEmSJEmSJK30/w9qTeVigNj8AAAA\nAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0xabb22d6c>"
       ]
      }
     ],
     "prompt_number": 13
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Although houses with EVs do report more charging, it doesn't seem signficant enough for an accurate classifier.\n",
      "Below is an attempt at classification using this statistic with 1D KNN, but it is not as accurate as the aggregate models used in part a."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from sklearn.neighbors import KNeighborsClassifier\n",
      "\n",
      "X = preds[train_classifier].apply(np.sum).values\n",
      "y = (when_charging[train_classifier].sum() > 0)\n",
      "\n",
      "clf = KNeighborsClassifier()\n",
      "clf.fit(np.array([X]).T, y)\n",
      "\n",
      "X = preds[test_classifier].apply(np.sum).values\n",
      "y = (when_charging[test_classifier].sum() > 0)\n",
      "pred = clf.predict(np.array([X]).T)\n",
      "\n",
      "accs = util.acc_info(pred, y)\n",
      "\n",
      "print \"acc: %.3f (%.3f), pos acc: %.3f (%.3f), neg acc: %.3f (%.3f)\" % (accs['acc'], accs['acc_unc'], accs['pos_acc'], accs['pos_acc_unc'], accs['neg_acc'], accs['neg_acc_unc'])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "acc: 0.642 (0.027), pos acc: 0.330 (0.048), neg acc: 0.772 (0.028)\n"
       ]
      }
     ],
     "prompt_number": 14
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Conclusions"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The most accurate model so far for predicting which houses own EVs has an accuracy of about 84% with an uncertainty of 2%.\n",
      "This was the 1D KNN with first derivative (weighted k=5).\n",
      "\n",
      "All of the most accurate 1D KNN or 2D KNN pairs seemed to rely on first or second derivatives with high k-values.\n",
      "This implies the tails of the distributions had the most significance.\n",
      "This isn't too surprising, since a house with an EV would have more sudden jumps, hence it would have many timeseries with large first or second derivatives.\n",
      "A different approach which extracts or fits the tails could result in a higher accuracy.\n",
      "\n",
      "Using the Markov Models were not that accurate, but the Markov Model was only used in one way (average time spent charging).\n",
      "It could be possible to use the results from the Markov Models in a different way, for example looking at distributions for how long each charge lasted (turn on to turn off), or the total number of charge groups, or the frequency of the charges.\n",
      "Real charges might have different behavior, which could still be encoded in the results of the Markov Model.\n",
      "\n",
      "Finally, combining the two approaches could result in a more accurate multi-dimensional classifier. For example, it might be worth re-running the first 2D KNN with the Markov Model's output as one of the statistics."
     ]
    }
   ],
   "metadata": {}
  }
 ]
}