Skip to content

Instantly share code, notes, and snippets.

@sallos-cyber
Created June 21, 2023 08:50
Show Gist options
  • Save sallos-cyber/06a9a07e59c9262e7b24ed0148233aa6 to your computer and use it in GitHub Desktop.
Save sallos-cyber/06a9a07e59c9262e7b24ed0148233aa6 to your computer and use it in GitHub Desktop.
Create a random variable with a Beta distribution using only draws from the Uniform distribution
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEICAYAAABRSj9aAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOy9eXxV1bn//36SkIQxYZIpwRDDECZBGRxanFrRaIFaRbyOrdbrLdpWv1bpvThc2/4utr3VerFWe6toq+DQa6UKWFRoxVYDyCSTTAFOCJAEEgKZz3l+f+x9wuHkJDknOdM+We/X67zOHtZe+9lhnQ9rP2ut5xFVxWAwGAyJS1KsDTAYDAZDZDFCbzAYDAmOEXqDwWBIcIzQGwwGQ4JjhN5gMBgSHCP0BoPBkOAYoe9EiMi/isjTEah3gIhsF5G0cNdtMAQiUm05iPv+SkTuifZ9O4oR+k6CiKQC84Ff2Pv9ROQTESkXkQoR+aeIXNzK9bNF5B8iUi0iq33PqeoRYBVwdwQfwWAAmrdl+9gLIrJTRDwickeAa+4XkcMiUikiL7bUKRGRVBF5S0SKRERF5FK/Ir8A/sO2wTEYoe88zAR2qGqxvX8S+A7QH+gNPAn8RURSWrj+GPA0sKCF868C/xo+cw2GFvFvywCbgO8Bn/sXFpHpwDzgCiAHyAX+s5X61wC3AIf9T6hqCbADmNFO22OCEfp2ICLnicgGEakSkTdF5HUR+amI9BaRd0WkVESO29tZPtettsv9Q0ROishfRKSviLwqIidEZK2I5PiUVxH5nojssu/1ExE5x+59nxCRN7w9i7buDVwN/M27o6q1qrpTVT2AAG4swe8T6JlV9QNVfQM41MKf5TMgV0TObtcf1RATEqEtA6jqs6r6IVAb4DFvB36vqltV9TjwE+COQH8PVa1X1adVdQ3WbyIQq4FrWv/LxhdG6EPEboxvA4uwRHEx8E37dBLwEnA2MBSoARb6VTEHuBUYApwD/NO+pg+wHXjMr/xVwPnABcBDwAvAzUA2MBa4Kch7jwN2BniezVg/jqXA/6rq0SD+DM1Q1UZgN3Bue643RJ9Ea8utMAarx+9lEzBARPqGUIcv23FYOzdCHzoXACnAM6raoKr/BxQCqGq5qv5JVatVtQr4GXCJ3/UvqeoeVa0ElgN77N5yI/AmMNGv/JOqekJVtwJfAH9V1b0+108M8t6ZQJX/w6jqeKAX8C9Yr6wdocq+j8EZJFRbboUeQKXPvne7Zwh1+OK4dt6SP9bQMoOBYj0zGtxBABHpBjyF1XPpbZ/rKSLJqup9DTzic11NgP0efvdrq/zAIO99nBYatqrWAovtmTMbVXVToHJB0BOoaOe1huiTcG25BU5idWa8eLdD+c/CF8e1c9OjD50SYIiIiM+xbPv7/wEjgamq2guYZh/3LRsp2rr3ZmBEG3V0wRqoChl7EDePM1+RDfFNIrdlX7ZypqvlXOCIqpa30758HNbOjdCHzj+xBmnuFZEUEZkJTLHP9cTqmVSISB+a+ygjSVv3XobP66+IXCAiX7Gnk3UVkYeBAViDqojIpSKiPuWTRSQd6y0wSUTSRaSLT/1TgCJV3R+RpzNEgoRoy9A0LTId6z+DLnb79OrbK8CdIjJaRHpjTc1c5HPtahF53Gc/za4LINWuy/c/uEuwXE2OwQh9iKhqPXAdcCfW69stwLtAHdb0w65AGfApsCKKprV1778Ao0RksL2fBjwLlAPFQAFwjap6Z9VkYwmBl1uxfnzPAV+1t3/nc/5m4LfhehhD5EmgtgzwV6w2eRHWIG8N9puAqq4Afo611mO//fH9zyMb+MRnf6d9/RDgfXv7bAARGQSMBv7cwWeMKmISj3QcEfkM+K2qvhRrW1pDRO4GRqvqD4Mo+7/Am6r6fhBlz8Ka7jbR9vcbHEoituU26snCaucXBln+v7EGnX/TkftGGyP07UBELsH6X7+M0z3ZXHsxhcHgGExb7hyYWTftYyTwBtasgj3A9eaHYXAopi13AkyP3mAwGBIcMxhrMBgMCU7cuW769eunOTk5sTbDkMCsX7++TFX7R/u+pm0bIklr7TruhD4nJ4d169bF2gxDAiMiMZnrb9q2IZK01q6N68ZgMBgSHCP0hoRnxYoVjBw5kry8PBYsaB5OX0SmicjnItIoItf7HL9MRDb6fGpFZJZ9bpGI7PM5NyGKj2QwhETcuW4MhnDidruZO3cuK1euJCsri8mTJwOk+xU7gBWf/EHfg6q6CpgAYC/F3421AtPLj1T1rUjZbjCEi04r9A0NDbhcLmprzULORCU9PR2Xy0VeXh65uVastjlz5rB58+YzQsyqahGAiHhaqe56YLmqVkfK3nBg2nXik56eTlZWFl26dGm7sE2nFXqXy0XPnj3JycnhzHhFhkRAVSkvL2fz5s1kZ2c3Hc/KygJoT77POcCv/I79TEQeBT4E5qlqXXvtDRemXSc23nbtcrkYNmxY0Nd1Wh99bW0tffv2NT+GBEVE6Nu3L/X19YFOh7RK0A5kNQ4rwJWXHwOjgMlYGZUebuHau0VknYisKy0tDeW27cK068TG265DfWPrtEIPmB9DgiMiDBw4kIMHDzYdc7lcAA0hVjUbeFtVm65T1RK1qMNKezcl0IWq+oKqTlLVSf37R2fqvmnXiU17/n07tdAbEp+xY8eya9cu9u3bR319PUuWLIHQswPdhJVPtQm7l48dp3wWVmo8gyEuMULvANatW8f3v//9WJsRVyxatIh7770XgNLSUqZOncrEiRP5+OOPzyiXkpLCwoULmT59Ovn5+cyePRugVkSeEJEZACIyWURcwA3A8yKy1Xu9iORgxSv/m58Jr4rIFmAL0A/4aSSeM5Ex7bo5wbbrUOm0g7FOYtKkSUyaNCnWZsQtH374IaNGjeLll18OeL6goICCgoKm/fnz56Oqj3r3VXUtkBXoWntGzpAAxy/voNmdHtOuW6etdh0KpkcfI06dOsU111zDueeey9ixY3n99dcBWLt2LRdddBHnnnsuU6ZMoaqqitWrV3PttdcC8Pjjj3Prrbdy+eWXM3z4cH73OyvJ06233so777zTVP/NN9/M0qVLz7inx+Phe9/7HmPGjOHaa6+loKCAt96ypoE/8cQTTJ48mbFjx3L33XfjjWp66aWXcv/99zNt2jTy8/NZu3Yt1113HcOHD2f+/PkAFBUVMWrUKO666y7Gjh3LzTffzAcffMDFF1/M8OHDKSwsBKCwsJCLLrqIiRMnctFFF7Fz585mf5fVq1czbdo0vvnNbzJ69GjuuecePB5r1uNLL73EiBEjuOSSS/jkEysh0MaNG3nooYdYtmwZEyZMoKamJjz/QIZ2Ydp1nLZrVY2rz/nnn6/RYNu2bad3lj2s+mJBeD/LHm71/m+99ZbeddddTfsVFRVaV1enw4YN08LCQlVVrays1IaGBl21apVec801qqr62GOP6fjx47W6ulpLS0s1KytLi4uLdfXq1Tpz5symunJycrShoeGMe7755pt69dVXq9vt1pKSEs3MzNQ333xTVVXLy8ubyt1yyy26dOlSVVW95JJL9KGHHlJV1aeffloHDRqkhw4d0traWh0yZIiWlZXpvn37NDk5WTdv3qxut1vPO+88/fa3v60ej0f//Oc/N9nlfR5V1ZUrV+p1113X7O+yatUqTUtL0z179mhjY6N+7Wtf0zfffFMPHTqk2dnZevToUa2rq9OLLrpI586dq6qqL730UtN2q//ONsA6TdC2bdq1adeBPqZHHyPGjRvHBx98wMMPP8zHH39MRkYGO3fuZNCgQd7Vm/Tq1YuUlObetZkzZ9K1a1f69evHZZddRmFhIZdccgm7d+/m6NGjLF68mG9961vNrl2zZg033HADSUlJDBw4kMsuu6zp3KpVq5g6dSrjxo3jo48+YuvWJjc1M2bMaLJ5zJgxDBo0iLS0NHJzc5tmtAwbNoxx48aRlJTEmDFjuOKKKxARxo0bR1FREQCVlZXccMMNjB07lvvvv/+Me/gyZcoUcnNzSU5O5qabbmLNmjV89tlnXHrppfTv35/U1FRuvPHG9v/xDRHDtOv4bNfGRw9wdfP4J5FmxIgRrF+/nmXLlvHjH/+YK6+8klmzZgU1dcq/jHf/1ltv5dVXX2XJkiW8+OKLza7TFpLM1NbW8r3vfY9169aRnZ3N448/fsY83bS0NACSkpKatr37jY2NZ5TxL+db5pFHHuGyyy7j7bffpqioiEsvvTSk5zPTBkPEtGvTrm1Mjz5GHDp0iG7dunHLLbfw4IMP8vnnnzNq1CgOHTrE2rVrAaiqqmpqTL6888471NbWUl5ezurVq5t6SnfccQdPP/00AGPGjGl23Ve+8hX+9Kc/4fF4OHLkCKtXrwZoavz9+vXj5MmTTf7NcFNZWcmQIda45qJFi1osV1hYyL59+/B4PLz++ut85StfYerUqaxevZry8nIaGhp48803I2KjoWOYdr2oxXKxbNemRx8jtmzZwo9+9COSkpLo0qULzz33HKmpqbz++uvcd9991NTU0LVrVz744INm106ZMoVrrrmGAwcO8MgjjzB48GAABgwYQH5+PrNmzQp4z29961t8+OGHjB07lhEjRjB16lQyMjLIzMzku9/9LuPGjSMnJ6fpBxZuHnroIW6//XZ+9atfcfnlLU9aufDCC5k3bx5btmxpGsBKSkri8ccf58ILL2TQoEGcd955uN3uiNhpaD+mXcdpu27JeR+rT0wGYx3EY489pr/4xS8Cnjt16pTm5uZqRUVFi9dXVVWpqmpZWZnm5uZqSUlJi2Vjge8AXTjo1IOxDsK069AItV2bHn2C8MEHH/Cd73yHBx54gIyMjBbLXXvttVRUVFBfX88jjzzCwIEDo2ilwRAapl2HB9EWBjLOKCRyFfBrIBn4X1Vd4Hf+HmAu4AZOAner6jb73I+BO+1z31dV38BQzZg0aZJGI93a9u3byc/Pj/h9DLEl0L+ziKxX1aiv1IlG2zbtunMQartuczBWRJKBZ4GrgdHATSIy2q/Ya6o6TlUnAD/HDudql5sDjAGuAn5j12cwGAyGKBHMrJspwG5V3auq9cASYKZvAVU94bPbndNhYGcCS1S1TlX3YWXoCRjlz2AwGAyRIRgf/RDgoM++C5jqX0hE5gIPYCV18A49DwE+9bu2WdwQEbkbuBtg6NChwdhtMBgMhiAJpkcfaDZ/M8e+qj6rqudgJWCYH+K1UY/ZbTAYDJ2FYHr0LqwwrV6ygEOtlF8CPNfOa2NGzrz3wlpf0YJrwlrf0qVL2bZtG/PmzQtrvcFSUFDAa6+9RmZmJs888wzPPfcc5513Hq+++mpM7DEEh2nXrdNZ2nUwQr8WGC4iw4BirMHVf/EtICLDVXWXvXsN4N1eCrwmIr8CBgPDgcJwGN7ZmDFjRlNsjliwbNmypu3f/OY3LF++POiclY2NjQFjmxgMpl1HhzZdN6raCNyLlS9zO/CGqm71TdwA3CsiW0VkI5af/nb72q3AG8A2YAUwV1XNckasEKhjx45t2v/lL3/J448/DsAzzzzD6NGjGT9+PHPmzAHOTEhwxx138P3vf5+LLrqI3NzcpqXdrYVr9eXSSy/FO82vrKyMnJycpntcd911XHXVVQwfPpyHHnqo6ZqcnBzKysq455572Lt3LzNmzOCpp57i2LFjzJo1i/Hjx3PBBRewefNmwAo7e/fdd3PllVdy2223sWjRImbNmsU3vvENhg0bxsKFC/nVr37FxIkTueCCCzh27FiLz25wDqZdx2e7Duq/I1VdBizzO+abuOEHrVz7M+Bn7TWwM7JgwQL27dtHWloaFRWBs96VlJSwZs0aduzYwYwZM7j++uv5v//7P4qKitiyZQtHjx4lPz+f73znOyHde+PGjWzYsIG0tDRGjhzJfffdR3b2ae/bb3/7W1asWMGqVavo168f9913HxMnTuTPf/4zH330EbfddhsbN24EYP369axZs4auXbuyaNEivvjiCzZs2EBtbS15eXk8+eSTbNiwgfvvv59XXnmFH/7wh0E9u8GZmHYdu3ZtgprFIePHj+fmm2/mj3/8Y4uvhrNmzSIpKYnRo0dz5MgRoPVwrcFyxRVXkJGRQXp6OqNHj2b//v2tll+zZg233norAJdffjnl5eVUVlYC1mt5165dm8pedtll9OzZk/79+5ORkcE3vvENgDNCvgbz7AZnYtp17Np1Qgt9uAeiwklKSkpThhngjPCp7733HnPnzmX9+vWcf/75ASP9+YZP9a5uDmaVs/+9fe/rX29ycnLAe/sS6J7esKvdu3dvse6WQr4G8+yG+MW06/hs1wkt9PHMgAEDOHr0KOXl5dTV1fHuu+8Clj/y4MGDXHbZZfz85z+noqKCkydPBlVnS+Fa/cnJyWH9+vUAHQ7dOm3atKYZCqtXr6Zfv3706tWrXXV15NkN8YFp182Jh3adsO/Gofbmwz1trC26dOnCo48+ytSpUxk2bBijRo0CwO12c8stt1BZWYmqcv/995OZmRlUnS2Fa/XnwQcfZPbs2fzhD39oNaxqMDz++ON8+9vfZvz48XTr1q1DiYw78uyGwJh23T4Srl23FNYyVp9whXI9++F39eyH323xvFPDubZFvIdrjTYmTHFiYNr1mZgwxZ0cE661OStWrOAHP/gBbrebu+66q9l5EZkGPA2MB+ao6ls+59zAFnv3gKrOsI8Pw1oc2Af4HLhVrVhQhghg2nXHSEihj+dB2EjTkv+ys+J2u5k7dy4rV64kKyvLm2Uo3a/YAeAO4MEAVdSoFZXVnyeBp1R1iYj8FisU93MByhnCgGnXHaNTD8ZqkKP5BmeiqmzZsoW8vDxyc3NJTU31LlbJ9CtXpKqbAU/AivwQa/rF5YC35/8yEDjPXQww7Tqxac+/b6cV+vT0dMrLy82PIkFRVcrLyzl27NgZC2OysrLAirAaLOkisk5EPhURr5j3BSrUWjUOLURlBSsyq339utLS0tAfJERMu05svO06Pd3/pbR1EtJ1EwxZWVm4XC6i8eMzxIb09HT69OkT6FQoKjhUVQ+JSC7wkYhsAU4EKBewTlV9AXgBrAxTIdy3XZh2nfikp6d7OyxB02mFvkuXLkEHLzI4l8OHD3Pw4Ol0Ci6XC6Ah2OtV9ZD9vVdEVgMTgT8BmSKSYvfq4yYqq2nXhkAkhOumMw++Glpn8uTJ7Nq1i3379lFfX8+SJUsAggo2IiK9RSTN3u4HXAxss6eyrQKut4veDrwTfusNhvCQEELfRNkuVj5yOZckbYq1JYY4ISUlhYULFzJ9+nTy8/OZPXs2QK1v9FURmSwiLuAG4HkR2Wpfng+sE5FNWMK+QO2k91gJdh4Qkd1YPvvfR/O5DIZQSCzXzcrH+HryesYn7eGiuv/BjclDbrCSSxQUFDTtz58/3z/66los98sZqOo/gHGB6lTVvZj8xwaHkDg9+vpq2PMhJdqHAVLBlKQdsbbIYDAY4oLEEfojW6GxlgUNc3CrMDVpe6wtMhgMhrggcYT+qOU6/VyH86VmMUH2xNgggyFxMBMenE0CCf126NINl/Znk+ccxiftIbTp0gaDwZCYJI7Ql++GvnkoSezUbPrISfoGXNdiMBgMnYvEEfoThyDDWuq+RwcDcI7ExRoWg8FgiCkJJPTF0MsS+D0eW+iTjNAbDKEQyBffln/e+O/jn4QQ+q7UQm1Fk9Afoi81mmp69AZDBzEinhgkhNAPlOPWRi8rgKCSRJEOJEcOx9Aqg8FgiA+CEnoRuUpEdorIbhGZF+D8AyKyTUQ2i8iHInK2zzm3iGy0P0vDabyXgXLM2rB79AAu7ccQKYvE7QwGg8FRtCn0IpIMPAtcDYwGbhKR0X7FNgCTVHU8VjKGn/ucq1HVCfZnRpjsPoOzsHv0PU+nF3Npf7KkDExcboMhahhXT3wSTI9+CrBbVffaOTGXADN9C6jqKlWttnc/JUDckEjSW05aG936Nh1zaT96So3luzcYDIZOTDBCPwQ46LPfYjYdmzuB5T77gTL0nEFHs/D0lipAID2j6Vix9rc2Kg4GvshgMISM6bE7k2CEXgIcC+gPEZFbgEnAL3wOD1XVScC/AE+LyDnNKlN9QVUnqeqk/v37B2HSmfTmJHTtDUmno1W6tJ+1UXEg5PoMBoMhkQhG6F1Ats9+wGw6IvI14D+AGapa5z3um6EHWI2VoSes9JEq6HZmyrhir9BXmh69wRAOfHvzpmfvLIIR+rXAcBEZJiKpwBzgjNkzIjIReB5L5I/6HA+YoSdcxnvJpAq6nin0x+nJKU0zrhuDoYNkUsXvuvySnWm381KXJ+Hk0bYvMsQVbQq9nRPzXuB9YDvwhqpu9c3Qg+Wq6QG86TeNsrUMPWGjj5w8YyDWQjisfaDKLJoyGNqL4OG3qU8zLWkLf3ZfzAVJ2+G1G0mhMdamGUIgqAxTqroMWOZ3zDdDz9dauK7FDD3hJDOA6wbgqPamfMs2ptwQaQsMhsTkG0n/5IKk7Tzc8F1ed1/Gas+5PHfo19yR/D5+k+8McUxCrIztQ2ChP0xvBnIsBhYZDAmAKvel/JltnrN5w30JAMs9U1ntPpfvpbwDDTUxNtAQLM4X+vpq0qUBuvZpNkB0RHszQCrMoimDoT0cLGR4UjGL3FeiPlLxm8YZlrt00+IYGmcIBecnB6+ttL595tB7Oaq9SZMGqDkesMdvMBhaYeOrnNI03nNfcMbhQh3FVs/ZjPn8FXLeGhAj4wyh4PwefV2V9R1A6I9ob2ujqiSKBhnijRUrVjBy5Ejy8vJYsGBBs/MiMk1EPheRRhG53uf4BBH5p4hsteM43ehzbpGI7POJ4zQhSo8THTwe+PJ9VnkmcIqufieFt91fgUMbGCbmt+UEEkbov714R7NTh43Qd3rcbjdz585l+fLlbNu2jcWLFwOk+xU7ANwBvOZ3vBq4TVXHAFdhLfjL9Dn/I584Thsj9Aix4fAmOHmYVe7Ay17+4r4QEGYmfxJduwztIgGE3koXeFL9ex1wBK/Qm3DFnZXCwkLy8vLIzc0lNTWVOXPmAPiKNapapKqbAY/f8S9VdZe9fQg4CoS+dNuJfPlXQFjtOTfg6SP0geypXJH0eXTtMrSLBBB6q0d/stnrJZSq/Xs+YXr0nZXi4mKys08v7M7KygJIDbUeEZliX7fH5/DPbJfOU96FgQGu61Acp5hR9DEMHEc5zV2iTYy4knFJRfT3Ro81xC0JI/RVAYS+jlSOaw/juunEaOAZVyFNwxKRQcAfgG+rqrfX/2NgFDAZ6AM83ML9OxTHKSa4G6B4PQy9sPVyw68E4NLkTVEwytAREkboA7luwBqQ/etnieU+NQRPVlYWBw+eDoPhcrkAGoK9XkR6Ae8B81X1U+9xVS1RizrgJaxw3onB4c3QUM3cNQFfUk4zYCyHtA+XJZnfV7yTMELffGaAhTWX3rxadlYmT57Mrl272LdvH/X19SxZsgQgqCQFdmynt4FXVPVNv3OD7G8BZgFfhNfyGHLgMwDWeUa0Xk6ET9xjuTBpmzVLxxC3JIDQn6BWu9DQwpKAUjLpLyb5SGclJSWFhQsXMn36dPLz85k9ezZArW+sJhGZLCIu4AbgeRHZal8+G5gG3BFgGuWrIrIF2AL0A34a1QeLJAc/hcyh1oBrG3ym+Vbin9LtUTDM0F6cv2CqrooqurV4ukwz6EeltTpWAoXWNyQ6BQUFFBQUNO3Pnz/fP1bTWgJkRVPVPwJ/DFSnql4eAVPjg5JNMHgiBDFZ7VNPvrWx/x8wYExk7TK0mwTo0VdR1YJ/HqBUM0gVt7U61mAwtEpPquF4EQwcH1R5l/anWPtC0ZrIGmboEAkh9IGmVnopVXt62CkHTW0zGGLEKLEzsgUp9CB85smndOtHJqZUHJMYQt9Kj77MOw/YJEswGNpkdNJ+a2Ng8NHFCz2j6C8n4NjeCFll6CiJIfSt9ujtRVMnj0TJIIPBuYyRIujWD3oODPqaDZ48a8O1LjJGGTpMAgj9iYCLpbyUaS9rw7huDIY2GZ20HwaND2niwi7NstJ2Fhuhj1cSQuhPteK6qaAHjZpkXDcGQ1t43AwXV8izZzwksUVzrdW0hrjE+UJfX82pZsEIT6MkWX56I/QGQ+tU7CdNGqHfyKCK+yb62ejJg8NbSA1+0bEhijhb6D1ucNdRra0v1S7TDDhlhN5gaJXSL63vfm2siA3ABs854K5ntOwPs1GGcOBsoa8/BUANrQt9qZoevcHQJmVeoR8e8qUb7QHZCUm7w2mRIUw4W+jt5MRtCX2ZEXqDoW3KvqRUe7Ur7eYR+kDPQUbo4xSHC73Vo2/TdUOGNevGBF4yGFqmbBd7dXD7rx88kbFSFDZzDOHD2UJfXw0E47rJBE8D1JrgZgZDi5R9yR5PB4R+4HhypaTJpWqIH4ISehG5SkR2ishuEZkX4PwDIrLNzrbzoYic7XPudhHZZX9uD6fxIbluwLhvDIaWOFUONcfYY/fofWfUBM2g8SSJwpGtbZc1RJU2hV5EkoFngauB0cBNIjLar9gGYJKqjgfeAn5uX9sHeAyYipWY4TER6R02623XTY22nhmu1BsGwcy8MRgCYw/E7umI68YbH6fEZJyKN4Lp0U8BdqvqXlWtB5YAM30LqOoqVa22dz/ldMjX6cBKVT2mqseBlcBV4TGdJtdNdTCzbsD06A2GlrCFfncHhD5nwSaOaQ8rQ5UhrghG6IcAB332XfaxlrgTWB7Kte1OoNwQnI/euG4MhjY4vo96TeaQ9utAJcJWTw6b130MtNP9Y4gIwQh9oKAXAeORisgtwCTgF6Fc2+4Eyl6hb2PWTSXdIamLcd0YDC1xbB8u7Y+ng/MztmoOI+WglWDcEDcE86/qArJ99rOAQ/6FRORrwH8AM+yEyUFf226CdN0oSdC9v+nRGwwtcbyIg3pWh6vZ5smxwiiU7gyDUYZwEYzQrwWGi8gwO1nyHGCpbwERmQg8jyXyvmr6PnCliPS2B2GvtI+Fh4bgVsYC0MMIvcHQIsf3sV8HdLiarWpPuDN++riiTaFX1UbgXiyB3g68oapbfZMrY7lqegBv2pai63EAACAASURBVAmUl9rXHgN+gvWfxVrgCftYeGioAYQ6urRdtscAE5PeYAhEzXGoreRAGHr0+3SQtYCxxAh9PBFUcnBVXQYs8zvmm1z5a61c+yLwYnsNbJX6akjtDrVBxM7ufhYc3hIRMwwGJ3PtE3/k3TTC4rrxkMR2Hcr5JZuAr3bcOENYcPbK2IZT0KVbcGV7nGXCIHRSVqxYwciRI8nLy2PBggXNzovINBH5XEQaReR6v3MBF/yJyPkissVeRPiMSAiZOuKMs8VyaYbDdQOw3TMUjm6lhTkbhhjgcKGv4UBVkI2px1ngabReUw2dBrfbzdy5c1m+fDnbtm1j8eLFQLMEBgeAO4DXfA+2seDvOeBuYLj9Cd/6kCgz1Bb6gxrCjLdW2KFDobaSgYTPS2voGM4W+vpTVLeSdOQMetivpcZP36koLCwkLy+P3NxcUlNTmTNnDkCmbxlVLVLVzYD/617ABX8iMgjopar/VFUFXgFmRf5pIkO2HIHu/TnVSkrOUNjhsSbajUo6EJb6DB3H2ULfUE0trYc/aKKH/Vpq5tJ3KoqLi8nOPj3DNysrCwi20bS44G+Ive1/vBntXgwYRc6Wo9A7J2z1fam20MvBNkoaooXDhb6mzRDFTXiF3kyx7FRYHe7mh4O8vKUFf0EvImz3YsAoMtRP6Du6ovUE3SEj2/To4whnC339qTYXSzXR3f6RGddNpyIrK4uDB0/3LF0uFxB0YtOWFvy5OB3Pyfe482isZ7CUQe9h4a33rNGmRx9HOFvoG6qpDVbo0zMgOc0IfSdj8uTJ7Nq1i3379lFfX8+SJUsAgk1MEHDBn6qWAFUicoE92+Y24J2IPECkqTxIsij0PrvtsqEwYAznyCForA9vvYZ24XChD8F1I2IvmopPP6khMqSkpLBw4UKmT59Ofn4+s2fPBqj1XfAnIpNFxAXcADwvIluhzQV//wb8L7Ab2MPpQH7OotIeasgcGt56B4yhi7hP56E1xJSgFkzFLaG4boANx1OZ2Nf06DsbBQUFFBQUNO3Pnz/ff8HfWs50xeBzLuCCP1VdB4wNv7VRxiv0vYYAO8JX74Ax1vfRbTDQ+X8mp+P4Hn3Qs26w49KbwViD4TQniq3vXq1FHm8HffOo12Q48kV46zW0C+cKvSq466gLQejLNMP46A0GXypdlGov6BLkepRgSe7Cbs2CI9vCW6+hXThX6BtrAahtI42gL6VkQnU5uBsjZZXB4CwqXZRo34hUvV2zTf7YOMHxQh9U5EobK6WgQnVZhIwyGBzGieIOZpVqmZ2ebKg6BNUmFEKscbDQW7lNQvPR2yvfjfvGYLDcn5UuSrRPRKrfofZMnqPGfRNrnCv0DTUA1GmoPXrMgKzBAFBbCfUnORQh180Ojy30xk8fc5wr9HaPPljXTc689yjFCL3B0IQ94yZSrpujZELXPmbmTRzgYKG3evShuG7K7B79k2/9PSImGQyOotIS+ki5bkCs+fTGdRNzHCz0ofXoAWpIh9Se9JdgV8AbDAlMpRWLJlKuG8AS+iPbTMKfGONcobd99KFMrwSgx1lG6A0GsFw3ksxRerddtr2cNdrKBFdRFLl7GNrEuULfjh49AD3Ooh8nImCQweAwKouh12A8EZCBplDHA+zwB2ZANqY4WOjtBVMh+OgB06M3GLycKA5/6AN/zhoFiFk4FWMcL/Sh9+gHGKE3GMDy0WdYQt/RZCMtktod+uSamTcxxvFCH6qP/hefVJAh1dBQGwmrDAZn4PHAiUOQETBoZ3gZMMb06GNMUEIvIleJyE4R2S0i8wKcnyYin4tIo4hc73fOLSIb7c/ScBnuFepQe/Rl3rn0p0xcekMnproM3PXQKxpCPxaO7YX6U5G/lyEgbQq9iCQDzwJXA6OBm0RktF+xA8AdwGsBqqhR1Qn2Z0YH7T1NO330ZnWswUDT1Eqv6yaiDBgDKBwNY7x7Q0gE06OfAuxW1b2qWg8sAWb6FlDVIlXdDERvsmw7ffQm3o3BQNNiqei4bux+ofHTx4xghH4I4Jvl12UfC5Z0EVknIp+KyKxABUTkbrvMutLSIF0qjbXUaQoa4jDD6R69EXpDJ6Yp4UgUhD4zB7p0NytkY0gwKikBjmkI9xiqqpOAfwGeFpFzmlWm+oKqTlLVSf379w+u1oba0GfcAOXGR28wQKWLWu0C3SIV/uA0Of++3OrVmwHZmBGM0LuAbJ/9LOBQsDdQ1UP2915gNTAxBPtaprF9Qt9ACse1h+nRGzo3lS6KtR9IoH5cBBgwxnLdaCh9REO4CEbo1wLDRWSYiKQCc4CgZs+ISG8RSbO3+wEXA+F5f2sMLY2gL6UmpaChs3OimBLtE7n58/4MGAs1x6GqJDr3M5xBm0Kvqo3AvcD7wHbgDVXdKiJPiMgMABGZLCIu4AbgeRHxvqPlA+tEZBOwCligqmES+pqQYtH7UqqZZtZNJ2LFihWMHDmSvLw8FixY0Oy8iKSJyOv29OHPRCTHPn6zz9TgjSLiEZEJ9rnV9pRj77mzovpQHSWCKQQDMmCM9W3cNzEhJZhCqroMWOZ37FGf7bVYLh3/6/4BjOugjYFprAs9/IFNKRlw0vQsOgNut5u5c+eycuVKsrKymDx5MoB/Juw7geOqmicic4AngRtV9VXgVQARGQe8o6obfa67WVXXReExwou7AaoOc4gLo3fPs3xm3gz/evTuawCcvDK2oaZdPnqw49KfPGL8hZ2AwsJC8vLyyM3NJTU1lTlz5gBk+hWbCbxsb78FXCHSzHl9E7A4stZGiaoSQCMbntimyTXUNROX9jM9+hjhXKFvrAs9RLHNEe0NDdVQZ6JYJjrFxcVkZ5+eS5CVlQU0exVsmkJsuyorAX8VvJHmQv+S7bZ5JMB/DEA7pw5HmkoXEMmEI4HZ4ck2Qh8jHCz07e/RH/E28BPGfZPoaOC3Nv+DrU4hFpGpQLWq+q74uVlVxwFftT+3tnD/0KcORxp7sVRxhFIItsQOHQplX0JjPRDBQGqGZjhY6Nvvoz+sdqKFqqBniRocSlZWFgcPnl7v53K5ABr8ijVNIRaRFCADOOZzfg5+vXlVLba/q7BCf0wJs+mR44S3Rx/FwVhgpycbPI2W2BuiioOFvn3z6AEOY3r0nYXJkyeza9cu9u3bR319PUuWLAHwj1O9FLjd3r4e+EjtVwERScKaTbbEW1hEUuzpwohIF+BawDnr+ytdkJ5BdbMx6ciyXYdaG8Z9E3WcK/QNtR3z0YPp0XcCUlJSWLhwIdOnTyc/P5/Zs2cD1PpODwZ+D/QVkd3AA4BvhNZpgMte8OclDXhfRDYDG4Fi4HcRf5hwUVkMGdltlwsz+3QQJKeamDcxIKjplXFJB3r0daRyXHvQ2/ToOwUFBQUUFBQ07c+fP99/enAtVq+9Gaq6GrjA79gp4PyIGBsNTrgin1kqAG6S+aJhMGNNjz7qOLdH31jbbh892H56s0rP0BmpdEUnPHEAduhQ47qJAc4UetUO9ejBnnlzwrhuDJ2M+morFEEMevQA2z3ZcPIwnCqLyf07K84UencDqIe6dvrowfToDZ0Ub3jiGPjowe7Rg+nVRxlnCn1TdqkO9OjpY8W7cfvPtDMYEphoZpYKwA6PJfRP/P7NmNy/s+JooW9v9ErwzqVXE8XS0LnwZpaKsuvGuziqnAxKNYNRciCq9+/sOFzo29+jP2xWxxo6IyeKAYFeg2NmwnbPUEYn7Y/Z/TsjDhX6OoAO+ejNXHpDp6TyIPQ4C1LSYmbCFzqMEXKQ1GYLlA2RwplC31ADdMxHb3r0hk5JZXHMZtx4+cKTQ6q4GSEH2y5sCAvOFHpvj74DPvpj9LRW6VUdMsGVDJ2HE8WQEYWE4K2wRYcBMDapKKZ2dCYcKvR2j74DrhsQ6DnQ9OgNnQdVe7FUbIX+oJ7FCe3GWNkXUzs6Ew4Vem+Pvv2uG8B6hTVz6Q2dhZrjVh6GGLtuQPjCk8PYJCP00cKZQt/ko+9Ijx5r5oF3AYnBkOg0LZaKbY8erAHZfDlo1rFECWcKfbh69BlZUOlC8ITBKIMhzrEzS8WF0HuGkSYNULoz1qZ0Chwq9PY8eu2o0GeDu55+VIbBKIMhzoknodcca6NkU0zt6Cw4Wug77LrJtJZjZ4kJsGToBFS6IKkLdD8r1pawTwdyUtOhZGOsTekUOFroO+66sQI7ZUmcJG02GCJJpQt6DSbn35fH2hKUJLZqjunRRwlnCn1DuHr0ltAPMT16Q2fgRGwyS7XEVk8OHN4CHnesTUl4ghJ6EblKRHaKyG4RmRfg/DQR+VxEGkXker9zt4vILvtzu/+17aKxFiQZN8kdqyetJ6RnGqE3dA7iYA69L194cqzpnmW7Ym1KwtOm0ItIMvAscDUwGrhJREb7FTsA3AG85ndtH+AxYCowBXhMRHp32OrGWujStcPVAJCZbYTekPh43HDiEP/zufU2HA+rwbdorrVxaENsDekEBNOjnwLsVtW9qloPLAFm+hZQ1SJV3QzN5ilOB1aq6jFVPQ6sBK7qsNWNtVb4gnCQebYRekPiU3UY1E2J9o21JU3s0cGQ2gOK18XalIQnGKEfAvhGH3LZx4IhqGtF5G4RWSci60pLgxgYbazlUHWQFrRFht2jVw1ThQZDHGJPrTwUR0LvIQmGnAcuI/SRJhihlwDHglXFoK5V1RdUdZKqTurfv3/btTbWdTDOjQ+Z2fSQWmt5uCEhWbFiBSNHjiQvL48FCxY0Oy8iaSLyuj0G9ZmI5NjHc0SkRkQ22p/f+lxzvohssa95RkQCtfX4wc4sVaz9YmyIH1mT4cgXTavdDZEhGKF3Ab5D9VlAsEHcO3Jty3QwMbgv//qXo9ZGhcl4k4i43W7mzp3L8uXL2bZtG4sXLwZI9yt2J3BcVfOAp4Anfc7tUdUJ9ucen+PPAXcDw+1Px12SkcQOf1DiDc8dLwyZBJ5GM80ywgQj9GuB4SIyTERSgTnA0iDrfx+4UkR624OwV9rHOkZDeIQ+Z957p3s4lSY2diJSWFhIXl4eubm5pKamMmfOHIBMv2IzgZft7beAK1rroYvIIKCXqv5TVRV4BZgVAfPDR6UL0jI4SbdYW3ImWZOsb9fa2NqR4LQp9KraCNyLJdDbgTdUdauIPCEiMwBEZLKIuIAbgOdFZKt97THgJ1j/WawFnrCPdYzG2g7FovelSegrjNAnIsXFxWRnn36pzMrKApo1nqaxJLu9VwJeZ/YwEdkgIn8Tka/6lHf5XN/iuFXI40+RIs6mVjbR4yxrhboR+oiSEkwhVV0GLPM79qjP9lost0yga18EXuyAjc1prOt4nBub4/SkSrvS83hRWOozxBcaeJDd/2BLY0klwFBVLReR84E/i8iYVsoHuv8LwAsAkyZNit2If+VBPiwJ07hWuMmaDAc+i7UVCY0zV8aG0UcPwn4dAMf2hqk+QzyRlZXFwYOn39ZcLhfQLFlp01iSiKQAGcAxVa1T1XIAVV0P7AFG2OV9OzbhGXuKJJXF8eef95I1GU644ERJXMzvT0QcKvR1YXPdABTpAPZ+uSVs9Rnih8mTJ7Nr1y727dtHfX09S5YsAajwK7YU8K7avh74SFVVRPrbCwYRkVysQde9qloCVInIBbYv/zbgnag8UHuoPwU1xzgUbzNuvAyx/fRmPn3EcKjQ14SxRw/7dQDZUgruxrDVaYgPUlJSWLhwIdOnTyc/P5/Zs2cD1PqOMQG/B/qKyG7gAcAb5mMasFlENmEN0t7jM8b0b8D/AruxevqxjxTWEpXWjJviOJpD7yVn3nswaDx1mmL89BEkKB993BFGHz1AkQ6ki7ih8gD0yQ1bvYb4oKCggIKCgqb9+fPn+48x1WJNJDgDVf0T8KdAdarqOmBs+K2NACesceN4WhV7BilpbNUczjvwGVa0FEO4cWiPPpw+eijyDLQ2jJ/ekIh4V8USp64boNCTD8XrSacu1qYkJA4V+vD76AE4ZpIVGxKQShcgHNaOxxOMFJ96RoGngYlJu2NtSkLiPKFXDXuPvpRMTmma6dEbEpPj+6HXEBrj2FO73jMSEKbIjlibkpA4T+jd9UAY8sWegbBfBxqhNyQmFfuh99mxtqJFcua9RxXdYOA4piZtj7U5CYnzhN4OftTh7FJ+FOkAKN8T1joNhrjg+H7IjF+hb+LsizkvaRc01pv59GHGeULfaA3WhNN1A9YUS44XmSmWhsSisQ6qSuK6R99EzsWkS4NJRBIBHCj0YUoM7sceHQyeBkvsDYZEoeIgoM7o0Q+9EICfP//7GBuSeDhQ6O0efVh99PClx17RXrrDvDYaEoeKIgBueD2+IzQA5PzkM770DDF++gjgQKGPTI9+t9rBB0vNqL8hgTi+H4CDGkRCnzjgE89YpiTtILVZOCJDR3Cw0Id3MLaadMgYaoTekFhU7IfkVI4Qv3PoffnYM46uUs/5SV/G2pSEwrFCH7ZUgr70H2mE3pBYHN8PGdmoQ37qn3pGU6/JTEvaHGtTEgpn/Ov7EqFZNwD0H0ltyQ6S8IS/boMhFsT5HHp/qknncx3BV5O2mLGyMOJAoY+Mjx6As/JJlway5Wj46zYYYoFT5tD78LF7HGOTiuhLZaxNSRgcKPSR7NGPAmCEuNooaDA4gLoqqDnmqB49WH56gIuTvoixJYmDA4Xe7tGHeXolYPnoMUJvSBDsGTdO69F/ocM4rj24JNn46cOF84TeDoEQ7lk3AKT1pMgzgNFJReGv22CINt7YTQ7LseAhib97xnNJ0iZy5/0l1uYkBM4Tett1UxsJ1w3wheYwTky4YkMCUG6H/O17TmztaAcr3efTT05wnphpluHAgUIfmXn0YEXR2+oZxtCkUqg5Hvb6DYaocmwP9BgAaT1jbUnIrPacS70m8/Xk9bE2JSFwoNDX4VGhgeSIVL9Fh1kbJcY/aHA45Xugj/N68wAn6canntFcmbTOykFh6BBBCb2IXCUiO0Vkt4jMC3A+TURet89/JiI59vEcEakRkY3257cdtrgp6Yh0uKpAbPXYA1clmyJSv8EQNcr3ONJt4+WvnkkMSzoCpTtjbYrjaVPoRSQZeBa4GhgN3CQio/2K3QkcV9U84CngSZ9ze1R1gv25p8MWhzm7lD/H6YVL+xmhNzib2ko4ddTRQr/Sfb61sdMsnOoowfTopwC7VXWvqtYDS4CZfmVmAi/b228BV4hIZLrcERZ6gC88w6BkY0TvYYgeK1asYOTIkeTl5bFgwYJm51t5I/26iKwXkS329+U+16y233K9b6tnRe2BgsGbRKdvXmzt6ABH6MNGzzmw9e1Ym+J4ghH6IcBBn32XfSxgGVVtBCqBvva5YSKyQUT+JiJfDXQDEblbRNaJyLrS0tLWrWmsi8wceh82ec6xZiycKo/ofQyRx+12M3fuXJYvX862bdtYvHgxQLpfsZbeSMuAb6jqOOB24A9+193s87YaX8upm6ZWOrdHD/CO+yI4vAWOmhhUHSEYoQ/UM/cfHWmpTAkwVFUnAg8Ar4lIr2YFVV9Q1UmqOql//zbCqTbWhj2NoD/rPCOsjYOfRfQ+hshTWFhIXl4eubm5pKamMmfOHIBMv2IB30hVdYOqegO5bwXSRSQtOpa3n5x570H5bjwq0GdYrM3pEO+6LwRJgi1vxtoURxOM0LuAbJ/9LMA/i0FTGRFJATKAY6pap6rlAKq6HtgDjOiQxY11EXfdbNZcSE6Fg59G9D6GyFNcXEx29unmm5WVBTTrKbT2RurlW8AGVa3zOfaS7bZ5JGKuyvZSvptD9IUuXWNtSYcoJROGXWIJvZl9026CEfq1wHARGSYiqcAcYKlfmaVYr7YA1wMfqaqKSH97MBcRyQWGA3s7ZHFjbWRWxfpQRyoMnggHjNA7HQ0sDsG+kVonRcZguXP+1ef8zbZL56v259ZANwrJLRlGtm0qZJfH38PqTB7cOcKKwulaG2tTHEubQm/3cO4F3ge2A2+o6lYReUJEZtjFfg/0FZHdWC4a7xTMacBmEdmE9Up8j6oe65DFDbXUa0qHqgiK7KlWkuKG2sjfyxAxsrKyOHjw9BCTy+UCmqUvCvhGau9nAW8Dt6nqHu8Fqlpsf1cBr2FNWmhGSG7JMJGMm3OkmJ2a3XZhB7DCPRlS0mHjq7E2xbEENY9eVZep6ghVPUdVf2Yfe1RVl9rbtap6g6rmqeoUVd1rH/+Tqo5R1XNV9TxV7XjgisYaqomCm3ToheCuh0OfR/5ehogxefJkdu3axb59+6ivr2fJkiUAFX7FWnojzQTeA36sqp94C4tIioj0s7e7ANcCcRNqMUcOkyaN7PLmQXY4J+kGY66DLW9ZETkNIeO8lbENNdRGQ+jPtgeB9qyK/L0MESMlJYWFCxcyffp08vPzmT17NkBtkG+k9wJ5wCN+0yjTgPdFZDOwESgGfhfN52qN4VIMwE7NSpzkHZO+A/UnzaBsO4mCDyTMNNRQE4k0gv507Q1DJsGeD+Hy/4j8/QwRo6CggIKCgqb9+fPno6qPevdVtRa4wf86Vf0p8NMWqj0/3HaGi5FyEI/K6YT3iUDWJLZ7hpK/7kU4/9sQZ2Pf8Y4De/TV1EShR58z7z2e2pcNxZ9DdceGFQyGaDI8ycV+PavpzTchevUivOq+wppT71oXa2schwOFvoaaCM+68fJ3z3hAYc9HUbmfwRAORoqLLxNkINaXt91fgfQM+MevzzieEP+RRRhnCb0qNFRHx0cPbNJzqNDusPvDpmOmURnimvpqcuUQOxJQ6E/RFSZ/F7a/C2W7Ym2Oo3CW0NvZpaLio8fKdPORZyLsXAZu/xl5BkMccngLyaJs8Tgrq1RbNHWwpt4DKWksefrB2BrkMJwp9FHq0QMsc0+F2grY+7eo3dNgaDd2ML7NCSb0XnJ+Wsgfar/CdckfQ8WBWJvjGBwm9NUAUfPRg52RPq0XbDMR9AwO4NAGjmomR+kda0sixm8aZ+IhCVb9f7E2xTE4TOitHn2tRq9HX0cqjLza8gs21rV9gcEQSw5tYLPH2YHM2qKEvixyT4dNS+Bw3KxTi2scJvTR79ED3LYux3LfbDcZ6Q1xTF0VlO7kC01soQd4rnEGldoVVj4KqJkk0QYOE/ro++jBct8c8PSH9Yuiel+DISQOFgLKek/HAsQ6gUp68HTjt2DPh1yTZMKJt4XDhN7u0Udp1o0XJYkl7suh6GNyxT9Cs8EQJ+z/BCS5Uwg9wMvu6TBoAo91eYVenDrjnOnhn4nDhN720Ue5Rw/whvtS6rQLdyWbBmSIU/b/AwZPoLpZAq3ExEMSfOPX9KWSR7v4J/8y+OJIoY+2jx6gjAxed1/K9cl/h0pX1O9vMLRKQw0Ur4ezL461JdFl8AQWumdZv8tNr8famrjFYULvdd3EJpvb843XWhkq1jwVk/sbDC2y/xMrrHZOwLTMCUvOvPd4pvE6Cj0j4b0HoHRnrE2KS5wl9PWWHy4q8egDUEx/lrgvg3UvwZFtMbHBYAjIl+9bY1fDOpfQA7hJ5gf191ppE1+9Hk7GV572eMBZQm8nHThJ7PJg/nfjDRz3dIXlD5Ez792Y2WEwNKEKX65gjWccOY90zgB8JfSFm16n5vhheO1GelAda5PiCocJ/QmqNQ03yTEzoYKe/LzxRij6mDuS34+ZHQZDE4e3QMUBPvRMjLUlMSVn4WG+33AvlGziD6kLoMY/kVjnxWFCXxXT3ryXxe7LYcRV/DjlNTi0MdbmGDo7mxZDUhcrt2onZ6VnEsx+hTGyDxZda+Lh2DhM6E9QpbEXehCY+RvKyIDXZsPx/bE2yNBZcTdQ9s8/wsirqaBnrK2JD/Kv5a6GB6FiP+VPXWQCEuI4oa+iKg569AB078sd9Q9DYy28MgPK98TaIkMn5P5HHqGfnICJt8balLghZ957/N1zLpefeJTj2tP6fS77EdSdjLVpMcNxQn8yLnr0VmPapVlwy9tQe4KyZy6FfX+PtVmGzoTHzb+lLGW7JxuGfz3W1sQde3UwM+p/ykuN06HwBVg4yQpj4m6MtWlRx1lCX3uCk3SLtRVnknU+3PlXKrU7vDwDls8zOWYN0WHdi4xIKuZ/Gr9pkmW3QDXp/Gfj7fCdv/J5ZXf4yw/gf86DT57pVL9TZwn9qVLKtVesrTiDnHnvQb/hXFv/M5h8J3z2WyqfHAMf/QwqDgYubzB0lKPb4YPH+dg9lmWeqbG2Jv4ZOpXr6v+T79Y/AL2GwMpH4L9Hwms3woZXE37ufVBCLyJXichOEdktIvMCnE8Tkdft85+JSI7PuR/bx3eKyPR2W+puhOpyawA0DqkhnZyPL2d63X9R6MnH87dfwK/Hw4tXwcf/zXU/fqppwZchuqxYsYKRI0eSl5fHggULmp1vT/tt6zcRUQ5vgT9cB6nd+VHDvwKmN98WVgdLWOmZRM6X98I9n8CkO+HIVnjne/DL4fDMefDOXPjsBWsA9+RRa41CAiDaxoOISDLwJfB1wAWsBW5S1W0+Zb4HjFfVe0RkDvBNVb1RREYDi4EpwGDgA2CEqrpbut+kSZN03bp1zU+cKIFfjWJ+w7f5ozv+/ZFZUsr1yX/jiqTPGZdUBIBbhSIdyDkjxvLqDg83X3khdOtrZbZPz4T0TK58di1/ffDrkNwFktOs75Q0SOpivZ5LEiD2tvmBt4Xb7WbEiBGsXLmSrKwsJk+ezObNm7eq6lhvmVDbr31Zq7+JQLTYtv1RBfWApxE8blA31J6wevHb34GNr0H3/lxV9gN26NDQ/yiGJor+q4Br/30h716rcOCfHNvxd/rI6UHbU5pG935D+eRoKhdPHMcLn1dx99cn2L/ZDEjrxU2vfMHie6ZBcqr1W01Os75T7N+vJLXxCc/vWETWq+qkQOdSgrh+CrBbVffalS0BZgK+jXom8Li9/RawUETEPr5EVeuAfSKy267vnyE/xe4PANijT/bnVwAABOVJREFUg0O+NBa4tD9PN17P01xP0X9M4q7/ep4xUsSopAOcc6qUq5L3wEcfNrvur2nA/zwc2s18xb/pPwHfY7EkRve/7N/honspLCwkLy+P3Fwrh+qcOXPYvHlzpl/pUNsvtP2bCJ5XZsKBzyxB9zRaIt8SXbpZM2wun8+On5g47B1GhC80Fy6+Bi7+PufNe5f+VDAiycUIcZElZdw5IJX00m1QtIZbkkth1Znu18WpwIs/66AdPsLv/c2c8dv1OfbDLdC9X0jVByP0QwBfZ7ML8HcKNpVR1UYRqQT62sc/9bt2iP8NRORu4G5796SItBSZqB88XBaEzXGFPEk/wHF22zjU9vv6wX1lQG+gl4h4Fzv0AQb5FW5P+23rNwGE2raD+TufAH5tf+ICh7YPy2550trxfgPsB3zfu+6KolFB0I/5/Vv6e5/d0kXBCH2gLpm/v6elMsFci6q+ALzQpiEi61p6NYlnnGo3ONd2r90icgMwXVXvso/fyuleeVPxAFW01n4DjW0F9IGath2fdDa7gxmMdQHZPvtZgH+apaYyIpICZADHgrzWYIgkkWi/pl0bHEUwQr8WGC4iw0QkFZgDLPUrsxS43d6+HvhIrVHepcAce1bDMGA4UBge0w2GoIhE+w2mToMhbmjTdWP7LO8F3geSgRdVdauIPAGsU9WlwO+BP9iDVcewGj52uTewBqkagbmtzbgJgjZfgeMUp9oNzrX9BYhc+w1UZzjsdSDG7ujSLrvbnF5pMBgMBmfjrJWxBoPBYAgZI/QGg8GQ4DhG6GO65DwERORFETkqIl/4HOsjIitFZJf93TuWNgZCRLJFZJWIbBeRrSLyA/t4XNsuIukiUigim2y7/9M+PswOZ7DLDm+QGmtbA+GUdg2mbUebcLZtRwi9HYbhWeBqYDRwk708PR5ZBFzld2we8KGqDgc+tPfjjUbg/6lqPnABMNf+G8e77XXA5ap6LjABuEpELgCeBJ6y7T4O3BlDGwPisHYNpm1Hm7C1bUcIPT5hGFS1HvAuOY87VPXvWDM3fJkJvGxvvwzMiqpRQaCqJar6ub1dBWzHWgUa17arhTc4SRf7o8DlWOEMIA7ttnFMuwbTtqNNONu2U4Q+UBiGZqEU4pgBqloCVqMDzoqxPa0iVvTGicBnOMB2EUkWkY3AUWAlsAeoUFVvhol4bS9Ob9fggPbhS2dt204R+qBCKRg6joj0AP4E/FBVT8TanmBQVbeqTsBaoToFyA9ULLpWBYVp11GkM7dtpwi905ecHxGRQQD2d1xmORCRLlg/hFdV9f/sw46wHUBVK4DVWH7Y/7+9O0ZpIAjDMPz+ldgtghfIATyBZapcwiLHEDxCrmFhlQMEjxCRgIWx9hApxmKmSLFJIZKdGd8HlixJ8xUfP8nC/BnKOgOoty+t9xoa6cd/73Yrg771I+fHR+wfgPWEWUZFRJBPiH6klFZHH1WdPSJuI2Io99fAnPwM9pW8zgAqzF203muovB9gtwFIKTVxAQvynz18AY9T5zmT8xn4Bg7kb2xL8srbDfBZXm+mzjmS+578E/AdeCvXovbswB2wLbl3wFN5f0beS7MHXoCrqbOeyN9Er0tWu33Z3H/WbVcgSFLnWnl0I0n6JQe9JHXOQS9JnXPQS1LnHPSS1DkHvSR1zkEvSZ37Aexd67T5z7wWAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 2 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAsIAAAF1CAYAAADiNYyJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdeVhUZf8G8PucmWFHUAFxBVwAEXHftcI9F7TUNDWXbDF3s90WyjL1Tc2l9c3SrKxfamnuS5pibrjjrggqKsgOss3MOb8/KF5RlG3gmeX+XBdXDHPmnFt6fbl7+M5zpN8jolUQEREREdkYWXQAIiIiIiIRWISJiIiIyCaxCBMRERGRTWIRJiIiIiKbxCJMRERERDaJRZiIiIiIbBKLMBFRGTw/uAtOHI6olGsdO7QXs998sVKudbcNvy7H91/MrfTrEhFVFhZhIqJKVJYC/cNX/8GgEeMLHr89eThG9WuNp3uGYNroPji4d/sDX3vq6H68PXk4hvcKwfODu9z3fPzN63h78nA81S0IE4d3L5StZ9jT+GvbOqSmJJYqLxGRpWARJiIyYxfPnkDWnQwEBLco+NpzU9/Bd78fwKptJzHhtdlY+MHLSE5MKPL19g5O6NZ3CEZPeLPI5xeET0V9/yCs3HQUI16YgbnvTERaShIAwM7eHi3bP4rdm38z/R+MiMgMsAgTEZXRxXMnMWlkT4zo3RyLZ7+KvNxcAMDhfTsxbUxfDO/dDK+PH4yYS2cBAAtnvYzE+Bv46PXnMaxHMNb++BUAYN7bEzEmrC2G9wrBWxOH4mr0hYJrHD3wF5o0b1four4NG0Oj1eY/kCQYjXokJtwsMqN/UDOE9n4C3rXq3vdc3NVoXL5wGk+PmwZ7ewd0fOxx+NQPwP6/thQcE9yiHSL37yr7N4mIyIyxCBMRldGebevx3vzl+PL/duPGtRj8umIpLp+PwtKP38CEVz/Cyo1H0WvA0/jojRegz8vF9HcWwKNGLcyc+1/8vD0KT47In/tt2f5RfPHzn1jxx2HU92+CBR9ML7hG7OXzqF2v/n3X/vC1cRjSNRCvvfAEgpu3R8PApqXOf+3KRXjXqgtHJ5eCr/k1bIyrVy4WPK7j07CgyBMRWRut6ABERJaqz6Bn4FmjFgBg8KgJ+O/C95GRnoKeA56Gf5PmAICujw/C6u8/x/nTxxHcol2R5+ne76mCz4c9Ow0jH2+OO5npcHapgjuZ6XB0cr7vNW/PWwaDQY8Th/ch7uplyHLp1zWys7Pg5Oxa6GtOzq5ISowveOzo5IysOxmlPjcRkSVgESYiKiMPr5oFn3vVqI2UxHgk3IrDn5vXYuOaFQXPGfR6JN9VLu9mNBrx49efYN+uzUhPTYYkSwCA9LQUOLtUgbOrG7Kz7hT5Wq1Wh1YdHsOGX5fDu7YP2nbuXqr8jo5OyMrKLPS1rKzMQsU7O+vOfWWZiMhasAgTEZXR3XO5t+NvoKpHDXh41cKQURMxZPTEIl8jSVKhx3u2r8ehiB344NOV8KpZB3cyMzDy8eaAqgIAfBsE4sa1Kw/NYTQacCsuttT56/o1QvyNq8jOyiwYj4i5dBaP9AgrOOZ67CX4Nmxc6nMTEVkCzggTEZXRprUrkZhwExnpqVi98nN07tYXPcOGYsu6n3Dh9HGoqoqc7CxE/v0nsv9ZeXWv6oFbN64VnCM7KxNanR1c3dyRm5ONH776T6FrtOrwGE4fO1jw+HrsZRzZvxu5uTkwGPTYvfV3nDlxuOANdfE3r2Ng5/qIv3kdAKAoCvJyc2EwGABVRV5uLvT6PABA7Xr14dcwCD9/uxh5ubk48NdWxFw+hw6P9i64XtSxQ2jZ/tGK+QYSEQnGFWEiojJ6pEcYwl8ejeTEeLTr3ANPjZ4EewdHTHxtNr5e+B5uXI+Bvb0DGjdtjSbN2wIABj3zEv67MBwrvpiDp0ZNQq+Bw3H80F48O7AjXKu4YfhzL2PL7z8WXKNBQDCcXFxx4fRx+DdpDlVV8fO3i3Dt3cmQNTJq1fHFK+8vRoOAYABAYsINeHrXRnXPGgCA08cP4Z0pwwvO91S3xmjSvB0+WroKADDj/cVY/NGrGPl4c3jUqIXXZ30Gt6rVAQB5ubk4emA35i9bVynfTyKiyib9HhGtig5BREQPduzQXmz+7Qe89fFXxR77f8uXws29GnoNHF7sscXZsHoFEhNuYsyEN8p9LiIic8QiTEREREQ2iTPCRERERGSTWISJiIiIyCaxCBMRERGRTWIRJiIiIiKbJGz7tHED2sDX11fU5YmIiIjIRly4FI2VG4/c93VhRdjX1xeRkZGiLk9ERERENqJhYNMiv87RCCIiIiKySSzCRERERGSTWISJiIiIyCYJmxEmIiIisiR6vR7Xr19HTk6O6Cj0AA4ODqhTpw50Ol2JjmcRJiIiIiqB69evw9XVFb6+vpAkSXQcuoeqqkhKSsL169fh5+dXotdwNIKIiIioBHJyclC9enWWYDMlSRKqV69eqhV7FmEiIiKiEmIJNm+l/ffDIkxERERkhdavX485c+ZU+nX79OmD1NTUSr9uWXBGmIiIiMgKhYWFISwsrNKvu2nTpkq/ZllxRZiIiIjIAsTExCA4OLjg8SeffILw8HAAwOLFixEUFISQkBAMGzYMALB8+XJMmjQJADBmzBhMmTIFHTt2RP369bF69WoAgKIomDBhApo0aYJ+/fqhT58+Bc/drajzZ2ZmYuzYsWjatClCQkKwZs0aAPl3D05MTERMTAwCAwMxevRohISEYPDgwcjKysLOnTvxxBNPFJx7+/btePLJJ03/DSsBrggTERERldL7f5zGmRvpJj1nUK0qeK9/kzK9ds6cObhy5Qrs7e0fOJZw8+ZNRERE4Ny5cwgLC8PgwYOxdu1axMTE4NSpU0hISEDjxo3x7LPPluj8s2bNgpubG06dOgUASElJue9158+fx7Jly9CpUyc8++yz+PzzzzFjxgxMnDgRt2/fhqenJ7777juMHTu2TH/u8uKKMBEREZGFCwkJwYgRI/DDDz9Aqy16nXPgwIGQZRlBQUGIj48HAERERGDIkCGQZRne3t4IDQ0t8fl37NiBiRMnFhxTtWrV+15Xt25ddOrUCQAwcuRIREREQJIkPPPMM/jhhx+QmpqK/fv34/HHHy/Xn7+suCJMREREVEplXbktD61WC0VRCh7fvU3Yxo0bsWfPHqxfvx6zZs3C6dOn73u9vb19weeqqhb6Z3GKOr+qqsXu0nDv8/8+Hjt2LPr37w8HBwcMGTLkgeW9onFFmMxfuFvhDyIiIhtUo0YNJCQkICkpCbm5udiwYQOA/Dnfa9euITQ0FPPmzUNqaioyMzNLdM7OnTtjzZo1UBQF8fHx2L17933HPOj8PXv2xNKlSwuOK2o04urVq9i/fz8AYNWqVejcuTMAoFatWqhVqxY+/PBDjBkzppTfCdNhESYiIiKyADqdDu+++y7atWuHfv36ITAwEABgNBoxcuRING3aFC1atMD06dPh7u5eonMOGjQIderUQXBwMF588UW0a9cObm6FF50edP63334bKSkpCA4ORrNmzbBr1677zt+4cWOsWLECISEhSE5OxksvvVTw3IgRI1C3bl0EBQWV47tSPtLvEdElWxM3sVlThyAyMlLEpcnS3LsKHJ4mJgcREdm0s2fPonHjxqJjmFxmZiZcXFyQlJSEtm3bYt++ffD29i73eWNiYtCvXz9ERUUV+fykSZPQokULjBs3rtzXultR/54aBjbF/GXr7zuWM8JERERENqxfv35ITU1FXl4e3nnnHZOU4OK0atUKzs7OmD9/foVf62FYhMn8FDcHXNTzXCUmIiIqk6Lmgk3B19f3gavBR44cqZBrlhaLMFknjlMQERFRMfhmOSIiIiKySSzCRERERGSTWISJiIiIyCaxCBMRERFZofXr12POnDnCrt+nTx+kpqYCABYvXozGjRtjxIgRwvIUhW+WIyIiIioLU9/t1MRv7A4LC0NYWJhJz1kamzZtKvj8888/x+bNm+Hn51ei1xoMhkq57TJXhImIiIgsQExMDIKDgwsef/LJJwgPDweQv+IaFBSEkJAQDBs2DACwfPlyTJo0CQAwZswYTJkyBR07dkT9+vWxevVqAPm3T54wYQKaNGmCfv36oU+fPgXP3e2xxx4ruBFaYmIifH19C67x5JNPonfv3mjUqBFee+21gtf4+voiMTER48ePR3R0NMLCwrBw4UIkJydj4MCBCAkJQfv27XHy5EkAQHh4OF544QX07NkTo0aNwvLlyzFw4ED0798ffn5+WLp0KRYsWIAWLVqgffv2SE5OLvf3lCvCZLZUFYhWa2K/EoQDShCi1ZrwklJRU0pGTSkJ3sj/Zz0pAT6iwxIREQk0Z84cXLlyBfb29gXjCPe6efMmIiIicO7cOYSFhWHw4MFYu3YtYmJicOrUKSQkJKBx48Z49tlnS3Xt48eP49ixY7C3t0dAQAAmT56MunXrFjz/5ZdfYsuWLdi1axc8PDwwefJktGjRAr///jv+/PNPjBo1CsePHweQv79wREQEHB0dsXz5ckRFReHYsWPIyclBw4YNMXfuXBw7dgzTp0/H999/j2nTppX9mwYWYTIzKXfysMUQWlB+E1AVAFADyQiUryJRdUOU4odEFP51VOsv/sa4zn7o2cQbGlmqmHDcm5iIiMxUSEgIRowYgYEDB2LgwIFFHjNw4EDIsoygoCDEx8cDACIiIjBkyBDIsgxvb2+EhoaW+trdunWDm1v+z8igoCDExsYWKsL3ioiIwJo1awAAXbt2RVJSEtLS8n+mhoWFwdHRseDY0NBQuLq6wtXVFW5ubujfvz8AoGnTpgUryeXBIkxmQVVV/H48Dh/8cQYphufhiRR0kM+gvXwWHeQz8JVuQbqr3+aqWsSr1XAT1XBKqY8VGc/jpR+Pom41R4zp6IenVEe4Stni/kBEREQmptVqoShKweOcnJyCzzdu3Ig9e/Zg/fr1mDVrFk6fPn3f6+3t7Qs+V1W10D9Lc+27r3vveTUaDQwGw0PPVdQ1pX9+yDs7Oz/w3LIsFzyWZbnY65QEZ4RJuOvvNcDYmbMx/ZcT8Mk+gz/sZuKQ/UQstvsMw7V/wk8uXIIBwF4yoJ6cgHbyOTyn3YTdr4Tiy5EtUcPVAbM2nEHH3CX4SD8cqapz0RclIiKyMDVq1EBCQgKSkpKQm5uLDRs2AMif87127RpCQ0Mxb948pKamIjMzs0Tn7Ny5M9asWQNFURAfH//A2y37+voW3Ba5qBni0njkkUfw448/Asi/vbOHhweqVKlSrnOWVbErwnm5uZg5aSj0eXkwGo3oGNobT4+bXuiYnZtWY8Xnc1DNowYAoO+gUejRf2jFJCarYVRUrNwfg3m58wAA72lXYJRmGzRSyf7r9G4aWULv4JroHVwTx6+lYtkX/8G3xsex0dgeS+yWoJWJsxMREVU2nU6Hd999F+3atYOfnx8CAwMBAEajESNHjkRaWhpUVcX06dPh7u5eonMOGjQIO3fuRHBwMPz9/dGuXbuCMYe7vfLKK3jqqaewcuVKdO3atVx/jvDwcIwdOxYhISFwcnLCihUrynW+8pB+j4h+aOtQVRU52VlwdHKGwaDHmy89heemvouA4BYFx+zctBqXz53CCy+/X+ILz5o6pODdh2R7Lt/OxKu/nsDRq6l4RD6Bj7TLUFdOLPsJ753XDXfDScUPk/RTEKd64JXeTfDiI/Uhl2d+mDPCREQ27ezZs2jcuLHoGCaXmZkJFxcXJCUloW3btti3bx+8vb1Fxyqzov49NQxsivnL1t93bLErwpIkwdEp/9fLRoMBRqOhYI6DqCwuxmfgqa/2QwWw4KlmeGLd8PtGH0whRL6CDXZv4U39c5i7RYP90UlY8FQzeLjYF/0CFl0iIrJB/fr1Q2pqKvLy8vDOO+9YdAkurRK9Wc5oNGLGuDDciovF40+MhH+T5vcds/+vLTh94hBq1fXDs5PfhmeNWiYPS5bvalIWRnxzEFqNjF9f7ABfD2fg/v9AM5kqUjaW6pagY5/n8f4fZ9Bn0V4sGtYCHRpUr7iLEhERWZAHzQXbghK9WU6j0eDT5Rvxzdq/cfHsScRGny/0fJtO3fD1r3uwaMVmNGvdCYs/erXI82xdtwozxoVhxrgw3L59u/zpyaLcSsvBiGUHkGdU8MO4dvkluBJIEjCinQ/WTewEFwctRnxzACv+jqmUaxMREZH5KtX2aS6uVRDcoh2OHdgDn/oBBV+v4la14PMe/Yfh+y/mFvn6XgOeRq8BTwPInxEm25GUmYuRyw4iOTMPPz3fHgHerpWeoXHNKvhjUmdM/+U43lt/Gtj0KkZrt1V6DiIislyqqnJE1IyVdDu4fxW7IpyWkoTMjHQAQG5uDk5E7kNtn/qFjklOTCj4/HDEDtTxaViqEGTd0nP0GP3dIVxLzsKyMW3QrG7J3slaEZzttfhsREv0kg/hPcMYfG/oISwLERFZFgcHByQlJZW6bFHlUFUVSUlJcHBwKPFril0RTklKwKKPXoWiGKEqKjp17YM2nbrhp28WomFgU7Tt3B0bVy/HoYid0Gg0cKnijikz/1OuPwhZj+w8I8YtP4xzNzPw31Gt0b6++NlcnUbGEt0STNRPxbuGsZCg4hntDtGxiIjIzNWpUwfXr1/neKcZc3BwQJ06dUp8fLFF2LdhYyz8bsN9Xx/+3P/2En5m/Gt4ZvxrJb4oWZGH7LSgqiomhc/GEaU5FuuWIPTng2azE4OdZMRnukWYoJ+GdwzPQoaCEdo/RcciIiIzptPp4OfnJzoGmRBvsUwV5qdDV7FTaYn3tCvQT3NQbJh7Czv+V4Zf0k/DTMNzkAAMZxkmIiKyGbzFMlWIa8lZmL3xLDrJURij2So6zgPZSwZ8ofsUofIxvGV4Dr8aHhEdiYiIiCoJizCZnKqqeGPtSQDAXN3XFXKzDFP6twx3lk/hLcNziFT8RUciIiKiSsAiTCb348Gr2HcpCTP7BqGOVI7bJpdGuFvhj1JykPT4TLcYtaVEjM+bhptqtQoISUREROaERZhM6lpyFmZvOosujTzwdNu6ouOUipt0B1/r5iMb9hifNx05eqPoSERERFSB+GY5MhlFlfDa6pOQJQlzBoVY5Ibj/nIcFui+wIv6lzHzvTfxie4rsx/tICIiorLhijCZzI/GbtgfnYSZfRujtruj6Dhl1ksTiamaNVijPIrvjL1FxyEiIqIKwiJMJnFN8cTHhuHo0sgDw9pY1khEUaZq16KHHImPDCOwz9hEdBwiIiKqACzCVG6qCrxmeAEyFIsdibiXLKlYoPsCftJNTNJPwTXFU3QkIiIiMjEWYSq3HUpL7Fea4A3tKoseibiXq5SN/+oWwAAZ4/XTkKdqREciIiIiE2IRpnJRVAnzDUPgK93CMM0u0XFMzk++hU90X+G06odPDYNFxyEiIiIT4q4RVC6blbY4p/rgU91n0EpK8S8owx6/ovXSROIpZRe+NPZHqOY42sjnRUciIiIiE+CKMJWZUVGxwDAYjaTr6C//LTpOhXpXuxK1pUS8rH8JmaqD6DhERERkAizCVGbrjsfhslobL2t/hUZSRcepUC5SDhbqPkec6oFZhmdExyEiIiITYBGmMtEbFSzaeRFBUgx6yZGi41SK1vIFjNf8gV+ModhmbCU6DhEREZUTizCVyZoj1xGblIUZ2l8hW/lq8N2maVcjSIrBm/rncDsjV3QcIiIiKgcWYWsV7lb4w4RyDUYs3nkRzeu6o6t8zKTnNnd2khGf6j5DBhzx5tqTUFXb+Y8AIiIia8MiTCX3T6n+JXwobqTlYMat12AF984oNX85Dm9of8aOswn45fA10XGIiIiojFiEqVRyVB2WGgairXQWneUo0XGEGaPZik4Nq2PWhjO4kZotOg4RERGVAYswlcoPxh5IQFXM0P1qk6vB/5IlFXOeDIFRVRG+/rToOERERFQGvKEGldgd1R6fG8LQRT6JdvK5og+ywBtmlFXdak6Y3t0fH28+h62nb6FXE2/RkYiIiKgUuCJMJbbW2AXJqIKp2rWio5iNZzv7IdDbFeHrTyMz1yA6DhEREZUCizCViKqq+MHYHcHSFbSSLoiOYzZ0GhkfPdEUt9JzsGAbvy9ERESWhEWYSuTQlWScV+thlGabTc8GF/LPLhqtvvPFCHk7lu+7jKi4NNGpiIiIqIRYhKlEvj8QCzdkor9mv+goZulV7S+ojjS8ufYUjAr3FiYiIrIELMJUrIT0HGyNuoWnNLvhKOWJjmOW3KQsvKtbiVNxafh+f4zoOERERFQCLMJUrJ8OXYVBUTFSs0N0FLPWTz6AR/098cnW87iZxr2FiYiIzB2LMD2U3qjgp4NX8ai/J3zkBNFxzJokAbNiRsCQl43wuXNsais5IiIiS8QiTA+17XQ8EjJyMaqDj+goFqGenIAp2rXYqrTFbmOI6DhERET0ECzC9FDf749BnaqOeCzAS3QUi/GcZhN8pVuYZXgGeqMiOg4RERE9AIswPdD5Wxk4eCUZI9v7QCNzz7SSspcMmKn9AZfV2li5P1Z0HCIiInoAFmF6oJUHYmCnlfFU67qio1ic7vJRdJZP4dMdF5B8hzttEBERmSMWYSpSRo4evx2NQ/+QWqjmbCc6jsWRJOAd7UrcyTNiwfbzouMQERFREbTFHZCXm4uZk4ZCn5cHo9GIjqG98fS46YWO0efl4tMPX8Hl81FwreKOVz5Ygho161RYaKp4vx2Lw508I98kVw4B8nWMaFcPPxyIxcj2Pgj0rnL/ThLhvBMdERGRKMWuCOvs7PDBoh/x6YpNWLh8A44e2IPzUccKHbN9w//BxbUKvvxlF8KGPovvv5hbYYGp4qmqiu/3x6JZHTc0q+suOo5Fm97dH64OOnzwxxmoKu84R0REZE6KLcKSJMHRyRkAYDQYYDQaIEmF3zh1KGIHQh8fBADo+NjjOHnkb/7Qt2AHryTjUkImnungKzqKxavqbIfp3Rvh78tJ2HYmXnQcIiIiukuxoxEAYDQaMWNcGG7FxeLxJ0bCv0nzQs8n346Hh1dNAIBGq4WTsysy0lJQxb1aoeO2rluFbetXAQD02emmyE8V4LejcXCx16JfSE3RUazCiPY++PHgVXy08SweU7WwlwyiIxERERFK+GY5jUaDT5dvxDdr/8bFsycRG134zT9Frv5K92+31WvA05i/bD3mL1sPT0/PsiWmCpWjN2JT1E30auINB51GdByroNPIeKdfEK4mZ+Fb4+Oi4xAREdE/SrVrhItrFQS3aIdjB/YU+np1L28kJtwEkD8+kXUnA65VOFtqiXafv42MHAMGNK8lOopVecTfE90CvbDUMBAJKm+9TEREZA6KLcJpKUnIzMgfY8jNzcGJyH2o7VO/0DFtO3XDrs1rAAB/796Mpi073DdHTJZh/Yk4eLjYoWOD6qKjWJ2ZfRsjFzosMgwSHYWIiIhQghnhlKQELProVSiKEaqiolPXPmjTqRt++mYhGgY2RdvO3dG931B8OutljB8aCtcqbpgRvrgyspOJpefoseNsAoa3rQethltMm1p9Txc8rfkTPxm74VnNZjSQb3I7NSIiIoGKLcK+DRtj4Xcb7vv68Of+t5ewnb09XvvwM9Mmo0q3NeoW8gwKxyIq0BTtWqw1dsF/DEPxpd2nouMQERHZNC77UYH1J26gXjUnNOfewRXGU0rH89qN2KK0xRGlkeg4RERENo1FmAAACRk52HcpEQOa1+J8dwV7XrMRHkjFHP3T4HbbRERE4pRoH2GyfhtO3ISiovBYxL3zq1R6RXwPnaVcTNOuwduGcdihtEQPzVEBwYiIiIgrwgQAWHfiBprUqoKGXq6io9iEoZrdqC/dwFzDMBhU/jUkIiISgT+BCTGJd3DiWirfJFeJdJIRr2l/wSW1DlYbHxEdh4iIyCaxCBPWHb8BSQL6N2MRrky95MNoKV3AQsNgZKt2ouMQERHZHBZhG6eqKtadiENb32qo6eYoOo5NkSTgTd0qxKMab71MREQkAIuwjTt9Ix3Rt+9gYIvaoqPYpDbyeXSXI/GFoT+SVM5nExERVSYWYRv3+7E46DQSHg/2Fh3FZr2h/RlZcMDnhgFlO0G4W+EPIiIiKhEWYRtmVFT8cfIGHvX3grsTZ1RFaSjfwJOavVhp7I6badmi4xAREdkMFmEbdvBKEuLTc7lbhBmYqlkDFTIWz5vJFV4iIqJKwiJswzaduglHnQbdG9cQHcXm1ZUT8bTmT/xqfBSxipfoOERERDaBRdhGqaqKHWcS8Ii/BxztNKLjEIBJ2t+hhRGfGgaJjkJERGQTeItlGxUVl45b6Tl4JSjgf1/kr+GF8pJSMVqzFV8b+2G88gcC5OuiIxEREVk1rgjbqO1nbkGWgK6B/DW8ORmv3QAX5GCBYbDoKERERFaPRdhGbTsTj9Y+1VDNmbtFmJOqUiae027EVqUtTij1RcchIiKyaizCNuhachbO3cpAjyC+Sc4cPavZgqrIwCeGp0RHISIismoswjZox9l4AEB3FmGz5Cpl4yXteuxVQnBACRQdh4iIyGqxCNugHWfj0dDLBX4ezqKj0AOM0mxDDSTjE/1QqKoqOg4REZFVYhG2MWnZehyMTubewWbOQdJjkvZ3RKoB2H3+tug4REREVolF2MbsPp8Ag6JyPtgCDNXsQl0pAQu2X+CqMBERUQVgEbYx28/Ew8PFDi3quouOQsWwk4yYrPkNp+LSsONsgug4REREVodF2IbkGRT8df42ugXWgCxLouNQCTyp2Quf6k5YsP0CFIWrwkRERKbEO8vZkINXkpCRa/jfWATvJGf2tJKCqd0a4eX/O4FtZ26hd3BN0ZGIiIisBleEbcj2M/Fw0Mno1NBDdBQqhbBmtVDf0xkLt1/kqjAREZEJsQjbCFVVseNMPLo08oSjnUZ0HCoFrUbG1G6NcD4+A5uiboqOQ0REZDVYhG3E6RvpuJGWw90iLFS/kFpo5OWCT3dchJGrwkRERCbBImwjtp+JhyQBXQO9REehMtDIEqZ198elhExsOHlDdBwiIiKrwCJsI3acjUerelXh4eLOe3MAACAASURBVGIvOgqV0ePB3gj0dsWiHRdhMCqi4xAREVk8FmEbEKdWx+kb6ejOsQiLJv+zKhydeAfrjnNVmIiIqLxYhG3ATmNLAOB8sBXo1aQGmtSqgsV/XoSeq8JERETlwn2EbcBOpQX8pJto8Flt0VGonCRJwvTu/nju+0isPXodQ9vUEx2JiIjIYnFF2MrlqDocVBrjUfmE6ChkIt0aeyGkjhuW/HmJq8JERETlUOyK8O34G1j04StITb4NSZLRM2wY+j81ttAxp44ewMdvvgCvmnUBAB0e7YWhY6dUTGIqlSOKP3Jgjy7yKdFRyET+XRUeu/xw/qqw6EBEREQWqtgirNFoMXbSW2gQEIzsrEzMeDYMzdt0Rl2/RoWOC2rWBm/PW1ZhQals9ihNoYMB7eUzoqOQCT0W4Ilm/6wKP6lqoJOMoiMRERFZnGJHI6p5eKFBQDAAwNHJBXV8GyIp8VaFByPT2KOEoJV8Hs5SrugoZEKSlL+DxPWUbKw1dhEdh4iIyCKV6s1y8TevI/rCafgHNb/vufNRxzBtdB9U86iBMRPfRL36/vcds3XdKmxbvwoAoM9OL2NkKqkE1Q1nVV+8plklOgpVgIJV4biBeFKzl6vCREREpVTiN8tlZ93B3JkTMG7qO3Bydi30XIOAJvh69V58umIT+gwehY/ferHIc/Qa8DTmL1uP+cvWw9PTs3zJqVj7lPyV/Ec4H2yVClaFVS+uChMREZVBiYqwwaDH3Lcn4NGeYejwaO/7nndydoWjkzMAoHWHUBgMBqSnJps2KT1cuFvhDwB7jCGojjQESbGCw1FFeSzAE82ky1hiHAi9qhEdh4iIyKIUW4RVVcXSj99AHZ8GGDDsuSKPSUm6DVVVAQAXzpyAqihwdatq2qRUKooqYa/SFJ3lKMiSKjoOVRBJkjBNu4arwkRERGVQ7Izw2ZOR2L31N/g0CMC0MX0BACNffAWJ8fm3eO09cAT+3r0ZW377ERqNBnb2Dnjl/cWQJKlik9NDnVPrIhHu6KI5KToKVbDH5ONoJl3CEuM/s8KiAxEREVmIYotwULM2+D0i+qHH9B00Cn0HjTJZKCq/PUoIAHD/YBsgScA07RqM1b+OtcYu3FeYiIiohHhnOSu1VwlBoHQVNaRU0VGoEjwmnyhYFebd5oiIiEqGRdgKZat2OKwEoIvMsQhb8e+q8HXVC2uPXhcdh4iIyCKwCFuhA0pj5EHHsQgbU7Aq/OclrgoTERGVAIuwFdqrNIU98tBWPic6ClWiglXhlGyuChMREZUAi7AV2quEoK18Dg6SXnQUqmSPySfy7zbHVWEiIqJisQhbmRtqNVxU6+ARzgfbJElC/t3mUrLx29E40XGIiIjMGouwlYkwNgUAFmEb9liAJ0LquGHJrotcFSYiInoIFmErs0cJgRdS4C9xRtRWSZKEad0b4VryQ1aFi7glNxERka1hEbYiRlVChBKMLvJJ8MZ+ti00wAshddywdBdnhYmIiB6ERdiKRKl+SIUrHuFtlW1buBuk990xNf5tXE3Owm/vDRSdiIiIyCyxCFuRvUr+fHBnOUpwEjIHXeVjaCpF4zPjABi4KkxERHQfFmErsscYgmDpCqpLGaKjkBmQJGCqdi1iVW/8dow7SBAREd2LRdhKZKt2OK42RCeuBtNduslHESxdwdJdl7gqTEREdA8WYStxVGmEPOjQXj4jOgqZkX/vNheblIXfj98QHYeIiMissAhbif1KEDQwoo18XnQUMjPd5KMIrl0FS/68yFVhIiKiu7AIW4kDShCaStFwkXJERyEzI0nA1G7+XBUmIiK6h1Z0ACq/rDwDTqgNME6zSXQUqggmuOFF98ZeaFIrf1V4YPNa/ItPREQErghbhSOxKdBDiw6cD6YHyL/bXP6qMHeQICIiyscibAX2X06CFga05nwwPUi4G7r/4p+/g8SabTCo/KtPRETEn4ZW4EB0EkKkaDhLuaKjkBkr2EFC9cZvxs6i4xAREQnHImzh7uQacPJ6GsciqES6yUfRVIrGEuMT0Ksa0XGIiIiEYhG2cIdjkmFQVO4fTCXy76rwVbUGV4WJiMjmsQhbuAPRydBpJLSSL4qOQhaiq3yMq8JERERgEbZ4+6OT0KyOO5w4H0wl9O+q8DXVi6vCRERk01iELVhGjh5RcWno0KC66ChkYbrKxxAiXeaqMBER2TQWYQsWGZMCo6KifX0WYSqdu1eF1xq7iI5DREQkBIuwBdsfnQQ7jYxWPlVFRyELFCof/2dVeCD0RkV0HCIiokrHImzBDkQnoXk9dzjo+KttKr1/V4Wvq15Yc+S66DhERESVjkXYQqX/Mx/MsQgqj1D5OJpLF7Hkz0vINRhFxyEiIqpULMIW6vCVZCgq0IFFmMpBkoCXtasRl5qN/4vkqjAREdkWFmELtf9yEuy0MlrUcxcdhSxcF/kUWvtUxWd/XkKOnqvCRERkO7SiA1AJhLvd96UD1TegJeeDyQQkCXi5hz+Gf3MQPx+6ijGd/ERHIiIiqhTFrgjfjr+BtycPx6QRPTB5ZC/88X/f3XeMqqr476fvY/zQUEwd/Tgun4+qkLCUL011xukb6ZwPJpPp0KA62vlVw2e7L3NVmIiIbEaxRVij0WLspLew9MftmPf1GmxeuxLXrhS+ne+RA7tx81oMvvj5T0x4dTa+/OSdCgtMwEElECrng8mEJEnC9B7+uJ2Rix8OxIqOQ0REVCmKLcLVPLzQICAYAODo5II6vg2RlHir0DGH9u7AY72fgCRJCAhugTuZ6UhOTKiYxIQDSmPYa2U053wwmVD7+tXRqWF1fPnXZWTlGUTHISIiqnClerNc/M3riL5wGv5BzQt9PTnxFjy8ahY8ru7ljeR7yjKZzn4lCK18qsJey/lgMq2Xe/gjMTMPK/dzVZiIiKxfid8sl511B3NnTsC4qe/Aydm10HOqqhbxCum+r2xdtwrb1q8CAOiz00uXlAAAaaoTzqn1MM2PYxFkeq18quFRf098+ddljGjvAxd7vp+WiIisV4lWhA0GPea+PQGP9gxDh0d73/d8dc+aSEy4WfA4KeEWqnnUuO+4XgOexvxl6zF/2Xp4enqWI7btOqL4Q4WMNn+Nyt9N4t8PIhOZ3sMfKVl6rPg7RnQUIiKiClVsEVZVFUs/fgN1fBpgwLDnijymbedu2L3lN6iqivNRx+Ds4opqHl4mD0vAISUQOhjQQr4kOgpZqeZ13dEt0Atf74lGeo5edBwiIqIKU+zvPc+ejMTurb/Bp0EApo3pCwAY+eIrSIy/AQDoPXAEWnUIxZH9uzF+aCjsHRww5a15FZvahh1WAhAsXYGjlCc6Clmx6T380W9JBJbtvYLpPfxFxyEiIqoQxRbhoGZt8HtE9EOPkSQJL874wGShqGg5qg4n1QZ4VrNZdBSycsG13dC7iTeWRVzBmI6+qOpsJzoSERGRyfEWyxbkhNoAemjRRj4vOgrZgJd7+uNOngFf/nVZdBQiIqIKwSJsQQ4rAQCA1izCVAn8a7hiYPPaWLE/BgnpOaLjEBERmRyLsAU5pATCX7oGd+mO6ChkI6Z1bwSDUcXSXXxzJhERWR8WYQthVCUcVRpxLIIqlU91ZwxpXRerDl3FteQs0XGIiIhMikXYQpxV6yETTmgrnxMdhazR3XtS37Mv9ZRuDSFJEhbvvCgoHBERUcVgEbYQh5VAAEAbFmGqZDXdHDGynQ/WHL2Oy7czRcchIiIyGRZhC3FYCURt3EYtKVl0FLJBE0IbwEGnwcLtF0RHISIiMhkWYQugqsAhJYBjESSMh4s9xnbyxYaTN3HmRrroOERERCbBImwBYlRvJMKd26aRUC90aQBXBy0WbOf/DomIyDqwCFuAf/cPbssiTAK5Oenw4iP1seNsAo5eTREdh4iIqNxYhC3AYTUAVZGBhlKc6Chk48Z28oOHix3mbTkHVVVFxyEiIioXFmELcFgJRGv5PCRJdBKydc72WkwKbYgD0cnYczFRdBwiIqJyYRE2cwkZOYhRvflGORLvnz2Gh29rhbpSAuZuPgdF4aowERFZLhZhM3f4Sv4sJu8oR+bCTjJihvZXnLmZjj9O3hAdh4iIqMxYhM3c4ZhkOCIHTaQY0VGICoTJfyPQ2xXzt11AnkHJ/+JD7k5HRERkjliEzdzhmGS0kC9BJxlFRyEqIEsqXu8diKvJWfjl8FXRcYiIiMqERdiMZeTocfZmOtpIHIsg8/NYgCfa+lXDop2XcCfXIDoOERFRqbEIm7EjsSlQVPCNcmSWJEnC670DkZiZi28jroiOQ0REVGoswmbscEwytLKEFvIl0VGIitTKpyp6BNXAV3uikay6io5DRERUKizCZuxwTAqa1HaDk5QrOgrRA73WKwBZeQZ8bggTHYWIiKhUWITNVK7BiOPXUtHGp6roKEQP1aiGKwa1rIPvjT0Rp1YXHYeIiKjEWITNVFRcGvIMClr7VhMdhahY03v4A1CxQD9YdBQiIqISYxE2U5Ex+TfSaMUVYbIAtdwdMVazFWuVLjit+IiOQ0REVCIswmYqMjYFvtWd4OlqLzoKUYlM0K6DG+5gtmEEVN55mYiILACLsBlSVRVHY1PQyodjEWQ53KQsTNWuxT4lGLuVZqLjEBERFYtF2AzFJGUh6U4eWvtyLIIsywjNDvhKt/CxYTgMRkV0HCIioodiETZDkTHJAIDWnA8mC2MnGfG6dhUuqHWx+sh10XGIiIgeSis6AN3vSGwK3Bx1aODpIjoK2apwtzK/tLd8GK2k85i/3R79m9WCsz3/b4aIiMwTf0KZocjYFLTyqQpZlkRHIXqwB5RlSQJm6n7EkxkB+HpP9D9bqxEREZkfjkaYmdSsPFxKyOS2aWTRWsqX0DekJr7eE42E9BzRcYiIiIrEImxmjsRy/2CyDq/3CoRBUbBg+wXRUYiIiIrE0QhzcNevmCP1Q6GVB6JZHXeBgYjKr151J4zq4Ivv9l3B2E5+CPB2FR2JiIioEK4Im5kjij+aqBfhOLtafkEux5uWiESb3LUhXOy1+GjTWai8ywYREZmZYovwktmvYXS/NpjyTO8inz919ACG9wrBtDF9MW1MX/zy3WKTh7QVeaoGJ9QGaC3zV8lkHdyd7DC1uz/2XLiNXecTRMchIiIqpNjRiK59BqPPoFFY9OErDzwmqFkbvD1vmUmD2aIo1Q+5sENr+bzoKEQmM6qDD346GIsPN5xF54aesNPyF1FERGQeiv2J1KR5W7hU4bxqZTii5G8z1Uq+KDgJkenoNDLe7heE6MQ7+H5/jOg4REREBUzyZrnzUccwbXQfVPOogTET30S9+kXvG7p13SpsW78KAKDPTjfFpa1KpOKPulICvKRU0VGITCo0wAuhAZ5YtOMiBraoDQ8Xe9GRiIiIyv9muQYBTfD16r34dMUm9Bk8Ch+/9eIDj+014GnMX7Ye85eth6enZ3kvbVVUNX9FuLXEsQiyTm/3C0K23oj52zgDT0RE5qHcRdjJ2RWOTs4AgNYdQmEwGJCemlzuYLbmquqFRLijFd8oR1aqgacLRnf0xc+Hr+L0jTTRcYiIiMpfhFOSbhdsi3ThzAmoigJXN94MorQi1QAA4I4RZNWmdGuEqk52eP+PM9xOjYiIhCt2Rnj+e1MQdfwg0lNTMO6Jjhg2biqMBgMAoPfAEfh792Zs+e1HaDQa2Nk74JX3F0OSpAoPbm0iFX+44g78peuioxBVGDdHHWb09MfM36KwOeoW+jStKToSERHZsGKL8Iz3H74vcN9Bo9B30CiTBbJVRxR/tJQvQpa4SkbWbVibeli5PxYfbTyLroFecNBpREciIiIbxQ09zUCa6owLal2ORZBN0MgS3uvfBHGp2fhmb7ToOEREZMNYhM3AUaURAKCVxCJMtqFDg+ro09Qbn+26jLjUbNFxiIjIRrEIm4FIxR8aGNFcviw6ClGlmdk3CKo+C7PmzQHC3fI/iIiIKhGLsBmIVPwRJMXCScoVHYWo0tR2d8Rk7W/YorTFbmOI6DhERGSDWIQF0xsVnFAbcP9gsknPazaivnQD7xnGIEfV/W9lmCvERERUCViEBTt9Ix05sEdrmXeUI9tjJxnxgXY5YlVvfG3sJzoOERHZGBZhwSJj8u/Cxx0jyFZ11kShr7wfnxkG4JrCW68TEVHlKXYfYapYR6+moDZuw1tKER2FyLSKGm0IL/rWyu/ofsDu3OYIN4zGMrtPKjgYERFRPq4IC6SqKiJjUrgaTDbPW0rBNO0a7FRaYruxpeg4RERkI1iEBbqeko2EjFzOBxMBGKPZCn/pGsL1o5Gt2omOQ0RENoBFWKDI2Pz5YO4YQQToJCM+0H2HOHjic8MA0XGIiMgGsAgLFBmTAhd7LQKka6KjEJmF9vI5DJQj8JWxHy4ptUTHISIiK8ciLNCR2BS0qOcOjaSKjkJkNmbqfoAjcvGm/jkoCv9uEBFRxWERFiQtW4/z8Rlo5VNVdBQis+IppWOm9kccVgOx6vBV0XGIiMiKsQgLcuxqClQVaO1TTXQUIrMzRPMXOsinMWfTOcSn54iOQ0REVopFWJAjsSmQJaB5PXfRUYjMjiQBs7XLkGdU8N6606LjEBGRlWIRFiQyJgWNa1aBiz3vaUJUFD/5FqZ2b4Qtp29h6+lbouMQEZEVYhEWQG9UcPxaKlpzPpjooZ7vUh+Na1bBu+uikJ6jFx2HiIisDIuwAGdvpiNbb0QrX84HEz2MTiNjzpNNcTsjF/O2nBMdh4iIrAyLsACRMSkAwBVhohJoVtcdYzr64YcDVxEZkyw6DhERWREOqApw5GoKark5oJa7o+goRJUr3K1ML5vR0x9bT9/CG2tPYeOUzrDXakwcjIiIbBFXhCuZqqo4EpPCsQiiUnC21+LDJ4JxKSETi3deFB2HiIisBFeEK1lcajZupedwLIKoJO5aQQ4FMETzAr7YDfQM8kazutx6kIiIyocrwpXsSGz+fDDvKEdUeu9oV6JGFQfM+PUEcvRG0XGIiMjCsQhXssiYFDjbaRDo7So6CpHFqSJlY+6gEFxKyMSC7RdExyEiIgvHIlzJImNT0KJeVWg1/NYTlcUj/p4Y3q4e/rs3GkdiuYsEERGVHWeEK1FGjh7nb6VjctdGoqMQWbS3+jTGngu38cqvJ7EpYwgcpbwHHxyeVnnBiIjIonBZshIdu5oKRQVa+3I+mKg8XOy1mDc4BFcS72CeYajoOEREZKFYhCtRZGwKZAlozne7E5VbxwYeGNPRF98ZH8cBJVB0HCIiskAcjahER2KTEeBdBa4OOtFRiKzCa70DsHv/AbyqH48tdq/DWcq9/6CibuLBcQkiIgJXhCuNwajg2NVU7h9MZEJOdlr8R/cVrqse+MAwSnQcIiKyMFwRriTnbmUgK8+YPx9cxtvMEtH92sjnMVGzDkuNT+AR+ST6ag6KjkRERBai2BXhJbNfw+h+bTDlmd5FPq+qKv776fsYPzQUU0c/jsvno0we0hpExuRv88QbaRCZ3lTtWjSXLuJN/XOIU6uLjkNERBai2CLctc9gvDv/uwc+f+TAbty8FoMvfv4TE16djS8/ecekAa3F4ZgU1HZ3RJ2qTqKjEFkdnWTEIt1nMELG9LwJMKqS6EhERGQBii3CTZq3hUuVB+9ycGjvDjzW+wlIkoSA4Ba4k5mO5MQEk4a0dKqq4uCVZLThtmlEFcZHTsAs3Xc4pDbGF8Yw0XGIiMgClHtGODnxFjy8ahY8ru7ljeTEW6jm4VXeU1uNmKQsJGbmoo1fNdFRiCzfQ2bsn5AjsFtuhoWGwegon0ZL+VLJzsFdJIiIbFK5i7CqqkV8tehfS25dtwrb1q8CAOiz08t7aYtx+Er+fHA7FmGiCiVJwIe6b3E0rxGm6idhk92bcJWyRcciIiIzVe7t06p71kRiws2Cx0kJt1DNo0aRx/Ya8DTmL1uP+cvWw9PTs7yXthgHrySjmrMdGni6iI5CZPWqSNlYpPsMcaoH3tOPER2HiIjMWLmLcNvO3bB7y29QVRXno47B2cWVYxH3OByTPx8sSXwDD1FlaCVfxBTtWqxVuuA3YyfRcYiIyEwVOxox/70piDp+EOmpKRj3REcMGzcVRoMBANB74Ai06hCKI/t3Y/zQUNg7OGDKW/MqPLQluZWWg6vJWRjVwUd0FCKbMknzO/42NsFb+nFoIsXAX44THYmIiMxMsUV4xvuLH/q8JEl4ccYHJgtkbQ79s39wW84HE1UqraRgqd0S9MmdjfH66Vhv9zZcpBzRsYiIyIzwFssV7PCVZDjbaRBUs4roKEQ2x0tKxRLdEsSo3nhd/zyKfG8vERHZLBbhCnboSjJa+lSFVsNvNZEIHTRn8Yr2F2xUOmCFsafoOEREZEbYzipQalYezsdnoK0vxyKIRBqv2YDu8hF8ZBiJo0pD0XGIiMhMsAhXoMiYFACcDyYSTZZUzNd9CW8pGRPzpiJZdRUdiYiIzEC5b6hBD3YoJhl20KPZigBA0ouOQ2TT3KQ7+EL3KZ7MC8dU/UQs182FRuLQMBGRLeOKcAU6dCUZzaTLcGAJJjILwXIM3teuwF4lBIsMg0THISIiwViEK0hWngFRcWloI58XHYWI7jJMswtDNLux2PgkNhrbiY5DREQCcTSighy7mgqDoqKN7pzoKER0F0kCPtR+i2ilJmbox8NHikfwvQeFu93zOK2y4hERUSXiinAFOXQlGbIEtJIviI5CRPewlwz40m4hqiEDz+fNQEIGb7RBRGSLWIQryOGYZDSuWQVVpGzRUYioCJ5SOr62m49UOGP8yiPINRhFRyIiokrGIlwB8gwKjl5N4bZpRGYuWI7FAt0XOHo1FW+tjYLKW88REdkUFuEKEHUjDTl6hTfSILIAj2sOY1r3Rlhz9Dq+2XtFdBwiIqpEfLNcBTh0JRkA0JpFmMgiTOnaCBfiMzB781k09HJBqOhARERUKbgiXAEOX0lGfU9neLrai45CRCUgyxI+GdIMQTWrYNJPRxGl+IqORERElYBF2MQURcXhmGSORRBZknA3OM2ujm+TRsI97xbG5L2Ga4qn6FRERFTBWIRN7Hx8BtJzDGjDIkxkcWpIqVhhNxd6aDFa/zqSVVfRkYiIqAKxCJvY4Zj8+WDuGEFkmRrKN7DM7hPEqR4Yl/cKslU70ZGIiKiCsAib2P7LSajt7og6VR1FRyGiMmotX8Ai3VIcVxtgsn4SDEZFdCQiIqoALMImZFRU/B11CZ0yNkF63/3+27QSkcXorYnE+9oV2KG0xjvrTnOPYSIiK8Tt00zo9I00pMEFneQo0VGIyARGabfjlloNnx8agJpuDpjSrZHoSEREZEJcETahfZeSAAAd5dOCkxCRqbyq/QVPtqyNBdsvYFkEb7hBRGRNuCJsQvsuJSJQugpPKV10FCIyEUkC5g0KQVauEbM2nIG9VsbI9j73H3jvKFR4WuUEJCKiMuOKsInk6I04HJPMsQgiK6TVyFj8dAt0DfTC279H4dfIa6IjERGRCbAIm8jR2BTkGhQWYSIrZaeV8fmIlujSyAOvrzmJdcfjREciIqJy4miEiURcSoRWltBWPic6ChFVEAedBl8/0xqjvzuEl38+Avs1o9Fbc1h0LCIiKiOuCJvIvkuJaFHPHS5SjugoRFSBHO00+HZMG4RI0Zisn4xdxuaiIxERURmxCJtAWpYeJ+PS0LGBh+goRFQJXOy1WG43FwHSNbyon4YdxpaiIxERURmwCJvA/ugkqCrQuRGLMJGtcJOy8IPdbARK1zBePw1/GNuLjkRERKXEImwC+y4lwtlOg+Z13UVHIaJK5C7dwY92s9FSuogp+kn4xfCY6EhERFQKLMImsO9SItr6VYNOw28nka1xlbKxwm4uusin8LrhBXxr6C06EhERlRCbWzndSM1GdOIddGrIsQgiW+Uo5eG/uvnoJR/CB4ZRWGoYAFVVRcciIqJisAiX075LiQA4H0xk6+wlAz7TLcaT8l58YhiKuVvOswwTEZk57iNcTvsuJcLDxQ4BNVxFRyEiwbSSgk90X8LRkIsv/wISM3Mx+4mmsNNyzYGIyByVqAgfPfAXvln0ARRFQY9+T2HQMy8Ven7nptVY8fkcVPOoAQDoO2gUevQfavq0ZkZVVey7nISODTwgSZLoOERkBmRJxYfab+EVOgELd1zAzbRsfD6iFdwcdaKjERHRPYotwkajEV8teA/vL/we1b288epzA9G2c3fU9WtU6LjOXfvihZffr7Cg5uhiQiZuZ+SiM+eDiegukgRM7d4Idao64o21JzHky7/x7Zg2qFPVSXQ0IiK6S7G/r7t49gRq1vGBd+160Ons0Ll7PxyM2F4Z2cxexMX8+eCODasLTkJE5mhQqzpYMbYtbqbl4InP/8ap62miIxER0V2KXRFOvn0LHl41Cx5X96yJi2eO33fc/r+24PSJQ6hV1w/PTn4bnjVq3XfM1nWrsG39KgCAPju9PLnNwt+XE+Fb3YmrPETWLtytzK/pCGCNUhtjNUvx1Ff7sXR4C3RrXMO0+YiIqEyKXREu8k3P98zDtunUDV//ugeLVmxGs9adsPijV4s8V68BT2P+svWYv2w9PD09yxTYXOiNCg5EJ3PbNCIqlr8ch98mdkRDLxc8/30k/rsnmjtKEBGZgWKLcHUvbyQm3Cx4nHT7Jqp5eBU6popbVejs7AEAPfoPw+Xzp0wc0/ycvJ6KzFwDOh19OX/l598PIqIieLk64JcX26NnkDc+2nQWk1cdw51cg+hYREQ2rdgi3CgwBDevxSD+xjXo9XmI2LEBbTt1L3RMcmJCweeHI3agjk9D0yc1MxEXkyBJQAf5jOgoRGQhnOy0+GJkS7zeOxCbTt3EE5/vw5XEO6JjERHZrGJnhDVaLZ5/ORzvvzwaRkVB975DUK++P376ZiEaBjZF287dsXH1chyK2AmNrzTtjwAAEi9JREFURgOXKu6YMvM/lZFdqN0XEtC0thuqJmaKjkJEluCf3xhJAF4C0FQbjMkZ7yFsSQQWDm2O7kGcGyYiqmzS7xHRQgbVZk0dgsjISBGX/v/27js8qjJ94/j3nCkJgSQkQAoCilJFWCtIEyH86IuislhXkUVxF2UVsaAoKquisqAsIthWrGBDASWIAoIoIi4lIEVADZAJhABpJJOZOb8/QgkYzEDKCZn7c11zZeacM+FOnkzycOY971tmngP5XPrUl4zq2Zx/LLnI7jgicpraMSKNYW+tImVnFnclNeWfSU0xTc1JLiJS3pq0aM2EVz/93XYtd3QKFmzwANCzVYLNSUTkdNYgJoIPhnXgmosa8MKXW/jra9+TnpVvdywRkZChRvgUzE/x0CSuFk3iatkdRUROZ2OjCf9XLM+mXMbTzums+nUfvSZ9TfJ6j93JRERCghrhk5SZ62XF9kx6ttJ4PhEpH4YB1zoXM/euTpwRU4Pb31zFgx+tI8973KwSxWeo0Uw1IiJlpkb4JC38KR1/wKJXq8TSDxYROQnn1KvFR3d05PYuZ/Peyt/oN3kZKTu1Gp2ISEVRI3ySFqz3cEbtGpx3RpTdUUSkuhkbjXtcDA+uuJS3nePILfAx4MVvmLLoZwr9AbvTiYhUO2qET0JOgY+vt2TQs1UChqEru0Wk4nRwbGD+iMvo3jKeZ5M3ceWUb0gJnGV3LBGRakWN8ElYvGk3Xl+AXitv0fg8EalwMTXdTL3xIl668UJ2ZxdwhfcJniq8lnzLZXc0EZFqQY3wSZif4qEuB7jI2Gx3FBEJIb3OS2Th3V0Y6FjCNH9/ennH862/pd2xREROe6WuLCdF8gv9LNq4m/6OH3AYtqxBIiKhpti7TtHA0y7oby7nAd9Qriscw6DAIkY92pC6Rlax5+jiOhGRYKkRDtI3P2eQ6/XT07XS7igiEsI6ODaQbN7PJN9VvOrvw2f+doxwfsRfHcm4Df/vn3D8EC41yiIiR2hoRJCS13uIDHPSwVxvdxQRCXE1DC8Put4j2X0/l5gbGee7kV7e8XzlPx/L0jtWIiLB0hnhIPj8Ab7YkE5SyzjcG0s44yIiYoNzzDRecz/HIv+feMJ3E7cW3keX11cypl9LmsRFlvykki7y1VliEQlROiMchO9/yWRfXiG9zkuwO4qIyO90dawh2X0/Y5wz+PG3ffSY+DWj3l9Damae3dFERKo0NcJBSE7xEO4yuaxZPbujiIiUyGX4GeKcz2JrCIPNz/hk1Ta6PfMFDxcOxmPF2B1PRKRK0tCIUgQCFsnr0+nSrB4Rbn27RKRqq2NkM8b1FkOd8/iP70pm+rsyy9+FmxwLucP56bEzTBymC+pEJETpjHAp1uzYjycrX8MiROS0kmDsY5zrdb5yj+QKx3Je9/eic8EkHi28mdSA3t0SEQE1wqVKXp+O0zTo1jze7igiIietobmHZ13T+cI9ij7mCt7xJ9HFO5Hh3jtZF2hsdzwREVupEf4Dhf4Ac9bsov05dYiO0JKmInL6OsdMY4J7GkvDRjDUMY/FgT/xZ++/uN47msX+NgQsw+6IIiKVTo3wH5i3No2d+w9yc/uz7I4iIlIuEox9POh6l+Vhd/Kg8x22BupzS+EDJHmfY7qvL5nWCaZdExGphnT11wlYlsVLS7bSNK4W3VrE2R1HRKRcRRkHud05l8GOz5kXuJR3fEk86buB53x/oc97/+P6dmdyyVkxGMYpnikuab7iY/brgjwRsZ8a4RNYsnkPGz3ZPHtNG0xTbxmKSPXkNvwMcHzDAMc3bAo04B1/Eh/91JfZq3fRNK4WAy9uQP8vk0gw9h19kppYEakmNDTiBF5aspWEqHCuOP8Mu6OIiFSK5uYOHnO9wYqHknjm6jbUDHPy5GcbaV8wmeu9o5nl60KWVcPumCIi5UaNcAlWp+7nu22ZDOnUGLdT3yIRCS0Rbid/uaQhs//RkUX3Xs5djo/ZZdXlPt/tXFwwlX+8/SOfrUsjt8Bnd1QRkTLR0IgSTFuylchwJ9d92R6+yrc7joiIbRrXrcndrg/5p/Uha6xzmO3vyJxt/Zm3Lg2306RTk7r0ODeepJbx1IsMszuuiMhJUSNc3NhotgcSmO99jjscc6jlUhMsIgJgGHC+sZXzza08PPp5Vv26jwUb0kle7+GrjbsxjHVc1CiGpJbxdG5al3MtA9Ow7I4tIvKH1AgfZ7q/Ly583OKcb3cUEZEqyekwaXd2HdqdXYeH+7ZkoyebBevTWbDBw/j5Gxk/H+rwIh3NFDqb6+jsWHfsxXYiIlWEGuFidlvRfOjvzNWOpcQZuipaRKQ0hmHQMjGKlolRjOjelN3Z+XzzcwZL3/8PXwda82mgI/jgHGMnbc1NXGxuoq2xkQaWdepTs4mIlBM1wsX819eLQpzc5phndxQREfuUNgfwH4iLDGfABQ0Y8MlULAs2Wg1ZGmjDd4GWzPO3411/NwDin/qSi8+K5ZIzY2jTsDbnJkYR7nKU11cgIhIUNcKH5BT4eNPfnV7mShqbHrvjiIic9gwDWhqptDRTuY15BCyDzVYDVgaas7Lxw6z8JZN5a9MAcJgGzeIjaXNGNG0aRtP6jGiaxUeqORaRCqVG+JB3V/xGNjUZ5pxjdxQRkWrJNCxaGKm0MFO56boPAUg7cJC1Ow6wbscB1u48wIINHmb+kHroeDizTk2ax0fSPOHorVFsBC6HprYUkbJTIwxs2JXFtK+3cqm5nj+Z2+yOIyJS/R0afpF46Nbz0Gp1lmWxY99BUnYeYKMnm02ebDanZ7Ngg4fAoUkonKZBo9gIGmcupbGRdujm4cyRX5EQFY5Dq4GKSJBCvhH+fF0a98xaQ1QNJ2OdM+yOIyJS9ZVhDHFpDMOgYWwEDWMj6N068cj2/EI/W9Jz2JSezfaMHLZn5LJtbx2+CbQin0PzFz/9FU7TILF2OA1qR9AwtgYNVk8kkb0kGJkkGpnEj15LZLjrxF+Llo8WCSlBNcI/freEV55/nEAgwP/1+wtX33THMfsLvQVMGncvWzelEBlVm3sfn0x8YoMKCVxeAgGLSWP+xgv+q7jA2MK0gonEmfvtjiUiIocVa1LDgdaHbkeEQcAy8BDDtkB9Uq167LDqkbq/Hjv21WPxtnrsZuBxn3MBNd0O4qPDiY8Mp473TuoaB6hnHKAuB6j7UzqxNd3ERLiJqekmKtyp2S1EqrFSG2G/38+0fz/KYxNnUCcugVF/u5K2nbrTsHHTI8d8MXcWtSKjeGnmIpYunMOMqeMZ9fjkCg1eFjkFPu6ZuZoF/qsY6FjMOOdrhBlaKlRExDaneJbZNCzqk0l9R2aJ+/MtFx4rFg8xpFuxeKxY0vyxpO+NYXdGbVJoTEYgihwiip7wxg/HPN9hGtSu4aJ2hIvaEUWNcdSW2UQaeUSRR5SRS+Sfn6RWmJOabicRYY6i+4ce13A7iHA7NKZZpIoqtRHe8tMaEhucScIZjQDo1L0fK5Z9cUwj/P2yhVx76wgAOlzem+kTx2JV0Tkif92by9AZP7B1Ty6POGcw2DGfKhhTRETKQbhRyFlGOmeR/ofH5Vsu9ljR7B26ir05BezLK2R/npd9ed4j9/fnFZKR42WbdQ7ZgQiyiMCHEz5OKTWH0zSo4XZQw+WghttBuNNBuMskzOkg7NDHcJeJ22kS5jRxOUzcjqLHriMfDZymictp4jINnI6j2xymgdM0cBZ7XHQD0yjaZppFjb3DMDCMov2mUbTfLHbfOLzNKNpmYIBRdPGiYRgYHD0OimYHMTAOfTx6zNH9+iMrVVepjXDmHg91446O06pTL5EtG1Yfd0z6kWMcTicRNSPJPrCPqNqx5Ry3bLakZzNw2rdYFrwxuC2d3r7O7kgiIlIFhBuFNDQyaPjqmb/fefy44bF9AbAsyMdN1shUcgp85Bb4yCnwkfffgeQSTo4VzkHCim6dHiTP6ye/0E+e10+Bz09+YYACn5/sfB8ZO36gADcFlgsvTgoj4vD6Anh9AXyB6rNU9eFmuei+Uex+sWM45kFJd094AuuY5wZx/KlQW182C0d2ITG6ht0xjii1EbZKev0d9xNllXRQCT91yZ+8y4JP3wVg56/baNKi9e+OqQxZ+zO55d1YoJEt/75Unqz9mVXuP2RSMVTr0FHptX7v+L9Vx/3tmHlJ6Z/jg4GlH3MCoXBVu3WC+3pdVz+dZ5W8vaJrvduzs8Ttpb6+6sQlkLE77cjjvXvSiK0bV+IxdeMS8ft85OVmExlV+3efq+cV19HzCvvPwo4c0p8Jr35qdwypBKp16FCtQ4dqHTpU69BhV61LHb3ftEUb0lJ/IX1XKoWFXpYtnEvbjt2POaZtxyQWfV40OfryxZ/T+sL2GhMkIiIiIlVaqY2ww+lk6D1jeeyemxl+Qw86dutLo7Ob8c4rE/l+2UIAuvcbRPaB/Qwb1JVPZ77KX4fdV+HBRURERETKIqihRxe378rF7bses+36v9195L47LIz7xk0p32QVqEd/+4dnSOVQrUOHah06VOvQoVqHDrtqbcxetq36XI4qIiIiIhIkzfAtIiIiIiGp2jbCP363hL9fl8SwQV358M2pv9tf6C3g2UfuZNigrowaOoD0tB02pJTyUFqtP3nvFYbf2IMRN/dmzIgbTjiFilR9pdX6sOWLPuPKTmfz88a1lZhOylMwtV725TyG39iDO2/syYSxIyo5oZSX0mq9x7OTh++8nrsH92PEzb354dtFNqSU8jD5yfu4ud8l3HVTrxL3W5bFy5MeY9igroy4uTdbN5W+WE1ZVctG+PCy0I889zqT30pm6cI5pG7fcswxxZeF7j/oVmZMHW9TWimLYGp9drNWTHjlE55/43M6XN6bN1582qa0UhbB1BrgYF4Ocz94g2bnnm9DSikPwdR6V+p2PnxrKk+/+D6T30pmyIgxNqWVsgim1rPemELHbn2Y+Ppc7h37AtMmPGJTWimrbn2u4ZEJr59w/6rvFpOW+gtT3/uKv496kpeeq/jXdbVshIsvC+1yuY8sC13c98sW0rX31UDRstBrVy0veWEQqdKCqXXrC9sTFl60ik3zVhewd4/HjqhSRsHUGuDtl//NgOtvw+UOsyGllIdgar1gzkz6XHUTtaKiAagdU9eOqFJGwdTaMAwO5uYAkJubTWzdeDuiSjlodX5bapWwzsRh3y9dyOW9BmAYBs3Pu4DcnCwyM3ZXaKZq2QiXtCx05p70444peVloOb0EU+viFs6dxYXtulRGNClnwdR62+b1ZOxO45KOSZUdT8pRMLXelbqdnanbeeCOgdx321X8+N2Syo4p5SCYWl976wgWL5jNkAEdeOLeWxn6z0crO6ZUksyM434e4hLIzKjYk1fVshEuz2WhpWoLptaHLU6ezc8b1zHg+qEVG0oqRGm1DgQCvPrCOAYPf6jyQkmFCOZ1HfD7SEv9hXGT32Hk2OeZMv5BcrKzKieglJtgar104ad0630Nr368nDHPvcakcSMJBAKVE1AqVcnvzFdsb1YtG+GTWRYa+MNloaVqC6bWAGtWLuODGVMYPX663jI/TZVW64N5Ofy2fTMP33kdQ6/pzOYN/+Nf99+mC+ZOQ0H9Dq+XQNvO3XE6XcTXb0j9Ro1J27G9sqNKGQVT64Vz36djtz4AtDjvQgoLCsg6kFmpOaVy1KmXeOzPw25PhQ+FqZaNsJaFDh3B1Hrb5vW8+OzDjH56usYRnsZKq3XNWlG8OW8VL3+wlJc/WEqzcy/gofHTadKijY2p5VQE87pu17kHKT9+B0DW/kx2pf5CfP1GdsSVMgim1vXi67N21XIAUn/5Ga+3gOjadeyIKxWsbackFs//GMuy2JTyP2rWiizx5FZ5CmpludNN8WWh/YEA3fsOPLIsdJMWrWnbqTvd+w1i0hP3MGxQVyKjohk59gW7Y8spCKbW/53yFPkHc3lmzHCg6JfqQ+Nftjm5nKxgai3VQzC1vqDdZaxeuZThN/bANE1u+fsDREXH2B1dTlIwtR48fDRTnhnNnJmvgWFw10PP6sTVaWrCo3eRsnoFWfv3MWRAB64dMgK/zwdArytv4KL2XVn17WKGDepKWHg4d41+psIzaWU5EREREQlJ1XJohIiIiIhIadQIi4iIiEhIUiMsIiIiIiFJjbCIiIiIhCQ1wiIiIiISktQIi4iIiEhIUiMsIiIiIiFJjbCIiIiIhKT/B63ipYJ2n5W6AAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 864x432 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"\"\"\"\n",
"Created on Tue Jun 20 07:35:08 2023\n",
"\n",
"@author: Irene Markelic\n",
"\"\"\"\n",
"\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"from scipy.stats import uniform, beta, gamma\n",
"import typing\n",
"\n",
"class Beta_Sim:\n",
" \"\"\"Create a random variable with a Beta distribution with a,b \n",
" positive integers, using Universality of the Uniform. \n",
"\n",
" This is an implementation of example 12.1.4 on p. 539 of the book\n",
" \"Introduction to Probability\" by Blitzstein and Hwang, 2nd edition\n",
"\n",
" The Beta distribution is a continuos function over the support [0,1] and\n",
" is a generalization of the Uniform distribution. A random variable(r.v.)\n",
" X is said to have the Beta(a,b) distribution if its probability density \n",
" function (pdf) is: f(x) = x^(a-1)*(1-x)^(b-1)*1/b(a,b),\n",
" where b(a,b) is a constant to make the integral of the pdf equal 1.\n",
" For a,b positive integers b(a,b) = (a-1)!(b-1)!/(a+b-1)!\n",
"\n",
" A Beta distribution can be described by two independent r.v.s X,Y with \n",
" the same rate lambda, that have a Gamma distribution (cf. p. 396, story \n",
" 8.5.1 \"Bank-post office\"), thus:\n",
" X / (X+Y) ~ Beta, with X~Gamma(a,1), Y~Gamma(b,1)\n",
" Note, here lambda=1.\n",
" We also know that the sum of independent Gammas with the same rate\n",
" lambda has a Gamma distribution, thus X+Y is Gamma(a+b,lambda).\n",
"\n",
" Page 539 gives the solution to the task of how to simulate an r.v.\n",
" with a Beta distribution:\n",
" \"Applying universality of the Uniform directly for the Beta is hard, so \n",
" let's first use the bank-post office story: if X~Gamma(a,1) and \n",
" Y~Gamma(b,1) are independent, then X/(X+Y)~Beta(a,b).\n",
" So if we can simulate the Gamma distribution, then we can simulate the\n",
" Beta distribution!\n",
"\n",
" To simulate X~Gamma(a,1), wen can use X1+X2+..Xa with the X_j i.i.d. \n",
" Expo(1) r.v.s. [...] Lastly, by taking the inverse of the Expo(1) CDF\n",
" and applying universality of the Uniform, -log(1-U)~Expo(1) for\n",
" U~Unif(0,1), so we can easily construct as many Expo(1) r.v.s as we \n",
" want.\"\n",
"\n",
" In other words, if we want a Beta r.v. we need 2 independent Gammas with\n",
" the same rate lambda (here we simply choose 1). And to get the Gammas we\n",
" use the fact that the sum of exponentially distributed r.v.s with the \n",
" same rate lambda is a Gamma r.v. We can create r.v.s that have an \n",
" exponential distribution (Expo(lambda), here lambda = 1), by using \n",
" universality of the Uniform. The pdf for a r.v. with an exponential \n",
" distribution is Expo(lambda) = lambda e^(-lambda*x).\n",
" Universality of the Uniform says that we can create a r.v. with any \n",
" desired distribution if we know the inverse of its Cumulative \n",
" Distribution Function (CDF), F^(-1) and plug in a Uniform r.v.\n",
" The CDF of Expo(lambda) = F(x) = 1-e^(-lambda x) and the inverse\n",
" F^(-1) = -log(1-U). U is a r.v. with the uniform distribution.\n",
"\n",
" So: First we draw many samples from the uniform, \n",
" then we plug in each of these values into F^(-1) to get our Expo(1) \n",
" r.v.s. For X we do this a times and for Y b times. X/(X+Y) is then\n",
" our beta r.v.\n",
"\n",
" \"\"\"\n",
"\n",
" def __init__(self, a: int, b: int, number_of_simulations: int) -> None:\n",
" \"\"\"\n",
"\n",
" Parameters\n",
" ----------\n",
" a : int\n",
" first parameter for Beta(a,b) distribution.\n",
" b : int\n",
" second parameter for Beta(a,b) distribution.\n",
" number_of_simulations : int\n",
" how often you want to sample.\n",
"\n",
" Returns\n",
" -------\n",
" None.\n",
"\n",
" \"\"\"\n",
" self.a = a\n",
" self.b = b\n",
" self.number_of_simulations = number_of_simulations\n",
"\n",
" def do_Beta_Universality_Uniform(self)->np.ndarray:\n",
" \"\"\"\n",
" Creates and array with samples from the desired beta(a,b) distribution.\n",
"\n",
" Returns\n",
" -------\n",
" beta_rv : numpy.ndarray\n",
" an array with samples from the desired beta(a,b) distribution.\n",
"\n",
" \"\"\"\n",
" # 1) Create \"a\" and \"b\" many Expos by drawing n samples\n",
" # from the Uniform distr. and plugging them into F^(1) of the\n",
" #Expo (-log(1-U))\n",
" n = self.number_of_simulations\n",
" X = []\n",
" for i in np.arange(self.a):\n",
" uniform_samples = uniform().rvs(size=n)\n",
" X.append((-1*np.log(1-uniform_samples)))\n",
" Y = []\n",
" for i in np.arange(self.b):\n",
" uniform_samples = uniform(0, 1).rvs(size=n)\n",
" Y.append(-1*np.log(1-uniform_samples))\n",
"\n",
" # 2)Create Gammas from Expos by summing all the uniform samples\n",
" X = np.array(X)\n",
" X = X.sum(axis=0)\n",
"\n",
" Y = np.array(Y)\n",
" Y = Y.sum(axis=0)\n",
"\n",
" def plot_gammas():\n",
" # plot them just to see that they are correct\n",
" fig, axs = plt.subplots(nrows=1, ncols=2, sharex=True)\n",
" axs[0].hist(X, density=True, bins=300,)\n",
" axs[0].set_title(f'gamma({self.a},1)')\n",
"\n",
" axs[1].hist(Y, density=True, bins=300)\n",
" axs[1].set_title(f'gamma({self.b},1)')\n",
"\n",
" # compare the plots to scipy.gamma\n",
" x_g = np.arange(0, 30, 0.01)\n",
" y_g1 = gamma(self.a).pdf(x_g)\n",
" axs[0].plot(x_g, y_g1)\n",
" axs[0].legend(['scipy gamma pdf', 'using uniforms'])\n",
"\n",
" y_g2 = gamma(self.b).pdf(x_g)\n",
" axs[1].plot(x_g, y_g2)\n",
" axs[1].legend(['scipy gamma pdf', 'using uniforms'])\n",
"\n",
" plot_gammas()\n",
"\n",
" # 3) The final Beta r.v. is X/(X+Y)\n",
" C = X+Y\n",
" beta_rv = X/C\n",
"\n",
" return beta_rv\n",
"\n",
" def do_Beta_Scipy(self) -> typing.Tuple[int, int]:\n",
" \"\"\"\n",
" compute the beta r.v. using scipy.stats.beta for comparing results\n",
"\n",
" Returns\n",
" -------\n",
" x : list[int]\n",
" x vals for plotting.\n",
" y : list[int]\n",
" y vals for plotting of beta.\n",
"\n",
" \"\"\"\n",
" x = np.arange(0, 1, 0.01)\n",
" y = beta(self.a, self.b).pdf(x)\n",
"\n",
" return x, y\n",
"\n",
"\n",
"def main():\n",
"\n",
" a = 3\n",
" b = 10\n",
" number_of_simulations = 10000\n",
" beta_custom = Beta_Sim(a, b, number_of_simulations)\n",
"\n",
" beta_scipy_x, beta_scipy_y = beta_custom.do_Beta_Scipy()\n",
" beta_uni_rv = beta_custom.do_Beta_Universality_Uniform()\n",
" #beta_mh_x, beta_mh_y = beta.do_Beta_MH()\n",
"\n",
" # plot the beta distributions on top of each other. They should look alike\n",
" fig, ax0 = plt.subplots(nrows=1, ncols=1, sharex=True,\n",
" figsize=(12, 6))\n",
" # plot beta_scipy\n",
" ax0.plot(beta_scipy_x, beta_scipy_y)\n",
" fig.set_facecolor('lightsteelblue')\n",
"\n",
" # plot beta_uni\n",
" ax0.hist(beta_uni_rv, density=True, bins=100)\n",
" ax0.legend(['using scipy', 'using uniform'])\n",
" ax0.set_title(f'beta({beta_custom.a},{beta_custom.b})')\n",
"\n",
"\n",
"if __name__ == \"__main__\":\n",
" main()\n",
"else:\n",
" print('__name__={}'.format(__name__))\n"
]
},
{
"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.6"
}
},
"nbformat": 4,
"nbformat_minor": 4
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment