{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Introduction\n", "\n", "This Jupyter Notebook introduces [Kernel Density Estimation](https://en.wikipedia.org/wiki/Kernel_density_estimation) (KDE) alongside with [KDEpy](https://kdepy.readthedocs.io/en/latest/).\n", "Kernel density estimation is an approach to solve the following problem.\n", "\n", "> **Problem.** Given a set of $N$ data points $\\{x_1, x_2, \\dots, x_N\\}$, estimate the probability density function from which the data is drawn.\n", "\n", "There are roughly two approaches to solving this problem:\n", "\n", "1. Assume a **parametric form**, e.g. a [normal distribution](https://en.wikipedia.org/wiki/Normal_distribution), and estimate the *parameters* $\\mu$ and $\\sigma$. These parameters uniquely determine the distribution, and are typically found using the [maximum likelihood principle](https://en.wikipedia.org/wiki/Maximum_likelihood_estimation). The advantage of this approach is that we only need to estimate a few parameters, while the disadvantage is that the chosen parametric form might not fit the data very well.\n", "2. Use **kernel density estimation**, which is a non-parametric method -- we let the data speak for itself. The idea is to place a *kernel function* $K$ on each data point $x_i$, and let the probability density function be given by the sum of the $N$ kernel functions.\n", "\n", "Assuming a parametric form is a perfectly valid approach, especially if there is evidence suggesting the presence of a theoretical distribution.\n", "However, for [exploratory data analysis](https://en.wikipedia.org/wiki/Exploratory_data_analysis) we might want to assume as little as possible when plotting a distribution, and this is where kernel density estimation and KDEpy comes to the rescue." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "from scipy.stats import norm\n", "import matplotlib.pyplot as plt\n", "from KDEpy import FFTKDE # Fastest 1D algorithm\n", "\n", "np.random.seed(123) # Seed generator for reproducible results\n", "\n", "distribution = norm() # Create normal distribution\n", "data = distribution.rvs(32) # Draw 32 random samples" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## The histogram\n", "\n", "A KDE may be thought of as an extension to the familiar [histogram](https://en.wikipedia.org/wiki/Histogram). \n", "The purpose of the KDE is to estimate an unknown probability density function $f(x)$ given data sampled from it. \n", "A natural first thought is to use a histogram – it’s well known, simple to understand and works reasonably well.\n", "\n", "To see how the histogram performs on the data generated above, we'll plot the true distribution alongside a histogram. \n", "As seen below, the histogram does a fairly poor job, since:\n", "- The location of the bins and the number of bins are both arbitrary.\n", "- The estimated distribution is discontinuous (not smooth), while the true distribution is continuous (smooth)." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xd8VFX6x/HPkxAIEJYgRRRcQUQUhEQJiHVpKkpRFH4CgouuIiqLdRUbKOqKoGBDEAVZC2KDBYVdLMBaEVGKAiqICFgQQpEWSvL8/jgJBgzJJJnJmbnzvF+veZGbe3PnmyF5cubcc88RVcUYY0ywJPgOYIwxJvysuBtjTABZcTfGmACy4m6MMQFkxd0YYwLIirsxxgSQFXdjjAkgK+7GGBNAVtyNMSaAyvl64ho1ami9evV8Pb0xxsSkzz//fKOq1izqOG/FvV69eixYsMDX0xtjTEwSkR9COc66ZYwxJoCsuBtjTABZcTfGmADy1udujIl9e/fuZd26dWRlZfmOEjjJycnUrVuXpKSkEn29FXdjTImtW7eOKlWqUK9ePUTEd5zAUFUyMzNZt24d9evXL9E5rFvGGFNiWVlZVK9e3Qp7mIkI1atXL9U7IivuxphSscIeGaV9Xa24G2NMAFmfu4lppw+bzY9bdoX1nHVSK/LRoLZhPaeJjMzMTNq1awfAL7/8QmJiIjVrups358+fT/ny5cs0z1133UWNGjW44YYbWLZsGT179iQhIYGpU6dS1nfkW3E3Me3HLbtYPaxjWM9Zb9CMsJ7PRE716tVZtGgRAPfccw8pKSnccsstBxyjqqgqCQll21ExZcoUunXrxt13312mz5vHumWMMYGzcuVKTjzxRPr378/JJ5/M2rVrSU1N3b9/8uTJXHnllQCsX7+eiy66iIyMDFq2bMm8efP+cL5nn32Wrl27cu6559KoUSPuv//+/fuGDh1Ko0aNOPvss1mxYgUA06dP58knn2Ts2LG0b98+wt9twazlbowJn9at//i5//s/uPZa2LkTzj//j/v79nWPjRuhW7cD982dW+Ioy5Yt47nnnmPs2LHs27fvkMcNHDiQW2+9lVatWrF69Wo6derEV1999Yfj5s+fz1dffUX58uVp0aIFnTp1Ys+ePbzxxhssWrSIPXv2kJ6ezqmnnkqXLl2YP3/+/i4aH6y4G2MCqUGDBrRo0aLI4959912++eab/dubN29m165dVKxY8YDjzj33XKpVqwbAhRdeyIcffkhWVhYXX3wxFStWpGLFinTu3Dm830QpWHE3xoRPYS3tSpUK31+jRqla6gerXLny/o8TEhJQ1f3b+cePq2pIF18PHpqYtx2tQ0Gtz90YE3gJCQlUq1aNFStWkJOTw9SpU/fva9++PaNHj96/nXeB9mBvv/02W7ZsYefOnUybNo3TTz+ds846iylTppCVlcVvv/3GW2+9FfHvJVTWcjfGxIWHHnqIDh068Oc//5nGjRuze/duAEaPHs0111zDc889x759+2jTps0BxT7PGWecQa9evfjuu+/o06cP6enpAHTt2pW0tDTq1avHWWedVabfU2Ek/1uVspSRkaG2WIcprXqDZkRkKGS4zxlUy5cv54QTTvAdI+KeffZZvvrqKx599NEyfd6CXl8R+VxVM4r6WuuWMcaYALJuGWOMKULemPhYYi13Y4wJICvuxhgTQCEVdxHpICLfiMhKERlUyHHdRERFpMjOfmOMMZFTZHEXkURgNHAe0BjoKSKNCziuCjAQ+DTcIY0xxhRPKBdUWwIrVXUVgIhMBi4Alh103H3AcOAWjDFxKdxTMIcy/XJKSgrbt2/fvz1x4kQWLFiwf+KuSpUqcdlllxX4tXPnzqV8+fKcdtppYcscLUIp7nWAtfm21wGn5D9ARE4CjlLVt0TEiruJD6rwxhuwdq2bFCs1FRo2hJYt3cdxKNxTMJd2+uX+/fsXun/u3LmkpKSEpbhnZ2eTmJhY6vOESyh97gVNnLD/zicRSQBGATcXeSKRfiKyQEQWbNiwIfSUxkSLHTtg/nz3sQjccAPcdBPcdRcMGADnnguXX+43o9nvnnvu4eGHHwbg8ccfp3HjxjRr1owePXqwevVqxo4dy6hRo0hPT+eDDz7ghx9+oF27djRr1ox27dqxZs0aAL777jtatWpFixYtGDx4MCkpKYD749CmTRt69epF06ZNATepWPPmzWnSpAnjxo3bnyUlJYXbbruN5s2b0759e+bPn0/r1q055phjmD59evi/+byJ7A/1AE4FZuXbvh24Pd92VWAjsDr3kQX8BGQUdt7mzZurMaV19G1vlc05s7NVn3pKtWZN1Ro1VHfudJ9fuVJ10ybVrCzVn39Wfecd1UWL3L5Vq1RbtFD9+OOwZ4wWy5YtO2A73P8foZwvISFB09LS9j+OOuoove6661RVdciQITpixAhVVT3iiCM0KytLVVU3b978h/2qqp06ddKJEyeqqur48eP1ggsuUFXVjh076qRJk1RVdcyYMVq5cmVVVZ0zZ45WqlRJV61atf8cmZmZqqq6c+dObdKkiW7cuFFVVQGdOXOmqqpeeOGFevbZZ+uePXt00aJFmpaWVuD3dvDrm3ueBVpE3VbVkFrunwENRaS+iJQHegD7/8yo6lZVraGq9VS1HjAP6KKqNreACYY1a6BtWzcn+YknwrRpkDcdbIMGUK0aVKgAtWtD+/aQlub2/fKLe5x+umvd79nj73sIsIoVK7Jo0aL9j6FDhxZ4XLNmzbj00kt58cUXKVeu4B7pTz75hF69egHQp08fPvzww/2f7969O8D+/XlatmxJ/fr1928//vjjpKWl0apVK9auXbt/AY/y5cvToUMHAJo2bcpf/vIXkpKSaNq0KatXry75C3AIRRZ3Vd0HDABmAcuBV1V1qYgMFZEuYU9kTDT55RfIyICFC2H8eHjvPQi1f/bUU2HpUvdHYdQoaNcO1q+PbF5zSDNmzOC6667j888/p3nz5oUu4JEnlOl8808tPHfuXN59910++eQTFi9ezEknnbR/euGkpKT950tISKBChQr7Pw4lS3GFNM5dVWeq6nGq2kBVH8j93GBV/UNHkaq2tla7CYzatWHgQNfPfsUVrp+9OKpUgSefhJdfhs8/hwcfjExOU6icnBzWrl1LmzZtGD58OFu2bGH79u1UqVKFbdu27T/utNNOY/LkyQC89NJLnHHGGQC0atWKN954A2D//oJs3bqVatWqUalSJb7++usCl+wrKza3jDEFmToVjjsOmjRxF0tLq0cPaNYMjj229OeKYnVSK4Z1gfE6qRWLPigE2dnZ9O7dm61bt6Kq3HjjjaSmptK5c2e6devGtGnTeOKJJ3j88ce54oorGDFiBDVr1uS5554D4NFHH6V379488sgjdOzYkapVqxb4PB06dGDs2LE0a9aMRo0a0apVq7DkL5FQOuYj8bALqiYcInFB9a/dhqgmJqp27Rr2c6uq6k8/uXP//HNkzl+GCrrgF0Q7duzQnJwcVVV9+eWXtUuXLmXyvKW5oGotd2PyW7iQ0dMechdF//WvyDzHjz/CrFnQsaNbVq5Klcg8jwmbzz//nAEDBqCqpKamMmHCBN+RimTF3Zg8P/4InTqxJbkKld98M3JFNyMDXn0VLrgA+vSBKVMgwebwi2Znnnkmixcv9h2jWOwnypg8//wnbN3K37oNhiOPjOxzdewIjzzihlWOGBHZ5zJxyYq7MXlGjoTZs/m6Vv2ijw2HgQPhkkvgxRchdz1PY8LFirsxixbBli3uRqSWLcvueUXg2Wfh44/dcxsTRlbcTXzbuhU6d3YtaB9SUlzf/o4dkbuAa+KSFXcT3266CX7+Ge6/32+OceOgb183y6QplsTERNLT02nSpAlpaWmMHDmSnJycQr9m9erVTJo0qYwS+mHF3cSvmTNhwgS49VZo0cJvlgEDoHlz6N8fNm70myXG5M0ts3TpUt555x1mzpzJvffeW+jXWHE3Jqh27HCFtHFjGDLEdxpISoLnnoPNm2HQIVeyDIbWrd0jAmrVqsW4ceN48sknUVVWr17NmWeeycknn8zJJ5/Mxx9/DMCgQYP44IMPSE9PZ9SoUYc8LpbZOHcTn7Ztc9MBDBoUPRczmzaFG2+Ehx92c8KffrrvRDHpmGOOIScnh19//ZVatWrxzjvvkJyczIoVK+jZsycLFixg2LBhPPzww7z11lsA7Ny5s8DjYpkVdxOfateG3F/sqDJkCPz0E9Ss6TtJTHN36cPevXsZMGAAixYtIjExkW+//bbA40M9LpZYcTfxRRUeeAB69YJjjvGd5o9SUuCll3yniGmrVq0iMTGRWrVqce+993L44YezePFicnJySE5OLvBrRo0aFdJxscT63E18efNNuPtud2doNFuzBq680o2/NyHbsGED/fv3Z8CAAYgIW7du5YgjjiAhIYEXXniB7OxsgD9M9Xuo42KZtdxN/Ni3D267DY4/3o1OiWabNrmRPIcdBsOH+04T1Xbt2kV6ejp79+6lXLly9OnTh5tuugmAa6+9losvvpjXXnuNNm3a7F9Yo1mzZpQrV460tDT69u17yONimRV3Ez8mToSvv3ZztScl+U5TuPR0+Otf4bHH4JproH4ZTYlQFubODevpCmtlN2zYkCVLluzffjB3sZSkpCTee++9A44t6LhYZt0yJj7s3OkuVp56qpuNMRbcfz8kJsIdd/hOYmKQFXcTH/buhW7d4KGHir9Uni916sDNN8PkyfDFF77TmBhj3TImPlSt6ro4Ys0tt8CuXZGfgrgUVDWkhaRN8eQN5ywpa7mb4Bs/HmbP9p2iZKpWdTc11a7tO0mBkpOTyczMLHUhMgdSVTIzM0s1JNNa7ibYNm6E66+HLl2gbVvfaUruo4/cTVdRdqGvbt26rFu3jg0bNviOEjjJycnUrVu3xF9vxd0E26OPuoupd97pO0npfPwxDBvmpic+7TTfafZLSkqifpBG8gSIdcuY4Nq8GR5/3F1IbdLEd5rSufZaqFXL3YBlTAisuJvgeuwxN0HYXXf5TlJ6lSu7G7Bmz3ateGOKYMXdBFft2nD11W72xyC4+mp3x+qwYb6TmBhgfe4muPr3950gvCpXhnvugawsNwGaDT80hbDiboJn507497+he/fon2aguP7+d98JTIywbhkTPBMnwqWXwqef+k4SGXv3wosvwg8/+E5iopgVdxMs2dkwciS0bBnclYx+/RWuuAJGjPCdxEQxK+4mWP79b/juO/jHP4LbJ12nDvTp4+68Xb/edxoTpay4m+BQda3ZY46Brl19p4msW291F1afesp3EhOlrLib4Ni0CbZvh5tuclPlBlmjRtCpE4wZ44q8MQex4m6Co3p1+PJLNx48Htx4o7trdc0a30lMFLKhkKZMnD5sNj9u2RW5J/j1V0hOhj/9CcrFyY91mzbw5Zec/tAcftyyIuynr5NakY8GxfBka3EuTn4LjG8/btnF6mEdw37eeoNmuA+GDHHL5/3wA1SoEPbniUq5F4y3rM9k9fUtoF69sJ5+/2trYlJI3TIi0kFEvhGRlSIyqID9/UXkSxFZJCIfikjj8Ec15hB++w1eeAHOOy9+CnseVaY/f5ObWMyYfIos7iKSCIwGzgMaAz0LKN6TVLWpqqYDw4GRYU9qzKG88ALs2OEWko43Ikxr/Bf4z39g+XLfaUwUCaXl3hJYqaqrVHUPMBk4YIVhVf0t32ZlwJZlMWVD1Y0Yad4cWrTwncaLl9Jz37HE4jKCJmJCKe51gLX5ttflfu4AInKdiHyHa7kPLOhEItJPRBaIyAJbucWEQ5P138HSpa7VHtSbloqQWTkVeveG55+HzEzfcUyUCKW4F/Qb84eWuaqOVtUGwG1AgRNoq+o4Vc1Q1YyaNWsWL6kxBVha+1hYuBB69vQdxa/rr3cLac+c6TuJiRKhjJZZBxyVb7su8FMhx08GxpQmlDHFkp7uO4F/TZvCN9/Accf5TmKiRCgt98+AhiJSX0TKAz2A6fkPEJGG+TY7AuEfdGvMwUaO5JEZI2HfPt9JokNeYVe75GVCaLmr6j4RGQDMAhKBCaq6VESGAgtUdTowQETaA3uBzcBfIxnaGLKz4bHHqC2p8XPTUihuuQVWrYIpU3wnMZ6F9FuhqjOBmQd9bnC+j68Pcy5jCvef/8CaNbx4QS8COrFvyaSk/D4zZoMGvtMYj2xuGRObxoyBI47gnYatfCeJLv36uUnTxthlr3hnxd3Enu+/dy33K69kX6J1yRzgyCPddMcTJrjlBk3csuJuYk9yspsR8aqrfCeJTtddB5s3w+TJvpMYj6zZY2LPEUfAI4/kbizxGiUqnXUWPPigmzXSxC0r7ia2vP8+7N4N7dvH7R2pRRKBQX+Y38/EGeuWMbHljjvcDIg2lrtoc+bYfDNxzIq7iR1LlsBHH0H//pBgP7pFmjLFrbW6caPvJMYD+w0xsWPMGHcxtW9f30liQ//+sGcPTJzoO4nxwIq7iQ3btsGLL8Ill7i1Uk3RmjSBM86Ap5+GnBzfaUwZs+JuYsM330CVKvG5IEdp9O8PK1e6/ncTV6y4m9iQkeHWR23Z0neS2HLxxZCW5sa9m7hiQyFN9MvMhKpVISnJd5LYk5zs5ru3YaNxx1ruJvoNHOjmbLfhjyUjAnv3wrff+k5iypAVdxPdfv0VXnsN2rWz1mdpXHopnHOOmyrZxAUr7ia6TZjgWp39+/tOEtu6d3fXLGbN8p3ElBEr7iZ6ZWe7YXytW8MJJ/hOE9suuABq1XKvp4kLVtxN9Jo9G1avtuGP4VC+PPztb/DWW7B2re80pgzYaBkTvdq1c90IrVv7ThIMV10Fw4bBK6+45fiKUL5cAvUGzQhrhDqpFfloUNuwntMUzIq7iV4JCe4ioAmP+vXhs8/gpJNCOnzPvhxWD+sY1gjh/mNhDs26ZUx0eughN22tDX8Mr+bNbdK1OGH/yyb67N7tFuP4+msb/hgJw4ZB796+U5gIs+Juos+UKbBhg11IjZSsLJg0ya1FawLLiruJPk89BQ0awNln+04STFde6d4RPfOM7yQmgqy4m+jy5Zfw4Ye2IEck1a0LnTrB+PFuvncTSPbbY6JLhQpw2WVw+eW+kwRb//5uaodp03wnMRFiQyFNdDnuOPjXv3ynCL5zzoHrr4dGjXwnMRFixd1Ej/ffd1P7pqX5ThJ8iYnw6KO+U5gIsuJuooMqDBjgis4XX9gQyLKyaBF8951b1MMEivW5m+jw8cfuYuo111hhL0sPPAD9+rnhkSZQrLib6DBmDPzpT9Crl+8k8aV/f9i0CV5/3XcSE2ZW3I1/Gza4BTkuuwxSUnyniS9t2kDDhjYVcABZcTf+LVzopqS1O1LLXkICXH21u7fgyy99pzFhZMXd+HfOObB+PTRu7DtJfOrbF2rWhKVLfScxYWSjZYxfv/0GVapApUq+k8Sv6tXhxx8hKcl3EhNG1nI3fvXqBR06+E5hkpLccNRff/WdxIRJSMVdRDqIyDcislJEBhWw/yYRWSYiS0TkPRE5OvxRTeB8/z3MnAmnnOI7iQE3DXCbNjaHfkAUWdxFJBEYDZwHNAZ6isjBnaMLgQxVbQa8DgwPd1ATQGPG/H5Bz/jXvj0sW+buFDYxL5SWe0tgpaquUtU9wGTggvwHqOocVd2ZuzkPqBvemCZwdu1ysxJeeCHUqeM7jQG45BKoVs1NuWxiXijFvQ6Qf7n0dbmfO5S/Af8pTSgTB15/3d08c911vpOYPJUqudk4p0yBn3/2ncaUUijFvaB7wQvslBOR3kAGMOIQ+/uJyAIRWbBhw4bQU5rg6d7dFfjWrX0nMfldcw3s2wcTJ/pOYkoplOK+Djgq33Zd4KeDDxKR9sCdQBdV3V3QiVR1nKpmqGpGzZo1S5LXBEVyspusyuaRiS7HHgtvvw033eQ7iSmlUIr7Z0BDEakvIuWBHsD0/AeIyEnA07jCbmOpTOHuvhtGj/adwhzK2We7RVNMTCuyuKvqPmAAMAtYDryqqktFZKiIdMk9bASQArwmIotEZPohTmfiXWYmPPyw3Q0Z7caPZ/hMm+89loV0h6qqzgRmHvS5wfk+bh/mXCaoJkxw08tee63vJKYw69fzf1++C19/Dccf7zuNKQG7Q9WUnexsN7b9rLPgxBN9pzGFufJK9iSUg7FjfScxJWRzy5iy89//urtShw3znaRQ5cslUG/QjLCft05qRT4a1Dbs542IWrWYefzpXDhxItx/v03FHIOsuJuyU62au1Gma1ffSQq1Z18Oq4d1DPt5I/EHI5L+dXJnLlz2P3j+eetGi0FW3E3ZOe009zAxYWGd42HgQJuKOUZZcTdlov2KT+HHdJtqINY89pjvBKaE7IKqibxNm3j8zeEweHDRx5ros24dPPus7xSmmKy4m8h79lkq7d0N11/vO4kpiRdegKuusnsTYowVdxNZe/fCE0/w0dHNoFkz32lMSfTr56aLsC6amGLF3UTW1Kmwbh0TMi4o+lgTnapXh8sucy34jRt9pzEhsuJuImv5cjjhBGY3aOE7iSmNgQPdncXPPOM7iQmRFXcTWUOGwMKFqNiPWkxr0gTOOw9++cV3EhMiGwppIuenn+DII22GwaB4801ITPSdwoTImlMmMtauhaOPhqef9p3EhEteYf/2W1tEOwZYcTeRMWqUKwAdOvhOYsJp5kxo1AjmzvWdxBTBirsJv02bYNw46NHDtd5NcLRtC4cfHvWTvxkr7iYSnnoKduyAW2/1ncSEW3Iy3HCDW4pv4ULfaUwhrLib8MrJcQtynHee3bQUVP37Q5UqMHy47ySmEDZaxoRXQgIsWABbtvhOYiIlNdUV+Kefhq1boWpV34lMAazlbsInJ8ddRD3sMDjmGN9pTCQNGuQWXrHCHrWsuJvwefVVaNHCjW83wXbYYe6hCnv2+E5jCmDF3YSHquuD3b4datf2ncaUhd274ZRT4L77fCcxBbDibsLjv/91oyf+8Q/X726Cr0IFt/jKk0/Cb7/5TmMOYr+FpvRU4d573Zj2Pn18pzFl6a673MXzJ57wncQcxIq7Kb1334VPP4U77oDy5X2nMWWpeXPo2BFGjoRt23ynMflYcTeld9ZZbmx7376+kxgfhgxxdyXbUnxRxca5m9KrUAEuv9x3CuNLixYwbRqce67vJCYfa7mbklN188e89JLvJMa3Ll3cH3mbLTJqWHE3JTdnDrzyit2NapyZMyEjww2HNd5Zt0wZOH3YbH7csius56yTWpGPBrUN6zmhGFlVee2l2zgq5TD+8v2R7B40I+xZgqZ8uQTqBfl1Ouww+OILGD0abrutwEMi9RpE6vchlllxLwM/btnF6mEdw3rOSBWJkLPOmAHDl8HYsXxzddciDw90UQvRnn05Yf85gCh6bVu1gvPPh4cegquvdnPQHCTwr0EUsW4ZU3w5OW7Y47HHwhVX+E5josk//wmbN7sCb7yy4m6KTwQefNDdmZiU5DuNiSZpadCrFzz2GKxf7ztNXLNuGVN8Iu7ttzEFuf9+N4qqVi3fSeKaFXdTPOPHw8qVMHSotdpNwerXdw/jlXXLmNBt2wZ33gmffALlrF1givDgg3ZNxiMr7iZ0Dz7o+lGHD3ddM8YUZtcueO45mDfPd5K4FFJxF5EOIvKNiKwUkUEF7D9LRL4QkX0i0i38MY1333/vJofq0wdatvSdxsSCW2+FI46A6693I6xMmSqyuItIIjAaOA9oDPQUkcYHHbYG6AtMCndAEyUGDYLERDfUzZhQpKS4d3vz59sUFR6E0nJvCaxU1VWqugeYDFyQ/wBVXa2qSwD78xxUd9/tZv2rW9d3EhNL+vRxE4vdcQfs3es7TVwJ5apYHWBtvu11wCmRiWOi1oknuocxxZGQAM884yYUs9FVZSqUlntBV85KNPWbiPQTkQUismDDhg0lOYUpa88848Ys22RQpqTS0iA9HYCkbGu9l5VQivs64Kh823WBEi1vr6rjVDVDVTNq1qxZklOYsrR+vbsotn49VK7sO42JdQMG8PSUB2xa4DISSnH/DGgoIvVFpDzQA5ge2VgmKtxyC+zYAWPG2NBHU3rHHkvbVQvg9dd9J4kLRRZ3Vd0HDABmAcuBV1V1qYgMFZEuACLSQkTWAd2Bp0VkaSRDmzLw3nvw4otulMzxx/tOY4Lg739nSe1jYeBA2LrVd5rAC2mcu6rOVNXjVLWBqj6Q+7nBqjo99+PPVLWuqlZW1eqq2iSSoU0ZuOsuaNAAbr/ddxITFImJ3HHuAPj1V/u5KgN2h6op2LRp8NprULGi7yQmQL7Ka7lPmgSZmb7jBJoVd3OAmts3QXa2m9HvpJN8xzFBdN99sHQpVK/uO0mgWXE3v9u9mxdfucvNx21MpKSkQJ06btTMJ5/4ThNYVtzN7+65h0Yb10Dfvr6TmHjw5JNw+unw/vu+kwSSFXfjzJsHw4czudk5cN55vtOYeHD55XDMMa4xsW2b7zSBY8XduGFpvXpB3bo80PZvvtOYeJGSAhMnwurVcMMNvtMEjhV3Az//7Ob9ePlltlWwO1FNGTrjDDep2IQJNnNkmNlyOsbdpLR0qVtdafoM32lMvLnnHvj2WzjySN9JAsVa7vFs6VK48UbIyrJl84w/5crBq69CmzZu2+aeCQsr7vFq82a48EJ4+WXYssV3GmNcUR88GK6+2gp8GFhxj0fZ2e4C6g8/wBtvQO3avhMZ4yany85200yPHu07Tcyz4h6P7rwT/vvf38cZGxMt7rsPOnd2o2dmz/adJqZZcY8369fDU09B//7Qr5/vNMYcKCHBzUbaqBF07w6rVvlOFLOsuMebww93Nyw99pjvJMYU7E9/gunT3fDcRYt8p4lZNkQiXixeDB9+CNddB40b+05jTOEaNIAVK6BKFbetagvGFJO13OPB6tVuSoFhw2yRBBM78gr7G2/ARRfBXlt/tTisuAfdDz+48cO7drmLqFWr+k5kTPFs3gz//jf07g379vlOEzOsWybI1q51hX3zZrdsXhNbIMvEoCuvdO84b7nFbb/0kt10FwJ7hYLsf/9zhf2dd6B5c99pjCm5m292/95yi+tXEQ0zAAAKz0lEQVR/nzTJCnwR7NUJoqwsSE52b2M7dIAaNXwnMqb0br7ZXVRduxYSE32niXrW5x40H3/sRhp8+KHbtsJuguSmm2DUKFfkly+HjRt9J4paVtyD5PXXoX17qFzZLWNmTFDt2+fuZD3tNFi50neaqGTFPQhyctyUAt27Q3o6fPQR1K/vO5UxkVOuHDz/PGRmQosWtP5uge9EUceKexC8+CL8859uVMGcOVCzpu9ExkTeaafBggVw9NFMeP1euP9+19AxgF1QjW2//eZu1e7dG6pVg06d7C4+E1/q14ePP2Zai450/eADmyo4H2u5x6I9e/jH//4Fxx3nlshLSHD9j1bYTTyqVIkbO93sbnRKTHSjaaZO9Z3KOyvusebDD+Gkk7hu3muupZ53i7Yx8UwEKlZ0H48Y4aYr6NHDNX7ilBX3WJGT46bpPfNM2L6dy7sNgWefdSvIG2N+98gjcO+9riV//PFu3YLsbN+pypwV92iX90OZkACbNrlxvkuXMqdBC7+5jIlWSUluub4vv4RTToG//x2GD/edqszZBdVolZUFEye6mRxnzXKLF0ye7Iq8MaZoDRu6350pU6B1a/e5efNg504351LAr1FZpYg2W7e6gl6vHlxzDRxxhCv0YIXdmOISgYsvhurV3fb990O7dtCiBbzySqBnmbRqEU327HGtjdtvdzcjzZnjphNIS/OdzJhgeP11GDcOtm1zF1yPPRbGj/edKiKsuPu0ciU89JAbpw5QvrzrG/ziCzf3euvWgX/raEyZSk6Gq66CZctcd03Dhm6tA4Dt2910wgFZ0Mb63MvaypWu7/yNN35fH/KUU2DLFkhNhb59vcYzJi4kJkLXru6Rd+PT1Klw2WXugmz79m7fOefA0Uf7zVpC1nKPpH37YMkSLl04062IBO6Czt13Q6VKMHKkWwJv3jxX2I0xZS/v3fGll7pu0Ouvh6+/hn793LWvNWvc/hUr3O9xjNwFG1LLXUQ6AI8BicCzqjrsoP0VgOeB5kAmcImqrg5v1CiXt4Dv+vWuq2X+fFi4EHbu5AGAt9Pd28ELL4RffoHDD/ed2BiTX0ICnHqqewwf7rpuPvkE/vxnt3/wYPeu+/DDoWVLd1H2pJPczYRRqMjiLiKJwGjgbGAd8JmITFfVZfkO+xuwWVWPFZEewEPAJZEI7E1e8c7Odv/Ba9e6x4oV7q/85Ze7GycqVICnn3YXQa+6CmbNYuHWHE668kp3npSU4t14lDeEa+7ccH9H4ZE/X1FZC9sf6r68dzhbthT+HHldXunpBZ83NdX1sZ5xRuh58j6XpzT/JwefvxT/z5MnDYJ5I0J/3Qt7DUuaqzjnDKfWrWHRIib/6SgY1jF85xVxy1LmX5ry9tvdTYTz58Onn8Kbb7qfr7zi3revu1B79NHuD8LRR7ubqE44IXy5iiGUlntLYKWqrgIQkcnABUD+4n4BcE/ux68DT4qIqHp+/6LqivHeve7jSpXc55ctczcE7djhfsF37HCLWpx/vtt/ww3w/fduIYDMTPdv587w3HPur/vVV7uvSU11V9tbt/69iKSmugsyeUuAtW7N7h2ZdmHUmFjXrJl7XHut296+HTZsOPCYZcvcYIidO912584wfbr7+OSTXR3q3BmGDo143FCKex1gbb7tdcAphzpGVfeJyFagOhCZZVIGDoT333dFe88e92/9+m7oIMDZZ/++P+/vyymnuL5tgJ49YcmSA8/Ztu3vxX3hQlegq1d3LfDq1d3XgyvSixe7t2aHaoHb2o7GBN/B78InTnT/qrpG4Zo1B9aCFi3gp5/cqLgyEEoVKqjJeXCLPJRjEJF+QD+AP+f1Y5VEjRruQkdS0u+PunV/33/RRZCRceD+I4/8ff8TT7g/CpUru0dKClSt+vv+//2v8Odv0KBYccuXS6TeoBnF+po8k1dlAtDjoK+vUqFcic9ZmGqVkop13vz5DpUVoGJSIvMK2V/Y1+bftyRrLwDNco+rmHTga5t3bOPc45Yd4rxLsvZSKUdZsCrzD/sOlTXv3HkKylqY/FkP/n4L+/6L8qoI8wr4PvIcfO6DX8NDHVvU/1l+xTlnOH9uJ6/KpHHWXsofFkVrqoq4GnXwEpdPP122MYrqORGRU4F7VPXc3O3bAVT1wXzHzMo95hMRKQf8AtQsrFsmIyNDFyyIg9VTStvfaH3u1ude3HMVtT+Afe6kp0fv70iYicjnqppR1HGhDIX8DGgoIvVFpDzQA5h+0DHTgb/mftwNmO29v90YY+JYkS13ABE5H3gUNxRygqo+ICJDgQWqOl1EkoEXgJOATUCPvAuwhxI3LXdjjAmjUFvuIV35U9WZwMyDPjc438dZQPfihjTGGBMZdoeqMcYEkBV3Y4wJICvuxhgTQFbcjTEmgKy4G2NMAFlxN8aYALLibowxAWTF3RhjAiikO1Qj8sQiG4AfDrG7BpGaUTL8LGvkxFLeWMoKsZU3lrJC5PMerao1izrIW3EvjIgsCOX22mhgWSMnlvLGUlaIrbyxlBWiJ691yxhjTABZcTfGmACK1uI+zneAYrCskRNLeWMpK8RW3ljKClGSNyr73I0xxpROtLbcjTHGlEJUF3cRuUVEVERqFH20PyJyn4gsEZFFIvK2iBxZ9Ff5ISIjROTr3LxTRSTVd6bCiEh3EVkqIjki4n0EQkFEpIOIfCMiK0VkkO88hRGRCSLyq4h85TtLUUTkKBGZIyLLc38Grved6VBEJFlE5ovI4tys9/rOFLXFXUSOAs4G1vjOEoIRqtpMVdOBt4DBRX2BR+8AJ6pqM+Bb4HbPeYryFXAR8L7vIAURkURgNHAe0BjoKSKN/aYq1ESgg+8QIdoH3KyqJwCtgOui+LXdDbRV1TQgHeggIq18Bora4g6MAm4Fov6igKr+lm+zMlGcWVXfVtV9uZvzgLo+8xRFVZer6je+cxSiJbBSVVep6h5gMnCB50yHpKrv45bCjHqq+rOqfpH78TZgOVDHb6qCqbM9dzMp9+G1DkRlcReRLsCPqrrYd5ZQicgDIrIWuJTobrnndwXwH98hYlwdYG2+7XVEaQGKZSJSD7dG86d+kxyaiCSKyCLgV+AdVfWaNaQ1VCNBRN4Fahew607gDuCcsk1UuMLyquo0Vb0TuFNEbgcGAEPKNGA+RWXNPeZO3Nvel8oyW0FCyRvFpIDPRe07t1gkIinAG8ANB71Ljiqqmg2k517HmioiJ6qqt2sb3oq7qrYv6PMi0hSoDywWEXDdBl+ISEtV/aUMIx7gUHkLMAmYgcfiXlRWEfkr0Alop1EwFrYYr200WgcclW+7LvCTpyyBIyJJuML+kqpO8Z0nFKq6RUTm4q5teCvuUdcto6pfqmotVa2nqvVwvzwn+yzsRRGRhvk2uwBf+8pSFBHpANwGdFHVnb7zBMBnQEMRqS8i5YEewHTPmQJBXOtuPLBcVUf6zlMYEamZN/JMRCoC7fFcB6KuuMeoYSLylYgswXUnRe2QLeBJoArwTu7QzbG+AxVGRLqKyDrgVGCGiMzynSm/3IvTA4BZuAt+r6rqUr+pDk1EXgY+ARqJyDoR+ZvvTIU4HegDtM39WV0kIuf7DnUIRwBzcmvAZ7g+97d8BrI7VI0xJoCs5W6MMQFkxd0YYwLIirsxxgSQFXdjjAkgK+7GGBNAVtyNMSaArLgbY0wAWXE3xpgA+n91tdS3Ku/MUQAAAABJRU5ErkJggg==\n", "text/plain": [ "