{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Estimation of rates of random events\n",
"\n",
"*Antonino Ingargiola* -- tritemio AT gmail.com
\n",
"**ORCID** [0000-0002-9348-1397](orcid.org/0000-0002-9348-1397)
\n",
"**twitter** [@tritemio_sc](https://twitter.com/tritemio_sc)
\n",
"\n",
"Date: April 2016\n",
"\n",
"\n",
"## Abstract\n",
"\n",
"
In this notebook I explore properties of different \n", "estimators of \"rate\" for random events generated by a Poisson process.\n", "The initial reason was computing the rate of photon detection\n", "events (for example in single-molecule spectroscopy experiments),\n", "but the treatment is general and valid for any Poisson process.
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# The Model\n", "\n", "Let's assume that $\\{d_i\\}$ is a set of $n$ samples distributed as $D \\sim\\operatorname{Exp}(\\lambda)$.\n", "\n", "Then, the can construct $m = n + 1$ times of random events $\\{t_i\\}$ as:\n", "\n", "$$t_0 = 0, \\quad t_i = \\sum_{j=0}^{i-1} d_j \\quad\\textrm{for}\\quad i > 1$$\n", "\n", "These events are stochastic, statistically independent and have a fixed rate $\\lambda$.\n", "\n", "# Problem: how to estimate $\\lambda$ from $\\{t_i\\}$?\n", "\n", "Given $\\{t_i\\}$ we can trivially compute $\\{d_i\\}$ as:\n", "\n", "$$d_i = t_{i+1} - t_i$$\n", "\n", "Instead of $\\lambda$, we can directly estimate $\\tau = 1/\\lambda$:\n", "\n", "$$\\hat{\\tau} = \\frac{1}{n}\\sum_{i=0}^{n-1} d_i = \\frac{T_n}{n}$$\n", "\n", "The estimator $\\hat{\\tau} \\sim \\frac{1}{n}\\operatorname{Erlang}(n, \\lambda)$.\n", "The Erlang distribution has mean $n/\\lambda$ and variance $n/\\lambda^2$.\n", "\n", "Note that:\n", "\n", "$$\\hat{\\tau} \\sim \\frac{1}{n}\\operatorname{Erlang}(n, \\lambda) = \\operatorname{Erlang}(n, n\\,\\lambda)$$\n", "\n", "thus, $\\hat{\\tau}$ has mean $1/\\lambda$ and variance $1/(n \\lambda)^2$, i.e.\n", "it is an unbiased estimator of $1/\\lambda$.\n", "\n", "However, the fact the $\\hat{\\tau}$ is an unbiased estimator of $1/\\lambda$ \n", "does not mean that $1/\\hat{\\tau}$ is an unbiased estimator of $\\lambda$.\n", "Indeed, as shown in the next section, the unbiased estimator of the rate is:\n", "\n", "$$\\hat\\lambda_u = \\frac{n-1}{T_n}$$ \n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It can be shown that $\\hat{\\lambda}_u$ has also the minimum variance among the\n", "unbiased estimators (i.e. it is a minimum variance unbiased estimator, MVUE)\n", "of the rate (see [Lehmann-Scheffé theorem](https://en.wikipedia.org/wiki/Lehmann%E2%80%93Scheff%C3%A9_theorem)\n", "and [Rao–Blackwell theorem](https://en.wikipedia.org/wiki/Rao%E2%80%93Blackwell_theorem))." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## MLE estimator of the rate\n", "\n", "For our sample of $n$ observations from an Exponential distribution, \n", "the likelihood function is\n", "\n", "$$L(\\lambda) = \\prod_{i=1}^n f(d_i;\\lambda) = \\prod_{i=1}^n \\lambda\\, e^{-\\lambda d_i}$$\n", "\n", "$$\\log L(\\lambda) = \\sum_{i=1}^n \\log (\\lambda\\, e^{-\\lambda d_i} ) = \n", "n\\log(\\lambda) - \\lambda \\sum_{i=1}^n d_i = \n", "n\\log(\\lambda) - \\lambda T_n$$\n", "\n", "Taking the first derivative and setting it to 0 we obtain:\n", "\n", "$$\\frac{\\rm d}{{\\rm d}\\lambda} \\log L(\\lambda) = \n", "\\frac{n}{\\lambda} - T_n = 0$$\n", "\n", "$$\\Rightarrow\\quad \\hat{\\lambda}_{MLE} = \\frac{n}{T_n} = \\frac{n}{\\sum_{i=1}^n d_i}$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Minimum MSE estimator\n", "\n", "Given an estimator $\\Lambda_n$ (a R.V.), we define the mean square error (MSE) as:\n", "\n", "$$\\epsilon(\\Lambda_n) = \\mathrm{E}\\{ (\\Lambda_n - \\lambda)^2 \\} = \n", "\\mathrm{Var}\\{ \\Lambda_n\\} + (\\mathrm{E}\\{\\Lambda_n\\} - \\lambda)^2$$\n", "\n", "In general the minimum MSE (MMSE) estimator is defined as:\n", "\n", "$$\\hat{\\lambda}_{MSE} = \\operatorname*{arg\\,min}_{\\Lambda_n} \\epsilon(\\Lambda_n)$$\n", "\n", "If $\\Lambda_n$ has the form:\n", "\n", "$$\\Lambda_{n,c} = \\frac{n - c}{T_n}$$\n", "\n", "we need to find $c$ that minimizes $\\epsilon(\\Lambda_{n,c})$.\n", "It can be shown that the minimum is achieved for $c=2$\n", "\n", "To derive this result use the expression for $\\mathrm{Var}\\{ \\Lambda_{n,c}\\}$\n", "derived in [Appendix 2](#Appendix-2.-Variance-of-the-rate-estimators)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Family of estimators\n", "\n", "We consider the family of estimators \n", "\n", "$$\\Lambda_{n,c} = \\frac{n - c}{T_n}$$\n", "\n", "where $T_n = \\sum d_i \\sim \\textrm{Erlang}(\\lambda, n)$ and $c < n$. We can write:\n", "\n", "$$\\textrm{Erlang PDF:}\\quad f(x \\:|\\: n, \\lambda) = \n", "\\frac{\\lambda^n x^{n-1} e^{-\\lambda x}}{(n-1)!\\,}$$\n", "\n", "$$\\mathrm{E}\\{\\Lambda_{n,c}\\} = \\int \\frac{n - c}{x} \\, f(x \\:|\\: n, \\lambda)\\, dx\n", "= (n - c) \\int_0^\\infty \\frac{1}{x}\\,f(x \\:|\\: n, \\lambda)\\, dx\n", "= (n - c) \\int_0^\\infty \\frac{\\lambda^n x^{n-2} e^{-\\lambda x}}{(n-1)!\\,}\\, dx \n", "= \\frac{(n - c)\\lambda^n}{(n-1)!} \\int_0^\\infty x^{n-2} e^{-\\lambda x}\\, dx $$\n", "\n", "Knowing that (see [Appendix 1](#Appendix-1.-Computation-of-integral-\"A\")):\n", "\n", "$$ \\int_0^{+\\infty} x^{n-2} e^{-\\lambda x}\\, dx = \\frac{(n-2)!}{\\lambda^{n-1}}\\qquad\\qquad \\textrm{(A)}$$\n", "\n", "we find that:\n", "\n", "$$\\mathrm{E}\\{\\Lambda_{n,c}\\} = (n - c)\\frac{\\lambda}{n-1} = \\overline{\\Lambda}_{n,c}$$\n", "\n", "Note that \n", "\n", "- for $c=1$ we obtain the unbiased estimator ($\\overline{\\Lambda}_{n,c} = \\lambda$). \n", "- for $c=0$ we obtain the MLE estimator derived in a previous section. \n", "- for $c=2$ we obtain the minimum MSE (MMSE) estimator introduced in a previous section. \n", "\n", "Furthermore, from simulations (see below), I found that for $c=1/3$ the estimator\n", "has median equal to the rate $\\lambda$. This result should not be difficult\n", "to prove analytically." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Simulation of rate estimators" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": true }, "outputs": [], "source": [ "%matplotlib inline\n", "import matplotlib.pylab as plt\n", "import numpy as np" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": true }, "outputs": [], "source": [ "n = 5 # number of delays\n", "num_times = n + 1 # number of timestamps\n", "λ = 1000 # simulated rate \n", "num_iter = 100 * 1000" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [], "source": [ "np.random.seed(5)\n", "d1 = np.random.exponential(scale=1/λ, size=(num_iter * n)).reshape(num_iter, n)\n", "t = np.cumsum(d1, axis=1)\n", "\n", "assert np.allclose(np.cumsum(d1, axis=1), t)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We $n$ delays, exponentially distributed with a fixed rate. Their cumulative sum (the times) will not have a sharp cut-off: " ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAEACAYAAABLfPrqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGSFJREFUeJzt3W2sXdWd3/HvD1JCZnKhMFPskQ0JETjBaSVjaZwZpS9u\nRMvDVAUUCeTMqEDjqJGAJmqkNjhSZXs00iRISZ1oZF4kzGBQIo9LNANREBCEzkiRILgFxiR24faF\nCTa5lwgHxjQzEQ//vjj7wrZ97sO+PvfJ9/uRtljnv/c6Z+3N5v5Ze+29dqoKSZK6OGOxGyBJWn5M\nHpKkzkwekqTOTB6SpM5MHpKkzkwekqTOZp08kpyR5OkkDzafz0vyaJLnkzyS5NzWtluTjCU5mOTK\nVnxjkv1JXkiysxU/K8meps4TSS4a1g5KkoavS8/ji8CB1uc7gMeq6qPA48BWgCTrgRuBy4BrgF1J\n0tS5C9hSVeuAdUmuauJbgKNVdSmwE7hzjvsjSVoAs0oeSdYCfwR8pxW+DtjdlHcD1zfla4E9VfVW\nVR0CxoBNSVYDI1W1r9nu3lad9nfdD1zRfVckSQtltj2P/wH8V6D9OPqqqpoAqKpx4IImvgZ4qbXd\nkSa2Bjjcih9uYsfVqaq3gdeSnD/73ZAkLaQZk0eSfwdMVNWzQKbZdJjznEz3O5KkRfa+WWzzSeDa\nJH8EfAAYSXIfMJ5kVVVNNJekXmm2PwJc2Kq/tolNFW/XeTnJmcA5VXX0xIYkcSIuSZqDqhrq/5TP\n2POoqq9U1UVV9RFgM/B4Vf0H4AfALc1mNwMPNOUHgc3NHVQXA5cATzWXtl5PsqkZQL/phDo3N+Ub\n6A/AT9UelyEt27ZtW/Q2nC6Lx9LjuZSX+TCbnsdUvgrsTfJZ4EX6d1hRVQeS7KV/Z9abwK31Xutv\nA+4BzgYeqqqHm/jdwH1JxoBX6ScpSdIS1Sl5VNXfAX/XlI8C/2aK7f4c+PMB8f8N/KsB8d/QJJ+Z\nfOYznzspdt555/D973+fV175+UnrVq36EOPjh2bz1ZKkWTqVnsei2LPnD06KnX32f+ef/mmcQWP2\nExOOvU9ldHR0sZtw2vBYDpfHc+nLfF0Pmw/9AfOT2zsysp5jxw4y+IavzNs1P0laDpJQCz1gvvy9\nnyQDl9WrP7zYjZOkZWkFJI/f0O+RnLxMTIybVCRpDpbdmMdwTSaW4zlOIknTWwE9D0nSsJk8JEmd\nmTwkSZ2ZPAbyDi1Jms4KHzCfyuCBdHAwXZLAnockaQ5MHpKkzkwekqTOTB6SpM5MHp15J5YkebdV\nZ96JJUn2PCRJnZk8JEmdmTwkSZ3NmDySvD/JT5I8k+S5JNua+LYkh5M83SxXt+psTTKW5GCSK1vx\njUn2J3khyc5W/Kwke5o6TyS5aNg7KkkanhmTR1X9BvhUVV0ObACuSbKpWf2NqtrYLA8DJLkMuBG4\nDLgG2JVkciT5LmBLVa0D1iW5qolvAY5W1aXATuDOIe2fJGkezOqyVVX9uim+n/4dWpO3Gw26veg6\nYE9VvVVVh4AxYFOS1cBIVe1rtrsXuL5VZ3dTvh+4ostOSJIW1qySR5IzkjwDjAM/aiWA25M8m+Q7\nSc5tYmuAl1rVjzSxNcDhVvxwEzuuTlW9DbyW5Py57JAkaf7NtufxTnPZai39XsR6YBfwkaraQD+p\nfH2I7fKBCUlawjo9JFhV/5CkB1xdVd9orfo28IOmfAS4sLVubRObKt6u83KSM4Fzquro4FZsb5VH\nm2Wp6D99fqJVqz7E+PihhW+OpBWp1+vR6/Xm9TdSNfhp6Xc3SH4XeLOqXk/yAeAR4KvA01U13mzz\nX4Dfr6o/bnol3wU+Qf9y1I+AS6uqkjwJfAHYB/wQ+FZVPZzkVuBfVtWtSTYD11fV5gFtqUFPd4+M\nrOfYsYMMfvI7U8SnWzeXOtN/30zHWZLmSxKqaqhXdGbT8/g9YHeSM+hf5vrrqnooyb1JNgDvAIeA\nzwNU1YEke4EDwJvArfXeX87bgHuAs4GHJu/QAu4G7ksyBrwKnJQ4JElLx4w9j6XEnockdTcfPQ+f\nMJckdWbykCR1ZvKQJHVm8pAkdWbykCR1ZvKQJHVm8pAkdWbykCR1ZvKQJHVm8pAkdWbykCR1ZvJY\nEP2p2gctq1d/eLEbJ0mddXqfh+bqN0w1meLEhO+9krT82POQJHVm8pAkdWbykCR1ZvKQJHVm8pAk\ndWbykCR1NmPySPL+JD9J8kyS55Jsa+LnJXk0yfNJHklybqvO1iRjSQ4mubIV35hkf5IXkuxsxc9K\nsqep80SSi4a9o5Kk4ZkxeVTVb4BPVdXlwAbgmiSbgDuAx6rqo8DjwFaAJOuBG4HLgGuAXUkmH2a4\nC9hSVeuAdUmuauJbgKNVdSmwE7hzWDsoSRq+WV22qqpfN8X303+wsIDrgN1NfDdwfVO+FthTVW9V\n1SFgDNiUZDUwUlX7mu3ubdVpf9f9wBVz2htJ0oKYVfJIckaSZ4Bx4EdNAlhVVRMAVTUOXNBsvgZ4\nqVX9SBNbAxxuxQ83sePqVNXbwGtJzp/THkmS5t2spiepqneAy5OcA/xNko9z8nwbg+ffmJtp5uzY\n3iqPNoskaVKv16PX683rb3Sa26qq/iFJD7gamEiyqqommktSrzSbHQEubFVb28SmirfrvJzkTOCc\nqjo6uBXbuzRZklac0dFRRkdH3/28Y8eOof/GbO62+t3JO6mSfAD4t8BB4EHglmazm4EHmvKDwObm\nDqqLgUuAp5pLW68n2dQMoN90Qp2bm/IN9AfgV4jBM+46266kpWw2PY/fA3YnOYN+svnrqnooyZPA\n3iSfBV6kf4cVVXUgyV7gAPAmcGtVTV7Sug24BzgbeKiqHm7idwP3JRkDXgU2D2XvloXBM+46266k\npSzv/V1f+pLUoD+0IyPrOXbsIIOHXTJFfLp1c6kz7O8Ly+nfjaSlKwlVNdT/I/UJc0lSZyYPSVJn\nJg9JUmcmD0lSZyYPSVJnJg9JUmcmD0lSZyYPSVJnJg9JUmcmD0lSZyYPSVJnJg9JUmcmD0lSZyYP\nSVJnJo8la/BLonxRlKSloNNraLWQBr8kCnxRlKTFZ89DktSZyUOS1JnJQ5LU2YzJI8naJI8n+VmS\n55L85ya+LcnhJE83y9WtOluTjCU5mOTKVnxjkv1JXkiysxU/K8meps4TSS4a9o5KkoZnNj2Pt4Av\nVdXHgT8Ebk/ysWbdN6pqY7M8DJDkMuBG4DLgGmBXkskR3ruALVW1DliX5KomvgU4WlWXAjuBO4ex\nc5Kk+TFj8qiq8ap6tim/ARwE1jSrB932cx2wp6reqqpDwBiwKclqYKSq9jXb3Qtc36qzuynfD1wx\nh32RJC2QTmMeST4MbAB+0oRuT/Jsku8kObeJrQFealU70sTWAIdb8cO8l4TerVNVbwOvJTm/S9sk\nSQtn1s95JPkg/V7BF6vqjSS7gD+tqkryZ8DXgc8NqV3TPMiwvVUebRZJ0qRer0ev15vX35hV8kjy\nPvqJ476qegCgqn7Z2uTbwA+a8hHgwta6tU1sqni7zstJzgTOqaqjg1uzfTZNlqQVa3R0lNHR0Xc/\n79ixY+i/MdvLVn8JHKiqb04GmjGMSZ8GftqUHwQ2N3dQXQxcAjxVVePA60k2NQPoNwEPtOrc3JRv\nAB6f095IkhbEjD2PJJ8E/gR4Lskz9OfM+Arwx0k2AO8Ah4DPA1TVgSR7gQPAm8CtVTU5z8ZtwD3A\n2cBDk3doAXcD9yUZA14FNg9l7yRJ8yLv/V1f+pLUoPmeRkbWc+zYQQbPBZUp4tOtm0udYX/f9HWW\n0783SYsrCVU11EnxfMJcktSZyUOS1JnJQ5LUmcljWfJFUZIWly+DWpZ8UZSkxWXPQ5LUmclDktSZ\nyUOS1JnJQ5LUmclDktSZyUOS1JnJQ5LUmclDktSZyUOS1JnJQ5LUmclDktSZyUOS1JnJQ5LUmcnj\ntDN4unanapc0TDMmjyRrkzye5GdJnkvyhSZ+XpJHkzyf5JEk57bqbE0yluRgkitb8Y1J9id5IcnO\nVvysJHuaOk8kuWjYO7pyTE7XfvwyMfHiorZK0ullNj2Pt4AvVdXHgT8EbkvyMeAO4LGq+ijwOLAV\nIMl64EbgMuAaYFeSyZdM3AVsqap1wLokVzXxLcDRqroU2AncOZS9kyTNixmTR1WNV9WzTfkN4CCw\nFrgO2N1sthu4vilfC+ypqreq6hAwBmxKshoYqap9zXb3tuq0v+t+4IpT2SlJ0vzqNOaR5MPABuBJ\nYFVVTUA/wQAXNJutAV5qVTvSxNYAh1vxw03suDpV9TbwWpLzu7RNkrRwZv0a2iQfpN8r+GJVvZHk\nxPegDn4v6txM8y7V7a3yaLNIkib1ej16vd68/saskkeS99FPHPdV1QNNeCLJqqqaaC5JvdLEjwAX\ntqqvbWJTxdt1Xk5yJnBOVR0d3Jrts2myJK1Yo6OjjI6Ovvt5x44dQ/+N2V62+kvgQFV9sxV7ELil\nKd8MPNCKb27uoLoYuAR4qrm09XqSTc0A+k0n1Lm5Kd9AfwBekrREzdjzSPJJ4E+A55I8Q//y1FeA\nrwF7k3wWeJH+HVZU1YEke4EDwJvArVU1eUnrNuAe4Gzgoap6uInfDdyXZAx4Fdg8nN2TJM2HvPd3\nfenrj7Oc3N6RkfUcO3aQwcMumSI+3bq51Bn29w2/Dcvp37Wk4UlCVU0zltydT5hLkjozeUiSOjN5\nSJI6M3lIkjozeUiSOjN5SJI6M3lIkjozeawYg18S5YuiJM3FrCdG1HI3+ZKok01MDPXZIUkrgD0P\nSVJnJg9JUmcmD0lSZyYPSVJnJg9JUmcmD0lSZyYPSVJnJg9JUmcmD0lSZyYPSVJnMyaPJHcnmUiy\nvxXbluRwkqeb5erWuq1JxpIcTHJlK74xyf4kLyTZ2YqflWRPU+eJJBcNcwclScM3m57HXwFXDYh/\no6o2NsvDAEkuA24ELgOuAXYlmZw46S5gS1WtA9YlmfzOLcDRqroU2AncOffdkSQthBmTR1X9GPjV\ngFWDZtO7DthTVW9V1SFgDNiUZDUwUlX7mu3uBa5v1dndlO8Hrph98zUcg2fcdbZdSVM5lTGP25M8\nm+Q7Sc5tYmuAl1rbHGlia4DDrfjhJnZcnap6G3gtyfmn0C51Njnj7vHLxMSLi9oqSUvXXKdk3wX8\naVVVkj8Dvg58bkhtmmF+8O2t8mizSJIm9Xo9er3evP7GnJJHVf2y9fHbwA+a8hHgwta6tU1sqni7\nzstJzgTOqaqjU//69rk0WZJWjNHRUUZHR9/9vGPHjqH/xmwvW4VWj6AZw5j0aeCnTflBYHNzB9XF\nwCXAU1U1DryeZFMzgH4T8ECrzs1N+Qbg8TntiSRpwczY80jyPfrXhn4nyc+BbcCnkmwA3gEOAZ8H\nqKoDSfYCB4A3gVuravL1dbcB9wBnAw9N3qEF3A3cl2QMeBXYPJQ9kyTNm7z3t33pS1KDXqU6MrKe\nY8cOMvg1q5kiPt26udQZ9vctjTYsp/ND0mBJqKqhvm/aJ8wlSZ2ZPCRJnZk8JEmdmTwkSZ2ZPCRJ\nnZk8JEmdmTwkSZ2ZPDSNwbPtOuOupLlOjKgVYXK23ZNNTAz1eSNJy4w9D0lSZyYPSVJnJg9JUmcm\nD0lSZyYPSVJnJg9JUmcmD0lSZyYPSVJnJg9JUmcmD0lSZzMmjyR3J5lIsr8VOy/Jo0meT/JIknNb\n67YmGUtyMMmVrfjGJPuTvJBkZyt+VpI9TZ0nklw0zB2UJA3fbHoefwVcdULsDuCxqvoo8DiwFSDJ\neuBG4DLgGmBXkslJkO4CtlTVOmBdksnv3AIcrapLgZ3AnaewP1owgydNdMJEaWWYMXlU1Y+BX50Q\nvg7Y3ZR3A9c35WuBPVX1VlUdAsaATUlWAyNVta/Z7t5WnfZ33Q9cMYf90IKbnDTx+GVi4sVFbZWk\nhTHXMY8LqmoCoKrGgQua+BrgpdZ2R5rYGuBwK364iR1Xp6reBl5Lcv4c2yVJWgDDmpJ98LzdczPD\nXN/bW+XRZpEkTer1evR6vXn9jbkmj4kkq6pqorkk9UoTPwJc2NpubRObKt6u83KSM4Fzquro1D+9\nfY5NlqSVYXR0lNHR0Xc/79ixY+i/MdvLVuH4HsGDwC1N+WbggVZ8c3MH1cXAJcBTzaWt15NsagbQ\nbzqhzs1N+Qb6A/CSpCVsxp5Hku/Rvzb0O0l+DmwDvgr8zySfBV6kf4cVVXUgyV7gAPAmcGtVTV7S\nug24BzgbeKiqHm7idwP3JRkDXgU2D2fXJEnzJe/9bV/6ktSg4ZWRkfUcO3aQwUMvmSI+3bq51Bn2\n9y3XNoTldE5JK0ESqmqo7472CXNJUmcmDw3Z4IcHfYBQOr0M61ZdqTH58ODJJiaG2muWtIjseUiS\nOjN5SJI6M3lIkjozeUiSOjN5SJI6M3lIkjozeUiSOjN5aAH5AKF0uvAhQS0gHyCUThf2PCRJnZk8\nJEmdmTwkSZ2ZPCRJnZk8JEmdmTwkSZ2dUvJIcijJ3yd5JslTTey8JI8meT7JI0nObW2/NclYkoNJ\nrmzFNybZn+SFJDtPpU1argY/A+LzH9LSdKo9j3eA0aq6vKo2NbE7gMeq6qPA48BWgCTrgRuBy4Br\ngF1JJm/uvwvYUlXrgHVJrjrFdmnZmXwG5PhlYuLFRW2VpMFONXlkwHdcB+xuyruB65vytcCeqnqr\nqg4BY8CmJKuBkara12x3b6uOJGkJOtXkUcCPkuxL8rkmtqqqJgCqahy4oImvAV5q1T3SxNYAh1vx\nw01MkrREner0JJ+sql8k+RfAo0me5+T5JwbPRyFJWrZOKXlU1S+af/4yyd8Cm4CJJKuqaqK5JPVK\ns/kR4MJW9bVNbKr4FLa3yqPNIkma1Ov16PV68/obqZpbxyDJbwFnVNUbSX4beBTYAVwBHK2qryX5\nMnBeVd3RDJh/F/gE/ctSPwIurapK8iTwBWAf8EPgW1X18IDfrEEdmZGR9Rw7dpDBnZxMEZ9u3Vzq\nDPv7lmsbht3us+kPpp9s1aoPMT5+aIrvkzQpCVU11NlHT6XnsQr4m/4fdN4HfLeqHk3yv4C9ST4L\nvEj/Diuq6kCSvcAB4E3g1novc90G3EP/L8VDgxKHVipn4pWWojn3PBaDPY/l0IaFbfdyOn+lxTIf\nPQ+fMJckdWbykCR1ZvLQMuaUJtJi8TW0WsYGD6Y7kC7NP3sekqTOTB6SpM5MHjoNDR4LcTxEGh7H\nPHQa8sFCab7Z85AkdWbykCR1ZvKQJHVm8tAK42C6NAwOmGuFcTBdGgZ7HpKkzkwe0rucK0uaLS9b\nSe9yrixptux5SDNykF06kclDmtFkj+TkZWJi3KSiFcnLVtIp8VKXVqYl0/NIcnWS/5PkhSRfXuz2\nSKfGS106vS2J5JHkDOAvgKuAjwOfSfKxxW3VStBb7AacRnonfO5+qcvE8p5er7fYTdAMlkTyADYB\nY1X1YlW9CewBrlvkNq0AvcVuwGmk12Fbx1BmYvJY+pZK8lgDvNT6fLiJSSvM4MQyXW/lzDN/24Sj\nBbfsBszPOeffnxT7x3/8+SK0RFpIU0+r8s47GbhuYuJsksED92ec8Vu8886vZx2f67pVqz7E+Pih\ngXW0vKVq8Am5oI1I/gDYXlVXN5/vAKqqvnbCdovfWElahqpqqLcALpXkcSbwPHAF8AvgKeAzVXVw\nURsmSRpoSVy2qqq3k9wOPEp/HOZuE4ckLV1LouchSVpeFu1uq9k8FJjkW0nGkjybZMNMdZOcl+TR\nJM8neSTJuQuxL0vBPB3PbUkOJ3m6Wa5eiH1ZCuZwPC9vxe9OMpFk/wnbr8jzc56Opedmx//Wk6xN\n8niSnyV5LskXWtt3PzerasEX+knr/wIfAv4Z8CzwsRO2uQb4YVP+BPDkTHWBrwH/rSl/GfjqYuzf\naXQ8twFfWuz9W07Hs/n8r4ENwP4T6qy483Mej6XnZvf/1lcDG5ryB+mPM8/5b+di9Txm81DgdcC9\nAFX1E+DcJKtmqHsdsLsp7waun9/dWDLm63gCrMRJmk7leFJVPwZ+NeB7V+L5OV/HEjw3Ox3Pqhqv\nqmeb+BvAQd57nq7zublYyWM2DwVOtc10dVdV1QRAVY0DFwyxzUvZfB1PgNubru93VsplFuZ2PI8M\n2OZEF6zA83O+jiV4bsIcj2eSD9Pv0T3ZhDqfm0vlCfPZmMv/ZXg3wNRmczx3AR+pqg3AOPCN+W3S\niuP5OXeem3OU5IPA/cAXq+r/TbHZjOfmYiWPI8BFrc9rm9iJ21w4YJvp6o5PdneTrAZeGWKbl7J5\nOZ5V9ctqLoIC3wZ+f4htXspO5XhOZ2IFnp/zciw9N9/V6XgmeR/9xHFfVT3Q2qbzublYyWMfcEmS\nDyU5C9gMPHjCNg8CN8G7T6C/1nSrpqv7IHBLU74ZeICVYV6OZ3MSTfo08NP53Y0l41SO56Rwcu9u\nJZ6f83IsPTfnfDz/EjhQVd8cUOeWpjy7c3MR7xq4mv5o/xhwRxP7PPCfWtv8Bf07C/4e2Dhd3SZ+\nPvBYs+5R4J8v1v6dJsfzXmA//Ts6/pb+mNKi7+syOJ7fA16mPyHVz4H/2MRX5Pk5T8fSc3P2x/Py\nJvZJ4O3mmD0DPA1cPddz04cEJUmdLacBc0nSEmHykCR1ZvKQJHVm8pAkdWbykCR1ZvKQJHVm8pAk\ndWbykCR19v8BkIHBfuoy0wIAAAAASUVORK5CYII=\n", "text/plain": [ "