{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "Sampling from linear model\n", "=========================" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A dataset consisting of three simultaneously obtained data sets (e.g., different photometric bands or RVs determined in individual spectral orders) is constructed and a model is set up with a single slope for all series but individual means. " ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pylab as plt\n", "from PyAstronomy import funcFit2 as fuf2\n", "import scipy.optimize as sco" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# Definition of model \n", "\n", "class LinMod(fuf2.MBOEv):\n", " \"\"\" Linear model \"\"\"\n", "\n", " def __init__(self, no):\n", " # no is number of measurement series\n", " self.no = no\n", " # 'parnames' specifies parameter names in the model\n", " parnames = [\"mean%d\" % (i+1) for i in range(no)]\n", " parnames.append(\"slope\")\n", " \n", " fuf2.MBOEv.__init__(self, pars=parnames, rootName=\"LinMod\")\n", "\n", " def evaluate(self, x, *args, **kwargs):\n", " \"\"\"\n", " Evaluate model at 'x' here and return y = const + x*slope\n", " \"\"\"\n", " r = np.zeros( (len(x), self.no) )\n", " for i in range(self.no):\n", " r[::,i] = self[\"slope\"] * x + self[\"mean%d\" % (i+1)]\n", " \n", " return r" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "# Generate mock data\n", "np.random.seed(1718)\n", "\n", "# Number of data points and orders\n", "nd = 10\n", "no = 3\n", "# Uncertainty of individual points\n", "s = 0.3\n", "\n", "# Mock data and time (x) axis\n", "d = np.zeros((nd,no))\n", "x = np.arange(nd)\n", "\n", "for i in range(no):\n", " offset = np.random.random() * 20\n", " d[::,i] = 0.6 * x + offset + np.random.normal(0, s, nd)\n", "\n", "# Uncertainty\n", "yerr = np.ones_like(d)*s" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Optimization terminated successfully.\n", " Current function value: 32.650283\n", " Iterations: 272\n", " Function evaluations: 470\n", "-------------------- Parameter summary ---------------------\n", " mean1 = 2.35174, free: T, restricted: F, related: F\n", " mean2 = 3.68232, free: T, restricted: F, related: F\n", " mean3 = 5.94872, free: T, restricted: F, related: F\n", " slope = 0.630558, free: T, restricted: F, related: F\n", "------------------------------------------------------------\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAEICAYAAABGaK+TAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABWtklEQVR4nO3dd1hURxfA4d/sUgVpAkpRAUWjYgV77z2W2GOiRk0s+dR0k2h6UdObJsaSoonGkkSNPUpUNEbsvaIUKSJFetmd74+L2CtlWZj3efZhuXvL2RUPw9yZM0JKiaIoimJ+dKYOQFEURXk4KoEriqKYKZXAFUVRzJRK4IqiKGZKJXBFURQzpRK4oiiKmVIJXFHuQAjxgxDiPRNc10cIIYUQFvex70ghxI7iiEspeVQCVwqNEOK8ECJDCJEihEgSQuwUQowTQtzXz9mDJC5FUVQCVwpfbylleaAqMAN4BZhv2pAUpXRSCVwpElLKZCnlKmAwMEIIEQAghOgphNgvhLgihIgQQrx13WHb8r4mCSFShRDNhRDVhBBbhBCXhRDxQojFQginO11XCPFF3nmvCCH2CiFaX/faW0KI34QQP+X9lXBUCBF03esNhRD78l5bCtjc5TojhRAhQojP8v7aOCeEaJG3PUIIESeEGHHd/o55170khLgghJh29S8TIYReCPFx3vs7B/S86VqOQoj5QohoIUSUEOI9IYT+3v8KSmmnErhSpKSU/wGRwNVEmgY8CTihJarxQoi+ea+1yfvqJKW0l1LuAgTwIeAJ1AIqA2/d5ZJ7gAaAC/ALsEwIcX0ifhRYknf9VcDXAEIIK+AP4Oe8Y5cBj93j7TUFDgEV8q61BGgMVAeGA18LIezz9v0KcAT8gLZ5n8GovNfGAr2AhkAQMOCm6/wA5OadtyHQBRhzj9iUskBKqR7qUSgP4DzQ6Tbb/wVev8MxnwOf5T33ASRgcZdr9AX2P0BMiUD9vOdvAZuve602kJH3vA1wERDXvb4TeO8O5x0JnL7u+7p5sVe8bttltF8meiAbqH3da88AwXnPtwDjrnuty9XPAagIZAG2170+FNh6XRw7TP1vrx6meaibRUpx8AISAIQQTdH6xgMAK8AarbV7W0KIisAXaC348mh/NSbeZf8XgdFoLXYJOACu1+0Sc93zdMAm76apJxAl87Jingv3eF+x1z3PAJBS3rzNPu/6ljed7wLa50LetSPucN2qecdGCyGubtPdtL9SRqkuFKVICSEaoyWqq0PdfkHruqgspXQEvkXrJgEt4d7sg7ztdaWUDmhdE+I2+5HX3/0yMAhwllI6Acl32v8m0YCXuC5LAlXu47j7EQ/koCXj688ddd21K9/huhFoLXBXKaVT3sNBSlmnkGJTzJhK4EqREEI4CCF6ofULL5JSHs57qTyQIKXMFEI0AYZdd9glwIjWT8x1+6cCyUIIL+Clu1y2PFpf8SXAQgjxBloL/H7syjt2khDCUgjRH2hyn8felZTSAPwGvC+EKC+EqAo8DyzK2+W3vOt6CyGcganXHRsNbAQ+yftMdXk3dtsWRmyKeVMJXClsq4UQKWgtx9eBT7l2sw5gAvBO3j5voCUvAKSU6cD7QEjeyI5mwNtAI7SW9F/AyrtcewOwHjiF1g2RyX12NUgps4H+aH3KCWijZ+52rQf1P7QbuOfQ/hr5BViQ99r3aLEfBPbd5rpPonU3HUPrPloOeBRibIqZEjd2+SmKoijmQrXAFUVRzJRK4IqiKGZKJXBFURQzpRK4oiiKmSrWiTyurq7Sx8enOC+pKIpi9vbu3RsvpXS7eXuxJnAfHx9CQ0OL85KKoihmTwhx21nBqgtFURTFTKkEriiKYqbumcCFEAvyahsfuW7bR0KIE0KIQ0KI3+9Wn1lRFEUpGvfTAv8B6HbTtk1AgJSyHtq05VcLOS5FURTlHu6ZwKWU28grBXrdto1Syty8b/8FvIsgNkVRFOUuCqMP/Clg3Z1eFEI8LYQIFUKEXrp0qRAupyiKokABE7gQ4nW0EpyL77SPlHKulDJIShnk5nbLMEZFURTlIT30OHAhxEi0dfw6SlXSUFEUpdg9VAIXQnRDW/mkbV4NZ0VRFOUORq3XSuIv7LawUM97P8MIf0VbraSmECJSCDEabSXv8sAmIcQBIcS3hRqVoiiKck/3bIFLKYfeZvP8IohFURRFeQBqJqaiKEoxiEqNuvdOD0glcEVRlCKwbt06Xn755fzvo9OiC/0aKoEriqIUAikloaGhGI1GAPbv38/SpUs5GXOSi6kXi+SaKoEriqIUglWrVtG4cWOCg4NJyU7Br48fNT6uwYANA7iYpiXwuj/Wpe6PdZl9YHahXLNYV6UPCgqSqh64oiilQXJyMsOGDWPgwIGMHDmSlLQUZi6dSWqVVHbE7CDLkIWPgw99qvdhS/gWDscf5vCIww91LSHEXill0M3bi3VBB0VRFHO2d+9eLl68SO/evXFwcCArO4toQzSz9sxi7bm1XNZfxineif7+/Xm02qPUqVAHIQQhUSFFEo9K4IqiKHeRmpqKvb09ANOnT+fMmTM06dCEdWHrcPifA78k/oLFCQvaebejd7XetPZqjaXe8pbzeNh5FHpsKoEriqLcwZw5c3jllVeIiIjAys6KAdMHEJIUQpcVXTBKI/Xc6jGt6TS6+nTFycbprufysvcq9PhUAlcURclz8eJFZsyYwcSJE6lZsyZNmzVl8EuDeS/0PbbHbictJw0POw/G1B1Db7/e+Dj6mDRelcAVRSnT4uLiSEtLw9fXF51Ox4IFC6jepDob0jaw5twaoqtEYxdrR5eqXehdrTeBFQPRiQcbwFfYNVCuUglcUZRS615FpAwGA/Xr16dt27Z8+8O3BCcF0+vXXnyf8D26IzqaezRnSqMptK/SHlsL2+IM/b6oBK4oSpny+eefs3nzZtasWYMRI5O+nMRpq9O0X9aeXGMu1Z2q80LgC/Tw64F7OXdTh3tXKoErilKqZV7JZM6cOYwdOxYLCwssLC2QlSTvhLzDpohNJGUlUcFYgWGPDKN3td7UdK6JEMLUYd8XlcAVRSl1cnJyMBgMAJyLOseS55fgXt2dhIoJbPDcQLh9OLFhsXSo3IHe1XrT3LM5FjrzS4fmF7GiKMpdxMXFUa9ePaZPn47B10C6SzqDlg/i7YtvIy9KAisGMipgFJ2rdqa8VXlTh1sgKoErimL25s6dS3Z2Ns8++yzu7u48NuIxTnud5mD8QQBSRSoTGkygl18vvMt7mzjawqMSuKIoZic3N5eDBw8SGBgIwPr160lLS6P3k71ZeGQhO2rvwJhszN8/PCWcbw58g1EamdBggqnCLnQqgSuKYnbefvttZs6cycWLF3F1deX1L1/nlzO/0Pv33ljqLBlYcyAj6ozgjZA3CI0NfegiUiWdSuCKopR4R44cYfTo0Xz33Xc0aNCAESNG0KhRI05lnOL1Ta+z8+JO7C3teSrgKYbXHo6rraupQy4WKoErilLiGI1GNm7ciLOzM02bNsXT0xODwUBSUhJGaSTcKpzlVss5tOUQLjYuTG40mcE1B9/2pmRRFJEqKVQCVxSlxLha+c9oNPLUU0/Rtm1bfv31V1xcXNj13y7WnltLvz/7cS75HF72XkxrOo0+1ftgY2Fzx3MWRRGpkkIlcEVRSoRx48axc+dODh48iIWFBRs3bsTf35+M3AxWnl7JD0d/ICYtBn9nf2a0nkFXn65mOXa7MJXtd68oiskcOHCAr7/+mm+++QZra2s6d+5M9erVMRgMWFhYUNm/MgtPLOSX47+QmJVII/dGTG82ndZere97pmRRFZEqKVQCVxSl2Jw5cwYXFxdcXFyIiYlhxYoVTJw4kYYNG/LYY48BEJsWy8/7f2bZqWWk56bTxrsNowNG06hiIxNHX/KoBK4oSrEIDw/H39+fjz76iBdffJHOnTsTHR2NjY3Wf30++Tw/HP2BVWdXYZAGuvl046mAp6jpUtPEkZdcKoErilJknnnmGRwcHPjoo4+oUqUKCxcupEuXLgDo9Xr0ej3HLh9j/uH5bLqwCUudJf39+zOizggql69s4uhLPpXAFUUpNGFhYWzbto3gisEA2Opt0ev1+a+PHDkSACkle2L2MP/I/DI7hrswqASuKEqBpKamYmdnhxCC+fPnM3PmTAb+MhDr8tbMnj37hn2N0sjWiK0sOLyAQ/H3HsOt3N2DrQukKIpynW3btlGpUiV2794NwKRJkwgLC8O6vPUN++UYc/jzzJ/0+7MfU7ZO4XLmZaY1ncaGxzYwpu4YlbwfkmqBK4py3zIyMpg1axZBQUH07NmThg0bMnz4cJydnQFwd89bweYIRKVGkZ6Tzu9nfldjuIuI+gQVRbmrtLQ0wsLCCAgIwNramh9//JHMzEx69uxJ+fLl+fbbb285JteYS3RaNF1XdCUpK+mhxnAr93bPBC6EWAD0AuKklAF521yApYAPcB4YJKVMLLowFUUxlQEDBnD27FlOnjyJTqfj8OHD2NnZ3bCPURo5fvk4O6J2EHIxhAOXDgBQz62eGsNdhISU8u47CNEGSAV+ui6BzwISpJQzhBBTAWcp5Sv3ulhQUJAMDQ0thLAVRSkqq1ev5s0332THjh2UK1eOHTt2IKWkVatWN7Se4zPi2XVxFzuidrDr4i4Ss+7chhtff3ypqsNd3IQQe6WUQTdvv2cLXEq5TQjhc9PmPkC7vOc/AsHAPRO4oiglT0ZGBitXrqRNmzZUrlwZR0dHHB0diY2NxdfXl1atWgHajciDsQcJuRhCSFQIxxOOA+Bi40JLr5a09GpJc4/mVLCtwKj1o0p1He6S4mH7wCtKKaPznscAFe+0oxDiaeBpgCpVqjzk5RRFKUxSStLS0rC3tyc2Npbhw4fz8ccf88ILL9CmTRu2bt0KaDciQ6K0hL07ZjdpOWnohZ76bvWZ1HASLb1a8ojLI+iEGtBmCgW+iSmllEKIO/bDSCnnAnNB60Ip6PUUpSQatX4UYPriSfcTh5SSJk2aULduXRYsWICPjw/79u2jfv36ZOZmEhobqiXtiyGEJYcBWk3t7r7daeXZiiYeTe5r2F9prsNdUjxsAo8VQnhIKaOFEB5AXGEGpShK4Vq3bh3BwcHMnDkTIQRDhgzB09MTKSVhyWEctjrMt39/y97YvWQZsrDWWxNUMYiBNQbS0qslvg6+Dzx6pDTX4S4pHjaBrwJGADPyvv5ZaBEpilJgUkpCQ0MJDAxEp9Oxb98+li5dyvTp05FWkvoD6rMjagddVnQhJi0GAD9HPwbWGEgrr1YEVgy86yIJSslwP8MIf0W7YekqhIgE3kRL3L8JIUYDF4BBRRmkopiDqNQoU4cAaHH88ccf9O/fn61bt9KmbRu6j+yOay9XJm6fyMFLBzFIA/aW9jTzaMbT9Z6mpWdLPO09Cy0GU3cllRX3Mwpl6B1e6ljIsSiKWYtOi773TkUkKSmJYcOGkVI7haS6STTu3ZgX5r/AauNq3vztTRIyEwCoXaE2TwU8RUuvltRzq4elztJkMZcpC3tqX0f9VainVTMxFaUQ5BhyAFh0bBFGaUQikVJixIiUEonEKI03vnbTcyNGkOQ/v/k44IZzXLp0ifSMdLy8vTBKI5caX8LgbACg++ruoAeXSy608GxxwxA/pfRQCVxRCmD2gdnMOTgn//uZe2be97ECgU7oEAiE0J5fHY6nEzp06BAi77Wrz68eIwSXL10mNzeXaMtoUrJSyPTLvOUag2oMYmLDiQV/o0qJpBK4ohTAY/6P8de5v4hKjcIgDewYsiM/KV9NtDckaXT52x/U7NmzmTp1KhERETg6OnL27FkqVKiAk5NT/j5qAk0JlnSh0E+pRt8rykOKS49j9MbRXM68jL+TPwCO1o6UtyqPvZU95SzLYWthi7XeGiu9FZY6S/Q6/X0n74sXLzJp0iROnToFQLNmzZgwYQK5ubkAVKtW7YbkrZRgUkJyRKGfViVwRXkIl9IvMXrDaC6lX+LbTt9ib2VfKBNX4uLiCAvTJs9cXSBh7969ADRq1IgZM2ZQocLd+7HVBJoS5vRmiD5QJKdWXSiK8oDiM+IZvXE0semxfNvpWxq4NwAKPnHFYDBQr149OnTowC+//IKHhwdxcXG3VP67FzWBpoSIOwG/DoHEsGvb3nLUvradCu1fLfAlVAJXlAdwOeMyYzaMISYthtkdZxe4TOqnn35KcHAwq1atQq/X8+2331KjRo381x80eSslQNplCP4QQhcgrewQXd6HE39B+E54K7lQL6USuKLcp4TMBMZsHENUahSzO80mqNK16p73O3ElPj6eZcuWMXbsWCwsLLC0tMTa2pqsrCysra3p27dvgWJUE2hMKDcb9nwP/8yErFQOWDbipTWX2Th1IuLkuiK5pOoDV5T7kJiZyJiNY4hMieTrjl/TuFLj+z42JyeHzExtiN+2bduYMGECISEhAPzvf/9j2bJlWFtb3+0USkkmJZz4i9yvGsOG1zB4NITxIRz3HY1n9XqkpqZq+zlWLvRLqwSuKPeQlJnE2I1jCb8Szlcdv6KpR9P7PjYmJgZvb28WLFgAQK9evTh06BBt27YtqnCVYpR1YS+5C3vCkmFk5eTSfXE6G93HgXsthg4dyo8//kj58nmVG52qFvr1VReKotxFclYyYzeNJSw5jK86fEUzj2b3PGbu3Lnk5OQwceJEKlasyLBhw6hbty4AVlZW+c9LtSKaOl5ipMaRvX46loeXkClssOj+EbaNRvD1Y5FUq1at2MJQCVxR7iA5K5mxG8dyLukcX3b4khZeLW67X25uLgcPHiQwMBDQSrdmZmYyceJEhBB89tlnxRm2UoS+/fpzqsdvopPVQaxyM9itD4K2L9O0aVd0cOfkXUS/yFQCV5TbuJJ9hWc2PcOZpDN80f4LWnq1vOO+b775Jh999BHR0dFUqFCBRYsWqdEjpURWVhYhISF0aN8ejv3BYzEzcbNIR/p2R3R+j6au1U0an+oDV5SbpGSn8MzGZziZeJLP2n1Ga+/WN7x++PBhmjZtyqFDhwAYOXIkv/32Gw4ODoAa+pevCKaOFwcptQJjAJ999hkvD+9M5pz2sGwkrp4+8OSfiKFLwMTJG1QCV5QbpGanMm7zOE4knuDTtp/StnJbDAYDGzZsYM+ePQB4eHhgMBhITNRWYff396dv375YWqrSrDcogqnjRe3EiRPUrVuXf/75B65cZFLlY4Q+bY91WiT0/gIxbgf4tTN1mPlUF4qi5EnLSWPc5nEciz/Gx+0+prGLNlRQSsmIESPo1KkTixYtwtXVldDQUBNHW8LlZmlfjUbQldx2otFo5I8//sDR0ZGOHTtSpUoVqnq6UTVsCYT8STljLrScgmj9Atg4mDrcW4irfyoUh6CgIKl+8JWSKC0njfGbx3Po0iE+bvsxS99byp49ezhw4AAAR44cwd/fX43XvpetH8I/M27d7tMa2r4MXoFgZdouJikl8fHxuLm5IaXkkUceoU6dOqxcvhwOL4O/34YrUVC7D3R6G1x8TRovgBBir5Qy6ObtqgWulHnpOemM+HMEJ1NPMrPVTDpV7URil0Rq1qxJbm4uFhYWBAQEmDpM82DvDkIHlnaQnQKBIyF8N5zfrj2EHirVhSrNoHITqNwMHIu3dsvEiRP566+/CAsLQ6fTsWHDBryJhvmdIGoveDSAx+ZB1duPOipJVAJXyqzTp09TzrEc0/ZN43TaaS7/eBmvhloyGThwoImjMzNSwpb3YPvHUKMbZCRDxC7o/YX2ekYiRIZC+L8QsRv2/QS7v9Vec6x8LZlXbgIVA0BfSKlpYU/2X0jmg6MezJ8/HwcHB/r160dAQAC5ublYpcfgE/oOHFkB5T2g7xyoN6REd/tcTyVwpUw6f/48NQNq0vGbjsRZxfFBqw/oOLgjNjZqJfYHZsiB1ZPhwGJo9CT0/Ax+6nPj1HFbZ/DvrD2uHhNzGCL+g4h/4cIuLYmC1nr3DtQSepWm4N0YbBwfKKRDhw7h6uqKJ5CVa2T79u2cOHGCJk2a0LlzZzq3aQbbZ8Cub7QD2r4CLSebvHvnQakErpQJUkqefvppXFxcmDlzJpW8K9H9u+5Eikjeb/U+Pf16mjpE85SVCstGwJnN0O5VLRFeXbDiblPH9Zbg1Uh7NBunbUuK0FrnEbu1lvr2j0EaAQHutbXW+dWuF2ffa9fJI6VECEFcXBwNGzZk6tSpvF8dmvo5ERkZiYWFBRgNcOAX2PIupMZC3UHQ6U1w9C6az6eIqQSulFrnzp0jJCSEJ554AiEEFhYW6HQ6MnMzmbRlEhEigvdbvU8vv16mDtU8pV6CXwZC9EGtqyRwZMHO51RZe9QdoH2flQpRoVorPfxfrYW+N6/aop37dQm9GePenk1mruSHH37A3d2dlStX0qpVK1j1JCI5XEveYdthw6tay9+7MQz5BbxvuS9oVtQoFKVUSU1Nxc7ODiEEr732Gh9//DGxsbE4OzsDkGXIYtKWSey6uIt3W75Ln+p9TByxmbp8FhY9BikxMPAHqNntxteLohaK0QCXTkD4vyQf3YQM/xcnozYWP0fqiZLu+LQedK0v3c5Vi+PCDnikF5xYAw7e0PltCHjslhZ8SXanUSgqgSulRnBwMD179mTr1q00adKE2NhYcnJy8PbW/jzONmQzeetkdkTt4J0W79DPv5+JIzZTUXth8SCte+PxZcXWik1ISMDR0RG9Xs/HH3/Ma6+9Rvix/6iUHX6t6+XiATDmaAdUqA7pCZCRoPWrt34Omj8LlrbFEm9hulMCN49brYpyGxkZGbz99tusXbsWgIYNG/LEE0/kL/RbsWLFG5L3c8HPsSNqB281f0sl74d1ehP80AusysHoTcWWvHfv3o2HhwebNm0CYOzYsURFRVGpegOo/Sh0fR/GbIZXI2DUevBtC5fPaMkbICdNGyWz4/NiibfYXJ33XxyPwMBAqZQeI9eNlCPXjSzWa6ampsrDhw9LKaU0GAzSx8dHvvbaa3c9Jjs3Wz67+VkZ8EOA/O3kb8URZum072cp33KWck4rKa/EFOmlsrKy5EsvvSQXLVp0w/enT59+sBMt6CHlmw5FEGHxAkLlbXKquompmJX+/ftz4cIFjh8/jk6n48iRI3ctHpVjyOGFf14gODKYaU2nMbCGGt/9wKSEbR/D1vfArz0M/hmsyxf6ZRITEzl58iTNmjXD0tKSLVu2oMsbj21lZcWsWbMK/ZrmTiVwpURbtWoVb7/9Njt27MDW1pZp06Yhrrv5dNfkbczhpW0vsTViK681fY3BjwwujpBLF6MB1r4IoQug3mB49GuwsCq008u8oX+gdYuEhIQQERGBhYUFu3btKpwCYUWwlFlJofrAlQKJSo0q1PNlZGSwePFioqK08zo4OODo6EhsbCwArVu3plWrVjck8dvJMebwyrZX+Dv8b6Y2mcrQR4YWapxlQk4G/PaklrxbPQf9vivU5L1mzRqqVavG5cuXAZg2bRp//fUXer0eoPCqOxbBUmYlhUrgSoFEp0UX+BxSyvyFX2NiYhg+fDjLli0DoF27dmzZsgUfH5/bHjtq/ShGrR91w7ZcYy5Tt01l04VNvNz4ZR6v9XiBYyxz0hO02ZQn/oLuH0Gntwo87C4tLY3vv/+ekydPAlClShVq166dX5a3QYMGNGrU6J6/nM3R4O92Mfi7XYV+XpXAlYeSkZvBxdSLAKw6u4q9sXuJSYvBKI0PdB4pJUFBQUyZMgUAX19f9u3bx6RJkx4qrlxjLq9tf42NFzbyYtCLPFH7iYc6T5mWFA4LumpD8gb+AE2ffuhTGQyG/BZ2RkYGEydOZMUKbcp8vXr1WLNmDdWrF+HCCKP+Kr3rclLAPnAhxHPAGEACh4FRUsrMwghMKbm+2vcVcw/Pzf/+9R2v5z+31Fniae+Jl73XtUd5L7ztvfGy98LJ2ol169axbds2ZsyYgRCCIUOGULnytX7Khg0bPlRcBqOB13e8zrrz63gu8DlG1Bnx8G+yrIo+BIsHQm4GPPE7+Nx5Kbn70bx5c7y9vVm5ciWurq4cP34cPz+/QgpWeegELoTwAiYBtaWUGUKI34AhwA+FFJtSAhmMBs5fOQ9AVYeqXLhygdV9VxOVGkVUahSRqZFEpWjPj10+RlJW0g3Hl7Moh2W6JZfFZexC7PBx9iFwcCBe9l6k5aRhZ/lwxYQMRgPTQqaxNmwtkxtN5qmApwr4Tsugc8GwZLi2cMFTG8C91gOfYvXq1fzxxx/Mnz8fgHHjxuWPy4e7LPpbBkQmphf6OQs6CsUCsBVC5ADlgIsFD0kpqaSUvPvvu/ndE8ERwVy4cgEfRx98HH1ue0xaTho//P4Dr854lefffh7biraEJ4cT7R7NmvNrSD9z4w+1k7XTbVvuXvZeeNp7YqW/9SZaVEoUb+x8gzXn1vC/hv9jTN0xRfDuS7nDy+H3ceDqD48vv+8a3UajkeDgYFq1aoWVlRXh4eGEhISQkJCAi4sLTz1Vtn+RpqamYm9vj9EoiUoq/M6Jh07gUsooIcTHQDiQAWyUUm68eT8hxNPA06DdtFDMk5SST/d+yorTK3i63tOMqDOC4IhgPOw8btk3MTGRYcOGMWzYMJ544glG9RqFbaotA9sNpHz58jecMykr6ZaWe1RqFCcTT7I1Yis5V6dFAwKBWzm3a0m9vBfxGfFEp0ez6uwqJjSYwNP1Hr6/tsza+RVsnAZVW2oFnmyd7nnI1eF/mzdvpmvXrqxcuZJ+/frxzDPPMGHChFJ5I/JBvf766/y49Hcmfb2C/RFJRXKNgnShOAN9AF8gCVgmhBgupVx0/X5SyrnAXNBqoTx8qIopzTs8jx+O/sDQR4bybINn87d72WsttdDQUGJjY+nZsydOTk5kZ2eTm5sLQLly5W7bEhNC4GzjjLONMwGut654Y5RG4tLj8pN6VEpeok+NYk/sHtacW4NE+5EaX3884+uPL4q3XnoZjVri/vcbqN1XGyZoefd66GlpafTo0YO+ffvy3HPP0aFDB5YuXUq3bloxKwuLsju15NSpU8yaNYuPPvqIuCwLzru3wmJAc2YHn8vfx2eqdkN1ckd/nutco8DXLMin3QkIk1JeAhBCrARaAIvuepRidpacWMKX+7+kl18vpjaZmt+6ysnIwdJWG6v72muvERkZSY8ePRBC8Pfffxf4ujqho5JdJSrZVSKwYuANr80+MJs5B+fkfz/n4BzmHJzD+PrjmdBgQoGvXerlZsHvz8DR36HpOOj64R1XodmxYwcXLlzg8ccfx87ODi8vr/x+bQsLCwYNGlSMgZcs4eHhWFpa4uHhQUpKCit3Hid63r8cvWzExlLHE818eKqVL1NXHGJ3WALnZxRy3fnbza+/nwfQFDiK1vctgB+B/93tGFULxfysPrtaBvwQIJ/9+1mZbcjO3/7VV19JBwcHeeXKFSmllKdPn5ZJSUnFHt/IdSNlwA8BxX5ds5aeKOXCnlqNkB2fS2k03rJLXFxc/vPBgwdLX19faTAYijHIki8pKUlaWVnJ5196Wf62J1x2+ewfWfWVNTLovU3yq79PyYTUrPx9B327U1Z9Zc1DX4s71EJ56HHgUsrdwHJgH9oQQh15XSVK6RAcEcy0HdNoUqkJz1d/nucnP8+ZM2cAaNasGePHjyc7OxuA6tWr4+j4YMteKSZw5SIs7KEtkND/e20ZsZv6q+fOnYunpyfR0dokrU8//ZQjR47k1yUpy6ZOncro0aMBkJa2jJj5C8H2HXhp+SGQ8NGAeux4pT3PdvDH2e7GG+5eToW/XF+BOqyklG8CbxZSLEoJsidmD88HP4+fvR9fdviSpLgkFixYQMuWLalevTpBQUEEBZWM1UxudyNVuY2447BoAGQmaXW8q7UH4OLFi0yfPp0JEyYQGBhI+/bteeONN7Cy0hKQp6enCYM2raioKNauXcvYsWMBrcsoy8qR6X8cYfneSDJybGjt78TY1n609ne9681bb+dyhR5f2b3joNzRkfgjPPv3s2RczMDqhBV2/e2w87IjLi7ursWjTOXqjVTlLi7sgl8Hg4UNjFpLeI4zKUePUqdOHezs7FizZg3t2rUjMDAQf39/pk+fbuqITSYjIwNLS0ssLCxYsWIFkydPpn379iRZuHCpZl9CjsWye084fRp4Maa1L49UcjBZrCqBK/k++eQTNu/fTEqvFJxtnJlQfQJNHm2S/3pJTN7KTW63lNmxVbBiDNKpMmL4CqRTVdr4+lKrVi3WrVuHo6MjkZGRhVc8yowdOnSINm3a8NNPP/Hoo48ybPgT2Ndqw4sbYjkQcRKncpZMbFedJ5tXxd3h/rtElj7TvEjiVQm8DIuPj2f58uWMHTsWvV5PhlUGMa1jcBSOfN/5eyo7lPwynAu7LTR1CCXbf9/D2pe4qPdi5FLY8GxVhBAsWLAAX1/f/N3KavI2GAx89NFHVKlShWHDhlGrVi2GDh2Km4c383eEsTAkjMjEDHwqlOPdPnV4LNCbclYlJ22WnEiUYpGTk4PBYMDGxoZ//vmH8ePHExAQQM1GNQnxDKFcdjm+72IeyVu5s5joaMLmP0Xz3J1QswfbDV3xrB5MWloa9vb2dOjQwdQhmkxmZiYnTpygQYMG6PV6li9fTsOGDRk2bBjx6blU6T2Jp/8KJyUzmiY+LkzvVZtOtSqi15W8yUlqUeMyJDo6mvr16/POO+8wbtw4srKyOH36NJX9KzNy/UiiUqOY12Ue9dzqmTpU5SFkZGSQu7Av5TMjibapjkdcMBHunaj8zFLQl8222tUSrtd3YYwYMYK1a9cSFRWFlZUV6enpnEvMYd72c6w5FI0EugdUYkxrPxpUdjJN4De506LGZfNftQz57rvvMBgMTJgwgUqVKjFs2DDq1KkDgLW1NX41/Ri7cSwXrlxgdqfZKnmbqaSkJHx9fXm9kysv1onD40okCQ3GU7nPhwWu423u4k7tp0GD8WzYsIGKFSsyefJknnjiCXQ6PX8fj+X77ef491wC9tYWjGjhw8gWPlR2KfwRI0VBJfBSJjc3l4MHDxIYqM1cXLt2Lbm5ufn1KT7//PP8fbMMWUzaOomjl4/ySbtPaObRzERRKw/jq6++IikpienTp+Pk6MiC53vSPXe99uKjX+HS6EnTBmgi2dnZrF69+lpDxd4JW1tbYmNjqVixIrXr1mflvii6frGds5fS8HC04bUejzCkSRUcbMzrXoDqQillXn31VT799FOio6NxcXEhLS3ttqNHco25vBD8AlsitvB+q/d5tNqjJohWeRCZmZns2rWL9u218dtPPfUUMTEx/DX7NcTy0ZAac+tBbadC+1eLOdLiJ6XM799PTEzEw8OD5557jnM+jxKZmE7I1I7Ep2bx864LLPr3ApfTsgnwcmBsaz961PXAUl+yJympLpRS6tChQ4wZM4YFCxYQEBDAU089RbNmzfKr/t0ueRulkTd3vsmWiC1MbTJVJe8S7GoDSwjBp59+yrRp0wgLC6Nq1ap898GLWG59D37oCeU9oc9s2L8YwkPgrWQTR168evTogV6vZ82aNTg7O7N7924CAgIYNu8/opIyeXXlIVbsiyI710inWu6Mae1HU18Xs6+aqBK4mTEYDGzatAlXV1eCgoLw8PDAaDTmL1vl7++Pv7//HY+XUjJrzyxWnV3FxAYT1XqRJdjx48cZOHAgs2fPpk2bNowYMYKgoCC8na1hzfNY7v0BLMtBh+nQbAJYlYMDv5g67GLxzz//8Ntvv/H1118jhKBfv343JGOXyv78sieSEzEpAKzcF8WAQG9Gt/Klmpu9qcIudCqBm4lrheGNjBgxgq5du/LTTz/h5ubGg3RLzTk4h8XHF/NE7Sd4pt4zRRix8qAMBgN//vknjo6OdOzYkapVq1KpUiWMRm2dUS83Z7xsDsLXo7QV44NGaV0k9m43nsixdA4BPXDgAI888gg2NjacPHmS33//nddffx1PT09GPDWa0POJvP/XMbaevMSZuNQbjs3KNfLL7nDc7K0LpYxrSaH6wB/A1dXPi3vyyJgxY9i3bx/79u0D4MiRI9SoUSO/VsX9+vnYz8zaM4u+1fvyTot3zP7Px9JASkl8fDxubm5IKalRowb169dn+fLl13YyGuDgr7DlfUi5CI/00laJd73NX1q3m4lZCoSEhNCqVSuWLl3KoEGDyMrKIjYlm+1nEwg+eYmQM/GkZxuw0uto6udCu5rutKvpxqsrDvHf+cTCL+NazFQfuBnZv38/33zzDbNnz8bKyorOnTtTq1YtDAYDer2egIBbFz+4l99P/86sPbPoXLUzbzZ/UyXvEmL8eG1429mzZ9HpdGzcuPHGlavO/A2b3oDYI+AVCAPmQ9UWpgu4mGRnZzNkyBBatWrF888/T7NmzZjz3Vwc/RvzwdrjBJ+M41Ss1sr2crKlfyMv2tVwp0X1CjfMlCztP+cqgZcQp06dws3NDWdnZy5evMjKlSuZNGkS9erVY/DgwQU696YLm3hr11u08GzBjNYzsNCpf3ZT2bdvHx9++CELFiygfPny9O/fn/r165Obm4uVldW16e0xR2DTdDi7BZyqwoAFUKf/vcd0m3HL+8iRI5w4cYIBAwZgZWWFTqcj1WjJr/+Fs/VEHCGRlUk7dwxLvaCJrwuDgirTrqYb1dzs75qoi6KMa0mh/ic/oKjUqEI/57lz56hZsyafffYZU6ZMoVu3bkRHR2NtbV3gc++8uJNXtr1CXde6fNbus9suCqwUrYMHD+Lm5oanpydZWVns2LGDEydO0LhxY7p06XLjzslRsPV97WakjSN0/QAajwGLgv8sFKfbzYC8nav3dgBmzZrF+g2b8Kjfhm1nLpPaahILYlJg5WE8HW3o09CLdjXcaFHdFXvr+09dRVHGtaRQCfwBRadFF/gcUkrGjh2Lq6srM2bMwM/PjwULFuSvK6jX69Hr9QW+zoG4A0zZOgVfR1++6fgN5SxL7w9ySSPzFv2NiYmhUaNGvPbaa7z77rs0a9aMiIiIW9eOzLwCIV/Arm9AGqDFs9D6BbB1Ns0bKAZLlixh1KhRhOw7yulUC2g1FjefYTw+fw8WOkFjHxde6/EI7Wq64+9+91Z2WaUS+H2SUpKekw5o46h14sEG/p87d46dO3cyfPhwhBDo9fobVjgZNWpUocZ7MuEkE/6egJutG991/g5Ha7VaTnEZPXo0UkoWLFhApUqVWLFiBa1btwa0PtkbkrchB/b+AMEzID0eAgZAx+ng7GOS2IvS5cuXmTFjBv0eG4DevToH8aHG/xbS/8fjAFRysOHRBp60reFOy+oVKF8IsyKLqoxrSaES+D3EpsUyLWQa/0b/m7+t/k/1AehStQtvtXiL8lblb3tsamoqdnZ2CCH4/vvv+fTTT+nVqxdOTk589913RRbzhSsXeGbTM9ha2PJ9l+9xtXUtsmspcOLECTZs2MDkyZMBbQWb60d39e3b99aDpISTa2HTm3D5NFRtBV3e0W5UlhKRiekkJiZy6dIlnCpVYcPxRH45b8uqP+PIkvFY6ASB1arS/hFtxEjNiuVVK/sBqWGEt5Gek86WiC2sOrOKf6P/RSJp4NaAxMxELqRcoLtvd0KiQriSfQULYUHDig1p49WGNt5t8HX0RQjBli1b6N27N//88w9BQUHExsaSm5uLl1fRrh4TkxbDiHUjyMjN4IduP+Dn5Fek1yurLl++jJOTE3q9nlmzZjFt2jTCw8OpVKnSvQ+O3Asbp0H4TnCtAZ3ehprdS1XRqUHf7uS/84lYnv4b4VWX7HLuALjZW9HhkYq0q+lGS39Xs6s9Yip3GkaoEngeozQSGhPKqrOr2HRhE+m56XjZe9G7Wm96+fWiqkNVRq0fRWhsKIdHHCbXmMuhS4fYFrmNbVHbOJ14GgAXvQtd/bsS6BTIis9X8NLzL911ZmRhSshMYOT6kcSlxzG/63zqVKhTLNcta3bt2kW7du1YvXo1Xbp0ISkpidyfBuJa3uruo0ASz8Pf78CRFWDnBu1ehUYjSkWpV6NRciImhc9/Xcv2kzFkV/DHILWVzh9xs6JXoC/tarhTy0O1sh+GGgd+B2HJYaw+u5o159YQnRaNnaUd3Xy70duvN40qNrqlr/vqAroWOgtq2NXA2saaKY9OIfJKJO1GtcOyoyW/n/6dXw2/YtPShi/Cv6CNUWudV7K7j9bZQ0rNTmX85vFcTL3It52+Vcm7EGVlZfHGG2/QoEEDhg4dSqNGjZg8eTJ+ftpfN05OTlD+LqN70hNg+yfw31wQemjzkrYavPXtu97MgZSSs5fS2HIkkhU7DhMrHUjKyAXKg0t5yGsXGoFjl7LpnGOktqfp1o4srcpkAk/OSmZ92HpWnV3FofhD6ISO5p7NmdJoCu2rtMfWwvaOx16/gG6/fv2Iiori6NGjeDt4c/Cng9jZ2ZGZm8memD1sj9rOtsht/BP5DwD+zv75XS313Oo99Hjsm2eEZuZm8uyWZzmVcIovOnxBUKWSsVq8OUtMTOTkyZM0a9YMKysrNm3ahF6vZ+jQoVhbWzNr1qx7nyQ3S0va2z6GzGRo+Di0fx0czG+VdyklEQkZ7DoXT8iZeHaevUx8ajYAuVdSaeZrxdBegTSvVgFPJ1sGf7eL3WEJZj8DsqQrMwk8x5DD9qjtrD67muDIYHKNufg7+/Ni0Iv08O2BWzm3e54jfFc4B389SMbBDGxtbZk+ffoNI0muVv6zsbChtXdrWnu35tUmrxKWHJbf1fLj0R+Zf2Q+DlYOtPRqSRvvNrTybIWTjdPDvS9jDi/88wL7Yvcxs81M2ni3eajzKNeG/oE2kuTff/8lPDwcCwsLdu/efe91I5MuXD2R1k3y9zvatmodofM7UOnBZ9A+qPsdf30/opMz2HX2MjvPXmbX2ctEJWUAIDOSqWafywv929Hcz4WE8FM0bNhQdY2YQKlO4FJKjl0+xqqzq1gXto7ErERcbFwY+shQHq32KDWda971hy4jI4MVK1bQoUMHPD09eb3t67y3/T1iY2Px8fGhTZt7J0shBH5Ofvg5+TEyYCQp2SnsvLiTbZHb2BG1g3Vh69AJHXVd69LGW2ud3yuuqwxGA69vf51tkduY3mw63X27P9Dno1yzZs0aJk+ezJ49e3BxcWH69OlIKfPH49/Xor/JEXBhp3aDMmovVKwLT/wO1cxj/cn41Cz+PXctYYfFpwFgqzdSSZfCO31a0NyvAj9+NYtWrVrSrYk25d/HtdFtz1eaZ0CWFGaRwB+0iFRMWgx/nfuL1WdXczb5LFY6K9pXac+j1R6luWdzLHV3/s94fWH4ixcv8sQTT/DFF18wadIkOnToUODFYMtblaerT1e6+nTFKI0cu3xMa51HbuOr/V/x1f6vcC/nTmuv1rTxbkMzj2a3nYATlRLF+7vfZ935dUxpNIVBNQcVKC6z9ZDFm9LS0li8eDHt2rWjRo0aeHt7U7t2bRITE3FxcaFhw4b3d6KMREg4B2mX8uLpDg5e0HcO1BsMuoJPyCoqyek5/BumJetdZy9zMlYrvWpvraeGs57He9aiebUKfP7WSxzYv5/hbw9Bp9Px3nvv3tf5S/MMyJLCLBL4/UjPSefv8L9ZdXYVu6N3I5E0dG/IG83foKtPVxys7n0DxWg0EhgYSOPGjZk7dy7VqlVj37591K9fv0hi1gkdAa4BBLgGMKHBBOIz4tkeuZ3tUdtZf349K06vwFJnSVDFoPzWeRUHrdUTnR7NslPLeCrgKUbXHV0k8ZU2BoOBpKQkKlSoQHp6OhMnTuTdd99l6tSpNGjQgNWrV996kJSQFq8l6esfiWHa14zEW4+5EgWJF0yWvCMT02+7PTUrlz3nE/K6ReI5evEKUoKNpY7GPi70aehJc78KBK/8ieefm8Ks06ep7unIF59/jq2treoiKYHMOoEbpZE9MXvyh/5l5GbgZe/FuPrj6OXXKz/Z3c1ff/1FSEgIH3zwATqdjqFDh1K1atX81++7JVYIXG1d6effj37+/cgx5LAvbl9+63zmnpnM3DMTHwcfkrKSABhYYyBTGk0ptvjMmZSSpk2b4uPjw/Lly3Fzc+P48eNUq1YNjEZtObKEc5AQdlOyDoPslGsnEjqt3raLn1ZcysVPe/wzC6L3l4iVcKKSMgHIzDGw90JifsI+GJmMwSix0utoWMWJyR39aVHNFSdjEmNGjWTkm2/SsEp1Kg0ehEelinh7ewNQrtyDt6RL+wzIksJsEvj1RaTOJZ/LH/oXkxaDvaU9PXx70Ltabxq5N7prS0FKSWhoKEFBQQghCA0NZcmSJUybNo1y5crx8ssvF8fbuSdLvSVNPZrS1KMpLzV+iYgrEby/+31CLobk77Ps1DKWnVrG+PrjmdBgggmjNbGrNw9vsmrVKlavXs3333+PkEZeGjsYT5ssCF0ACeeonhAGG/OSdG7GtQN1FloFQBc/qNL8WpJ28QOnKmBxmyGDu74pojf3YDKyDQAMmbuLfReSyDYY0esE9bwdGdfWjxbVXGlY2YndO7cjxEWa+NYgM7Mc2dnZpKdrLXcPDw+GDBliyreh3CezmMhzdQLNa01fY/XZ1RyOP4xO6Gjh2YJHqz1K+8rtsbG4vxsmy5YtY9CgQWzfvp1WrVqRkZGBtbX1DaNJSroR60awL24fh0ccNnUoprewJ1zYAW8lY8zJYvfG5TSuVgGL5HAObF1J0rl9tK7jhf5KJBiyrx2ntwYX3+uS83XPHbwffHLNwp7aL5LnjhTu+7tPb686ysKd52/Z3qe+J+/3r4u9tUV+5T8pJXXq1MHb25uNGzcWf7DKAzPriTxXW98f7P6AGs41HmjoX0JCAo8//jjDhw/n8ccfp0ePHsybNy+/X9vW9s5jvkuqBy2kVWqlXoIrF7XnXzSAxAs0xwh7tE31Le2goT/CxQ9qP3pjki7vCYX9S9up6r33KWRxKZl8s+UMv/wXjrWFDhc7K6KTM28Zf/3aa6+xaNEiwsLC0Ov1rFy58oauQsU8legEPvvAbOYcnHPDtlOJp0jLSbtr8t6zZw+XLl2iR48eODs7k5mZSU5ODqCN1R492vxv+l2dEVrmZF6BE2vgn5na1PSrEsPQAVH29XHt8w7WlWoh7N1LVX2R6yVn5PDdP2dZGHKebIORQUGVmdzRn8lL9hOdnMnJkyf56KOP+Pjjj3FycqJdu3bY2dmRnZ2Nra0tjzzyiKnfglIIzKoL5W5dBtcXhu/cuTPR0dEcOWKaP2eLmqnW5jSZnEw4vQEOL8d4ch06Y47W2q07gJc+mM1HrTNKxM3D4pCRbWDhzjC+DT7LlcxcHq3vyXOda+DraseFCxf439JDXDaW44tOjnTo0IFVq1bRrl07U4etFFCRdKEIIZyAeUAAWvWDp6SUuwpyzofx5ZdfMn36dKKiorC3t2fOnDm4ud27e0UpwQy5EBYMh1dgPL4KXXYq2Lmz/pIH8/9LZvmuAwidjo/Cd2t94KVcdq6RpXvC+XLLGS6lZNHhEXde7FIzv75IUlIS/v7+VOswmHr9JhAYGEhMTMxDjSBRzEdBu1C+ANZLKQcIIayAIvtpub7LIDIykpkzZ/Lcc8/h5+dH8+bNmThxItnZ2k2q6tWrF1UYJUKpbXlLCRH/weFlcOwPSLtElrDh1wOp9Jn2M84NHqVBbBw/Ozkhru+/dqxsspCLmsEo+fNAFJ9tPkVEQgZNfFyY/XgjGvu48PLLL5OQkMC8efNwcnLixx9/ZHGYtvSaEEIl7zLgoRO4EMIRaAOMBJBSZgPZdzumIFxyXDh//jw+Pj4ALFiwgNatW+Pn50fjxo1p3LhxUV1aKUpSQuxROLIcDq+A5HCyjTrSvNvg3Gss4VQlUv8nRp+2oLfA0/M2haBMcPOwqEkp2XQslk82nuJkbAq1PRz4qHclko7tpLGPNsba0tISCwuL/BouQ4cOZaiJ41aKV0Fa4L7AJWChEKI+sBeYLKVMu34nIcTTwNMAVarce2LN7RgNRv4Y9wdWfaz46aef8Pb2Ji4uLr94lGKGEsKuJe1Lx5FCj6jWgbRmzxEw4BXemzWIx2v1wh+YNq2uqaMtVjvPxvPRhpPsD0/Cp0I5vhhcn971vfjyyy947rnn6NSpI76+vrz//vumDlUxsYe+iSmECAL+BVpKKXcLIb4Arkgpp9/pmIIs6LBy5Upq1apFrVq1Hup4pQg8aB2SlBg4+jscXg5R2s+BrNKc15fsJ8y2Pr/+qY1JzsnJub/iUaXMocgkPtpwku2n4/FwtOGxmra8+1RPflm8iJ49e5KUlERCQkJ+HXKl7CiKm5iRQKSUcnfe98uBqQU4313179+/qE6tFKWMJDi+WuvXPr8dpJE4XSWWnnDn2W//RjhVoWPVv7Up7XnKWvI+E5fCxxtOsf5oDDYil56euXwyvht6jJwdOiR/SruTk5O2eISi5HnoBC6ljBFCRAghakopTwIdgWOFF5pitrLT4dR6rSb26Y1gyCZBOOPQfAoWDYewccMe9p7bTLplBeyAjh07mjpik4hISOfTjSf482A05awsmNLJnx9fHUFuo3rYWD4G6JkzZ849z6OUXQUdhfI/YHHeCJRzwKiCh6SYlat1SAw5cHar1q994i/ITsVoXxFd47H8l1GFpv3Hs3FjCzq71WT48JoMHz7ctHGb0KWULL7ZeobFuy+Qm5tL9tFN7Fo0g4pOdjzd4m81ekS5bwVK4FLKA4Bav6usklJbxGDNc3D0D8hIABsnsmr0ZsDbS2g5dARTu71OkNHI2bNdSnXf7f2shJOckcO0n7ey+lQqektrBjWuTHu3THSd++DmoJV0UMlbeRAleiq9UsJkpUL0QW21mahQiMorOnJwCSdFNQ7o2jL4xblYW1jR9IQ3bdtri1/odLpSnbzvJjs7m+W/r+Ks3pvlR6+QnJGDTcIZPh/diW4t65k6PMXMqQSu3J4hFy4dh8jQvIS9T/teGm/dNyedmhxm78VkpN4SAUybNq3YQy4prq7qZG1bjoXbz/De9gz09pdpV9ONF7vUJMBLLfSrFA6VwM3RQy4jdkdSan3ZVxN11F64eOBajWxbZ/AKRNbqBV6BCK9A3vvsW954YzrGNxzgrWRycnIYVsZGj9zs6ko43bp3J8OtDtTtRXhCOgE+lXizXwOaVVPlHZTCpRJ4WZSecC1RX32kx2uv6a3Boz4EjgTvIPBqBM6+HDt+nIEDB/Ldd21oVcOVkSNH0jRpFXASKHtD/673zz//sG/J1+Q0HcnmY7EkNH6GSzlW1LK2YOHIxrSr6aaWI1OKhErgpV1OBsQcvjFZJ5zLe1GAW02o0RW8ArVHxTqgt8RgMPDHH3/gfOk8HTr4UbVqVSpVqoTBoK344u3tjXcdN0i6/fqLpd2BAweoVasWegtLgg+e4XImOABjfgrFp4IT07vUpFddD3Q6lbiVoqMSuLm63TJiRiPEn7ouWYdqdUaMudrr5T21FnWjJ7Vk7dEAbK4t9iylJD4+Hjc3N3Q6HS+//DKBgYF06NABOzs7/v7771uvWQrrkNxN3JVMFm/YyTvf/ETDTv0JT7cg11gJh7Yj8/c5fzmds3GpKnkrRU4lcHOVHKGtRnM1WUeGav3WVxfgtSoPXg2hxf/AK68rxOE2haCu88wzz7B582bOnDmDTqdj8+bND12/pjTIzDFw9OIVQs9dYs6yDRicKpNi1LqKnJsNwNbenuH13GhYxYkFO8I4GJl8y0o4ilKUVAI3N/GnIV7rd+bTvLowOguoGAD1Bmkta+8gqOB/zyXD9u3bx4wZM1iwYAH29vY89thjNGzYEIPBgE6nw9fX9+6xFNZN1AK4n/HX90NKSXhCOgciktgfnkTIyYucS8jCILVWtJWtO1WtMpncsR4NqzhTx9MBG0t9/vG/7A4v0PUV5WGoBG4urlyEX4ZAzMFbX2s5BTresYbYDQ4cOEDFihXx8PAgIyOD7du3c/LkSQIDA+natWvhxlyCpWTmcCgymf3hiewPT2J/RBIJaVo1ZFtLPZapF8k8uZd5M14j0KcC7g73XjTby+n+FtZWlMKiEnhJl54AIZ/D7u/AaICm47QRJJH/3fcyYlfrRUdHR9OoUSOmTZvGO++8Q4sWLYiIiMDConT/GBiMktNxKRwIT8pL1omcjkvlaiHO6u72+Fqlcn7VQlZ+/ymt6voRFRmBvX0fKlSocN/X8XZWsyiV4lW6/+eas+w02P0t7PgCsq5A/SHQ7lVwrnptHPh9GDVqFDqdjvnz5+Ph4cGKFSto27YtoK3aUhqS99Xx11ddSsniQEQSByK01vXBiCTSsrXRM07lLGlY2Yn21Rw5tn0tY/p1om2LJpw9e5ZvYqtTw70cFnqdWrFdMQvm/7+3tDHkwL4f4Z9ZkBoLNbpr3SMV69y43x2WETt+/DibNm1i0qRJAHh5eaG7ri+8X79+RRa6KUgpiUrKZMGOMK3/OiKRiARtApKFTlDLw4HHAr1pUNkJPweBrSGVmjVrkpKSQuUnZ9Hc14m2LZpQrVo1Pv3004eOo6B98IryMMxiVfoywWiEoythy3uQGAZVmkOnt6BKs1v3vWkm5uXLl3FyckKv1zNjxgzefPNNwsPDqVixYvHFX4wMRkno+QT+OhzNr/+Fk2PQfoY9HG1oWMWJhpWdaVDFiQBPR2yttBuNUkpq165N5cqV2bhRWzgiPT1dFY9SzEKRrEqvFAIp4czf8Pdb2oQb9zow7Dfw7wL3MXsvJCSEDh06sGbNGjp37sy4ceMYM2YMrq6uRR97Mbo+aa87EsOllKxb9olOzmSQe3nGttEKZ82bN4+ff/6Z4OBghBB8+umneHhcWxxbJW/F3KkEbkoRe2DzW3BhhzYhpv/3EDDgrsP/srKymHasNoGBgQwBAgMDmTx5cv6Qv9K0YsvtkraNpY72Nd3pUdeDDo+489QPe9gdlsD5GT1JTk5myZIlpKV5YWdnh62tLY6OjiQnJ+Pk5ET37t1N/ZYUpVCpBG4Kcce1rpITa8DODXp8DI1GgIXVbXdPTEzk5MmTNGvWDCsrKzZt2oSVlRVDhgzBxsaGWbNmFfMbKDr3k7TtrK/92BqNhvznBw8eZNy4cbi5udG/f38ef/xxHn/8cVO8DUUpFiqBF6ekcAieAQd/BSt7aD8Nmo0Ha/tbdr069A+0kSR79uwhPDwcvV7Pnj17SlXxqAdN2lelpaWx9vUBePZ7CehJ69at2b9/P/Xr1y/+N6EoJqASeHFIi4ftn8CeeYCAZhOg1fNgd/sxxqtXr2bKlCmEhobi7OzMm2++CYBer92QKw3J22CU7DmfwNoHTNqLFi0iKiqKV155BTs7O6o264abndblJISgQYMGxfxOFMV0VAIvSlkpsGs27PwKctKgwTBtLLej9w27paamsnjxYjp06IC/vz9eXl7UqVOHxMREnJ2dadiwoYnewN096DT2uyXtnvU8aF/z1qRtMBjYt28fjRs3BiA4OJjDhw/z8ssvI4Sgbp9nCvdNKYoZUcMIi0JuFoQuhG0faXW2a/WGDtO10q15DAYDSUlJVKhQgdjYWLy8vPjggw94+eWXTRj4g7mfBP4wSft6n332Gc8//zxnz57Fz8+P9PR0bG1tVX1tpUxRwwgLw71WwjEa4PAy2Pq+1t/t01oby+194+cupaRx48ZUq1aNZcuWUbFiRU6cOEG1atWKNv5iUpCkfe7cOUaMGME777xD+/btGTx4MJ6ennh6apUU1dA/RblGJfDCICWcWg9/vwNxx7QVbXp/AX7t88dyr1q1ijVr1jB37lyEEEyYMOGGOhvVq1c3VfQFcnUa+8MmbSklwcHB6HQ62rZti4eHB9nZ2aSlpQHg6enJ4MGDi/U9KYq5UAm8oC7s1MZyR+wGl2owYCHU7osRCN66ldatW2NpaUlYWBg7duwgKSkJJycnxowZY+rIC+zqNPY3/jzywC3t1NRU7O210Tfjx4/Hx8eHtm3bYmtry+7du4vzbSiK2VIJ/EFdXQkn5gj8/Tac3gj2laDX59BwOFJngRCC9WvX0rNnT1atWkXv3r2ZOHEikyZNKjV9t9tPX+JAhFYN8bfQiPvu0waYOnUqS5Ys4dy5c+h0On7//Xd8fHyKIWpFKV3uXvFfuVVyBKwYC9+20lrdnd6CSftJfWQgbTt04ptvvgGgc+fOLF26lE6dOgFgYWFRKpK3wSgZ9O1Onpj/H9kGIwCZOUbWHYnhdGzqbZP3iRMnGDNmDMnJWsJv164dY8eOJStLmw5fq1YtbG1ti+9NKEopoVrg9ys1DhLOas+Pr4ZWU9gpmnAhIpmhrcphbwUeHh44OGhrTFpaWjJo0CATBlz44q5kMmnJfv47n8iAQG/Ox6cReiHxtsuInT9/HhsbGypVqsSVK1f47bffePLJJ2nTpg3dunWjW7duJngHilK6qGGE95KZrK2EE77zlpd+i/Pj9U2pnDp1qlS0ru8m5Ew8k5fsJzUrl3f7BDAwqDKDv9uVX4fkeomJibi7u/Piiy/y4YcfIqUkIyNDjSBRlIekhhE+qJxMbebk9k8gIwHq9OOXdbsYViWGuAmncXd3p2VUFAc/di7VydtglHz592m+3HKaam72/DK2GTUqls9//eoyYi+99BLJycnMnTsXZ2dnfvzxR1q2bAloMyRV8laUwqcS+M2MBq1WydYP4UokRzLcsew+l5rtBhN4ui1kx+SvYuPl5WXiYItWXEomU5YcYOfZy/Rv5MV7fQMoZ6W994iICM7tWIVfq0cBrY/fwsIiv4bLsGHDTBm6opQJKoFfJSWcXEv2+ulYJZ0Fz0akdJxJh85P8UUDIzWBmh72kFQZFxcXU0db5HaeiWfSkgOkZuUwa0A9BgZ6k5GRgUEv0Ov1LFu2jNBFM6j4iDbF/cMPPzRxxIpS9qgEDsjzOxCb34bI/4hI0vHb5Vq8+uYWygtBVFTUjcWjnEr3WokGo+TrLWf44u9T+LrasXhMU2pWKs+BAwdo27YtS5YsoXv37owaNYp+/frl1yFXFKX4FTiBCyH0QCgQJaXsVfCQilHMYc59PwI/w1lkeQ9E7y+4kOjJ0Oo18mdQlobKf/frUkoWzy09wI4z8fSp74Hr+c0c+CeWmoMHU7t2bYYOHZrfbeTs7Iyzs7OJI1aUsq3Ao1CEEM8DQYDDvRJ4SRiFEh0dzZ8/fMHY6nHoj64kS2fDn5er0eudVZRzvH1517Jg19nLTPp1H8kZ2bzbty6DgirTqFEjmjVrxpw5c0wdnqKUaUUyCkUI4Q30BN4Hni/IuYpSRkYGBoMBe9LJ+XMKT6Wvg+PW0GoK1i0nM8i27LYkjUbJN1vP8NnmU1hlXyFl7Sf0f2sXQghCQkLU6BFFKcEK2oXyOfAyUP5OOwghngaeBqhSpUoBL/fgEhISqP+IL0smtaSlOEDl3ExSag3Aodf74OBx7xOUYn+u38ILyw+T61qdPg08GVDVjazuM/IXjlDJW1FKtodO4EKIXkCclHKvEKLdnfaTUs4F5oLWhfKw13sQX375Jampqbz20vO4nPyVkxNsKZcbAnX6IdpPw8HVPCv/FVR2djarVq2iXr16JFpU4O3d2eS6+PBs0wq80LdBqR7PriilUUFa4C2BR4UQPQAbwEEIsUhKObxwQrt/mZmZ7Nq1i/bt2wOwN3QPDXUn4etfIDmCcn7toeMb4NWouEMzOSklqamplC9fnpSUFIYNe5weL37OIVmFqhXsWDS2BbU9HUwdpqIoD+GhE7iU8lXgVYC8FviLxZm8r958FUIwa9Ys3nrrLcIvXMA77RA/NDmJiD8Jdg2hz9fg1664wipxunbtip2dHb///jtY29Nz1l/sj8mid31PPuxfF/t7VA5UFKXkMo9qhAt7XlsNBzh69CgBAQH8+++/gLZq+3/LPsNrwyhYMgwhDTDwRxi7tcwl7+DgYCZNmpT/C65///706tWLPecT6PnlDo7G5/B+vwC+HNJAJW9FMXOFksCllMFFOQbcYJQsD41m69atAFStWpVKlSqRk5MDMUeovO15gg6/gUiO1FbCmbAb6vTNH8tdmkkp2b9/f35p1mPHjrF8+XLi4uIAePrpZ8iq1pYhc//FxlLHyvEteLxpVdXfrSilgFlUIzTO70G1qcE07fgoS5Ys0TYmhMHWD7Q1KG0coNVz0OQZsCpbIyf++ecf2rVrx4oVK+jfvz9ZWVno9XosLCxISMvm+d8OEHzyEj3reTCjf13K25SdiUmKUlqYdTVCnU6wb3JlHF5drNXl3vaRtuq7Tg8tJ0OrKVBGxnJnZWUxZMgQ2rdvz6RJk2jVqhVz586lXbt2AFhbWwMQej6B//26n8up2bzbN4DhTauoVreilDJmkcABnHNj4J8ZsOsbyM2ERk9A21fAwdPUoRW5Q4cOcfr0aR577LH8BH31Lye9Xs/YsWPz9zUaJXO3n+OjDSfxcrJl5YQWBHg5miRuRVGKlnkk8JQY7eu2WVC7L3SYDqV8LPf1i/7OmDGDLVu20LdvX/R6PVbdXmYnMPmmYxLTsnlh2UG2nIijR91KzHisHg6qy0RRSq2SPQpl64fwliMknLm27dgfWr93KbZ48WLc3d2Jjo4GtFKtR48ezZ8heTt7LyTS88vt7Dgdz9uP1uGbYY1U8laUUq5kJ/D2r8JbyVBFW9mFt5K1R/tXTRtXIYuPj+fFF19k7969ADRt2pTx48fnv161alUqVLh9oS0pJd9vO8fg73ah1wuWj2/OiBY+qr9bUcoA8+hCKYXJKCEhgfj4eGrUqIGVlRXz5s2jevXqBAYGUr16dT755JO7Hh+ZmE5SejYvLjvI5uNxdKtTiZkD6uFoq1rdilJWmEcCB3CsbOoIGPzdLgCWPtO8QOeRUtK8eXN8fX1Zv349Dg4OXLx48YGKR0UlZdLzyx3EpWTyVu/aqtWtKGWQ+SRwM18JZ968eSxevJgtW7YghODzzz/Hw+NaNcQ7JW8pJalZucReySIuJZO4K1mEJ6QD2h8my8e1oH5lp+J4C4qilDDmk8DNTFJSEkuXLuWJJ56gXLlyWFtb4+DgwJUrV3B0dKRbt26kZuVyJi6VuCuZxKVkEXvT16vb07MNt71GZGIGfb4JYXJHf57rXKOY36GiKKZmHgl81F+mjiBfZGL6HV8zGAxkZmZi0Fuxdsc+pnwwmwhdJSrXCCDWKZDKg+sw9tdjxKVkEnsli4ycWxOzraWeig7WuDvYEODlSEUHG9zLW+d/dXew4ZUVB9l7IYnzM3reJgpFUcoK80jgJYRRSqKSMtl5Jp7YvO6M2CtZxKZkEpOYTuix01iWr0Au2nC/ikM/YNFZ4OwJylnpqehgg1t5a+p6O9GxvLWWqMvb4O5wLUHbW1vcsy/bQleyBw8pilI8VAK/D5dSsvh513n2hycBMGze7vzXLDBQ2dUB9/LW+DgIanhb0rhujbyEfC05F3blPy8nm0I9n6Io5kcl8Ls4HZvCvO1hLNsbgfE2Nb+qpJ0gJ3Q5W3buzGs1F2x0yoPwdi5bRbsURbmVSuA3kVKy8+xlvt9+juCTl7C20DG0SRVGt/Kl//hXSarejeBxdfDx8SE9vT22ts+r4XuKopiESuB5cgxG1hy6yPfbwjgWfQVXeytGBbmy4ZvpdOn8On5udakc1JmkJKhUqRJgukV/CzoOXVGU0qHMJ/DkjBx+/S+cH0LOE3MlEw87wei61rw0qD3GnCzWfpRMWloaAOWc3fAiHRsb1f+sKIrpldkEHpGQzoKQMH7bE0FatoEW1SrwQb8AnunThl3Vq2HzeCewLMfu3btvOE71PSuKUlKUuQR+ICKJ77efY93haHRC4GWIIW3jQhbt/QedTseff/yOr6+vqcNUFEW5pzKRwA1Gyebjsczbfo495xOxMGYzooUfT7evwYGdwex37E52djY2NjbUrl37judRfc+KopQkpTqBZ2QbWL43gm+3nibqSjZeTrY8Wbccn08aRYfeK/BwtMWje3e6d+9u6lAVRVEeWKlM4HEpmfy86wI//3uBpPQcsqNP0dELFr40Gb1OMLXfeZONIFEURSksZpHA77eM66nYFOZtP8fyPeEYEXSpU4mxbfw4uSOJ1q1bYaHXpqCr5K0oSmlgFgn8bqSU/LHrOF9vPs7ZdBtsLHX4iRi8004z98lZADT2edzEUSqKohQ+s03gSSmpbDqZwIKQCxyPvoIhNYOnO1RlYtd6uNhZmTo8RVGUImc2CfxqGdfk9Bw+/mMXP+68gM7eBX93e97oVo2WXpbU9K9m4igVRVGKj1kkcKMhl6ikTB7/5A/2J9uQnm3AzdrA/9q7MaJLY1WLRFGUMsksEnhkcjYAO+P09GtUidGtfanj6WjiqBRFUUyrRCfwzzad4ou/T+d/L4WOlfujqOxSTiVwRVHKPCHlbQpdF5GgoCAZGhr6wMcN/m4Xu8MS1BJiiqKUSUKIvVLKoJu3q7W5FEVRzNRDJ3AhRGUhxFYhxDEhxFEhxOTCDOxmagkxRVGUGxWkDzwXeEFKuU8IUR7YK4TYJKU8Vkix3UCVcVUURbnRQ7fApZTRUsp9ec9TgOOAV2EFpiiKotxdoYxCEUL4AA2B3ffY9aGoMq6Koii3KvBNTCGEPbACmCKlvHKb158WQoQKIUIvXbpU0MspiqIoeQqUwIUQlmjJe7GUcuXt9pFSzpVSBkkpg9zc3ApyOUVRFOU6BRmFIoD5wHEp5aeFF5KiKIpyPwrSAm8JPAF0EEIcyHv0KKS4FEVRlHt46JuYUsodgKoipSiKYiJqJqaiKIqZUglcURTFTKkEriiKYqZUAlcURTFTxVpOVghxCbjwkIe7AvGFGI65U5/HNeqzuJH6PG5UGj6PqlLKWybSFGsCLwghROjt6uGWVerzuEZ9FjdSn8eNSvPnobpQFEVRzJRK4IqiKGbKnBL4XFMHUMKoz+Ma9VncSH0eNyq1n4fZ9IEriqIoNzKnFriiKIpyHZXAFUVRzJRZJHAhRDchxEkhxBkhxFRTx2Mqxb2QtLkQQuiFEPuFEGtMHYupCSGchBDLhRAnhBDHhRBldjkrIcRzef9PjgghfhVClLqV0Ut8AhdC6IFvgO5AbWCoEKK2aaMymasLSdcGmgETy/Bncb3JaGuyKvAFsF5K+QhQnzL6uQghvIBJQJCUMgDQA0NMG1XhK/EJHGgCnJFSnpNSZgNLgD4mjskk1ELStxJCeAM9gXmmjsXUhBCOQBu0hVaQUmZLKZNMGpRpWQC2QggLoBxw0cTxFDpzSOBeQMR130dSxpMWFP1C0mbkc+BlwGjiOEoCX+ASsDCvS2meEMLO1EGZgpQyCvgYCAeigWQp5UbTRlX4zCGBKze510LSZYUQohcQJ6Xca+pYSggLoBEwR0rZEEgDyuQ9IyGEM9pf6r6AJ2AnhBhu2qgKnzkk8Cig8nXfe+dtK5PuZyHpMqQl8KgQ4jxa11oHIcQi04ZkUpFApJTy6l9ly9ESelnUCQiTUl6SUuYAK4EWJo6p0JlDAt8D+AshfIUQVmg3IlaZOCaTUAtJ30hK+aqU0ltK6YP2c7FFSlnqWln3S0oZA0QIIWrmbeoIHDNhSKYUDjQTQpTL+3/TkVJ4Q/eh18QsLlLKXCHEs8AGtDvJC6SUR00clqlcXUj6sBDiQN6216SUa00XklLC/A9YnNfYOQeMMnE8JiGl3C2EWA7sQxu9tZ9SOKVeTaVXFEUxU+bQhaIoiqLchkrgiqIoZkolcEVRFDOlEriiKIqZUglcURTFTKkEriiKYqZUAlcURTFT/wfwT047NBeQdAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Get best-fit parameters\n", "\n", "# Instantiate model\n", "lm = LinMod(no)\n", "# Guess slope\n", "lm[\"slope\"] = 0.5\n", "# Thaw parameters\n", "lm.thaw([\"mean%d\" % (i+1) for i in range(no)])\n", "lm.thaw([\"slope\"])\n", "\n", "# Because ln inherited from MBOEv, it has a default chi-square objective function\n", "# Find best-fit parameters\n", "fr = sco.fmin(lm.chisqr, x0=lm.freeParamVals(), args=(x,d,yerr))\n", "# Set model parameters to best-fit values and display\n", "lm.setFreeParamVals(fr)\n", "lm.parameterSummary()\n", "\n", "# Plot best-fit model\n", "plt.title(\"Data and model\")\n", "model = lm.evaluate(x)\n", "for i in range(no):\n", " plt.errorbar(x, d[::,i], yerr=s, fmt='+-')\n", " plt.plot(x, model[::,i], 'k:')\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "EMCEE progress: 100% ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| ETA: 0:00:00\r" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Parameter mean1, mean = 2.38839, median = 2.39174, std = 0.132071, 95% HPD = 2.12382 - 2.62801\n", "Parameter mean2, mean = 3.77283, median = 3.77296, std = 0.133714, 95% HPD = 3.50435 - 4.02043\n", "Parameter mean3, mean = 5.94215, median = 5.94154, std = 0.128204, 95% HPD = 5.68971 - 6.19776\n", "Parameter slope, mean = 0.624081, median = 0.623544, std = 0.0189487, 95% HPD = 0.588183 - 0.661993\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD4CAYAAADhNOGaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA7PElEQVR4nO2daZgU1dWA3wMMIIKCLAZBHYzgChIdcUdcUBQjJG4YSTBRSKIYYxRB44po8DOJcYsGIy5RASUGMSEirhgRZXGQTWVEjKMEcUAUFBzwfD+qmqnp6aW6u7q7avq8z9NPVd26y6mamnvueo6oKoZhGEbp0qTYAhiGYRjFxRSBYRhGiWOKwDAMo8QxRWAYhlHimCIwDMMocZoVW4B4OnTooOXl5cUWw2jELFiw4DNV7Vjocu3bNvJJLt916BRBeXk58+fPL7YYRiNGRD4sRrn2bRv5JJfv2tfQkIgMEJF3RaRKRMYkuH+7iFS6v/dE5HPPvWEissL9DctWUMMwDCM/pFUEItIUuAc4BdgfOFdE9vfGUdXLVLW3qvYG7gKectPuAlwPHAb0Aa4XkXaBPoFhZE5PEVnsNlzmg/Otisgst8EyK/adisOdbiPobRE5OJaJNXKMxoKfHkEfoEpVV6rqN8BkYFCK+OcCk9zzk4FZqrpOVdcDs4ABuQhsGAFxnNt4qXCvxwAvqGp34AX3GpwGUHf3NwK4F6yRYzQu/MwRdAE+8lxX43z8DRCRPYFuwIsp0nZJkG4Ezj8Ze+yxhw+RDCNwBgH93POHgZeB0W74I+rYYpkrIm1FpLMbd5aqrgMQkVgjZxJGA2pra6murmbz5s3FFiXytGzZkq5du1JWVhZYnkFPFg8BpqrqtkwSqeoEYAJARUWFGT8yCsFzIqLAX9zvb1dVXe3e+x+wq3uerDFjjZwMqK6upk2bNpSXlyMixRYnsqgqNTU1VFdX061bt8Dy9TM09DGwu+e6qxuWiCHUbxFlktYwCsU7qnowzrDPxSLS13vTbf0H0iBR1QmqWqGqFR07FnzFamjYvHkz7du3NyWQIyJC+/btA+9Z+VEE84DuItJNRJrjVPbTEwi4L9AOeN0TPBM4SUTaueOnJ7lhhlFMagFU9VPgHzhj/GvcIR/c46du3GSNGWvkZIgpgWDIx3tMqwhUdSswEqcCXw48oapLRWSsiJzuiToEmKweu9bu+OlNOMpkHjA2NqZqlAaPPgobNxZbijo2bdoE7ncvIjviNE6W4DRuYit/hgFPu+fTgZ+4q4cOBza4Q0jWyDEaDb7mCFR1BjAjLuy6uOsbkqSdCEzMUj4jwsydCz/+MfzkJ/Dww8WWxmHNmjUA+4rIIpzv/3FVfVZE5gFPiMgFwIfA2W6SGcCpQBXwFfBTcBo5IhJr5IA1cjKifMy/As1v1fiBWaXr168fv//976moqEgfOQtuuOEGWrduzRVXXNHg3pFHHsmcOXOSpr3lllu4+uqr8yJXPKHbWWw0Hr780jl+8klx5fCy1157ASzzLBsFQFVrgBPi47s93IsT5WWNnPDwdvXn2897dW1bNDkyIZUSgMIqAjM6ZxhGo2fTpk0MHDiQgw46iAMPPJApU6Y0iDNp0iR69uzJgQceyOjRo7eHt27dmssuu4wDDjiAE044gbVr1wLw/vvvM2DAAA455BCOOeYY3nnnnYRlL1u2jH79+rHXXntx55131ssXYPXq1fTt25fevXtz4IEH8uqrrzJmzBi+/vprevfuzXnnnRfkq0iIKQIj75g3VKPYPPvss+y2224sWrSIJUuWMGBA/X2tn3zyCaNHj+bFF1+ksrKSefPmMW3aNMBRIhUVFSxdupRjjz2WG2+8EYARI0Zw1113sWDBAn7/+99z0UUXJSz7nXfeYebMmbz55pvceOON1NbW1rv/+OOPc/LJJ1NZWcmiRYvo3bs348ePZ4cddqCyspLHHnss+BcShw0NGYbR6OnZsyeXX345o0eP5rTTTuOYY46pd3/evHn069eP2BLf8847j9mzZzN48GCaNGnCOeecA8DQoUP54Q9/yMaNG5kzZw5nnXXW9jy2bNmSsOyBAwfSokULWrRoQadOnVizZg1du3bdfv/QQw/lZz/7GbW1tQwePJjevXsH/PTpsR6BkXds1aBRbHr06MHChQvp2bMn11xzDWPHjs06LxHh22+/pW3btlRWVm7/LV++PGH8Fi1abD9v2rQpW7durXe/b9++zJ49my5dunD++efzyCOPZC1btpgiKAALF8L3vgfOysXSI9ehoTBNNhvR5JNPPqFVq1YMHTqUUaNGsXDhwnr3+/TpwyuvvMJnn33Gtm3bmDRpEsceeywA3377LVOnTgWcYZyjjz6anXbaiW7duvHkk08Czo7fRYsWZSXbhx9+yK677srw4cO58MILt8tWVlbWYBgpX9jQUAEYNQoqK+H11+HEExvef+cdZ619nlawFY0gegKPPuosQX3tNTjyyNzzM8KBn+We3pVAubJ48WJGjRpFkyZNKCsr49577613v3PnzowfP57jjjsOVWXgwIEMGuTY1txxxx158803GTduHJ06ddo+0fzYY4/xy1/+knHjxlFbW8uQIUM46KCDMpbt5Zdf5rbbbqOsrIzWrVtv7xGMGDGCXr16cfDBB+d9nkA0ZDN5FRUV2ticd5x4IrzwAjz3HPTvD19/7YTdcYdT+ccqzJD9KXLm+eed5z3hBOc8G37+c5gwAe67zzkPAhFZEL98tBA0xm/bL8uXL2e//fbLKE22iiDo5aOtW7dmY5h2RZL4febyXdvQUAGIr+jfegvmzIFf/ap4MhWCxqbYDKOxYoqgAMQPkcSuv/668LIYhpEZYesN5ANTBAUk1kKOKYLKyqKJUhCCXC1kvYvoE7Zh6KiSj/doiqAAxA8NNbG37htbeto4aNmyJTU1NaYMciTmj6Bly5aB5murhgpAvCIotcrN/veNrl27Ul1dvd08gx/WrM9u7HT5lztklS4qxDyUBYkpggKQbI7AMEqFsrKyjD1qnZKlhdJsLZGWMjZIUUBKtUcQxPNar8Iw8ocpggJQ6nMEuVTipaY0DaMYlFiVlB/WrHEqrFmzEt8v1aGhUnlOw4g6vhSBiAwQkXdFpEpExiSJc7aILBORpSLyuCf8/9yw5SJypzRCx6Vvvukc77gjdbxSHRoyDCPcpJ0sFpGmwD1Af6AamCci01V1mSdOd+Aq4ChVXS8indzwI4GjgF5u1P8AxwIvB/kQxeTZZ+Hbb1PHKdVVQ0GO69scgWHkDz89gj5AlaquVNVvgMnAoLg4w4F7VHU9gKp+6oYr0BJoDrQAyoA1QQgeBmbMgFNOgVtvda6TVVZ+5wj+979g5csXc+bAtm2FKatUlKZhFBM/iqAL8JHnutoN89ID6CEir4nIXBEZAKCqrwMvAavd30xVbWC0W0RGiMh8EZmfyTrjYrN6tXP84APnmE4RJLuOkcB7XsG58kqYPTv5/f/8B446Cm6+OX1ehajE//Mfp5wPPoAStedmGDkT1D6CZkB3oB/QFZgtIj2BDsB+bhjALBE5RlVf9SZW1QnABHAsNAYkU97JdLgi3dBQGIY/brvN+SWT5eOPnePSpYWTKRFr1jimqZ95xrk+6yxYsMAZqjv55OLKZtSnPMv9AEbh8NMj+BjY3XPd1Q3zUg1MV9VaVf0AeA9HMfwAmKuqG1V1I/Bv4IjcxQ4H8UM9foeGwjDc8eqr6ec2wkT8uz35ZDjjDPjyS+d6yRLnuHJlYeUyjMaAH0UwD+guIt1EpDkwBJgeF2caTm8AEemAM1S0EvgvcKyINBORMpyJ4sT+3DLgnXec1l9YiHnQSuagKL7iTzZHEKvsliyB//43GNkSMWMG9O0Ld96ZvzKCIpnSjA3HFWquwjAaM2kVgapuBUYCM3Eq8SdUdamIjBWR091oM4EaEVmGMycwSlVrgKnA+8BiYBGwSFWfyVXo/fYLhzev+FZqbM4gGbEWeLoeQc+esOee2cuVjg8/dI7vvANffAErViSPu3w5HHooVFfnTx7DMIqLrzkCVZ0BzIgLu85zrsBv3J83zjYgIL9SxWXrVvjRj+Cqqxz/w5kQ8ztQU+McwzRHcMwx8Pbbycu+915nEnbaNBg5Mrsy8vlcYRpuM4yoYjuLfbJiBTz5JJx7buZpO3VyjjHLsZlWWnfd5aRJN6a/eTMccUTdBjc/vP126vuxMrOpzAtZOZtCMIzsMUXgk1z8CsenSTdHEM8VVzjHb75JXc6iRTB3Lhx2WN3kaSIyeZYgnzufNFa/z4ZRCEwR+CRRRZNtpZNpqzVZ/NWrYcuWxPF69sxcrkzKzjWPzZvhn/9Mn/b113Mv3zCM1Jgi8EkuLU6/q4aSkahMVdhtN2feIlk5YSGR/JddBt//PsyblzrtW28lzyNReFjfgWGEGVMEPkm0V2DuXH9pg6qsEqV76qnM882m/KDnCKqqnOPnnzvHDRtSb1K75BIYMiR53rW1mctnGIaDKQKfxCoe74TtQw9ll1eysf50lW26+/GV42WXJd/b4Jegxt63bYO993Ym3BPld+KJcOCBqfPwmuCIT2/7CQwje0wR+CTIoaEHH0wc7513/JedSI74cv70JzjhBF8iJiWooZaNG+H99+HCCxPnb3aCDKN4RFoRpFtFEyRBrkpplmT3xhdfpE6XaY8A6vYuZJNfqrix3dTZEtTqHlslZBi5E2lFUEhbOUEqgjZtskuXjSJIxf33p48T21F8++31w//zn+zMYKRTdn7J9e8hIk1F5C0R+ad73U1E3nCdL01xzakgIi3c6yr3frknj6vc8HdFxEzdGZEl0opg9OjsK4L3389MkQSpCLJdPprp0FAQrFrlHD/6qOG9a67xl0ehZM2QS6lv9+pW4HZV3RtYD1zghl8ArHfDb3fjISL749jdOgAYAPzZdeJkGJEj0orgzjtT28CprYULLqizrRNj2TJn4nL8+ORpP/ywzuQy1FVciSrE+DLjiR+eSbWhLJVzmqCGQTIxcZGqwt60Kbvyv/gi+bMUyHpoGTAQ+CuA6z71eBzbWAAPA4Pd80HuNe79E9z4g4DJqrrFtbhbhePEyTAiR6QVAaSuHF98ESZOhBEj6ofHhjReeSV52vJy6NoVHnnEqQw3bPAnT6KK7N//rn/dNEm7UTW1TEFsZsuF+OWdfmWIVyYPPJD83poc/df57GnsDlwJxPqE7YHPXQOLUN/50nbHTO79DW58Pw6bDCMSNGpFkIxMhiX++EfnGN+ryEaemPO1VBvKEk2ApxqW8j5LrsohXfp0yzuDLCtf+f7T2c68VVULYsg8qt73jNKipBWBn7SxOMla8Zlw+eX1y0/EtdemlyVXMlGEqeL6lUc1eUs/m7mCbdvqemiZLhh47bXXANqKyCoc/9vHA3e4YbH1XF7nS9sdM7n3dwZq8OewCVWdoKoVqlrRsWPHzIQ1jAJR0orATyWSyOH8X/8Kf/975uXGSNUjSNTziJmx9s4f5KNFnY3rzS++gAEDEs/VeCv555+vO6+qgpdeyk5GgJtuqjufMaP+vXSK5Xe/+x3A26pajjPZ+6KqnofjR+NMN9ow4Gn3fLp7jXv/Rdfs+nRgiLuqqBuOR74M7L4aRngIymdx0UhVeW3cmDj8xRed4wsv+C/HW8EMH56dPDEytT4a4/nnYd99My/7669hhx2c89Wr05ue9pKqYl2wAHbe2Tm/6Sb4y1/8yfPnP/svPxGpDNHloCBHA5NFZBzwFhCbyXgA+JuIVAHrcJQHrnOmJ4BlwFbgYtf/hmFEjsj3CFJxptu+e+45uP76uvBMzC4k6hHkSpB5+eHXv64733NPuPvuxPEyrUTTeS2LzXe8/HKwPZgghqucuPqyqp7mnq9U1T6qureqnqWqW9zwze713u79lZ70N6vqd1V1H1X9d7JyDCPs+KqSRGSAu2mmSkTGJIlztogsE5GlIvK4J3wPEXlORJa798sDkh3w/48/dqxX1szzzyTNf/+b2h9Asp3F6Ug3MZzsXXiHmzI1zpbJcy9fXmdEDvyvtDIMo7ikrZLcTTL3AP1xlsjNE5HpqrrME6c7cBVwlKquF5FOniweAW5W1Vki0pq6JXuRIJuWbMzfcLK02bqqzPeS0VzyF4H994d99kluMylZuiAxkxOGkTl+egR9gCq36/wNzkqLQXFxhgP3qOp6AFX9FLbvvmymqrPc8I2q+lVg0lP/H//jj/0t85w5M/Ny/FZYiSqidu385bV8eeLwQlBZGUw+776bODyTFUap+Oyz1OaqC2l2xDAaC34UgZ+NMz2AHiLymojMFZEBnvDPReQp167LbYm24eey1tpbcXTt6mwES0cmJouXuf2eXOzdH3RQ3fmdd8L55yeOl0lLOuiW7/e+F1pTEPWoqYH165Pftx6BYWROUNOWzXCWz/UDzgXuF5G2bvgxwBXAocBewPnxiXNZa53JP34mFW225aSLd+ml2cuQbdm5Dvlki7fcZPnEh/sxZBc25WQYUcePIvCzcaYamK6qta7dlfdwFEM1UOkOK20FpgEH5yy1h0wqudiGrkzYay/n2Lp15mmDJpGJiSBbwH7y8u589lbI6Spnv0rq3HMdb2Sp8glq1ZBhGA5+FME8oLtrprc5zjrq6XFxpuH0BhCRDjhDQivdtG1FJNbMPx5n3XVReOklZzNYJsQqnSi3QoOUPVklHWQFnGx5Kzjmr1OxZUvdBjzDMPyRVhG4LfmRwEwcs71PuJtpxorI6W60mUCNiCzD2aE5SlVr3A02VwAviMhiQAAfVvD9k0kF9PXXqTeD5Zp/WFF1fldemT5eOrJdErp5c+LwCRMyy2fWrNSK7coroVWrzPI0jFLH14p2VZ0BzIgLu85zrsBv3F982llAr9zETCVbvnKun382q4Yefhh6985+JcvDD8Mee6QvJ1VYjE8+gdtuy1wGv5vf0r2fsrLE4U89lTzNl18mLico5zaGYTg0ahMTQea/bl3maZOtDipU+hhBTXTHI5K9Kep04ZDYFEaUh+gMI6xE3sREVVVim0J+hzCefDL1/ZiHrpi5inSEdSjJTwX6cQPbmanJZANcNpvoMlnmaxhG9kROEcRPJA4aBKec0jDewT7XJp19Nrz3Xt11+/aOV7N4tmzxl1+8NcwgSeeYJtVOZj+K4NlnE6cNglSyJduElkwRZGuiwzCMxERKEaxYkXjVSqKVJJm4PPS6XFy3zvFqli1jElpiCoZLL4XHHstf/onmMlIpgnTLR709tVQt/4ULM5Nnl12S52UYRuZEShFcdVV+8o2SWYKhQ5PfS1bZPvecs9omHZkqAi+JJna9LkKzGRpKJo/NExhGsERGEXzzTWpnMIcdln3eUVIEMTKdi/Ca4U5GLu/h0UezS5eqUk80NGRKwDCCJzKKIN2Qy5s5+IYK6wRvpqQyHOenAi2GQtyyJbFsV1xRf8guhkjuDu4Nw6hPZKbdVqzIX96NpUfw85/nlmem7yGfrfM//AFOOy1/+RuGUUdkegT5rHS2bWv8Zgk++CB9nFdfbRg2e3bwsvil0J7cDKNUsX814IwzHLMEkyYVW5Likot11lyoqkocbnMEhlEYTBFQN+Z8773FlSMT8jGvUawhsmuvTRyeSJ6//S2/shhGKWKKwEOioZFSIpO9F5C6dZ6hf6GERHHuxjCiSGQUgQ0JFJbnnksfZ+vW5PdGjap/nU0PprGs5jKMsBMZRVAIEpmqCCv5riSTDdf4JQg7QWZryDAKgykCD1GxYRNvEyjRrt6w4WfVUjw2NGQYhSEyiqAQQ0NNm+a/jCDwGskDf35+MyWXDXpQ36UlwB13ZJ7H7runj2MYRu5ERhEUgmnTii2Bf7xDQ9n4Ssg38YogG2L+og3DyC++FIGIDBCRd0WkSkQSGnsQkbNFZJmILBWRx+Pu7SQi1SKSwhutkS2jRxdbgoYEMaxjk8WGURjSjoqLSFPgHqA/UA3ME5HpqrrME6c7cBVwlKquF5FOcdncBBRxj2p4+eyzzNPU1NTfhFVbG5w8XhI5/PFLEIpg8eLc8zAMIz1+egR9gCpVXamq3wCTgUFxcYYD96jqegBV/TR2Q0QOAXYFfCxILD169848zdixUFERuCgNaNMm+7RBKIJU/owNwwgOP4qgC/CR57raDfPSA+ghIq+JyFwRGQAgIk2APwBXpCpAREaIyHwRmb82yU6kxrqPIFP3kIkI4xCKrfgxjOgQ1GRxM6A70A84F7hfRNoCFwEzVLU6VWJVnaCqFapa0bFjx4BEKh3yoQiOPz639KYIDCM6+Fk5/zHgXcjX1Q3zUg28oaq1wAci8h6OYjgCOEZELgJaA81FZKOq5tGhY3qH9I2NZK4ecyHXHpgpAsOIDn56BPOA7iLSTUSaA0OA6XFxpuH0BhCRDjhDRStV9TxV3UNVy3GGhx7JVglkUjGdfXY2JRhBYruCDSM6pFUEqroVGAnMBJYDT6jqUhEZKyKnu9FmAjUisgx4CRilqjX5EtpoSNAt8FyHm6xHYBjRwdccgarOUNUeqvpdVb3ZDbtOVae756qqv1HV/VW1p6pOTpDHQ6o6MltBwzghGibCVvGGTZ4YmzdvBthPRBa5e15uBHB7vG+4e2WmuL1fRKSFe13l3i+P5SUiV7nh74rIyUV5IMMIgMjsLDZFkJqwvZ9UlkmLSYsWLQDeVdWDgN7AABE5HLgVuF1V9wbWAxe4SS4A1rvht7vxEJH9cYZJDwAGAH9299wYRuSIjCKYM6fYEoSboBVBrpPF+drkliviPFisv1Lm/hQ4Hpjqhj8MDHbPB7nXuPdPECeTQcBkVd2iqh8AVTh7bgwjckRGEZx3XrElCDdhG4oJ+2SxiFQCnwKzgPeBz935MKi/V2b7Phr3/gagPf721/jaI2MYxSYyiuCII4otQbgJ29BQ2BWBqvbGWQrdB9g3j+XYHhkj9ERGEYStxRs2wrZqKOyKAEBVP8dZ5XYE0FZEYvtqvHtltu+jce/vDNTgb3+NYUSCiLhiMUWQDusR+MMdnmkKICI74BhTvBVHIZyJY0trGPC0m2S6e/26e/9FVVURmQ48LiJ/BHbD2UCZoxcHIwjKx/wrq3Srxg8MWJLoEBlFsGlTsSUIN6+8Emx+uU4WO6s0w8fq1asB9hGRt3F6xE+o6j/dPTCTRWQc8BbwgJvkAeBvIlIFrMNZKYS7l+YJYBmwFbhYVUOq/gwjNZFRBFOmFFuCcBOm3dQ1NbBqVbGlSEyvXr0AlqlqPfutqrqSBKt+VHUzcFaivNw9NTfnQUzDKCiRmSMI61BDWPjqq2DzW7YsfZxkdOgQnByGYeSfyCgCZx+QUSicERTDMEqByCiCk20Dv2EYRl6IjCKwHoFhGEZ+iIwisJ3FhmEY+SEyimDnnYstgWEYRuMkMorAMAzDyA+mCAzDMEocUwSGYRglji9FICIDXC9MVSKS0OewiJwtIstcr0+Pu2G9ReR1N+xtETknSOENwzCM3ElrYsL1unQPjnGuamCeiExX1WWeON2Bq4CjVHW9iHRyb30F/ERVV4jIbsACEZnpWn00DMMwQoCfHkEfoEpVV6rqNzjWGQfFxRkO3KOq6wFU9VP3+J6qrnDPP8FxBGJG2Q3DMEKEH0XgxxNTD6CHiLwmInNFZEB8JiLSB2iO4w0q/p55cTIMwygSQU0WN8Oxx94POBe4X0Taxm6KSGfgb8BPVbWBZwG/XpyGDQtI2hJhzz2LLYFhGFHAjyLw44mpGpiuqrWuI+/3cBQDIrIT8C/gt6o6NxdhTzghl9Slx9VXF1sCwzCigB9FMA/oLiLdRKQ5jmOO6XFxpuH0BhCRDjhDRSvd+P8AHlHVqbkK26NHrjkYhmEY8aRVBKq6FRgJzASW43h0WioiY0XkdDfaTKDG9fL0EjBKVWuAs4G+wPkiUun+emcr7GGHZZvSMAzDSIYvD2WqOgOYERd2nedcgd+4P2+cR4FHcxfTyJQnn4R164otRfHYtAl23LHYUhhGNLCdxY2UQYPC59C+kDxqzQ/D8I0pgkZKWRl822B9VulQykrQMDLFFEGJ0qeBm/bGhSkCw/CPKYJGTKrKUKRwcviloqLYEhhGaWKKoBFz/PHJ7331lf98Dj88d1n80L59cHlZj8Aw/GOKoBGz777J7y1e7D+fE0/MXRY/BNlLMUVgGP4xRWCk5dpriy2BYRj5JNKKoJmvXRBGrjRrBq+8kv9yMu0RPP988nvWIzAM/0RaEcyZU2wJokvbtnXn++yTPn7fvnkTJWtS2Z4yI7aG4Z9IK4IWLYotQXTxLh895JDiyeElkx5Bu3ap748dm5sshlFKRFoR9OxZbAkaB+kq4DAuNf3xj4stgWE0HiKtCMJYQUWFML67TGQKo/yGEVUirQiixtKlhS/zgAMSh3srUqtUDaO0ibwi+O1viy2Bf/bfH774orBlJqvkM5n8LZSisB6BYRSHyCuCvffOLt2yZcHKEVYSVZh//3t9Jz9hqVRNERhGcYisIog5qdlhh+zSt2kTnCyZUOgKzE95UaxUs5X5o48+AughIstEZKmIXOrkJ7uIyCwRWeEe27nhIiJ3ikiViLwtIgfXySDD3PgrRMQ8ahuRJbKKoKzMOZ5xBlxzTebp81H5hWWYKp0nt/hnL6YiGDo0u3TZytzM2YVYrar7A4cDF4vI/sAY4AVV7Q684F4DnILjf7s7MAK41ylfdgGuBw4D+gDXx5SHYUQNX4pARAaIyLtuq2hMkjhne1pZj3vC89pqatYMbrop6Fz9M3Jk3Xl5edHEqEcTz181UYW5bVvhZElHtgbtslUEnTt3BvgKQFW/xHG/2gUYBDzsRnsYGOyeD8Lxua2qOhdoKyKdgZOBWaq6TlXXA7OAAdlJZRjFJa0iEJGmwD04LaP9gXPdFpQ3TnfgKuAoVT0A+LUbHniraebMXFIHzw9/WHcelg1u6VYExZtmKGaP4IwziieHiJQD3wPeAHZV1dXurf8Bu7rnXYCPPMmq3bBk4fFljBCR+SIyf61tdzZCih9rPX2AKlVdCSAik3FaSd7p1uHAPW7LCFX91A3f3mpy08ZaTZOyFThW2eZaaQRV6TRtWnceRo9giZ5z48a68x/8oLiKIF3ZS5bAgQdmni59udIa+Dvwa1X9QjwZqqqKSCDWilR1AjABoKKiItIWkMrH/KvYIhh5ws/QkJ+WTw+cCbjXRGSuiAzIIG1WraYwTnCGRRE0yXDmJyyKIJEcBxwArVqlTpdNsThK4DFVfcoNW+MO+eAeY42Zj4HdPWm7umHJwg0jcgQ1WdwMZzKtH3AucL+ItPWbWFUnqGqFqlZ07NgxIJEKg7dCGjy4aGLUI13l2rNnOBRpdbU/Oby9rlxRxyzpnsByVf2j59Z0IDaHNQx42hP+E3f10OHABncIaSZwkoi0c4c7T3LDDCNy+Bka8tPyqQbeUNVa4AMReQ9HMXyMoxy8aV/OVlgIzrxwPirCdIbQCkVsRRUkfs7ddqt/XSyl0KUL1NSklyPTHk4qXnvtNYD2wPEiUukGXw2MB54QkQuAD4Gz3XszgFOBKpxJ5p8CqOo6EbkJmOfGGxsbAjWMqOFHEcwDuotIN5yKfQjwo7g403B6Ag+KSAecoaKVwPvALZ4J4pNwJpWzJqYIwtCiDSu7e9R2ovf0/e/Diy8WTp5E/PSnztHP3zGRIsj273/00UcDLFDVRB6SGxi2VqcLcXGivFR1IjAxO0kMIzykbWup6lZgJE63dznwhKouFZGxInK6G20mUCMiy4CXgFGqWuO2kGKtpnkE2GoKiyIIQo5ULiXzgdcXQbGY6FaffmweBTk0ZBhGQ3z5+FLVGThdZG/YdZ5zBX7j/uLTBtpq6t3bOY5JuJshNeecA1OmOOfeSqd9+/pDFIUmvgKsqID587PPz1kqnzhvP+WHhS7usoJEisA8kBlGcERuZ3G7dk4lcPLJ9cNPPz1xfC8nnQStWzcMr6zMTab1651fUPz5z7mlP/jguvNYJZ9q93X79rmVFxTxCunMM52jd2jI3JMaRvBEThEko3lzf/EStSS7dvVfTrzCEXGGWvwOtyQa7/5R3IyL32eJ57rr4NFH6ypQLwMHJk93/fVw993ZlZkJo0dnFj+mGPbaqy4sfhf5e+/l1nsyDKMRKQI/BGGD/+mnc6s0W7WCP/4R/vIX5/rtt4OzUdS8OZx3XubPucMOcHHC6dBg6dUr9f1kso4a5RxPPbXhve7dw7NayzCiSqNRBH4r9jCMh192GYwY4fRO4tf0X3xxsOPfsbzDMKae6N2n+nvE7sXmCETq5nhiR8MwcqfRKAI/fOc7wecZtGJp2TL7SjtVRRsGRZAtJ53k+Cj+859h+XIn7MMP6+6HQbkbRpQpGUUwbhycckrddS6VRzrrnpnMOcQjEqwiiBGf56ZNztFboQbJV181DFuXYOFwqmGs2HXz5vDII7DHHonfjSkCw8iNRqMI0lUGAwI0EJyqrKVL4ZVXcss7Hz2CeGIbyhYu9J9/RaItWElI5DDojjsahnl3QWeLKQLDyI1GowjSUajKYv/9Yeed666rqwtTLmQ2NBSrlN94w3/+3vX8J56YmWwAK1Y0DPMalMvWYY4pAsPIjUajCDKtDHKpPNKtyvGGdWlgazU93g1hyejQAe66K3m58WGqsM8+deE77eSE9enjXy6vMtlzT//psiXR8+yyi3O0/QSGERyNRhFcdJFzjC3LLCa5TMzuuCP4McBaVlbfOxqknyN49VWYPTt72eLzC5r4vQ6Jnie21PbHPw6+fMMoVRqNIujb16mcRowotiTZceWVztGvl7Nshk06doRjjslMLi9efwuJFMH3vpd93gDDfDgybdnSOdpwkGEER6NRBJmSz1VD+ayk5s5Nfi9RuYWsOPv1Cza/eHPZkHjOw5SCYeRGySiCWGUR9UojfmmqqrNBDRI/28MPw9VXw1FH5V52uuGgRJ7EcuGSSxqGxZRwlPdFGEbYKBlFEGTFka41migsE/MK8bJ26uRPnkTldu4MN98cjHOXIN7hfvulvu81gJfI6mjsGcPiFtQwGgMlowjiKWTPYM4ceOklf3FVG1a4r76aPg3k55m6das7D6LyffPN1Etq0y23TdQjsN6BYeSGLcIrAEcckT5Osh7BZZdBjx6p0+ZTEXjz/Oab3PNr3bqhKfDbboO993bOY3MayTjwQOd43HG5y2IYhoMpAh8sXFjfxn82Q0N+ie8RxPI67TR44YXUafOtCM44A5YtC76MK67wH/fgg2HNmvpLbKM+72MYxcbX0JCIDBCRd0WkSkQa+AYTkfNFZK2IVLq/Cz33/k9ElorIchG5U6Q4/7aZ7Fp13No6piJUUy+L7Nkzd9lSyRMLf+YZx35P584wdChMm1YXJ59DI165vBPOicpMJ0cQ5iTAmTOxyt8wgiNtj0BEmgL3AP2BamCeiExX1fi24RRVHRmX9kjgKCA2Vfof4Fjg5RzlDoSWLbPbIRur8EaMSD+UkU3eqXocTZrA3/6WWJ589whibkINw2hc+OkR9AGqVHWlqn4DTAYG+cxfgZZAc6AFUAasyUbQTPDuUE21WmfTpsRDHWGyceOnjEIpgl13bVhmKuIVVhArl/xyyCGFK8swoo6ff80uwEee62o3LJ4zRORtEZkqIrsDqOrrwEvAavc3U1WXxycUkREiMl9E5q9duzbjh4jn5z93jqedlnrIpUkT/5XT4sV1RtOOPNI5Dh6cPO9cyHSzVKEUQaIyUzF0aP1NZoUczunQoXBlGUbUCaqN9gxQrqq9gFnAwwAisjewH9AVR3kcLyINjByo6gRVrVDVio5+DO1kQboNZRdf7PgjXrTIK1fd+YEH1q1s6dXLGa/Pl2nrKCiCROGJlMPzz9ed9+8fjEx+ZLE5BMPwjx9F8DGwu+e6qxu2HVWtUdUt7uVfgVjH/AfAXFXdqKobgX8DPhZT5kb//k5r/U9/angvWUv27rsdf8TpNn7FSGRvPwjSzREkS+M3birmzWv4znLN07sp7NBDc8srEwo5DGUYUcfPv8s8oLuIdBOR5sAQYLo3goh4DSefDsSGf/4LHCsizUSkDGeiuMHQUNC0bAn/+Ad897vJ4/ip4LKtBLNJl0nLO542bZxjrsqpogIuvdRf+dlsLstXKz2Wr9d8t/UIDMM/aVcNqepWERkJzASaAhNVdamIjAXmq+p04FcicjqwFVgHnO8mnwocDyzGmTh+VlWfCf4x0pNNq7kYO1az6RHccAO0bZsf08xRqFBtaMgwcsPXhjJVnQHMiAu7znN+FXBVgnTbgJ/nKGPBKUYlkqxMP0McrVo5huXyQZDvYqedgssrEWaR1DCyo2RGUrNxg1jIHsGOOzrHVq3qD/EUu0KLL/+pp+DBBzOTK+al7eSTg5PLSyLT1IkM1hmGkZiSNTGRqpLPtfLNJv1ll8G2bY7p5bIyuP56uPHG8CmCH/zAOc6blzj+I49AeXn9sJh56nzPEXgxRWAY/ilZRZCKYswNtGhR54bRK0PYFEEqVBPPUxTqfXrLMZ/GhuGfkhkaiqfYFWw6YtY1g/b6lSmZbCirrMwur1xp29Y5xhz0gCkCw8iERv/vkqw1GrahoXj69YPNm/37MM4XmTxL0B7K/NKyZd3f8yp3yYLtIzAM/5TMv0vYewCJKLYSgOTvbeedG4bFdl7HEybHMT/72c8ADhKRJbEwEdlFRGaJyAr32M4NF9dibpVrPuVgT5phbvwVIjKs8E9iGMHR6HsEuZBtBRZFpZMp3s1b6SjGfEeyss4//3wefPDBFXHBY4AXVHW8a2Z9DDAaOAXo7v4OA+4FDhORXYDrgQqc/TELXIu86/PxLEZhKB/zr6zSrRo/MH2kkFMyPYIYhaiMwtQCzpVM5gjSvdswKMi+ffuCs/HRyyBc+1jucbAn/BF1mAu0dXfRnwzMUtV1buU/CwjQ8pRhFJaSUQTZVM65VlzFGjPPlV12qTsPovK+6CLnmCd7ggnJUO5dVXW1e/4/IGZwO5nlXb8WeQ0jEpSMIsiGUh0a+uijOquh6WwgXXIJ3HJL6vwuv9x5l/neWewl27+BqirOcE9AcgRrYt0w8kHJKIL4iiGfq4b8lBFmWrWq2+mcybsIgwJs3tw5Hn98RsnWxAwnusdP3fBklnfTWuSNUQgT64aRKyWjCApJGCrEXEk3wRtWJXf44c5x991Tx4tjOhBb+TMMeNoT/hN39dDhwAZ3CGkmcJKItHNXGJ3khhlGJCnZVUOFsDUU1srSD35X+oRN6aWT+9xzzwXYF2d1aDXO6p/xwBMicgHwIXC2G30GcCpQBXwF/NQpQ9eJyE04JtoBxqrqusAfxjAKRMkqgnxW0mGrHLPBryIIq7JLJvekSZOYPHny26paEXfrhPi47nzBxYnyUdWJwMQcxTSMUNDoh4ZyqagaQ4WeLWGxdZQpYVVMhhFmGr0iiJFNhZZtpVJW5hwvTtiWjAZRHRqKEVa5DCOMlNzQUD5dVMZo1gy2bo22vZtMegRhaoWHSRbDiAq+qioRGSAi77o2V8YkuH++iKwVkUr3d6Hn3h4i8pyILBeRZSJSHqD8oaVp02i3SrOZIwjD80Z1SMswiknaHoGINAXuAfrj7KCc59pVWRYXdYqqjkyQxSPAzao6S0RaA1m4Pc8daylmRtQr1KjKbRjFwE+PoA9QpaorVfUbYDKODZa0iMj+QDNVnQWgqhtV9auspTUKTlTnCAzD8I8fReDXrsoZrqneqSIS287TA/hcRJ4SkbdE5Da3h1GPQmzDj1VYf/qTs3O2TZu8FNNoOPhg6NDBcZeZjjD1tsIki2FEhaCmM58BylW1F44lxpglx2bAMcAVwKHAXsD58YkLuQ1/2DDYtMmfB6tSrlR22gnWroVjj/WfJgy9g6gPaRlGMfCjCNLaVVHVGlXd4l7+FTjEPa8GKt1hpa3ANOBgQo5VIunp1cs5xkw6hA37GxqGf/wognlAdxHpJiLNgSE4Nli2EzPY5XI6sNyTtq2IxJr5xwPxk8yho2mDwSsjnuOOgw8/hHPPhb32csK++93iygSl3YszjGxJO0CiqltFZCSOUa2mwERVXSoiY4H5qjod+JWInI7j8GMd7vCPqm4TkSuAF0REgAXA/fl5lGTyZ55m4kS49dbMhkVKkT32cI5DhkCXLnDMMcWVx4v1CAzDP742lKnqDBwDXN6w6zznVwFXJUk7C+iVg4yBkEnF0LUr3HVX/mRpbIiA4/jLMIwoEuG9r4bREJssNozMKTkTE0b4efVVx0RHNpgiMIzMMUXgoba2lurqajZv3lxsURo1LVu2pGvXrpTFrPPFcfTRBRbIMEqcRq8IzjwTliyBXXdNH7e6upo2bdpQXl6OWJMyL6gqNTU1VFdX061btzzkH3iWhtHoafRzBNdeC59/Dp06pY+7efNm2rdvb0ogj4gI7du3z3uvy/6EhuGfRq8ImjSBnXf2H9+UQP7J5zu2OQLDyJxGrwiM0sIUgWFkjimCCNCvXz/mz5+ft/xffvll5syZk7f8DcMIN41+sthIz8svv0zr1q058sgjG9zbunUrzfxY6AsJNlnsj/Ix/yq2CEaIiM5/eIH59a+hsjLYPHv3dsxgJ2PTpk2cffbZVFdXs23bNq699lrOOeecenEmTZrELbfcgqoycOBAbr31VgBat27N8OHDee655/jOd77D5MmT6dixI++//z4XX3wxa9eupVWrVtx///3su+++2/NbtWoV9913H02bNuXRRx/lrrvu4oEHHqBly5a89dZbHHXUUQwZMoRLL72UzZs3s8MOO/Dggw+yzz77sG3bNkaPHs2zzz5LkyZNGD58OJdccgkLFizgN7/5DRs3bqRDhw489NBDdO7cmUJiQ0OG4R9TBCHi2WefZbfdduNf/3Jaaxs2bKh3/5NPPmH06NEsWLCAdu3acdJJJzFt2jQGDx7Mpk2bqKio4Pbbb2fs2LHceOON3H333YwYMYL77ruP7t2788Ybb3DRRRfx4osvbs+zvLycX/ziF7Ru3ZorrrgCgAceeIDq6mrmzJlD06ZN+eKLL3j11Vdp1qwZzz//PFdffTV///vfmTBhAqtWraKyspJmzZqxbt06amtrueSSS3j66afp2LEjU6ZM4be//S0TJ04syDu0HoFhZI4pgiSkarnni549e3L55ZczevRoTjvtNI6Js+I2b948+vXrR8xnw3nnncfs2bMZPHgwTZo02d57GDp0KD/84Q/ZuHEjc+bM4ayzztqex5YtW/DDWWedRVPXDOuGDRsYNmwYK1asQESora0F4Pnnn+cXv/jF9qGjXXbZhSVLlrBkyRL69+8PwLZt2wraG7DJYsPIHFMEIaJHjx4sXLiQGTNmcM0113DCCSdw3XXXpU+YABHh22+/pW3btlRmMca14447bj+/9tprOe644/jHP/7BqlWr6NevX9J0qsoBBxzA66+/noXUwWGKwDD8Y6uGQsQnn3xCq1atGDp0KKNGjWLhwoX17vfp04dXXnmFzz77jG3btjFp0iSOdW1lf/vtt0ydOhWAxx9/nKOPPpqddtqJbt268eSTTwJOJb1o0aIG5bZp04Yvv/wyqVwbNmygSxfHO+lDDz20Pbx///785S9/YatrGGjdunXss88+rF27drsiqK2tZenSpVm+EcMwCoH1CELE4sWLGTVqFE2aNKGsrIx777233v3OnTszfvx4jjvuuO2TxYMGDQKcFvybb77JuHHj6NSpE1OmTAHgscce45e//CXjxo2jtraWIUOGcNBBB9XL9/vf/z5nnnkmTz/9NHclsL995ZVXMmzYMMaNG8fAgQO3h1944YW899579OrVi7KyMoYPH87IkSOZOnUqv/rVr9iwYQNbt27l17/+NQcccEDQryshO+zgHJtYE8coENmswFo1fmD6SAVENGSzaxUVFZrPNfOpWL58Ofvtt19Rys6V1q1bs3HjxmKL4Zt8veuPPoIHHoDrr08+PCQiC1S1IvDC01DMbzseWz5aXPKhCHL5rq1HYDQqdt8dbrih2FIYRrTw1YEWkQEi8q6IVInImAT3zxeRtSJS6f4ujLu/k4hUi8jdQQlu1CdKvQHDMMJF2h6BiDQF7gH6A9XAPBGZrqrxTuinqOrIJNncBMzOSdICoapmeC7PhG040jBKHT89gj5AlaquVNVvgMnAIL8FiMghwK7Ac9mJWDhatmxJTU2NVVR5JOaPoGXLlsUWxTAMFz9zBF2AjzzX1cBhCeKdISJ9gfeAy1T1IxFpAvwBGAqcmKwAERkBjADYY489fIoePF27dqW6upq1a9cWTYZSIOahzDCMcBDUZPEzwCRV3SIiPwceBo4HLgJmqGp1quEWVZ0ATABnZUVAMmVMWVlZXrxmGYZhhBk/Q0MfA7t7rru6YdtR1RpVjdku+CtwiHt+BDBSRFYBvwd+IiLjc5LYMEJCukUUhhEV/PQI5gHdRaQbjgIYAvzIG0FEOqvqavfydGA5gKqe54lzPlChqvYPY0SeDBZR5BXbD2AEQVpFoKpbRWQkMBNoCkxU1aUiMhaYr6rTgV+JyOnAVmAdcH4eZTaMMLB9EQWAiMQWUWSlCKxCLy2y/Xvna0dy6HYWi8ha4MMktzsAnxVQnFyJkrylJOueqtoxFwFE5ExggKpe6F7/GDgsfgm1dyEEsA/wbgbFROlv4iWKckdRZqgvd9bfdeh2Fqd6EBGZXwzTANkSJXlN1vzgXQiRKVF6Ti9RlDuKMkNwcptpLsPIjrSLKAwjKpgiMIzs2L6IQkSa4yyimF5kmQwjK0I3NJSGrLrYRSRK8pqsGZBsEUXAxRT9ObMkinJHUWYISO7QTRYbhmEYhcWGhgzDMEocUwSGYRglTmQUQVi284vIKhFZ7PpdmO+G7SIis0RkhXts54aLiNzpyvy2iBzsyWeYG3+FiAwLSLaJIvKpiCzxhAUmm4gc4j57lZs2a3vdSWS9QUQ+9vi1ONVz7yq33HdF5GRPeMLvwp3EfcMNn+JO6IaGXHx8iMg2T3jBJqj9/A+KyNkiskxElorI457wwL93v+QodyjftYjc7pHrPRH53HMv83etqqH/4UzGvQ/sBTQHFgH7F0mWVUCHuLD/A8a452OAW93zU4F/AwIcDrzhhu8CrHSP7dzzdgHI1hc4GFiSD9mAN9244qY9JWBZbwCuSBB3f/dv3gLo5n4LTVN9F8ATwBD3/D7gl8X+jjP5nnF259+dJP3GkMrcHXjL8710yuf3nm+5w/yu4+JfgrNYIet3HZUeQU4+EQrAIByLq7jHwZ7wR9RhLtBWRDoDJwOzVHWdqq4HZgEDchVCVWfjmPgIXDb33k6qOledL+4RT15ByZqMQcBkVd2iqh8AVTjfRMLvwu2pHA9MTfDcYSDs33Mi/Mg8HLjH/W5Q1U/d8Lx87wWQu1hk+n2cC0xyz7N611FRBIl8InQpkiwKPCciC8QxHwCwq9YZ3fsfjiMeSC53IZ8nKNm6uOfx4UEz0h2qmhgbxspC1vbA56q6Nc+yZovfv/8Z7ruYKiLezWstRWS+iMwVkcH5FNSDH5l7AD1E5DVXtgEZpM0XucgN4X3XAIjInji95BczTeslavsIwsDRqvqxiHQCZonIO96bqqoiEso1uWGWzeVeHLem6h7/APysqBIVj2Q+PsCxKfOxiOwFvCgii1X1/aJJWkcznGGWfjg7rWeLSM+iSuSPhHKr6ueE913HGAJMVdVtuWQSlR5BaLbzq+rH7vFT4B843bg17tAJ7jHWtUwmdyGfJyjZPnbP8yazqq5R1W2q+i1wP867zUbWGpyhrmZx4WEhFx8f3m9wJfAy8L18Cuvi55utBqaraq07hPceTgVbzP/fXOQO87uOMYS6YaFM09ZR6ImQbH44GnslThcoNnlyQBHk2BFo4zmfgzP+dhv1J2T/zz0fSP0J2Te1bkLnA5zJnHbu+S4ByVhO/QnYwGSj4WTxqQHL2tlzfhnOvADAAdSfLF6JM6GW9LsAnqT+ZPFFxf6OM/me497FD4C57nk7oIV73gFYQQEWTviUeQDwsEe2j3CG6fL2vedZ7tC+azfevjiLV8QTltW7Lvo/RQYv51QcTf0+8NsiybCX+0dZBCyNyeF+NC+4H8rz1FWcguO85H1gMY5jnlheP8OZ9KwCfhqQfJOA1UAtTivngiBlAyqAJW6au70fYECy/s2V5W0cuz3eyvC3brnv4lmtlOy7cP9Wb7rP8GTsHzosv0RyA2OB093z37nf2CLgJWBfN/xI9x0tco8XhEhmAf6I45NhMa4iztf3nm+5w/yu3esbgPEJ0mb8rs3EhGEYRokTlTkCwzAMI0+YIjAMwyhxTBEYhmGUOKYIDMMwShxTBIZhGCWOKQLDMIwSxxSBYRhGifP/WrfZWJNDabQAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Sample from the posterior\n", "\n", "fuf2.sampleEMCEE2(lm, pargs=(x, d, yerr), dbfile=\"chain1.emcee\", sampleArgs={\"burn\":300, \"iters\":1000})\n", "\n", "\n", "ta = fuf2.TraceAnalysis2(\"chain1.emcee\")\n", "\n", "# Calculate mean, median, standard deviation, and\n", "# credibility interval for the available parameters\n", "for p in ta.availableParameters():\n", " hpd = ta.hpd(p, cred=0.95)\n", " print(\"Parameter %5s, mean = % g, median = % g, std = % g, 95%% HPD = % g - % g\" \\\n", " % (p, ta.mean(p), ta.median(p), ta.std(p), hpd[0], hpd[1]))\n", "\n", "ta.plotTraceHist(\"slope\")\n", "ta.show()" ] } ], "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.9.9" } }, "nbformat": 4, "nbformat_minor": 4 }