Skip to content

Instantly share code, notes, and snippets.

@phobson
Last active September 17, 2020 16:35
Show Gist options
  • Star 2 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save phobson/82f223f4aae24787bc18 to your computer and use it in GitHub Desktop.
Save phobson/82f223f4aae24787bc18 to your computer and use it in GitHub Desktop.
Lognormal python ref
Display the source blob
Display the rendered blob
Raw
{
"metadata": {
"name": "",
"signature": "sha256:a80a44b573fab4a827df2ed5c88c14a378339d1662d63ed1cd9a17526ade1180"
},
"nbformat": 3,
"nbformat_minor": 0,
"worksheets": [
{
"cells": [
{
"cell_type": "heading",
"level": 1,
"metadata": {},
"source": [
"Sorting out the ambiguities with lognormal distributions in Python"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Notes:\n",
"\n",
" + $\\log(x)$ denotes the natural logarith of $x$\n",
" + This notebook refers to the `scipy.stats` modules as just `stats`\n",
" + See here for more info: http://en.wikipedia.org/wiki/Log-normal_distribution"
]
},
{
"cell_type": "heading",
"level": 1,
"metadata": {},
"source": [
"TL; DR"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Using scipy, you can easily use maximum likelihood estimation (MLE)\n",
"do determine the parameters describing the distribution of a sample\n",
"of data. However, you can have to get the terminology straight\n",
"for each type of distribution. This notebook aims to clear that\n",
"up for lognormal distributions.\n",
"\n",
"If you want to estimate the location $(\\mu)$ and scale $(\\sigma)$\n",
"(as defined by Wikipedia) that describe a lognormally distributed\n",
"sample $(X)$ with `stats.lognorm.fit`, do the following:\n",
"\n",
" + fit the distribution, making sure to feed a guess that scipy's `loc` is 0\n",
" + set $\\sigma$ to scipy's shape parameter\n",
" + set $\\mu$ to the natural log of scipy's scale parameter\n",
" + check to confirm that scipy's loc parameter is close to 0\n",
" \n",
"If you already know $\\mu$ and $\\sigma$ of your distribution, \n",
"then you can create that distribution in the manner below:\n",
"\n",
"`myDist = stats.lognorm(sigma, loc=0, scale=numpy.exp(mu))`"
]
},
{
"cell_type": "heading",
"level": 2,
"metadata": {},
"source": [
"More details"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"import numpy\n",
"from matplotlib import pyplot\n",
"import scipy\n",
"from scipy import stats\n",
"import seaborn\n",
"\n",
"seaborn.set()\n",
"%matplotlib inline\n",
"\n",
"print('numpy: {}'.format(numpy.version.full_version))\n",
"print('scipy: {}'.format(scipy.version.full_version))"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"numpy: 1.9.1\n",
"scipy: 0.14.0\n"
]
}
],
"prompt_number": 1
},
{
"cell_type": "heading",
"level": 2,
"metadata": {},
"source": [
"Basic Concepts"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"A (perhaps over-)simplification of the lognormal distirbution would state that\n",
"a sample $(X)$ is lognormally distributed if $\\log(X)$ is normally\n",
"distributed. This implies that the value of a sample from a populate\n",
"that is *truly* distributed lognormally can never be zero. You can\n",
"get really dang close, but not actually there.\n",
"\n",
"So let's define the parameters that describe the distirbution of \n",
"$\\log(X)$ as $\\mu$ (kind of the mean) and $\\sigma$ (kind of the\n",
"standard deviation). In other words:\n",
"\n",
"\\begin{equation*} X = e^{(\\mu + \\sigma \\: Z)} \\end{equation*}\n",
"\n",
"Where $Z$ is a standard normal variable.\n",
"\n",
"Wikipedia calls these the location and scale.\n",
"For our purposes, let's change that to `wikiLoc` and `wikiScale`.\n",
"That's a pain in the ass, but it gets worse due to things out of\n",
"my control, so please bare with me."
]
},
{
"cell_type": "heading",
"level": 2,
"metadata": {},
"source": [
"Sample statistics"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In addition to $\\mu$ and $\\sigma$, we can also talk about the sample mean\n",
"and the sample standard deviation. Let's call them `wikiM` and `wikiS`. \n",
"\n",
"In python you can compute those very easily with e.g., `numpy.mean(X)`"
]
},
{
"cell_type": "heading",
"level": 2,
"metadata": {},
"source": [
"Generate a random sample from `numpy`"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You can generate a random sample from a lognormal distribution\n",
"using numpy and `wikiLoc` and `wikiScale`. The wording of the \n",
"arguments (`mean` and `sigma`) is unfortunate, but I don't \n",
"expect that to change anytime soon."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"numpy.random.seed(0)\n",
"N = 37000\n",
"wikiLoc = 2.25\n",
"wikiScale = 0.875\n",
"\n",
"# _n denotes numpy\n",
"X_n = numpy.random.lognormal(mean=wikiLoc, sigma=wikiScale, size=N)"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 2
},
{
"cell_type": "heading",
"level": 3,
"metadata": {},
"source": [
"Plot the histogram"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"bins = numpy.logspace(-1, 3, num=250)\n",
"\n",
"fig, ax = pyplot.subplots(figsize=(8, 5))\n",
"_ = ax.hist(X_n, bins=bins, cumulative=False, histtype='step', \n",
" normed=True, lw=1.25, label='Hist of random numpy sample')\n",
"\n",
"ax.set_xscale('log')\n",
"ax.set_xlabel(r'$X$')\n",
"ax.set_ylabel('PDF')\n",
"ax.legend()"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 3,
"text": [
"<matplotlib.legend.Legend at 0x91d6e80>"
]
},
{
"metadata": {},
"output_type": "display_data",
"png": "iVBORw0KGgoAAAANSUhEUgAAAf0AAAFPCAYAAABHzhf2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8XFX9//HXTJbJMkmatJN0b7qeLtAWaKW0UECpiFIs\nIiquoAUpIPhFEVBBvn7FBQREsCpQ4OvCV38oIFYBEVkjtuwtS04XaEOXtGnT7Mkks/z+mGmYliZN\n2ty5M5n38/FQcufcO/POnE4+c89djicajSIiIiKDn9ftACIiIpIcKvoiIiIZQkVfREQkQ6joi4iI\nZAgVfRERkQyhoi8iIpIhsp16YmOMF1gOzASCwFJr7caE9sXANUAIuNtae5cxJhe4C5gEdAGXWmtf\ncyqjiIhIJnFyT38JkGutnQ9cBdy0t8EYkwPcDCwCTgQuMMaUA+cDbfFtzgfudjCfiIhIRnGy6C8A\nHgWw1q4C5iS0TQM2WGsbrbVdwHPAQmB6wjbrgFHGmGIHM4qIiGQMJ4t+MdCUsByOD/nvbWtMaGsG\nSoBXgdMBjDHzgABQ6GBGERGRjOHYMX1iBb8oYdlrrY3Ef27cr60Y2AP8BZhmjHkWqALWAfW9vUg0\nGo16PJ4BCy0iIpLiDrnoOVn0q4DFwP3xvfY1CW3VwGRjTCnQCpwA3AB8APiXtfZyY8wc4APW2mBv\nL+LxeKira3bkFxDnBQJF6r80pb5Lb+q/9BUIFB18pR44WfQfBBYZY6riy+cZY84B/NbaO40xlwOP\nETvEsMJau90YEwT+aIz5NtBB7GQ+ERERGQCeQTDLXlTfVtOX9jbSl/ouvan/0lcgUHTIw/u6OY+I\niEiGUNEXERHJECr6IiIiGUJFX0REJEM4efa+iEjaCEci7Gnq9QrhQ1Ja7CPLq/0rSQ0q+iIiwJ6m\nIN/61fMD/rw3XHgcw4bkD/jzihwKFX0RkQRXfe5oyop8h/089c1Bfvz7l3td5+WXX+Qvf3mA//7v\nH3Y/9stf3kZl5XgmTZpCVdUznHvu0gNu+9prr+D3FzFx4qSDZtm6dQtXXHEZRxwxk29/+3v9+0X6\n4Itf/DS/+c0fB/x5B6sVK37N0KHDWLLkrKS/toq+iEiCsiJf0vbMD3QL8b2PTZ48hcmTp/S47cqV\nf+GUU07tU9Ffs+ZV5s8/gUsu+fqhh5UB4+at41X0RURccqCbo+19LHEU4Ic//G+2bt1CMBjk7LM/\nQ2XlBFavfp7169dRWTmeiorh3dvfdtstrF37GgCLFn2EhQtP4re/vYdgMMjo0aNZsuST3eueddbp\njBs3nvHjx/Oxj32c22+/hXA4QmNjA9/85lUcccRMPvOZM5k5czY1NZspLS3j+utvoLOzk+9//xoa\nGxsYNWo0kUhsWpV166r52c9+itfrJTfXx5VXfodIJMK1115NRcVwamu386EPfZh33tnIunWW445b\nwFe/evE+v/+BXu/RR/9GTc1mLrzwEoLBIJ///Nncf//DXHLJBUyebHj77Y0UFOQzc+ZRrF79PC0t\nzdx88y949tmnWLXq3zQ0NNLY2MCXv3wB48dP4Pvfv4Y77/xfAK699mrOOefzTJs2A4BgMMi1115F\na2srwWAHF1xwEXPnzuPPf/4jzzzzFO3t7QwZMoQf/vCn/OMfj1BV9QydnZ3s3r2Ls88+h2effZq3\n397IJZdcxvHHn8iXvnQOY8aMYceOWiZNmsK3vvWdfX7fX/3qdtaseZVIJMKnP/1ZTj75lMP9Z9Ur\nFX0RERe9/PKLfO1rX+1e3rZtK0uXXti9N9jW1sZrr73CHXfcC8Dq1f/BmKkce+x8Tjnl1H0KflXV\ns9TWbuOOO+4lFApx0UVLOeaYOXzhC+dRU7N5n4IPUFe3k3vuuY/i4mKeeOJxLrnk60yYMInHH3+U\nv/3trxxxxEy2b9/Gbbf9mkCgnGXLvsJbb73J2rWvMn78BM4/fxk1NZu44orYCMJPfnI9V199LZMm\nTea5557mtttu4ZJLvs727du49dbldHR0cPbZZ/DQQ4/i8/n45CcXv6/oH+j1etoz9ng8TJ8+g8su\n+wbf+Mal5Ofnccstv+D666/j1VdfwuPxEIlEufXW5ezevYuvfvU8/vjHh/D5fGza9A5lZWVs376t\nu+BD7FBIU1MjN910G3v27KGmZjPRaJSmpiZ+9rPleDweLr/8a7z11ht4PB7a29u5+ebbeeKJf/DH\nP97HHXfcy8svv8j99/+B448/kdrabdx8820MHTqM7373Sp555qnu13r++Sq2b9/G8uV3EQwGufDC\n85g7dx5+v7///5D6SEVfRMRFRx89Z59j+r/61e37tBcUFHDppd/gJz+5ntbWVk499bQen2vz5k3M\nmnUUANnZ2cyYcSTvvPMOcOBRhZKSIRQXFwMwbFiAe+9dgc/no62tlcJCf/c6gUA5AOXlFXR2Bqmp\n2cxxxy0AYOzYSoYMKQVg9+5dTJo0GYCZM4/q/l1GjhxFQUEhWVnZlJUNpagoNmHMgWr5gV5vX/v+\nHsZMBcDv91NZOQGAoqIiOjs7ATjmmLkADB06DL+/iMbGBs4440z+/ve/UlExnI985KP7PN+ECRM5\n44xPcN113yEUCvHJT34Gj8dDdnY21133bfLzC6ir20EoFAJg8mQDQGGhn8rK8e97/crKCQwdOiz+\nnszi3Xc3d7/WO+9sxNrq7i994XCY2trt3e+hE1T0RUQS1DcPzGV7A/U8u3fvwtq3+OEPbyQYDHLW\nWadz6qkfxePxEA6H91m3snI8f//7w3zqU58lFArx+uuv8dGPns6GDesP+Nxe73tV99Zbf8r3vvcD\nxo2rZMWKX1Nbux04cGGurJzA2rWvccIJJ7F16xYaGxuA2BeHjRs3MHHiJF599WXGjBkXf46+H8M+\n0Kq5ubns3r0LAGur99+i1+errn4TOIv6+t10dLRTWlrGSSd9iPvu+y0lJUP4wQ9+ss/6b7+9gba2\nNm644Wfs2rWLZcu+QkXFcJ599mnuuONeOjo6WLr0C91fog72u23ZUkNLSwt+v5+1a9dw2mmnxzPF\nvjAdffQxfOtbsS8Yv/3tPYwcOarX5ztcKvoiIgkOdsb9QPJ4PL0WDY/Hw9Chw6iv382yZV/G683i\ns5/9AllZWUyffgS/+tXtjBo1irFjKwGYP/94XnnlJS688Mt0dXXxoQ8tYvJkw4YN63t4nfceO/XU\n07jmmispL69g6tTp3UX2QJmWLDmLH/3o+yxb9hVGjBhJUVFstODKK7/DLbfcQDQaJTs7m6uuuoZo\nNLrfa/f084Ef83g8HHvsfB588E9cdNFSjJnWPQrRF1u2vMtll11EW1sL3/zm1Xg8HnJzc5k9+2ga\nGxu6Rx32Gj16LHfffSdPPvlPIpEI559/IaNHjyY/P5+LLz6fkpIhTJkylV27dnXnS/zve7lj/83N\n9fGDH1xLfX09M2fOZv7846mujh2yOP74hbzyyktcfPH5tLe3sXDhyRQUFPT5dzsUmmVPXKWZvtLX\nYOu7TLs5z2DrvwN55JGVNDQ0cM45n39f2y233MCJJ36Qo4+e42gGJy5nPJxZ9rSnLyICZHm9uonO\nIHSgAY7LL7+EIUNKHS/4sdd37/K8A9GevrgqE/Y2Biv1XXpT/6Wvw9nTT70xJxEREXGEir6IiEiG\nUNEXERHJECr6IiIiGUJFX0REJEOo6IuIiGQIFX0REZEM4djNeYwxXmA5MBMIAkuttRsT2hcD1wAh\n4G5r7V3xbe4CpgAR4HxrrXUqo4iISCZxck9/CZBrrZ0PXAXctLfBGJMD3AwsAk4ELjDGlAMfBgqt\ntccD3weudzCfiOta2rt4bHUNj62u4elXt7odR0QGOSdvw7sAeBTAWrvKGJN4v8NpwAZrbSOAMeY5\nYCFQB5QYYzxACdDpYD4R1zW1dvLHf21giD8Xr9fDibOdnWFLRDKbk3v6xUBTwnI4Pny/t60xoa2Z\nWJF/DsgDqoFfA7c5mE8kZSyaO8btCCKSAZzc028CEucs9FprI/GfG/drKwIagCuBKmvtd4wxo4F/\nGWOOsNb2uscfCBT11iwpLpP7ryP+ifAX+vB6vWn3XqRbXtmX+i/zOFn0q4DFwP3GmHnAmoS2amCy\nMaYUaCU2tP9T4GjeGx3YA+QAWQd7IU0akb4yfdKP+vpWAFpag0QikbR6LzK979Kd+i99Hc6XNSeL\n/oPAImNMVXz5PGPMOYDfWnunMeZy4DFihxhWWGu3GWNuBO4xxjxLrOBfba1tdzCjiIhIxnCs6Ftr\no8Cy/R5el9C+Eli53zYNwJlOZRIREclkujmPiIhIhlDRFxERyRAq+iIiIhlCRV9ERCRDqOiLiIhk\nCBV9ERGRDKGiLyIikiFU9EVERDKEir6IiEiGUNEXERHJECr6IiIiGUJFX0REJEOo6IuIiGQIFX0R\nEZEMoaIvIiKSIVT0RUREMkS22wFEZOAEO8N0dIYAyMnOoiBPH3EReY/+IoikiW27Wlm/pQGA8tIC\npo0rfd86T7+2jT88sR6AhbNGcO5p05KaUURSm4q+SJpY924Dv398HTnZXmZPHnbAog8wvKyA8tL8\nJKcTkXSgoi+SRoYNyWfK6BK6wpEe18nJ9lLg00dbRN5PJ/KJiIhkCO0OiKSgPzyxntr6NgA+sXAC\nYyuKXE4kIoOB9vRFUtCGrY20B0Os2biblvYut+OIyCDh2J6+McYLLAdmAkFgqbV2Y0L7YuAaIATc\nba29yxhzLvCl+Cr5wCygwlrb5FROkVQ1c+JQNmxtdDuGiAwiTg7vLwFyrbXzjTHHAjfFH8MYkwPc\nDMwB2oAqY8zD1tp7gXvj69wO3KWCLyIiMjCcLPoLgEcBrLWrjDFzEtqmARustY0AxpjngIXAn+LL\nc4AZ1tpLHMwnkhZ21LfxcNUmGlqCeDwet+OISBpz8ph+MZC4lx6OD/nvbUsct2wGShKWvw1c52A2\nkbTR0RVm3bsNzDHlnHLMaLfjiEgac3JPvwlIPOXYa63de3Fx435tRcAeAGPMEGCKtfbpvr5QIKAz\nm9NZJvdfR/wT4S/04fV6eb56J6ter2X77jaOn+3DE28D+NLiGRTk5fDzP75CVihywPfN7/eRne3F\n58vBl5vl+HubyX03GKj/Mo+TRb8KWAzcb4yZB6xJaKsGJhtjSoFWYkP7N8bbFgJP9OeF6uqaDz+t\nuCIQKMro/quvbwWgpTVIJBLh7ZoGdta3ccoxoxlZmk803gawa1cL+b5sOjq66ApHDvi+tbQECYUi\nBINdRCNhR9/bTO+7dKf+S1+H82XNyaL/ILDIGFMVXz7PGHMO4LfW3mmMuRx4jNghhhXW2u3x9aYA\nG9//dCKDW0cwzLt1LYwaVsiZCyf0uq6taeD2B9Zy/MwRzJ40LEkJRSTdOVb0rbVRYNl+D69LaF8J\nrDzAdj91KpNIqirwZVNa5KO1vYuxFf5e151WWUpuThar3tyBGTMkSQlFZDDQHflEUsDsyQFOnD2q\nT+vOmz6cedOHY2saHE4lIoONir6IS6o376Guod3tGCKSQVT0RVxy65/W0NkVRpfei0iyqOiLuOhr\nn5x50BPxdjV0JCmNiAx2KvoiSfaHJ9bz7s4WOkPhPq3/5Ctbe2zb0xxkc20zXaEI9/1zHY2tnfjz\ncwYqqogMMir6IklWs6OZcCTKx44bR0Vpfq/r/vqbJ3X/nOV9/3GAR1fX8OjqGs4+eSKbaps568QJ\nlBT6eHNT/UDHFpFBQEVfxGHtwRANLbEb7PhysgCYNq6UJSf0fi0+QHZWz3fKvvrzR9MVjvD1nz/X\n/dipHxhLdpZXRV9EDkhFX8Rhr79Tzy8feh2AIycMHbDnzfdlkxOOHHxFEZE4JyfcEZG4wrxsTj6q\nb9fhi4g4RUVfJAk8Hg++3Cy3Y4hIhtPwvsggsG1Xa49twc4w9zzyFgBej4cLzpiRrFgikmK0py+S\n5ny5Wbxo63ocSQhHIqx+aye7Gjt4oXpnktOJSCrRnr5IGsvO8vLLy0/s07rHmACbazWVqkgmU9EX\nGUB3PPxG9+VyXzptKkdNDricSETkPRreFxlAbcEQ40cUEwxFCIWjrmbZUd/Ov1/fTm39e5P6RKJR\n1r69m41bG11MJiJuUdEX6afXNuziheqdvFC9k/Zg6H3tI4YWkpPl5bHVNfziwbU0tXYmPWNhXg47\nG9q5+2/V2Jo9AHjwkOX1cuv9a/jDv9YnPZOIuE9FX6SffvcPy4q/vckvH3qd+ubgAdc5YdYIKkoL\neMnW9fke+wPpcx+ewk0XLyAwJK/7semVpdxxxUksOWF80vOISGpQ0Rc5BB8/vvfCefZJkzhTxVVE\nUoyKvoiISIZQ0RdxWO3utu6ft9S1sL2+rZe1RUSco0v2RByUneXh+Tdq8efnMKbcT31TBwCjAn6X\nk4lIJlLRFzmIprbO7jPwC/Ny+rzdsCH53HHFyfs8dtyM4QOaTUSkP1T0RQ7iuTXb+dNTGwH44NGa\nKU9E0pdjRd8Y4wWWAzOBILDUWrsxoX0xcA0QAu621t4Vf/xqYDGQA9xurf1fpzKK9NXYCj+lfl+P\n7c+u2UbV2lq21rUwcmhhEpMdur89v4m3tzUBsGjOGKaOK3U3kIg4zskT+ZYAudba+cBVwE17G4wx\nOcDNwCLgROACY0y5MeYk4Lj4NicBExzMJ9JnWV4v2dk9f1x2N3awu7GDU+aMYcb4siQmO3Rvb2ti\nV2MHb27aw54e7jcgIoOLk8P7C4BHAay1q4wxcxLapgEbrLWNAMaY54CFwNHAWmPMQ0AxcIWD+UQG\n1LCSvINev59qzJghtLR3uR1DRJLEyaJfDDQlLIeNMV5rbSTelnjz72agBBgGjAM+Rmwv/2FgqoMZ\nRQ7LDfe9TH5uNvNmVLgdpUcPPPO22xFEJEU4WfSbgKKE5b0FH2IFP7GtCGgAdgPV1toQsM4Y02GM\nGWat3dXbCwUCRb01S4pL9f4rLPSRk+PF58smPz8Xb5aXwFA/3/js0Wzc2sgjz2+KrZOblXK/y2Wf\nOZqOztj8AFPHD6UwPyeWNTsr9vsU5OL1eigqzjuk7Kn2+0r/qP8yj5NFv4rYCXn3G2PmAWsS2qqB\nycaYUqCV2ND+jUAHcBlwszFmJFBI7ItAr+rqNEd4ugoEilK+/1pbg3R1RQgGQ7RndxIJRwi2dzJn\n0ghCnSGi0Whsnc5wyv0uw0t8QOwExLaWDtpaOmJZQ+HY79PWSSQSpbmpo9/Z06HvpGfqv/R1OF/W\nnCz6DwKLjDFV8eXzjDHnAH5r7Z3GmMuBx4idTLjCWrsd+JsxZqExZnX88Yuste7OTyoiIjJIOFb0\n48V62X4Pr0toXwmsPMB2VzqVScQJ4XCUNzfvIcvjcTuKiEivdHMekcNQVJDLtPj17aPL0/fWuuu3\nNBAlyuxJAQry9GdBZLDSp1ukH17bsJumtvcucZswspjLPz3bxUSHr6Qwl1c27OKpV7dx/fnFKvoi\ng5g+3SJ9NHvSMIaV5AEwclh63HWvJ6FQlK5Q7GKaa8+dS0NLkMtvryLYFaY9GCLfpz8NIoORPtki\nPWjrCLG9vpXd8ZnxFhw5wuVEA2fzjthZ28PLCvZ5/Pv3vkhebhbLLz/RjVgi4jAVfZEebK5t4sY/\nvArA+BHFLqcZOCcdNYq5U8sByI8P5fvzc/jeuXOxNXt46Ll33IwnIg5S0Rc5iF9940QG04n5/vwc\n/Pn7ThGcneVl3PAi6ps7XEolIsmgoi/SCw+Qm5PldgwRkQHh5Cx7IiIikkJU9EVERDKEir6IiEiG\n0DF9EXmfLTtb2La7FYhd1je2QrOxiQwGKvoiso+OzjDX3r0agCyvhw9/YIyKvsggoaIvAnSFIkSi\nsQkds7wesrMy88jXhBHFXHrWTAD8BTk8rGv2RQYVFX0R4N5H3uL5N3YA8IVTDScfNcrlRO4o8fuY\nPdnndgwRcUhm7s6IHMDsScMYWqyCJyKDl4q+SJy/IIecbN2IZ38tbV1s29VKezDkdhQROUwq+iLS\nq2fXbOe7d63ijXfq3Y4iIodJx/RFpEdf/fgMQuEo37njP25HEZEBoKIvsr9olHAkQiTqdhD3FebF\nJuYZTBMOiWQyFX2R/XR0hTn/hqeA2IQ7IiKDhYq+SA/OPW0qw8sK3I6RMt6q2UN7Z4i5U8vJy9Wf\nDpF0pBP5RHowtsLPlDFD3I6REoYW5/Hq+l3c8/dqWtt1Fr9IulLRF9nPo6tq3I6Qcq778ge48nNH\nux1DRA6TY2N0xhgvsByYCQSBpdbajQnti4FrgBBwt7X2rvjjLwON8dXettZ+xamMIvs7Y0ElbfHr\n0cuK8lxOIyIysJw8MLcEyLXWzjfGHAvcFH8MY0wOcDMwB2gDqowxfwGaAay1JzuYS6RH82YMdzuC\niIhjnBzeXwA8CmCtXUWswO81DdhgrW201nYBzwEnArOAAmPMY8aYJ+JfFkRERGQAOFn0i4GmhOVw\nfMh/b1tjQlszUAK0Ajdaa08FLgR+n7CNiIiIHAYnh/ebgMRJuL3W2kj858b92oqAPcA6YAOAtXa9\nMWY3MALY2tsLBQKa6zudpUL/+fJyyMnypkSWVBX2xr5/lw0tJFAau5RR71d6U/9lHieLfhWwGLjf\nGDMPWJPQVg1MNsaUEtu7XwjcCJxH7MS/i40xI4mNCGw/2AvV1TUPcHRJlkCgKCX6L9jRRSjLmxJZ\nUtXuhnYA6ne34gmFU6bv5NCo/9LX4XxZc3Lo/EGgwxhTRewkvv8yxpxjjDk/fhz/cuAx4N/ACmvt\ndmAFUGyMeQb4A3BewuiAiKSAq+94nm8ur3I7hogcAsf29K21UWDZfg+vS2hfCazcb5sQ8AWnMonI\noSspzOW/PjULW9PAc2u2uR1HRA6B7qUpIn3iy8niyAlDaQ/qjnwi6UpnxouIiGQI7elLRqpv6uDp\nV2ND1MWFuS6nERFJjh739I0x+x+PFxk06puD/PXfm/jPm7U8/WqvV4SKiAwavQ3vX7D3B2PM00nI\nIpJ0Jx01yu0IIiJJ09dj+sWOphARERHH6UQ+ERGRDNHbiXx+Y8xCwLPfz1EAa+0zScgnIiIiA6S3\nor8V+O8D/LyXpr+VtBONRmlo6aS5rdPtKCIiSddj0bfWnpTEHCJJEYlG+cYvdAtZEclMvV6nb4yZ\nBiwFpgJtwJvE7pNfk4RsIo65YPF0KkcUs2bjbrejiIgkTW/X6Z8GPAvkAX8D/gWUAy8aY05KSjoR\nh5QV5zG8rACPx+0k6amjK8xv/v4mVWsPOgmmiKSQ3vb0/wc41Vr7UuKDxph7iM2ad4KTwUQkNRUX\n5DJpVAlPvbyFEWUFLDhyhNuRRKSPertkz7d/wQew1q4GCp2LJCKpbOq4Ur75maNYMHOk21FEpJ96\nK/q9TaWlQVERYXNtM8sfXMsL1TvdjiIifdDX6/Qhfn1+fNnvaCqRJGrtCLGzoZ0RQzWA1R/TKsto\nbO7g5XV1bKpt5qlXtvKJEycwcWSJ29FEpAd9vU5/f1scyCKSdHk5WXg8sKc5yORRQ9yOk1bmzxzJ\n5BFFjKsoYldjO39/vobW9i63Y4lIL3or+p8DbgcmA1XAVdbaPUlJJZIkJ8wayQmzdGz6cCyMv3//\nWP2uy0lE5GB6O6Z/D/AWcAXgA25OSiIRERFxRG97+iOttd8GMMb8E3gtOZFEBlZXKEIkGjslRWeg\nikgm663od9+c3FrbZYwJJiGPyID7/ePreOa1bQCcffJEl9OIiLint6KvnSIZNGaML2NXQ7vbMURE\nXNVb0Z9hjHknYXlkwnLUWjvBwVwiA8qfn0Nzq2bWE5HM1lvRn3I4T2yM8QLLgZlAEFhqrd2Y0L4Y\nuIbYTYDuttbeldBWDrwEfMhau+5wcoiIiEhMb1PrbjrM514C5Fpr5xtjjiV2v/4lAMaYHGJXA8wh\nNntflTHmYWvtznjbr4HWw3x9EUmyXY0dbNnZwohhBWR5e7s4SETc4OSncgHwKIC1dhWxAr/XNGCD\ntbbRWtsFPAcsjLfdCPwS0PRdImnmd/9Yx7V3r6a1o7e7eIuIW5ws+sVAU8JyOD7kv7etMaGtGSgx\nxpwL1Flr/xF/XCcTiqSJn168gKs/fzQAV/7yea65a5XLiURkf70d0z9cTUBRwrLXWhuJ/9y4X1sR\n0ABcCkSNMacAs4H/NcZ83Fq7o7cXCgSKemuWFOdU/z350rvc/fAbtHV0Me/IEWRnZ+Ev9AEwZEiB\n/t0MgMT3MABUlBfx3fM+wGsbdrHq9e16j1Oc+ifzOFn0q4DFwP3GmHnAmoS2amCyMaaU2LH7hcCN\n1to/713BGPMk8NWDFXyAurrmAQ0uyRMIFDnWf7vrW4lEo3z5Y9MoK87jd1stLa2x2000NLRRV5fr\nyOtmip76bkKFn207mgmHo/pspjAnP3virMP5suZk0X8QWGSMqYovn2eMOQfwW2vvNMZcDjxG7BDD\nCmutjuHLgMvLzeID0yrcjiEikhIcK/rW2iiwbL+H1yW0rwRW9rL9yQ5FE5Ek+eu/N/HCW7HBukVz\nx3DCTE1uJOImXVMjIo5paAni8XhoC4ZoadO0uyJuc3J4XyRpbv5/r/LOttjFIhefeSRTx5W6nEj2\nqijNZ3eTLsQRSQXa05dBoSMYZsb4MtqCIcKRqNtxRERSkoq+DBqjAn68Hu1Rioj0REVfREQkQ6jo\ni4iIZAgVfRERkQyhoi8iIpIhVPRFxBEtHV3Ymga3Y4hIAl2nLyIDLjAkj1kThwIwYUQxu5s6XE4k\nIqCiLyIOMGNLMWPfu0HSqrcOOm+WiCSBhvclozz07DtuRxARcY329CVjfOqDk2gPhgAYMbTA5TQi\nIsmnoi8ZY3plmdsRRERcpeF9ERGRDKGiLyIikiE0vC+DTn1zB79/fB1b61rcjiIiklK0py+DTkt7\nF0+8tIXcnCxmjNdxfBGRvbSnL4NO1dpaAC5YPJ2CvByX04iIpA4VfRlUPnj0aKLRKNPHlZKVpYEs\nEZFEKvoyqJxzymS3I4iIpCztComIiGQIFX0REZEM4djwvjHGCywHZgJBYKm1dmNC+2LgGiAE3G2t\nvcsYkwUcMBY3AAAXfklEQVTcCUwBosCF1to3nMooIsn37s4WfvHAWgDKS/O5/NOzXU4kkjmc3NNf\nAuRaa+cDVwE37W0wxuQANwOLgBOBC4wx5cBiIGKtPR74LnC9g/lEJMl+/LuXWLHyTXY2tDOmwq8p\nd0WSzMmivwB4FMBauwqYk9A2DdhgrW201nYBzwELrbUPAV+Nr1MJ7HEwn4gk2botjYwp97Pk+PFM\nHj3E7TgiGcfJol8MNCUsh+ND/nvbGhPamoESAGtt2BhzL/Bz4D4H84lIEj1ctQmA+UeO4Izjx5Od\n5XE3kEgGcvKSvSagKGHZa62NxH9u3K+tiIS9emvtucaYK4FVxphp1tr23l4oECjqrVlS3ED0X05O\nFoWFufq3kGR9fb/PW3wELW1dABwxaSilRXn4/XlkZ3vVZy7Se595nCz6VcSO0d9vjJkHrEloqwYm\nG2NKgVZgIXCjMeYLwGhr7Y+AdiAS/1+v6uqaBzq7JEkgUDQg/dfVFaa1tVP/FpKoP303uiwfyvIB\nCHV0UdfRRUtLB6FQRH3mkoH67EnyHc6XNSeL/oPAImNMVXz5PGPMOYDfWnunMeZy4DFihxhWWGu3\nG2P+BNxrjHkayAEus9YGHcwoaeo/b9Ry58o3AThqcsDlNCIi6cGxom+tjQLL9nt4XUL7SmDlftu0\nA592KpMMHlGgwJeNGVtKNBp1O46ISFrQzXkkbeVkeykr9rkdQw5DfXOQ2x9Yy+Mvvut2FJGMoHvv\ni4grRgf8HH/kCN7cVE9ebpbbcUQygvb0Je21B0N0doXdjiH9NGXMED63aAoTRha7HUUkY6joi6s2\nbW/C1uzB1uyhqbXzkJ6juqaBmp0tA5xMRGTw0fC+uOo3f3+TF97cAcBXz5jBsdMr+rX9JxZO4PT5\nlQD4sjVELCLSG+3pi+s+PHcMxQU5h7RtXm42xQW5FBfk4tNx4bRVXbOH2/68hjUbd7kdRWRQ056+\npAaPbsmaqaZXlpGfm83zb9QyvbLM7Tgig5qKvqSsto4Qwa4w//fPdexsaCfL6+HskyYxcVQJ4bCu\nzR8sjpsxnONmDOetzZpfS8RpKvqScppaO9la18IDz77Nxq2xOZv8+Tl0dIZ5p7aJG/7vFQCG+HPd\njCkiknZU9CXlrN/SwC8efB2AEUMLuPSsmRTkZXPVr5/vXudrZx1JaZFuzCMi0h8q+pKSCvOyue3r\nC9/3+CP/qQGgcnixir6ISD+p6Eva+MyHJhOOxI7lF/j0T3cwCoUj3PHXN7qXl54+Ha9O8hQZMPrL\nKWnjhJkj3Y4gDotEo/znjR2MH1HEO9ubWXr6dLcjiQwquk5fRFLOnKnlbkcQGZRU9EVERDKEir6I\niEiGUNEXERHJECr6IiIiGUJFX0RSxubaZrcjiAxqKvqSMkLhCDU7mqlr6HA7iriguDCX9VsaKS3y\nke3VnyYRJ+g6fUkZTW2dXHfPC0DsjnySWa4456jun21NbPKddTUN5PuyGTe8yK1YIoOK/rJKyvne\nuXMJDMl3O4akgBv+7xXGVRTxvfPmuh1FZFDQGJqknHxfFgXa089oU8YM4a4rT+bskye6HUVkUHHs\nL6sxxgssB2YCQWCptXZjQvti4BogBNxtrb3LGJMD3A2MA3zAD6y1f3Uqo4ikJo/HgweI/b+IDBQn\n9/SXALnW2vnAVcBNexvixf1mYBFwInCBMaYc+BxQZ61dCHwEuN3BfCIiIhnFyaK/AHgUwFq7CpiT\n0DYN2GCtbbTWdgHPAQuB+4FrE7KFHMwnLqpau50b/+8VqjftcTuKiEjGcLLoFwNNCcvh+JD/3rbG\nhLZmoMRa22qtbTHGFBH7AvAdB/OJi3Y3drBjTxsfnV/J9Moyt+OIiGQEJ8+WagISr7PxWmsj8Z8b\n92srAvYAGGPGAA8Av7DW/qEvLxQI6HKedFNQ6GPEMD+fP20aAPc+Wo2/0AdAWZmfwLBCN+NJHzn9\n2fP7fWTnePUZd4je18zjZNGvAhYD9xtj5gFrEtqqgcnGmFKgldjQ/o3GmArgH8BF1ton+/pCdXW6\ni1e6aWsN0tUZO3pTV9dMJBKlpTUIQH19C9nRSG+bSwoIBIoc/+y1tAQJdUX0GXdAMvpPnHE4X9ac\nLPoPAouMMVXx5fOMMecAfmvtncaYy4HHiB1iWGGt3W6MuRUoAa41xuw9tn+atVa3aBskOrvCNLV1\n0hbU6RoiIsnmWNG31kaBZfs9vC6hfSWwcr9tLgMucyqTuK+6poGf3f8aAFNGl7icRkQks+gOKJJ0\n2Vkerj9/HtlZujeUHNzupg7ufaQaM2YIkWiU9vgo0bHTKygqyHU5nUh6UdGXpPN4PLrNrvTJiKEF\nzBhfxjOvbeOZ17bt02bGlqroi/STir6IpKxZk4Yxa9Iwjp1WQW19GwDjRxTxk/tecTmZSHpS0ReR\nlDd78rDun1vauwC44b6X8eVm8dOLFrgVSyTt6KCqiKQVX04WXz1jBifOHkVTa6fbcUTSioq+iKSV\nnGwvx06vYOrYIW5HEUk7KvoiIiIZQkVfREQkQ6joS1KEwhHWb2lg265Wt6PIIBIKR7n01me7b/gk\nIr3T2fuSFO3BED/63ctA7JisyOEaU+5n2ZIjeMnupL456HYckbSgoi9J9b1z5zIqoBn05PCV+H3M\nnVpO7e5WFX2RPtIulziusytMZ1ds1rysLI9uvysDbuPWRpbd/DT3P7XB7SgiKU17+uK47961il2N\nmihRnDF3WgVjyov46783EQpF3Y4jktJU9GXAdXSGaA+GAfDlxPbqT59fyexJw3TPfRlww8sKGF5W\n8L5784vI+6noy4B78pWt3P/kRgBOPmoUAIEheUwYWexmLBGRjKeiL44YFSiktMjndgwREUmgoi+H\n5Nk129i4tQmAedMrmDqudJ/2nCwvebn65yXJ9a+Xt/Dc2m2ce9o05k4tdzuOSMrRadRySNbVNPDm\npnpWv7WDbbt1wx1x38fmj+OiJUfgwUN7MER7MEQkqhP7RBKp6MshM2OGaAhfUsbEkSUcNSVATo6X\nex+p5uJbnqFuT7vbsURSioq+iAwql39qNl8/e5bbMURSkg66yoBpautkXU0DW3a2dD+2c08bwa6w\ni6kk04wp9+PPzwFg++42usIRRgf8LqcSSQ3a05cBs62uleUPvc6ajbvJ92WTl5vF5h0tRKOxE/tE\nku3nf17TPeeDiGhPXw4iGo1yxS//3b38w/PnkZuT1es2t319odOxRHpV4s/l55edwJqNu/j94+vd\njiOSMhwv+sYYL7AcmAkEgaXW2o0J7YuBa4AQcLe19q6EtmOBH1trT3Y6p/SsvinI3KnlvFC9k988\nZikrznPkdXTSlQwUr8eDPz9Hl42K7CcZY65LgFxr7XzgKuCmvQ3GmBzgZmARcCJwgTGmPN72LeBO\nQKeHp4Ajxpdx3IwKtu9u5eV1dY68xlOv6jaqMvCi0SitHV20B0NuRxFxXTK+Bi8AHgWw1q4yxsxJ\naJsGbLDWNgIYY54DFgJ/AjYAnwB+m4SMchDlpfmcMGskT7y0hSdf2Trgz/+zrx0/4M8pAtDRGeZr\nP3uW6ZWlfPMzR7kdR8RVydjTLwaaEpbD8SH/vW2NCW3NQAmAtfYBYkP+IiKHZOrYUq47by4nzBzh\ndhSRlJCMPf0moChh2WutjcR/btyvrQjY098XCASKDr6SHJJo/I5mJUMKCASKKPL7yMry4suLXRKV\nne3F788jECiitjEI9L8/1H/pKx36btyYUt6oaaC5PZQWeZNJ70fmSUbRrwIWA/cbY+YBaxLaqoHJ\nxphSoJXY0P6N/X2BurrmgcgpB7C36Dc2tFFX10xzS5BwOEKwowuAUChCS0sHdXXNNDS0Af3rj0Cg\nSP2XptKp79raOunsCqVN3mRIp/6TfR3Ol7VkFP0HgUXGmKr48nnGmHMAv7X2TmPM5cBjxA41rLDW\nbt9ve908W0QO24YtjVyx/N98YHo5p84dS74vm5xs3T9CMovjRd9aGwWW7ffwuoT2lcDKHrbdBMx3\nLJyIZITZk4dRVpzHbx+zPPKfGh75Tw0Xn3kkx5iA29FEkkoXsUqPqjfvobXj/edS1u5uY1dDu6Yu\nlbQxcWQJE0eWcOSEMoJdEa7/zYs88MxGHnuhhkvOPJLiwly3I4okhYq+9OiPT25ga10r+b5sPB4P\nANMrS/nKx6YBEBiSz9vbm3p7CpGUMqwkH4AzFoxnT3OQx198l1A4cpCtRAYPFX3p1ScWTuAjx47t\nXh4xtJARQwtdTCRy+D5y7Fh27Gnj8RffdTuKSFLpLBYRyWjfv/cFrrt7tdsxRJJCRV8OWyQSJRSO\ndF/eJ5IOigtyuWDxdI6ZWk5dYwdVa7fz1uZ+3yZEJK1oeD/DXXrrs93HNG9YNr97HvL+uO+f67nv\nn+s5++SJAx1PxDH5vmzmzRhObk4Wr23Yxe8eX8e0saVMG1fqdjQRx2hPP8N1dIY5ZkqAjs7wIW1/\nweIZfOMzswc4lUjyHD0lwE8vWsDCmSPp6Ayxta6FprZOt2OJOEJFXxgZOPQT88YNL2LSyJIBTCPi\nnuqaBq5ZsZqnHZhUSiQVaHg/Q63ZuJvmts4BPQ5fs6NlwJ5LJNnOXDiej80fx633rzn4yiJpSkU/\nQ/3t+U1sqWthiD+XnKzDHPDxQElhLtWb91Cim5xImsrLzSYvF7KyPESicO2K1XSFI2Rnebj6c8dQ\nkKc/l5L+9K94EHq46h227IztdZ967Fgm9jD8fuoHxnLGgvHd6x4qX04Wt3zt+MN6DpFU8pfn3gFg\nbIWfmh0tRHRligwSKvr7CYUj7KiPzRbn9XrS8kY0699toLm9i227Wjl2+vB+b//7f6zj6de2EQpH\nOHZahQMJRVLXFz5s6OiM3X46Eonyk/tecTmRyMBR0d/PnuYg16yI3aijqCCHWy89weVEh+bICUPZ\n1dBxSNuGo1EmjCjiY/MrGV5WMMDJRFLbmHJ/98/vxkfBqtZup8Sfy7xD+BItkkpU9HtwxoJKnszg\nM3hL/D6OnDDU7RgirsrO8lBRVsAj/9mM1+uhorSAwrxsykv1ZVjSky7Z64EvN6vHtrqGdn78u5f4\n8e9eYvlDrycxlYgk04ihhfzognl84sSJNLR08j//+yIPPPO227FEDpn29A9BZ1eYdVsamVFZyjvb\nNMucyGA3/4jhfGBaOff9cz0v2Z18+cf/AiA7y8s3PzObKWOGuJxQpG9U9A/D9PFl1NZn7iEAkUyR\nneUlO8vLcTOGM7bcTyQSpSsc4c9Pv605JyStqOiLiPTRtHHv3Zs/Eo3y56ff5o6/vonX4yEciVBe\nWsBVnzva5ZQiPVPRTwHPv1HLP16Izet9zJQAp8+vdDeQiByUB/j8h6fg9XrI8nqo3tzAS3YnP/jN\ni0SjUU6aPYqJo0oYOSz9LvuVwUtFPwU0tXbS2BKkqCCX3U29X2a3bVcrbR2xa4jLy/IpLtAd8ETc\n4PF4+ODRo7uXRwf8DB9awBMvvovH6+GeR6qB2CWA5y+ezuiAv6enEkkaFf0UUVyY26dr4v/01EZe\n3bBrn8dK/Lm0B0MsnDWShuYgDa2aIUwk2caPKGb8iGIWx0fqqtZup7G1kz89tZGuUMTdcCJxKvpp\n6EPHjGaOCdAWDPHGO/U0tHTSHgzxxEtbyMvNxpfjJS9+yeEvHlwLwNBiH81tXcybUcGuxg72NAfd\n/BVEBr0FR44g2BnmT09t5PEX3qW4MJczFlTSGh+py8vNokgjdZJkKvppyOvxYMbGTiY6anKgx/Uq\nhxfT3NbJ29ua2NMcJBKN8mJ1HSX+XPJys993aODldXU88dIW6ps6mF5Z5ujvIJIJvF6YNXEob23e\nQ2NrZ/e5OwDHzajA6/HETg4AzjxhAmXFeS4llUzhWNE3xniB5cBMIAgstdZuTGhfDFwDhIC7rbV3\nHWwb6Z8Z42OFe96Mvt06tLWji/qmDhYvGM/wsnwno4lkhJzsLC47exZNrZ3YdxsACIUirHprB8+/\nsWOfdV9/u54RQwv47KIpbK1rBWB4WQHjhhclPbcMXk7u6S8Bcq21840xxwI3xR/DGJMD3AzMAdqA\nKmPMw8DxgO9A20hy5Puy+fDcMW7HEBlUigtzmTu1vHt5+vgygl1hIDZL5TOvbaN2dxv/ebOWa+Nz\nf+xV4MtmZKCQI+Kjbzv2tDFr0jDaOkKMH1FMcWEupUW+5P0yktacLPoLgEcBrLWrjDFzEtqmARus\ntY0AxpjngIXAccAjPWwjIjIolBTue2ht8fxKttS1kOeLnYtTMSSfkYFC6puCPPjM20QjUd7YVM/m\n2mY6Q5H3jRIA5GZ7mTCyGI/Hg9cTu7rA6/UQjcbuIjpr0jCysjw0tXYyeXQJgcYOmhrb8XhixxfK\ninyU+H3EF4kdedi7EDsKsXddSV9OFv1iIPEetWFjjNdaG4m3NSa0NQMlB9nmgL543aNEIgN3R6xw\nwnO1tHXxX7c91+s69c0dB1ynPzq6wlSUxobTn3+9llfX7+px3daOLk4+anSP7Yfj4apNFOXnOPLc\nItK70QE/X/iwed/jC2eNfN9jXaEwoXCUTdubCEWivFC9k6HFeUSjUSLRKNFo7OZB0QisemsHpUU+\nVr21g9r6NoKd4QHJ6+n+v9iXg8QvC92PvnfKQvyLgydhGwh2hYlGYzOaOi1pX1eS8MXo998/7ZC3\ndbLoNwGJB6MSi3fjfm1FQMNBtjmg31z3Ecfe4S+efsSArJPqAoEi/nrTx119fUlP6jt3jR0dO6H3\ng8dW9rjOxUnKIunByVn2qoCPAhhj5gFrEtqqgcnGmFJjTC6xof1/H2QbEREROQwepyaLMMZ4eO9M\nfIDzgGMAv7X2TmPM6cC1xL54rLDW/vJA21hr1zkSUEREJMM4VvRFREQktTg5vC8iIiIpREVfREQk\nQ6joi4iIZAgVfRERkQwxaCfcMcZ8EDjHWnu+21mkb4wx84EL4ouX7b1jo6QPfe7SkzHmQ8CngQLg\nBmutLpdOI8aYY4BLiN2D6FvW2p09rTso9/SNMROB2YCmrEov5xMr+iuI/QGSNKLPXVrLt9ZeAPwU\n+LDbYaTffMDXgb8Ru519jwZl0bfWbrTW3ux2Dum3LGttJ7AdGOF2GOkffe7Sl7V2pTGmELgUuNfl\nONJP1tp/A9OBbwKv9rZu2gzvx2fd+7G19uSepuA1xvwPMAlYZq1tcDGu7Kcv/Qe0xe/QOBKodS+t\n7K+P/ScpqI9/O4cBNwDXWmt7nvxDkq6P/TcXeBE4DfgecFlPz5cWe/rGmG8BdxIbwoCEaXuBq4hN\nwYu19hpr7Tkq+Kmlr/0H3AH8mtgw/2+TnVMOrB/9JymmH313E1AB/MgYc1bSg8oB9aP//MDdwI3A\n73t7znTZ098AfIL3CsHx9Dxtbzdr7ReSE08Ook/9Z619mdjtmiW19Ovzp89dSunrZ+9L7sSTg+hr\n/z0JPNmXJ0yLPX1r7QNAKOGhIg4wBW9yU0lfqf/Sm/ovfanv0psT/Zeund3vKXglpaj/0pv6L32p\n79LbYfdfuhZ9TcGb3tR/6U39l77Ud+ntsPsvXY7p77V3SsAHgUXGmKr4so4Dpwf1X3pT/6Uv9V16\nG7D+09S6IiIiGSJdh/dFRESkn1T0RUREMoSKvoiISIZQ0RcREckQKvoiIiIZQkVfREQkQ6joi4iI\nZAgVfRERkQyhoi8iIpIhVPRFREQyRLrde19EXGSMWQLcAGwEvgj4gKeAR4CfWWs3updORA5G994X\nkX4xxnwZ+Li19uPGmPmAx1pbdbDtRMR9Kvoi0i/GGD9QA5wN5FtrV7ocSUT6SEVfRPrNGHMnkGOt\nPdftLCLSdzqRT0T6xRhTANQDC9zOIiL9o6IvIn1mjPEC5wPfBVqMMSe7HElE+kFFX0T6xBjjAS4E\n7rXWdgF3AUvdTSUi/aGiLyIHZYz5KPAoMN9a2xh/eCRwpjHmv9xLJiL9oRP5REREMoT29EVERDKE\nir6IiEiGUNEXERHJECr6IiIiGUJFX0REJEOo6IuIiGQIFX0REZEM8f8BPKGumRV9+5YAAAAASUVO\nRK5CYII=\n",
"text": [
"<matplotlib.figure.Figure at 0x91d6e48>"
]
}
],
"prompt_number": 3
},
{
"cell_type": "heading",
"level": 2,
"metadata": {},
"source": [
"Enter ~~scipy~~ the confusion"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can fit distributions to samples using scipy. The problem is that \n",
"the terminology is pretty ambiguous. Scipy needs three parameters to\n",
"fully define a lognormal distribution. Here they are with `sp` as a \n",
"prefix to differentiate them from Wikipedia's parameters:\n",
"\n",
" + `spShape`\n",
" + `spLoc`\n",
" + `spScale`.\n",
"\n",
"You might be thinking: `spShape` is weird, but the rest make sense."
]
},
{
"cell_type": "heading",
"level": 2,
"metadata": {},
"source": [
"__They don't make sense__."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
" + `spShape` actually corresponds to `wikiScale`.\n",
" + Remember when we said that lognormal samples can't \n",
" take values at or below 0? Scipy disagrees. `spLoc` is \n",
" used to translate the PDF of $X$ anywhere along the x-axis.\n",
" In other words, `spLoc` (almost) *always* should be set to 0\n",
" + `spScale` doesn't directly coorespond to `wikiLoc`. It's\n",
" actually $\\log(\\mathtt{wikiLoc})$\n",
" \n",
"The take away here is that if you have real-world sample that you\n",
"think may be best described by a lognormal distribution, don't\n",
"just throw it into `scipy.stats.lognorm.fit` and you can't just\n",
"take the returned parameters at face value."
]
},
{
"cell_type": "heading",
"level": 2,
"metadata": {},
"source": [
"Fitting lognormal distributions with scipy"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"When using the `.fit` method of `scipy.stats` distributions, you\n",
"can guide the parameters by specifying a start/guess value. For us,\n",
"this means doing:\n",
"\n",
"`stats.lognorm.fit(mydata, loc=0)`"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"spShape, spLoc, spScale = stats.lognorm.fit(X_n, loc=0)\n",
"ln_dist = stats.lognorm(spShape, spLoc, spScale)\n",
"print(\"\"\"\n",
"spShape = {:.3f}\n",
"spLoc = {:.3f}\n",
"spScale = {:.3f}\n",
"\"\"\".format(spShape, spLoc, spScale))"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n",
"spShape = 0.866\n",
"spLoc = -0.044\n",
"spScale = 9.505\n",
"\n"
]
}
],
"prompt_number": 4
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"So as we expected, this doesn't really makes \n",
"sense at facevalue. But information we can use \n",
"is in there somewhere."
]
},
{
"cell_type": "heading",
"level": 3,
"metadata": {},
"source": [
"Another graph, this one with the theoretical distribution we just fit"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"fig, ax = pyplot.subplots(figsize=(8, 5))\n",
"_ = ax.hist(X_n, bins=bins, cumulative=False, histtype='step', \n",
" normed=True, lw=1.25, label='Numpy random number')\n",
"\n",
"# theoretical PDF\n",
"x_hat = numpy.logspace(-1, 2, num=1000)\n",
"y_hat = ln_dist.pdf(x_hat)\n",
"\n",
"ax.plot(x_hat, y_hat, '-', lw=1.25, label='Fit with scipy')\n",
"ax.set_xscale('log')\n",
"ax.set_xlabel(r'$X$')\n",
"ax.set_ylabel('PDF')\n",
"ax.legend()"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 5,
"text": [
"<matplotlib.legend.Legend at 0xa766358>"
]
},
{
"metadata": {},
"output_type": "display_data",
"png": "iVBORw0KGgoAAAANSUhEUgAAAf0AAAFPCAYAAABHzhf2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xd4HNW9//H3bFVb9VW1VSzLR7JsYxsDLtjGYEpITMlN\nyCUBQhLKJZUEkgskhJCQQghwc4PJJZQAuVzID0IJppiEYMA2zdjGVceWq6zee9n2+0OykZvkot3Z\n1X5fz6PH2jk7Ox9rtPruOTNzxggEAgghhBBi7LOYHUAIIYQQoSFFXwghhIgSUvSFEEKIKCFFXwgh\nhIgSUvSFEEKIKCFFXwghhIgStmC9sFLKAjwITAP6gGu01juGtC8Bbge8wGNa60cG13kEmAT4gWu1\n1jpYGYUQQohoEsye/iWAQ2s9F7gFuHd/g1LKDtwHnAssBK5TSmUA5wHxWuszgZ8DvwxiPiGEECKq\nBLPozwNeB9BafwDMGtJWClRordu01h5gJbAA6AGSlFIGkAT0BzGfEEIIEVWCWfQTgfYhj32Dw/f7\n29qGtHUwUORXAjFAOfAQ8Icg5hNCCCGiSjCLfjvgGrotrbV/8Pu2Q9pcQCvwn8AqrbUCpgNPKKUc\nw20kMDCPsHzJV0R+VdZ1BJbc9FLg+be2B772izdMzyNf8iVfEfF1woJ2Ih+wClgCPKuUmg1sGNJW\nDhQrpVKALgaG9n8HzOTT0YEWwA5Yh9uIYRg0NHSMcnQRKm63K6r3X3NzFwCdXX34/f6I+llE+76L\ndLL/Ipfb7Rr5SUcRzKL/AnCuUmrV4OOvKaUuBxK01g8rpX4ALGdgtOFRrXW1Uuoe4M9KqXcZKPi3\naq17gphRCCGEiBpBK/pa6wBwwyGLtw1pXwYsO2SdVuDSYGUSQgghoplMziOEEEJECSn6QgghRJSQ\noi+EEEJECSn6QgghRJSQoi+EEEJECSn6QgghRJSQoi+EEFGipqaa885byHe+cz1XXnkl3/nO9Tz+\n+CNs376Nxx9/BIC3336LxsbGY37N5uYm7r33bgDWr1/Ljh0VAFx00fmjkvl///dxtm7dPCqvJYI7\nOY8QQohhdPZ30efrG7XXc1qdJDjih31OYeEE/vCHhw6bka+4eBIAzz33DIWFhUD6MW0zNTWNm276\nTwBeeeXvLF58PkVFEzGME/s/HOqKK64enRcSgBR9IYQwhcfn4fbVv6Lf7xm113RY7Px2/s+wW+3H\ntd7atWt46aXnueCCC9m+fRt33fUzHnzwEWy2gRLxjW9cyb33/oGEhAQuvPAcli79E8XFiq9//Qp+\n9rNfctddd3DTTbfwwQfvsX37NgoKCunv93DnnT+hrq6WpKQkfvGLuw+8HsDzzz/L66+/gsVioaRk\nMjfeeDOVlXu5++678Hq9OJ0x3Hnnr1i69L9YvPh8mpubeP/9VbS2ttHW1srXv34dhYUT+PnPb+fh\nh58A4Kc/vZXLL7+C0tKy0fqRjjlS9IUQwgR2q51fzL1t1Hv6IxX83bt38p3vXI/dbsXj8XHHHXcd\naJsz50yKiyfxwx/edlCBnj9/IR98sBq3O4OcnFw++ugDbDY748fn4XA4MAwDpUqYPXsuixefT2Zm\nFj093Vx//bfJysriO9+5nu3b9UHF+LXXXuamm26lpKSUF198Dp/Px9Kl/8VVV32d00+fzcqV77B9\neznGkCEDvz/A73//IE1NjVx//df4619fxOl0snv3LlJTU6mpqZaCPwIp+kIIYZIERzwJDD8cP9oK\nCg4f3t+7d8+w6yxYsIgnnniUrKxsrrvumzz33DP4/X4WLTrnqOskJiaRlZUFDBwC6O3tPaj91lvv\n4Jln/pfq6iqmTJlGIBCgsnIvU6ZMBeDMMxcA8I9/LD+wzqmnngZAWlo6CQku2tpaueiiS3n11ZfJ\nzMziggsuPM6fRvSRE/mEEEIcYLFY8Pv9By2bMKGI6uoqysu3MGfOPLq7u1m58h3mzJnHwN3NBxiG\ngc/nG/x++O28/PKL3HzzrTzwwJ/Ytk2zadMG8vML2bJl4KS9f/zjdf72t/930Drl5VuAgZMHe3t7\nSElJ5ayzzuHDD9/nnXdWcN55UvRHIj19IYSIIsYRqrFhGAeWT5kyjbvuuoP771+Ky/XpLVxnzpxF\nbW01hmEwY8ap7N69C6cz5qDXnDx5Cg89tJScnFzAOGwbQxUVFfGtb11DXFw8bncGZWVT+da3Mvjt\nb3/FE088SmxsLD/5yc/ReuuBdfbtq+R73/sm3d2d3HzzrRiGgcPhYPr0mbS1tR6UVxyZMfRTWoQK\nyD2hI1e039O7urGLnzzyAV9cVMSbH+/jd9+cZ3akYxbt+y7SRdr+e+21ZbS2tnL55Vcc1nb//b9l\n4cKzmTlzlgnJQs/tdp3wtREyvC+EECIiHOmQwQ9+8G06OjqipuCfLBneF0IIEfY+85nPHXH5ffc9\nEOIkkU16+kIIIUSUkKIvhBBCRAkZ3hdCCBP4/H5a2kdvYp79UhKdWC3SnxNHJkVfiDGkqrGL3TXt\nAGSmxjExN8nkROJoWtr7+NH/vDfqr/vb/5hDenLsqL+uGBuk6AsRIbZVtvLh1joACrMTmTc1+7Dn\nbN7VzLNvVWC1GMwuy5SiHwFu+cpMUl3Ok36d5o4+fvPU2mGfs3btGm677WaefPKvuN0D17T/8Y9/\noKCg8KgnyoW7q676Ek8++VezYxxm7do1vP32v/j+939kdpSDSNEXIkJUN3axamMtCbF2uvu8Ryz6\nADnp8eSmh3ZqV3HiUl3OkPbM7XYHv/rVnTz11F+AI0/WI05euP5cpegLEUGSXU4mjUvC4/OP/GQh\nDmEYxuD17AGeeuopzjvvogNttbU13HHHbTz00J8BuP76r3Hnnb/ilVf+TnX1Plpb22hvb+Xzn7+M\nFSvepLJyLz/+8Z2kpqbyy1/+jNjYWJqaGpk7dz5f//p1fPnL/8bDDz+Jy+XihReeo6enmy9/+aoD\n27vyysvIy8vHZrPz7W/fyO9+92v6+/tpamrk2mtvYP78s/jqV/+dGTNOpaJiO4Zh8Jvf3EtsbBz3\n3PNrduzYTkZGJl1dXQDU1FTz61///MAUwjfe+EMmTizmS1+6hKlTT6Gyci+nnnoaXV2dbNmymby8\nfG6//ecH/XyOtD2ty3nppee5885fAXDxxefz0kvL+eUvf4bNZqeurob+/n4WLz6PVavepa6ull//\n+l4Atm3T3HjjN+nq6uTSS7/IhRcuYceOCn7/+98RCARISkri1lt/itbl/PGPf8DhcHDRRZdy/vnB\nm05YzvYQQogosX8G1ptuuoXHH3+cqqp9I65jGAZOZwz33vvfLFx4Nu+9t4q7776fK664mjffXI5h\nGNTV1fKLX9zNww8/yUcffUBFxTbOPfcC/vnPgZvlvPHGa3zmM0sOet3e3l6uvvpa7rzzV+zZs5t/\n//cruP/+pfzoRz/m+eefBaC7u5vFiy/ggQf+hNudwfvvr+bdd1fQ19fLn/70ODfffCtdXZ0ALF36\nX1x22Zd54IE/8b3v3cxvfvMLYODDzHXXfZOlSx/muef+yuc/fxkPP/wEGzZ8cmDd/Y60vaP12A3D\nICcnh/vue4CCgkJqamq4557fs3Dh2axa9S4wcB+D++9fygMP/Im//OXPtLa2cvfdd3HTTbfwhz88\nxOzZ83jqqScxDAOPx8PSpQ8HteCD9PSFCEsfbq2jrasfgFMnuUlNjDE5kRhLEhOTuO2227jrrjuY\nOvWUIz5n6BTtkyaVAJCQ4KKwcMKB7/v7B35HJ0+eQkxMzIHvKyv38tnPXszPfnYbp5wyg9TUVFJS\nUg7bRl5ePjBwF74nn3yMZcteOuimPQPbVgBkZGTS399PdXUVJSWTAUhOTiY/vxCAPXt2M336TACK\niydRXz9w/ktSUjIZGZkAxMbGkJ9fMJg/nv7+fuIPORJ26PYO/7kMfe6nP5f9r+tyJdLfP3BVxrRp\n0w98aCooKKS2tpq9e3fzu9/9GgCv18v48XkH/SyCTYq+EGHojY8qaWzrpb2rn9z0eCn6Y1hzx+hc\ntne8r7No0SJefvlVXnttGd/85ndxOBy0tDTj9/vp6uqipqb6CGsFONL9Wnbs2I7X68UwDLZu3cxF\nF11KVlYWCQkJPPnkY3zuc5ccMYNl8NLCRx/9H5YsuZTZs+fyyit/57XXlh14zqE97YKCQv7xj+Vc\ndtnltLe3U1m5F4D8/ELWr1/LmWcuYPt2TVpa2uD6x/VjOWx7DoeTpqZGYGDUoL297Zhfa+vWzQQC\nAXp6etizZze5ueMZP37gsEJGRibr16+lra3tiNsNFin6QoSpc2eN4/l3dpodQwTZSGfcj6ahd9MD\n+N73buLjjz8CBnrbp512BtdccxW5ueMYN278QesNfnfg+4F/Pl3+ox99n/b2NhYvPu/AaMCSJZfy\n+9//jjvuuOtIaQ58t2jRYpYu/S+effYZysqm0NHRftT/w/z5Z7F27cdce+1XSU93k5o6UNy//e0b\nufvuu3jmmf/F6/Vyyy0/PWw7h975b7hM+5WUlOJyubjuuqspKCgcvIPg4LOHGfrf/++NN36L7u5O\nrr32BlwuFzfffCu/+MVP8fl8WCwWbrnldhoa6kNW9IN2lz2llAV4EJgG9AHXaK13DGlfAtwOeIHH\ntNaPKKWuBr46+JRY4BQgU2t99N8AucteRIu0O32NtqPdZe+uJ9cwozid59/ZyU1fms7kglRWrKti\n+UeVB07ku25J2WGv98ZHlazaWENuejwOu4WrP1MatOzRvu9OltmT84zW/qupqeb+++/ht7+9/7C2\nt976Jzt37uAb37j+pLcjPnUyd9kLZk//EsChtZ6rlDoDuHdwGUopO3AfMAvoBlYppf6utX4ceHzw\nOQ8Aj4xQ8IUYU1r72lhTt57mtLWs7PXgmNzP63W7aLaX4QlkmB1PjCKrxTImJtEZGD04fPlDDy1l\n/fqPufvu/wp9KHFUwSz684DXAbTWHyilht73sBSo0Fq3ASilVgILgOcGH88CyrTW3w5iPiHCRn+g\nl/7MT/jp6mWkOJOw+lLIsWbS1NxIX5LB/21aBhYfMSmT8HMaIw9TChEaWVnZ3H334b3866//lglp\nxEiCecleIjC0l+4bHPLf3zb0bIgOYOjUYbcBPwtiNiHChiWhhbd7nsYf38C1U6/kjjk/Iqn1VCbY\np+OrncB05zn0rj+LBemLsWTsYavj7/QZMgAmhDh+wezptwOuIY8tWuv9M4q0HdLmAloAlFLJwCSt\n9dvHuqH900mKyBTN+++92vU4Sj5knHMqNdsK6Cl0s/SdzdQ2dxMf78QAEuKdELDwjQUX8lVjMTe/\n+HsqYl+lxZjIpPQJB71eQoITm82C02nH6bAG/WcbzftuLJD9F32CWfRXAUuAZ5VSs4ENQ9rKgWKl\nVArQxcDQ/j2DbQuAN49nQ3IyUeSK5pPBPq5bz1Nbnsa7bxKlpfOp9u1j+54Wqhu7OGt6LllJMQSA\nzq6Bk70aGzuJddoo6j2HXdbV3LXiD3x3xrXkJ356lnVnZx9er5++Pg8Bvy+oP9to3ndjgey/yHUy\nH9aCObz/AtCrlFrFwEl831dKXa6UulZr7QF+ACwHVgOPaq1rBtebBOw44isKMUbo5goe3/IMZ2cv\nxls7MLmIzx+g1+NlXHo8XziriJL8wyczATCwYK+dTqa1kP9e9wgN3U2hjC6EiGBB6+lrrQPADYcs\n3jakfRmw7JB2tNa/C1YmIcJBY08Tj276XxaOm8vslLks4wMA2jr7WbWxllMnuYddPzHeQU+ll8rV\neeSe0cmfNj7BTad+ixjbyd+pTQgxtsnkPEKEkMfv5eGNf2G8K5dLiz5LXXMvADOL3cwcLPZOu3XY\n1/i3hUX828Iifvroh5wWewEfev/G0/pvfK3sy0HPL4SIbHLDHSFC6JWdb9DW187VZZdjtXxa3GOc\nNjJT4shMiSM54dh77A4jhm9MuYK19Rv4uG59MCILIcYQKfpChEhF6y7+ufdtvlL6BVyOBJ5dUcGr\n7+856dfNSxzHZwrO4Rn9Ar3+rlFIKoQYq6ToCxECPr+PZ/TznJ41k6npA3cI+9fHVdQ0dTF9Yjp2\n68lNtnN+/tmkxaSwqW/laMQVQoxRUvSFCIEV+1bR2tfOpRM/e9DyJfMK+e4XphEXYz/qurVN3SO+\nvtVi5TJ1KVXebXicDSedVwgxNsmJfEIEWWtfG6/seoNLij47MKz/VgWVDZ30e30jrmu1WFi9qRar\n5cgjAT19Xlo6+mjt7OPpfzRgTcijM20dgf485DO9EOJQUvSFCLJlO98gIzadM3PPAGBXTTteX4AL\nTs8jY4Qbrvzph2cN2/7iyl28uHIXX1xUxI7qdi48cxErPU/T5N9GDpNH678ghBgjpCsgRBDVdtXx\nfs0aLp54IRbj07fb5IIUvrhoIjnp8Sf82t/9t6n8/BunH7TskjmlXFB4NnWODfjwnvBrCyHGJunp\nCxFEL+9cTk5MHs+82ILBh0wclzTySscoPTkWr89/2PKzxs3jtYoV1Fu3AlNHbXtCiMgnPX0hgmRP\neyXrGzYxPeFM6lp6iIux0dDaG/TtxticZHimUmPdSI+3J+jbE0JEDin6QgTJa7vfZGp6KW5HNg6b\nlYLsxJBtO82jsGBlReXqkG1TCBH+pOgLEQTVnbVsbNzC+flnm7J9C1ayfFNYsW8l/T6PKRmEEOFH\njukLEQRv7FlBcfIECpPyaaypD/r2tu5pOWyZ2zeJ5sAm3qn8gOWvDny+t1oMfvMfc4KeRwgRnqSn\nL8Qoa+xp5uP69ZyXvygk2yvKSaSn10tRTiLGkMv5rdhZOG4uK6repam9m8KcRJrag39OgRAifElP\nX4hR9ubed8iNz6I0dVLQt2WzWvjxVbOO2r5w3Dz+sfdtrKm1FGYXs26bzNYnRDSTnr4Qo+gv/9zE\nu5Uf0L5n/BGH3EMtwRHPqe4ZWDP3mh1FCBEGpOgLMYoqejZjBOw0V6bQ1Wvu5Dg+X4Defi+zM8/A\n6mqlzSe9fCGinRR9IUaJP+CnLWYb462TcVjtVNZ3sHlXM/2ekefYD4ZVm2r55n3vUFHhx9eeym7v\nBnz+ADfc+zb3PL3OlExCCHPJMX0hjtNdT66hpaMPgB9ePoOs1DgAtjZvw2PpJJvJ1FjaefW9vSxb\nvYcvLioKecZLFkzg3NPGs/SFjQB46/KoTtrIdZeey8ZtHdS3jnznPiHE2CM9fSGOU1tnH8Xjkmjp\n6MPnDxxY/va+1ST05+Ekjvu/cya/uX62aRkzkmMpzE7Ebh14i/tbMoizxdMeU0F2WpxpuYQQ5pKi\nL8QJyM9yHfS4saeZLU2a5N7gn7F/YizMSDmVVVUfEAgERn66EGJMkqIvxCh4r+YjsuIziPG6zY5y\nVFOSp9HU20Kzv9rsKEIIk0jRF+Ik+QN+3q9Zw9zs0zAwDmvfV9954Ps9te3sa+gKZbwDXPZEJqcp\n9vm2mrJ9IYT5pOgLcZK2Nm+jo7+T07JmHrTcMAziY2xs2NFEfIyN4twkTi/NZPbkTCbkhO7mO0PN\nzT6NWv8O/Ea/KdsXQphLzt4X4iStrv6IqemTcTkSDlqelhTDH25ccNCyGZPMHf6fkl6KFTs9sXuB\nuaZmEUKEnhR9IUaw/MO9vPDuTgDOmp57UFuXt4uNjVu4ftpXzYh23GwWG7lWRVX8brOjCCFMIEVf\niBH4/AFSEpy44hx4ff6D2ja3bMDlSCDBm8N7m2pp6egjJy3epKTHZpy1hN2OT1izuwKHLwWA8RkJ\npLicJicTQgRb0Iq+UsoCPAhMA/qAa7TWO4a0LwFuB7zAY1rrRwaX3wosAezAA1rrJ4KVUYhjFRdj\nJynBcdjyT1rWMzvrVNZta+TV9/fginMQ67SakPDYuSxp2PtTeWnLu1RvzCcQgGs/N5k5U7LMjiaE\nCLJgnsh3CeDQWs8FbgHu3d+glLID9wHnAguB65RSGUqps4A5g+ucBUwIYj4hTooR20FDbz2nD57A\nV5STxL3fmseSeYUmJxtZbHce7Y49nD0zV3r4QkSRYBb9ecDrAFrrD4Ch9/8sBSq01m1aaw+wElgA\nnAdsVEq9CLwM/D2I+YQ4Kda0atIdWRj9CSM/2URN7b2HLYvtGY/X6KGdGhMSCSHMEsxj+olA+5DH\nPqWURWvtH2xrG9LWASQB6UA+8FkGevl/B0pG2pDb7RrpKSKMhfv+i493YrdbcDptxMY6sFgtxMc5\nsKbVUL29gDs+/IjPnzURu8Madv8Xq9XCv9ZWAZCSEo/b7SI+3kmMJYEEfzYt9t1YLJNwJcacUPZw\n+/+K4yP7L/oEs+i3A0N/o/YXfBgo+EPbXEAr0ASUa629wDalVK9SKl1r3TjchhoaOkYxtgglt9sV\n9vuvq6sPj8dPX5+XHls/fp+fhr5KrM4+Pj9jPn97s2rgOf2+sPu/3HH1LPbPumuzGTQ0dAxk9fpI\n6M2n1vIxtsAEOtp7jzt7JOw7cXSy/yLXyXxYC+bw/irgQgCl1Gxgw5C2cqBYKZWilHIwMLS/moFh\n/gsG18kB4hn4ICBEWNnbr5mUUkRGQorZUYZlt1lx2Ae+LMbBswW6+sfjx4c/vt6kdEKIUAtmT/8F\n4Fyl1KrBx19TSl0OJGitH1ZK/QBYzsAHj0e11jXAK0qpBUqpDweXf1NrLXcHEWElgJ+9fdu4rOAi\niOA71FpxkEoeLUn7zI4ihAiRoBX9wWJ9wyGLtw1pXwYsO8J6/xmsTEKMBl9CHd6Al+kZU9i2uwu/\nH5rb+8yOdULc/iKaEv6Fxy/T8goRDWTufSGOky9xH7mOCcTaYgHw+vys3BiZZ8GnMA4CFqr6d5od\nRQgRAjIjnxDHob2nG19CHfkxnwFgalEaD//oLIAj3mEv3FmwYenM4t32tby3ysqPrzqVzJQ4s2MJ\nIYJEevpCHIe1tVsJBCDHUQCAxTCwWixYLRYslsgr+gDnFM3CltJIZ18Pfr+cQiPEWCZFX4hj9PUL\nSzl9tp9pGaXMnZw78gphrL2rn6a2gUl7Liw7beBDS1Ijm3c1s257g8nphBDBIkVfiKNo6ejj3Q3V\n7K4ZmGPKaguwtbWc07OnY7eF9/z6w4mLsREIQHefl4RYOw6rHZU8Cae7nr+9vZOHX95idkQhRJDI\nMX0hjqK2qYs/v1pOelIM4zMS2NqsCQQCTE5VZkc7KWfPHMfZM8cdtOz0nOnotv/Hl8uK+fMr201K\nJoQINin6QgzDAH57w1wAHt/8NGVpihjb2LtBzeRURSAQYF/vbrOjCCGCSIb3hTgGHp+HjY1bmOGe\nanaUoIixOSlLU+zq0WZHEUIEkRR9IY5Bect2vAEfU9JLzY4SNDPcU9nbswMMn9lRhBBBIkVfiGOw\nrn4jk1MVMbYYs6METVl6Kf6AD1zD3t9KCBHBpOgLMQKv38uGxs3MyBibQ/v7xdpiyIkpgKRas6MI\nIYJEir4QI9jWsoN+n4epY3hof7/82ImQWI8/4B/5yUKIiCNn7wsxgg2NW5iUUnRgrv2xLC+2CGzL\neWzFarztyVgMg1MmpnPG5EyzowkhRoH09IUYVoCNjVuYlj7Z7CAhEWeNJ9CVxIdVG/lwaz26spW9\n9R1mxxJCjBIp+kIA/R4fPX1eevq8eH2fDm0b8e209rUxNUqK/rSiND5XNpucog7+56aF5KbHmx1J\nCDGKZHhfCOCJ18t5b3MdAFeer1g0Y2BufWtyPXmuXFJiks2MFzJWi4VTs6bw6p7ltPQ3mR1HCDHK\npKcvxKDpE9NJSzz4kjxLSh3T0stMSmSOzLgMMmLT2dAoc/ALMdZI0RdiUEKcHbvt07dEu6cVS1wn\n09zRVfQNw2Bq+mQp+kKMQVL0hTiKnd3b8ffFkhOfZXaUkJvmLmNX2x68Ri/vflLDTx/9gM27ms2O\nJYQ4SVL0hTiKXV3b8bdkYBiG2VFCbkJSPnH2WPKKu1kyr4D6lh56+rxmxxJCnCQp+kIcQaeni+re\nSvwtGWZHMYXFsDAlrZQ2ayXnzhp/0GEPIUTkkneyEIfw+wO8vnUNNhz4O1LMjmOaae4ytjZvo9/X\nb3YUIcQokUv2hDhEv9fHP3etgUAaNlv0vkVKUycRCPjZ1rLD7ChCiFESvX/RhDgKX8CLJamRi/Mv\n4fwvnmF2HNM4rQ6KU4rY1FQOpPFJRSPN7b3MPyWHWKf86RAiEsnwvhCHaPTtA8NPYfwEs6OYriyt\nhE2NW8l1x7Ozpp1n/lVBd6+c0CdEpJKiL8QhPqnfgr89DYfVaXYU001JK6Wlr5WrLs7le188xew4\nQoiTJEVfiCHmTc3CSKpHJU8iMc5hdhzTuePSyIxzs7mp3OwoQohRELQDc0opC/AgMA3oA67RWu8Y\n0r4EuB3wAo9prR8ZXL4WaBt82k6t9TeClVGIQ02f6uTVD7u4cs4CUmNjRl4hCpSllbCpaSunJEXv\n+Q1CjBXBPBvnEsChtZ6rlDoDuHdwGUopO3AfMAvoBlYppV4COgC01ouCmEuIo9rcWE5WXAbpsalm\nRwkbZWklrNi3ih5vj9lRhBAnKZjD+/OA1wG01h8wUOD3KwUqtNZtWmsPsBJYCJwCxCmlliul3hz8\nsCBEyGxqKqcsvcTsGGFlYnIhDoudHe1y6Z4QkS6YPf1EoH3IY59SyqK19g+2tQ1p6wCSgHLgHq31\no0qpYuA1pdSkwXWOyu12jXJ0EUrhsP+cMXYMq4dd7Xv4yoyLwyJTOJmWXUpl307ATWpaPO6UOCA8\n9p04cbL/ok8wi347MPQ3yjKkeLcd0uYCWoBtQAWA1nq7UqoJyAaqhttQQ0PHaGUWIeZ2u8Ji//X1\nemiz78YR6yCdzLDIFE6KE4p5oeJVIJ3mpi4Mry9s9p04MbL/ItfJfFgL5vD+KuBCAKXUbGDDkLZy\noFgplaKUcgALgPeArzFw7B+lVA4DIwI1QcwoxAFt1n2UphZjtVjNjhJ2ytJK6PZ2YcS3jfxkIUTY\nCmbRfwHoVUqtYqCQf18pdblS6trB4/g/AJYDq4FHtdY1wKNAolLqHeAZ4GsjDe0LMRoC+Gm37KMs\nTY7nH0mAtJUaAAAgAElEQVSS00V2XA7W5AZueeg9blq6yuxIQogTELThfa11ALjhkMXbhrQvA5Yd\nso4XuDJYmYQ4mm5LE16jj8lS9I9qanoJHu9mpuXl8e4n1WbHEUKcAJmcRwig3bqPOH86SU45selo\nprkn09hfR1qa2UmEECdKir4QDBT9JP84s2OEtfzEcSTY46nq22V2FCHECZKiL6JeW187PdZmknxS\n9IdjMSyUpZWwT4q+EBFL7o8potKO6jbueXodACn59dhSYogLpJucKvyVpZWwtu5ZAsYks6MIIU6A\n9PRFVAoEoN/j5/SSTHqd1bi8uRgYZscKe6WpxXgDHgJxzWZHEUKcACn6IqplpjnxxNaTKEP7xyTO\nHofbnkMgod7sKEKIEyBFX0S1Zn8NAcOHy5djdpSIMS6mEL9Lir4QkUiKvohq9d7d2HrTsOEwO0rE\nyHUWQkwHDV1NZkcRQhwnKfoiqtX59uDoyTI7RkRJsbnBE8O6ms1mRxFCHCcp+iJqGc5uugKt2Lul\n6B8PwzAwOjJYV7PJ7ChCiOMkl+yJqOLz+7nn6fX09HmxJDcQayRg9SSaHSviWDoz2FS3Ac8kL3aL\n/BkRIlJIT19EnW2VreRnusgu6CTTWiCX6p0AozMdr99LRetOs6MIIY6DFH0Rlc6YkkabUUOWvcDs\nKBHJ8NspcU9kc1O52VGEEMdBir6ISpXdezCAdJtcn38iOno87Nzq5L29G8yOIoQ4DlL0RVTa3V1B\ncUoRNsNudpSIk5/l4orzFBm2fHqNdhq65dI9ISKFFH0RhQLs6tpBWVqJ2UEiUmZKHItm5KIy87D5\n4tjcLEP8QkQKKfoi6hixnXR625mSVmp2lIhmGAb+Njevb1nDx7rB7DhCiGMgRV9EHWtyA6mONNJj\nUwHo9/rp6vWanCryTMpLoTBhIh2WGpa9X8Fflmv2NXSaHUsIMQwp+iLqWJIaKIibeOBxfUsPG3bI\ncenjNX96Lt86dxGGBWLS2lixvorm9l6zYwkhhiGzaoio0u3pxuJqpSC+CICFp+Qwp2xgRj6bVa7X\nP15Oq4OS1IlkjOtlV3mM2XGEECOQnr6IKuUt28FvISd2PAAOu5WEWDsJsXZiHPIZ+ERMTlNsbiwH\nAmZHEUKMQIq+iCqbm8rxt6VjNaxmRxkzytJKaOxtBmeX2VGEECM4atFXSt0QyiBCBMveug7WbWtg\nra5nc6PG1+o2O9KYkhGbTnpsGkaSnMEvRLgbrqd/3f5vlFJvhyCLEEHxr7VVPPjiJh58YyXdvm58\nbVL0R5NhGJSllWAk1psdRQgxgmMd3pfbkImINqskg5TcNpIsbvA4zY4z5pSllUBCCx5/v9lRhBDD\nkGP6Imr0x9aSYc03O8aYVJw8AYCq3j0mJxFCDGe405UTlFILAOOQ7wMAWut3QpBPiFHhMXrwOlsG\ni76ccDbaHFY7dKSxN3EnMMfsOEKIoxiu6FcBdx7h+/0WDffCSikL8CAwDegDrtFa7xjSvgS4HfAC\nj2mtHxnSlgF8DJyjtd52bP8VIY6uw7oPw+cg2ZIB7DI7zpgUaHNTmbqTQCCAYcicB0KEo6MWfa31\nWSf52pcADq31XKXUGcC9g8tQStmB+4BZQDewSin1d611/WDbQ0h3TIyiDlsVjp4sjEQ5ohUsgfYM\nOn1buPXJN/jJZWeTECt3MBQi3Aw7G4lSqhS4BihhoDhvAR7VWu89hteeB7wOoLX+QCk1a0hbKVCh\ntW4b3M5KYAHwHHAP8Efg1uP7rwhxZH78dFirie+ZaXaUMe3fF5zCK83raQ7s459rKkl2OTlreq7Z\nsYQQQxy16CulPgP8Bfgr8AoDx/KnAWuUUpdprVeM8NqJQPuQxz6llEVr7R9saxvS1gEkKaWuBhq0\n1m8opW5l4ByCEbndrmN5mghTwdp/DS097Kltp8FTjd/hJa4/h4T4gTP3k5Pj5PdmFAz9GX7x3BIa\nVk9jdV8572+pO7BMhC95D0Sf4Xr6vwDO11p/PHShUurPDAzVzx/htduBob9R+ws+DBT8oW0uoBX4\nLhBQSi0GpgNPKKUu1lrXDbehhoaOEaKIcOV2u4K2/95eX8UTr2uc+dtxxbrxe6x0dvUB0NraTUOD\nIyjbjRZH2ndT0xTv7FvJkjOyeW11tbw3w1gw33siuE7mw9pwBzidhxZ8AK31h0D8Mbz2KuBCAKXU\nbGDDkLZyoFgplaKUcjAwtL9aa71Qa32W1noRsB64aqSCL8RwMlJiyZnQzeLiU82OEhWKkguxWWzU\neY7lCKAQItSGK/rD3WD8WIbdXwB6lVKrGBgZ+L5S6nKl1LVaaw/wA2A5sJqB8wRqjjW0EMfKb+ui\ntquOsnQZZg4Fu8VGSUox1f1yhYQQ4ehYr9OHT2+hZQAJI72w1joAHDp//7Yh7cuAZcOsP+wlgUIc\nC298HcnOJHLis4BKs+NEhbI0xYvNy4GJvPLebj7aOjA97+JZ4zlzWrap2YSIdsd6nf6h9gUhixCj\nzhtfOzAvvFw3HjJlaSU8rZ8nxtlGc0c6/gB093no6JYpeoUw23BF/yvAA0AxA8fnb9Fat4QklRDH\n6W9v76CmqRuAi+YVkJfpwhvw4ottoCztApPTRZeUmGSSrGl0J9QBRWSlxtLULh+6hAgHwx3T/zOw\nFfgh4GRgMh0hwpLe20pbZx/rtjXQ0e0BoMEzMJyvUiaaGS0q5TgK8cXLObhChJvhin6O1vo2rfVr\nwLXAGSHKJMQJmTYxHYvl0x5ljWc31p50YmxyV71Qy3EU4o9txkuf2VGEEEMMV/QPHIAbPNte3r0i\nYgQCAWo8u7F1Z5odJSql27LBb6PVqDI7ihBiiOGKvhyEExGrrrueLn87tq4ss6NEJYthxdqVQash\nV0wIEU6GO5GvTCk19GLbnCGPA1rrCUHMJcRJ2dRUToIlGcMz4tWlIkisXZm0uraSceBqXyGE2YYr\n+pNClkKIUba5sZxsewG1ZgeJYpauDPqNdfRamziGqT2EECEw3K11d4cwhxCjps/fS0XbLs5MuEiK\nvoks3ljiAml02KqxSh9CiLAgNxcXY86+nt3YLDbcNrmtq5l6+71YOzPotMnJfEKECyn6YszZ1b2d\nyamTsBrDHb0SwRQfYyPF5cToyKDH2ojfIhf/CBEO5K+iGGMC7O3eyRfGL6Ff5oYxzYxJbmZMcuPz\n+7hl5Yf0xdQC482OJUTUk56+GFOMhBZ6/T1MSTvyXfW2V7aFOFF0s1qslKZOoi9GbqIpRDiQoi/G\nFEtyPVnOXFyOw88Wz0yNo7Gtl3HuBJx2qwnpolNZWgl9MXUEAn6zowgR9WR4X4wplqR6CuJmHbHt\nhkumhDiNAJicpghY+2n11wOFZscRIqpJT1+MGV3+VozYLgri5AY74cTlSMDel0K9b4/ZUYSIelL0\nxZhR799NoDeWFHu62VHEIRy92VL0hQgDUvTFmNHg342/LQPDkNtGhBtnbzZtgQba+jrMjiJEVJOi\nL8YEv9FPi78Gf2smHT39vLWuCl3ZanYsMcjen4KDWLY0a7OjCBHVpOiLMaE3pgYrdgKdKbR09PGX\n5Zqd1e24k2PNjiYAAwO3dTybm8rNjiJEVJOz98WY0BtTRbolj70BCzVN3QD89KuziIuxm5xM7Jdh\nzae8+V18fh9Wi1wyKYQZpKcvIp7P76M3ppYMawE2q4X3N9dhs8qvdrhxW/Po9faxs01O6BPCLNLT\nFxGvonUXAcOL25LHH2868kx8wnwOI4bCpDw2N5VTnDLB7DhCRCXpDomIt7FxC45+N3YjxuwoYgRl\naSVyMp8QJpKiLyJaIBBgY+MWYnpzzI4ijsHkNEVVZw0tvXJlhRBmkKIvIlp1Vy2Nvc3E9kjRjwTj\nEnJIdLjY0iS9fSHMIMf0RURb37CJ3IRs8LnMjiKOgcWwkB9XxP/7eDXLlvlxp8Tyg8ummx1LiKgR\ntKKvlLIADwLTgD7gGq31jiHtS4DbAS/wmNb6EaWUFXgYmAQEgP/QWm8OVkYR+T5p2MQp6WV8Ih3H\niPDYq1tp8rjwxG0k1x1LTVOv2ZGEiCrBHN6/BHBorecCtwD37m9QStmB+4BzgYXAdUqpDGAJ4Nda\nnwn8BPhlEPOJCNfQ3URVZw3TM6aaHUUco5UbarD3ZmBY/CRldZodR4ioE8yiPw94HUBr/QEw9H6n\npUCF1rpNa+0BVgILtNYvAtcPPqcAaAliPhHhPmncRHpMKjnxWWZHEcfgkx1NACyZPYmS1InU+3eb\nG0iIKBTMY/qJQPuQxz6llEVr7R9saxvS1gEkAWitfUqpx4FLgS8cy4bcbjmeG8lOdP9t+WQrs/Nn\nkpGRiN1uJT7eIb8LIXasP++pE91UNXRyakkG43OSmJs0k2c+eYVYW7HsMxPJzz76BLPotwNDf6P2\nF3wYKPhD21wM6dVrra9WSv0n8IFSqlRr3TPchhoa5M5dkcrtdp3Q/mvra0c37eRzBRfQ0NCBx+Oj\nq6tffhdC6Hj23SXzCg56PCGmiG5fJxZrk+wzk5zoe0+Y72Q+rAVzeH8VcCGAUmo2sGFIWzlQrJRK\nUUo5gAXAe0qpK5VStw4+pwfwD34JcZAPqjYQZ42nryWRfQ1ybDjSpMQkk2LNxBNfbXYUIaJKMHv6\nLwDnKqVWDT7+mlLqciBBa/2wUuoHwHIGPng8qrWuUUo9BzyulHobsAPf01r3BTGjiFAfVn9Ce3Uq\n97y3nhnF6WbHEScgxz6BrfGbzI4hRFQJWtHXWgeAGw5ZvG1I+zJg2SHr9ABfClYmMTZ0e7qp7a8k\ntvd0Fs4aR1ObXPYViXLsRWx2vEddVz2Z8RlmxxEiKsiMfCLibGoqx27YMbqkhx/JXJYUAr3x/Omd\nf/GvtfvMjiNEVJCiLyLO+oZN5DonYMivb0TLTU8gxz6BhsAudla3j7yCEOKkyV9NEVH6fP1sadKM\nd040O4o4SSX5KXzl9AX4YlrwGN1mxxEiKkjRFxFlU+NWDMMgx1lwYNnWPS3srZNLjyJRfuJ47IFY\nWq17zY4iRFSQG+4IUz3y0iY2bK8H4NIFE5hSmDbs89fWb2Ba+mRs2AE4dZKbjORYAIpyk4IbVow6\ni2EhyZdHmxR9IUJCir4wVXVjJzarhb11nXT1eId9bq+3l81NW/la2VfoHvicgMpLQeWlhCCpCJZk\nXx4Vjn/S1ttJgiMOq0UGIIUIFnl3CdMVZifidFhHfN6mxq1YDSuTUyeFIJUIlQR/FgGflZv/8gIr\n1slkPUIEk/T0RcRYW7+Bqell2K12s6OIUXTx3CL6d5SiexvMjiLEmCdFX4Stzh4P3X1eHnl5C1XN\nrVC2lfMzL6Wv34fHK7MzjxUZKXEsLJjFlrYn8ARkAk4hgkmKvgg73b1emtt7eWNNJSs31ACQlt9E\nl9+CpcvNDfe9DUBygsPMmGIUlaQWY/it1Hh2AkVmxxFizJKiL8LO1j3NLH1hYE72SeOSuHZJGc/s\n/D+27MjCmjRw7P/GL04jOcFpZkwximwWG/buHPY5K4BzzY4jxJglRV+EpVinjbv/Yw5Wi0HA4mFb\n63ZoncnH3QPHfcdnuEhxSdEfSxydudS73qfH20OsLdbsOEKMSVL0RViyGJAQO3DC3vs167FbHUzJ\nLsHnNUhPisFukwtPxhp7TyZW7Pxj2xpKE6cCMDE3CcMwTE4mxNghRV+EvY/rPuEUdxlXLphqdhQR\nRAYWMm2FvLr1fV7c7gHgkf9chJR8IUaPdJdEWGvv76C8ZTunZc4wO4oIgWxrEZakRi5akGt2FCHG\nJCn6Iqx9XPcJLnsCk1LkjO5o4LaOB7+VBv8es6MIMSZJ0Rdh7aPadczKmo7FkF/VaGAxrPhaMqnx\nVZgdRYgxSf6SirBV193Ano5Kzsg61ewoIoR8zVnU+/aCtd/sKEKMOVL0Rdj6qHYdOfFZ5CZkmx1F\nhJC/LQ07TqypdWZHEWLMkaIvwlKAAB/VruX0rJlmRxEhZyHHNhFrmtx8R4jRJkVfhKe4Vpp6W5iV\nOd3sJCKEvL4AALm2SVgTW2jubTE5kRBjixR9EZYCyfsoTp5ASkyy2VFECL3wzk4Aki2Z+HvjWFO3\n3uREQowtMjmPCBsd3f3c+9f1tHX1wrhqTsu6xOxIIoSuv7jswN0TO7r78W3PZk3dOs7PXySz8gkx\nSqSnL8KGx+dn865mUnNbsdgCzMiYYnYkEULj3AkUZidSmJ2I027F15RDTVcdj6/40OxoQowZUvRF\n2Amk7eW0rFPkpitRLDnBycLSSdj7U9nSusnsOEKMGVL0RXix97KjfQdzsk8zO4kwUWZqHFedr5iU\nUEZv/F78Ab/ZkYQYE4J2TF8pZQEeBKYBfcA1WusdQ9qXALcDXuAxrfUjSik78BiQDziBu7TWLwcr\nowg/tvQqUpwpTEwuNDuKCAO5tmI2W99le8tOVOpEs+MIEfGC2dO/BHBorecCtwD37m8YLO73AecC\nC4HrlFIZwFeABq31AuAC4IEg5hNhJhAIYHXvY3rqDDlxSwDgNOJw9GTzXs0as6MIMSYEs+jPA14H\n0Fp/AMwa0lYKVGit27TWHmAlsAB4FvjpkGzeIOYTJvqkopEnXy9nd037gWVN/moMZw+npJ5iYjIR\nbmI6C1jfsIFuT4/ZUYSIeMEs+olA+5DHvsEh//1tbUPaOoAkrXWX1rpTKeVi4APAj4OYT5hoT20H\na3QDJfmp5GUmAFDp2YK/zU2iI9HkdCKcOHtycFqdfFwv1+wLcbKCeZ1+O+Aa8tiitd5/Nk7bIW0u\noAVAKTUeeB5YqrV+5lg25Ha7Rn6SCCtx8U7ysxP50ZUDA0DPvqup9e3A2zCV1NQE3OnxJicUxyLY\n772EBCd2u40zCmfzUcNaPj/9vKBuL9rI387oE8yivwpYAjyrlJoNbBjSVg4UK6VSgC4GhvbvUUpl\nAm8A39Rav3WsG2po6Bi91CIkurv68PQPHL1paOjA66rEhgN/awbNzZ3Y5GztsOd2u4L+3uvs7MPr\n8XNK8iks2/Ym63dtkxswjZJQ7D8RHCfzYS2Yw/svAL1KqVUMnMT3faXU5UqpaweP4/8AWA6sBh7V\nWtcAtwFJwE+VUm8NfsUEMaMIA4FAAH/qbsbbSyEgV5GKw+UkZFGQmMd7NR+ZHUWIiBa0nr7WOgDc\ncMjibUPalwHLDlnne8D3gpVJmK+prZeNu5rYXftpD6OidRc4O8i3lbGJBhPTiXA2J3sWf9/5OpcU\nXYjNIjOIC3EipFslQqqqsYsnX9fUNHWRkjgwiPNu1XsYnZnEWeQEPnF0p2aeQr/Pw8bGrWZHESJi\nycdlEXJ2m4VfXz8HgNaeNtY1bMTSfBpkmRxMhKX61m7++7kNTJmQSra1iBc2r6ByewILp+eSFO8w\nO54QEUV6+sJUb+5cRYozGaMzw+woIgzlZ7k4c2oONU1dvLRyF9vWJ9MUqOSlDzfR3tVvdjwhIo70\n9IVpfH4f/9yxkvm5c3jlY5mBTxyuND+F0vwUSvKSqWnuJhAYz0e+XVRnVJodTYiIJEVfmGZj01ba\n+zqYk3Mar7DO7DgijM2Y5GbG4Pexu+fwXPcr/PWtchLjYrluSZmp2YSIJDK8L0zzr73vcmb+6STY\nZSIecexmZc7AZjXojt3LmvJ6s+MIEVGk6AtT7GmvZEfbLj6nzjE7iogwibGxLBh/Ov7U3WZHESLi\nSNEXpnhz7zuUpk5ifFKO2VFEBJqfO5uG/lqMuFazowgRUaToi9Cz97CuYSPn5C0wO4mIUJnxGeTG\n5GPJ2Gt2FCEiihR9EXKWjN1kxWVQklJsdhQRwaYlzsJIqaa+o4V+j8/sOEJEBCn6IiS8Pj+bdjWx\no6YRI72Ss8fPxzDkMj1x4vJjiwj0xfGTl/7K7/4qt90V4lhI0Rch0dPn5b6/fsLrFavAb2VW1oyR\nVxJiGEW5yZxfuICY7H34Da/ZcYSICFL0RegYPpIm7OPSknOwyw1TxEmKi7FxoToTi2GhO3aX2XGE\niAhS9EXIWN378AV8zM+dbXYUMUY4rHbyrGW0xZZz/7PreHt9ldmRhAhr0t0SQffk6+U0dXRjy97F\n6emzibHFmB1JjCHzx81l55717GvaQVZjgtlxhAhr0tMXQbdpVzMt9p3YHT7m58w1O44YY04rGs/c\nnFmQsYMAAbPjCBHWpKcvRt32fa1s3tUMQEF2IgH89CZrzh0/n5yUZJPTibHo3LyzWFX1Ia2BKmCS\n2XGECFtS9MWoq6hqY/lHlditFk4r8eBN3Iff382i8fPNjibGKHdcGon9hVQ61hIInCWXgwpxFDK8\nL4IiOzWOkvwU/PjwpJVTHHsKCQ65sY4InrSeMtqpZ0PdNjxev9lxhAhL0tMXJ+Qfayop39MCwMLp\nuUwrSjvi8+otmgAeSmNnhTKeiEIOfxLe5iz++P4LXF18NbPLssyOJETYkZ6+OCF7azuobe5mW2Ur\njW09R3yOHw9VlvXYm4txWOSMfRFcX15czH+ccTFWVzN1/fvMjiNEWJKiL07YhOxEEuMdR21vcpQD\nYG+ZEKpIIoqlJ8cyI28CRkc2r+99k2/et4L6lm6zYwkRVqToi6DwW/podG4i1z8DIyBHkUTofKnk\nQqyuZvpj6gnIFXxCHESKvgiKzqTN2PyxZPjl8ikRWvNLFDPSZ2AfrwkE5IQ+IYaSLpgYNRVVbTz0\n0ia6aQW1k7yeRbyvG+j3ym1PRWidnXM2H9evZ+mKN8iyFPOtS6eaHUmIsCBFXwwrEAhQvrf1wGM1\nPhmL5cjXQHu9fpraeymctwubtYCLJ55OU34fABPHJYUkrxAAmQmp5FmmUJu4iebNqWbHESJsSNEX\nI7rn6XUHvv/jTQtxWqxHfa4luYF67z5um/l9suPlj60wR1yMje/M/zw/WflrSNttdhwhwkbQi75S\nygI8CEwD+oBrtNY7hrQvAW4HvMBjWutHhrSdAfxGa70o2DnF8L64qIhn39rBL55YQ2ZKLHHOw391\nvH4P9rxy5ufOITs+87i3UbGvbTSiCgFAvD2OGYnz+CDrHTr6O3E55GY8QoTiRL5LAIfWei5wC3Dv\n/gallB24DzgXWAhcp5TKGGz7EfAw4AxBRjGCtMQYvry4mJy0OOpajnxd/kctqzEsPj5XeN5xv35W\nSix1LT3kpMdjtcj5pWJ0TE6YTqA/hr/pV+nq9ZgdRwjThWJ4fx7wOoDW+gOl1NCp2UqBCq11G4BS\naiWwAHgOqAA+D/wlBBnFCJLiHZxemolhGFQ3HX7tc3VnLWtb36d/zzTi7LHH/fq3XHHqaMQU4iAW\nw0r/nhI+jFlDlU7jx184x+xIQpgqFF2qRKB9yGPf4JD//rahY7odQBKA1vp5Bob8RZgLBAI8rf/G\n+LhC/C3HP6wvRLCU5qdw52WfJcWfT5PrY/xyCZ+IcqHo6bcDriGPLVrr/e+8tkPaXEDL8W7A7XaN\n/CRxQgKDs5skJcfhdrtwJTixWi04Y+wA2GwWqoxyqjpquK7025Sv2HLc+0P2X+SKhH2XB8zcsZC3\nOp9iY8cGFhfJ3R73i4T9J0ZXKIr+KmAJ8KxSajawYUhbOVCslEoBuhgY2r/neDfQ0NAxGjnFEewv\n+m2t3TQ0dNDR2YfP56dv8PhoH+2saX2LS4svhJ6B0y+OZ3+43S7ZfxEqovZdbywpnVN5ct3z5DsL\nSXbKJaQRtf/EQU7mw1oohvdfAHqVUqsYOInv+0qpy5VS12qtPcAPgOXAauBRrXXNIevLRJphKoCf\n7sw1pNqyWThurtlxhBhWoLGQeEsSj2947sCHWSGiTdB7+lrrAHDDIYu3DWlfBiw7yrq7AakmYarO\ntgmftYNZcRdjMeSMexG+nHYrHV1eOtdOpGnKah55900+VzaH7LR4s6MJEVLyl1ocN6/XT7Ovnhr7\neuIaTyHOIscFRXj73NwCfv/d+Swum0xCRwnret6ioq7e7FhChJwUfXFU73xSzcurdh+2vL6jnZ2O\nt0j2FWDvHB/6YEKcoMsXF/PLi67E8MTxVPmzfP03b9Lc3mt2LCFCRoq+OKq31lWxcmMNxeOSiB2c\ngW/ulExmnFNFZlIC3597BQZHnodfiHBls9i4qvTfcSS3Ys3ca3YcIUJK5t4Xwzp75jguOCPvwONV\ndavZ2b6DH532XdLiZVpTEZnOKCqi2f8ZXva/Qm13HamJ+WZHEiIkpKcvjtnmJs1LO17j8pJ/O6G5\n9YUIJzPSZuJvdfM/nzzB/X9bY3YcIUJCir44JjVddTy26SnOGb+A07Nmmh1HiJOWEOvg/OwlWC1W\ndljfpruvn36Pz+xYQgSVDO9Huf/75zb8/oFrlr941kT+f3t3Hl1lnd9x/H3vzd3IHhISiCwBzI9t\nAIFatnFkGMcqVrTjVJ22OrY6pzgzjm2tp3849nQ650yL1c7pIrUutTPFWu04xwV16owICrKoyDby\nM8QQIAFCSMh+c5fn6R8JGZSEJIZwc3M/r3NywrP8nnzDN8/9PuvvFwycO2xua7SNdbv/A5M/jeun\n/c45yzfvrsUePs3sMg2lK6kjK+znxmXljNl9Ay/W/Rf3/s/TzBmzmHtumpvs0ESGjc7009xbu2qx\nR07z5gc1xBLn9ksed6M8uucpwhkhbpt1yznv418+s5jJxdnsPFCnUcwkJS0rL2f1pBvxl1bS4teD\nfTK6qegLS+eU9DrfJcHOyKu0x9q5e96fEMo4d5Tj1cvL+MZXyoc7RJFhkxX2c7W5nMnuImoz32Hd\n/21m7yenkh2WyLBQ0ZdeJZwETYXbaHUb+e78u8gNqgMeGd0WFy0lPzadffyCfTXVyQ5HZFio6Kep\nX753hGd/VUHCOfeSfsyJ8+T+9cSC9SwOrWZseGD36vdW6uxIUtcV80r5m6vvIBQv4r3YK5zqaEh2\nSCIXnB7kS1M7D9TR3B5j3rRCssOBnvnRRIzH9/6EmtZa8k+sIDun/4Lv9dLzEJ8e5pNU5vP6yG9c\nQvteOvUAAAzuSURBVOeEbTy46Z/g4BKIhVm7ZgmZ3cNJi6QyFf1R6NVt1dScbAXgK4smUjY+p9f1\nlswu5vplZRyt61q3Pd7BEx89Q31HA/cuWMO6g1UD+nn+DB9/cfP8CxO8SJJ53Qxqds4kUP4+4Rk7\naN6zAA3KJ6OFiv4o9NGhBpraYhw71caC8nGUje+/jSfYzqP7HsPn9TI7sYp//GkFp5ojMFOd8Eh6\nuf1qQySWIOpcxguHnyUyczsn2ueTFS5NdmgiQ6Z7+p9Rf7qDNY9sYs0jm7h/3dZkh/O5zZs+lqD/\n3Hfue3Ok7TDBWe+SF8jlvoV3E+sIEQpm8I2rypk1JX+YIxUZWUqLspg2IZeZlxRx69Q/xG3P4dF9\nj7HvxMFkhyYyZCr6n+ECndEEi0wRnaO8dy7XdXmj+i3WV/0niYYS7phxG+GMMADj8sJcOb+UScV6\nal/Sl9/rJ1oxn9Zjhazb+wQ7j+9KdkgiQ6Ki34cJhZl9LuuMJrCHG7GHG6msbbqIUV04jqeTXfHX\neLXqDa67ZDWx6tn4vAO7MiCSLkoKwvzwziVce8l1JGrLeXr/s6zdtJ6EM7pPCGT00j39z6G+qYO/\nf6briH9sToiH7l6a5IgGZ8/J/Rwf9zphN8z9l99Doi2T59iR7LBERhx/ho8JhZnMm16I41zJezXj\nOFK0jR+8/c9ckX8tmb4cZk0pICcz0P/GREYAFf0huOnKaWz8oCbZYQxYa7SN5yte5P0Tu8lsNywe\n+0XGZxZztK012aGJjGhTSnKYUpJDy+tRNu2B2LQ9/KzzaWLVs7gve5WKvqQMFf0h8HiSHcHAJJwE\nm2veZUPVG+QFc/jLRd/hmZfq8BUq/SKDccvKS7npyuk47lfZVLuZDd432HAsQlHR18gN5uI44M/Q\nXVMZufSpP4q5uBzq+JgNO7bSHG1lVdlVXFG6pPvefV2ywxNJOQG/j0B3Hz3XTF3JCxtaqZiynwfe\nWUu8Zhrj4rP54Z1LkhukyHmo6I8A2/Yf55fvHwXgsksLWbVkypC25+JS51ZC+XbeamhjeeliVpVd\nRVag74cTRWRwPMB3r12Oz7Mc27aHt70baYge5YGf1ZITncLKBZOYWJzFuLxwskMV6aGiPwI0tUU5\n1RwhOxygvily3nUdx8Wlq3swj8eD96x7DB3xDrYf/4CjBRtJOG3QOpGvl/0BV5ipwxq/SDryeDzM\nn14IwBf4MjOPzObVyjep8u/kVOd+/m3TVNzT47ls+jhuWjGN4vwxSY5YREV/xMjNDFBS0P+Hwr+8\nsJcPD9b3TBfkBJlUFqfeW0FToIqA18+YyBTm5Cxga+1pMufrPXuRi2H2xGJmT7yVps5V/KJ6I++E\nt+NzK9l9pJQVbeNU9GVEUNFPQQtnFOELN9Pgq6I2UcGBjA6clnwS1bNoaSwG18vCJZnAaTbuOsrG\nXUeZPaWAuOMyd+pYwkG9jy8yXHKDOfx++WpWlV3FpsPv8krsLdZ9/GNMneFLky5nVoFRnxiSNCr6\nKaI11saBhgqOZW6jI3CMmKeDydkTubF4JQvGzSU/lEdDc4TG1k4A8rOCtLRHOX6qncraZn59qBGA\nn2/+pGebnx0R71RThINHm6hv6iAc0J+GyFBk+sewcuKV/O/zLt7cevYW1nKg8Sd4XT/ZiVJm5c/k\n6pkLCHiDXeuHM/B59eS/DC99so9ACSfBifaTVDVV80lzNVVNhznRXkfIFyLDU8RkdxG3L1tOQejT\n/eIX5IQoyAn1TH/zmpk9/44nHGJxh7rGDto743iAwtzQp9r/urqB5zdWEgz4WFheNKy/o0g68Pu9\nrF2zjEPHWvjvX1XgPx2n0XsYCuvZ2vw6W7a/itOaj9OSz40LFjG3ZBqlY/OSHbaMYsNW9I0xXuBR\nYC7QCdxpra08a/nvAt8H4sBT1ton+msz2kQTURoijZyIH6I9p5ojoTaivtP8+eZm4k6cvGAuZTmT\nWDbhcspyJzM5+xL+9YX9FOWFzyn4/cnwecnweZlccv57/IW5IdauSa0eBkVGKq/HQ2FumMLcMItm\njOuZ/9KWKg4cOUWr7zhZhU1UeQ/x8rHnePk4jPUX4XRkE07kMym3lBUzZlCUnUtIV9/kAhjOv6Ib\ngIC1dqkx5reBh7vnYYzxA48Ai4B2YIsx5iVgORDsrU0qSTgJWmNtNEdbaYm20BJtpbn7e1O0mVMd\nDdRHGmiJdvWE58WHd0w2fgoY60zlujlfoDRrPPkhHfGLjEbXLyvjesqArjdyGloi7D1Ux/otOziR\n2YRnTAve8AlqI++zffeLuIkMQm42fiebzpYgU4tKiLb7KR9fTOGYPC4tGUdeZgh/hp4VkPMbzqK/\nDHgdwFq73Riz6KxlM4GD1tomAGPMO8AVwBLgtT7a9MpxHRJOAsd1cHBxXAfXdXBcFwene7pr/pl5\nPctdp3vaJeE6xJ049c1tePNOUBN3cPLqeLtmG3EnTsyJdX+P09jajr/sKO9HKohMbOah9z6gIx4h\nEo/QkYgQTUR74vN5fGQHsnq+cgM5zB47k7HhfArDBYwNFfDuh41sqz5BScEYwsEM5hTOuGBJEJGR\nzevtuhowf2oJId8XAcjO9FNamEVbNMKLO/aSyGgj4m3mZOcpPJlNfNReg9ffSVVjAhqBGnDjfnxO\nEI/jx+P48Tp+PG4ATyIDEn4mFubh9/rp7IRJRXnMmzCHjtZY12u/nq5xRM7cHvRwpsdRT0/Po12r\nnZk4s06KdEsqPYaz6OcAzWdNJ4wxXmut073s7OHpWoDcftr06pbnvn0BQwYcL4GpXvZ3+kkUuTy7\npwKP6wPHC64PXC+u4wWPj4AnTLTZ4UhrCBK54PjxJDLIcPx44gFIBCHhpx0P7cCJT/2g9u6vo0Ri\nCYrzuzrweHffcT6sqD8nrDPaIjFWXHbJhf2du7205RDZYf+wbFtEzi8/O8iSOSWfnkeQb1+zvM82\nxxqbaIw0s6uqBjI6iTjtRJ1Oom4nUaeTmNvJyZYWvBlxqqO1xNwYDnEOnnR4c0cdiZZ8XNclnnCH\nFHv38UHPQcHZBwpntp01DJ8tw3HMMSyHMRc40PU/uOZztx3Oot8MnH0D+ezi3fSZZdnA6X7a9Oq5\nm9fpUHOIioqyefnh1Un9+ZKalLvkOvP//6W5s5IciaSK4Xw/ZAtwLYAxZjGw56xlB4BLjTH5xpgA\nXZf2t/bTRkRERIbA47pDu6zTF2OMh988iQ9wB7AQyLLWPm6MuQ54kK4Djyettet6a2Ot/XhYAhQR\nEUkzw1b0RUREZGRR908iIiJpQkVfREQkTajoi4iIpAkVfRERkTQxajtzNsZ8GbjVWntXsmORgTHG\nLAW+1T35vTM9Nkrq0H6XmowxK4GbgTHAWmutXpdOIcaYhcB36Opb6H5rbV1f647KM31jzDRgPhDq\nb10ZUe6iq+g/SdcHkKQQ7XcpLWyt/RbwD8BXkx2MDFoQuBfYQFd39n0alUXfWltprX0k2XHIoPms\ntVHgGDA+2cHI4Gi/S13W2leMMZnAPcDTSQ5HBslauxWYBdwHfHi+dVPm8n73qHt/Z61d0dcQvMaY\nvwWmA2ustaeTGK58xkDyB7R399A4ATievGjlswaYPxmBBvjZWQisBR601vY9+IdcdAPM328B7wHX\nAH8NfK+v7aXEmb4x5n7gcbouYcBZw/YCf0XXELxYa79vrb1VBX9kGWj+gH8HHqPrMv9PL3ac0rtB\n5E9GmEHk7mGgGPiRMeZrFz1Q6dUg8pcFPAU8BKw/3zZT5Uz/IPB7/KYQLKfvYXt7WGv/6OKEJ/0Y\nUP6stR/Q1V2zjCyD2v+0340oA933bk9OeNKPgeZvI7BxIBtMiTN9a+0LQPysWdn0MgTvxY1KBkr5\nS23KX+pS7lLbcOQvVZM96CF4ZURR/lKb8pe6lLvUNuT8pWrR1xC8qU35S23KX+pS7lLbkPOXKvf0\nzzgzJODPgauMMVu6p3UfODUof6lN+Utdyl1qu2D509C6IiIiaSJVL++LiIjIIKnoi4iIpAkVfRER\nkTShoi8iIpImVPRFRETShIq+iIhImlDRFxERSRMq+iIiImlCRV9ERCRNqOiLiIikiVTre19EksgY\ncwOwFqgEbgOCwFvAa8CPrbWVyYtORPqjvvdFZFCMMX8MrLbWrjbGLAU81tot/bUTkeRT0ReRQTHG\nZAGHga8DYWvtK0kOSUQGSEVfRAbNGPM44LfWfjPZsYjIwOlBPhEZFGPMGKABWJbsWERkcFT0RWTA\njDFe4C7gAaDVGLMiySGJyCCo6IvIgBhjPMCfAk9ba2PAE8CdyY1KRAZDRV9E+mWMuRZ4HVhqrW3q\nnj0BuNEY82fJi0xEBkMP8omIiKQJnemLiIikCRV9ERGRNKGiLyIikiZU9EVERNKEir6IiEiaUNEX\nERFJEyr6IiIiaeL/AXOILqYEpPiHAAAAAElFTkSuQmCC\n",
"text": [
"<matplotlib.figure.Figure at 0xa26dcc0>"
]
}
],
"prompt_number": 5
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"See? Buried in all that nonsense is really good info. Let's get it."
]
},
{
"cell_type": "heading",
"level": 2,
"metadata": {},
"source": [
"Estimating sample statistics"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"So let's estimate our sample statistics both directly\n",
"and from the fit theoretical formulas in the wiki article.\n",
"\n",
"\\begin{equation*} m = e^{\\mu + \\frac{1}{2}\\sigma^2} \\end{equation*}\n",
"\n",
"\\begin{equation*} s = \\sqrt{\\left(e^{\\sigma^2} - 1 \\right) e^{2\\mu + \\sigma^2}} \\end{equation*}"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"wikiM = numpy.exp(wikiLoc + 0.5*wikiScale**2)\n",
"wikiM_est = X_n.mean()\n",
"\n",
"wikiS = numpy.sqrt((numpy.exp(wikiScale**2) - 1) * numpy.exp(2*wikiLoc + wikiScale**2))\n",
"wikiS_est = X_n.std()\n",
"\n",
"summary = \"\"\"\n",
"\\tTheory\\tSample\n",
"wikiM\\t{0:.3f}\\t{1:.3f}\n",
"wikiS\\t{2:.3f}\\t{3:.3f}\n",
"\"\"\"\n",
"print(summary.format(wikiM, wikiM_est, wikiS, wikiS_est))"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n",
"\tTheory\tSample\n",
"wikiM\t13.913\t13.804\n",
"wikiS\t14.922\t14.748\n",
"\n"
]
}
],
"prompt_number": 6
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Everything is pretty close. Cool."
]
},
{
"cell_type": "heading",
"level": 2,
"metadata": {},
"source": [
"Now let's compute `wikiLoc` and `wikiScale` from our sample estimates"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"sampleLoc = numpy.log((wikiM_est**2) / numpy.sqrt(wikiS_est**2 + wikiM_est**2))\n",
"sampleScale = numpy.sqrt(numpy.log(1 + (wikiS_est/wikiM_est)**2))\n",
"\n",
"summary = \"\"\"\n",
" \\tTheory\\tSample\n",
"wikiLoc \\t{0:.3f}\\t{1:.3f}\n",
"wikiScale\\t{2:.3f}\\t{3:.3f}\n",
"\"\"\"\n",
"\n",
"print(summary.format(wikiLoc, sampleLoc, wikiScale, sampleScale))"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n",
" \tTheory\tSample\n",
"wikiLoc \t2.250\t2.244\n",
"wikiScale\t0.875\t0.873\n",
"\n"
]
}
],
"prompt_number": 7
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Again. Pretty close so we're cool."
]
},
{
"cell_type": "heading",
"level": 2,
"metadata": {},
"source": [
"The crux - making sense of the results of scipy's `.fit` method"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Recall what scipy returned:"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"print(\"\"\"\n",
"spShape = {:.3f}\n",
"spLoc = {:.3f}\n",
"spScale = {:.3f}\n",
"\"\"\".format(spShape, spLoc, spScale))"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n",
"spShape = 0.866\n",
"spLoc = -0.044\n",
"spScale = 9.505\n",
"\n"
]
}
],
"prompt_number": 8
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"But we want the estimates of `wikiLoc` (known = 2.25) and `wikiScale` (0.875) from the scipy fit."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"sp_wikiLoc = numpy.log(spScale)\n",
"sp_wikiScale = spShape\n",
"\n",
"summary = \"\"\"\n",
" \\tTheory\\tSample\\tScipy\n",
"wikiLoc \\t{0:.3f}\\t{1:.3f}\\t{2:.3f}\n",
"wikiScale\\t{3:.3f}\\t{4:.3f}\\t{5:.3f}\n",
"\"\"\"\n",
"print(summary.format(wikiLoc, sampleLoc, sp_wikiLoc, wikiScale, sampleScale, sp_wikiScale))"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n",
" \tTheory\tSample\tScipy\n",
"wikiLoc \t2.250\t2.244\t2.252\n",
"wikiScale\t0.875\t0.873\t0.866\n",
"\n"
]
}
],
"prompt_number": 9
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Phew. It all makes sense."
]
},
{
"cell_type": "heading",
"level": 2,
"metadata": {},
"source": [
"Now work backwards"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"So now we can use our original `wikiLoc` and `wikiScale` \n",
"values to freeze a `stats.lognorm` object."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"frozen_dist = stats.lognorm(wikiScale, loc=0, scale=numpy.exp(wikiLoc))"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 10
},
{
"cell_type": "heading",
"level": 3,
"metadata": {},
"source": [
"And then plot it"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"But the theoretical and fit distributions are so close we can't really see the difference."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"fig, ax = pyplot.subplots(figsize=(8, 5))\n",
"_ = ax.hist(X_n, bins=bins, cumulative=False, histtype='step', \n",
" normed=True, lw=1.25, alpha=0.5,\n",
" label='Numpy random number')\n",
"\n",
"# theoretical PDF\n",
"x_hat = numpy.logspace(-1, 2, num=1000)\n",
"y_hat_1 = ln_dist.pdf(x_hat)\n",
"y_hat_2 = frozen_dist.pdf(x_hat)\n",
"\n",
"ax.plot(x_hat, y_hat_1, '-', lw=1.25, label='Fit with scipy')\n",
"ax.plot(x_hat, y_hat_2, '--', lw=1.25, label='Theoretical')\n",
"ax.set_xscale('log')\n",
"ax.set_xlabel(r'$X$')\n",
"ax.set_ylabel('PDF')\n",
"ax.legend()"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 11,
"text": [
"<matplotlib.legend.Legend at 0xa71b4e0>"
]
},
{
"metadata": {},
"output_type": "display_data",
"png": "iVBORw0KGgoAAAANSUhEUgAAAf0AAAFPCAYAAABHzhf2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl83FX97/HXzGQm62Tf22xt09M23Xfa0lKgZVEQ9Iri\nD1SQRUV/oPDzB/iriOKCCuhl8SKLiHLFK4JAkU122tJC9y2nSdosbZZmnSyT2ef+kbSkW9Ilk+9M\n5vN8PPJoZs73O99355vkM+e7nGMKBoMIIYQQYvQzGx1ACCGEECNDir4QQggRJaToCyGEEFFCir4Q\nQggRJaToCyGEEFFCir4QQggRJWJC9cJKKTPwCDAdcAPXaa2rBrRfAqwCfMCTWuvH+9d5HJgIBIDr\ntdY6VBmFEEKIaBLKnv5lgE1rvQi4HbjvUINSygrcD6wAlgE3KKWygZVAotZ6CfAT4GchzCeEEEJE\nlVAW/cXAawBa6/XA3AFtk4FKrbVDa+0FPgSWAr1AilLKBKQAnhDmE0IIIaJKKIt+MtA54LG///D9\noTbHgLYu+or8h0AcUA48CjwYwnxCCCFEVAnZOX36Cr59wGOz1jrQ/73jqDY70AH8N7BGa/1DpdRY\n4G2l1FSt9Ql7/MFgMGgymYY5uhAjw9Ht5pU1+5g5MYs9tR1ctmy80ZGEEOHvtIteKIv+GuAS4O9K\nqYXAtgFt5UCpUioN6KHv0P5vgNl8enSgHbAClsE2YjKZaG7uGuboYqRkZdmjev91Oj309Ljp6HDS\n0+2KqPci2vddpJP9F7mysuxDL3QCoSz6LwArlFJr+h9fo5S6EkjSWj+mlPo+8Dp9pxie0FrXK6V+\nDfxRKfUBfQX/Dq11bwgzCiGEEFEjZEVfax0EvnXU03sGtK8GVh+1TgdweagyCSGEENFMBucRQggh\nooQUfSGEECJKSNEXQgghooQUfSGEECJKSNEXQgghooQUfSGEECJKSNEXQogo0dBQz8qVy/jud2/k\n6quv5rvfvZGnnnqcioo9PPXU4wC89947tLS0nPRrtrW1ct999wKwZcsmqqoqAbj00guGJfNf/vIU\nu3fvHJbXEqEdnEcIIcQguj09uP3uYXu9WEssSbbEQZcpKRnHgw8+esyIfKWlEwF47rlnKSkpATJP\napvp6Rnceut/A/DKKy9x/vkXMH78BIZrdPSrrvr68LyQAKToCyGEIbx+L6vW/hxPwDtsr2kzW/nV\n2T/GarGe0nqbNn3Ciy8+z4UXXkxFxR7uuefHPPLI48TE9JWIb3zjau6770GSkpK4+OLzePjhP1Ba\nqrj22qv48Y9/xj333MWtt97O+vXrqKjYQ3FxCR6Pl7vv/h+amhpJSUnhpz+99/DrATz//N957bVX\nMJvNTJo0hVtuuY26ulruvfcefD4fsbFx3H33z3n44d9y/vkX0NbWykcfraGjw4HD0cG1195ASck4\nfvKTVTz22J8A+NGP7uDKK69i8uSy4XpLRx0p+kIIYQCrxcpPF9057D39oQp+dfVevvvdG7FaLXi9\nfu66657DbWedtYTS0on813/deUSBPvvsZaxfv5asrGzy88fw8cfriYmxUlBQiM1mw2QyodQkFi5c\nxPnnX0BOTi69vU5uvPE75Obm8t3v3khFhT6iGL/66svceusdTJo0mX/+8zn8fj8PP/xbvvrVa5k/\nfyEffvg+FRXlDJxQLRAI8rvfPUJraws33ngNf/vbP4mNjaW6eh/p6ek0NNRLwR+CFH0hhDBIki2R\nJAY/HD/ciouPPbxfW1sz6DpLly7nT396gtzcPG644ds899yzBAIBli8/74TrJCenkJubC/SdAnC5\nXEe033HHXTz77F+orz/A1KnTCQaD1NXVMnXqNACWLFkKwJtvvn54nTlz5gGQkZFJUpIdh6ODSy+9\nnH/962VycnK58MKLT/HdiD5yIZ8QQojDzGYzgUDgiOfGjRtPff0Byst3cdZZi3E6nXz44fucddZi\ngsHg4eVMJhN+v7//+8G38/LL/+S22+7goYf+wJ49mh07tlFUVMKuXX0X7b355mv84x//74h1yst3\nAX0XD7pcvaSlpXPOOeexYcNHvP/+u6xcKUV/KNLTF0KIKGI6TjU2mUyHn586dTr33HMXDzzwMHb7\np1O4zp49l8bGekwmE7NmzaG6eh+xsXFHvOaUKVN59NGHyc8fw9FTvh+93fHjx3PTTdeRkJBIVlY2\nZWXTuOmmbH71q5/zpz89QXx8PP/zPz9B692H19m/v46bb/42Tmc3t912ByaTCZvNxsyZs3E4Oo7I\nK47PNPBTWoQKypzQkSva5/TudHp4e+N+ykrS2Vffycr5hUZHOmnRvu8iXaTtv1dfXU1HRwdXXnnV\nMW0PPPArli07l9mz5xqQbORlZdlP+94IObwvhBAiIhzvlMH3v/8durq6oqbgnyk5vC+EECLsXXTR\nZ4/7/P33PzTCSSKb9PSFEEKIKCFFXwghhIgSUvSFEEKIKCHn9IUYRTqdHjq6+kZ4S4q3kp4cZ3Ai\nIUQ4kaIvRIRodbg40NINQGpSLIU5x96T3Nzey859bZjMJsZmJUrRF0d46KHfovVu2tpa8Xo95OTk\nsW/fXubMmcfdd/88pNveu7eSrq4uZsyYxV133cmqVT85YqjfoVx66QW89NLrQy8oBiVFX4gI0en0\nUNvUjc1qwesLHLfoA9gTbdjjT23CFWEMv7OHgNN5zPPmhAQsCccOz3uqyx/tO9+5Bei7572lpYGr\nr76ezZs38s9//uM00p+ad955i4yMTGbMmHVaHzCGa9a+aCdFX4gIEmezkJESRyAQ8YNqCaD9zTdo\ne/nFY55Pv+RzZH7u8jNefjCHBmYLBoPs31/Hbbf9J+3t7SxefDbXXnsDVVWV/O53vyEYDJKSksId\nd/yIxMQkHnzwAbZv3wrAihUX8sUvfpmf/ezHdHY66Ozs5Fe/+i3PPPMntm3bQiAQ4Etf+grTps3g\n1VdXY7PZUGoSq1bdzl//+jyNjQ3HzKrX1tbCQw/9Fr8/gMPRwW233c7UqdNP6f8mTkyKvhBCGCRt\nxUpSFi855nlzQsKwLH+yPB43v/zl/fj9Pr7whc9y7bU3cO+99/DDH/6YoqJiVq9+kWeeeZqysmk0\nNtbzhz88hc/n49vfvo45c+ZiMpmYM2c+V1xxJevWraGhoZ5HHnkct9vNN795DQ8++AcuvvgSMjIy\nmTy5DJPJRDAYPO6sep2dXXznO7cwbtwE3nzzNV555WUp+sNIir4QYehAczcub9/EJXnpiSTEya/q\naGRJSDypw/Knu/zJGjduPDExMcTExGCxWACoqdnHb37zCwB8Ph8FBYXU1FQzY8YsAGJiYigrm8a+\nffsAKCwsAvrO3Wtdzne/eyMAfr+fhob64273eLPqbd26haeeeoLY2Ficzh4SE5OG/f8bzeQviRBh\nqKq+E6fLh8vjw55gk6IvQuzYE+aFhcWsWvUTsrNz2LJlEw6HA6vVyr/+9RJXXPEVfD4fO3Zs5aKL\nPsP69Z9OqFNUVMLs2XP4wQ9+iM/n489//iNjxow97ux9h2bVmzt3Pm+++RqdnZ288spL3HXXPRQV\nFfPEE4/S2NgwIu9AtJC/JEKEqXH5yeyuaTc6hhilDhXpgTPs9bcAcNttd/DTn/4Iv9+PyWTijjt+\nxNixBWzevJFvfvNavF4v5523gokTJx3xekuWLGXz5o3cdNP19PY6Wbp0OQkJCSg1iYcf/t8UFRUD\nfdu86aabj5hVb9Wqn+LzeVm16r/Jzs5h0qQptLa2HJFLnJmQzbKnlDIDjwDTATdwnda6akD7JcAq\nwAc8qbV+XCn1deBr/YvEAzOAHK115yCbkln2IlikzfQ13E40y977W+vJTU9gd007Z03NJTs1nn0N\nnVQdcBy+kG+Oyj7m9aoOOKg92I093orFYmJWaVbIskf7vot0sv8i15nMshfKnv5lgE1rvUgptQC4\nr/85lFJW4H5gLuAE1iilXtJaPwU81b/MQ8DjQxR8IUaVtpb96LWv0bG3irfnZeL2Bti/N4sZOZPI\nDI4zOp4QIsKFsugvBl4D0FqvV0oNnPdwMlCptXYAKKU+BJYCz/U/nguUaa2/E8J8QoQNV3sjprf/\nRNNjDSTGWvBlpTEpcQqNvl4IOPmHfhU/PkpjZ3FWcDEmLEZHFkJEoFAW/WRgYC/dr5Qya60D/W2O\nAW1dQMqAx3cCPw5hNiHChmXvh/DRuyTarQS++nnKFl3Mmh0HmZ2ewO7udqbkppHnbCYutxG9+RW2\n79lE8fyvAcce3hdCiMGEsuh3AgOHDDtU8KGv4A9sswPtAEqpVGCi1vq9k91QVtbxRyYTkSGa99/W\nPeuoCm4idm4pcfOvIn9CHrv2dxMwmUhNTSCxuYfU1ATsSQl8ccEF1Hebqf394+xofYi2slWorPFH\nvF6r00u700tyUiwxZlPI39to3nejgey/6BPKor8GuAT4u1JqIbBtQFs5UKqUSgN66Du0/+v+tqXA\nW6eyIbkYJXJF88VEG5u28NTOZxlfvIySCcvZV99JXX0HB9t7yctMxOz30+P00NHhpKfHTXNzNwmz\nF+O6wsSUfzzBe7/6CZ03/BfFqUWHX7Oj3UlXtxv8ASwWU0jf22jed8MhEAzicvuG/XXjYmMwn8SY\ntbL/IteZfFgLZdF/AVihlFrT//gapdSVQJLW+jGl1PeB1+mb3vcJrfWhmzEnAlXHvpwQo4duq+Sp\nXc9yYeFKvA3FAASC4PMHSU60UVacfsJ1g8WldFx6LerlJ1n/1H0k3nAXWQkZI5RcDBeX28cbH9cN\n++uunFdAQpzMvSCOL2RFX2sdBL511NN7BrSvBlYfZ73fhCqTEOGgpbeVJ3b8hWVjF3F2/hLebtgP\ngMvjo7api/zMwUdci7VaaM0uxLnsCqa/8yzPrPkD31z+feJiYkcivhhmZ0/PJz72zC/M7HX7+WDb\n8Ue+O2TTpk+4887bePrpvx3uLf7+9w9SXFzCRRd99owzGOGrX/0STz/9N6NjHGPTpk947723+d73\nfmB0lCPI4DxCjCBn4wH++dGTFBSN4fLxn6HH1T/UbkYieRl9xd5iHvzQ7JTidKYUp/OO2YR/2h20\nO57nr/ofXFP2lZDnF8MvPtYyoj1zq9XGz39+N88882eAowbmEcMlXN9XKfpCjJCg30/FQ7+hxOJi\n/sXfw2K2AH1FP8ZiIs526r+O8ak5fKPgKn698SGmZ04hlZJhTi1GE5PJxOzZc4EgzzzzDCtXXnq4\nrbGxgbvuupNHH/0jADfeeA133/1zXnnlJerr99PR4aCzs4PPf/4K3n33LerqavnhD+8mPT2dn/3s\nx8THx9Pa2sKiRX2z9H3lK1/gsceexm6388ILz9Hb6+QrX/nq4e1dffUVFBYWERNj5TvfuYXf/OYX\neDweWltbuP76b3H22efwta99mVmz5lBZWYHJZOKXv7yP+PgEfv3rX1BVVUF2dg49PT0ANDTU84tf\n/OTwUL+33PJfTJhQype+dBnTps2grq6WOXPm0dPTza5dOyksLGLVqp8c8f4cb3tal/Pii88fng74\nc5+7gBdffJ2f/ezHxMRYaWpqwOPxcP75K1mz5gOamhr5xS/uA2DPHs0tt3ybnp5uLr/8i1x88SXH\nnb1Q63J+//sHsdlsXHrp5VxwwcUh+xkwh+yVhRBHqHzx/2Jq7SD7a9dgtyWxs7qNijrH0CsOoTB5\nLBcVn8ez+gW6fXJhljixQyOw3nrr7Tz11FMcOLB/yHVMJhOxsXHcd9//Ztmyc1m3bg333vsAV131\ndd5663VMJhNNTY389Kf38thjT/Pxx+uprNzDihUX8u9/vw7AG2+8ykUXXXLE67pcLr7+9eu5++6f\nU1NTzZe/fBUPPPAwP/jBD3n++b8D4HQ6Of/8C3nooT+QlZXNRx+t5YMP3sXtdvGHPzzFbbfdQU9P\nNwAPP/xbrrjiKzz00B+4+ebb+OUvfwr0fZi54YZv8/DDj/Hcc3/j85+/gsce+xPbtm09vO4hx9ve\niXrsJpOJ/Px87r//IYqLS2hoaODXv/4dy5ady5o1HwBgNpt54IGHeeihP/DnP/+Rjo4O7r33Hm69\n9XYefPBRFi5czDPPPI3JZMLr9fLww4+FtOCD9PSFGBG99Qfwvf42B86fykXj5gOwr74Te4KV3IyE\nIQ/pD+WConPZ1ryT9w6+xbSY84cjshjFkpNTuPPOO7nnnruYNm3GcZcZOET7ofH1k5LslJSMO/y9\nx+MBYMqUqcTFxR3+vq6uls985nP8+Md3MmPGLNLT00lLSztmG4dm5ktPz+Dpp59k9eoXMZlM+P3+\nAdtWAGRn5+DxeKivP8CkSVMASE1Npaio7+hWTU01M2fOBqC0dCIHDzYBkJKSSnZ2DgDx8XH9Y/9D\nUlIiHo+HxKMuoTl6e8e+LwOX/fR9OfS6dnsyHo8bgOnTZx7+0FRcXEJjYz21tdXHzF448L0INenp\nCxFiwWCQiicepD47lsWX3nBEmypMY+GUXKwxJ76Qq7vXO+Q2LGYLV6jL6a7cTLNn+K8IF6HT6/bj\ndHnP+KvX7R96YwMsX76cwsIiXn2173pqm81Ge3sbgUCArq6uE0yHG+R487VUVVXg8/nw+/3s3r2T\nceMmkJubS1JSEk8//SSf/exlx81gNveVoCee+D9ceOFnWLXqJ8yaNeeI2fiO7mkXF5ewY0ffHeCd\nnZ3U1dUCfTP2bdmyCYCKCk1GRkb/+qfwphxnezZb7OFJfxobG+jsPPmjc7t37yQYDOJ0OqmpqWbM\nmAIKCvpOKzz44KPceONNLF689LjbDRXp6QsRYm1NNXhbmom79oskx9rZua+Nzh4P/sDQk12ZTVB3\nsBvzCY4E+PwBet0+XB4/9bvcXPZuJx/1vEhg8X9ikV/viDDUFffD6egZ9W6++VY2bvwY6Ottz5u3\ngOuu+ypjxoxl7NiCI9br/27A7Hx9jw89/4MffI/OTgfnn7/y8NGASy65nN/97jfcddc9x0tz+Lvl\ny8/n4Yd/y9///ixlZVPp6jrxlCtnn30OmzZt5Prrv0ZmZhbp6X3F/TvfuYV7772HZ5/9Cz6fj9tv\n/9Ex2xl6pr5j2ydNmozdbueGG75OcXEJ+fljPl16kEP/h/695ZabcDq7uf76b2G324+YvdBsNnP7\n7atobj44YkU/ZLPsjSCZZS+CRcMAIX/Z/Xf2O/bzgwU3YzaZ+XBbA4FgkIzkOApykkhOsJ3W676z\naT+Onr7Dj2Ul6ezc18aEfWvwbniHXVd9gSkZ82WWvTA2WgbnaWio54EHfs2vfvXAMW3vvPNv9u6t\n4hvfuPGMtyM+Fa6z7AkR9Rp7mvio4RNumvkNzKZPz6ZlpcYzuejYc5ynYsGUHHz+IG9v+vRirElf\n/CJ7Nq2Fj/+Nb+WsM3p9EVpmk2lUDKLTd/Tg2OcfffRhtmzZyL33/nbkQ4kTkqIvRAi9vPd1ipOK\naahOoLHmABnJwzeATkKclcBRpwjMsbFkXXI505/7K9tmf8Q89flh254Qx5Obm8e99x7by7/xxpsM\nSCOGIhfyCREiNZ11bGnewbKcc+lx+bBazPS4hv9w7tEyl51LID4R88fv0uvrDfn2hBCRQ4q+ECEQ\n8Hp4be+bTMucTH7CGGLMJtLsIzNMrikmBvcXbmDbtAzerVs7ItsUQkQGKfpChEDty89R+v/WcUHR\nuYZs35yRy5S0Bby7/0M8/qFv+RNCRAc5py/EMAt4PPS88w6OOWMoSSniQHP30CudoZYO1zHPqcRp\nbO/+iDUHNuBt6hsAxGSCFXMLjllWCBEdpKcvxDCrf+8N/H4f6oIvjMj20pPj8PoDpCfHHXGbsdVs\nY9nYRbxd9z7dvW7S7LH0jsA1BUKI8CVFX4hhFAwEaHvzNWqmZDE5d2rIt2c2m1g6I//w19H3Zy8b\nu5hubzdNwSpSk05vPAAhxOghRV+IYbT19XewdnTTPeVsmh3HHnIfaUm2RFa2ZeM/sM7oKEKIMCDn\n9IUYRlW9tfROTcUeOwWv99TGQh9uwWDfML2lbVYSKvZzcH4jMDJ3EAghwpP09IUYJoFggI8zGvAu\nWY7NEkNnj4eD7c6TGmM/FGqbuli9tprAnHPJa/WxV79FIBhk9dpq1mxvMCSTEMJYUvSFOEXvb63n\n9Q21vL6h9ogZ8Ha37aHb38kUe990mhX7Hazd0YjbgB7/5OI0ls0cQ0JsDKasXNry8kjdvJ3pKpm8\njETDPogIIYwlRV+IU+Ry+8hIjqPX7SMwYMKq9/avpSiulMSYJC5cUMj5c8caljExzkqaPfbw7Hyd\nk5cyodbF3uYN2BMif7x3IcTpkaIvxGlIOepK+JbeNna1aibGzzQo0eCcY0oJJCbQtO69486HLoSI\nDlL0hRgGH9V9RG5iNlnWfKOjHJ/ZTOJNN7OuJMiB3lqj0wghDCJFX4gz5KzeS9Fvn2NJyjRMx5lj\ntLN/znuAjm43XU7PMcuMhLTcYqZkTmJn51ZDti+EMJ4UfSHOUO0bL9KcFsOccYuOajFhs1poau/F\nGmMmPTmWMZlJjMlKIn2EJt852qK8eVR0l+MJuA3ZvhDCWHKfvhBnIOhxE9iyg67lk7DbkoBPx9lP\niIvh4oVFRyyfl5E4wgmPNDVzMjazlVqPBkoMzSKEGHlS9IUYQuUBB7tr2gEoybUf0daxeT2BQIAJ\nSy40ItopizHHMMk+jX1dO4HIyCyEGD5yeF+IIQQCQeJtFlISbUfcogfQtvZtakrs5NtLqTvYjctj\n7Ch8J6MsZgLZFdWUH6yhqc1JU5uTXrdMxCNENAhZT18pZQYeAaYDbuA6rXXVgPZLgFWAD3hSa/14\n//N3AJcAVuAhrfWfQpVRiJNljTETZ7Mc+WQgQHugm9gFc2lq7WXP/g5irRaslvD+LJ3hsXHBR12s\nKXiL9MTlBINB5qhsCrKTjI4mhAixUP51ugywaa0XAbcD9x1qUEpZgfuBFcAy4AalVLZS6hzgrP51\nzgHGhTCfEGeky9TOPxbHMW1e32HydHscF8wvRBWmGZxsCDn59KQlE1e+jeJcO/GxcpZPiGgRyqK/\nGHgNQGu9Hpg7oG0yUKm1dmitvcCHwFJgJbBdKfVP4GXgpRDmE+KMNAYqyIvPI9Ec3kX+eIfu3ZPm\nMGFfFw0uuWdfiGgSyo/4yUDngMd+pZRZax3ob3MMaOsCUoBMoAj4DH29/JeASUNtKCvLPtQiIoyF\n+/472OWh2+MnIc5KfGwMXe4AySlxNAUrGR+cxYY9zUwpzsAdCL//S1JSHE0ON4mJsWRkJJGWHMfB\nLg/OOUtIWPcOlXXrKSz6PGlpCaeVPdz+v+LUyP6LPqEs+p3AwJ+oQwUf+gr+wDY70AG0AuVaax+w\nRynlUkplaq1bBttQc3PXMMYWIykryx72+6+jw0lXlwufx4fHZqGn20XFwYO4gj0sLZzD3lrP4WXC\n7f8yb2IGh6499Lo8NLu9dHQ46TTF48vPIXb7ThzpF9Le7iTJemoH/iJh34kTk/0Xuc7kw1ooD++v\nAS4GUEotBLYNaCsHSpVSaUopG32H9tfSd5j/wv518oFE+j4ICBFWdnXuYGLaeDISUo2OMiiL2UyM\npe/r6NECfUsuYee4OA76awxKJ4QYaaHs6b8ArFBKrel/fI1S6kogSWv9mFLq+8Dr9H3weEJr3QC8\nopRaqpTa0P/8t7XWMjuICCvWxhrGbH+f9K99zegoZ8Q0tpSEzskc8OzhyEtuhBCjVciKfn+x/tZR\nT+8Z0L4aWH2c9f47VJmEGA7mirUk9PiYmTsdR2cQgsGIvc99QuIU/u18CbffDcgte0KMduF9Q7EQ\n4SYQILW6kg5VRHxMPAD+QJCapsg8N1oQPw4zZqq6KoyOIoQYAXKDrhCnwF1Vjt3lJXnGYgBy0uK5\ndEnkjmEfY4ohJ6aE9fVbCLTlcfaMfJLirUbHEkKEiPT0hTgF3Vs/YH+OjeKc6QCYTCbMA74i0aKx\ns+h2VeP0uAgG5RIaIUYzKfpCnKSZEzLJaKyjd8o4inPC+6r9obg9fpz91yHM7DLx1RebaPPt42BH\nLw2tPQanE0KEihR9IU6g1+2jprGLju7+uefNAf7ymXTGLF2BxRy5vzrWmL7sPl+AWKuZZDUFS8CE\npX4Lu6vb2bSn2eCEQohQkXP6QpxAd6+XzRXNJMRZSUm0sbtN44kxMSV3qtHRzkhJXjIleclHPBec\nNI7s6hqmX5DK9soOg5IJIUItcrsrQowAk8nEynkFLJiSw+aD2ynLUMTFxBoda9jlLDyH4v0u9rZp\no6MIIUJIir4QJ8Hr97K9ZRezsqYZHSUk0mbNxRI00bx9zdALCyEilhR9IU5CeXsFvqCfqZmTjY4S\nEua4OHzzp9PUvR9/0G90HCFEiEjRF2IIwWCQ6rX/Zqp9AnExcUbHCZkJX/sme4piaQnUGR1FCBEi\nUvSFGIKzugr10iZmJo43OkpIxcfEMc4+nib/XqOjCCFCRIq+EEOoW/c2B9OtTJ0w3+goITcxWXHQ\nX0MgGBh6YSFExJGiL8QQ3Fu345iYf3is/dGsNHkiHnp5W+/gk/KDfFJ+kP3N3UbHEkIMEyn6QgzC\n2tFCfGsXqbNHfy8fIMmaRIoph82NO9jf3E1rpwtHj8foWEKIYSJFXwjA5w/g9fV9BQKfjj9vq91I\nu93C5MmLDEw3cnLSE7jYYSe3egefXVRMcoLN6EhCiGEkI/IJAWytbKHuYN9h7BkTMg+PWLc/PUDM\n/AIWxKcZGW/EmE0mCvyJmHYcpNXVYnQcIcQwk56+EP1yMxJIiDtyWtndOZ1kLFxiUCJj5Mw/m6wO\nH7sq1hsdRQgxzKToC9EvNsaCxfzp9Ljt7na6g21MzyozMNXIix1bgCclkY5NG4yOIoQYZlL0hTiB\nckc5cdjJT8w1OsqIMplM2GZMI7WiEVfASW1jF+9s2s/BdqfR0YQQZ0iKvhAnoDvKyTYXYzKZhl54\nlBmzYBn5zV48gSomFqbS7fLh9QeHXlEIEdak6AtxHN3eHmq6a8g2lxgdxRAJEyay7eolVJsOMD4/\n5YjTHkIcJau5AAAgAElEQVSIyCVFX4ijBIOw+6lHWbS9l1RTntFxDGGyWCidspDd7RV4/HKfvhCj\nhRR9IY7i8/ux7NBYbBnEmKP3V2Ry+kSCwQB72quMjiKEGCbR+xdNiBPwNzeQ1OOmYMFiLlkcnYf3\nAWItNkrTxrOjtRyApjYnlQcceH0yLr8QkUqKvhBHcexaQ0eShZKSeUZHMVxZxiR2tOzGHm+lvcvN\njr2tUvSFiGBS9IU4il/vonFMBnExcUZHMVxZ+iSCzS2MH29mYVmO0XGEEGdIir4QA4xNjyWpzYFF\nlRFnsxgdx3CpXV6+troNvUdG5xNiNAjZ2PtKKTPwCDAdcAPXaa2rBrRfAqwCfMCTWuvH+5/fBDj6\nF9urtf5GqDIKcbSkLB8PXZbBj866mPhYmZrCmpOLOzURx9aNMOVio+MIIc5QKP+qXQbYtNaLlFIL\ngPv6n0MpZQXuB+YCTmCNUupFoAtAa708hLmEOKGdLeVk23PISsoyOkpYMJlM2MqmYK/aitPXa3Qc\nIcQZCuXh/cXAawBa6/X0FfhDJgOVWmuH1toLfAgsA2YACUqp15VSb/V/WBBixOxoLacsc5LRMcJK\n3rwl5Dd5qGjcaXQUIcQZCmVPPxnoHPDYr5Qya60D/W2OAW1dQApQDvxaa/2EUqoUeFUpNbF/nRPK\nyrIPc3QxksJh/yU3dOEJuth3oIb/mPW5sMgULjKWLODAwxY6dn9M4pgvkJmZRGJ832yE8j5FNtl/\n0SeURb8TGPgTZR5QvB1HtdmBdmAPUAmgta5QSrUCecCBwTbU3Nw1XJnFCMvKsofF/uvs7GVfbzk2\ns41McsIiUzhxL5jBPsd+8lNctLR044yLCZt9J06P7L/IdSYf1kJ5eH8NcDGAUmohsG1AWzlQqpRK\nU0rZgKXAOuAa+s79o5TKp++IQEMIMwpxmKviY6YlFGMxy1X7Ryu58hp25UNn8KDRUYQQZyCURf8F\nwKWUWkNfIf+eUupKpdT1/efxvw+8DqwFntBaNwBPAMlKqfeBZ4Frhjq0L8RwCPR2M/eN7ZR50o2O\nEpZSYu2MScynOVjLvz+p4/UNtUZHEkKchpAd3tdaB4FvHfX0ngHtq4HVR63jA64OVSYhTqS36mPi\nLCZKZ51jdJSwNTVzMtsO7mJCXgo1TXJYWIhIJIPzCAFQtZ2mMcmkJqQanSRsTcucTENvA2ab2+go\nQojTJEVfRL1gIEBqbQO9xdE7uc7JKEoeS5I1kaquSqOjCCFOkxR9EfVaK3YS5/ITXzrf6ChhzWwy\ns6LKSte2dUZHEUKcJhlnVESltk4Xa3Y0AtDTVY5vRgpT0iYYnCr85XVb8NbuI1DgNzqKEOI0SE9f\nRC2/P8DYzET2JrTRPHcWJpPJ6EhhL2f2QsY2uGj3DTp0hhAiTEnRF1EtPt5Mk7eWsXHjjI4SEdKm\nz8HmC+Jt2Db0wkKIsCNFX0S1+t46/PgYE1dsdJSIYElKwpmbQVJd1dALCyHCjhR9EdX29VSSYckj\n1hxndJSIETOxjPz6Tpp7Wo2OIoQ4RVL0RVTb111JnlVu1TsVKWet5M2luWxukFn3hIg0UvRF1DI3\n7Wb5a1XkWgqNjhJRzCmpxKVOYHPDDqOjCCFOkdyyJ6JKIBhkzfYGfP4gtrqt2AIWkmMyjY4VcbIt\nhexoehvvRB9Ws/wZESJSSE9fRJ1Wh4vURBsZDfvpKSmQW/VOQ6alAF/AR2XHXqOjCCFOgRR9EZWy\nY5yktDmJnzTL6CgRyWqyMSlrAjtby42OIoQ4BVL0RVRq3Po+PXFmcornGh0lInm8AeJ6c9m9d6vR\nUYQQp0CKvohKrsrddJRkYY2JNTpKxElNimX6+Ayml7dy7ps1NDvl1j0hIoUUfRF1gsEg/5qfgPnS\nC4yOEpGS4q2U5CWTOW0eWR0+dldvNDqSEOIkSdEXUaebNjr8nZSNnWl0lIhmHluEO87K/vXrqG/p\nMTqOEOIkSNEXUaclUEtWXBaZ8ekA+AJBvL6AwakiT2ZaAr1FxaTWNrCrrpmtlS109niMjiWEGIQU\nfRF1WoK1lKaUHn7c0+ulsc1pYKLIVJSbzIQlyyhodNMVOEB1Yxe9bp/RsYQQg5BRNURUcXqdOIKN\nTEzpO59flGNnbFYSABaz3K9/qlKmzWRPlh2oxWweY3QcIcQQpKcvooretY7kHihKKgIgxmIm1moh\n1mohxiK/DqfKkpSE8/ovssVdQzAYNDqOEGII8ldORBXf6tdYvCuIxWQxOsqoUZYxiRZXGz3BDqOj\nCCGGcMKir5T61kgGESJUHN1uGlp7OFDfhr22he4xE42ONKpkx2eSGZ9BS6DW6ChCiCEM1tO/4dA3\nSqn3RiCLECGxr6GLj3cf5JN33wTAnDfb4ESji8lkoixjEs2BGqOjCCGGcLKH95NDmkKIEMvPTCSx\nfhfNeUlYbalGxxl1yjIm0Raox+N3Gx1FCDEIOacvokMwSFpdPe5x44xOMiqVkM7ZG7vY56g0OooQ\nYhCD3bKXpJRaCpiO+j4IoLV+fwTyCTEsnJ5ODhTEUDhlAU6X0WlGn9jYeGbs6WHHnk9gwjyj4wgh\nTmCwon8AuPs43x+yfLAXVkqZgUeA6YAbuE5rXTWg/RJgFeADntRaPz6gLRvYCJyntd5zcv8VIU7s\ngL+WDXOymZ07meZqucp8uFmSkujMTIc9muBFQUwmGfNAiHB0wqKvtT7nDF/7MsCmtV6klFoA3Nf/\nHEopK3A/MBdwAmuUUi9prQ/2tz0KyGDeYtjsd+0j11qE2SRntELFOVaRt+9jXvh4G5+dNRWbVW6L\nFCLcDDoin1JqMnAdMIm+4rwLeEJrfTL35iwGXgPQWq9XSg2cuHwyUKm1dvRv50NgKfAc8Gvg98Ad\np/ZfEeL4AkE/B9zVzE441+goo1rhwsX4Nq9jW8sO9tYXEhdroThXrgEWIpycsOgrpS4C/gz8DXiF\nvnP504FPlFJXaK3fHeK1k4HOAY/9Simz1jrQ3+YY0NYFpCilvg40a63fUErdQd81BEPKyrKfzGIi\nTIVq//X0enF0u2kONOILehmfOpHU1AQSE3tIz0gkKy0hJNuNJgP3XebKBbz/dBzJTXtoK1kBTi/z\npsnQvOFM/nZGn8F6+j8FLtBaHzFZtlLqj/Qdqj97iNfuBAb+RB0q+NBX8Ae22YEO4D+BoFLqfGAm\n8Cel1Oe01k2Dbai5uWuIKCJcZWXZQ7b/qhs72VLRQkVgB9m2fDxOEx0dTnp63LS19mDy+UOy3Whx\nvH1nvvqLbD34GvNTTOw74JLfzTAWyt89EVpn8mFtsBOcsUcXfACt9QYg8SReew1wMYBSaiGwbUBb\nOVCqlEpTStnoO7S/Vmu9TGt9jtZ6ObAF+OpQBV+IwSTGxTDv3Y9YapEe50gYP3sZ3sRYqnv2GR1F\nCHEcg/X0B5sj82QOu78ArFBKrel/fI1S6kogSWv9mFLq+8Dr9H3weEJr3XBSiYU4Bb7mSorquogf\nW0Zli9FpRj+rOYZJaaVUdVUyjjyj4wghjnKy9+lD//35/Y+ThnphrXUQOHr8/j0D2lcDqwdZf9Bb\nAoU4GaaaTbSn2igdO4nKlnqj40SFsgzFy1VvUBK3hD11HRxo6bsRZ3x+MoU5cg5ZCCOd7H36R9sf\ngixCDDt7bTWuiYVy3/gIKsuYxF/183QGWuh1J0AwiNcXwO2VayiEMNpgRf8/gIeAUvrOz9+utW4f\nkVRCnKJd1W1093oBUAWppCTF4unuIKOlG+/lc4dYWwyntLhUsqyZNLuqmEghifFWet2DnS0UQoyU\nwS7k+yOwG/gvIJa+wXSECEstDhcuj5+GVufhHmXr7nV4Y0xMmLHU4HTR53Or95OstxgdQwhxlMGK\nfr7W+k6t9avA9cCCEcokxGnJSYs/4grTnWPgX5dMJT5W7scfaabCErL3N+Py9xodRQgxwGBF33Po\nG621l77x84WICMFgkL09laSkTDQ6SlSyl81j7EEP9V0y654Q4WSwoi9XPomI1eQ8iMPrICem2Ogo\nUSlm/CRMQRPOfZuMjiKEGGCwC/nKlFIDR9jIH/A4qLWWiclF2NrRWk6aLZ0kS6rRUaKSyRaLIy+H\nhH3VBMuCQ68ghBgRgxV9OS4qItbOlnLGJU0YfIgpEVLukhlQ/x6t3iYSyTQ6jhCCwafWrR7BHEIM\nm96uVlprK5hRNh+vY+jlRWi4y5awefwBCtz7mGiToi9EOJDJxcWo07Th31z+Vhtj4guMjhLVfP4A\nOTFF7HfLOPxChAsp+mLU8ZXvpHNcDtYYm9FRopY1xky8zUK2pYgWbwOugNPoSEIIpOiL0cbrJqWu\nlaTpM41OEtXyMhJZPnssn587h4SYeJq8tUZHEkIgRV+MMsHGLZiDQSbMO++47a2drhFOFN0sZguT\n0yfS4JVD/EKEAyn6YlSJrdtJS56dlORjLxxLirfS6/KRkmjDYpFhKEbK9J4UUst3EQgGjI4iRNQb\n7JY9ISJOo93F2NJpx22bNyl7hNMIgIIeK5YtHTQuPYAi3eg4QkQ16emLUaPD08bHykTugpVGRxED\npM+cR6IrQHO1jM4nhNGk6ItRY29PJfEkkxWXZXQUMYA1PZ2u9GSo2GV0FCGinhR9MWrsc1aQZS7C\nZJLz9eHGVVRKVm0rDneX0VGEiGpS9MWo4Am4qO+tI9tcjMcbYF9DJ60OuVI/bJTMIbfVx+7azUYn\nESKqSdEXo0KDtxqr2UaaKY9ej4+tlS20dblJjJNrVcOBL7eAreeXsatbbt0TwkjyF1GMCsnr3+Gs\n9ETMRRa6nV4AzpmZjzXGYnAyAYDZTPz0RWxoeR1/wI/FLPtFCCNIT19EPJ/XQ8nuA2RbszGbTexv\n7sZilvP64aYosQSXz81eR43RUYSIWtLTFxGvassH2LwB0srO5pyyYqPjiBOItyRQklLIztZyStPG\nGR1HiKgkPX0R8Vo2rqMlK4k4u0zfGu7KMiaxq00bHUOIqCVFX0S0QCBAnK6hq1h6jpFgSoaiqaOe\nNmeb0VGEiEpS9EVEO1CzE3uXF8uE+UZHESchz5zKDc+3ULHtA6OjCBGVpOiLiLaNRv515WTiMqSn\nHwmsiUm4MpJp/Gg9//6kjnU7Go2OJERUCdmFfEopM/AIMB1wA9dprasGtF8CrAJ8wJNa68eVUhbg\nMWAiEAS+qbXeGaqMIvJtbd7BjJJZIAO9RYTNe5rpKhpPtt4GCRacvT6jIwkRVULZ078MsGmtFwG3\nA/cdalBKWYH7gRXAMuAGpVQ2cAkQ0FovAf4H+FkI84kI1+xs5UB3AzOzjz+rngg/NU1dWMYvIMPh\nw90lF/QJMdJCWfQXA68BaK3XA3MHtE0GKrXWDq21F/gQWKq1/idwY/8yxUB7CPOJCLe1ZQeZcenk\nJ+YaHUWchMY2JwBq9jScybE4d20wOJEQ0SeU9+knA50DHvuVUmatdaC/zTGgrQtIAdBa+5VSTwGX\nA//rZDaUlWUflsDCGKe7/3Zt3c3CotlkZydjt7eTmpogPwsj7GTf75KCNLp6PKSnJpCfm0LTvGl0\ntu6lJNEm+8xA8t5Hn1AW/U5g4E/UoYIPfQV/YJudAb16rfXXlVL/DaxXSk3WWvcOtqHmZjmhG6my\nsuyntf86ulro2KOZWLiS5uYuurpcdMSY5GdhBJ3KvivMSICMBAB8bi9jL7+K/7P251gcdTQ3Z4cy\npjiB0/3dE8Y7kw9roTy8vwa4GEAptRDYNqCtHChVSqUppWzAUmCdUupqpdQd/cv0AoH+LyGOsGvN\na3zhrQ6SvGl09niMjiNOUVpcKrlx+dT7qoZeWAgxbELZ038BWKGUWtP/+Bql1JVAktb6MaXU94HX\n6fvg8YTWukEp9RzwlFLqPcAK3Ky1docwo4hQvVu34M7LoEK3k5chRT8SldoVW1q3GB1DiKgSsqKv\ntQ4C3zrq6T0D2lcDq49apxf4UqgyidGhx9VFZm0r9fPOZvyYFJwuue0rEk2wKz5ofoemnoPkJMoh\nfiFGggzOIyKO3vg2sZ4gpuIFRkcRZyDdlkGiKZVXdq1nX0Pn0CsIIc6YFH0RcRyfbKA9Pw3i5crj\nSJacYGNOZya2jWtp65SzeEKMBCn6IqK4/R72JPQQmL/Q6CjiDGWmxjM3OY8ZWxtxujuMjiNEVJCi\nLyLKjpbdlE9MIm/OysPPtTh6ccgV/BGpaMG5xPiDdFSuGXphIcQZC+XV+0IMaVP5QfbW9U2zOrko\njey0hMGXP7iN6ZlTsJqtAORlJJAQ1/djnG6PC21YMexiEhJpG5OJrUKm2BBiJEhPXxiq0+nBbDLR\n5fTi8Q0+JIPL52Jn625mZ884/FxmSjzj81MYn59Cmj021HFFCPgmTiO3toVOZyeBYNDoOEKMalL0\nheHS7LHEWIb+UdzRshuLycKU9IkjkEqMlMTJS4j1BHnlpReobpAR4oQIJSn6ImJsOriNaZllWC1W\no6OIYTRJFbL18/PYl9VtdBQhRj0p+iJsebx+elxe3t9azz/XaCY8u4ZJTfH4/AECATkMPFokxVsp\nm3cuDaY63H6X0XGEGNXkQj4Rdrw+P063n70HHNQ09R3u9Ti2MvagB3/GJFavrQYgziY/vqPFpPRS\nLMRQ1b2HycjofEKEivT0Rdhp7nDxzqb91DR1kZESx8p5hSTv30F7WiIxGTkAnFWWy1llOQYnFcMl\nxhxDvnU8umu30VGEGNWkqyTCkjXGzPlzCzCbTHgCLtL0ATonzMDd6gQgOdFGfKz8+I4mY2ylrO95\nhV5fL/Ex8UbHEWJUkp6+CEsmk4lYqwVrjJld298lrctP+pzzSIiNYWxWEhazyeiIYpjlxBSQ4LXw\n0a61tDpctDpcBOUWPiGGlXSVRNhr2LkRS14qcxaWGR1FhJDZZOHSdb20xb3GB4tKAPjckhKDUwkx\nukhPX4S1Tk8Xbxb0kPDdG42OIkbC5OkU1LYwPl8GWhIiFKToi7C2sWkrdmsSE7OU0VHECEifeS42\nb5CD5e8bHUWIUUmKvghrHzduZm7uTMwm+VGNBpYEO0156bBjs9FRhBiV5C+pCFtNzmZquupYkDvH\n6ChiBDlKppJd3YzPJyP0CTHcpOiLsPVx42byE3MZk5RndBQxgoIFC9lXmEibe4/RUYQYdaToi7AU\nDAZpXvMOi03FRkcRI80WR+MFZ1MTW2d0EiFGHSn6Iiw5vPXMWduIciYZHUWMoENzKkyyT6Uj2ECb\nq93gREKMLlL0RVgK1KzFGjSRveBso6OIEbS7pq/I58blE08ynzRtMTiREKOLDM4jwobb62ftjgZ6\nPV6yKqvwThmPJV6GY40Wc1U2/v6evsfrJ89cysdNm1lZtByTSUZgFGI4SE9fhI1AIMjB9l56PHso\nrO+l4JwLjY4kRlByoo00eyxp9lgsFhO55lIaepp4X8skPEIMFyn6IuwEy9fiT4gldeoso6MIg8TZ\nYijLL+TCT3x0vvui0XGEGDWk6Iuw4gp282F2F7avXYnJLD+e0Sop3srMCZkkp48lR+/FH/AbHUmI\nUSFk5/SVUmbgEWA64Aau01pXDWi/BFgF+IAntdaPK6WswJNAERAL3KO1fjlUGUX4qQ9o4lMzKZ21\nzOgoIgykzV1J3EfbqNz6AWrWOUbHESLihbIrdRlg01ovAm4H7jvU0F/c7wdWAMuAG5RS2cB/AM1a\n66XAhcBDIcwnwkwwGKQ+oJmXPUcu3BIAJGQWcDDXTtMHbxkdRYhRIZRFfzHwGoDWej0wd0DbZKBS\na+3QWnuBD4GlwN+BHw3I5gthPmGgxjYnWypb6OhyH35uf28tvXQxO3OmgclEuOmeOIPU8v30ODuN\njiJExAtl0U8GBv6W+vsP+R9qcwxo6wJStNY9WutupZSdvg8APwxhPmEgR7ebhpYeslLjSU60AbDT\nsZVMUwEpsSkGpxPhJHbiMnxWM9t2vWd0FCEiXijv0+8E7AMem7XWgf7vHUe12YF2AKVUAfA88LDW\n+tmT2VBWln3ohURYaex0k+8PsnhGPgAVTe049mxjXP5nyMxMwp5gMzihOBmh/t1r6fbS4Uqn6dYv\nsrutikvld31Yyd/O6BPKor8GuAT4u1JqIbBtQFs5UKqUSgN66Du0/2ulVA7wBvBtrfU7J7uh5uau\n4UstRoSjw0lnlwvo23/VTR/zv15vofY8Dy0t3bjirQYnFEPJyrKH/Hevo8NJd5eLGaUzWV35Nlv2\n7ZEJmIbJSOw/ERpn8mEtlIf3XwBcSqk19F3E9z2l1JVKqev7z+N/H3gdWAs8obVuAO4EUoAfKaXe\n6f+KC2FGEQaCwSBx+iM8ifE488YbHUeEofykXIqTC1nX8LHRUYSIaCHr6Wutg8C3jnp6z4D21cDq\no9a5Gbg5VJmE8ZwuHwc7nLR3f3oBX0X7Xkor22D2UpB788UJnJU3l5f2vsZl4y8mxiwjiAtxOuQv\nrBhRXU4PWypa6HZ6ibf1/eHevv5fJPf4iZsj9+aLE5uTMwOP38v2FhmWV4jTJUVfjDiL2cT5cwuY\nOymbjl4Htk920jG2AFNKqtHRRBjqcXn5aFcjjc1e5jtyqP/r/0XXtuPyyB29QpwqKfrCUG/tXUNz\ncRq9M1YYHUWEoZQkG4U5drqdXspr24l1j2XC9ga279iG2xsY+gWEEEeQoi8M4w/4+XfVhxQsvQB3\nXrHRcUQYykqNZ9q4DMpK0hmfn8LUBefQmRqLreJdo6MJEZGk6AvDbG/dTae7i7Py5xkdRYS5vIxE\nJhakoorSMS2Yw9jKWrZV1LNRHzQ6mhARRYq+MMzbtR+wpGg+SdZEo6OICKJWfIFYT4Cu8rc40NJj\ndBwhIooUfWGIms46qhz7+Kw6z+goIsLYUzPonFZMbNV2o6MIEXGk6AtDvKf/TZl9PAUp+UZHERGo\n5MpreGGBFUdADu8LcSqk6IsR1xvsIuHtDZz3TpPRUUSEys0opDh5HLV+6e0LcSqk6IsRd6B3E1P2\nuchbKrfpidM3P3shDYFKWp0d+Pxy+54QJ0OKvhgRgUCQg+1Omjo7Sd2zGXNcPPa5ctW+OH0Tk0uJ\nJ5knN7zKuh2NRscRIiJI0RcjwusPsHZHI+/tXcOM8m4yVlyIKUbGTxenLyM5nvOLltLAbnxBr9Fx\nhIgIUvTFiPEHfSRUf0SC30z6crlqX5wZa4yF5cULmLKnm64d/zI6jhARQYq+GDEHAuW0JZvJ/NKX\nsSTIvfnizNksVkoC2eR/spG12+qobuw0OpIQYU2Kvgi5LZUtfKwbqQ5spnT6uWQuk16+GD4lF36Z\n+F4/7VvfoMsph/mFGIwUfRFyB9t7qXLuJGjysnTMIqPjiFFmXHExbTOKydz2McGAXMUvxGCk6Ith\n1+pwsbumnd017TS2OQkEA+xxb+TcwiWkJ9qNjidGoYmX/gfJHS7adr9ndBQhwpoUfTHs2rpcVB1w\nUN3QSVObk/0+TY+vm+UFZxsdTYxSOWMncGBCLj36I4LBoNFxhAhbUvRFSCQlWMlMjccf9NPQ8iEL\n4qaRZJOL90ToBFZ8hVdnW9jZUoFfDvMLcVxS9MVpqap3sH5XE+t3NdHU5jzhcru7t7JoQxOzN8jg\nKSK0kuOyyTGN59kd/6Kh5cQ/k0JEMyn64rQ4uj1093pp7XThdPuOu4w34KG+4h0KGt3Yll00wglF\ntJk+LoMvTb2I9mA9dT01RscRIixJ0RenLc0eS6zVcsL23T2bmLutna5xUzDnymx6IrQS4qxMzCog\nzzKef1W/yeq1++julVv4hBhIir4ICXegl7bKD8hvdNE551yj44goctnEC2kPNtDkkd6+EEeToi9C\nYqdzHYu2dOGbNhdveo7RcUQUmZpXxMruMUx69Vn8/uOfehIiWknRF8OmrdPFGxtq2Vizl72eHXSs\n+Cy1ZUvpPcE5fyFCZc7sz5LicPPR80+zYXeT0XGECBsyzZkYVDAYpMXhOvw4IyUOs8l03GUDgSA9\nLi91sRsYZxvH/Ann0uvxA5CeHDcieYUAyM4tYsfsieRs2MD+orNhshxtEgKkpy9OwprtDYe/AoHB\nBz5pCdZQ66zmyimXkZWWQGGOncIcO0nx1hFKK0TfDHyL/uNbBMwQs+VFo+MIETZC3tNXSpmBR4Dp\ngBu4TmtdNaD9EmAV4AOe1Fo/PqBtAfBLrfXyUOcUgysrSWfnvjbe21JPUrwVa8yxnxe9AS/av5az\nx55FXuKp96zaOt3DEVUIAOyJqTjOWUjh62voaKwlNbfQ6EhCGG4kevqXATat9SLgduC+Qw1KKStw\nP7ACWAbcoJTK7m/7AfAYEDsCGcUQEmJjmDY+A3uC9YS3Qb3X8B5+fHy2ZOUpv35ifAzdvV7sCbYT\nnj4Q4lRNWPJFdk1M4/W97+Dx+o2OI4ThRqLoLwZeA9BarwfmDmibDFRqrR1aay//v707D6+qvvM4\n/r77kvVmARKJBAF/oBBQQBG1uGut1K3FMh21rdVn9LHTdp4Zp390eeaZ9en2dKYdO3V7tNplRqut\ng1tbpWMBFREFRHJYgyJJyHqTm+Su58wfN0iERBJDcnPJ5/U8POGe3zkn3+Sbe77n/M65vx+sAz7R\n37YbuAFQBZgAAj4Ps6pLqCgJDdp+MNbE3rfXcv0GCLn9I97/hXXVXLp4Opcunk44qEdN5MTwenxE\nl1zHy/HtvLBtW67DEcm58Sj6xUDXgNeZ/i7/w23RAW3dQAmAZVlPku3ylwnOcRx+veMJrticIOgu\nwuUZesAekfFUURLis0vPozY0h02xtdiOxuSXyW08Lqm6gIHzqbotyzr8zose1VYEdIz0G1RWarrW\nseI4DgUFAcrKC6ksC9PRl6a1O0FxcfaKP+3AzsQ2Im/sorgnw3srPzXifCh/+Stfcveprmv42Y4f\ns617K5fN0myPh+VL/uTEGY+ivx5YCTxujFkGbB3QVg/MMcZEgB6yXfvfG+k3aGnpPhFxyiAcx6Gn\nJ4iprkwAABEXSURBVEF7Wwx3JkNHRy/dsQTe/psujdFmNjX9L5/f0kN45WfodvwjykdlZZHyl6fy\nKnfxAAvCy/n5m08yIzCT0kBJriPKubzKn3zIaE7WxqN7/ykgboxZT/Yhvq8bY1YbY27vv4//N8AL\nwAbgQcuyGo/aXpNjT1C2Y7Op9wWuer2XgtNmE1imKyiZuGq9dVSnCtl4//ewM3qoTyanMb/StyzL\nAe48avHOAe1rgDVDbNsALB+z4GRUtnVvpCfVTmltHVNXXkuXW8M+yMTk9bhJpR2mJ5dQtf1RXnns\nPuo+extF4ZE/dCqSz3SUlhGzHYemvoNs7lrPgsKLCV99E/4pU3IdlsiQTq8p5ZPnzuCMuQvZtayO\nsvWv0bhzR67DEhl3KvoypP1N3ex8r/OY5Z29MV5se5qZYUONz+QgMpGPZ8Fp5az8y7t5v6qItl/+\njN+u3am5IWRSUdGXIe1r6mJ/c4zy4uAHI/DVTCmgpfhVikNBbl90Ey4NpCN5xufxUXPzXYR64oQ2\n/0+uwxEZVyr68pFOqyrmwoXVlBRmB0b8v4Pr2NW5mzvqbqEoGM5xdCIfj5k9j9iNV9Hka+Jg7Ohn\nh0VOXir6Mmzb2yz+uG0Nt79sUx7Thyokvy1cdh0NC+Zz37ZHefnthlyHIzIuVPRlWBp7mnl4y2N8\n7jWbQm8Yb1l5rkMSGZWAz8MNM6/H43bzp7ZnSabTpDMasU9Obir6k9y2vW1s3dPK1j2tQx7wYske\nfvrWQ1zzVpqinjTVd92N23dkqtz9Td28Xn+InrgeiJL84fd5WFA7lVW1N9GaeZ+frHuKTdahXIcl\nMqZU9Ce5hsYuWqNx9h7sImMf22WftBPcu/Uh5lsxqna2UH3XV/CWRj5oP6WigJJCP++3xEimNeCJ\n5J95VaeyataN7LU38V7X27kOR2RMqegLNVMKB12ecdI8/f4TeFo7WPRqI9Nu/SKhWbM/tM7cGREW\nnKaufslfAZ+HC2aczbLQMmY/9kvW/+bXNHf05joskTGhOUxlUBk7w6ux5+ly2rhnxd0ULugjcMop\nuQ5LZMxcUHspWxbu4PTnn+dAYSlTr7wq1yGJnHC60p+k9hyM8vbeNgbp0Sdlp3lw+y9oTR/khlNW\nUx4qG1bBb27vG4NIRcZHbVUx19z291gLqgk9+d80bduU65BETjgV/UnqYEsPTe29TC0LEfB5Plie\nzKS4b+sjNET3c3HxjZQHKo+7L5cLpkRCuN3ZryL5yuP2UHLhHew/vZLWe+/l2ade4JlXGkim9LyK\nnBzUvX8S2vleJ929SQBOqy4hUhQYdL2aKYWYUyNEe7Lr9qb6eKT+Mdp62vja4jt5Z2diWN/P43az\nfH7ViQleJMe8bj+9S75EQ+o+mttepKBcf9ty8lDRPwm1RvtIJDN096WoKi8YsugP1Ot08R9bfkNV\nYy9Xru9ge+p9ugPFUDEOAYtMIAtnV5DOlJGo+w6P7nqMd2K/5fzeSmpKVPwl/6l7/yi98RRrNjSw\nZkMDv3/9vVyH87FNLQvjcQ9vXPyG7gY2pp/kjIYEFz2/H9ssxBUpo25WOZWlwTGOVGRiKS7wU1Yc\npCpSwpfmfYEiVwX//tZ/Ud+2N9ehiYyaiv5RHCCdsamuKCBjn9yjczmOwx/2/4kHtj/IRdtcLHxp\nD5WrPkfqkpUUhALUTiv+YMx9kcnI5/ZR57mCMqeW/9xyP683vZnrkERGRUV/CEVh35Bt6YxNa7SP\n1mgf7V3xcYzqxEnYcdY0PsGze3/PlzeHMfWNVN5xF5HLrsh1aCITRmHIy2VLTuWGmdcyx3MOz657\nlA3//A0SnR25Dk3kY9E9/Y+hN55m3dbszFzhoI8rltbkOKKR2dqynec6HyfsC3PPOV/F7d3NG/N8\nnLbo7FyHJjKheNxuisN+MmUFXM5F1CdDOJ2/of7b9+BbdSvFpy+isjREwO85/s5EJgAV/VE4o7aM\nhqbuXIcxbLFkD4/v+h1vNG/BBM/mkuqLqSqYQvTsCMnNB3IdnsiEFSkKECkKkEgtof7qCsIbf8Xp\njzzAjnlzCdxyJ5UVJbkOUWRYVPRHwTW85+RyLmNnePn9V3hm3x8oDRTzd0vu5t19HrzuoW9hiMix\n5s8s48zaCPYFC3nlxV9Q/szLNPzgm7i/9Q9EQhEcx8Hj1l1TmbhU9E9ijuNQH93BA7vXkoxG+fyu\nEDMXL6O0uIZ3OZjr8ETyjtdzuKB7WHHlrfzKO4uW1pf4+avfZ5Z7MfOLlnDZ4hk5jVHko6joTwAH\nDsXY29gFwLSyMKfXlI5qf47jsDtWz4bEy6R2t/PpAxGq32zFV15BoOL4I+yJyPBccvZi4Czean+T\n3x/4Awe6dtC48Txqg3OZVV1KSaGfgqB61GTiUNGfAOKpDL3xNAGfm97jzElvO0cGy3cBrgH3GPrS\nfbzWtJkX2v9MX6qTi3eEmWt14QtnKF+1muLzL8Dl0QNHIieCy+ViWlkYgKvKLmRB+Xye37eWjZ0v\nsi32GivWhbBnX8LUmmmcUVtGYUjFX3JPRX+CCPo9wzoobNzRTFPbkWk/QwEvTqidnT3b2Rt/B7/b\nT23gTJZWL8X7558TvPJaaq66DLdPBxyRsXRKJMJtkRuIJi7jpa3PULnrjxRu2Upj7XTar76RgrqF\nHzpJF8kFFf08VF1RQLfTwt5eC6v7Hfq6u4i4q5nrXsEUZuJOeSjwFbHvqpvpLQxwoL6VKZEQtu0w\ntSw84L6kiJxoJYFirl+6mljdSja/9Diu9a8R//GP2FxRSPCTlzP3wmvwuNXjJrmhop8nYqke6tt3\nsb7zLZpS+/F3RDm7LcwtTRnKzlzK9OtW0RtPE09mbw8E/V4SqQyxvhQdXQlaOrPT3u7Yf2RQkaNn\nxOtLpGnvitMXT+P16sRAZDQKA4Usv+JWni64kIPR7RTtfJWGvWt50H6D6cHTmF8+j3NrFhD0ZEe9\n9PncuNUTIGNMRX8CytgZmntb2Bfdz96u/eyLvktz7yGmxjycW29z+aFegp3d+Co8FNTVUbQgO6hO\nOOglHDyS0rPmHHloz7YdMrZDTzxFKm3j6l9/oJbOPrbva8frcVNVXjAuP6vIyczjdnHVOTPoiFXx\n9tSzmONK4O/bTVtqP2sOPM3TB56i1DWNut1pZk2r4dT5S6monanbADJmxqzoG2PcwL1AHZAAvmxZ\n1p4B7SuBbwFp4CHLsh443jYnm2QmSXu8g72x99gbf5/eeAfRZAuP72yhMwSlgRJmFp/K+dXnMLNk\nBlXdbvZs/xXuxecyY8Uy/FVVwz44uN0u3G4XpccZS78g6OPyPBthUGSicrlchIM+wkEfp1QcPpGe\ng/VuB82dMRoT79LtaSLcug7vW9vp+N1zNAU8RCtLSE6Zimv5CubPnktpsAifV7cEZPTG8kr/OsBv\nWdZyY8y5wA/6l2GM8QE/BJYAvcB6Y8zTwAVAYLBt8knGzhBL9dCVjNGd7KY7GaOr/2s02UVbXzut\n8XbiPV0ssvoo6XGY2+OipCdDqCeBq7CAKf/yT0SCR310rwRi195CQdBHoLo8Nz+ciIyaOTWCOTUC\n1GA7DvEzPs2Btiivv/kq3pY9FLQ2UvReA79tbuWJdhde/BR7Syn0lOB3ijhrdwt2sJjiadMoLpvG\n1KpTKIiU4vWq81Y+2lj+hZwPPA9gWdZrxpglA9rmAbsty4oCGGPWAZ8AzgOeG2KbQdmOTcbOYDs2\nNg62Y+M4NrbjYGP3v3b617Ox7XS23baxMxlsJ4Nj29gBH2k7TSwe55DdhLerkYPJbjZsrSeTSpJJ\nJcikU9ipJPFEku0BDw0HobOvhy2bXPSl4yQTfSzY3IwnkcGfsvGnHQIp8OLmjevnUeQvpMRfzJnl\n8ygPRSj3FuHa+BjpwjK6S4pwl1eQjFRwxjnz8B5d8EXkpOTu7w2oqSgldM5FwEUEfG6Kwn7qUkne\n2L+XWKaLrkwnrb2txNKHcNe/Q3FPkoK+7EygzYDtgodXz8HrCeJzBfC7/PhcQbz4mLdpD/6CEK5A\nkIzXT6iygkWXfpruzgS4sh//DQW8hAJHSsLhTkQXH/znSNsH6+g2RL4Zy6JfDHQNeJ0xxrgty7L7\n26ID2rqBkuNsM6i1q1bhcgCn/3PrjoPjcvHTVccOQuO2Hb7y65Zjlmdc8JPVU7Lr4MaNl10tXki6\nuOiRfcesn/a42PK55YQ8RaSdEIG+IEUuP163l9rOl3H8AZxAGFdRGHwhnECIRZnzIO6CePYnbCf7\njytuI52xKQz5KAz5aGzv5aDV3f8rOVYynWFm1dh8/M56t5OAT12IIrkQCnipmVJ4zLLLzpx/7MoX\nZb+0R6NE25s5ePAAmd5Ozp0WIZGJk8gkiGfiJOwE3bFWAh3NeJrTeFMZfCmb9hIvayLTCSUqcMg+\n8zMah4v/h04B+l8c3veYHFvy5JzDdYIDvfmaMz/2tmNZ9LuAogGvBxbv6FFtRUDncbYZ1KWPPzHo\nb/PioTZYPfjiTwy1/l8MvnjFUOvf8JmhWiasykqYXZu72wWVlUXHX0kmJOUut7K//+nA4o9ecdV4\nRCP5YCw/l7UeuBrAGLMM2DqgrR6YY4yJGGP8ZGvuhuNsIyIiIqPgcpzRdesMxRjj4siT+ABfJHs6\nWmhZ1v3GmGuAb5M98XjQsqyfDraNZVk7xyRAERGRSWbMir6IiIhMLBp2TUREZJJQ0RcREZkkVPRF\nREQmCRV9ERGRSeKkHbPRGHMJsNqyrNtzHYsMjzFmOXBH/8uvHh6xUfKH3nf5yRhzKXATEAa+a1mW\nPi6dR4wxi4G7yQ5XdI9lWYeGWvekvNI3xswCFgHBXMciI3I72aL/INkDkOQRve/yWsiyrDuA7wNX\n5DoYGbEA8DXgGbLD2Q/ppCz6lmXtsSzrh7mOQ0bMY1lWEmgEqnIdjIyM3nf5y7KsNcaYAuCvgYdz\nHI6MkGVZG4AzgL8F3vqodfOme79/1r1/syzr4qGm4DXG/CMwG7jTsqzOHIYrRxlO/oDe/hEaq4Gm\n3EUrRxtm/mQCGuaxswL4LvBty7JacxiuHGWY+VsKbAI+CXwH+OpQ+8uLK31jzD3A/WS7MGDAtL3A\nN8hOwYtlWd+yLGu1Cv7EMtz8AfcBPyPbzf/oeMcpgxtB/mSCGUHufgBMBf7VGHPjuAcqgxpB/gqB\nh4DvAb/4qH3my5X+buAGjhSCCxh62t4PWJZ18/iEJ8cxrPxZlrWZ7HDNMrGM6P2n992EMtz33q25\nCU+OY7j5WwusHc4O8+JK37KsJ4H0gEVFDDIF7/hGJcOl/OU35S9/KXf5bSzyl6/JHvEUvDKhKH/5\nTfnLX8pdfht1/vK16GsK3vym/OU35S9/KXf5bdT5y5d7+ocdnhLwKeByY8z6/te6D5wflL/8pvzl\nL+Uuv52w/GlqXRERkUkiX7v3RUREZIRU9EVERCYJFX0REZFJQkVfRERkklDRFxERmSRU9EVERCYJ\nFX0REZFJQkVfRERkklDRFxERmSRU9EVERCaJfBt7X0RyyBhzHfBdYA9wCxAA/gQ8B/zIsqw9uYtO\nRI5HY++LyIgYY74EXGtZ1rXGmOWAy7Ks9cfbTkRyT0VfREbEGFMIvAt8FghZlrUmxyGJyDCp6IvI\niBlj7gd8lmV9IdexiMjw6UE+ERkRY0wYaAfOz3UsIjIyKvoiMmzGGDdwO/BNIGaMuTjHIYnICKjo\ni8iwGGNcwF8BD1uWlQIeAL6c26hEZCRU9EXkuIwxVwPPA8sty4r2L64GrjfGfD13kYnISOhBPhER\nkUlCV/oiIiKThIq+iIjIJKGiLyIiMkmo6IuIiEwSKvoiIiKThIq+iIjIJKGiLyIiMkn8P537VmkW\nhmS8AAAAAElFTkSuQmCC\n",
"text": [
"<matplotlib.figure.Figure at 0xa922438>"
]
}
],
"prompt_number": 11
}
],
"metadata": {}
}
]
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment