{ "cells": [ { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# 7: Random number generation\n", "\n", "\n", " #### [Back to main page](https://petrosyan.page/fall2020math3215)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Sampling from the logistic function" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The logistic function has the following form\n", "\n", "$$F(X)=\\frac{e^x}{1+e^x}\\quad \\Rightarrow \\quad f(x)=\\frac{e^x}{(1+e^x)^2}.$$\n" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6IAAAE/CAYAAABcn34zAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdd3hc9Z32//sz6s2SrGZbzbbcexEGTO8mQIAECJBGGiEbQrJ50n6bbLK7efI8yY/NbpINCSGEVBJaqInBOKEYMOAqyw0LuarYsmRb1aoz3+cPSY4QMpbt0ZyZ0ft1XVzSnDk6c8NlM7rnnPP5mnNOAAAAAACEis/rAAAAAACA0YUiCgAAAAAIKYooAAAAACCkKKIAAAAAgJCiiAIAAAAAQooiCgAAAAAIKYooECHMzJnZlL7vk8zsGTNrMrNHvc4GAMBoZGYT+96fY/se55nZKjNrMbMfep0PCGexXgcAcEpukJQnKcs51+N1GAAAIEm6XVKDpDHOOed1GCCccUYUiEzFkioooQAAhJViSdsoocCJUUQBD5nZHjP7/8xsm5kdMbNfm1li33NfNbP9ZlZrZp8c8DP/Lunbkj5kZq1m9imv8gMAEK3MrNDMHjezejM7ZGY/NbMYM/tPM2sws12Srhqw/28kfVzS1/reny/1KjsQCbg0F/DehyVdIalN0jOSvmVmr0r6iqRLJO2W9Mv+nZ1z3zEzJ2mKc+4jHuQFACCqmVmMpL9IekHSRyX5JZVK+oykqyUtVO/79p/7f8Y5d5uZSVK1c+5boc4MRBrOiALe+6lzrso5d1jS9yTdIukmSb92zm1xzrVJ+jcvAwIAMMoskTRB0ledc23OuQ7n3KvqfX/+0YD37f/raUogglFEAe9VDfh+r3rf+CYMsR0AAIRGoaS9Q8xi4P0ZCBKKKOC9wgHfF0mqlbR/iO0AACA0qiQV9S/LMgDvz0CQUEQB733ezArMbKykf5H0sKRHJN1mZrPMLFnSdzxNCADA6LJGvaXz+2aWYmaJZnaOet+f7+p7386U9A1PUwIRjCIKeO+Pkp6XtKvvn//tnHtW0o/UOyShsu8rAAAIAeecX9I1kqZI2iepWtKH1Ds8cIWkTZI2SHrcq4xApDOWOQK8Y2Z7JH3aOfc3r7MAAAAAocIZUQAAAABASA2riJrZMjPbYWaVZjbktfBmdqGZlZnZVjN7ObgxAQAAAADR4oSX5vYt6Fsh6TL1Xh+/VtItzrltA/bJkLRa0jLn3D4zy3XOHRy52AAAAACASDWcM6JLJFU653Y557okPSTp2kH73CrpcefcPkmihAIAAAAAjmc4RTRf71y4t7pv20DTJGWa2Utmtt7MPhasgAAAAACA6DJ4kd6h2BDbBl/PGytpsaRLJCVJet3M3nDOVbzjQGa3S7pdkmbNmrV469atJ58YAIChDfV+hZOwbNky99xzz3kdAwAQPY773jycM6LVkgoHPC6QVDvEPs8559qccw2SVkmaP/hAzrn7nHOlzrnSpKSkYbw0AAAIlYaGBq8jAABGieEU0bWSpprZJDOLl3SzpKcH7fOUpPPMLNbMkiWdKWl7cKMCAAAAAKLBCS/Ndc71mNmdklZIipH0gHNuq5nd0ff8vc657Wb2nKRySQFJ9zvntoxkcAAAAABAZBrOPaJyzi2XtHzQtnsHPb5b0t3BiwYAAAAAiEbDuTQXAAAAAICgoYgCAAAAAEKKIgoAAAAACCmKKAAAAAAgpCiiAAAAAICQoogCAAAAAEKKIgoAAAAACCmKKAAAAAAgpCiiAAAAAICQoogCAAAAAEKKIgoAAAAACCmKKAAAAAAgpCiiAAAAAICQoogCAAAAAEKKIgoAQIQxswfM7KCZbTnO82ZmPzGzSjMrN7NFoc4IAMB7oYgCABB5fiNp2Xs8f6WkqX3/3C7p5yHIBADAsMV6HQAAAJwc59wqM5v4HrtcK+l3zjkn6Q0zyzCz8c65/SEJCAAYUc45Odf3ff9jacC2Ac8P2ta//z9+tu+bQfsMfr13ZXjXPu/OmZOWcNx/B4ooAADRJ19S1YDH1X3bKKIARp2unoCOdvXoaJdf7d1+tfd9PdrV/32P2rt69+nyB9Td49TtD6jbH+h9PGDbscf+vsc9AfUEnPwBp4Dr/yoFAk5+17vt2PcBDdjnePsOKpF6d8mMJHu+f9Vxn6OIAgAQfWyIbUP+CmNmt6v38l0VFRWNZCYAOG1tnT2qb+lUQ2vnP762dqnpaJeaO3rU1N6t5vbu3q8dvV87ugMn/ToxPlNcjCkuxqf4GJ/iYnyKix30uO/5xDiffGaK8ZlizOTzmXzWewyf2bHner+qd1v/vqZj38f4TGYms97/iVvf/8lN/9jWv7H/eev73/3AnzH7x1tA/z7/ONa7t/XuZ8feOP5xrHe+ldgQ7yzv2jTUTsdBEQUAIPpUSyoc8LhAUu1QOzrn7pN0nySVlpZG4OftAKKFc071LZ2qOnJUVYfbVX3kqKqPtKuq7+vB5k61d/vf9XNm0pjEOI1JitWYxDilJ8WpJCdV6Un/2JaSEKuk+Bglx8coMa73a1JcTN+22N7v42KUENdbMmN8wy9UODUUUQAAos/Tku40s4cknSmpiftDAYSTw21d2lLTpLcPturtupZjX5s7et6xX3ZqggoykzQnP13jZyYqOy1BOakJA77Ga2xyvGJjmMEaaSiiAABEGDP7k6QLJWWbWbWk70iKkyTn3L2Slkt6n6RKSUclfcKbpAAgdfb4VV7dpE1VjSqratSm6kZVHW4/9vzYlHhNzU3V+xdM0JScVBVnp6gwM0n5GclKio/xMDlGEkUUAIAI45y75QTPO0mfD1EcAHiHQMBp2/5mvVbZoFcrG7R2z+Fj92nmZyRpfmG6PnJmseYWpGt6XpqyUo8/WRXRiyIKAAAA4LR09QT0xq5Dem7rAT2/tU4NrZ2SpKm5qbr5jCKdXZKlhUUZyk1L9DgpwgVFFAAAAMBJc85pze7Demx9tVZsPaDmjh4lx8foohm5umRGrs6Zkq28MRRPDI0iCgAAAGDYDrZ06NF11Xp0XZX2HDqq1IRYXTF7nK6cM07nTs1WYhz3deLEKKIAAAAATujtuhbd/8puPbGxRl3+gM6aPFZ3XTJVy+aMU3I8tQInhz8xAAAAAI5rc3WT/vtvFXrhrYNKiPXppjMK9MlzJmlyTqrX0RDBKKIAAAAA3qXyYKt++PwOPbvlgDKS4/TPl07TR84qYsotgoIiCgAAAOCYpvZu/dfzO/T7N/YqKS5Gd10yVZ85b5LSEuO8joYoQhEFAAAAIOecniyr0ff++pYOt3XqI2cV64uXTOUMKEYERRQAAAAY5fY3tetrj5XrlbcbNL8wQ7/5xBmak5/udSxEsWEVUTNbJunHkmIk3e+c+/6g5y+U9JSk3X2bHnfO/UcQcwIAAAAYAU+V1ehfn9yibr/Td6+drQ+fWSyfz7yOhSh3wiJqZjGS7pF0maRqSWvN7Gnn3LZBu77inLt6BDICAAAACLKObr++9eQWPba+WguLMvTfNy3QxOwUr2NhlBjOGdElkiqdc7skycweknStpMFFFAAAAEAE2HfoqO74w3pt29+suy6eorsumarYGJ/XsTCKDKeI5kuqGvC4WtKZQ+x3tpltklQr6SvOua1ByAeMWvub2vVUWa3au/xaNmecZo4f43UkAAAQBVZXNuhzD26Qc06/vu0MXTQj1+tIGIWGU0SHukDcDXq8QVKxc67VzN4n6UlJU991ILPbJd0uSUVFRScZFRg9Xq6o1+cf3KDWzh6ZpJ/8/W1986qZ+vR5k72OBgAAItgTG6v1tcfKNSk7Rfd/7AwVZSV7HQmj1HDOv1dLKhzwuEC9Zz2Pcc41O+da+75fLinOzLIHH8g5d59zrtQ5V5qTk3MasYHotX1/sz77+3XKTI7T966box99aIEWFWfqf/91u54qq/E6HgAAiEDOOd3zYqX++eFNKi0eq0fvWEoJhaeGU0TXSppqZpPMLF7SzZKeHriDmY0zM+v7fknfcQ8FOywQ7Xr8Af3zw2VKiI3Rly6dprwxiUpJiNVnzpukqbmp+pfHN+tgc4fXMQEAQARxzukHz+3Q3St26LoFE/TbTy5RelKc17Ewyp2wiDrneiTdKWmFpO2SHnHObTWzO8zsjr7dbpC0pe8e0Z9Iutk5N/jyXQAn8PjGGr11oEW3LCl8xxtErM+n25ZOVGdPQHev2OFhQgAAEEmcc/o/y7fr3pd36qNnFeu/blqg+FiGEsF7w1pHtO9y2+WDtt074PufSvppcKMBo0u3P6Af/+1tFWcla3FR5ruezxuTqAun5+jxjTX60mXTlJ+R5EFKAAAQKZxz+u5ftuuB13brtqUT9Z1rZqnvIkbAc3wcAoSJZzbVqqaxXdfOn3DcN4nLZ42TJN3/yq5QRgMAABHox39/Ww+8tlufPGcSJRRhhyIKhIkH39ynvDEJmpufftx9xqbEa8nEsXpoTZVaO3tCmA4AAESSP7yxVz/629u6cXGB/vXqmZRQhB2KKBAGKupatH7vEZ0/NeeEbxQXTMtRe7dfz27eH6J0AAAgkjy7eb/+9aktumRGrv7vB+ZSQhGWKKJAGHhkbZVifaalJVkn3LckJ0V5YxL02PrqECQDAACRpLy6UV96uEwLCzP001sXKTaGX/cRnviTCXjMOaflm/dr1vgxSks88Sh1M9PSkmy9ufuwqg4fDUFCAAAQCQ42d+j2361XdmqCfvmxUiXFx3gdCTguiijgsfLqJtU2dWhx8bsn5R7PGRN7912x9cBIxQIAABGks8evO/6wXk3t3frlx0qVlZrgdSTgPVFEAY89u+WAYnym+YUZw/6Z3LREFWYmUUQBAIAk6d+e3qoN+xr1w5vma9aEMV7HAU6IIgp4qP+y3Bl5aUpNGNayvscsLMrUuj1HVN/SOULpAABAJHh6U63+tKZKn7uwRO+bO97rOMCwUEQBD1XUtWrf4aNadBKX5fZbWJQhJ+nv2+uCHwwAAESEvYfa9C+Pb9bi4kx9+bJpXscBho0iCnjo5YqDkvSea4ceT0FGkjKT47Tq7fpgxwIAABGgqyegL/xpo3wm/fjmBYpjQi4iCH9aAQ+tqmjQhIxEjU2JP+mfNTPNGj9Gr77dIH/AjUA6AAAQzv5rZYXKq5v0/98wXwWZyV7HAU4KRRTwSHuXX2t2H9bs8Sd/NrTf7Anpau7o0abqxiAmAwAA4W7DviO6b9VO3XxGoZbNGed1HOCkUUQBj6zZc1hd/oBmn8Zku1njx8gkvVLRELxgAAAgrHV0+/XVRzdpfHqSvnnVTK/jAKeEIgp4ZFVFveJiTFPzUk/5GKmJsZqYnayXK7hPFACA0eK/V1ZoZ32bvv/BuUpLjPM6DnBKKKKAR159u0FTclOVEBtzWseZOX6MNlU1qrWzJ0jJAABAuNq474h++cou3bKkSOdNzfE6DnDKKKKABxqPdqmirkXT89JO+1jT89Lkd04b9h4JQjIAABCuevwBffOJLcpNS9S/vG+G13GA00IRBTywbs8ROUlTc0+/iJbkpMpn0prdh08/GAAACFu/f2Ovtu1v1revmcUluYh4FFHAA2v2HFaszzQpO+W0j5UYF6PirGS9uftQEJIBAIBwdLC5Qz98vkLnT8vRlUzJRRSgiAIeeHP3IU3KTlF8bHD+Ck7NTdOmqiZ1dPuDcjwA4c/MlpnZDjOrNLNvDPF8upk9Y2abzGyrmX3Ci5wAguN7y7eryx/Qf7x/tszM6zjAaaOIAiF2tKtHW2qaNTX31KflDjYtL01d/oDKq5uCdkwA4cvMYiTdI+lKSbMk3WJmswbt9nlJ25xz8yVdKOmHZhYf0qAAguKNXYf0VFmtPndBiSYG4WoqIBxQRIEQ27ivUf6A07QgDCrqNyWnt9Su3cN9osAosURSpXNul3OuS9JDkq4dtI+TlGa9p05SJR2WxHhtIMIEAk7f++t2TUhP1OcuLPE6DhA0FFEgxN7cfVg+6x0yFCypibHKz0jSm7u4TxQYJfIlVQ14XN23baCfSpopqVbSZklfdM4FQhMPQLA8talGm2ua9NVl05UYd3pLvgHhhCIKhNja3YdVODZZSfHBfTOZlpeqdXuPqMfP75nAKDDUDWJu0OMrJJVJmiBpgaSfmtmYdx3I7HYzW2dm6+rr64OfFMAp6+j26+7ndmheQbqunT/4syYgslFEgRAKBJw2VTdq8gjc31GSk6qjXX5V1LUG/dgAwk61pMIBjwvUe+ZzoE9Ietz1qpS0W9K7Fh50zt3nnCt1zpXm5OSMWGAAJ+9Xr+5WbVOHvvm+mfL5GFCE6EIRBUJoV0Objnb5NTEr+EV0ck7vMcuqGoN+bABhZ62kqWY2qW8A0c2Snh60zz5Jl0iSmeVJmi5pV0hTAjhlh9u69POXduryWXk6c3KW13GAoKOIAiG0uaa3JI5EEc1JTVBaYqzKqo4E/dgAwotzrkfSnZJWSNou6RHn3FYzu8PM7ujb7buSlprZZkl/l/R151yDN4kBnKxfvLxTR7t69LVl77qQAYgKsV4HAEaT8uomJcT6ND49MejHNjNNzErRhn2cEQVGA+fccknLB227d8D3tZIuD3UuAKfvYEuHfvv6Hl23IF9TgrjcGxBOOCMKhFBZVaOKxiaP2H0ek3NStLO+VS0d3SNyfAAAMPJ+9uJOdfud7rpkqtdRgBFDEQVCpMcf0LbaZhVnJY/Ya0zOTpFzvWdeAQBA5Nnf1K4/vrlPNywq0MQRGG4IhAuKKBAilfWt6uwJjMj9of0m9b1hbdzHfaIAAESin75QKSenL1wyxesowIiiiAIh0n+WciSLaHJ8rCZlpzA5FwCACFR1+KgeWVelm88oUkHmyF1BBYQDiigQIuXVjUqM8yl3TMKIvs6CwnSVVTXKucFr2wMAgHB236pdMpk+fxFnQxH9KKJAiGyqalLx2BT5bGQXpJ5XkKGG1i5VH2kf0dcBAADBU9/SqUfWVemDi/M1bgSm6wPhZlhF1MyWmdkOM6s0s2+8x35nmJnfzG4IXkQg8nX1BLTjQIsmjuCgon7zC9MlSRu5PBcAgIjxwGu71e0P6PbzS7yOAoTECYuomcVIukfSlZJmSbrFzGYdZ78fqHdxbQADVNS1qMsfUPEI3h/ab1pemhJifdpEEQUAICI0d3TrD6/v1ZVzxx8bPAhEu+GcEV0iqdI5t8s51yXpIUnXDrHfFyT9WdLBIOYDosLmmr5BRdkjf0Y0LsanmePHHHtNAAAQ3v7wxl61dPbocxdwNhSjx3CKaL6kqgGPq/u2HWNm+ZKul3Tvex3IzG43s3Vmtq6+vv5kswIRq7yqUSnxMcpJHdlBRf3mFaRra02TAgEGFgEAEM46uv164NXdOn9ajubkp3sdBwiZ4RTRoSarDP7t9keSvu6c87/XgZxz9znnSp1zpTk5OcPNCES8supGFWelyEZ4UFG/Ofnpauvya1dDW0heDwAAnJpH11erobWLs6EYdYZTRKslFQ54XCCpdtA+pZIeMrM9km6Q9DMzuy4oCYEI19HtV0Vda0gGFfWbV9D7ieoWLs8FACBsBQJOD7y6W/ML0nXW5LFexwFCajhFdK2kqWY2ycziJd0s6emBOzjnJjnnJjrnJkp6TNI/OeeeDHpaIAK9daBF/oALyaCiflNyUpUY51N5NUUUAIBw9VLFQe1uaNMnz50UsqumgHARe6IdnHM9Znaneqfhxkh6wDm31czu6Hv+Pe8LBUa7zdW902tDeUY0tm9gEWdEAQAIX79+bY/yxiTofXPHex0FCLkTFlFJcs4tl7R80LYhC6hz7rbTjwVEj03VTUpLjNXYlPiQvu68/HQ9tr5a/oBTjI9PWQEACCcVdS165e0GffWK6YqLGc5FikB04U89MMI2VTWqOCs55Jfc9A8s2t3QGtLXBQAAJ/br1/YoIdanW5YUeR0F8ARFFBhB7V1+7axv1cQQ3h/ab15BhiSxnigAAGHmSFuXHt9QrQ8syg/5FVNAuKCIAiNo2/4mBZw8KaIlOSkMLAIAIAz9cc0+dfYEdNvSSV5HATxDEQVGUH8JDOWgon6xMT7NnpDOwCIAAMJItz+g37++V+dOydb0cWlexwE8QxEFRlB5dZPSk+KUkezNZTdz89O1paZZ/oDz5PUBAMA7Pb+1TgeaO/SJcyZ6HQXwFEUUGEGbqho9ORvab05+utq7/dpVz8AiAADCwYNv7lV+RpIunJ7rdRTAUxRRYIS0dvZod0ObJ/eH9ptXkC6JgUUAAISDnfWtWr3zkG49s4il1TDqUUSBEbK1pklO0sRs74poSU6qkuJiGFgEAEAY+NOb+xTrM91YWuB1FMBzFFFghPSXv+Kx3l2aG+MzzZ4whoFFAAB4rKPbr0fXV+uK2eOUm5bodRzAcxRRYIRsqm5UVkq8xiTFeZpjTn66ttYysAgAAC/9tXy/mtq79eGziryOAoQFiigwQsqrm1Ts4aCifnP7BhbtZGARAACeefDNvZqcnaKzJ2d5HQUICxRRYAQ0He3WvsNHPR1U1O/YwCLuEwUAwBPbapu1YV+jbj2zSGYMKQIkiigwIrbU9t0fGgZnRCfnpCo5PobJuQAAeOTBN/cqPtanGxYzpAjoRxEFRsCxQUVhcEa0f2ARRRQAgNA72tWjp8pqdfXc8cpIjvc6DhA2KKLACCivblROWoJSE2K9jiKpf2BRk3r8Aa+jAAAwqjy7+YBaO3v0oTMKvY4ChBWKKDACyqsbPV22ZbB5Benq6A5oZ32b11EAABhVfr6iTDmJTksmjfU6ChBWKKJAkB1u61JNY0dYDCrqNze/b2ARl+cCABAyuxvaVNlsOivPMaQIGIQiCgRZf9mbmB0+Z0QnZacqJT5Gm6sbvY4CAMCo8dj6Kpmk8woSvI4ChJ3wuIENiCL9Za94bPicEe0dWJTOGVEAAELEH3B6bH21Lpyeo5uuWeJ1HCDscEYUCLLy6iaNG5OopPgYr6O8w5z8dG3b38zAIgAAQmDV2/Wqa+7UTaUMKQKGQhEFgmxTdWNYrB86WP/Aosr6Vq+jAAAQ9R5dV6WxKfFKb9unNWvWeB0HCDsUUSCIDjZ3qK65M6wGFfWb0z+wqJrLc4FoYGbLzGyHmVWa2TeOs8+FZlZmZlvN7OVQZwRGq0OtnVq5rU7XL8zXofqDqqur8zoSEHYookAQHRtUFIZnRCdnp/QOLOI+USDimVmMpHskXSlplqRbzGzWoH0yJP1M0vudc7Ml3RjyoMAo9WRZrbr9jstygfdAEQWCqLy6ST6TCsNoDdF+Pp9pdj4Di4AosURSpXNul3OuS9JDkq4dtM+tkh53zu2TJOfcwRBnBEYl55weXVel+QXpmj4uzes4QNiiiAJBtLm6UePTk5QYF16DivrNy0/XtloGFgFRIF9S1YDH1X3bBpomKdPMXjKz9Wb2sZClA0axzTVNeutAi27kbCjwniiiQJA457SppiksBxX1m1uQrs6egN4+yMAiIMLZENvcoMexkhZLukrSFZL+1cymvetAZreb2TozW1dfXx/8pMAo8/iGGsXH+nTN/AmSpJSUFKWkhN/sCMBrrCMKBMmB5g4dau3SxFnh+2Yzt39gUU2TZo4f43EaAKehWtLA0y0FkmqH2KfBOdcmqc3MVkmaL6li4E7Oufsk3SdJpaWlg8ssgJPQ7Q/omU21umxmntKT4iRJF198scepgPDEGVEgSMr7ptGG8xnRiVkpSk2IZXIuEPnWSppqZpPMLF7SzZKeHrTPU5LOM7NYM0uWdKak7SHOCYwqqyrqdaitS9cvHHylPIDBOCMKBMnm6ibF+EyFmeFbRH0+0+wJYxhYBEQ451yPmd0paYWkGEkPOOe2mtkdfc/f65zbbmbPSSqXFJB0v3Nui3epgej3+MYaZSbH6fxpOce2rV69WpK0dOlSr2IBYYkiCgRJeXWj8jOSFB8b3hcazCtI129f36tuf0BxMeGdFcDxOeeWS1o+aNu9gx7fLenuUOYCRqvmjm6t3Fanm88ofMfvAocOHfIwFRC++C0UCALnnMprmlQUhsu2DDYnP11dPQG9XcfAIgAAguXZzfvV1RPgslxgmCiiQBBUHW5X49FuTcoO30FF/eYVZEiStnB5LgAAQfPExhpNyk7RgsIMr6MAEWFYRdTMlpnZDjOrNLNvDPH8tWZWbmZlfSPgzw1+VCB8bapulCRNDONBRf2KxyYrLSFW5TWNXkcBACAq1DS2641dh3X9wnyZDbW6EoDBTniPqJnFSLpH0mXqHQW/1syeds5tG7Db3yU97ZxzZjZP0iOSZoxEYCAclVc3Ki7GlJ+Z5HWUE/L5THPy07W5ptnrKAAARIUnN9ZIkq5b8O7LctPT00MdB4gIwxlWtERSpXNulySZ2UOSrpV0rIg65wbebJaidy+qDUS1TdW994fG+iLjave5Ben6zeo9DCwCAOA0Oef0xMYalRZnqmiIK6POP/98D1IB4W84v4HmS6oa8Li6b9s7mNn1ZvaWpL9K+mRw4gHhzx9w2lzdpOKs8L8/tN/cvoFFFXUtXkcBACCibalpVuXBVl2/iCFFwMkYThEd6kL3d53xdM494ZybIek6Sd8d8kBmt/fdQ7quvr7+5JICYaryYKvau/2aFGFFVGJgEQAAp+uJjTWKj/Hp6rkThnx+1apVWrVqVYhTAeFvOEW0WlLhgMcFkmqPt7NzbpWkEjPLHuK5+5xzpc650pycnCF+Gog8xwYVZYf/oKJ+xVnJSkuMVXk1RRQAgFPV4w/o6U21unhGrtKT44bcp6mpSU1NvN8Cgw2niK6VNNXMJplZvKSbJT09cAczm2J9I8LMbJGkeEms3otRoby6UUlxMcobk+h1lGEzM83NT+eMKAAAp+GVygY1tHZyWS5wCk5YRJ1zPZLulLRC0nZJjzjntprZHWZ2R99uH5S0xczK1Dth90POOQYWYVTYVNWkidnJ8kXYuPa5+enavr9FXT0Br6MAABCRnthQo4zkOF00PdfrKEDEGc7UXDnnlktaPmjbvQO+/4GkHwQ3GhD+Onv82r6/WZfOzPM6ykmbW5CuLn/vwKI5+YyWB+1Esm8AACAASURBVADgZLR29uj5bQf0wUUFio9lAj1wsvhbA5yG7ftb1BNwmpQdOYOK+s0vyJAklVU1epwEAIDI8+zm/eroDugDJ7gsNysrS1lZWSFKBUQOiihwGsr7BxUNsW5YuCvITFJ2aoI27qOIAgBwsp4sq1FxVrIWFWW+535Lly7V0qVLQ5QKiBwUUeA0bKpqUnpSnMamxHsd5aSZmRYUZmhj1RGvowAAEFH2N7Vr9c5Dum5BvizCZkQA4YIiCpyGsqojKs5Kjtg3oYVFGdpV36amo91eRwEAIGI8VVYr56TrF554Wu4LL7ygF154IQSpgMhCEQVOUWtnj3bVt2lSVuTdH9pvYVHffaLVXJ4LAMBwOOf0xIYaLSzK0MRhzIhoa2tTW1tbCJIBkYUiCpyi8upGOWlYb0Lhal5Bhsykjfu4PBcAgOHYtr9ZO+pa9IFhnA0FcHwUUeAU9Q/5mRzBRTQ1IVbT89IYWAQAwDA9ubFGcTGmq+dN8DoKENEoosAp2rD3iManJyolYVjL8YatBYUZKqtqlHPO6ygAAIQ1f8DpqbJaXTg9V5kROKgQCCcUUeAUOOe0ft+RiD4b2m9hUYaa2ru1u4H7VwAAeC+vVTboYEvnSV2Wm5eXp7y8vBFMBUSmyD6VA3hkd0ObGo92qyQ31esop21h3/pnG/c1anJO5P/7AAAwUp7YWKMxibG6aEbusH9myZIlI5gIiFycEQVOwYa+eypLoqC4leSkKjUhlvVEAQB4D22dPXpuywFdNW+8EuNivI4DRDyKKHAKNuw7ouT4GI1PT/Q6ymmL8ZnmF6arrIqBRQAAHM+KrQfU3u3X9QsLTurnVq5cqZUrV45QKiByUUSBU7Bh7xFNyk6Rz8zrKEGxsDBT2/e3qL3L73UUAADC0hMba1SQmaTS4syT+rmOjg51dHSMUCogclFEgZPU0tGtirqWqBhU1G9BYYb8AafNNU1eRwEAIOzUNXfotcoGXb8wXz5fdHwIDXiNIgqcpE1VTQq46Lg/tN+CogxJvZccAwCAd3qqrEYBJ11/EtNyAbw3iihwkjbsOyKTNDknes6IZqcmaFJ2itbtoYgCADDY4xtqNL8wg+nyQBBRRIGTtGHvEU3ISFJyfHStflRanKl1ew8rEHBeRwEAIGxs39+stw60nNTaoQPl5+crP58zqcBgFFHgJAQCThuqjkTV/aH9zpg0Vo1Hu1VZ3+p1FAAAwsYTG2sU6zNdM3/CKf38okWLtGjRoiCnAiIfRRQ4CTvrW9Xc3qOS3Oi7NGfJxLGSpLV7DnucBACA8OAPOD1VVqMLp+dqbEq813GAqEIRBU7Cm7t7S9q0vOgrosVZycpOTdDa3RRRAAAkafXOBtU1d+oDi0790tpnn31Wzz77bBBTAdEhum5yA0bYmt2HlZEcp5zUBK+jBJ2ZacmkTK1lYBEAAJKkJzbUKC0xVhfPyD3lY/T09AQxERA9OCMKDJNzTm/uPqSpuakyi841xEqLx6qmsV21je1eRwEAwFNHu3r03NYDunreeCXGxXgdB4g6FFFgmKqPtKuuuVPTctO8jjJilkziPlEAACRpxdYDOtrl1/ULC7yOAkQliigwTP33h06NwvtD+80Yl6bUhFiKKABg1Ht8Q40KMpNUWpzpdRQgKnGPKDBMa3YfUkpCjCZkJHkdZcTExvi0sChDa3dznygAYPSqa+7Qa5UN+vxFU+Tznd7tOMXFxUFKBUQXzogCw7Rm92FNzUmTL0rvD+23ZOJY7ahrUdPRbq+jAADgiafKahRw0vULT31abr958+Zp3rx5QUgFRBeKKDAMB5s7tOfQ0ai+LLdfKeuJAhHBzJaZ2Q4zqzSzb7zHfmeYmd/MbghlPiCSPb6hRvMLMzQ5J/rf9wGvUESBYVizJ/rvD+23sChD8bE+vb7rkNdRAByHmcVIukfSlZJmSbrFzGYdZ78fSFoR2oRA5Nq+v1lvHWjRB4JwNlSSnnnmGT3zzDNBORYQTSiiwDC8seuQEmJ9Khqb7HWUEZcYF6PS4kyt3kkRBcLYEkmVzrldzrkuSQ9JunaI/b4g6c+SDoYyHBDJHl1XrfgYn94/f4LXUYCoRhEFhuG1ykOalpemWN/o+CuztCRL2/c363Bbl9dRAAwtX1LVgMfVfduOMbN8SddLujeEuYCI1tUT0JNlNbp0Vq4yU+K9jgNEtdHxWzVwGvY3tWt3Q5tmjo/e9UMHO7skW1LvmWAAYWmoqWlu0OMfSfq6c87/ngcyu93M1pnZuvr6+qAFBCLRC2/V6XBbl25cXOh1FCDqUUSBE3itsreMzRw3xuMkoTOvIF0p8TFavbPB6ygAhlYtaeBvygWSagftUyrpITPbI+kGST8zs+sGH8g5d59zrtQ5V5qTkzNSeYGI8Oi6auWmJei8qdleRwGi3rCK6Ikm85nZh82svO+f1WY2P/hRAW+srmxQWmKs8jOjd/3QweJifFoyaSz3iQLha62kqWY2ycziJd0s6emBOzjnJjnnJjrnJkp6TNI/OeeeDH1UIDIcbOnQSxX1+sCiAsXGBO9cTUlJiUpKSoJ2PCBanPBv2TAn8+2WdIFzbp6k70q6L9hBAS845/RKZYOm50X/+qGDnTMlW7vq23SgqcPrKAAGcc71SLpTvdNwt0t6xDm31czuMLM7vE0HRKYnNtTIH3C6sbQgqMedNWuWZs1611BrYNSLHcY+xybzSZKZ9U/m29a/g3Nu9YD931DvJUJAxNtZ36r6lk4tmz3O6yghd3ZJliRp9c4GfWARf6WBcOOcWy5p+aBtQw4mcs7dFopMQKRyzunR9dVaXJypkiCvHdrT0yNJio0dzq/dwOgxnOsOTjiZb5BPSXr2dEIB4eLY/aGjaFBRv5njxigjOe7YfwMAAKJVWVWjKg+26sbFwf/g9dlnn9Wzz/KrMTDYcD6aGc5kvt4dzS5SbxE99zjP3y7pdkkqKioaZkTAO69VNig7NV45qQleRwk5n8+0tCRLr1U2yDknG2WXJgMARo9H11crMc6nq+aN9zoKMGoM54zocCbzyczmSbpf0rXOuSFPoTCZD5Gkxx/QG7sOaca4MaO2hF0wLUcHmju0o67F6ygAAIyI9i6/nimr1fvmjldaYpzXcYBRYzhF9IST+cysSNLjkj7qnKsIfkwg9DZWNaq5o0dzJoyeZVsGu2BariTp5R2sLQgAiE4rth5QS2cPa4cCIXbCIjrMyXzflpSl3jXKysxs3YglBkLkpR0H5TNp1iguouPSEzVjXJperqCIAgCi08Nrq1Q4NklnThrrdRRgVBnW+K4TTeZzzn1a0qeDGw3w1ks76jUlN1XJ8aN7yt0F03L0wGu71dbZo5SE0f3fAgAQXXbVt+r1XYf01Sumy+cbmdtwpk+fPiLHBSJd8FbrBaLIweYOba1t1pwJ6V5H8dwF03LU7Xd6fSfTcwEA0eVPa/Yp1mdBXzt0oGnTpmnatGkjdnwgUlFEgSG81Hcp6tx8imjpxLFKjo/RSxUHvY4CAEDQdHT79dj6al02K0+5aYkj9zodHero6Bix4wORiiIKDOHlHfXKTI5TQWaS11E8Fx/r09KSbL20o17ODblyEwAAEWfF1gM6crRbt545sksKrly5UitXrhzR1wAiEUUUGKTHH9Cqt+s1e0L6qF22ZbALpueo+ki7dta3eR0FAICg+OOb+1Q0NlnnlGR7HQUYlSiiwCAb9jWqpaOHy3IHuHhG7zIuf9te53ESAABOX+XBVr25+7BuWVI0YkOKALw3iigwyN+21ynWZ5o1fvQu2zJYfkaS5uSP0cptFFEAQOTrH1J0w+KRG1IE4L1RRIEBnHN6bssBzRifpqT4GK/jhJXLZo7Thn1HdLCFgQsAgMjV0e3XnzdU64rZ45STluB1HGDUoogCA+yoa9G+w0e1sDDT6yhh5/LZeXJO+vt2pucCACLXX8v3qzEEQ4r6zZo1S7NmzQrJawGRhCIKDLBiS51M0oLCDK+jhJ0Z49JUODaJy3MBABHLOaffvr5HU3JTtbQkKySvWVJSopKSkpC8FhBJKKLAACu2HtCU3FSlJ8V5HSXsmJkumzlOr1Y2qK2zx+s4AACctI1VjSqvbtLHzy4O2WT81tZWtba2huS1gEhCEQX6VB0+qm37mzkb+h4un52nrp6AVlXUex0FAICT9pvX9igtIVYfWBS6IUUvvviiXnzxxZC9HhApKKJAnxVbD0iSFhZRRI+ntDhTmclxeq7vvxUAAJHiYHOHlm/erxtLC5WSEOt1HGDUo4gCfZ7bckAFmUnKTUv0OkrYio3xadmccfrbtjq1d/m9jgMAwLA9+OY++Z3Tx84u9joKAFFEAUlSbWO71u09osVFTMs9kWvmTVBbl18v7mB6LgAgMnT1BPTgm/t00fRcTcxO8ToOAFFEAUnSX8prJUlLJo31OEn4O3NylrJTE/TMplqvowAAMCzLN+9XQ2unPr50otdRAPThAnlA0lNltZqUnaK8MVyWeyIxPtPV88brT2v2qaWjW2mJTBgGAIQv55x++couleSk6Lwp2SF//Xnz5oX8NYFIwBlRjHqVB1u1tbZZZ0zkstzhumb+eHX2BPS37awpCgAIb6t3HtLW2mZ95rzJ8vlCs2TLQMXFxSou5r5UYDCKKEa9pzfVyiQtmchlucO1sDBT+RlJembTfq+jAADwnu59eady0hJ03cJ8T16/sbFRjY2Nnrw2EM4oohjVnHN6uqxG08elKSM53us4EcPXd3nuqop6HW7r8joOAABD2lbbrFfebtBtSycqMS7GkwyvvPKKXnnlFU9eGwhnFFGMauXVTdpz6ChDik7BBxYVqCfg9OTGGq+jAAAwpF++skvJ8TH6yJlcGguEG4ooRrVH1lUpPsZUWsz9oSdr+rg0zStI1yPrquSc8zoOAADvUNPYrqc31ermM4qUnsxgPSDcUEQxarV3+fVUWa0WF49VcjwDpE/FjaWFeutAi7bWNnsdBQCAd3jg1d2SpE+dN8njJACGQhHFqLV88361dvboXA9GuUeL98+foIRYnx5ZV+V1FAAAjjnU2qk/vrlP758/QfkZSV7HATAEiihGrYfXVilvTIKm5aV6HSVipSfF6YrZ4/RUWa06uv1exwEAQJL0y1d2q6PHr89fNMXrKFq0aJEWLVrkdQwg7FBEMSrtqm/Vmj2HdU5JtsxCv6ZYNLmptFBN7d16fhtrigIAvHe4rUu/e32Prpk3QVNyvf+wOT8/X/n53iwdA4QziihGpUfWVctn0tKSLK+jRLylJVkqHJukP7yx1+soAADo/ld2qb3br7su8f5sqCQdOnRIhw4d8joGEHYoohh1Orr9emjtPs0vyGDt0CDw+UwfPatYa3Yf1vb9DC0CAHjnSFuXfrt6j66aO15TctO8jiNJWr16tVavXu11DCDsUEQx6jxVVqPGo926ZGau11Gixk2lhUqM8+m3q/d4HQUAMIr96tXdOtrt112XTPU6CoAToIhiVHHO6YHX9qgwM0nT88Ljk9JokJEcr+sX5uvJsho1Hu3yOg4AYBQ60tal36zeo/fNGa9pvMcDYY8iilHl9V2HtONAiy6ZkceQoiD7+NKJ6ugO6OG1LOUCAAi9e16s1NGuHn3xUs6GApGAIopR5dev7VFaYqzOnDzW6yhRZ8a4MTpz0lj97vW96vEHvI4DRD0zW2ZmO8ys0sy+McTzHzaz8r5/VpvZfC9yAqFQfeSofvf6Xt2wuICzoUCEoIhi1NjT0Ka/bavT+VNzFBfDH/2R8OnzJqumsV3PlNd6HQWIamYWI+keSVdKmiXpFjObNWi33ZIucM7Nk/RdSfeFNiUQOv+1skJm0pcuneZ1lHdZsmSJlixZ4nUMIOwM67fxYXzqOsPMXjezTjP7SvBjAqfv5y/tVFyM6eIZDCkaKZfMyNX0vDT9/KWdCgSc13GAaLZEUqVzbpdzrkvSQ5KuHbiDc261c+5I38M3JBWEOCMQEtv3N+uJjTW67ZyJmpCR5HWcd8nLy1NeXp7XMYCwc8IiOsxPXQ9LukvSfwY9IRAENY3t+vOGap07JUfpSXFex4laPp/pny4qUUVdq/62vc7rOEA0y5c08Ibs6r5tx/MpSc+OaCLAIz947i2lJcTqny4Ij3VDB6urq1NdHe+JwGDDOSM6nE9dDzrn1krqHoGMwGn7xcs7JUlXzOYTyZF21dzxKhqbrHte2innOCsKjJChpq0N+RfOzC5SbxH9+nGev93M1pnZuvr6+iBGBEbei28d1Es76nXnxVOUnhyeHzSvWbNGa9as8ToGEHaGU0RP9lNXIKwcbO7QQ2uqdPbkLGWlJngdJ+rFxvj02Qsma1NVo16tbPA6DhCtqiUVDnhcIOldN2eb2TxJ90u61jl3aKgDOefuc86VOudKc3JyRiQsMBK6egL6j79s0+TsFN22dJLXcQCcpOEU0WF/6nrCA/GpKzxw78u71BMI6Mq547yOMmp8cFGBJqQn6j+fr+CsKDAy1kqaamaTzCxe0s2Snh64g5kVSXpc0kedcxUeZARG1AOv7dbuhjZ9+5pZio9lCCEQaYbzt3ZYn7oOB5+6ItSqDh/V79/Yo6Ul2cpNS/Q6zqiRGBejL102TZuqGrVi6wGv4wBRxznXI+lOSSskbZf0iHNuq5ndYWZ39O32bUlZkn5mZmVmts6juEDQ1TV36H/+/rYunZmnC6czhBCIRLHD2OfYp66SatT7qeutI5oKCJIfPr9DJun98yd4HWXU+cDCfN23apfuXrFDl87MUyxL5gBB5ZxbLmn5oG33Dvj+05I+HepcQCj83+Xb1e13+terZ3odBcApOuFvhsP51NXMxplZtaQvS/qWmVWb2ZiRDA6cyJaaJj1ZVqtLZuZpbEq813FGndgYn75y+XTtrG/T4xtqvI4DAIgSL1fU68myWt1xwWQVZ6V4HeeEli5dqqVLl3odAwg7wzkjOpxPXQ+I9ckQRpxz+v6zbyk1IVZXzuHeUK9cMTtPCwoz9MOVO3TVvPFKSRjW/3IAABhSW2eP/uXxzSrJSdHnLw7P5VoGy8rK8joCEJa4Vg5RaeW2Or1a2aCr5o5Xcjzlxytmpm9fM0t1zZ36yQtvex0HABDh/mtlhWoa2/X9D85TQmyM13GGpaamRjU1XBkEDEYRRdQ52tWjf3t6qwoyknTRDIZieW1RUaZuXFygX72yW5UHW72OAwCIUGVVjfr1a7v1kbOKdMbEsV7HGbYNGzZow4YNXscAwg5FFFHnf16oVG1Thz58ZpFiffwRDwdfv3KGkuJj9G9Pb2U5FwDASWvv8uvLj5QpNy1RX1s2w+s4AIKA39IRVSoPtuiXq3ZpaUmWpualeR0HfbJTE/S/LpumVysb9PSmU1r9CQAwiv2f5du1q75NP7xpvsYkxnkdB0AQUEQRNXr8AX3l0XIlxPl0wyJmZ4Wbj549UQuLMvSdp7fqYEuH13EAABHihbfq9Ps39uoz503SOVOyvY4DIEgooogav1i1S2VVjfrwkmKNSeLT0nAT4zPdfcN8He3y65tPbOESXQDACTW0duprj5Vrxrg0feWK6V7HARBEFFFEhW21zfrvlRU6Y2KmlkyKnAEGo82U3FR95fJpWrmtTk9sZIIgAOD4evwB3fWnjWru6NGPb14YMVNyBzvvvPN03nnneR0DCDsUUUS89i6//vnhMqUmxOrDS4q9joMT+NS5k1VanKlvP7VVuxvavI4DAAhTP1xZodU7D+l7183R9HGRO/chIyNDGRkZXscAwg5FFBHv209tUUVdi25bOlGpiawZGu5ifKYf37JQsTGmzz+4QR3dfq8jAQDCzIqtB/Tzl3bq1jOLdGNpoddxTsvevXu1d+9er2MAYYciioj28Np9enR9ta6eN15z8tO9joNhys9I0g9vnK9t+5v13b9s8zoOACCMVB5s1Vce2aT5Ben6zjWzvI5z2srLy1VeXu51DCDsUEQRsTZXN+lfn9yqWePTdM28CV7HwUm6ZGaePnv+ZD345j49tGaf13EAAGGgvqVTt/16jRLifPrZRxZH7H2hAE6MIoqIVNPYrk/+Zq3GJMXq0+dOls9nXkfCKfjqFdN1/rQcfevJLVpd2eB1HACAh9q7/Pr0b9eqobVTv/r4GcrPSPI6EoARRBFFxGnu6NYnfr1GbV09+sLFU1mqJYLFxvj001sXalJ2iu74w3rtrG/1OhIAwAP+gNMXH9qo8pom/eTmhZpfyHAfINpRRBFROrr9+uzv1mtnfZs+d0EJn5ZGgTGJcXrgtjMUF+PTR+9/UzWN7V5HAgCEUCDg9NXHNun5bXX69tWzdPnscV5HAhACFFFEjM4ev+74/Xq9seuQbjt7omaOH+N1JARJ4dhk/faTS9TS2aMP//INHWzu8DoSACAEnHP65pNb9PiGGn35smn6xDmTvI4UdBdddJEuuugir2MAYYciiojQ1RPQnX/cqJcq6vXRs4p1dkmW15EQZHPy0/WbTyzRwZZOffj+NymjABDlnHP692e26U9r9umfLizRFy6e4nWkEZGamqrU1FSvYwBhhyKKsNfW2aNP/XatVm6r061LinT+tByvI2GELC7O1K8+foZqGtt14y9eV9Xho15HAgCMgB5/QF99rFy/Wb1Hnz53kr56xXSZRefgwZ07d2rnzp1exwDCDkUUYe1wW5du/eUbeq2yQbedPVEXz8j1OhJG2NklWXrw02eq8Wi3brh3td460Ox1JABAEHV0+/W5BzfosfXV+tKlU/XNq2ZGbQmVpG3btmnbNtbMBgajiCJs7TjQomt/+qq27W/W5y4o0blTs72OhBBZWJSphz97lpyTPviz1Xp+6wGvIwEAguBQa6c+9qs1WrmtTv/+/tn60qXTorqEAjg+iijC0nNb9uu6n72mlo4efeXy6VpYlOl1JITYjHFj9PSd52pKbqpu//16/fSFt+Wc8zoWAOAUbd/frGvveU2bqhv1k1sW6uNLJ3odCYCHKKIIK509fn3vr9t0xx82aPyYRH3rqpkqyeEG/9FqXHqiHv7s2bp2wQT95/MV+sRv1qq+pdPrWACAk7R883598Oer1e0P6JHPnq33z5/gdSQAHov1OgDQ760Dzfrin8q0o65FF03P0U2lhYqL4bOS0S4xLkY/+tACLS7O1Pf+ul1X/niV7r5hvi7ifmEACHsd3X79x1+26Y9v7tOCwgz94qOLlTcm0etYAMIARRSe6+j2696Xd+qeFyuVHB+ruy6eonkFGV7HQhgxM33s7Ik6a3KW7vrTRn3iN2t1/cJ8ffOqmcpOTfA6HgBgCNtqm/Wlhzeqoq5Vnz1/sv7X5dMVHzv6PmC+7LLLvI4AhCWKKDz1ckW9vv3UFu09dFRnTMzUrUuKlJYY53UshKlpeWl68vPn6GcvVurnL+/UC28d1NeWTdeHSgsVy9lzAAgLHd1+/c8Lb+sXL+9SRnK8fvvJJbpgFC+9lpjIGWBgKBRReGJTVaPuXrFDr1Y2aNyYRH350mmaNWGM17EQARLjYvTly6fr/Qvy9a0nN+ubT2zRr1/bo69cPl1XzM5j+iIAeGh1ZYO+9eQW7Wpo042LC/TNq2YqIzne61ieqqiokCRNmzbN4yRAeKGIIqS21DTpf154Wyu21ik1IVY3lRbooum53AuKkzYlN1V/+sxZWrG1TneveEt3/GG9FhZl6M6Lpuii6bny+SikABAqlQdb9P1n39Lfth9U0dhk/eFTZ7LsWp8dO3ZIoogCg1FEMeICAacX3jqo+1/ZpTd2H1ZSXIyunT9Bl87MU1J8jNfxEMHMTMvmjNOlM3P12Ppq/c8LlfrUb9dpSm6qPnPeJF27IF+JcfwZA4CRUn3kqH720k49vLZKyXEx+vqyGfrEORP5fy+AE6KIYsRUHT6qP2+o1mPrq1V9pF1jU+J14+ICnTc1W8nx/NFD8MTG+HTzkiJ9cHGB/lq+X79YtUtf//Nmfe+v23XdwnzdVFqoOfnpXscEgKixq75VP3tpp57cWCMz6cNnFumLl0xVFgPkAAwTbQBBVdfcoZXb6vTX8v16fdchmaQZ49N0+3mTtag4Q7E+LsHFyImL8em6hfm6dsEEvb7rkB5eW6WH1lbpd6/v1Yxxabp63ngtmzNeU3JZmxYATpY/4PRyxUH9/vW9eqmiXvExPn3krGLdfv5kTchI8joegAhDEcVp8f+/9u49Rq7qPuD498xjZ/Zpr9frXRsb24CxMWmBxAFKEym1REuRKtqqkSqkhFSpUKXS/lG1SVWpEooqlZBIVVWhJpRWKuojaStaUIt4tEJKI7XBboAQ28EGY/D6uezTO7s7szNz+sed3Z192QbWs9719yMd3XPPPffOGUtzz/7u/d3rauTImVG+d6yflw6d4/WTwwD0dOR44PYt3HNDl1dH1XAhBO65cSP33LiRr41P8ewbp/jX107xzZeO8s2XjnLTpjbu3dvDZ27ayKe2d5pCJkkXcXJwnGdfP8U/vnqSU8MTdLfn+J2fu4kv/MwOutud4yV9NCHGuCIfvG/fvnjw4MEV+Wx9dKVylaPnLvDqu4P8z/EBfnB8gNHJMgA7NrZw+9b13HF9J1vW5X176Qr57M0byWUMrBZzdmSSFw+d5YUfn+XAiUHK1Uguk+LOnRu4+4Yu7ti2np/aus7/Qmj18qTzMTk3a9rp4Qle+PFZnnvj9MxF5rtv2MAX7t7Bz9/a40sGP4RyOfk7KZPx/o+uSUvOzf4itKSxYpnj/WMcOTPKj/pGePPUCEfOjDJVSS5ebGrP8dNb17Ont509ve3X/OvZdfXrXZfnoXt28NA9Oxgrlnn13QG+f2yA77/dzzdeTN5qGALs2tTG7dvWs6e3g9297ezqaaO7LefFFUlr1lSlymvvD/PKW+d55Sfn+cnZCwDcsrmDr963h1+6bTNbO1tWeJSrkwGotLjL+mWEdNeiaQAACuJJREFUEO4D/hxIA0/FGB+btz3Utt8PjANfijH+cJnHqmUWY2R4fIrTIxOcHp6kb2ic4/0F3ukf4+3zY5y/UJzp29KU5voNLezfs4ntG1q5sbvVlFutam25DPv39LB/Tw8Aw+MlXj85PFNePnyOfzrYN9O/syXLzT3t7NzYyrYNLUnpbGbbhha6WpsMUiWtKhOlCq+dHOLAu0McODHID98fYrxUIZMK7NvRyR/dv4f9e3p8pn4ZHD58GIC9e/eu8Eikq8slA9EQQhp4ArgX6AMOhBCeizEeruv2i8CuWrkL+MvaUg0WY2SsWGaoMMVAochgocRAocRgocRQrX5udJJTQxOcGZlkYqoyZ/+WpjS96/LctKmNz9y0kc3r8mxZ30x3e46Uf2hf1XKlU2w//wTZ996E7rvglq9A67aVHtaqsb6lic/t3sTndm8Ckt9S/1iRo2fHOHruAsfOX+Ctsxd4+fA5BgqlOfs2Z9P0dOTobs+xqT1Pd3tutrTlWNeSZV3zbDGlTVIjDRVKHDkzyqHToxw+M8qh0yO801+gUo2EAHt6O/j8p7Zy9w1d/OyujXT4eMKyeueddwADUWm+y7kjeifwdozxOEAI4TvAA0B9IPoA8HRMHjj93xDC+hDC5hjjmWUf8SpSrUbK1Ug1JstKrZSrVSrVSKlcpViuUpyqUixXKJartbbKbHulSnFqer3CWLHCWHGKsWI5qU9OcWGyTKFYrrWVqS7x2G82HWjPZ1nfnKWztYkbu9vobM3S1ZqbWXbkM97ZWYVypVPc/dZ+0pUCKcow/Dq8+/dw/xsGox9RCIFN7Xk2tecX/KfshWKZvqEJTg6O8/7gOH1DE5y/MEn/hSJHzo7yvaNFLhTLSx67pSk9E5R25LO05NK0NKXJZ5NlS1OG5mya5qZkvTmbtDVlUmTTgaZ0qlZPSlMm1NVry3TSN50K/qbXKLOVNC3GyGChxJmRSd4bGOfEQIHj/QVODBQ48UFhzsWz3o48t27p4Bdu7eWT13fyye2drGs28JTUeJcTiF4HnKxb72Ph3c7F+lwHLBmIvjcwzsNPHyQCyfuSIjHCdAwVY5zZFmvrzPSFON0/1tWnDz6vrf5Ytc0Q44J9L3XsaoxUKpFKnA0qK9XZ9XI1Uq1rX+7XQAUgn03T3JQin02Tz6TJZ1O05TN0t+eStmyKlmyG9nxS2nIZ2vNZ2vMZcpmUf5CuUdvPPzEbhALEKSiPwZHHYd9frOzg1qDWXIbdve3s7m1fss9EqcIHY0XOXygyOjHFSK3U10cmphidnGKwUOLUUIXxUoWJqQrjpTKTU9VlG28IkA6BVCqQDklwmgrUlvPaUwv7hhBIpyAQCKH21oHauSTUqskyzKwn25INs31q+9dtn6kv2Hd+22zfme8173s++cV9y/ZvdrUzW2ntizEyOllmqFBiaDwpg4UphsdL9I8VOTsyyZmRSc6OTHJ2dJJSee45o6cjx46uVu7d28POja3s3dLB3s0dPlbTaIWTfGLy26yvHIWDL5mtJNW5nEB0schlfox1OX0IITwMPAyQ79nJgSPvznStP8BsPS44eJh32LDI/ov1XfqYS3/2gvYQSRHJEmkCUkRCraRq/WfWU7H2h1es6zd3PU0kHaqkSUqKOFNPh7igPUVyTKZq5SKmgMFa0dq3e/sLpJrn3YGLU/S98QxP/XvXygxKF9VaK1sW25iGmIIyqaTEFGXSVAhUCVRJUY1hznqFQDUm6xVStfakLdbOhjGGpFSZbatbJn1rdQKVedth7ok9Ljizwvyz9OyL2eefvRc/5lLHXfwsP9+1E4hittKKijFSqlRnMptK5Xn1SmUmq6m+vViu1DKYKowXyxRKZQrFykxW03gpqY9OTjE8PkV5iRSnpnSK3nV5ejvy3L5tPZvX5WfWr+9qYUdXK605X5Cz4gon4fnb2D41SooKHDthtpJU53LOUn1A/a9lK3D6I/Qhxvgk8CTUXhH/p5//UIOVtISDQ3Ds28md0Gkhy9bbfpVHv/zoig1L0hVzRbKV+oYm+IN/fgNYmCnETL0uS2lOn9oy1mUE1TKLqPWbn3k0/xjU913kePPbqWtfcNy6z67G5HGZ+RlN1Tg3u6laZTbDqbZt+jGbZHuyXI7/+a6lKU1rLkPr9DKXYWNbE9u7WmjPZ9nQmqWzpYkNrU10tjTR2dpEZ0vyaE17zsdoVoUjj0N5LAlCwWwlaZ7LCUQPALtCCDuBU8CvAw/O6/Mc8EjtiuxdwIhXXKUGuuUryVXW8lgy0YUsZNqSdklr0RXJVsr13MB//N/xRe5fT6dIL8xUmru+WCbSpbcn64tlPC29feFnLLU9aU+F2ayluRlKST1bq6fC0n1CmM12ms5cmslmWmQ9RTXJfppehipZKmSozs1wKiwYOoVa6Vu4SavEb257hq3N81LYzFbSNebRRx9dctslA9EYYzmE8AjwIskLEf4mxngohPBbte3fAp4neRnC2yQvRPiNjz9sSZetdVuS6nPkcRh4Fbru9DkUaW27ctlKj/3K8o5UulYdHDBbSbqIy3qAIMb4PEmwWd/2rbp6BH57eYcm6UNp3Waqj3TtMFtJutqZrSRdlE+yS5K0ypitJK0CZitJF2UgKknSKmS2krQKmK0kLSm10gOQJEmSJF1bDEQlSZIkSQ1lICpJkiRJaigDUUmSJElSQxmISpIkSZIaykBUkiRJktRQBqKSJEmSpIYyEJUkSZIkNZSBqCRJkiSpoQxEJUmSJEkNZSAqSZIkSWooA1FJkiRJUkMZiEqSJEmSGspAVJIkSZLUUAaikiRJkqSGMhCVJEmSJDWUgagkSZIkqaEMRCVJkiRJDRVijCvzwSH0A++tyIdLa9tG4IOVHoS0Aj6IMd630oNYzZybpSvGuVnXqiXn5hULRCVdGSGEgzHGfSs9DkmSlHBulhYyNVeSJEmS1FAGopIkSZKkhjIQldaeJ1d6AJIkaQ7nZmkenxGVJEmSJDWUd0QlSZIkSQ1lICqtISGEnSGEH4QQjoUQvhtCaFqi3/UhhJdCCEdCCIdDCDsaO1JJkq4Nzs3S4gxEpbXl68CfxRh3AUPAl5fo9zTwjRjjLcCdwPkGjU+SpGuNc7O0CANR6SoUQvh0COFHIYR8CKE1hHAohPCJS+wTgP3Av9Sa/hb45UX67QUyMcaXAWKMYzHG8WX+CpIkrSnOzdLyyqz0ACQtFGM8EEJ4DvgToBn4O+C9EMLrS+zyIMmV0+EYY7nW1gdct0jfm4HhEMIzwE7gP4E/jDFWlvM7SJK0ljg3S8vLQFS6en0NOABMAr9bm4xuX6pzCKF7kebFXoudAT4L3AG8D3wX+BLw1x9zvJIkrXXOzdIyMRCVrl4bgDYgC+RDCCngv5fo+yBwBFgfQsjUrrxuBU4v0rcPeC3GeBwghPBvwN042UmSdCnOzdIyMRCVrl5PAn9MkqLz9RjjI1zkqitACOEV4NeA7wAPAc8u0u0A0BlC6I4x9pM8u3JwOQcuSdIa5dwsLRNfViRdhUIIXwTKMcZ/AB4DPh1C2H8Zu34V+L0QwttAF7UrqSGEfSGEpwBqaUS/D/xXCOFNIAB/dQW+hiRJa4Zzs7S8QoyLpalLkiRJknRleEdUkiRJktRQBqKSJEmSpIYyEJUkSZIkNZSBqCRJkiSpoQxEJUmSJEkNZSAqSZIkSWooA1FJkiRJUkMZiEqSJEmSGur/ATUMKhslmzVWAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# nbi:hide_in\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "from ipywidgets import interact, FloatSlider\n", "plt.rcParams['figure.figsize'] = (16, 5)\n", "\n", "def cdf_func(xdata):\n", " val = np.divide(np.exp(xdata), 1+np.exp(xdata))\n", " return val\n", "\n", "def pdf_func(xdata):\n", " val = np.divide(np.exp(xdata), (1+np.exp(xdata)**2))\n", " return val\n", "\n", "def logistic_cdf_pdf(x):\n", " xdata = np.linspace(-10, 10, 1000)\n", " \n", " fig, [ax1, ax2] = plt.subplots(1, 2)\n", " \n", " ax1.plot(xdata, pdf_func(xdata))\n", " xshade = xdata[xdata<=x]\n", " ax1.fill_between(xshade, pdf_func(xshade), alpha=0.3, zorder=1)\n", " ax1.scatter(x,0, s=30, zorder=2, c=\"orange\")\n", " ax1.axhline(y=0, color='k', linewidth=0.5, zorder=1)\n", " ax1.set_xlim(-10, 10)\n", " ax1.set_ylim(-0.06,0.6)\n", " ax1.set_xticks([x])\n", " ax1.set_xticklabels([\"x={}\".format(x)])\n", " ax1.set_title(\"pdf\")\n", " ax1.spines['top'].set_visible(False)\n", " ax1.spines['right'].set_visible(False)\n", " \n", " ax2.plot(xdata, cdf_func(xdata))\n", " ax2.vlines(x, 0, cdf_func(x), linestyle=\"dashed\", alpha=0.4)\n", " ax2.scatter(x,0, s=30, zorder=2, c=\"orange\")\n", " ax2.axhline(y=0, color='k', linewidth=0.5, zorder=1)\n", " ax2.set_xlim(-10, 10)\n", " ax2.set_ylim(-0.1,1.1)\n", " ax2.set_xticks([x])\n", " ax2.set_xticklabels([\"x={}\".format(x)])\n", " ax2.set_title(\"cdf\")\n", " ax2.spines['top'].set_visible(False)\n", " ax2.spines['right'].set_visible(False)\n", " \n", " plt.show();\n", " \n", "logistic_cdf_pdf(0.6); " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To apply the inverse transform method we note that \n", "\n", "* $y=\\frac{e^x}{1+e^x}\\Rightarrow e^x=\\frac{1-y}{y}.$ \n", "* $F^{-1}(y)=\\log(1-y)-\\log(y).$" ] }, { "cell_type": "code", "execution_count": 33, "metadata": {}, "outputs": [], "source": [ "# The following code generates samples from the logistic distribution\n", "import numpy as np\n", "\n", "def logistic_inv(y):\n", " return np.log(1-y)-np.log(y)\n", " \n", "# random sample from uniform distribution on [0,1)\n", "y = np.random.sample(20) \n", "\n", "# random sample from logistic distribution\n", "x = logistic_inv(y)" ] }, { "cell_type": "code", "execution_count": 34, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAasAAAFWCAYAAADXHbKqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXyV5YH28d+dk40lAbIRSNiCCEkgBIiAOygiLhUXrFq76OhQ2rG1+rbU6XRedXzbUavTdiqVsa1Up53SVltFZVzAoohSNkGWsCaBhED2ELKcJCfnfv9IwBAScoAkz1mu7+eTT3LO85yTK5DkynOf57lvY61FRETEn4U5HUBERKQ7KisREfF7KisREfF7KisREfF7KisREfF74U594nnz5tm3337bqU8vIiL+yXR2p2NHVuXl5U59ahERCTAaBhQREb+nshIREb+nshIREb/n2AkWnWlubqaoqAi32+10FAGio6NJTU0lIiLC6SgiEuL8qqyKioqIiYlh9OjRGNPpCSHSR6y1VFRUUFRUxJgxY5yOIyIhzq+GAd1uN/Hx8SoqP2CMIT4+Xke5IuIX/KqsABWVH9H/hYj4C78rKxERkY5UVu1UVFSQnZ1NdnY2ycnJpKSknLzd1NTUp1lWrVrFzTffDEBDQwNXXXUV2dnZvPLKK32aQ0TEH/jVCRZOi4+PZ+vWrQA89thjDBw4kO9+97un7GOtxVpLWFjf9fzmzZsxxpzMJiISarr9jWuMedEYU2qM2dHFdmOM+U9jzH5jzGfGmKk9H9NZ+/fvZ+LEiSxatIipU6dSWFjI4MGDT25fvnw5999/PwAlJSXceuut5OTkMH36dNavX3/a83k8Hh566CEmTpxIVlYWv/zlLwF46623GD9+PJdddhmvv/46AMXFxdxzzz1s2rSJ7OxsCgoKev8LFhHxM74cWf0WeA54uYvt1wHj2t5mAM+3vT8vj7+xk13FNef7NKfIGB7Lo1/IPKfH7tq1i2XLlrF06VI8Hk+X+337299m8eLFzJw5k4KCAm688UZ27Di1559//nmKi4vZtm0bLpeLyspK6uvr+frXv84HH3xAWloaCxYsAGD48OEsXbqU5557jtdee+2csouIBLpuy8pa+6ExZvQZdpkPvGyttcB6Y8xgY8wwa+2RHsroF8aOHctFF13U7X6rVq1iz549J29XVVXR0NBAv379TtnnO9/5Di6XC4C4uDg2bdrEhRdeyNixYwG4++67efnlrv4+EJFAY62lxWtpaXvv8Vq8Hd63eNttsxZPS/vHeGnxgsfrxVqwFrzWYjnx8gRYbNv9bfedsq1t/7bHtWb6/L4T22l7Hm/b5zjxnPbzL6T13ak3iR8YyY1Zw3vt368nXrNKAQrb3S5qu++0sjLGLAQWAowcOfKMT3quR0C9ZcCAASc/DgsLw9qT/3WnXItkrWXDhg1ERkZ2+VzW2k5PC9ep4iK9J+6pOKrcVU7HOCujGt50OoLPslIH+X1ZdfYb1nZyH9baF4AXAHJycjrdJxCEhYUxZMgQ9u3bx9ixY/nrX/9KYmIiAHPmzGHJkiU89NBDAGzdupXs7OxTHj937lyef/55Lr/88pPDgBkZGezdu5f8/HxGjx7NH/7whz7/ukSCRYvXUlHXSGlNI2XHGyk97qbKXcW/Zm+noq6JY/XNVDc0cayhmer6Zo67ux7aby/SFUb/KBf9I1xER7qIdIURFeEiyhVGZHjbW/uP225HRYQR5Qoj3BWGK8zgCjOEt70/+Wba7ncZwowhPCwMVxi4wsIIDzOEtT0mzBjCDISFGQxgTOsfuq0ft24zmLb7P/84zJx4D7S/r+05wtr+WG6/b/vnaXvUyX1ab7e9b8vem3qirIqAEe1upwLFPfC8fu2pp55i3rx5jBw5koyMDBobGwFYsmQJ3/jGN1i2bBkej4fZs2ezZMmSUx779a9/nX379pGVlUV4eDjf+MY3WLRoEUuXLuW6664jISGBSy+99JThRBH5nLWW0uONFJTXUVjVQGFlfetbVT2FlQ2UHnfj7fjncD94fWsx8QMjGdwvgqSYaC5MiiG2XwSD+0cwqO19TFQEA6LCGRDlon9kOP0jXQyIDKdfpIvIcF3t4xTTfjiry51aX7N601o7sZNtNwAPANfTemLFf1prp3f3nDk5OXbTpk2n3Jebm0t6erpPwaVv6P9EnFZd38SOwzXsKTnO/tLj7C2pZV/JcWraHQ0ZA8Nio0mN68+IIf0ZPjiapJgoEmOiSYqNIikmihH/OQD7aMAO6ISSTg/Ruj2yMsb8AZgFJBhjioBHgQgAa+1SYCWtRbUfqAfu7Zm8IhJqmjxeth8+xtbCaj4rqmZbYTUFFfUntw/pH8G4oTHclD2ccUkxjEkYwMi4/gwf3E9HPUHOl7MB7+pmuwX+qccSiUjIaPFadhXX8PGBcj4+UMHGgkrqm1oASI6NZvKIQdyeM4Ks1EFMSI4lYWCkTkQKUZrBQkT6VENTC2v3lfHurhLe311KZV3rVGbjkgZy+7RULh4bz5SRQxgaG+1wUvEnKisR6XVNHi/v7y7ltU8Ps2ZvKe5mLzHR4Vw1IYmrJiRxcVo8SSonOQOVlYj0mh2Hj/HK5iJe33qYqvpmEgZGcfu0EczNHMqMMfF6nUl8prISkR7lafHy7q4Slq3LZ2NBFZHhYVyTMZQFU1O5fFwC4S4VlJw9fdd04HK5Ti4Lkp2dzZNPPtkjz3v99ddTXV191tvO5Le//S0PPPDAOeVp/9iysjJmzJjBlClTWLt27Tk9n4i7uYXfrsvnyp+s4Zu/38LRGjf/emMGG38whyVfmsrsCUkqKjlnOrLqoF+/fr2yFMfKlStPu+/EciOdbetLq1evZsKECbz00kuO5pDA1Ohp4Y8bC/nl3w5wtMbN9NFxPPqFDK5OH9rrsxpI6NCfOT4aPXo0P/jBD7j44ovJyclhy5YtXHvttYwdO5alS5cCsGbNGq644gpuueUWMjIyWLRoEV6v9+Tjy8vLKSgoID09nW9+85snlxs5sQ3g5ZdfJisri8mTJ/OVr3wFgDfeeOPkkc+cOXMoKSk5Y9ba2lruvfdeJk2aRFZWFq+++ioAy5Yt48ILL+TKK69k3bp1QOt0UIsXL2blypVkZ2fT0NDQK/9+EnystazYVszsn6zh/76+k5Fx/fmff5zBnxZdzNzMZBWV9CgdWXXQ0NBwylx+//zP/8wdd9wBwIgRI/jkk0946KGHuOeee1i3bh1ut5vMzEwWLVoEwIYNG9i1axejRo1i3rx5/OUvfzm53McJe/bsYdmyZSfXsTph586d/OhHP2LdunUkJCRQWVkJwGWXXcb69esxxvDrX/+ap59+mmeffbbLr+GJJ55g0KBBbN++HWid+f3IkSM8+uijbN68mUGDBjF79mymTJlCdnY2//Zv/8amTZt47rnnzv8fUELCjsPHePyNnWwsqGJiSixPL5jMpRfE6xoo6TV+XVbm8Z7/xu9uupUzDQPedNNNAEyaNIna2lpiYmKIiYkhOjr65GtO06dPJy0tDYC77rqLjz766LSyGjVqFDNnzjzt+d9//30WLFhAQkIC0Lp0CEBRURF33HEHR44coampiTFjxpzxa1i1ahXLly8/eXvIkCG89tprzJo16+SEu3fccQd79+494/OIdNTQ1MIz7+7hxXX5xPWP5MlbJ3F7zggdRUmv8+uy8rd5vKKiooDWWddPfHzi9okFGTv+ZdnZX5rtlxtpr6ulQ771rW/x8MMPc9NNN7FmzRoee+yxM+bUEiTSUY8vjxENBV64639b3wLBkOghTkeQ8+DXZRWINmzYQH5+PqNGjeKPf/wjCxcu9PmxV199NbfccgsPPfQQ8fHxVFZWEhcXx7Fjx0hJSQHw6SSIuXPn8txzz/Gzn/0MaB0GnDFjBg8++CAVFRXExsby5z//mcmTJ5/bFykBp8pddc5//DW3eHnm3T288GEeqUP68dRtWVwyNqGHE4qcmU6w6ODEa1Yn3h555JGzevzFF1/MI488wsSJExkzZgy33HKLz4/NzMzkX/7lX7jyyiuZPHkyDz/8MACPPfYYt99+O5dffvnJIcIz+eEPf0hVVRUTJ05k8uTJ/O1vf2PYsGE89thjXHzxxcyZM4epU6ee1dcloam4uoE7X1jPf32Qx50XjeTtB69QUYkjfFoipDcE4xIha9as4ZlnnuHNNwNndc/uBPr/ibQyj5uzPrL6cG8ZDy7/lCaPlx/fOon52Sm9lE7kFOe2RIiIhJ6XPi7g8Td2cuHQGJbcPZWxiQOdjiQhTmXVg2bNmsWsWbOcjiFyzjwtXp54cxcvfXKQOelD+fmd2QyI0q8JcZ7ffRd2dSab9D2nhojFGe7mFr71h095b1cJ9182hn++Pl2npIvf8Kuyio6OpqKigvh4XVzoNGstFRUVREdr2YZQUN/kYeHLm/lofzmPfSGDey4987V8In3Nr8oqNTWVoqIiysrKnI4itP7xkJqa6nQM6WU17mbuXbaRTw9V8ZMFWdyeM8LpSCKn8auyioiI6HZ2BhHpOfVNHv5h2Ua2FVbzi7umckPWMKcjiXTKr8pKRPqOu7mFhS9vZsuhKp770lSun6SiEv+lshIJQZ4WL9/6w6d8tL+cZ26frKISv6cZLERCjLWWR1fs5L1dJTx+UyYLpul1SfF/KiuREPObj/L5/d8PsejKsXztktFOxxHxicpKJIS8u/MoP1qZy3UTk1l87Xin44j4TGUlEuTinopjSPQQ9hw9zoPLt5KVOpj/+GI2YbrgVwKITrAQCXJV7ipqvt/E/OfWMSAqnF99ZRr9Il1OxxI5KyorkRDw/Vc/42BlPb+/fwZJsZqVRAKPhgFFQsDK7Uf53rXjmZkW73QUkXOishIJYjsOHwPgmoyhfP2KNIfTiJw7lZVIkHI3t/CdP24F4OnbsjQ5tAQ0lZVIkHrq7d3sL60FYMiASIfTiJwflZVIEFq7r4xl6wr42sWjnI4i0iNUViJBprbRw/df+YyxiQN45Lp0p+OI9Aidui4SZJ59dw9Haty8sugSXU8lQUNHViJBZGthNb/9uIAvzxjFtFFDnI4j0mNUViJBornFyyOvfkZSTBTfm6d5/yS4aBhQJEj85qN8dh89zn99ZRqx0RFOxxHpUTqyEgkCJTVu/nP1PuakJ3FtZrLTcUR6nMpKJAg8/fYePC2WH96Q4XQUkV6hYUCRALe1sJpXtxRROvBuxiw5dtr2IdE60UICn8pKJIB5vZbHVuwkMSaKg55j2Eet05FEeoWGAUUC2IptxWwtrOb78yY4HUWkV6msRAJUk8fLM+/uIXN4LLdOSXE6jkivUlmJBKg/bDhEUVUDi+dN0BL1EvRUViIBqK7Rwy/e38fMtDiuGJfgdByRXqeyEglAL36UT3ltE4vnTdA6VRISVFYiAaaqrokXPsxjbsZQpo7UaekSGlRWIgHm1x/lUdvk4bvXav4/CR0qK5EAcqy+mZc+Psj1E4dx4dAYp+OI9BmVlUgAWfZxPrWNHh646gKno4j0KZ/Kyhgzzxizxxiz3xjzSCfbBxlj3jDGbDPG7DTG3NvzUUVC23F3My9+lM81GUNJHxbrdByRPtVtWRljXMAS4DogA7jLGNNxtsx/AnZZaycDs4BnjTGRPZxVJKS9/MlBatwevn3VOKejiPQ5X46spgP7rbV51tomYDkwv8M+FogxrefQDgQqAU+PJhUJYfVNHn7zUT6zxicyKXWQ03FE+pwvZZUCFLa7XdR2X3vPAelAMbAdeNBa6+2RhCLCK5uLqKxr4p9m67UqCU2+lFVnVxx2nNr5WmArMBzIBp4zxpw2qG6MWWiM2WSM2VRWVnbWYUVCUYvX8puP8skeMZicUZ1fVxX3VJyWApGg5ktZFQEj2t1OpfUIqr17gb/YVvuBfOC0aaCttS9Ya3OstTmJiYnnmlkkpLy3q4SDFfX84+VpXc5WUeWuovL7lX2cTKTv+FJWG4FxxpgxbSdN3Ams6LDPIeBqAGPMUGA8kNeTQUVC1a/W5pE6pB/XZg51OoqIY7otK2utB3gAeAfIBf5krd1pjFlkjFnUttsTwCXGmO3AauD71try3gotEiq2HKpi88Eq7rtsDOEuXRYpocunlYKttSuBlR3uW9ru42Jgbs9GE5Ffr80jNjqcL+aM6H5nkSCmP9VE/FRhZT1v7zjKl2aMYkCUT39XigQtlZWIn/rv9QcxxnDPJaOdjiLiOJWViB9yN7fwp02FXJs5lORB0U7HEXGcykrED72xrZjq+ma+PHOU01FE/ILKSsQP/W79QS5IGsjFafFORxHxCyorET+zrbCabUXH+MrMUVqyXqSNykrEz/xu/UH6R7q4ZWrHKThFQpfKSsSPVNU1sWJbMTdPSSE2OsLpOCJ+Q2Ul4kde2VxEo8fLV3RihcgpdKWhiJ+w1vLHTYVMGTm4y5WA456Ko8pdddr9mnFdgp3KSsRPbDlUzf7SWp68dVKX+1S5q7CPdlyhRyT4aRhQxE/8eVMh/SJc3Dh5uNNRRPyOykrED9Q1enhjWzE3ZA1joOYBFDmNykrED6zcfoS6phbuuEizq4t0RmUl4gf+tKmQtIQBXS5bLxLqVFYiDssrq2VjQRW354zQjBUiXVBZiTjsT5uKcIUZbtOMFSJdUlmJOMjT4uXVLUXMHp9IUqyWAhHpispKxEEf7S+n7HgjC6alOh1FxK+prEQc9PrWYmKjw5k9IcnpKCJ+TWUl4pC6Rg9v7zjKDVnDiQp3OR1HxK+prEQc8t6uEhqaW7g5WzNWiHRHZSXikNe2HiZlcD8uGh3ndBQRv6eyEnFA2fFG1u4rZ372cMLCdG2VSHc0CZmIA978rJgWr+XmKWe+tqrjkiBaCkRClcpKxAGvbS0mY1gsFw6NOeN+WhJEpJWGAUX6WF5ZLdsKq7mlm6MqEfmcykqkj722tRhj4Atat0rEZyorkT5kreXNz4qZOSae5EGaXknEVyorkT60++hx8srquHHyMKejiAQUlZVIH3rrsyOEGZiXmex0FJGAorIS6SMnhgAvGZtA/MAop+OIBBSVlUgf2VlcQ0FFPTdkaQhQ5GyprET6yFvbj+AKM1yrIUCRs6ayEukD1lre+uwIl16QQNyASKfjiAQclZVIH9hxuIZDlfXcOElDgCLnQmUl0gfe3F5MeJhhbuZQp6OIBCSVlUgvOzEEeNm4BAb31xCgyLlQWYn0sm1FxyiqauAGDQGKnDOVlUgvW7n9CBEuw9yMszsLMO6pOC0JItJGS4SI9CJrLW/vOMolYxMY1D/irB6r5UFEPqcjK5FetPvocQ5V1jNvoq6tEjkfKiuRXvT2jqMYA3PSdRagyPlQWYn0ond2HuWiUXEkxmguQJHzobIS6SUHK+rYffS4rq0S6QEqK5Fe8s7OowCaC1CkB6isRHrJOztLyBwey4i4/k5HEQl4KiuRXlBa42bzwSodVYn0EJWVSC94d1cJgE5ZF+khKiuRXvDOzqOMSRjAuKSBTkcRCQoqK5Eedqy+mU8OVHBtZjLGGKfjiAQFn8rKGDPPGLPHGLPfGPNIF/vMMsZsNcbsNMZ80LMxRQLH6t0leLyWa3XKukiP6XZuQGOMC1gCXAMUARuNMSustbva7TMY+CUwz1p7yBiT1FuBRfzdOzuPkhwbzeTUwU5HEQkavhxZTQf2W2vzrLVNwHJgfod9vgT8xVp7CMBaW9qzMUUCQ0NTCx/sLWNu5lDCwjQEKNJTfJl1PQUobHe7CJjRYZ8LgQhjzBogBvi5tfbljk9kjFkILAQYOXLkueQV8Wtr95XhbvaeXA4k7qk4qtxV5/RcWh5E5HO+lFVnfx52XLcgHJgGXA30Az4xxqy31u495UHWvgC8AJCTk6O1DyTorM4tJSYqnOlj4gAt8yHSU3wpqyJgRLvbqUBxJ/uUW2vrgDpjzIfAZGAvIiHC67Ws3l3KFeMTiQzXibYiPcmXn6iNwDhjzBhjTCRwJ7Ciwz6vA5cbY8KNMf1pHSbM7dmoIv5tW1E15bWNXKPlQER6XLdHVtZajzHmAeAdwAW8aK3daYxZ1LZ9qbU21xjzNvAZ4AV+ba3d0ZvBRfzN6txSXGGGWeMTnY4iEnR8WtbeWrsSWNnhvqUdbv8E+EnPRRMJLKtyS5g2agiD+0c6HUUk6GhgXaQHFFXVs/voceak6xJDkd6gshLpAatzWy8t1PL1Ir1DZSXSA1bllpCWMIC0RE1cK9IbVFYi56m20cPf8yq5WkOAIr1GZSVyntbuLaOpxcvVGgIU6TUqK5Hz9F5uCYP6RZAzStMjifQWlZXIeWjxWtbsKWP2+ETCXfpxEukt+ukSOQ+fHqqisq5JQ4AivUxlJXIeVuWWEh5muFKzVoj0KpWVyHlYlVvCjLQ4YqMjnI4iEtRUViLn6GBFHftLa7l6goYARXqbykrkHK3SrBUifUZlJXKOVueWMC5pICPj+zsdRSToqaxEzsGxhmY25FcyJ0NHVSJ9QWUlcg4+2FuGx2s1y7pIH1FZiZyD1bklxA2IJHuEZq0Q6QsqK5Gz1Nzi5W+7S5k9PglXmHE6jkhIUFmJnKVNBVXUuD1ck6EhQJG+orISOUurc0uIdIVx+TjNWiHSV1RWImfBWsuq3BJmjo1nQFS403FEQobKSuQsHCiro6CiXmcBivQxlZXIWVidWwKgWdZF+pjKSuQsrM4tJX1YLCmD+zkdRSSkqKxEfFRV18Smg5UaAhRxgF4hFvHR3/aU4rW+DQHGPRVHlbuKIdG6aFikJ6isRHy0OreUxJgoslIGdbtvlbsK+6jtg1QioUHDgCI+aPJ4+WBvGVdPSCJMs1aI9DmVlYgPNuRXUtvo0VmAIg5RWYn4YFVuCVHhYVx2QYLTUURCkspKpBsnZq247IIE+kW6nI4jEpJUViLd2FtSS1FVg4YARRykshLpxqqTs1bo+ioRp6isRLqxKreErNRBDI2NdjqKSMhSWYmcQdnxRrYWVnP1BA0BijhJZSVyBn/bXYq1MEcLLYo4SmUlcgarcksYPiiajGGxTkcRCWkqK5EuuJtbWLuvnKvTh2KMZq0QcZLKSqQLnxyooKG5RWcBivgBTWQr0oVVuSX0j3QxMy3+5CzqvtJs6yI9S2Ul0glrLatzS7liXCLRES7Noi7iMA0DinRiZ3ENR2vcGgIU8RMqK5FOrMotwRi4aoLKSsQfqKxEOrEqt4SpI4cQPzDK6SgigspK5DRHjjWw43CNhgBF/IjKSqSD1bmlAFyjWdZF/IbKSqSD1bkljIzrzwVJA52OIiJtVFYi7dQ3eVh3oIKr05M0a4WIH1FZibSzdl85TR6vhgBF/IzKSqSd1bklxESHc9GYOKejiEg7KiuRNi1ey/u7S7nywkQiXPrREPEnPv1EGmPmGWP2GGP2G2MeOcN+FxljWowxC3ouokjf2HKoivLaJq7NTHY6ioh00G1ZGWNcwBLgOiADuMsYk9HFfk8B7/R0SJG+8O7Oo0S6wpg1PtHpKCLSgS9HVtOB/dbaPGttE7AcmN/Jft8CXgVKezCfSJ+w1vLOzhIuuSCemOgIp+OISAe+lFUKUNjudlHbfScZY1KAW4ClZ3oiY8xCY8wmY8ymsrKys80q0mv2lBznUGU9czNOHQKMeyoO87jRkh8iDvNliZDOLjbpuFbCz4DvW2tbznRtirX2BeAFgJycHK23IH7jnR2tE9fOyTh1iiUtDSLiH3wpqyJgRLvbqUBxh31ygOVtRZUAXG+M8VhrX+uRlCK97N1dR5k6cghJMdFORxGRTvgyDLgRGGeMGWOMiQTuBFa038FaO8ZaO9paOxp4BfimikoCRWFlPTuLa5iboQuBRfxVt0dW1lqPMeYBWs/ycwEvWmt3GmMWtW0/4+tUIv7uvV0lAMzVKesifsunZe2ttSuBlR3u67SkrLX3nH8skb7zzs6jXDh0IGMSBjgdRUS6oMv0JaRV1jWxsaBSFwKL+DmVlYS0VbkleC2nnbIuIv5FZSUh7d2dJQwfFM3ElFino4jIGaisJGTVN3lYu6+MuZnJWrtKxM+prCRkfbi3jEaPl7mZOmVdxN+prCRk/e+OowzuH8H00Vq7SsTfqawkJLmbW1idW8q1GcmEa+0qEb+nn1IJSR/uLaO20cP1WcOcjiIiPlBZSUhauf0Ig/tHcMnYeKejiIgPfJrBQiSYNHpaWJVbyvWTkk9bvj7uqTiq3FUnb2tpEBH/oLKSkLN2bzm1jR6um3T6EKCWBBHxTxoGlJCzcvsRYqPDuXRsgtNRRMRHKisJKY2eFt7bVcLczGQiw/XtLxIo9NMqIeWjfeUcb/RwQydDgCLiv1RWElJWbj9KTHQ4l16gIUCRQKKykpDR5PHy3q6jXJMxVEOAIgFGP7ESMtbtL6fGrSFAkUCkspKQ8ca2YmKjw7lsnIYARQKNykpCQkNTC+/sPMr1k4YRFe5yOo6InCWVlYSE1btLqGtq4abs4U5HEZFzoLKSkPD61mKGxkYxY4zmAhQJRCorCXrH6ptZs6eUL2QNxxWmFYFFApHKSoLe/+44QnOLZX52itNRROQcqawk6L2+tZi0hAFMTIl1OoqInCPNui5B7egxN+vzK/j2VeMw5tQhwI7LgYCWBBHxVyorCWpvflaMtXR6FqCWAxEJHBoGlKD2+tZiJqUMYmziQKejiMh5UFlJ0NpXcpzth48xX9dWiQQ8lZUErVe2FBEeZrh5is4CFAl0KisJSp4WL3/dcphZ45NIGBjldBwROU8qKwlKa/eXU3q8kQXTUp2OIiI9QGUlQemVzUUM6R/BVROSnI4iIj1AZSVB51h9M+/tLGF+dooWWRQJEvpJlqCz4rNimlq8GgIUCSIqKwk6r2wuYkJyDJnDNb2SSLBQWUlQ2VtynG2F1SyYlnra9EoiErhUVhJU/ufvh4h0hXHrVA0BigQTlZUEDXdzC3/ZUsS1E5OJGxDpdBwR6UEqKwkab352hBq3hy9NH+l0FBHpYZp1XYLGHzYcIi1xADPT4k7e19kyICdoORCRwKGykqCw5+hxNh+s4oc3pJ9yYoWWAREJDhoGlKDwP38/SGR4GLfpxAqRoB+VAv0AABA3SURBVKSykoDX0NTCXz49zPUTkxmiEytEgpLKSgLeG9uKOe72cJdOrBAJWiorCWjWWl5cl8+E5Bimj4nr/gEiEpBUVhLQ1udVsvvoce69dLRmrBAJYiorCWjL1uUzpH8E87O1GrBIMFNZScAqrKznvdwSvjRjJNERLqfjiEgvUllJwHrp4wJcxvCVmaOdjiIivUxlJQGpttHDHzcVct2kYSQPinY6joj0Mp/Kyhgzzxizxxiz3xjzSCfb7zbGfNb29rExZnLPRxX53Kubizju9nDvpaOdjiIifaDbsjLGuIAlwHVABnCXMSajw275wJXW2izgCeCFng4qcoKnxcuv1uYxZeRgpo7U/H4iocCXI6vpwH5rbZ61tglYDsxvv4O19mNr7YnZQtcDmvNGes1b249QVNXAN2dd4HQUEekjvkxkmwIUtrtdBMw4w/73Af/b2QZjzEJgIcDIkZptQM6etZbn1xxgXNJArp6QdNr2jrOsa2Z1keDgS1l1dqVlp9NYG2Nm01pWl3W23Vr7Am1DhDk5OZoKW87amj1l7D56nGdvn0xY2OnfmpplXSQ4+VJWRcCIdrdTgeKOOxljsoBfA9dZayt6Jp7IqZ5fc4Dhg6K5KXu401FEpA/58prVRmCcMWaMMSYSuBNY0X4HY8xI4C/AV6y1e3s+pghsKqhkQ0El/3hFGhEuXXUhEkq6PbKy1nqMMQ8A7wAu4EVr7U5jzKK27UuB/wvEA79sm5/NY63N6b3YEop+tmof8QMiueOiEd3vLCJBxaeVgq21K4GVHe5b2u7j+4H7ezaayOf+nlfBR/vL+eEN6fSP1ALXIqFGYykSEH66ai+JMVHcPWOU01FExAEqK/F7Hx8oZ31eJd+cNZZ+kZqwViQUqazEr1lr+el7e0mOjdZKwCIhTGUlfm3N3jI2FlTxT7PHahkQkRCmshK/5Wnx8u8rcxkd3587LtJRlUgoU1mJ33plcxF7S2r5/rwJRIbrW1UklOk3gPilukYPz763l2mjhjBvYrLTcUTEYSor8UsvfJhH2fFGfnB9Om0XmotICFNZid85cqyBFz7M44ZJw5g2SrOmi4iPM1iI9KUn3tyF11oeuW7Cads6LgHSkZYEEQlOKivxKx/uLWPl9qM8fM2FjIjrf9p2LQEiEpo0DCh+o9HTwqMrdjI6vj8Lr0hzOo6I+BEdWYnfeOGDPPLL63jpH6brAmAROYWOrMQvHCir5bm/7ee6iclceWGi03FExM+orMRxLV7L9/68jegIF4/flOl0HBHxQxoGFMctW5fPlkPV/PSOySTFRjsdR0T8kI6sxFF5ZbX85J09zElP4ubsFKfjiIifUlmJY5pbvPyfP28jKjyMH98ySTNViEiXNAwojvnpe3v59FA1v7hriob/ROSMdGQljvhoXznPf3CAO3JG8IXJw52OIyJ+TmUlfa7seCMP/WkrYxMH8uhNGU7HEZEAoGFA6VOeFi8PLv+UYw3N/Pd90+kfqW9BEemejqykT/145W4+PlDB/7t5IhOSY52OIyIBQmUlfeaVzUW8uC6fey4ZzRdzRjgdR0QCiMZgpE9sOVTFD/66nYvT4vmXG9J9ekxny4FoCRCR0KSykl6XV1bLfb/dSHJsNEvunkqEy7cDei0HIiInaBhQelVpjZuvvriBMGN46R+mEzcg0ulIIhKAdGQlvabG3czXlm2ksq6J5QtnMiZhgNORRCRA6chKekWNu5mv/mYD+0qO8/yXp5GVOtjpSCISwFRW0uNOFNWOw8f45d1TtT6ViJw3DQNKjzpW38zXln1eVHMzk52OJCJBQGUlPaa4uoGvvbiBgxX1KioR6VEqK+kRe44e52svbqCu0cNL/zCdi8fGOx1JRIKIykrO2/u7S3jwD1vpH+XiT4suJn2YplESkZ6lspJz5vValvxtP/+xai/pybG88NVppA7p73QsEQlCKis5J1V1TXz/1c94d1cJ87OH8+StWfSLdDkdS0SClMpKztrH+8t5+E/bqKhr5Ic3pHPfZWO0JL2I9CqVlfjM3dzCf7y3l1+tzWNMwgB+9dVLmZQ6yOlYIhICVFbikw/3lvHD13ZwqLKeu6aP5F9vTNfCiSLSZ/TbRs5oyJNxVDe2W6ajHzy5vfWt1z+3lgMRkTYqK+lUjbuZpWsOUN1YxbimlXxj1li+MWss0RE6iUJE+p7KSk5R2+jh9+sP8l8f5lFZ1wT9YPX/uZIRcTolXUSco7ISoPVU9GUfF/DbdfnUuD1cPi6BxddOIOs3qKhExHEqqxC3rbCa3//9ICu2FeNu9nJt5lC+OesCJo/Qkh4i4j9UViGosq6Jt7Yf4c+bCvms6Bj9IlzcMiWFey8dw4VDY5yOJyJyGpVViDjW0Mz7u0tYsbWYtfvK8Xgt44fG8G/zM7l5Sgqx0RFORxQR6ZLKKkhZazlQVsv7u0t5f3cpmwqq8HgtKYP7cf/laczPHs6E5BjNPCEiAUFlFSRay6mODfmVbMivYEN+JcXH3ABMSI5h4RVpXJ2exJQRQwgLU0GJSGBRWQUgay2HqxvYcbiGncXH2HH4GJ8VHaOirgmAhIFRzEiL45tp8cyekETK4H4OJxYROT8qKz/mafFSWNVAXlkteWV15JXXcqCsjr0lx6mubwbAFWa4IHEgV45PZProOGakxTM6vr+G90QkqKisHNLitVTXN1Fe20TxsQaKqxs4Uu2m+Njn74urG2husScfM6R/BGmJA5mXmczElEFMTBnEhOQYzSohIkHPp7IyxswDfg64gF9ba5/ssN20bb8eqAfusdZu6eGsfsdaS31TC7WNHo67mznu9lDb6KHW7eF4o4fjbg/HGpqprGuksq6Jitqm1vd1TVTXN+G1pz5fmIGhsdEMH9yPSSmDuH7SMNISBpCWOIC0hIEMGRDpzBcqIuKwbsvKGOMClgDXAEXARmPMCmvtrna7XQeMa3ubATzf9r5X7S+tpbbRQ4vXi6fF4vG2vbV4295bPG3bWryWZq+Xlvb3n/zY0uTx0uhpwd3spbG5hUaPF3dzC25PC43NXtxt29wntjW1UNfkOa1wOjIGBveLIG5AJPEDohibOJCLxkQSPyCy9b6BUQwf1FpQSTFRhLvCevufTUQk4PhyZDUd2G+tzQMwxiwH5gPty2o+8LK11gLrjTGDjTHDrLVHejxxO9/546fsOFzTI88V4TJEh7uIinARHRFGVHgY0REuoiNcRIWHEdsvou1+18n3MdHhDIwKZ2B0ODHREcS0fTwwKpyY6HBioiIYEOVSAYmInCdfyioFKGx3u4jTj5o62ycFOKWsjDELgYUA0dHR5OTknG3eU9Q1eoiyFjAYA6b1c2AADJh293e8bdpuc2J7O+62t2NAWVkZiYmJ55Wzvc3Fm3vsufqCK8xFzhvn9/90Jj3979sXAi2z8vauQMsL/p158+bNb1tr53W835ey6uy0so6DX77sg7X2BeAFgJycHLtp0yYfPr2zcnJyCIScJyhv7wu0zMrbuwItL/h95tOKCsCX8akiYES726lA8TnsIyIick58KauNwDhjzBhjTCRwJ7Ciwz4rgK+aVjOBY739epWIiISObocBrbUeY8wDwDu0nrr+orV2pzFmUdv2pcBKWk9b30/rqev39l7kvrVw4UKnI5wV5e19gZZZeXtXoOWFwMxsWk/g63uB8pqViIj0qU6n39E51SIi4vdUViIi4vdUVj76xS9+wfjx48nMzGTx4sVOx/HJM888gzGG8vJyp6Oc0fe+9z0mTJhAVlYWt9xyC9XV1U5H6tTbb7/N+PHjueCCC3jyySe7f4CDCgsLmT17Nunp6WRmZvLzn//c6Ug+aWlpYcqUKdx4441OR/FJdXU1CxYsYMKECaSnp/PJJ584HemMfvrTn5KZmcnEiRO56667cLvdTkfynbXWkbdp06bZQPH+++/bq6++2rrdbmuttSUlJQ4n6t6hQ4fs3Llz7ciRI21ZWZnTcc7onXfesc3NzdZaaxcvXmwXL17scKLTeTwem5aWZg8cOGAbGxttVlaW3blzp9OxulRcXGw3b95srbW2pqbGjhs3zq/znvDss8/au+66y95www1OR/HJV7/6VfurX/3KWmttY2OjraqqcjhR14qKiuzo0aNtfX29tdba22+/3S5btszZUJ3rtDN0ZOWD559/nkceeYSoqCgAkpKSHE7UvYceeoinn346IJYKmTt3LuHhrSemzpw5k6KiIocTnW7Dhg1ccMEFpKWlERkZyZ133snrr7/udKwuDRs2jKlTpwIQExNDeno6hw8fdjjVmRUVFfHWW29x//33Ox3FJzU1NXz44Yfcd999AERGRjJ48GCHU52Zx+OhoaEBj8dDfX09w4cPdzqSz1RWPti7dy9r165lxowZXHnllWzcuNHpSGe0YsUKUlJSmDx5stNRztqLL77Idddd53SM0xw+fJgRIz6/7j01NdXvf/mfUFBQwKeffsqMGb0+t/R5+c53vsPTTz9NWFhg/FrKy8sjMTGRe++9lylTpnD//fdTV1fndKwupaSk8N3vfpeRI0cybNgwBg0axNy5c52O5TOtZ9Vmzpw5HD169LT7f/SjH+HxeKiqqmL9+vVs3LiRL37xi+Tl5Tl61HKmvD/+8Y959913HUjVtTPlnT9//smPw8PDufvuu/s6XrdsJ5d4BMJRa21tLbfddhs/+9nPiI2NdTpOl958802SkpKYNm0aa9ascTqOTzweD1u2bOEXv/gFM2bM4MEHH+TJJ5/kiSeecDpap6qqqnj99dfJz89n8ODB3H777fzud7/jy1/+stPRfKKyarNq1aoutz3//PPceuutGGOYPn06YWFhlJeXOzoRZFd5t2/fTn5+/smjqqKiIqZOncqGDRtITk7uy4inONO/L8BLL73Em2++yerVq/2yBFJTUyks/Hyu5qKiIr8fQmlubua2227j7rvv5tZbb3U6zhmtW7eOFStWsHLlStxuNzU1NXz5y1/md7/7ndPRupSamkpqaurJI9YFCxb49Yk3q1atYsyYMSd/b9166618/PHHAVNWgXG87bCbb76Z999/H2gdEmxqaiIhIcHhVJ2bNGkSpaWlFBQUUFBQQGpqKlu2bHG0qLrz9ttv89RTT7FixQr69+/vdJxOXXTRRezbt4/8/HyamppYvnw5N910k9OxumSt5b777iM9PZ2HH37Y6Tjd+vd//3eKioooKChg+fLlXHXVVX5dVADJycmMGDGCPXv2ALB69WoyMjIcTtW1kSNHsn79eurr67HWsnr1atLT052O5TPHZrAwxnQ6Dbw/apsT8UUgG2gCvmutfd/ZVL4xxhQAOdZavz1/3RizH4gCKtruWm+tXeRgpE4ZY64Hfsbn0479yOFIXTLGXAasBbYD3ra7f2CtXelcKt8YY2bR+jPm9+evG2OygV8DkUAecK+1tsrZVF0zxjwO3AF4gE+B+621jc6m8o1jZSUiIuIrDQOKiIjfU1mJiIjfU1mJiIjfU1mJiIjfU1mJiIjfU1mJiIjfU1mJiIjf+/+zHbIJGiiJjQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# nbi:hide_in\n", "# Compute the empirical cdf and compare against the true cdf\n", "\n", "# plot setup\n", "plt.figure(figsize=(7,6)) \n", "plt.gca().spines['top'].set_visible(False)\n", "plt.gca().spines['right'].set_visible(False)\n", "plt.gca().spines['bottom'].set_position(('data',0))\n", "\n", "def ecdf(x):\n", " erange_x, counts = np.unique(x, return_counts=True)\n", " cdf_emp = np.cumsum(counts)/x.size\n", " return cdf_emp, erange_x\n", "\n", "# add padding \n", "def padding(cdf_values, range_x):\n", " edge = (range_x[0]-2, range_x[-1]+2)\n", " range_x_padded = np.pad(range_x, (1,1), 'constant', constant_values=edge)\n", " cdf_values_padded = np.pad(cdf_values, (1,1), 'constant', constant_values=(0, 1))\n", " return cdf_values_padded, range_x_padded\n", "\n", "# sort the data\n", "x = np.sort(x)\n", "\n", "# compute empirical cdf\n", "ecdf_values, erange_x = ecdf(x)\n", "ecdf_values_padded, erange_x_padded = padding(ecdf_values, erange_x)\n", "\n", "# the true cdf\n", "N=100\n", "delta = 1/N\n", "xval = np.linspace(np.min(x)-4, np.max(x)+4, 1000)\n", " \n", "# plot\n", "plt.plot(xval, cdf_func(xval), label ='True cdf')\n", "plt.step(erange_x_padded, ecdf_values_padded, where='post', \n", " color=\"green\", linewidth=1, label='Empirical cdf')\n", "plt.legend()\n", "plt.show();\n", "\n" ] } ], "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 }