diff --git a/unit8/SupplementaryExercises83.ipynb b/unit8/SupplementaryExercises83.ipynb index bb5ac373..43c97b8a 100644 --- a/unit8/SupplementaryExercises83.ipynb +++ b/unit8/SupplementaryExercises83.ipynb @@ -50,21 +50,397 @@ "```{admonition} Solution\n", ":class: tip, dropdown\n", "\n", - "Solution to be added. Feel free to share yours on Ed Discussion!\n", + "Code in hidden cell below!\n", "```" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 1, "id": "a6d402b7", "metadata": { "tags": [ "hide-cell" ] }, - "outputs": [], - "source": [] + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Initializing NUTS using jitter+adapt_diag...\n", + "Multiprocess sampling (4 chains in 4 jobs)\n", + "NUTS: [alpha, beta0, beta1]\n" + ] + }, + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "7b603102c2a740ecbd24730a25603948", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "Output()" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "
\n"
+      ],
+      "text/plain": []
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    },
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "Sampling 4 chains for 1_000 tune and 3_000 draw iterations (4_000 + 12_000 draws total) took 2 seconds.\n"
+     ]
+    },
+    {
+     "data": {
+      "application/vnd.microsoft.datawrangler.viewer.v0+json": {
+       "columns": [
+        {
+         "name": "index",
+         "rawType": "object",
+         "type": "string"
+        },
+        {
+         "name": "mean",
+         "rawType": "float64",
+         "type": "float"
+        },
+        {
+         "name": "sd",
+         "rawType": "float64",
+         "type": "float"
+        },
+        {
+         "name": "hdi_3%",
+         "rawType": "float64",
+         "type": "float"
+        },
+        {
+         "name": "hdi_97%",
+         "rawType": "float64",
+         "type": "float"
+        },
+        {
+         "name": "mcse_mean",
+         "rawType": "float64",
+         "type": "float"
+        },
+        {
+         "name": "mcse_sd",
+         "rawType": "float64",
+         "type": "float"
+        },
+        {
+         "name": "ess_bulk",
+         "rawType": "float64",
+         "type": "float"
+        },
+        {
+         "name": "ess_tail",
+         "rawType": "float64",
+         "type": "float"
+        },
+        {
+         "name": "r_hat",
+         "rawType": "float64",
+         "type": "float"
+        }
+       ],
+       "ref": "8708461e-2f0c-44c9-8f7f-26346d5eb096",
+       "rows": [
+        [
+         "alpha",
+         "1.37",
+         "0.232",
+         "0.936",
+         "1.801",
+         "0.004",
+         "0.003",
+         "3480.0",
+         "4173.0",
+         "1.0"
+        ],
+        [
+         "beta0",
+         "-4.45",
+         "0.749",
+         "-5.943",
+         "-3.153",
+         "0.013",
+         "0.009",
+         "3552.0",
+         "3973.0",
+         "1.0"
+        ],
+        [
+         "beta1",
+         "-0.686",
+         "0.45",
+         "-1.557",
+         "0.137",
+         "0.006",
+         "0.005",
+         "5719.0",
+         "5149.0",
+         "1.0"
+        ],
+        [
+         "mean0",
+         "24.849",
+         "7.131",
+         "14.939",
+         "37.348",
+         "0.089",
+         "0.263",
+         "7245.0",
+         "6449.0",
+         "1.0"
+        ],
+        [
+         "mean1",
+         "41.615",
+         "13.881",
+         "23.289",
+         "65.1",
+         "0.175",
+         "0.377",
+         "7367.0",
+         "7045.0",
+         "1.0"
+        ],
+        [
+         "median0",
+         "18.541",
+         "4.882",
+         "10.894",
+         "26.944",
+         "0.059",
+         "0.159",
+         "7516.0",
+         "7403.0",
+         "1.0"
+        ],
+        [
+         "median1",
+         "31.003",
+         "9.423",
+         "17.703",
+         "47.228",
+         "0.115",
+         "0.217",
+         "7949.0",
+         "7223.0",
+         "1.0"
+        ]
+       ],
+       "shape": {
+        "columns": 9,
+        "rows": 7
+       }
+      },
+      "text/html": [
+       "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
meansdhdi_3%hdi_97%mcse_meanmcse_sdess_bulkess_tailr_hat
alpha1.3700.2320.9361.8010.0040.0033480.04173.01.0
beta0-4.4500.749-5.943-3.1530.0130.0093552.03973.01.0
beta1-0.6860.450-1.5570.1370.0060.0055719.05149.01.0
mean024.8497.13114.93937.3480.0890.2637245.06449.01.0
mean141.61513.88123.28965.1000.1750.3777367.07045.01.0
median018.5414.88210.89426.9440.0590.1597516.07403.01.0
median131.0039.42317.70347.2280.1150.2177949.07223.01.0
\n", + "
" + ], + "text/plain": [ + " mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk \\\n", + "alpha 1.370 0.232 0.936 1.801 0.004 0.003 3480.0 \n", + "beta0 -4.450 0.749 -5.943 -3.153 0.013 0.009 3552.0 \n", + "beta1 -0.686 0.450 -1.557 0.137 0.006 0.005 5719.0 \n", + "mean0 24.849 7.131 14.939 37.348 0.089 0.263 7245.0 \n", + "mean1 41.615 13.881 23.289 65.100 0.175 0.377 7367.0 \n", + "median0 18.541 4.882 10.894 26.944 0.059 0.159 7516.0 \n", + "median1 31.003 9.423 17.703 47.228 0.115 0.217 7949.0 \n", + "\n", + " ess_tail r_hat \n", + "alpha 4173.0 1.0 \n", + "beta0 3973.0 1.0 \n", + "beta1 5149.0 1.0 \n", + "mean0 6449.0 1.0 \n", + "mean1 7045.0 1.0 \n", + "median0 7403.0 1.0 \n", + "median1 7223.0 1.0 " + ] + }, + "execution_count": 1, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "import pymc as pm\n", + "import arviz as az\n", + "import pandas as pd\n", + "import numpy as np\n", + "import pytensor.tensor as pt\n", + "\n", + "data = pd.read_csv(\"../data/dukes.csv\")\n", + "\n", + "censored = data[\"censored\"].to_numpy()\n", + "treatment = data[\"treatment\"].to_numpy()\n", + "months = data[\"months\"].to_numpy()\n", + "\n", + "with pm.Model():\n", + " alpha = pm.Exponential(\"alpha\", 1)\n", + " beta0 = pm.Normal(\"beta0\", 0, 100)\n", + " beta1 = pm.Normal(\"beta1\", 0, 100)\n", + "\n", + " beta = pm.math.exp(beta0 + beta1 * treatment) ** (-1 / alpha)\n", + "\n", + " obs_latent = pm.Weibull.dist(alpha=alpha, beta=beta)\n", + " pm.Censored(\"lik\", obs_latent, lower=None, upper=censored, observed=months)\n", + "\n", + " # https://docs.pymc.io/en/latest/api/distributions/generated/pymc.Weibull.html\n", + " pm.Deterministic(\n", + " \"mean0\",\n", + " pt.gamma(1 + 1 / alpha) * (pt.exp(beta0 + beta1 * 0) ** (-1 / alpha)),\n", + " )\n", + " pm.Deterministic(\n", + " \"mean1\",\n", + " pt.gamma(1 + 1 / alpha) * (pt.exp(beta0 + beta1 * 1) ** (-1 / alpha)),\n", + " )\n", + "\n", + " pm.Deterministic(\n", + " \"median0\",\n", + " (np.log(2) * (pt.exp(beta0 + beta1 * 0)) ** (-1 / alpha)),\n", + " )\n", + " pm.Deterministic(\n", + " \"median1\",\n", + " (np.log(2) * (pt.exp(beta0 + beta1 * 1)) ** (-1 / alpha)),\n", + " )\n", + "\n", + " trace = pm.sample(3000, target_accept = .9)\n", + "\n", + "az.summary(trace)" + ] }, { "cell_type": "markdown", @@ -98,7 +474,14 @@ "```{admonition} Solution\n", ":class: tip, dropdown\n", "\n", - "Solution to be added. Feel free to share yours on Ed Discussion!\n", + "A custom distribution is implemented in PyMC using the Rayleigh survival model:\n", + "\n", + "$$f(c,t) = \\begin{cases} \n", + " \\frac{t}{\\lambda} \\exp(\\frac{-t^2}{2\\lambda}) & \\text{if } c = 1 \\\\\n", + " \\exp(\\frac{-t^2}{2\\lambda}) & \\text{if } c = 0 \\\\\n", + " \\end{cases}$$\n", + "\n", + "Censored points are labeled with $c=0$:\n", "\n", "```\n", "\n", @@ -108,15 +491,307 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 8, "id": "019cb623", "metadata": { "tags": [ "hide-cell" ] }, - "outputs": [], - "source": [] + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Initializing NUTS using jitter+adapt_diag...\n", + "Multiprocess sampling (4 chains in 4 jobs)\n", + "NUTS: [λ]\n", + "/opt/anaconda3/envs/pymc_env/lib/python3.13/multiprocessing/popen_fork.py:67: RuntimeWarning: os.fork() was called. os.fork() is incompatible with multithreaded code, and JAX is multithreaded, so this will likely lead to a deadlock.\n", + " self.pid = os.fork()\n" + ] + }, + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "e2fca449758e4e22891d17d698108d77", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "Output()" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/opt/anaconda3/envs/pymc_env/lib/python3.13/multiprocessing/popen_fork.py:67: RuntimeWarning: os.fork() was called. os.fork() is incompatible with multithreaded code, and JAX is multithreaded, so this will likely lead to a deadlock.\n", + " self.pid = os.fork()\n" + ] + }, + { + "data": { + "text/html": [ + "
\n"
+      ],
+      "text/plain": []
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    },
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "Sampling 4 chains for 1_000 tune and 5_000 draw iterations (4_000 + 20_000 draws total) took 1 seconds.\n"
+     ]
+    },
+    {
+     "data": {
+      "application/vnd.microsoft.datawrangler.viewer.v0+json": {
+       "columns": [
+        {
+         "name": "index",
+         "rawType": "object",
+         "type": "string"
+        },
+        {
+         "name": "mean",
+         "rawType": "float64",
+         "type": "float"
+        },
+        {
+         "name": "sd",
+         "rawType": "float64",
+         "type": "float"
+        },
+        {
+         "name": "hdi_2.5%",
+         "rawType": "float64",
+         "type": "float"
+        },
+        {
+         "name": "hdi_97.5%",
+         "rawType": "float64",
+         "type": "float"
+        },
+        {
+         "name": "mcse_mean",
+         "rawType": "float64",
+         "type": "float"
+        },
+        {
+         "name": "mcse_sd",
+         "rawType": "float64",
+         "type": "float"
+        },
+        {
+         "name": "ess_bulk",
+         "rawType": "float64",
+         "type": "float"
+        },
+        {
+         "name": "ess_tail",
+         "rawType": "float64",
+         "type": "float"
+        },
+        {
+         "name": "r_hat",
+         "rawType": "float64",
+         "type": "float"
+        }
+       ],
+       "ref": "51f52fa5-b51e-4c03-a4d7-3f1ed3dd5666",
+       "rows": [
+        [
+         "S_new",
+         "0.492",
+         "0.092",
+         "0.31",
+         "0.669",
+         "0.001",
+         "0.001",
+         "8524.0",
+         "13589.0",
+         "1.0"
+        ],
+        [
+         "λ",
+         "2663.574",
+         "760.235",
+         "1423.04",
+         "4205.367",
+         "8.133",
+         "7.232",
+         "8524.0",
+         "13589.0",
+         "1.0"
+        ]
+       ],
+       "shape": {
+        "columns": 9,
+        "rows": 2
+       }
+      },
+      "text/html": [
+       "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
meansdhdi_2.5%hdi_97.5%mcse_meanmcse_sdess_bulkess_tailr_hat
S_new0.4920.0920.310.6690.0010.0018524.013589.01.0
λ2663.574760.2351423.044205.3678.1337.2328524.013589.01.0
\n", + "
" + ], + "text/plain": [ + " mean sd hdi_2.5% hdi_97.5% mcse_mean mcse_sd ess_bulk \\\n", + "S_new 0.492 0.092 0.31 0.669 0.001 0.001 8524.0 \n", + "λ 2663.574 760.235 1423.04 4205.367 8.133 7.232 8524.0 \n", + "\n", + " ess_tail r_hat \n", + "S_new 13589.0 1.0 \n", + "λ 13589.0 1.0 " + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "import pymc as pm\n", + "import arviz as az\n", + "import pandas as pd\n", + "import numpy as np\n", + "from pymc.math import log, exp\n", + "\n", + "t = np.array([23, 40, 41, 67, 69, 72, 84, 88, 100, 100])\n", + "c = np.array([ 1, 1, 1, 1, 1, 1, 1, 1, 0, 0])\n", + "\n", + "# CustomDist requires all data to be together in a matrix\n", + "val = np.concat([[t],[c]])\n", + "\n", + "# log-probability of custom distribution\n", + "def logp(value, lam):\n", + " t = value[0,:]\n", + " c = value[1,:]\n", + " return (c * log(t/lam) - ((t**2) / (2*lam))).sum()\n", + "\n", + "with pm.Model() as m:\n", + " λ = pm.Gamma(\"λ\", .001, .001)\n", + " \n", + " rayleigh_surv = pm.CustomDist('rayleigh_surv', λ, logp=logp, observed=val)\n", + "\n", + " #mean prediction\n", + " S_new = pm.Deterministic(\"S_new\",exp(-(60**2) / (2*λ)))\n", + "\n", + " trace = pm.sample(5000)\n", + "\n", + "az.summary(trace,hdi_prob=.95)\n" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "7a0cfedd", + "metadata": { + "tags": [ + "hide-cell" + ] + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "MLE version of λ: 3315.25\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAHFCAYAAAAOmtghAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjMsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvZiW1igAAAAlwSFlzAAAPYQAAD2EBqD+naQAAhAVJREFUeJzs3Xd4FNXXwPHvbnqH9ARCEnroEIr03quoVKUJyqtIEwVEqgV/ooioFKkiCAgC0kR6B+k1kRoIJSEQIIX03Xn/WFhdE0I22U09n+eZh9nZO3fOZDbsycwtKkVRFIQQQgghCgl1XgcghBBCCGFKktwIIYQQolCR5EYIIYQQhYokN0IIIYQoVCS5EUIIIUShIsmNEEIIIQoVSW6EEEIIUahIciOEEEKIQkWSGyGEEEIUKpLcCPEvS5cuRaVS6RdLS0t8fHzo1asXV65cMfvxVSoVU6ZMyfV9AwIC6NSpU7b2BYiOjmb8+PFUqlQJBwcHXFxcqFixIm+88Qbnzp3Ldr05cePGDVQqFUuXLjXbMfbu3YtKpWLv3r1ZKpfR8uqrr5otvqzYunXrcz83AQEBDBgwIFfjEcIULPM6ACHyoyVLllCxYkWSkpI4dOgQn332GXv27OHvv/+mePHieR1eho4cOULJkiVz/bjx8fG89NJLxMfH88EHH1C9enUSExO5fPky69at48yZM1SrVi3X4/Lx8eHIkSOUKVMm14/9PJ9//jnNmzc32Obm5pZH0ehs3bqVH374IcMEZ/369Tg7O+d+UELkkCQ3QmSgSpUq1K5dG4BmzZqh0WiYPHkyGzZsYODAgXkcXcZeeumlPDnumjVruHr1Krt37073xT169Gi0Wq1JjqMoCklJSdjZ2WWpvI2NTZ79TJ6nXLly+S6mzNSsWTOvQxAiW+SxlBBZ8CzRuXfvnn5bUlIS77//PjVq1MDFxQVXV1fq16/P77//brBvy5YtqVixIv+do1ZRFMqWLUvHjh0zPXZkZCRvv/02JUuWxNramsDAQKZOnUpaWppBuYweSx08eJD69etja2tLiRIlmDhxIgsXLkSlUnHjxo10x9q2bRu1atXCzs6OihUrsnjx4hf9aIiOjgZ0d0oyolb/89/MgAEDCAgISFdmypQpqFSqdOczbNgw5s2bR1BQEDY2NixcuBBPT0/eeOONdHU8fvwYOzs7Ro8eDaR/LLVhwwZUKhW7du1Kt+/cuXNRqVT6R2gnTpygV69eBAQEYGdnR0BAAL179+bmzZsv/Hlk1/MeATVr1oxmzZrpXz97xLVy5UomTJiAr68vzs7OtGrVikuXLqXbf9u2bbRs2RIXFxfs7e0JCgpi+vTpgO56/PDDDwAGj8qefTYyiik8PJzXX38dT09PbGxsCAoK4uuvvzZIYp/97L/66itmzpxJYGAgjo6O1K9fn6NHjxrUd/36dXr16oWvry82NjZ4eXnRsmVLzpw5Y/wPUYin5M6NEFkQFhYGQPny5fXbkpOTefjwIWPGjKFEiRKkpKSwc+dOunfvzpIlS+jXrx8AI0aMoGvXruzatYtWrVrp9//jjz+4du0as2fPfu5xIyMjqVu3Lmq1mkmTJlGmTBmOHDnCp59+yo0bN1iyZMlz9z137hytW7emfPny/PTTT9jb2zNv3jyWL1+eYfmzZ8/y/vvvM27cOLy8vFi4cCFvvvkmZcuWpUmTJs89Tv369QHo168fH330EY0bNzbZo5YNGzZw4MABJk2ahLe3N56enoSFhTFv3jx++OEHg0cmK1euJCkp6bl31jp16oSnpydLliyhZcuWBu8tXbqUWrVq6R+f3bhxgwoVKtCrVy9cXV2JiIhg7ty51KlTh5CQENzd3bN1PlqtNl1SammZvf+GP/roIxo2bMjChQuJjY1l7NixdO7cmdDQUCwsLABYtGgRQ4YMoWnTpsybNw9PT08uX77MhQsXAJg4cSJPnjxh7dq1HDlyRF/38xLV+/fv06BBA1JSUvjkk08ICAhg8+bNjBkzhmvXrjFnzhyD8j/88AMVK1Zk1qxZ+uN16NCBsLAwXFxcAOjQoQMajYYvv/ySUqVK8eDBAw4fPszjx4+z9XMRAgBFCKG3ZMkSBVCOHj2qpKamKnFxccq2bdsUb29vpUmTJkpqaupz901LS1NSU1OVN998U6lZs6Z+u0ajUUqXLq107drVoHz79u2VMmXKKFqtVr8NUCZPnqx//fbbbyuOjo7KzZs3Dfb96quvFEC5ePHic/d97bXXFAcHB+X+/fsGsVSqVEkBlLCwMP12f39/xdbW1uA4iYmJiqurq/L2228/95yfmTZtmmJtba0ACqAEBgYqQ4cOVc6ePWtQrn///oq/v3+6/SdPnqz8978jQHFxcVEePnxosP3cuXMKoPz4448G2+vWrasEBwfrX4eFhSmAsmTJEv220aNHK3Z2dsrjx4/120JCQhRA+e677557fmlpaUp8fLzi4OCgfPvtt/rte/bsUQBlz549z9333+UyWq5cuaIoiu4a9O/fP92+TZs2VZo2bZqurg4dOhiU+/XXXxVAOXLkiKIoihIXF6c4OzsrjRo1MviM/de7776b7mf/zH9jGjdunAIof/31l0G5//u//1NUKpVy6dIlRVH++dlXrVpVSUtL05c7duyYAigrV65UFEVRHjx4oADKrFmznhufENkhj6WEyMBLL72ElZUVTk5OtGvXjuLFi/P777+n+yt7zZo1NGzYEEdHRywtLbGysmLRokWEhobqy6jVaoYNG8bmzZsJDw8H4Nq1a2zbto133nkn3eOYf9u8eTPNmzfH19eXtLQ0/dK+fXsA9u3b99x99+3bR4sWLQzuMqjVanr06JFh+Ro1alCqVCn9a1tbW8qXL5+lRzETJ04kPDycxYsX8/bbb+Po6Mi8efMIDg5m5cqVL9z/eVq0aJGuAXfVqlUJDg42uGsVGhrKsWPHGDRoUKb1DRo0iMTERFavXq3ftmTJEmxsbOjTp49+W3x8PGPHjqVs2bJYWlpiaWmJo6MjT548Mbi2xvrf//7H8ePHDRY/P79s1dWlSxeD18/uOj27XocPHyY2NvaFnzFj7N69m0qVKlG3bl2D7QMGDEBRFHbv3m2wvWPHjvq7SBnF6OrqSpkyZZgxYwYzZ87k9OnTJmujJYo2SW6EyMCyZcs4fvw4u3fv5u233yY0NJTevXsblFm3bh09evSgRIkSLF++nCNHjnD8+HEGDRpEUlKSQdlBgwZhZ2fHvHnzAN3tejs7uxd+Gd+7d49NmzZhZWVlsFSuXBmABw8ePHff6OhovLy80m3PaBtk3GvHxsaGxMTETGP8d70DBw5k3rx5nDt3jn379mFtbc2IESOytH9Gnvd4ZNCgQRw5coS///4b+CdB+e81+q/KlStTp04dfWKk0WhYvnw5Xbt2xdXVVV+uT58+fP/99wwePJg///yTY8eOcfz4cTw8PLL888hI6dKlqV27tsFiY2OTrbr+e72e1fMsvvv37wOYtAdddHR0htfE19dX/74xMT5rA9W2bVu+/PJLatWqhYeHB8OHDycuLs5kcYuiR9rcCJGBoKAgfSPi5s2bo9FoWLhwIWvXrtWPS7J8+XICAwNZvXq1wV/GycnJ6epzcXGhf//+LFy4kDFjxrBkyRL69OlDsWLFMo3D3d2datWq8dlnn2X4/rMvlYy4ubkZNIB+JjIyMtNjmkqTJk1o06YNGzZsICoqCk9PT2xtbTP8+TwvSXveHYfevXszevRoli5dymeffcbPP/9Mt27dstRNf+DAgbzzzjuEhoZy/fp1IiIiDNrpxMTEsHnzZiZPnsy4ceP025+1sTKXzH422Wnj4+HhAcDt27dzHNszbm5uREREpNt+9+5dgGzF6e/vz6JFiwC4fPkyv/76K1OmTCElJUX/x4AQxpI7N0JkwZdffknx4sWZNGmS/ra5SqXC2tra4As4MjIyXW+pZ4YPH86DBw949dVXefz4McOGDXvhcTt16sSFCxcoU6ZMur/4a9eunWly07RpU3bv3m2QOGi1WtasWZPV086Se/fuZfgoQaPRcOXKFezt7fVJXEBAAFFRUQZJV0pKCn/++adRxyxevDjdunVj2bJlbN68mcjIyBfeBXumd+/e2NrasnTpUpYuXUqJEiVo06aN/n2VSoWiKOnuqCxcuBCNRmNUnMYICAhIN+Dh5cuXM+wBlRUNGjTAxcWFefPmpeup92//vZuSmZYtWxISEsKpU6cMti9btgyVSpVuKABjlS9fno8//piqVaumO4YQxpDkRogsKF68OOPHjyc0NJRffvkF0CUely5d4p133mH37t389NNPNGrU6LmPUsqXL0+7du3Ys2cPjRo1onr16i887rRp07CysqJBgwbMnTuX3bt3s3XrVubMmUOnTp0y/at8woQJaDQaWrZsya+//sqmTZvo3LkzT548AQy7aOfEzz//TIUKFZg8eTKbN2/mwIEDrFy5kjZt2nDx4kXGjBmDtbU1AD179sTCwoJevXqxdetW1q1bR5s2bbKVNAwaNIiIiAiGDRtGyZIlDXqiZaZYsWK8/PLLLF26lI0bN9K/f3+Dn4WzszNNmjRhxowZLFy4kJ07dzJx4kQ+++yzF95py4k33niDkJAQ3nnnHXbt2sXixYvp0qWL/g6MsRwdHfn666/Zv38/rVq1YtWqVezZs4cFCxYYJNZVq1YFdO2B/vrrL06cOEFKSkqGdY4aNYoSJUrQsWNHFixYwPbt2xkxYgRz5szh//7v/wx6E2bFuXPnaNKkCd999x3btm1j9+7dfPzxx/qefkJkWx43aBYiX3nWW+r48ePp3ktMTFRKlSqllCtXTt8D5IsvvlACAgIUGxsbJSgoSFmwYEGGPX+eWbp0qQIoq1atyvB9/tPjSVEU5f79+8rw4cOVwMBAxcrKSnF1dVWCg4OVCRMmKPHx8Znue+DAAaVevXqKjY2N4u3trXzwwQfK//73PwUw6DHk7++vdOzYMV08/+2pk5GQkBDl/fffV2rXrq14eHgolpaWSvHixZWmTZsqP//8c7ryW7duVWrUqKHY2dkppUuXVr7//vvn9pZ69913n3tcjUaj+Pn5KYAyYcKEdO9n1Fvqme3bt+t7K12+fDnd+7dv31ZeeeUVpXjx4oqTk5PSrl075cKFC+l6DxnbW2rNmjXPLaPVapUvv/xSKV26tGJra6vUrl1b2b1793N7S/23rued79atW5WmTZsqDg4Oir29vVKpUiXlf//7n/795ORkZfDgwYqHh4eiUqkMetJl1IPr5s2bSp8+fRQ3NzfFyspKqVChgjJjxgxFo9Gki2XGjBnpzvPfn9N79+4pAwYMUCpWrKg4ODgojo6OSrVq1ZRvvvnGoJeVEMZSKUom9yuFECb1yiuvcPToUW7cuIGVlVWexNCmTRtu3LjB5cuX8+T4QghhbtKgWAgzS05O5tSpUxw7doz169czc+bMXEtsRo8eTc2aNfHz8+Phw4esWLGCHTt26BtwCiFEYSTJjRBmFhERQYMGDXB2dubtt9/mvffey7VjazQaJk2aRGRkJCqVikqVKvHzzz/z+uuv51oMQgiR2+SxlBBCCCEKFektJYQQQohCRZIbIYQQQhQqktwIIYQQolApcg2KtVotd+/excnJyWSTyQkhhBDCvBRFIS4uDl9f3xcOQlrkkpu7d+9mexZeIYQQQuStW7duvXBC2CKX3Dg5OQG6H46zs3MeRyOEEEKIrIiNjcXPz0//PZ6ZIpfcPHsU5ezsLMmNEEIIUcBkpUmJNCgWQgghRKEiyY0QQgghChVJboQQQghRqBS5NjdCCCGeT6PRkJqamtdhiCLK2tr6hd28s0KSGyGEECiKQmRkJI8fP87rUEQRplarCQwMxNraOkf1SHIjhBBCn9h4enpib28vg5yKXPdskN2IiAhKlSqVo8+gJDdCCFHEaTQafWLj5uaW1+GIIszDw4O7d++SlpaGlZVVtuuRBsVCCFHEPWtjY29vn8eRiKLu2eMojUaTo3okuRFCCAFkbXA0IczJVJ9BSW6EEEIIUajkaXKzf/9+OnfujK+vLyqVig0bNrxwn3379hEcHIytrS2lS5dm3rx55g9UCCGEyOcCAgKYNWtWXoeRL+RpcvPkyROqV6/O999/n6XyYWFhdOjQgcaNG3P69Gk++ugjhg8fzm+//WbmSIUQQuQ3AwYMQKVS6Rc3NzfatWvHuXPn8jq0TN24ccMg7n8vR48efeH+S5cupVixYum2Hz9+nLfeessMERsqCElUnvaWat++Pe3bt89y+Xnz5lGqVCn9DzUoKIgTJ07w1Vdf8corr5gpyixSFG7+fRI3r5I4FvMEEwxCJIQQInPt2rVjyZIlgK47+8cff0ynTp0IDw/P48hebOfOnVSuXNlgW056q3l4eOQ0pFyVkpKS4/FsnqdAfQMfOXKENm3aGGxr27YtJ06ceO6ImsnJycTGxhos5hAXE43/6pY4zq5A6lR3oqcGcvOzWlz6ui0hS9/jzt5FaO6cgdQksxxfCCGKIhsbG7y9vfH29qZGjRqMHTuWW7ducf/+fX2ZsWPHUr58eezt7SldujQTJ07Uf2fcuHEDtVrNiRMnDOr97rvv8Pf3R1EUAEJCQujQoQOOjo54eXnxxhtv8ODBA335tWvXUrVqVezs7HBzc6NVq1Y8efIk09jd3Nz0sT9bnnV/Pnv2LM2bN8fJyQlnZ2eCg4M5ceIEe/fuZeDAgcTExOjv9kyZMgVIf0dFpVIxf/58OnXqhL29PUFBQRw5coSrV6/SrFkzHBwcqF+/PteuXdPvc+3aNbp27YqXlxeOjo7UqVOHnTt36t9v1qwZN2/eZNSoUfrjP/Pbb79RuXJlbGxsCAgI4OuvvzY434CAAD799FMGDBiAi4sLQ4YMyfTnkxMFKrmJjIzEy8vLYJuXlxdpaWkGH7J/mz59Oi4uLvrFz8/PLLHFPrxPDA4AWKk0uCkP8U+9RoW4o1S6sYwSe0djsaApms98uD+jDo/WfwhXd0JKglniEUKI7FIUhYSUtDxZniUT2REfH8+KFSsoW7aswR0QJycnli5dSkhICN9++y0LFizgm2++AXRfuK1atdLf/XlmyZIl+sdeERERNG3alBo1anDixAm2bdvGvXv36NGjBwARERH07t2bQYMGERoayt69e+nevXuOzqVv376ULFmS48ePc/LkScaNG4eVlRUNGjRg1qxZODs7ExERQUREBGPGjHluPZ988gn9+vXjzJkzVKxYkT59+vD2228zfvx4fUI3bNgwg59hhw4d2LlzJ6dPn6Zt27Z07txZfyds3bp1lCxZkmnTpumPD3Dy5El69OhBr169OH/+PFOmTGHixIksXbrUIJ4ZM2ZQpUoVTp48ycSJE7P983mRAjeI33+7iT378Dyv+9j48eMZPXq0/nVsbKxZEpwSpYNgyl2eJCQQFXmHR/dvEx8dQeKDcJSoUFzjLlOOcIqr4vF4chnOXoaz80lTWZHiWwf7ai9Dle7g4G7y2IQQwhiJqRoqTfozT44dMq0t9tZZ/2ravHkzjo6OgK4dp4+PD5s3bzaYn+jjjz/WrwcEBPD++++zevVqPvzwQwAGDx7M0KFDmTlzJjY2Npw9e5YzZ86wbt06AObOnUutWrX4/PPP9fUsXrwYPz8/Ll++THx8PGlpaXTv3h1/f38Aqlat+sLYGzRokG4epZiYGCwsLAgPD+eDDz6gYsWKAJQrV05fxsXFBZVKhbe39wuPMXDgQH0SNnbsWOrXr8/EiRNp27YtACNGjGDgwIH68tWrV6d69er6159++inr169n48aNDBs2DFdXVywsLHBycjI4/syZM2nZsqU+YSlfvjwhISHMmDGDAQMG6Mu1aNEi02TMVApUcuPt7U1kZKTBtqioKCwtLZ/7nNLGxgYbG5vcCA8AB3t7AkuXI7B0OYPtGq3C5chYdl++RNSFPbhHHaaB6jwliMbyzmG4cxjttnFQpjnqaj2hYkewdsi1uIUQoiBq3rw5c+fOBeDhw4fMmTOH9u3bc+zYMX2isXbtWmbNmsXVq1f1iYizs7O+jm7dujFs2DDWr19Pr169WLx4Mc2bNycgIADQ3ZXYs2ePPon6t2vXrtGmTRtatmxJ1apVadu2LW3atOHVV1+lePHimca+evVqgoKCDLZZWFgAMHr0aAYPHszPP/9Mq1ateO211yhTpozRP59q1arp1589+fh34uXl5UVSUhKxsbE4Ozvz5MkTpk6dyubNm/UjBScmJr6wDVNoaChdu3Y12NawYUNmzZqFRqPRn1ft2rWNPofsKFDJTf369dm0aZPBtu3bt1O7du0cDdOcGyzUKoJ8XQjyrQvN6hKTkMqOkEhOnj6BY/guOqoOUUN9Xfeo6upOFGtHVLX6Qb23oXhAXocvhChC7KwsCJnWNs+ObQwHBwfKli2rfx0cHIyLiwsLFizg008/5ejRo/Tq1YupU6fStm1bXFxcWLVqlUF7EGtra9544w2WLFlC9+7d+eWXXwzarmi1Wjp37sz//ve/dMf38fHBwsKCHTt2cPjwYbZv3853333HhAkT+OuvvwgMDHxu7H5+fgax/9uUKVPo06cPW7Zs4Y8//mDy5MmsWrWKl19+2aifz7+/G5894chom1arBeCDDz7gzz//5KuvvqJs2bLY2dnx6quvkpKSkulxFEV57pOVf3NwyJ0/2vM0uYmPj+fq1av612FhYZw5cwZXV1dKlSrF+PHjuXPnDsuWLQNg6NChfP/994wePZohQ4Zw5MgRFi1axMqVK/PqFLLNxd6KV2v78WptP6Li2rP8yE2mHD1Ks5S9dFUfIjDlHhydg/LXPFQVO0L9YeBXD2QEUSGEmalUKqMeDeUnKpUKtVpNYmIiAIcOHcLf358JEyboy9y8eTPdfoMHD6ZKlSrMmTOH1NRUunfvrn+vVq1a/PbbbwQEBGBpmfHPRaVS0bBhQxo2bMikSZPw9/dn/fr1Bs0ijFW+fHnKly/PqFGj6N27N0uWLOHll1/G2to6x9MTPM+BAwcYMGCAPomKj4/nxo0bBmUyOn6lSpU4ePCgwbbDhw9Tvnx5/V2b3JSnDYpPnDhBzZo1qVmzJqC7DVezZk0mTZoE6Bpp/ftWWGBgIFu3bmXv3r3UqFGDTz75hNmzZ+d9N/Ac8nSyZXSbCqwa/zpeXaYyxOVH+qWMZZ+mGipFC6GbYHFbWNIBbh7J63CFECLfSE5OJjIyksjISEJDQ3nvvfeIj4+nc+fOAJQtW5bw8HBWrVrFtWvXmD17NuvXr09XT1BQEC+99BJjx46ld+/e2NnZ6d979913efjwIb179+bYsWNcv36d7du3M2jQIDQaDX/99Reff/45J06cIDw8nHXr1nH//v10j5z+Kzo6Wh/7syUpKYnExESGDRvG3r17uXnzJocOHeL48eP6+gICAoiPj2fXrl08ePCAhATTdUwpW7Ys69at48yZM5w9e5Y+ffro7+o8ExAQwP79+7lz546+M8/777/Prl27+OSTT7h8+TI//fQT33//fa60r8mQUsTExMQogBITE5PXoTyXRqNVfjt5S6n/+U6l1bh5yi8TuirJk90UZbKzbln+qqLcPZvXYQohConExEQlJCRESUxMzOtQjNK/f38F0C9OTk5KnTp1lLVr1xqU++CDDxQ3NzfF0dFR6dmzp/LNN98oLi4u6epbtGiRAijHjh1L997ly5eVl19+WSlWrJhiZ2enVKxYURk5cqSi1WqVkJAQpW3btoqHh4diY2OjlC9fXvnuu++eG3dYWJhB3P9eVq5cqSQnJyu9evVS/Pz8FGtra8XX11cZNmyYwfUZOnSo4ubmpgDK5MmTFUVRFH9/f+Wbb77RlwGU9evXpzvu6dOn9dv27NmjAMqjR4/0ZZo3b67Y2dkpfn5+yvfff680bdpUGTFihH6fI0eOKNWqVVNsbGyUf6cRa9euVSpVqqRYWVkppUqVUmbMmGFw3v+NLyOZfRaN+f5WPf0BFBmxsbG4uLgQExNj0KAsP0pK1bD08A1+2HMV+6QoRliuo6flPix4ejuwyivQehq4lMzbQIUQBVpSUhJhYWEEBgZia2ub1+Hkmc8++4xVq1Zx/vz5vA6lyMrss2jM93eBGuemqLG1smBo0zLs+6A5LerW4KO0wbRM/pLt6ka6Ahd+g+/rwqHZoMl4EEMhhBCZi4+P5/jx43z33XcMHz48r8MRJiDJTQHg6mDN9O5V+WVIPbTFy/BWwjt0TP6c63ZVIfUJ7JgI85tIexwhhMiGYcOG0ahRI5o2bcqgQYPyOhxhAvJYqoBJSEnj6+2XWXwoDBQtbzv/xRjVciyTH+kK1OoHbT8HG6e8DVQIUWDIYymRX8hjqSLK3tqSiZ0qsXZoA3yLOTAvtj4N4r/kSsmnPcZOLYO5DeHm4bwNVAghhMgjktwUUMH+xdkyvBGtgjyJ0jjQ+uorzPL7Fq1LKXh8U9dtfMckSEvO61CFEEKIXCXJTQFWzN6aBf1q81GHilioVcy64sHL2i+Jr9QLUODQt7CgBdy/lNehCiGEELlGkpsCTqVS8VaTMvz69kt4Odtw9r6Wppde5WqL+WDvBvcuwI/N4fzavA5VCCGEyBWS3BQSwf6u/P5uIyr7OhP9JIUO2134s+kGCGyi61H125uw9QNIy3x+ECGEEKKgk+SmEPF2seXXt+vTKsiLlDQtb6+/xfclZqA0el9X4NiPsLQDxNzO20CFEEIIM5LkppBxsLFk/hvBDG6km4n2q53XmBj/Mtpeq8DWBW4fh/lN4daxPI5UCCGEMA9JbgohC7WKjztV4pOulVGpYPnRcMac9SFt8D7wrgoJD2BpJzi3Jq9DFUKIHBkwYAAqlYqhQ4eme++dd95BpVIxYMAAfdlu3bo9t66AgABUKlW65YsvvjBT9MJcJLkpxN6oH8CsnjWwUKtYd/oOw/54SHK/LVChI2iSYd1g2DMditY4jkKIQsbPz49Vq1aRmJio35aUlMTKlSspVaqUUXVNmzaNiIgIg+W9994zdcjCzCS5KeS61ijBvNeDsbZQs+1iJENWXSLx5aXQ4On8Kfu+gN8GQ2pSnsYphBDZVatWLUqVKsW6dev029atW4efnx81a9Y0qi4nJye8vb0NFgcHB1OHLMxMkpsioHUlLxYPqIOdlQX7L9+n/08nSWg2GTrPBrUlXFgLy7tDUkxehyqEyA8UBVKe5M2SzTvJAwcOZMmSJfrXixcvlnmiijDLvA5A5I5G5dxZPrguA5Yc51jYQ4YsO8Gi/q9jWzwAVvWFm4dgaUd4fR04euZ1uEKIvJSaAJ/75s2xP7oL1sbfKXnjjTcYP348N27cQKVScejQIVatWsXevXuNqmfs2LF8/PHHBts2b95Ms2bNjI5J5B25c1OEBPu78tOgujhYW3DoajTvrDhFSqnGMHALOHhA5HlY3BYe3cjrUIUQwiju7u507NiRn376iSVLltCxY0fc3d2NrueDDz7gzJkzBku9evXMELEwJ7lzU8TUKlWcxQPq0H/JMXb/HcXwlaf5vk9NLAf9CT93g4fXYVFbeGMdeFXO63CFEHnByl53ByWvjp1NgwYNYtiwYQD88MMP2arD3d2dsmXLZjsGkT/InZsiqF5pNxb0q421pa6R8ehfz6IpXhoGbQfPShAfCUvaw51TeR2qECIvqFS6R0N5sahU2Q67Xbt2pKSkkJKSQtu2bU34AxEFjdy5KaIal/Ngbt9avP3zSTaevYujrSWfdauCauBWWNEDbh+DZd3gjfVQMjivwxVCiBeysLAgNDRUv56RmJgYzpw5Y7DN1dVV32U8Li6OyMhIg/ft7e1xdnY2fcDCbOTOTRHWMsiL2b1rolLBL3+F88Oeq2BXXPdIqlR9SI7RPaq6dTyvQxVCiCxxdnbONBHZu3cvNWvWNFgmTZqkf3/SpEn4+PgYLB9++GFuhC5MSKUoRWsEt9jYWFxcXIiJiZFM/KmfDt9g8saLAMx4tRqv1faD5Hj4pYeuF5W1E7z+G5SSRnVCFEZJSUmEhYURGBiIra1tXocjirDMPovGfH/LnRtB/wYB/F+zMgCMW3eePZeiwMYR+q6BgMaQEqcbByf8aB5HKoQQQryYJDcCgA/bVqB7zRJotArvrjjFuduPdY37+vwKgU0gJR5WvAZ3T+d1qEIIIUSmJLkRALrJ4V6pRqOy7iSkaBi09Di3HyWAtT30Xg3+DSE5Fn7uDlGheR2uEEII8VyS3Ag9a0s1c1+vRZCPMw/iUxj80wmeJKc9TXBWgW8tSHyo60X18HpehyuEEEJkSJIbYcDJ1oqF/Wvj7mjD35FxjFp9Bq1WAVtnXaPiZ+Pg/NQVYu7kdbhCCBMqYv1LRD5kqs+gJDcinRLF7Jj/hm4m8e0h9/h6xyXdG/au8MYGcC0DMeGwrCs8eZCnsQohcs7KygqAhISEPI5EFHUpKSnA88cpyioZxE9kKNi/OP97tSqjVp/lhz3XKOfpRLeaJcDJC/r9rhvBOPqKrrt4/03ZmuhOCJE/WFhYUKxYMaKiogDdoHWqHIwULER2aLVa7t+/j729PZaWOUtPJLkRz/VyzZJcvhfP3L3X+PC3c/i72VOzVHEo5qebPXxxG7hzEtYMhF6/gIV8nIQoqLy9vQH0CY4QeUGtVlOqVKkcJ9cyiJ/IlFar8NbPJ9kZeg9vZ1s2D2+Eu6ON7s3wv2BZF0hLgpqvQ5fvczQvjBAi72k0GlJTU/M6DFFEWVtbo1Zn3GLGmO9vSW7EC8Unp9H1+4Ncu/+E+qXd+PnNulhaPP3w/b0VVvcFRQtNPoAWH+dtsEIIIQolGaFYmJSjjSXz3wjGwdqCI9ejmbH90j9vVuwAnb7Rre+fAccX5U2QQgghxFOS3IgsKevpxIzXqgMwf991/jgf8c+bwQOg6Tjd+tYP4OrO3A9QCCGEeEqSG5FlHar68FaT0gB8sPYcV6Pi/3mz2Tio3gcUja6BsYxiLIQQIo9IciOM8mHbCrxU2pX45DSGLj+pG8EYdA2JO8+CUg100zSs6AHx0utCCCFE7pPkRhjF0kLNd71r4eVsw9WoeCb9fvFfb9pArxXgWlo3yN+qPpCamHfBCiGEKJIkuRFG83CyYXavmqhV8Nup26w9efufN+1ddTOJ27rA7eOw4R3QavMuWCGEEEWOJDciW+qVdmNkq/IATNxwwbD9jXs56Lkc1JZwcR0c+CqPohRCCFEUSXIjsu3d5mVpUMaNxFQNw345RVKq5p83A5tAx5m69T2fwd9b8iZIIYQQRY4kNyLbLNQqZvWsgZuDNX9HxvHJ5hDDAsH9oe5buvV1b0kPKiGEELlCkhuRI57OtnzTswYAK/4KZ8u5CMMCbT+HgMaQEg8re0PCw9wPUgghRJEiyY3IsSblPfi/ZmUAGL/uHHcf/6uHlIUVvPYTuJSCR2GwdhBo0vIoUiGEEEWBJDfCJEa3Lk91v2LEJqUx+tczaLT/mrLMwQ16/wJW9nB9D+ycnHeBCiGEKPQkuREmYWWhZlbPGthbW3D0+kMWHLhuWMC7KnSbq1s/8j1cXJ/7QQohhCgSJLkRJhPo7sCkTpUA+Hr7JS7ciTEsULkbNBypW9/wLty/hBBCCGFqktwIk+pZx482lbxI1SiMXH2GxBSNYYEWE3XdxFOfwOrXITkubwIVQghRaElyI0xKpVLxxSvV8HDSTc8w/Y//dP+2sIRXFoOTLzy4DL+/C4qScWVCCCFENkhyI0zO1cGar16rDsCyIzfZe+k/E2g6ekCPZaC2gpDfdW1whBBCCBOR5EaYRdPyHgxoEADA2N/OEZOQaljArw60/0K3vmMy3DiYuwEKIYQotCS5EWYztl1FAt0duBebzNRNF9MXqP0mVOsFigbWvgnx93M/SCGEEIWOJDfCbOysLfjqteqoVbDu9B3+vBhpWEClgk4zwaMixEfCuiGg1WRcmRBCCJFFktwIswr2L85bTXSjF09Yf57o+GTDAtYOuhGMnw3wd+DrPIhSCCFEYSLJjTC7Ua3LUd7LkQfxKXy84QLKf3tHeVaEjk+Tmr3TIWx/7gcphBCi0JDkRpidjaUFM3vUwFKt4o8LkWw8ezd9oRp9oMbroGjht8EQH5W+jBBCCJEFktyIXFGlhAvDWpQFYPLGi9yPS05fqMMM8AiC+Hu6BEfa3wghhMiGPE9u5syZQ2BgILa2tgQHB3PgwIFMy69YsYLq1atjb2+Pj48PAwcOJDo6OpeiFTnxbvOyVPJx5nFCKpM3XkhfwNoeejxtfxO2Dw5+k/tBCiGEKPDyNLlZvXo1I0eOZMKECZw+fZrGjRvTvn17wsPDMyx/8OBB+vXrx5tvvsnFixdZs2YNx48fZ/DgwbkcucgOKws1X75aDQu1iq3nI9l2ISJ9IY8Kujs4AHs+h1vHcjdIIYQQBV6eJjczZ87kzTffZPDgwQQFBTFr1iz8/PyYO3duhuWPHj1KQEAAw4cPJzAwkEaNGvH2229z4sSJXI5cZFeVEi4MbVoagI83XORxQkr6QjX6QpVX/xn/JvFx7gYphBCiQMuz5CYlJYWTJ0/Spk0bg+1t2rTh8OHDGe7ToEEDbt++zdatW1EUhXv37rF27Vo6duz43OMkJycTGxtrsIi89V6LcpTxcOBBfDLTNoekL/Bs/Jti/hATDptGyPxTQgghsizPkpsHDx6g0Wjw8vIy2O7l5UVkZGSG+zRo0IAVK1bQs2dPrK2t8fb2plixYnz33XfPPc706dNxcXHRL35+fiY9D2E8WysLvny1OioVrDt1hz3/nXsKwNYFXl0CaksI2QCnluV6nEIIIQqmPG9QrFKpDF4ripJu2zMhISEMHz6cSZMmcfLkSbZt20ZYWBhDhw59bv3jx48nJiZGv9y6dcuk8YvsCfYvzsAGgQB8tO48cUmp6QuVDIYWE3Xrf4yFqL9zMUIhhBAFVZ4lN+7u7lhYWKS7SxMVFZXubs4z06dPp2HDhnzwwQdUq1aNtm3bMmfOHBYvXkxERAaNUwEbGxucnZ0NFpE/jGlbnlKu9kTEJDHjz0sZF2owHEo3h7REXffwtAy6kAshhBD/kmfJjbW1NcHBwezYscNg+44dO2jQoEGG+yQkJKBWG4ZsYWEBkH7UW5Hv2Vtb8kX3qgD8fPQmJ28+Sl9IrYaX54O9G9w7D7s/yeUohRBCFDR5+lhq9OjRLFy4kMWLFxMaGsqoUaMIDw/XP2YaP348/fr105fv3Lkz69atY+7cuVy/fp1Dhw4xfPhw6tati6+vb16dhsiBBmXdeTW4JIoC49edIyVNm76Qkxd0+V63fvg7uL43V2MUQghRsORpctOzZ09mzZrFtGnTqFGjBvv372fr1q34+/sDEBERYTDmzYABA5g5cybff/89VapU4bXXXqNChQqsW7cur05BmMCEDkG4OVhz+V488/ddy7hQxQ4QPFC3vv7/IOFh7gUohBCiQFEpRex5TmxsLC4uLsTExEj7m3zk9zN3GLHqDNYWav4Y2ZgyHo7pC6U8gflNIPoqBHWBHst03caFEEIUesZ8f+d5bykhALpU96VJeQ9SNFo+WncerTaDnNvaAV5ZqOseHroRTi/P/UCFEELke5LciHxBpVLxWbcq2FlZ8FfYQ9acfE6Xfd+a0OJj3fofY+Hh9dwLUgghRIEgyY3IN/xc7RndujwAn20J5UH8c7p9NxgO/o0g9Ymu/Y3MHi6EEOJfJLkR+crAhgFU9nUmNimNz7eEZlxIbQEvzwVrJ7h1FA59m7tBCiGEyNckuRH5iqWFms9erqqbmuH0HQ5fe5BxwWKloMOXuvU9n0PEudwLUgghRL4myY3Id2r4FeP1errhAD5ef4HktOc8dqreGyp2Am0qrHsLUpNyMUohhBD5lSQ3Il/6oF0FPJxsuP7gCfP3PafRsEoFnb8FB0+4HyqjFwshhAAkuRH5lLOtFRM7VQLg+z1XufHgScYFHdyhy9NZ4Y/8AGEHcilCIYQQ+ZUkNyLf6lzNh8bl3ElJ0zLx9wvPnz+sQjuo1R9QYMM7kByXq3EKIYTIXyS5EfmWSqXik65VsLZUc+DKAzady3jmdwDafgbF/CEmHP6ckHtBCiGEyHckuRH5WoC7A8OalwXg080hxCWlZlzQxgm6zdGtn/oJruzMpQiFEELkN5LciHzvrSalCXCzJyoumVk7rzy/YEAjeOkd3frGYZD4KHcCFEIIka9IciPyPVsrC6Z0qQzA0sM3+Dsy9vmFW04Ct7IQF6GbnkEIIUSRI8mNKBCaVfCkXWVvNFqFSRsuPr9xsZUddJsHKjWcWw0hG3M3UCGEEHlOkhtRYEzqXAk7KwuO3XjIulN3nl/Qrw40HKlb3zwKnjxnlGMhhBCFkiQ3osDwLWbH8JblAJj+Rygxic9pXAzQbBx4VoKEB7Dl/VyKUAghRH4gyY0oUN5sFEhZT0cexKfw9fZLzy9oaaPrPaWygJANcHF9rsUohBAib0lyIwoUa0s107rqGhcvP3qTC3dinl/YtyY0Hq1b3/I+xN/PhQiFEELkNUluRIHToIw7nav7olVg8sZMGhcDNPkQPCtDQjRsHZN7QQohhMgzktyIAumjDhWxt7bg5M1HrD+dSeNiS2t5PCWEEEWMJDeiQPJxseO9FrrGxZ9v/fv5IxcD+NaQx1NCCFGESHIjCqxBjQIo7e7Ag/hkvs1s5GL4z+Mp6T0lhBCFmSQ3osCysbRg8tORi5ccvsHle5nMBm7weOp33SKEEKJQkuRGFGhNy3vQtrIXGq3C5N9f0LjYtwY0GqVb3/I+JDzMlRiFEELkLkluRIH3ccdK2FiqOXI9mi3nIzIv3PRD8KgIT+7L3FNCCFFISXIjCjw/V3veaVYWgM+3hJKQkvb8wpY20PUH3dxT53+FS3/kUpRCCCFyiyQ3olB4u2lpShSz425MEnP3Xsu8cMnaUP9d3frmUZD42OzxCSGEyD2S3IhCwdbKgomdggCYv/864dEJme/QfAK4loG4CNg+IRciFEIIkVuMTm4CAgKYNm0a4eHh5ohHiGxrW9mbhmXdSEnT8umWkMwLW9npHk+hgtPL4equXIlRCCGE+Rmd3Lz//vv8/vvvlC5dmtatW7Nq1SqSk5PNEZsQRlGpVEzpXBkLtYrtIffYf/kFg/X514e6b+nWN42E5HizxyiEEML8jE5u3nvvPU6ePMnJkyepVKkSw4cPx8fHh2HDhnHq1ClzxChElpXzcqJ//QAApm66SKpGm/kOLSeBSymICYdd08wfoBBCCLPLdpub6tWr8+2333Lnzh0mT57MwoULqVOnDtWrV2fx4sWZjzcihBmNaFUONwdrrt1/wk+Hb2Re2MYROs/SrR/7EcKPmjs8IYQQZpbt5CY1NZVff/2VLl268P7771O7dm0WLlxIjx49mDBhAn379jVlnEJkmYudFR+2qwDAtzuvcD/uBY9Ny7aEGq8DCvw+DFKTzB+kEEIIszE6uTl16hTvvfcePj4+vPfee1SuXJkLFy5w8OBBBg4cyIQJE9i4cSPr18vsyyLvvBbsR9USLsQlp/HVn5devEPbT8HRC6KvwL7/mT9AIYQQZmN0clOnTh2uXLnC3LlzuX37Nl999RUVK1Y0KFOpUiV69eplsiCFMJZarWJKl0oA/HryFudvx2S+g11x6DhTt37oW4g4a+YIhRBCmItKMbJxzM2bN/H39zdXPGYXGxuLi4sLMTExODs753U4wsxGrjrNhjN3qe1fnDVD66NSqTLf4df+ELIBvKvCkD1gYZUrcQohhMicMd/fRt+5ad68OdHR0em2P378mNKlSxtbnRBmNbZ9ReysLDhx8xEbz9598Q4dvtLdxYk8D4e/M3+AQgghTM7o5ObGjRtoNJp025OTk7lz545JghLCVHxc7Hi3eRkAvvjj78znnQJw9IC203Xre7+AB1fNHKEQQghTs8xqwY0bN+rX//zzT1xcXPSvNRoNu3btIiAgwKTBCWEKgxuXZtXxW9x+lMi8vdcY3aZC5jtU76WbVPPabtg0HPpvBrXMVCKEEAVFltvcqJ/+565SqdKNYWNlZUVAQABff/01nTp1Mn2UJiRtboqmbRciGLr8FDaWanaOboqfq33mOzy6CXPqQ+oTXUPjOm/mTqBCCCEyZJY2N1qtFq1WS6lSpYiKitK/1mq1JCcnc+nSpXyf2Iiiq21lb+qXdiM5Tcv0P0JfvENxf93oxQA7JkOMPHIVQoiCwuh77WFhYbi7u5sjFiHMRqVSMblLJdQq2Ho+kqPX0zeKT6fuEChZB1LiYMv7IKNuCyFEgZClNjezZ8/mrbfewtbWltmzZ2dadvjw4SYJTAhTq+jtTJ96pVh+NJxpm0LY9F4jLNSZdA1XW0CX72BeY7j8B1xcB1Veyb2AhRBCZEuW2twEBgZy4sQJ3NzcCAwMfH5lKhXXr183aYCmJm1uiraHT1JoNmMPsUlpTO9eld51S714p71fwN7pYO8Ow46Dvav5AxVCCGHAmO9vowfxK+gkuRGLD4YxbXMIbg7W7PmgGc62LxioLy0F5jeB+6FQvQ+8PDd3AhVCCKFn1kH8hCjo3qjvTxkPB6KfpPDdrisv3sHSGrrMBlRw9hddF3EhhBD5Vpbu3IwePTrLFc6cOTNHAZmb3LkRAHsvRTFgyXEs1Sq2j2pCaQ/HF++09UM4Nh+KlYJ3joK1g/kDFUIIARj3/Z2lBsWnT5/O0oFfOG+PEPlEswqeNK/gwZ5L9/lsSyiLBtR58U4tJ8KlrfA4HPZ8Dm0/M3+gQgghjCZtbkSRde1+PG2/2U+aVuGnQXVpWt7jxTtd2QErXgWVGgbvhBLB5g9UCCGEtLkRIivKeDjSv0EAAJ9sDiFVo33xTuVaQ9XXQNHCxuGgSTVvkEIIIYyWpcdS3bt3Z+nSpTg7O9O9e/dMy65bt84kgQmRG4a3KMe6U7e5GhXPiqM3GdDw+UMd6LX7Aq7ugnsX4PBsaPy++QMVQgiRZVm6c+Pi4qJvT+Pi4pLpIkRB4mJvpZ9I85udV3j0JOXFOzm46xIcgL3/k5nDhRAin5E2N6LIS9No6Tj7IJfuxdG/vj9Tu1Z58U6KAstfgWu7wL8R9N8kM4cLIYQZ5Uqbm6ioKA4cOMDBgweJiorKbjVC5DlLCzWTOlcCYPlf4Vy+F/finVQq6PQNWNnDzYNwepmZoxRCCJFVRic3sbGxvPHGG5QoUYKmTZvSpEkTSpQoweuvv05MTIw5YhTC7BqWdad1JS80WoVPNoeQpRuaxf2hxUTd+vZJEBth3iCFEEJkidHJzeDBg/nrr7/YvHkzjx8/JiYmhs2bN3PixAmGDBlijhiFyBUTOgRhZaHiwJUH7LmUxbuR9d4G31qQHAN/fGDeAIUQQmSJ0cnNli1bWLx4MW3btsXZ2RknJyfatm3LggUL2LJli9EBzJkzh8DAQGxtbQkODubAgQOZlk9OTmbChAn4+/tjY2NDmTJlWLx4sdHHFeK/AtwdGPS0t9Snm0NJSctC1/BnM4erLSF0k24RQgiRp4xObtzc3DLsFeXi4kLx4sWNqmv16tWMHDmSCRMmcPr0aRo3bkz79u0JDw9/7j49evRg165dLFq0iEuXLrFy5UoqVqxo7GkIkaFhLcri7mjN9QdPWHbkRtZ28q4CDUfo1reMgcTH5gpPCCFEFhjdW+rHH39kzZo1LFu2DB8fHwAiIyPp378/3bt35+23385yXfXq1aNWrVrMnfvPLMtBQUF069aN6dOnpyu/bds2evXqxfXr13F1dTUmbD3pLSVeZNWxcMatO4+TrSV7xzTDzdHmxTulJsG8hhB9FYIHQudZZo9TCCGKEpP3lqpZsya1atWiVq1azJs3j6NHj+Lv70/ZsmUpW7YspUqV4vDhw8yfPz/LQaakpHDy5EnatGljsL1NmzYcPnw4w302btxI7dq1+fLLLylRogTly5dnzJgxJCYmPvc4ycnJxMbGGixCZOa12n4E+TgTl5TGzB2Xs7aTlS10/la3fnIJ3DhkvgCFEEJkKksjFHfr1s3kB37w4AEajQYvLy+D7V5eXkRGRma4z/Xr1zl48CC2trasX7+eBw8e8M477/Dw4cPntruZPn06U6dONXn8ovCyUKuY3LkSvX48yspj4bz+kj9BPlm4yxfQCIIHwMmlsGk4DD2kS3qEEELkqiwlN5MnTzZbAP+dSVxRlOfOLq7ValGpVKxYsULf7mfmzJm8+uqr/PDDD9jZ2aXbZ/z48YwePVr/OjY2Fj8/PxOegSiMXirtRvsq3vxxIZJPt4Sw/M16WZv1vtVUuLRN93jqwFfQ4mPzByuEEMJAng2p6u7ujoWFRbq7NFFRUenu5jzj4+NDiRIlDBo0BwUFoSgKt2/fznAfGxsbnJ2dDRYhsuKjDkFYW6o5dDWaHSH3sraTXTHoMEO3fvAbuHfRbPEJIYTImNHJjUaj4auvvqJu3bp4e3vj6upqsGSVtbU1wcHB7Nixw2D7jh07aNCgQYb7NGzYkLt37xIfH6/fdvnyZdRqNSVLljT2VITIlJ+rPYMb6bqGf7Y1lOQ0TdZ2rNQFKnYCbRpsfA+0WdxPCCGESRid3EydOpWZM2fSo0cPYmJiGD16NN27d0etVjNlyhSj6ho9ejQLFy5k8eLFhIaGMmrUKMLDwxk6dCige6TUr18/ffk+ffrg5ubGwIEDCQkJYf/+/XzwwQcMGjQow0dSQuTUO83L4uFkw83oBJYeupH1HTt8BTbOcOckHPvRbPEJIYRIz+jkZsWKFSxYsIAxY8ZgaWlJ7969WbhwIZMmTeLo0aNG1dWzZ09mzZrFtGnTqFGjBvv372fr1q34+/sDEBERYTDmjaOjIzt27ODx48fUrl2bvn370rlzZ2bPnm3saQiRJY42lnzQVjdr+He7r3I/LjlrOzr7QOtpuvVdn8Cjm2aKUAghxH8ZPc6Ng4MDoaGhlCpVCh8fH7Zs2UKtWrW4fv06NWvWzPfzS8k4N8JYWq1C1x8Ocf5ODL3q+PHFK9WyuiP81AluHoIyLeH133QTbgohhDCaWWcFL1myJBERugkCy5Yty/bt2wE4fvw4NjZZGOxMiAJG/bRrOMDqE7e4cCeLCbxarRv7xsIGru2Cc7+aMUohhBDPGJ3cvPzyy+zatQuAESNGMHHiRMqVK0e/fv0YNGiQyQMUIj+oHeBK5+q+KApM25TFWcMB3MtB0w9169vGwZMH5gtSCCEEkI3HUv919OhRDh8+TNmyZenSpYup4jIbeSwlsuvO40Rafr2XpFQtP/SpRcdqPlnbUZMKPzaDexegag94ZYFZ4xRCiMLImO/vHCc3BY0kNyInZu64zOxdVyhRzI5d7zfF1soiazveOQkLW4Gihb5roVxr8wYqhBCFjFnb3ABcunSJYcOG0bJlS1q1asWwYcO4dOlStoIVoiAZ2rQ0Pi623HmcyMID17O+Y4lgeOkd3fqmkZAcZ5b4hBBCZCO5Wbt2LVWqVOHkyZNUr16datWqcerUKapUqcKaNWvMEaMQ+Ya9tSXj2lcE4Ic914iMScr6zs0/gmL+EHtb1z1cCCGEWRj9WKp06dK8/vrrTJs2zWD75MmT+fnnn7l+3Yi/ZvOAPJYSOaUoCq/MPcyp8Md0r1mCmT1rZH3na3vg526ACt7cDn51zRSlEEIULmZ9LBUZGWkwavAzr7/++nNn8xaiMFGpVEzpUhmAdafvcCr8UdZ3LtMcavQFFPh9GKRlcVBAIYQQWWZ0ctOsWTMOHDiQbvvBgwdp3LixSYISIr+rVrIYrwXr5jObuikErdaIG6BtPgUHD3hwCQ7MNFOEQghRdFlmpdDGjRv16126dGHs2LGcPHmSl156CdB1B1+zZg1Tp041T5RC5EMftKvA1vMRnL31mPWn7/BKcBYnb7V3hfZfwtqBcOBrqNwNPIPMGqsQQhQlWWpzo1Zn7QaPSqVCo8nfMyBLmxthSvP2XeOLP/7G08mG3WOa4WiTpb8XQFFgVR+4tBVK1Na1v1FnsVu5EEIUQSZvc6PVarO05PfERghTG9gwAH83e6Likpmz52rWd1SpoOPXT2cOPyEzhwshhAlla5wbIYSOjaUFH3fUzTu18EAYN6OfZH1nZ99/zRw+DR7dMH2AQghRBGUrudm3bx+dO3embNmylCtXji5dumTYyFiIoqBVkCeNy7mTotHy2ZZQ43au1R8CGkNqAmwaoXtcJYQQIkeMTm6WL19Oq1atsLe3Z/jw4QwbNgw7OztatmzJL7/8Yo4YhcjXVCoVkzpVwkKtYnvIPQ5eMWJyzGczh1vawvW9cGaF2eIUQoiiwuhB/IKCgnjrrbcYNWqUwfaZM2eyYMECQkON/Ms1l0mDYmEuUzZeZOnhG5TzdGTriMZYWRjxt8Oh2bBjIti6wLvHwMnbfIEKIUQBZNZB/K5fv07nzp3Tbe/SpQthYWHGVidEoTGqVXlcHay5EhXP8qM3jdv5pXfAtyYkxcDWMeYJUAghigijkxs/Pz927dqVbvuuXbvw8/MzSVBCFEQu9laMaVMBgG92XCY63ojRhy0socv3oLaE0E1wcYN5ghRCiCIgi4Ny/OP9999n+PDhnDlzhgYNGqBSqTh48CBLly7l22+/NUeMQhQYPev4seKvm1y8G8tX2y8zvXvVrO/sXQUajYb9X8LWDyCwiW7APyGEEEYxus0NwPr16/n666/17WuCgoL44IMP6Nq1q8kDNDVpcyPM7VjYQ3rMP4JKBZuGNaJKCZes75yWDPObwP2/oVov6D7ffIEKIUQBYrY2N2lpaUydOpXatWtz8OBBoqOjiY6O5uDBgwUisREiN9QNdKVzdV8UBaZuuohRfz9Y2ugeT6GCc6vg8nazxSmEEIWVUcmNpaUlM2bMkJGIhXiB8e0rYmdlwfEbj9h49q5xO/vVgfrv6tY3j4SkWJPHJ4QQhZnRDYpbtWrF3r17zRCKEIWHbzE73m1eBoDPt4byJDnNuAqaT4DigRB7B3ZMMkOEQghReBndoLh9+/aMHz+eCxcuEBwcjIODg8H7Xbp0MVlwQhRkgxuX5tcTtwl/mMD3e64ytl3FrO9sbQ9dvoOfOsHJJVClu66BsRBCiBcyukFxZjOEy6zgQhjaEXKPIctOYGWhYvuopgS6O7x4p3/bPApOLIbiAfB/h8HayP2FEKKQMOsgfjIruBBZ1yrIk6blPUjVKEzbdDEbFUwF55K6STV3fWLy+IQQojAyKrm5efMmCxYsYO7cuYSEhJgrJiEKDZVKxeTOlbCyULHn0n12hd4zrgJbZ+jydPyov+ZB+FHTBymEEIVMlpOb/fv3U7lyZd5++23effddatSowcqVK80ZmxCFQmkPRwY1CgRg2uYQklKNvMNZthXUeB1Q4Pd3ITXR9EEKIUQhkuXkZuLEiTRv3pzbt28THR3NoEGD+PDDD80ZmxCFxnstyuHpZMPN6AQWHczGHGxtPwMnH4i+Cns+M32AQghRiGQ5uTl//jzTp0/H19eX4sWL8/XXX3P37l0ePXpkzviEKBQcbSz5qEMQAN/tvsKdx0befbErBp2fPp468gPcOm7aAIUQohDJcnLz+PFjPD099a8dHBywt7fn8ePH5ohLiEKnaw1f6ga4kpSq5bMt2WizVr6tbkoGRQu/vwOpSaYPUgghCgGjGhSHhIRw7tw5/aIoCqGhoQbbhBAZU6lUTO1aGQu1iq3nIzl45YHxlbSbDo5e8OAy7PvC9EEKIUQhkOVxbtRqNSqVKsN5cp5tl3FuhHixKRsvsvTwDUp7OLBtRBOsLY0ckSF0M6zuCyo1vLkTSgabJ1AhhMhHjPn+zvIIxWFh2WgEKYRIZ1Tr8mw+d5fr95+w+FAYQ5uWMa6CoE5Q9TU4vwY2/B+8vR+sbM0TrBBCFEBGj1Bc0MmdG5EfrD15mzFrzmJvbcGu95vi42JnXAUJD2HOSxB/DxqOgNbTzBOoEELkE2YdoVgIkXPda5agVqliJKRo+GxLqPEV2LtCp1m69cPfSe8pIYT4F0luhMgDarWKaV2roFbB5nMRHL6ajcbFFTv803tqw//J4H5CCPGUJDdC5JEqJVx4/SV/ACb+foGUNK3xlbT/Ahy9IfqKDO4nhBBPSXIjRB56v00F3B2tuXb/CQsPXje+Arvi/wzud/h7mXtKCCGQ5EaIPOViZ/XPyMW7rnL7UYLxlVRoB9X7AIru8VTKE9MGKYQQBUyWuoLXrFkTlUqVpQpPnTqVo4CEKGperlmCVcdvcSzsIdM2hfBjv9rGV9JuOoTtg4fXYecU6DDD5HEKIURBkaXkplu3bmYOQ4iiS6VS8UnXKnScfYDtIffY83cUzSt6vnjHf7MrBl2+g+Xd4diPULEjlG5mjnCFECLfk3FuhMgnPt8ayo/7r1PK1Z7to5pga2VhfCWbR8OJReDiB/93CGxdTB+oEELkARnnRogCaETLcng72xL+MIE5e69lr5LW06B4AMTcgm0fmTQ+IYQoKIxObjQaDV999RV169bF29sbV1dXg0UIkT0ONpZM6lwJgHl7r3H9frzxldg4Qre5gArOLIdLf5g2SCGEKACMTm6mTp3KzJkz6dGjBzExMYwePZru3bujVquZMmWKGUIUouhoX8WbpuU9SNFomfj7hQwnqn0h/wZQ/13d+sbh8CTatEEKIUQ+Z3Rys2LFChYsWMCYMWOwtLSkd+/eLFy4kEmTJnH0qIyxIUROPGtcbGOp5tDVaDaevZu9ilpMBI+K8CQKNo+EotW0TghRxBmd3ERGRlK1alUAHB0diYmJAaBTp05s2bLFtNEJUQSVcrNneMtyAHyyOYSYhFTjK7Gyhe4/gtoSQjfCudUmjlIIIfIvo5ObkiVLEhERAUDZsmXZvn07AMePH8fGxsa00QlRRA1pXJqyno48iE/hyz//zl4lPtWh2Tjd+tYP4PEt0wUohBD5mNHJzcsvv8yuXbsAGDFiBBMnTqRcuXL069ePQYMGmTxAIYoia0s1n3arAsAvx8I5Ff4oexU1HAUl60ByrG70Ym025q8SQogCJsfj3Pz1118cOnSIsmXL0qVLF1PFZTYyzo0oSN7/9Sy/nbpNkI8zm4Y1xNIiG6M3RF+DeY0gNQHafv5PY2MhhChAjPn+Njq5SUhIwN7ePkcB5iVJbkRBEh2fTIuv9xGTmMqEDkEMaVI6exWdWAybR4GFDby9DzyDTBuoEEKYmVkH8fP09OT111/nzz//RCu3uIUwKzdHGyY8nVhz5o7L3HqYjYk1AYIHQrk2oEmGdUMgLdmEUQohRP5idHKzbNkykpOTefnll/H19WXEiBEcP37cHLEJIYDXapfkpdKuJKZq+HhDNse+Uamgy/dg7waR52H3p6YPVAgh8gmjk5vu3buzZs0a7t27x/Tp0wkNDaVBgwaUL1+eadOmmSNGIYo0lUrFZy9XxdpCzb7L99l0LiJ7FTl56SbXBDj8HYTtN12QQgiRj2R7biknJycGDhzI9u3bOXv2LA4ODkydOtWUsQkhnirj4ciwFmUBmLbpIo8TUrJXUcWOUKs/oMD6oZCYzV5YQgiRj2U7uUlKSuLXX3+lW7du1KpVi+joaMaMGWN0PXPmzCEwMBBbW1uCg4M5cOBAlvY7dOgQlpaW1KhRw+hjClEQDW1ahnJPx7754o9sjn0Duh5TrqUh9o5uFnEZvVgIUcgYndxs376d/v374+XlxdChQ/H09OTPP/8kPDyc//3vf0bVtXr1akaOHMmECRM4ffo0jRs3pn379oSHh2e6X0xMDP369aNly5bGhi9EgWVtqebz7rrRwVcdv8XR69mcM8rGEbovBJUFXFwH5341YZRCCJH3jO4Kbm9vT8eOHenbty8dO3bEysoq2wevV68etWrVYu7cufptQUFBdOvWjenTpz93v169elGuXDksLCzYsGEDZ86cyfIxpSu4KOg+Wn+eX/4Kp7S7A1tHNMbWyiJ7Fe37EvZ8BjbOMPQAFA8waZxCCGFKZu0KHhkZyZo1a+jWrVuOEpuUlBROnjxJmzZtDLa3adOGw4cPP3e/JUuWcO3aNSZPnpztYwtRkI1tVxFPJxuuP3jCd7uvZL+iRqPBr55u9OLfhoAmzXRBCiFEHspSchMbG5vu9fOWrHrw4AEajQYvLy+D7V5eXkRGRma4z5UrVxg3bhwrVqzA0tIyS8dJTk7OdoxC5EcudlZ88nRqhnn7rnPxbkz2KrKwhO4LwMYFbh+DfcY9VhZCiPwqS8lN8eLFiYqKAqBYsWIUL1483fJsu7FUKpXBa0VR0m0D0Gg09OnTh6lTp1K+fPks1z99+nRcXFz0i5+fn9ExCpHftK3sTYeq3mi0CmN/O0eaJpsDahb3h04zdesHvoIbh0wXpBBC5JEs3f7YvXs3rq6u+vWMkg9jubu7Y2Fhke4uTVRUVLq7OQBxcXGcOHGC06dPM2zYMAC0Wi2KomBpacn27dtp0aJFuv3Gjx/P6NGj9a9jY2MlwRGFwpQulTl0NZoLd2JZdDCMt5uWyV5FVV+Fa7vhzApY9xb830GwM/4PFSGEyC9yPHFmTtSrV4/g4GDmzJmj31apUiW6du2arkGxVqslJCTEYNucOXPYvXs3a9euJTAwEAcHhxceUxoUi8JkzYlbfLD2HDaWaraNbEKg+4t/BzKUHA/zm8DDa1CpK7z2k25UYyGEyCfM2qC4dOnSTJw4kUuXLmU7wGdGjx7NwoULWbx4MaGhoYwaNYrw8HCGDh0K6O669OvXTxeoWk2VKlUMFk9PT2xtbalSpUqWEhshCptXg0vSuJw7yWlaxq87l72pGUDXPfyVhaC2gpDf4dRPpg1UCCFykdHJzbBhw9i2bRtBQUEEBwcza9YsIiKyNxx8z549mTVrFtOmTaNGjRrs37+frVu34u/vD0BERMQLx7wRoihTqVR8/nJV7KwsOHr9Ib8cy8HvS4la0HKibv2PcRAVapoghRAil2X7sdTly5dZsWIFq1at4vr16zRv3pzXX39df6clv5LHUqIwWnQwjE82h+BgbcGfo5pQsrh99irSamHFK7o2OB5BMGQ3WGezLiGEMCGzPpZ6pnz58kydOpVLly5x4MAB7t+/z8CBA7NbnRAiBwY0CKC2f3GepGgYv+589h9PqdXw8nxw9IL7ofDneNMGKoQQuSDbyQ3AsWPHGDlyJC+//DKXLl3i1VdfNVVcQggjWKhVfPlqNWws1Ry48oDVx29lvzJHT12CgwpOLoUL60wVphBC5Aqjk5vLly8zefJkypUrR8OGDQkJCeGLL77g3r17rF692hwxCiGyoLSHI2PaVADgsy2h3H2cmP3KyjSHxk+HUNg0Ah6GmSBCIYTIHUYnNxUrVuSPP/7g3Xff5datW/qJNJ2cnMwRnxDCCIMaBVKzVDHiktNy9ngKoNlH/0zPsHYQpKWYLlAhhDAjo5IbjUbDvHnz2LZtGyNHjsTb29tccQkhssFCrWLGq9WxtlSz7/J91py8nYPKLOGVRWBbDO6egp0yn5sQomAwKrmxsLBg+PDhxMRkcy4bIYTZlfV0ZHRr3RQln2wKydnjqWJ+0G2ubv3oHAjdZIIIhRDCvIx+LFW1alWuX79ujliEECYy+F+Ppz5cm4PB/QAqdoD6uilP2PCutL8RQuR7Ric3n332GWPGjGHz5s1ERETIjNtC5EOWFmq+fq06tlZqDl59wPKjN3NWYaspULIuJMfAmgGQlmyKMIUQwiyMHsRPrf4nH/r3BJrPZvPWaDSmi84MZBA/UZQsORTG1E0h2FlZ8MeIxgRkd+4pgJjbMK8xJD6EOkOg41emC1QIIV7AmO/vLM0K/m979uzJdmBCiNzVv34A2y/e48j1aMasOcvqt+tjoc7mhJguJaH7j7DiVTi+APzrQ5VXTBuwEEKYQJ7OCp4X5M6NKGpuPUyg/bcHiE9OY3z7irzdtEzOKtw5FQ7OBGtH3fQMHhVME6gQQmTCmO9vo5Ob/fv3Z/p+kyZNjKku10lyI4qi1cfDGfvbeawt1Gx6rxEVvHMwLpUmDX7uBjcOgHsFXYJj42iyWIUQIiNmTW7+3eZGX8m/2t5Imxsh8h9FUXjzpxPs/juKIB9nNrzbABtLi+xXGB8F85tAXARU7g6vLgZVNh93CSFEFph14sxHjx4ZLFFRUWzbto06deqwffv2bActhDAflUrFF69UxdXBmtCIWGZuv5yzCh094bWloLaEi+vgr3kmiVMIIUzB6OTGxcXFYHF3d6d169Z8+eWXfPjhh+aIUQhhAp5OtkzvXhWAHw9c5/C1BzmrsNRL0OZT3fr2jyH8aA4jFEII08jRrOD/5uHhwaVLl0xVnRDCDNpW9qZXHT8UBcb8epaYxNScVVhvqO6xlDZNN/5N3D2TxCmEEDlhdFfwc+fOGbxWFIWIiAi++OILqlevbrLAhBDmMbFTJY5cj+ZmdAKTfr/At71qZr8ylQq6fAdRIXD/b1jTH/ptBEtr0wUshBBGylaDYpVKlW4495deeonFixdTsWJFkwZoatKgWAg4Ff6I1+YdQaNV+LZXDbrWKJGzCh9cgQUtdDOI130LOswwTaBCCPGUWQfxCwsznFdGrVbj4eGBra2tsVUJIfJIrVLFea9FWWbtvMLHGy5Qq1Rx/Fzts1+hezndAH8re8GxH8GnBtTsa7J4hRDCGEa3ufH39zdY/Pz8JLERogAa1rysbnLNpDRGrj5DmkabswortIdm43Xrm0fBnVM5D1IIIbIhy8nNX3/9xR9//GGwbdmyZQQGBuLp6clbb71FcrJMpidEQWFpoWZ2r5o42Vhy8uYjZu+6kvNKm3wI5duDJhlWvwHx93NepxBCGCnLyc2UKVMMGhOfP3+eN998k1atWjFu3Dg2bdrE9OnTzRKkEMI8/Fzt+fxp9/Dv9lzlyLXonFWoVkP3+eBWFmJv6xoYp6WYIFIhhMi6LCc3Z86coWXLlvrXq1atol69eixYsIDRo0cze/Zsfv31V7MEKYQwn87VfelRuySKAqNWn+HRkxwmI7Yu0OsXsHaCm4dg21jTBCqEEFmU5eTm0aNHeHl56V/v27ePdu3a6V/XqVOHW7dumTY6IUSumNKlMqU9HIiMTWLsb+fS9YY0mkcFeGUhoIITi+HYApPEKYQQWZHl5MbLy0vfUyolJYVTp05Rv359/ftxcXFYWVmZPkIhhNnZW1syu1dNrC3UbA+5x89Hb+a80grtoNVk3fofY+H6vpzXKYQQWZDl5KZdu3aMGzeOAwcOMH78eOzt7WncuLH+/XPnzlGmTBmzBCmEML8qJVwY1143TtWnm0O5cCcm55U2HAlVe4Ci0bW/eRj2wl2EECKnspzcfPrpp1hYWNC0aVMWLFjAggULsLb+ZxTSxYsX06ZNG7MEKYTIHQMbBtAqyIsUjZZ3fzlFbFIOp2dQqaDLbPCtBYmPYGVvSIo1TbBCCPEcRo9QHBMTg6OjIxYWFgbbHz58iKOjo0HCkx/JCMVCZO5xQgodZx/kzuNEOlb14fs+NVGpVDmrNPYu/Ngc4iOhXBvotRIsjB5DVAhRhBnz/Z2tWcH/m9gAuLq65vvERgjxYsXsrfm+T02sLFRsOR/BclO0v3H2hd6/gKUtXNkO2yfkvE4hhHgOk80KLoQoPGqWKs649kEAfLI5lPO3TdD+pkQwvDxft/7XPOlBJYQwG0luhBAZGtQwgNaV/ml/E5OYw/Y3AJW7QctJuvU/xsKVnTmvUwgh/kOSGyFEhlQqFV+9Wp2Sxe0If5jA+7+eQavN4fg3AI1GQ/U+T3tQDYB7ITmvUwgh/kWSGyHEc7nYWzGnby2sLdXsDI1i7r5rOa9UpYLO34J/Q0iJg196QlxkzusVQoinJLkRQmSqWsliTOtSGYCvt1/iwBUTTIZpaQ09l4NrGYgJh196QHJ8zusVQggkuRFCZEGvuqXoWdsPrQLDV57mzuPEnFdq7wp914C9G0SchbUDQZOW83qFEEWeJDdCiCyZ2rUyVUo48yghlXeWnyQ5TZPzSt3KQJ9fwdJO10V8y2jI6bxWQogiT5IbIUSW2FpZMLdvMMXsrTh7O4YpGy/mfIJNgJK14dVFgApO/QQHvs55nUKIIk2SGyFElvm52jOrZw1UKlh57BbL/wo3TcUVO0L7L3Xruz+BMytNU68QokiS5EYIYZRmFTz5sK1ugs2pGy/y1/Vo01Rc7y1o8J5u/fd34fJ209QrhChyJLkRQhhtaNPSdK7uS5pW4Z0Vp0zTwBig1TSo1vOfWcRvHTdNvUKIIkWSGyGE0VQqFV++Uo3Kvs5EP0nhrWUnSEwxQQNjtRq6/gBlW0FqAvzyGty/lPN6hRBFiiQ3QohssbO2YP4bwbg6WHPxbixjfztnmgbGFlbQY5luLqrER/Bzd4i5nfN6hRBFhiQ3QohsK1ncnjl9a2GpVrHx7F3m7DXBCMYA1g7QZw24lYPY27oE54mJ2vYIIQo9SW6EEDnyUmk3pjwdwXjGn5fYej7CNBU7uMEb68HJFx5cguXdISnWNHULIQo1SW6EEDn2+kv+DGwYAMDoX89w9tZj01RczA/6bXg6ivEZWNkLUhJMU7cQotCS5EYIYRIfd6xE8woeJKVqGbzsBHdN1YPKowK8vg5snOHmIfj1DUhLMU3dQohCSZIbIYRJWKhVzO5dkwpeTtyPS2bwTyd4kmyiuaJ8a/wzTcPVnbBuCGhN0DtLCFEoSXIjhDAZJ1srFvavjbujNSERsYxYdRqN1kRzRfnXh14rwMIaQjbA78NAqzVN3UKIQkWSGyGESfm52jP/jdpYW6rZGRrF5I0XTNNFHKBsS3h1Mags4OwvsHmEJDhCiHQkuRFCmFywf3G+fToH1fKj4czbd910lQd1hlcWgEoNp5bB1vdlJnEhhAFJboQQZtG+qg8TO1YC4H/b/ub3M3dMV3mVV6DbPEAFJxbDHx9KgiOE0JPkRghhNoMaBfJmo0AAxqw5y+FrD0xXefWeuqkaUMGxH+HPCZLgCCEASW6EEGY2oUMQHav6kKpRePvnk/wdacKB+Gr2hc7f6taP/gB/fiQJjhBCkhshhHmp1Sq+7lGdugGuxCWl0W/RMcKjTTgQX3B/6PSNbv3oHHlEJYSQ5EYIYX62VhYs6Febit5ORMUl88biv4iKSzLdAWoPgi7foX9EtXmU9KISogiT5EYIkStc7K1YNqgufq523IxOoP/i48QkppruALX6Qbc5gApOLoFNwyXBEaKIkuRGCJFrPJ1t+XlQPdwdbQiNiGXITydISjXhSMM1+kD3H3XdxE//DBuGgsZEoyQLIQqMPE9u5syZQ2BgILa2tgQHB3PgwIHnll23bh2tW7fGw8MDZ2dn6tevz59//pmL0QohcirA3YFlg+riZGvJsRsPeWfFKVLSTHiHpVoPeGWhbqC/c6thTX9ISzZd/UKIfC9Pk5vVq1czcuRIJkyYwOnTp2ncuDHt27cnPDw8w/L79++ndevWbN26lZMnT9K8eXM6d+7M6dOnczlyIUROVPJ1ZlH/OthYqtn9dxQjVp0mTWPCBKfKK9BzOVjYwN+b4ZcekBxvuvqFEPmaSjHZuOjGq1evHrVq1WLu3Ln6bUFBQXTr1o3p06dnqY7KlSvTs2dPJk2alKXysbGxuLi4EBMTg7Ozc7biFkKYxt5LUby17CQpGi1da/gys0cNLNQq0x3g+l5Y2QdSn0DJutD3V7Arbrr6hRC5xpjv7zy7c5OSksLJkydp06aNwfY2bdpw+PDhLNWh1WqJi4vD1dX1uWWSk5OJjY01WIQQ+UOzCp780LcWlmoVv5+5y/h159CaaqJNgNLNoP9GsC0Gt4/B0s4Qd8909Qsh8qU8S24ePHiARqPBy8vLYLuXlxeRkZFZquPrr7/myZMn9OjR47llpk+fjouLi37x8/PLUdxCCNNqXcmLb3vVRK2CX0/cZvLGi6abaBOgZG0YsAUcPOHeeVjUGqKvma5+IUS+k+cNilUqw1vQiqKk25aRlStXMmXKFFavXo2np+dzy40fP56YmBj9cuvWrRzHLIQwrY7VfPi6R3VUKvj56E2mbgoxbYLjXQXe/BOKB8Ljm7oE585J09UvhMhX8iy5cXd3x8LCIt1dmqioqHR3c/5r9erVvPnmm/z666+0atUq07I2NjY4OzsbLEKI/OflmiX5ontVAJYevsGk3y+a9hGVa2l4czv4VIeEaN0jqqs7TVe/ECLfyLPkxtramuDgYHbs2GGwfceOHTRo0OC5+61cuZIBAwbwyy+/0LFjR3OHKYTIRT3rlOLLV6rp7+B8/PsF0yY4jp66R1Slm+saGf/SE86uMl39Qoh8IU8fS40ePZqFCxeyePFiQkNDGTVqFOHh4QwdOhTQPVLq16+fvvzKlSvp168fX3/9NS+99BKRkZFERkYSExOTV6cghDCxHnX8mPGq7hHVL3+F89H686ZNcGycoM+vUPU10KbB+rdh7xcyH5UQhUieJjc9e/Zk1qxZTJs2jRo1arB//362bt2Kv78/ABEREQZj3syfP5+0tDTeffddfHx89MuIESPy6hSEEGbwanBJZvaojloFq47f4sPfzqExZYJjaQ0v/wgNR+pe750O64fKYH9CFBJ5Os5NXpBxboQoOH4/c4dRq8+gVaBjVR++6VkDa0sT/012cilsHg2KBvwbQc+fwf75w0sIIfJGgRjnRgghXqRrjRLM6VsLKwsVW85HMHjZCRJSTDxXVPAA6LsGrJ3g5kHpKi5EISDJjRAiX2tXxYdF/etgZ2XB/sv36bfomGlnEwco21LXVdy5JERfhQUt4Noe0x5DCJFrJLkRQuR7Tcp7sHywbrLNEzcf0fvHozyIN3H7GK/KMGQXlKgNSY9h+StwdJ40NBaiAJLkRghRIAT7u7LqrZdwd7QmJCKWV+Ye5saDJ6Y9iJO3rqt49d66NjjbxsLG9yAtxbTHEUKYlSQ3QogCo7KvC7++XZ+Sxe24GZ3AK3MPc+bWY9MexMoWus2FNp+CSg2nf4afOkNc1qaFEULkPUluhBAFSmkPR9a904AqJZyJfpJC7x+PsvtvE0+GqVJBg/d04+HYuMCtozC/CdzM2qS+Qoi8JcmNEKLA8XSyZdVb9WlS3oPEVA1Dlp1k1bHwF+9orHKt4a094FkJ4u/p7uAcnSvtcITI5yS5EUIUSI42lizqX5tXg0ui0SqMW3ee6X+EmnY0YwC3MjB4J1R5RTei8bZx8NtgSDFxex8hhMlIciOEKLCsLNTMeLUaw1uWA2D+vuu8vfwkT5JNPBaOtQO8sgjafQEqC7iwFn5sBvcumvY4QgiTkORGCFGgqVQqRrcuz6ynoxfvCLnHa/OOcPdxoqkPBC/9HwzYDE4+8OCybjyckz/JYyoh8hlJboQQhUK3miVYOeSfruJdfzjE6fBHpj+QfwMYehDKtIS0JNg0HNYNgeQ40x9LCJEtktwIIQqNYP/ibHi3IRW8nLgfl0zP+UdZfdwMDY0d3KHvWmg5WfeY6vwaXW+qOydNfywhhNEkuRFCFColi9vz2zsNaF3JixSNlrG/nWfC+vOkpGlNeyC1GhqPhoFbwbkEPLwOi9rA/hmg1Zj2WEIIo0hyI4QodBxtLJn/ejDvty6PSgUr/gqn94Kj3ItNMv3BSr2ke0xVqZuuN9XuT2FpJ3hshjtGQogskeRGCFEoqdUq3mtZjsX96+Bka8nJm4/o9N1BjlyLNv3B7F3htaW6kY2tHSH8MMxtCGdWSmNjIfKAJDdCiEKteUVPNg1rpG+H03fhUb7bdQWNqcfDUamgRh/dXZySdSE5FjYMhVV9IM7EIygLITIlyY0QotALcHdg/bsNeDW4JFoFvt5xmQFLjpl+ZnEA10AY+Ae0nARqK7i0FebUg/Nr5S6OELlEpShF67ctNjYWFxcXYmJicHZ2zutwhBC5bO3J20zccIHEVA2eTjbM6lWDBmXczXOwexdhw/9BxFnd64qdoOPXutnHhRBGMeb7W+7cCCGKlFeDS7JxWEPKeToSFZdM34V/MX1rKMlpZujh5FUZBu+CZh+B2hL+3gzf14UTS0Br4t5bQgg9uXMjhCiSElLSmLYphFXHbwFQyceZb3vVoJyXk3kOGHkeNg6Hu6d0r0s1gM7fgkd58xxPiELGmO9vSW6EEEXanxcjGffbOR4lpGJjqWZ8+4r0qx+AWq0y/cG0Gvhrvq67eOoTsLCGRqOh0UiwsjP98YQoRCS5yYQkN0KI/4qKTeKDtefYd/k+AA3KuPG/V6rh52pvngM+DofNo+HqDt3r4gHQ/kso39Y8xxOiEJDkJhOS3AghMqIoCj8dvsEX2/4mKVWLvbUF49pX5PV6/ua5i6MoELIBtn0EcXd12yp0gHbTdcmOEMKAJDeZkORGCJGZGw+e8OFv5zgW9hCAeoGufPlqNfzdHMxzwOR42Pc/ODpHN8KxpS3UHwaNRoGNo3mOKUQBJMlNJiS5EUK8iFar8PPRm3zxx98kpmqwsVQzvGU5hjQujbWlmTqZRv0NW8fAjQO6147eurFyqvfWzWMlRBEnyU0mJLkRQmRVeHQC49ad4/DTKRvKezny2ctVqRPgap4DKgqEboIdE+HRDd02n+rQ9nMIaGSeYwpRQEhykwlJboQQxlAUhfWn7/DZllCin6QA0LO2H2PbV8TVwdo8B01L1vWq2j9DN40DQNnW0GoyeFc1zzGFyOckucmEJDdCiOx4nJDCF3/8rR8Xx8XOilGtyvH6S/5YWpjpsVH8fdg7HU79pGuPgwqq9YDmH0mjY1HkSHKTCUluhBA5cfzGQyb9fpHQCN0dlfJejkzuXJmGZc00hQNA9DXd2DgX1+leq62gVj9o/D64lDDfcYXIRyS5yYQkN0KInNJoFVYeC+fr7Zd4lJAKQOtKXoxtV5Gynmbs4XT3NOycCtf36F5bWEPwAN1AgM4+5juuEPmAJDeZkORGCGEqjxNSmLXzCj8fvYlGq2ChVtGjth+jWpXD09nWfAe+cRD2TIebB3WvLWwguD80GA7F/Mx3XCHykCQ3mZDkRghhalfuxfG/bX+zMzQKADsrCwY3DuStJqVxsrUy34HD9uuSnPDDutdqS6jWExqOlDmrRKEjyU0mJLkRQpjLsbCHfL41lDO3HgO6RsdvNSlN/wYBONpYmuegigJh++DATN2/AKggqBM0GAF+dcxzXCFymSQ3mZDkRghhToqisO1CJDO2X+L6/ScAFLe34q0mZehX3x8HcyU5ALdP6JKcS1v+2eZXD+q/CxU7gdrCfMcWwswkucmEJDdCiNyg0SpsOnuXb3ddIeyBLslxdbBmYIMA+tUPwMXejI+rokLh0Gw4vwa0ugbPFCsF9YZCjT5gV9x8xxbCTCS5yYQkN0KI3JSm0fL7mbvM3n2Fm9EJADjaWNK3XinebBRo3obHcZFwfCEcXwSJurmysLTTjZVTd4gMCCgKFEluMiHJjRAiL6RptGw5H8Hcvdf4OzIOAGsLNS/XLMHARgFU9Dbj/0cpCXBuFRxbAFEh/2z3ewlqD4SgLmBtb77jC2ECktxkQpIbIUReUhSF3X9HMWfvNU7efKTf3qisO4MaBdCsvCdqtcpcB4ebh+H4At0cVto03XYbF6j2GtTqDz7VzHNsIXJIkptMSHIjhMgvTtx4yOJDYWy7EIn26f/Ege4O9K1XildqlaS4ueauAoiNgNM/65bH4f9s96mum4m8yqvg6GG+4wthJEluMiHJjRAiv7n9KIFlR26y8lg4cUm6uynWlmo6VfWh70ulqFWqOCqVme7maLW6LuSnlsHfm0GjmxwUlQWUbQXVe0KFDmBlZ57jC5FFktxkQpIbIUR+9SQ5jY1n77L86E0u3o3Vby/n6chrtUvSrWYJPJ3M2AD5STRc+E3XPufOyX+2WztChfZQ5RUo0wIsbcwXgxDPIclNJiS5EULkd4qicPZ2DCuO3mTTubskpWoBsFCraFreg9eCS9K8oie2VmYct+b+ZTi3Gs79CjH/emxl46IbIDCoM5RuDlZmTLaE+BdJbjIhyY0QoiCJTUpl89kI1p68xanwx/rtTjaWtKnsTZcavjQs44alhdo8ASiKbnDAC7/BxfUQH/nPe9aOUK61boDAcm3AVv5PFeYjyU0mJLkRQhRU1+7Hs/bkbX4/fYe7MUn67W4O1rSr4k37Kj7UK+2KlbkSHa0Gwo9AyEZd+5zYO/+8p7YE/wZQvh2UawvuZc0TgyiyJLnJhCQ3QoiCTqtVOHHzERvP3mHr+UgePknRv+diZ0WrIC/aVfGmUVl37KzN9OhKUeDuKQjdrEt0Hlw2fN+1tK5BcpkWENAYbBzNE4coMiS5yYQkN0KIwiRVo+XwtWi2XYhg+8V7RP8r0bGxVNOgjBstKnrSvKInJYubcaC+6GtwZTtc3gY3Dv0z7QOA2gpKvQSlm0JAEyhRCyzMOP2EKJQkucmEJDdCiMJKo1U4ceMh2y5Gsv3iPe48TjR4v7yXI43LedConDv1Al2xtzbTJJ7JcXB9H1zbBVd3weObhu9bOYB/fd0dHf8G4FMDLM04po8oFCS5yYQkN0KIokBRFC7fi2fX3/fY83cUJ28+0g8UCLqpH4L9i9OgjBsvlXGjWkkXbCzN9Ajr4XVdkhO2H24c/Geeq2cs7aBkbd3dnVIvQYlgmdxTpCPJTSYkuRFCFEWPE1I4dDWaA1fuc+DKg3R3dWyt1NQqVZx6gW7UDihOdb9iONqY4c6OVgtRFyHsgC7RCT+SPtkBcK8AJevokp4StcAjSO7uFHGS3GRCkhshRFGnKAphD55w6OoDjl5/yNHr0QZtdQDUKgjycaa2f3FqlCpGtZLFCHRzMP28V4qia4wcfgRuHoHbx3R3ev7Lwhq8qoBvTfCtoZvR3CNIxtkpQiS5yYQkN0IIYUhRFK7dj+fItWiO33jEyZuP0t3ZAXCytaRaSReqlSxGZV9nKvu64O9qb/qE58kD3dg6t4/rlogzkBSTvpzKAjwq6BIdr8rgWQk8g8C5BJhrugqRZyS5yYQkN0II8WIRMYmcvKlLdM7djuHCnRiS07TpyjnaWBLk40RFb2cqeDtRwduJ8l5OuNiZsDeUosCjMLh7WrdEnIPI8xk/zgKwcQaPiuBRHtzL6x5xuZeDYv5gYaZG1MLsJLnJhCQ3QghhvFSNlsv34jh3O4Zzt2MIiYjl74jYDBMeAB8XW8p6OlLGw5Eyno6U8XCgjIcjnk42ppkEVFEg9i7cuwCR5yAqVLc8uAzatIz3UVtB8QBwKwtuZXRj8bgG6ra5+En39HxOkptMSHIjhBCmkabRcv3BEy7ejeHvyDguR8ZxKTLOYPTk/7K3tsDfzYEAN3sC3B0o5WqPX3F7/Fzt8C1ml/PRldNSIPoq3A+FB1d0yc79yxB9BdKeHxcqNbiU1N3dKVZKt7j4QTE/3WMu5xLSviePSXKTCUluhBDCvGISU7kaFce1qCdcvR/Ptah4rt6P59bDBIPu6P+lVoGPix2+xWzxLWanX3ycbfF20S2u9tbZa+Oj1eqmi4i+Cg+vQfR13b+PbuiWzBKfZ+zdwdlXl+g4eevWnbzByQccvXSLgzuozTihaREmyU0mJLkRQoi8kZKm5fajBG5EPyHsQQI3o59w62ECtx4lcuthwnMfcf2blYUKTydbPJxs8HSyefqvLe5O1rg52ODuaI27ow2ujtY42Vhm7RGYokD8PXgYBo/DdbOgPw6Hx7d0/8behbT0DawzpFLrkiBHT12i4+ChW+zddK/t3f5Z7Fx14/lIO6AsKVDJzZw5c5gxYwYRERFUrlyZWbNm0bhx4+eW37dvH6NHj+bixYv4+vry4YcfMnTo0CwfT5IbIYTIfxRF4X5cMrcfJ3L36XLnUSJ3HidxLzaJiJgkop8kY8w3lqVaRXEHa1ztrSnuYEUxO2uK2VvhYv/PurOtFc52lk//tcLJ1hInW0vDAQ0VBRIf6e78xNyBuLsQF6lLeuIiIS4C4qPgyX0gG1+pNs5gVwxsi/3zr63L03UXsHHRzbhu4/z0XyfdjOw2zro5uyxti0TvMGO+v/M0XVy9ejUjR45kzpw5NGzYkPnz59O+fXtCQkIoVapUuvJhYWF06NCBIUOGsHz5cg4dOsQ777yDh4cHr7zySh6cgRBCCFNQqVR4Otvi6WxLrVIZj06cqtFyLzaJqLhkomKTuR+fzP2nr6OfpPAgPpno+BSi45N5kqIhTatLmO7HJRsdj7WFGkdbSxxtLHGwscTRxgIHG0scrD1xsPHB3toSOxsLHJwssLO2xN7aAnsLBSclFmfNQxxTHmKX9gi7lEfYJEdjnRyNZdIjLJIeokp8iCohGpIe6w6WHKtbCM/eD09tCdYOuoTH2vHpugNY2YO1vW66C2t73Wsre13bISt7XVJkZff0X1vdSNFWtrrXlrZgaQMWNrrBEy1tdWMNFZAkKk/v3NSrV49atWoxd+5c/bagoCC6devG9OnT05UfO3YsGzduJDQ0VL9t6NChnD17liNHjmTpmHLnRgghCr+kVA2PElJ4+CSFR09SeZiQQkxCCo8TUnmcmMrjhFRiElOITUwjNimV2MRUYpPSiE9+Tk8rE7O2UGNnCW4WibhZPMHNIgE3VTwu6gRcVE9wJgFnnuBEPA5KAvZKAnbaJ9hpn2CrTcBW+wQbbRYflZmQVm2NYqFbsLBCsbBBUVuBhTWKhRWorVHUligO7tj1XWHSYxeIOzcpKSmcPHmScePGGWxv06YNhw8fznCfI0eO0KZNG4Ntbdu2ZdGiRaSmpmJllb4bX3JyMsnJ/2TtsbGxJoheCCFEfmZrZYGPix0+LnZG7afVKsSnpBH/NNGJS0olPllDQrLu9ZPkNJ6kaEhISSMhRUNiioYnT/9NStWQmPrPelKqhqQ0rX79342pUzRaUjQQgw3XsQFcjT5HFVocSMKRROxVydiThAPJOKgSsScZe1WS7l+SsVMlY0cytqRgp0rBlmTsSMGWFGxVKdiSig0p2Kie/ksaNqRio0o1OKZamwLaFEh9TlBPReGKcT9508qz5ObBgwdoNBq8vLwMtnt5eREZGZnhPpGRkRmWT0tL48GDB/j4+KTbZ/r06UydOtV0gQshhCi01GqVrv2NrWnHvFEUhTStQnKaluRUje7fNC0pzxaNRv86VaOQqnm2XUva09e6RSFNoyVVq/s3Tat7T6PV1Z/2tLxGUUjQKsQ93a55umgVRf++Vmv4b5pGQVFAq+heKwpoNRosScVCm4qlkoJam4oVqVgpqVgoqVjqlzQs0OjX7e1s+dykP0Hj5HkT7f+2ZFcUJdPW7RmVz2j7M+PHj2f06NH617Gxsfj5+WU3XCGEEMJoKpUKKwsVVhZq80xIKgzk2U/Y3d0dCwuLdHdpoqKi0t2decbb2zvD8paWlri5uWW4j42NDTY2NqYJWgghhBD5Xg6Hgsw+a2trgoOD2bFjh8H2HTt20KBBgwz3qV+/frry27dvp3bt2hm2txFCCCFE0ZNnyQ3A6NGjWbhwIYsXLyY0NJRRo0YRHh6uH7dm/Pjx9OvXT19+6NCh3Lx5k9GjRxMaGsrixYtZtGgRY8aMyatTEEIIIUQ+k6cP/nr27El0dDTTpk0jIiKCKlWqsHXrVvz9/QGIiIggPPyffv+BgYFs3bqVUaNG8cMPP+Dr68vs2bNljBshhBBC6OX5CMW5Tca5EUIIIQoeY76/8/SxlBBCCCGEqUlyI4QQQohCRZIbIYQQQhQqktwIIYQQolCR5EYIIYQQhYokN0IIIYQoVCS5EUIIIUShIsmNEEIIIQoVSW6EEEIIUagUuXnXnw3IHBsbm8eRCCGEECKrnn1vZ2VihSKX3MTFxQHg5+eXx5EIIYQQwlhxcXG4uLhkWqbIzS2l1Wq5e/cuTk5OqFQqk9QZGxuLn58ft27dKpTzVRX28wM5x8KgsJ8fyDkWBoX9/MB856goCnFxcfj6+qJWZ96qpsjduVGr1ZQsWdIsdTs7OxfaDysU/vMDOcfCoLCfH8g5FgaF/fzAPOf4ojs2z0iDYiGEEEIUKpLcCCGEEKJQkeTGBGxsbJg8eTI2NjZ5HYpZFPbzAznHwqCwnx/IORYGhf38IH+cY5FrUCyEEEKIwk3u3AghhBCiUJHkRgghhBCFiiQ3QgghhChUJLkRQgghRKEiyU0OzZkzh8DAQGxtbQkODubAgQN5HVK2TJ8+nTp16uDk5ISnpyfdunXj0qVLBmUGDBiASqUyWF566aU8ith4U6ZMSRe/t7e3/n1FUZgyZQq+vr7Y2dnRrFkzLl68mIcRGy8gICDdOapUKt59912gYF7D/fv307lzZ3x9fVGpVGzYsMHg/axct+TkZN577z3c3d1xcHCgS5cu3L59OxfP4vkyO7/U1FTGjh1L1apVcXBwwNfXl379+nH37l2DOpo1a5buuvbq1SuXz+T5XnQNs/K5LKjXEMjwd1KlUjFjxgx9mfx+DbPyHZGffhclucmB1atXM3LkSCZMmMDp06dp3Lgx7du3Jzw8PK9DM9q+fft49913OXr0KDt27CAtLY02bdrw5MkTg3Lt2rUjIiJCv2zdujWPIs6eypUrG8R//vx5/XtffvklM2fO5Pvvv+f48eN4e3vTunVr/XxkBcHx48cNzm/Hjh0AvPbaa/oyBe0aPnnyhOrVq/P9999n+H5WrtvIkSNZv349q1at4uDBg8THx9OpUyc0Gk1uncZzZXZ+CQkJnDp1iokTJ3Lq1CnWrVvH5cuX6dKlS7qyQ4YMMbiu8+fPz43ws+RF1xBe/LksqNcQMDiviIgIFi9ejEql4pVXXjEol5+vYVa+I/LV76Iisq1u3brK0KFDDbZVrFhRGTduXB5FZDpRUVEKoOzbt0+/rX///krXrl3zLqgcmjx5slK9evUM39NqtYq3t7fyxRdf6LclJSUpLi4uyrx583IpQtMbMWKEUqZMGUWr1SqKUvCvIaCsX79e/zor1+3x48eKlZWVsmrVKn2ZO3fuKGq1Wtm2bVuuxZ4V/z2/jBw7dkwBlJs3b+q3NW3aVBkxYoR5gzORjM7xRZ/LwnYNu3btqrRo0cJgW0G6hoqS/jsiv/0uyp2bbEpJSeHkyZO0adPGYHubNm04fPhwHkVlOjExMQC4uroabN+7dy+enp6UL1+eIUOGEBUVlRfhZduVK1fw9fUlMDCQXr16cf36dQDCwsKIjIw0uJ42NjY0bdq0wF7PlJQUli9fzqBBgwwmiS3o1/DfsnLdTp48SWpqqkEZX19fqlSpUiCvbUxMDCqVimLFihlsX7FiBe7u7lSuXJkxY8YUqDuOkPnnsjBdw3v37rFlyxbefPPNdO8VpGv43++I/Pa7WOQmzjSVBw8eoNFo8PLyMtju5eVFZGRkHkVlGoqiMHr0aBo1akSVKlX029u3b89rr72Gv78/YWFhTJw4kRYtWnDy5MkCMdpmvXr1WLZsGeXLl+fevXt8+umnNGjQgIsXL+qvWUbX8+bNm3kRbo5t2LCBx48fM2DAAP22gn4N/ysr1y0yMhJra2uKFy+erkxB+11NSkpi3Lhx9OnTx2BCwr59+xIYGIi3tzcXLlxg/PjxnD17Vv9YMr970eeyMF3Dn376CScnJ7p3726wvSBdw4y+I/Lb76IkNzn077+IQXfR/7utoBk2bBjnzp3j4MGDBtt79uypX69SpQq1a9fG39+fLVu2pPtFzY/at2+vX69atSr169enTJky/PTTT/rGi4Xpei5atIj27dvj6+ur31bQr+HzZOe6FbRrm5qaSq9evdBqtcyZM8fgvSFDhujXq1SpQrly5ahduzanTp2iVq1auR2q0bL7uSxo1xBg8eLF9O3bF1tbW4PtBekaPu87AvLP76I8lsomd3d3LCws0mWbUVFR6TLXguS9995j48aN7Nmzh5IlS2Za1sfHB39/f65cuZJL0ZmWg4MDVatW5cqVK/peU4Xlet68eZOdO3cyePDgTMsV9GuYlevm7e1NSkoKjx49em6Z/C41NZUePXoQFhbGjh07DO7aZKRWrVpYWVkV2Ov6389lYbiGAAcOHODSpUsv/L2E/HsNn/cdkd9+FyW5ySZra2uCg4PT3TLcsWMHDRo0yKOosk9RFIYNG8a6devYvXs3gYGBL9wnOjqaW7du4ePjkwsRml5ycjKhoaH4+Pjobwf/+3qmpKSwb9++Ank9lyxZgqenJx07dsy0XEG/hlm5bsHBwVhZWRmUiYiI4MKFCwXi2j5LbK5cucLOnTtxc3N74T4XL14kNTW1wF7X/34uC/o1fGbRokUEBwdTvXr1F5bNb9fwRd8R+e530aTNk4uYVatWKVZWVsqiRYuUkJAQZeTIkYqDg4Ny48aNvA7NaP/3f/+nuLi4KHv37lUiIiL0S0JCgqIoihIXF6e8//77yuHDh5WwsDBlz549Sv369ZUSJUoosbGxeRx91rz//vvK3r17levXrytHjx5VOnXqpDg5Oemv1xdffKG4uLgo69atU86fP6/07t1b8fHxKTDn94xGo1FKlSqljB071mB7Qb2GcXFxyunTp5XTp08rgDJz5kzl9OnT+t5CWbluQ4cOVUqWLKns3LlTOXXqlNKiRQulevXqSlpaWl6dll5m55eamqp06dJFKVmypHLmzBmD383k5GRFURTl6tWrytSpU5Xjx48rYWFhypYtW5SKFSsqNWvWzBfnpyiZn2NWP5cF9Ro+ExMTo9jb2ytz585Nt39BuIYv+o5QlPz1uyjJTQ798MMPir+/v2Jtba3UqlXLoOt0QQJkuCxZskRRFEVJSEhQ2rRpo3h4eChWVlZKqVKllP79+yvh4eF5G7gRevbsqfj4+ChWVlaKr6+v0r17d+XixYv697VarTJ58mTF29tbsbGxUZo0aaKcP38+DyPOnj///FMBlEuXLhlsL6jXcM+ePRl+Nvv3768oStauW2JiojJs2DDF1dVVsbOzUzp16pRvzjuz8wsLC3vu7+aePXsURVGU8PBwpUmTJoqrq6tibW2tlClTRhk+fLgSHR2dtyf2L5mdY1Y/lwX1Gj4zf/58xc7OTnn8+HG6/QvCNXzRd4Si5K/fRdXToIUQQgghCgVpcyOEEEKIQkWSGyGEEEIUKpLcCCGEEKJQkeRGCCGEEIWKJDdCCCGEKFQkuRFCCCFEoSLJjRBCCCEKFUluhBBCCFGoSHIjhMixKVOmUKNGjTw7/sSJE3nrrbf0r5s1a8bIkSPzLJ7M1KlTh3Xr1uV1GEIUapLcCCEypVKpMl0GDBjAmDFj2LVrV57Ed+/ePb799ls++uijPDm+sSZOnMi4cePQarV5HYoQhZYkN0KITEVEROiXWbNm4ezsbLDt22+/xdHRMUuzVZvDokWLqF+/PgEBAXly/H9LSUl5YZmOHTsSExPDn3/+mQsRCVE0SXIjhMiUt7e3fnFxcUGlUqXb9t/HUgMGDKBbt258/vnneHl5UaxYMaZOnUpaWhoffPABrq6ulCxZksWLFxsc686dO/Ts2ZPixYvj5uZG165duXHjRqbxrVq1ii5duqTbrtVq+fDDD3F1dcXb25spU6YYvB8eHk7Xrl1xdHTE2dmZHj16cO/evXTn8G8jR46kWbNm+tfNmjVj2LBhjB49Gnd3d1q3bg3oHtOVKlUKGxsbfH19GT58uH4fCwsLOnTowMqVKzM9LyFE9klyI4Qwi927d3P37l3279/PzJkzmTJlCp06daJ48eL89ddfDB06lKFDh3Lr1i0AEhISaN68OY6Ojuzfv5+DBw/i6OhIu3btnntH5NGjR1y4cIHatWune++nn37CwcGBv/76iy+//JJp06axY8cOABRFoVu3bjx8+JB9+/axY8cOrl27Rs+ePY0+z59++glLS0sOHTrE/PnzWbt2Ld988w3z58/nypUrbNiwgapVqxrsU7duXQ4cOGD0sYQQWWOZ1wEIIQonV1dXZs+ejVqtpkKFCnz55ZckJCTo28aMHz+eL774gkOHDtGrVy9WrVqFWq1m4cKFqFQqAJYsWUKxYsXYu3cvbdq0SXeMmzdvoigKvr6+6d6rVq0akydPBqBcuXJ8//337Nq1i9atW7Nz507OnTtHWFgYfn5+APz8889UrlyZ48ePU6dOnSyfZ9myZfnyyy/1r7du3Yq3tzetWrXCysqKUqVKUbduXYN9SpQoQXh4OFqtFrVa/sYUwtTkt0oIYRaVK1c2+OL28vIyuINhYWGBm5sbUVFRAJw8eZKrV6/i5OSEo6Mjjo6OuLq6kpSUxLVr1zI8RmJiIgC2trbp3qtWrZrBax8fH/2xQkND8fPz0yc2AJUqVaJYsWKEhoYadZ7/vWv02muvkZiYSOnSpRkyZAjr168nLS3NoIydnR1arZbk5GSjjiWEyBq5cyOEMAsrKyuD1yqVKsNtz3oNabVagoODWbFiRbq6PDw8MjyGu7s7oHs89d8ymR1LURT93aF/+/d2tVqNoigG76empqbbx8HBweC1n58fly5dYseOHezcuZN33nmHGTNmsG/fPn1MDx8+xN7eHjs7uwzPSwiRM3LnRgiRL9SqVYsrV67g6elJ2bJlDRYXF5cM9ylTpgzOzs6EhIQYdaxKlSoRHh6ub+8DEBISQkxMDEFBQYAuoYqIiDDY78yZM1mq387Oji5dujB79mz27t3LkSNHOH/+vP79CxcuUKtWLaNiFkJknSQ3Qoh8oW/fvri7u9O1a1cOHDhAWFgY+/btY8SIEdy+fTvDfdRqNa1ateLgwYNGHatVq1ZUq1aNvn37curUKY4dO0a/fv1o2rSp/jFTixYtOHHiBMuWLePKlStMnjyZCxcuvLDupUuXsmjRIi5cuMD169f5+eefsbOzw9///9u5exWFgSgMw9+0FkEbyyykS2MnBASrQPp0aQKxSG2jnZV3EW8gdoK1jYHcSS4gpEiZ7Sz2B1dYUIf3qWeGU36cmTkftzVVVf34hgjA/yDcAHgJo9FI1+tVrusqjmP5vq/VaqW+7+U4zq/78jxXWZYPDcUzxuh0OmkymWi5XCoMQ3mep+PxeFsTRZF2u522263m87m6rlOapnfPHo/HOhwOWiwWms1mulwuOp/PtzlATdOormtlWfbnegE8xgxfL5UB4I0Mw6AgCLRer5UkybPLuWuz2ahtWxVF8exSAGvRuQHw1owxKori24+kVzWdTrXf759dBmA1OjcAAMAqdG4AAIBVCDcAAMAqhBsAAGAVwg0AALAK4QYAAFiFcAMAAKxCuAEAAFYh3AAAAKsQbgAAgFU+AR148SJa9yhbAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "#visualize survival function with mean lambda\n", + "\n", + "import matplotlib.pyplot as plt\n", + "\n", + "# MLE\n", + "λ_mle = np.sum(t**2) / (2*sum(c))\n", + "print(str(\"MLE version of λ:\") + \" \" + str(λ_mle))\n", + "\n", + "\n", + "x = np.linspace(1,200,100)\n", + "y1 = np.exp(-(x**2)/(2*az.summary(trace)['mean']['λ']))\n", + "y2 = np.exp(-(x**2)/(2*λ_mle))\n", + "\n", + "plt.plot(x,y1,label=\"Bayes Estimator\")\n", + "plt.plot(x,y2,label=\"MLE\")\n", + "plt.title(\"Rayleigh Survival Functions\")\n", + "plt.xlabel(\"Time (hours)\")\n", + "plt.ylabel(\"Survival Probability\")\n", + "plt.legend()\n", + "plt.show()" + ] }, { "cell_type": "markdown", @@ -146,22 +821,305 @@ "```{admonition} Solution\n", ":class: tip, dropdown\n", "\n", - "Solution to be added. Feel free to share yours on Ed Discussion!\n", + "Solution in hidden cell below!\n", "\n", "```" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 6, "id": "95c9385d", "metadata": { "tags": [ "hide-cell" ] }, - "outputs": [], - "source": [] + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/opt/anaconda3/envs/pymc_env/lib/python3.13/site-packages/pymc/model/core.py:1311: ImputationWarning: Data in x_impute contains missing values and will be automatically imputed from the sampling distribution.\n", + " warnings.warn(impute_message, ImputationWarning)\n", + "/opt/anaconda3/envs/pymc_env/lib/python3.13/site-packages/pymc/model/core.py:1311: ImputationWarning: Data in likelihood contains missing values and will be automatically imputed from the sampling distribution.\n", + " warnings.warn(impute_message, ImputationWarning)\n" + ] + }, + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "c2567e1bc3974f07a3d0267bee913a4e", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + " 0%| | 0/6000 [00:00\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
meansdhdi_2.5%hdi_97.5%mcse_meanmcse_sdess_bulkess_tailr_hat
alpha0.5500.0150.5210.5780.00.05980.08982.01.0
beta1-0.4160.016-0.447-0.3840.00.06393.09438.01.0
beta2-0.5890.022-0.631-0.5450.00.011754.012590.01.0
\n", + "" + ], + "text/plain": [ + " mean sd hdi_2.5% hdi_97.5% mcse_mean mcse_sd ess_bulk \\\n", + "alpha 0.550 0.015 0.521 0.578 0.0 0.0 5980.0 \n", + "beta1 -0.416 0.016 -0.447 -0.384 0.0 0.0 6393.0 \n", + "beta2 -0.589 0.022 -0.631 -0.545 0.0 0.0 11754.0 \n", + "\n", + " ess_tail r_hat \n", + "alpha 8982.0 1.0 \n", + "beta1 9438.0 1.0 \n", + "beta2 12590.0 1.0 " + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "import pymc as pm\n", + "import arviz as az\n", + "import pandas as pd\n", + "import numpy as np\n", + "from pymc.math import switch, ge\n", + "\n", + "data = pd.read_csv(\"../data/stagnant_water.csv\")\n", + "\n", + "x = data[\"x\"].to_numpy()\n", + "y = data[\"y\"].to_numpy()\n", + "\n", + "with pm.Model() as m:\n", + "\n", + " tau = pm.Gamma(\"tau\", .001, .001)\n", + " alpha = pm.Normal(\"alpha\", 0, sigma = 1000)\n", + " beta1 = pm.Normal(\"beta1\", 0, sigma = 1000)\n", + " beta2 = pm.Normal(\"beta2\", 0, sigma = 1000)\n", + " theta = pm.Uniform(\"theta\",-1.3,1.1)\n", + "\n", + " x_impute = pm.Uniform(\"x_impute\", -5,5,observed=x)\n", + "\n", + " mu = alpha + beta1 * x_impute + beta2 * (x_impute-theta) * switch(ge(x_impute - theta, 0), 1, 0)\n", + "\n", + " likelihood = pm.Normal(\"likelihood\", mu=mu, tau=tau, observed=y)\n", + "\n", + " trace = pm.sample(5000, nuts_sampler = 'numpyro', target_accept=.8)\n", + " # trace = pm.sample(5000)\n", + "\n", + "az.summary(trace, hdi_prob=.95, var_names=['alpha','beta1','beta2'])" + ] } ], "metadata": { @@ -170,7 +1128,7 @@ "formats": "ipynb,py,md" }, "kernelspec": { - "display_name": "Python 3 (ipykernel)", + "display_name": "pymc_env", "language": "python", "name": "python3" }, @@ -184,7 +1142,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.1.-1" + "version": "3.13.3" } }, "nbformat": 4,