segasai/untitled8.ipynb Secret

## untitled8.ipynb
{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "Untitled8.ipynb",
      "version": "0.3.2",
      "provenance": [],
      "collapsed_sections": [],
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/gist/segasai/2a5bbd6f5bedbac9e7d2c841d0796de9/untitled8.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "16lCF-ZPfCUE",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "import numpy as np\n",
        "import scipy.stats\n",
        "\n",
        "\n",
        "class LogNorm:\n",
        "    def __init__(self, mu, s):\n",
        "        self.m1 = 0\n",
        "        self.m2 = np.inf\n",
        "        self.d = scipy.stats.lognorm(s=s, scale=mu)\n",
        "\n",
        "    def pdf(self, x):\n",
        "        return self.d.pdf(x)\n",
        "\n",
        "    def cdf(self, x):\n",
        "        return self.d.cdf(x)\n",
        "\n",
        "\n",
        "class TruncLogNorm:\n",
        "    def __init__(self, mu, s, m1, m2):\n",
        "        self.m1 = m1\n",
        "        self.m2 = m2\n",
        "        self.d = scipy.stats.lognorm(s=s, scale=mu)\n",
        "        self.norm = self.d.cdf(self.m2) - self.d.cdf(self.m1)\n",
        "\n",
        "    def pdf(self, x):\n",
        "        return self.d.pdf(x) * (x >= self.m1) * (x <= self.m2) / self.norm\n",
        "\n",
        "    def cdf(self, x):\n",
        "        return (self.d.cdf(x) - self.d.cdf(self.m1)) / self.norm\n",
        "\n",
        "\n",
        "class PowerL:\n",
        "    def __init__(self, slope, m1, m2):\n",
        "        self.slope = slope\n",
        "        self.m1 = m1\n",
        "        self.m2 = m2\n",
        "\n",
        "    def pdf(self, x):\n",
        "        return x**self.slope * (self.slope + 1) / (\n",
        "            self.m2**(self.slope + 1) -\n",
        "            self.m1**(self.slope + 1)) * (x >= self.m1) * (x <= self.m2)\n",
        "\n",
        "    def cdf(self, x):\n",
        "        return (np.clip(x, self.m1, self.m2)**(self.slope + 1) -\n",
        "                (self.m1**(self.slope + 1))) / (self.m2**(self.slope + 1) -\n",
        "                                                self.m1**(self.slope + 1))\n",
        "\n",
        "\n",
        "class BrokenPowerL:\n",
        "    def __init__(self, slopes, breaks):\n",
        "        assert (len(slopes) == len(breaks) - 1)\n",
        "        nsegm = len(slopes)\n",
        "        pows = []\n",
        "        for i in range(nsegm):\n",
        "            pows.append(PowerL(slopes[i], breaks[i], breaks[i + 1]))\n",
        "        weights = [1]\n",
        "        for i in range(1, nsegm):\n",
        "            rat = pows[i].pdf(breaks[i]) / pows[i - 1].pdf(breaks[i])\n",
        "            weights.append(weights[-1] / rat)\n",
        "        weights = np.array(weights)\n",
        "        self.slopes = slopes\n",
        "        self.breaks = breaks\n",
        "        self.pows = pows\n",
        "        self.weights = weights / np.sum(weights)\n",
        "        self.nsegm = nsegm\n",
        "        self.m1 = breaks[0]\n",
        "        self.m2 = breaks[-1]\n",
        "\n",
        "    def pdf(self, x):\n",
        "        x1 = np.asarray(x)\n",
        "        ret = x1 * 0.\n",
        "        for i in range(self.nsegm):\n",
        "            xind = (x1 < self.breaks[i + 1]) & (x1 > self.breaks[i])\n",
        "            if xind.sum() > 0:\n",
        "                ret[xind] = self.weights[i] * self.pows[i].pdf(x1[xind])\n",
        "        return ret\n",
        "\n",
        "    def cdf(self, x):\n",
        "        x1 = np.asarray(x)\n",
        "        ret = x1 * 0.\n",
        "        cums = np.r_[[0], np.cumsum(self.weights)]\n",
        "        for i in range(self.nsegm):\n",
        "            xind = (x1 < self.breaks[i + 1]) & (x1 > self.breaks[i])\n",
        "            if xind.sum() > 0:\n",
        "                ret[xind] = cums[i] + self.weights[i] * self.pows[i].cdf(\n",
        "                    x1[xind])\n",
        "        return ret\n",
        "\n",
        "\n",
        "class CombinedDistr:\n",
        "    def __init__(self, distrs):\n",
        "        nsegm = len(distrs)\n",
        "        self.distrs = distrs\n",
        "        weights = [1]\n",
        "        breaks = [_.m1 for _ in self.distrs] + [self.distrs[-1].m2]\n",
        "        for i in range(nsegm-1):\n",
        "            assert(self.distrs[i].m2==self.distrs[i+1].m1)\n",
        "        for i in range(1, nsegm):\n",
        "            rat = distrs[i].pdf(breaks[i]) / distrs[i - 1].pdf(breaks[i])\n",
        "            weights.append(weights[-1] / rat)\n",
        "        weights = np.array(weights)\n",
        "        self.breaks = breaks\n",
        "        self.weights = weights / np.sum(weights)\n",
        "        self.nsegm = nsegm\n",
        "\n",
        "    def pdf(self, x):\n",
        "        x1 = np.asarray(x)\n",
        "        ret = x1 * 0.\n",
        "        for i in range(self.nsegm):\n",
        "            xind = (x1 < self.breaks[i + 1]) & (x1 > self.breaks[i])\n",
        "            if xind.sum() > 0:\n",
        "                ret[xind] = self.weights[i] * self.distrs[i].pdf(x1[xind])\n",
        "        return ret\n",
        "\n",
        "    def cdf(self, x):\n",
        "        x1 = np.asarray(x)\n",
        "        ret = x1 * 0.\n",
        "        cums = np.r_[[0], np.cumsum(self.weights)]\n",
        "        for i in range(self.nsegm):\n",
        "            xind = (x1 < self.breaks[i + 1]) & (x1 > self.breaks[i])\n",
        "            if xind.sum() > 0:\n",
        "                ret[xind] = cums[i] + self.weights[i] * self.distrs[i].cdf(\n",
        "                    x1[xind])\n",
        "        return ret"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "8ku2dVXkfhco",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 286
        },
        "outputId": "b5ae960f-9ee6-48ce-dc4b-077b2b8d683d"
      },
      "source": [
        "import matplotlib.pyplot as plt\n",
        "xgrid=np.linspace(0.100001,9.9999,100000)\n",
        "plt.plot(xgrid,CombinedDistr([PowerL(-1.1,0.1,2),PowerL(-2.1,2,4),TruncLogNorm(5,.1,4,10)]).pdf(xgrid));\n",
        "plt.plot(xgrid,CombinedDistr([PowerL(-1.1,0.1,2),PowerL(-2.1,2,4),TruncLogNorm(5,.1,4,10)]).cdf(xgrid))"
      ],
      "execution_count": 7,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "[<matplotlib.lines.Line2D at 0x7f420b02ef60>]"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 7
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xt81PWd7/HXJzOT+41cCBAIQUAQ\nRLwAotbWavVYq2Jba+3F2p5tOdV2bbvt7lb3bG17tqfddbdnq25rrXVrW1er1lrXUqyr1MsqKggi\ndwIiBAIJScj9OvM9f/wGDJBJAkwy+f3yfj4e85jf/OY7v/mMhne++f6+8/2Zcw4REQmWtFQXICIi\nyadwFxEJIIW7iEgAKdxFRAJI4S4iEkAKdxGRAFK4i4gEkMJdRCSAFO4iIgEUTtUbl5SUuMrKylS9\nvYiIL61evfqAc650sHYpC/fKykpWrVqVqrcXEfElM3tnKO00LCMiEkAKdxGRAFK4i4gEkMJdRCSA\nFO4iIgGkcBcRCaBBw93MppjZCjPbaGYbzOwr/bS5yMyazGxt/Pat4SlXRESGYijz3HuBrzvn3jCz\nPGC1mT3jnNt4VLsXnXNXJr/EI23Z18JT6/Zy4/mVlORmDPfbiYj40qA9d+dcjXPujfh2C7AJKB/u\nwhLZXtfKXc9VUd/anaoSRERGveMaczezSuAs4NV+nj7PzN40sz+a2dwEr19qZqvMbFVdXd1xFwsQ\nSjMAeqKxE3q9iMhYMORwN7Nc4LfAV51zzUc9/QYw1Tk3H7gLeKK/Yzjn7nXOLXDOLSgtHXRphH6F\n4+EejbkTer2IyFgwpHA3swhesD/onHv86Oedc83Oudb49jIgYmYlSa00LhzySu5VuIuIJDSU2TIG\n/BzY5Jz7YYI2E+LtMLNF8ePWJ7PQQw713Hs1LCMiktBQZstcANwAvGVma+P7bgMqAJxz9wDXAjeZ\nWS/QAVzvnBuWrrWGZUREBjdouDvnXgJskDZ3A3cnq6iBhEPxE6oKdxGRhHz3DdVQmldyNKZhGRGR\nRHwX7u+OuavnLiKSiP/CPT4so9kyIiKJ+S/c0zQVUkRkMD4Md02FFBEZjO/C/dDyA+q5i4gk5rtw\nj4QOzZZRuIuIJOK7cA9pWEZEZFC+C/eIZsuIiAzKd+Ee0jx3EZFB+S7cI1oVUkRkUL4L99DhhcM0\n5i4ikojvwj18+EpM6rmLiCTiu3A3M0JppqmQIiID8F24gzc006NhGRGRhHwZ7pE0I6phGRGRhHwZ\n7qE002wZEZEB+DLcw6E0ejUsIyKSkD/DXSdURUQG5Ntw11RIEZHE/BnuoTT13EVEBuDPcNcJVRGR\nAfky3ENppiV/RUQG4Mtw92bLqOcuIpKIP8NdPXcRkQH5M9xDGnMXERmIP8Nd89xFRAbky3D3Tqgq\n3EVEEvFluEe0/ICIyIB8Ge5aOExEZGCDhruZTTGzFWa20cw2mNlX+mljZnanmVWZ2TozO3t4yvVE\nQml096rnLiKSSHgIbXqBrzvn3jCzPGC1mT3jnNvYp80HgZnx27nAT+L3wyI9lEaPpkKKiCQ0aM/d\nOVfjnHsjvt0CbALKj2q2BPil86wECs1sYtKrjUsPp2nhMBGRARzXmLuZVQJnAa8e9VQ5sLvP42qO\n/QWQNJGQaVhGRGQAQw53M8sFfgt81TnXfCJvZmZLzWyVma2qq6s7kUMA3pi7hmVERBIbUribWQQv\n2B90zj3eT5M9wJQ+jyfH9x3BOXevc26Bc25BaWnpidQLeMMy3Qp3EZGEhjJbxoCfA5uccz9M0OxJ\n4DPxWTOLgSbnXE0S6zxCumbLiIgMaCizZS4AbgDeMrO18X23ARUAzrl7gGXAFUAV0A58Lvmlvss7\noapwFxFJZNBwd869BNggbRzwpWQVNZhIKI2Yg2jMEUobsDQRkTHJl99QjYS8sjU0IyLSP1+Ge3o4\nHu4amhER6Zc/wz3kDcWo5y4i0j9fhvuhYRmdVBUR6Z8vw/3QsIzCXUSkf74Md51QFREZmL/DXT13\nEZF++TLcMw4Py2hlSBGR/vgy3DUsIyIyMJ+GuzcVUidURUT658tw15eYREQG5stw17CMiMjAfBnu\nGZrnLiIyIF+Gu3ruIiID82e4q+cuIjIgX4Z7+uEvMWmeu4hIf/wd7hqWERHply/DPRLWPHcRkYH4\nMtwP9dy7ehTuIiL98WW4h0NphNOMzt5oqksRERmVfBnuAJmREJ09CncRkf74ONzT6NSwjIhIv3wb\n7hnhEF3quYuI9Mu34Z4ZSdOYu4hIAj4O95CGZUREEvB5uKvnLiLSHx+He5rCXUQkAf+GezhEl5Yf\nEBHpl3/DXcMyIiIJ+TbcMzTPXUQkoUHD3czuN7NaM1uf4PmLzKzJzNbGb99KfpnHyoyE6NJUSBGR\nfoWH0OYXwN3ALwdo86Jz7sqkVDREmWFNhRQRSWTQnrtz7gWgYQRqOS6aLSMikliyxtzPM7M3zeyP\nZjY3ScccUFYkRG/M0as13UVEjpGMcH8DmOqcmw/cBTyRqKGZLTWzVWa2qq6u7qTeNDMSAqBT0yFF\nRI5x0uHunGt2zrXGt5cBETMrSdD2XufcAufcgtLS0pN638yIV7qGZkREjnXS4W5mE8zM4tuL4ses\nP9njDibjUM9d4S4icoxBZ8uY2UPARUCJmVUDtwMRAOfcPcC1wE1m1gt0ANc759ywVRx3eFhGM2ZE\nRI4xaLg75z4xyPN3402VHFGZYQ3LiIgk4ttvqGalez33DoW7iMgxfBvu2eneHx3t3Qp3EZGj+Tbc\nczK8nntbV2+KKxERGX38G+7xnrvCXUTkWP4N9wwNy4iIJOLbcM+On1BtVc9dROQYvg33jHAaoTSj\nvVvhLiJyNN+Gu5mRkx6irUvDMiIiR/NtuIM37q4TqiIix/J1uGenh3RCVUSkH74O99yMMG0acxcR\nOYavwz07XcMyIiL98XW452TohKqISH98He7Z6WFNhRQR6Yevwz0nI0yreu4iIsfwd7inh9RzFxHp\nh6/DPTsjTHt3lGhs2C/8JCLiK74O9/xMb/EwrS8jInIkf4d7VgSA5o6eFFciIjK6+DvcM71wb1K4\ni4gcwdfhXqCeu4hIv3wd7vlZ3ph7c6fCXUSkL1+H+6Geu4ZlRESO5Otwf/eEqmbLiIj05etwz00P\nk2bquYuIHM3X4Z6WZuRlRjTmLiJyFF+HO3jj7potIyJyJN+He35WWMMyIiJH8X24F2RFaO7UCVUR\nkb4CEe7quYuIHMn34V6YnU5jW3eqyxARGVUGDXczu9/Mas1sfYLnzczuNLMqM1tnZmcnv8zEinPS\naWzv1rK/IiJ9DKXn/gvg8gGe/yAwM35bCvzk5MsauuKcdGIODrar9y4ickh4sAbOuRfMrHKAJkuA\nXzrnHLDSzArNbKJzriZJNQ6oKDcDgPq2borj2yK+5hz0dEDnQehug95O6O3y7ns64487IdoDLgou\nBrH4/RHb0SO3+32vIe/06hrqQRK2FQAqFsP09w/rWwwa7kNQDuzu87g6vu+YcDezpXi9eyoqKpLw\n1lCSkw5AfWs3lCXlkCLDp6sVDmyFxrehuQaa90LzHmipgfYG6GzyQj2qv0QD7YKv+iLch8w5dy9w\nL8CCBQuS8qu9+HDPvSsZhxNJntY6qH7du9W86YV60+4j20SyIb8c8ifChNMhswAyCyGr0NtOz4Vw\nJkSyIJzhbR+6D0XAQmBpkBY6ajvt2P1mCQrtZ//xtE3UPuExZCQkI9z3AFP6PJ4c3zciiuI99wbN\nmJFU62qFt1+Aqmdg+wqvdw6QFobxp0HFeVB6I5ScCsUzvFDPLFAIyrBIRrg/CXzZzB4GzgWaRmq8\nHWBcdgQzONCqcJcU6GyGzX+A9Y/Bjuch1uP1tqe9Fxb+BUxeCBPnez1vkRE0aLib2UPARUCJmVUD\ntwMRAOfcPcAy4AqgCmgHPjdcxfYnHEqjMCtCfauGZWSEOAdvPw+r/h22LvdObhZUwOIvwszLYMpi\nCKenukoZ44YyW+YTgzzvgC8lraITUJyboWEZGX6dzfDmQ/D6fd74eVYRnP0ZmPcxr4eu4RUZRUb0\nhOpwKc5J54B67jJcOpvh1Xvg5buhqwnKz4Fr7oG5H4ZIZqqrE+lXIMK9LD+TtbsPproMCZrudnj9\nZ/DSv0JHA8z6EFz4dZh8TqorExlUIMJ9YkEmy9d34pzD9KexnKzeLlj9C3jhn6GtFmZ8AN5/m9dj\nF/GJQIT7hIJMuqMxGvQtVTkZ0R5Y+yA8fwc0V8PU98B1v4Sp56W6MpHjFohwn1jgjXvWNHUq3OX4\nxaLw1qPw5+9D407v5Og1/wbT3qeTpOJbgQj3CQXeHOJ9TZ2cXl6Q4mrEN2Ix2PiEF+oHtsKEM+CT\nj3jTGRXq4nOBCPfDPffmzhRXIr7gHGxZBiu+D/vfgtLZcN2vYPaVkOb7SxyIAAEJ95LcDEJpxr6m\njlSXIqOZc7Dlj15Pfd86KDoFPnIfnP4Rbw0WkQAJRLiH0oyyvAxqmtRzl37Eol6ov/BP3gJe46bB\nNT+BeddBKBD/BESOEZif7ImFWexpVM9d+uhq9Wa/rPyJt4jXuEpY8mM44+MKdQm8wPyEVxRl8+qO\n+lSXIanmHOx+zVsmYP3j3jdKJy+CS74Fp12tUJcxIzA/6VOLs3li7R46e6JkRjR+OqZEe70107c9\nDRt/Dw07vHXST7sKFi2FyQtSXaHIiAtMuFcW5+AcVDe2M2N8XqrLkeHU3e6NnVe/5oX6zpego9Fb\nN33qBXDhN2DO1ZChnwMZuwIT7hXF2QDsPKBw9yXnoKvFW8Olo9G7tR/aPugtA3BgG9RXHXk1o3GV\ncOoH4dTLYPrF3sUvRCQ44V5ZnAPAOw3tKa5kjIv2etcAPRTKh4J6oH2HQjzRRZwBMvK9qYsVi6Hk\nM1A21/smae74kftsIj4SmHAflx0hLzPMO/VtqS4lWGIxrzfdVvfurbXPdtsBr1fdVgdt9dDdMvDx\nMvLj1wcthKxx3uXnsoogu8h7nDXOe3xoO7vIa6uLX4gcl8CEu5lRWZzDznr13IestwtaaqB575G3\nlr7b+/rvUVsIckogp9S7H1cJ2cV9grnw3YA+FOSZBZqtIjJCAvUvbcb4XFZqOuS7nPN61I07oeFt\n777vrWXvsa9Jz4X8Sd7tlIsgbwLklsWDfHw8zEu9sNZX9UVGrUCF+6llefxuzR6a2nsoyI6kupyR\nE+2Fg+9A3eb4bat3f2Ab9Bw1TJU3CYqmwfT3Q2EFFEyGvImQX+4FemZ+aj6DiCRVoMJ99gRvlszW\n2hYWVhaluJph0t0G+9Z7UwFr1nr3B7ZCtM81ZPPLoeRUOPsGKJruDZmMq/TCXJeFExkTAhXus+Lh\nvnlfQMLdOa/3vetl2LUS9q7xgtzFvOezS2DifG8KYOlsKJ3lhbp63yJjXqDCfWJBJnmZYbbsa051\nKScm2gv73oR3XoFdr3iB3n7Aey67xJv6N+caL9AnnekNp2jdcRHpR6DC3cyYVZbHln2DTMcbTRp2\nwPbnYPsKePtFby0U8IZRZl7mzeueej4Uz1CQi8iQBSrcAeZOyufR1dVEY45Q2igMw+42L8i3P+uF\neuNOb3/BFJi7xLu029TzvZObIiInKHDhflbFOB545R227m/htImjZOy5tdZbT3zLMi/Yo13elMNp\n74XFX/LGzIunq2cuIkkTuHA/c0ohAGt2HUxtuDfvhfW/hY1Peotb4aCgAhZ8DmZd4fXOQ2NouqaI\njKjAhfvU4mzGZUdYu7uRT55bMbJv3t4Am56Etx7zVirEeRddvuhWmH0FlJ2u3rmIjIjAhbuZceaU\nQtbsOjgybxiLwY7nYPUvYMtyiPV4Jz8v+iacfi2UzBiZOkRE+ghcuAMsqCxixZYt1Ld2UZybMTxv\n0rIP1vwa3ngADu7y1lVZtBTOuM6bqqgeuoikUCDD/YIZJdzx9BZe3l7PVfOTPOtk31vw8l3eeHqs\n1zsp+oFvw+wrITxMv0hERI7TkMLdzC4HfgSEgPuccz846vnPAncAe+K77nbO3ZfEOo/LvPIC8jLD\n/HfVgeSEu3PetMWX74IdKyCS4/XSF37em+UiIjLKDBruZhYC/g24FKgGXjezJ51zG49q+hvn3JeH\nocbjFkozzjulmBe3HcA5h53MEMmO5+G5f/Au6ZY7AS653ZvxkjUueQWLiCTZUHrui4Aq59wOADN7\nGFgCHB3uo8p7Zpbwp4372XGgjemlucd/gOrV8Oy34e0XvJUUP/RDOOvTGnoREV8YyoLc5UCfi1ZS\nHd93tI+a2Toze8zMpiSlupNwyWllADy9Yd/xvbBlP/zuJrjvYqjdBJf/AG5ZAwv/QsEuIr6RrKst\n/CdQ6Zw7A3gGeKC/Rma21MxWmdmqurq6JL11/8oLs5g/uYDl64cY7rEovPJjuHsBvPUoXPBVL9QX\n36RlckXEd4YS7nuAvj3xybx74hQA51y9c64r/vA+4Jz+DuScu9c5t8A5t6C0tPRE6j0ul58+kXXV\nTVQ3DnLpvQNVcP/l8PStMGUR3LwSLv0OZOQNe40iIsNhKOH+OjDTzKaZWTpwPfBk3wZmNrHPw6uB\nTckr8cRdMW8CAE+tq+m/gXOw8h645wJvnfSP/Aw+9Zi+eCQivjdouDvneoEvA0/jhfYjzrkNZvZd\nM7s63uwWM9tgZm8CtwCfHa6Cj8fU4hwWVo7jkdd345w78snOZnjkBlj+t95KjDev9L6ApC8fiUgA\nDGmeu3NuGbDsqH3f6rN9K3BrcktLjo8vrOAbj77Ja283cO4pxd7O2k3w8Ke85XYv+x6c9yWFuogE\nSuAvX/+heRPJywjzH6/t8na8/QL8/H9Adyt89ik4/8sKdhEJnMCHe1Z6iI8tmMIf1tVQv/I/4Ncf\nhfyJ8PlnvWV3RUQCKPDhDvD5C6exJO1FipbfDOUL4H8uh8KUT8UXERk2YyLcJ1Uv447wPbzq5lB7\nzYNaOkBEAi/44b59Bfz2C3RPWsQXer7B//vznsFfIyLic8EO97ot8MiNUDqLzM88ysfOm81vXt/F\n5n3Nqa5MRGRYBTfcOw7Cf3wcwunwyd9AZj63XDKDvMwIf//EemIxN/gxRER8Kpjh7hw89TXvCkkf\nfxAKvWupFman8/dXzuH1nY38+8s7U1ujiMgwCma4r/k1bHgcLv47qDj3iKc+enY5l8wezz8t38zW\n/S0pKlBEZHgFL9yba2D5rVB5obey41HMjO9/ZB75WRH+169W09TRk4IiRUSGV/DC/elbIdoNV98J\naaF+m4zPz+THnzqb3Q3tfPXhNfREYyNcpIjI8ApWuFc9Cxt+B+/9ayg6ZcCmCyuL+O6S01mxpY6/\neWydTrCKSKAMaeEwX4jF4L++DYVT4YJbhvSST55bQWN7N3c8vYWMcBrf+/A8QmlaZ0ZE/C844b7p\nSdi3Dj780+O6HN7NF02nqyfKnc9V0dDWzZ2fOIvMSP/DOSIifhGMYZlYDP78fSidDfM+dlwvNTP+\n6rJZfOfquTyzaT/X/fQVdjcMcuUmEZFRLhjhvv05qNsM7/lawpOog7nx/Ep++ulzeLuujSvveok/\nHe+FtUVERpFgDMu8+hPILYO5Hzmpw1w2dwJP3ZLHzQ++wdJfreaq+ZO4/ao5lOQOfZhH/KOpo4fN\nNc3UtXbR3hUlMz1ESW46p5bl6f+5+J7/w/1AFVT9F7z/f3tLDZykqcU5PH7z+dzz5x3cvWIbL26r\n4y8vnsmnF1eQEdZYvN9tr2vl92v3snx9DVv3tyZsN6Uoi4tOHc8V8yay+JQiTBd0EZ+xY64tOkIW\nLFjgVq1adfIHeu4f4MV/gb/aBHkTTv54fWzb38J3/nMjL1UdoLwwi698YCZLzpykkPehN3Y18uMV\n2/mvTfsxg3OnFfGeGSXMm1zIxIJMstNDdPZEqWnqZMu+FlbuqOfl7fW0d0epLM7m04un8slzK8hO\n939/SPzNzFY75xYM2s7X4e4c/Gg+FE+HG36XnML68dK2A/zj8s28taeJ8XkZ3Hh+JZ86t4LC7JP/\nS0GG196DHXxv2Sb+sK6GwuwInz2/kk8uqmB8fuagr+3sibLsrRoeem0Xr+9spDgnnc9feAqfPb+S\nrHT9gpfUGBvhvutVuP8yb/rj/OuTU1gCzjle3HaAn724gxe3HSA9lMalc8q49pzJXDizhHAoGOem\ngyIWc/z8pbf5l2e24BzcdNF0lr73lBPuea9+p4E7n63i+a11TMjP5G8/OIsl88tJ0/ciZISNjXB/\n+u/gtXvhb3ZARl5yChuCzfuaefi13fx+7R4a23soyU3n0jllXDZnAudNL9Y8+RSrbmznG4++ycod\nDXzgtPF8++q5TB6XnZRjv76zgf/z1EbWVTcxf3IB37pqDudMLUrKsUWGYmyE+90LIb8cPvNEcoo6\nTt29MZ7bXMt/vrmXP2+ppa07SnZ6iAtmlHD+9GLOn17CqWW5Ohk3Qpxz/G7NHm7//QZiznH71XP5\n2DmTk/7fPxZzPLF2D/+4fDP7m7u4ev4kvvnB2UwqzErq+4j0J/jh3vA23HkmXP4DWHxT8go7QV29\nUV7eXs8zG/fz4rY6djd0AFCck865pxQxf3IhZ0wuZN7kAnIzdFIu2Zrae7jtibf4w7oaFlaO44fX\nncmUouT01hNp7+7lnud38NPnt2MGN71vBkvfe4rG42VYBT/cX/sZLPsG/OUb3gnVUWZ3Qzuv7Khn\n5fZ6XtvZQHWjF/ZmML00lzPKC5hZlsepZbmcWpZHeWGWxm9P0GtvN/C136xlf3MnX7v0VL74vukj\nukZQdWM73//jZv6wroZJBZl884rTuOqMifqLTYZF8MP9kRthz2r42vrkFTWM6lu7WLeniXW7m3iz\n+iAb9jaxv7nr8PPZ6SFmjM9l5vg8Zk3IZWZZHrPK8phYkKmQSKA3GuPO56q4+7ltTCnK5kfXn8WZ\nUwpTVs+rO+r57lMb2bC3mQVTx3H7VXOZN7kgZfVIMAU73J2DH54GUy+Aa3+e3MJGUFN7D9tqW9i6\nv5VttS1s29/K1v0t1La8G/p5GWFmluUyr7yAMyYXMn9KAaeU5I75Xn5VbQt//dg61uw6yEfPnsx3\nlswdFcNd0Zjj0VW7uePpLTS0d/PRsydzy8UzqSge3iEiGTuCHe6N78CPzoAr/hkWfSG5hY0CjW3d\nbN3fwtbaVrbua2HzvmY27G2mvTsKQG5GmDkT85lZlsvM8V4vf+b4XErzMgLfy++Jxvjp89u589kq\nsjNCfOfquSw5szzVZR2jubOHu57dxgOvvEM05lgyfxI3v386M8aP3KwuCaZgh/u6R+DxL8AXX4IJ\n85Jb2CgVjTmqaltZV32QddVNbNjbRFVtK82dvYfb5GeGmVqcw5SiLCaPy2bKuPh9URaTCrN8/e1K\n5xzPb63j/y7bxNb9rXzojIl8+6q5lOaN7jVg9jd38rMXdvDgq7vo7I3yvlNL+fS5U3n/7PG6doCc\nkGCH+/LbYNX9cNueE14FMgicc9S1dFFV28q2Wm9oZ1dDB9WN7VQ3dtDde+TlA3MzwpTmZVCal8H4\nvAzG52UyPt/bHpeTTmFWhHHZ6YzLTicvMzwqhn6cc6x6p5E7n93Gi9sOUFmczd99aA6XzilLdWnH\npb61iwdeeYeHX9tFbUsXkwoyWXJWOVecPpHTy/MD/xeXJE9Sw93MLgd+BISA+5xzPzjq+Qzgl8A5\nQD3wcefczoGOeVLh/sDV0NUCS1ec2OvHgFjMUdfaRXVjO7sbOtjb1EFtcxd1rV3UNXdR29JJbUvX\n4aGeo6UZFGRFKMxOpzA7QmFWhPysCDkZYfIywuTEb7kZIXIzIuRkhMjNCJObGSYnPUxuRpjMSIiM\ncNoJ/ZJo7+7lTxv288ArO1mz6yDjsiN8+eKZ3LB4Kulh/34buCca49lNtTz02i5eqjpANOaYUpTF\nJbPLOG96MYunFVOQHUl1mTKKJS3czSwEbAUuBaqB14FPOOc29mlzM3CGc+6LZnY98GHn3McHOu4J\nh7tzcMd0mHUFLLn7+F8vR2jt6qWupYuGtm6aOrppbOvhYEcPTe3dNLZ72wfbu2ls76als5e2rl5a\nu3rp7Bn6RcXTw2lkRUJkRtLIjITIDPfZ7mf/gdZuVmyppb07ytTibD7/nmlce86UwM0fb2zr5pmN\n+1m2voaVO+rp7IlhBrPK8ji9vIC5k/KZMzGfaaU5lOYG/3yKDM1Qw30og7CLgCrn3I74gR8GlgAb\n+7RZAnw7vv0YcLeZmRuOMZ/WWmivh7LTk37osSg3w+tlTyvJOa7X9UZjtHVFae32Av9Q8Ld19dLS\n1Ut7Vy+dvTE6uqN09kbp6onR2ROloydKZ0+Uzvjjgx09dDZ5bQ7tz80Is+TMcq45cxILK4tGxfDQ\ncBiXk851C6dw3cIpdPVGeXN3E/9ddYA1uw/y5y21PLa6+nDbzEgaFUXZTBmXTUluBkW56RTnpFOU\n4w2jZaWHyEkPk5UeIjt+y4yECKcZoTTTL4YxaCjhXg7s7vO4Gjg3URvnXK+ZNQHFwIFkFHmE/fF5\n7WVzk35oGbpwKI2C7DQNISRJRjjEomlFLJr27jo1tc2dbKxpZldDO7vq29nV0M7uxg7e2tNEQ1s3\nvbGh950Ohfyh+0gojdCh4Idjwv/QQzPwWvTZ16f94VfZu9v6RTK46xdO4fMXnjKs7zGi0yfMbCmw\nFKCiouLEDhLJ9oZkFO4ScOPzMxMuTeyco7mjl/q2Lg529NDRHaW9O0p7d+/h7Y6eKNGYozfmiMZi\n9MYcvVEX3xfz7qMOd/iY8XscHN5+9/0OPXb9PHf410xq5mf4zkhc6Wso4b4HmNLn8eT4vv7aVJtZ\nGCjAO7F6BOfcvcC94I25n0jBTD3Pu4mMYWZGQXZEfzlJQkOZdvA6MNPMpplZOnA98ORRbZ4Eboxv\nXws8Nyzj7SIiMiSD9tzjY+hfBp7Gmwp5v3Nug5l9F1jlnHsS+DnwKzOrAhrwfgGIiEiKDGnM3Tm3\nDFh21L5v9dnuBD6W3NJERORE+ffbICIikpDCXUQkgBTuIiIBpHAXEQkghbuISAClbMlfM6sD3jmO\nl5QwHMsZjH5j8XOPxc8MY/NmmT3wAAADQElEQVRzj8XPDCf3uac650oHa5SycD9eZrZqKCuhBc1Y\n/Nxj8TPD2PzcY/Ezw8h8bg3LiIgEkMJdRCSA/BTu96a6gBQZi597LH5mGJufeyx+ZhiBz+2bMXcR\nERk6P/XcRURkiHwR7mZ2uZltMbMqM/tmqusZbmY2xcxWmNlGM9tgZl9JdU0jycxCZrbGzJ5KdS0j\nwcwKzewxM9tsZpvMbExcsMDMvhb/+V5vZg+ZWf9XJvE5M7vfzGrNbH2ffUVm9oyZbYvfj0v2+476\ncI9foPvfgA8Cc4BPmNmc1FY17HqBrzvn5gCLgS+Ngc/c11eATakuYgT9CFjunJsNzGcMfHYzKwdu\nARY4507HW048qEuF/wK4/Kh93wSedc7NBJ6NP06qUR/u9LlAt3OuGzh0ge7Acs7VOOfeiG+34P1j\nL09tVSPDzCYDHwLuS3UtI8HMCoD34l0TAedct3PuYGqrGjFhICt+9bZsYG+K6xkWzrkX8K5z0dcS\n4IH49gPANcl+Xz+Ee38X6B4TQQdgZpXAWcCrqa1kxPwr8DdALNWFjJBpQB3w7/GhqPvMLCfVRQ03\n59we4J+BXUAN0OSc+1NqqxpRZc65mvj2PqAs2W/gh3Afs8wsF/gt8FXnXHOq6xluZnYlUOucW53q\nWkZQGDgb+Ilz7iygjWH4E320iY8xL8H75TYJyDGzT6e2qtSIX5I06dMW/RDuQ7lAd+CYWQQv2B90\nzj2e6npGyAXA1Wa2E2/47WIz+3VqSxp21UC1c+7QX2aP4YV90H0AeNs5V+ec6wEeB85PcU0jab+Z\nTQSI39cm+w38EO5DuUB3oJiZ4Y3BbnLO/TDV9YwU59ytzrnJzrlKvP/PzznnAt2bc87tA3ab2az4\nrkuAjSksaaTsAhabWXb85/0SxsCJ5D6eBG6Mb98I/D7ZbzCka6imUqILdKe4rOF2AXAD8JaZrY3v\nuy1+LVsJnr8EHox3XnYAn0txPcPOOfeqmT0GvIE3O2wNAf22qpk9BFwElJhZNXA78APgETP7C7zV\nca9L+vvqG6oiIsHjh2EZERE5Tgp3EZEAUriLiASQwl1EJIAU7iIiAaRwFxEJIIW7iEgAKdxFRALo\n/wMOLzL3Fkf0MgAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": []
          }
        }
      ]
    }
  ]
}
	{
	"nbformat": 4,
	"nbformat_minor": 0,
	"metadata": {
	"colab": {
	"name": "Untitled8.ipynb",
	"version": "0.3.2",
	"provenance": [],
	"collapsed_sections": [],
	"include_colab_link": true
	},
	"kernelspec": {
	"name": "python3",
	"display_name": "Python 3"
	}
	},
	"cells": [
	{
	"cell_type": "markdown",
	"metadata": {
	"id": "view-in-github",
	"colab_type": "text"
	},
	"source": [
	"<a href=\"https://colab.research.google.com/gist/segasai/2a5bbd6f5bedbac9e7d2c841d0796de9/untitled8.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
	]
	},
	{
	"cell_type": "code",
	"metadata": {
	"id": "16lCF-ZPfCUE",
	"colab_type": "code",
	"colab": {}
	},
	"source": [
	"import numpy as np\n",
	"import scipy.stats\n",
	"\n",
	"\n",
	"class LogNorm:\n",
	" def __init__(self, mu, s):\n",
	" self.m1 = 0\n",
	" self.m2 = np.inf\n",
	" self.d = scipy.stats.lognorm(s=s, scale=mu)\n",
	"\n",
	" def pdf(self, x):\n",
	" return self.d.pdf(x)\n",
	"\n",
	" def cdf(self, x):\n",
	" return self.d.cdf(x)\n",
	"\n",
	"\n",
	"class TruncLogNorm:\n",
	" def __init__(self, mu, s, m1, m2):\n",
	" self.m1 = m1\n",
	" self.m2 = m2\n",
	" self.d = scipy.stats.lognorm(s=s, scale=mu)\n",
	" self.norm = self.d.cdf(self.m2) - self.d.cdf(self.m1)\n",
	"\n",
	" def pdf(self, x):\n",
	" return self.d.pdf(x) * (x >= self.m1) * (x <= self.m2) / self.norm\n",
	"\n",
	" def cdf(self, x):\n",
	" return (self.d.cdf(x) - self.d.cdf(self.m1)) / self.norm\n",
	"\n",
	"\n",
	"class PowerL:\n",
	" def __init__(self, slope, m1, m2):\n",
	" self.slope = slope\n",
	" self.m1 = m1\n",
	" self.m2 = m2\n",
	"\n",
	" def pdf(self, x):\n",
	" return x*self.slope (self.slope + 1) / (\n",
	" self.m2**(self.slope + 1) -\n",
	" self.m1*(self.slope + 1)) (x >= self.m1) * (x <= self.m2)\n",
	"\n",
	" def cdf(self, x):\n",
	" return (np.clip(x, self.m1, self.m2)**(self.slope + 1) -\n",
	" (self.m1(self.slope + 1))) / (self.m2(self.slope + 1) -\n",
	" self.m1**(self.slope + 1))\n",
	"\n",
	"\n",
	"class BrokenPowerL:\n",
	" def __init__(self, slopes, breaks):\n",
	" assert (len(slopes) == len(breaks) - 1)\n",
	" nsegm = len(slopes)\n",
	" pows = []\n",
	" for i in range(nsegm):\n",
	" pows.append(PowerL(slopes[i], breaks[i], breaks[i + 1]))\n",
	" weights = [1]\n",
	" for i in range(1, nsegm):\n",
	" rat = pows[i].pdf(breaks[i]) / pows[i - 1].pdf(breaks[i])\n",
	" weights.append(weights[-1] / rat)\n",
	" weights = np.array(weights)\n",
	" self.slopes = slopes\n",
	" self.breaks = breaks\n",
	" self.pows = pows\n",
	" self.weights = weights / np.sum(weights)\n",
	" self.nsegm = nsegm\n",
	" self.m1 = breaks[0]\n",
	" self.m2 = breaks[-1]\n",
	"\n",
	" def pdf(self, x):\n",
	" x1 = np.asarray(x)\n",
	" ret = x1 * 0.\n",
	" for i in range(self.nsegm):\n",
	" xind = (x1 < self.breaks[i + 1]) & (x1 > self.breaks[i])\n",
	" if xind.sum() > 0:\n",
	" ret[xind] = self.weights[i] * self.pows[i].pdf(x1[xind])\n",
	" return ret\n",
	"\n",
	" def cdf(self, x):\n",
	" x1 = np.asarray(x)\n",
	" ret = x1 * 0.\n",
	" cums = np.r_[[0], np.cumsum(self.weights)]\n",
	" for i in range(self.nsegm):\n",
	" xind = (x1 < self.breaks[i + 1]) & (x1 > self.breaks[i])\n",
	" if xind.sum() > 0:\n",
	" ret[xind] = cums[i] + self.weights[i] * self.pows[i].cdf(\n",
	" x1[xind])\n",
	" return ret\n",
	"\n",
	"\n",
	"class CombinedDistr:\n",
	" def __init__(self, distrs):\n",
	" nsegm = len(distrs)\n",
	" self.distrs = distrs\n",
	" weights = [1]\n",
	" breaks = [_.m1 for _ in self.distrs] + [self.distrs[-1].m2]\n",
	" for i in range(nsegm-1):\n",
	" assert(self.distrs[i].m2==self.distrs[i+1].m1)\n",
	" for i in range(1, nsegm):\n",
	" rat = distrs[i].pdf(breaks[i]) / distrs[i - 1].pdf(breaks[i])\n",
	" weights.append(weights[-1] / rat)\n",
	" weights = np.array(weights)\n",
	" self.breaks = breaks\n",
	" self.weights = weights / np.sum(weights)\n",
	" self.nsegm = nsegm\n",
	"\n",
	" def pdf(self, x):\n",
	" x1 = np.asarray(x)\n",
	" ret = x1 * 0.\n",
	" for i in range(self.nsegm):\n",
	" xind = (x1 < self.breaks[i + 1]) & (x1 > self.breaks[i])\n",
	" if xind.sum() > 0:\n",
	" ret[xind] = self.weights[i] * self.distrs[i].pdf(x1[xind])\n",
	" return ret\n",
	"\n",
	" def cdf(self, x):\n",
	" x1 = np.asarray(x)\n",
	" ret = x1 * 0.\n",
	" cums = np.r_[[0], np.cumsum(self.weights)]\n",
	" for i in range(self.nsegm):\n",
	" xind = (x1 < self.breaks[i + 1]) & (x1 > self.breaks[i])\n",
	" if xind.sum() > 0:\n",
	" ret[xind] = cums[i] + self.weights[i] * self.distrs[i].cdf(\n",
	" x1[xind])\n",
	" return ret"
	],
	"execution_count": 0,
	"outputs": []
	},
	{
	"cell_type": "code",
	"metadata": {
	"id": "8ku2dVXkfhco",
	"colab_type": "code",
	"colab": {
	"base_uri": "https://localhost:8080/",
	"height": 286
	},
	"outputId": "b5ae960f-9ee6-48ce-dc4b-077b2b8d683d"
	},
	"source": [
	"import matplotlib.pyplot as plt\n",
	"xgrid=np.linspace(0.100001,9.9999,100000)\n",
	"plt.plot(xgrid,CombinedDistr([PowerL(-1.1,0.1,2),PowerL(-2.1,2,4),TruncLogNorm(5,.1,4,10)]).pdf(xgrid));\n",
	"plt.plot(xgrid,CombinedDistr([PowerL(-1.1,0.1,2),PowerL(-2.1,2,4),TruncLogNorm(5,.1,4,10)]).cdf(xgrid))"
	],
	"execution_count": 7,
	"outputs": [
	{
	"output_type": "execute_result",
	"data": {
	"text/plain": [
	"[<matplotlib.lines.Line2D at 0x7f420b02ef60>]"
	]
	},
	"metadata": {
	"tags": []
	},
	"execution_count": 7
	},
	{
	"output_type": "display_data",
	"data": {
	"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xt81PWd7/HXJzOT+41cCBAIQUAQ\nRLwAotbWavVYq2Jba+3F2p5tOdV2bbvt7lb3bG17tqfddbdnq25rrXVrW1er1lrXUqyr1MsqKggi\ndwIiBAIJScj9OvM9f/wGDJBJAkwy+f3yfj4e85jf/OY7v/mMhne++f6+8/2Zcw4REQmWtFQXICIi\nyadwFxEJIIW7iEgAKdxFRAJI4S4iEkAKdxGRAFK4i4gEkMJdRCSAFO4iIgEUTtUbl5SUuMrKylS9\nvYiIL61evfqAc650sHYpC/fKykpWrVqVqrcXEfElM3tnKO00LCMiEkAKdxGRAFK4i4gEkMJdRCSA\nFO4iIgGkcBcRCaBBw93MppjZCjPbaGYbzOwr/bS5yMyazGxt/Pat4SlXRESGYijz3HuBrzvn3jCz\nPGC1mT3jnNt4VLsXnXNXJr/EI23Z18JT6/Zy4/mVlORmDPfbiYj40qA9d+dcjXPujfh2C7AJKB/u\nwhLZXtfKXc9VUd/anaoSRERGveMaczezSuAs4NV+nj7PzN40sz+a2dwEr19qZqvMbFVdXd1xFwsQ\nSjMAeqKxE3q9iMhYMORwN7Nc4LfAV51zzUc9/QYw1Tk3H7gLeKK/Yzjn7nXOLXDOLSgtHXRphH6F\n4+EejbkTer2IyFgwpHA3swhesD/onHv86Oedc83Oudb49jIgYmYlSa00LhzySu5VuIuIJDSU2TIG\n/BzY5Jz7YYI2E+LtMLNF8ePWJ7PQQw713Hs1LCMiktBQZstcANwAvGVma+P7bgMqAJxz9wDXAjeZ\nWS/QAVzvnBuWrrWGZUREBjdouDvnXgJskDZ3A3cnq6iBhEPxE6oKdxGRhHz3DdVQmldyNKZhGRGR\nRHwX7u+OuavnLiKSiP/CPT4so9kyIiKJ+S/c0zQVUkRkMD4Md02FFBEZjO/C/dDyA+q5i4gk5rtw\nj4QOzZZRuIuIJOK7cA9pWEZEZFC+C/eIZsuIiAzKd+Ee0jx3EZFB+S7cI1oVUkRkUL4L99DhhcM0\n5i4ikojvwj18+EpM6rmLiCTiu3A3M0JppqmQIiID8F24gzc006NhGRGRhHwZ7pE0I6phGRGRhHwZ\n7qE002wZEZEB+DLcw6E0ejUsIyKSkD/DXSdURUQG5Ntw11RIEZHE/BnuoTT13EVEBuDPcNcJVRGR\nAfky3ENppiV/RUQG4Mtw92bLqOcuIpKIP8NdPXcRkQH5M9xDGnMXERmIP8Nd89xFRAbky3D3Tqgq\n3EVEEvFluEe0/ICIyIB8Ge5aOExEZGCDhruZTTGzFWa20cw2mNlX+mljZnanmVWZ2TozO3t4yvVE\nQml096rnLiKSSHgIbXqBrzvn3jCzPGC1mT3jnNvYp80HgZnx27nAT+L3wyI9lEaPpkKKiCQ0aM/d\nOVfjnHsjvt0CbALKj2q2BPil86wECs1sYtKrjUsPp2nhMBGRARzXmLuZVQJnAa8e9VQ5sLvP42qO\n/QWQNJGQaVhGRGQAQw53M8sFfgt81TnXfCJvZmZLzWyVma2qq6s7kUMA3pi7hmVERBIbUribWQQv\n2B90zj3eT5M9wJQ+jyfH9x3BOXevc26Bc25BaWnpidQLeMMy3Qp3EZGEhjJbxoCfA5uccz9M0OxJ\n4DPxWTOLgSbnXE0S6zxCumbLiIgMaCizZS4AbgDeMrO18X23ARUAzrl7gGXAFUAV0A58Lvmlvss7\noapwFxFJZNBwd869BNggbRzwpWQVNZhIKI2Yg2jMEUobsDQRkTHJl99QjYS8sjU0IyLSP1+Ge3o4\nHu4amhER6Zc/wz3kDcWo5y4i0j9fhvuhYRmdVBUR6Z8vw/3QsIzCXUSkf74Md51QFREZmL/DXT13\nEZF++TLcMw4Py2hlSBGR/vgy3DUsIyIyMJ+GuzcVUidURUT658tw15eYREQG5stw17CMiMjAfBnu\nGZrnLiIyIF+Gu3ruIiID82e4q+cuIjIgX4Z7+uEvMWmeu4hIf/wd7hqWERHply/DPRLWPHcRkYH4\nMtwP9dy7ehTuIiL98WW4h0NphNOMzt5oqksRERmVfBnuAJmREJ09CncRkf74ONzT6NSwjIhIv3wb\n7hnhEF3quYuI9Mu34Z4ZSdOYu4hIAj4O95CGZUREEvB5uKvnLiLSHx+He5rCXUQkAf+GezhEl5Yf\nEBHpl3/DXcMyIiIJ+TbcMzTPXUQkoUHD3czuN7NaM1uf4PmLzKzJzNbGb99KfpnHyoyE6NJUSBGR\nfoWH0OYXwN3ALwdo86Jz7sqkVDREmWFNhRQRSWTQnrtz7gWgYQRqOS6aLSMikliyxtzPM7M3zeyP\nZjY3ScccUFYkRG/M0as13UVEjpGMcH8DmOqcmw/cBTyRqKGZLTWzVWa2qq6u7qTeNDMSAqBT0yFF\nRI5x0uHunGt2zrXGt5cBETMrSdD2XufcAufcgtLS0pN638yIV7qGZkREjnXS4W5mE8zM4tuL4ses\nP9njDibjUM9d4S4icoxBZ8uY2UPARUCJmVUDtwMRAOfcPcC1wE1m1gt0ANc759ywVRx3eFhGM2ZE\nRI4xaLg75z4xyPN3402VHFGZYQ3LiIgk4ttvqGalez33DoW7iMgxfBvu2eneHx3t3Qp3EZGj+Tbc\nczK8nntbV2+KKxERGX38G+7xnrvCXUTkWP4N9wwNy4iIJOLbcM+On1BtVc9dROQYvg33jHAaoTSj\nvVvhLiJyNN+Gu5mRkx6irUvDMiIiR/NtuIM37q4TqiIix/J1uGenh3RCVUSkH74O99yMMG0acxcR\nOYavwz07XcMyIiL98XW452TohKqISH98He7Z6WFNhRQR6Yevwz0nI0yreu4iIsfwd7inh9RzFxHp\nh6/DPTsjTHt3lGhs2C/8JCLiK74O9/xMb/EwrS8jInIkf4d7VgSA5o6eFFciIjK6+DvcM71wb1K4\ni4gcwdfhXqCeu4hIv3wd7vlZ3ph7c6fCXUSkL1+H+6Geu4ZlRESO5Otwf/eEqmbLiIj05etwz00P\nk2bquYuIHM3X4Z6WZuRlRjTmLiJyFF+HO3jj7potIyJyJN+He35WWMMyIiJH8X24F2RFaO7UCVUR\nkb4CEe7quYuIHMn34V6YnU5jW3eqyxARGVUGDXczu9/Mas1sfYLnzczuNLMqM1tnZmcnv8zEinPS\naWzv1rK/IiJ9DKXn/gvg8gGe/yAwM35bCvzk5MsauuKcdGIODrar9y4ickh4sAbOuRfMrHKAJkuA\nXzrnHLDSzArNbKJzriZJNQ6oKDcDgPq2borj2yK+5hz0dEDnQehug95O6O3y7ns64487IdoDLgou\nBrH4/RHb0SO3+32vIe/06hrqQRK2FQAqFsP09w/rWwwa7kNQDuzu87g6vu+YcDezpXi9eyoqKpLw\n1lCSkw5AfWs3lCXlkCLDp6sVDmyFxrehuQaa90LzHmipgfYG6GzyQj2qv0QD7YKv+iLch8w5dy9w\nL8CCBQuS8qu9+HDPvSsZhxNJntY6qH7du9W86YV60+4j20SyIb8c8ifChNMhswAyCyGr0NtOz4Vw\nJkSyIJzhbR+6D0XAQmBpkBY6ajvt2P1mCQrtZ//xtE3UPuExZCQkI9z3AFP6PJ4c3zciiuI99wbN\nmJFU62qFt1+Aqmdg+wqvdw6QFobxp0HFeVB6I5ScCsUzvFDPLFAIyrBIRrg/CXzZzB4GzgWaRmq8\nHWBcdgQzONCqcJcU6GyGzX+A9Y/Bjuch1uP1tqe9Fxb+BUxeCBPnez1vkRE0aLib2UPARUCJmVUD\ntwMRAOfcPcAy4AqgCmgHPjdcxfYnHEqjMCtCfauGZWSEOAdvPw+r/h22LvdObhZUwOIvwszLYMpi\nCKenukoZ44YyW+YTgzzvgC8lraITUJyboWEZGX6dzfDmQ/D6fd74eVYRnP0ZmPcxr4eu4RUZRUb0\nhOpwKc5J54B67jJcOpvh1Xvg5buhqwnKz4Fr7oG5H4ZIZqqrE+lXIMK9LD+TtbsPproMCZrudnj9\nZ/DSv0JHA8z6EFz4dZh8TqorExlUIMJ9YkEmy9d34pzD9KexnKzeLlj9C3jhn6GtFmZ8AN5/m9dj\nF/GJQIT7hIJMuqMxGvQtVTkZ0R5Y+yA8fwc0V8PU98B1v4Sp56W6MpHjFohwn1jgjXvWNHUq3OX4\nxaLw1qPw5+9D407v5Og1/wbT3qeTpOJbgQj3CQXeHOJ9TZ2cXl6Q4mrEN2Ix2PiEF+oHtsKEM+CT\nj3jTGRXq4nOBCPfDPffmzhRXIr7gHGxZBiu+D/vfgtLZcN2vYPaVkOb7SxyIAAEJ95LcDEJpxr6m\njlSXIqOZc7Dlj15Pfd86KDoFPnIfnP4Rbw0WkQAJRLiH0oyyvAxqmtRzl37Eol6ov/BP3gJe46bB\nNT+BeddBKBD/BESOEZif7ImFWexpVM9d+uhq9Wa/rPyJt4jXuEpY8mM44+MKdQm8wPyEVxRl8+qO\n+lSXIanmHOx+zVsmYP3j3jdKJy+CS74Fp12tUJcxIzA/6VOLs3li7R46e6JkRjR+OqZEe70107c9\nDRt/Dw07vHXST7sKFi2FyQtSXaHIiAtMuFcW5+AcVDe2M2N8XqrLkeHU3e6NnVe/5oX6zpego9Fb\nN33qBXDhN2DO1ZChnwMZuwIT7hXF2QDsPKBw9yXnoKvFW8Olo9G7tR/aPugtA3BgG9RXHXk1o3GV\ncOoH4dTLYPrF3sUvRCQ44V5ZnAPAOw3tKa5kjIv2etcAPRTKh4J6oH2HQjzRRZwBMvK9qYsVi6Hk\nM1A21/smae74kftsIj4SmHAflx0hLzPMO/VtqS4lWGIxrzfdVvfurbXPdtsBr1fdVgdt9dDdMvDx\nMvLj1wcthKxx3uXnsoogu8h7nDXOe3xoO7vIa6uLX4gcl8CEu5lRWZzDznr13IestwtaaqB575G3\nlr7b+/rvUVsIckogp9S7H1cJ2cV9grnw3YA+FOSZBZqtIjJCAvUvbcb4XFZqOuS7nPN61I07oeFt\n777vrWXvsa9Jz4X8Sd7tlIsgbwLklsWDfHw8zEu9sNZX9UVGrUCF+6llefxuzR6a2nsoyI6kupyR\nE+2Fg+9A3eb4bat3f2Ab9Bw1TJU3CYqmwfT3Q2EFFEyGvImQX+4FemZ+aj6DiCRVoMJ99gRvlszW\n2hYWVhaluJph0t0G+9Z7UwFr1nr3B7ZCtM81ZPPLoeRUOPsGKJruDZmMq/TCXJeFExkTAhXus+Lh\nvnlfQMLdOa/3vetl2LUS9q7xgtzFvOezS2DifG8KYOlsKJ3lhbp63yJjXqDCfWJBJnmZYbbsa051\nKScm2gv73oR3XoFdr3iB3n7Aey67xJv6N+caL9AnnekNp2jdcRHpR6DC3cyYVZbHln2DTMcbTRp2\nwPbnYPsKePtFby0U8IZRZl7mzeueej4Uz1CQi8iQBSrcAeZOyufR1dVEY45Q2igMw+42L8i3P+uF\neuNOb3/BFJi7xLu029TzvZObIiInKHDhflbFOB545R227m/htImjZOy5tdZbT3zLMi/Yo13elMNp\n74XFX/LGzIunq2cuIkkTuHA/c0ohAGt2HUxtuDfvhfW/hY1Peotb4aCgAhZ8DmZd4fXOQ2NouqaI\njKjAhfvU4mzGZUdYu7uRT55bMbJv3t4Am56Etx7zVirEeRddvuhWmH0FlJ2u3rmIjIjAhbuZceaU\nQtbsOjgybxiLwY7nYPUvYMtyiPV4Jz8v+iacfi2UzBiZOkRE+ghcuAMsqCxixZYt1Ld2UZybMTxv\n0rIP1vwa3ngADu7y1lVZtBTOuM6bqqgeuoikUCDD/YIZJdzx9BZe3l7PVfOTPOtk31vw8l3eeHqs\n1zsp+oFvw+wrITxMv0hERI7TkMLdzC4HfgSEgPuccz846vnPAncAe+K77nbO3ZfEOo/LvPIC8jLD\n/HfVgeSEu3PetMWX74IdKyCS4/XSF37em+UiIjLKDBruZhYC/g24FKgGXjezJ51zG49q+hvn3JeH\nocbjFkozzjulmBe3HcA5h53MEMmO5+G5f/Au6ZY7AS653ZvxkjUueQWLiCTZUHrui4Aq59wOADN7\nGFgCHB3uo8p7Zpbwp4372XGgjemlucd/gOrV8Oy34e0XvJUUP/RDOOvTGnoREV8YyoLc5UCfi1ZS\nHd93tI+a2Toze8zMpiSlupNwyWllADy9Yd/xvbBlP/zuJrjvYqjdBJf/AG5ZAwv/QsEuIr6RrKst\n/CdQ6Zw7A3gGeKC/Rma21MxWmdmqurq6JL11/8oLs5g/uYDl64cY7rEovPJjuHsBvPUoXPBVL9QX\n36RlckXEd4YS7nuAvj3xybx74hQA51y9c64r/vA+4Jz+DuScu9c5t8A5t6C0tPRE6j0ul58+kXXV\nTVQ3DnLpvQNVcP/l8PStMGUR3LwSLv0OZOQNe40iIsNhKOH+OjDTzKaZWTpwPfBk3wZmNrHPw6uB\nTckr8cRdMW8CAE+tq+m/gXOw8h645wJvnfSP/Aw+9Zi+eCQivjdouDvneoEvA0/jhfYjzrkNZvZd\nM7s63uwWM9tgZm8CtwCfHa6Cj8fU4hwWVo7jkdd345w78snOZnjkBlj+t95KjDev9L6ApC8fiUgA\nDGmeu3NuGbDsqH3f6rN9K3BrcktLjo8vrOAbj77Ja283cO4pxd7O2k3w8Ke85XYv+x6c9yWFuogE\nSuAvX/+heRPJywjzH6/t8na8/QL8/H9Adyt89ik4/8sKdhEJnMCHe1Z6iI8tmMIf1tVQv/I/4Ncf\nhfyJ8PlnvWV3RUQCKPDhDvD5C6exJO1FipbfDOUL4H8uh8KUT8UXERk2YyLcJ1Uv447wPbzq5lB7\nzYNaOkBEAi/44b59Bfz2C3RPWsQXer7B//vznsFfIyLic8EO97ot8MiNUDqLzM88ysfOm81vXt/F\n5n3Nqa5MRGRYBTfcOw7Cf3wcwunwyd9AZj63XDKDvMwIf//EemIxN/gxRER8Kpjh7hw89TXvCkkf\nfxAKvWupFman8/dXzuH1nY38+8s7U1ujiMgwCma4r/k1bHgcLv47qDj3iKc+enY5l8wezz8t38zW\n/S0pKlBEZHgFL9yba2D5rVB5obey41HMjO9/ZB75WRH+169W09TRk4IiRUSGV/DC/elbIdoNV98J\naaF+m4zPz+THnzqb3Q3tfPXhNfREYyNcpIjI8ApWuFc9Cxt+B+/9ayg6ZcCmCyuL+O6S01mxpY6/\neWydTrCKSKAMaeEwX4jF4L++DYVT4YJbhvSST55bQWN7N3c8vYWMcBrf+/A8QmlaZ0ZE/C844b7p\nSdi3Dj780+O6HN7NF02nqyfKnc9V0dDWzZ2fOIvMSP/DOSIifhGMYZlYDP78fSidDfM+dlwvNTP+\n6rJZfOfquTyzaT/X/fQVdjcMcuUmEZFRLhjhvv05qNsM7/lawpOog7nx/Ep++ulzeLuujSvveok/\nHe+FtUVERpFgDMu8+hPILYO5Hzmpw1w2dwJP3ZLHzQ++wdJfreaq+ZO4/ao5lOQOfZhH/KOpo4fN\nNc3UtXbR3hUlMz1ESW46p5bl6f+5+J7/w/1AFVT9F7z/f3tLDZykqcU5PH7z+dzz5x3cvWIbL26r\n4y8vnsmnF1eQEdZYvN9tr2vl92v3snx9DVv3tyZsN6Uoi4tOHc8V8yay+JQiTBd0EZ+xY64tOkIW\nLFjgVq1adfIHeu4f4MV/gb/aBHkTTv54fWzb38J3/nMjL1UdoLwwi698YCZLzpykkPehN3Y18uMV\n2/mvTfsxg3OnFfGeGSXMm1zIxIJMstNDdPZEqWnqZMu+FlbuqOfl7fW0d0epLM7m04un8slzK8hO\n939/SPzNzFY75xYM2s7X4e4c/Gg+FE+HG36XnML68dK2A/zj8s28taeJ8XkZ3Hh+JZ86t4LC7JP/\nS0GG196DHXxv2Sb+sK6GwuwInz2/kk8uqmB8fuagr+3sibLsrRoeem0Xr+9spDgnnc9feAqfPb+S\nrHT9gpfUGBvhvutVuP8yb/rj/OuTU1gCzjle3HaAn724gxe3HSA9lMalc8q49pzJXDizhHAoGOem\ngyIWc/z8pbf5l2e24BzcdNF0lr73lBPuea9+p4E7n63i+a11TMjP5G8/OIsl88tJ0/ciZISNjXB/\n+u/gtXvhb3ZARl5yChuCzfuaefi13fx+7R4a23soyU3n0jllXDZnAudNL9Y8+RSrbmznG4++ycod\nDXzgtPF8++q5TB6XnZRjv76zgf/z1EbWVTcxf3IB37pqDudMLUrKsUWGYmyE+90LIb8cPvNEcoo6\nTt29MZ7bXMt/vrmXP2+ppa07SnZ6iAtmlHD+9GLOn17CqWW5Ohk3Qpxz/G7NHm7//QZiznH71XP5\n2DmTk/7fPxZzPLF2D/+4fDP7m7u4ev4kvvnB2UwqzErq+4j0J/jh3vA23HkmXP4DWHxT8go7QV29\nUV7eXs8zG/fz4rY6djd0AFCck865pxQxf3IhZ0wuZN7kAnIzdFIu2Zrae7jtibf4w7oaFlaO44fX\nncmUouT01hNp7+7lnud38NPnt2MGN71vBkvfe4rG42VYBT/cX/sZLPsG/OUb3gnVUWZ3Qzuv7Khn\n5fZ6XtvZQHWjF/ZmML00lzPKC5hZlsepZbmcWpZHeWGWxm9P0GtvN/C136xlf3MnX7v0VL74vukj\nukZQdWM73//jZv6wroZJBZl884rTuOqMifqLTYZF8MP9kRthz2r42vrkFTWM6lu7WLeniXW7m3iz\n+iAb9jaxv7nr8PPZ6SFmjM9l5vg8Zk3IZWZZHrPK8phYkKmQSKA3GuPO56q4+7ltTCnK5kfXn8WZ\nUwpTVs+rO+r57lMb2bC3mQVTx3H7VXOZN7kgZfVIMAU73J2DH54GUy+Aa3+e3MJGUFN7D9tqW9i6\nv5VttS1s29/K1v0t1La8G/p5GWFmluUyr7yAMyYXMn9KAaeU5I75Xn5VbQt//dg61uw6yEfPnsx3\nlswdFcNd0Zjj0VW7uePpLTS0d/PRsydzy8UzqSge3iEiGTuCHe6N78CPzoAr/hkWfSG5hY0CjW3d\nbN3fwtbaVrbua2HzvmY27G2mvTsKQG5GmDkT85lZlsvM8V4vf+b4XErzMgLfy++Jxvjp89u589kq\nsjNCfOfquSw5szzVZR2jubOHu57dxgOvvEM05lgyfxI3v386M8aP3KwuCaZgh/u6R+DxL8AXX4IJ\n85Jb2CgVjTmqaltZV32QddVNbNjbRFVtK82dvYfb5GeGmVqcw5SiLCaPy2bKuPh9URaTCrN8/e1K\n5xzPb63j/y7bxNb9rXzojIl8+6q5lOaN7jVg9jd38rMXdvDgq7vo7I3yvlNL+fS5U3n/7PG6doCc\nkGCH+/LbYNX9cNueE14FMgicc9S1dFFV28q2Wm9oZ1dDB9WN7VQ3dtDde+TlA3MzwpTmZVCal8H4\nvAzG52UyPt/bHpeTTmFWhHHZ6YzLTicvMzwqhn6cc6x6p5E7n93Gi9sOUFmczd99aA6XzilLdWnH\npb61iwdeeYeHX9tFbUsXkwoyWXJWOVecPpHTy/MD/xeXJE9Sw93MLgd+BISA+5xzPzjq+Qzgl8A5\nQD3wcefczoGOeVLh/sDV0NUCS1ec2OvHgFjMUdfaRXVjO7sbOtjb1EFtcxd1rV3UNXdR29JJbUvX\n4aGeo6UZFGRFKMxOpzA7QmFWhPysCDkZYfIywuTEb7kZIXIzIuRkhMjNCJObGSYnPUxuRpjMSIiM\ncNoJ/ZJo7+7lTxv288ArO1mz6yDjsiN8+eKZ3LB4Kulh/34buCca49lNtTz02i5eqjpANOaYUpTF\nJbPLOG96MYunFVOQHUl1mTKKJS3czSwEbAUuBaqB14FPOOc29mlzM3CGc+6LZnY98GHn3McHOu4J\nh7tzcMd0mHUFLLn7+F8vR2jt6qWupYuGtm6aOrppbOvhYEcPTe3dNLZ72wfbu2ls76als5e2rl5a\nu3rp7Bn6RcXTw2lkRUJkRtLIjITIDPfZ7mf/gdZuVmyppb07ytTibD7/nmlce86UwM0fb2zr5pmN\n+1m2voaVO+rp7IlhBrPK8ji9vIC5k/KZMzGfaaU5lOYG/3yKDM1Qw30og7CLgCrn3I74gR8GlgAb\n+7RZAnw7vv0YcLeZmRuOMZ/WWmivh7LTk37osSg3w+tlTyvJOa7X9UZjtHVFae32Av9Q8Ld19dLS\n1Ut7Vy+dvTE6uqN09kbp6onR2ROloydKZ0+Uzvjjgx09dDZ5bQ7tz80Is+TMcq45cxILK4tGxfDQ\ncBiXk851C6dw3cIpdPVGeXN3E/9ddYA1uw/y5y21PLa6+nDbzEgaFUXZTBmXTUluBkW56RTnpFOU\n4w2jZaWHyEkPk5UeIjt+y4yECKcZoTTTL4YxaCjhXg7s7vO4Gjg3URvnXK+ZNQHFwIFkFHmE/fF5\n7WVzk35oGbpwKI2C7DQNISRJRjjEomlFLJr27jo1tc2dbKxpZldDO7vq29nV0M7uxg7e2tNEQ1s3\nvbGh950Ohfyh+0gojdCh4Idjwv/QQzPwWvTZ16f94VfZu9v6RTK46xdO4fMXnjKs7zGi0yfMbCmw\nFKCiouLEDhLJ9oZkFO4ScOPzMxMuTeyco7mjl/q2Lg529NDRHaW9O0p7d+/h7Y6eKNGYozfmiMZi\n9MYcvVEX3xfz7qMOd/iY8XscHN5+9/0OPXb9PHf410xq5mf4zkhc6Wso4b4HmNLn8eT4vv7aVJtZ\nGCjAO7F6BOfcvcC94I25n0jBTD3Pu4mMYWZGQXZEfzlJQkOZdvA6MNPMpplZOnA98ORRbZ4Eboxv\nXws8Nyzj7SIiMiSD9tzjY+hfBp7Gmwp5v3Nug5l9F1jlnHsS+DnwKzOrAhrwfgGIiEiKDGnM3Tm3\nDFh21L5v9dnuBD6W3NJERORE+ffbICIikpDCXUQkgBTuIiIBpHAXEQkghbuISAClbMlfM6sD3jmO\nl5QwHMsZjH5j8XOPxc8MY/NmmT3wAAADQElEQVRzj8XPDCf3uac650oHa5SycD9eZrZqKCuhBc1Y\n/Nxj8TPD2PzcY/Ezw8h8bg3LiIgEkMJdRCSA/BTu96a6gBQZi597LH5mGJufeyx+ZhiBz+2bMXcR\nERk6P/XcRURkiHwR7mZ2uZltMbMqM/tmqusZbmY2xcxWmNlGM9tgZl9JdU0jycxCZrbGzJ5KdS0j\nwcwKzewxM9tsZpvMbExcsMDMvhb/+V5vZg+ZWf9XJvE5M7vfzGrNbH2ffUVm9oyZbYvfj0v2+476\ncI9foPvfgA8Cc4BPmNmc1FY17HqBrzvn5gCLgS+Ngc/c11eATakuYgT9CFjunJsNzGcMfHYzKwdu\nARY4507HW048qEuF/wK4/Kh93wSedc7NBJ6NP06qUR/u9LlAt3OuGzh0ge7Acs7VOOfeiG+34P1j\nL09tVSPDzCYDHwLuS3UtI8HMCoD34l0TAedct3PuYGqrGjFhICt+9bZsYG+K6xkWzrkX8K5z0dcS\n4IH49gPANcl+Xz+Ee38X6B4TQQdgZpXAWcCrqa1kxPwr8DdALNWFjJBpQB3w7/GhqPvMLCfVRQ03\n59we4J+BXUAN0OSc+1NqqxpRZc65mvj2PqAs2W/gh3Afs8wsF/gt8FXnXHOq6xluZnYlUOucW53q\nWkZQGDgb+Ilz7iygjWH4E320iY8xL8H75TYJyDGzT6e2qtSIX5I06dMW/RDuQ7lAd+CYWQQv2B90\nzj2e6npGyAXA1Wa2E2/47WIz+3VqSxp21UC1c+7QX2aP4YV90H0AeNs5V+ec6wEeB85PcU0jab+Z\nTQSI39cm+w38EO5DuUB3oJiZ4Y3BbnLO/TDV9YwU59ytzrnJzrlKvP/PzznnAt2bc87tA3ab2az4\nrkuAjSksaaTsAhabWXb85/0SxsCJ5D6eBG6Mb98I/D7ZbzCka6imUqILdKe4rOF2AXAD8JaZrY3v\nuy1+LVsJnr8EHox3XnYAn0txPcPOOfeqmT0GvIE3O2wNAf22qpk9BFwElJhZNXA78APgETP7C7zV\nca9L+vvqG6oiIsHjh2EZERE5Tgp3EZEAUriLiASQwl1EJIAU7iIiAaRwFxEJIIW7iEgAKdxFRALo\n/wMOLzL3Fkf0MgAAAABJRU5ErkJggg==\n",
	"text/plain": [
	"<Figure size 432x288 with 1 Axes>"
	]
	},
	"metadata": {
	"tags": []
	}
	}
	]
	}
	]
	}