From 1ac70f1c869f4b3ca15fc1111ba306d169948de4 Mon Sep 17 00:00:00 2001 From: mbi6245 Date: Tue, 13 Aug 2024 13:01:03 -0700 Subject: [PATCH] probably fully functional MVP ensemble fitter with L2 norm; analytical solutions to all distributions except fisk, weibull still broken --- simulations.ipynb | 136 ++++++++++++++++++++++++++++++++- src/ensemble/distributions.py | 41 ++++------ src/ensemble/ensemble_model.py | 42 +++++++--- tests/test_distributions.py | 5 +- 4 files changed, 182 insertions(+), 42 deletions(-) diff --git a/simulations.ipynb b/simulations.ipynb index 2c38e9a..349b892 100644 --- a/simulations.ipynb +++ b/simulations.ipynb @@ -4,12 +4,64 @@ "cell_type": "code", "execution_count": 1, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[-1.88951825 -1.67562088 -1.65214474 -1.62064284 -1.562736 -1.50039701\n", + " -1.46558166 -1.41107123 -1.33930541 -1.33348901 -1.22021545 -1.16896446\n", + " -1.11823493 -1.11364725 -1.0355334 -0.99155624 -0.9322802 -0.91862877\n", + " -0.89289741 -0.87632799 -0.85337447 -0.82952583 -0.80436763 -0.77718609\n", + " -0.77606566 -0.74560126 -0.6891231 -0.66738905 -0.65431026 -0.60114386\n", + " -0.51794973 -0.51688887 -0.4953134 -0.44496799 -0.40446697 -0.39829545\n", + " -0.39799426 -0.31320473 -0.284331 -0.28041655 -0.22160632 -0.17987713\n", + " -0.17326124 -0.11955245 -0.08257138 -0.00575042 0.01167684 0.03075504\n", + " 0.05233829 0.06073081 0.08723611 0.09763968 0.10284089 0.13299994\n", + " 0.1459378 0.25129287 0.2637708 0.26529311 0.28513058 0.3164859\n", + " 0.35035412 0.35802117 0.36623514 0.3697894 0.40023227 0.41085433\n", + " 0.48258617 0.48964917 0.50151693 0.50449019 0.51650857 0.53792918\n", + " 0.56122631 0.57075674 0.57988314 0.5942923 0.62610387 0.6539129\n", + " 0.66528933 0.73613501 0.74289003 0.80039489 0.80555165 0.81712376\n", + " 0.94335029 0.99278054 1.01699842 1.02001017 1.05444052 1.07084284\n", + " 1.09992702 1.16661601 1.17901527 1.20899573 1.2415302 1.3356225\n", + " 1.4081054 1.61840104 1.83143626 1.85685823]\n", + "\n" + ] + } + ], "source": [ "import numpy as np\n", "import scipy.stats\n", + "import matplotlib.pyplot as plt\n", "from ensemble.distributions import distribution_dict\n", - "from ensemble.ensemble_model import EnsembleFitter" + "from ensemble.ensemble_model import EnsembleFitter\n", + "\n", + "\n", + "std_norm = distribution_dict[\"normal\"](0, 1).rvs(100)\n", + "model = EnsembleFitter([\"normal\", \"gumbel\"], None)\n", + "res = model.fit(std_norm)\n", + "print(res.weights)\n", + "\n", + "\n", + "\n", + "# mean = 1\n", + "# var = 2\n", + "# x = np.linspace(0, 5, 100)\n", + "# y = np.linspace(-4, 5, 100)\n", + "# std_norm = distribution_dict[\"normal\"](0, 1)\n", + "# q = std_norm.ppf(0.05)\n", + "# p = std_norm.cdf(q)\n", + "# print(q, p)\n", + "\n", + "# a = 0.3\n", + "# b = 0.7\n", + "# exp = distribution_dict[\"exponential\"](mean, var)\n", + "# fisk = distribution_dict[\"fisk\"](mean, var)\n", + "# # print(a * exp.(x) + b * norm.pdf(x))\n", + "\n", + "# temp = EnsembleFitter([\"exponential\", \"fisk\"], None)\n", + "# temp.fit(x)\n" ] }, { @@ -18,6 +70,86 @@ "source": [ "working on finding some way to take linear combination of all PDFS and then take draws from that distribution, cant really do that with scipy.stats.rv_continuous.rvs unfortunately" ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# How does Scipy work?" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "# helper funcs\n", + "def reverse_z(z, mean, variance):\n", + " return z * np.sqrt(variance) + mean" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "(4.0, 1.0)\n", + "0.022750131948179195\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAA94AAAH9CAYAAADsyqSaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABGm0lEQVR4nO3dfVRVdd7//9cR5KCMnBSSI5eA1OWkhhmCKVBpl4mZmtU0Uk5Yk5Yub5LoTsaaQedK0iklb9Cc1ZLMS+WaVZhNNokz5c1ojaI0aV3dzGCYcYbRb4FYgcH+/cHPUydu9CCfc0Cfj7X2WpzP/uzNe2/xs3mx72yWZVkCAAAAAABGdPJ3AQAAAAAAXMgI3gAAAAAAGETwBgAAAADAIII3AAAAAAAGEbwBAAAAADCI4A0AAAAAgEEEbwAAAAAADCJ4AwAAAABgEMEbAAAAAACDCN4AAAAAABhE8AYAAC3auXOnxo8fr8jISNlsNm3evLlRnw8//FC33HKLHA6HunXrpmHDhqmsrMw9v6amRrNnz1Z4eLhCQkJ0yy236PPPP/fhVgAA4D+B/i6grdTX1+uLL75Qt27dZLPZ/F0OAOAiZ1mWTp48qcjISHXq1LH/zn3q1CkNGjRIv/zlL/Wzn/2s0fx//OMfuvbaazVlyhTNnz9fDodDH374oYKDg919MjIy9Nprr2nTpk0KCwvTww8/rHHjxqm4uFgBAQHnVAfHegBAe3Oux3ubZVmWD+sy5vPPP1dUVJS/ywAAwMPRo0fVu3dvf5fRZmw2mwoLC3Xrrbe62+6880517txZL730UpPLVFZW6tJLL9VLL72ktLQ0SdIXX3yhqKgobd26VaNHjz6n782xHgDQXp3teH/BnPHu1q2bpIYNDg0N9XM1aJdOnZIiIxu+/uILKSTEv/UAuKBVVVUpKirKfXy6UNXX1+v111/XY489ptGjR+vgwYOKjY1VVlaWO5wXFxfr9OnTSk1NdS8XGRmpuLg47dmzp9ngXVNTo5qaGvfnM+cKPI71jO0AAD861+P9BRO8z1xyFhoaSvBG0354KWNoKL+cAfCJC/2S6IqKClVXV+vpp5/Wf//3f2vRokX605/+pNtvv11vvfWWhg8fLpfLpaCgIHXv3t1j2YiICLlcrmbXnZOTo/nz5zdq9zjWM7YDANqBsx3vO/ZNZwAAwK/q6+slSRMmTNBDDz2kq6++WnPnztW4ceO0evXqFpe1LKvFX1SysrJUWVnpno4ePdqmtQMA4CsEbwAA0Grh4eEKDAzUgAEDPNr79+/vfqq50+lUbW2tvvzyS48+FRUVioiIaHbddrvdfXabK9oAAB0ZwRsAALRaUFCQhgwZoo8++sij/eOPP1ZMTIwkKSEhQZ07d1ZRUZF7fnl5uQ4dOqTk5GSf1gsAgD9cMPd4AwAAM6qrq/Xpp5+6P5eWlqqkpEQ9evRQdHS0Hn30UaWlpen666/XDTfcoD/96U967bXX9Pbbb0uSHA6HpkyZoocfflhhYWHq0aOHHnnkEQ0cOFA33nijn7YKAADfIXgDAIAW7d+/XzfccIP7c2ZmpiTpnnvuUX5+vm677TatXr1aOTk5evDBB3XFFVfo5Zdf1rXXXuteZunSpQoMDNTEiRP1zTffaOTIkcrPzz/nd3gDANCRXTDv8a6qqpLD4VBlZSX3gKFpp05JP/lJw9fV1Tz5FoBRHJfaXpP7lLEdAOBH53q85x5vAAAAAAAMIngDAAAAAGAQwRsAAAAAAIMI3gAAAAAAGNSq4J2Xl6fY2FgFBwcrISFBu3btOqfl/vrXvyowMFBXX311o3kvv/yyBgwYILvdrgEDBqiwsLA1pQEAAAAA0K54HbwLCgqUkZGhefPm6eDBg7ruuus0ZswYlZWVtbhcZWWlJk+erJEjRzaat3fvXqWlpSk9PV3vvfee0tPTNXHiRL377rvelgcAAAAAQLvi9evEhg4dqsGDB2vVqlXutv79++vWW29VTk5Os8vdeeed6tu3rwICArR582aVlJS456WlpamqqkpvvPGGu+2mm25S9+7dtXHjxnOqi9e24Kx45QwAH+K41PZ4nRgAoL0x8jqx2tpaFRcXKzU11aM9NTVVe/bsaXa5tWvX6h//+Id+85vfNDl/7969jdY5evToFtdZU1OjqqoqjwkAAAAAgPbGq+B9/Phx1dXVKSIiwqM9IiJCLperyWU++eQTzZ07V//zP/+jwMDAJvu4XC6v1ilJOTk5cjgc7ikqKsqbTQEAAAAAwCda9XA1m83m8dmyrEZtklRXV6dJkyZp/vz5+ulPf9om6zwjKytLlZWV7uno0aNebAEAAAAAAL7R9CnoZoSHhysgIKDRmeiKiopGZ6wl6eTJk9q/f78OHjyoWbNmSZLq6+tlWZYCAwO1bds2/dd//ZecTuc5r/MMu90uu93uTfkAAAAAAPicV2e8g4KClJCQoKKiIo/2oqIiJScnN+ofGhqq999/XyUlJe5p+vTpuuKKK1RSUqKhQ4dKkpKSkhqtc9u2bU2uEwAAAACAjsSrM96SlJmZqfT0dCUmJiopKUlr1qxRWVmZpk+fLqnhEvBjx45p3bp16tSpk+Li4jyW79mzp4KDgz3a58yZo+uvv16LFi3ShAkT9Oqrr2r79u3avXv3eW4eAAAAAAD+5XXwTktL04kTJ7RgwQKVl5crLi5OW7duVUxMjCSpvLz8rO/0/rHk5GRt2rRJTzzxhJ588kldfvnlKigocJ8RB9B6fea+7vUyR54ea6ASAABgirfHe471gG95/R7v9or3peKsLtJ3vRK8Af/guNT2eI830DyCN+AfRt7jDQAAAAAAvEPwBgAAAADAIII3AAAAAAAGEbwBAAAAADCI4A0AAAAAgEEEbwAAAAAADCJ4AwAAAABgEMEbAAAAAACDCN4AAAAAABhE8AYAAAAAwCCCNwAAAAAABhG8AQAAAAAwiOANAAAAAIBBBG8AAAAAAAwieAMAAAAAYBDBGwAAAAAAgwjeAAAAAAAYRPAGAAAAAMAggjcAAAAAAAYRvAEAAAAAMIjgDQAAAACAQQRvAAAAAAAMIngDAAAAAGAQwRsAAAAAAIMI3gAAAAAAGETwBgAAAADAIII3AAAAAAAGBfq7AADtT5+5r3u9zJGnxxqoBAAAmMCxHvAtzngDAAAAAGAQwRsAAAAAAIO41BzoQFpzWRgAAOhYON4DFx7OeAMAAAAAYBDBGwAAtGjnzp0aP368IiMjZbPZtHnz5mb7Tps2TTabTbm5uR7tNTU1mj17tsLDwxUSEqJbbrlFn3/+udnCAQBoJwjeAACgRadOndKgQYO0YsWKFvtt3rxZ7777riIjIxvNy8jIUGFhoTZt2qTdu3erurpa48aNU11dnamyAQBoN7jHGwAAtGjMmDEaM2ZMi32OHTumWbNm6c0339TYsZ6vHKqsrNQLL7ygl156STfeeKMkaf369YqKitL27ds1evRoY7UDANAecMYbAACcl/r6eqWnp+vRRx/VlVde2Wh+cXGxTp8+rdTUVHdbZGSk4uLitGfPnmbXW1NTo6qqKo8JAICOiOANAADOy6JFixQYGKgHH3ywyfkul0tBQUHq3r27R3tERIRcLlez683JyZHD4XBPUVFRbVo3AAC+QvAGAACtVlxcrOeee075+fmy2WxeLWtZVovLZGVlqbKy0j0dPXr0fMsFAMAvCN4AAKDVdu3apYqKCkVHRyswMFCBgYH67LPP9PDDD6tPnz6SJKfTqdraWn355Zcey1ZUVCgiIqLZddvtdoWGhnpMAAB0RARvAADQaunp6fr73/+ukpIS9xQZGalHH31Ub775piQpISFBnTt3VlFRkXu58vJyHTp0SMnJyf4qHQAAn+Gp5gAAoEXV1dX69NNP3Z9LS0tVUlKiHj16KDo6WmFhYR79O3fuLKfTqSuuuEKS5HA4NGXKFD388MMKCwtTjx499Mgjj2jgwIHup5wDAHAhI3gDAIAW7d+/XzfccIP7c2ZmpiTpnnvuUX5+/jmtY+nSpQoMDNTEiRP1zTffaOTIkcrPz1dAQICJkgEAaFcI3gAAoEUjRoyQZVnn3P/IkSON2oKDg7V8+XItX768DSsDAKBj4B5vAAAAAAAMIngDAAAAAGBQq4J3Xl6eYmNjFRwcrISEBO3atavZvrt371ZKSorCwsLUpUsX9evXT0uXLvXoc+bdnz+evv3229aUBwAAAABAu+H1Pd4FBQXKyMhQXl6eUlJS9Pzzz2vMmDH64IMPFB0d3ah/SEiIZs2apauuukohISHavXu3pk2bppCQED3wwAPufqGhofroo488lg0ODm7FJgEAAAAA0H54HbyXLFmiKVOmaOrUqZKk3Nxcvfnmm1q1apVycnIa9Y+Pj1d8fLz7c58+ffTKK69o165dHsHbZrPJ6XS2ZhsAAAAAAGi3vLrUvLa2VsXFxUpNTfVoT01N1Z49e85pHQcPHtSePXs0fPhwj/bq6mrFxMSod+/eGjdunA4ePNjiempqalRVVeUxAQAAAADQ3ngVvI8fP666ujpFRER4tEdERMjlcrW4bO/evWW325WYmKiZM2e6z5hLUr9+/ZSfn68tW7Zo48aNCg4OVkpKij755JNm15eTkyOHw+GeoqKivNkUAAAAAAB8olXv8bbZbB6fLctq1PZju3btUnV1td555x3NnTtX//mf/6m77rpLkjRs2DANGzbM3TclJUWDBw/W8uXLtWzZsibXl5WVpczMTPfnqqoqwjcAAAAAoN3xKniHh4crICCg0dntioqKRmfBfyw2NlaSNHDgQP3rX/9Sdna2O3j/WKdOnTRkyJAWz3jb7XbZ7XZvygcAAAAAwOe8utQ8KChICQkJKioq8mgvKipScnLyOa/HsizV1NS0OL+kpES9evXypjwAAAAAANodry81z8zMVHp6uhITE5WUlKQ1a9aorKxM06dPl9RwCfixY8e0bt06SdLKlSsVHR2tfv36SWp4r/czzzyj2bNnu9c5f/58DRs2TH379lVVVZWWLVumkpISrVy5si22EQAAAAAAv/E6eKelpenEiRNasGCBysvLFRcXp61btyomJkaSVF5errKyMnf/+vp6ZWVlqbS0VIGBgbr88sv19NNPa9q0ae4+X331lR544AG5XC45HA7Fx8dr586duuaaa9pgEwEAAAAA8B+bZVmWv4toC1VVVXI4HKqsrFRoaKi/y0F7dOqU9JOfNHxdXS2FhPi3nlboM/d1f5fQrCNPj/V3CUC7wnGp7TW5Ty+AsR34sfZ6vOdYDzR2rsd7r+7xBgAAAAAA3iF4AwAAAABgEMEbAAAAAACDCN4AAAAAABhE8AYAAAAAwCCCNwAAAAAABhG8AQAAAAAwiOANAAAAAIBBBG8AAAAAAAwieAMAAAAAYBDBGwAAAAAAgwjeAAAAAAAYRPAGAAAAAMAggjcAAAAAAAYRvAEAAAAAMIjgDQAAAACAQQRvAAAAAAAMIngDAAAAAGAQwRsAAAAAAIMI3gAAAAAAGETwBgAAAADAIII3AAAAAAAGEbwBAAAAADCI4A0AAAAAgEEEbwAAAAAADCJ4AwAAAABgEMEbAAAAAACDCN4AAAAAABhE8AYAAAAAwCCCNwAAAAAABhG8AQAAAAAwiOANAAAAAIBBBG8AAAAAAAwieAMAgBbt3LlT48ePV2RkpGw2mzZv3uyed/r0aT3++OMaOHCgQkJCFBkZqcmTJ+uLL77wWEdNTY1mz56t8PBwhYSE6JZbbtHnn3/u4y0BAMA/CN4AAKBFp06d0qBBg7RixYpG877++msdOHBATz75pA4cOKBXXnlFH3/8sW655RaPfhkZGSosLNSmTZu0e/duVVdXa9y4caqrq/PVZgAA4DeB/i4AAAC0b2PGjNGYMWOanOdwOFRUVOTRtnz5cl1zzTUqKytTdHS0Kisr9cILL+ill17SjTfeKElav369oqKitH37do0ePdr4NgAA4E+c8QYAAG2qsrJSNptNl1xyiSSpuLhYp0+fVmpqqrtPZGSk4uLitGfPnmbXU1NTo6qqKo8JAICOiOANAADazLfffqu5c+dq0qRJCg0NlSS5XC4FBQWpe/fuHn0jIiLkcrmaXVdOTo4cDod7ioqKMlo7AACmELwBAECbOH36tO68807V19crLy/vrP0ty5LNZmt2flZWliorK93T0aNH27JcAAB8huANAADO2+nTpzVx4kSVlpaqqKjIfbZbkpxOp2pra/Xll196LFNRUaGIiIhm12m32xUaGuoxAQDQERG8AQDAeTkTuj/55BNt375dYWFhHvMTEhLUuXNnj4ewlZeX69ChQ0pOTvZ1uQAA+BxPNQcAAC2qrq7Wp59+6v5cWlqqkpIS9ejRQ5GRkbrjjjt04MAB/fGPf1RdXZ37vu0ePXooKChIDodDU6ZM0cMPP6ywsDD16NFDjzzyiAYOHOh+yjkAABcygjcAAGjR/v37dcMNN7g/Z2ZmSpLuueceZWdna8uWLZKkq6++2mO5t956SyNGjJAkLV26VIGBgZo4caK++eYbjRw5Uvn5+QoICPDJNgAA4E8EbwAA0KIRI0bIsqxm57c074zg4GAtX75cy5cvb8vSAADoEFp1j3deXp5iY2MVHByshIQE7dq1q9m+u3fvVkpKisLCwtSlSxf169dPS5cubdTv5Zdf1oABA2S32zVgwAAVFha2pjQAAAAAANoVr4N3QUGBMjIyNG/ePB08eFDXXXedxowZo7Kysib7h4SEaNasWdq5c6c+/PBDPfHEE3riiSe0Zs0ad5+9e/cqLS1N6enpeu+995Senq6JEyfq3Xffbf2WAQAAAADQDtisc7k+7AeGDh2qwYMHa9WqVe62/v3769Zbb1VOTs45reP2229XSEiIXnrpJUlSWlqaqqqq9MYbb7j73HTTTerevbs2btx4TuusqqqSw+FQZWUlrxtB006dkn7yk4avq6ulkBD/1tMKfea+7u8SmnXk6bH+LgFoVzgutb0m9+kFMLYDP9Zej/cc64HGzvV479UZ79raWhUXFys1NdWjPTU1VXv27DmndRw8eFB79uzR8OHD3W179+5ttM7Ro0e3uM6amhpVVVV5TAAAAAAAtDdeBe/jx4+rrq5OERERHu0RERHuV4c0p3fv3rLb7UpMTNTMmTM1depU9zyXy+X1OnNycuRwONxTVFSUN5sCAAAAAIBPtOrhajabzeOzZVmN2n5s165d2r9/v1avXq3c3NxGl5B7u86srCxVVla6p6NHj3q5FQAAAAAAmOfV68TCw8MVEBDQ6Ex0RUVFozPWPxYbGytJGjhwoP71r38pOztbd911lyTJ6XR6vU673S673e5N+QAAAAAA+JxXZ7yDgoKUkJCgoqIij/aioiIlJyef83osy1JNTY37c1JSUqN1btu2zat1AgAAAADQHnl1xluSMjMzlZ6ersTERCUlJWnNmjUqKyvT9OnTJTVcAn7s2DGtW7dOkrRy5UpFR0erX79+khre6/3MM89o9uzZ7nXOmTNH119/vRYtWqQJEybo1Vdf1fbt27V79+622EYAAAAAAPzG6+CdlpamEydOaMGCBSovL1dcXJy2bt2qmJgYSVJ5ebnHO73r6+uVlZWl0tJSBQYG6vLLL9fTTz+tadOmufskJydr06ZNeuKJJ/Tkk0/q8ssvV0FBgYYOHdoGmwgAAAAAgP94/R7v9or3peKs2uG7Xtvrezp9hfeB4kLGcant8R5vdEQc6znW48Jm5D3eAAAAAADAOwRvAAAAAAAMIngDAAAAAGAQwRsAAAAAAIMI3gAAAAAAGETwBgAAAADAIII3AAAAAAAGEbwBAAAAADCI4A0AAAAAgEEEbwAAAAAADCJ4AwAAAABgEMEbAAAAAACDCN4AAAAAABhE8AYAAAAAwCCCNwAAAAAABhG8AQAAAAAwiOANAAAAAIBBBG8AAAAAAAwieAMAAAAAYBDBGwAAAAAAgwjeAAAAAAAYRPAGAAAAAMAggjcAAAAAAAYRvAEAAAAAMIjgDQAAAACAQQRvAAAAAAAMIngDAAAAAGAQwRsAAAAAAIMI3gAAAAAAGETwBgAAAADAIII3AAAAAAAGEbwBAAAAADCI4A0AAAAAgEEEbwAAAAAADCJ4AwAAAABgEMEbAAAAAACDCN4AAKBFO3fu1Pjx4xUZGSmbzabNmzd7zLcsS9nZ2YqMjFSXLl00YsQIHT582KNPTU2NZs+erfDwcIWEhOiWW27R559/7sOtAADAfwjeAACgRadOndKgQYO0YsWKJucvXrxYS5Ys0YoVK7Rv3z45nU6NGjVKJ0+edPfJyMhQYWGhNm3apN27d6u6ulrjxo1TXV2drzYDAAC/CfR3AQAAoH0bM2aMxowZ0+Q8y7KUm5urefPm6fbbb5ckvfjii4qIiNCGDRs0bdo0VVZW6oUXXtBLL72kG2+8UZK0fv16RUVFafv27Ro9erTPtgUAAH/gjDcAAGi10tJSuVwupaamutvsdruGDx+uPXv2SJKKi4t1+vRpjz6RkZGKi4tz92lKTU2NqqqqPCYAADoigjcAAGg1l8slSYqIiPBoj4iIcM9zuVwKCgpS9+7dm+3TlJycHDkcDvcUFRXVxtUDAOAbBG8AAHDebDabx2fLshq1/djZ+mRlZamystI9HT16tE1qBQDA1wjeAACg1ZxOpyQ1OnNdUVHhPgvudDpVW1urL7/8stk+TbHb7QoNDfWYAADoiAjeAACg1WJjY+V0OlVUVORuq62t1Y4dO5ScnCxJSkhIUOfOnT36lJeX69ChQ+4+AABcyHiqOQAAaFF1dbU+/fRT9+fS0lKVlJSoR48eio6OVkZGhhYuXKi+ffuqb9++Wrhwobp27apJkyZJkhwOh6ZMmaKHH35YYWFh6tGjhx555BENHDjQ/ZRzAAAuZK06452Xl6fY2FgFBwcrISFBu3btarbvK6+8olGjRunSSy9VaGiokpKS9Oabb3r0yc/Pl81mazR9++23rSkPAAC0of379ys+Pl7x8fGSpMzMTMXHx+vXv/61JOmxxx5TRkaGZsyYocTERB07dkzbtm1Tt27d3OtYunSpbr31Vk2cOFEpKSnq2rWrXnvtNQUEBPhlmwAA8CWvz3gXFBQoIyNDeXl5SklJ0fPPP68xY8bogw8+UHR0dKP+O3fu1KhRo7Rw4UJdcsklWrt2rcaPH693333XfQCXpNDQUH300UceywYHB7dikwAAQFsaMWKELMtqdr7NZlN2drays7Ob7RMcHKzly5dr+fLlBioEAKB98zp4L1myRFOmTNHUqVMlSbm5uXrzzTe1atUq5eTkNOqfm5vr8XnhwoV69dVX9dprr3kEb5vN5n5ACwAAAAAAFwqvLjWvra1VcXGxUlNTPdpTU1O1Z8+ec1pHfX29Tp48qR49eni0V1dXKyYmRr1799a4ceN08ODBFtdTU1OjqqoqjwkAAAAAgPbGq+B9/Phx1dXVNXr1R0RERKPXiDTn2Wef1alTpzRx4kR3W79+/ZSfn68tW7Zo48aNCg4OVkpKij755JNm15OTkyOHw+GeoqKivNkUAAAAAAB8olUPV7PZbB6fLctq1NaUjRs3Kjs7WwUFBerZs6e7fdiwYbr77rs1aNAgXXfddfrf//1f/fSnP23xPrCsrCxVVla6p6NHj7ZmUwAAAAAAMMqre7zDw8MVEBDQ6Ox2RUVFo7PgP1ZQUKApU6boD3/4w1lfHdKpUycNGTKkxTPedrtddrv93IsHAAAAAMAPvDrjHRQUpISEBBUVFXm0FxUVKTk5udnlNm7cqHvvvVcbNmzQ2LFjz/p9LMtSSUmJevXq5U15AAAAAAC0O14/1TwzM1Pp6elKTExUUlKS1qxZo7KyMk2fPl1SwyXgx44d07p16yQ1hO7Jkyfrueee07Bhw9xny7t06SKHwyFJmj9/voYNG6a+ffuqqqpKy5YtU0lJiVauXNlW2wkAAAAAgF94HbzT0tJ04sQJLViwQOXl5YqLi9PWrVsVExMjSSovL1dZWZm7//PPP6/vvvtOM2fO1MyZM93t99xzj/Lz8yVJX331lR544AG5XC45HA7Fx8dr586duuaaa85z8wAAAAAA8C+bZVmWv4toC1VVVXI4HKqsrFRoaKi/y0F7dOqU9JOfNHxdXS2FhPi3Hkl95r7u7xL86sjTZ7/1BOioOC61vSb3aTsc24Ef4ljPsR4XtnM93rfqqeYAAAAAAODcELwBAAAAADCI4A0AAAAAgEEEbwAAAAAADCJ4AwAAAABgEMEbAAAAAACDCN4AAAAAABhE8AYAAAAAwCCCNwAAAAAABhG8AQAAAAAwiOANAAAAAIBBBG8AAAAAAAwieAMAAAAAYBDBGwAAAAAAgwjeAAAAAAAYRPAGAAAAAMAggjcAAAAAAAYRvAEAAAAAMIjgDQAAAACAQQRvAAAAAAAMIngDAAAAAGAQwRsAAAAAAIMI3gAAAAAAGETwBgAAAADAIII3AAAAAAAGEbwBAAAAADCI4A0AAAAAgEEEbwAAAAAADCJ4AwAAAABgEMEbAAAAAACDCN4AAAAAABhE8AYAAAAAwCCCNwAAAAAABhG8AQAAAAAwiOANAAAAAIBBBG8AAAAAAAwieAMAAAAAYBDBGwAAnLfvvvtOTzzxhGJjY9WlSxdddtllWrBggerr6919LMtSdna2IiMj1aVLF40YMUKHDx/2Y9UAAPgGwRsAAJy3RYsWafXq1VqxYoU+/PBDLV68WL/73e+0fPlyd5/FixdryZIlWrFihfbt2yen06lRo0bp5MmTfqwcAADzCN4AAOC87d27VxMmTNDYsWPVp08f3XHHHUpNTdX+/fslNZztzs3N1bx583T77bcrLi5OL774or7++mtt2LDBz9UDAGAWwRsAAJy3a6+9Vn/+85/18ccfS5Lee+897d69WzfffLMkqbS0VC6XS6mpqe5l7Ha7hg8frj179jS5zpqaGlVVVXlMAAB0RIH+LgAAAHR8jz/+uCorK9WvXz8FBASorq5OTz31lO666y5JksvlkiRFRER4LBcREaHPPvusyXXm5ORo/vz5ZgsHAMAHOOMNAADOW0FBgdavX68NGzbowIEDevHFF/XMM8/oxRdf9Ohns9k8PluW1ajtjKysLFVWVrqno0ePGqsfAACTOOMNAADO26OPPqq5c+fqzjvvlCQNHDhQn332mXJycnTPPffI6XRKajjz3atXL/dyFRUVjc6Cn2G322W3280XDwCAYQRvoI30mfu6v0vocFqzz448PdZAJQDO19dff61OnTwvpAsICHC/Tiw2NlZOp1NFRUWKj4+XJNXW1mrHjh1atGiRz+sFWoNjvfc41gMNWnWpeV5enmJjYxUcHKyEhATt2rWr2b6vvPKKRo0apUsvvVShoaFKSkrSm2++2ajfyy+/rAEDBshut2vAgAEqLCxsTWkAAMAPxo8fr6eeekqvv/66jhw5osLCQi1ZskS33XabpIZLzDMyMrRw4UIVFhbq0KFDuvfee9W1a1dNmjTJz9UDAGCW18G7oKBAGRkZmjdvng4ePKjrrrtOY8aMUVlZWZP9d+7cqVGjRmnr1q0qLi7WDTfcoPHjx+vgwYPuPnv37lVaWprS09P13nvvKT09XRMnTtS7777b+i0DAAA+s3z5ct1xxx2aMWOG+vfvr0ceeUTTpk3Tb3/7W3efxx57TBkZGZoxY4YSExN17Ngxbdu2Td26dfNj5QAAmGezLMvyZoGhQ4dq8ODBWrVqlbutf//+uvXWW5WTk3NO67jyyiuVlpamX//615KktLQ0VVVV6Y033nD3uemmm9S9e3dt3LjxnNZZVVUlh8OhyspKhYaGerFFuGicOiX95CcNX1dXSyEhbbp6Lj/zDS4/Q0fBcantNblPDY/twA9xrPcNjvXoSM71eO/VGe/a2loVFxd7vINTklJTU5t9B+eP1dfX6+TJk+rRo4e7be/evY3WOXr06BbXybs9AQAAAAAdgVfB+/jx46qrq2vyHZxn3s95Ns8++6xOnTqliRMnuttcLpfX68zJyZHD4XBPUVFRXmwJAAAAAAC+0aqHq3nzDs4f2rhxo7Kzs1VQUKCePXue1zp5tycAAAAAoCPw6nVi4eHhCggIaHQmuqV3cJ5RUFCgKVOm6A9/+INuvPFGj3lOp9PrdfJuTwAAAABAR+DVGe+goCAlJCSoqKjIo72oqEjJycnNLrdx40bde++92rBhg8aObfywhKSkpEbr3LZtW4vrBAAAAACgI/DqjLckZWZmKj09XYmJiUpKStKaNWtUVlam6dOnS2q4BPzYsWNat26dpIbQPXnyZD333HMaNmyY+8x2ly5d5HA4JElz5szR9ddfr0WLFmnChAl69dVXtX37du3evbutthMAAAAAAL/w+h7vtLQ05ebmasGCBbr66qu1c+dObd26VTExMZKk8vJyj3d6P//88/ruu+80c+ZM9erVyz3NmTPH3Sc5OVmbNm3S2rVrddVVVyk/P18FBQUaOnRoG2wiAAAAAAD+4/UZb0maMWOGZsyY0eS8/Px8j89vv/32Oa3zjjvu0B133NGacgAAAAAAaLda9VRzAAAAAABwbgjeAAAAAAAYRPAGAAAAAMAggjcAAAAAAAYRvAEAAAAAMIjgDQAAAACAQQRvAAAAAAAMIngDAAAAAGAQwRsAAAAAAIMI3gAAAAAAGETwBgAAAADAIII3AAAAAAAGEbwBAAAAADCI4A0AAAAAgEEEbwAAAAAADCJ4AwAAAABgEMEbAAAAAACDCN4AAAAAABhE8AYAAAAAwCCCNwAAAAAABhG8AQAAAAAwiOANAAAAAIBBBG8AAAAAAAwieAMAAAAAYBDBGwAAAAAAgwjeAAAAAAAYRPAGAAAAAMAggjcAAAAAAAYRvAEAAAAAMIjgDQAAAACAQQRvAAAAAAAMIngDAAAAAGAQwRsAAAAAAIMI3gAAAAAAGETwBgAAAADAIII3AAAAAAAGEbwBAAAAADCI4A0AAAAAgEEEbwAAAAAADCJ4AwCANnHs2DHdfffdCgsLU9euXXX11VeruLjYPd+yLGVnZysyMlJdunTRiBEjdPjwYT9WDACAbxC8AQDAefvyyy+VkpKizp0764033tAHH3ygZ599Vpdccom7z+LFi7VkyRKtWLFC+/btk9Pp1KhRo3Ty5En/FQ4AgA8E+rsAAADQ8S1atEhRUVFau3atu61Pnz7ury3LUm5urubNm6fbb79dkvTiiy8qIiJCGzZs0LRp03xdMgAAPsMZbwAAcN62bNmixMRE/fznP1fPnj0VHx+v3//+9+75paWlcrlcSk1NdbfZ7XYNHz5ce/bsaXKdNTU1qqqq8pgAAOiICN4AAOC8/fOf/9SqVavUt29fvfnmm5o+fboefPBBrVu3TpLkcrkkSRERER7LRUREuOf9WE5OjhwOh3uKiooyuxEAABjSquCdl5en2NhYBQcHKyEhQbt27Wq2b3l5uSZNmqQrrrhCnTp1UkZGRqM++fn5stlsjaZvv/22NeUBAAAfq6+v1+DBg7Vw4ULFx8dr2rRpuv/++7Vq1SqPfjabzeOzZVmN2s7IyspSZWWlezp69Kix+gEAMMnr4F1QUKCMjAzNmzdPBw8e1HXXXacxY8aorKysyf41NTW69NJLNW/ePA0aNKjZ9YaGhqq8vNxjCg4O9rY8AADgB7169dKAAQM82vr37+/+/cDpdEpSo7PbFRUVjc6Cn2G32xUaGuoxAQDQEXkdvJcsWaIpU6Zo6tSp6t+/v3JzcxUVFdXoL9pn9OnTR88995wmT54sh8PR7HptNpucTqfHBAAAOoaUlBR99NFHHm0ff/yxYmJiJEmxsbFyOp0qKipyz6+trdWOHTuUnJzs01oBAPA1r4J3bW2tiouLPR6MIkmpqanNPhjlXFVXVysmJka9e/fWuHHjdPDgwRb788AVAADaj4ceekjvvPOOFi5cqE8//VQbNmzQmjVrNHPmTEkNf2DPyMjQwoULVVhYqEOHDunee+9V165dNWnSJD9XDwCAWV4F7+PHj6uurs6rB6Oci379+ik/P19btmzRxo0bFRwcrJSUFH3yySfNLsMDVwAAaD+GDBmiwsJCbdy4UXFxcfrtb3+r3Nxc/eIXv3D3eeyxx5SRkaEZM2YoMTFRx44d07Zt29StWzc/Vg4AgHmteo+3Nw9GORfDhg3TsGHD3J9TUlI0ePBgLV++XMuWLWtymaysLGVmZro/V1VVEb4BAPCjcePGady4cc3Ot9lsys7OVnZ2tu+KAgCgHfAqeIeHhysgIMCrB6O0RqdOnTRkyJAWz3jb7XbZ7fY2+54AAAAAAJjg1aXmQUFBSkhI8HgwiiQVFRW16YNRLMtSSUmJevXq1WbrBAAAAADAH7y+1DwzM1Pp6elKTExUUlKS1qxZo7KyMk2fPl1SwyXgx44d07p169zLlJSUSGp4gNq///1vlZSUKCgoyP3akfnz52vYsGHq27evqqqqtGzZMpWUlGjlypVtsIkAAAAAAPiP18E7LS1NJ06c0IIFC1ReXq64uDht3brV/bqQ8vLyRu/0jo+Pd39dXFysDRs2KCYmRkeOHJEkffXVV3rggQfkcrnkcDgUHx+vnTt36pprrjmPTQMAAAAAwP9a9XC1GTNmaMaMGU3Oy8/Pb9RmWVaL61u6dKmWLl3amlIAAAAAAGjXvLrHGwAAAAAAeIfgDQAAAACAQQRvAAAAAAAMIngDAAAAAGAQwRsAAAAAAIMI3gAAAAAAGETwBgAAAADAIII3AAAAAAAGEbwBAAAAADCI4A0AAAAAgEEEbwAAAAAADCJ4AwAAAABgEMEbAAAAAACDCN4AAAAAABhE8AYAAAAAwCCCNwAAAAAABhG8AQAAAAAwiOANAAAAAIBBBG8AAAAAAAwieAMAAAAAYBDBGwAAAAAAgwjeAAAAAAAYRPAGAAAAAMAggjcAAAAAAAYRvAEAAAAAMCjQ3wUA7VGfua/7uwQ0w9t/myNPjzVUCQCgI+NY33615t+G4z3aO854AwAAAABgEMEbAAAAAACDCN4AAAAAABhE8AYAAAAAwCCCNwAAAAAABhG8AQAAAAAwiOANAAAAAIBBBG8AAAAAAAwieAMAAAAAYBDBGwAAAAAAgwjeAAAAAAAYRPAGAAAAAMAggjcAAAAAAAYRvAEAAAAAMIjgDQAAAACAQQRvAADQpnJycmSz2ZSRkeFusyxL2dnZioyMVJcuXTRixAgdPnzYf0UCAOBDBG8AANBm9u3bpzVr1uiqq67yaF+8eLGWLFmiFStWaN++fXI6nRo1apROnjzpp0oBAPAdgjcAAGgT1dXV+sUvfqHf//736t69u7vdsizl5uZq3rx5uv322xUXF6cXX3xRX3/9tTZs2ODHigEA8A2CNwAAaBMzZ87U2LFjdeONN3q0l5aWyuVyKTU11d1mt9s1fPhw7dmzp9n11dTUqKqqymMCAKAjalXwzsvLU2xsrIKDg5WQkKBdu3Y127e8vFyTJk3SFVdcoU6dOnnc7/VDL7/8sgYMGCC73a4BAwaosLCwNaUBAAA/2LRpkw4cOKCcnJxG81wulyQpIiLCoz0iIsI9ryk5OTlyOBzuKSoqqm2LBgDAR7wO3gUFBcrIyNC8efN08OBBXXfddRozZozKysqa7F9TU6NLL71U8+bN06BBg5rss3fvXqWlpSk9PV3vvfee0tPTNXHiRL377rvelgcAAHzs6NGjmjNnjtavX6/g4OBm+9lsNo/PlmU1avuhrKwsVVZWuqejR4+2Wc0AAPiS18F7yZIlmjJliqZOnar+/fsrNzdXUVFRWrVqVZP9+/Tpo+eee06TJ0+Ww+Fosk9ubq5GjRqlrKws9evXT1lZWRo5cqRyc3O9LQ8AAPhYcXGxKioqlJCQoMDAQAUGBmrHjh1atmyZAgMD3We6f3x2u6KiotFZ8B+y2+0KDQ31mAAA6Ii8Ct61tbUqLi72uEdLklJTU1u8R+ts9u7d22ido0eP5r4vAAA6gJEjR+r9999XSUmJe0pMTNQvfvELlZSU6LLLLpPT6VRRUZF7mdraWu3YsUPJycl+rBwAAN8I9Kbz8ePHVVdX5/U9Wmfjcrladd/X/PnzW/09AQBA2+jWrZvi4uI82kJCQhQWFuZuz8jI0MKFC9W3b1/17dtXCxcuVNeuXTVp0iR/lAwAgE95FbzP8PYeLRPrzMrKUmZmpvtzVVUVD10BAKCdeuyxx/TNN99oxowZ+vLLLzV06FBt27ZN3bp183dpAAAY51XwDg8PV0BAgNf3aJ2N0+ls1X1fdru91d8TAACY8/bbb3t8ttlsys7OVnZ2tl/qAQDAn7y6xzsoKEgJCQke92hJUlFR0Xndo5WUlNRondu2beO+LwAAAABAh+f1peaZmZlKT09XYmKikpKStGbNGpWVlWn69OmSGi4BP3bsmNatW+depqSkRJJUXV2tf//73yopKVFQUJAGDBggSZozZ46uv/56LVq0SBMmTNCrr76q7du3a/fu3W2wiQAAAAAA+I/XwTstLU0nTpzQggULVF5erri4OG3dulUxMTGSpPLy8kbv9I6Pj3d/XVxcrA0bNigmJkZHjhyRJCUnJ2vTpk164okn9OSTT+ryyy9XQUGBhg4deh6bBgAAAACA/7Xq4WozZszQjBkzmpyXn5/fqM2yrLOu84477tAdd9zRmnIAAAAAAGi3vLrHGwAAAAAAeIfgDQAAAACAQQRvAAAAAAAMIngDAAAAAGAQwRsAAAAAAIMI3gAAAAAAGETwBgAAAADAIII3AAAAAAAGEbwBAAAAADCI4A0AAAAAgEEEbwAAAAAADCJ4AwAAAABgEMEbAAAAAACDCN4AAAAAABhE8AYAAAAAwCCCNwAAAAAABhG8AQAAAAAwiOANAAAAAIBBBG8AAAAAAAwieAMAAAAAYBDBGwAAAAAAgwjeAAAAAAAYRPAGAAAAAMAggjcAAAAAAAYRvAEAAAAAMIjgDQAAAACAQQRvAAAAAAAMIngDAAAAAGAQwRsAAAAAAIMI3gAAAAAAGETwBgAAAADAIII3AAAAAAAGEbwBAAAAADCI4A0AAAAAgEEEbwAAAAAADAr0dwGAaX3mvi5J6lL7rT78/9v6P/knfRMU7L+i4DNn/v29ceTpsQYqAQCY0pqxHhcWb38GONbD1zjjDQAAAACAQQRvAAAAAAAMIngDAAAAAGAQwRsAAAAAAIMI3gAAAAAAGETwBgAAAADAIII3AAAAAAAGEbwBAMB5y8nJ0ZAhQ9StWzf17NlTt956qz766COPPpZlKTs7W5GRkerSpYtGjBihw4cP+6liAAB8p1XBOy8vT7GxsQoODlZCQoJ27drVYv8dO3YoISFBwcHBuuyyy7R69WqP+fn5+bLZbI2mb7/9tjXlAQAAH9uxY4dmzpypd955R0VFRfruu++UmpqqU6dOufssXrxYS5Ys0YoVK7Rv3z45nU6NGjVKJ0+e9GPlAACYF+jtAgUFBcrIyFBeXp5SUlL0/PPPa8yYMfrggw8UHR3dqH9paaluvvlm3X///Vq/fr3++te/asaMGbr00kv1s5/9zN0vNDS00V/Gg4ODW7FJAADA1/70pz95fF67dq169uyp4uJiXX/99bIsS7m5uZo3b55uv/12SdKLL76oiIgIbdiwQdOmTfNH2QAA+ITXZ7yXLFmiKVOmaOrUqerfv79yc3MVFRWlVatWNdl/9erVio6OVm5urvr376+pU6fqvvvu0zPPPOPRz2azyel0ekwAAKBjqqyslCT16NFDUsMf4l0ul1JTU9197Ha7hg8frj179jS5jpqaGlVVVXlMAAB0RF4F79raWhUXF3scNCUpNTW12YPm3r17G/UfPXq09u/fr9OnT7vbqqurFRMTo969e2vcuHE6ePCgN6UBAIB2wrIsZWZm6tprr1VcXJwkyeVySZIiIiI8+kZERLjn/VhOTo4cDod7ioqKMls4AACGeBW8jx8/rrq6Oq8Omi6Xq8n+3333nY4fPy5J6tevn/Lz87VlyxZt3LhRwcHBSklJ0SeffNJsLfwVHACA9mnWrFn6+9//ro0bNzaaZ7PZPD5bltWo7YysrCxVVla6p6NHjxqpFwAA07y+x1vy7qDZXP8ftg8bNkzDhg1zz09JSdHgwYO1fPlyLVu2rMl15uTkaP78+a0pHwAAGDJ79mxt2bJFO3fuVO/evd3tZ24hc7lc6tWrl7u9oqKi0R/oz7Db7bLb7WYLBgDAB7w64x0eHq6AgIBGZ7dbOmg6nc4m+wcGBiosLKzpojp10pAhQ1o8481fwQEAaD8sy9KsWbP0yiuv6C9/+YtiY2M95sfGxsrpdKqoqMjdVltbqx07dig5OdnX5QIA4FNeBe+goCAlJCR4HDQlqaioqNmDZlJSUqP+27ZtU2Jiojp37tzkMpZlqaSkxOMv4j9mt9sVGhrqMQEAAP+YOXOm1q9frw0bNqhbt25yuVxyuVz65ptvJDVc5ZaRkaGFCxeqsLBQhw4d0r333quuXbtq0qRJfq4eAACzvL7UPDMzU+np6UpMTFRSUpLWrFmjsrIyTZ8+XVLDmehjx45p3bp1kqTp06drxYoVyszM1P3336+9e/fqhRde8Ljva/78+Ro2bJj69u2rqqoqLVu2TCUlJVq5cmUbbSYAADDpzNtNRowY4dG+du1a3XvvvZKkxx57TN98841mzJihL7/8UkOHDtW2bdvUrVs3H1cLAIBveR2809LSdOLECS1YsEDl5eWKi4vT1q1bFRMTI0kqLy9XWVmZu39sbKy2bt2qhx56SCtXrlRkZKSWLVvm8Q7vr776Sg888IBcLpccDofi4+O1c+dOXXPNNW2wiQAAwLQzz29pic1mU3Z2trKzs80XBABAO9Kqh6vNmDFDM2bMaHJefn5+o7bhw4frwIEDza5v6dKlWrp0aWtKAQAAAACgXfPqHm8AAAAAAOAdgjcAAAAAAAYRvAEAAAAAMIjgDQAAAACAQQRvAAAAAAAMIngDAAAAAGAQwRsAAAAAAIMI3gAAAAAAGETwBgAAAADAIII3AAAAAAAGEbwBAAAAADCI4A0AAAAAgEEEbwAAAAAADCJ4AwAAAABgEMEbAAAAAACDCN4AAAAAABhE8AYAAAAAwCCCNwAAAAAABgX6uwDAG33mvu7vEnARaM3P2ZGnxxqoBAAuThzvYRrHevgaZ7wBAAAAADCI4A0AAAAAgEEEbwAAAAAADCJ4AwAAAABgEMEbAAAAAACDCN4AAAAAABhE8AYAAAAAwCCCNwAAAAAABhG8AQAAAAAwiOANAAAAAIBBBG8AAAAAAAwieAMAAAAAYBDBGwAAAAAAgwjeAAAAAAAYRPAGAAAAAMAggjcAAAAAAAYRvAEAAAAAMIjgDQAAAACAQQRvAAAAAAAMCvR3Abi49Zn7ur9LANpEa36Wjzw91kAlANC+cKzHhYJjPc4HZ7wBAAAAADCI4A0AAAAAgEEEbwAAAAAADCJ4AwAAAABgEMEbAAAAAACDCN4AAAAAABjUquCdl5en2NhYBQcHKyEhQbt27Wqx/44dO5SQkKDg4GBddtllWr16daM+L7/8sgYMGCC73a4BAwaosLCwNaUBAIB2zNvfIQAAuBB4/R7vgoICZWRkKC8vTykpKXr++ec1ZswYffDBB4qOjm7Uv7S0VDfffLPuv/9+rV+/Xn/96181Y8YMXXrppfrZz34mSdq7d6/S0tL029/+VrfddpsKCws1ceJE7d69W0OHDj3/rWwF3tPnPd7TCXiHccY89nH74u3vEKbx8+E9jvWAdxhnfKMj7Gevz3gvWbJEU6ZM0dSpU9W/f3/l5uYqKipKq1atarL/6tWrFR0drdzcXPXv319Tp07Vfffdp2eeecbdJzc3V6NGjVJWVpb69eunrKwsjRw5Urm5ua3eMAAA0L54+zsEAAAXCq/OeNfW1qq4uFhz5871aE9NTdWePXuaXGbv3r1KTU31aBs9erReeOEFnT59Wp07d9bevXv10EMPNerTUvCuqalRTU2N+3NlZaUkqaqqyptNalZ9zddeL9NW37ujas0+86W62m915l+oruZr1Vv1fq0HaI2LfZzxlj/H8jPrsSyrTdbX0bXmd4hzOtafOvX911VVUl3dOdfEsd577f1YD1wILvZxpjU6wvHeq+B9/Phx1dXVKSIiwqM9IiJCLperyWVcLleT/b/77jsdP35cvXr1arZPc+uUpJycHM2fP79Re1RU1LluTptz5PrtW+McOc58kTfZn2UArcY4Y15b7+OTJ0/K4XCcveMFrjW/Q3h9rI+MPO86z4b/gwBMY5zxDV8f772+x1uSbDabx2fLshq1na3/j9u9XWdWVpYyMzPdn+vr6/X//t//U1hYmMdyVVVVioqK0tGjRxUaGtrCVl3Y2A8N2A/fY180YD80YD80aMv9YFmWTp48qUgfhMGOxJvjPcd677EvGrAfGrAfGrAfGrAfvueP471XwTs8PFwBAQGN/jJdUVHR6C/YZzidzib7BwYGKiwsrMU+za1Tkux2u+x2u0fbJZdc0mz/0NDQi/4HTGI/nMF++B77ogH7oQH7oUFb7QfOdH+vNb9DcKxvPfZFA/ZDA/ZDA/ZDA/bD93x5vPfq4WpBQUFKSEhQUVGRR3tRUZGSk5ObXCYpKalR/23btikxMVGdO3dusU9z6wQAAB1La36HAADgQuH1peaZmZlKT09XYmKikpKStGbNGpWVlWn69OmSGi4LO3bsmNatWydJmj59ulasWKHMzEzdf//92rt3r1544QVt3LjRvc45c+bo+uuv16JFizRhwgS9+uqr2r59u3bv3t1GmwkAAPztbL9DAABwofI6eKelpenEiRNasGCBysvLFRcXp61btyomJkaSVF5errKyMnf/2NhYbd26VQ899JBWrlypyMhILVu2zP0Ob0lKTk7Wpk2b9MQTT+jJJ5/U5ZdfroKCgjZ5h7fdbtdvfvObRpeqXWzYDw3YD99jXzRgPzRgPzRgP5h1tt8hWot/t++xLxqwHxqwHxqwHxqwH77nj31hs3jPCQAAAAAAxnh1jzcAAAAAAPAOwRsAAAAAAIMI3gAAAAAAGETwBgAAAADAoIsmeB85ckRTpkxRbGysunTpossvv1y/+c1vVFtb6+/SfCIvL0+xsbEKDg5WQkKCdu3a5e+SfConJ0dDhgxRt27d1LNnT91666366KOP/F2W3+Xk5MhmsykjI8PfpfjcsWPHdPfddyssLExdu3bV1VdfreLiYn+X5XPfffednnjiCffYeNlll2nBggWqr6/3d2lG7dy5U+PHj1dkZKRsNps2b97sMd+yLGVnZysyMlJdunTRiBEjdPjwYf8Ui2YxtjeNsZ2xnbH94h7bW9oPp0+f1uOPP66BAwcqJCREkZGRmjx5sr744gv/FWzI2X4efmjatGmy2WzKzc01Vs9FE7z/7//+T/X19Xr++ed1+PBhLV26VKtXr9avfvUrf5dmXEFBgTIyMjRv3jwdPHhQ1113ncaMGePx2rcL3Y4dOzRz5ky98847Kioq0nfffafU1FSdOnXK36X5zb59+7RmzRpdddVV/i7F57788kulpKSoc+fOeuONN/TBBx/o2Wef1SWXXOLv0nxu0aJFWr16tVasWKEPP/xQixcv1u9+9zstX77c36UZderUKQ0aNEgrVqxocv7ixYu1ZMkSrVixQvv27ZPT6dSoUaN08uRJH1eKljC2N8bYztguMbZf7GN7S/vh66+/1oEDB/Tkk0/qwIEDeuWVV/Txxx/rlltu8UOlZp3t5+GMzZs3691331VkZKTZgqyL2OLFi63Y2Fh/l2HcNddcY02fPt2jrV+/ftbcuXP9VJH/VVRUWJKsHTt2+LsUvzh58qTVt29fq6ioyBo+fLg1Z84cf5fkU48//rh17bXX+ruMdmHs2LHWfffd59F2++23W3fffbefKvI9SVZhYaH7c319veV0Oq2nn37a3fbtt99aDofDWr16tR8qxLlibGdsZ2xvwNjO2H7Gj/dDU/72t79ZkqzPPvvMN0X5QXP74fPPP7f+4z/+wzp06JAVExNjLV261FgNF80Z76ZUVlaqR48e/i7DqNraWhUXFys1NdWjPTU1VXv27PFTVf5XWVkpSRf8v39zZs6cqbFjx+rGG2/0dyl+sWXLFiUmJurnP/+5evbsqfj4eP3+97/3d1l+ce211+rPf/6zPv74Y0nSe++9p927d+vmm2/2c2X+U1paKpfL5TFu2u12DR8+/KIeNzsCxnbGdsb2BoztjTG2N6+yslI2m+2iuzqkvr5e6enpevTRR3XllVca/36Bxr9DO/WPf/xDy5cv17PPPuvvUow6fvy46urqFBER4dEeEREhl8vlp6r8y7IsZWZm6tprr1VcXJy/y/G5TZs26cCBA9q3b5+/S/Gbf/7zn1q1apUyMzP1q1/9Sn/729/04IMPym63a/Lkyf4uz6cef/xxVVZWql+/fgoICFBdXZ2eeuop3XXXXf4uzW/OjI1NjZufffaZP0rCOWBsZ2xnbP8eY3tjjO1N+/bbbzV37lxNmjRJoaGh/i7HpxYtWqTAwEA9+OCDPvl+HT54Z2dna/78+S322bdvnxITE92fv/jiC9100036+c9/rqlTp5ousV2w2Wweny3LatR2sZg1a5b+/ve/a/fu3f4uxeeOHj2qOXPmaNu2bQoODvZ3OX5TX1+vxMRELVy4UJIUHx+vw4cPa9WqVRfdL2cFBQVav369NmzYoCuvvFIlJSXKyMhQZGSk7rnnHn+X51eMmx0LYztjO2P79xjbm8fY/r3Tp0/rzjvvVH19vfLy8vxdjk8VFxfrueee04EDB3z279/hg/esWbN05513ttinT58+7q+/+OIL3XDDDUpKStKaNWsMV+d/4eHhCggIaHR2u6KiotFf/C4Gs2fP1pYtW7Rz50717t3b3+X4XHFxsSoqKpSQkOBuq6ur086dO7VixQrV1NQoICDAjxX6Rq9evTRgwACPtv79++vll1/2U0X+8+ijj2ru3LnucXTgwIH67LPPlJOTc9H+cuZ0OiU1nB3p1auXu/1iHTc7AsZ2xnaJsf2HGNsbY2z3dPr0aU2cOFGlpaX6y1/+ctGd7d61a5cqKioUHR3tbqurq9PDDz+s3NxcHTlypM2/Z4cP3uHh4QoPDz+nvseOHdMNN9yghIQErV27Vp06Xfi3uAcFBSkhIUFFRUW67bbb3O1FRUWaMGGCHyvzLcuyNHv2bBUWFurtt99WbGysv0vyi5EjR+r999/3aPvlL3+pfv366fHHH78ofjGTpJSUlEavHPr4448VExPjp4r85+uvv240FgYEBFzwr5xpSWxsrJxOp4qKihQfHy+p4XkZO3bs0KJFi/xcHX6Isb0BY3sDxvbvMbY3xtj+vTOh+5NPPtFbb72lsLAwf5fkc+np6Y2ehzF69Gilp6frl7/8pZHv2eGD97n64osvNGLECEVHR+uZZ57Rv//9b/e8M38Bu1BlZmYqPT1diYmJ7jP9ZWVlmj59ur9L85mZM2dqw4YNevXVV9WtWzf3FQAOh0NdunTxc3W+061bt0b3PoaEhCgsLOyiuifyoYceUnJyshYuXKiJEyfqb3/7m9asWXNRXAXzY+PHj9dTTz2l6OhoXXnllTp48KCWLFmi++67z9+lGVVdXa1PP/3U/bm0tFQlJSXq0aOHoqOjlZGRoYULF6pv377q27evFi5cqK5du2rSpEl+rBo/xtjegLG9AWP79xjbG1ysY3tL+yEyMlJ33HGHDhw4oD/+8Y+qq6tzj509evRQUFCQv8puc2f7efjxHxw6d+4sp9OpK664wkxBxp6X3s6sXbvWktTkdDFYuXKlFRMTYwUFBVmDBw++6F610ty//dq1a/1dmt9djK+csSzLeu2116y4uDjLbrdb/fr1s9asWePvkvyiqqrKmjNnjhUdHW0FBwdbl112mTVv3jyrpqbG36UZ9dZbbzU5Jtxzzz2WZTW8duY3v/mN5XQ6Lbvdbl1//fXW+++/79+i0Qhje/MY2xnbGdsv3rG9pf1QWlra7Nj51ltv+bv0NnW2n4cfM/06MZtlWZaZSA8AAAAAAC78m5wBAAAAAPAjgjcAAAAAAAYRvAEAAAAAMIjgDQAAAACAQQRvAAAAAAAMIngDAAAAAGAQwRsAAAAAAIMI3gAAAAAAGETwBgAAAADAIII3AAAAAAAGEbwBAAAAADCI4A0AAAAAgEH/H3tU04rxJfbZAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "# initialize Normal(0, 1)\n", + "MEAN = 4\n", + "VARIANCE = 1\n", + "std_norm = distribution_dict[\"normal\"](MEAN, VARIANCE)\n", + "print(std_norm.stats(moments=\"mv\"))\n", + "support_01 = np.linspace(0.00000001, 0.999999999, 1000)\n", + "quantiles = std_norm.ppf(support_01)\n", + "vals = reverse_z(quantiles, MEAN, VARIANCE)\n", + "\n", + "fig, ax = plt.subplots(1, 2, figsize=(12, 6))\n", + "ax[0].hist(quantiles, density=True, bins=30)\n", + "quantile = 2\n", + "print(std_norm.cdf(quantile))\n", + "ax[0].axvline(quantile, color=\"red\")\n", + "\n", + "ax[1].hist(vals, bins=30)\n", + "ax[1].axvline(quantile, color=\"red\")\n", + "plt.show()\n", + "\n", + "\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "* this shows that cdf and pdf take values from the real line, as if cdf were to take " + ] } ], "metadata": { diff --git a/src/ensemble/distributions.py b/src/ensemble/distributions.py index da3d242..fa96ad8 100644 --- a/src/ensemble/distributions.py +++ b/src/ensemble/distributions.py @@ -62,32 +62,18 @@ def _create_scipy_dist(self) -> None: self._scipy_dist = scipy.stats.gamma(a=alpha, scale=1 / beta) -# TODO: INVESTIGATE analytic SOLUTION +# analytic sol class InvGamma(Distribution): support = "strictly positive" def _create_scipy_dist(self) -> None: strict_positive_support(self.mean) - optim_params = scipy.optimize.minimize( - fun=self._shape_scale, - # a *good* friend told me that this is a good initial guess and it works so far??? - # alpha = 3 is because alpha > 2 must be true due to variance formula - # beta = mean * (alpha - 1) after isolating beta from formula for mean - x0=[3, self.mean * 2], - args=(self.mean, self.variance), - ) - shape, scale = np.abs(optim_params.x) - self._scipy_dist = scipy.stats.invgamma(a=shape, scale=scale) - - def _shape_scale(self, x: list, samp_mean: float, samp_var: float) -> None: - alpha = x[0] - beta = x[1] - mean_guess = beta / (alpha - 1) - variance_guess = beta**2 / ((alpha - 1) ** 2 * (alpha - 2)) - return (mean_guess - samp_mean) ** 2 + (variance_guess - samp_var) ** 2 + alpha = self.mean**2 / self.variance + 2 + beta = self.mean * (self.mean**2 / self.variance + 1) + self._scipy_dist = scipy.stats.invgamma(a=alpha, scale=beta) -# TODO: investigate analytic solution +# numerical sol class Fisk(Distribution): support = "positive" @@ -102,6 +88,7 @@ def _create_scipy_dist(self): ) alpha, beta = np.abs(optim_params.x) # parameterization notes: numpy's c is wikipedia's beta, numpy's scale is wikipedia's alpha + # additional note: analytical solution doesn't work b/c dependent on derivative self._scipy_dist = scipy.stats.fisk(c=beta, scale=alpha) def _shape_scale(self, x: list, samp_mean: float, samp_var: float) -> None: @@ -153,23 +140,23 @@ def _shape_scale(self, x: list, samp_mean: float, samp_var: float) -> float: return (mean_guess - samp_mean) ** 2 + (variance_guess - samp_var) ** 2 -# somewhat broken +# analytic sol (M.O.M. estimators) class LogNormal(Distribution): support = "strictly positive" def _create_scipy_dist(self) -> None: strict_positive_support(self.mean) - # using method of moments gets close, but not quite there - loc = np.log(self.mean / np.sqrt(1 + (self.variance / self.mean**2))) - scale = np.sqrt(np.log(1 + (self.variance / self.mean**2))) - # loc = np.log(self.mean**2 / np.sqrt(self.mean**2 + self.variance)) - # scale = np.log(1 + self.variance / self.mean**2) - self._scipy_dist = scipy.stats.lognorm(loc=loc, s=scale) + mu = np.log(self.mean / np.sqrt(1 + (self.variance / self.mean**2))) + sigma = np.sqrt(np.log(1 + (self.variance / self.mean**2))) + # scipy multiplies in the argument passed to `scale` so in the exponentiated space, + # you're essentially adding `mu` within the exponentiated expression within the + # lognormal's PDF; hence, scale is with exponentiation instead of loc + self._scipy_dist = scipy.stats.lognorm(scale=np.exp(mu), s=sigma) # analytic sol class Normal(Distribution): - support = "positive" + support = "real line" def _create_scipy_dist(self) -> None: self._scipy_dist = scipy.stats.norm( diff --git a/src/ensemble/ensemble_model.py b/src/ensemble/ensemble_model.py index 26b14fe..ebbc6e5 100644 --- a/src/ensemble/ensemble_model.py +++ b/src/ensemble/ensemble_model.py @@ -12,31 +12,43 @@ class EnsembleResult: weights: npt.NDArray ensemble_density: npt.NDArray + def __init__(self, weights, ensemble_density) -> None: + self.weights = weights + self.ensemble_density = ensemble_density + + +# class EnsembleModel: +# def __init__(self, distributions, weights): +# self.distributions = + class EnsembleFitter: def __init__(self, distributions: List[str], objective): self.supports = set() for distribution in distributions: - self.supports.add(distribution_dict[distribution]) + self.supports.add(distribution_dict[distribution].support) # TODO: HOW SHOULD WE TELL THE USER WHICH DISTRIBUTION IS THE "TROUBLEMAKER"? if len(self.supports) != 1: raise ValueError( "the provided list of distributions do not all have the same support: " - + self.supports + + str(self.supports) ) self.distributions = distributions + def objective(self, vec): + return scipy.linalg.norm(vec, 2) + def ensemble_obj(self, weights): # return data - F @ weights pass def get_ensemble_density( - self, cdfs: np.ndarray, fitted_weights: np.ndarray + self, pdfs: np.ndarray, fitted_weights: np.ndarray ): - return cdfs @ fitted_weights + return pdfs @ fitted_weights def ensemble_func(self, weights: list, ecdf: np.ndarray, cdfs: np.ndarray): - return ecdf - cdfs @ weights + return self.objective(ecdf - cdfs @ weights) def fit(self, data: npt.ArrayLike) -> EnsembleResult: # TODO: SWITCH CASE STATEMENT FOR BOUNDS OF DATA NOT MATCHING THE ELEMENT OF SELF.SUPPORTS @@ -46,30 +58,38 @@ def fit(self, data: npt.ArrayLike) -> EnsembleResult: ecdf = scipy.stats.ecdf(data).cdf.probabilities # may be able to condense into 1 line if ub and lb are not used elsewhere - support_lb = np.min(data) - support_ub = np.max(data) - support = np.linspace(support_lb, support_ub, len(data)) + # support_lb = np.min(data) + # support_ub = np.max(data) + # support = np.linspace(support_lb, support_ub, len(data)) + equantiles = scipy.stats.ecdf(data).cdf.quantiles # fill matrix with cdf values over support of data num_distributions = len(self.distributions) cdfs = np.zeros((len(data), num_distributions)) + pdfs = np.zeros((len(data), num_distributions)) for i in range(num_distributions): curr_dist = distribution_dict[self.distributions[i]]( sample_mean, sample_variance ) - cdfs[:, i] = curr_dist.cdf(support) + cdfs[:, i] = curr_dist.cdf(equantiles) + pdfs[:, i] = curr_dist.pdf(equantiles) # initialize equal weights for all dists and optimize initial_guess = np.zeros(num_distributions) + 1 / num_distributions + bounds = tuple((0, 1) for i in range(num_distributions)) minimizer_result = scipy.optimize.minimize( fun=self.ensemble_func, x0=initial_guess, args=(ecdf, cdfs), + bounds=bounds, ) + fitted_weights = minimizer_result.x + + # return minimizer_result.x res = EnsembleResult( - weights=minimizer_result.x, - ensemble_density=self.ensemble_density(cdfs, minimizer_result.x), + weights=fitted_weights, + ensemble_density=self.get_ensemble_density(cdfs, fitted_weights), ) return res diff --git a/tests/test_distributions.py b/tests/test_distributions.py index dd4bad7..b2c3ccd 100644 --- a/tests/test_distributions.py +++ b/tests/test_distributions.py @@ -22,8 +22,8 @@ NEG_MEAN = -2 BETA_MEAN = 0.5 BETA_VARIANCE = 0.2 -MEAN = 2 -VARIANCE = 8 +MEAN = 1 +VARIANCE = 2 def test_exp(): @@ -48,6 +48,7 @@ def test_invgamma(): assert np.isclose(res[1], VARIANCE) +# TODO: WRITE ADDITIONAL TESTS DUE TO NUMERICAL SOLUTION def test_fisk(): fisk = Fisk(MEAN, VARIANCE) res = fisk.stats(moments="mv")