{ "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": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.hist(t.ravel(), bins=50);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The delays are exponential by construction:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEPCAYAAABP1MOPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAH2pJREFUeJzt3XmcVNWZ//HP0yASlaC4RGVzAWTCjJEZ3MWQQQVX1DAq\naqI4Ucdfgo7GH5hxTHfHmKCTZHRGMK+ZIOOG/ETciCRBf9pGXNG4RUVQUDaDwYUgGKPwzB/nthZF\nVfXt7nur6lZ9369Xvarrnrs8VV23nnvuOfdcc3dERKR+NVQ6ABERqSwlAhGROqdEICJS55QIRETq\nnBKBiEidUyIQEalzSgQiInVOiUBEpM7VRSIws8Fm9jMzG97O5Y4xs6vMbLe0YhMRqTSrhyuLzewl\n4BB3X9eBZQcBN7j7yOQjExGpvLqoEQDbdiQJALj7ImDPhOMRKcjMppvZDyodh9SXekkEVcXM3jSz\nP5tZr7zpz5nZJjPrV6nY2sPMvm1mC6L3cmOM+Xcws7vN7EMzW2pm42KWdTOzX0Sf21oz+52Zje5E\n3EW3FZXfYmZvm9kHZrbQzP6xHes+zMwei5ZdY2aPmtnf5ZQvNbO/70Tsb5rZhuhzeM/M5pvZ+WZm\nHV1n2vQ9qX71kgiK7iRmdlT0T96rjPE4sBTI/RL/NfCFqCwrVgJXAtNizj8V+DOwM3AmcIOZ/VWM\nsq7AMmC4u/cErgDu6ETCLLUtgB8De7r79sAJwA/NbGhbKzWzHsAc4DpgB6A30Ax83ME4C3Hg2Ohz\n6A9MBiYR/39QCfqeVDt3r+kH0AV4rUS5AS+0sY4lCce0FPgX4Omcaf8GfA/YCPSLpu0G3Am8A7wB\nTMiZfxLwOvAn4PfAiQW28V3gBeB94HagW0qf8ZXAjW3Msw3hB3HvnGk3AT8qVVZifS8AJ7X1ObUn\njiLz7wOsAsbG+Bz+DnivRPnN0f93ffR/uxTYD3gWWAvMjP5PP2jju/P3edP2j9b75VKfBzARmJW3\n7HXAtWl8L/Q9yc6jHmoE5wMzSpT/LfC7NtbxWzM7PrmQAHgS6GFm+5hZA3AqcCtR7SWq6s8BniN8\ngUcCF5nZkdHyrwOHuvsXCUedt5rZl/K28Q/AUYQ2jq8AZxcKxMzmmNn70amG/Of7Enq/g4BP3P2N\nnGkvAEPaKCsU75eiZX4f43NqTxy525hiZuuBVwk7+NwY73ERsNHM/sfMRpvZ9rmF7v5NwhHrcdH/\n7TrgHsIPTC9gFvD1GNvZjLsvAFYAw9v4PGYCR5vZttF7bCB8R26Lsx19TxL7nlSdmk4EZnYpsIu7\nN+dNH2hmV5rZ0YTq44Ol1uPuZwNjzGxswiHeApwFHMnnX6RW+wM7uftV7r7R3d8EfgGcFsU0291X\nR3/PAhYDB+St/zp3X+3uHxB2gv0KBeHux7v7Du7eq8DzCQm91+0IR8G5/gT0aKNsM2bWlZAwp7t7\n63su+jm1M47PuPu3o3kPA+4ixukdDx0SDgM2Af8FvGNm95rZLvlvI3o+GOjq7v8RxT4bWNDWdopY\nRUgmRb837r6M8EN4UrTMSGB9lEjapO9JMt+TalTTicDdfwK8a2ZXtU4zs20IR14/cfdfEf6x/z+n\n/GAzOzh3PVED173ufmfCId4KnE44Ur85r6w/0Ds62nrPzN4nnDraJYrpmxYal9+PyoYAO+WtY3XO\n3xsIX9hK+RD4Yt60nsC6qKxnkbLPREd1txJ2tgnR5H6U/pxON7N1ZvYnM7s/7rYAPHgc6AtcEOdN\nuvtr7n6Ou/cD/hrYHfj3IrPvRjh/nuutONspoDfwHm18bwi149a2qXGUri1XQl18T6pN10oHUAZT\ngNeAy6PXJwMvuftaM+sG9HD3P7TO7O5PFFjHCHc/J+nA3H2ZmS0Fjgby17+c0DaxT/5yUePXfwFf\na43XzJ6jRKN4KWY2FxhO4YbqR9392I6sN88ioKuZ7Z1T3f4K8HIbZbmmEZLdMe6+MZpW9HMCcPcZ\n5PzYRQcCXWJsK1dXYO84bzJv24vM7H+A83In5/z9NuEHPFc/wmm/2Mxsf0LCmU84sCn6eRAdBJlZ\nb0LN4KB2bEffkxS+J9WgpmsEAO6+idBg3Gpn4Pno75HAU2Y2CsDMhpnZj6MjinI5h9D491He9KeB\ndWY20cy6m1kXMxtiZsOAbQmnH9aYWYOZjSccfXaIux/j7j3c/YsFHkV37iim7oTPt6uZbW1mXQrN\n6+4bCFXnH5jZNmZ2GHA8cHNUNrtA2S052/o5MBg4wd3/EvNzak8ct0Tb2dnMTjWzbaPPdhTh9EHJ\n04fRsvuY2SXRjyxm1pdw1J17cPEHoLWH2hPAp2Y2wcy6mtnJbHl6r9T2epjZcYQG5lvc/eW2Pg93\nXwM8Akwn/DC+Fnd7+p4k8z2pSl7h1upyPMjp9QPsSqiqHw2Mj/7+h6hsd2BKqeWTioe8nh/R9C5s\n3mtoV8JRytvAu8DjrcsRemC8S+gB8RPgYeCcYtsAGgk7U5Lvo5GQkDbmPL6fUz4XuCzn9Q7A3YRq\n95vAqTHL+kXb2UComq8jnK8d19bnVCTuUtvaCWghnGb5gNBAeE7Mz2N34P8RGm7XEY5CpwLb5cxz\nAuH0z3vAJXzeWWEt4Qd9s15DBT7DpYReR2sJvcEeA/6JaJSAOJ8HoSvkRuCSAu9hs+3pe5L896Qa\nH6kMMREdUV9JONe3wN1vaWORVJnZUndv8+pgM9uDqBeFu6/KmR5reRGRLErr1NAYoA/wF8LRUaWt\nN7P8hp9CdiYcbX2WHc1sH8LRgIhITYqVCMxsmpmtNrMX86aPtnBp9SIzm5RTtA/wmLtfCvyfBOPt\nqJOAy83skFIzufsCd5/q7m9DGH2UcPro1DLEKCJSEbFODUUNJR8SzjHvG01rILTijyT0YV5A6Ku8\n0MzOAD529zvNbKa7F+urKyIiFRarRuDu8wkNU7kOABa7+1vu/gnhqsUxUdldwGgzu47QQ0FERKpU\nZ64j6E3oFdFqBVHXNw9dIb/ViXWLiEiZVOyCMjPL0iibIiJVw90TvdapM72GVhL67bbqw5aXy5fU\n2NjIww8/XPE+tB15NDY2VjwGxV/5OBR/9h5Zjf3hhx+msbGxEz/ZxbUnERibD2GwABhgZv2joRpO\nA9o1AmFTUxMjRoxozyIiInVpxIgRNDU1pbLuuN1HZxCuwhtkZsvMbLyHMTwmAPMI42/MdPdX27Px\npqYmWlpa2hmyiEj9aWlpSS0RVOzm9Wbmldp2ElpaWjJdm1H8laX4KyfLsQOYGZ5wG0FFE0FjYyMj\nRozI9D9FRKQcWlpaaGlpobm5ubYSQZZrBCIilZBGjaDmh6EWEZHSKpoI1FgsIhKPGotFRATQqSER\nEUmBTg2JiGSATg2JiAigU0MiIpICnRoSEckAnRoSERFAp4ZERCQFSgQiInVOiUBEpM6psVhEJAPU\nWCwiIoAai0VEJAVKBCIidU6JQESkzikRiIjUOSUCEZE6p+6jIiIZoO6jIiICqPuoiIikQIlARKTO\nKRGIiNS5rmms1My+ClwJvAzc7u6/LTTf3nsXX8dZZ8H3v59GdCIikiuVRAA4sA7YGlhRbKZ58wpP\nb2mB2bPTCEtERPLFSgRmNg04Dljt7vvmTB8NXEs4xTTN3a8GiGoAvzWzXYCfAWcWWm+xGsFrr8V/\nAyIi0jlx2wimA6NyJ5hZA3B9NH0IMM7MBuct9wHQrbNBiohIemLVCNx9vpn1z5t8ALDY3d8CMLOZ\nwBhgoZmdREgQPQnJQkREqlRn2gh6A8tzXq8gJAfc/W7g7o6uuFs3ePppGDWq+DyTJ8PQoR3dgoiI\ntEqrsTiW3MulR4wYwYgRIwAYORLuuAM++aTwclOnwqOPKhGISO1raWlJfSie2ENMRKeG5rQ2FpvZ\nQUCTu4+OXl8GeGuDcYz1dXiIiQsvhAEDwrOISD2p9BATFj1aLQAGmFl/M+sGnAbc156Na9A5EZF4\n0hx0LlYiMLMZwOPAIDNbZmbj3X0jMAGYR7hwbKa7v5pKlCIikppMjj564YWwZg0ceWTh8p494aST\nwBKtPImIVF4ap4YymQgeewx+8Yvi5bNnw7PPwsCBHQxORKRKpZEIKt5rKLe3UFyHHhoexTz2GOhW\nByJSS9LsPZTJGkFbBg2CX/4yPIuI1JJK9xpKnHoNiYjEo1tVttOgQXDppdC7d+HyAQNgn31S2bSI\nSKpqro0gLWedBffeW7hswwZYtUojnIqItMpkY3FbLr+8eNmSJXDEEYluTkQkdWosTlBrIliypOyb\nFhHptJprLBYRkcqryTaCtnz6KSxfXrx8111hq63KF4+ISCXVZBtBKb16wfbbwyGHFC5ftw4mTIAr\nryxbSCIibVIbQRlNnQq//314FhGpNmojEBGRxCkRiIjUOQ0xISKSARpiooxuvhkuvhj69i1c3rUr\n3H67hrgWkcrQEBNlcOaZsO++xcsvughefVWJQERqhxJBnoYG2G+/4uU9e5YvFhGRclBjsYhInVMi\nEBGpc2osbqcTT4S1a4s3JvftCz/8IViiTTkiIkHN3by+sbGx7ENMdNbrr8PjjxcvHz8ePvoIunUr\nX0wiUvtah5hobm6urUSQxRpBW7p1gw8/VCIQkXRoiAkREUmcEoGISJ3TdQQpuOuu4vczGDoU9tqr\nvPGIiJSSWhuBmW0DPAI0uvvcAuU12UYwaVJoUC7knXege3d44IHyxiQitSNTvYbMrBlYB7xST4mg\nlIceCl1LH3qo0pGISFZVrLHYzKaZ2WozezFv+mgzW2hmi8xsUs70I4BXgD8C6lEvIlLF4jYWTwdG\n5U4wswbg+mj6EGCcmQ2OikcABwKnA99KJFIREUlFrMZid59vZv3zJh8ALHb3twDMbCYwBljo7v8a\nTfsmsCbBeDPvL3+BNSU+kR131FXJIlJenek11BtYnvN6BSE5fMbdb+7E+mtO376wdCkMHly4fP16\nuOEGOPvssoYlInWuot1Hc++2k7WhJjpi4EBYubJ4+cSJoWeRiEir1qEl0hS711B0amiOu+8bvT4I\naHL30dHrywB396tjrq/ueg21ZeJE2Gmn8CwiUkilh5gwNu8BtAAYYGb9zawbcBpwX3s2rnsWi4jE\nU/F7FpvZDEJPoB2B1YSLxKab2dHAtYSEMs3dJ8fesGoEW5g4Ed57D049tXD5dtvBwQeXNyYRqS6Z\nuqCszQ1ndBjqND3wAFxzTfHyJ56AZ54p3tgsIrVLw1ALAPvuC7feGp5FpD5Vuo0gcWojEBGJp+Jt\nBKlsWDWCdlONQERUIxARqVOqEQgQagLHHQe9excu//KX4WtfK29MIlJeNddrSImgfW67DR5/vHDZ\n+vUwbx6sWlXemESkvNJIBBUfYkLdR+M744zwKGTVqpAIRKQ2pTnUhGoENWLVKhg2TDUCkVpXc43F\nIiJSebp5fQ1ZuxYuvrh4+SmnaIgKEdmS2ghqxG67wc9/XvymN089Bf/930oEIlmlNgLptBtvhPnz\nw7OIZJfaCEREJHFKBCIidU5tBHVk6VK4887CZV26wPHHQ1d1HxCpSmojkE57/XX43veg2Ef+5JNw\n3XXw9a+XNy4RaZ+au7JYymfAAJg1q3j5KafAxo3li0dEqofaCERE6pwSgYhInVMbgQBw6qnhYrS9\n9ipcvttu0NQEDTp0EKkoDUMtqVm8GEp1SPjOd0Ki6NGjbCGJSAE111is7qPVY+DA8CjmkkvKF4uI\nbEndR6XievQIQ1yrRiBSWRpiQkREEqfrCCS2WbPgC18oXDZ0KAweXN54RCQZOjUksTQ3w8KFhcv+\n+Ef4+GN49NHyxiRSjzLTa8jMBgMXATsCD7n7zwvMo0RQIx57DCZODM8ikq7MtBG4+0J3vwA4FTgk\njW2IiEgyYiUCM5tmZqvN7MW86aPNbKGZLTKzSXllxwO/BOYmF66IiCQtbo1gOjAqd4KZNQDXR9OH\nAOOiU0IAuPscdz8WODOhWEVEJAWxeg25+3wz6583+QBgsbu/BWBmM4ExwEIz+ypwMrA1cH+C8UoV\n6tIlNCSXGsL6sstg//3LF5OIxNeZ7qO9geU5r1cQkgPu/gjwSCfWLRly4IFwyy3w0UeFy2+/HX79\nayUCkWpV8SEmWmmoiewyg2OOKV7+/PPli0Wk1qQ5tESr2N1Ho1NDc9x93+j1QUCTu4+OXl8GuLtf\nHXN96j5aJ664Arp1C88i0jmV7j5q0aPVAmCAmfU3s27AacB97dl4U1NT6plORKQWtLS0bHYWJUmx\nagRmNgMYQbhAbDXQ6O7Tzexo4FpCQpnm7pNjb1g1grpxxRVhiOtDDy1cvs024YK07t3LGpZIJlVs\nGGp3P73I9F8Bv+roxjUMdX0491zYdtvi5ddfD0cfrcZkkVI0DLXUtP33h6lTlQhE4qh0G0Hi1EYg\nIhJPxdsIUtmwagQS2X9/+Pa3YciQwuW77gp9+5Y3JpFqpVtVSk065hiYMqVw2aefhnslL19euFyk\nXqiNQOrWu+/CoEHhWURqsI1AREQqT43FIiIZoMZiqVs6NSSyOZ0aEhGRxFW015BIW7p2hQ0b4OCD\ni89z3nkwfnz5YhKpNeo+KlWtZ0947jl4//3C5b/5DTz4oBKB1D51HxUp4rbbYO7c8CxSD9RGICIi\niVMbgWTe+vWwbFnhsoYG6NOnvPGIZI1ODUmmPfMMjB0LmzYVLl+zBmbNgmOPLW9cImnRWEMieYYN\ngzffLF7+jW/Ae++VLRyR1KixWKSDvvENOOqo8CxSC9RYLCIiiVMiEBGpc+o1JDXvmWdgu+0Kl+28\nMxx2WHnjEak2aiOQmnbPPXDTTcXL586FVatgxx3LF5NIZ6jXkEg7nXhieBSz887Fu56KVBP1GhJJ\nyc47wyuvhGeRLFCvIRERSZwai6Xu/fSnsO22hcsOOgiOPLK88YiUm04NSV27/fZwaqiQt98OQ2A/\n+2x5YxIpJY1TQ0oEIkU8+2y46Y0SgVSTTPUaMrMxwLFAD+BGd38grW2JiEjHpZYI3P1e4F4z2x74\nN0CJQDJn9Wr42c+Kl59yioa5luyLnQjMbBpwHLDa3ffNmT4auJbQA2mau1+dt+i/AlMSiFWkrIYM\nCbfAXLGicPlTT8G778JVV5U3LpGkxW4jMLPDgA+Bm1sTgZk1AIuAkcAqYAFwmrsvjMonA/Pc/aEC\n61MbgWTaVVfBhg1KBFJeFb2OwN3nA/m3ED8AWOzub7n7J8BMYAyAmU0gJIixZnZeQvGKiEjCOttG\n0BtYnvN6BSE54O7/CfxnqYWbmpo++1tDTYiIbCnNoSVaVXysIRERKS7/ILm5uTnxbXR2iImVQL+c\n132iabE0NTWlnulERGpBS0tLagfP7bqgzMz2AOa4+99Er7sArxHaAt4GngbGufurMdalxmLJNDUW\nSyVU9IIyM5sBjAB2NLNlQKO7T48ahefxeffRNpNAKw1DLVnWqxdcfTXcemvh8u7dYd486N+/vHFJ\nbdIw1CJVaNMmWL68ePnYseFitOHDyxeT1L5MDTERh2oEkmUNDaWP9rt3L18sUvtUIxDJoOHD4Uc/\nUo1AkqUb04iISOIqmgjUfVREJJ6q6T6a6IZ1akhq3PDhcPrpsN9+hct33RX23LO8MUn26cY0Ihky\neTLce2/hso0bYckSWLOmvDFJ9tVcImhsbFSvIalL69fDLruEZ5E4WnsNNTc311YiUI1A6tWf/ww7\n7QT9+hWf59xz4eKLyxeTZEPN1QiUCKSerVwJa9cWLps3Dx59FGbPLm9MUv10QZlIDendOzwKeTX2\nQC1SL3RBmUidmT0bZsxQjUC2pAvKREQkcUoEIiJ1TolARKTOaYgJEZEM0BATInVm9my45hr47ncL\nl5vBCSfA1luXNy6pvJrrPioihR14IOy1F9x5Z+HyBQtg3To455zyxiW1SYlApAr16QO33168/Nxz\nw3hFIklQY7GISJ1TIhARqXMaYkIko155BR54oHDZDjvAsGHljUfSpSEmRGQzd90FN9xQvHz+fFi6\nNNz8pr0++ABmzSpe3q0bnHEGdFULY0Wo15CIAHDyyeFRTN++8MknHVv3jTfC9Olw0EGFy+fODYPl\nHXFEx9Yv1UeJQEQ24w5HHQU//Wnh8iOPDPNI7VAiEJF2mzULXnyxcNmee5aurUj1USIQqVFz5oS7\noLXX88+H22gWc+ml4cY5q1ZtWbZpE0ycqGscsiaVxmIz2xO4HPiiu59SZB41Fouk5Oqr4dlnO778\nd74Dhx/e/uU2bgyNyUoE6cncrSrN7A4lApH6oUSQvordmMbMppnZajN7MW/6aDNbaGaLzGxSkoGJ\niEh5xL2yeDowKneCmTUA10fThwDjzGxw3nKJZi0REUlerETg7vOB9/MmHwAsdve33P0TYCYwBsDM\nepnZDcB+qimIiFS3zvQa6g0sz3m9gpAccPf3gAs6sW4RESmTio811EpjDomIbCnNMYZadSYRrAT6\n5bzuE01rFyUAEZHiWn8jq2LQOTPbA5jj7n8Tve4CvAaMBN4GngbGufurMden7qMiNUbdR9NXye6j\nM4DHgUFmtszMxrv7RmACMA94GZgZNwm00s3rRUTi0c3rRSQTVCNIX8VqBGlRjUBEJB7VCEQkE1Qj\nSF/N1QhERKTydGpIRCQDdGpIRDJBp4bSp3sWi0jVa2iA444rXj52LJx9dtnCkRgqPsSEriwWqR1d\nusCTT8LbbxcuX7Ag3OZSiaD9quLK4sQ3rFNDInXn/vth6tTwLB2jXkMiIpI49RoSEckA9RoSkZqg\nU0Odp1NDIiKSOCUCEZE6p+sIRKRqvPEG3HNP8fLtt4fx48O1Cvnc4cYb4YMPii8/ZgwMGND5OGuN\nriMQkapx7bXw8sswdGjh8iuvhMMPh4EDtyxbsgQuuQS+9a3Cyz73XJhnypTk4i0nXUcgIjWhrcbi\nCRNg0KDwXMjAgTB3buFE8MYbcNRR4bmQKVPglVeymwhaqbFYREQSp0QgIlLnlAhEROqcEoGISJ1T\n91ERyYyddoJDD4WuBX65Pv0U9tyz9LLTp8Pddxcu32qr0BA9ZEgysWaJuo+KSGY8+CCsXVu8vGfP\n4mWnnALDhxcvP+us0L20WhNBmt1HK54IRETi2nbb8OgIM9h99+Ll3bt3bL3l0nrQ3NzcnPi61UYg\nIlLnlAhEROqcEoGISJ1LpY3AzLYBpgIfA4+4+4w0tiMiIp2XVo3gZGCWu58PnJDSNkREJAGxEoGZ\nTTOz1Wb2Yt700Wa20MwWmdmknKI+wPLo740JxVpVsn6LTcVfWYq/crIce1ri1gimA6NyJ5hZA3B9\nNH0IMM7MBkfFywnJACDRUfKqRda/TIq/shR/5WQ59rTESgTuPh94P2/yAcBid3/L3T8BZgJjorK7\ngbFmNgWYk1SwIiKSvM40Fvfm89M/ACsIyQF33wCc04l1i0gN2morWLAAjj++cPlLL8Fll5U3prj+\n+Z+L3+sA4IQT4NxzyxdPkmLfmMbM+gNz3H3f6PXXgVHufl70+kzgAHe/MOb6dFcaEZEOSPrGNJ2p\nEawE+uW87hNNiyXpNyIiIh3Tnu6jxuYNvwuAAWbW38y6AacB9yUZnIiIpC9u99EZwOPAIDNbZmbj\n3X0jMAGYB7wMzHT3V9MLVUREUuHuiTyA0cBCYBEwqcg8/wEsBp4H9mtrWWAHQqJ5DfgN0DOpeMsU\n/zXAq9H8s4EvZiX2nPLvApuAXln67KOyCdHn/xIwOUvxA18BngCeA54GhlVR/ENzpk8DVgMv5s1f\nzftunPjLsu+mFX9Oeaz9N6k30gC8DvQHtoqCHZw3z9HA/dHfBwJPtrUscDUwMfp7Ulo7c4rxHwE0\nRH9PBn6cldij8j7Ar4GlbX2Rqi1+YAThh6hr9HqnjMX/G+ConOUfrrb4o9eHAfux5Q9p1e+7bcSf\n+r6bZvxRWez9N6khJkpdU9BqDHAzgLs/BfQ0sy+1sewY4Kbo75uAExOKtyzxu/uD7r4pWv5JPr/I\nrupjj/w78H9TiLkc8V9A+PH5NFpuTcbi3wS03mZle9rREaOM8eOFrzFqXaba992i8Zdp300t/kjs\n/TepRFDomoLeMecpteyX3H01gLv/AdgloXjzpRV/rnOAX3U60i2lEruZnQAsd/eXkg44Zmxx5im1\n7CDgcDN70sweNrNhiUbddmxx5im17MXAT8xsGeE0xfcSjDlObKXmWVlgnny7VPG+Gyf+XGntu5BS\n/O3dfyt5h7KOdB+tpmsPYsdvZpcDn3j1jMJaMnYz+wLwL8CRcZcpszixdAV2cPeDzGx/4A5gr3TD\nii1O/BcAF7n7PWY2FriRzf8fWVNN+25sVbjvtqkj+29SNYI41xSsBPoWmKfUsn9orQKZ2a7AOwnF\nmy+t+DGzs4FjgNOTC3eLuJKOfW9gD+AFM1saTX/WzNI4qkvrs18B3AXg7guATWa2Y3JhbxZbGvGf\n5e73ALj7nURX7aegM/GXsjoD+25JZdh3IZ3427//JtTg0YXPGzy6ERo8/ipvnmP4vMHjID5vMCu6\nLKHBaZKn3+CUVvyjCV1rd0wj7jRjz1t+KeHoOjPxA+cDzdHfg4C3MhJ/a2Pxy8BXo79HAguqLf6c\n8j2Al/KmVf2+20b8qe+7acafV97m/pvkGxpN6Cq2GLgsmnY+cF7OPNdHb/oF4G9LLRtN7wU8GJXN\nA7ZP8R+SRvyLgbeA30WPqVmJPW/9S0i/+2jSn/1WwC2ErqPPEP2oZij+Q6K4nyN0Ix1apfHPAFYR\nbkK1DBgfTc/Kvlss/rLsu2nFn7f+Nvff2GMNiYhIbdI9i0VE6pwSgYhInVMiEBGpc0oEIiJ1TolA\nRKTOKRGIiNQ5JQIRkTqnRCAiUuf+F1iaVFCeYB/rAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.hist(d1.ravel(), bins=40, histtype='step');\n", "plt.yscale('log')\n", "plt.title('$\\{d_i\\}\\quad$ Mean = %.3e Std.Dev. = %.3e ' % (d1.mean(), d1.std()));" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here we compare a set of estimators:\n", "\n", "$$\\hat{\\Lambda}_{n,c} = \\frac{n - c}{T_n}$$" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false, "scrolled": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " c n Mean Std MSE Median % > λ\n", " -1.00 5 1.49712 85.90% 99.247% 1.282 71.34%\n", " 0.00 5 1.24760 71.58% 75.744% 1.068 55.67%\n", " 0.33 5 1.16442 66.81% 68.804% 0.997 49.76%\n", " 1.00 5 0.99808 57.27% 57.267% 0.854 36.93%\n", " 2.00 5 0.74856 42.95% 49.769% 0.641 18.41%\n", " 3.00 5 0.49904 28.63% 57.702% 0.427 5.27%\n", " median 5 1.92114 181.83% 203.832% 1.440 73.49%\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAEPCAYAAACneLThAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xd4FNX6wPHv7KZX0gglxVBFWujSBAVpgkhRpFjRi+Wi\niHotoICK1x9yEbEhIhcQRLmAICpF0EgHAQNIDS0JCTW91/n9McmW9LKb+n6eZx/3zJw5cxLiu7Nn\nzpxXUVUVIYQQdZOuujsghBDCeiTICyFEHSZBXggh6jAJ8kIIUYdJkBdCiDpMgrwQQtRhpQZ5RVH8\nFEX5TVGUk4qinFAU5YUi6vRTFCVeUZSjea+Z1umuEEKI8rApQ51sYLqqqqGKorgARxRF2a6q6pkC\n9Xapqnq/5bsohBCiokq9kldV9ZqqqqF575OB00DTIqoqFu6bEEKISirXmLyiKLcBwcDBInb3VBQl\nVFGUnxVFucMCfRNCCFFJZRmuASBvqGYd8GLeFb2pI0CAqqqpiqIMBTYCrSzXTSGEEBWhlGXtGkVR\nbICfgC2qqn5chvqXgC6qqsYW2C4L5QghRAWoqlqhIfGyDtcsA04VF+AVRfE1ed8d7cMjtqi6qqrW\n+NesWbOqvQ9lecnvs/71szb0Ufpp+VdllDpcoyhKb2AicEJRlL8AFXgTCNRijLoEGKsoyrNAFpAG\njKtUr4QQQlhEqUFeVdW9gL6UOp8Bn1mqU0IIISxDnngtQv/+/au7C3VKbfl91oZ+1oY+gvSzJinT\njVeLnUxR1Ko8X12nKEqlx+uEEDVf3v/rFbrxWuYplEKI+um2224jPDy8urtRLwQGBnL58mWLtlmn\nruSzsuLIyrpVaLu9fVP0eiernbe6yJW8qAryd1Z1ivtdy5V8nmvXlnPhwvRC24ODQ2jQoF819EgI\nIapXnQry+Vq1WoJe70RKyikiIt6v7u4IIUS1qZOzaxo2fAhf34l4eAys7q4IIUS1qpNX8kII65k2\nDUJDrXuO4GBYuNC656gv6uSVvBDCekJDrRvkrd1+ZWVlZfHggw8SFBSETqdj165dJdaPi4tj1KhR\nuLi4EBQUxJo1a6qopxq5khdClFtwMISEWKft2vB8Ut++fXnppZd48MEHS6373HPP4eDgwM2bNzl6\n9Cj33XcfwcHBtGnTpgp6Wk+u5M+fn05o6EDD6+bNH6q7S0IIC7ly5QpjxoyhYcOG+Pj48MILhTKU\nWpStrS0vvPACvXr1QqcrOYSmpqayYcMG3nvvPRwdHenduzcjR47km2++sWofTdXpIG9j446bW290\nOkdyc9PJyUkiPn4nGRmR1d01IYQF5ObmMnz4cIKCgoiIiCAqKoqHH364yLpr1qzBw8MDT09PPDw8\nzN57enpy5coVi/fv3Llz2Nra0rx5c8O2jh07cvLkSYufqzh1erjG1bUznTvvMZSzsmLZu9erGnsk\nhLCkQ4cOcfXqVebNm2e4qu7Vq1eRdcePH8/48eOrsnskJyfj5uZmts3NzY2kpKQq60OdvpIXQtRt\nkZGRBAYGljpsUpn2XV1dcXV1LRSsy8LFxYXExESzbQkJCbi6ulqqi6WSIC+EqLX8/f2JiIggNze3\n1LrffvutIVibvvK3FTVc4+/vT1JSEklJSYWCdVm0atWK7OxsLly4YNh27Ngx2rZtW+62KqpOD9cI\nIawjNNR6s2BCQ7XZO2XRvXt3GjduzOuvv87s2bPR6/UcOXKkyCGbCRMmMGHCBIv0MTMz0/DBkpGR\nQUZGBvb29oXqOTk5MXr0aN5++22++uorjh49yubNm9m3b59F+lEWciUvhCiX4OCyB2Frt6/T6di8\neTNhYWEEBATg7+/P2rVrrde5PK1bt8bZ2Zno6GiGDBmCk5MTERERAPz73//mvvvuM9T97LPPSE1N\npWHDhkyaNInFixdX2fRJqGOrUEZGfsSFC9Pp0yceGxv3Qvvzb7y2aPExfn7WnWZVFWR1QFEV5O+s\n6lhjFUq5khdCiDpMgrwQQtRhEuSFEKIOkyAvzMTExLB06dLq7oYQwkJkCqUws2bNGqZNm8aoUaPw\n8qqfTwd//PHHHD58mH79+tGxY0d27txJTk4Op0+f5r///S+2trZWOW9ycjLvv/8+iqLg5OSEnZ0d\nr776qlXOJeoPuZIXZjZu3IiqqmzevLm6u1Itzp8/T5MmTRg3bhxvvPEGx48f5/XXX2fGjBmcOnWK\nLVu2WOW8iYmJ9OvXj+DgYObOnUvnzp15++23iYuLs8r5RP0hQV4YJCQkGNbG3rhxYzX3pnocOXKE\nQYMGceLECVq1asXkyZMN++Li4lCUCs1iK9WLL76In58fDz30EADu7u689NJLeHh4WOV8ov6Q4Rph\n8PPPP9O6dWtOnjzJjh07SE9Px8HBobq7VaXGjRsHwO7duxkxYoRhe3R0NJGRkXTp0qXI49atW8f6\n9euL/RBQVRVFURg2bBiTJk0y23ft2jVWrVrFpk2bDNt69epV7EJbQpSHBHlhsHHjRj799FMGDRpE\nWloa27dv5/7776/ublU5VVXZt28fb775pmHb5s2badu2LU2aNCnymLFjxzJ27NgKne/QoUPk5uZy\n1113Vej4qjZt6zRCr1k3dVNwo2AWDpH8f5YgwzUC0NbiOH36NP369aN/3qIk9XXIJjQ0lIyMDLp3\n727Y9v333xuC+Oeff27R8+Xm5uLq6oqLi4vZ9vDwcNLT0y16LksIvRZq1SBv7fYr4rXXXsPb2xsf\nHx9ef/31Euvu3LmTNm3a4OLiwoABAwzLHVSXWn0lHxHxIXFx2w3l9PTwauxN7bZjxw7uueceAEaO\nHMmvv/7KTz/9ZBhmqE/27NlD586dsbOzM2z7888/ee+999i8eTO9e/e26Pn69euHXq/n+vXr+Pr6\nAnDmzBmWLVvGvHnzLHouSwluFEzI4yFWabv/8v5WabeivvzyS3788UdOnDgBwMCBA2nWrBn/+Mc/\nCtWNiYlhzJgxLFu2jOHDhzNz5kzGjRvH/v37q7rbBrX6Sj4l5STx8bvIyUklJycVW1sf3Nx6Afrq\n7lqts2nTJkaOHAloQV5RFGJiYtizZ08pR9Y9kZGRhYZennnmGb7//nuuX79Ox44dLXo+Dw8PNmzY\nwPTp05k7dy7vvPMOe/bsqbEBvqaxdvq/lStX8vLLL9O4cWMaN27MK6+8wvLly4usu2HDBtq1a8fo\n0aOxs7Nj9uzZHDt2jHPnzlm0T+VRq6/kAezsGtO5897q7katpqoqu3btMgxDNG3alC5dunDkyBE2\nbdpE3759q7mHhS1btoxvvvmm3N8yVFVFr9ezdOlSbrvttiLrFBVcP/zww4p0s8z69etHv379rHqO\nuig//d/AgQNZvXo1Op2Ow4cPF1l3zZo1PPfcc2aLgOW/VxSF48eP4+fnV+i4kydPmn2wl5S+r2Bd\nJycnWrRowcmTJ2nVqlVlftQKq/VBXlTegQMH6NKlC3q98RvQyJEjOXz4MJs2bWL+/PnV2LuiPfnk\nkzz55JPV3Q1Rzaoi/V9ycjLu7sZVbd3c3EhOTi62bsOGDc22VXW6v4Jq9XCNsIyNGzcahmryPfDA\nAwBcvHixSpMOC1Ee1k7/B4VT+CUkJBS6SV5c3fz6VZnuryAJ8oItW7YwdOhQs21t27Y1ZJg3nb8t\nRE1i7fR/oP2/cOzYMUM5NDS02PR9bdu2JTTUODMoJSWFCxcuVGm6v4JkuKaeO3PmDE2bNi3yymTk\nyJEsWLCATZs2mc0Zr8t0Ol2FZxNZYyZSfps5OTkWbbeyQq+FWm0WTOi1UIIblS01VFWk/3v00UdZ\nsGABQ4cORVVVFixYwLRp04qsO2rUKP71r3/xww8/MGzYMObMmUNwcHC1jceDBPl6r6ihmnwPPPAA\nCxYs4MiRI0RHRxd6EMiSC3mVd3Gur7/+mlWrVlX4xutXX31FUFBQof1luSKs78oagCvTflnPkZ/+\nb+rUqQQEBKDT6ZgwYYJFnxaeMmUKly5don379iiKwtNPP83TTz9t2N+uXTtmzJjB+PHj8fb2Zv36\n9Tz//PNMmjSJHj168N1331msLxWiqmqJL8AP+A04CZwAXiim3iIgDAgFgoupo1rSqVOPqfv2BZa5\nfmZmjPr776iRkR9btB/VxRK/zx49eqhRUVFF7svNzVUbNmyo6nQ69YsvvjDbFxYWpq5du1bdvHmz\n6u3trS5dutSwr1OnTuqmTZvK3IeEhAS1c+fO6vfff6+qqqr+8ssvqoODgxobG1uBn0hYmqX/vxXF\nK+53nbe91Hhd1KssY/LZwHRVVdsCPYHnFUW53bSCoihDgeaqqrYEpgCLK/vhY02Rkf/hyJEehtfl\ny+9Vd5eqxdq1azl06BCxsbFF7lcUhZ49e6KqKqtWrTK7yrXkQl61cXGuGTNm4Orqyn333cejjz7K\noEGD0Ol0dO3alUcffZRx48bh6+vLgAEDqruror4r76cCsBEYUGDbYmCcSfk04FvEsZb6wFNVtfxX\n8llZieqxY0PMXiEh9urp009atF9VpaK/z1OnTqm33367qtPpVJ1Op9rb26sTJ040q/Puu++qbdq0\nMdTR6XSqn5+f+uCDD5rVGzp0qPrvf//bUI6KilL1en2x3w4Kunr1qmpjY6P+/PPPFfpZqsugQYPU\n69evG8rr1q1TdTqdevbsWcO28PBw9YknnqiO7lmUpf+/FcUr7ndNJa7kyzUmryjKbUAwcLDArqZA\npEk5Km/b9fK0b202Nq506GC+Hvj+/f7V1Jvq06ZNG06fPl1inZkzZzJz5swS66gVWMiroNq2OBdo\nsysmTJhgNh/6jz/+oGHDhmY32Hx9fWnZsmV1dFEIgzIHeUVRXIB1wIuqqhb9JEAZzJ492/C+f//+\nhsWwRO1TloW8nnvuuRLbKGlxLl9f3xq51PHOnTt59tlnzbbt3r270JPBV69epXXr1lXZNVFHhISE\nEBISYpnGynK5j/ZhsBUtwBe1v+BwzRlq4HBNUfbt86t3wzWWsmjRIrVXr15m21xcXNS9e/eqP/74\noxoaGqqqqqouXbpU9ff3V2NiYgq1ERsbq3p6eqrXrl0zbDt9+rT66quvWrfzFhQfH6/q9Xp10aJF\n1d0Vq6juv7P6pLjfNVUwXLMMOKWq6sfF7P8ReB74XlGUO4F4VVVr1FCNsLySFvJq3769IemGqqpk\nZGRw4MABhg0bZlbfdHGuO+64g5ycHJo0aVKrFufas2cPqqrWqiEnUX8oat5CPcVWUJTewC606ZNq\n3utNIBDt02VJXr1PgSFACvCEqqpHi2hLLe185XH69OPEx4fQs+flCrexf78/Hh6DuP32ry3Wr6pi\nutBSTZeUlMTOnTsNyyXUJa+//jpLliwpdpZSbVeb/s5qu+J+13nbK/SkXalX8qqq7qUMa/eqqvrP\ninRA1A9//PFHsanzartdu3ZZfI15ISxFnngVVpeVlcWxY8cYPnx4dXfF4tLS0jhy5AjvvVePnrWY\nNg1CrZy5KTgYFkr6P0uQBcqE1UVGRhaZRacuOHDgANnZ2fVrLfjQUOsGeWu3XwFlTf93+vRpunXr\nhqenJ15eXgwaNKjU6crWJlfywuqaNWtW3V2wmpCQEBwdHYscijp06BB79uwhIyOD1q1bk5ycTEJC\nAlOnTq2GnlpYcDBYaopfQTVsWnV50v81bdqUtWvXEhQUhKqqfPrppzz88MNmq1hWNbmSF6Kcrly5\nwsSJExk1ahQfffQRWVlZjB07lkceecQsmUR8fDxBQUGcP3+e0aNHM3z4cHbt2lWNPa+balL6Pzc3\nN8PCdzk5Oeh0Oi5cuGDR/pSXXMkLUU5+fn6sXr261HqDBg3i+eefNyxvu2PHDvr06WPt7tUrNS39\nXz4PDw9SUlLIzc3l3XffrcRPWHlyJS+EFe3evdsQ2FevXs348ePZtm1bNfeq7jBN/+fg4ICdnV2J\n6f/i4uKIjY0lLi7O7H1sbGyRAR7Kl/4vX1xcHAkJCXz66acWT/xeXnIlL4SVxMbG0qJFC+zt7QEI\nDAxk+/btDBo0qJp7VnfUtPR/phwdHZkyZQo+Pj6cOXMGb29vq/WxJHIlL4SVeHp6smHDBkN50aJF\nTJo0qVCiZ1FxNS39X0E5OTmkpqYSFRVVth/ICuRKXghRfqGh1psFExqqzd4pg5qW/m/Hjh14e3vT\noUMHkpOTmTlzJp6enrRp06bc57UUuZIXQpRPcHCZg7C1289P/xcWFkZAQAD+/v6sXbvWot2ZMmUK\nI0aMoH379nTs2JH777+/UPq/NWvWANqMqvHjx9OgQQNatmzJpUuX2Lp1K3Z2dhbtU3mUunaNRU8m\na9dYlKwpIqqC/J1VHWusXSNX8kIIUYdJkBdCiDpMgrwQQtRhEuSFEKIOkyAvhBB1mAR5IYSowyTI\nCyFEHSZBXggh6jBZ1gC4eXM9iYn7DWUnpzto125dNfZIiJprWlgYoaWswlhZwS4uLGzZ0qrnqC/q\n/ZW8l9dwPD0H4ezcDmfndmRnx5Gefqm6uyVEjRWanGzVIG/t9ssrJCSEe+65hwYNGpQpy9nOnTtp\n06YNLi4uDBgwgIiIiCroZfHq/ZV8q1ZfmJVPnBhBRkZ0NfVGiNoh2MWFkE6drNJ2/7/+skq7FeXs\n7MzkyZOZMGEC77//fol1Y2JiGDNmDMuWLWP48OHMnDmTcePGsX///hKPs6Z6fyUvhKjdrJ3+r1u3\nbkycONGQ1q8kGzZsoF27dowePRo7Oztmz57NsWPHOHfunEX7VB4S5IUQtVZ++r+goCAiIiKIiori\n4YcfLrLumjVr8PDwwNPTEw8PD7P3np6exa4nXx4FUwU6OTnRokWLUtMFWlO9H64RQtRepun/8rND\nlZT+b/z48VbtT3JycqGkMG5ubiQlJVn1vCWRK3khRK1VFen/yqNgqkDQ0gW6urpWU48kyAsharGq\nSP9XHm3btiU0NNRQTklJ4cKFC2VOF2gNMlwjhCi30ORkq82CCU1OJrgMibKhatL/qapKZmYmmZmZ\n5ObmkpGRgU6nw9bWtlDdUaNG8a9//YsffviBYcOGMWfOHIKDg2nVqlW5z2spciUvhCiXYBeXMgdh\na7dfFen/du3ahaOjI8OHDycyMhInJycGDx5s2G+a/s/b25v169fz5ptv4unpyeHDh/nuu+8s2p/y\nqvfp/wo6HjqcjMxoOrU5ZLbdxrXmfemRtGyiKsjfWdWxRvq/mhe5qtitzbfITTOO58WejwXvZPb0\n2GOspIf+2f2rvnNCCFFJ9T7Ihz0XRsaVDOOGudp/ms9vDkDMlhjiQ+KroWdCCFF59T7IA3iP8ibo\nPe1ptnPxDchS0/Dv7Q9AVlyWBHkhRK0lQR6w8bDB+Q5n7f0JPTkZcj9aCFE3SDQTQog6rNQgryjK\n14qiXFcU5Xgx+/spihKvKMrRvNdMy3dTCCFERZRluOa/wCfAyhLq7FJV9X7LdKkGyoE/g/8029R2\nXVucWjhVU4eEEKJsSg3yqqruURQlsJRqFZq/WRs4t3HGa6SXoZwZlUnS4SRy00t/jFoIIaqbpW68\n9lQUJRSIAl5VVfWUhdqtdr4TffGd6Gso31h3g1MP1pkfT4hyC5sWRnKodTM3uQS70HKhpP+zBEvc\neD0CBKiqGgx8Cmy0QJtCiBoqOTTZqkHe2u2X1/z582nfvj1ubm40b96c+fPnl1i/zqX/U1U12eT9\nFkVRPlcUxVNV1dii6s+ePdvwvn///vTv37+yXRBCVDGXYBc6hVgn/d9f/WtW+j+Ab775hg4dOnD+\n/HkGDRpEQEAADz30UKF6lkr/FxISQkhIiEX6XtYgr1DMuLuiKL6qql7Pe98dbT2cIgM8mAf5qpaT\nksOlmeZJurPisqqpN0IIS7hy5Qovvvgiu3fvRlVVxo8fz6JFiyzW/iuvvGJ436pVK0aOHMnevXuL\nDPKm6f9Ai3fe3t6cO3euXCtRFrwAnjNnToX7X2qQVxTlW6A/4KUoSgQwC7ADVFVVlwBjFUV5FsgC\n0oBxFe6NleVm5HJl4RUUewWdvTZSpegVdI7yuIAQtVF++r+BAweyevVqdDodhw8fLrLumjVreO65\n58wWAct/rygKx48fx8/Pr9Rz7t69m2eeeabIfSWl/6uu5YbLMrumxAWYVVX9DPjMYj2qAs3nNcfv\nhdL/MYUQNVtVp/+bNWsWqqryxBNPFLm/Jqb/k2UNhBC1VlWm//v0009ZtWoVe/bsKTJhCEj6PyGE\nsKiqSv+3bNky5s2bx2+//Ubjxo2LrSfp/yopMvI/ZGREG8pJSX+WUFsIYS3JoclWmwWTHJqMS3DN\nSf+3evVqZsyYQUhICIGBJT8XKun/Kun69VVcufIxV68u4erVJWRkRKDXO1d3t4SoV1yCXcochK3d\nflWk/3vrrbeIjY2lW7duhqv+5557zrC/pqf/q1VX8gBeXvfRvv2m6u6GEPVWTXsS1c/Pjx9++MFq\n7V+8eLHE/X///bdZ+Z577uH06dNW60951bogXxVSU8/y1199Tbbo6dQppLq6I4QQFSZBvgBn5w7k\n5KQaymlpF8jIKHxD5taGWyQdNE6Lcu/nLqtSCiFqHAnyBTRrNtesfPHiTCIiPihU7/Ksy2blNqva\nSJAXQtQ4EuTLyWuYF3dG3Gkop19KJ7RfaAlHCCFE9ZEgX056Jz16J72hLOvKCyFqslo1hVIIIUT5\nSJAXQog6TIK8EELUYTImL4Qol7CwaSQnW3eygYtLMC1bLrTqOeoLuZIXQpRLcnKoVYO8tdsvr4UL\nF9K8eXPc3d3x8/Pj5ZdfLnFBtDqX/k9obq67Seo540NU7r3d8RzkWY09EsJ6XFyCrfYU+F9/9bdK\nuxU1cuRIHnvsMTw8PIiPj2fMmDEsWrSIadOmFaprqfR/liRX8hZya+Mtwt8JN7zidsRVd5eEqBeu\nXLnCmDFjaNiwIT4+PrzwwgsWbT8oKAgPDw8AcnJy0Ol0nD9/vsi6pun/7OzsmD17NseOHePcuXMW\n7VN51Okr+cwbmSTsTjCUs5OyLX4Op5ZO9Ff7m23b5bTL4ucRQhRWVen/1qxZwzPPPENSUhI+Pj4s\nWLCgyHq1Mv1fbZZyMoWTY09WdzeEEFZSVen/8o+9cOECK1euxNfXt8h6kv6vmrRe1hrXrsb0W3aN\n7aqxN0IIS6nK9H8AzZs354477uDZZ59l/fr1hfbXxPR/9SLIOzZzxKW9dZIcpObkcLDAPyqoxGVl\nWeV8Qggj0/R/pQX6b7/9lilTpqAoitn2/OGaU6dOFTtcYyorK6vYNebbtm3LihUrDGVJ/1cHRGVk\ncM+xY2bbtuTC+mvX+CPSmLUq0N6esQW+xglRWyUnh1ptFkxyciguLsFlqlsV6f++/vpr7r//fnx8\nfDh16hQffPABQ4cOLbJuTUz/J0G+nI4kJbE42phnNiFbu5k7IyCAez21KZNphJKtwisXLhjq3evh\nIUFe1AllDcCVab+s58hP/zd16lQCAgLQ6XRMmDCh2HH5iti7dy8zZswgJSUFHx8fHnroId555x3D\n/nbt2jFjxgzGjx9vSP/3/PPPM2nSJHr06CHp/2qbS2lpLL16FR9bW2zzvvY1sbOjs6sr/Ro0AOAP\nnY4X/RrzYZ8gAAYWuNIXojaraU+iWjv937Jly0rcL+n/6qjfOnaknUvR4/wKYKfT4Wqj/Xr1BcYA\nhRCiqsjDUEIIUYdJkBdCiDpMhmtKkZGbA6icubEDgOTEBDpwBTWnDWCdaZlCCGEpEuRLcSAxkSBy\nuXbqXgBuAz4GcjJ7AkU/9VaUX+PisP3jD0O5oa0tURacASCEEEWRIF+KaOexfJrYnFf8/AFwSNuP\nR8xHNLS1LXMbjzVqxN15M28AtsTGEpWRYfG+CiFEQRLkS5FmG8gxFMa36A/AjRsKp2I+ws2m5CAf\n+Z9Irnx8BYDbgZ4DPejwcwcAYrOz2XDzpjW7LYQQgAR5q/B/xR81SzWUr628ZlYWQoiqIkHeCoLe\nCTIrx/8RX009EcLypk2bRmiodTM3BQcHs3Bh9T50NWfOHM6fP88333xDZGQkbdu2JSEhodDaNzWd\nBHkhRLmEhoYSGhpKcLB1ljew9gdIeeQHdH9//0KrS9YWEuSFEOUWHBxMSEiIVdru37+/Vdqtr+Rh\nqOpw5AgpycnMf+kls1fKjRvV3TMhap2goCDmz59Px44dcXV15emnn+bGjRsMGzYMNzc3Bg0aREKC\nliHuwIED9O7dGw8PDzp16sQfJtOaL1++TP/+/XF3d2fw4MHcunXLsC88PBydTmdI4L18+XLuuOMO\n3NzcaNGiBUuWLDHU/eOPP/D392fBggX4+vrStGlTli9fXjW/jCKUeiWvKMrXwHDguqqqHYqpswgY\nCqQAj6uqWnO+b1WHDRsgJ8dYjmkAWSazcVJSSGncmFdHjjQ7bFJmJs4IIcprw4YN7Ny5k6ysLIKD\ng/nrr79YtmwZt99+O0OHDmXRokVMnjyZ4cOHs3r1agYPHszOnTsZM2YMZ8+excvLiwkTJtC7d29+\n/fVXDhw4wH333ccDDzxgOIfpWLyvry+//PILt912G7t372bIkCF0797dMIR17do1kpKSiI6OZvv2\n7YwdO5ZRo0bh7u5e5b+bsgzX/Bf4BFhZ1E5FUYYCzVVVbakoSg9gMXCn5bpYC02aBGlpJhs+AU8n\noD8AH/31Fx++8w7kJR5Yun49L/n6Qni4+XFNmoCzhH0hSjN16lS8vb0B6Nu3L76+vnTooF2Tjho1\nip07d2Jvb899993H4MGDARgwYABdu3bll19+oX///hw+fJidO3dia2tL3759GTFiRLHnM11Pvm/f\nvgwaNIjdu3cbgrydnR1vvfUWOp2OoUOH4uLiwtmzZ+nevbu1fgXFKnW4RlXVPUBcCVVGkvcBoKrq\nQcBdUZSyPwpaVz3xBPz9t/ZycobEROjdG3r3xmH9elwyMnCxscHFxgaHvKTCjBgBrVoZX7t3V+/P\nIEQtYZpz1dHRsVA5OTmZ8PBw1q5di6enJ56ennh4eLB3716uXr1KdHQ0Hh4eODo6Go4LDAws9nxb\ntmyhZ887MKGKAAAgAElEQVSeeHl54eHhwZYtW8yGd7y8vMwyVTk5OZGcnGypH7dcLHHjtSkQaVKO\nytt23QJtV7mTKSnEmqTui6jok6ne3pCf8svtDGSq4OSkldu3B5MnYLn9doiPh88+g9xcCAuDOXMq\n+BMIIQpSFIWAgAAeffRRvvzyy0L7IyIiiIuLIy0tzRDoIyIiikwpmJmZydixY1m1ahUjR45Ep9Mx\natQoVLVmPgsjs2sKeP3iRX6KiTHbpq9so0FB4KKH7U8WuVv18YH4eHJHjyLXzg4OHEA3Z46We7Ky\n5xbCCkJDQ602C8Za0zMnTZpE165dGTNmDAMHDiQzM5ODBw/SsmVLAgIC6Nq1K7NmzWLu3LkcPHiQ\nzZs3M9Lkvll+EM/MzCQzMxNvb290Oh1btmxh+/bttG/f3uJ9tgRLBPkowN+k7Je3rUizZ882vO/f\nv3+NnC7VwtGRxSY5GUsMtKqqvQpuK4fdEbuBRjTdvQPUXPS54LZpE+/Fnee5crUkhPVZa368afvl\nOUfBh5OKe1ipadOm/Pjjj7z66quMHz8eGxsbunfvzhdffAHA6tWreeyxx/Dy8qJnz5489thjxMfH\nF2rXxcWFRYsW8eCDD5KZmcmIESPMPgzK0sfShISEWGyKqlKWrxiKotwGbFZVtdBHlaIow4DnVVW9\nT1GUO4GFqqoWeeNVURS1Ml9pDh/uhL19AO3bbypT/bjf4zh2zzGCQ4Jp0K9B6QcAI06cIDojgyNd\nuxa5/8aNdZw69SBdu57AxaWdNrRSVJLeV1+FefMA+KHtDyTqE9nw5gbDbnd7d1aO0u5lv3JwBf8J\nO0Knxp2w09uREpvC356t+Pjq37ww/p/F9lVRlBr7FVHUHfJ3VnWK+13nba/QF/uyTKH8Fm1aiJei\nKBHALMAOUFVVXaKq6i+KogxTFOU82hTKJyrSkVpv9Gjo2NFYNllGODEjkSaXm/DEs9qvJic3hyTn\nJBil7W+hT4MLn/DLA1dp5NKI1T9/wyTgYuxFVh1fZWinU6NOtG3Ytip+GiFEHVFqkFdVdUIZ6hR/\nuVlfjB4NEycWuetMxzPE+cUxpMUQAM5tO4dbnBtBH2tr3CRmFP24tM/2PWQu/chQPjXhcdq++l8L\nd1wIUZfJjdcq8MewP3Cxc2HaI9MACHswDMdfHbkr8C6zeg42DgD4ezaDjCw6Znky+GZjyM3B9uoN\nQu68VuV9F0LUbhLkK+jcuX+g17tCairMgwC7k3iU8djW3q25aX+TFQ+sKHK/3R3t4K+/sPlmNbZe\nXqQlxmLr7mW5zgsh6g0J8uVka+uJm9udgEpOTiI5ufGkdINGsca78L+E/UJOrnFZg9i0WFzsyp8P\nduGVK6y/dYvszAxsXn6ZFq5pec/MCiFE2UiQLycP5z54BP1iKKeG/c6hrDFmdcauHUtadprZtgD3\ngDKfw1ano4mdHSdSUjiRkkKumsu14cMZGPJfvv6kiJk8QghRDAny5fXjj/Dgg8ZyU2BV4WqPdHiE\nF3u8aCi72ruW+RRdXF3NknwnJMbR4OgxvJ288WmiTe3cfG5zubsuREUEBgbWukQZtVVJSylUlAT5\ninrtNWjcGOxuAnO1p1pNNHJpRJcmXSxyKju9HQDBSU68dkK7kv/koB0vWKR1IUp2+fLl6u6CqAQJ\n8hU1aRK0awepYXBorrZipLWdPw/ffQfAVJAgL4QolQT5SopP1xbo/OH0D9w8oy0dnJWbVdIh5eeg\nTa3k//4P1qwB4FJDO7hp4fMIIeqceh/kF0dFkWSS4ON8WhpORaw8V5yEdC3jzPrT69lxY73F+1ea\nQ3cbb8RmN2pIrzV7qrwPQoiaq94H+bkREVwpsJxwZ5fyT3d8psszbO36qaFc2o2qrJtZ7Guyz2xb\n19Cu2DW0K9P5bnk7ws0sfE5cAMAnKZfrnuHl7LUQoq6rU0E+ZmsMN9YY86RmXs0s03GP+Pryhcki\nYxVJfKsoCnpd2RYldu/jjppjXIQo5VQKiXsTUXO1bdHRsHq1sX6WAnSFiAggbybmrf8kwDCFw19o\n30I85/oTcFVyxAohzNWpIJ96OpXrK69jH2CPotOupB1uc0DnUHLYtlUUnPWVXjW+zHwn+uI70Zi5\nJmpxFIl7jevXRETAv/5lcoAdsA2+T47mbGgsAPsua7seekj77/d+0DhLyzti6plnoAp/NCFEDVOn\ngny+bse7YeNe+3+0H36Ae++F+BQFv+/dyHSFUylatvi05kmAll0Q4NwIQJ/JP686mbWxun8C9jbG\nJOKLF0Pr1lXSfSFEDVD7I2Ed5uCg5fHW6XS4v9UZgNS8fbrPDpKLMcOgY7MmNN5/g/37mwNwJe4G\nSeoNlrtmk5trS0wMnDwJSUlV/3MIIaqPBPnSvPUWzJ1rLFdD8gRHRy0FrKlWByHMpNysVRcIu8ad\n57Vhn/TYRBySYVKaiq0D/PSTlidcCFG/SJAvTX56v7feMt/u42OxU5z8Wxt2f2QSpDpAXFwFGvn8\nc7Pigafvpf/SHTT4vwbk6HXk5gIz4XTCXrpimSdxhRA1nwT5stDr4Z13rNZ8bCw0QptVk543pN6l\nC7i5VbxN/7wF0V7s8SK5Nnp2nzzPvvj1fPaZypavjfWefhruvrvi5xFC1GwS5Mtp2/ltTNs2zVBu\noE/h37dbpu3vvgO/DpZpq7lHMwDeH/A+2NryfvxP7ItfT3S09qGSlQWXL8M990iQF6IukyBfTkmZ\nSZy5dYbBzQfj7uCOuy4JiMTLqexJPY4cgYsXjeWICO1KviIWRkYa3nva2vJoo6Jb6tABOAkbNkDX\nJnDlCvj7V/CkQohaQ4J8Bc0fNJ92DduRmhrGoUNbcMvcxZkzTxr2e3oOoWHDh4o8dulSbSpjvhFA\ndyDxx5tcCzX+k3gM8MC+qX2J/XjpwgXD+zZOToWDfEIC2Npik5SCWzqQnV3WH1EIUQdIkK8kRbHB\n3t6f9PRLpKdfAiAjIxJbW+9igzyApyfs2qW9T/sekt+FG2+dx/SZ1fZb2hcb5A936YI7ENe7NwCP\nnznDubS0whXzbhAPARKAz/y/5GiP49rN3S6w+/RdBGw3jjf5+0ObNmX84YUQNZ4E+UpydAyiZ88I\ns227djkVU9vIxsY4xz37FV+yHvc07Ev+K5mTY0+WeLybjfZP18BWe9DJruCiagMHgpOxHxH7thCw\n7lf+G7qcIzeXaxtHwMofv2LlAmOQ/+c/4ZNPSu2+EKKWkCBfA9i42WDjZvynyLxRtjV3StSjh/bK\n06hZIKz7lS0TfyGrc0cux0TTe0U3XnsN7vfT6gwZAsePmy+N0LgxjB5d+e4IIaqHBPl6Ij+7lI+z\nD7g2IVfVlkdo0QJ6aQ/TYmurDSHlDyMB9OkjQV6I2kyCfBXYsAGuXTOW89ebKc7hk4cZylB09+vM\nlsR87733eOWVV6zTSeDsWfMHekeP1rY9+qhxm4MDLFlitS4IISxMgnwVmD8f9u8339awYfH1VVUl\niyyGNR9Gs9bNyM7JZvFPi0mOSiYrq+hsUFkpKVy+fp17/vjDbPu2Pn2wNV2GcskS+Okn3DISmbUf\nvFqEQ96VvLe3eZuBgRAVBXvy8pDExGgfAhLkhag9JMhXkbvugv/9z1guJacIAP3O9KP7me5kkMFi\nFjNn4RzmLJxTZN0zc+aQ9ssv/F5ge1Z6unmQ/+orANyA2cAbjRfzYfpOw+6xd4xles/pAKxaZd7W\n9Ona9E8hRO0hQb6K2NsXf/X++++/s2LFCkP5xjVtImXLL1rStU9XMpIzmNxzMh6DPGhwVwMAFixY\nQGxsLC+++CIA2efO4evry4wZMwD4bONGzv72m/Ekw4ebjcXEnDuGV+tgfJx8cLHTMmH9fvl3ujQu\neV2b/KV8CirLh5YQoupJkK8Ghw4d4vz584byjh07WLFiBX5+fujzrroDAwPxbOOJSzsXHNMcmcQk\nmt3TjIDXtDVpfvjhB2JjY1m5cqWhndatWzN16lQAfoyI4OxvvxEREYGTvXGuva+vL/b29oYndKf3\nnM70R54CwGte6U/tJidDwdmaFy9CUFAFfhFCCKuTIF/Q0qXawut5orf+j4ZqDo3maQPWGTkZxR1Z\nrJSU42RnJ7F3r1b+8MMP2bRpU6F6f/75J42KWZagoMOHD6MoCnGlLFnZxiStIcDBQ4fo3q1b2Tpe\nwODB4O5uLIeGwsaNFWpKCFFFJMgXtGkT/PwzuLoC4JmTQaIdPNzuYbNqno6eRR0NQGoqfPghfPll\n/pZngP306WNyvKcn+wvcjfUueOezElrffTc7TMdQLlyADRt4+NFH8XFzg0xtLv6b33zDyLybuU/s\nTyfp1k7muhjXz/d39+fRjtr0msGDtVe+FSskyAtR00mQL0qnTtoqYsB7v83kgz0fkD3s03I10aQJ\nzJqlvf/6a7Cx6cySJR8Y9tvZ2dGqwFW2JY3v3Zsm+Y/UArv372frlSs0cnPD3caG9IQEQoCbJhPj\n5wOfdD/NC84zDcf1CehjCPLFeeQRs4drmTUL8lZbEEJUs3oV5NNycphnsmojQIKVFuxq0yaU3r21\nuYY//XQdV1d37r33Xqucqyi93d3pbTK28pOzM1sbNWJR5850dXPjSng4/rfdxvYRI0jsot1sVT/4\ngA6NB5E5cx0Ag1YNIju3+N+Pjw/07Am5udpYfUqK9sTss89a92cTQpRd/QryubnMvny5Ss7l4/Mr\n5879CkB6OtjaFj+8Uy3ybvD+b/Nm/rd5s2HzPy9e4l69th6OTtEVeWi+YcO0V75jxyA42PJdFUJU\nXL0K8vkWtmjBC02bWq39yZMv8OSTKjPzRj10OusNy1RUkyZNSEhIMNsW2KABK8+eZUuLFgBEJ0Xj\nfJszPFEdPRRCWEKNDvLx8bvNyjk5yRZpVwEUK07svn59BuvXZ3JJW3mYiIgM2rQpeV34qqbT6XAr\nkF9wnJ0dyU2aGJbHXP9bBKlX45mwfoKhjqONI1+P/JqSxMbC1avGsouL4T62EKKKlSnIK4oyBFiI\ntpLK16qq/l+B/f2ATUB+vqMNqqq+V9nOhYbeVWibk9MdlW3W6nJyvuPiRTsyM7XZMq6uenx8Kh/k\n0y+nk7AvofSKJVgUFUXjmzcN5TE+PnTPC/aLnZ3h0iXyP52uAn8kw9on1gJoi5rZwdcxJQf5p54y\nL8+aBbNnV6rbQogKKjXIK4qiAz4FBgDRwJ+KomxSVfVMgaq7VFW939Id9PWdRKNGjxvKtraWm2Zo\nTR06/INDh+YBcPRoL/R6l0q3Gb04mujF0RU6VqcoOOh0/C8vwKuqSoaq0tLR0RDk+fxzw9RKgOGv\nvEJLnQ5GjgRg/bb13Iq+RVJSklnbzs7O6HQ6/PzMM14BPPNMhborhLCQslzJdwfCVFUNB1AU5Ttg\nJFAwyFtl/MPBoRkeHgOs0XStobPT0WGbeYbvM0+c0T5yy2iYlxdpdxm/GV1JT8f/wAHzSuPGmRVf\nWr4cDh/WlqIEMpPT+W8mhYZ5Ll68SFBQEF5eMGWKeZPPPAPbt2uzb/K1awePP172vgshKq4sQb4p\nYDrv8Apa4C+op6IooUAU8Kqqqqcs0L8a7+bNm8ybN6/AVgsk/TCh6BU8B5nPztE764upXT6X09M5\nYnJl7mNrS4CDg1bo0EHLCZurrT0/OjOb5npYO0H7wIm7EEfkvkhupd4iiKLXNXB21qZVHj+ulVNT\ntS8GEuSFqBqWuvF6BAhQVTVVUZShwEagyCkls00GZ/v370///v0t1IXqER8fz/z587Gzs8PGJv/X\naY8+bxpiTTc3IoK5Ecb0hf9s2pRPWrbUCh99ZFbXd1xfXvxxH6t7aE/IJmYmArBq5SoO+h801Bs6\ndCjNmzcHzK/gATp2tPRPIETdExISQkhIiEXaKkuQjwICTMp+edsMVFVNNnm/RVGUzxVF8VRVNbZg\nY7Pr2B24W7e0//buvYyWLScC2hOuffua18vMvM6NG+sMZZ3OHm/vEVXVzUK8bG35sV07s20TTp8u\n8ZhuTbuB7TFOPa99SXs65mmWspRF8xaZ1Vuzdo0hyAshyq/gBfCcOUUvMV4WZQnyfwItFEUJRJtw\n8TAw3rSCoii+qqpez3vfHVCKCvB1Uf6V6oEDkB8jfXwKTxlMSTnOqVMPGsq2tg3x9r5e6fMfaG4c\nV3ds4UjHbWW7VHbU6xlRYK0cu3JOK+0+uDtLk00WmL8BrIDl65eTFG8cAurcuTNduhiXMN682Xyh\ns8BA43COEMKySg3yqqrmKIryT2A7ximUpxVFmaLtVpcAYxVFeRbIAtKAccW3aDlqboGFzXOr4qxF\ne/xxbXJKUdq0+YacnDRDOTz8HeLj/yi6chl5j/SG+eDeS4uW8bviyYgs/wqZ5ZaUBHkrZT6h5jAx\nx42fVr5FSmMvDh05xOIVi9n2/Ta2fb/N7LDbbrtN62d8Lm3aDGHgQG31tq1btSaFENZRpjF5VVW3\nAq0LbPvS5P1nwGeW7Vrpzjx5husrKn81bOpi3CV08dd4Z9OTABy5esRs/7Fjx1i+fLmhHB5e8lK/\nAI6O5kMXtralr9temuYfNof50OabNgCcfOgkKX+nVLrdEvXsqd05zWNz9iw2ISE81GYMBAXR3rs9\ni6cv5qsRXzG05VAApkyZgpeX8edduXIl8fFLyMwMAbTct5mZTYHfEEJYXo1+4rUs9O56/Kf7m21T\n7Msx7BAVBTk5hmJS/DVy0mLYcXGHYZu/u7H9sLAwFi5ciLOzM3q9Pm/iiVutudFamnU3b3LU5NK6\ng4sLX+Svlvngg9or34oVEBICCQkQF4dDShoNbEHvmIPipv0bLFmzBFc7V1zttfGrRo0aEWFyo/en\nnw6QkXGJhx4y78fy5eYrWwohKqbWB3mbBjbc9vZtFW/gzjvhyhVDsSNwOsCRiJciij8GOHDgAO3a\ntSMsDFq10pqp7fo1aECSyQfekaQkdGUZp+/UCYB2QBww6vIzPNnG+BTUrH6zmN1/NgD/939mD0vT\nvv1jnD27i7//1sr59zWys8F0Ov4bb0Brs++SQoiyqPVB3iL69IEntFW4Pj74MeH6ZBZUc5eqw4YC\ns20GhIaSWVRC13zdusHChYZi2uXzOC78lGe7PsPQu7XAP+WnKcUdDUDnztrN61N5T1UsXQrvvqs9\ng6XTQVoa3LgBkydLkBeiIiTIg3Yp/qQ2Br/D8Qeikyq2dEB1S7uQxp/t/zTb1uVoF3S2JS8ZXJKE\n7Gx2x8cbynpFoVf+1Jg77tBeeRyPHYOFnzIo0QduBQKw/jwc4lsmxJwz1OvapCvTe043lK9fv87I\nvKUTQFuueM2aNTg5ObFzJwwcWOHuC1HvSZCvI9x6uKHmGK+6U8+kknoqtYQjyuZESgp3hYYayq56\nPYkFHwIo6N13DW+3AZ8MieETz8MAXI6/TFp2miHIBwQE0Lp1a8M4/Y0bN4iOjmby5MnY2dkZVrOM\niHiTAvf+hRBlUKeDfExWFocSEw1l0/Hm6pabm86tW+bJvD08BqPXO1SoPf+X/fF/2XiDOHxuOJdm\nXqpUHz9s3pw4k8xZX0RFsb2kxOEtWsC+febbevViao+pTJ06G4COizty7Noxpv4yVdvfE+6/537m\n3K097LF06VLmzp3Lgbx1dRIS0oDrPP/8U7z2Wmtyc1NISppFr17Qvr3xNN27d+ehgndvhRB1O8gf\nT05m2IkTFm3zet6MzXXrYP9+Y7m8cnIS+fvvB8y29ex5Fb2+USV7aDmdCzzR9XNMTMkHODtr0yxL\n4G7vzpXEK3z797cAJGcm08S1iSHIP/XUUzxlslbxsmU7mTx5IH36aNPzExLSWbfuP+zcacvevXYA\npKSkMHnyZAnyQhShTgf5fJ+3bGkWsALty762+9WrV7luEskPHNCujivxlDEBAW/SuPE/DOWbN/9H\nRMS/K95gLbLriV1m5cc2Psau8F3F1NaehgWIj38Dvd6LlBRt8Tdf3//Qq5f2beDnn/0oZUUGIeqt\nehHk73B2pkeB5XHL6rPPPmPu3LmFtn/1FQwdaix7eJS9TQcHfxwcjEMriYkHS6hdOX+P+htFZ5wG\n2fw/zXFqWXsmoLu4uBAcHExKSgopKSnk5ICDQzCK4m2YkZOWBpcuXWbt2h8MxzVo0IBBg+6upl4L\nUXPUiyBfWYqisGHDBkAbcv7wQ2jXLgArpomtNLtGdrh0ciEzWrvyzU7IJv1iOoFvB1aq3fTcXF46\nf95s28zAQLxsS3gYLDHR+CzCjh3g6GjY1f3PS9gkFn+DuEePHvz1118l9kmvh6tXdzJu3E6TrV1Y\nvvywoaQo8OijJTYjRJ0kQb4MdDodDzygjZ/n34d0qXyiJ6tqPLkxjSc3NpRv/XSLv0f8Xak2HXQ6\nHHU6luVNeclUVdJzc3mhadOSg/xHHxVatjjf80Ckpw0s0cpXEq9wKa7wDeMOvh1wd3AvtB3g9de3\nk5SUZSh/8slU4DiPPz7GsE1R9Dz66NqSf0Ah6qB6H+SvpVzn97+/45UFWwGISY2hbcO21dKXjIwI\ncnONC5nZ2fmi11t2aCX5aDK5qcaV3BxbOmLfuGz3KN5v1oz3mzUzlFdcu8bjZ87w3Y0b+JgE+YEe\nHtyWf7X+1VfmjagqpKcbJr8ffWwQnmFRDF2tjX1tPb+1yHPvenwXfQOLnro5d6553t+kpPbs2xcH\naHPzIyKukZ4eX8SRQtR99T7I56q52OpsGNZimGFbgHtACUdYz9GjPczK7dtvwctriEXPcW7KObNy\nq69a0eSpJpVq881L5lfeG9q2NQb5glm9C9B5eGGnv0FsmrYydbcm3QiLDWPFAytwtnXm6NWj/GvH\nv8rVn//+13ytvN69Z7Bv37/5vMAyoQ8//DCenuYZt4Soa+pXkM/JgfBws036HBVfZ1++ul+74vzP\nf/7DX3/9xaSvJgEQavIgkLV4eNzN7bcvN5TT0s4THv6eRc/hdqcbHXca15rPupnFqYcrl6FxrI8P\nAxo0MJRPpaYyuJwLwwc36giuCRx8yno3nzUqzz//vNkWT09PgoK0tIVpaWl4eXnhbbLGvrOzc6F8\ntkLUNvUryCckQIGMRb7kf6nXhISEsH37dvz9tdkv2dnQpElz9uzR9ltjqp6TU2ucnIxPcyYkHLB4\nkLfztsPuHjtDOf1KeqXbdNbrcdYbc83G5N2wiM7MJMxkSWJPW9uSx+zj46HAwmVMmwblmOpakl69\n3mDfvhfzl8EnPX0T8fH/YPz48SUeN378eLMF1RwcHPDx8bFIn4SoKvUryOebOBEGDQJg2tZpNOzY\nCtPR3nbt2nHkiLaO/MyZMHdu4XR+dUV2fDYZV43JRvROemzcK/dn8c+wMLPyrMBAZgcVnegb0IL8\n66+bb/PyAicnfG/8zfjjsHXfSkKvGb9V3RV4Fx0blS0LVs+eLkyebLxTfvr0cPbt+4U339QesLp2\nLZ4dO0IJDDTeb/jf/55hzZo1rFmzxrBt8ODBbN1a9D0DIWqq+hnku3c3zKdbHzODQc18S6yuKLB9\nu/m2vERHtd7FVy9y8dWLhnLTfzal5SctK9RWgL09q9q0Mds26fRpvr95kxMpxoQm3d3ceC0g777H\nV1/B4sXGAxYt0gL+008D2vLF3wJ93Zeyx2T258y+M3GyNd6UdrN3w9el6H/H0aO1V761axuzb19j\n3n/ftNZ4Dh0yLdszZEgOY8dqpaKelRCiNqifQb6cdLq6txKiTQMbWi1uZbbt/PTzxdQuGw9bWyb6\nmgfaD/IWHjuXps0aCktNNc/SaGdnVp9//ANGjTIUs/fsxmbyU/w0fjPZvXsSmxZLq09b8d7u93hv\nt3FIa3KnySy9fyll0bs3bDJfNoiMDLj9dvDz08qNGz9Ohw7aEscAXxWcJSRELSFBvgaLilrErVsb\nDWUfn1F4eg62SNs2LjY0mWI+q+bimxeLqV1xJ7p1Myt3/PPPYmrm8fAwe3zYJjISAPeQ/XA1Frec\nLPbaTCGmbTMSm2lPo03dMpVj14/xwZ4PDMc1cmnE48GPF3mKpk0p9UG2onKlbNu2Ddu8ewvZ2dk0\natSIyLz+nT59moMHD6LTmS/rPHjwYJrW5KfmRJ0nQb4G0ulssbX1JSnpKElJRwGVrKwbODoGWSzI\nFyfu9zjOPnPWUHZs4UjAK5adUrorPp5+Jk+xNrKz4/u2pTybkDe2Ygv0yt+Wl0WkR2wS391+mDcG\nGJ9w7dK4S7FBvqyuX4f8yVUDBjxO584DDJ8/77//PteuXTME/eJs27ZNgryoVhLkayBX1y707n3N\nUM7JSWP3buuvN2PX0I6sW1nc2ngLgOzYbNx6uFk0yHdzdcXTJDCeTkkhIiOj+AN69YILF4zlpCSY\nPh1MZrk03xTO6y0m8cqMLwAY9f0otp3fhsv7xputNjob4l8v3wNRK1ZoL42WznDGDK00blxvYmOP\nctddxvoxMTE8/fTTuLi4cPjwYcaMGcOQIUNQFIVcLRkwwcHBvPrqq4ZjAgMD6d27d7n6JUR51Psg\nH/9DPH/o/uCRdY8AcPToURo1qjnL/ZqKivq80Br0ltT9dHezcuiAUJIOJ/FXP+NVt95JT4ctHSp8\njqW3325Wfuz0aTbHxPCJSZ5dgClNmmCn02nr3Jg8ZQvAzp1mRcXPD5tDf2Lz6hsAfHgjm2dd7mLX\nmK4A7I7YzaGoQzy87mGz4+YPmo+fm1+R/fz+e/KStGvybxN8kDcilJMzDBjG/v3mx73xBjRsqL2f\nOXOmYXtqaioLFiwgNDSUiRMnGrZPmDBBgrywqnoV5FOzUnECfjr3E/t3assHpxxPITwjHDVSy6qU\nne1ASoof77yjHbOr+FVwq4yi6PDwuNdsW0LCHquf16WjC6Z3SVPDUslJtHzilbjsbF4osOjZY40a\naUG+LFxdITISli0DoF1KCu2Cg7n/y/kAfLT/IxLSEwxTMBMzErmafJW3+71dbJP3329eLpjqNj+f\nQHynCVoAAB7ESURBVL7QUPjtN2O9gIAA3jXJkJWbm8uUKeb5bu+991527NjBPffcQ25uLrm5uaiq\nanh/4MABXFxccHZ2RjG5SbB161Y6dizb9FEh6lSQT8vJIcUk+1OCSVYjgIzsDJyAXy/8yuf7fgdA\nVVWa9WnG2V+1cegRI+Cnn2DWLONxJs/7VAudzp6OHc3ncB482AoI48IF4yP/trYNCQh4xWLnbbGg\nhVn5/PTzXF161WLtA3zWsiULWhjP83lUFG9fvkzjffswvfe5rWNHersXvUBZoSfUhg+HiAgt8AMv\n+Y3lJb+xkPeA29qTaxm3bhwHrhzgWrJxWOwOnzto5FK2b3Fjx2KYXgnwxRdakO/WTft7yc2F6Ggt\nP7z2jIUOaEXnzpB/+6FHjx5ER0eTnZ2NTqdDr9ej0+kMryFDhvDbb78xbtw49Ho94eHhbNu2jewC\nf9dClKROBfnF0dFMNx2/LcboNqP5+K31APgv8adPQB+z/Z06weHDRR1Zc+h0WprAqKhPAcjNzcDJ\nqbVFg3xRcjNyubLIfGil0WONKvwAlYuNDaYLet7VoAEv+xmHUMIzMlh38ya5BS+lS3PiBAQUuJfQ\nrh0AgzISOJ4AD8RO5qLJ0jWrR69mQvsJ5fwJNK1amS9lvH+/9rT0V18VXqPNuFzOWsaOhS+/LNs5\nNm/ezLZt29i9ezfXrmkfTqqq0r59ewIK/KxKUdODRL1Up4J8vnnNmuFo8lW/hcn65WWhKNrc+Jqs\nW7fjgMJdd2nLB5w8+RApKZVbSrgs1EyV8y+aD614jfCq9FOy+fo1aEA/k/VwdsbFse7mTaaGhdHA\nRjtHYk4ODWxsGGWyzswdzs4MyJ/68sILZnPt+f57bT5+3jIJzhGXaR8Wyf8Gfklqi0CuJFxh8uan\n+McPT/Lsz88CkJKZgp+bn9mQTlPXpgxuUfTspgEDtFe+9HSIijKWMzO1b4emjxGsXAnJyeX69QDw\n0ksvlVpnyBDjwnbt27dn3rx55T+RqBPqZJD/R5MmuNvYQFaWlrAiORmSk1FKSkItShX0XhCBM42P\nnd5Yc4Owf4aRsCeBtAvGJZJdOrhg19CuqCbKrYGNDf1MhmlisrP5O+/p2d/jjbNlJjdqZAzyeUtW\nGHdONivarl0L48bR+X7jGPnDwDf/upuj3bUbwwsPLiQ8IZzJP5of+/kw40qWjraOxU7TdHAotEwS\nawssZ791K2zZog3xgDae37y59plUlL59+3LI5LHchIQEvvvuO8M6SwCzZ8+mW7duxMZqq3qeOnWK\nVJN1hDZu3MiuIm40zZw5U1bkrKNqVZC/OOMi6eHGhbUS9yWWfMCePXDPPYaiHfAGELHnb7a+oc3E\niI+vO+uMq2ouOTnmWZZ0OkeLfXXXO+nROxlvUOhdtPdnHj1jXs9Fj+d9xoDhM9aHhmMbVuicXVxd\nCenUyVDOVVXiCoxJtyvtAauCOnQA0yvbmzfhww95JNaPR65oA+b//vVOMrsGG6qsOr6KI67JPMdz\nhm2+zr6Vmovfv782bp/vl1/gyBE4eVIrqyqcOqV9KdG+xDQAuvHBB8YJRwMLPIo9y/RmEnD33Xdz\n6NAhOnTQZkSdyEts75qX8zgzM5OMjAw8PDzw8vIqtq/dunWjW4EH20TtUKuCfOwvsaSeS8Wucd5V\nog4c/B1KP3DqVGjZkpiEm3zw1rsoh8OwPbbAsNvGplb9GoqVlnaW3budzbbddVcmilLyAzsV5TnE\nk+DdxkCYeiaVCy9dwK6xHcmh2jhE2tk0Mq9lknrS+OHjdLsTDcdVLOjrFKXQipZ6YNm1a3xjknC9\ns6sr+zt3LrqR22/XXvnOndNyOi5dqr0AB8Bh7wHtkhx4Lj2djIF38/7L2oJlr/76Kmv+XkOPpVoO\nAFVVOXnzJB8M+AAvJ2OwvPu2u2nsaszQZargWP38+XDggLG8fr2WyPz0aW0IMSlJy6LYqpX52kl9\n+xqeCyukb9++eJg8QdyiRQucnJxYtWoVAMuXL+eJJ57g7beLn2kExm8IovZR1PLe0KrMyRRFLc/5\nQkIUAgPfJihoDgCHOx3GPsCe9pvaF1n/o8hIpl+4QHyfPtpwze+/a1fyISHQrx8Xoi7Qwq8Fo18a\nzfoF64tsY8QI7eoqbxHKGk1RFPJ/n7dubSI11fikalzcTuLitqPXu5tdyXfs+Buurp0KtWUtIUpI\noW3eD3jT7od2FjvHv8PDSTSZVbX2xg3isrMZajL8kJKTw9cF5ui76/XY6HTasJ7JB4SBqyvkDxXd\neef/t3fu4VVVVwL/rXPzIA8IJOGhRaE8lAoooqJIO6BSRGVqxerU6XyFdupM1aqdGa0zvuvXqaXW\nqYpl8FWx7Wd1tFYBtRUt+EClogQQIcgrQAiEkAd53Jvc3LPmj33OfeUdSXJN9+/77nfO2Wffc9dd\nyV1777XXXttYWM/XvaV8C8XhMh77J5OQra0drZ674jnOGx3bUDw7PZus9K7NEfmsWNEytBOMn3/h\nQnMeDkNhYWK97GxoKwloMBiktra2zc9UVUaMGMHMmTOZNWsWAD/+8Y8ZN25cNIWD67rs2LEjYZHX\n7t27uf766xkwYABlZWXs3LmzRZrms846iwlJfxNL63i/9W4NyftHF/Yz8OyzEJcgkb17/aHx54vC\nwksTrnNzp5CdHfsBhUK7OXJkBQcOPEJmZmyZfUHBvB41+jPdmQnX609bT9XqKtafGQtfyjw+k8nL\nW2+4O8N/jUrcnPxIOMzq6mr+6hmvHV5ytJfWrk2oN3rAACZkx1YS/+iEEzgvrtebwNChxsh7qYYn\nVlUxUYT506cDEIycTn1THVXXfAfNNTtaXfWHq7jiuSsSHrNo9iJ+NKNrO135zJ4djQoFzPm555r2\n6YEHTJm/ePjWWxPf66vIdU275XltyMrKIqudwAS/E/Hmm2/y5ptvRst37NjBN79pFpc988wziAgP\nP2wivYKevp9//vl2v88tt9ySkNM/Pz8/Or9QWlqaMJfgM3589zKk/i3zOTRnx5abbjK/3XjaGuV/\nnsjPn0N+fmwCsqJiJUeOrKCsLDFeLyPjuB418snzAYNnDia4OzZJW1dUR+P+RqreSJwUz5uZh5PW\nvRCnR5N8F+/V1LA+rrf6elUVm+rrGZqeTmU4zK5QiIpwmD9VVjIzL48IEFHFVY2eR+68k+mDBrHU\nf/b3v29iHz0feJb3KizaDsOHMzIc5JOa2AriCC4fHFgPL97Ck9wCgAJNZ53Ogl+Zxmd58XL2Hd2H\nxK0QaIw0cuPZN+KIA2lQMBwy0zJxxGHkyJaLtF580fTmfZYsMZGk/p/hqadM45CcKvudd2IZOOOJ\nHy22RXzOfYBdu3axbVviPE11dTXjx48nPz+fyspKpk2bxqJFixI2ZRk1ahSPePGkl112WbSxiOeX\n3obwqhqVK/583rx5dnSQRL828iXVJYwCfvbOzyiuXsbRqtYnaq+80rhkfTrIOfW5pKDgEmbOjE1Y\nNjbu5/33R1NV9TquGxvK5OScxpAhs3pMjuRc9VsXbOXQbw6xcfbGhPL0oekEcgOoa368E5+dyMCz\nBsYqSOdjwafn5TE9LkLn+iRrtjcU4puffIKDMbzpIgxwHAKYOYCACO/U1HBcKG43rSVLzMvnd78z\nBn/bNiguJhv4UpxxVFVGhvK8TzDkHj7KX45u4ILvm9HEz16HkxvBTfpaeVffRiRuQd69F9zL1OOm\nGuOGMXCuuihK4EvKuCFfZNIw4w67InEgwZw5iXsj7NgBa9caN2Vmpmk0KipgwQIzMeyHE/tH/zw/\nv2X0kM+YMWMYk5yKIo6RI0fywgsvJJTNnz+fkpKShNBPo1Yzd3D77bezZ8+eDsNHXdcl7LVyruuS\nn58fTVPywx/+kF27EjOtOo7DU7EERYAZUTiOQ3NzM3V1dQkL1eKPn5e1CP3aJ//8z3/EH265j7fH\npVM1MB1tVoKbg1x181U8/fOnAbMIcs4ceOKJz/79epvO9LLaIhTaz/vvn9CifMSIhYweHVuO7zgZ\nZGR0b5K0MzQUN9BU3hS9rny1koZtDdHInerV1TTub5nALL0wneOviaVKzp2ay9Cv99zWfOd8+CHF\nwSCTc8zEtgIloRBPnHwygvlbCMRe3vXknByGtNFrqDlxOHn7yluUN19yEQChrZvJ3bWfn6/+bzQt\nwIrtK1i7b22L+skMzR7K/1wYCywYljOMOWPntFp31Sq4//7Y9aZNUNbJRc2+58TfCMx3c/rBT8uX\nm2mN/fth9+7ENQKhkJks9tansWnThzhOE/Fb6ubm5jJ5svmtNzQ00Oj5okQkamD946effsqZZ57Z\nKbmnTTMjrL8m7hIT5eqrr2bYsGGsXbuWNWvWdOqZixcvbvf+tdde2yINdVf4LD75lDHy9fVbKCo6\nL6EsHD7crpHfUFvLvrgMhi8dPsyvDx2i+swzyUtL4/G7b+TqRY8yPH8weYUxQ3XbbbfxbW954t+q\nkVd1aW6uSShbt24Mzc0tQ0rT0sziJNcN47r1nHFG4nLgnJxJOM6x2Y81meDuIId+G5sUDR8JU/qQ\nt8rI/5f3VJD3lVhvPfvkbE5+rI2Qk25w3fbtbI3zEa/uZOjt4LQ0zvMWd9VHIpSHw1zsTwjv28dp\n4TBXxmdCy801jnaAn/wE7rjDhHwGAjSL0uA2geOgAQcCAQJ19QQnjKPm/HMRhG2/vI2ipMwMewZD\nwQ23RK93Vu3k+mnXR6+rglWcMvQUstKzaKgXPt7sAILrCuqaZiuDHDIkG1W4995ohogob78dW/G7\neLGJBOoq3/tey4ijzlJbW8trSdu3PfPMM0yZMiWhbOzYsdG5hNWrV7PFj1cFli5dGr0OBAJEvMn8\nuXPncv755xOJRHBdN3p0XZd7/CRXHVBdXU0gEGD37t28/vrrZGRktNpY+WWTJk3i3HOjSbV73siL\nyFzgAUwCjidUdVErdR4CLgLqgYWqWtRKnTaNfF3dJtavP40hQy4kKys21MvPv5jCwnlASyO/cOtW\nnmolKqJ63jzy6utZDZwPPHbHD/jePYsJheDRRxPr3nknXH75356Rb42DB3+L68YazfLy35OTMxHf\nmpaWPtTq+44//jqysmLhGwUFXyM7u/cmyLZcuYVwRcwRXbexDgQKL42tiA0fDjPuAZMjp2FrAzXv\n1kRHCwBuyGXI+UPInhibiA1kBQjktJ64aFcwyMGmJhTPJ4zJ5eafb6ir479LSjg+MxMBgq4bnQBO\n837Qzd7fbqCXHCmsSsh1WeB3eUtKYNcunli1ikBzM0Qi5uW65rh3b2Ia5jjU61aL162u8drg9Ahk\nN8Ot3vKRL+8FFaiMm3stbIDXxsKHcZGfJYPh7HON70dRhucM56ZzTQqNldtXUttYa+YNgOoaKC2v\n49LjryXda/zr6mDSSQMpHJJOKGQSvMXPcd98s3EF+dkZNm40CUj9QVA4DMFgoocsnvJyk4okfvfJ\n5maTomTePJNPyDHtYvTlu54+C5WVldE00q3xi1/8ImHeobPccMMNPPjgg9HrHjXyIuIA24ELgAPA\nB8A3VXVbXJ2LgB+o6iUicjbwoKqe08qzOjTyEyf+gaFD57dapzUj/+eqKl72hnTb77uXA48vY/C0\nU5H0NNbt28Uj67byq8V3ce0P7qay0uwPncx3v5to5NesWRMNF0tljrWR74j6+i0EgzGf5sGDT1FR\n0Xooqoj5cas2UlQEM2ZM8q5dVCOMHRv7x3ecbPLzv9rqc7rD1oVbqX4j1tNuzd3TGkUUMYXEnt8J\nN5kuq6oSORphxIIRRudK9BV/nfflPJwMB3UVbW75t5GAIAHzW/1pSQkVcbOkv9y/n5GZmfjNSok3\nSj0jN5c0ERTY8u67TJ4xw0QFuS5VNTVEVPmKt43iU8EgM3JzTdwkwJYtlDY1cc2ePeb6pZc4PHgw\n523YQFokQsB1CbguBITGIQMZsr+cgOuSGQ7jxP1v1WQ5PDK3kIgjHGg4RLMDEYGIA80OOAqvjDdH\nR6FpL2SPhIAbK/v6SX/PGcNPR1wl1NTAoNEnM36C6a3efgc0VA4mKxyL/IofHTzwxAHIaenaomIC\nNCeulfGDhVqZt22FNVx44SwA/vxn41IaONA0AK5rXE0nnwxeEBX19eb6q3H/rgUFbc9PvPfee6xN\niuqqq6vj8ssvZ5iXkzp+AhlgwoQJLFy48JgZ+c5MvE4DPlXVEu/DngEuBeKnzy8FfuMJuk5E8kRk\nuKq2EnzceXbftTvhurGskcwTE90CmSJM9VbvPbFmBUvKyuClRKfihvIhrFwZG0L+9KcQn/U1eZvR\nz4uR721yciZ6PXtDfv5cXDc2GdnUVEZp6ZJo8jSAsrJH2L59BLNnm/1kKyrMhNvHH3894dkFBV/D\n9CcEcLzz+GOs3HWbyM+fw+DBsWQxgUA26en5qCoTnjzJ+8G4qLrUbTxK/ZY6HI310jWsFMwrIG1w\nGuGqMOW/L2f5quV842KTWnLH9TuQNKF0iXENuQ2mt1b2WPsOa8kQ0gvSaSprarPO0CvM3MHXXKVx\nXyPH/bPpMt/ESUTqItERyPLDFfzx8GHEixf4qK6OtL98yMFJUznQ2MjBpiaaPMOwsin2eZ/U1HCc\nZ+HKPNfQK/4eCee06Hu1SoYqc/DarnXrzFEE13FwRSgaN47vvfxytP4HEyZwjxdi+dapp1K8/w1m\nFcYazLKCAiat/ICRB15AVHFcF0dfpjTzftIiyg1AesTlxPJyAi5kRiAtEuHFd5pQUR4Z34xv4vym\nZ1wlrJgNFYPSKKht5qJP4dN8KMg23zlS18y+wRlUTPwuTsRlc9P/MUzHIa4SiCjDNlfxSu1Ovli9\nEIB/PRuCNfuoGHslWQMHUry/lLxR6zgaGcZrxRBpVprcEC98dCq33++7lgVUuPhi49Z65WVzbUYg\nQrjpRGAUCxYIAUeikVNLH8snQDoHDphFb6ecEhtV1NU5dHURd3t0xsh/AYiLzmU/xvC3V6fUK/tM\nRr7knpKE6xAhiqsbWextHbdy8UMEt2/nyo3GM1TcYP6xr7nzJrKGHsfevfD8c/D4Ly7l8bhWPScn\nPhOgpbs4TjqOE5tUTEsbyPjxDyTUGTv2Z7z++t1MmnQ3APX1WxMahj177iYUKqGxcW+CYW7taLZB\nrCQSOUp5+dNdE/YEM3dgPI7G71m6T2Cf+XE2Tz9C7Sd1BC+MAMLIT/0GRxBxaDzQhDYqiGN+qNGj\naXwaS5oIhreTdXgWEs4mB6HpcJi0gQEGnJiF2+hy+I+lZBYOpvbQFwAhtMPoofbHXn/Js2I77zPH\n4SreflQxltXWsfDd/V64i3htoiCOkUXEO3dM7y+iih5xCVw8CMShKhI2bhVHUDUf6YpvwKE6HGF3\nY5BB6emogCKkhU9j2G6XkinGHXOgOcLYT6Bavg4CYYRR22DdiFNwHUg/DDkVu9g4dD5OmhJyHFTg\n8b+bZZ4pxBlsiX71aFlcnWHVR0zIj0JU3YC4LiogaoZQNcCvJzioKLWDIqau6yKAinmzw0QqxI1+\nRtkYpWTDW1xQOCRatm3smcxe/y659UFmFIBE0tFABa4TW2ynw9dSP2ktKPxp+gymFhdHp4du+Af4\neMxYTtu8JlrfzTmBzNIPQJWKwpEcGvoFBtUdRVQYmudw++XDqRgY4cSDJpb7nuIGmspj7/+s9HoI\nZXZ2ByOOpquIjxdrToOwb0eCQXgLSIpjf9vviafBCGDZvVcTDJ8Uvb90KZxxRqx+8qSRpffIyflS\nwvXkyV3b6SoSCXHo0O8Syg4ffp5Bg6bhG+SYITfHioo/kp4+DMfJAPyYar9BUUKhXYRCe2huzuDQ\nod94jYp6dePO090WZdF+5WhzaBwf6936RKe3p0DnHEftsAxY2NJFpknHZHwTNaiN+z4FQBueB7qy\nCeSyRlh4/u87rtjHLNsJl8Vll54PLbuw7TCTZSR5+LjMFegoA0Q7npeT9qQz9vg2b3eZzvjkzwHu\nVtW53vV/Aho/+SoiS4HVqvqsd70NmJnsrhGR3nMgWywWSz+iJ33yHwDjRGQUUIbJynpVUp3lwHXA\ns16jUN2aP767QlosFoule3Ro5FU1IiI/AF4jFkK5VUT+1dzWR1X1FRG5WER2YEIov9OzYlssFoul\nM/TqYiiLxWKx9C49ssmdiMwVkW0isl1EbmmjzkMi8qmIFInIlNbq9DQdySkiM0WkWkQ+8l6394GM\nT4jIIRHZ1E6dVNBlu3KmiC5HishfRGSLiGwWkRvaqNen+uyMnCmiz0wRWSciGzw572qjXl/rs0M5\nU0GfnhyO9/nL27jfdV36GdyO1QvTcOwARgHpQBEwIanORcDL3vnZwPvHWo5jJOdMYHlvy5Ykw5cx\n8/eb2rjf57rspJypoMsRwBTvPBcoTtH/zc7I2ef69OTI9o4B4H1gWqrps5Nypoo+/w34XWuydFeX\nPdGTjy6eUtUw4C+eiidh8RSQJyLD6V06IyfEMqT0Car6DtDe5rSpoMvOyAl9r8uD6qXbUNU6YCtm\nPUc8fa7PTsoJfaxPAFX1E/pkYub4kv2/fa5P77M7khP6WJ8iMhK4GHi8jSrd0mVPGPnWFk8l/4O2\ntXiqN+mMnADTvaHRyyJySu+I1iVSQZedJWV0KSKjMSOPdUm3Ukqf7cgJKaBPz72wATgIrFLV5LWa\nKaHPTsgJfa/PXwI30/Zyh27pskd88v2ID4ETVXUK8DDwYh/L83kmZXQpIrnA88CNXk85JelAzpTQ\np6q6qno6MBI4u68b77bohJx9qk8RuQQ45I3g/IzVx4SeMPKlJC6OG+mVJdc5oYM6PU2HcqpqnT/M\nU9VXgXQRSbWECKmgyw5JFV2KSBrGcP5WVVtbbpsS+uxIzlTRZ5w8R4HVwNykWymhT5+25EwBfc4A\nviYiu4DfA+eJyG+S6nRLlz1h5KOLp0QkA7N4KnmmeDnwbYiuqG118VQP06Gc8f4uEZmGCTmt7F0x\nzcfTdsueCrr0aVPOFNLlr4FPVPXBNu6nij7blTMV9CkihSKS551nAV8lMXEhpIA+OyNnX+tTVW9V\n1RNVdQzGFv1FVb+dVK1bujzmuWv0c7J4qjNyAt8QkWuAMBAE/qG35RSRp4FZQIGI7AXuAjJIIV12\nRk5SQ5czgG8Bmz3/rAK3YiKsUkafnZGTFNAncBzwlMQSBj3r6S+lfuudkZPU0GcLjoUu7WIoi8Vi\n6cfYiVeLxWLpx1gjb7FYLP0Ya+QtFoulH2ONvMVisfRjrJG3WCyWfow18haLxdKPsUbe0iuIiCsi\n98Vd/4eI3NnLMjwpIvO988dEZMJnfN4oEdncRnmDlzL2YxFZJiKB1p6R9J7kHdcsls+MNfKW3qIR\nmN/dpeIdGcmuoqpXq2ry6sxuPaqN8h2qOhU4FbMU/coOnvNF4B87qGOxdBlr5C29RTPwKPDvyTe8\nXuwbXgbAVV7KVb/n/b8i8h6wSETu8nrFb4nIbhG5TEQWicgmEXnFbwhE5A4xm0RsErPJfAtEZLWI\nTBWRvxezmcRHYjaQ2endP0NE1ojIByLyqr/s3Ssv8laiXtfRl1ZVF/grXrZA77u+JSLrvdc5XtV7\ngS97ctwoJmviz73vUSQiV3dN3RaLwRp5S2+hwK+Ab4nIwKR7i4EnvQyAT3vXPl9Q1emqepN3PQaT\nPuFSzOYKb6jqqUAIuMR/nqqe7ZVni8nw17pQqitU9XSv170RuE9McrCHgMtV9SzgSeCn3lt+DVzn\nZTRsDwEQkQGYDR7+5JUfAmar6pmYHCX+d/1P4G1Vnerlq/lnTG6SszF7H/yLiIzq4DMtlhYc89w1\nFktbqGqdiDwF3IjJD+IzHbjMO/8tsCju3nNJj3lVVV3PF+6o6mte+WZgtHd+gYjcDGQDQ4CPgZfb\nk01EfgQ0qOpSEZkITAJWiYhgOkMHvCRXeaq6Nk7W5KyLPmNF5CNMo7RSVT/2yjOAh8Vs3RYBxrfx\n/jnAZBG5wrse5NUtae97WCzJWCNv6W0eBD7C9I592kugVJ903QgmY5OIhOPKXSBNRDIxI4apqnpA\nzH6eA9oTSERmA5cDX/GLgI9VdUZSvbz2npPEDlWdKiIFwFoRmaeqKzHbux1U1VM991KwjfcLcL2q\nrurCZ1osLbDuGktvIQCqWgX8H8Yd4fMu4EeW/BPwdleemcQATKNxRMymG99o9wHGBfIwcIWqNnnF\nxcBQ318uImkicoqq1gDVInKuV+9bHcmmqkcwrphbvfI8oMw7/zZmz1GAWiDejfVn4FrPdYSIjBeT\nJtdi6RLWyFt6i/je+v1AQVzZDcB3RKQIYzhvbOU9HT3TFBhD/DiwBXgVM+nZWn3/fAGQD7zoTcCu\n9Pb8vQIz2VsEbMC4lAC+CyzxXDGdkk1VXwSyvBTCS4CF3sTtScRGKpsA15PhRlV9DPgE+MhzTS3F\njrwt3cCmGrZYLJZ+jO3JWywWSz/GGnmLxWLpx1gjb7FYLP0Ya+QtFoulH2ONvMVisfRjrJG3WCyW\nfow18haLxdKPsUbeYrFY+jH/D7KVju+SIODeAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "kws = dict(bins=np.arange(0, 4, 0.05), histtype='step', lw=1.5, normed=True)\n", "print('{:>8} {:>6} {:>8} {:>8} {:>8} {:>8} {:>8}'\n", " .format('c', 'n', 'Mean', 'Std', 'MSE', 'Median', '% > λ'))\n", "for c in (-1, 0, 1/3, 1, 2, 3):\n", " r_hc = (n - c) / d1.sum(axis=1) / λ\n", " r_hc_err_rms = np.sqrt(np.mean(((r_hc - 1)**2)))\n", " plt.hist(r_hc, **kws, label = 'c = %.1f' % c);\n", " print('%8.2f %6d %8.5f %8.2f%% %8.3f%% %8.3f %8.2f%%' % \n", " (c, n, r_hc.mean(), r_hc.std()*100, r_hc_err_rms*100, np.median(r_hc),\n", " (r_hc > 1).sum() * 100 / num_iter))\n", " \n", "r_hm = 1/np.median(d1, axis=1) / λ\n", "r_hm_err_rms = np.sqrt(np.mean(((r_hm - 1)**2)))\n", "plt.hist(r_hm, **kws, label = 'median');\n", "print('%8s %6d %8.5f %8.2f%% %8.3f%% %8.3f %8.2f%%' % \n", " ('median', n, r_hm.mean(), r_hm.std()*100, r_hm_err_rms*100, np.median(r_hm),\n", " (r_hm > 1).sum() * 100 / num_iter))\n", "plt.xlabel('Normalized Rate')\n", "plt.axvline(1, color='k');\n", "plt.legend()\n", "plt.text(0.35, 0.75, r'$\\Lambda_{n,c} = \\frac{n - c}{T_n}$',\n", " fontsize=24, transform=plt.gcf().transFigure);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can observe that:\n", "\n", "- with $c = 2$ we obtain the mean square error (MSE) in estimating the rate,\n", "- with $c = 1$ we obtained the unbiased estimator of the rate (but higher MSE error),\n", "- with $c = 0$, we obtain the MLE estimator,\n", "- with $c = 1/3$, the estimator has median equal to the true rate ($\\lambda$)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "By contrast, when estimating $\\tau = \\lambda^{-1}$,\n", "the MVUE is the empirical mean of the delays ($\\hat{\\tau}$).\n", "The distribution of $\\hat{\\tau}$ is shown below:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " n Mean Std Err.RMS Median % > 1/λ\n", " 5 1.00177 44.74% 44.742% 0.936 44.33%\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAEPCAYAAACneLThAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHbdJREFUeJzt3XmUFOW9//H3dwZRQYIKXheQISaKqFHcEH8aGTQq7ite\ngZi45XijRI3iVWNUvNHkovG4JlFUXPAal2iuYFDhXJ0YFwQXxCCIIouIonFDRHSE7++PqqYXZunp\nqe6qrvm8zqlDP0/XVH/n0flMTS1PmbsjIiLpVBN3ASIiUj4KeRGRFFPIi4ikmEJeRCTFFPIiIimm\nkBcRSbFWQ97M7jCzZWY2q4V1bjSzt8xsppkNiLZEEREpVTF78ncCBzf3ppkdAnzP3bcFzgBuiag2\nERFpp1ZD3t2fBT5tYZWjgHvCdV8EupvZ5tGUJyIi7RHFMflewLs57ffCPhERiZlOvIqIpFinCLbx\nHrB1Trt32LcOM9NEOSIiJXB3K+Xrit2Tt3BpykTgJwBmNgj4zN2XNbchd0/8cvnll8deQzGLxrPj\n1VkNNarO6Jf2aHVP3szuA+qBHma2GLgc6Bzki49z98lmdqiZvQ18CZzSropERCQyrYa8u48oYp1R\n0ZQjIiJR0onXJtTX18ddQqpUy3hWQ53VUCOoziSx9h7vadOHmXklP69aPPUU3HXXuv2jR8POOzf/\ndWbW7uN1IpJ84c96SSdeo7i6RtqosRHWrMm2Z8+GCRNg662hUydYtQrefx9GjGg55EVEWqOQr4Br\nr4UFC7LtP/yh6fWmT4cttoBp02DvvStTm4ikm0K+Ah5+GGbMgO7dg3b37vD55/Db3+avt9FGla9N\nRNJNIV8hQ4bAlClxVyFSfn379mXRokVxl1E16urqWLhwYdm2r5Avgy++yG+vXh1PHSJxWLRokS4I\naAOzks6nFk0hXwabbw5ffZXfd+CB8dQiIh2bQr5MBg+GI47Ituvq4qtFRDouhXyZDBwI55/fvm3c\neCP89a/Z9jHHwNCh7dumiHQsCvkEWm+94FLKV18NFndYtgy22UYhLyJto5BPoN13D26GyvjqK+jS\nJb56RKR48+bNY7vttou7jLUU8iJSEeeeCzNnxl1FvgED4Prro9veuHHjuPXWWznzzDM57bTTottw\nO2iCsnZauBD69ctfCq+sEZEg4JMU8lHXc9ttt7F8+XJefvllPvzwQ8aPHx/dxttBe/Lt1NgI8+bB\noEHQt2/Qt9tusMsusZYlkkgDBkBDQ9xVBKKegHLw4MFrD9NcfPHFvPnmm9F+QIkU8hEZNQpGjoy7\nChGJS+Fx+H79+sVUST4drhERSTGFvIhIO8ybN49BgwZRW1tLTU1N3tK5c2fez71ULgY6XCMiUqLl\ny5dz/vnnc/XVV9O3b18mTJjAjjvuyG677QbAhhtuyGabbRZrjQp5EZESvfDCC4wbN44tt9wSgOnT\np3PuuefStWvXmCvLUsiLiJTo4IMPXvt65cqVLF68OFEBDwr5qrJwITz3XNxViEhTpkyZQv/+/eMu\nYx0K+TaaPBnOPDPbbmys3GffckuwiEjyTJw4ce2x+CRRyLfRypWwaBEceSRsskm2f5ttyveZnTuv\n+1Spk0+GpUvL95kiUjx357HHHuPHP/5x3KWsQyFfoquugp12qsxn1dau+9CRhB32E+nQli5dSteu\nXRk4cGDcpaxDIS8iFTNzZvTTCZRq5sxgmoUo9OrViwULFkSzsYgp5EWkIqIK1KgMGJC8mspBIS8i\nFRHllL5SPIV8lTv++Ozr3r31gyQi+RTyVWrbbeGtt2Du3KC9eHEQ8iIiuTRBWZX629+Cf//5z2DR\ns19FpCkKeRGRFFPIi4ikmEJeRCTFFPIpMmcOmOUvlZxbR0SSR1fXpMSwYbDDDtn2M8/A00/HV4+I\nJENRIW9mQ4HrCfb873D3sQXvfwe4F+gD1ALXuvtd0ZYqLRk2LFgyrrpKIS8iRYS8mdUANwMHAEuB\nGWb2qLvPzVntLGC2ux9pZj2BN83sXnf/tixVi0hi1dXVYWZxl1E16urqyrr9YvbkBwJvufsiADO7\nHzgKyA15B7qFr7sBHyvgRTqmhQsXxl2C5CjmxGsv4N2c9pKwL9fNwA5mthR4DTgnmvLid8UVwbS+\nmWXkyLgrEhEpXlQnXg8GXnX3/c3se8BUM9vZ3VcUrjhmzJi1r+vr66lPyryjzWhsDB4UMnp0fn/P\nnvHUIyLp19DQQENDQyTbKibk3yM4oZrRO+zLdQrwOwB3n29mC4DtgZcKN5Yb8tWithauuSbuKkSk\noyjcAb7iiitK3lYxh2tmAN83szoz6wycCEwsWGcR8CMAM9sc2A54p+SqREQkEq3uybv7ajMbBUwh\newnlHDM7I3jbxwFXAneZ2azwy/7T3T8pW9UiIlKUoo7Ju/sTQL+CvltzXr9PcFxeREQSRNMaiIik\nmEJeRCTFFPIiIimmkBcRSTHNQplyo0cH1/ln/PKXsPXW8dUjIpWlkE+pzp2hWze4886g3dgIq1bB\niBEKeZGORIdrUuqCC2D58uzy0ENxVyQicVDIi4ikmEJeRCTFFPIiIimmkBcRSTGFvIhIiukSygKL\nF8OKnEedfPRRfLWIiLSXQr7AWWfBY4/l9+XeTCQiUk0U8k3o2xfGjs229eB5EalWCvkmbLopnHBC\n3FWIiLSfQr6DOe002GijbPuii+CII+KrR0TKSyHfQfTsCT/6Uba9ahU8+ywsWxZfTSJSfgr5DmLQ\nIJg6NdteskQTlYl0BLpOXkQkxRTyIiIpppAXEUkxhbyISIop5EVEUkwhLyKSYgp5EZEUU8iLiKSY\nQl5EJMUU8iIiKaaQFxFJMYW8iEiKKeRFRFJMId/BXXklDByYXX7/+7grEpEoaarhDmqDDeCQQ/L7\npk6FvfaKpx4RKQ+FfAfVsydMnpzf16NHPLWISPkUFfJmNhS4nuDwzh3uPraJdeqB64D1gI/cfUiE\ndZbNz38OH3+cbb/0Emy1VXz1iIhEqdWQN7Ma4GbgAGApMMPMHnX3uTnrdAf+ABzk7u+ZWc9yFRy1\nxx6DL77IBvsmm8B3vxtvTSIiUSlmT34g8Ja7LwIws/uBo4C5OeuMAB529/cA3P1fURdaTscdB3fc\nEXcVIiLRK+bqml7AuzntJWFfru2ATc3saTObYWYnRVWgiIiULqoTr52A3YD9ga7AC2b2gru/HdH2\nRUSkBMWE/HtAn5x277Av1xLgX+6+ClhlZs8AuwDrhPyYMWPWvq6vr6e+vr5tFYuIpFxDQwMNDQ2R\nbMvcveUVzGqBNwlOvL4PTAeGu/ucnHW2B24ChgLrAy8C/+7ubxRsy1v7vErbems46KDqPCZvZkQ5\nnj16wIgRcNNNkW1SRCIQ/qxbKV/b6p68u682s1HAFLKXUM4xszOCt32cu881syeBWcBqYFxhwIuI\nSOW1uicf6YdpTz5S5diT/+QTqK3N9u27L0T0V6OIlKise/LScZx3Hnz1VbZ9772wenV89YhI+ynk\nZa1LLslvv/ACfPNNPLWISDQ0C6WISIop5EVEUkyHa6RFK1bAa69l2zU18IMfxFePiLSNQl5aNHMm\nDBiQbXfrBsuXx1ePiLSNQl6adcUVMGpUtj1+PPz97/HVIyJtp5CXZu27b377H/9QyItUG514FRFJ\nMYW8iEiKKeRFRFJMIS8ikmIKeRGRFFPIi4ikmEJeRCTFFPIiIimmkBcRSTGFvIhIiinkRURSrEPN\nXfPFFzBsWH7fRx/FU4uISCV0qJBvbIQnn4Q+fWDLLYO+AQNgm23irUtEpFw6VMhnnH8+nH123FWI\niJRfhwx5Kd2XX8KBB+b33XNP9i8jEUkWhbwUrW9fGDQIVq4M2h9+CG+/DatWxVqWiLRAV9dI0c4+\nG557Lrv8+tdxVyQirVHIi4ikmEJeRCTFFPIiIimmkBcRSTGFvIhIiinkRURSTNfJS7tNnQqbb55t\nDxyom6NEkkIhL+12xhn57UcegWOOiacWEcmnkJeSHX44vPJKtv3mmzB8eHz1iMi6FPJSsh49giWj\nRmd4RBJHP5YiIilWVMib2VAzm2tm88zswhbW29PMGs3s2OhKFBGRUrUa8mZWA9wMHAzsCAw3s+2b\nWe+/gSejLlJEREpTzJ78QOAtd1/k7o3A/cBRTaz3C+AvwIcR1iciIu1QTMj3At7NaS8J+9Yys62A\no939T4BFV56IiLRHVFfXXA/kHqtPRNDPng1XXZVtf/11fLWIiMShmJB/D+iT0+4d9uXaA7jfzAzo\nCRxiZo3uPrFwY2PGjFn7ur6+nvr6+jaWXLwPP4Q//xl69YIuXYK+bbeFjTcu20eKiLRbQ0MDDQ0N\nkWzL3L3lFcxqgTeBA4D3genAcHef08z6dwKT3P2RJt7z1j4vSk8/DfvvDw0NMHhwxT62YsyMSo5n\na157DQYM0B2vIlELf9ZLOkLS6p68u682s1HAFIJj+He4+xwzOyN428cVfkkphYiISPSKOibv7k8A\n/Qr6bm1m3VMjqEuq2HHHgeXsc1x2GVx+eXz1iHRkmtZAIvNv/7buw72vvBISdERJpMNRyEtkttwS\nfvOb/L4rr4ynFhEJaO4aEZEUU8iLiKSYQl5EJMV0TF7K7qmn4Ntvs+3+/WHkyPjqEelIFPJSVrW1\n8PzzwQKwejUcfbRCXqRSFPJSVrl78AC77BJPHSIdlY7Ji4ikmEJeRCTFFPIiIimmkBcRSTGFvIhI\niinkRURSTCEvIpJiCnkRkRTTzVBScatWwbJl2XZNDWy2WXz1iKSZQl4q7oknYIstsu2+fWHBgtjK\nEUk1hbxU1CWXwMcfZ9t3352/Vy8i0VLIS0WdcEJ+e9o0hbxIOSnkJXbffguLF+f39e4dHKsXkfZJ\nVch/9hl88km2vXRpfLVI8ZYsgbq6/L7PP4fvfCeeekTSJFUhf+edcN55cVchbXHaaTB4cLb95JPw\n4IPx1SOSNqkK+Yw//Qm6dMm2+/WLrxZp2X77BUvGZ58p5EWilMqQHz4cunePuwoRkfjp1JaISIop\n5EVEUkwhLyKSYgp5EZEUS+WJV6l+N9wAG2yQbQ8fHtwgJSJto5CXRLrssvz2oEEKeZFS6HCNJMov\nfgErVmSXSZPirkikumlPXhJlvfWCJWPDDeOrRSQNtCcvIpJiCnkRkRRTyIuIpFhRIW9mQ81srpnN\nM7MLm3h/hJm9Fi7PmtkPoi9VRETaqtUTr2ZWA9wMHAAsBWaY2aPuPjdntXeA/dz9czMbCtwGDCpH\nwdIxXXwx9OyZbf/sZ3DYYfHVI1Itirm6ZiDwlrsvAjCz+4GjgLUh7+7TctafBvSKskjpuDbaCHbZ\nJXtJ5bffwuzZcMABcVcmUh2KCflewLs57SUEwd+c04HH21OUSMZee8HMmdn2xx/n79GLSMsivU7e\nzIYApwD7NrfOmDFj1r6ur6+nvr6+5M8bPx6efz7bfv31kjclIpIYDQ0NNDQ0RLItc/eWVzAbBIxx\n96Fh+yLA3X1swXo7Aw8DQ919fjPb8tY+ry1OPhkmTIAtt8zvf+ONjvF8UDMjyvGsBpk9+RtvDO6O\nFekIwp91K+Vri9mTnwF838zqgPeBE4HhBQX0IQj4k5oL+HLZemtYuLCSnygiUj1aDXl3X21mo4Ap\nBJdc3uHuc8zsjOBtHwdcCmwK/NHMDGh095aO24u0yxNPwOefZ9sDBsDhh8dXj0hStXq4JtIPK8Ph\nmoaGjrsn35EP1xQ67TS4/fbK1yNSCeU+XCOSGJtuCo2N+X11dfHUIlINFPJSVcygU6d1+0SkaZq7\nRkQkxRTyIiIpppAXEUkxHZOXVHj0UZg1K9vu3x/uvju+ekSSQiEvVW/IkODSyowZM2DNmvjqEUkS\nhbxUvQkT8tuHHw4ffBBPLSJJo2PyIiIppj15SaWVK/OnKIZg6gORjkYhL6k0Zw7sumu23anTunfK\ninQECnlJnUsugdNPz7bvvTe4+kakI1LIS+rsvXd+e8aMeOoQSYKqCvm334ZVq7LtTz+NrxYRkWpQ\nVSE/bNi6J9M0A6EU49tvgweC53rkEfje9+KpR6RSqirkAfbYAy68MNvu0iW+WqQ67LADHH10tr1k\nCbz0EowdC5ttlu0/4YR1fxGIVLuqemjIrrtCnz46iZbRER8aEoWHHoIRI7Jtd1i9Gv7nf/L7RZKi\nPQ8N0c1Q0uEMGxZcTplZ3ngj7opEykchLyKSYgp5kdCpp0K3btnl0kvjrkik/aruxKtI1DbZBM47\nL7/vhhvgm2/iqUckSgp56fA22wyuvTa/749/jKcWkagp5EWaMX8+TJ6cbffsCQMHxlePSCkU8iLN\nePjhYMk4+GB44on46hEphUJepAnPPBNcP59x6qkwbRrss0+2b9NNYdKkytcm0hYKeZEm7LlnfnvI\nEJg7N9ueMwfmzYPFi/PX69ULamvLX59IsRTyIkW46ab89n/8B9x667pzJ33wAWy+eeXqEmlNokM+\nd88J8megFInTSSfln4R9+ulg3vrPPoPOnbP9G20E661X+fpEMhId8v37r9u33XaVr0Ok0D775B+f\n//rrIOS33z5/vSefhIMOqmxtIrkSPUGZGRx7bDDXSEavXvDDH5ahuCqkCcqSY9asYG8+Y8GC4Iaq\njTeG9dfP9l96KZx1VuXrk+rWngnKEr0nD7DTTnDiiXFXIdKynXcOloz58/MPLzY2wvjx8OWXla9N\nOrbE78lfdhlccUUZi6pi2pOvHitXQteu0KNHcOllxhFHrHu3rUihVO/Ji6RBbS0MH57fN2kSPPBA\n8FjLjD591r2SR6Q9FPIiFbD++nDfffl9Rx4J776bvdZ+wYLgypzCSzB/9Suo0XyxUiKFvEhMJk7M\nbx9/fDCNQuEUxwsXQqecn9SLL9azjaV4RYW8mQ0FrieYf/4Odx/bxDo3AocAXwInu/vMwnVaMn8+\nnH56W75CJF0eeCB/KoXrrgueQztpUnB+atmyoP/BB/OfRXvBBXDooZWtVapHqydezawGmAccACwF\nZgAnuvvcnHUOAUa5+2Fmthdwg7sPamJbzZ54nTUr+B+3f//8hyv/9KfBvCGV1NDQQH19fWU/tATV\ncuK1WsYz6XU+9xyceWYD3bvXYxZcvTN9Ouy1F3z3u8FzatesCf7NLJMnB9fzW84pu65dyz/RWtLH\nMqNa6iz3ideBwFvuvij8sPuBo4Dc+1GPAu4BcPcXzay7mW3u7svaWtCVVwbXxsepWv7DV4tqGc+k\n17nPPnDMMQ2MGVMPwNKlQd+//gWffBKc3K2tDY7fZ17vuWfwy2H//YNtzJsXTL0wcmTLn3XeebD7\n7qXXmvSxzKiWOtujmJDvBbyb015CEPwtrfNe2NdsyM+alX9M8oMPiqhERNbaaqvgZG1bXHUV3HVX\n8BdAUzJX+tx3H+yxR/AXQFNLTU3z75nBO+/A889n2wsWBDcy7rhjsP3ly2HFCthvv/zt5W537txg\nYrim3qupCaaQ6N8/+MukVJ9/XvrXVouKn3jt1i34d8WKSn+yiFxySbA0Z/58OOec4NxAc8uaNa33\nffNNEOTu8NJLwfsffACvvRaE85o1wec98kjL9V5/fXTfe3Nuuy34N/fIZzGvS/maYr6+X7/g/OTo\n0a3XXoxijskPAsa4+9CwfRHguSdfzewW4Gl3fyBszwUGFx6uMbPkH0AWEUmgch6TnwF838zqgPeB\nE4GC2zqYCJwFPBD+UvisqePxpRYpIiKlaTXk3X21mY0CppC9hHKOmZ0RvO3j3H2ymR1qZm8TXEJ5\nSnnLFhGRYlR07hoREamsstwsbWZDzWyumc0zswubWedGM3vLzGaa2YBy1NGa1uo0s8Fm9pmZvRIu\nv46hxjvMbJmZzWphnSSMZYt1JmQse5vZU2Y228xeN7Ozm1kv1vEsps6EjOf6Zvaimb0a1nl5M+vF\nPZ6t1pmE8QzrqAk/f2Iz77d9LN090oXgF8fbQB2wHjAT2L5gnUOAv4Wv9wKmRV1HRHUOBiZWuraC\nGvYFBgCzmnk/9rEsss4kjOUWwIDw9UbAmwn9f7OYOmMfz7COLuG/tcA0YGDSxrPIOpMynr8E7m2q\nllLHshx78mtvnnL3RiBz81SuvJungO5mVuknYxZTJ0CsJ4vd/Vng0xZWScJYFlMnxD+WH3g43Ya7\nrwDmENzPkSv28SyyToh5PAHcfWX4cn2Cc3yFx39jH8/ws1urE2IeTzPrDRwK3N7MKiWNZTlCvqmb\npwr/B23u5qlKKqZOgL3DP43+ZmY7VKa0NknCWBYrMWNpZn0J/vJ4seCtRI1nC3VCAsYzPLzwKvAB\nMNXdZxSskojxLKJOiH88rwMuoOlfQFDiWGoC05a9DPRx9wHAzcD/xlxPNUvMWJrZRsBfgHPCPeVE\naqXORIynu69x912B3sBecf/ybk4RdcY6nmZ2GLAs/AvOiPCvinKE/HtAn5x277CvcJ2tW1mn3Fqt\n091XZP7Mc/fHgfXMbFOSJQlj2aqkjKWZdSIIzgnu/mgTqyRiPFurMynjmVPPcuBpYGjBW4kYz4zm\n6kzAeO4DHGlm7wB/BoaY2T0F65Q0luUI+bU3T5lZZ4KbpwrPFE8EfgJr76ht8uapMmu1ztzjXWY2\nkOCS008qW2bw8TT/mz0JY5nRbJ0JGsvxwBvufkMz7ydlPFusMwnjaWY9zax7+HpD4EDyJy6EBIxn\nMXXGPZ7u/it37+Pu2xBk0VPu/pOC1Uoay8jnrvEquXmqmDqB483s50Aj8BXw75Wu08zuA+qBHma2\nGLgc6EyCxrKYOknGWO4DjAReD4/POvArgiusEjOexdRJAsYT2BK424LpyGuAB8LxS9TPejF1kozx\nXEcUY6mboUREUkwnXkVEUkwhLyKSYgp5EZEUU8iLiKSYQl5EJMUU8iIiKaaQl8iZ2Rozuyanfb6Z\nXVbhGu40s2PD17eZ2fbt3F6dmb3eTP9KM3vZzN4ws2lm9tMitjfYzCa1pyaRYlT8Qd7SIXwNHGtm\nvyvlrkEzq3X31VEV4+4/i2pTzfS/7e67w9oJxf5qZrj73SVuTyQy2pOXcvgWGAecV/hGuOf7f+Fs\nf1PD6VUze95/MrMXgLFmdrmZ3WVmz5jZAjM7xszGmtksM5tsZrXh111qwQMhZlnwQPl1mNnTZrab\nmR1hwYMjXrHgYTHzw/d3N7MGM5thZo9nbnEP+2eGd52eVcw37u4Lw+/7nHAbXSx4oMq0cG//iCbq\n29PMng/ff9bMtg37/25mO+es9w8z+4GZ7ZfzfbxsZl2LqU06JoW8lIMDfwBGmlm3gvduAu4MZ/u7\nL2xn9HL3vd19dNjehmCqhKMIHqTwf+6+M7AKOCyzPXffK+zvYsFsfk0X5T7J3Xd1992A14BrLJgI\n7EbgOHffE7gT+G34JeOBs8LZC9viFaBf+PqSsO5BwP7A78P5U3LNAfYN/xq4HPhd2H874a3rZrYd\nsL67vw6MBs4Mv48fEtyGL9IkhbyURTg17t2Ee7Q59iaYZQ9gAsHsexkPFaz7uLuvAV4Hatx9Stj/\nOtA3fH1AuJc8CxgC7NhabWb2n8BKd7+FIIx3AqaGe+yXAFuFE1p1d/fncmotVu4kbQcBF4XbbiCY\nz6dPwfobA38Jj/lfB2Smwf0LcFj4V8spwF1h/3PAdWb2C2CTcIxEmqRj8lJONxDs1d6Z09fScegv\nC9pfQzA7k5k15vSvATqZ2foEfzHs5u5LLXh25wYtFWRmPwKOI9gDhiCQ/+nu+xSs172l7bRiN4K9\n88z2j3P3twq2v0VO8zcEsw4ea2Z1BFPh4u5fmdlU4GhgGLB72D/WzB4j+GvmOTM7yN3ntaNeSTHt\nyUs5GIC7fwo8CJyW897zwPDw9Y+Bf7RlmwU2IPil8bEFD9g4vsUNBAF6MzDM3b8Ju98ENrNg6lbM\nrJOZ7eDunwOfmdn/C9cbWUxt4YnXawgOAQE8CZyd835TD1/uTnZe8MKZBe8ItzU9rAkz28bdZ7v7\n1QRTZrfryiFJN4W8lEPu3vq1QI+cvrOBU8xsJkFwntPE17S2zaAjCL3bgdnA48D0ZtbPvP4psCnw\nv+GJy8c8eL7vMIKTvTOBVwkOKQGcCvzRzF5ppbZtMpdQEjwr+Hp3zzzw4TcED6CYFR6O+a8mvv5q\n4L/N7GUKfibd/RVgOdlDNQDnmtnrYb3fhN+7SJM01bBIgpnZVgSHcrS3LiXRnrxIQpnZScALBA8M\nESmJ9uRFRFJMe/IiIimmkBcRSTGFvIhIiinkRURSTCEvIpJiCnkRkRT7/7xcAG9ocBzCAAAAAElF\nTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "kws = dict(bins=np.arange(0, 4, 0.05), histtype='step', lw=1.5, normed=True)\n", "print('{:>6} {:>8} {:>8} {:>8} {:>8} {:>8}'\n", " .format('n', 'Mean', 'Std', 'Err.RMS', 'Median', '% > 1/λ'))\n", "d_h = np.mean(d1, axis=1) * λ\n", "d_h_err_rms = np.sqrt(np.mean(((d_h - 1)**2)))\n", "plt.hist(d_h, **kws, label = r'$\\hat\\tau$');\n", "print('%6d %8.5f %8.2f%% %8.3f%% %8.3f %8.2f%%' % \n", " (n, d_h.mean(), d_h.std()*100, d_h_err_rms*100, np.median(d_h),\n", " (d_h > 1).sum() * 100 / num_iter))\n", "plt.legend(fontsize=18)\n", "plt.xlabel('Normalized Delays')\n", "plt.axvline(1, color='k');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Since both $\\hat\\tau = T_n / n$ and $\\hat\\lambda_u = (n-1)/T_n$\n", "are unbiased estimators their empirical mean \n", "computed on $N$ (`num_iter`) sets of `n` delays, will converge\n", "to $\\lambda^{-1}$ and $\\lambda$ respectively. In symbols:\n", "\n", "$$\\textrm{E}[\\hat\\lambda_u] = \\lambda$$\n", "\n", "$$\\textrm{E}[\\hat\\tau] = \\lambda^{-1} \\quad\\Rightarrow\\quad \\frac{1}{\\textrm{E}[\\hat\\tau]} = \\lambda$$\n", "\n", "Below we plot the convergence of the expected value,\n", "approximated by an empirical average with number of \"trials\" $N \\to \\infty$.\n", "\n", "Let's start plotting $1/\\textrm{E}[\\hat\\tau]$ for $N \\to \\infty$ first:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZEAAAEACAYAAAByG0uxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmUXGWd//H3J4QdDWFLlLDIIpAcFlFCHB3tI4wJi7Ko\nbCo7okME+c0ogXNmklE8QziKJkZEfiwT+InAMCBwRAgMNAxKMJiFwGQjSEiCCUgWSNhC+vv747ll\nVXfX1tXV3VXdn9c5lbr3uc+997m3K/Wt57nPc68iAjMzs1oM6usCmJlZ83IQMTOzmjmImJlZzRxE\nzMysZg4iZmZWMwcRMzOrWcUgIulGSaslPVsmz1RJSyTNlXRYljZC0qOSnpc0X9LFBfmHSpohaZGk\nhyQNqc/hmJlZb6qmJnIzMLbUQknHAPtGxP7AhcB12aL3gf8TEaOATwIXSTowWzYBeCQiDgAeBS6v\nsfxmZtaHKgaRiHgSWFsmywnALVnep4EhkoZFxKqImJulbwAWALsXrDM9m54OnFhb8c3MrC/V45rI\n7sDygvmV5IMFAJL2Bg4DZmZJu0XEaoCIWAXsVodymJlZL+vxC+uSdgDuAi6JiI0lsvneK2ZmTWhw\nHbaxEtijYH5EloakwaQAcmtE3FuQZ3XW5LVa0nDg1VIbl+QAY2ZWg4hQT++j2pqIslcx9wFnAkga\nA6zLNVUBNwH/GxFTiqxzdjZ9FnAvZUSEXxFMnDixz8vQKC+fC58Ln4vyr95SsSYi6TagBdhZ0svA\nRGArICLi+oh4QNKxkl4ANpIFB0mfAr4KzJc0h9RkdUVEPAhMBu6UdC6wDDil7kdmZmY9rmIQiYgz\nqsgzvkja74EtSuRfAxxdTQHNzKxxecR6E2lpaenrIjQMn4s8n4s8n4vep95sO6uFpGj0MpqZNRpJ\nRANdWDczM+vEQcTMzGrmIGJmZjVzEDEzs5o5iJiZWc0cRMzMrGYOImZmVjMHETMzq5mDiJmZ1cxB\nxMzMauYgYmZmNXMQMTOzmjmImJlZzRxEzMysZg4iZmZWMwcRMzOrmYOImZnVrGIQkXSjpNWSni2T\nZ6qkJZLmSvpYpXUlTZS0QtLs7DWue4dhZmZ9oZqayM3A2FILJR0D7BsR+wMXAr+oct1rIuLw7PVg\ntQU2M7PGUTGIRMSTwNoyWU4AbsnyPg0MkTSsinV7/Nm/ZmbWs+pxTWR3YHnB/MosrZLxWfPXDZKG\nlMv45pvwy192p4hmZtYT+urC+rXAPhFxGLAKuKZc5o99DL75zV4pl5mZdcHgOmxjJbBHwfyILK2k\niHitYPb/AveX3cHKSQBMmgQtLS20tLTUUEwzs/6rtbWV1tbWXt+vIqJyJmlv4P6IOLjIsmOBiyLi\nOEljgJ9GxJhy60oaHhGrsulLgSMi4owS+46ddgrWrIEqimpmZoAkIqLHrz1XrIlIug1oAXaW9DIw\nEdgKiIi4PiIekHSspBeAjcA55daNiJuBqyUdBrQBL5F6dZXk4GFm1piqqon0JUmx447BunUOJmZm\n1eqtmohHrJuZWc2aIoi4BmJm1piaojkLUhkbvKhmZg3DzVlmZtbwHETMzKxmDiJmZlYzBxEzM6uZ\ng4iZmdXMQcTMzGrmIGJmZjVzEDEzs5o5iJiZWc0cRMzMrGYOImZmVjMHETMzq5mDiJmZ1cxBxMzM\nauYgYmZmNXMQMTOzmlUMIpJulLRa0rNl8kyVtETSXEkfq7SupKGSZkhaJOkhSUO6dxhmZtYXqqmJ\n3AyMLbVQ0jHAvhGxP3Ah8Isq1p0APBIRBwCPApdXXWIzM2sYFYNIRDwJrC2T5QTglizv08AQScMq\nrHsCMD2bng6c2IUym5lZg6jHNZHdgeUF8yuztHJ2i4jVABGxCtitDuUwM7NeNrivC5CJ8osnpX8n\nQUtLCy0tLT1dHjOzptLa2kpra2uv71cRFb6/AUl7AfdHxCFFll0HPBYRd2TzC4HP5moaxdaVtABo\niYjVkoZn6x9UYt+RizFVFNXMzABJRIR6ej/VNmcpexVzH3AmgKQxwLpcACmz7n3A2dn0WcC9VZbD\nzMwaSMWaiKTbgBZgZ2A1MBHYCoiIuD7LMw0YB2wEzomI2aXWjYibJe0E3AnsASwDTomIdSX275qI\nmVkX9VZNpKrmrL7kIGJm1nWN1pxlZmbWiYOImZnVzEHEzMxq5iBiZmY1cxAxM7OaNVUQ2bixr0tg\nZmaFmiqIvPpqX5fAzMwKNVUQMTOzxuIgYmZmNWuqIKIeH3tpZmZd4SBiZmY1cxAxM7OaOYiYmVnN\nHETMzKxmDiJmZlYzBxEzM6tZUwURP5TKzKyxOIiYmVnNHETMzKxmFYOIpBslrZb0bJk8UyUtkTRX\n0mEF6eMkLZS0WNJlBekTJa2QNDt7jev+oZiZWW+rpiZyMzC21EJJxwD7RsT+wIXAdVn6IGBatu4o\n4HRJBxasek1EHJ69HqymsK6JmJk1lopBJCKeBNaWyXICcEuW92lgiKRhwGhgSUQsi4hNwO1Z3pwu\n97VyEDEzayz1uCayO7C8YH5FllYqPWd81vx1g6Qh1ezIQcTMrLEM7oFtVlPDuBb4fkSEpCuBa4Dz\nSmefBMBPfgInndRCS0tLtwtpZtaftLa20tra2uv7VVTx817SXsD9EXFIkWXXAY9FxB3Z/ELgs8BH\ngEkRMS5LnwBEREyudtvZ8oBUxqVLYZ99unB0ZmYDlCQioseHaFfbnCVK1zDuA84EkDQGWBcRq4FZ\nwH6S9pK0FXBalhdJwwvWPxl4rppCPFdVLjMz6y0Vm7Mk3Qa0ADtLehmYCGxFqlVcHxEPSDpW0gvA\nRuAc0sLNksYDM0jB6saIWJBt9uqsK3Ab8BKpV1dFt9wCX/xiVw7PzMx6UlXNWX2psDkLfHHdzKwa\njdacZWZm1omDiJmZ1cxBxMzMauYgYmZmNXMQMTOzmjmImJlZzRxEzMysZg4iZmZWMwcRMzOrWVMG\nkXXr+roEZmYGTRpEhg6F9ev7uhRmZtZUQWSHHfLT773Xd+UwM7OkqYJIIfX4bcXMzKySpgoihXfw\ndRAxM+t7TRtEzMys7zVVEDEzs8bSdEHkqafSu5uzzMz6XtMFkdmz07uDiJlZ32uqIBIBg7ISe8Ch\nmVnfqxhEJN0oabWkZ8vkmSppiaS5kg4rSB8naaGkxZIuK0gfKmmGpEWSHpI0pOoCZyW++eZq1zAz\ns55STU3kZmBsqYWSjgH2jYj9gQuB67L0QcC0bN1RwOmSDsxWmwA8EhEHAI8Cl1dT2MLeWT/4QTVr\nmJlZT6oYRCLiSWBtmSwnALdkeZ8GhkgaBowGlkTEsojYBNye5c2tMz2bng6cWE1hP/MZ2H77anKa\nmVlvGFyHbewOLC+YX5GlFUsfnU0Pi4jVABGxStJulXZy/vlwxBHw0Y+m+cH1KLmZmXVLT3wV19Jv\nqsIwQnHDDXDDDfmU9993Dy0zs75WjyCyEtijYH5ElrYVsGeRdIBVkoZFxGpJw4FXy+8i+MY34GMf\ngyOPhMMPhw9+0HfyNTMrRb30K7vaLr6idA3jPuBMAEljgHVZU9UsYD9Je0naCjgty5tb5+xs+izg\n3ooFULqwnjsvb7xRZcnrQIJPf7r39mdm1iwq1kQk3Qa0ADtLehmYSKplRERcHxEPSDpW0gvARuAc\n0sLNksYDM0jB6saIWJBtdjJwp6RzgWXAKZXLkYJIW1uXj7Eu/vznvtmvmVkjqxhEIuKMKvKML5H+\nIHBAkfQ1wNHVFDAnF0QKu/nOnw8HH9yVrXRdLnjst1/P7sfqa8IEmDw5TT/+eOrZV0/vvgtbb13f\nbZo1o6YZsV6sJvL88z2/3332Se9PPNHz+7LqzZ2bb9pcurTz8lwAAfjsZ2HDhuq3/WqFK3QPPADb\nbJP2//rrMHOmH5JmA1dTBRFoXxPpq6Yt61nXXAOLFsFf/5pP+/rX02dAguXLUycLgN/9LtUSjzwy\nLRs+vHivvQ98IP9DpFBum4WvYcPguOPyed59N6W/9VZ6L1y2yy7wyU+mWsm9Fa/sdXbFFfDKK+Xz\nLC/oKD9sGHzve13fj/WcTZvSC9KPiSlTUqefU07J/9AAeO65vitjj4qIhn4BAREXXRQxdWrEzJm5\nRq2IW2+NHpdvROv5fQ10b77Z/nxDxNNPt/+bd+UVEfHaa53TV62KaGurvP6oURG77FJ82dy55ddd\nsKBz2lNPRaxcWTz/5Zd3Ph/vvptf/qUvtc//4ou9+7cZCDZtinj22fLLly1L06+/HnHLLbV9LiHi\n/vvTdtraIkaOjDjnnPofT/p674Xv6N7YSbcKmAWR8eMjpkyJ+MMf8n+IqVO7e5rLe//9zl9KVl+5\nc/ub31T+j3frrdX9B/2P/0j/OXM2b65uvSeeSO9bblk+X86776Yg1dYWsf32tX+h5F677JK2u3hx\ndfkvvjgFky98oXf/Zv1RqR8VW24ZMX9+xJ13Vv93vPvu/PQ773T9czBsWH76ySfbl3Pz5ojf/z7i\n0kvT8jvuKH1MDiK5AmZB5Nvfjvj4xyM+85n8Cf7Xfy33sei+jn/cd95pv3zmzPQlYrWp9J/pyivz\n05demtZ58cWI5cvT9Nix8bcv9d/8JmLt2vL7e/nlzvtYty7iu9+N+O1v2+d99NF8nr/+NT9dSi5Q\nHXBAPu8ZZ0RMmxbx5S9H3HZb533nvPJK6XOw226p5pQLjhERhx1WOv9NN1V//pvJ2rXp+CZMaH+8\nxx2X3t94o/I2nnkmBdzuBPpRo9rPV7Pfjgo/T5BqJd0pU+61fn36YbNkScTxx4eDyN8KmAWRiy/u\nfNK+//2u/wG7Iref3/0uP51LnzatfVp3vf9++1/P/d177xX/j/D44+3z1fMcR6TmirVrIz70oYhf\n/rJ+263W6tXF0//0p/bn4amnym9n48bqvlg2b66+bDfdFDFiRMQee0Q8/3z57ZY6jq56++0IKQXK\njRtTc19ExIYNtX+ZFgbyrrzOOy/irbdSU1XOD34QccghEWPGRPzoR/U55nLa2vKBac2a4uXcZpv0\n96lcwyaiN76je2Mn3SpgFkQuuaTzSfrhD7vz56ost5+XXspPd2ziqscX3Ntvp+3st1/nZaeemmo8\n/UVbW3rlquO5wNnWln79d7RhQzo/A8EnPxlx0knV53/nnYgVK9L0xo2pXb3UF0ol55/f9S/dXHPy\nH/8Yce+96f/GnDnt84wcmd4HDUrvgwfX9gUPESefHLHvvhHPPZe+QOfNizj22Nyv7sqvk05KPxzm\nzMkf91tvpXPX6N57Lx1vR5s3p/Oeq61CCnrpehoRDiL5IPKd73T+UEyZ0tU/RfVmzcrvZ8WK8h/O\n7iq1renT67ePRlDsl5PV3+bNxZvItt02dVAp1PFawCOPpOuPkJoTn3km1Q5eeSV/MfiRR/L5jz66\n9qCQe33ve8XTK9XGimlrS/9v5s2LmDFjYH/GeiuIKO2rcUkKCC69FH7yk/bLZs5MXTt7Zr/56Yjy\nN3vszil8+23Ybrv229q4EZ58EsaNq88+GsVvfwvHH5+fb2mBxx7rs+IMCBH5B7nl3HwznH12mi78\nXL//PmyxRe3bBdhhh9Sldaut0mf7+edh5Mj0vv32aeDnN7/Zed+FvvxlGDMG/vmfqyuLFSeJiOjx\nG2g13TiR3nDyyfnpn/0svbe0dM739tvpfdGi4tuZOhX23BP+8z9L76swgBxwQBoUt8MO7QMIpOPv\njcGVPakwgIADSG8oNjbmnHPyY2Jynnuu+gCS225uHM+6dfn6w5tvpgACsO228IlPpM/4EUekYPKt\nb3Xed0d33eUA0kyapibyT/8EP/5x+2VPPZV+sdR/n/np3On505/Sf4hChTWUYqex2HZK5VmwAA46\nqPPyu+9uH9TKbauR/f73+ZtYLl4MCxfCF77Qt2UaiN59N422L9SMnyerzDWRDor9cumJD//GjcXT\nd9qp/HovvNC9/e68c+e0ZcvgpJM6pzfLzSDfey/93b773XwAaWuD/fd3AOkrW2+d/t/ccAO0tjqA\nWPc1dRC5vKons3fNDjsUT//IR/LTS5d2vhfT/vu3n58zp/383/0dXHttehbKpk2weXO+OWDpUth1\n1/b577svNYUBnHde+2UPPVT5OBpB7gaFP/pRPs0PEmsM552X7ilm1l1N0Zx1zjnBrrvC1Vd3Xr7l\nlqk9N/fY3O7vLz+9aVP7x/AWa7oq1WTVlS/Ltrb27cTr16eHbhV67bX0fuihMGkSfOMb1W+/rxSe\ng112yR+DmfU8N2d1UOpLedMmmD27fvsprHFU8xz3Vas6p916a366mhidO7bDDkvvHQMIpJrKrrvC\nX/4CF17YszefPOQQ+P73u7eN3DH94Q/pHDiAmPVPTRFE1q5tf2vvjjo2HdVq3rz89YZf/rJ4ni9+\nsf38sGFw8cVpOncnzzPPbJ/nqquq2//s2dUHhy22qP8X82OPpZ408+fDxIm1b6fwqZM91QXbzBpD\nUwSRmTPLL//4x+uznyuvzE8Xay6KKH6771zeUrchv+yy0vvcvDk/XanrI8BLL+Wnd9st9fnvriee\ngDVr4HOfy3dbrkUErF4NQ4ak+YceKj6WwMz6j6b4L75mTfnl9Xpg1F13pfeuPk991Kj0fuyx7dP/\n7d/y02+/3bkWA13/kt1rr/bz556bAs9//3fXtvPii6lb8W23pQusxXqHddWgQSmQ5nz+893fppk1\ntqa4sA7lyzhiRPsH99S+r/R+9dWpW2ot6+Ycc0x6Al5Hmzenrq/bbltbGXPefz91Kuio2j9nNRf+\nV65M5Rw6tHSetrbUnPfzn7dPf++94uUzs97RUBfWJY2TtFDSYkmdGmck7SjpbknzJM2UNLJg2SWS\n5mevSwrSJ0paIWl29hrXcbtVH0Qd6lOFT9H7zne6v71iAQTStYzuBhBIF/2L1Ww6Wreu/XzuCYAd\nfeAD6RpLRP66zO67p/ExxUY952yxRecAAg4gZgNFxa9fSYOAacBYYBRwuqQDO2S7ApgTEYcCZwFT\ns3VHAecBnwAOA46XtE/BetdExOHZ68FaD6Irt2sopTCINMsX4L33pi/3adPgwx/Op+dG0l92WapF\nPPRQaqqT4I9/7LydiHQxfJdd0nyxIHPbbalnWDn/8z/p+eSF59LM+rdqfsOPBpZExLKI2ATcDpzQ\nIc9I4FGAiFgE7C1pV+Ag4OmIeDciNgOPA4U38ahLVaseQaQ7vZEg3zMLen8U8EUXpZ5lkAJArmaW\nG1czbhx85Sud1zv33HRblXLuuCO9f+1rKVBNmZJflhvdf8896Zg//enUDbke11fMrDlUE0R2Bwqv\nOKzI0grNIwsOkkYDewIjgOeAv5c0VNJ2wLHAHgXrjZc0V9INkobUeAzdvuVIBNx5Z5rOfWl2VW5M\nyY47dq8stSo2tqSUJ59M15BuvLH4bVUA3nkn1TxOOaV9+ne+k5ZBfnS/b2FiNnBVMZyuKlcBUyTN\nBuYDc4DNEbFQ0mTgYWBDLj1b51rg+xERkq4EriE1fRUxqWC6JXvVz4IF+emOX5pd0Zd9FHJ3Ti2W\n/t576VrIPffAs8/Cpz5VeXtbb53vaTVjRvueVi0t+Y4HRxxRn5qgmXVPa2srra2tvb7fir2zJI0B\nJkXEuGx+AulhJyWH/0n6M3BwRGzokP5DYHlEXNchfS/g/og4pMi2AoKddsp39R0+vPNI8e58gZ96\nar4m0uCd1cqaNy8/6v3Xv0739zr//Ppt/9574cQT26ctW5a/x5eZNY5G6p01C9hP0l6StgJOA+4r\nzCBpiKQts+kLgMdzASS7NoKkPYGTgNuy+YIRBZxMavoqqXCsyNixnZfXehuQZ57JB5BLL61tG43i\n0ENTU9P48ekaSD0DCMAJJ8BPf5qfP+88BxCzga6qcSJZ99sppKBzY0RcJelCUo3k+qy2Mh1oA54H\nzouI9dm6TwA7AZuASyOiNUu/hdRjqw14CbgwIlYX2XencSKbNnXuQfXaa/neRV1xyin5h0Y9/DAc\nfXTXtzHQlHuGipk1ht6qiTTlYMPcU9QKx4esWpXuY9X17bffrplZf9BIzVkNqeNYhmruuNuRg4aZ\nWfc0bRDpqJaHHfneTmZm3dO0zVlpWT6tlmsiHQNPg58KM7OquTmriwpvqW5mZr2j3wSR7j7pb/r0\n+pTDzGwgcRABRo/u/DRCMzOrrF63Pek1ufs1ddTV5qzcY3Cvuab44EUzM6usaWoiBx+c3kvVOLpa\nE8ndjfbEE2HkyPJ5zcysuKYJIl/6Unov7EE1bVq+ZlJtEPna19KjanNBZP78+pXRzGygaZogklMY\nLC66CPbJHnFVTXPWb38Lv/oV/Mu/wHbbpTTfxtzMrHZNHUQKPfxw5XUnZ/cd/vGPU01k9OjaBima\nmVnSNEHkgAPSe+ETBAtddFHlbRQ+3vWCC4o/KtbMzKrXNEFkzJji6V2pSXT3CYhmZtZe0wSRSrck\n+ehHyy/v7mBEMzPrrGmCSCk775zet9++fL6rruqcdsEF9S+PmdlA0vRB5K67qst3zz2d03IByMzM\natP0QWToUBgyBE4+uXy+Z57pnLbttj1TJjOzgaJpgki5ayLr16exH9U46qj89KhR3SuTmdlA1zRB\npJxhw+DUU6vLO2RI/jYnuVHwZmZWm6qCiKRxkhZKWizpsiLLd5R0t6R5kmZKGlmw7BJJ87PXxQXp\nQyXNkLRI0kOShtR6EMceCwceWF3eu++Gr3611j2ZmVmhik82lDQIWAwcBbwCzAJOi4iFBXmuBt6M\niB9IOgD4eUQcLWkU8GvgCOB94EHgwoh4UdJk4PWIuDoLTEMjYkKR/QcES5fCvvumtI5Fzo0VKXco\nuTxTp8K3vgVr1sBuu5U9dDOzptVITzYcDSyJiGURsQm4HTihQ56RwKMAEbEI2FvSrsBBwNMR8W5E\nbAYeB3KXwE8Aco+Cmg6c2K0jIQWGShYvhsGDHUDMzOqhmiCyO7C8YH5FllZoHllwkDQa2BMYATwH\n/H3WdLUdcCywR7bOsIhYDRARq4CyX+u5WsZXvtJ52cSJ6b1Ul91LLslPT5tWbi9mZtYV9Xoo1VXA\nFEmzgfnAHGBzRCzMmq0eBjbk0ktso0xj1KS/ffm/+moL0NJu6a67li/c1Kn56fHjy+c1M2tGra2t\ntLa29vp+q7kmMgaYFBHjsvkJQETE5DLr/Bk4OCI2dEj/IbA8Iq6TtABoiYjVkoYDj0XEQUW2FRBE\nwH77wb//e+fayHXXpescUPy6SOH9te66y72yzKz/661rItXURGYB+0naC/gLcBpwemGGrGfVWxGx\nSdIFwOO5ACJp14h4TdKewElA7laK9wFnA5OBs4B7KxWk1A0UB1dZn7r/fjjuuOrymplZZRW/fiNi\ns6TxwAzSNZQbI2KBpAvT4riedAF9uqQ24HngvIJN/JeknYBNwD9GxBtZ+mTgTknnAsuAU2o9iGpv\nrnj88bXuwczMiqnYnNXXCpuzSvnEJ+BPf0rT5ZqzGvxQzczqpreas/pFECm85uEgYmbWWONEzMzM\niuoXQeS00/q6BGZmA1O/CCIf/GB++vXXYYstYO3afNoRR8DTT/d+uczM+rt+EUTOPDM/vcsuqbfW\nrFn5tFmzYOnS3i+XmVl/1y+CyAEHdE574IH28zfd1DtlMTMbSPpFENlii85pr7/efv7b3+6dspiZ\nDST9oovvG2+kh011FJFegwalJi71eGc3M7PG0Ei3PWlajz8OO+2Uph1AzMzqr180Z229dfH0lhZ4\n8cVeLYqZ2YDSr4MIwIIFvVcOM7OBpl9cE0n5yi9v8MM0M6sr3/bEzMwaXr8JIi+8AGec0delMDMb\nWPpNENl339LPWDczs57Rb66JAKxaBR/6UPFlDX6YZmZ15eeJZLoSRFL+4ukNfphmZnXlC+t1NGNG\nX5fAzKx/6rdB5B//MT/9uc/1XTnMzPqzqoKIpHGSFkpaLOmyIst3lHS3pHmSZkoaWbDsUknPSXpW\n0q8kbZWlT5S0QtLs7DWufocFP/95frrYDRrNzKz7KgYRSYOAacBYYBRwuqQDO2S7ApgTEYcCZwFT\ns3U/DHwbODwiDiHdq6vwOYTXRMTh2evBbh+NmZn1qmpqIqOBJRGxLCI2AbcDJ3TIMxJ4FCAiFgF7\nS9o1W7YFsL2kwcB2wCsF6/m2iGZmTayaILI7sLxgfkWWVmgecDKApNHAnsCIiHgF+DHwMrASWBcR\njxSsN17SXEk3SCpyM3czM2tk9boV/FXAFEmzgfnAHGCzpB1JtZa9gPXAXZLOiIjbgGuB70dESLoS\nuAY4r/jmJzFpUppqaWmhpaWlTsU2M+sfWltbaW1t7fX9VhwnImkMMCkixmXzE4CIiMll1nkROAQY\nB4yNiAuy9K8DR0bE+A759wLuz66bdNxWl8aJXHYZvPUW/OxncOCBsGiRx4iY2cDTSA+lmgXsl33R\n/4V0Yfz0wgxZU9RbEbFJ0gXAExGxQdLLwBhJ2wDvAkdl20PS8IhYlW3iZOC5ehzQ5ILQNmtWeqKh\nmZn1jIpBJCI2SxoPzCBdQ7kxIhZIujAtjuuBg4DpktqA58mapSLij5LuIjVvbcrer882fbWkw4A2\n4CXgwroeGfCBD9R7i2ZmVqjf3fbEzMx82xMzM2sCDiJmZlYzBxEzM6uZg4iZmdWsKYLIuLremtHM\nzOqlKYJIqQdNmZlZ32qKIGJmZo2pKYKIayJmZo2pKYKImZk1JgcRMzOrmYOImZnVrCmCiK+JmJk1\npqYIImZm1piaIoi4JmJm1piaIoiYmVljaoog4pqImVljaoogYmZmjclBxMzMalZVEJE0TtJCSYsl\nXVZk+Y6S7pY0T9JMSSMLll0q6TlJz0r6laStsvShkmZIWiTpIUlD6ndYZmbWGyoGEUmDgGnAWGAU\ncLqkAztkuwKYExGHAmcBU7N1Pwx8Gzg8Ig4BBgOnZetMAB6JiAOAR4HLS5ehK4fUf7W2tvZ1ERqG\nz0Wez0Wez0Xvq6YmMhpYEhHLImITcDtwQoc8I0mBgIhYBOwtadds2RbA9pIGA9sBK7P0E4Dp2fR0\n4MSaj2J9p5RWAAAFGUlEQVSA8H+QPJ+LPJ+LPJ+L3ldNENkdWF4wvyJLKzQPOBlA0mhgT2BERLwC\n/Bh4mRQ81kXEf2fr7BYRqwEiYhWwW6kCuCZiZtaY6nVh/SpgqKTZwEXAHGCzpB1JNY69gA8DO0g6\no8Q2otTGo+QSMzPrS4oK39CSxgCTImJcNj8BiIiYXGadF4FDgHHA2Ii4IEv/OnBkRIyXtABoiYjV\nkoYDj0XEQUW25RBiZlaDiOjxdpzBVeSZBewnaS/gL6QL46cXZsh6Vr0VEZskXQA8EREbJL0MjJG0\nDfAucFS2PYD7gLOByaSL8fcW23lvnAQzM6tNxZoIpC6+wBRS89eNEXGVpAtJNZLrs9rKdKANeB44\nLyLWZ+tOJAWeTaRmrvOzYLMTcCewB7AMOCUi1tX9CM3MrMdUFUTMzMyKadgR65UGODYrSSMkPSrp\neUnzJV2cpZccfCnpcklLJC2Q9PmC9MOzQZyLJf20IH0rSbdn6zwlac/ePcrqSRokabak+7L5AXke\nIDULS/rP7Piel3TkQDwfxQYoD6TzIOlGSaslPVuQ1ivHL+msLP8iSWdWVeCIaLgXKbi9QOrVtSUw\nFziwr8tVp2MbDhyWTe8ALAIOJF0b+l6WfhlwVTY9ktQMOBjYOzsvuRrk08AR2fQDpE4MAN8Crs2m\nTwVu7+vjLnM+LgX+H3BfNj8gz0NWxv8AzsmmBwNDBtr5IPXifBHYKpu/g3TNdMCcB+DTwGHAswVp\nPX78wFBgafa52zE3XbG8fX3CSpzEMcDvCuYnAJf1dbl66Fh/AxwNLASGZWnDgYXFjh34HXBklud/\nC9JPA36RTT9I6gUHabDna319nCWOfQTwMNBCPogMuPOQle+DwNIi6QPqfJCCyLLsC20wqQPOgPv/\nQfoBXRhEevL4X+2YJ5v/BXBqpbI2anNWNQMcm56kvUm/OGaSPiDFBl92PBcrs7TdSeclp/Ac/W2d\niNgMrMs6MjSanwDfpf0YoYF4HgA+AvxV0s1Z8971krZjgJ2P6DxAeX1EPMIAOw9FlBqcXY/jX58d\nf6ltldWoQaTfk7QDcBdwSURsoPNgy3r2eGi4btKSjgNWR8RcypevX5+HAoOBw4GfR8ThwEbSr8yB\n9rnoOEB5e0lfZYCdhyo0zPE3ahBZSbp1Ss4I8vfcanpK9xG7C7g1InLjY1ZLGpYtHw68mqWvJHWD\nzsmdi1Lp7daRtAXwwYhY0wOH0h2fAr6oNDD118DnJN0KrBpg5yFnBbA8Ip7J5v+LFFQG2ufiaODF\niFiT/Uq+B/g7Bt556Kg3jr+m791GDSJ/G+CodOv400hto/3FTaT2yikFabnBl9B+8OV9wGlZj4qP\nAPsBf8yqtOsljZYk4MwO65yVTX+F7OaYjSQiroiIPSNiH9Lf99GI+DpwPwPoPORkTRXLJX00SzqK\nNOZqQH0uSM1YYyRtk5X/KOB/GXjnQbSvIfTG8T8E/INSL8GhwD9kaeX19QWkMheWxpF6Li0BJvR1\neep4XJ8CNpN6nM0BZmfHuhPwSHbMM4AdC9a5nNTrYgHw+YL0jwPzs3M0pSB9a9JAziWk6y179/Vx\nVzgnnyV/YX0gn4dDST+g5gJ3k3rJDLjzAUzMjulZ0iDmLQfSeQBuA14h3eXjZeAcUkeDHj9+UqBa\nAiwGzqymvB5saGZmNWvU5iwzM2sCDiJmZlYzBxEzM6uZg4iZmdXMQcTMzGrmIGJmZjVzEDEzs5o5\niJiZWc3+Pwz3xkQoGKAyAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "rate_t = 1 / ((d1.mean(axis=1)).cumsum() / np.arange(1, num_iter+1)) / λ\n", "\n", "plt.plot(rate_t[100:])\n", "plt.ylim(0.98, 1.02)\n", "plt.axhline(1, color='k');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And now, let's plot the same for $\\textrm{E}[\\hat\\lambda_u]$:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZEAAAEACAYAAAByG0uxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XucHFWd9/HPNwlBIRrut4QEIcgdAnKJ6z4yKkpA9gHx\npQIqyCLioyjiolzcNVEQgd1FQBYQhQgKC66yAosCBhjQfW0wGEiCEggBEggSRAiQcAuZ3/PHqdq+\nTN+m09PTnfm+X6+mq0+dqjpVGepX51SdU4oIzMzMmjFiqAtgZmbdy0HEzMya5iBiZmZNcxAxM7Om\nOYiYmVnTHETMzKxpdYOIpCskLZM0r0aeiyQtlPSApMlZ2nhJd0r6o6T5kr5clH9DSbdLeljSbZLG\ntmZ3zMysnRqpicwADqw2U9JBwHYRsT1wAnBZNutN4KsRsQvwbuCLknbM5p0GzIyIHYA7gdObLL+Z\nmQ2hukEkIn4HvFAjy6HA1Vnee4GxkjaPiGci4oEsfQXwEDCuaJmrsumrgMOaK76ZmQ2lVtwTGQc8\nWfR7KYVgAYCkbYDJwKwsabOIWAYQEc8Am7WgHGZm1maDfmNd0hjg58BJEbGySjaPvWJm1oVGtWAd\nS4Gti36Pz9KQNIoUQH4SETcW5VmWNXktk7QF8Gy1lUtygDEza0JEaLC30WhNRNmnkpuAowEkTQGW\n501VwJXAnyLiwgrLfCabPga4kZoCCCKG92fatGlDXoZO+fhY+Fj4WNT+tEvdmoika4EeYGNJS4Bp\nwGggIuLyiPiVpIMlPQqsJAsOkt4DfBKYL+l+UiQ4IyJuBc4Ffibp74HFwMdbvmdmZjbo6gaRiDiq\ngTwnVkj7b2BklfzPAwc0UkAzM+tc7rHeRXp6eoa6CB3Dx6LAx6LAx6L91M62s2akG+upjB1eVDOz\njiGJ6KAb62ZmZv04iJiZWdMcRMzMrGkOImZm1jQHETMza5qDiJmZNc1BxMzMmuYgYmZmTXMQMTOz\npjmImJlZ0xxEzMysaQ4iZmbWNAcRMzNrmoOImZk1zUHEzMya5iBiZmZNcxAxM7Om1Q0ikq6QtEzS\nvBp5LpK0UNIDkvast6ykaZKekjQn+0xds90wM7Oh0EhNZAZwYLWZkg4CtouI7YETgEsbXPb8iNgr\n+9zaaIHNzKxz1A0iEfE74IUaWQ4Frs7y3guMlbR5A8sO+rt/zcxscLXinsg44Mmi30uztHpOzJq/\nfiRpbAvKYWZmbTZUN9YvAbaNiMnAM8D5Q1QOMzNbA6NasI6lwNZFv8dnaVVFxF+Kfv4QuLn2Jqan\n/06Hnp4eenp6BlxIM7O1WW9vL729vW3friKifiZpG+DmiNitwryDgS9GxIclTQEuiIgptZaVtEVE\nPJNNnwzsExFHVdl2QCpjA0U1MzNAEhEx6Pee69ZEJF0L9AAbS1oCTANGAxERl0fEryQdLOlRYCVw\nbK1lI2IGcJ6kyUAf8ATpqS4zM+syDdVEhpJrImZmA9eumoh7rJuZWdMcRMzMrGkOImZm1jQHETMz\na5qDiJmZNc1BxMzMmuYgYmZmTXMQMTOzpjmImJlZ0xxEzMysaQ4iZmbWNAcRMzNrmoOImZk1zUHE\nzMya5iBiZmZNcxAxM7OmOYiYmVnTHETMzKxpDiJmZta0ukFE0hWSlkmaVyPPRZIWSnpA0p71lpW0\noaTbJT0s6TZJY9dsN8zMbCg0UhOZARxYbaakg4DtImJ74ATg0gaWPQ2YGRE7AHcCpzdcYjMz6xh1\ng0hE/A54oUaWQ4Grs7z3AmMlbV5n2UOBq7Lpq4DDBlBmMzPrEK24JzIOeLLo99IsrZbNImIZQEQ8\nA2zWgnKYmVmbjRrqAmSi9uzpjBgB06dDT08PPT097SiTmVnX6O3tpbe3t+3bVUSd8zcgaSJwc0Ts\nXmHeZcBdEXF99nsBsH9e06i0rKSHgJ6IWCZpi2z5napsOyBYZx14440m9tDMbBiSRERosLfTaHOW\nsk8lNwFHA0iaAizPA0iNZW8CPpNNHwPcWK8ADcQ6MzNrs7o1EUnXAj3AxsAyYBowGoiIuDzLczEw\nFVgJHBsRc6otGxEzJG0E/AzYGlgMfDwillfZfkAwciS8+eaa7ayZ2XDRrppIQ81ZQykPIiNGwOrV\nQ10aM7Pu0GnNWWZmZv10TRDp8AqTmdmw5CBiZmZN65ogYmZmnaergohrI2ZmncVBxMzMmtY1QUSD\n/qCamZkNVNcEEXBNxMys03RNEJEcRMzMOo2DiJmZNc1BxMzMmuYgYmZmTeuqIGJmZp2la4IIuCZi\nZtZpuiaIuDnLzKzzOIiYmVnTHETMzKxpXRVEzMyss3RNEHnlFddEzMw6Td0gIukKScskzauR5yJJ\nCyU9IGlyUfpUSQskPSLp1KL0aZKekjQn+0xtpLD3399ILjMza5dGaiIzgAOrzZR0ELBdRGwPnABc\nlqWPAC7Olt0FOFLSjkWLnh8Re2WfWxsp7CuvNJLLzMzapW4QiYjfAS/UyHIocHWW915grKTNgX2B\nhRGxOCJWAddleXMDvsvh5iwzs87Sinsi44Ani34/laVVS8+dmDV//UjS2EY21Ne3pkU1M7NWGjUI\n62ykhnEJ8O2ICElnAecDx1XPPh2Aa66B9dbroaenZ40LaWa2Nunt7aW3t7ft21U00EYkaSJwc0Ts\nXmHeZcBdEXF99nsBsD/wDmB6REzN0k8DIiLObXTd2fyAVMabb4ZDDhnA3pmZDVOSiIhB7xzRaHOW\nqF7DuAk4GkDSFGB5RCwDZgOTJE2UNBo4IsuLpC2Klj8ceLCRQvieiJlZZ6nbnCXpWqAH2FjSEmAa\nMJpUq7g8In4l6WBJjwIrgWNJM1dLOhG4nRSsroiIh7LVnpc9CtwHPEF6qquu226Dv/u7geyemZkN\npoaas4ZScXMWuDZiZtaITmvOGlJnnDHUJTAzs0q6IoiYmVlnchAxM7OmOYiYmVnTHETMzKxpDiJm\nZtY0BxEzM2uag4iZmTXNQcTMzJrWFUHEvdTNzDpTVwQRMzPrTA4iZmbWtK4IIhr0IcTMzKwZXRFE\nfE/EzKwzdUUQMTOzzuQgYmZmTXMQMTOzpnVFEPE9ETOzztQVQeS559L3fvu1f9vf/356OuzGG9u/\nbTOzTlc3iEi6QtIySfNq5LlI0kJJD0iaXJQ+VdICSY9IOrUofUNJt0t6WNJtksbWKsPcuY3uTut9\n+cvp+4c/HLoymJl1qkZqIjOAA6vNlHQQsF1EbA+cAFyWpY8ALs6W3QU4UtKO2WKnATMjYgfgTuD0\nWgXohH4it9wy1CUwM+s8dYNIRPwOeKFGlkOBq7O89wJjJW0O7AssjIjFEbEKuC7Lmy9zVTZ9FXBY\n7TLUK6WZmQ2FVtwTGQc8WfT7qSytWjrA5hGxDCAingE2q7UBBxEzs840ahDW2UzjU80wcd99aZX3\n3ju0TVud0KxmZtZJWhFElgJbF/0en6WNBiZUSAd4RtLmEbFM0hbAs7U28K53BWPGwOGHF250t8vk\nyYUb+7Nnwz77FOa5hmRmnUptuupttDlLVK9h3AQcDSBpCrA8a6qaDUySNFHSaOCILG++zGey6WOA\nmg/QRsCuu7a/JrBqVemTYeuuWzr/vvvaWx4zs05TtyYi6VqgB9hY0hJgGqmWERFxeUT8StLBkh4F\nVgLHkmaulnQicDspWF0REQ9lqz0X+JmkvwcWAx+vX472Xvn/5jfwoQ+VpvX1lf5evrx95TEz60R1\ng0hEHNVAnhOrpN8K7FAh/XnggEYKmGt3ECkOINtsA088AStWlOb54AfhhRdggw3aVy4zs07SFT3W\nYWhvaj/8MEyYAJ/9bP95G27Y/vKYmXWKrggieQ1kqG5kjx4NS5bAggVDs30z6x6vvAK//W268JXg\nmGPSuUuCO+/sn//Pf4aXXmp/OVula4JIO5qz3ngDXnxxcLdhZmuPhQsLwSL/rL8+vPe9hTxXXw0j\nsjPtBz6Q8nzjG+n3Y4/BVlvB2LGNn9/237+wrT/8obX704yuCCLQniCy7rrp/sbddxfStt66NM/f\n/m0qx+uvD25ZzKwzvfIKHHhgOie9853V891/f/+HcXJnn52W3267QtqIEf0DUh50pLTdiRPhnnsK\ny+y9d/VttEtXBpFvfQtuuql2/mbWn+vpKUwffXRpvt/9Ln2PHg2LFrW2DGY2dGbPhvnz4a9/TSfn\nb34TXnstNTctWpQuLvOaxu23p2W22CLlj0ife+4pTE+eXDhvLV8OzzyTpl9+uXS7t94Kxx9fuUzF\nzV/rr5+a1QHefLOwnpEj2//gUTFFh/eYkxSTJwfvfz9suSWccko6YPvsA7//fSu3Uzn9kUdg++1L\n5+eHrK8PRo2CxYtTU1jxVUX5ejv8MJsNaxGFJqdGLVgAO/R79rQ1ZcmDwumnpxHEt98+jdhx7bVw\nyCHwtrelvGedBf/0T9XWJCJi0B9J6sqaCLTnpLxkSfrHg3QVUW7EiFSOCRNg0qTSeflVi3WPvFng\nH/+x/9Xi4sWF6SeeSP+2Tz6JtVlfX+1Ovo8/njoJ5yJSM9B55xWahZ5/Ps3Lb2bPmtU/gNx1F/zi\nF/3Xf++9cOqpqRyDEUCgcN6Q4JxzUk1n1qy0L0ceWQggkP5WIypfwLbLYIyd1XKVbqy3I4i85S2F\n6fwf7ogjGlu2/CmMvr6BX+lY+zz6aOGCAeA730n/ZhGpuQDgjjtSG3VuwoSUxxcLg+OVV9LFW7UT\n5GGHwUEHwec+l07ojzzS2Ho33rj6vDwAjcrOjOVPhkqw776NbaedHn00fT/7LLz6avrbbNf5pmtO\na0MRRDbZpDD91rem7wPKukhWuxoZM6b0d3lHResMixalt1YWB5DciBGFAAKlAaQ4z9ixsGzZ4JWx\nWH41nYtI7fbdKiL9P7V6dWl6fu+h1hX2L38JJ5yQ8tYLIJ/7XGpy/uUvK8//9rdTWUaNKgSQ8vJ0\nw8XCZpulm+/tLGtXBJFKNZFWPpFQLSAV/0NI8KlP9R8Kpfgkc8YZ6XvLLVM1s1j5DXprj5NOqn3B\nMWlSuqIF2HNPuP761NZdfhGQ2yx7aUHx399LL6UbrPnfaXlT2JpasSKd/Mr/HqUUxN76VrjyytZu\ns9gdd8AnPwmXXJI61776au38edmKj1FEKuecOemm8IoVhUdf77gDPvax9Kjr176WaoHlzj67cMM6\nInUA3n//0jx5zTH/rFyZmr4i4Ac/gHXWgUMPLc3z2mspgFW/r2B1RURHf4DYbbeIU0+N+O53IyLS\nP/8ee0TLvPFG4c9qt90K040o/ZOMePHF/mn55/XXW1dmq+2WW0qPfV9fxMyZpXnK/32K3XBDIX3V\nqohPfKKwntx551X+d27l32a1v6V6n+XL13zbCxbU3kZE2s4BB0TsvXf/+WPHRvzzPze/D8uWrfk+\nDGfp9D745+iuqInkoqw565ZbKjcxDNTKlen71Vfhqqv6b2sgxtZ4W3zebmmD78MfLv09YkRqNpk9\nO90Yf+ON0vmHlb1b8yMfKZzORo2C664r1DRyX/ta+psprxXPnQtf/SqMH1+oQed/Y41aubL/1fFH\nP1oo0403ppoxVB4IdIMN4N3vHtg2y/so7Lhj/zz77Vea/0tfgpkzK9/sfvHFdIxqmTGjcvnnzy/U\n+qzDtSNSrcmHrCZy2mkR3/lOHmEjtt46Yscd+19BNmPcuMJ6+voi/vu/G1+23tXUL34Rccophd+H\nHFJ9PdOmrfGuDCsvvlg5/ZVXqv97TJ9e+vvoo1tTlhdeSOv7/vfrX73nKqX39fVf5vnn629/1ar0\nvXp1/e1Wcs89lZfbd9/+eV9+uX++vMaxenUh31lnpbS//CXi0ksjlixJ6Y89Vr88tuZoU01kyINE\n3QJmQeT009MfZWSl3myzxv8HqWdN1lMviERE3Hpr7f+h//zn1u1LNyk+cc2enU7+228fsWhRabNR\nROEkDREvvRQxa1bhpAUR//ZvEVdeGfHaa6XHsq8v4thjq//7vPlm6/dr8eLq29t888rpV1yRlt1n\nn9L0l19urgwvvVS6nm99q3K+8qD17/8eceONET/9aWrmrSVf5q67miujDS4HkaIgsuuuEWecEXHm\nmfnBidh009adeL/+9TUPIosW9T8xzJmT8jz3XO0g0ttbfd7ChRF77dVc2TrRo49GPPts5SvuSp/H\nHov40Y+qz691wi7fbl5zHTMmbb88UA2GvCxPPVW9jPn0xz5WmG5F2X7844jPfrb68SmvsT311MC3\n0Yp7LzY42hVEuuaeyGA+nbXjjvCZzzS37De+kYY32HbbQtoLL6Sy7rln+r3xxqkNO1f+uG/xMCuz\nZ5fu6w03pCdaVq1KT7V0s5deSk9DbbYZnHlmY8tsu23lIfhzEydWTv/tb0t/b7cd/OlPaXrFivY9\nstnXlz7jxvXvnPjcc+n75z9P3//xH+l71qzWlO2YY1Jv5yeeqDx/vfXS91ZbpeF8xo0b+DZq3QO0\n4aErgkhE4cSaD8deHET+8pc1W//q1c13zDnzzHSSL1bpJVWbb16YrtVnJO/gOH9++j711PQ9enR6\nPLGb5SdNgGnT0ndEaR+BGTNK32OfO+CA1HO3uOf4l79cmD7xRLjmmjT9ve+lgTLLSemE2s5+FcXB\navz40rpA3untox8t5D/ooNKb160wcWLaXl8fXHpp6bxrroGlS+E972ntNm346IogUvyUS35VVRxE\nHntszda/fHnzV/nFJ4nXXqv+DH3x+v/619Lli+VP3ORXzcV+9atUy+lWH/tY6e+8o2Y+fExEqhHm\nV+bF/8a//jVstFHqifvqq+npqryPzgUXwPe/D0cdldbxla9UL8PEiWm05k6zfHnqa9LqgUWLSfD5\nz5cGsqPqvrfUrLauGIDxC18INtkk/U/w7nfD1Knw9renq6df/zpdkV544ZpsI30P5qHo60sdEL/7\n3fT7K19JV8z5tr/+9TS+TyM6/J+sqvKA2ch+1Pq36etLtbq3v33Ny2a2tpE6aABGSVMlLZD0iKRT\nK8zfQNINkuZKmiVp56J5J0man31OKkqfJukpSXOyz9Rq2y9uzspPJn19aVgEgN12a3R3+6s0sOJg\nGDEi9brNXXBB6Uk1Dy6VLFxY+vsPf0jNEOec09oytlpfH/T2pum8X0Zxz+NG1Mo7YoQDiNlQqxtE\nJI0ALgYOBHYBjpRU3g3pDOD+iNgDOAa4KFt2F+A4YG9gMnCIpKJb0JwfEXtln1trlyO13R50UPod\nUWjO2HXXentR3eOPN79sq0yZkk6I1cYJmjQp1bhye++dmr1OP7012x87No1OmjfFPfZY5ea0gRo5\nEt73vjRERd6EVHwfw8y6XyM1kX2BhRGxOCJWAdcB5bd4dwbuBIiIh4FtJG0K7ATcGxGvR8Rq4G7g\n8KLlBlTVKu7Z2tdXCCI//vFA1lKqfIyrofA//5O+b7gBbr4Zjj22MC+/Cp86dXCasSZNSk9NTZmS\nntbJ37a2yy5rtt5vfaswXdzzOa89mtnaoZEgMg4ofjjxqSyt2Fyy4CBpX2ACMB54EPg/kjaUtB5w\nMFD8wtkTJT0g6UeSaj4sWN6evnp1ekELpMdim1U+ZPtgiyh9T0HxyXb33dMLZy69NA15cWvNulmy\nJo+CPvfc4Lyd8bnnYPr0/uk//Wnrt2VmQ6tV7xM5B7hQ0hxgPnA/sDoiFkg6F/gNsCJPz5a5BPh2\nRISks4DzSU1f/cyePZ0xY/J3C/cAPSVX5eXDSHe6w4vqYt/8Zv/5666bhqau5B3v6N8Et3hx9f4S\ntVxySe35K1emmsOzzw5sHKM77ihMP/007LFHegz74x8feBnNrDG9vb305jch26ju01mSpgDTI2Jq\n9vs0Uk/Ic2ss8ziwW0SsKEv/DvBkRFxWlj4RuDkidq+wrvj854Nx46oP17zbbjBvXs3dqCq/kn/5\n5erDfw+Gvr70qO+mmw582aefhgcfhAMPTL8nTep/872RdeSdy5YsSX0znnmmtGZz332w9dapj0sj\nL19644100/9v/ib9Xrmy0KHNzNqrk57Omg1MkjRR0mjgCKDkaXZJYyWtk00fD9ydB5Ds3giSJgAf\nAa7Nfm9RtIrDSU1fVb34YvV5jfTxWLkyNXs9WLSV888vTLczgEC6kd5MAIHUw/hDHyrUYuqNDvzm\nm6VPoUWU9k7eeuvKT6ntvXehk+Tdd6fHaV95pfp2fvazQgABBxCz4aBuEMluiJ8I3A78EbguIh6S\ndIKkz2XZdgIelPQQ6Smuk4pW8QtJDwI3Al+IiOzNxpwnaZ6kB4D9gZNrleNf/qX6vEaCyJgx6bWW\nxY8D/8M/1F+ukxXfT5k/v/pQMJddll6U9fjjqbNZce/8j3ykNG9fX+X7JKeckl4RvP761YPWH/5Q\nmM5HFjCztVtD90Syx293KEv7QdH0rPL5RfPeWyW9Ze/622uv+nkOOQT+678qz7voolaVpP3yd7fv\nnjUEVmqdzIdlKR7fK3fDDaW/pZSv/H3ixQHiq19NT8S9/nqqFQH853+mvi+5aq8NNrO1S1cMe1LP\n9deXDiVSSa3mqnyokW7UyNNZM2ZUTr/tturLvP/9hY5+5YPs3XxzGvcpDyBQWqN5b8XLBjNbG60V\nQQTqvzmu+Aq9/Gp91arWl2eoPPlk/zfUVVP+vvhqli1LPeQbGTl50aJ0/8TMhoeuCCIDGWOpmuKh\nwWfOLJ3X7a/hLH4F6YQJlfNcfHG6J3LbbfCudw1s/euumwbqk0oDyZVXpu+Xsrtcr75aucnMzNZe\nXRFEasl7rdcLIk8/XZh+/vnC9Dvf2foytdt55xXeDV/N+PGwzTap9nHffc33fs8DyapVqWd9RLrh\nHlEYxt7Mho+uDyJ5R8OB9NwuHq6922shuaPLHlP43vdS5768R/7ee7duWxKMalU3VTPral0RRBq5\nah7IS6U22qgwveGGAy9Pp1q0KL3gaNNN00uarr8+DYBY3i/EzKxVuuJ9Ip/7XLDTTnByjZ4kS5eW\nPi3Ufz2lv/Ph5efNW7Oh5M3MOlEn9VgfclJhsMVqBhoL82awThgK3sysW3VFEIkovTFeyUAHYcyH\n75g0qbkymZlZlwQRqH/Po5E+DMXyN+LtvHPtfGZmVl1XBJGINQ8iO+3UuvKYmVnSFUEE6geRes1Z\na9NTWGZmnaJrgkjeqbCaejWRvr70xkAzM2udrgkia1oTWb0a9tyzdeUxM7MuCiIrVtSe30hNpF5t\nxszMBqZrgki9R3wbCSIjRpTm+/Of17xcZmbDWVcEkUY6EtZrzsqDiFRo1ho9es3LZmY2nHVFEGlE\nvXeCzJ1beI1uPnjgOusMbpnMzNZ2a00QKX6VazWvv56+8+DhIGJmtmYaCiKSpkpaIOkRSadWmL+B\npBskzZU0S9LORfNOkjQ/+3y5KH1DSbdLeljSbZLGlq8310hzVq0b7/nyeU3kHe9I337/hZnZmqkb\nRCSNAC4GDgR2AY6UtGNZtjOA+yNiD+AY4KJs2V2A44C9gcnA30nK3313GjAzInYA7gROr1WOL32p\ndjkPPrj6vPxmel7z2H772usyM7PGNFIT2RdYGBGLI2IVcB1waFmenUmBgIh4GNhG0qbATsC9EfF6\nRKwG7gYOz5Y5FMjfx3cVcFitQnziE7ULedxx1eflQeTd707fZ5yR3kVuZmZrppEgMg4oPuU+laUV\nm0sWHCTtC0wAxgMPAv8na7paDzgY2DpbZvOIWAYQEc8ANd8xuCadDfv60pNY+TtF1lknvS7WzMzW\nTKtecnoOcKGkOcB84H5gdUQskHQu8BtgRZ5eZR1V73zMmTO96FdP9ilVL4gM5M2HZmbdpre3l97e\n3rZvt+6bDSVNAaZHxNTs92lARMS5NZZ5HNgtIlaUpX8HeDIiLpP0ENATEcskbQHcFRH9xtqVFMcd\nFxx/PEyZUr2c11wDRx1Ved7KlTBmzMBfXGVm1q066c2Gs4FJkiZKGg0cAdxUnEHSWEnrZNPHA3fn\nASS7N4KkCcBHgPwdhTcBn8mmjwFurFaA/FW2/1voEYV7JB/+cPquVROZNavuPpqZWRPqNmdFxGpJ\nJwK3k4LOFRHxkKQT0uy4nHQD/SpJfcAfSU9k5X4haSNgFfCFiHgpSz8X+JmkvwcWAx+vVY7i5qi+\nPvjhD+Hcc2GjjeCgg2oHkZUr6+2lmZk1o6F7IhFxK7BDWdoPiqZnlc8vmvfeKunPAwc0tv3SmgjA\n296WPgDvfGftILLllo1sxczMBqprbjcXB5HyMa9GjqwdROr1MTEzs+Z0RRCRSoPIQw+Vzh85svYo\nvvfeOzjlMjMb7roiiJS/Y33bbUvnv/EGvPpqe8tkZmat6ycy6MrviRSbMSN9n3xye8piZmZJV9RE\noHYQMTOzodEVNZHip7N6evrPf+tbCyPzmplZ+3RFEIHCPZGzz+4/75RTSt+f/tOfppF699uvkObH\nfM3MWq9rgkheE6k0dEn5I76f/nT/vJ/61OCVzcxsuFor7omMGlV44VSx/JW5kyfDkUcOTrnMzIaz\nrgkitUbhfe45mDevf/pZZ6Vvj+JrZjY4uubUWvwukHIXXAC33NI//fzz07eDiJnZ4OiKU2sEvPBC\nmt577/7z8zcW/ulPpen5e9cffLD2sChmZtacrggiAJtskr4r3Rv54hfT9667Vn9niIc+MTNrva4J\nIttuC88+W3nezJnpOwIefbRyHr9T3cys9bomiABsumnl9OJRfddbr3Keww9vfXnMzIa7rgoi1bzl\nLYXpURV6vowcCbvt1r7ymJkNF10RROq9G33ddesv77G3zMxaryuCSD1PP12YrvQUloOImdngWCuC\nSHFNpDyIXHRR//eRmJlZazR0apU0VdICSY9IOrXC/A0k3SBprqRZknYumneypAclzZN0jaTRWfo0\nSU9JmpN9plbbfr3mrD32KEy/733pOx+o8aST8nI0sqdmZjYQdYOIpBHAxcCBwC7AkZJ2LMt2BnB/\nROwBHANclC27FfAlYK+I2J004OMRRcudHxF7ZZ9bm92J4iCyaBFssw2cdlqzazMzs0Y1UhPZF1gY\nEYsjYhUkSoRNAAAIbElEQVRwHXBoWZ6dgTsBIuJhYBtJ+QO5I4H1JY0C1gOK7mDQkvpB+RNZI0e6\n5mFm1g6NBJFxQHFXvaeytGJzgcMBJO0LTADGR8TTwL8CS4ClwPKImFm03ImSHpD0I0ljm9wHxpYt\n6fsfZmbt0ar3iZwDXChpDjAfuB9YLWkDUq1lIvAi8HNJR0XEtcAlwLcjIiSdBZwPHFdp5XPnTmf6\n9DTd09NDT9nrDXfaqUV7YWbWpXp7e+nt7W37dhV17lpLmgJMj4ip2e/TgIiIc2ss8xiwOzAVODAi\njs/SPw3sFxEnluWfCNyc3TcpX1d8+tPB1VdXL2Olp6/KH+utd3PezGxtIomIGPSG/UYafmYDkyRN\nzJ6sOgK4qTiDpLGS1smmjwfuiYgVpGasKZLeIknAB4CHsnxbFK3icODBagWoFwDq3f/YsfwxADMz\na4m6zVkRsVrSicDtpKBzRUQ8JOmENDsuB3YCrpLUB/yRrFkqIn4v6eek5q1V2ffl2arPkzQZ6AOe\nAE5o6Z4VKR8i3szMWqNuc9ZQkxSf+lTwk5/Uy1f6OwK23BKeecZNWWY2/HRSc9aQazYIjBzZ2nKY\nmVmprggijdhll/5pW23V/nKYmQ0nXRFEGuk4mL9PvdjZZ8OZZ7a+PGZmlrSqn8igaqQ5a511+qcd\ncED6mJnZ4OiKmkgj3OHQzKz9uiKINFIT2WyzwS+HmZmV6oog0giPl2Vm1n4+9ZqZWdMcRMzMrGkO\nImZm1rS1MohssslQl8DMbHjoiiAy0GFPVq8enHKYmVmprggiA+UgYmbWHmtVEFm6NH339Q1tOczM\nhou1ZtgTSAMuPvhg5SFQzMys9boiiAxEpdF8zcxscHRFc1Yjo/iamVn7dUUQ8ZsJzcw6U0NBRNJU\nSQskPSLp1ArzN5B0g6S5kmZJ2rlo3smSHpQ0T9I1kkZn6RtKul3Sw5JukzS2dbtlZmbtUDeISBoB\nXAwcCOwCHClpx7JsZwD3R8QewDHARdmyWwFfAvaKiN1J92COyJY5DZgZETsAdwKnr/nurN16e3uH\nuggdw8eiwMeiwMei/RqpiewLLIyIxRGxCrgOOLQsz86kQEBEPAxsI2nTbN5IYH1Jo4D1gOxBXA4F\nrsqmrwIOq1YAN2cl/h+kwMeiwMeiwMei/RoJIuOAJ4t+P5WlFZsLHA4gaV9gAjA+Ip4G/hVYQgoe\nyyPijmyZzSJiGUBEPAP4jSBmZl2mVTfWzwE2lDQH+CJwP7Ba0gakGsdEYCtgjKSjqqyjan3DNREz\ns86kqHOGljQFmB4RU7PfpwEREefWWOYxYHdgKnBgRByfpX8a2C8iTpT0ENATEcskbQHcFRH9XnIr\nySHEzKwJETHoHSQa6Ww4G5gkaSLwZ9KN8SOLM2RPVr0SEaskHQ/cExErJC0Bpkh6C/A68IFsfQA3\nAZ8BziXdjL+x0sbbcRDMzKw5dWsikB7xBS4kNX9dERHnSDqBVCO5PKutXAX0AX8EjouIF7Nlp5EC\nzypSM9dns2CzEfAzYGtgMfDxiFje8j00M7NB01AQMTMzq6Rje6zX6+DYrSSNl3SnpD9Kmi/py1l6\n1c6Xkk6XtFDSQ5I+VJS+V9aJ8xFJFxSlj5Z0XbbM/0ia0N69bJykEZLmSLop+z0sjwOkZmFJ/5Ht\n3x8l7Tccj0elDsrD6ThIukLSMknzitLasv+SjsnyPyzp6IYKHBEd9yEFt0dJT3WtAzwA7DjU5WrR\nvm0BTM6mxwAPAzuS7g19PUs/FTgnm96Z1Aw4CtgmOy55DfJeYJ9s+lekhxgA/h9wSTb9CeC6od7v\nGsfjZOCnwE3Z72F5HLIy/hg4NpseBYwdbseD9BTnY8Do7Pf1pHumw+Y4AH8LTAbmFaUN+v4DGwKL\nsr+7DfLpuuUd6gNW5SBOAX5d9Ps04NShLtcg7esvgQOABcDmWdoWwIJK+w78Gtgvy/OnovQjgEuz\n6VtJT8FB6uz5l6Hezyr7Ph74DdBDIYgMu+OQle/twKIK6cPqeJCCyOLshDaK9ADOsPv/g3QBXRxE\nBnP/ny3Pk/2+FPhEvbJ2anNWIx0cu56kbUhXHLNIfyCVOl+WH4ulWdo40nHJFR+j/10mIlYDy7MH\nGTrN94CvUdpHaDgeB4B3AM9JmpE1710uaT2G2fGI/h2UX4yImQyz41BBtc7Zrdj/F7P9r7aumjo1\niKz1JI0Bfg6cFBEr6N/ZspVPPHTcY9KSPgwsi4gHqF2+tfo4FBkF7AX8W0TsBawkXWUOt7+L8g7K\n60v6JMPsODSgY/a/U4PIUtLQKbnxFMbc6npK44j9HPhJROT9Y5ZJ2jybvwXwbJa+lPQYdC4/FtXS\nS5aRNBJ4e0Q8Pwi7sibeA/xfpY6p/w68X9JPgGeG2XHIPQU8GRH3Zb9/QQoqw+3v4gDgsYh4PrtK\n/k/gbxh+x6FcO/a/qfNupwaR/+3gqDR0/BGkttG1xZWk9soLi9LyzpdQ2vnyJuCI7ImKdwCTgN9n\nVdoXJe0rScDRZcsck01/jGxwzE4SEWdExISI2Jb073tnRHwauJlhdBxyWVPFk5LemSV9gNTnalj9\nXZCasaZIektW/g8Af2L4HQdRWkNox/7fBnxQ6SnBDYEPZmm1DfUNpBo3lqaSnlxaCJw21OVp4X69\nB1hNeuLsfmBOtq8bATOzfb4d2KBomdNJT108BHyoKP1dwPzsGF1YlL4uqSPnQtL9lm2Ger/rHJP9\nKdxYH87HYQ/SBdQDwA2kp2SG3fEApmX7NI/UiXmd4XQcgGuBp0mjfCwBjiU9aDDo+08KVAuBR4Cj\nGymvOxuamVnTOrU5y8zMuoCDiJmZNc1BxMzMmuYgYmZmTXMQMTOzpjmImJlZ0xxEzMysaQ4iZmbW\ntP8PMKB1gBVIVs4AAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "rate_t2 = ((n - 1) / d1.sum(axis=1)).cumsum() / np.arange(1, num_iter+1) / λ\n", "\n", "plt.plot(rate_t2[100:])\n", "plt.ylim(0.98, 1.02)\n", "plt.axhline(1, color='k');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In both cases there is convergence to $\\lambda$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Appendix 1. Computation of integral \"A\"\n", "\n", "Here we want to compute the integral:\n", "\n", "$$ \\int_0^{\\infty} x^{n-2} e^{-\\lambda x}\\, dx $$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The $\\Gamma$ function is defined as:\n", "\n", "$$\\Gamma(s) = \\int_0^{\\infty} t^{s-1}\\,e^{-t}\\,{\\rm d}t$$\n", "\n", "Therefore, with the variable substitution $t = \\lambda x$:\n", "\n", "$$ {\\rm d}t = \\lambda {\\rm d} x$$\n", "\n", "$$\\frac{1}{\\lambda^{n-2}} \\int_0^\\infty t^{n-2} e^{-t}\\,\\frac{{\\rm d}t}{\\lambda}\n", "= \\frac{1}{\\lambda^{n-1}} \\int_0^\\infty t^{n-2} e^{-t}\\,{\\rm d}t\n", "= \\frac{1}{\\lambda^{n-1}} \\Gamma(n-1) = \\frac{(n-2)!}{\\lambda^{n-1}}$$\n", "\n", "So we obtain:\n", "\n", "$$ \\int_0^{\\infty} x^{n-2} e^{-\\lambda x}\\, dx = \\frac{(n-2)!}{\\lambda^{n-1}}\\qquad\\qquad \\textrm{(A)}$$\n" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "Finally, note that $n$ is the number of delays, while the number of times/events is $m = n+1$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Appendix 2. Variance of the rate estimators\n", "\n", "We want to compute the variance for the family of rate estimators $\\Lambda_{n,c}$ defined as\n", "\n", "$$\\Lambda_{n,c} = \\frac{n-c}{T_n}$$\n", "\n", "where $T_n = \\sum d_i$ is the sum of the $n$ delays, and $c$ is a real parameter with $c < n$.\n", "\n", "We start by writing:\n", "\n", "$$\\mathrm{Var}\\{\\Lambda_{n,c}\\} = \\mathrm{E}\\left\\{(\\Lambda_{n,c} - \\mathrm{E}\\{\\Lambda_{n,c}\\})^2 \\right\\} \n", "= \\int \\left( \\Lambda_{n,c} - \\overline{\\Lambda}_{n,c} \\right)^2 f(x \\:|\\: n, \\lambda)\\, dx\n", "= \\int \\left(\\frac{n - c}{x} - (n - c)\\frac{\\lambda}{n-1} \\right)^2 f(x \\:|\\: n, \\lambda)\\, dx =$$\n", "\n", "$$= (n-c)^2 \\int \\left(\\frac{1}{x} - \\frac{\\lambda}{n-1} \\right)^2 f(x \\:|\\: n, \\lambda)\\, dx$$\n", "\n", "For brevity, let's define:\n", "\n", "$$\\left(\\frac{1}{x} - \\frac{\\lambda}{n-1} \\right)^2 = \\left(\\frac{1}{x^2} + \\left(\\frac{\\lambda}{n-1}\\right)^2 \n", "-2 \\frac{\\lambda}{x(n-1)})\\right)\n", "= \\left(\\frac{1}{x^2} + a + \\frac{b}{x} \\right)$$\n", "\n", "From **term 2**:\n", "\n", "$$(n-c)^2 \\int a \\, f(x \\:|\\: n, \\lambda)\\, dx = a (n-c)^2 = (n-c)^2\\left( \\frac{\\lambda}{n-1} \\right)^2$$\n", "\n", "From **term 1**:\n", "\n", "$$(n-c)^2 \\int \\frac{1}{x^2} \\, f(x \\:|\\: n, \\lambda)\\, dx\n", "= (n - c)^2 \\int \\frac{\\lambda^n x^{n-3} e^{-\\lambda x}}{(n-1)!\\,}\\, dx \n", "= \\frac{(n - c)^2\\lambda^n}{(n-1)!} \\int x^{n-3} e^{-\\lambda x}\\, dx \n", "= \\frac{(n - c)^2\\lambda^n}{(n-1)!} \\frac{(n-3)!}{\\lambda^{n-2}} =$$\n", "\n", "$$= \\frac{(n - c)^2\\lambda^2}{(n-1)(n-2)}$$\n", "\n", "And from **term 3**:\n", "\n", "$$(n-c)^2 \\int \\frac{b}{x} \\, f(x \\:|\\: n, \\lambda)\\, dx\n", "= b(n - c)^2 \\int \\frac{\\lambda^n x^{n-2} e^{-\\lambda x}}{(n-1)!\\,}\\, dx \n", "= \\frac{b(n - c)^2\\lambda^n}{(n-1)!} \\int x^{n-2} e^{-\\lambda x}\\, dx \n", "= \\frac{b(n - c)^2\\lambda^n}{(n-1)!} \\frac{(n-2)!}{\\lambda^{n-1}} =$$\n", "\n", "$$= \\frac{b(n - c)^2\\lambda}{(n-1)} = -\\frac{2\\lambda}{n-1} \\frac{(n - c)^2\\lambda}{(n-1)}\n", "= -2 \\frac{(n - c)^2\\lambda^2}{(n-1)^2}$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Summing these 3 terms we obtain the variance:\n", "\n", "$$\\mathrm{Var}\\{\\Lambda_{n,c}\\} = \\lambda^2(n-c)^2\\left[ \\frac{1}{n-1} + \\frac{1}{n-2} - 2 \\frac{1}{n-1}\\right]\n", "= \\lambda^2(n-c)^2\\left[\\frac{1}{n-2} - \\frac{1}{n-1}\\right]$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# References\n", "\n", "- [Random: Chapter 6. Point Estimation](http://www.math.uah.edu/stat/point/index.html)\n", "- http://www.amstat.org/publications/jse/v9n1/elfessi.html\n", "- Any statistics book explaining RV parameter's estimators, MLE, MMSE and MVUE." ] } ], "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.5.1" } }, "nbformat": 4, "nbformat_minor": 0 }