{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# How to generate random real numbers with a power-law distribution?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Notes from Appendix D of:
\n", "Clauset, A., Shalizi, C. R., & Newman, M. E. J. (2007). Power-law distributions in empirical data. SIAM Review, 51(4), 43. https://doi.org/10.1109/ICPC.2008.18" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First, generate set of random real numbers $r$ that are uniformaly distributed in the interval $[0, 1)$." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from math import floor, ceil" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "r = np.random.uniform(low=0.0, high=1.0, size=1000)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's plot to see if they're uniformly distributed." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "NumBins = 10\n", "MyBins = np.linspace(0, 1, NumBins+1)\n", "\n", "r_hist, r_bins = np.histogram(r, bins=MyBins, density=True)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAewAAAFWCAYAAACre65zAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAATVUlEQVR4nO3df6xkd3nf8c+DXYfwy1Tspq38g3VS07JyIzlaOUSkxSmktUGyVRVSu6WNU5eNkppUMkW4SetSRyoEGqVNawKmRE6ixj9CVboNSx0lMQLRGHnBYMV23W6ME28d4g0xphUxxuLpHzPA9eWud+7du3Pmu/f1kq40Z+a7c589vuu3zrkzZ6q7AwCstudMPQAAcHyCDQADEGwAGIBgA8AABBsABiDYADCA06ce4GTYtWtX79mzZ+oxAGBTPvWpT/1xd+/e6LFTMth79uzJoUOHph4DADalqn7/WI85JQ4AAxBsABiAYAPAAAQbAAYg2AAwAMEGgAEINgAMQLABYACCDQADEGwAGIBgA8AABBsABnBKfvgHwCrbc92Hpx7hGR5+5+umHoEFOMIGgAEINgAMQLABYAB+hz04vwsD2BkcYQPAAAQbAAYg2AAwAMEGgAEINgAMQLABYACCDQAD8D7sBazSe529z3nnWKWfu8TPHkzNETYADMARNrAlzgDAcjnCBoABCDYADECwAWAAgg0AAxBsABiAYAPAAAQbAAYg2AAwAMEGgAEINgAMQLABYACCDQADmDTYVfWLVfVYVf3uMR6vqvr5qjpcVfdW1fcse0YAWAVTH2HfnOSSZ3n80iTnz7/2J/mFJcwEACtn0mB398eS/MmzLLk8yS/3zF1JXlxVf2E50wHA6pj6CPt4zkryyJrtI/P7AGBHWfVg1wb39YYLq/ZX1aGqOnT06NGTPBYALNeqB/tIknPWbJ+d5NGNFnb3Td29r7v37d69eynDAcCyrHqwDyT5B/NXi78iyRPd/YdTDwUAy3b6lN+8qm5JcnGSXVV1JMm/TPJnkqS735vkYJLXJjmc5MtJfmSaSQFgWpMGu7uvPM7jneQfL2kcAFhZq35KHACIYAPAECY9Jc6pb891H556hGd4+J2vm3oEgC1xhA0AAxBsABiAYAPAAAQbAAYg2AAwAMEGgAEINgAMQLABYACCDQADEGwAGIBgA8AABBsABiDYADAAwQaAAfh4TXa0Vfr4Tx/9ub38t926Vdp3yXj772RxhA0AAxBsABiAYAPAAAQbAAYg2AAwAMEGgAEINgAMQLABYACCDQADEGwAGIBgA8AABBsABiDYADAAwQaAAQg2AAxAsAFgAIINAAMQbAAYwOlTDwAAm7Hnug9PPcI3PPzO1y3teznCBoABCDYADECwAWAAkwe7qi6pqger6nBVXbfB4+dW1Z1VdU9V3VtVr51iTgCY0qTBrqrTktyY5NIke5NcWVV71y3750lu7+4Lk1yR5D3LnRIApjf1EfZFSQ5390Pd/VSSW5Ncvm5NJ3nR/PaZSR5d4nwAsBKmflvXWUkeWbN9JMn3rlvz9iS/UVVvTvL8JK9ZzmgAsDqmPsKuDe7rddtXJrm5u89O8tokv1JV3zJ3Ve2vqkNVdejo0aMnYVQAmM7UwT6S5Jw122fnW095X53k9iTp7t9J8twku9Y/UXff1N37unvf7t27T9K4ADCNqYN9d5Lzq+q8qjojsxeVHVi35g+SvDpJqurlmQXbITQAO8qkwe7up5Nck+SOJA9k9mrw+6rqhqq6bL7sLUneVFWfTXJLkqu6e/1pcwA4pU39orN098EkB9fdd/2a2/cneeWy5wKAVTL1KXEAYAGCDQADEGwAGIBgA8AABBsABiDYADAAwQaAAQg2AAxAsAFgAIINAAMQbAAYgGADwAAEGwAGINgAMADBBoABCDYADECwAWAAgg0AAxBsABiAYAPAAAQbAAYg2AAwAMEGgAEINgAMQLABYACCDQADEGwAGIBgA8AABBsABiDYADAAwQaAAQg2AAzg9M0srqpK8pokP5jkryU5N8muJH+a5LEkn0ny20kOdPf/2d5RAWDnWijYVfW8JD+R5Eczi3TNH3oys1B/e5LvTPJdSf52kn9XVf8tyc929//Y7qEBYKc57inxqvqRJP87yb/O7Ej6X2V2hP3i7n5ed5/d3S/JLP57k/zDJP85yaVJPl5Vt1XVuSfrLwAAO8EiR9gfSPKhJO/o7ruPtai7O8n/nH/dXFUvSvLDSa5LclWSG054WgDYoRYJ9r7u/vRmn7i7v5Tk31fV+5Ps2eyfBwC+6bjB3kqs1/35JzM76gYAtmjTb+uqqoeq6saTMQwAsLGtvA97d5IntnsQAODYthLs+zJ7+9a2qKpLqurBqjpcVdcdY80PVdX9VXVfVf3qdn1vABjFpi6cMvfzSf5jVX13d997It+8qk5LcmNmbxM7kuTuqjrQ3fevWXN+kn+W5JXd/XhVfceJfE8AGNFWgn0kyW8m+URVvS/J3Uk+n6TXL+zujx3nuS5Kcri7H0qSqro1yeVJ7l+z5k1Jbuzux+fP+dgWZgaAoW0l2B/NLM6V5NpsEOo1TjvOc52V5JE120eSfO+6NS9Lkqr6xPz53t7d/339E1XV/iT7k+Tcc12nBYBTy1aCfUOePdKbURvct/65T09yfpKLk5yd2dXTLujuLz7jD3XflOSmJNm3b992zQcAK2HTwe7ut2/j9z+S5Jw122cneXSDNXd191eTfK6qHsws4Me86hoAnGqm/njNu5OcX1XnVdUZSa5IcmDdmg8l+YEkqapdmZ0if2ipUwLAxCYNdnc/neSaJHckeSDJ7d19X1XdUFWXzZfdkeQLVXV/kjuTvLW7vzDNxAAwjeOeEp9/TOa/6O7PbPbJq+rbkvx4kie7+xc2WtPdB5McXHff9Wtud2Yvbrt2s98fAE4Vixxh/+Ukn6qqj1TV36mq5x7vD1TVy6vqHUk+l+RnkvzfE5wTAHa0RV50tjfJP0nyk0n+RpKnqurTSQ4l+cMkjyd5bpKXZBb3V2T2dq1K8htJ/ml3/+72jw4AO8cin9b11ST/pqrek+TvJbk6swuefN/Xl+SZb886muT9Sd5zoldCAwBmFn5bV3d/ObMQv7+qXpRZsM/N7Mj6T5M8luTe7r7vZAwKADvZVi6cku7+Umav3gYAlmBLwU6SqnpBkr+V5MIkZ2b2kZv3JPkv3f3/tmc8ACDZYrCr6g1J3pvkxXnm7687yb+tqh/t7g9uw3wAQLYQ7Kr6wSS3JPlakl/O7MNAPp/kz2d2RbK/m+SWqvpid//m9o0KADvXVo6wr0/ylSR/tbs/ve6xX6qq/5DkY/N1gg0A22Arlya9MMltG8Q6SdLdh5LcnuR7TmQwAOCbthLsr2R2wZRn8+h8HQCwDbYS7I8n+f7jrHllZqfFAYBtsJVgvy3JX6mqd1bV89c+UFXPr6p3JbkgyXXbMSAAsLUXnb0tyb1J3ppk//y64n+U5M9l9nvrMzM7un5b1dp3fKW7++oTGxcAdqatBPuqNbdfnOSvb7DmVfOvtTqz65ADAJu0lWCft+1TAADPatPB7u7fPxmDAADHtpUXnQEASybYADAAwQaAAQg2AAxAsAFgAIINAAMQbAAYgGADwAAEGwAGINgAMADBBoABCDYADECwAWAAgg0AAxBsABiAYAPAAAQbAAYg2AAwAMEGgAEINgAMQLABYACCDQADmDzYVXVJVT1YVYer6rpnWff6quqq2rfM+QBgFUwa7Ko6LcmNSS5NsjfJlVW1d4N1L0zyE0k+udwJAWA1TH2EfVGSw939UHc/leTWJJdvsO6nk7wryZPLHA4AVsXUwT4rySNrto/M7/uGqrowyTnd/evLHAwAVsnUwa4N7utvPFj1nCQ/l+Qtx32iqv1VdaiqDh09enQbRwSA6U0d7CNJzlmzfXaSR9dsvzDJBUk+WlUPJ3lFkgMbvfCsu2/q7n3dvW/37t0ncWQAWL6pg313kvOr6ryqOiPJFUkOfP3B7n6iu3d1957u3pPkriSXdfehacYFgGlMGuzufjrJNUnuSPJAktu7+76quqGqLptyNgBYJadPPUB3H0xycN191x9j7cXLmAkAVs3Up8QBgAUINgAMQLABYACCDQADEGwAGIBgA8AABBsABiDYADAAwQaAAQg2AAxAsAFgAIINAAMQbAAYgGADwAAEGwAGINgAMADBBoABCDYADECwAWAAgg0AAxBsABiAYAPAAAQbAAYg2AAwAMEGgAEINgAMQLABYACCDQADEGwAGIBgA8AABBsABiDYADAAwQaAAQg2AAxAsAFgAIINAAMQbAAYgGADwAAEGwAGINgAMIDJg11Vl1TVg1V1uKqu2+Dxa6vq/qq6t6p+q6peOsWcADClSYNdVacluTHJpUn2JrmyqvauW3ZPkn3d/d1JPpjkXcudEgCmN/UR9kVJDnf3Q939VJJbk1y+dkF339ndX55v3pXk7CXPCACTmzrYZyV5ZM32kfl9x3J1ko9s9EBV7a+qQ1V16OjRo9s4IgBMb+pg1wb39YYLq96YZF+Sd2/0eHff1N37unvf7t27t3FEAJje6RN//yNJzlmzfXaSR9cvqqrXJPmpJK/q7q8saTYAWBlTH2HfneT8qjqvqs5IckWSA2sXVNWFSd6X5LLufmyCGQFgcpMGu7ufTnJNkjuSPJDk9u6+r6puqKrL5sveneQFSX6tqj5TVQeO8XQAcMqa+pR4uvtgkoPr7rt+ze3XLH0oAFgxU58SBwAWINgAMADBBoABCDYADECwAWAAgg0AAxBsABiAYAPAAAQbAAYg2AAwAMEGgAEINgAMQLABYACCDQADEGwAGIBgA8AABBsABiDYADAAwQaAAQg2AAxAsAFgAIINAAMQbAAYgGADwAAEGwAGINgAMADBBoABCDYADECwAWAAgg0AAxBsABiAYAPAAAQbAAYg2AAwAMEGgAEINgAMQLABYACCDQADEGwAGMDkwa6qS6rqwao6XFXXbfD4t1XVbfPHP1lVe5Y/JQBMa9JgV9VpSW5McmmSvUmurKq965ZdneTx7v6LSX4uyc8sd0oAmN7UR9gXJTnc3Q9191NJbk1y+bo1lyf5pfntDyZ5dVXVEmcEgMlNHeyzkjyyZvvI/L4N13T300meSPKSpUwHACuiunu6b171hiR/s7v/0Xz77ye5qLvfvGbNffM1R+bbvzdf84V1z7U/yf755l9K8uAS/gqbtSvJH089xKDsu62z706M/bd19t3mvbS7d2/0wOnLnmSdI0nOWbN9dpJHj7HmSFWdnuTMJH+y/om6+6YkN52kObdFVR3q7n1TzzEi+27r7LsTY/9tnX23vaY+JX53kvOr6ryqOiPJFUkOrFtzIMkPz2+/Pslv95SnBQBgApMeYXf301V1TZI7kpyW5Be7+76quiHJoe4+kOQDSX6lqg5ndmR9xXQTA8A0pj4lnu4+mOTguvuuX3P7ySRvWPZcJ8lKn7Jfcfbd1tl3J8b+2zr7bhtN+qIzAGAxU/8OGwBYgGBvM5daPTEL7L9rq+r+qrq3qn6rql46xZyr6Hj7bs2611dVV5VX784tsu+q6ofmP3v3VdWvLnvGVbbAv9tzq+rOqrpn/m/3tVPMObzu9rVNX5m9cO73knxnkjOSfDbJ3nVrfjzJe+e3r0hy29Rzr8rXgvvvB5I8b377x+y/xffdfN0Lk3wsyV1J9k099yp8Lfhzd36Se5L82fn2d0w996p8Lbj/bkryY/Pbe5M8PPXcI345wt5eLrV6Yo67/7r7zu7+8nzzrszeu89iP3tJ8tNJ3pXkyWUOt+IW2XdvSnJjdz+eJN392JJnXGWL7L9O8qL57TPzrdfbYAGCvb1cavXELLL/1ro6yUdO6kTjOO6+q6oLk5zT3b++zMEGsMjP3cuSvKyqPlFVd1XVJUubbvUtsv/enuSNVXUks3cFvTls2uRv6zrFbHSkvP5l+Ius2akW3jdV9cYk+5K86qRONI5n3XdV9ZzMPu3uqmUNNJBFfu5Oz+y0+MWZndX5eFVd0N1fPMmzjWCR/Xdlkpu7+2er6vsyu7bGBd39tZM/3qnDEfb22sylVvNsl1rdoRbZf6mq1yT5qSSXdfdXljTbqjvevnthkguSfLSqHk7yiiQHvPAsyeL/bv9rd3+1uz+X2WcVnL+k+VbdIvvv6iS3J0l3/06S52Z2nXE2QbC3l0utnpjj7r/5ad33ZRZrv0f8pmfdd939RHfv6u493b0ns9//X9bdh6YZd6Us8u/2Q5m94DFVtSuzU+QPLXXK1bXI/vuDJK9Okqp6eWbBPrrUKU8Bgr2N5r+T/vqlVh9IcnvPL7VaVZfNl30gyUvml1q9Nskx336z0yy4/96d5AVJfq2qPlNV6//HsCMtuO/YwIL77o4kX6iq+5PcmeStve4TA3eqBfffW5K8qao+m+SWJFc5UNk8VzoDgAE4wgaAAQg2AAxAsAFgAIINAAMQbAAYgGADwAAEGwAGINgAMADBBoABCDawkKraU1VdVTdX1cuq6raqeqyqvlZVF089H5zqfLwmsFnfleSTSf5Xkv+U5NuTfGnSiWAHEGxgs74/yTu6+yenHgR2Eh/+ASykqvYk+VySP0ryUp9FDsvld9jAZn1WrGH5BBvYrM9PPQDsRIINbJbfo8EEBBsABiDYADAAwQaAAQg2AAzA+7ABYACOsAFgAIINAAMQbAAYgGADwAAEGwAGINgAMADBBoABCDYADECwAWAAgg0AA/j/qGDMkq4xKH4AAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "fig = plt.figure()\n", "ax = fig.add_axes([0, 0, 1, 1])\n", "\n", "ax.bar( r_bins[:-1], r_hist, width=0.9*(MyBins[1] - MyBins[0]) )\n", "\n", "ax.set_xlabel(\"r\", fontsize=20)\n", "ax.set_ylabel(\"p(r)\", fontsize=20)\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let $p(r)$ be the uniform probability density function, from which we have drawn the above random real numbers." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, suppose $p(x)$ is a continuous probability density from which we want to draw\n", "random real numbers $x \\geq x_{min}$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Probability densities $p(r)$ and $p(x)$ are related as follows:\n", "\n", "$$p(x)dx = p(r)dr$$\n", "\n", "Note that $p(r) = 1$ for $r \\in [0, 1)$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\\Rightarrow p(x) = \\frac{dr}{dx}$$\n", "\n", "Integrate both sides wrt $x$:\n", "\n", "$$\\Rightarrow \\int_{x}^{\\infty} p(x') dx' = \\int_{r}^{1} dr'$$\n", "$$\\Rightarrow P(x) = 1 - r$$\n", "$$\\Rightarrow x = P^{-1}(1 - r)$$\n", "\n", "Where, $P$ is the CDF, and $P^{-1}$ is the functional inverse of the CDF." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Recall that for power-law distribution, PDF is:
\n", "$$p(x) = \\frac{\\alpha - 1}{x_{min}} \\left(\\frac{x}{x_{min}}\\right)^{-\\alpha}$$\n", "\n", "So, CDF is:\n", "$$P(x) = \\int_{x}^{\\infty} p(x')dx'\n", "= \\int_{x}^{\\infty} \\frac{\\alpha - 1}{x_{min}} \\left(\\frac{x'}{x_{min}}\\right)^{-\\alpha} dx'\n", "= \\left(\\frac{x}{x_{min}}\\right)^{-\\alpha+1}$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So:\n", "$$ P(x) = 1 - r $$\n", "$$ \\Rightarrow \\left( \\frac{x}{x_{min}}\\right)^{-\\alpha+1}\n", "= 1 - r $$\n", "$$ \\Rightarrow x = x_{min} (1 - r)^{1/(-\\alpha+1)}$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let $x_{min} = 7$, $\\alpha = 2.5$" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "xmin = 7.\n", "alpha = 2.5" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "x = xmin * ((1-r)**(1/(-alpha+1)))" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "7.000574541399359 2987.649692743292\n" ] } ], "source": [ "print(min(x), max(x))" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "7 2988\n" ] } ], "source": [ "print(floor(min(x)), ceil(max(x)))" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "MyBins_x = np.arange(floor(min(x)), ceil(max(x))+1, 1)\n", "\n", "x_hist, x_bins = np.histogram(x, bins=MyBins_x, density=True)" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "px_theory = (alpha-1)/(xmin) * (x_bins[:-1]/xmin)**(-alpha)" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "PlotYlim_bottom = (sorted(set(x_hist))[1] - sorted(set(x_hist))[0])/2" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfUAAAFaCAYAAAAHAsQqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deZxf873H8ddnZhK7SUQUQYitqbRFgtBbpZpIQnSvopulunDbWq6lyq2llipdbJVLr95LaasLEsG1NUoGM60llhRhiC3ECIpMJvO9f/wymkxmz2/m/H6/83o+Hnkkc37nd857Jkfezu+c8/1GSglJklT+qrIOIEmSisNSlySpQljqkiRVCEtdkqQKYalLklQhLHVJkipETdYBsrL++uunzTffPOsYkiT1SkNDw6sppeEdvZbbUt98882pr6/POoYkSb0SEY2dvebH75IkVQhLXZKkCmGpS5JUISx1SZIqRO5KPSKmRsS0RYsWZR1FkqSiyl2pp5RuSCkdXltbm3UUSZKKKnelLklSpbLUJUmqEJa6JEkVwlKXJKlCWOqSJFUIS12SpAphqUuSVCEsdUmSKoSlLklShbDUJUmqELkrdcd+lyRVqtyVumO/S5IqVe5KXZKkSmWpS5JUISx1SZIqhKUuSVKFsNQlSaoQlrokSRXCUpckqUJY6pIkVQhLXZKkCmGpS5JUISz1Ykgp6wSSJFnqq6zxHrh4V3jzpayTSJJyzlJfRXMWrU7Lq0/x2h+PzTqKJCnnLPVV0NDYxOd+/zIXLpnKek/fwD/uuT7rSJKkHLPUV0HdvIU0t7RySctUnml9HxvcdRIseTfrWJKknMpdqUfE1IiYtmjRolXe1vhRwxhcU0VLDOa0dChD3nkW7v55EVJKktR7kXJ65/a4ceNSfX39Km+nobGJunkLGT9qGGPvOxoenwHfng3DtixCSkmSVhQRDSmlcR29lrsz9WIbO3IoR+y5FWNHDoW9z4TqwXDjsT7mJkkacJZ6Ma27Eex1Mjx1Ozzyp6zTSJJyxlIvtp0Og40+DDedCO++kXUaSVKOWOrFVlUN+/4U3noZ7jgz6zSSpByx1PvDiLEw7hC471J48cGs00iScsJS7y97nQJrDoPpR0Hr0qzTSJJywFLvL2sMKdwN/3wDNFyRdRpJUg5Y6v3pg5+HLXaH206FtxZknUaSVOEs9f4UAVPOg+a34ZaTs04jSapwlnp/G74NfOS78NA18PSsrNNIkiqYpT4Qdj8WhoyEGcdAS3PWaSRJFcpSHwiD1oApP4FX/wH3/CLrNJKkCmWpD5RtJsLo/WDWudD0TNZpJEkVyFIfSJPOhqoauPE4J3yRJBWdpT6QakfAHifCEzfD49OzTiNJqjCW+kDb5ZvwvjEw83hY/FbWaSRJFcRSHwANjU1cdMeTNDQ2QXUN7HM+vPE83HlW1tEkSRWkJusAla6hsYmDLqujuaWVwTVVXHXYeMaO3AV2/CrUXQIfPgA2HJN1TElSBfBMvZ/VzVtIc0srrQmWtLRSN29h4YVP/LAwPvyMo6G1NcuIkqQKYan3s/GjhjG4porqgEE1VYwfNazwwprrwYTT4bl74YErsw0pSaoIkXL6aNW4ceNSfX39gOyrobGJunkLGT9qGGNHDv3XCynBf0+BVx6DIxtgrWEDkkeSVL4ioiGlNK6j1zxTHwBjRw7liD23WrHQoTDhy77nw+I34f9OySacJKliWOpZ22A07Hpk4SP4xtlZp5EklbGKKPWIGBURl0fEtVln6ZOPHQe1mxVumlu6JOs0kqQylXmpR8SvImJBRMxpt3xSRMyNiCcj4oSutpFSmpdSOrR/k/ajwWvB5HNgwaNQd3HWaSRJZSrzUgeuACYtvyAiqoGLgMnAB4ADIuIDEfHBiJje7tcGAx+5H7x/Cmw7Be48G15/Lus0kqQylHmpp5RmAa+1W7wz8OSyM/Bm4Brgkymlh1NK+7b7tWDAQ/eXyecUfp95fLY5JEllKfNS78QIYPnT1fnLlnUoIoZFxC+BHSLixC7WOzwi6iOi/pVXXile2mIZshl87HiYOwPmzsw6jSSpzJRqqUcHyzp9oD6ltDCl9M2U0pYppU4HVE8pTUspjUspjRs+fHhRghbDCmPD73oEDB9dmJ61+Z9ZR5MklZFSLfX5wKbLfb0J8EJGWfpV29jw590yl4Muq6Nh/luFZ9cXPQuzzs06niSpjJRqqd8PbB0RW0TEYOCLwPUZZ+oXHY4NP3I32P4guOcCWPBY1hElSWUi81KPiKuB2cC2ETE/Ig5NKbUARwI3A48Bv0spPZJlzv7S6djwE06DwWvDjGMKw8lKktSN3I39HhFTgalbbbXV15944oms4wBdjA3f8Gu44TvwqUtg+wOzCyhJKhldjf2eu1JvM5ATuvRZaytvXbIX1a/PY+7n72T7bbbIOpEkKWNO6FKmGp5bxIEv7c+g5jeYe+UxhbvjJUnqhKVewurmLWROy6b899JJ7F91G08/cGfWkSRJJcxSL2FtN9FdsPSzvJTWY0rjObC0JetYkqQSlbtSj4ipETFt0aJFWUfp1tiRQ7nqsPF8Y+L2/PPjZ7Lma4/BfZdmHUuSVKK8Ua5cpAS/+QI03gNH3Ae1nY6aK0mqYN4oVwkiYPKPobUFbu50eHtJUo5Z6uVkvS1g92Ph0evgiVuzTiNJKjGWernZ7TswbGu48RhY8k7WaSRJJSR3pV5ON8p1qGY12Oc8aHoG7jov6zSSpBKSu1JPKd2QUjq8trY26yh9N+pj8KH94a8/g1dLY6hbSVL2clfqleLBDxzLu1Wr88YfvuOEL5IkwFIvSw2NTex/1VP86J3Pse6L9zDvjl9nHUmSVAIs9TLUNgf7VUv34sHWUWw4+zR45/WsY0mSMmapl6G24WMjqjg1fZ01Wprg9tOzjiVJylhN1gHUe23DxxbmYN+NeLQR7r20MOf6iLFZx5MkZSR3w8RGxFRg6lZbbfX1J56okDvH330DLtwJ1tkQvn47VFVnnUiS1E8cJnY5FfFIW3urrwuTzoQXH2DWb8523nVJyqnclXqlalh7T+5OH2SHJy7ge5fNtNglKYcs9QpR9/RrnLzkawymheP4H+rmLcw6kiRpgFnqFWL8qGG8UD2CS5bux9Tq2UxY7dGsI0mSBpilXiHa7ohffc9jeHedzdmm/oew5N2sY0mSBpClXkHGjhzKN/fajtU/9VN47Sm4++dZR5IkDaDclXrZz9LWE1t+HLb7TGEWt4VPZZ1GkjRAclfqFflIW0f2PhOqB8ONxzrhiyTlRO5KPTfW3Qj2Ohmeuh0e+VPWaSRJA8BSr2Q7HQYbfRhuOrEw6pwkqaJZ6pWsqhr2/Sm89TLccWbWaSRJ/cxSr3QjxsK4Q+C+S+HFB7NOI0nqR5Z6Hux1Cqw5DKYfBa1Ls04jSeonlnoerDGkcDf88w3QcEXWaSRJ/cRSz4sPfh622B1uOxXeWpB1GklSP7DU8yICppwHzW/DLSdnnUaS1A9yV+q5GFGuM8O3gY98Fx66Bp6+K+s0kqQiy12p52ZEueU0NDZx0R1PFuZY3/1YGDISZhwNLc1ZR5MkFVHuSj1vGhqbOOiyOs67ZS4HXVZHwwvvwpSfwKv/gHt+kXU8SVIRWeoVrm7eQppbWmlNsKSllbp5C2GbiTB6P5h1LjQ9k3VESVKRWOoVbvyoYQyuqaI6YFBNFeNHDSu8MOlsqKqBG49zwhdJqhA1WQdQ/xo7cihXHTaeunkLGT9qGGNHDi28UDsC9jgRbjkJHp8Oo6eu8L6GxqaV3yNJKmmRcnqWNm7cuFRfX591jGwtbYFpH4N3muCI+2C1tYF/XYdvbmllcE0VVx023mKXpBIREQ0ppXEdvebH73lWXQP7nA9vPA93nvXe4g6vw0uSSp6lnneb7QI7fhXqLoGX5gBdXIeXJJU0P37PoZWul7/9Glw4DoZtBQffBFVVXlOXpBLV1cfv3iiXMx1fL18PJpwO130bHrgSdvwKY0cOtcwlqcz48XvOdHa9vGHoZF6o3YGWm0+Gf3oNXZLKUe5KPddjv9Px9fKGxiYOuvxeDn7lANK7b/Lqn4/POqYkqQ9yV+p5HPt9eW3PrR89cdv3HlVrO3uf27oJly+dwvpP/B4aZ6/wvhXGj5cklSSvqedQ++vlbWfvS1pauTQ+y8Fr/43VZhwN35gF1YN8bl2SykTuztS1suXP3i87bA9W2/cnsOBRqLsY8Ll1SSoXnqkLaH/2PgW2nQJ3ng3bfWaFM3mfW5ek0uWZujo2+ZzC7zed0OF1eElS6fFMXR0bshl87Hi49T9h7kzGbjvZMpekEueZujq36xEwfHRhetbmf2adRpLUDUtdnaseBPueD4uehVnnZp1GktQNS11dG7kbbH8Q3HMBLHg86zSSpC5Y6urehNNg8Now42jI6QRAklQOLHV1b631C8XeeDc8eHXWaSRJnbDU1TM7fBk22Rlu+UFhqlZJUsnpU6lHxDoRMSkiDo+I70fEURFxUERsV+yAKhFVVbDvT+Gd1+G2U7NOI0nqQI+fU4+INYADgEOBnfnX/xDEst/TsvVeBf4AXJJSerh4UZW5DcfA+G/B7Ath+y/BpjtlnUiStJxI3dz4FBE1wHeAk4ChwLtAA3A/8BLwGrAGMAx4PzAe2IxCyd8KHJNSmtNP+fts3Lhxqb6+PusY5Wfxm3DhzrDmMDj8Tqh2/CJJGkgR0ZBSGtfRaz35F/lxYAvgJuDXwHUppcXd7HAb4GvAV4C/R8ShKaX/6VVqlabV1oHJZ8PvvgL3XVoYoEaSVBJ6ck39UWDHlNI+KaXfdVfoACmlf6SUvg+MAv6Dwpl8SYiIqRExbdGiRVlHKV+j92PRJnvSfOsZPPToo1mnkSQt022pp5T2Syk92JeNp5SaU0o/Syld2pf394eU0g0ppcNra2uzjlK2Gp59nc82forWliW88Nvv0dDYlHUkSRI+0qY+qJu3kHktw7mg5dNMint5of6GrCNJkuhDqUfER3u43nd7H0floG1+9V+17sO8tDETnzkXlryTdSxJyr2+nKnfHhE/6OzFiBgaEdcB5/c9lkpZ2/zqR04cQ8ukn7Dam8/CXedlHUuScq8vzyM9CZwaEXsAX0opvdT2QkTsBlwNbAr8qSgJVZLGjhy6bH71rVj45KcZctdPeWz9SYz5UIdPWUiSBkBfztTHAlcCH6fwuNoEgIg4EbgT2AA4MqX02WKFVOlqaGxiv3/szVutg3nzD9+h4RmHkJWkrPS61FNKb6eUvkrhOfS1gZkR8RBwBvAUsEtK6eKiplTJqpu3kBdb1uXclv3ZNR5hYd1VWUeSpNzq893vywaTOWHZNsYArwK7p5QeKlI2lYG2m+Z+27oXD6Ut2bPxZ4Xx4SVJA66vE7pURcSPgJ8DbwH3AMOBOyPig0XMpxLXdtPc9yaOpma/nzPo3dfg9jOyjiVJudSXR9o2BWYBJwIPA+NSSv9GYWz4rYF7I+LbRU2pkjZ25FCO2HMrPjD2o7Dz4XD/ZfB8Q9axJCl3+nKm/gCwG3AJMD6l9A+AlNJZwB7AK8AFEfHHYoVUGdnzJFj7fTD9aGhdmnUaScqVvpR6FfC5lNIR7ceBTyndA3wYuB74ZBHyqdysvi5MOhNefADuvzzrNJKUK30p9R1SSp2ehaeUXk8pfRpwRLm82u4zMGpPuP10ePOl7teXJBVFXx5pe6aH613Y6zSqDBGwz3nQshhuPinrNJKUG07oov4xbEv46NEw51p46vas00hSLnRb6hFxQ0Rs35eNR8RqEXFURHyrL+9XmfvI92C9UTDjWFjybtZpJKni9eRM/f1AQ0TMjIj9I2L17t4QEaMj4izgaeAc4M1VzKlyNGj1wsfwrz0Fd/886zSSVPF6MqHLByjc9PZ9YCLQHBF/A+qBF4EmYHVgGIX/ARgPjAACuAU4NqU0p/jRVRa2/Hjhxrm7zoMPfq7wsbwkqV9ESqlnK0asCRwEHEphUpfqZS8lCgXe5hUKM7RdXMpDxo4bNy7V19dnHSMf3ngRLtwJNt0JvvTHwo10kqQ+iYiGlFKHU2L2eOrVlNLbwH8B/xUR6wK7AptROEN/B1gAPJRSemTVI6uirLsR7HUyzDwOHvkTjPlM1okkqSL1ZT51UkpvADcXOYsq2U6HwQNXwU0nwlafKAxSI0kqqj4/0hYRa0fElyPi/Ii4fNnvX46ItYsZUBWiqhr2/Sm89TLccWbWaSSpIvXpTD0iPg/8EhjCitfTE/CziPhGSunaIuRTJRkxFsYdAvddCtsfABt9OOtEklRR+jJL2wTgamAd4H+Ag4HJy37/32XLr46ITxQxpyrFXqfAmsNg+lFO+CJJRdaXj99PARZTmKHt4JTSr1NKNy/7/WsUZnBbsmw9aUVrDIGJPypMzdpwRdZpJKmi9GlCF+C3KaW/dfRiSqke+B2w46oE642I+FRE/FdEXBcREwdqv+qjD30BNv8o3HYqvLUg6zSSVDH6UuqLKQw605UXlq3XrYj4VUQsiIg57ZZPioi5EfFkRJzQ1TZSSn9OKX0d+Bqwf0/2q/7T0NjERXc8SUNj0wp/fk8E7HM+NL/Nwj8dv/Lr7bYhSeqZvtwodxfwb92s8xFgVg+3dwVwIYXr8wBERDVwETABmA/cHxHXUxjw5qx27z8kpdR2uveDZe9TRhoamzjosjqaW1qpqQqIoGVpK4NrqrjqsPGMHTm0sOLwbXhxzDfY6KEL+etjH+CC6jHvvb78NlZ6nySpU305Uz8e+GBEnB0Ray3/QkSsFRE/BsYAXZ5dt0kpzQJea7d4Z+DJlNK8lFIzcA3wyZTSwymlfdv9WhAF5wAzO7sssCzf4RFRHxH1r7zySi++ZfVU3byFNLe00ppgydLEkrY/t7RSN2/hCutet+4BPJuGc1rNf0NL83uvr7CNDt4nSepYX87UjwceAv4DOHzZOPAvA++jcB29lsJZ+vGx4nCgKaV0aA/3MQJ4brmv5wO7dLH+vwOfAGojYquU0i87WimlNA2YBoVhYnuYRb0wftQwBtdUsaSlleplZ+pLl7YyqKaK8aOGrbDuTluP4Iy/HMK06nM4fNCNjB+1+0rb6Oh9kqSO9Xjs9/feENHax32llFJ1Ry9ExObA9JTSmGVffx7YO6V02LKvvwzsnFL69z7ueyWO/d5/GhqbqJu38L0ybvtzRx+hNzQ2UXvDoYxq+itVR94HQzdfaRt+9C5J/1KUsd+Xs8Uq5umJ+cCmy329CYWb71QGxo4cukIRd1XKY0cOhS9fABftzOt/OIqrRv2Y8Vuuv9I2LHlJ6l6vSz2l1NgfQdq5H9g6IrYAnge+CBw4APtVFmpH8NyHv8em95/Bw09fxQV37LLCzXHeOCdJPdPnsd+LJSKuBmYD20bE/Ig4NKXUAhxJYdKYx4DfFWv2t4iYGhHTFi1aVIzNqUhuWGMqj7Vuxik1v2ZQy9sr3BznjXOS1DOZl3pK6YCU0kYppUEppU1SSpcvW35jSmmblNKWKaUfFXF/N6SUDq+trS3WJlUEu2z5Pk5Lh7JxvMb3Bv1xhZvj2m6cqw68cU6SutCnCV2kYhs7cijHHvZVHrnpAQ556XpiteeBoe+9dtVh472mLknd6PXd75XCu9+z1emNb2+/BheOg2FbwcE3QVXmHyZJUknp6u53/8XUgGu78e28W+Zy0GV1Kw4Fu+Z6MOF0eO5eeODK7EJKUhnKXal7o1z2ur3xbfsDYbPd4P9OgX96U5wk9VTuSt0b5bLX7Y1vEbDv+bD4TbjVGXwlqae8UU4Drkc3vm0wGnY9Eu7+GWz/JRi568AHlaQy441yKl3N/4SLxsNqa8M3ZkH1oKwTSVLmvFFOJaHXc6QPXgsmnwMLHoW6i/s3nCRVAD9+14Do81Cv758C206BO8+G7T4DQzbt/j2SlFO5O1P37vdsrNJQr5PPKfx+0wn9E06SKkTuSt2737OxSkO9DtkMPnY8PD4d5s7sv5CSVOa8UU4DZpWmT126BH750cLNc0fUFa63S1IOeaOcSsLYkUM5Ys+t+jZ2e/WgwrPri56FWecWP5wkVQBLXSVthTvmR+4G2x8E91wACx7vel1JyiHvflfJ6vCO+QmnweMzYMYx8LXphdHnOlvX2dwk5UzuztS9+718dHjH/Frrw4RTofGv8OA1Xa8rSTmTu1L37vfy0ekd8zt8BTbZGW45qTBVa1frSlKOePe7Slqnd8y/NAcu3R12/DJM/XnX60pSBenq7nevqaukjR05tOOC3nAMjP8WzL6wMOHLpjt1vq4k5UTuPn5XBdnjBFhnY5h+FCxtyTqNJGXOUlf5Wm0dmHw2vPww3Dct6zSSlDlLXeVt9H4s2mRPmm89nYcefXSFlzp6bt1n2SVVMq+pq6w1PPs6xzd+iulVf+WF336PJYf8lrEjh3b43Drgs+ySKlruztR9Tr2y1M1byLyW4VzQ8mkmxb28UH/De8vbP7fus+ySKl3uSt3n1CtL2/Ppv2rdh3lpYyY+cy4seafD59Z9ll1SpfM5dZW9tufTJ6w+l21uPhB2Pw4+flKHz637LLukctfVc+qWuirLHw+HOX+Eb8+G9bfOOo0kFZ1Tryo/Jp4Bg9eEGUdDTv+HVVJ+WeqqLGtvAHudAk/PgoevzTqNJA0oS12VZ+zBsPGOcPP34Z3Xs04jSQPGUlfZ6mwgmYbn3uC3Gx5NevtVuP2MjNJJ0sBz8BmVpY4Gl1lx0JkqFg+ayJfvv4zY/gAYMTbryJLU73J3pu7gM5Whs4Fkll9+3pLP8fbgYTD9aGhdmnFiSep/uSt1B5+pDJ0NJLP88sU1a/Pyrv8JLz4A91+ecWJJ6n8+p66y1dlAMiss32wI/O+n4fkGOPJ+WGfDDBNL0qpz8JkOWOo5svApuHhXGD0VPucZu6Ty5uAzyrdhW8JHj4Y518JTd2SdRpL6jaWufPjI92C9UTDjGFjybtZpJKlfWOrKh0Grwz7nwWtPwd0/zzqNJPULS135seXHYbvPwF3nFa6zS1KFsdSVL3ufCdWD4cb/cMIXSRXHUle+rLsRfPwH8NRt8Oifs04jSUXlMLEqO509n97d+kPXHEzT282M3/yzjN3oNzDzBP4+eCz3zG/u0bZW2s6yAW+6ytI+a2+zS1JvWOoqK52N+d7d+ouXtJKAqoDBNVX86ZOn8f7pn+Kh/z2O85Z8udttdbSdmqqACFqWdpylfdZT9t2O06Y/0uPsktRbufv43bHfy1tnY753t37b1fO2993+5ibM2fizfCluYjTPdLutDrezNLGkiyzts86c82KvsktSb+Wu1B37vbx1NuZ7d+u3HehVy71v6R4n08S6nDHoVwyuocttdbid6mBQF1naZ508ZqNeZZek3nKYWJWdVb6mvtz7nr79V2wx6ygad/0RI/c+stfbAa+pSxpYjv3eAUtdQOGxtl9PhZcegiPrYe0Nsk4kSV1y7HepMxGwz/nQ/DbccnLWaSRplVjq0vBt4CPfhYeugafvyjqNJPWZpS4B7H4sDBlZmPClpTnrNJLUJ5a6BDBoDZjyE3h1Lsy+IOs0ktQnlrrUZpuJMHo/+MuPoemZrNNIUq9Z6tLyJp0NVTVw43FO+CKp7Fjq0vJqR8AeJ8ITN8Pj07NOI0m9YqlL7e3yTXjfGJh5PCx+K+s0ktRjlrrUXnVN4dn1N56Hv5yddRpJ6jFLXerIZrvAjl+F2RfDy49knUaSesRSlzrziR/CGkNg+lHQ2pp1GknqlvOpK5c6mlhlpWVrrgcTTofrvs3t15xP7UcO7XYSlq4mj+lu/5K0qix15U5DYxMHXVZHc0srg2uquOqw8QArLRs7cigNQyfTmt7PDnN/yuTHNuWiwyZ2WsJt2128pDDvelWwwra62r/FLqkYcvfxe0RMjYhpixYtyjqKMlI3byHNLa20JljS0krdvIUdLgOoe/o1frDkENbmHY7mN+8t72q7bU+3t99WV/uXpGLIXamnlG5IKR1eW1ubdRRlZPyoYQyuqaI6YFBNFeNHDetwWdu6jdWbcfnSKXyh+k72Wmtet9tt+4+qqt22utq/JBWD86krl3p0TX25deufmM/XHtif1dashW/MgupBXW7Xa+qS+ktX86lb6lJPPX4jXHMATDitMFWrJGWgq1LP3cfvUp+9fwpsOwXuPBtefy7rNJK0Ektd6o3J5xR+v+mEbHNIUgcsdak3hmwGHzuuMNnL3JlZp5GkFVjqUm/teiQMH12YnrX5n1mnkaT3WOpSb1UPgn3Ph0XPwqxzs04jSe+x1KW+GLkbbH8Q3HMBLHg86zSSBFjqUt9NOA0Grw0zjoGcPhoqqbRY6lJfrbU+TDgVGv8KD16TdRpJstSlVbLDV2CTneGWk+Dt17JOIynnLHVpVVRVwb4/hXdeh9tOzTqNpJyz1KVVteEYGP8taLgCnrs/6zSScsxSl4phjxNgnY1h+lGwtCXrNJJyylKXimG1dWDy2fDyw3DftKzTSMopS10qltH7wdYT4Y4fwRsvZJ1GUg5Z6lKxRMDkH0NrixO+SMqEpS4V03pbwO7HwqPXwRO3Zp1GUs5Y6lKx7fYdGLY13HgMLHkn6zSScsRSl4qtZjXY5zxoegbuOj/rNJJyxFKX+sOoj8GH9oe7fwavPpF1Gkk5YalL/WXiGTBoDZhxtBO+SBoQlrrUX9beAPY6BZ6eBQ9fm3UaSTlgqUv9aezBsPGOcPP3C+PDS1I/qohSj4jREfHLiLg2Ir6VdR7pPVXVhQlf3n4Vbj8j6zSSKlzmpR4Rv4qIBRExp93ySRExNyKejIguR/JIKT2WUvom8AVgXH/mlXpt4+1h58Ph/svg+b9lnUZSBavJOgBwBXAh8D9tCyKiGrgImADMB+6PiOuBauCsdu8/JKW0ICL2A05Yti2pzxoam/jD3+YTwHYb19L0djPjRw1j7Mih3a7/mR036Xi9PU+CR/5cmPDl67dDVfV773v1zcUMX2e19/Y1dM3BK+2zobGJunkLO83R9vrQNQcz54VFK2UHVnh/V9tr//13tL3Ofhad+c29zzJzzotMHrMRB+6yWac/x66+x1LVUe7llwFl+X11JWs8KoIAAAkNSURBVKu/q3I8RgY6c+alnlKaFRGbt1u8M/BkSmkeQERcA3wypXQWsG8n27keuD4iZgC/6b/EqmQNjU0cMG02zUv/dbd6VcDgmiquOmx8hwW4/Pq/b5jP1V9feT1WXxcmnQnXHgL3X07Dhp9faT8AAaR2+wQ46LI6mltaO8zR0NjEQZfVsXhJK+3vsa8KqKkKiKBlaeH9p+y7HadNf6TD7XX0/bffXmc/i8785t5n+f6fHgbgrideBVip2Nu+h86+x1LVUW74199X+599uXxfXcnq76ocj5EsMmf+8XsnRgDPLff1/GXLOhQRe0TELyLiUuDGLtY7PCLqI6L+lVdeKV5aVYy6eQtZ0q7QWhMsaWmlbt7CbtfvbD0AtvsMjNoTbj+dhx57fKX9AO+V8vL7rJu3kOaW1k5ztL3eUQ23JliyNLFkuffPnPNip9vr6Pvv6c+iMzPnvNjl18t/D33ZfpY6yr3CsnY/+3L5vrqS1d9VOR4jWWQu1VKPDpZ1+i9NSunOlNJ3UkrfSCld1MV601JK41JK44YPH16UoKos40cNY1D1iodfVcCgmqr3Pkrtav3O1gMKE77scx60LOaTCy5eaT/wrwN/+X2OHzWMwTVVVHeSo+31jv5jrgoYVB0MWu79k8ds1On2Ovr+e/qz6MzkMRt1+fXy30Nn32Op6ij3Csva/ezL5fvqSlZ/V+V4jGSROVIJDIqx7OP36SmlMcu+3hX4YUpp72Vfnwiw7OP3ohg3blyqr68v1uZUQfrlmvry7jgL/nI2/9j7Sq54aXOvqbf7Hsrpeil4Td1r6l3rj8wR0ZBS6vCm8FIt9RrgH8BewPPA/cCBKaVHirVPS12ZWfIuXLIrEPCte2DQ6lknklRGuir1zD9+j4irgdnAthExPyIOTSm1AEcCNwOPAb8rVqFHxNSImLZo0aJibE7qvUGrFz6Gf+0puPvnWaeRVEFK4kw9C56pK3O/PxgenwHfng3Dtsw6jaQyUdJn6lJu7X0mVA+GG//DCV8kFYWlLmVl3Y3g4z+Ap26DR/+cdRpJFcBSl7K002Gw0YfhphPh3TeyTiOpzOWu1L1RTiWluqYw4cubL8EdZ2adRlKZy12pp5RuSCkdXltbm3UUqWDEWBh3CNx3Kbz4YNZpJJWx3JW6VJL2OgXWHAbTj4bW1qzTSCpTlrpUCtYYAhN/BM/Xw9+uyDqNpDJlqUul4kNfgM0/Crf+EN5ywiFJvWepS6UiAvY5H5rfhlt+kHUaSWUod6Xu3e8qacO3gY98Fx66Bp6+K+s0kspM7krdu99V8nY/FoaMhBnHQEtz1mkklZHclbpU8gatAVN+Aq/OhdkXZJ1GUhmx1KVStM1EGL0f/OVcaHom6zSSyoSlLpWqSWdDVTXceJwTvkjqEUtdKlW1I2CPE+GJm+Hx6VmnkVQGclfq3v2usrLLN+F9Y2Dm8bD4razTSCpxuSt1735XWamuKTy7/sbz8Jezs04jqcTlrtSlsrPZLrDjV2H2xfDyI1mnkVTCLHWpHHzih4Xx4acf5YQvkjplqUvlYM31YMLp8Ny98MCVWaeRVKIsdalcbH8gbLYb/N8p8M+FWaeRVIIsdalcRMC+58PiN+HWU7JOI6kEWepSOdlgNOx6JPz9SmicnXUaSSUmUs5GqoqIqcBUYH/giW5WrwV68kB7T9brbp2uXl8feLUHOUpJT392pbKfVdlOb99brOPKY6r099PXbWV1THW3jsdU9vupBYaklIZ3+GpKyV+d/AKmFWu97tbp6nWgPuufRX/97EplP6uynd6+t1jHlcdU6e+nr9vK6pjqbh2Pqez30922/Pi9azcUcb3u1unpvsrFQH0/xdrPqmynt+8t1nHlMVX6++nrtrI6prpbx2Mq+/10ua3cffxejiKiPqU0LuscqhweUyo2j6nS4Jl6eZiWdQBVHI8pFZvHVAnwTF2SpArhmbokSRXCUpckqUJY6pIkVQhLvcxExKiIuDwirs06iypHRHwqIv4rIq6LiIlZ51H5i4jREfHLiLg2Ir6VdZ68sNRLQET8KiIWRMScdssnRcTciHgyIk4ASCnNSykdmk1SlZNeHld/Til9HfgahdEWpZX08ph6LKX0TeALgI+6DRBLvTRcAUxafkFEVAMXAZOBDwAHRMQHBj6aytgV9P64+sGy16WOXEEvjqmI2A/4K3DbwMbML0u9BKSUZgGvtVu8M/DksjPzZuAa4JMDHk5lqzfHVRScA8xMKf1toLOqPPT236qU0vUppd2AgwY2aX5Z6qVrBPDccl/PB0ZExLCI+CWwQ0ScmE00lbEOjyvg34FPAJ+LiG9mEUxlq7N/q/aIiF9ExKXAjdlEy5+arAOoU9HBspRSWgj4j676qrPj6hfALwY6jCpCZ8fUncCdAxtFnqmXrvnApst9vQnwQkZZVDk8rlRsHlMlxFIvXfcDW0fEFhExGPgicH3GmVT+PK5UbB5TJcRSLwERcTUwG9g2IuZHxKEppRbgSOBm4DHgdymlR7LMqfLicaVi85gqfU7oIklShfBMXZKkCmGpS5JUISx1SZIqhKUuSVKFsNQlSaoQlrokSRXCUpckqUJY6pIkVQhLXZKkCmGpS5JUISx1SaskIv4cESki/r2D105f9tplWWST8sax3yWtkohYD/g78D5g15TS35ct3wu4BXgc2Cml9HZ2KaV8sNQlrbKI2A34C/A0sCOwJvAgUEuh0J21SxoAfvwuaZWllO4BTga2Bi4FrgQ2BL5joUsDxzN1SUUREQHMBPZetujqlNKBGUaScsczdUlFkQpnCH9abtHPssoi5ZVn6pKKIiK2Bv4GLKFwLf0RYOeU0ruZBpNyxDN1SassIlYDfgusBXwROAv4IJ6tSwPKUpdUDD8BdgB+nFK6BfhP4G7gGxHxhUyTSTnix++SVklEfIrCtfR7gX9LKbUsW74p8ABQA+yQUpqXXUopHyx1SX0WEZtRKO4qCsX9dLvXPwn8GbifQuE3D3xKKT8sdUmSKoTX1CVJqhCWuiRJFcJSlySpQljqkiRVCEtdkqQKYalLklQhLHVJkiqEpS5JUoWw1CVJqhCWuiRJFeL/AUIokz4zVfOJAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "fig = plt.figure()\n", "ax = fig.add_axes([0, 0, 1, 1])\n", "\n", "ax.plot( x_bins[:-1], x_hist, \".\")\n", "ax.plot( x_bins[:-1], px_theory, \"-\")\n", "\n", "ax.set_xscale('log')\n", "ax.set_yscale('log')\n", "\n", "ax.set_xlabel(\"x\", fontsize=20)\n", "ax.set_ylabel(\"p(x)\", fontsize=20)\n", "\n", "ax.set_ylim(bottom = PlotYlim_bottom)\n", "\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "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.7.7" } }, "nbformat": 4, "nbformat_minor": 4 }