From 06fbfdd88369ab3e714672ad350fa4e3bc62a384 Mon Sep 17 00:00:00 2001 From: beckynevin Date: Thu, 25 Jan 2024 21:22:17 -0600 Subject: [PATCH] added an inference fxn and also the notebook now loads and makes a corner plot --- notebooks/train_SBI.ipynb | 226 +++++++++++++++++++++++++++++++------- src/scripts/evaluate.py | 4 + 2 files changed, 190 insertions(+), 40 deletions(-) diff --git a/notebooks/train_SBI.ipynb b/notebooks/train_SBI.ipynb index 057830d..760f93f 100644 --- a/notebooks/train_SBI.ipynb +++ b/notebooks/train_SBI.ipynb @@ -11,32 +11,36 @@ }, { "cell_type": "code", - "execution_count": 1, + "execution_count": 9, + "id": "ebb1cf48-c144-4a56-b46a-edef83f443fa", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "import matplotlib\n", + "# remove top and right axis from plots\n", + "matplotlib.rcParams[\"axes.spines.right\"] = False\n", + "matplotlib.rcParams[\"axes.spines.top\"] = False" + ] + }, + { + "cell_type": "code", + "execution_count": 18, "id": "486dda47-bf7b-45ea-88fe-55960d81c4bb", "metadata": {}, - "outputs": [ - { - "ename": "ModuleNotFoundError", - "evalue": "No module named 'sbi'", - "output_type": "error", - "traceback": [ - "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", - "\u001b[0;31mModuleNotFoundError\u001b[0m Traceback (most recent call last)", - "Cell \u001b[0;32mIn[1], line 1\u001b[0m\n\u001b[0;32m----> 1\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01msbi\u001b[39;00m\n\u001b[1;32m 2\u001b[0m \u001b[38;5;66;03m# from sbi import inference\u001b[39;00m\n\u001b[1;32m 3\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01msbi\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01minference\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m SNPE\n", - "\u001b[0;31mModuleNotFoundError\u001b[0m: No module named 'sbi'" - ] - } - ], + "outputs": [], "source": [ "import sbi\n", "from sbi.inference import SNPE\n", "from sbi.inference.base import infer\n", + "from sbi.analysis import pairplot\n", "import torch" ] }, { "cell_type": "code", - "execution_count": 2, + "execution_count": 1, "id": "f64b72b1-3c46-45af-932e-59512b2adbc8", "metadata": {}, "outputs": [], @@ -49,7 +53,7 @@ }, { "cell_type": "code", - "execution_count": 3, + "execution_count": 6, "id": "cd5034fb-94da-4b3d-b5ca-89f16ed98ff9", "metadata": {}, "outputs": [], @@ -68,22 +72,10 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": 9, "id": "9fe446e0-e80e-4c6a-a67e-8a8bd19d2787", "metadata": {}, - "outputs": [ - { - "ename": "NameError", - "evalue": "name 'torch' is not defined", - "output_type": "error", - "traceback": [ - "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", - "\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)", - "Cell \u001b[0;32mIn[4], line 3\u001b[0m\n\u001b[1;32m 1\u001b[0m num_dim \u001b[38;5;241m=\u001b[39m \u001b[38;5;241m2\u001b[39m\n\u001b[0;32m----> 3\u001b[0m low_bounds \u001b[38;5;241m=\u001b[39m \u001b[43mtorch\u001b[49m\u001b[38;5;241m.\u001b[39mtensor([\u001b[38;5;241m0\u001b[39m, \u001b[38;5;241m-\u001b[39m\u001b[38;5;241m10\u001b[39m])\n\u001b[1;32m 4\u001b[0m high_bounds \u001b[38;5;241m=\u001b[39m torch\u001b[38;5;241m.\u001b[39mtensor([\u001b[38;5;241m10\u001b[39m, \u001b[38;5;241m10\u001b[39m])\n\u001b[1;32m 6\u001b[0m prior \u001b[38;5;241m=\u001b[39m sbi\u001b[38;5;241m.\u001b[39mutils\u001b[38;5;241m.\u001b[39mBoxUniform(low \u001b[38;5;241m=\u001b[39m low_bounds, high \u001b[38;5;241m=\u001b[39m high_bounds)\n", - "\u001b[0;31mNameError\u001b[0m: name 'torch' is not defined" - ] - } - ], + "outputs": [], "source": [ "num_dim = 2\n", "\n", @@ -95,38 +87,192 @@ }, { "cell_type": "code", - "execution_count": 5, + "execution_count": 12, "id": "b4d1c9af-cedc-483b-92ba-8bb1d195bccf", "metadata": {}, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "7d46e37a2a3c473f9fe9bc96c50b46c2", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "Running 10000 simulations.: 0%| | 0/10000 [00:00 1\u001b[0m posterior \u001b[38;5;241m=\u001b[39m \u001b[43minfer\u001b[49m(simulator, prior, \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mSNPE\u001b[39m\u001b[38;5;124m\"\u001b[39m, num_simulations\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m10000\u001b[39m)\n", - "\u001b[0;31mNameError\u001b[0m: name 'infer' is not defined" + "Cell \u001b[0;32mIn[2], line 4\u001b[0m\n\u001b[1;32m 2\u001b[0m path \u001b[38;5;241m=\u001b[39m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m../savedmodels/sbi/\u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m 3\u001b[0m model_name \u001b[38;5;241m=\u001b[39m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124msbi_linear\u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[0;32m----> 4\u001b[0m inference_model\u001b[38;5;241m.\u001b[39msave_model_pkl(path, model_name, \u001b[43mposterior\u001b[49m)\n", + "\u001b[0;31mNameError\u001b[0m: name 'posterior' is not defined" ] } ], "source": [ - "posterior = infer(simulator, prior, \"SNPE\", num_simulations=10000)" + "inference_model = evaluate.InferenceModel()\n", + "path = \"../savedmodels/sbi/\"\n", + "model_name = \"sbi_linear\"\n", + "inference_model.save_model_pkl(path, model_name, posterior)" + ] + }, + { + "cell_type": "markdown", + "id": "3dcfdf32-7ced-4380-8047-46a507eb0de6", + "metadata": {}, + "source": [ + "Test that this worked." ] }, { "cell_type": "code", - "execution_count": null, - "id": "07042d00-c9b1-494b-8870-d93ad18e2e11", + "execution_count": 13, + "id": "c4d1dbcb-2baf-478d-a3f2-d1c31228b7ab", "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "../savedmodels/sbi/\n" + ] + } + ], "source": [ "inference_model = evaluate.InferenceModel()\n", - "path = \"saved_models/sbi/\"\n", + "path = \"../savedmodels/sbi/\"\n", "model_name = \"sbi_linear\"\n", - "inference_model.save_model_pkl(self, path, model_name, posterior)" + "posterior = inference_model.load_model_pkl(path, model_name)" + ] + }, + { + "cell_type": "markdown", + "id": "687c0227-aed3-410a-85a3-42757bc7dfe4", + "metadata": {}, + "source": [ + "Run inference on the posterior." + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "faa33718-1c47-428b-9e94-23d4ce59c5ca", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjsAAAGwCAYAAABPSaTdAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA/nElEQVR4nO3df3RU9Z3/8deQSKAiSYMhAWYwtGWLrdZaUYqalqw5R/vDxobUI6BS14NHF2oCbqvUVWErjVvPYoJa3Xq66qkkiiT+KP11bJhoLCm/LFZri3QbNAYSqCwJqCBMPt8//M6UCTPJTObOnXvvPB/nzDmbO3dmPnOrzms/9/15f3zGGCMAAACPGpXpAQAAAKQTYQcAAHgaYQcAAHgaYQcAAHgaYQcAAHgaYQcAAHgaYQcAAHgaYUeSMUb9/f2i5RAAAN5D2JF06NAh5efn69ChQ5keCgAAsBhhBwAAeBphBwAAeBphBwAAeBphBwAAeBphBwAAeBphBwAAeBphBwAAeBphBwAAeBphBwAAeBphBwAAeBphBwAAeFpupgcAAAC8KRQKqb29XXv37tWkSZNUVlamnJwc28dB2AEAAJZraWlRTU2N3nnnncgxv9+vhoYGVVVV2ToWbmMBAABLtbS0qLq6OiroSFJ3d7eqq6vV0tJi63h8xhhj6yc6UH9/v/Lz89XX16fx48dnejgAALhWKBRSaWnpSUEnzOfzye/3q7Oz07ZbWszsAAAAy7S3t8cNOpJkjFFXV5fa29ttGxNhBwAAWGbv3r2WnmcFwg4AALDMpEmTLD3PCoQdAABgmbKyMvn9fvl8vpjP+3w+BQIBlZWV2TYmwg4AALBMTk6OGhoaJOmkwBP+u76+3tZ+O4QdAABgqaqqKq1fv15TpkyJOu73+7V+/Xrb++yw9FwsPQcAIB2c0kGZsCPCDgAAXsZtLAAA4GmEHQAA4GmEHQAA4GmEHQAA4GmEHQAA4GkZDTsvvfSSLr/8ck2ePFk+n0/PPvts5Lljx47p1ltv1dlnn61TTz1VkydP1rXXXqs9e/ZEvceBAwe0YMECjR8/XgUFBbr++ut1+PBhm78JAABwqoyGnffee0/nnHOOHnzwwZOee//99/XKK6/ojjvu0CuvvKKWlhbt3LlT3/jGN6LOW7Bggf70pz/phRde0IYNG/TSSy/phhtusOsrAAAAh3NMnx2fz6dnnnlGV1xxRdxztm7dqgsuuEBvvfWWpk6dqj//+c/6zGc+o61bt2rmzJmSpF//+tf66le/qnfeeUeTJ09O6LPpswMAgHe5qmanr69PPp9PBQUFkqSOjg4VFBREgo4kVVRUaNSoUdq8eXPc9zl69Kj6+/ujHgAAwJtcE3aOHDmiW2+9VfPmzYvMvvT09GjixIlR5+Xm5qqwsFA9PT1x36uurk75+fmRRyAQSOvYAQBA5rgi7Bw7dkxXXnmljDF66KGHUn6/5cuXq6+vL/Lo6uqyYJQAAMCJcjM9gOGEg85bb72ljRs3RtXUlJSUaN++fVHnHz9+XAcOHFBJSUnc98zLy1NeXl7axgwAAJzD0TM74aCza9cu/fa3v9WECROinp89e7YOHjyo7du3R45t3LhRAwMDmjVrlt3DBQAADpTRmZ3Dhw/rr3/9a+Tvzs5O7dixQ4WFhZo0aZKqq6v1yiuvaMOGDQqFQpE6nMLCQo0ePVpnnnmmLrvsMi1atEgPP/ywjh07piVLluiqq65KeCUWAADwtowuPW9ra1N5eflJxxcuXKgVK1Zo2rRpMV8XDAY1Z84cSR81FVyyZIl+/vOfa9SoUZo7d67WrFmjcePGJTwOlp4DAOBdjumzk0mEHQAAvMvRNTsAAACpIuwAAABPI+wAAABPc3yfHQAAkJhQKKT29nbt3btXkyZNUllZmXJycjI9rIwj7AAA4AEtLS2qqanRO++8Eznm9/vV0NCgqqqqDI4s87iNBQCAy7W0tKi6ujoq6EhSd3e3qqur1dLSkqGROQNLz8XScwBA+qT71lIoFFJpaelJQSfM5/PJ7/ers7Mza29pMbMDAECatLS0qLS0VOXl5Zo/f77Ky8tVWlpq6UxLe3t73KAjScYYdXV1qb293bLPdBvCDgAAaWDXraW9e/daep4XEXYAALBYKBRSTU2NYlWKhI/V1tYqFAql/FmTJk2y9DwvIuwAAGAxO28tlZWVye/3y+fzxXze5/MpEAiorKws5c9yK8IOAAAWs/PWUk5OjhoaGiTppMAT/ru+vt7y4uRQKKS2tjY1NTWpra3NklmqdCHsAABgMbtvLVVVVWn9+vWaMmVK1HG/36/169db3mfHjsJrK7H0XCw9BwBYK7wcvLu7O2bdTrqWg9vRQTlceD34e4VnkdIRrlJF2BFhBwBgvXAokBQVDJwcCobj1p4+3MYCACAN7L61ZAe39vRhbywAANKkqqpKlZWVntmc0609fQg7AACkUU5OjubMmZPpYVjCrT19uI0FAAAS4taePoQdAACQkEz19EkVYQcAACTMjYXXLD0XS88BAEiWHT19rELYEWEHAAAv4zYWAADwNMIOAADwNMIOAADwNMIOAADwNMIOAADwNMIOAADwNMIOAADwNMIOAADwNMIOAADwtNxMDwAAAAzPTdszOA1hBwAAh2tpaVFNTY3eeeedyDG/36+GhgZHbrzpNOyNJfbGAgBYIx2zLy0tLaqurtbgn2ufzydJjt1p3EkIOyLsAABSl47Zl1AopNLS0qj3PJHP55Pf71dnZ2fcUMXtLwqUAQBIWXj2ZXAo6e7uVnV1tVpaWkb0vu3t7XGDjiQZY9TV1aX29va44yotLVV5ebnmz5+v8vJylZaWjng8bkXYAQAgBaFQSDU1NSfdZpIUOVZbW6tQKJT0e+/du3fE56UrgLkRYQcAgBSkOvsylEmTJo3ovHQGMDci7AAAkIJUZl+GU1ZWJr/fHylGHszn8ykQCKisrCzqeDoDmBsRdgAASMFIZ18SkZOTo4aGBkk6KfCE/66vrz+p4DidAcyNMhp2XnrpJV1++eWaPHmyfD6fnn322ajnjTG68847NWnSJI0dO1YVFRXatWtX1DkHDhzQggULNH78eBUUFOj666/X4cOHbfwWAIBsNtLZl0RVVVVp/fr1mjJlStRxv98fd9l5OgOYG2U07Lz33ns655xz9OCDD8Z8/kc/+pHWrFmjhx9+WJs3b9app56qSy+9VEeOHImcs2DBAv3pT3/SCy+8oA0bNuill17SDTfcYNdXAABkuZHOviSjqqpKu3fvVjAYVGNjo4LBoDo7O+MuaU93AHMd4xCSzDPPPBP5e2BgwJSUlJh77703cuzgwYMmLy/PNDU1GWOMeeONN4wks3Xr1sg5v/rVr4zP5zPd3d0Jf3ZfX5+RZPr6+lL/IgCArNTc3Gz8fr+RFHkEAgHT3NycsfH4fD7j8/mixhQ+lqlxZYJja3Y6OzvV09OjioqKyLH8/HzNmjVLHR0dkqSOjg4VFBRo5syZkXMqKio0atQobd68Oe57Hz16VP39/VEPAABSkezsy4lCoZDa2trU1NSktrY2S1ZJjeT2l1c5dm+snp4eSVJxcXHU8eLi4shzPT09mjhxYtTzubm5KiwsjJwTS11dnVauXGnxiAEA2S4nJ0dz5sxJ6jXp3PeqqqpKlZWVUR2UL7zwQm3atElNTU2Wd1R2ardmx87spNPy5cvV19cXeXR1dWV6SACALGRH479wAJs3b54OHDigT37yk2npqOzkbs2ODTslJSWSpN7e3qjjvb29kedKSkq0b9++qOePHz+uAwcORM6JJS8vT+PHj496AABgJ7sb/6UzWDm9W7Njw860adNUUlKi1tbWyLH+/n5t3rxZs2fPliTNnj1bBw8e1Pbt2yPnbNy4UQMDA5o1a5btYwYAIFF2Nv5LZ7ByQ7fmjIadw4cPa8eOHdqxY4ekj4qSd+zYobfffls+n0+1tbW6++679fzzz+u1117Ttddeq8mTJ+uKK66QJJ155pm67LLLtGjRIm3ZskW/+93vtGTJEl111VWaPHly5r4YAADDsLPxXzqDlRu6NWe0QHnbtm0qLy+P/L1s2TJJ0sKFC/XYY4/pe9/7nt577z3dcMMNOnjwoC6++GL9+te/1pgxYyKvWbt2rZYsWaJLLrlEo0aN0ty5c7VmzRrbvwsAAMmwovFfogXB6QxWbujW7DOx5p2yTH9/v/Lz89XX10f9DgDAFqFQSKWlperu7o55C8jn88nv96uzszNmgElmFVdbW1vU5EI8wWAw6dVk6Xxvqzi2ZgcAAC9LpfNysgXB6eyo7IZuzYQdAAAyZCSN/0ZSEJzOLS3s2C4jVdzGErexAACZlUwzvlRuG8W69RUIBFRfX59yA8N0vneqCDsi7AAA3KOpqUnz588f9rzGxkbNmzfvpOPp7HLs1A7Kjt0uAgAAnCzVVVwj2dIiUel871RQswMAgIu4oSDYaQg7AAC4iBsKgp2GsAMAgI1CoZDa2trU1NSktra2EW2jMJJVXNmMAmVRoAwAsEcyjQAT4dSCYKch7IiwAwBeFisQSLI9JIQbAQ7+2Q3femJGJn0IOyLsAIBXxZpJmTBhgiTp3XffjRxLZXYlEeGtIeJtmDnc1hBIDTU7AABPirelwrvvvhsVdKT42yxYxQ07g3sZYQcA4DlDbakQS7xtFqzihp3BvYywAwDwnOFmUmJJ5+xKqo0AkRrCDgDAc1KZIUnH7AqNADOLsAMA8JxUZkh6e3stv5XltkaAVvQCchLCDgDAc4abSRnK0qVLVVpaanmxcrKNADMVOFpaWlRaWqry8nLNnz9f5eXlabkedmLpuVh6DgBeFF6NJSnhQuWwdPa+SaQRoNXNBxPl1V5AhB0RdgDAqxLtsxNLIr1v0tHBOFOBw8u9gAg7IuwAgJfF66B8//33a+nSpcO+PhgMas6cOScdT8fsSyYDR1tbm8rLy4c9L971cLLcTA8AAIB0ysnJifnjXFxcnNDrY63Oijf7Em5OONLZl2SaD1odOLzcC4gCZQBAVkq29024YHjt2rW68cYbY9YBpdqcMJOBw8u9gAg7AICslEzvmxNXKF199dXav39/3PdNpTlhJgOHl3sBEXYAAFkp0d43zz33XMw9toYzktmXTAYOt/UCSgZhBwCQtYbrfVNZWZnUHlsnGsnsS6YDR7K9gNyC1VhiNRYAZLt4S8gTXaF0IitWTMVa6RUIBFRfX29L4EjHkvpMIuyIsAMAiK2pqUnz589P+Hwre+F4LXBkEkvPAQCII9lbUX6/37LZl3hL5pE8wg4AAHGEC4a7u7vj1u0UFRXpvvvu05QpU5h9cSgKlAEAiGO4gmGfz6eHH35YCxYs0Jw5cwg6DkXYAQBgCF5doZRNKFAWBcoAgOFle8Gwm78/YUeEHQDwEjf/KDtVOjY9tRO3sQAAnnHitg7z589XeXm5SktL1dLSkumhuVZ409PBHaTDm5664doysyNmdgDAC+LtRG5l75tsEwqFVFpaGnerDCsaKNqBmR0AgOuFQqG42zqkuhN5Nmtvbx9yT7BUNj21E2EHAOB6XvlRdppENzMdyaandiLsAABczys/yk6TaAfpkWx6aifCDgDA9bzyo+w04Q7Sgxsqhvl8PgUCAZWVldk8suQQdgAArueVH2WnGa6DtCTV19c7ujhZIuwAADzAKz/KTuSFDtKODjuhUEh33HGHpk2bprFjx+qTn/ykfvCDH0RV2xtjdOedd2rSpEkaO3asKioqtGvXrgyOGgCySygUUltbm5qamtTW1paxFU9e+FF2qqqqKu3evVvBYFCNjY0KBoPq7Ox0zTV1dJ+dH/7wh1q9erUef/xxffazn9W2bdt03XXXadWqVbr55pslSf/5n/+puro6Pf7445o2bZruuOMOvfbaa3rjjTc0ZsyYhD6HPjsAMDJO7KxLB2UM5uiw8/Wvf13FxcX66U9/Gjk2d+5cjR07Vk888YSMMZo8ebJuueUW/du//Zskqa+vT8XFxXrsscd01VVXxXzfo0eP6ujRo5G/+/v7FQgECDsAkASa+MEtHH0b68ILL1Rra6vefPNNSdKrr76ql19+WV/5ylckSZ2dnerp6VFFRUXkNfn5+Zo1a5Y6Ojrivm9dXZ3y8/Mjj0AgkN4vAgAeQxM/uElupgcwlNtuu039/f2aMWOGcnJyFAqFtGrVKi1YsECS1NPTI0kqLi6Oel1xcXHkuViWL1+uZcuWRf4Oz+wAgFek+1ZOMk385syZY9nnDsYtKyTC0WFn3bp1Wrt2rRobG/XZz35WO3bsUG1trSZPnqyFCxeO+H3z8vKUl5dn4UgBwDnsqKNxQhM/J9YLwZkcfRvru9/9rm677TZdddVVOvvss3XNNddo6dKlqqurkySVlJRIknp7e6Ne19vbG3kOALKJXTtUZ7qJnxd24oZ9HB123n//fY0aFT3EnJwcDQwMSJKmTZumkpIStba2Rp7v7+/X5s2bNXv2bFvHCgCZZmcdTSab+FEvZD2ntA9IF0eHncsvv1yrVq3SL37xC+3evVvPPPOMVq9erW9+85uSPvqXqba2Vnfffbeef/55vfbaa7r22ms1efJkXXHFFZkdPAAMw+ofGDs3w8xkEz82/bRWS0uLSktLVV5ervnz56u8vFylpaXemh0zDtbf329qamrM1KlTzZgxY8wnPvEJc/vtt5ujR49GzhkYGDB33HGHKS4uNnl5eeaSSy4xO3fuTOpz+vr6jCTT19dn9VcAgJiam5uN3+83kiIPv99vmpubR/yejY2NUe8X79HY2JjW7xEIBFL6HsPJxPf0qubmZuPz+U66dj6fz/h8vrT+72gnR/fZsQtNBQHYKV39adra2lReXj7secFg0NIVUnaviMrU9/SaUCik0tLSuLNkPp9Pfr9fnZ2drl/hRtgRYQeAfdL5AxN+7+7u7pj1LF758cqW75lu2RQaHV2zAwBek856k2zZDNOu7+n1ol0ntA+wC2EHAGyU7h+YbNkMM93fMxuKdjPdPsBO3MYSt7EA2MeuWwfZ0lk4Hd8zW/b8yqbbgYQdEXYA2CebfmDcKJuKdqV/BDtJUf88ei3YcRsLAGyULXU1bpVtPXyy5bYnYQcAbJYtPzBulE1Fu2FVVVXavXu3gsGgGhsbFQwG1dnZ6al/Dh29ESgAeFVVVZUqKyuzoq4mUU6oM8qmot0T5eTkuH55+VCo2RE1OwAgZTZsOGUHc2qqvInbWACAjC61dtIO5tRUeRMzO2JmB0B2G2qptTFGK1eu1PTp09My2+PU1U+xZpoCgYDq6+s9VcuSLQg7IuwAcLdUbj8NFzYGs/rWkpO3LHBCDRGsQYEyALhYqrUuwy21Hix8a8mqVWNOXv3k9aLdbELNDgC4lBW1LsmGiPDNgNraWkv2isrW1U+wF2EHAFwoFAqppqYm5oqhZALJSEKElY31ysrK5Pf7TyoGDvP5fAoEAiorK0v5s5C9CDsA4EJWdfodLmwMxYpbS8mufvL6TuRID8IOALhQsrUu8ULCUGFjOKncWjpxPIWFhVq3bt2wHaWzYSdypAerscRqLADuk8wqpgMHDgxbxByr0DmeVJeDxyuqXr16tYqKimKufsqWnciRHoQdEXYApM7uZcqJdvpdvXq1rrzyyoRCwonfYdeuXVqxYoWk5HfDHupajCS0OLUXD1zEwPT19RlJpq+vL9NDARDD8ePHTTAYNI2NjSYYDJrjx4876vObm5uN3+83kiIPv99vmpub0zqu5uZm4/P5jM/ni/rs8LF169adNK7B5wUCgbjXM9b3CgQCQ36voa7F8ePHRzSeYDAY9zUnPoLBoJWXFx5C2DGEHcDJMhUkEv38cOCI9cPt8/lsCTzxAokVISGZoDnctVi5cuWIxtPY2JjQ6xobGy26qvAawo4h7ABO5YQgMdTnpzpzYpV4gcTOkJDIrE1hYeGIxsPMDlJFzY6o2QGcKNN1Gol8/umnn679+/cP+17hrQ7sruuxcyuGRD8rEYPHw07kSBVLzwE4klV9ZNL5+YkEHemj5d+ZWDZtZ8O+RJfCFxYWJj0ediJHqgg7ABwp03smWfm+u3btSnlbh5GwMyQk2nOnpqZmROOpqqrS+vXrh+3FA8RC2AHgSMnumWR1Z91EP7+oqGjImQq/369HHnkk5W0dRsqukJDoLNLtt98+4vFUVVVp9+7dCgaDamxsVDAYVGdnJ0EHw8tgvZBjUKAMOE+44DVWgbAGFf+mY8VWop//9NNPD7n8e6QrkKxmx/L94ZbCn/i/R6bbCSC7EHYMYQdwqkR+PNO5YivRH++hln9n27LpkfTmAdKN1VhiNRbgZLG2FggEAqqvr1dlZWXaV2wN9fkn3j6Jt9LKzhVRdkhkRZndq86A4RB2RNgBnC7TQSKVH28vLZuOt6fViXtsxUMAQiblZnoAADCcnJycmGHFrhVb8T4/0dc2NDSourpaPp8v5j5Tblg2HW9Pq/CKsqGKi0cakghIsAqrsQC4VrIrtjLF7cumQ6GQampqRrSiLBySkl12n4m+RPAubmOJ21iAW7ntFpFbZypGertwpF2wR7IzOjAUZnYAuJbbOuuGb4fNmzdPc+bMccy4hjPS24Uj6YKdyiwSEA9hB4Cruf0WkRuM9HbhSEJSprcJgTdRoAzA9aqqqlRZWXnSLSLpo1swbrttdCIn3PoKd0ce7nbh4D2tRhKSMr1NCLyJsAPAEwavmEplmXQmnRhudu3apUceeSTj32GkK8pGEpLcUnQOl7G3h6Ez0UEZ8JZ4XZXDj9raWkduURCr+/DghxWdoa0c33DdkZPZQsKY5LYJARLFaiyxGgvwkuFWAJ3ISTM98VYgxZLJVWYjua2WaBfqE8+vrq6WpJizSNRiIVlJh52FCxfq+uuv15e+9KV0jcl2hB3AOxJdJi2N/MfT6jqaZALaidyyxYSU/DVLNiABQ0m6Zqevr08VFRU644wzdN1112nhwoUnrYIAnMQJBZ6wTzKFq8YY+Xw+1dbWqrKyMqF/LtJRCzTcCqR43FSkm2wX6nhF5/y7i5FIeun5s88+q+7ubt1000166qmnVFpaqq985Stav369jh07ZvkAu7u7dfXVV2vChAkaO3aszj77bG3bti3yvDFGd955pyZNmqSxY8eqoqJCu3btsnwccCe6sGafZAtXTRJLmUfaDXg4Iw0tXi/SdWtfIjjPiPrsFBUVadmyZXr11Ve1efNmfepTn9I111yjyZMna+nSpZaFjf/7v//TRRddpFNOOUW/+tWv9MYbb+i//uu/9PGPfzxyzo9+9COtWbNGDz/8sDZv3qxTTz1Vl156qY4cOWLJGOBe6fphgrOFVwANbjI4nOECRzqb3SUbWnw+nwKBwElLvQHEkUp18549e8w999xjPv3pT5tTTz3VXHvtteaSSy4xubm5ZvXq1SlVThtjzK233mouvvjiuM8PDAyYkpISc++990aOHTx40OTl5Zmmpqa4rzty5Ijp6+uLPLq6uliN5THhFR0aYkULKzq8K94KoKEewWBwyPcMBoOWvE8sw61AGvzPbqZWYwFulfTMzrFjx9Tc3Kyvf/3rOuOMM/T000+rtrZWe/bs0eOPP67f/va3Wrdunf7jP/4j+eQ1yPPPP6+ZM2fqW9/6liZOnKhzzz1XjzzySOT5zs5O9fT0qKKiInIsPz9fs2bNUkdHR9z3raurU35+fuQRCARSHiuchS6s2S1eV+V4ioqK1N3drba2trgzM+lsdjfUtheD0RkaSF7SYWfSpElatGiRzjjjDG3ZskXbtm3TjTfeGLWKqby8XAUFBSkP7m9/+5seeughTZ8+Xb/5zW9000036eabb9bjjz8uSerp6ZEkFRcXR72uuLg48lwsy5cvV19fX+TR1dWV8ljhLHRhRVVVlXbv3q1gMKja2lpJ8YPE/v37dfXVVw9Z05XuZndDbXuxcuVKNTY2KhgMqrOzk6ADJCnppec/+9nP9K1vfUtjxoxJ15giRo8erZkzZ2rTpk2RYzfffLO2bt2qjo4Obdq0SRdddJH27NkT9R+YK6+8Uj6fT0899VRCn8PSc+8Z6S7N8K5Yq6hiibcc3a4d1hNZPcgKQyA5Sc/sXHPNNbYEHemj/w/pM5/5TNSxM888U2+//bYkqaSkRJLU29sbdU5vb2/kOWSn4YpUKfDMPifO9DzxxBMqKiqKeZ6JU2xs1w7rw61AYoUhkDxH73p+0UUXaefOnVHH3nzzTZ1xxhmSpGnTpqmkpEStra2R5/v7+7V582bNnj3b1rHCWez6YYK7hIPElClTtH///rjnxavpyvQO66wwBEYok9XRw9myZYvJzc01q1atMrt27TJr1641H/vYx8wTTzwROeeee+4xBQUF5rnnnjN//OMfTWVlpZk2bZr54IMPEv4c9sbyrpHs5QPva2xsTGhlVWNjY8zXHz9+3ASDQdPY2GjbHlusMARGzvF7Y23YsEHLly/Xrl27NG3aNC1btkyLFi2KPG+M0V133aWf/OQnOnjwoC6++GL9+Mc/1j/90z8l/BnU7Hgb9Q0YzI01XW4cM+AUjg87diDsANnFrmJjKzU1NWn+/PnDntfY2Kh58+bZMCLAPRxdswMA6eDGmq50L30HvIywAyArZbrYOFmsMARGjttY4jYWkM3SVdOVjvcNr8aSFHX7LV5voEyMEXAiwo4IOwCsFauBod/vV0NDQ8ozRrHeOxAIqL6+Pqn3TucYAach7IiwA6RLNs4chGdfBv+ndaSzL7Gkel3tGCPgJIQdEXaAdMjGmYPwKq94W1I4YZWXG8YIWI0CZQCWy9ZOv+3t7UPuvWXidGa2kxvGCFiNsAPAUqFQSDU1NTH715g4+055xd69ey09Lx3cMEbAaoQdAJZy+8xBKBRSW1ubmpqa1NbWllQoc0MvHDeMEbAaYQeApdw8c5DqjuJu6IXjhjECViPsALCUW2cOrKgzckNnZjeMEbAaYQeApdw4c2BlnZEbOjO7YYyAlVh6LpaeA4lIpreL1Z1+Ux3PcNKxo7gbegy5YYyAFXIzPQAAzpdsz5zwzEGs1yTb6deK8QwnHXVGOTk5CQejTHHDGAErMLMjZnaAoaTSbTed+0NZ2f03HTM7AJyDsCPCDhCP07rtpms84fft7u6OWbdDV2HA3ShQBhCXHT1zkulrk67xsEIJ8DZqdgCbuakoNN09c+LV3qxevVpFRUUnXaN0jifddUYAMoewA9jIbZtjprNnTrzam3feeUdXXnll1LHwNUp3D5+qqipVVla6JowCSAw1O6JmB/ZIR2FtuqWrlmW42ptYnyNJTz31lJYtW0ZtDYCkULMD2MCtm2Omq5ZluNqbwcLX6JZbbtF9991n+XgAeBthB7CBmzfHTEe33ZHU1ISv0emnn073XwBJoWYHsIGbN8eUrK9lSWVfrL1792revHnU1gBIGGEHsIFbN8c8kZXddsP7Z8WrvRlK+BrR/RdAoriNBdjAjZtjptNQtUDxZNs1AmAdwg5gA5rWnSxeLVAs2XqNAFiDpedi6TnsE6vPTiAQyOqmdYObLP7973/X0qVLuUYALEPYEWEH9nJTB+VM4RoBsBJhR4QdAAC8jJodAADgaYQdAADgaYQdAADgaYQdAADgaYQdAADgaYQdAADgaYQdAADgaYQdAADgaex6DrgEXYUBYGQIO4ALxNpTy+/3q6Ghgf2iAGAY3MYCHK6lpUXV1dVRQUeSuru7VV1drZaWlgyNDADcwVVh55577pHP51NtbW3k2JEjR7R48WJNmDBB48aN09y5c9Xb25u5QQIWCoVCqqmpUawt7MLHamtrFQqF7B4aALiGa8LO1q1b9d///d/63Oc+F3V86dKl+vnPf66nn35aL774ovbs2cO0Pjyjvb39pBmdExlj1NXVpfb2dhtHBQDu4oqwc/jwYS1YsECPPPKIPv7xj0eO9/X16ac//alWr16tf/7nf9Z5552nRx99VJs2bdLvf//7DI4YsMbevXstPQ8AspErws7ixYv1ta99TRUVFVHHt2/frmPHjkUdnzFjhqZOnaqOjo6473f06FH19/dHPQAnmjRpkqXnWSkUCqmtrU1NTU1qa2vjVhoAx3J82HnyySf1yiuvqK6u7qTnenp6NHr0aBUUFEQdLy4uVk9PT9z3rKurU35+fuQRCASsHjZgibKyMvn9fvl8vpjP+3w+BQIBlZWV2TqulpYWlZaWqry8XPPnz1d5eblKS0splgbgSI4OO11dXaqpqdHatWs1ZswYy953+fLl6uvrizy6urose2/ASjk5OWpoaJCkkwJP+O/6+npb++2kujqMGSEAdnN02Nm+fbv27dunL3zhC8rNzVVubq5efPFFrVmzRrm5uSouLtaHH36ogwcPRr2ut7dXJSUlcd83Ly9P48ePj3oAmTRUAKiqqtL69es1ZcqUqNf4/X6tX7/e1oL8VFeHMSMEIBN8JtZ/tRzi0KFDeuutt6KOXXfddZoxY4ZuvfVWBQIBFRUVqampSXPnzpUk7dy5UzNmzFBHR4e++MUvJvQ5/f39ys/PV19fH8EHtku0YaATOii3tbWpvLx82POCwaDmzJkTdSw8IzT4PznhGSq7gxuA7OHoDsqnnXaazjrrrKhjp556qiZMmBA5fv3112vZsmUqLCzU+PHj9Z3vfEezZ89OOOgAmRQvAIRvCZ0YAHJyck4KEHZLdNVXc3OzJEUC2XAzQuH+WZWVlWyBAcByjr6NlYj77rtPX//61zV37lx96UtfUklJCVPicAUnNwyMd1st0VVfDzzwQNQtKvoFAcgkR9/Gsgu3sZAJqdwSSqehbqtVVlaqtLRU3d3dMUPaYOFbVDU1Naqvrx/2/MbGRs2bN2/EYweAWFw/swO4lRMbBg630uq5556LuzoslnAgWrt2bUKfn4l+QQC8j7ADZIgVDQOtXMad6G21ysrKmKvD4jHGaP/+/SoqKnJcvyAA2YGwA2RIqg0DrV7GnUxdTVVVlXbv3q1gMKglS5Yk9P4LFiyQ5Jx+QQCyB2EHSEEiMyvxzkmlYWCqjf1iSfa2Wnh1WLjtw3DizQhlol8QgCxjYPr6+owk09fXl+mhwEWam5uN3+83kiIPv99vmpubUz4nEAhEnXOi48ePn3T+iQ+fz2cCgYA5fvx4Ut8nGAzGfc8TH8FgMOZ4fD5fQuM5fvy4CQaDprGx0QSDwaTHCQDJYjWWWI2VjVJt0JdIgzxJCTfRS2Y8qa7iivdZoVBoyJVWPp9Pfr9fnZ2dJ40tfD0kRb2WhoEAHCGTScspmNnJLonMtgwlkZkVv9+fltkXY4xpbGxMaAamsbEx6e/e3NxsfD7fSbM04WNDXaNkZ6gAwC6EHUPYySbhH/NY4WO4H/OwRG/3JPIYfEsoESO93ZTod08ltHCLCoATcRtL3MbKFuHbNPFWHA11m+ZETU1Nmj9/viVjGkkTvZHcbkr2uzthHy4AsAqrsZA1rNqywMrGdyN5r2RXcYVCId1///1JfffwSqt58+Zpzpw5BB0ArkbYQdZItWNxeAl5d3f3sA3y/H5/Sj10hlNVVZXQMu5wL56lS5cm9L52dmsGALs4etdzwEqpdCyOtV9ULOFwE555qa6uls/ni7lCKdUmelVVVaqsrIx7uyneirGhsF0DAC+iZkfU7GSLkS6tTiY0BAIB1dfXR82sDA5Jg89Jh+FqdAZLtF4JANyIsCPCTjZJth9MIqGhqKhI9913n6ZMmRKzkNfOYt/wZ7W2turuu+9O6DXhmaeVK1dq+vTpFCQD8BzCjgg72SaZ2ZZUG/jZKdFbbYNNmDBBkvTuu+9Gjvn9fjU0NNAIEIAnUKCMrHPiJpaNjY0KBoPq7OyM+cOealGzXeLtlTWc6667TgcOHIgKOlJqe2wBgNMwsyNmdhCfG2Z2kq3PkT66dRVeyZVq3yEAcDpmdoAhlJWVpXUJuRWG6x80WPi7LFq0yJK+QwDgdIQduFq4901TU5Pa2toUCoUsff9kG/hlQrK30MK9eKZPn56W9wcApyHswLXCDfPKy8s1f/58lZeXq7S01PI6k0Qb+GVKor1x/v3f/z2qPimVvkMA4CbU7IiaHTeK1/sm3hJyKzh1v6iR9g8a6esAwG0IOyLsuI1VG3p6SbL9g1J9HQC4Cbex4DqJbui5YsWKtNTxONFIb7U5/RYdAFiBmR0xs+MmoVBIK1asSLg7sJRdDfJGeqvNqbfoAMAKhB0RdtxipB2CE7klw489AHgXYUeEHTcYyQ7eJxqqjidWiMqm2SAA8DrCjgg7TjeSDsHxDO50PNSqLjbHBABvyM30AIDhJNsheCgnNsgLhUKqqamJOVsUPnbXXXdFjjHbAwDuxGqsDEp391+vsLKD74kN8pINUWyOCQDuRNjJELu6/3pBoh18i4qKktrDKtkQFZ7tqa2tJZgCgIsQdjIgXCcyeFaBmYPYEt2M88c//nHk78HPSyfvYTWSbRDYHBMA3IewY7NE6kSYOYiW6Gac1dXVSTXIGy5EDSXdm2NyixMArEPYsVmi3X+ZOYiWaKffqqoq7d69W8FgUI2NjVEbXw42VIgaTjo3x+QWJwBYi6XnsnfpeVNTk+bPnz/seY2NjZo3b15ax+JG6Wj+l0yzwnTvu5WJDU4BwOsIO7I37LS1tam8vHzY8wb3g0F6nRiidu3apRUrVkiyd3NMNjgFgPQg7MjesBP+Qevu7o5Zt8MPmjPEmu0JBAKqr69P28wKQRgA0oOmgjYL14lUV1dHuvSGxVs1BPtVVVWpsrLS1v2yEi16TndxNAB4DWEnA8LFtrH2Y0rnzAGSk5OTY+sMSqJFz+ksjgYAL+I2ljK3NxY7beNE3OIEgPRgZieD7J45gLNxixMA0oM+O4CDJNpPCACQOEeHnbq6Op1//vk67bTTNHHiRF1xxRXauXNn1DlHjhzR4sWLNWHCBI0bN05z585Vb29vhkYMpC6ZxogAgOE5umbnsssu01VXXaXzzz9fx48f1/e//329/vrreuONN3TqqadKkm666Sb94he/0GOPPab8/HwtWbJEo0aN0u9+97uEPydTNTsYGjVNAAArODrsDLZ//35NnDhRL774or70pS+pr69PRUVFamxsVHV1tSTpL3/5i84880x1dHToi1/8YkLvS9hxnlh9bvx+vxoaGpjhAAAkxdG3sQbr6+uTJBUWFkqStm/frmPHjqmioiJyzowZMzR16lR1dHTEfZ+jR4+qv78/6uFU2bghJLvCAwCs5JqwMzAwoNraWl100UU666yzJEk9PT0aPXq0CgoKos4tLi5WT09P3Peqq6tTfn5+5BEIBNI59BHLxg0h2RUeAGA114SdxYsX6/XXX9eTTz6Z8nstX75cfX19kUdXV5cFI7RWts5usCs8AMBqruizs2TJEm3YsEEvvfSS/H5/5HhJSYk+/PBDHTx4MGp2p7e3VyUlJXHfLy8vT3l5eekcckqGm93w+Xyqra1VZWWlpQW7TigIZssEAIDVHD2zY4zRkiVL9Mwzz2jjxo2aNm1a1PPnnXeeTjnlFLW2tkaO7dy5U2+//bZmz55t93Atk4nZDafcMmPLBACA1Rw9s7N48WI1Njbqueee02mnnRapw8nPz9fYsWOVn5+v66+/XsuWLVNhYaHGjx+v73znO5o9e3bCK7GcyO7ZjfAts8EzSeFbZnY2sysrK5Pf7x92y4SysjJbxgMAcD9HLz0Pt8gf7NFHH9W3v/1tSR81FbzlllvU1NSko0eP6tJLL9WPf/zjIW9jDea0pedtbW0qLy8f9rxgMJjydhPh/ZjizSQlsh/T4NtfF154oTZt2jTi22Hh8CUp5pYJdBIGACTD0WHHLk4LO3ZtCBkKhXT//fdr6dKlw54bL1jF6oeTk5MTtVpqJP1xYr1vIBBgV3gAQNIIO3Je2JGGnt0wxmjlypWaPn36iAuJY4WJoTQ2NmrevHkxxzjcP0IjnZFxQsE0AMD9CDtyZtiRYgeSCRMmSJLefffdyLFkZ04SDSknGjyzM9ztr8Gsmo0CACBZhB05N+xI0bMbu3bt0ooVK04KKcnMnCQbUiSpqKhI9913n6ZMmRKZXUm0rmgwK+qMAABIhqNXY+Gj+pc5c+ZEQkqqvXeGW9Yey/79+3X11VdL+scs0tGjR5P7Iv8f/XEAAHZzdJ8d/INVvXdSDRvh5ei7du0a0evpjwMAsBthxyWs6r2TaNjIz8+PeTw8s/TII4/I7/fHbQ8wmM/nUyAQoD8OAMB2hB2XsKqzcLhpX7yQ4vP5VFRUFNlhPhZjjN555x0tWrQo8pqhhJ+vr6+nOBkAYDvCjkskElKGmjkJhUJqa2vTunXr4oaU8N8LFixIaEzTp0/X+vXrNWXKlKjjgwON3++nESAAIGNYjSVnr8Y60Ug7Cye6hD3ctK+wsDCpDs5Wd1AGAMBKhB25J+xIyXcWjtdTZ6jmhHZ0cKZhIADALoQduSvsSIkHhVT2vUrn/lSxAttItpQAACARhB25L+wkKtUNRdOxP9VQM02S+zf5ZMYKAJyHpoIeEO8HNtXl6lVVVaqsrLTsxzsUCqmmpiblxohOxYwVADgTYcflhvqBtWK5eriDsxWSaYzoti0l4s1YhZswun3GCgDcjKXnLhb+gR0cIMI/sPv3709pubrVEp1pam5uVltbm0KhUJpHZI3hZqwkqba21jXfBwC8hrDjUon8wN5yyy267777JMXvqWNno79EZ5oeeOABlZeXq7S0VC0tLWkeVeqs2soDAJAehB2XSvQH9tVXX9WKFStOavyXiUZ/wzVGHCw8Q+X0wGPVVh4AgPSgZselEv3hvPvuuyVJU6ZMidlTx045OTlqaGhQdXV1pM/PUNxStGzVVh4AgPRgZselkv3h3LNnj1asWKG8vDzNmTMnY8Ghqqoq5hYT8bjhFlCqW3kAANKLsONSyd4SclKhbFVVlXbv3q1gMKglS5Yk9Bon3wIKz1hJzqiNAgBEI+y41FA/sPFkcpYkvBFpU1OT2traJElz5szR3LlzE3q9028BxZuxYhNUAMg8OijL3R2UY/XZGU5jY6PmzZuXxlFFG6oXUGVlZdr34bITHZQBwHkIO3J32JH+8QPb2toaKUgeSrztIdIhke0hJKVtHy4AAAg7cn/YCbNjt/KRjCeRjUife+45y/fhAgBAIuxIsi/s2HGLI527lScr2Y1IuQUEAEgH+uzYxK5NIsOFsrE+y+5ZkmSb7Vm5DxcAAGHM7Cj9MzuJ1K1YHUKcMEuS7MwOAADpQNhResNOMnUrXrtl47QaIgBAdqLPTppl8yaRNNsDADgBYSfNsn2TSJrtAQAyjQLlNGOTyI8CT2VlZcZriAAA2YmaHdlTs0PdCgAAmcFtrDSjbgUAgMwi7KTJiRtfFhYWat26ddStAACQAdTspEG8BoKrV69WUVERdSsAANiImh1ZW7OTiQaCAAAgPsKOrAs72dxAEAAAp6Jmx0LZ3EAQAACnIuxYKNsbCAIA4ESEHQvRQBAAAOch7FiorKxMfr//pH46YT6fT4FAQGVlZTaPDACA7OWZsPPggw+qtLRUY8aM0axZs7Rlyxbbx0ADQQAAnMcTYeepp57SsmXLdNddd+mVV17ROeeco0svvVT79u2zfSxsfAkAgLN4Yun5rFmzdP755+uBBx6QJA0MDCgQCOg73/mObrvttpPOP3r0qI4ePRr5u7+/X4FAwNK9sUKhEBtfAgDgAK7voPzhhx9q+/btWr58eeTYqFGjVFFRoY6Ojpivqaur08qVK9M6rpycHM2ZMyetnwEAAIbn+ttYf//73xUKhVRcXBx1vLi4WD09PTFfs3z5cvX19UUeXV1ddgwVAABkgOtndkYiLy9PeXl5mR4GAACwgevDzumnn66cnBz19vZGHe/t7VVJSUmGRjVy1PoAAGAt19/GGj16tM477zy1trZGjg0MDKi1tVWzZ8/O4MiS19LSotLSUpWXl2v+/PkqLy9XaWmpWlpaMj00AABcy/UzO5K0bNkyLVy4UDNnztQFF1yg+vp6vffee7ruuusyPbSExdstvbu7W9XV1Y5Zts7MEwDAbTyx9FySHnjgAd17773q6enR5z//ea1Zs0azZs1K6LVW7Xo+Um7ZLb2lpUU1NTVR4/T7/WpoaHBEEAMAIBbPhJ1UZDrstLW1qby8fNjzgsFgxpazx5t5CneGdsrMEwAAg7m+ZscLnL5beigUUk1NzUlBR1LkWG1trUKhkN1DAwBgWIQdB3D6bunt7e1xb7FJHwWerq4utbe32zgqAAASQ9hxAKfvlu70mScAAIZC2HEAp++W7vSZJwAAhkLYcQgn75bu9JknAACGwmosZX411omc2scmvBpLUlShMquxAABOR9iRs8KOk8XqsxMIBFRfX0/QAQA4FmFHhJ1kOHXmCQCAeAg7IuwAAOBlFCgDAABPI+wAAABPI+wAAABPI+wAAABPI+wAAABPI+wAAABPI+wAAABPI+wAAABPI+wAAABPy830AJwg3ES6v78/wyMBAADJOu200yIbU8dC2JF06NAhSR9tagkAANxluO2e2BtL0sDAgPbs2TNsMkxWf3+/AoGAurq62HMrjbjO9uFa24PrbA+usz3suM7M7CRg1KhR8vv9aXv/8ePH8y+SDbjO9uFa24PrbA+usz0yeZ0pUAYAAJ5G2AEAAJ5G2EmjvLw83XXXXcrLy8v0UDyN62wfrrU9uM724DrbwwnXmQJlAADgaczsAAAATyPsAAAATyPsAAAATyPsAAAATyPspNGDDz6o0tJSjRkzRrNmzdKWLVsyPSRXq6ur0/nnn6/TTjtNEydO1BVXXKGdO3dGnXPkyBEtXrxYEyZM0Lhx4zR37lz19vZmaMTud88998jn86m2tjZyjGtsne7ubl199dWaMGGCxo4dq7PPPlvbtm2LPG+M0Z133qlJkyZp7Nixqqio0K5duzI4YvcJhUK64447NG3aNI0dO1af/OQn9YMf/EAnrs3hOifvpZde0uWXX67JkyfL5/Pp2WefjXo+kWt64MABLViwQOPHj1dBQYGuv/56HT58OD0DNkiLJ5980owePdr8z//8j/nTn/5kFi1aZAoKCkxvb2+mh+Zal156qXn00UfN66+/bnbs2GG++tWvmqlTp5rDhw9HzrnxxhtNIBAwra2tZtu2beaLX/yiufDCCzM4avfasmWLKS0tNZ/73OdMTU1N5DjX2BoHDhwwZ5xxhvn2t79tNm/ebP72t7+Z3/zmN+avf/1r5Jx77rnH5Ofnm2effda8+uqr5hvf+IaZNm2a+eCDDzI4cndZtWqVmTBhgtmwYYPp7Ow0Tz/9tBk3bpxpaGiInMN1Tt4vf/lLc/vtt5uWlhYjyTzzzDNRzydyTS+77DJzzjnnmN///vemvb3dfOpTnzLz5s1Ly3gJO2lywQUXmMWLF0f+DoVCZvLkyaauri6Do/KWffv2GUnmxRdfNMYYc/DgQXPKKaeYp59+OnLOn//8ZyPJdHR0ZGqYrnTo0CEzffp088ILL5gvf/nLkbDDNbbOrbfeai6++OK4zw8MDJiSkhJz7733Ro4dPHjQ5OXlmaamJjuG6Alf+9rXzL/8y79EHauqqjILFiwwxnCdrTA47CRyTd944w0jyWzdujVyzq9+9Svj8/lMd3e35WPkNlYafPjhh9q+fbsqKioix0aNGqWKigp1dHRkcGTe0tfXJ0kqLCyUJG3fvl3Hjh2Luu4zZszQ1KlTue5JWrx4sb72ta9FXUuJa2yl559/XjNnztS3vvUtTZw4Ueeee64eeeSRyPOdnZ3q6emJutb5+fmaNWsW1zoJF154oVpbW/Xmm29Kkl599VW9/PLL+spXviKJ65wOiVzTjo4OFRQUaObMmZFzKioqNGrUKG3evNnyMbERaBr8/e9/VygUUnFxcdTx4uJi/eUvf8nQqLxlYGBAtbW1uuiii3TWWWdJknp6ejR69GgVFBREnVtcXKyenp4MjNKdnnzySb3yyivaunXrSc9xja3zt7/9TQ899JCWLVum73//+9q6datuvvlmjR49WgsXLoxcz1j/HeFaJ+62225Tf3+/ZsyYoZycHIVCIa1atUoLFiyQJK5zGiRyTXt6ejRx4sSo53Nzc1VYWJiW607YgSstXrxYr7/+ul5++eVMD8VTurq6VFNToxdeeEFjxozJ9HA8bWBgQDNnztQPf/hDSdK5556r119/XQ8//LAWLlyY4dF5x7p167R27Vo1Njbqs5/9rHbs2KHa2lpNnjyZ65xFuI2VBqeffrpycnJOWqHS29urkpKSDI3KO5YsWaINGzYoGAzK7/dHjpeUlOjDDz/UwYMHo87nuidu+/bt2rdvn77whS8oNzdXubm5evHFF7VmzRrl5uaquLiYa2yRSZMm6TOf+UzUsTPPPFNvv/22JEWuJ/8dSc13v/td3Xbbbbrqqqt09tln65prrtHSpUtVV1cnieucDolc05KSEu3bty/q+ePHj+vAgQNpue6EnTQYPXq0zjvvPLW2tkaODQwMqLW1VbNnz87gyNzNGKMlS5bomWee0caNGzVt2rSo58877zydcsopUdd9586devvtt7nuCbrkkkv02muvaceOHZHHzJkztWDBgsj/zTW2xkUXXXRS64Q333xTZ5xxhiRp2rRpKikpibrW/f392rx5M9c6Ce+//75GjYr+qcvJydHAwIAkrnM6JHJNZ8+erYMHD2r79u2RczZu3KiBgQHNmjXL+kFZXvIMY8xHS8/z8vLMY489Zt544w1zww03mIKCAtPT05PpobnWTTfdZPLz801bW5vZu3dv5PH+++9HzrnxxhvN1KlTzcaNG822bdvM7NmzzezZszM4avc7cTWWMVxjq2zZssXk5uaaVatWmV27dpm1a9eaj33sY+aJJ56InHPPPfeYgoIC89xzz5k//vGPprKykiXRSVq4cKGZMmVKZOl5S0uLOf300833vve9yDlc5+QdOnTI/OEPfzB/+MMfjCSzevVq84c//MG89dZbxpjErulll11mzj33XLN582bz8ssvm+nTp7P03I3uv/9+M3XqVDN69GhzwQUXmN///veZHpKrSYr5ePTRRyPnfPDBB+Zf//Vfzcc//nHzsY99zHzzm980e/fuzdygPWBw2OEaW+fnP/+5Oeuss0xeXp6ZMWOG+clPfhL1/MDAgLnjjjtMcXGxycvLM5dcconZuXNnhkbrTv39/aampsZMnTrVjBkzxnziE58wt99+uzl69GjkHK5z8oLBYMz/Hi9cuNAYk9g1fffdd828efPMuHHjzPjx4811111nDh06lJbx+ow5oY0kAACAx1CzAwAAPI2wAwAAPI2wAwAAPI2wAwAAPI2wAwAAPI2wAwAAPI2wAwAAPI2wAwAAPI2wAwAAPI2wAwAAPI2wAwAAPI2wA8Bz9u/fr5KSEv3whz+MHNu0aZNGjx6t1tbWDI4MQCawESgAT/rlL3+pK664Qps2bdKnP/1pff7zn1dlZaVWr16d6aEBsBlhB4BnLV68WL/97W81c+ZMvfbaa9q6davy8vIyPSwANiPsAPCsDz74QGeddZa6urq0fft2nX322ZkeEoAMoGYHgGf97//+r/bs2aOBgQHt3r0708MBkCHM7ADwpA8//FAXXHCBPv/5z+vTn/606uvr9dprr2nixImZHhoAmxF2AHjSd7/7Xa1fv16vvvqqxo0bpy9/+cvKz8/Xhg0bMj00ADbjNhYAz2lra1N9fb1+9rOfafz48Ro1apR+9rOfqb29XQ899FCmhwfAZszsAAAAT2NmBwAAeBphBwAAeBphBwAAeBphBwAAeBphBwAAeBphBwAAeBphBwAAeBphBwAAeBphBwAAeBphBwAAeBphBwAAeNr/A0YK3s+AH/S3AAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "# generate a true dataset\n", + "theta_true = [1, 5]\n", + "y_true = simulator(theta_true)\n", + "\n", + "# and visualize it\n", + "plt.clf()\n", + "plt.scatter(np.linspace(0, 100, 101),\n", + " np.array(y_true), color = 'black')\n", + "plt.xlabel('x')\n", + "plt.ylabel('y')\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "id": "3b68d4b4-91e2-4120-8cc3-8bc1eafbf883", + "metadata": {}, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "ab092726a226474a80e7032465acf102", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "Drawing 10000 posterior samples: 0%| | 0/10000 [00:00]" + ] + }, + "execution_count": 19, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAHRCAYAAACmZ/R8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAogElEQVR4nO3df4xb1d3n8a/t8cx4ZjLjSWYySYb8ggAtDUkQkBAWRHiULWq7WaFqq4h2IeJpi2iXiDSpVqgFglS13We1oFQtLVJbBO2qVfpIbVVt22g3eQqFKpCSNO1DoZSQkJBJMj8Szy9nEs/Y3j9S+55zZnzHHh/7+vq+XxKSPdc/7jjEn9zzPed7QtlsNisAAFgU9voEAAD1h3ABAFhHuAAArCNcAADWES4AAOsIFwCAdYQLAMA6wgUAYB3hAgCwrsHrEwCC5j+GP1Xx9whFG/O3s5Mp669p83XhP/8v86+zPoYrFwCAdYQLAMA6hsWAOlSJIat6GQarxJAhpuPKBQBgHeECALCOYTHAp2wM7wRxBlgQfsdaQLjUiL7hCUkknf/pO1sbpTce8/CMAGDuCJca0Dc8IZufflkmJtP5n8WiEdm/6y4CBoAvES41IJFMycRkWvZsXSerFrbJsYFx2bH3qCSSKcIFgC8RLjVk1cI2Wd3b4fVpwCcqUTswazCVfj+39zffjynE/kK4VBF1FQBBQbhUiVtdBQDqDeFSJW51FaCawq2t+duZZFI75vXQk9t71tJQWBCncJeKcKky6ioAgoAV+gAA67hysYyiPQAQLlZRtIdtNmogao1lNm7vYb6OWa9Redmaphr1EGossyNcLKJoDwBXEC4VQNEeQNBR0AcAWMeVC1DD1LH9udYSiq2NiIiEGqMFn1fKmhi38y70uNle08axmY6jMrhyAQBYR7gAAKxjWMxjxwbGvT4F1IFIZ6d2PzPuDGGVMgyUTU3mb7sNmU173hyHqUoZsirl/ed6boUeN9tjMR3h4pHO1kaJRSOyY+9REbmyHqaztfDYNAD4CeHikd54TPbvuiu/BoaV/ADqCeHiod54jEABUJcIF8An3Mb81RqL+Vi32oFZxwi1OS1esuOFpzCL6PUZt3MtZSqyWjtKJxIF36+UWsls71ns47zejsBvCJcaliv2M2QGwG8IlyoodUbYTMX+/bvuImAA+AbhUkFznRGmFvvV5peECwopZQqvet8c2speGC74vHBPt3Y/0z/oHJsf146l+wcKnpvaXdl8j2wqNePjRPQhu4xynuUw38Nt6I2hsNIQLhVUzowwiv0A/IxwqTBCAkAQ0f4FAGAdVy5AHVJrCWY3Y3W6r1lzUWsg4c64/qJmfUa5H44a9Rm1rtJm1E6Ux6aVus1sMkodx601jfk7hY33TydSBR8717qK2xRqN/XcYoZw8Ql1xhlTkwHUOsKlxpkzzkSYmgyg9hEuZeobnsjPBqtEh2NzxhlTkwH4AeFShr7hCdn89MsyMZnO/6wS3Y2ZcYZSqXUWt3YoZj1ClU1e1H9grnsx1rYUZNQ1Mt3znXNLDOuvqdR5zHrMbG1dCp6XS5ua6b+/c9+tdmKuj1Ef69pup45rLCbCpQyJZEomJtOyZ+s6WbWwTUSohwCACOFixaqFbbK6t8Pr0wCAmkG4AFVma2jEbbqx9n4uLVbcOhubw0shY7pxRhnSChnTlrXpx+Z7nPhAOTfzs1CmN5tTmFtbZnycyezmbE6FNoe0NMpn5TYMZ553yKVtjPp+bp/3tPfw+RAaiygBANYRLgAA6xgWm4Pc9ONKTD0uFnu9AKhlhEuJzOnHlZh67Ia9XvzP1ti523Rj9T3c2qGYtZqGq3oLvp9Z5zBrGdpj1enGk1Pasanutvzt6KDxD7ThUef2wgX6seSEc9v4ndLXXpW/HTnZr5+LWw3GrIEo9yPmFgPGtOliudXD3JSyE2Yt7pJJuJTInH5c7SsH9noB4AeEyxx5Of2YRZUAah3hAvhUsUMh5lRk19Xkyqp8deqviIjE27W7mZam/O3ISNI45rznZLxNO9Z0xhkKG71RH/pqHI3nb0eHL+vnvcR5/0hSH85qGBrL355Yu0w71vzH97T76tBXqLOj4DFzGMxtOE39jKdP/S5umrKplI7NtTIUpmK2GADAOsIFAGAd4QIAsI6aC+BTxY6zu7UccWsNY3ZBzkbnaffDF52aSLpDrx1ERpxpw5nGiHZsYpnzOuHJrHYs1e48tumM0U1ZeZ3QZEZ/3pJ4/nbs2JB+3uaUZkVaqRtNY3ZlVusxxvRm9c9C3ZVSxP3zL2WacrHtfmoFVy4AAOsIFwCAdYQLAMA6ai5AHapEO5DQ4AX9vrYORq9dpJQ1KZmo/m/YpqFLBd9j7GqnrpC8Rl9XMxVzXqf1tP68rPIeifWLtGPzjuv1ictdzfnbLQeP6S+k1Gey11+tH1N+/8jSJfr7J0akELWuVcqfhbkGqRJ1lkq2jeHKBQBgHeECALCOYTGgDrm2gylyh8XMhWHtmLkzpdoqJtulT1NufMsZtxq+Wx9eGlzrvEfPG/oQWUb5RkrN0//tO3xtWHlcs3asccyZmjzVHNKOpVv1KdXqsJzZ/kWdmqxOtRbRd+LMRvWvzmltZBRZZUqz6y6Y4t4qxu3Pba5DWpVsG8OVCwDAOsIFAGAd4QIAsI6aC1DnXMfqjdYk6nRXs43JNMpz1Zb3IiIJpc5i1kCi407LF7N1fmqeU/NIXqUdksZh5Twb9NdMXOd8lXWcSOvvN3hRu6+2n0ncoG8H0PmWsx3ApeVx7Vjzn51dMs1WOCqzHiMrlzq3T3ygHQp36u8xdbqv4OuqNa90/0DBx5m82qWSKxcAgHWECwDAOobFitA3PCGJ5JXLyWMD47M8GqgtbrsfqlOPRUTSCeexWodkmaGDcm+Pc9t4z/Z3nGGygdv0lfbtJ6fyt8+v0YelGsecV7o8rv/bNzboHBtdaQ61ObcT1+pdmC92Fx7ei7+nD8tNLHa2D2+Y0Dsvj9+2wjnP0SkppMHYQTMyOJy/nW3TpyKnjc7LDVf1FjzmNhTmNvTl1S6VhMss+oYnZPPTL8vEpDOOG4tGpLO18Dg2AAQd4TKLRDIlE5Np2bN1naxaeOVfWZ2tjdIbj83yTAAILsKlSKsWtsnq3sKrcAEADsIF8AlzSrGNsfR0IlHwPabtoGjsvqhVPeJ6XSXT4tRnel7R30OdCtw4rldr1BYv89/WpxSPrHRqKam4Xg9pP+GcTeJDYtDrM2GlXKLufCmiT3E2d8nMRJ1juSnU4WxGVo+dlAUTo3IhOk/ebFsmqU59VCP8vjMV2tzBcloda9L5zMNGfUYa4/mbZv1lrp2XK4lwAYA5uOPCW/LFk7+V7pSz/mUw2i7fuf4/ySsLV3t4ZrWBqcgAUKI7LrwlT767V7qUYBERWTA5Kk+9+RO5c+BNj86sdhAuAFCCcDYjXzz5WxExB9yufKFmReS/vftrCWcz5lMDhWExwCdKaaPv9thptZRCr2nWA1xqMJM36L1amk46uzamlsS1Y1Mtzr9pzdYwLWedOsfYUr0eklzh1GCyEb0eMnST85rZJv1LfeG6Ie3+2XPOupdMg76Dprpe5lKncW6DV1537eAJbSjMFBaRhZdH5LrQaTnSs0oaViwp+NiQ2dY/OeHc0Xb6FMkY6160Yy67VKrrmsw/Q1ruA0CN6LpUOFhUCybGZn9QHePKpc6o3QRYjwPYN9TcPvuDROR8rHBzyyAgXOpAriXN+WRKHv7x4Xw3gVg0Ivt33UXABIDbroWltAPRprTOMnymto5Rd54UEckqU5MbhvXdJtuUTsgXl+nTbUeXOV9JGePbafE1zrDQmdPz9fdrcIbCllx1QTs2emCR/kIfcYaimm/Xh8wmDnUp768Pi12OXxmme+maVXLujQ5ZmByZcegnIyL98+LyhzWrJBMOS8OE8zvGzk7oDx7Wr4LUqcpmx+RihzOnnY8yZDbbTpja+5U5ZEa4+Fhna6PEohHZsfdo/mexaERe/Of1kkimZMfeo5JIpggXwKJMOCz/a+O98j/3vygZ0WsLuYj7l033SiYc7KoD4eJjvfGY7N91V34YTMQZCnuzb8TDMwPq27+tXCP/ffM2+fLBX8qipPN3rb8tLv9y971y4No1Hp5dbSBcCsjVLmq9C3JvPMaVCeCBf1u5Rl5avlpuOndc5k+NyWBruxxZcrWkG4N9xZJDuMzA7IRMF2TUK3MKszo+77aDpSnr1hqmW6+PTMWd6b9mW/uGCWeK8aUuveYxOOy054/F9TrOdd1OPebv+6/RX/M/6O1nOiPOlOb2Zn0q8EiHcz7puN5WP3bS+TymYrlzi8jr11wrzYkrz2scEWm4pLetUdvzT8b1qc8iPdq98F+PSyHqTpTm5+1G/TM1pyxXcpdKwmUGZidkZl0BQGkIFxd0QgaAuSFcgDpU7Gp+tynMs63QV3eqNKfNivLcdIv+OpmoU5OIXNSHnhb8xRmmmpynrxPJNjvvPz6sr16/NN/5Kvun/3xYO/Z/37teu3/rig/yt++Z/+/asadGtjivOa4PYUWUkbhkr94hIN3s/E5Ro9PzyMrm/O3uP+uft7rzpYhITK7O3zabx0RO9ksh6rRws/Oy9rgKdNYuhMoTAMA6wgUAYB3hAgCwjpoLEDDqOHspY/BuNZi00bE30tPt3DZqBQ1DTr0kG9W/goY2Os+72KPXLnr+t9NSZurDesfkM0q7mXiT3mIlk9Yfe1XMmZr8vZObpKAR/fdNKm1j2g/r9Zis8ms0J/TzznVTFhEZX6K/ZsuAXnMKTTqPDZsdk6PKc81dKpU/i+m7W3qzMyXhUufURaBMqQZQLYRLnSrUd4xGlgCqgXD5B7VVfa23fCmG2Xfs2MA4jSwBVA3hItPbvYjUR8sX+o5hNm7rXKY91ljnorYjEeNYJjHsvKZRH1B3VMxet0I7ptYgGi7ptZLYOacGEV2qr3MZ63NqLv2tF7VjqxYPaPfPp5w2Mr/+8L9qx+7803917vQW/kfmiv+ibzHw1lmnjUvoD23asfh7zmdjtvFXd+UUEUm3OvUSs+aiChm7VIr6ESf0prXqO6jrYURE0onC9TfVXOo2hItMb/ciQn0CAMpBuCho9wIAdhAuAIpiDo1kLgznb4eNoS+tHYwxNVZ73IjeqqRFGQpqbtGn+6a6nZGEJmO6r7pbfd95/R+Izc36kN3fTzk7U36s88/ascR5p+VMdkofwlqx3BnO+8vRldqx6Ljz2El9VEzGlzhfs2F95rEseGNYuz+xzHn/xjP6g7OTynRjKZ46FJZOJAo+zvaUZRZRAgCsI1wAANYRLgAA66i5AMhzaw3j9lh1SquISEhp+262I1HrMZmWZikkfEZvKXPhjlUFH7vkd87tsaV60aPhn4a0+2tWOtOIP9Gi72i5XWn5snbdCe3YX05c5byGcex4wtltc2xIrz9NXOVMqW59X59efXmJUaBRjy3Xd/BsPKN8XU/q9Rj1cwxNGtPCldqYyW2LBXPXylJx5QIAsI5wAQBYx7AYgBmVMjU10tmp3c9ow2L68NrU6T7nWH/hrsyXP36rdmzRq87K80uL9RXq0VFnKGiqWR9q6+9v1+6PjDrDVjfs/6J2bM3m9/K3r2oZ1o4djzvDVIMT+tDXvGZnCvXEUFw7tvig0/nj7EZ9CvXFhfpXcMMl53jr3/RhQW3oK6l3ftauEozV+2Gle0LGaDqgDoWVOwxm4soFAGAdVy4Bk2vKSXsbAJVEuASE2YKf9vsAKolwCQi1BT/t92GDOo3Vra2IOU1Zq8+4tIZpea/wazaf1e8P3ObUVcydIBfv17/mLnzYOe9UR0Y7pk4pfvP0Eu1YZshpR3Prhje1Y797fbXz/qvGtGOJ8865tfZph6T5gj6lWK0d9X1Cf/+uN51p0w1RfUpzZHBYClI+Y9edRo2p5+XWYwiXAKEFP4BqoaAPALAu0Fcuud0n62HnSQCoJYENF3P3yXrYeRKopmLXwZhj+Wp9JtKzUH+wuibj+Cn9dZQagNniv+cVZy1JyGiN0vdx/T1azjk1mUtd+uDNuLKj5eJr9HUm4284rfp/13G9dkxtuZ96X2/p0qqcTiquHZq2zqVtyjm3lkG9HhRJOp+N2RpH3SXU3KUyq6w5MtcjudXKVLO1AppJYMPF3H2SqbkAYE9gwyWH3ScBwL7AhwuA4phDI6VMay0k3T+g3VeHbcwuvdrrG0M/oWFn+u/la3u0Y4v+YEwNvsEZtmoe0qctT7Y6w2SjBxZpx1JdzmOjffrvN9mbKnxMGSVbeFjvWKzuUikiEp50hsIaR9PascjJ/vxt8/eXqPO62eRF7ZA6ZGZOKQ636sOLKroiAwBqDuECALCOcAEAWEfNJcDU9T3MloNNxe5oaY75Z1x2sNQe169PxZ1c70wNbjp5QT+XqP4113EsohzT/32daXDa2l/q0t8zcsmZbhw1lsZ1vu2ca+OYXisZXea8/6X5etuWlkH9sUOrnfdX272IiKRucHbCbDwzqh0LKXUVWbhAO6b/hoU/b3N6t/o82r+gKGYTSxEaWQKwi3AJILWJpYjQyBKAdYRLQNHEEqUqZWfKYp9XynCL2+s0vnXaeVxc33lShvUhpPTyeP525KK+mr/7T8403ssL9OG8oTXO12XrWX31vGoqVriU3TKgv9/oMn3or+2sM0zWMHxZOzYVb5JCspPKVOTTwwUfZ9I6H4/rfxbFDm0WQkEfAGAd4QIAsI5wAQBYR80FQEW5tY0pZSxffaxZf1E7/2YvDOvPM6Y0R5VaRuTd0/pjlbYq0UG9JtkwEc/fbhyc0I6NXj8vf7ulXz+3xlFn+rH5vMb5+ldw7JxzbqFJfZqyOsXabPGiCl21WLufPe1s2xlq1D9vs85SyFzqbVy5AACsI1wAANYRLgAA66i5IC/XDoZWMLDJHK93q5241WDUx05rG5MsvNvitLrCoX933u/qFfp7KK1iJrv1HSVjbzm1i2yr+ffDqbmoO0aK6OtVUt3688z6jCp80Wj/siSevx19t3DNJf3OMe2++nmYn4Vaj1Jb89tAuGBaO5hYNCLP3X+zLPjHts+EDYBSES7Q2sGcT6bk4R8flm3PH8ofp+8YgFIFKlz6hie0flpwqO1g6DuGSnKb1lrs7pZubWOyKfdps+rrpD84ox1TOwNHWvR2K5n4PCmk/R1nt8t0qz71ufGMM4U4anRhNqdCZ1YscV6nQx/6iw4631np5fpum+oulSb185g2hbvIaeFzaf8SmHDpG56QzU+/LBPK3PFYNCKdraV/aPWOvmMAyhWYcEkkUzIxmZY9W9fJqoVXCnXUEgCgMgITLjmrFrbJ6t4Or08DAOpa4MIFgD/MtW2MWY9xe50GYyqyWoMJ/f19/XWUtipqSxURkZByu2HlUuP9nCm+4Yvu030jI0obm8RI4ccN63WddP+Ac8xlKrb5WWgt9+ew26QbFlECAKyr+yuX3AwxZocBQPXUdbiYM8SYHQb4RylTlotd2W8+dur4+wWPmcNEYZfuwqFOp46bfue4dixz0/XO4/70jv68+XH9sf2Dhd9D6disDrXNxu1zVKdeh4wV+urz5tIVua7DxZwhxuwwAKiOug6XHGaIAUB1UdAHAFgXiCsXAMhxm9Js7lqp0jovG/UQddqwWZ/QWry0Gd2cjV0zw1cvc17HmO6stXEx6iNql2iz/Y36O5q/XzqRkErhygUAYB3hAgCwjmExFIWNxACUgnCBq5k2EmNvF1RKKetVij1WCrWW4faabrUK83dQH2vuoGnWQDLHTxV+XeWx5vOKPW9bn1Mx6i5c2LPFLnUjMfZ2AVCsugoX9mypDPZ3AVCqugoX9mwB/K2awzY230/tRKx2IRYpbRfHYs/HtaWLMfSmTqGea5sc2r/8AyvyAcBbTEUGAFhHuAAArKvLYTEAmIu51hnMOovKrXW9uWukFFkfcT1mtKYRlx0mKzltmSsXAIB1XLmgZOr6IWbjAZgJ4YKimav1RVixj/oy16GguT7PbaX/3F+zMjt4lopwQdHU1foiwop9AAXVRbjkWr7Q7qXyWK0PoBi+Dxez5QvtXqqPjskATL4PF7PlC19w1UPHZKA0Zo2j2u1uTJV8f9+HSw4tX6qPjskACqmbcIE3qMEAmAmLKAEA1vnyyoUNwWoXCyyBwqpRY3Gr61Sz5uO7cGFDsNrEAksAKt+FCxuC1SYWWAJQ+SZczIWSzA6rPTMV91kDA1RXJTsdl8IX4cJCSf9hDQwQbL4IFxZK+g9rYIBg80W45DAU5i+sgQGCy1fhAn9jmjIQHIQLKq7QNOXn7r9ZFrQ2EjRAHaqpcFEXR6pYKOlv5jTl88mUPPzjw7Lt+UMiQrEfqEeehosaJrkvHHVxpIoZYv5m1l/MYv8fT1yQBOuWgDlRV97Pdbqx7dX7noVLoZX2L/7zelkwQ4jwhVNfcmEz25CZCH/2gB95Fi6stIfI7ENmItPDRsX/M0BtKjpcBkYvycDYZWtvzEp75BQaMhOZOWxUbsGj4v8xoLpC2Ww26/VJAADqC/u5AACsI1wAANYRLgAA6wgXAIB1hAsAwLqipiJns1kZGxur9LkAFTVv3jwJhUJenwYQCEWFy9jYmHR0sE4A/jYyMiLt7e1enwYQCEWtcynnymV0dFSWLl0qH3zwAX+xLeJzLR1XLkD1FHXlEgqFyv4Ca29v50uwAvhcAdQiCvoAAOsIFwCAdRUPl6amJtm9e7c0NTVV+q0Chc8VQC2jcSUAwDqGxQAA1hEuAADrCBcAgHWECwDAOivh8uyzz8qKFSukublZNmzYIIcOzbwlbc6ePXvk+uuvl1gsJkuXLpUvfelLcunSJRunUvd+//vfy5YtW2TJkiUSCoXkl7/8pdenBADTlB0ue/fulZ07d8ru3bvlyJEjsnbtWrnnnntkYGBgxsf/5Cc/kccee0x2794tb7/9tvzwhz+UvXv3yle+8pVyTyUQksmkrF27Vp599lmvTwUACip7KvKGDRvk1ltvle985zsiIpLJZGTp0qWyfft2eeyxx6Y9/pFHHpG3335bDhw4kP/Zrl275PXXX5dXX321nFMJnFAoJL/4xS/k3nvv9fpUAEBT1pVLKpWSw4cPy+bNm50XDIdl8+bNcvDgwRmfc/vtt8vhw4fzQ2fHjx+X3/zmN/Lxj3+8nFMBANSQohpXFjI0NCTpdFp6enq0n/f09Mjf/va3GZ/z6U9/WoaGhuSOO+6QbDYrU1NT8vDDDzMsBgB1pOqzxV566SX5xje+Id/97nflyJEj8vOf/1x+/etfy9e+9rVqnwoAoELKunLp6uqSSCQi/f392s/7+/tl0aJFMz7niSeekPvvv18+97nPiYjIjTfeKMlkUh566CH56le/KuEws6MBwO/K+iZvbGyUm2++WSvOZzIZOXDggGzcuHHG51y8eHFagEQiERG5sikZAMD/yrpyERHZuXOnbNu2TW655RZZv3697NmzR5LJpDz44IMiIvLAAw9Ib2+vfPOb3xQRkS1btsgzzzwjN910k2zYsEGOHTsmTzzxhGzZsiUfMihsfHxcjh07lr9/4sQJOXr0qMyfP1+WLVvm4ZkBgKPscNm6dasMDg7Kk08+KefOnZN169bJvn378kX+U6dOaVcqjz/+uIRCIXn88celr69Puru7ZcuWLfL1r3+93FMJhDfeeEPuvvvu/P2dO3eKiMi2bdvkhRde8OisAEBHy30AgHVUzwEA1hEuAADrCBcAgHWECwDAOsIFAGAd4QIAsI5wAQBYR7gAAKwjXAAA1hEuAADrCBcAgHWEiwc2bdok27dvlx07dkhnZ6f09PTI97///Xw36Xnz5smqVavkt7/9rdenCgBzQrh45MUXX5Suri45dOiQbN++Xb7whS/Ipz71Kbn99tvlyJEj8tGPflTuv/9+uXjxotenCgAloyuyBzZt2iTpdFpeeeUVERFJp9PS0dEhn/zkJ+VHP/qRiIicO3dOFi9eLAcPHpTbbrvNy9MFgJJx5eKRNWvW5G9HIhFZsGCB3Hjjjfmf5fbDGRgYqPq5AUC5yt4sDHMTjUa1+6FQSPtZKBQSkSvbRgOojL7hCUkkU/n7na2N0huPeXhG9YNwARBIfcMTsvnpl2ViMp3/WSwakf277iJgLCBcAARSIpmSicm07Nm6TlYtbJNjA+OyY+9RSSRThIsFhAuAQFu1sE1W93Z4fRp1h3DxwEsvvTTtZ++///60nzGRD4BfMVsMAGAd4QIAsI5wAQBYR7gAAKyjoA8gMNRFk8cGxj0+m/pGuAAIhEKLJjtbGz08q/pFuAAIBHPRpAjtXiqJcAEQKCyarA4K+gAA6wgXAIB1hAsAwDpqLgAwA/Z6KQ/hAgAG9nopH+ECAIrc4kr2eikP4QIAcmXYKxaNyI69R0XkypXKrSvnEyZzRLgAgIj0xmOyf9dd+ToLNZbyEC4A8A+98RiBYglTkQEA1hEuAADrCBcAgHWECwDAOgr6AFAiVu/PjnABgBKwer84hAuAulLJqwpW7xePcAFQNwpdVTx3/81a4JSK1fulI1wA1A1zK+PzyZQ8/OPDsu35QyJyJRQ6WxtLfl1W75eOcAFQd9StjG2FAqv3S0O4AKhrhII3WOcCALCOcAEAWMewGABYkJumTLH/CsIFAMow0zRlFlQSLgBQFnWaMgsqHYQLAJSJGWnTES4AfC/X8iVX94D3CBcAvma2fJnrKnzYRbgA8DWz5QuztWoD4QKgLqgtX+A9FlECAKwjXAAA1hEuAADrCBcAgHWECwDAOsIFAGAd4QIAsI5wAQBYR7gAAKxjhT4AWKY20AxqOxrCBQAsMTcOEwnu5mGECwBYom4cJiKB3jyMcAEAi9g47AoK+gAA67hyAeA7uZ0nRYTdJ2sU4QLAV8ydJ0XYfbIWES4AfMXceVIkuNN9axnhAsCX2HmythEuAHwhV2ehxuIPhAuAmmfWWaix1D7CBUDNM+ss1FhqH+ECwDeos/gHiygBANYRLgAA6wgXAIB1hAsAwDrCBQBgHeECALCOcAEAWEe4AACsI1wAANaxQh8AKizXbDNIbWsIFwCokM7WRolFI7Jj71ERudJwc/+uuwIRMIQLAFRIbzwm+3fdld8qYMfeo5JIpggXAEB5euOxQISJiYI+AMA6wgUAYB3DYgBqUm5bYxFha2MfIlwA1BxzW2MRtjb2G8IFQM0xtzUWCdYakXpAuACoWWxr7F+ECwBUkVo/querMcIFAKrAXK0vUt8r9gkXAKgCdbW+iNT9in3CBQCqJEir9VlECQCwjnABAFhHuAAArKPmAqBm5Fq+0O7F/wgXADXBbPlCuxd/I1wA1ASz5Us9LzAMAsIFQE2h5Ut9oKAPALCOcAEAWEe4AACsI1wAANZR0AcAD+XW9NTb7DjCBQA8YLbgr7f2+4QLAHhAbcFfj+33CRcA8Eg9t+CnoA8AsI5wAQBYR7gAAKwjXAAA1hEuAADrCBcAgHWECwDAOta5APBMbltjEWFrY9E/A7+3gyFcAHjC3NZYJLhbG5utYET83w6GcAFQVbmrlWMD49q2xiL+/9f6XKmtYESkLtrBEC4Aqsa8WolFI3Lryvm+/QK1qd5awRAuAKomkUxpVytBvVIJAsIFQNWtWtgmq3s7vD4NVBBTkQEA1hEuAADrCBcAgHWECwDAOgr6ACqKVfhzl/u8/DirjnABUDGswp8bc8W+H1frEy4AKsZc1yLiz3+FV5u6Yt+vq/UJFwAVx7qW0vl9xT4FfQCAdYQLAMA6wgUAYB3hAgCwjnABAFhHuAAArCNcAADWsc4FAHxAbZ3jh4WohAsA1DCzFYyIP9rBEC4AUMPUVjAi4pt2MIQLANQ4P7aCIVwAWJdrs0+L/eAiXABYZbbZp8V+MBEuAKwy2+z7YWYT7CNcAFQEbfaDjXABUDa2MoaJcAFQFrYyxkwIFwBlYStjzIRwAWAFNRaoaFwJALCOKxcA8KHcxIlaHYIkXADAR8xGlrXaxJJwAQAfURtZ1nITS8IFAHzGD40sKegDAKzjygUAfK4Wd6kkXAAURW3xoqLdi3dqeZdKwgXArGZq8aKi3Ys3anmXSsIFwKxmavGiqpWhmCCq1eI+4QKgaLR48YdaWGBJuABAnailBZaECwDUiVpaYEm4ACgoN0OMGWH+USs1GMIFwIzMGWLMCEMpCBcAMzJniDEjDKUgXAC4YoYY5oJwAZCnrsKnzoJyEC4ARGTmVfjUWfzPq75jhAsAEZl5FT51Fv/yuu8Y4QJAQ42lPhTqO/bHExckUYUJGoQLEDB0Nw4Odc1LtVfvEy5AgNDdOLhmWr2fu4oRsT8ESrgAAUJ342DLXckUqsc8d//NsqC10cr/B6FsNpst83wB1Di1jcuOvUfl/2y/g7pKwKnDo+eTKXn4x4e1bgzlDplx5QLUOdq4YCZmDzLbDS8JF6AOFCrSi1wp1NPGBbMxw6bc9TGEC1BlZhC4/cV1C40cc0hjJrFoRG5dOZ9Qwaxmq8eISFFDqtRcgCr78BP7pq2CV//i5hQTGrO9Rg5XKyiFWz1GROT9//GJWV+DcAEAWBf2+gQAAPWHcAEAWEe4AACsI1wAANYRLgAA61jnAlRRNpuVsbExr08DKNu8efMkFAoVPE64AFU0NjYmHR309IL/jYyMSHt7e8HjrHMBqigIVy6jo6OydOlS+eCDD1y/fOpF0H7fHK5cgBoSCoUC8wXU3t4emN9VJHi/72wo6AMArCNcAADWES4ArGpqapLdu3dLU1OT16dSFUH7fYtFQR8AYB1XLgAA6wgXAIB1hAsAwDrCBQBgHeECwJpnn31WVqxYIc3NzbJhwwY5dOiQ16dUMU899ZSEQiHtvw996ENen1bNIFwAWLF3717ZuXOn7N69W44cOSJr166Ve+65RwYGBrw+tYr5yEc+ImfPns3/9+qrr3p9SjWDcAFgxTPPPCOf//zn5cEHH5QbbrhBnnvuOWlpaZHnn3/e61OrmIaGBlm0aFH+v66uLq9PqWYQLgDKlkql5PDhw7J58+b8z8LhsGzevFkOHjzo4ZlV1rvvvitLliyRq6++Wj7zmc/IqVOnvD6lmkG4ACjb0NCQpNNp6enp0X7e09Mj586d8+isKmvDhg3ywgsvyL59++R73/uenDhxQu68886673pdLLoiA8AcfOxjH8vfXrNmjWzYsEGWL18uP/vZz+Szn/2sh2dWG7hyAVC2rq4uiUQi0t/fr/28v79fFi1a5NFZVVc8HpfrrrtOjh075vWp1ATCBUDZGhsb5eabb5YDBw7kf5bJZOTAgQOyceNGD8+sesbHx+W9996TxYsXe30qNYFhMQBW7Ny5U7Zt2ya33HKLrF+/Xvbs2SPJZFIefPBBr0+tIr785S/Lli1bZPny5XLmzBnZvXu3RCIRue+++7w+tZpAuACwYuvWrTI4OChPPvmknDt3TtatWyf79u2bVuSvF6dPn5b77rtPzp8/L93d3XLHHXfIa6+9Jt3d3V6fWk2g5T4AwDpqLgAA6wgXAIB1hAsAwDrCBQBgHeECALCOcAEAWEe4AACsI1wA+N6mTZtkx44dXp8GFIQLAMA6wgUAYB3hAqAuTE1NySOPPCIdHR3S1dUlTzzxhNDdyjuEC4C68OKLL0pDQ4McOnRIvvWtb8kzzzwjP/jBD7w+rcCicSUA39u0aZMMDAzIX//6VwmFQiIi8thjj8mvfvUreeuttzw+u2DiygVAXbjtttvywSIisnHjRnn33XclnU57eFbBRbgAAKwjXADUhddff127/9prr8m1114rkUjEozMKNsIFQF04deqU7Ny5U9555x356U9/Kt/+9rfl0Ucf9fq0AottjgHUhQceeEAmJiZk/fr1EolE5NFHH5WHHnrI69MKLGaLAQCsY1gMAGAd4QIAsI5wAQBYR7gAAKwjXAAA1hEuAADrCBcAgHWECwDAOsIFAGAd4QIAsI5wAQBYR7gAAKz7/0zhCg8nJIBGAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "# sample from the posterior\n", + "posterior_samples_1 = posterior.sample((10000,), x = y_true)\n", + "# that last little part is conditioning on a data value\n", + "# plot posterior samples\n", + "fig, axes = sbi.analysis.pairplot(\n", + " posterior_samples_1, \n", + " labels = ['m', 'b'],\n", + " #limits = [[0,10],[-10,10],[0,10]],\n", + " truths = theta_true,\n", + " figsize=(5, 5)\n", + ")\n", + "axes[0, 1].plot([theta_true[1]], [theta_true[0]], marker=\"o\", color=\"red\")" ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "639399bb-378f-4944-af34-4342ff569bc0", + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { diff --git a/src/scripts/evaluate.py b/src/scripts/evaluate.py index 16e54ea..a62c307 100644 --- a/src/scripts/evaluate.py +++ b/src/scripts/evaluate.py @@ -28,10 +28,14 @@ def load_model_pkl(self, path, model_name): :param model_name: Name of the model :return: Loaded model object that can be used with the predict function """ + print(path) with open(path + model_name + ".pkl", 'rb') as file: posterior = pickle.load(file) return posterior + def infer_sbi(self, posterior, n_samples, y_true): + return posterior.sample((n_samples,), x=y_true) + def predict(input, model): """