added example MC_3_1

3f7e71b8 · Mirko Birbaumer · b66aa7a5 · 3f7e71b8 · 3f7e71b8
Commit 3f7e71b8 authored 2 weeks ago by Mirko Birbaumer
--- a/notebooks/Model Comparison/Daten/dummy.csv
+++ b/notebooks/Model Comparison/Daten/dummy.csv
+-1.081 9.357
+-0.888 8.322
+-0.623 8.745
+-0.480 7.346
+-0.505 7.028
+-0.345 4.032
+0.052 1.490
+0.061 3.265
+0.215 4.525
+0.395 3.339
+0.643 2.045
+0.871 -1.336
+0.776 1.264
+0.833 -0.813
+1.224 -0.335
+1.300 0.152
+1.563 -1.797
+1.669 -4.772
+1.967 -3.043
+1.950 -5.293
+2.030 -4.401
+2.284 -5.294
+2.534 -2.912
+2.531 -3.680
+2.765 -3.975
+2.913 -7.202
+3.168 -4.317
+3.119 -7.854
+3.385 -5.167
+3.488 -5.135
+3.538 -6.754
+3.900 -6.662
+3.970 -6.080
--- a/notebooks/Model Comparison/MC_3_1.ipynb
+++ b/notebooks/Model Comparison/MC_3_1.ipynb
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "id": "2bc64f1b-8f55-40b3-bfe2-38e144d482f1",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import arviz as az\n",
+    "import matplotlib.pyplot as plt\n",
+    "from matplotlib.ticker import FormatStrFormatter\n",
+    "import numpy as np\n",
+    "import pandas as pd\n",
+    "import pymc as pm"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 13,
+   "id": "3677ec51-54fd-48cd-b379-ce7161c0c1f9",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/html": [
+       "<div>\n",
+       "<style scoped>\n",
+       "    .dataframe tbody tr th:only-of-type {\n",
+       "        vertical-align: middle;\n",
+       "    }\n",
+       "\n",
+       "    .dataframe tbody tr th {\n",
+       "        vertical-align: top;\n",
+       "    }\n",
+       "\n",
+       "    .dataframe thead th {\n",
+       "        text-align: right;\n",
+       "    }\n",
+       "</style>\n",
+       "<table border=\"1\" class=\"dataframe\">\n",
+       "  <thead>\n",
+       "    <tr style=\"text-align: right;\">\n",
+       "      <th></th>\n",
+       "      <th>0</th>\n",
+       "    </tr>\n",
+       "  </thead>\n",
+       "  <tbody>\n",
+       "    <tr>\n",
+       "      <th>0</th>\n",
+       "      <td>51.06</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>1</th>\n",
+       "      <td>55.12</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>2</th>\n",
+       "      <td>53.73</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>3</th>\n",
+       "      <td>50.24</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>4</th>\n",
+       "      <td>52.05</td>\n",
+       "    </tr>\n",
+       "  </tbody>\n",
+       "</table>\n",
+       "</div>"
+      ],
+      "text/plain": [
+       "       0\n",
+       "0  51.06\n",
+       "1  55.12\n",
+       "2  53.73\n",
+       "3  50.24\n",
+       "4  52.05"
+      ]
+     },
+     "execution_count": 13,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "df = pd.read_csv(\"./Daten/chemical_shifts.csv\",header=None)\n",
+    "#df = df.iloc[:,0]\n",
+    "\n",
+    "az.style.use(\"arviz-darkgrid\")\n",
+    "\n",
+    "df.head()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 24,
+   "id": "2db6d34d-e1f4-4e8b-b001-fac7c7ab8910",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "Auto-assigning NUTS sampler...\n",
+      "Initializing NUTS using jitter+adapt_diag...\n",
+      "Multiprocess sampling (4 chains in 4 jobs)\n",
+      "NUTS: [μ, σ]\n"
+     ]
+    },
+    {
+     "data": {
+      "text/html": [
+       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\">/opt/conda/lib/python3.10/site-packages/rich/live.py:231: UserWarning: install \"ipywidgets\" for Jupyter support\n",
+       "  warnings.warn('install \"ipywidgets\" for Jupyter support')\n",
+       "</pre>\n"
+      ],
+      "text/plain": [
+       "/opt/conda/lib/python3.10/site-packages/rich/live.py:231: UserWarning: install \"ipywidgets\" for Jupyter support\n",
+       "  warnings.warn('install \"ipywidgets\" for Jupyter support')\n"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    },
+    {
+     "data": {
+      "text/html": [
+       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"></pre>\n"
+      ],
+      "text/plain": []
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    },
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 3 seconds.\n",
+      "Sampling: [y]\n"
+     ]
+    },
+    {
+     "data": {
+      "text/html": [
+       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\">/opt/conda/lib/python3.10/site-packages/rich/live.py:231: UserWarning: install \"ipywidgets\" for Jupyter support\n",
+       "  warnings.warn('install \"ipywidgets\" for Jupyter support')\n",
+       "</pre>\n"
+      ],
+      "text/plain": [
+       "/opt/conda/lib/python3.10/site-packages/rich/live.py:231: UserWarning: install \"ipywidgets\" for Jupyter support\n",
+       "  warnings.warn('install \"ipywidgets\" for Jupyter support')\n"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    },
+    {
+     "data": {
+      "text/html": [
+       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"></pre>\n"
+      ],
+      "text/plain": []
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    },
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "Auto-assigning NUTS sampler...\n",
+      "Initializing NUTS using jitter+adapt_diag...\n",
+      "Multiprocess sampling (4 chains in 4 jobs)\n",
+      "NUTS: [μ, σ, ν]\n"
+     ]
+    },
+    {
+     "data": {
+      "text/html": [
+       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\">/opt/conda/lib/python3.10/site-packages/rich/live.py:231: UserWarning: install \"ipywidgets\" for Jupyter support\n",
+       "  warnings.warn('install \"ipywidgets\" for Jupyter support')\n",
+       "</pre>\n"
+      ],
+      "text/plain": [
+       "/opt/conda/lib/python3.10/site-packages/rich/live.py:231: UserWarning: install \"ipywidgets\" for Jupyter support\n",
+       "  warnings.warn('install \"ipywidgets\" for Jupyter support')\n"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    },
+    {
+     "data": {
+      "text/html": [
+       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"></pre>\n"
+      ],
+      "text/plain": []
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    },
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 4 seconds.\n",
+      "Sampling: [y]\n"
+     ]
+    },
+    {
+     "data": {
+      "text/html": [
+       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\">/opt/conda/lib/python3.10/site-packages/rich/live.py:231: UserWarning: install \"ipywidgets\" for Jupyter support\n",
+       "  warnings.warn('install \"ipywidgets\" for Jupyter support')\n",
+       "</pre>\n"
+      ],
+      "text/plain": [
+       "/opt/conda/lib/python3.10/site-packages/rich/live.py:231: UserWarning: install \"ipywidgets\" for Jupyter support\n",
+       "  warnings.warn('install \"ipywidgets\" for Jupyter support')\n"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    },
+    {
+     "data": {
+      "text/html": [
+       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"></pre>\n"
+      ],
+      "text/plain": []
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "with pm.Model() as model_g:\n",
+    "    μ = pm.Uniform('μ', lower=40, upper=70)\n",
+    "    σ = pm.HalfNormal('σ', sigma=10)\n",
+    "    y = pm.Normal('y', mu=μ, sigma=σ, observed=df)\n",
+    "    trace_g = pm.sample(2000, idata_kwargs={\"log_likelihood\": True}, random_seed=4591)\n",
+    "    trace_g.extend(pm.sample_posterior_predictive(trace_g, random_seed=4591))\n",
+    "    \n",
+    "with pm.Model() as model_t:\n",
+    "    μ = pm.Uniform('μ', 40, 75)\n",
+    "    σ = pm.HalfNormal('σ', sigma=10)\n",
+    "    ν = pm.Exponential('ν', 1/30)\n",
+    "    y = pm.StudentT('y', mu=μ, sigma=σ, nu=ν, observed=df)\n",
+    "    trace_t = pm.sample(2000, idata_kwargs={\"log_likelihood\": True}, random_seed=4591)\n",
+    "    trace_t.extend(pm.sample_posterior_predictive(trace_t, random_seed=4591))\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 26,
+   "id": "0e05f801-5548-4533-9c46-7ffc290582e9",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "/opt/conda/lib/python3.10/site-packages/arviz/stats/stats.py:792: UserWarning: Estimated shape parameter of Pareto distribution is greater than 0.70 for one or more samples. You should consider using a more robust model, this is because importance sampling is less likely to work well if the marginal posterior and LOO posterior are very different. This is more likely to happen with a non-robust model and highly influential observations.\n",
+      "  warnings.warn(\n",
+      "/opt/conda/lib/python3.10/site-packages/arviz/stats/stats.py:792: UserWarning: Estimated shape parameter of Pareto distribution is greater than 0.70 for one or more samples. You should consider using a more robust model, this is because importance sampling is less likely to work well if the marginal posterior and LOO posterior are very different. This is more likely to happen with a non-robust model and highly influential observations.\n",
+      "  warnings.warn(\n"
+     ]
+    },
+    {
+     "data": {
+      "text/html": [
+       "<div>\n",
+       "<style scoped>\n",
+       "    .dataframe tbody tr th:only-of-type {\n",
+       "        vertical-align: middle;\n",
+       "    }\n",
+       "\n",
+       "    .dataframe tbody tr th {\n",
+       "        vertical-align: top;\n",
+       "    }\n",
+       "\n",
+       "    .dataframe thead th {\n",
+       "        text-align: right;\n",
+       "    }\n",
+       "</style>\n",
+       "<table border=\"1\" class=\"dataframe\">\n",
+       "  <thead>\n",
+       "    <tr style=\"text-align: right;\">\n",
+       "      <th></th>\n",
+       "      <th>rank</th>\n",
+       "      <th>elpd_loo</th>\n",
+       "      <th>p_loo</th>\n",
+       "      <th>elpd_diff</th>\n",
+       "      <th>weight</th>\n",
+       "      <th>se</th>\n",
+       "      <th>dse</th>\n",
+       "      <th>warning</th>\n",
+       "      <th>scale</th>\n",
+       "    </tr>\n",
+       "  </thead>\n",
+       "  <tbody>\n",
+       "    <tr>\n",
+       "      <th>model_t</th>\n",
+       "      <td>0</td>\n",
+       "      <td>-122.479620</td>\n",
+       "      <td>3.940166</td>\n",
+       "      <td>0.000000</td>\n",
+       "      <td>1.000000e+00</td>\n",
+       "      <td>9.110407</td>\n",
+       "      <td>0.000000</td>\n",
+       "      <td>True</td>\n",
+       "      <td>log</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>model_g</th>\n",
+       "      <td>1</td>\n",
+       "      <td>-131.694721</td>\n",
+       "      <td>5.685827</td>\n",
+       "      <td>9.215101</td>\n",
+       "      <td>1.243450e-14</td>\n",
+       "      <td>12.510035</td>\n",
+       "      <td>5.110462</td>\n",
+       "      <td>True</td>\n",
+       "      <td>log</td>\n",
+       "    </tr>\n",
+       "  </tbody>\n",
+       "</table>\n",
+       "</div>"
+      ],
+      "text/plain": [
+       "         rank    elpd_loo     p_loo  elpd_diff        weight         se  \\\n",
+       "model_t     0 -122.479620  3.940166   0.000000  1.000000e+00   9.110407   \n",
+       "model_g     1 -131.694721  5.685827   9.215101  1.243450e-14  12.510035   \n",
+       "\n",
+       "              dse  warning scale  \n",
+       "model_t  0.000000     True   log  \n",
+       "model_g  5.110462     True   log  "
+      ]
+     },
+     "execution_count": 26,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "cmp_df = az.compare({\"model_g\": trace_g, \"model_t\": trace_t})\n",
+    "# cmp_df.to_markdown()\n",
+    "cmp_df"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 27,
+   "id": "82618723-2f74-4c98-8a0e-011f9889cd6e",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "<Axes: title={'center': 'Model comparison\\nhigher is better'}, xlabel='elpd_loo (log)', ylabel='ranked models'>"
+      ]
+     },
+     "execution_count": 27,
+     "metadata": {},
+     "output_type": "execute_result"
+    },
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmMAAADTCAYAAADNnRQhAAAAP3RFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMS5wb3N0MSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8kixA/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAA6NklEQVR4nO3deVhUdf8//ucMMoPJIrmQaBmi4RaBbCIaCOIaKvhR1MIlUTTN3TT3tNtb7xS3O3NLuzMT8w4s1+4QIjVQUMxU1CBR2VUURRhmYM7vD3+er5NsAwMH8Pm4rrkuzvY+zzfDwItz3uccmSAIAoiIiIhIEnKpAxARERG9yFiMEREREUmIxRgRERGRhFiMEREREUmIxRgRERGRhFiMEREREUmIxRgRERGRhFiMEREREUmIxRgRERGRhFiMEdVx+/btQ1BQUKXW3bx5M2bNmlXDiWrf0qVLsWnTJqljEBHVCBZjRAYSFBQEOzs7nDt3Tmf+P/7xD9jZ2WHfvn0SJav/VqxYgenTp0sdg4ioRrAYIzKg119/HQcPHhSnNRoNjh07hrZt20oXqh4TBAElJSVSxyAiqlEsxogMyM/PD//73/+gUqkAADExMbCzs4OVlZXOehEREejfvz+cnZ0xZswYpKSkiMsyMjIwduxYODo6YuTIkcjIyNDZNjU1FcHBwXBzc0OfPn2wd+/eSueLiYlBQEAAnJyc4OnpifDwcACAWq3GmjVr8Pbbb6NHjx5YsGABHj16JG739Mhe//794eDggKVLlyI3NxeTJk0Sc2ZlZemsv2fPHvj6+sLV1RWLFy+GWq0GAOTn52Py5Mlwd3eHi4sLJk6cqNPHoKAghIaGIigoCA4ODrh48SIWLFiAtWvXAgDu37+PKVOmwMXFBS4uLhg+fDhyc3MBAHfu3MGHH34INzc3eHt7Y+vWrdBqtQCAM2fOwMPDA3v27EHPnj3h7u6OnTt3Vvp7R0RUU1iMERlQs2bN4OjoiMjISABAeHg4/P39ddY5c+YMVq1ahdWrV+O3336Du7s7QkJCxGJl9uzZsLGxQVxcHBYvXoz//ve/4raFhYUYN24cvL29cerUKWzfvh07duzA6dOnK8x26dIlzJw5E9OnT8fZs2cRERGBjh07AgC2bduGuLg4HDhwAD/99BMePHiA5cuX62x/4sQJ7N+/H0eOHMHx48cxYcIETJs2DWfPnsXLL7+Mzz//XGf9I0eOICwsDMeOHUNSUhK++OILAIBWq8XQoUMRFRWFX375Baampvjkk090tg0PD8fChQuRmJiILl266CzbtWsXBEHAr7/+iri4OCxbtgxKpVL83pmamuKXX37B7t278f333+t8/+7fv4/MzExERUVh27Zt2LBhA27evFnh946IqCaxGCMysICAAERERCA3NxeJiYnw9fXVWf7jjz/C398fDg4OUCgUCAkJgUqlwrlz55CRkYELFy5gzpw5UCqV6Nq1K/z8/MRto6Oj0aJFC4wePRrGxsZo164dhg8fjsOHD1eY67vvvoO/vz+8vLxgZGSEl19+GZ07dxYzTZ06FVZWVjAzM8O8efNw7NgxsUAEgIkTJ8LCwgKtW7dGt27d8Oabb8Le3h7GxsYYMGAALl++rLO/SZMmoVmzZmjWrBmmTJkiZjQ3N0f//v3RuHFjNGnSBCEhITh79qzOtkOHDkWnTp0gl8uhUCh0lhkbG+PBgwe4efMmjIyM0LVrVzRp0gRZWVmIj4/HggUL0LhxY7Rt2xbvv/++zmljuVyOmTNnQqFQwN7eHjY2NkhKSqrwe0dEVJMaSR2AqKHx8vLC8uXLsWPHDvj6+opHbZ7Kzs5Gr169xGm5XA5ra2tkZ2ejcePGMDMzg5mZmbjc2tpaLBjS09ORlJQEZ2dncXlJSYnOdFkyMzPh4eFR6rLs7Gy0adNGnG7Tpg1KSkpw9+5dWFtbAwCaN28uLm/cuDFatGihM11QUKDT5tPtAKB169bIzs4G8OTo3j//+U+cPHkSeXl5AICCggKo1Wqx8Hp227+bMGECVCoVpk+fjoKCAgwePBizZs1CdnY2zMzMYGFhUep+AcDCwkKnuCstNxFRbWMxRmRgCoUCAwcOxO7du0u9gtLKygrp6enitFarRWZmJqysrNCyZUs8evQI+fn5MDU1BfCkiHrK2toajo6O2LNnj965WrVqhVu3bpW6zMrKCmlpaeJpy7S0NMjlcp0CTF8ZGRliexkZGeK4uV27diE5ORn79+9Hy5YtcfXqVQwZMgSCIIjbymSyMttt0qQJPvroI3z00Ue4desWJk6cCBsbG/Tq1QuPHj3Cw4cPYW5uDuBJ8fr38XpERHUNT1MS1YApU6Zg9+7dcHR0fG6Zn58fIiIicPHiRWg0GuzYsQPGxsZwcnKCtbU1HBwcEBoaCrVajStXruDQoUPitl5eXkhPT8eBAwegVqtRXFyMa9eu4eLFixVmGj58OA4ePIiYmBiUlJQgNzdXPOLm5+eHLVu2ICcnB/n5+Vi3bh0GDhz43ClCfezcuRO5ubnIzc3F1q1bMWjQIADA48ePYWJiAnNzc+Tl5WHLli16tRsdHY0bN25Aq9XC1NQUjRo1gpGREV555RU4OztjzZo1UKlUuHXrFnbv3o0hQ4ZUuQ9ERLWBxRhRDWjWrBnc3d1LXda9e3fMnz8f8+bNg7u7uzgQ/2nhs3btWvz5559wc3PDihUrMGzYMHHbJk2aYNeuXYiKioKnpyfc3d2xdOlSPH78uMJMb775JtauXYv169fD2dkZAQEBuHr1KgBg8uTJcHV1xbBhw+Dr6wtTU1MsW7asWt+DAQMGIDAwEP369UOHDh0wZcoUAMDYsWOh0Wjg7u6OESNGoEePHnq1e/PmTQQHB6Nbt24YPHgwPDw8xIJr3bp1ePDgATw9PTF27FgMGTIEw4cPr1Y/iIhqmkx49twAEZEB2NnZ4ejRo7C1tZU6ChFRnccjY0REREQSYjFGREREJCGepiQiIiKSEI+MEREREUmIxRhRLfL29savv/5a6rKMjAw4Ojrq3PW+PHZ2djrPtDS04OBgHDhwQO/tajoXEVFDw2KMqI6wtrZGYmJite7tZUg7d+6s9dtCBAUFPXej3PIKWCKihoDFGNELThAElJSUSB2jxhQXF0sdgYioXCzGiGpZcnIyAgIC0K1bN0ycOFF8PmNaWhrs7OxQVFQE4MlpyzFjxsDR0RGjRo1CaGgogoKCdNpKSEhA//794eTkhLlz5+qc4jx58iQCAgLg7OwMf39/JCQkiMuCgoLE9hwcHEq9g/+zR6lu376NMWPGwMnJCW5ubnj33XfL7eNvv/0GX19fuLq6YvHixZXK9dlnnyEhIQGrVq2Co6Mj5s2bh9mzZyMjIwPTpk2Do6Mj1q9fDwBITU1FcHAw3Nzc0KdPH+zdu1dsf/PmzZg2bRo+/vhjODs748svv6z4TSEikpJARLWmd+/ewpAhQ4T09HQhPz9fCAwMFDZs2CAIgiDcvn1beOONNwSVSiUIgiAEBgYKn3zyiaBSqYQrV64IHh4ewnvvvSe29cYbbwjjx48XcnNzhbt37wo+Pj7CgQMHBEEQhKSkJMHV1VWIj48XSkpKhKioKMHV1VW4d++eIAiC8N577wkeHh7ClStXhJKSEqGoqOi5rO+9957w7bffCoIgCLNmzRKWLFkiqNVqQa1WC/Hx8WX28Y033hACAwOFu3fvCnfv3hUCAgLEPlYm19N9Pvs9i4mJEacLCgoET09PYe/evYJarRZSUlIET09P4dSpU4IgCMKmTZuEzp07C4cPHxZKSkqEwsJCPd4hIqLaxyNjRLVszJgxsLa2RpMmTdCvXz9cuXLluXUyMjJw4cIFzJ49G0qlEp06dYKfn99z64WEhMDS0hLNmjWDl5eX2FZYWBj+7//+D87OzpDL5ejduzc6duyoM/Zq6NCh6NSpE+RyeYXj1IyNjXHnzh1kZGTA2NgYzs7O5a4/adIkNGvWDM2aNcOUKVNw+PDhSueqSHR0NFq0aIHRo0fD2NgY7dq1w/Dhw8V9AECXLl0waNAgyOVymJiYVLptIiIpNJI6ANGLpnnz5uLXJiYmKCgoeG6dnJwcmJmZwdTUVJz3yiuv4NKlS+W2dffuXQBAeno6zp49i/3794vLi4uL4eHhIU5bW1tXOvNHH32ETZs2YcyYMWjUqBECAwMxadKkMtd/tu3WrVsjOzu70rkqkp6ejqSkJJ2CsKSkRGdan74REUmNxRhRHdSyZUs8evQI+fn5YkGWlZVV6e1btWqF4OBgfPjhh2WuI5PJKt1es2bN8MknnwAAkpKSMG7cOLz55ptlPgw9IyMDHTt2FL+2srKqdK6KWFtbw9HREXv27ClzHX36RkQkNZ6mJKqDrK2t4eDggA0bNkCtVuPq1as4dOhQpbcPDAzE/v37kZCQAK1WC5VKhbi4OL0KumcdPXoUmZmZAAAzMzPI5XLI5WX/+ti5cydyc3ORm5uLrVu3YtCgQZXK1bx5c9y6dUunrebNm+P27dvitJeXF9LT03HgwAGo1WoUFxfj2rVrpV6EQERUH7AYI6qj1q5di2vXrsHNzQ3Lly+Hn59fpe9B1qVLF6xZswafffYZ3Nzc0Lt3b+zevRtarbZKWS5fvozAwEA4ODjgvffew9ixY+Hm5lbm+gMGDEBgYCD69euHDh06YMqUKZXKNWbMGERFRcHFxQXz588H8GT82c6dO+Hs7IwNGzagSZMm2LVrF6KiouDp6Ql3d3csXboUjx8/rlLfiIikxmdTEtUTn376KQoLC/GPf/xD6ihERGRAPDJGVEddunQJqampEAQBCQkJOHjwIPr27St1LCIiMjAO4Ceqo+7du4cPP/wQubm5aN68OaZOnQpPT0+pYxERkYHxNCURERGRhHiakoiIiEhCLMaIiIiIJMRijIiIiEhCLMaIiIiIJMSrKQ3k/v375S63sLBAXl5eLaWpXQ21b+xX/VKb/SoqKkJiYiIcHR2hVCprfH98z2qepaWl1BHoBcYjY7WkvEfH1HcNtW/sV/1Sm/3SaDRISkqCRqOplf3xPSNq2HhkjIhIT6amphg/frzUMYiogeC/JUREeiopKUFubi5KSkqkjkJEDQCLMSIiPRUUFOCbb75BQUGB1FGIqAFgMUZEREQkIRZjRERERBJiMUZEREQkIRZjRER6UiqV8PHxqZV7jBFRw8dbWxAR6UmhUKBLly5SxyCiBoJHxoiI9FRYWIjjx4+jsLBQ6ihE1ACwGCMi0lNxcTGuX7+O4uJiqaNQPXP58mVERkbiypUrUkehOoSnKYmIiGrYyZMnMXv2bCQkJMDExAQqlQouLi5Yt24devXqZbD9BAUF4cKFC2jU6P/9ebexsUF4eDjs7Oxw9OhR2Nra6myTlpYGHx8fvPTSSwAAc3NzDBo0CHPmzIGRkZFOm3K5HNbW1vDx8UFwcDBMTU0Nlv1FxiNjRERENejkyZPo06cPHB0dkZycjMLCQiQnJ8PBwQF9+vTByZMnDbq/hQsXIjExUXyFh4dXaru4uDgkJiZi586dOHjwIA4cOPBcm/Hx8VizZg0uXryIUaNG8VS9gbAYIyLSk0wmQ5MmTSCTyaSOQvXAnDlzMHbsWGzfvl08KmVra4vt27dj7NixmDt3rsQJdXXo0AHOzs64fv36c8vkcjk6d+6MjRs3Ijc3t9KFHpWPxRgRkZ5MTU0xYcIEnqKhMuXl5eH27duIjIxEfHw85s+fX+p68+fPx9mzZ3HixAnk5eXVcsrSXb9+HQkJCejcuXOZ65iZmaFHjx6Ij4+vxWQNF4sxIiI9lZSU4MGDB3xQOJVp/fr1eO211+Dr6wtjY+Pnxmk9ZWtrC2NjY/Tp0wfr1683yL5Xr14NZ2dn8VVWIfh3Hh4ecHFxwbRp0zBq1CgMGzas3PVbtmxZZwrI+o4D+ImI9FRQUICvv/4a48ePh5mZmdRxqA6aNWsWJkyYgOvXr6NPnz5ISUkptSBLSUmBRqNBZGQknJ2dDbLvBQsWYNSoUXpvd/r0ab1uZJyTkwMLCwu990PPq/aRsdu3byM9Pd0QWYiIiBoECwsLvPrqq/Dx8YGLiwvWrFlT6npr1qyBq6srfHx86lVhk5+fj9jYWLi4uEgdpUHQuxibN28eEhMTAQDff/89BgwYgAEDBnAQHxERUSnWrVuH//znP5g0aRJSUlIAPDkiNmnSJPznP//B2rVray2LRqNBUVGR+NL3XnmCIODq1auYOXMmmjZtioCAgBpK+mLRuxg7ffq0+BiQr776Crt27cL+/fuxbds2g4cjIiKq73r16oXIyEhcuHAB7du3h4mJCdq3b4/ff/8dkZGRBr3PGACsWrUKjo6O4svDw0NcNmTIENjb24uvzZs369Wms7Mz5s2bhy5dumDfvn1o3LixQbO/qPQeM1ZUVASFQoHs7Gzcu3cPrq6uAIA7d+4YPBwRUV2kVCrh7e3NB4VTpfXq1Qtnz57FlStXkJGRAWtr63KvVqyqPXv2lLns2rVrVVpWXptkGHoXY7a2tti2bRvS09PFajs3N5fVMRG9MBQKBbp27Sp1DKqHOnfuXCNFGNVvep+mXLp0KaKjo3Hjxg1Mnz4dAPDrr7/qHAYlImrICgsL8dNPP/Hu40RkEHofGevatSvCwsJ05g0dOhRDhw41VCYiojqtuLgY165dQ48ePaSOQkQNQKWKsVOnTlWqsZ49e1YrDBEREdGLplLF2NKlSytcRyaT4cSJE9UORERU112/fh0pKSm4fv06nJycpI5DRPVcpYqxqKioms5BRFTnxcbGYvHixUhMTIRSqcSePXvQrVs3rFy5Eu7u7lLHI6J6qkp34C8uLkZCQgKOHj0K4MmjQQoKCgwarLomTJhQqRvRhoeHY8KECeWuk5SUhB9++MFQ0YioHoqNjYW/vz+cnJyQnJwMlUqF5ORkdOvWDf7+/oiNjZU6IhHVU3oP4E9OTsbkyZMhk8lw9+5dDBw4EGfOnMGPP/5osIec1jVJSUk4cuQIhgwZInUUIpLIkiVLMGbMGOzYsUOcZ2trK04vXboUP//8s1TxiKge07sYW758OSZNmoQRI0aIz6RydXWt1Liy0tjZ2WHGjBmIiYlBVlYWpk2bBmNjY4SFhSE7OxshISEYOXIkgCdF0aeffoqHDx9CJpNh8uTJGDhwIADg999/x/Lly6HRaGBrawuVSiXu4/Hjx1i7di3++OMPaDQa2NjYYPny5WjatGmF+bKzs7Fp0yY8fPgQQ4YMQZs2bfD5559Xqa9E1fXw4UM8fPiwWm3k5+fjwYMHhglUh9Rkv5KTk3H+/Hl89913pS5fsGAB2rdvj5iYmFIfBl1dfM/KZ25uDnNz8+oHIpKKoCdnZ2dBq9UKgiAILi4uOvOr4o033hC2bNkiCIIgJCUlCV27dhU2b94sCIIg/PXXX4K9vb2gUqmEoqIioXfv3sLx48cFQRCE1NRUwcXFRUhOThaKiooET09P4eeffxYEQRB+//13oVOnTsL3338vCIIgLFmyRPjmm2/EfW7evFlYvHixIAiC8P333wvvv/9+uRkrs05JSUkVek+kn2XLlgkA+JLgZWxsXO57Y2xsLHnGF/W1bNmy2vkAUrUkJycLb7zxhjg9YcIE4bvvvhOnw8LCBA8PD8HBwUG4deuWcOPGDcHf319wcHAQ/v3vf0sRudbofWTslVdeQUpKCtq3by/Ou3r1Klq3bq1vUyI/Pz8AQMeOHSEIAgYMGAAAsLGxgUKhQFZWFlQqFYqKitCvXz8AQNu2bdGjRw/ExsZCo9FAo9GgT58+AAB7e3vx+ZkAEBkZiXPnzon/1Wo0GrRo0aLKeUuTl5dX7nJLS0vcv3/foPusKxpq3+piv8aPH49hw4ZVq42mTZs2yKMsNdmvlJQU+Pv7IyUlpdQjXykpKdBoNIiIiKiRI2N8z8pnbm5e7c+qpaVltXNURmZmJg4dOgQ/Pz+0atXK4O0HBQXhwoULaNTo//15t7GxEcdQ29nZ4ejRo8/9nKalpcHHxwcvvfQSgCff00GDBmHOnDkwMjLSaVcul8Pa2ho+Pj4IDg6GqalplbLu3LlT/Fqj0WDVqlXYs2cP7O3tAQCLFy/GW2+9Vanx3/Wd3sVYcHAwQkJCMHHiRBQXF+PgwYPYvn07pk2bVuUQzz7fTSaT6UzL5XIUFxdDJpOVuX1py56dp9VqsW7dOnTs2LHKGYnqAkOcjrG0tKzyL8+6rCb71aZNG3Tr1g2rV6/WGTP21OrVq+Hk5ARPT88a2T/fs4YjOjoamZmZiI6OxujRo2tkHwsXLsSoUaOqtG1cXByUSiX+/PNPjB07Fq+99po4VOhpu1qtFlevXsXatWsxatQofPfdd9V+JOK9e/egUqlgZ2cnzktLSxMPwOiruLhYpyCt6/S+mnLIkCFYtGgRoqKi0KpVKxw+fBhz5swRx27VFBsbGyiVSnGA7O3btxEbGwt3d3fY2NjA2NhYvAXHpUuXcOnSJXHbvn37Yvv27VCr1QCePMqkvIei/p2ZmRny8/MN2Bsiqm9WrlyJr7/+GhMnTkRKSgqAJ0fEJk6ciK+//horVqyQOCHVdZmZmfjzzz8BAH/++ScyMzMlTlS2Dh06wNnZGdevX39umVwuR+fOnbFx40bk5uaWeeSqqKgIixYtgqurK/r27fvcFcdBQUHYt28fUlJS0L9/fwBA9+7dERAQgHfffRdnzpzBqlWr4OjoiIsXL0KtViM0NBTe3t5wc3PDnDlzxLNSaWlpsLOzQ3h4OLy9vTF48GAAwMmTJxEQEABnZ2f4+/sjISFBZ/8bNmxAUFAQHB0dMXLkSKSnp4vL//rrL0ycOBFubm5wc3PT+YyX125VVOnWFt7e3ti+fTuOHDmCnTt3wsfHp1ohKsPY2Biff/45du/eDT8/P0ydOhXLli2Dra0tFAoFNm7ciM2bN+Odd97Bjh074OjoKG67YMECNG/eHMOGDYOfnx9GjBiBy5cvV3rf7u7uMDIyEvdLRC8ed3d3RERE4Pz582jfvj1MTEzQvn17JCYmIiIigvcZIx0qlQp5eXk6r//973/iWRuZTIb//e9/OsufvfBMatevX0dCQkK5DzU3MzNDjx49EB8fX+ryLVu24Pr16zh27Bj27duHQ4cOlbqera0tDh8+DODJkbnw8HDs3bsXzs7OWLhwIRITE2Fvb4/Q0FBcvnwZBw4cQExMDIyNjZ/7J+jkyZP48ccfER4ejqtXr2Lu3LlYuHAhzp49i+nTp2Pq1KnIzc0V1//hhx+wZMkSnDlzBi1btsTGjRsBPLnwb/z48ejWrRtiYmIQExMjHnSqTLv6qtQxvIMHD1aqsao8n/LvR6j++OMPnekzZ86IX3fq1Anffvttqe289dZbiIiIKHXZSy+9hIULF5a6LCAgAAEBAeVmNDU1LXO/RPTicHd3R2RkJH7//XdERkaiT58+eOutt6SORXVQXFwcYmJiylwuCAJSU1OxYcMGcZ6npye8vLyqve/Vq1dj3bp14rSPjw/WrFlTqW09PDwgk8lgaWmJUaNGVThGtWXLlrhy5Uqpyw4fPoxFixahWbNmAICJEydW+YCGIAgICwvD999/L7Y3Y8YM+Pr64l//+pe43rRp08RT32FhYfi///s/ODs7AwB69+6Njh074tdffxXrlYCAALzxxhsAgHfeeQebNm0CAPzyyy8wMzPDlClTxLaftlOZdvVVqWJs3759OtOXLl2Cubk5rKyskJ2djUePHqFr1658WDgRvRDatWsHY2NjtGvXTuooVEd1795d5wzNwYMHcfPmTQiCIM6TyWRo27at+Lfz2fHS1bFgwYIqjxk7ffq0XjlycnJgYWFR5jJra2txujoX+uXm5qKwsBCBgYE685/e87S0faSnp+Ps2bPYv3+/OK+4uBgeHh7idPPmzcWvTUxMxBvYZ2Rk4LXXXis1S2Xa1VelirFnd7hmzRp4eXkhJCQEcrkcWq0W27dvr/BqwrouJiYGoaGhz8339fWt1sUJRET04jExMYGJiQmAJ2PFUlNTn1vn6dGxgoKCGrmysqbl5+cjNjZW5+jRs1q2bImMjAzx4rmMjIwq78vS0hImJiY4ePAg2rRp89zytLQ0ALoX77Vq1QrBwcH48MMP9d6ftbU1fvzxx1KXVafdsug9Ziw8PByTJk2CXP5kU7lcjuDg4Hp/6amnpyd++OGH514sxIiIqDqio6PLvCOATCZDdHR0rebRaDQoKioSX8XFxXptLwgCrl69ipkzZ6Jp06ZlDvUZOHAgtm3bhtzcXOTm5pZ6JXJlyeVyBAYG4p///CdycnIAPLkCMzIyssxtAgMDsX//fiQkJECr1UKlUiEuLg5ZWVkV7s/T0xMPHjzA9u3boVKpoFKpxEH61Wm3zP7pu4GZmdlzg/XOnTv3wl2eTEQvLplMBhMTk3JvuUMEPClccnJydE5P6rO8Kp5egfj09ffTZ0OGDIG9vb342rx5s17tOjs7Y968eejSpQv27dtX5m0tpk6dinbt2qFfv34YOXIk3nnnnWr1a+7cuejYsSPeffdd8erHv48zf1aXLl2wZs0afPbZZ3Bzc0Pv3r2xe/duaLXaCvdlamqK3bt3Iy4uDm+//Ta8vLxw/PjxardbFpmg50/AkSNHsGjRIrz99tuwtrZGRkYGTp48iU8//RSDBg2qcpD6rqIbDtbFG4gaSkPtG/tVvzTUfgENt291qV81edPXpzctL4tSqRRPadKLSe87og0aNAgdO3bE8ePHkZOTIz5bsibuOk1EVBdptVo8fvwYTZo0EYdsEJXl2fFjRKWp0u1pbW1teb8tInphPX78GLt378b48eNhZmYmdRwiquf0/pdOo9Fgw4YN6N27N7p27YrevXtjw4YN4t3tiYiIiKjy9D4ytm7dOpw/fx4rVqxA69atkZaWhi1btqCwsBAff/xxTWQkIiIiarD0LsaOHz+ucwfcdu3aoXPnzggICGAxRkRERKQnvU9TFhcXP3d3XqVSiZKSEoOFIiKqy5RKJTw9PQ12x3QierHpXYy9/fbbmDVrFq5fv478/Hxcu3YNc+fOhaenZ03kIyKqcxQKBd566y0oFAqpoxBRA6B3MbZw4UI0a9YMw4YNg4uLC4YPHw5LS8syH8RNRNTQqFQqREZGQqVSSR2FiBoAvceMmZqaYvXq1Vi1ahXu378PS0tL3meHiF4oGo0GV65cgZubG+8fRUTVVqX7jKnVaty6dQsFBQVIT08X59vb2xssGBEREdGLQO9i7OjRo1i2bBlUKpXOf4QymQxnz541aDgiIiKihk7vYuyf//wnVq5cif79+9dEHiKiOk8mk0GpVPJB4URkEHoXY4IgwNfXtyayEBHVC6ampggJCZE6BhE1EHqPvB83bhy2bt0KQRBqIg8RUZ2n1WqRn58PrVYrdRQiagD0PjLWr18/jB8/Hl9++SUsLS11lp04ccJgwYiI6io+KJyIDEnvYmzGjBl48803MXDgQF7STURERFRNehdjN27cwIEDB2BkZFQTeYiIiIheKHqPGevRowcuX75cE1mIiIiIXjh6Hxlr0aIFJk6cCG9vb7Ro0UJn2ezZsw0WjIiorlIqlejVqxcfFE5EBqF3MVZUVARvb28AwJ07dwweiIiorlMoFHB0dJQ6BhE1EFW66SsR0YtMpVLht99+Q48ePXghExFVG5/wTUSkJ41Gg0uXLkGj0UgdhYgaABZjRERERBJiMUZEREQkIRZjRERVoFAopI5ARA1EpQbwx8fHV6oxFxeXaoUhIqoPzMzMMHnyZKljEFEDUalibNasWTrTDx48QElJCUxNTZGfnw8jIyM0bdoUp06dqpGQRER1iVarRWFhIRo3bgy5nCcYiKh6KlWMPVtk7d69G6mpqZg7dy7MzMzw8OFDhIaG4vXXX6+pjEREdQofFE5EhqT3fcZ27tyJ6OhocbyEubk5Pv74Y3h7e2PcuHGGzkdERETUoOl9fN3IyAh//fWXzrwbN27wweFEREREVaD3kbH3338f48aNg7+/P6ytrZGRkYGDBw8iJCSkJvIRERERNWh6F2Pjxo1Dhw4dcOTIEfz5559o2bIlPvvsM/Ts2bMm8hER1TkKhQI9e/bk7S2IyCD0LsYAwMPDAx4eHobOQkRlyMnJwYkTJ+Dj44OWLVtKHeeFp1Qq0a1bN6ljEFEDofeYMbVajQ0bNqBPnz5wcnICAJw8eRLffPONwcMR0RNxcXG4c+cO4uLipI5CePKg8OjoaKhUKqmjEFEDoHcxtmbNGly6dAmrV6+GTCYDALRv3x5hYWEGD0dET46KpaamAgBSU1ORk5MjbSCCRqPBH3/8wQeFE5FB6H2a8qeffsLx48dhamoq3uywVatWyMrKMng4In0UFRVBrVYbrD25XI5Hjx4ZrL2qOnXqFGQyGQRBgEwmw6lTp+Dr61vl9upKv56lUCigVCqljkFEJAm9izEjI6Pn7jj98OFDmJubGyxUfWRhYVHhnbgtLS1rKU3tqwt9++WXXxATEyN1jBolCALS0tKwe/duqaMYlKenJ7y8vKrdTm39HD79rFtYWMDCwqJW9lkXPmM1oaH2i0gfehdjnp6e+PTTT7F48WIAQElJCUJDQ+Ht7W3wcPVJXl5eucstLS1x//79WkpTu+pK3zp27Ih27doZrD0LC4sK39ea9vPPPyM9PR2CIIjzZDIZWrduXeWjY3WhX3+nUCiq/TNUmz+Hjx49QqNGjZCXlwetVlvj+6srnzFDq0v9YlFIUtK7GPvoo4/w8ccfw9XVFSUlJXBwcICXlxdWr15dE/mIKk2pVBr0VJeFhUWt/KEtS05ODtLS0p6b//ToWGFhYZWurJS6Xw2BmZkZPvjgA6ljEFEDoXcxJpPJsHnzZuTm5iItLQ3W1tZo3rw5MjIy0KRJk5rISPRCiouLE8eK/Z1MJkNcXBwGDx4sQTLSarVQqVQwMTHhg8KJqNr0/i3ywQcfQKPR4OWXX4a9vb1YiI0ZM6Ym8hG9kARBwL1790otxCqznGrW48ePsXPnTjx+/FjqKETUAOh9ZKxdu3aYNWsWNm/eDJlMhrS0NIwZMwZjx46tiXxELySZTIbRo0eXe3WoQqEQby9DRET1l95HxpYuXYpGjRph8eLFYiE2btw4FmNEBqZUKmFmZlbmi7eCICJqGPQuxmQyGT777DNkZmbCz88P77//Pk9REhEREVVRpU5Tzp49+7nTIQqFAo0bN0ZiYiISExMBAOvWrTN8QiKiOkahUMDDw4MPCicig6hUMVbWvZu6du1q0DBERPWBUqkUn81LRFRdlSrGpk2bVtM5iKgcSUlJyM7OxiuvvIKOHTtKHeeFp1KpEBcXh+7du8PExETqOERUz+l9NSUA3L59G0lJSSgoKNCZP3ToUENkIqL/X2xsLBYvXozExESYmJhApVKhW7duWLlyJdzd3aWO98LSaDS4ePEinJycWIwRUbXpXYx9+eWXWL9+Pdq1a4fGjRuL82UyGYsxIgOKjY2Fv78/xowZgwMHDsDW1hYpKSlYvXo1/P39ERERwYKMiKgB0LsY27VrF7799lvY29vXRJ5qmzBhAgYNGoSAgIBy1wsPD8eRI0fw5Zdf1lIyIv0sWbIEY8aMwY4dO8R5tra24vTSpUvx888/SxWPiIgMpEqPQ+rcuXNNZHnhPXz4EA8fPpQ6ht7y8/Px4MEDqWMYnJT9Sk5Oxvnz5/Hdd9+VunzBggVo3749YmJiYGtrq1fbfL+qr6CgAHK5HJmZmbXy0PXK9M3c3Bzm5uY1noWIDE/vYmzy5MnYuHEjZs6cCSMjoyrv2M7ODjNmzEBMTAyysrIwbdo0GBsbIywsDNnZ2QgJCcHIkSMBPBm8/Omnn+Lhw4eQyWSYPHkyBg4cCAD4/fffsXz5cmg0Gtja2kKlUon7ePz4MdauXYs//vgDGo0GNjY2WL58OZo2bVqpjI8fP8bSpUtx6dIlWFpaomvXrrh27Rr27Nnz3LoWFhYVPqPO0tKy3OUbN27EJ598Uqls1PAZGxuXWWjZ2trC2NgY/v7+tZyKnrV06VKpI4iWLVuG5cuXSx1DbxX9XiR6EehdjG3duhX379/HV199BQsLC51lp06d0qstmUyG/fv34+rVqxg+fDhCQkIQFhaGGzduYOjQofD394dMJsPUqVMxf/589OvXDzdv3sTw4cNhZ2eHV199FTNmzMDixYvRp08fXLx4ESNHjsSwYcMAAGvWrIGdnR2WLVsGAPj3v/+NdevWYeXKlZXKt2XLFhQXF+PYsWPQaDQYN24cGjUq/VtW0X/HlpaWuH//frnrjB8/XsxenzRt2rRBHmmRsl8pKSnw9/dHSkpKqQVZSkoKNBoNIiIi9D4yxver+gRBgEajgbGxca08kqoyfTM3N6/wd0xdU5nfi7WFRSFJSe9iLDQ01GA79/PzAwB07NgRgiBgwIABAAAbGxsoFApkZWVBpVKhqKgI/fr1AwC0bdsWPXr0QGxsLDQaDTQaDfr06QMAsLe3R5cuXcT2IyMjce7cOfFUj0ajQYsWLSqd77fffsOsWbMgl8uhVCoxePBgHD161CB9L019Pc1gaWkJU1NTqWMYnJT9atOmDbp164bVq1frjBl7avXq1XBycoKnp6febfP9qr5Hjx5h9+7dGD9+PMzMzGp8fw31PSOiJ/QuxlxdXQ2282efrSeTyXSm5XI5iouLy/2vs7Rlz87TarVYt25dte7L9Gx7fCgz1aaVK1eKpyEXLFigczXl119/jYiICIkTEhGRIej9bEoAuHTpEr766iusX78eoaGh4qsm2NjYQKlUileN3b59G7GxsXB3d4eNjQ2MjY0RFRUl5rp06ZK4bd++fbF9+3ao1WoAQGFhIa5du1bpfbu7uyMiIgJarRZqtRqHDh0yYM+Iyvf05+/8+fNo3749TExM0L59eyQmJvK2FkREDYjeR8b27t2Lzz77DD179kRMTAw8PT1x+vRpeHt710Q+GBsb4/PPP8fKlSuxadMmyGQyLFu2TBwns3HjRixfvhyhoaGwtbWFo6OjuO2CBQuwYcMGnXFY48ePh52dXaX2/cEHH2DJkiUYMGAAmjZtCjs7O9y9e9ewHSQqh7u7OyIjI3H16lVkZWXxDvxERA2QTBAEQZ8NfH19xfEqLi4uiI+Px6lTp3Do0CGsWbOmpnJKori4GMXFxTAxMYFarcaHH34IFxcXBAcHP7duRYNQ69JAVUNrqH1jv+qX2uxXUVERLl68CHt7e53hFTWF71nN4wB+kpLeR8bu3bsnPiBXLpdDq9WiZ8+emD17tsHDSa2goADjxo1DSUkJioqK4OrqiqCgIKljEZHElEolXFxcpI5BRA2E3sWYtbU1bt++jVdffRWvv/46fvrpJ1haWsLY2Lgm8tW4mJiYUse7+fr6Ytq0aQgPD5cgFRHVZSqVCmfPnoWrqyufTUlE1aZ3MRYcHIwbN27g1VdfxQcffIAZM2ZAo9Fg0aJFNZGvxnl6elbp9gBE9OLSaDS4cOECHB0dWYwRUbXpXYz17NlTPLfu6emJs2fPQqPRoEmTJgYPR0RERNTQ6XVrC0EQ4OPjozNPoVCwECMiIiKqIr2KMZlMBltbW2RlZdVUHiKieqGiZ9ESEVWW3qcpBw4ciMmTJyMoKAitWrXSuSt9z549DRqOiKguMjMzw7Rp06SOQUQNhN7F2LfffgvgyQPDnyWTyXDixAnDpCIiqsMEQYBarYZCoeBj0oio2vQuxp4+eoiI6EWVn59fqw8KJ6KGjYMeiIiIiCTEYoyIiIhIQizGiIiIiCTEYoyISE8KhQLdu3eHQqGQOgoRNQB6D+AnInrRKZVKuLq6Sh2DiBoIHhkjItJTUVERTp48iaKiIqmjEFEDwGKMiEhParUaiYmJUKvVUkchogaAxRgRERGRhFiMEREREUmIxRgRERGRhHg1JRGRnszMzDB9+nSpYxBRA8EjY0REehIEARqNBoIgSB2FiBoAFmNERHrKz8/HF198gfz8fKmjEFEDwGKMiIiISEIsxoiIiIgkxGKMiIiISEIygSNQiYiIiCTDI2NEREREEmIxRkRERCQhFmNEREREEmIxRkRERCQhPg7JwO7du4c5c+YgKSkJBQUF+OOPP0pdLywsDMuWLUNoaCgGDRoEAFiwYAFOnTqFZs2aAQBatmyJHTt21Fr28lSnXw8ePMDHH3+MlJQUyGQyTJ06FYMHD67N+GWqqF/vv/8+7ty5AyMjIzRq1AgzZ85Ez549AQCbN2/GN998g1deeUVc/4cffqjV/OWpTt/UajWWLVuG+Ph4AMDo0aPx/vvv13ofSlNev7Kzs7Fw4UKkpaVBqVSiefPmWLJkCWxsbADU389YRf2qy58xoOKfxY8//hhxcXHIyMjA/v374eDgIC6ry+8ZkcEIZFCPHj0SYmNjhaSkJKFr166lrpOcnCyMGDFCGDFihHD48GFx/vz584Vt27bVVlS9VKdfS5YsEVatWiUIgiDcvn1bcHV1FTIyMmold0Uq6ldeXp749eXLlwUHBwdBrVYLgiAImzZtEpYsWVJrWfVVnb5t3bpV+PDDDwWtVis8ePBA6N27t5CYmFhb0ctVXr/u3LkjxMfHi9NfffWVMGLECHG6vn7GKupXXf6MCULFP4uxsbFCVlZWqT9ndfk9IzIUnqY0MFNTU3Tv3h2mpqalLler1Vi4cCFWrFgBhUJRy+mqrjr9OnToEN577z0AQJs2bfD222/j2LFjNZ65Mirql7m5ufh1Xl4e5HJ5vXkeYXX6dvjwYYwePRoymQwWFhYYPHgwDh06VCu5K1Jev5o3bw5nZ2dxulu3bkhPT6/NeFVWnX7V5c8YUPHPYvfu3WFlZVXLqYjqDp6mrGXr169H//79YWdnV+ryffv24dChQ2jatCmCg4Ph6elZywmrpqx+PXjwAIWFhbC2thbntWnTBhkZGbUdscoWLVqE2NhY5Ofn49///rdOsRkVFYULFy7gpZdewsiRIzF06FDpglZBWX1LT09H69atxfXatGmDyMhIqWJW2VdffQVfX1+defX1M/asZ/vVED5jFWkI7xlReViM6WnKlCk4f/58qcumT5+Od999t8xtf/vtNyQlJeGjjz4qdfnMmTPRokULGBkZ4ffff0dISAh2796NTp06GSR7eWqyX1KqTr+e+sc//gEAOHPmDFatWoW9e/fC1NQUgYGBCAkJgUKhwF9//YXg4GBYWFigd+/eBu1DWWqyb1IyRL8AIDQ0FOnp6fj000/FefX1M/as0volNUP1rTRSvmdEtYXFmJ6++OKLKm8bGxuL5ORk+Pj4AADu3LmDP//8E+np6Zg0aZLOQPC33noLLi4uOH/+fK380qnJfjVu3BgZGRl49dVXAQBpaWm19ou0Ov36Ozc3NxgZGSEhIQFeXl5o2bKluKxdu3bo27cv4uPja60Yq8m+WVtbIz09Xec9e/ZIWU0yRL/Wrl2LhIQE7Ny5E40bNxbn19fP2FOl9atp06aSfsYAw/4s/p2U7xlRbeGYsVo0Z84cnDp1ClFRUYiKioKDgwOWLFmCSZMmAQAyMzPFdW/duoXExER06dJFqriVVlG/3nnnHXzzzTcAnvyR+PXXX9G/f38pI1fKw4cPkZ2dLU5fuXIFt27dgq2tLQDd9+vevXv49ddf0bVr11rPWRUV9e2dd97Bt99+C0EQkJeXhx9//FG8OrauW716NRITE/Hll18+d5Svvn7GgPL7VV8/Y5VRn98zosrisylrgLe3N9RqNe7cuQMrKyt06NABX3755XPrBQUFYeTIkeIfueHDh0OlUkEul0Mul2PcuHEYMmRIbccvU1X7df/+fZ3L7j/44IM6NbaqrH6lp6djxowZUKlUMDIygomJCaZMmQIvLy8AwLRp05CamgojIyNotVoMGzYM48aNk7Qvf1fVvhUVFWHp0qU4d+4cBEHAqFGjEBwcLG1nnlFWv86dO4fRo0fj9ddfh4mJibh+eHg4jIyM6u1nrKJ+1fXPGFD+74/Zs2cjISEBd+/ehYWFBYyNjXH06FGYmprW+feMyBBYjBERERFJiKcpiYiIiCTEYoyIiIhIQizGiIiIiCTEYoyIiIhIQizGiIiIiCTEYoyIiIhIQizGiIiIiCTEYoyojgsPD8eIESOqvL2HhwfOnDlTo/uoyJgxY3Du3DmD7yskJASnTp0ySFtERFJhMUZENerUqVPQarVwcnIyeNshISEIDQ01eLtERLWJxRgR1ahvv/22xh5f061bNzx69Ah//PFHjbRPRFQbWIwR1QF37tzBzJkz0aNHD3h6emLz5s3QarWlrmtnZ4c9e/bA19cXrq6uWLx4MdRqtbj8q6++Qq9eveDu7o5du3ZVKc/FixcRGBgIJycn+Pn5ITo6WlymVquxZs0avP322+jRowcWLFiAR48eldqORqPB6dOn0b179yrtq6ioCIsWLYKrqyv69u2LvXv3ws7OTmd7Nzc3nW2IiOobFmNEEtNqtZgyZQratm2L6OhofPfddzhx4gT++9//lrnNkSNHEBYWhmPHjiEpKQlffPEFAOD06dP44osvsHXrVvzyyy9ITU3F/fv39cqTl5eH4OBgDBs2DGfOnMH8+fMxa9YspKSkAAC2bduGuLg4HDhwAD/99BMePHiA5cuXl9rWzZs3UVJSgldffbVK+9qyZQuSk5Nx/PhxhIWF4ejRo8+1YWtri6SkJL36SERUl7AYI5LYpUuXkJWVhZkzZ0KpVMLKygrjxo3D4cOHy9xm0qRJaNasGZo1a4YpU6aI6x4+fBj+/v7o0qULlEol5s6dW+YRtrL88ssvsLa2xogRI9CoUSP07NkTvXv3xqFDhwAAP/74I6ZOnQorKyuYmZlh3rx5OHbsmM7Ruafy8vLQpEmTKu/r8OHDmDx5Ml5++WW8/PLLCA4Ofq6NJk2a4OHDh3r1kYioLmkkdQCiF11aWhpyc3Ph4uIiztNqtWjVqlWZ21hbW4tft27dGtnZ2QCAnJwcdOzYUVxmbm4OU1NTvfJkZ2ejdevWOvOe3Ud2djbatGkjLmvTpg1KSkpw9+5dnVwAYGFhgcePH1d5Xzk5OTrfh9K+J48fP4a5uXkle0dEVPewGCOSmLW1NV555RVERUVVepuMjAyx6MrIyICVlRUAoGXLlsjMzBTXe/jwIfLz8/XKY2VlhfT0dJ156enpaNu2rbg8LS1N3H9aWhrkcjmaN2/+XFuvvfYajIyMcPv27VJPVVa0r6f9ebqvZ/v2VEpKCjp16qRXH4mI6hKepiSS2JtvvglLS0t8/vnnKCgogFarRWpqKs6ePVvmNjt37kRubi5yc3OxdetWDBo0CAAwcOBAREREICkpCUVFRQgNDYVcrt/H3NPTE+np6fjvf/+L4uJi/Pbbb4iOjoafnx8AwM/PD1u2bEFOTg7y8/Oxbt06DBw4EAqF4rm2FAoFevToUeZ9zira18CBA7Ft2zaxr6VdkBAfHw8vLy+9+khEVJewGCOSmJGREbZu3YqbN2+ib9++cHFxwcyZM3Hnzp0ytxkwYAACAwPRr18/dOjQAVOmTAEA9OrVCyEhIZg0aRK8vLzw2muvwdLSUq88TZs2xfbt27F//364ublh1apVWLt2LWxtbQEAkydPhqurK4YNGwZfX1+Ymppi2bJlZbY3evRo/PDDD1Xa19SpU2FjY4N+/fph5MiR8PHxgbGxsbj9hQsX8NJLL8He3l6vPhIR1SUyQRAEqUMQUeXZ2dnh6NGjYsFSHwQFBWHmzJnVvvHrzz//jH/961/4+eefATwpDN9991306tXLEDGJiCTBMWNEVOP27NlTpe3u3LmD1NRUODk5ITMzE1u2bIGvr6+4fOvWrYaKSEQkGRZjRC+I4OBg8fmQz/Lz88OKFSskSFSxkpISrFixAmlpaXjppZfg7e2NqVOnSh2LiMigeJqSiIiISEIcwE9EREQkIRZjRERERBJiMUZEREQkIRZjRERERBJiMUZEREQkIRZjRERERBJiMUZEREQkIRZjRERERBL6/wB7spVteld8bQAAAABJRU5ErkJggg==",
+      "text/plain": [
+       "<Figure size 600x200 with 1 Axes>"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "az.plot_compare(cmp_df)"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.10.11"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
+%% Cell type:code id:2bc64f1b-8f55-40b3-bfe2-38e144d482f1 tags:
+
+``` python
+import arviz as az
+import matplotlib.pyplot as plt
+from matplotlib.ticker import FormatStrFormatter
+import numpy as np
+import pandas as pd
+import pymc as pm
+```
+
+%% Cell type:code id:3677ec51-54fd-48cd-b379-ce7161c0c1f9 tags:
+
+``` python
+df = pd.read_csv("./Daten/chemical_shifts.csv",header=None)
+#df = df.iloc[:,0]
+
+az.style.use("arviz-darkgrid")
+
+df.head()
+```
+
+%% Output
+
+           0
+    0  51.06
+    1  55.12
+    2  53.73
+    3  50.24
+    4  52.05
+
+%% Cell type:code id:2db6d34d-e1f4-4e8b-b001-fac7c7ab8910 tags:
+
+``` python
+with pm.Model() as model_g:
+    μ = pm.Uniform('μ', lower=40, upper=70)
+    σ = pm.HalfNormal('σ', sigma=10)
+    y = pm.Normal('y', mu=μ, sigma=σ, observed=df)
+    trace_g = pm.sample(2000, idata_kwargs={"log_likelihood": True}, random_seed=4591)
+    trace_g.extend(pm.sample_posterior_predictive(trace_g, random_seed=4591))
+
+with pm.Model() as model_t:
+    μ = pm.Uniform('μ', 40, 75)
+    σ = pm.HalfNormal('σ', sigma=10)
+    ν = pm.Exponential('ν', 1/30)
+    y = pm.StudentT('y', mu=μ, sigma=σ, nu=ν, observed=df)
+    trace_t = pm.sample(2000, idata_kwargs={"log_likelihood": True}, random_seed=4591)
+    trace_t.extend(pm.sample_posterior_predictive(trace_t, random_seed=4591))
+```
+
+%% Output
+
+    Auto-assigning NUTS sampler...
+    Initializing NUTS using jitter+adapt_diag...
+    Multiprocess sampling (4 chains in 4 jobs)
+    NUTS: [μ, σ]
+
+
+
+    Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 3 seconds.
+    Sampling: [y]
+
+
+
+    Auto-assigning NUTS sampler...
+    Initializing NUTS using jitter+adapt_diag...
+    Multiprocess sampling (4 chains in 4 jobs)
+    NUTS: [μ, σ, ν]
+
+
+
+    Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 4 seconds.
+    Sampling: [y]
+
+
+
+%% Cell type:code id:0e05f801-5548-4533-9c46-7ffc290582e9 tags:
+
+``` python
+cmp_df = az.compare({"model_g": trace_g, "model_t": trace_t})
+# cmp_df.to_markdown()
+cmp_df
+```
+
+%% Output
+
+    /opt/conda/lib/python3.10/site-packages/arviz/stats/stats.py:792: UserWarning: Estimated shape parameter of Pareto distribution is greater than 0.70 for one or more samples. You should consider using a more robust model, this is because importance sampling is less likely to work well if the marginal posterior and LOO posterior are very different. This is more likely to happen with a non-robust model and highly influential observations.
+      warnings.warn(
+    /opt/conda/lib/python3.10/site-packages/arviz/stats/stats.py:792: UserWarning: Estimated shape parameter of Pareto distribution is greater than 0.70 for one or more samples. You should consider using a more robust model, this is because importance sampling is less likely to work well if the marginal posterior and LOO posterior are very different. This is more likely to happen with a non-robust model and highly influential observations.
+      warnings.warn(
+
+             rank    elpd_loo     p_loo  elpd_diff        weight         se  \
+    model_t     0 -122.479620  3.940166   0.000000  1.000000e+00   9.110407
+    model_g     1 -131.694721  5.685827   9.215101  1.243450e-14  12.510035
+    
+                  dse  warning scale
+    model_t  0.000000     True   log
+    model_g  5.110462     True   log
+
+%% Cell type:code id:82618723-2f74-4c98-8a0e-011f9889cd6e tags:
+
+``` python
+az.plot_compare(cmp_df)
+```
+
+%% Output
+
+    <Axes: title={'center': 'Model comparison\nhigher is better'}, xlabel='elpd_loo (log)', ylabel='ranked models'>
+
+