{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Class 6: Bayesian Methods\n", "\n", "In this class we will review Bayesian likelihood methods for solving statistical problems, determining the posterior probabilities of model parameters, and selecting between two models. At the end of these activities you should be able to:\n", "\n", "- Understand the application of Bayes' theorem in model-fitting and the role of priors\n", "- Obtain parameter values and confidence ranges via likelihood methods\n", "- Search parameter space with MCMC algorithms\n", "- Apply model selection using the Bayes factor or Akaike information criteria\n", "\n", "## Bayesian methods\n", "\n", "Bayesian statistics is a framework that allows us to **assign probabilities to a model**. It makes use of conditional probabilities $P(A|B)$, meaning \"the probability of $A$ on the condition that $B$ has occurred\".\n", "\n", "An important role in Bayesian statistics is played by **Bayes' theorem**, which can be derived from elementary probability:\n", "\n", "$P(A|B) = \\frac{P(B|A) \\, P(A)}{P(B)}$\n", "\n", "Bayes' theorem can be re-written for science in the form:\n", "\n", "$P(model|data) = \\frac{P(data|model) \\, P(model)}{P(data)}$\n", "\n", "- $P(model|data)$ is the **posterior** probability of the model in light of the data\n", "- $P(data|model)$ is the **likelihood** function of the data given the model\n", "- $P(model)$ is the **prior** probability of the model\n", "- $P(data)$ is the **evidence**, which can typically be absorbed into the normalization of the posterior\n", "\n", "### The likelihood\n", "\n", "For Gaussian variables, the **likelihood** can be written as\n", "\n", "Likelihood $\\propto e^{-\\chi^2/2}$\n", "\n", "### The prior\n", "\n", "Bayesian statistics cannot determine probabilities of a model without assigning a **prior probability**. The importance of the prior probability is both the strong and weak point of Bayesian statistics. Adopting a **uniform** (or constant) prior is effectively equivalent to **fitting range** of a parameter:\n", "\n", "$P(a|data) \\propto P(data|a) \\, P(a) \\propto e^{-\\chi^2/2}$\n", "\n", "## Posteriors and confidence limits\n", "\n", "We can use the posterior probability distribution $P(a)$ to determine **summary statistics** and **confidence intervals** for a parameter $a$:\n", "\n", "- Mean $= \\mu_a = \\int_{-\\infty}^\\infty a \\, P(a) \\, da$\n", "- Variance $= \\sigma_a^2 = \\int_{-\\infty}^\\infty (a - \\mu_a)^2 \\, P(a) \\, da$\n", "\n", "Only if the probability distribution is Gaussian is the mean equal to the best-fitting value, and the standard deviation equal to the 68% confidence region. For a general probability distribution, we can determine the confidence regions by numerically determining the range $a_{bot} < a < a_{top}$ which encloses 68% of the probability, such that\n", "\n", "- $\\int_{-\\infty}^{a_{bot}} a \\, P(a) \\, da = 0.16$\n", "- $\\int_{a_{bot}}^{a_{top}} a \\, P(a) \\, da = 0.68$\n", "- $\\int_{a_{top}}^\\infty a \\, P(a) \\, da = 0.16$\n", "\n", "## Marginalization\n", "\n", "Now suppose we have determined the 2D posterior probability distribution of a 2-parameter fit, $P_{2D}(a,b)$.\n", "\n", "What is the probability distribution for parameter $a$, considering all possible values of parameter $b$? This is known as **marginalization** of parameter $b$.\n", "\n", "Marginalization can be performed by **summing (integrating) over one axis of the probability distribution**, for example:\n", "\n", "$P_{1D}(a) = \\sum_b P_{2D}(a,b)$\n", "\n", "If $P_{2D}(a,b)$ is normalized, then $P_{1D}(a)$ will also be normalized.\n", "\n", "## Worked example\n", "\n", "Let's solve our example of fitting a straight line $y = ax + b$ to some data, using a Bayesian likelihood approach. Here is the best-fit from before:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(-0.5, 4.0)" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEKCAYAAAA4t9PUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3de3yU5Zn/8c9lQNDQQhvQIketWJe6FTQlba1bigeIUuiuWqXaqgWp+1vU7tZutd2VqYftam1XSw+WDQiKjVq0lVJodBVitTUSji2gFVEUSCUiHhJRILl+f9wTEoYcSCaZ55nM9/165cXMM09mLvNCvrnv+3nuy9wdERGRjjos6gJERCS7KUhERCQtChIREUmLgkRERNKiIBERkbQoSEREJC2xCRIzyzOz1Wa2uJnXepnZA2a2ycwqzGx45isUEZHmxCZIgGuAjS28NhXY5e7HA/8D3JqxqkREpFWxCBIzGwycC5S0cMpkYH7y8ULgDDOzTNQmIiKt6xF1AUl3AP8OfKCF1wcBrwK4+z4zewsoAF5vepKZTQemA+Tn55964okndlnBIiLd0cqVK1939wHt+Z7Ig8TMJgI73H2lmY1N573cfTYwG6CwsNArKys7oUIRkdxhZlva+z1xmNo6DZhkZi8D9wPjzGxByjnbgCEAZtYD6AvszGSRIiLSvMiDxN2vd/fB7j4cuAh4wt0vSTltEXBp8vH5yXO026SISAxEPrXVEjO7Eah090XAHOBeM9sEvEEIHBERiYFYBYm7LweWJx/f0OT4e8AF0VQlIiKtiXxqS0REspuCRERE0qIgERGRtChIREQkLQoSERFJi4JERETSoiAREZG0KEhERCQtChIREUmLgkRERNKiIBERkbQoSEREJC0KEhERSYuCRERE0qIgERGRtEQeJGbW28yeNbO1ZrbezL7XzDmXmVm1ma1Jfk2LolYRETlYHBpbvQ+Mc/caM+sJPGVmS939mZTzHnD3GRHUJyIirYg8SJK912uST3smv9SPXUQkS0Q+tQVgZnlmtgbYATzm7hXNnHaema0zs4VmNiTDJYqISAtiESTuXufuo4DBwBgzOynllN8Cw939E8BjwPzm3sfMpptZpZlVVldXd23RIiIxkkgkMLP9X4lEImOfbWFmKT7M7AbgXXe/vYXX84A33L1va+9TWFjolZWVXVGiiEgsjR07FoDly5d3+D3MbKW7F7bneyIfkZjZADPrl3x8BHAW8FzKOQObPJ0EbMxchSIi0prIgwQYCCwzs3XACsIayWIzu9HMJiXPuTp5afBa4GrgsohqFRE5QJRTSnERu6mtzqKpLRHJlM6YUopLHVk5tSUiItlNQSIiImlRkIhI1tL6RDxEfme7iEhHJRKJ/esBUa9P5DKNSEREJC0KEhERSYuCRERE0qIgERGRtChIREQkLQoSERFJi4JERETSoiAREcl1e/fCb34DEyd26Nt1Q6KISK564QWYMwfmzYPXXoOBA9v8luZoRCIikkt274YFC2DsWDjhBLj9dvjUp+C3v4VXXunQW2pEIiKSC1avDqOPBQvgrbfgox+F738fLr20wyORBgoSEZHu6q234Je/hJISWLUKevWC88+HadPgH/4BDuucSanIg8TMegNPAr0I9Sx095kp5/QC7gFOBXYCF7r7yxkuVUQkturq6ti5cyc177zDH2+9lU+tX89hCxeGqayTT4ZZs+Dii+FDH+r0z448SID3gXHuXmNmPYGnzGypuz/T5JypwC53P97MLgJuBS6MolgRkbipq6vjwrFjKV6/nqnufOy666jNy+OIadM47Ior4JRTwKzLPj/yxXYPapJPeya/Uvv/TgbmJx8vBM4w68KfikiMqQeH7FdXB0uWsOP00yl96iluc2cHcClwXO/eLJk4EU49tUtDBOIxIsHM8oCVwPHAT929IuWUQcCrAO6+z8zeAgqA11PeZzowHWDo0KFdXbZIJNSDQ3j5ZZg7F+6+G7Zu5YP5+dwJlADPJ0+xd99lzZo1TOzgvSHtEfmIBMDd69x9FDAYGGNmJ3XwfWa7e6G7Fw4YMKBzixQRaUbD2sSWLVtYvHgxdXV1XfNB778PDz4IZ58Nxx0HN98MJ50ECxey/N57+V6fPvtDBCA/P59Ro0Z1TS0pYhEkDdz9TWAZMCHlpW3AEAAz6wH0JSy6i4hEpq6ujvHjx7NhwwZefvllpkyZwvjx4zs3TNavh3/7Nxg0CC68EJ5/HhKJMCpZuhTOO48JkyZRVFTEYcmrsPr06UNRURHFxcWdV0crIg8SMxtgZv2Sj48AzgKeSzltEWHaD+B84Al3T11HERHJqKVLl1JRUUF9fT0ANTU1VFRUsHTp0vTeuKYmTF195jNh1PGTn8C4cVBWBps3ww03QJPp+7y8PMrKyhg5ciTDhw+ntLSUsrIy8vLy0qvjEEUeJMBAYJmZrQNWAI+5+2Izu9HMJiXPmQMUmNkm4N+A6yKqVURiJGPTSi1YvXo1tbW1Bxyrra1lzZo17X8zd3j2WZg+PdwgOHUqvPkm/PCHsG1b47RWC+GQl5dHQUEBw4YNY+LEiRkLEYjBYru7rwNGN3P8hiaP3wMuyGRdIhJvTaeV6uvrmTJlCkVFRRn9TXz06NHk5+dTU1Oz/1i71yZ27gx3m8+ZA3/+Mxx5ZJjCmjYNPv3pLr/iqjPEYUQiItJuXTat1A7FxcUdW5uor4fHH4cpU+CYY+Ab34AjjoDZs6GqqnFaKwtCBBQkIpKlOnVaqYPavTaxbRvccgscfzyceWZY87jySli7Fioq4Ior4IMfzFj9nSXyqS0RkY7olGmlTtCwNlFQUND8PRt798KSJWG/qyVLwmhk3LgQKP/4j9C7d0br7QoKEhHJSg3TSsuWLaO+vj7jl7y2qaHXx/z58Le/hQX0666Dr30t7LzbjShIRCQrNUwrjRo1ipqaGmbNmkVxcXFGr1Y6yO7d8NBDYfRRXh6usDr33LBwXlwMPbrnP7nd879KRHJCm9NKGXJ8TQ3nVlWFUUdDr4//+q/Q6+OYYyKrK1MUJCIiHfHWW1BaCiUllKxcyR4z+PKXO73XRzZQkIiIHCp3eOqpMHX1q1+FqaxPfII7jz+ex446isULFkRdYSRyJzJFRDrqtdfgBz+AE08Mo41f/xq++lVYsQLWrOHXgwZR07Nn1FVGRiMSEZHm1NXBo4+G0ceiRbBvH3z2s3D99XDBBZCfH3WFsaEgERFp6uWXQ5+PuXNh61YYMCDceT51ahiRyEEUJCIi778fRh0lJfDYY+HY+PFwxx3whS/A4YdHW1/MKUhEJHdt2BBuGrznHnj99bA1+8yZcPnlB2zTLq1TkIhIbqmpCVuyl5TAn/4EPXvC5Mnhst0zz2xxm3Zpma7aEskyUffgyErN9frYtQtuvz2sg/zqV2EqSyHSIZGPSMxsCHAPcDTgwGx3vzPlnLHAI8BLyUMPu/uNmaxTJA7i0IMjq7zxRuj1UVKStb0+skHkQQLsA77p7qvM7APASjN7zN03pJz3B3ePbg8EkRhorQdHlFuExEp9PSxfHsLj4YfDQvonPwm/+AVcdFFWbtMed5EHibtXAVXJx++Y2UZgEJAaJCI5r7UeHDkfJNu2wbx5YfH8pZegX78wlTV1Kpx8ctTVdWuRB0lTZjac0Ha3opmXP21ma4HtwLXuvr6Z758OTAcYqisupBuKSw+O2Ni7l9Nefz1smDh0aBiNfP7zcPPNodfHEUdEXWHGJBIJysvLATAzZs6cSSKRyMhnm7tn5IPaYmZ9gHLgFnd/OOW1DwL17l5jZucAd7r7iNber7Cw0CsrK7uuYJEINKyRpPbgyLk1khdeCDcMzpsHf/sbrx9+OP2vvTayXh9jx44FYPny5Rn/7M5mZivdvbA93xOLEYmZ9QQeAu5LDREAd3+7yeMlZvYzM+vv7q9nsk6RqMWyB0em7N4d1jxKSsIaSLLXx/UvvsizBQU8fsstUVeYsyK//NfMDJgDbHT3H7VwzkeS52FmYwh178xclSLx0dCDY9iwYUycOLH7h8iaNTBjRujrcckl8OqrodfHK6/AI4/wp/79qdPVV5GKw4jkNOArwJ/NbE3y2HeAoQDufhdwPvDPZrYP2A1c5HGZkxORztek1wcrV0KvXnDeeeGy3c99Lla9PqJcm4iLyIPE3Z8CWv11wt1/AvwkMxWJSCTc4emnQ3g8+OD+Xh/MmhUaRn34w1FX2KxEIpFzwZEq8iARkRy3Y0fY66qkBJ5/Hj7wgdDrY9o0OPVU3TSYBeIzPhSR3FFXB7//PZx/PgwaBN/6FvTvH7Zvr6qCu+6CwsI2Q6RhWqm8vBwzy/mRQVQ0IhGRzNmyJVy2e/fdYdG8f3+45ppw0+Df/V27307TSvGgIBGRrtVcr4+zz4Yf/QgmTVKvj25AQSIiXSO118eQIXDDDaHXx7BhUVcnnUhBIiKdp6YmbMleUgJ//KN6feQIBYmIpMcdVqwI4XH//fDOO6G3+e23w1e+AkcdFXWF0sUUJCLSMc31+vjSl8Lo4zOf0WW7OURBIiKHrqVeH3fdFXp99O0bdYUSAQWJiLRt+/bGXh+bN6vXhxxAQSIizdu7F5YsCeHxu9819vq46aac6/UhrVOQiMiBNm0K4ZHs9cFHPgLf/nbo9XH88VFXJzGkIBGR5nt9nHNOWDg/5xzooX8qpGX62yGSy9auDeGxYAG8+SYcd1zo9XHppaH/h8ghUJCI5Jq3327s9VFZGeteH5IdIv8bY2ZDzGyZmW0ws/Vmdk0z55iZ/djMNpnZOjM7JYpaRbJWQ6+Pyy+HgQPhyivDpbs//nG4Iuu++8JCukJEOiAOI5J9wDfdfZWZfQBYaWaPufuGJucUAyOSX0XAz5N/ikhrUnt99OkT2tVOm3ZI27SLHIrIf/1w9yp3X5V8/A6wERiUctpk4B4PngH6mdnADJcqEgtt9uBortdHQUHYvr2qCn7xi3AToUJEOonFqfW5mQ0HngROcve3mxxfDPx3si0vZvY48G13r0z5/unAdIChQ4eeumXLlgxVLhIDW7aEPh9z5zb2+rj00g73+pDcZGYr3b2wPd8Th6ktAMysD/AQ8I2mIdIe7j4bmA1QWFgYn4QU6Sp79jT2+nj00XBMvT4kw2IRJGbWkxAi97n7w82csg0Y0uT54OQxkdykXh8SI5EHiZkZMAfY6O4/auG0RcAMM7ufsMj+lrtXZapGkViorYUHH2zs9dGjR2Ovj7POUq8PiUzki+3AacBXgHFmtib5dY6ZXWlmVybPWQJsBjYB/wv8v4hqlRyWSCQws/1fGekV3tDr4+tfD5ftfu1rYfv222+Hbdtg4UKYMEEhIpGK1WJ7ZyosLPTKysq2TxRph7FjxwKwfPnyrv2gN94I93aUlMC6dWGDxAsvVK8P6XJZvdgukvOa6/VRWKheHxJ7cZjaEmlTJNNKmbJ9e9jfasQIOOMMWLoUrrgCVq9unNZSiEiMaUQiWSGRSOyfTuryaaVM2Lcv9PooKVGvD8l6ChKRTNq0KdwwOG9euMtcvT6kG1CQiHS1995r7PWxbFnYGPHcc8PCeXEx9OwZdYUiaWkzSMzsMeBad1+bgXpEuo+1a8NNgwsWwK5dodfHLbeEbUsGpW4nJ5K9DmWx/dvAHWZ2tzZKFGnD2283boo4alR4PGECPP44vPACfOc73SJEuvXFD9JubY5Ikjvzft7MzgN+b2YPA7e5++4ur04kG7iHO81LSsKd5+++C3//93DnnXDxxWHn3W6m2138IGk5pMt/k9uYPE/oA3IV8IKZfaUrC5N40G+erdixA374Qxg5Ej772XCX+SWXwLPPhmmtq6/uliEikupQ1kieBo4F1gPPAJcBzwHXmNnp7j69SyuUSOk3zwMd5k7hrl1wwQXwyCOwd2+403zu3HCsT5+oSxTJuEO5ams6sMEP3kvlKjPb2AU1icRPstdHaUUFR7//frh096qrQq+PkSOjrk4kUoeyRrK+lZfP7cRaROKlmV4fW/r142cf/SjfW7UKevWKuECReEhrixR339xZhYjExsaNcO214eqqCy4IvT9uuAFeeol//8QnKB8wQCEi0oRuSBSBxl4fc+bA00+r14dIOyhIJHe5Q2VlmLoqLYV33oGPfQx+8AP46lfhqKOirlAkK8QiSMxsLjAR2OHuJzXz+ljgEeCl5KGH3f3GzFUo3UpzvT6+9KUw+jjtNPX6EGmnWAQJMA/4CXBPK+f8wd0nZqYc6Xbq66G8PITHQw+FXh+nnqpeHyKdIBZB4u5PmtnwqOuQbmj79rDT7ty58OKL0K9f6PUxdWrYwkRE0pZNja0+bWZrzWypmX28uRPMbLqZVZpZZXV1dYc/SHdzZ7l9+8Jlu5MmwdCh8N3vwpAhYfPE7dth1iyFiEgnisWI5BCsAoa5e42ZnQP8BhiRepK7zwZmQ+jZ3tEP093cWerFF8NVV017fXzrW6HXx4iD/rqISCfJihGJu7/t7jXJx0uAnmbWP+KyJA7eew9++UsYNy40hrr11rD28ZvfwCuvwPe/32khUldXx86dO9myZQuLFy+mrq6uU95XJNtlxYjEzD4CvObubmZjCAG4M+KyJErr1oWF84ZeH8ce26W9Purq6hg/fjwbNmygvr6eKVOmUFRURFlZGXm6x0RyXCyCxMxKgbFAfzPbCswEegK4+13A+cA/m9k+YDdwUTN7f0k3VldXx3s7djBuxw7ePOEE+r3wAhx+OJx3Xrhsd+zY0HmwiyxdupSKigrq6+sBqKmpoaKigqVLlzJxoi4mlNwWiyBx9yltvP4TwuXBkmvcqfvDH3jsoot4vKqKfGD9rl0sOOEE/vkPfyAvQzcNrl69mtra2gOO1dbWsmbNGgWJ5LysWCORHFRdvb/XR97nPsdpVVXcB4wBTqqv5/rt21n67LMZK2f06NHk5+cfcCw/P59RuvpLREEiMVJXB2VlYaPEQYPCxokf/jCLvvhFjgG+DqxIntowGsiU4uJiioqKOCw5fdanTx+KioooLi7OWA0icaUgkei98gp873tw3HGhv/myZaHXx/r18PTTHDZ16kENozI9GsjLy6OsrIyRI0cyfPhwSktLtdAukqQgkWjs2RNa006YAMOHhyA58cSwA++2bY0tbInPaCAvL4+CggKGDRvGxIkTczpEdCm0NBWLxXbJIRs3hpsG77knrIMMHgz/+Z9w+eUhUJrRMBoYNWoUNTU1zJo1i+Li4pz+hzxKuhRaUilIpOvV1sKvfhXu+2jo9TFpUrhs9+yzD6nXR8NooKCgQFdJRUyXQksqTW1Jqzo8hdHQ6+PKK2HgwDDiqK6G226DrVvDDrzFxWoYlYVauxRacpNGJNKiDk1hqNdHt9dwKXRNTc3+Y7oUOrdpRCItam0K4wD19eFKq4svhmOOgauvhp494ec/D5snzpsHn/2sQqSbiMvFDxIfGpFIi9q8m7shJObMCTvv9u0bRh5Tp8Lo0dEULV1OFz9IKgWJtKi5KYwPHnkkZ+/ZA5Mnw+9+F24i/NznIJEI+14dcUR0BUvG6OIHaUpTW9KiplMYxwG39ezJpr17GXPTTVBREXp9/PWvsHw5XHKJQkQkR2lEIi3K27uXRy+7jGefeYZP1dbidXXhSqsrroBzzgnrICKS8xQkcrAmvT4O27WLob17UzJ8ONOeeqpLen2ISHbT1JYEb78Ns2fDmDFw8snwi1/A+PHwf//HxWPGsGDYMIWIiDQrFkFiZnPNbIeZ/aWF183Mfmxmm8xsnZmd0pX15Mw+Qu7wxz+GnuYDB8LXvw67d8Odd8L27VBaCmecgeuyXRFpRVymtuYRGlfd08LrxcCI5FcR8PPkn50uJ/YRqq6Ge+8N01cbN4addS++OFy6+8lP6n4PEWmXWIxI3P1J4I1WTpkM3OPBM0A/MxvYFbUc8k142Sa118c3vwn9+oV7QKqqGqe1FCIi0k6xCJJDMAh4tcnzrcljBzCz6WZWaWaV1dXVHfqgbrePUHO9PmbMgL/8pXFaK6XXh4hIe8RlaqtTuPtsYDZAYWGhd+Q9usU+Qnv2wG9/G6auysrCWshZZ8EPfhBuJOzVK+oKRaQbyZYRyTZgSJPng5PHOl1W7yP03HPhJsHBg+H888Oo4z/+AzZvhkcfDZsnKkREpJNly4hkETDDzO4nLLK/5e5VXfFBWbePUCf0+hARSUcsgsTMSoGxQH8z2wrMBHoCuPtdwBLgHGAT8C5weVfWE/t9hNxh5coQHr/8JbzzDpxwQuj18dWvwtFHR12hiOSQWASJu09p43UH/iVD5cTXrl2NvT7Wrg17W11wQRh9aJv2LpdIJCgvLwfAzJg5cyaJRCLaokRiIFvWSHJSIpHAzBhrxgIz9g0YAFddFaavfvazcNPg/Plw+ukKkQxIJBK4+/4vhYhIEIsRiTSjqopE795MO+IIBu/eHXp9XHKJen2ISOwoSOJk3z5YujRMXSV7fbzety/zhw3ju6tW5fQ27ZpWEokvTW3FwYsvwne/C0OHhiuuKirg2mvh+ef5xqhRPHb00TkdIqBpJZE404gkKu+9B7/+dRh9PPEEHHZY6PExdSqce656fYhI1lCQZNq6dWF/q3vvDVdhHXss3HwzXHaZtmkXkaykIMmEt9+G++8Po48VK+Dww+Gf/ilctvv5z4fRiEgW0ZqVNKUg6Sru8Kc/hfB44AF49134+MfhjjvC1VcFBVFXKNJhiURCwSH7KUg6W2qvj/x8+PKXw+hD27SLSDekOZXOUF/fuCliQ6+Pvn1DmFRVwf/+LxQVZWWINExhlJeXY2b6LVREDqIRSTpeeQXuvjt8bdkSpqtmzAhXXn3841FX1yk0hSEibVGQtFdLvT5uu029PkQkJylIDtVzz4XLdufPD+sggwaFXh+XXx4u4RURyVEKklb0rqsLwVFSAk89FTZL/MIXwsL5+PHq9SEigoLkYO6wahX/+te/csaOHSFA1OtDRKRFsbhqy8wmmNnzZrbJzK5r5vXLzKzazNYkv6Z1ehG7dsFPfwqnnAKFhYx/7TWe6t8fnnyysYWtQkRE5CCRj0jMLA/4KXAWsBVYYWaL3H1DyqkPuPuMTv1wdygvD1NXDz0U9r865RT42c84/777qOnRgwmnn96pHyki0t1EHiTAGGCTu28GSPZlnwykBknnqaoKax9z5sCmTeGej699LVy2e8opANQ88ECXfbyISHcSh6mtQcCrTZ5vTR5LdZ6ZrTOzhWY2pLk3MrPpZlZpZpXV1dUHvrhvHyxeDF/8IgwZAtdfH668uvfe0GmwYVoL3YQnItIecRiRHIrfAqXu/r6ZfR2YD4xLPcndZwOzAQoLCx2AzZth7txw0+D27WGd49prwwjkhBOa/TDdhCcicujiECTbgKYjjMHJY/u5+84mT0uA29p81zfegDPOaOz1UVwcRh3q9SEi0qniMLW1AhhhZsea2eHARcCipieY2cAmTycBG9t815deCqORm24K25c0TGtlUYjU1dWxc+dOtmzZwuLFi6mrq4u6JBGRg0Q+InH3fWY2AygD8oC57r7ezG4EKt19EXC1mU0C9gFvAJe1+cYjRoTLdrO010ddXR3jx49nw4YN1NfXM2XKFIqKiigrKyNPN0KKSIyYu0ddQ5coLCz0ysrKqMvosMWLFzNlyhRqamr2H+vTpw+lpaVMnDgxwspEpDszs5XuXtie78nOX9dzwOrVq6mtrT3gWG1tLWvWrImoIhGR5ilIYmr06NHk5+cfcCw/P59Ro0ZFVJGISPMUJDFVXFxMUVERhyXXePr06UNRURHFxcURVyYiciAFSUzl5eVRVlbGyJEjGT58OKWlpVpoF5FYivyqLWlZXl4eBQUFFBQUaIFdRGJLIxIREUmLgkRERNKiIBERkbQoSEREJC0KEhERSYuCRERE0qIgERGRtChIREQkLQoSERFJi4JERETSEosgMbMJZva8mW0ys+uaeb2XmT2QfL3CzIZnvkoREWlO5EFiZnnAT4FiYCQwxcxGppw2Fdjl7scD/wPcmtkqRUSkJZEHCTAG2OTum919D3A/MDnlnMnA/OTjhcAZZmYZrFFERFoQhyAZBLza5PnW5LFmz3H3fcBbQEFGqhMRkVbFIUg6jZlNN7NKM6usrq6OuhwRkZwQhyDZBgxp8nxw8liz55hZD6AvsDP1jdx9trsXunvhgAEDuqhcERFpKg5BsgIYYWbHmtnhwEXAopRzFgGXJh+fDzzh7p7BGkVEpAWRd0h0931mNgMoA/KAue6+3sxuBCrdfREwB7jXzDYBbxDCRkREYiDyIAFw9yXAkpRjNzR5/B5wQabrEhGRtsVhaktakEgkKC8vp7y8HDMjkUhEXZKIyEGsuy41FBYWemVlZdRliIhkFTNb6e6F7fkejUhERCQtChIREUmLgkRERNKiIBERkbQoSEREJC0KEhERSYuCRERE0qIgERGRtChIREQkLQoSERFJi4JERETSoiAREZG0KEhERCQtChIREUlLpEFiZh82s8fM7IXknx9q4bw6M1uT/EptwysiIhGKekRyHfC4u48AHk8+b85udx+V/JqUufJERKQtUQfJZGB+8vF84IsR1iIiIh0Qdc/2o929Kvn4b8DRLZzX28wqgX3Af7v7b5o7ycymA9OTT983s790arXZqz/wetRFxIR+Fo30s2ikn0Wjj7X3G7o8SMzs/4CPNPPSd5s+cXc3s5b6/g5z921mdhzwhJn92d1fTD3J3WcDs5OfW9nedpHdlX4WjfSzaKSfRSP9LBolf2lvly4PEnc/s6XXzOw1Mxvo7lVmNhDY0cJ7bEv+udnMlgOjgYOCREREMi/qNZJFwKXJx5cCj6SeYGYfMrNeycf9gdOADRmrUEREWhV1kPw3cJaZvQCcmXyOmRWaWUnynL8DKs1sLbCMsEZyKEEyuysKzlL6WTTSz6KRfhaN9LNo1O6fhbm3tCwhIiLStqhHJCIikuUUJCIikpZuGSRmNsHMnjezTWbW0t3y3Z6ZDTGzZWa2wczWm9k1UdcUNTPLM7PVZrY46lqiZGb9zGyhmT1nZhvN7NNR1xQVM/vX5P8ffzGzUjPrHXVNmdtIsyYAAAMcSURBVGJmc81sR9N77g5166qmul2QmFke8FOgGBgJTDGzkdFWFZl9wDfdfSTwKeBfcvhn0eAaYGPURcTAncDv3f1E4GRy9GdiZoOAq4FCdz8JyAMuiraqjJoHTEg5dqhbV+3X7YIEGANscvfN7r4HuJ+wFUvOcfcqd1+VfPwO4R+LQdFWFR0zGwycC5S0dW53ZmZ9gX8A5gC4+x53fzPaqiLVAzjCzHoARwLbI64nY9z9SeCNlMPt3rqqOwbJIODVJs+3ksP/eDYws+GEGzkroq0kUncA/w7UR11IxI4FqoG7k9N8JWaWH3VRUUje7Hw78ApQBbzl7o9GW1XkDnXrqv26Y5BICjPrAzwEfMPd3466niiY2URgh7uvjLqWGOgBnAL83N1HA7UcwvRFd5Sc/59MCNdjgHwzuyTaquLDw/0hbd4j0h2DZBswpMnzwcljOcnMehJC5D53fzjqeiJ0GjDJzF4mTHeOM7MF0ZYUma3AVndvGJ0uJARLLjoTeMndq919L/Aw8JmIa4raa8ktq2ht66qmumOQrABGmNmxZnY4YeEsJ5thmZkR5sE3uvuPoq4nSu5+vbsPdvfhhL8TT7h7Tv7m6e5/A141s4ZdXs8gd7cdegX4lJkdmfz/5Qxy9MKDJtrcuipV1NvIdzp332dmM4AywhUYc919fcRlReU04CvAn81sTfLYd9x9SYQ1STxcBdyX/GVrM3B5xPVEwt0rzGwhsIpwleNqcmi7FDMrBcYC/c1sKzCTsFXVg2Y2FdgCfKnN99EWKSIiko7uOLUlIiIZpCAREZG0KEhERCQtChIREUmLgkRERNKiIBERkbQoSEREJC0KEpEMSvaHOSv5+GYzmxV1TSLp6nZ3tovE3EzgRjM7irAb86SI6xFJm+5sF8kwMysH+gBjk31iRLKaprZEMsjM/h4YCOxRiEh3oSARyZDkltz3Efpf1JhZaotTkaykIBHJADM7ktDr4pvuvhG4ibBeIpL1tEYiIiJp0YhERETSoiAREZG0KEhERCQtChIREUmLgkRERNKiIBERkbQoSEREJC3/HwVG2YbXHJmKAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from scipy import stats\n", "from scipy.optimize import minimize\n", "%matplotlib inline\n", "\n", "# Function which evaluates the chi-squared value of a model y=a*x+b for the (x,y) dataset\n", "def chisqstraightline2(p,x,y,yerr):\n", " a,b = p[0],p[1]\n", " return np.sum(((y - (a*x + b))/yerr)**2)\n", "\n", "# Read in the example dataset sampled from a straight line\n", "# You will need to change the file path to the location where you've saved the data\n", "data = np.loadtxt('../datasets/straightline.dat')\n", "x,y,yerr = data[:,0],data[:,1],data[:,2]\n", "\n", "# Find the best-fitting model parameters by minimising chi-squared (same code as in Class 3)\n", "p0 = [0.2,1.]\n", "result = minimize(chisqstraightline2,p0,args=(x,y,yerr))\n", "a,b = result['x']\n", "\n", "# Plot the dataset and best fit\n", "xmin,xmax = 0.,10.\n", "ymin,ymax = -0.5,4.\n", "plt.plot(x,y,marker='o',markersize=5,color='black',linestyle='None')\n", "plt.errorbar(x,y,yerr=yerr,color='black',capsize=2.,linestyle='None')\n", "plt.plot([xmin,xmax],[a*xmin+b,a*xmax+b],color='red')\n", "plt.xlabel(r'$x$')\n", "plt.ylabel(r'$y$')\n", "plt.xlim(xmin,xmax)\n", "plt.ylim(ymin,ymax)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We now compute $\\chi^2$ over a grid and convert this to a 2D probability distribution:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Text(0, 0.5, '$b$')" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAEKCAYAAAA1qaOTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAamElEQVR4nO3dfaxlVXnH8d/DDC8yA4LMRHmZcbSdluJLVW5Ba2JsxQaoYWzBFppWMJpJrcS++Ee1NNjYpFWbWG0g0gkQoS+itaaOcZSCotRYlIEMwkDVkdgwhMrw4vA64MDTP+7hznOXd6+77r7r7r3Pud9PMnHve/bZZ92di+v89rPW2ubuAgBgsQ7quwEAgMlAhwIAqIIOBQBQBR0KAKAKOhQAQBV0KACAKnrvUMxsnZndYGZ3mtlOM/vjOY4xM/sHM9tlZt81s9f00VYAQLOVfTdA0n5J73P3W83sCEm3mNl17n5nOOYMSRtH/06V9MnR/wIABqL3hOLu97n7raPtRyXdJen45LBNkq72aTdJOsrMju24qQCAjCEklBlmtkHSqyV9O3npeEn3hP3do5/dl7x/s6TNkrRq1aqTTzzxxKVqKgBMpFtuueUBd1/b5r2D6VDMbLWkf5f0J+7+SJtzuPsWSVskaWpqyrdv316xhQAw+czsf9u+t/dbXpJkZgdrujP5F3f//ByH3CtpXdg/YfQzAMBA9N6hmJlJukLSXe7+sYbDtkp6+2i012sl7XX3+xqOBQD0YAi3vF4v6Q8k3W5mO0Y/+wtJ6yXJ3S+TtE3SmZJ2SXpC0jt6aCcAIKP3DsXdvynJ5jnGJb2nmxYBANro/ZYXAGAy0KEAAKqgQwEAVEGHAgCogg4FAFAFHQoAoAo6FABAFXQoAIAq6FAAAFXQoQAAqqBDAQBUQYcCAKiCDgUAUAUdCgCgit6Xr8fCTK/kvzDTzzADgKVFQgEAVEFCGaA2KaTG+UgyABaDhAIAqIKEMhC1U0ntNpBeAMyHhAIAqIKEgiJN6YXkAuA5JBQAQBUklAlRowbTJm1QdwHwHBIKAKAKOhQAQBXc8hpjQ58Aye0wYHkhoQAAqhhEh2JmV5rZ/WZ2R8PrbzSzvWa2Y/Tv4q7bOM7cfdH/2pyvtE0AJsNQbnl9StIlkq7OHPNf7v6WbpoDAFioQXQo7n6jmW3oux3joPQbfe3jSsXaSDx3rmZCrQWYDIO45VXodWZ2m5l92cxe1ndjAACzDSKhFLhV0ovd/TEzO1PSf0jamB5kZpslbZak9evXd9vCnuW+5feVVkqVphfSCjBsY5FQ3P0Rd39stL1N0sFmtmaO47a4+5S7T61du7bzdgLAcjYWCcXMXiTpx+7uZnaKpjvCB3tu1lhqm2RKF4esnSiaajLztQNA9wbRoZjZpyW9UdIaM9st6YOSDpYkd79M0jmS3m1m+yU9KelcZ7wpAAzKIDoUdz9vntcv0fSw4onV5pt4jT616Rxt6y6lo7zaJApqLcCwjUUNBQAwfHQoAIAqBnHLC8MQbxu1vZ1WWryP+7Vvh+VuwwFYOiQUAEAVJJQJkRu+m/t5UyrJHVejfaUJpU2SKb0WJBegLhIKAKAKEsoA1R4a3EbbJNMkly4OOuigxtdyx5VqSiLUWoC6SCgAgCpIKGOsTZIprS+k4nHPPvts0ftLJzbG80mzk0jcTj8rvtYmySz1sjHAckNCAQBUQUIZuLaJYrFKayhpusi9FuUSRVNCSY+Ln5VLJbmRYk3HkVaAhSOhAACqoEMBAFTBLa8JUXpbpu3qwE3nSG9rPfPMM3O+tpBnmaxYsWLebWn2ba74Wq5436Zgz/BioAwJBQBQBQllzCx20mPtIn+aUOL+/v37G4/LfW5T8kgTysEHHzzn+dLjmgYU5AYDRCzlApQhoQAAqiChTKjSb8uLXYgxFZNIrKfE7XQ/V6OIqWHlytl/rvEcMa3Ebak55ZTWdUqHJFNrwXJHQgEAVEFCGWM1FpFsk1DajChLE8pPf/rTme1Ya0nfl0socf+QQw6Z2U4TyqGHHjrna2mtJSpNJbmfMzkSyw0JBQBQBQllQrT9dlx6XJsl5UsTStxOj82NyoppIyaUmEjS88XX4nvSz2qbXiKWcsFyQ0IBAFRBQlnmmr4tly7EWLqwYyqOBksTSq6+EsUUEdPKYYcdNuu45z3veXOeL01NTbWWVOm1iBgNhuWAhAIAqIIOBQBQBbe8JtRiC8Klw4bT2zzxNlTTExVT6a2np59+emb7ySefnNlOb41FscAeb3FJ0r59+2a2Dz/88Jnt9HZabEe8bZbeokqHLzdps5QLt78wzkgoAIAqBpFQzOxKSW+RdL+7v3yO103SJySdKekJSRe4+63dtnJ5Sb8pNz3bPZdQcgs75r6Jx+QQ00VMK+lxsR3pcOBVq1bNbB955JEz2zEJpedrWoZfml28rzHUmII9JsVQEsqnJJ2eef0MSRtH/zZL+mQHbQIALMAgEoq732hmGzKHbJJ0tU9/fbvJzI4ys2Pd/b5OGjjmanwDbqqbpN/K47f5WGtIh+GWLoHSlFbS/aa0Is2uhzz22GMz20cfffSs45omW6Y1nvg7pkOUSyxkWRfqKxgnQ0ko8zle0j1hf/foZ7OY2WYz225m2/fs2dNZ4wAAA0kotbj7FklbJGlqampxT46aUG2/AbcZ5RUTSjoyKrfcfNyPn5WOyoo1lZhW0ppHbNOjjz4657YkrVmzZmY71ldyCSVuxxFkUvmCnW3qK6QVDNG4JJR7Ja0L+yeMfgYAGIhxSShbJV1oZtdIOlXSXuon3WpKJem38Kbl5tMUEkdipaOymtJL+q08poOYVp566qlZx8U2xjbFeookPfHEE3Nup6PBYn0ltiG9Ful8mBKMBsM4G0SHYmaflvRGSWvMbLekD0o6WJLc/TJJ2zQ9ZHiXpocNv6OflgIAmgyiQ3H38+Z53SW9p6PmLCuLvUefW0SyafFGaXYqSUdKxf1Yl0iXpY9tjEkhTSixvhLbl9ZQYmKJCSUdXRYTS1M9Jf2sUul7mkbAMRoMQzQuNRQAwMDRoQAAqhjELS8MQ5uib+kSLemw4dwTFmMxOy6bErfT98XPTYcXx1tWaYE9evzxx2e2c0X+uJ8bXhx//9LbX7mFJ3MTQBlejCEgoQAAqiChYE65om/uuKaEkn67jkX6NKHEQnxMA3FhR0l65JFH5txOBwDEYnlMIWnyiL9LU4Feak48aTLKFexLxcSSm1AaMbwYfSGhAACqIKFgwdoMKc49pCpNKLEWsXr16pntNFEcddRRM9txCHBMK+lnxXbEtJKeP54jHV5cWkOJ+22GEC8ED/PCEJBQAABVkFBQpM0oolxCidKRTTGxxLpETCuSdMwxx8xsx7SRLqkSE0Y8d1rXiLWS+Llp4mlaliX3SOHciK/a6YXFJtEXEgoAoAoSChaszSii3LfmNKHEczbVU6TZ6SCmlTShxP2YZNLFG/fu3TuzHeehpLWWOMorN8elaZRXmoyaHg1QA6PB0CUSCgCgCjoUAEAV3PLCotSeACk1DylOh+XGiY7xGSXp6sBNzzlJb2XF4nu8TZYOV25aliW95RXbm5vYWHpbKndrrO3q0E3vp2CPNkgoAIAqSChYMjWeqRKXUUmfmxITQByyG9OKNDuVxHSRLqnS9Lz5+B5pdiqJn5V+bu5pjlH8/dNr0fRaelzp4p1N517I+4AmJBQAQBUkFFS12AmQOemij3HYb0wDaVI47rjjZrZjukhrLfG13CKScT9upxMbS+tJceHM9FrEelJ8LVdDKR2GXPoa9RSUIqEAAKogoWDJtJ0417RMe24CZEwouUUam9KFNHuUVtxOj4sJKFe7iaPGckvPx+SV/o5xPyaZ9HEATemlNMnk6lhMgEQpEgoAoAoSCjqxkHkOpTWVJumcj6bRYOm8kbgf00buuHi+uHRL+rlx1FjucchxW5qdXnLHNSWZxV7L9BzMV0EOCQUAUAUJBb0ovS9fep8/fkNPF31sWjo+V2tp2k73czPgYw0lviddvHLPnj0z2+lcm6ZUkj6ULO6X1lpKR43l5tAwGgwRCQUAUAUdCgCgCm55oXdtC71NBef0Fk28BdY01Djdj8X29JZX7tkmUWx7XOYl96TIePtLkg4//PCZ7Xg7LL2tF1+Lt8bSyaClt7ya8HwV5AwioZjZ6Wb2PTPbZWbvn+P1C8xsj5ntGP17Vx/tBAA06z2hmNkKSZdKerOk3ZJuNrOt7n5ncuhn3P3CzhuIXrVJKGkhOmparkVqTi/pt+2mhJJbNuWBBx6Y2U4Xm4znT1976KGHZrbjcv0xuUizf6+4nSaU0uHFbZZyaXo/lo/ihGJmbzCzG81sp5n9q5mdUqkNp0ja5e53u/vTkq6RtKnSuQEAHVlIQrlS0rsl7ZB0sqSPm9nH3f2zi2zD8ZLuCfu7JZ06x3Fnm9kbJH1f0p+6+z3pAWa2WdJmSVq/fv0im4W+dLXAZPotP1dfaTouSpNR3I+fG1OH9LNLu0SxvhLft2rVqlnHxd8lbqe1llhfaUorUvlSLqV1GIYULw8LqaE84O7Xufsed/+KpN+QdPEStSv1RUkb3P2Vkq6TdNVcB7n7FnefcveptWvXdtQ0AIBUkFDM7GpJt0r6ppldLOlv3H2/pKck7cu+ucy9ktaF/RNGP5vh7g+G3cslfbTC52IMtBk5VGO5kVKxDbmEEpNBmhoefvjhme00rcTfJb6WTo6Mkyhjklm9evWs42J6aUoradtZDh+lSv7Lu0LSs5JeoOnaxi4zu17S/0i6oUIbbpa00cxeYmaHSDpX0tZ4gJkdG3bPknRXhc8FAFQ0b0Jx929I+sZz+2a2UtKJkl4l6ZcX2wB3329mF0q6VtIKSVe6+04z+5Ck7e6+VdJ7zewsSfslPSTpgsV+LsZP7fkqteUenBUTQLqwY6yHpAtMxkUq4/nTJBMfCNaUVtLPapq7kra3dLHJ3DI5TccxX2WyLHjY8Oh21x2jf/9coxHuvk3StuRnF4ftD0j6QI3PAgAsjUFMbAQAjL/eJzYCbS3l8OI2bchNbMzd8oorBR9xxBGzXovPUYm/YzoAIC4VEydHpsX7n/zkJzPbsWCfrl4chzkzARKlSCgAgCpIKJgIQxhenJv01/TkRWl2OkgnW8YieizEp79j/Ky4mOW+fbNH9sfEEtNP+rmxTbln3jMBEhEJBQBQBQkFE6ev4cW5px42DSGWZqeBdNJjTA4xXaTPuY/iZ6VL78f3xaHGcYmXtB2xfWnb4/5SPgFyrnNieEgoAIAqSCiYeF2NBsvVDXLPec8llLgfJz2miSJOgIyfmyaK+PvHtJLWWuJ+3E7rP02fm0sobUeAUV8ZPhIKAKAKEgqWlaUcDVb6rTw3XyVNAHF5lDjiKy6vIs2eexLrJulnxXQUl+iPSSPdj0kmrd2UzleJ+7mHkpE8xhsJBQBQBQkFy1aXo8FK52ikNY+mB4LFGe/S7BFgseaRjvJqmr2fPlAsvi/Ows8lmTYP7EpTYpsRYKSa4SChAACqoEMBAFTBLS9gpPbw4tKhsrkhxU0F+3SplFikjxMW0yJ6/L2abn+lmm5/pftxu+0TIHO3vCImQA4TCQUAUAUJBZhDl4tN5oYUlz71MU6OjMvSpxMWYxE9/l6x+J9+bpQW7+N+TDLpYICmocLp+dpMgMRwkFAAAFWQUIB5LPXw4tIhxbkFJmNiiWklffZ8rKmk6SDK1XWieI6mtJI7Lndt29RT0veRcrpFQgEAVEFCARZosaPBFrJwYlN9JU0NcT/WQ3IJJdZTcrWM3ITFpt85TRfx/LkUknuNtDF8JBQAQBUkFGARasyHyNVaYjrIJZmmhJKrtTTNIZGak1d6vjYJJbc4ZGl6yWGJlv6QUAAAVdChAACq4JYXUEnb4aulr+UmQJYW72MhPnfLq6mInhvWnBt40GQhRfnF3r6iyL/0SCgAgCoGkVDM7HRJn5C0QtLl7v7h5PVDJV0t6WRJD0r6XXf/UdftBBaixmKTpYXoWIhvSg3S7MSSW/QxvpabANmUjHK/Y9thwzUmPWJp9Z5QzGyFpEslnSHpJEnnmdlJyWHvlPSwu/+8pL+X9JFuWwkAmM8QEsopkna5+92SZGbXSNok6c5wzCZJfzXa/pykS8zMvPRrCtCz0m/sbSdHNqWIXM0jppC01tKUUEr/kyt9YmNbNYYDM6S4vt4TiqTjJd0T9nePfjbnMe6+X9JeScekJzKzzWa23cy279mzZ4maCwCYyxASSjXuvkXSFkmampoivWAs1J4cWbrYZOmy+bFNaRJqSiyln7uQZNCUKHLXjBTSrSEklHslrQv7J4x+NucxZrZS0vM1XZwHAAzEEBLKzZI2mtlLNN1xnCvp95Jjtko6X9J/SzpH0teon2AS1ZjLkvtPo+m1NKE01U1KE0qufUudFNqkEpJMHb13KO6+38wulHStpocNX+nuO83sQ5K2u/tWSVdI+icz2yXpIU13OgCAAem9Q5Ekd98maVvys4vD9j5Jb+u6XUDfFlsPyM1riefILV8f37OQxRxL5Jbrr42Z8ktvCDUUAMAEoEMBAFQxiFteAOZXWngvvZWTK5S3ecJiyc9zbZhrH+OFhAIAqIKEAoyhGku55MTie5sFG2skjRrpZbFDiNt+7nJFQgEAVEFCASZMjaVcms7Rdkn5GkOKMXwkFABAFSQUYILllnIpPW4pl3xpO1qN9DJMJBQAQBUkFGAZafpm39fIpqUYDVb7/CwcWY6EAgCogg4FAFAFt7wAtC7eAxEJBQBQBQkFwM8oLd6XvGchmAA53kgoAIAqSCgAipXWWuZ7X9P7SRvjjYQCAKiChAKgtTa1lpL3L+acbc+PxSOhAACqIKEAqK70AWC1zolhIKEAAKogoQDo1FKkFwwDCQUAUAUdCgCgCm55ARiMId4OYzBAuV4Tipm9wMyuM7MfjP736IbjnjGzHaN/W7tuJwBgfn3f8nq/pK+6+0ZJXx3tz+VJd3/V6N9Z3TUPwFCY2YL/1Tg3yvXdoWySdNVo+ypJb+2xLQCARei7Q3mhu9832v4/SS9sOO4wM9tuZjeZGZ0OgCK1kwzylrwob2bXS3rRHC9dFHfc3c2sqer2Yne/18xeKulrZna7u/9wjs/aLGmzJK1fv36RLQcALMSSdyjuflrTa2b2YzM71t3vM7NjJd3fcI57R/97t5l9XdKrJf1Mh+LuWyRtkaSpqSlmSAFAh/q+5bVV0vmj7fMlfSE9wMyONrNDR9trJL1e0p2dtRAAUKTvDuXDkt5sZj+QdNpoX2Y2ZWaXj475JUnbzew2STdI+rC706EAwMD0OrHR3R+U9KY5fr5d0rtG29+S9IqOmwYAWKC+EwoAYELQoQAAqqBDAQBUQYcCAKiCDgUAUAUdCgCgCjoUAEAVdCgAgCroUAAAVdChAACqoEMBAFRBhwIAqIIOBQBQBR0KAKAKOhQAQBV0KACAKuhQAABV0KEAAKqgQwEAVEGHAgCogg4FAFAFHQoAoAo6FABAFXQoAIAq6FAAAFXQoQAAqqBDAQBUQYcCAKii1w7FzN5mZjvN7Fkzm8ocd7qZfc/MdpnZ+7tsIwCgTN8J5Q5Jvy3pxqYDzGyFpEslnSHpJEnnmdlJ3TQPAFBqZZ8f7u53SZKZ5Q47RdIud797dOw1kjZJunPJGwgAKNZrh1LoeEn3hP3dkk6d60Az2yxp82j3KTO7Y4nbNi7WSHqg70YMBNfiAK7FAVyLA36x7RuXvEMxs+slvWiOly5y9y/U/Cx33yJpy+hzt7t7Y11mOeFaHMC1OIBrcQDX4gAz2972vUveobj7aYs8xb2S1oX9E0Y/AwAMSN9F+RI3S9poZi8xs0MknStpa89tAgAk+h42/FtmtlvS6yR9ycyuHf38ODPbJknuvl/ShZKulXSXpM+6+86C029ZomaPI67FAVyLA7gWB3AtDmh9LczdazYEALBMjcMtLwDAGKBDAQBUMfYdynzLspjZoWb2mdHr3zazDd23shsF1+INZnarme03s3P6aGNXCq7Fn5nZnWb2XTP7qpm9uI92dqHgWvyhmd1uZjvM7JuTvBJF6TJOZna2mXluSahxV/B3cYGZ7Rn9Xewws3fNe1J3H9t/klZI+qGkl0o6RNJtkk5KjvkjSZeNts+V9Jm+293jtdgg6ZWSrpZ0Tt9t7vla/Jqkw0fb717mfxdHhu2zJH2l73b3dS1Gxx2h6eWgbpI01Xe7e/y7uEDSJQs577gnlJllWdz9aUnPLcsSbZJ01Wj7c5LeZPOs9TKm5r0W7v4jd/+upGf7aGCHSq7FDe7+xGj3Jk3Pb5pEJdfikbC7StKkjtQp+f8LSfprSR+RtK/LxnWs9FosyLh3KHMty3J80zE+PQR5r6RjOmldt0quxXKx0GvxTklfXtIW9afoWpjZe8zsh5I+Kum9HbWta/NeCzN7jaR17v6lLhvWg9L/Rs4e3Rb+nJmtm+P1Wca9QwEWxcx+X9KUpL/ruy19cvdL3f3nJP25pL/suz19MLODJH1M0vv6bstAfFHSBnd/paTrdOBOT6Nx71BKlmWZOcbMVkp6vqQHO2ldt1ii5oCia2Fmp0m6SNJZ7v5UR23r2kL/Lq6R9NYlbVF/5rsWR0h6uaSvm9mPJL1W0tYJLczP+3fh7g+G/y4ul3TyfCcd9w6lZFmWrZLOH22fI+lrPqo4TRiWqDlg3mthZq+W9I+a7kzu76GNXSm5FhvD7m9K+kGH7etS9lq4+153X+PuG9x9g6Zra2e5e+vFEges5O/i2LB7lqZXKsnre7RBhdEKZ0r6vqZHLFw0+tmHNP2HIEmHSfo3SbskfUfSS/tuc4/X4lc0fa/0cU2ntJ19t7nHa3G9pB9L2jH6t7XvNvd4LT4haefoOtwg6WV9t7mva5Ec+3VN6Civwr+Lvx39Xdw2+rs4cb5zsvQKAKCKcb/lBQAYCDoUAEAVdCgAgCroUAAAVdChAACqoEMBAFRBhwIAqIIOBeiQmZ1jZjeZ2W2jZ4+s7btNQC1MbAQ6ZGbHuPuDo+0PSnrA3S/tuVlAFSQUoFsXmNl3zOw2TT/8bZKfuYFlZmXfDQCWCzN7u6YfbPTr7v6Ymd2o6bWSgIlAQgG68wpJ3xp1JmdL+lVJt/fcJqAaaihAR8zsZZI+r+mnhv6npN9x91/ot1VAPXQoAIAquOUFAKiCDgUAUAUdCgCgCjoUAEAVdCgAgCroUAAAVdChAACq+H8mMJWxE5cXXAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Determine the chi-squared values over a grid of model parameters a,b\n", "# Set the range of parameters a,b and number of grid points\n", "amin,amax,na = 0.,0.5,100\n", "bmin,bmax,nb = -1.,2.,100\n", "da,db = (amax-amin)/na,(bmax-bmin)/nb\n", "alst = np.linspace(amin+0.5*da,amax-0.5*da,na)\n", "blst = np.linspace(bmin+0.5*db,bmax-0.5*db,nb)\n", "# Evaluate the chi-squared value at each grid point (i.e. for each set of model parameters)\n", "chisq2 = np.empty((na,nb))\n", "for i in range(na):\n", " for j in range(nb):\n", " chisq2[i,j] = chisqstraightline2([alst[i],blst[j]],x,y,yerr)\n", "# Find the minimum value of chi-squared across the grid\n", "minchisq = np.amin(chisq2)\n", "# Convert the chi-squared values into a probability, using P \\propto exp(-chisq/2)\n", "# We subtract the minimum chi-squared value in the exponent, to avoid numerical issues of exp(-large number)\n", "prob2 = np.exp(-0.5*(chisq2-minchisq))\n", "# Normalise the probability grid such that \\sum P(a,b) = 1\n", "prob2 /= np.sum(prob2)\n", "# Plot the probability grid as a greyscale\n", "plt.imshow(prob2,aspect='auto',origin='lower',interpolation='nearest',clim=(0.,np.amax(prob2)),cmap='Greys',extent=[amin,amax,bmin,bmax])\n", "plt.xlabel(r'$a$')\n", "plt.ylabel(r'$b$')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then we marginalize to obtain the posterior probability distributions for each parameter, $P(a)$ and $P(b)$:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deZhU9ZX4//fpvYFmX2QT6BUa2TdRRDRqEBVNXKPOT0kMMYmjScxk4kQddX7JJCYmwyTqiHEXo0GjYkJcYlBUQGSHpml6YV9kkU22prvO94+qwqbt5XZ13bq3qs7reerprqpb956CPnXqfu5nEVXFGGOM8ZsUrwMwxhhjGmIFyhhjjC9ZgTLGGONLVqCMMcb4khUoY4wxvpTmdQB1de3aVfv37+91GCbJLV26dI+qdvM6DjdYjhk/cJpjvipQ/fv3Z8mSJV6HYZKciGzyOga3WI4ZP3CaY9bEZ4wxxpesQBljjPElK1DGGGN8yQqUMcYYX7ICZYwxxpesQBljjPElK1DGGGN8yQqUMcYYX7ICZYwxxpesQBljjPElK1DGGGN8yQqUMcYYX7ICZYwxxpesQBljjPElK1DGGGN8yQqUMcYYX7ICZYwxxpesQBljjPElK1BJ4MSJE+zcuZNAIOB1KMbEnV27dnH48GGvw0hKVqAS3K9+9Suys7Pp2bMnZ599Nvv37/c6JGPigqpy++2306NHD9q1a8fVV19NbW2t12ElFStQCWzWrFn89Kc/5eKLL+aBBx5g6dKlXHjhhRw6dMjr0Izxvfvvv5/f//73TJs2jdtuu42XX36Ze++91+uwkkqa1wEYd2zatIlvfvObTJo0iVdeeYWMjAyGDRvG5Zdfzv/8z/9wzz33eB2iMb61Zs0a7r//fm666SaeeOIJAI4dO8YvfvELLrjgAs477zyPI0wOdgaVoGbMmEEgEODZZ58lIyMDgKlTp3LZZZfxu9/9joMHD3ocoTH+9dBDD9GmTRseeughRAQR4fe//z3du3fnt7/9rdfhJQ0rUAlo//79PP7441x77bX07dv3lOfuuece9u3bxyOPPOJRdMb4244dO5g1axbTpk2jS5cuJx/Pyspi+vTp/O1vf2PDhg0eRpg8XCtQIlIkIivq3A6KyA/cOp75wsyZM/n888+58847v/TcmDFjmDx5MjNmzLALvsY04JFHHqGmpoYf/ODLH1ff+c53SElJ4dFHH/UgsuTjWoFS1TJVHa6qw4FRwBHgVbeOZ77w9NNPM3HiREaMGNHg89OmTWPnzp189NFHMY7MGH9TVV588UUuvPBC8vPzv/R8nz59uOKKK3jyySepqanxIMLkEqsmvq8Alaq6KUbHS1plZWWUlpZy9dVXN7rNlClTyMrKYvbs2TGMzLSEiDwpIrtEZE0jz98gIqtEZLWILBCRYbGOMRGVlZVRUVHB5Zdf3ug21157LXv37mXRokUxjCw5xapAXQf8qaEnRGS6iCwRkSW7d++OUTiJ6/XXXwdoMsHatWvHlClTeOWVV2zwrn89DUxu4vkNwLmqOgT4L2BmLIJKdG+88QYAl112WaPbXHjhhaSmpjJ37txYhZW0XC9QIpIBTAUa/LquqjNVdbSqju7WrZvb4SS81157jVGjRn2pc0R9V111FTt27GDhwoUxisy0hKrOBz5r4vkFqrovdHcR0CcmgSW4OXPmMGLEiCbzp2PHjkyYMMEKVAzE4gzqYmCZqn4ag2MltR07drBo0SKuuOKKZre99NJLycjIOHnGZeLat4C/N/aktVI4s3v3bhYsWNDk2VPYlClTWLlyJdu2bYtBZMkrFgXqGzTSvGei680330RVmTp1arPb5uTkMG7cOObNmxeDyIxbROQ8ggXq3xvbxlopnHn33XcJBAJccsklzW47ZcoUADuLcpmrBUpE2gIXAn9x8zgmaP78+XTt2pUhQ4Y42v68885j2bJlHDhwwOXIjBtEZCjwR+ByVd3rdTzx7qOPPqJt27aMHDmy2W0HDx5M7969effdd2MQWfJytUCp6mFV7aKq9gkYA++//z4TJ05ERBxtP2nSJAKBAB988IHLkZloE5HTCX7x+xdVXe91PIngww8/5MwzzyQtrfkZ4ESEs88+267huqzZAiUi/yoinWIRjIncli1b2LBhAxMnTnT8mvHjx5OZmWnNfC6KNH9E5E/AQqBIRLaKyLdE5FYRuTW0yb1AF+CR0ED4JVEMO+kcPHiQVatWMWHCBMevOeuss9i8eTNbt251MbLk5mSy2B7AJyKyDHgSeEtV1d2wTEuFz4JaUqCysrIYP3487733nktRGSLMH1X9RjPP3wLcEp0QzaJFiwgEAi0uUAALFy5sctyhiVyzZ1CqejdQADwB3AyUi8gvRCTP5dhMC8yfP5/27dszdOjQFr1u0qRJLF++3NaJconlT3z48MMPSUlJYdy4cY5fM3z4cLKzs1mwYIGLkSU3R9egQt/4doZuNUAn4GURedDF2EwLzJ8/nwkTJpCamtqi102YMAFV5ZNPPnEpMmP5438ffvghw4cPJycnx/Fr0tPTGTNmjBUoFzm5BnWHiCwFHgQ+Aoao6ncJzq93pcvxGQf2799PaWkpZ599dotfO2rUKAAWL14c7bAMlj/xIBAI8MknnzB+/PgWv/ass85i2bJlHD161IXIjJNrUJ2Br9efR09VAyJyqTthmZZYtmwZEJypvKU6duxIYWGhnUG5x/LH58rLy/n8889PfllriTPPPJOamhqWL19+8pqUiR4nTXy59ZNLRJ4DUNVSV6IyLbJkSbADl5PxGw0ZM2aMFSj3WP743PLly4HI8if8mhUrVkQ1JhPkpEANrntHRFIJNk8Yn1iyZAkDBgw4ZXG1lhgzZgzbt2+3aVvcYfnjc8uWLSMjI4Pi4uIWv7ZPnz506dLlZJEz0dVogRKRu0TkEDA0tNjgwdD9XYBN4OYjS5YsYfTo0RG/fuzYsQB2FhVFlj/xY9myZQwZMoT09PQWv1ZEGD58uJ1BuaTRAqWq/62qOcCvVbV96JYTmhnirhjGaJqwd+9eNmzY0KoCNXz4cNLS0qxARZHlT3xQVZYvX97o4p5ODB8+nNWrV3PixIkoRmagiU4SIjJQVdcBs0XkS42zqrrM1ciMI0uXLgVoVYHKzs7mjDPOsAIVRZY/8WHz5s189tlnEV+/BRgxYgTHjx+nrKyMM844I4rRmaZ68d0JfBt4qIHnFDjflYhMi7S2g0TYiBEjbGbm6LL8iQPha0etPYOCYEcJK1DR1WiBUtVvh36eF7twTEutWLGC3NxcOnbs2Kr9DB06lKeeeopPP/2UHj16RCm65GX5Ex+WL19OSkpKi2dgqauoqIisrCyWL1/OjTfeGMXoTFNNfF9v6oWqakto+MCqVatalVxh4X2sWrWKCy+8sNX7S3aWP/FhzZo15OXl0aZNm4j3kZaWxpAhQ6wnnwuaauJrallJxdZ48tzRo0cpLy/n2muvbfW+rEBFneVPHCgpKYlKs9zQoUOZM2dOFCIydTXVxDctloGYllu7di2BQMDxAoVN6dq1K7169WLlypVRiMxY/vjfsWPHqKioiMpM5IMHD+aJJ55g9+7d2KrF0dNUE9+Nqvq8iPyooedV9bfuhWWcWL16NUBUChQEvwWuWrUqKvtKdpY//ldWVkZtbS2DBw9ufuNmhAf5rl27lnPPPbfV+zNBTc0k0Tb0M6eRm/HY6tWrycrKIj8/Pyr7Gzp0KGvXrrXxHNFh+eNzJSUlAFEvUCZ6mmrieyz08/7YhWNaYtWqVQwePLjFS2w0ZtiwYZw4cYJ169ZF7awsWVn++F9JSQmpqakUFha2el99+vQhJyfnZNEz0eFkuY1cEXlDRHaLyC4ReV1EcmMRnGna6tWro1pI6naUMNFh+eNfJSUlFBYWkpmZ2ep9iQjFxcVWoKLMyWSxLwB/BnoCvYDZwJ+c7FxEOorIyyKyTkRKRaTlC66YBu3atYtPP/00Kl3MwwoLC0lNTbVmiuiKOH+Mu9asWROV5r2wwYMHW+5EmZMC1UZVn1PVmtDteSDL4f5nAG+q6kBgGGDLC0TJmjVrAKI6cj0jI4OCggJKS+2/KYpakz/GJUeOHKGqqiqqBaq4uJhdu3axZ8+eqO0z2TU1m3lnEekM/F1Efioi/UWkn4j8BGh2ThwR6QBMBJ4AUNVqVd0frcCTXbiIRLJEQFOKi4vtW2AUtDZ/jLvKyspQ1ajmT7jYWf5ET1MDdZcSHFAoofvfqfOcAs3NyDwA2A08JSLDQvu7Q1UP191IRKYD0wFOP/1055EnudLSUnJycujVq1dU9zto0CBef/11jh8/HpW2+STW2vwxLiorKwNg4MCBUdtnuNiVlJQwceLEqO03mTXVi29AFPY9EvhXVf1YRGYAPwXuqXecmcBMgNGjR2srj5k0SktLGTRoECLS/MYtUFxcTG1tLeXl5TbxZStEIX+Mi8rKyhARCgoKorbPPn360KZNm5PFz7ReU2dQJ4nIGUAxddrOVfXZZl62Fdiqqh+H7r9MsECZKCgtLeWiiy6K+n7D3wJLS0utQEVJhPljXFRWVsbpp59OdnZ21PaZkpJCYWGhFagoarZAich/ApMIJthc4GLgQ6DJBFPVnSKyRUSKVLUM+ApgjbNRcODAAXbs2BH1608QnJlZRKwdPUoizR/jrrKyMoqKiqK+36KiIltXLYqc9OK7imBx2RmaX2wY0MHh/v8VmCUiq4DhwC8iitKcItxBYtCgQVHfd3Z2Nrm5uVagoifi/BGRJ0Njp9Y08ryIyP+KSIWIrGpoYUTzZapKWVlZVK8/hRUWFrJx40aOHz8e9X0nIycF6qiqBoAaEWkP7AL6Otm5qq5Q1dGqOlRVr1DVfa0J1gSFi4cbBSq8XytQURNx/gBPA5ObeP5ioCB0mw482oo4k8a2bds4fPiwa2dQgUCAioqKqO87GTkpUEtEpCPwOMGeScuAha5GZZpUWlpKZmYmAwa4cx2+uLiY9evXU1NT48r+k0zE+aOq84HPmtjkcuBZDVoEdBSRnq0NONGFrxG5VaAA1q9fH/V9J6Nmr0Gp6vdCv/6fiLwJtFdVmwvHQ6WlpSdnfXBDUVER1dXVbNq0iby8PFeOkSxczp/ewJY697eGHttRdyMbynEqNwtUeF4/6ygRHU7OoBCRr4vIbwleU7JPLI+tW7fOteY9+CJxLcmiw+v8UdWZoab20bZWUfDvum3btvTu3Tvq+27fvj2nnXaa5U6UOJks9hHgVmA1sAb4jog87HZgpmHV1dVs2LDBlW9/YVagosfl/NnGqdez+oQeM00oKyujsLAw6mMIw4qKiqyJL0qcjIM6HxikqgogIs8ANmWvRyorKwkEAlFZIqAxXbp0oVOnTpZk0eFm/swBbhORF4FxwAFV3dHMa5JeeXk5Y8aMcW3/RUVFvPLKK67tP5k4aeKrAOo2XPcNPWY8EC4abp5BiQhFRUV2BhUdEeePiPyJYIeKIhHZKiLfEpFbReTW0CZzgarQ/h4HvtfIrkxIdXU1GzdujOoMEvUVFRWxd+9e9u7d69oxkkVTS76/QXDOsBygVEQWh54aCyxu7HXGXeEC5WaCQTDJ3nnnHVePkciikT+q+o1mnlfg+62JM9ls3LiRQCDgav6EV7iurKykS5curh0nGTTVxPebmEVhHCsrK6N79+507NjR1eMUFRXxzDPPcOjQIXJybIXyCFj++FB5eTnwRRFxQ7j4lZeXM3bsWNeOkwyamiz2/fDvItIDCDfaLlbVXW4HZhq2fv16V5v3wsLXuMrLyxk50iYoaCnLH38KD6B18wxqwIABiIgN1o0CJ734riHYJHE1cA3wsYhc5XZgpmHhHkhus5580WH54y/l5eW0b9+erl27unaMrKws+vbtawUqCpz04vsZMCb8rU9EugH/IDg7uYmh/fv3s2vXrpgUqPz8fETEClTrWf74SEVFBQUFBa51MQ/Lz88/2ZxoIuekF19KvSaJvQ5fZ6Is/Acfiya+rKws+vXrZ13NW8/yx0fKy8td72AEwSZEO4NqPSeJ8qaIvCUiN4vIzcDfsCWrPRE+m4nFGRQEk8y+Bbaa5Y9PhLuYu9lBIiw/P5+9e/eyb5/Nj90aTRYoCZ4H/y/wGDA0dJupqv8eg9hMPeXl5aSkpJCbmxuT44ULVGiMqWkhyx9/iUUX87C6Xc1N5Jq8BqWqKiJzVXUI8JcYxWQaUV5ezumnn05mZmZMjldQUMCBAwfYs2cPNodby1n++EssupiH1e1qPnr0aNePl6icNPEtExH35gUxjsWq/TysbpKZiFn++ET4mlAsClS4lcOuQ7WOkwI1DlgkIpWhVTtXh1bINTGkqlag4pPlj09UVlbSrl27mLQGZGdn06dPH8udVnLSzfyrrkdhmrV3714OHDgQ0wI1YMAAUlNT7Vtg61j++ERlZSV5eXmudzEPy8/Pt2tQrdTsGZSqbgK6EFy9cyrQJfSYiaFYtp+Hpaen079/f/sW2AqWP/4RLlCxkpubS1VVVcyOl4iczCRxL/AMwSTrCjwlIne7HZg5VSymaGmIdTVvHcsff6itrWXDhg0xLVB5eXns3LmTw4cPx+yYicZJE98NwDBVPQYgIr8EVgD/f3MvFJGNwCGgFqhRVevOEqFwF/MBAwbE9Lj5+fl89NFHqGrMmkYSTMT5Y6Jn27ZtVFdXx7xAAVRVVTFkyJCYHTeROOkksR3IqnM/k5at2nmeqg634tQ65eXl9OvXj4yMjJget6CggEOHDrFrl81vGqHW5o+JgvC1IC8KlF2HipyTM6gDQImIvENwfZsLgcUi8r8Aqnq7i/GZkFj34Aur25OvR48eMT9+ArD88QEvClS4q7kVqMg5KVCvhm5h77Vg/wq8LSIKPKaqM+tvICLTgekAp59+ev2nDcEu5hUVFZx55pkxP3a4QFVWVjJhwoSYHz8BtCZ/TJRUVlaSlpZG3759Y3bMzp0707FjRytQrdBsgVLVZ1qx/wmquk1EugPviMg6VZ1fb/8zgZkAo0ePtjl1GhDuYh7LHnxh/fr1s67mrdDK/DFRUllZSf/+/UlLc/KdPHry8vKsQLWCq7Mqq+q20M9dBL9F2vKSEYjlCPj60tPT6devnxUoE9di3cU8LC8vz7qat4JrBUpE2opITvh34CJgjVvHS2ReFqjwca1AmXilqp4WqI0bN1JTUxPzYyeC5mYzTxWR30S47x7AhyKykuCKon9T1Tcj3FdSq6ioQERi3sU8LLz4ms1q3jKtzB8TJfv27ePAgQMxWwWgrtzcXGpqatiyZUvMj50ImpvNvFZEIroyrqpVwLCIojKnqKysjOks5vXl5+dz4MABPvvsM7p06eJJDPGoNfljoifcxObVGRQEc9irL5jxzMkVw+UiMgeYDZwcEq2qtnxAjFRUVHiSXGF117axAtVilj8eCxcoL86g6g7WNS3n5BpUFsFlqs8HLgvdLnUzKHOqiooKz64/wRcFyq5DRcTyx2PhXnRenMH07t2b9PR0K1ARctLNfFosAjENCy8Y6GWBGjBgACJiBSoCrckfEZkMzABSgT+q6i/rPX86wXn+Ooa2+amq2nLy9VRVVdG9e3dycnJifuzU1FT69+9vBSpCTiaL7SMir4rIrtDtFRHpE4vgzBff/rwsUFlZWfTp08cKVAQizR8RSQUeBi4GioFviEhxvc3uBv6sqiOA64BHoh1/IqiqqvKkeS8sNzfXxkJFyEkT31PAHKBX6PZG6DETA153MQ8L9+QzLRZp/owFKlS1SlWrgRcJLtlRlwLtQ793IDjvn6nH6wJlY6Ei56RAdVPVp1S1JnR7GnB/SUoDfFGgvEwwsMXXWiHS/OkN1O2bvDX0WF33ATeKyFZgLvCvUYg3oVRXV7N582ZPOxnl5uayf/9+9u3b51kM8cpJgdorIjeGxnSkisiNBC/6mhioqKigZ8+etG3b1tM48vLy2L17NwcPHvQ0jjjkZv58A3haVfsAU4DnRORLOS0i00VkiYgs2b17d5QOHR82b95MIBDwvIkPrCdfJJwUqG8C1wA7gR3AVYB1nIgRr0bA11e3q7lpkUjzZxtQd2bTPnx5mY5vAX8GUNWFBHsMdq2/I1WdqaqjVXV0t27J1fjhZRfzMCtQkXO05LuqTlXVbqraXVWvUNXNsQjOeN/FPMzWtolMK/LnE6BARAaISAbBThBz6m2zGfgKgIgMIligkusUqRl+KlCWOy3XaDdzEfmJqj4oIr8neDH2FLaOjfuOHDnC9u3bfVWgrCefM63NH1WtEZHbgLcIdiF/UlVLROQBYImqzgHuBB4XkR+GjnGz2nxUp6isrCQzM5NevXp5FkNOTg7dunWzM6gINDUOqjT0c0ksAjFf5uUULfXl5OTQvXt3+xboXKvzJzSmaW69x+6t8/ta4OxI958Mqqqq6N+/Pykpri7c0Kzc3FwrUBFotECp6huhsRhDVPXHMYzJhPhhDFRdNqu5c5Y//lBVVeWLL3i5ubksXLjQ6zDiTpNfK1S1FvuG5plwMfBDgoEtvtZSlj/eUlXPx0CF5ebmsmXLFk6cOOF1KHHFyWSxK2yyS29UVlbSuXNnOnXq5HUoQPAM6vnnn+fYsWNkZWV5HU68sPzxyGeffcbBgwd9U6Bqa2s9H5MVb5wUqLqTXYYpYAnmMq9nMa8vLy8PVWXDhg0MGjTI63DiheWPR/zQgy+s7qzmfsppv7PJYn2soqKCM8880+swTqo7q7kVKGcsf7wTbo72Q0GwsVCRcTJZbKGIvCsia0L3h4rI3e6Hltyqq6vZtGmTL5IrzJbdaDnLH++Ei4EfFgrs1asXGRkZVqBayEnfy8eBu4ATAKq6iuCgQeOiTZs2EQgEfNODD6Bz58506NDBOkq0jOWPR6qqqujRo4fn04SBLbsRKScFqo2qLq73WI0bwZgv+Kl5IkxErKt5y1n+eMQvPfjCbCxUyzkpUHtEJI/QaHgRuYrgnGLGRX5ZZqM+62reYpY/HvFbgQrnjk324ZyTAvV94DFgoIhsA34A3Or0AKEZnJeLyF8jjDEpVVZW0rZtW3r06OF1KKfIz89n48aNNp7DuVblj4lMdXU1W7Zs8VWBys3N5cCBA7bsRgs4KVCqqhcQXMNmoKpOcPi6sDv4YtoX41B4klgR8TqUU+Tn51NTU8PmzTZfsEOtzR8TAT8ss1Gf9eRrOSeJ8gqAqh5W1UOhx152svPQ0taXAH+MLLzk5ZdZzOuzSWNbLOL8MZELN0NbgYpvTc1mPhAYDHQQka/Xeao9wcGHTvwP8BMgJ+IIk1BtbS1VVVVMnTrV61C+xNaFciZK+WMi5KeJlsPC3d0td5xraqBuEXAp0BG4rM7jh4BvN7djEbkU2KWqS0VkUhPbTQemA5x++ukOQk58W7dupbq62pdnUD179iQ7O9vOoJrXqvwxrVNZWUlWVhY9e/b0OpSTbNmNlmtqNvPXgddFZHxotc6WOhuYKiJTCH5jbC8iz6vqjfWOMxOYCTB69Gjr3oJ/e/BBsKt5Xl6eFahmRCF/TCtUVlYyYMAAz5fZqM96wbaMk/+9r4lIexFJD42I3y0iNzb3IlW9S1X7qGp/ggMT/1m/OJmG+W2ZjfpsLFSLRJQ/pnX8OuedFaiWcVKgLlLVgwSbKzYC+cC/uRlUsquoqCAzM5PevXt7HUqD8vPzqaqqIhAIeB1KPLD8iTFVpbKy0rcFasuWLVRXV3sdSlxwUqDSQz8vAWar6oGWHkRV31PVS1v6umRVUVFBbm6u75onwvLy8jh+/Dhbt271OpR40Or8MS2ze/duDh8+7KsefGG5ubmoKhs3bvQ6lLjg5BPwDRFZB4wC3hWRbsAxd8NKbn7tYh5WUFAAWFdzhyx/YsyP04SFhWOyZj5nmi1QqvpT4CxgtKqeILjo2uVuB5asws0Tfi5QNqu5c5Y/sRcPBcp68jnT7HpQIpIO3AhMDM1q8D7wfy7HlbS2b9/OkSNHTp6l+FHfvn3JzMykvLzc61B8z/In9qqqqhAR+vfv73UoX3LaaafRpk0bO4NyyMmKuo8SbEd/JHT/X0KP3eJWUMks/KHv5zOolJQUcnNz7QzKGcufGKusrKR3795kZflvPLSIkJubawXKIScFaoyqDqtz/58istKtgJJd+EPfz2dQEIzPCpQjlj8x5rdZzOuzZTecc9JJoja0XAAAIpIL1LoXUnIrLy8nIyODvn37eh1Kk8JjoayrebMsf2LMr13Mw/Ly8qiqqrJlNxxwcgb1b8A8EakCBOgHTHM1qiQW7mKemprqdShNKigo4NixY2zfvp0+ffp4HY6fWf7E0OHDh9mxY4evm8jz8vI4cuQIO3fu9NVUTH7UbIFS1XdFpIDg3GIAZap63N2wkld5ebnvm/fgi2tk5eXlVqCa0Jr8EZHJwAwgFfijqv6ygW2uAe4juCDiSlW9PiqBxym/z8ICp3Y1twLVtGab+EQki+Cia/cB/wl8N/SYibJAIOD7MVBhNhbKmUjzR0RSgYeBi4Fi4BsiUlxvmwLgLuBsVR1McDHEpBYPBcqGaTjn5BrUswSXDfg98IfQ78+5GVSy2rFjB0ePHo2LM6g+ffqQkZFhXc2bF2n+jAUqVLVKVauBF/ny+KlvAw+r6j4AVd0VtajjVPhD38/XoPr160dqaqoVKAecXIM6Q1XrfnObJyJr3QoomcVDF/Ow1NRU8vLyrEA1L9L86Q1sqXN/KzCu3jaFACLyEcFmwPtU9c36O0qmJW0qKiro2rUrHTp08DqURqWnp9O/f38rUA44OYNaJiJnhu+IyDhgiXshJa/wh308nEFBME4rUM1yM3/SgAJgEvAN4HER6Vh/I1WdqaqjVXV0t27donRof/L7LCxhtiKAM04K1ChggYhsFJGNwEJgjIisFpFVrkaXZCoqKuKii3lYYWGhdTVvXqT5sw2o+4fQJ/RYXVuBOap6QlU3AOsJFqykVVFR4evmvbBwgbKu5k1z0sQ32fUoDADr168nPz/f913MwwoLCzl+/DhbtmyhX79+XofjV5HmzydAgYgMIFiYrgPq99B7jeCZ01Mi0pVgk1/SjgA9fvw4mzdvjpszqL7R8c8AACAASURBVAMHDrB37166du3qdTi+5aSb+aZYBGKCBaqwsNDrMBwLx7p+/XorUI2INH9UtUZEbgPeInh96UlVLRGRB4Alqjon9NxFoWtatcC/qereaMUebzZs2ICqxk2BgmCTpBWoxvlzwaEkVFtbS0VFRdwWKBN9qjpXVQtVNU9Vfx567N5QcUKDfqSqxao6RFVf9DZib8VDF/Mw62rujBUon9i0aRPV1dVxVaBOO+002rVrZwXK+EI8dDEPGzBgACJiBaoZVqB8IvwhX1RU1MyW/iEiFBQUWIEyvlBeXk779u3josksMzOTvn37WoFqRqPXoETkEMHpUxqkqu1diShJhT/k4+kMCoLxfvLJJ16H4TuWP7EXvoYbWnfL9/Lz822YRjMaPYNS1ZxQEs0Afkpw4GAf4N+B/4lNeMlj/fr1dOjQgXgbp1JYWMjGjRs5ftymZ6zL8if24mUey7DCwkIrUM1w0sQ3VVUfUdVDqnpQVR/FlqyOunj79hdWWFhIIBCw9W0aZ/kTA8ePH2fTpk1x1QJRUFDAZ599xt69SdvxsllOCtRhEblBRFJFJEVEbgAON/ciEckSkcUislJESkTk/taHm7jirYt5mPXka1ZE+WNaprKyElWNuzMosNxpipMCdT1wDfBp6HY1Xx4w2JDjwPmh1USHA5PrTvlivnD06FE2b94cVx0kwsJJVlZW5nEkvhVp/pgWiLdpwuCL3LFmvsY5Gai7kQiaJDQ4h8fnobvpoZvN69GA8JQn8ZRcYR07duS0005j3bp1XofiS5Hmj2mZ8FlIPOXQgAEDSE1NtTOoJjhZD6pQRN4VkTWh+0NF5G4nOw81a6wAdgHvqOrHDWwzXUSWiMiS3bt3tzT+hBD+cB80aJDHkURm4MCBVqAa0Zr8Mc6Vl5fTtWtXOnXq5HUojqWnpzNgwAArUE1w0sT3OMFF0U4AqOoqgvOCNUtVa1V1OMHeS2NF5IwGtkmamZYbs27dupNjiuJRuEDZxJcNijh/jHPl5eVxew3Xmvga56RAtVHVxfUeq2nJQVR1PzAPm3i2QaWlpfTr1482bdp4HUpEioqK2LdvH3v27PE6FD9qdf6Y5q1fvz4uv+CFB7rbl7uGOSlQe0Qkj9D1IxG5CtjR3ItEpFt4bRoRyQYuBKwdqAHr1q1j4MCBXocRsXDs1szXoIjyxzh3+PBhtm/fHpcFqrCwkCNHjrB9+3avQ/ElJwXq+8BjwEAR2Qb8APiug9f1JLh66CqCSwe8o6p/jTjSBBUIBCgrK7MClbgizR/jULzOwgLWk685TnrxVQEXiEhbIEVVDznZcaitfUQr40t4W7du5ciRI3HbQQKCy4hnZWVZgWpApPljnAsPcYjHL3nhs76ysjImTZrkbTA+5KQXX62I/BI4Ek4uEVnmemRJIvyhHo/JFZaSkkJRUZEVqAZY/rgvnjsZ9e3bl+zsbBtH2AgnTXwloe3eFpHOocfiaz4eHystLQXiu0BBMH5LsgZZ/rhs3bp1DBgwgKysLK9DaTH7ctc0JwWqRlV/AvwR+EBERmEDbqNm3bp1dOrUKe4mia1v4MCBVFVVcfToUa9D8RvLH5eVlZXF5SwsYTaOsHFOCpQAqOpLwLXAU0Cum0Elk3Xr1lFUVBR3k8TWV1xcjKraWdSXWf64KFE6GW3cuNG+3DXASYG6JfyLqq4BzgFudy2iJLN27VqKi4u9DqPVBg8eDATfjzmF5Y+Ltm7dytGjR+P+DEpVrSdfA5pasPB8Vf0n0E9E+tV7+vOGXmNaZs+ePezatevkh3s8KygoIDU11QpUiOVPbCRCJ6O6wzSGDh3qcTT+0lQ383OBfwKXNfCcAn9xJaIkUlJSApAQBSojI4OCgoKT78lY/sRCIhSogoICRMSuQzWg0QKlqv8Z+jktduEkl0QqUBB8H6tXr/Y6DF+w/ImNsrIyOnToQPfu3b0OJWJt2rShX79+VqAa0FQT34+aeqGq/jb64SSXtWvX0r59e3r37u11KFFRXFzMq6++yvHjx8nMzPQ6HE9Z/sRGonQysp58DWuqk0ROMzfTSiUlJRQXF8d9coUVFxcTCARs+YAgy58YKCkpSYgWiPA4wkAg4HUovtJUE58t0e6ykpISpk6d6nUYURP+oCgpKWHIkCEeR+Mtyx/37d27l08//TQhCtSgQYM4cuQImzZtYsCAAV6H4xvNzsUnIlnAt4DBwMmh2qr6TRfjSni7d+9m9+7dCZFcYYWFhaSkpFhHiTpakz8iMhmYAaQCf1TVXzay3ZXAy8AYVV0SjbjjQfjvLJGGaZSUlFiBqsPJOKjngNOArwLvE1x80Ca8bKVwd+xEKlCZmZkUFhayZs0ar0Pxk4jyR0RSgYeBi4Fi4Bsi8qVPYhHJAe4AvrRadaJLpByqW6DMF5wUqHxVvQc4rKrPAJcA49wNK/GFP8QTIbnqGjJkiPXkO1Wk+TMWqFDVKlWtBl4ELm9gu/8CfgUci1bA8aKkpIR27drRt29fr0NptY4dO9K7d28rUPU4KVAnQj/3h5Zs7wDEb59On1i1ahWdO3emV69eXocSVUOHDqWyspLPP7exqCGR5k9vYEud+1tDj50kIiOBvqr6t6Z2JCLTRWSJiCzZvXu388h9LtE6GQ0ePNgKVD1OCtRMEekE3APMAdYCD7oaVRJYtWoVQ4cOTZjkCgt3jrBEO8mV/BGRFOC3wJ3NbauqM1V1tKqOjvdJietau3ZtQrVADB48mLVr11JbW+t1KL7RbIFS1T+q6j5VfV9Vc1W1u6r+XyyCS1SBQIDVq1cn5LQm4fe0atUqjyPxh1bkzzagbttVn9BjYTnAGcB7IrIROBOYIyKjoxW7nyVSD76wwYMHc+zYMTZs2OB1KL7hpBdfJnAl0L/u9qr6gHthJbaqqioOHz6ckAWqX79+tGvXzgpUSCvy5xOgQEQGECxM1wHX13n9AaBrneO8B/w4WXrxJVIPvrC6HSXy8/M9jsYfnDTxvU7w4mwNcLjOzUQo/OGdiAUqJSXFOkqcKqL8UdUa4DbgLaAU+LOqlojIAyKSOIPnIpSInYzCxdaax7/Q7BkU0EdVJ7seSRJZtWoVKSkpCZVcdQ0ZMoTZs2ejqgl3jS0CEeePqs4F5tZ77N5Gtp0UyTHi1apVq+jYsWNC9OALa9++PaeffroN06jDyRnUAhFp8bQAItJXROaJyFoRKRGROyKILyGtWrWKgoIC2rRp43Uorhg6dCj79u1j27ZtzW+c+CLKH9O0lStXMmzYsIT7AjR06FBWrlzpdRi+4aRATQCWikiZiKwSkdUi4uQCQw1wp6oWE7yA+/2GBhomo5UrVyZk817YsGHDAFixYoXHkfhCpPljGpHInYyGDRtGWVkZx44l3bC2Bjlp4rs4kh2r6g5gR+j3QyJSSnAcR1KvaHfw4EGqqqqYNi1xV2EIf7Ndvnw5l156qdfheC2i/DGNC3cyCn8RSiTDhg2jtraWkpISRo0a5XU4nmv0DEpE2od+PdTIzTER6Q+MoIHpWBJ1EGFjwmcVI0eO9DgS9+Tk5FBQUMCyZcu8DsUz0cwfc6pE7mQULrrWzBfU1BnUC8ClwFKCK4DWbexVINfJAUSkHfAK8ANVPVj/eVWdCcwEGD16tDoLO34tXboUSOwCBcH3t2DBAq/D8FJU8sd82cqVKxO2k1FeXh5t27a1AhXS1HIbl4Z+Rjy1roikEyxOs1TVlrgGli1bRq9evTjttNO8DsVVI0eO5MUXX2Tv3r106dLF63BiLhr5YxqWyJ2MUlNTGTJkiBWoECcDdRv6qn8A2BQaq9HY6wR4Aii11UO/sHTp0oQ/e4IvzhCXL1/OBRdc4HE03ok0f0zjVq5cyZgxY7wOwzXDhg3jpZdesmEaOOvF9wiwiGAz3OOh32cDZSJyUROvOxv4F+B8EVkRuk1pbcDx7PDhw6xbty4pLn6OGDECIKmvQ4VEmj+mAfv372fDhg0J2UEibNiwYezfv58tW7Y0v3GCc1KgtgMjQpNNjgKGA1XAhTQx6aWqfqiqoqpDVXV46Da3se2TwcqVK1HVpDiD6ty5M/369WP58uVeh+K1iPLHNCz8hWf06MSdcnD48OGAfbkDZwWqUFVPzr2hqmuBgapa5V5YiSncQSIZzqAg2My3ZElSTA3XFMufKEqGHBo2bBipqakn32syc1Kg1orIoyJybuj2SOixTL5Y68Y4sHTpUrp3755wa0A1ZuzYsVRUVPDZZ595HYqXLH+iaOnSpfTr1y+hO960adOGwYMH25c7nBWom4AK4AehWxVwM8HkOs+1yBLQxx9/zNixY5Pmwue4ccGFYxcvXuxxJJ6y/ImiJUuWJPTZU9jo0aNZsmQJqgk/8qZJTRYoEUkF5qrqQ6r6tdDtN6p6RFUDqmrLpjq0f/9+1q1bx5lnnul1KDEzevRoRISPP/7S+OykYPkTXfv376eysjJpCtSePXvYvHmz16F4qskCpaq1QEBEOsQonoQVPotIpgKVk5NDcXFx0hYoy5/oSoYOEmHhIpzszXxO5uL7HFgtIu9QZx0bVb3dtagS0KJFixCRhB6/0ZBx48bx+uuvJ/OYDsufKEmGDhJhQ4cOJS0tjSVLlnDllVd6HY5nnBSov4RuphUWLVpEcXEx7du3b37jBDJu3DiefPJJqqqqyMvL8zocL1j+RMnixYvp379/QneQCMvKymLIkCF2BtXcBqr6TCwCSWSqyscff8zXvvY1r0OJubFjxwLBDiLJWKAsf6JDVVmwYAHnnnuu16HEzJgxY3jppZcIBAKkpDjpz5Z4mn3XIlIgIi+HFh6sCt9iEVyiCHe1DvdqSyZnnHEG7dq146OPPvI6FE9Y/kTHli1b2L59O+PHj/c6lJg566yzOHDgAKWlpV6H4hknZfkp4FGCCxCeBzwLPO9mUIkm/OGcTMkVlpaWxvjx45k/f77XoXjF8icKFi5cCCRXDp111lkASfvlDpwVqGxVfRcQVd2kqvcBl7gbVmKZP38+nTt3prg4ORcUPuecc1izZk2yDti1/ImChQsXkp2dndBz8NWXn59P165dk3rZGicF6riIpADlInKbiHwNaOdyXAll/vz5TJw4MWnbkc855xwgab8JWv5EwYIFCxgzZgzp6elehxIzIsJZZ51lBaoZdwBtgNuBUQRnKL/JzaASybZt26isrGTixIleh+KZcePGkZ6ezgcffOB1KF6w/Gmlo0ePsnz58qRq3gs7++yzKS8vJxlWG2+Ik158n4R+/RyY5m44iSd87SWZC1R2djZjxoxJygJl+dN6n3zyCTU1NSevySST8HtesGABl19+ucfRxF6jBUpE5jT1QlWdGv1wEs/8+fPJyclJqrbzhpxzzjk89NBDHD58mLZt23odjussf6LnvffeQ0RONhUnk1GjRpGRkcEHH3xgBaqe8cAW4E/Ax0BSTgPQWu+//z4TJkwgLc3JmOjEdd555/GrX/2KDz74gMmTJ3sdTixY/kTJvHnzGD58OJ06dfI6lJjLzs7mzDPPZN68eV6H4ommrkGdBvwHcAYwg+ACa3tU9X1VfT8WwcW7bdu2UVpaynnn2aTV55xzDhkZGbzzzjtehxIrrc4fEZksImUiUiEiP23g+R+FxletEpF3RaRfVN+BDxw7doyFCxcmdQ6df/75LF++nH379nkdSsw1WqBUtVZV31TVm4AzCS4Z8J6I3Baz6OLc22+/DcBXv/pVjyPxXps2bZgwYULSFKjW5k9oJvSHgYuBYuAbIlJ/nMJyYLSqDgVeJgFX6F20aBHHjx9P+gKlqrz/fvKdFzS33EamiHyd4MDC7wP/C7wai8ASwdtvv81pp53GkCFDvA7FFy688EJWr17Nzp07vQ4lJlqZP2OBClWtUtVq4EXglIsQqjpPVY+E7i4C+kQncv+YN28eKSkpSXn9KWzcuHFkZ2fzz3/+0+tQYq7RAiUizwILgZHA/ao6RlX/S1W3xSy6OBYIBHjnnXe46KKLknUW7y+58MILAXj33Xc9jsR9Ucif3gSvYYVtDT3WmG8Bf28klukiskRElsRbd+V58+YxcuRIOnRI3hVLMjIyOOecc6xA1XMjUEBwHMcCETkYuh0SkYPN7VhEnhSRXSKyJlrBxpPly5ezd+9eLrroIq9D8Y0RI0bQpUuXk02fCa5V+dMSInIjMBr4dUPPq+pMVR2tqqO7desWzUO76sCBAyxcuJALLrjA61A8d/7551NSUpI0rQ9hTV2DSlHVnNCtfZ1bjqo6WTPiaSApums15M033wSw5KojJSWFyZMnM3fuXGpra70Ox1VRyJ9tQN869/uEHjuFiFwA/AyYqqrHoxO9P/zjH/+gpqaGSy6xmaHC17H//vcGT5ITlmtz76jqfCApJ18DeP311xk7diw9evTwOhRfmTp1Knv27Dk5+adp1CdAgYgMEJEM4DrglLFVIjICeIxgcdrlQYyumjt3Lh07dkyqVagbM2zYMHr37s3f/vY3r0OJqeScHM5lW7du5ZNPPuGKK67wOhTfmTx5Munp6cyZ0+Q41qSnqjXAbcBbQCnwZ1UtEZEHRCQ8yPfXBOf1my0iK5obHBxPVJW5c+dy0UUXJf0YQgjOy3fJJZfw9ttvU11d7XU4MeN5gYrnC7iNCX/4WoH6svbt23Peeefx+uuvex2K76nqXFUtVNU8Vf156LF7VXVO6PcLVLWHqg4P3RJmdooVK1awc+dOpkyZ4nUovnHppZdy6NChpJoyzPMCFa8XcJvy2muvUVRUxKBBg7wOxZemTp3K+vXrWbdundehGJ964403EJFkmXXEkfPPP5/MzEz++te/eh1KzHheoBLNvn37mDdvnp09NeHyyy9HRJg9e7bXoRifmj17NhMmTLBruHW0bduWr3zlK7z22muoqtfhxIRrBUpE/kRwHEiRiGwVkW+5dSw/efnll6mpqeGqq67yOhTf6tOnDxMnTmTWrFlJk2jGudLSUtasWcPVV1/tdSi+c80117Bx40YWL17sdSgx4WYvvm+oak9VTVfVPqr6hFvH8pPnn3+egQMHMmrUKK9D8bXrr7+esrIyVqxY4XUoxmdmz56NiHDllVd6HYrvXHHFFWRkZPDiiy96HUpMWBNfFG3cuJH58+dz44032uwRzbjyyitJS0vjhRde8DoU4zPh5r1evXp5HYrvdOjQgcmTJzN79mwCgYDX4bjOClQUhT9sb7jhBo8j8b8uXbowefJkXnjhBWpqarwOx/jEypUrWbNmDddcc43XofjWtddey7Zt2/jwww+9DsV1VqCiJBAI8NRTT3HOOefQv39/r8OJC9/85jfZvn170g0+NI174oknyMzM5Prrr/c6FN+aOnUq7dq148knn/Q6FNdZgYqSt99+m4qKCr773e96HUrcuOyyy+jduzePPvqo16EYHzh69CjPPfccX//61+ncubPX4fhWu3btuOGGG3jppZcSfo0oK1BR8oc//IEePXrYhd0WSEtL49vf/jZvvfUWVVVVXodjPPbqq6+yf/9+brnlFq9D8b3p06dz7NgxZs2a5XUorrICFQVVVVXMnTuX6dOnk5GR4XU4ceWWW24hNTWVhx9+2OtQjMceffRRcnNzmTRpkteh+N7IkSMZNWoUjz32WEIP1bACFQUPPfQQaWlpfOc73/E6lLjTu3dvrrvuOh577DH27NnjdTjGIwsXLuTDDz/kjjvuICXFPpac+N73vseaNWsSepVq+0topa1bt/LHP/6RadOm0bt3U+vJmcbcddddHD58mBkzZngdivHIr3/9azp37sy3vpUU4/mj4oYbbqBXr17893//t9ehuMYKVCs9+OCDBAIB7rrrLq9DiVuDBw/ma1/7Gr///e/Zv3+/1+GYGFu3bh2vvfYa3//+92nbtq3X4cSNzMxM7rzzTt577z0WLVrkdTiusALVChs2bGDmzJncdNNN1rW8le69914OHjzIz3/+c69DMTH2s5/9jLZt23Lbbbd5HUrcmT59Op07d+bee+9NyGtRVqBa4c477yQtLY3777/f61Di3vDhw7n55puZMWMGFRUVXodjYuSjjz7iL3/5Cz/5yU/o3r271+HEnXbt2nHPPffwzjvvnFzFO5FYgYrQO++8w6uvvsrdd99t156i5Oc//zkZGRn88Ic/TMhvg+ZUgUCAH//4x/Ts2ZMf/ehHXocTt773ve9RUFDAnXfemXCzsliBisDBgweZPn06+fn5/PCHP/Q6nITRs2dP7r//fv7617/y3HPPeR2Ocdkf/vAHFi1axC9/+Uu79tQKGRkZ/OY3v6G0tJQHH3zQ63CiS1V9cxs1apTGg5tuuklTUlJ0wYIFXoeScGpqavScc87R9u3b66ZNmzyJAViiPsgHN25+ybGKigpt06aNXnzxxRoIBLwOJyFce+21mp6eritWrPA6lGY5zTE7g2qhJ554gmeeeYa7776b8ePHex1OwklNTeXpp59GVfn617/OkSNHvA7JRNnRo0e59tprSUtLY+bMmTbzf5Q8/PDDdO7cmeuvv55Dhw55HU5UWIFqgX/+85/ceuutXHTRRdxzzz1eh5OwcnNzeeGFF1i2bBnf/OY3k2JZgWShqnznO99h6dKlPP/88/Tp08frkBJGly5dmDVrFmVlZVx//fXU1tZ6HVKrWYFyaP78+Vx++eUUFRXx5z//mbS0NK9DSmiXXnopv/rVr3jppZe45ZZbrEglAFXlxz/+Mc899xwPPPAAl112mdchJZyvfOUrzJgxg7/+9a9897vfjf+8cdIOGKubX9rH65szZ45mZ2frwIEDdevWrV6Hk1Tuu+8+BfSaa67Rw4cPx+SY2DWoqKupqdHbb79dAb399tvtupPL/uM//kMBvfnmm7W6utrrcL7EaY55njB1b34rUNXV1Xr33XcroCNHjtRPP/3U65CS0oMPPqgiosOHD9e1a9e6fjwrUNG1e/duveiiixTQH/7wh1acYiAQCOj999+vgJ599tm6bds2r0M6hRWoVnrvvff0jDPOUECnTZumR44c8TqkpDZ37lzt3LmzZmRk6H333aeHDh1y7VhWoKKjtrZWn3/+ee3WrZump6fr448/HrNjm6AXXnhB27Rpox07dtTHHntMa2pqvA5JVa1ARaS6ulpfffVVnTRpkgLat29fff311z2NyXxh586des011yig3bp10wceeMCVb4ZWoFrn6NGjOmvWLB06dKgCOm7cOF29erXrxzUNW7dunZ577rkK6MCBA/Wpp56KWXN5Y3xRoIDJQBlQAfy0ue1jXaACgYBu2bJF//SnP+m0adO0W7duCmivXr30d7/7nef/iaZhCxcu1K9+9asKaEpKip577rn64IMP6uLFi/X48eOt3r9fClRz+QNkAi+Fnv8Y6N/cPt3Kse3bt+tLL72kN910k3bu3FkBLSws1Oeff94339qTWSAQ0D//+c8nvzTk5OToDTfcoLNmzdJNmzbFvNnVaY5JcNvoE5FUYD1wIbAV+AT4hqqubew1o0eP1iVLlrT4WKpKTU0NJ06c4MSJE1RXV3Ps2DGOHDnC4cOHOXDgAPv27WPPnj1s376dLVu2UFlZSUlJyck1iDp06MCUKVO47rrrmDJlivXSiwMVFRU8++yzvPrqq6xZswaA9PR0Bg0aRGFhIf3796dXr150796drl270qFDB9q1a0fbtm3Jzs4mMzOTjIwM0tLSSEtLIyUlBRFBRJaq6mgv35uT/BGR7wFDVfVWEbkO+JqqXtvUfpvLsUAgQG1t7ck8On78+Mk8Onjw4Mk82rFjB5s3b6aiooKSkhK2b98OQMeOHbnkkku4+eabOf/8821tJ59RVT744AOefvpp3njjjZOff927d+eMM84gPz+f008/nZ49e9K1a1c6depE+/btadeuHdnZ2WRlZZGRkUF6evopOdNSTnPMzQI1HrhPVb8aun8XgKo2unhJamqqZmdnE9ru5OPh3+tW1kAgcMrPFsRFjx49yMvLY+DAgQwbNozx48czfPhwK0pxbMeOHXzwwQcsX76c1atXU15ezubNmzl27FiL9nPHHXcwY8YMPxSoZvNHRN4KbbNQRNKAnUA3bSIhevTooYcOHfrSN9VwYWqJTp06kZuby+DBgxk2bBgTJkxgxIgRpKenR/COTazV1tayYsUKFixYwIoVK1izZg0bNmxg9+7dLdpPSkrKyVu4WIW+6J1SvOr+fvjwYc8L1FXAZFW9JXT/X4Bxqnpbve2mA9NDd88A1rgSkDNdAS+XdfX6+H6IwevjAxSpao6XATjJHxFZE9pma+h+ZWibPfX25accc4sf/m7ckKjvy1GOeX7KoKozgZkAIrLEy2+uyX58P8Tg9fHDMXh5/GjzU465xd5XfHGaY242EG8D+ta53yf0mDGmeU7y5+Q2oSa+DsDemERnTAy4WaA+AQpEZICIZADXAXNcPJ4xicRJ/swBbgr9fhXwz6auPxkTb1xr4lPVGhG5DXgLSAWeVNWSZl420614HEr244P3MXh9fPBBDI3lj4g8QLCL7hzgCeA5EakAPiNYxJrj+Xtzib2v+OLofbnWScIYY4xpDRukYIwxxpesQBljjPGlmBcoEZksImUiUiEiP23g+UwReSn0/Mci0t+DGCaKyDIRqQmNR4n18X8kImtFZJWIvCsi/TyI4VYRWS0iK0TkQxEpjuXx62x3pYioiES1q62D93+ziOwOvf8VInJLNI/vJRG5WkRKRCQQ7X9XLzj9W4onIvKkiOwKjXVLGCLSV0TmhT7fSkTkjiZf4GQ+pGjdCF7srQRygQxgJVBcb5vvAf8X+v064CUPYugPDAWeBa7y4PjnAW1Cv3/Xo3+D9nV+nwq8Gcvjh7bLAeYDi4DRMX7/NwN/iOa/u19uwCCgCHgvmv+uHr0XR39L8XYDJgIjgTVexxLl99UTGBn6PYfgdF6N/n/FnfiPiAAAA7VJREFU+gxqLFChqlWqWg28CFxeb5vLgWdCv78MfEUimeypFTGo6kZVXQW4sRylk+PPU9UjobuLCI6BiXUMB+vcbQtEszeNk78DgP8CfgW0bL6i6B0/IalqqaqWeR1HlCTk/6WqzifYMzOhqOoOVV0W+v0QUAr0bmz7WBeo3sCWOve38uXgTm6jqjXAAaBLjGNwU0uP/y3g717EICLfD02f8yBweyyPLyIjgb6q+rcoHtfx8UOuDDWzviwifRt43njP63w2EQpdvhlBcCb+BlknCR8TkRuB0cCvvTi+qj6sqnnAvwN3x+q4IpIC/Ba4M1bHbMAbBJevGAq8wxdn9XFBRP4hImsauMX92YWJfyLSDngF+EG91ppTxHouvpZM37LVpelbvJ6CydHxReQC4GfAuap63IsY6ngReDSGx88hOKnpe6HW3dOAOSIyVVWjMU9es+9fVev+zf2R4Flk3FDVC7yOIUa8zmfTQiKSTrA4zVLVvzS1bazPoPwwfYvXUzA1e3wRGQE8BkxV1V0exVBQ5+4lQHmsjq+qB1S1q6r2V9X+BK/DRas4NXt8ABHpWefuVIJt5cZ/vM5n0wKh/gRPAKWq+ttmX+BBL44pBHtuVAI/Cz32AMEPIIAsYDbBVUIXA7kexDCGYFv2YYJnbyUxPv4/gE+BFaHbHA/+DWYAJaHjzwMGx/L49bZ9jyj3NnPw/v879P5Xht7/wGj/H3h1A74W+vs+Hvo7e8vrmKL9fxnvN+BPwA7gROj/6ltexxSl9zWBYIerVXU+36Y0tr1NdWSMMcaXrJOEMcYYX7ICZYwxxpesQBljjPElK1DGGGN8yQqUMcYYX7ICZYwxxpesQBljjPElK1AJQkSuEpFFIrIytH5TN69jMiaRiMgkEXnO6ziSiRWoxDFPVc9U1WEEJze9xuuAjEkww4DlXgeRTKxAJY6bRWSxiKwkuOhjtNdQMibZDQd6h1b6rhKRSV4HlOisQCUAEfn/CC7cdn7oDKqM4DxyxpjoGQYcUtVxwK0EF9Q0LrIClRiGAAtU9XMRuRI4C1jtcUzGJIzQEhFdgV+EHloRum9cZAUqMTwNfE9EFhNcobJKVQ97G5IxCWUgwaXlq0P3RxKc6d64yGYzN8aYZojIvxBcjmUgkE6wI9IPVXWRp4EluFivqGuMMfFoGPAXYAGQDfyXFSf32RmUMcYYX7JrUMYYY3zJCpQxxhhfsgJljDHGl6xAGWOM8SUrUMYYY3zJCpQxxhhfsgJljDHGl/4fCZaRBREomvMAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Marginalise the 2D probability distribution to obtain the 1D posterior probabilities of the 2 parameters\n", "# Marginalise parameter b, by summing over the probabilities corresponding to all values of b\n", "proba = np.sum(prob2,axis=1)\n", "# Marginalise parameter a, by summing over the probabilities corresponding to all values of a\n", "probb = np.sum(prob2,axis=0)\n", "# Normalise the probabilities such that \\int P(a) da = 1 (currently, \\sum P(a) = 1)\n", "proba /= da\n", "probb /= db\n", "# Plot the posterior probability distributions for a and b\n", "fig = plt.figure()\n", "plt.subplot(121)\n", "ymin,ymax = 0.,1.1*np.amax(proba)\n", "plt.plot(alst,proba,color='black')\n", "plt.xlabel(r'$a$')\n", "plt.ylabel('Marginalized posterior probability')\n", "plt.xlim(amin,amax)\n", "plt.ylim(ymin,ymax)\n", "plt.subplot(122)\n", "ymin,ymax = 0.,1.1*np.amax(probb)\n", "plt.plot(blst,probb,color='black')\n", "plt.xlabel(r'$b$')\n", "plt.ylabel('Marginalized posterior probability')\n", "plt.xlim(bmin,bmax)\n", "plt.ylim(ymin,ymax)\n", "fig.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "By integrating under these distributions, we identify the 68% confidence regions:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Lower limit of 68% confidence region for a = 0.21537749641344503\n", "Upper limit of 68% confidence region for a = 0.3254708870429598\n", "Lower limit of 68% confidence region for b = 0.12026579012928584\n", "Upper limit of 68% confidence region for b = 0.7555475906345718\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOzdd3gc5dXw4d9Rl3s3xl3NvWAbMMVgakwzBAimfgFCzAskQEJCQl4gkAKEhASSAC8mAYcWwBBqTMdgwDbuTZZkFXcbXHBDVt/z/TG7sixLu7Nldma1z31dc620u5o5hjl7dp55iqgqhmEYhuE1KW4HYBiGYRgtMQXKMAzD8CRToAzDMAxPMgXKMAzD8CRToAzDMAxPSnM7gKZ69OihgwYNcjsM7ygpsR6HDHE3jgCvxeOQJUuW7FTVnm7H4YSkyzGvnLNeicMj7OaYpwrUoEGDWLx4sdtheMfkydbjJ5+4GcVBXovHISKywe0YnJJ0OeaVc9YrcXiE3RwzTXyGYRiGJ3nqCspo5oYb3I7gUF6LxzBC8co565U4EowpUF42bZrbERzKa/EYRiheOWe9EkeCMU18XrZpk7V5hdfiMYxQvHLOeiWOBGOuoLzsqqusR6/cWPVaPIYRilfOWa/EkWDMFZRhGIbhSaZAGYZhGJ5kCpRhGIbhSaZAGYZhGJ5kOkl42W23uR3BobwWj2GE4pVz1itxJBhToLzsvPPcjuBQXovHMELxyjnrlTgSjGni87KSkoOTTHqB1+IxjFC8cs56JY4EY66gvOz6661Hr4yd8Fo8hhGKV85Zr8SRYMwVlGEYhuFJpkAZhmEYnmQKlGEYhuFJpkAZhmEYnmQ6SXjZnXe6HcGhvBaPYYTilXPWK3EkGFOgvOz002Oym7q6Onbt2kWvXr1ISYniojlG8RhG3MTgnN2+fTvt27enffv2rsaRjEwTn5ctX25tUfjDH/5AdnY2ffr04YQRI9izZ4+r8RhGXEVxzqoqN998M71796ZDhw5874wzaGhoiHscycxcQXnZrbdajxGOnXj++ef55S9/ybmpqRyTmclvi4s5Y/hwPi4poWPHjnGPxzDiLopz9t577+Vvf/sb16Sm0j4lhb9/+CF3X3opv581K65xJDNzBdVGbdiwgWuvvZbJGRm82tDAXTU1vAIs3raNhy+91O3wDMPTVq9ezb333sv309L4Z0MDf62r4zrgvldeYc7MmW6HlzRMgWqjHnnkEXz19TyjSob/uanAecBfZs9m37ZtLkZnGN720EMP0S41lYd8PgQQ4G9AL+DPP/2pu8ElEVOg2qA9e/bw5JNPMk2E/nV1h7x2F7AbeOwHP3AlNsPwum3btvH8889zDdDd52t8PguYDvx3927WffyxW+ElFccKlIgMEZHlTbZ9InKrU8czDpoxYwbffvstt6WnH/ba0cAU4JF336Xh22/jHptheN1jjz1GfV0dt6amHvba9Vgfmo+b5TPiwrFOEqpaAowFEJFUYAvwmlPHa5Puuy+iP5s5cyYnpaVxVHV1i69fA0xT5Yu//IWT7rrL8XgMwzVhnrOqyosvvsgZqank1dYe9no/4ALgqeXLuW/PHtK6dHEkDsMSrya+04ByVd0Qp+O1Dccfb21hKCkpoaioiO+18O0v4Gys5opZTz/teDxGZETkKRHZLiKrW3n9ChFZKSKrRGSeiIyJd4wJIcxztqSkhLKyMs5vofUhYBqwC1jwt785FodhiVeBuhT4d0sviMh0EVksIot37NgRp3ASxLx51haGN954A4Dzm917aqoDVpF6dd06fAcOOBqPEbGZWK2xrVkHnKyqo4DfAjPiEVTCCfOcfeuttwA4r6am1fecAaQCs//d4kdaTOIwLKKqzh5AJAPYCoxQ1a+DvXfChAm6ePFiR+NJKJMnW49hjJ04/vjjqV20iMX19UHf92/gcuDzBx7ghF/8wrF4EpGILFHVCR6IYxDwtqqODPG+rsBqVe0bap9Jl2NhnrOTJk2i8ssvWRrkCx7AZGBPSgrL6+rAzuwsSZI7dtnNsXhcQZ0FLA1VnIzobdu2jQULFnCBSMj3ngtkAG889ZTjcRmO+wHwTmsvmlYKe3bs2MG8efM4z8aX9rOBFT4fW957z/nAklg8CtRltNK8Z8TWu+++i6oy1UaCdQSOBeaUlYHDV9GGc0TkFKwC1eplsKrOUNUJqjqhZ8+e8QsuwXz00Uf4fD7OsXFFdLb/cfZjjzkbVJJztECJSHusJtv/OHkcwzJ37lx6pKUxKkTzXsApwFKfj72LFjkbmOEIERkN/AM4X1V3uR1Povviiy9on5LCuBZ67zU3AugLfDR3ruNxJTNHC5SqVqpqd1Xd6+RxDMunn3zCSQ0NhG7gs0wGfMBn4fbmM1wnIgOwvvhdpapr3Y6nLfj888+ZiL2xNwKcAMzftw+qqpwNLImF/H8hIj8GnlPV3XGIx2jq4Ydtv3XTpk2sW7+eWzIzIUgPpKaOAzKBOe+8w7kxjsewRJo/IvJvrO8QPURkM/BrIB1AVf8PuBvoDjwm1j3Hei907PAcm+fsvn37WLlyJXelpkKT2SOCOR54Gdj8zjv0u/DCmMRhHMrOl4XewCIRWQo8BbynTnf9Myxjx9p+62effQbASWEsB5CFVaQ+2bTJug8VqnNFGPEYjSLKH1W9LMTr1wHXxSbENszmObtgwQJ8Ph8nBhn/1FxgVNP8V1/le6EKlMmdiIRs4lPVO4F84J/A1UCpiNwnIrkOx2Z8+KG12TB37lw6paYy2ub9p4DJwDKfjz3z58c0HsNi8sdlNs/Zzz//nBTgWJutD2BNk5MNzLNzH8rkTkRsTXWkqioiXwFfAfVAV+AVEflAVW93MsCk9rvfWY82VuOcO3cuJ6rS+vwRLTsRUGDRv//NGaFGuocRj3GQyR8X2TxnP//8c8ampdExjC946VhzW87bujV0C4TJnYiEvIISkVtEZAnwIPAFMEpVbwDGAxc5HJ9hw549eygqKuKEtPCnVhzvf1xoBhA6wuSP9/l8PhYtXMhxNu89NXU8Vk/YqtUtzkhlRMnOJ1o34MLm8+ipqk9EbN1bN5y1dOlSwPo2F64uQAGwqLw8liEZB5n88bjS0lK+raxkfFYWtDLBcmsmYl0SL3vpJY4fNcqR+JKZnW7mOc2TS0SeBVDVIkeiMsISmLrGzviNlhwNLKqqArP8hhNM/njcsmXLgMjyZ5z/cblpgXCEnQI1oukv/qUzxrfyXsMFixcvZnBaGt0j/PujsSZLNNO2OMLkj8ctXbqUDBGGR9DE1w+rr/+yIvNdwwmtNvGJyB3Ar4BsEdkXeBqoxcycHB9PPGHrbYsXLWJCGN3LmzvG/7jozTfpe1GQ2yI24zFM/niGjXN26dKljEpJIT2CHBKs3nzLd++2xk+1Nk2SyZ2ItHoFpar3q2pH4I+q2sm/dfTPDHFHHGNMXkOGWFsQu3btYt369UzIyIj4MGOxvqks+uKLqOMxLCZ/PCLEOauqLFu6lKOiGNo5FlilSt3aIBN6mNyJSLArqKGqWgzMEpFxzV9X1aWORmaAf20azjuv1bcsWbIEgAlRJFg2MBJYtHFj1PEYFpM/HhHinN24cSPf7N7NuAg6SAQcBdQAJe+8w8ihQyOKw2hZsF58twE/BB5q4TUFTnUkIuOgh/z/6YOc1NF2kAg4CphdVwe7d0PXrhHHYzQy+eMFIc7ZQAeJo0Ks/xRMYI6I5XPmMPInP4koDqNlrRYoVf2h//GU+IVjhGv58uXkpKXRJcwZJJobDTwNfP3pp/S+4IKYxJbMTP4khmXLlpECjI7iHu4QrGnDli1ZwpWxCswAgjfxBZ1cSlXNEhoesHLFiqiSK2B0YH8ffMAZpkBFzeRPYli9ejW5aWm0i+ILXhowCli2fXvM4jIswZr4gl2LKmaNJ9dVVVVRWlbGtPR0iLKJr7FAzZ/PGdGHZpj8SQiFq1czMkZf8N6sr4edO6FHj+gDM4DgTXzXxDMQI3xr1qzB5/MxysYKoKH0AI4EVpgZJWLC5I/3VVdXU1Zezvdi8AVvBNZswDvmzaPn1Kkxic8I3sR3pao+JyI/bel1Vf2zc2EZADz7bNCXV61aBcCoCHsfNTcaWLlvHzQ0QGoL086GiMc4yOSPRwQ5Z0tKSmhoaGBEGEtstGa4/3HNnDmc3FKBMrkTkWBNfO39jx3jEYjRgv79g768atUqslJSyItgBHxLRgMfAXXFxaSPGHH4G0LEYxzC5I8XBDlnCwsLARgRxhIbrWksUIsWcXKYcRitC9bE94T/8d74hWMc4qWXrMdp01p8eeXKlYwQCXuJjdaMAeqA4nffZVRLBSpEPMZBJn88Isg5W1hYSKoIBTFYf7Uf1jeRwtYG65rciYid5TZyROQtEdkhIttF5A0RyYlHcEnv8cetrRWrVq5kVIyunqBJR4lPP40oHuNwJn9cFuScLSwspCA1lcwYHEawrqIKv/km7DiM1tm5u/4C8DLQB+s++izg33Z2LiJdROQVESkWkSIROS7yUI2mtm/fztfbtzM6iimOmisAUoE1K1fGbJ9G5PljOGv1qlWMiEEPvoARwJqGBti1K2b7THZ2ClQ7VX1WVev923NY49LseAR4V1WHYrUgmSl/Y2S1f4G0kTG8gsrAWpu86KuvYrZPI6r8MRxy4MABKtati0kHiYDhwHZg57x5Mdtnsmu1QIlINxHpBrwjIr8UkUEiMlBEbgdmh9qxiHQGTsLqfYmq1qrqnlgFnuyK/NP7D49iipaWDAfW1NRADG4cJ7No88dwVklJCarK8BgM0QgI3LVdY9aGiplgvfiWYA0oFP/v1zd5TYFQMzIPBnYAT4vIGP/+blHVyqZvEpHpwHSAAQMG2I88yRUVFdExJYUjY3gFBTAMeAOoKSwkc9xhc5wa9kWbP4aDSkpKABgaoyEacLAnX+HChZwUs70mt2C9+AbHYN/jgB+r6pci8gjwS+CuZseZgX99nAkTJkTfnaYteeWVVl8qKipiGAc//WJlONAAlH78MSObF6gg8RiHikH+GLHQyjlbUlKCYDVpx0o/oB1QUlpqOw4juGBXUI1EZCTWZ1dj27mqPhPizzYDm1X1S//vr2AVKMOuIFOmFK1Zw5kOHDLwLbBowQJGhhGP0boI88eIhVbO2ZKSEgakppIdw04SKVgdjUpa6slnciciIQuUiPwamIyVYLOBs4DPgaAJpqpficgmERmiqiXAacCaqCNOJjNnWo9XX33I03v37mXbV18xPIo1bFozBOuqbM2KFbbjMVoXaf4YMdLKOVtSUsKQGDePg5U/i+rqoLIS2rc/+ILJnYjYuUN4MVZx+co/v9gYoLPN/f8YeF5EVmItm3JfRFEmq5kzD57YTQQ6SAyLwQDD5rKBHGDN1q224zGCijh/ROQp/9ip1a28LiLyVxEpE5GVLS2MmPRaOGdVlZLiYoa2NJ1XlAqA9Vj3cEPFYYRmp0BVqaoPqBeRTlg9KW3N26Gqy1V1gqqOVtULVHV3NMEaljVrrAvRYQ71tBsGrDlwAKJcY8oAosgfYCYwJcjrZ2HdRsnH6mhkRoLasGXLFioPHGBImq07HGEZAviAss8+i/m+k5GdArVYRLoAT2L1TFoKzHc0KiOooqIiMkVw6i78cGAtUO/v6WREJeL8UdW5QCtTEwBwPvCMWhYAXUSkT7QBt3WBHnxDYnj/KWCI/3Gtf6VrIzohv0Ko6o3+H/9PRN4FOqmqmWrARUVFRRSkpJDqQIKBlWS1wIbPPiO3pTn5DNsczp++wKYmv2/2P7et6ZvMUI5DNRaoGI8hBKuJD6CkpXu4RthsjVITkQtF5M9Y95RynQ3JCKW4qIhhDtzgDQh8CyxZuNCxYyQTt/NHVWf4m9on9OzZM96H95ySkhLap6TQ14F9dwKOAEo2b3Zg78nHTi++x4A8Ds4fdr2InK6qNzkamQGzD59woLa2lnXr13N5DBZZa01jgVq+nLNDxGME53D+bOHQ+1n9/M8ZAS2csyUlJRSIxHwMYcAQYO2334IqiLQahxGanbuEpwLDVK0uYyLyL6Aw+J8YMdGu3WFPlZeX4/P5KIjhFC3NdQe6Ams3bAgZjxGSk/nzJvAjEXkROBbYq6rbQvxNcmnhnC1du5ajHW6BeFUVvv4ajjii1TiM0Ox8ypUBTRuu+/ufM5z22GPW1sRa/3ozQxy6egJrHNQQoGTv3pDxGCFFnD8i8m+sDhVDRGSziPxARP5HRP7H/5bZQIV/f08CN7ayq+TV7Jytra1l/YYN5MdwktjmhgC7gF1Nm8hN7kQk2JLvb2HNGdYRKBKRwH/tYwBzcyIeXn7Zerzx4OdOoEDlO/gNEKwk+6ChAXbvhq5dW43HaFks8kdVLwvxugKmqT2YZufs+vXr8fl85DvYApHnfyz/8ku6B5Z/N7kTkWBNfH+KWxSGbSUlJfRKTaWLQz34AoYA/wL2L11Kx9NOc/RYbZTJHw8q9c+Tl+dAD76AwPx+pcuXc4xjR0kOwSaLbVxWVUR6A0f7f12oqtudDsxo2dq1axniwAwSzQW6y5Z+/jnjTIEKm8kfbyors1pX8x38gjcYq5m8rLjYsWMkCztLvl+C1STxPeAS4EsRudjpwIyWlRQVNRYPJzX25FuyJA5Ha7tM/nhLaWkpnVJScHLq1iz8Nxq//trBoyQHO734/hc4OvCtT0R6Ah9izU5uxNGePXvYvnMnBQ5MEttcHta3wJLmc4oZ4TL54yFlZWXkO9jFPCAPKD1w4NCu5kbY7BSolGZNEruwOcDXiFKzlTkD7efxaOLLAgYCa5su/25WCo2EyR83Nc+htWs5xuEORmDdh3pFFb76Cvr0MbkTITsF6l0ReY+DAw2nYZasdkVgipaCOC3Hng+UVlWBzwcO9npq40z+eESgi7mTg9wD8rC+iexetoyufcz0iJEK+qkjIgL8FXgCGO3fZqjqL+IQm/GnP1mbX2lpKSlYy2HEQz5QqooGpm1pFo8RnMkfD2hyzsaji3lAY1fzBQsOi8OwL+gVlKqqiMxW1VHAf+IUkxHw9tvW489+BlgFakBqKpkOdzEPyAf2AjsXL6bngAGHxWMEZ/LHA5qcs/HoYh7Q2NV82TImNIvDsM/OV4mlInJ06LcZTistLXV8gG5TjUkW+BZoRMLkj0cEupjnxeELXqCVw3Q1j46dAnUssEBEyv2rdq7yr5BrxJGqUrp2LfkOrALamqYDDo2ImfzxiPLycjqkpBCP+dyzsWbuLW3aycgIm51OEt9xPAojpF27drF33z7ys7LittLtYCAVKPNPr2RExOSPR5SXl5Mbhy7mAXlAeWWl1dXciIidBQs3iMg44ESsucW+UNWljkdmQHZ244+N7edxPNnTgUFA6fbth8Vj2GPyx2VNztnysjJGxOn+LVjNfLNVYedOkzsRsrMe1N1Yo+ADN3mfFpFZqvo7RyMz4J13Gn9snKIlTl3MAxq7mjc0HBKPYY/JH5f5z9mGhgbWrVvH1PR0iEMnCbBWpvwKqFy9mvYmdyJip4nvCmCMqlYDiMgDwHIgZIKJyHpgP9AA1KvqhMhDTW6BLuaD43zcPOALQDdtQgYNivPR24SI88eInS1btlBbV0duVlZcCxRAxZdfMuqUU+JyzLbGToHaijWxQGBunUzCW7XzFFXdGW5gBvDb31qPd91FaWkpA1NTyYhjEwVYV1D7ge2LFtH72Wcb4zFsizZ/jGj4c6j8xBMByI1j/gQKVPmKFYxqksuGfXYK1F6gUEQ+wGpDPwNYKCJ/BVDVmx2ML7l99JH16C9Q8exiHtC0q3nvwMSxJsnCYfLHTf4cKvfP5pAbp6snONjVvLywELb5Fzo2uRMWOwXqNf8W8EkY+1fgfRFR4AlVndH8DSIyHZgOMGDAgOYvG1hdzMtKS5mYmhq3HnwBgQJVvmIFJ8b1yG1GNPljxEh5eTlpIvSPYyejbkAXoHzLFujWLW7HbUvs9OL7VxT7P1FVt4hIL+ADESlW1bnN9j8DmAEwYcIE0x+zBYEu5nlx7GIeMBB/V/OyMjD3oMIWZf4YMVJeXs6glBTS4txEnguU79sX12O2JY5OSqWqW/yP27G+RZoFJiPQOALehfEU6VhFqmy7WWPPSFzl5eXkutBEngtU1NdbEy4bYXOsQIlIexHpGPgZOBNY7dTx2qTu3aF794MFKs5dzAPygLKqKquZont3V2IwjIh074527055WRm5cZyFJSAXWA/40tNN7kQg1GzmqSIS6RS8vYHPRWQF1oqi/1XVdyPcV3J69VV49VXKysoQ4t/FPCAPKAX04YetmAxboswfIxZefZXdTz7J3n37yEmzc8s9tnKAemDDtGkmdyIQajbzBhGJ6N64qlYAYyKKyjhEeXl5XGcxby4PqyvaN0uX0t10ZLEtmvwxYqeiogKAXBeayBu7mi9b5toXzERm5yvFMhF5E5gFVAaeVFWzfIDT7rgDsO5BudF+HhBY26bmwQfhyy/h/vtdiyUBmfxx0x130L6oCIAcF5rIAwWq+zvvWPlscicsdgpUFtbikKc2eU4x69s4b/58wJqs9cKUFGu6IRcEClRqSQlkZLgSQwIz+eOm+fNpt3Ej4E4TeV+sjkadd+1qzGfDPjvdzK+JRyBGy+rr69n5zTdWF3OXCtRgQICqAwdcOX4iiyZ/RGQK8AhWT/9/qOoDzV4fAPwLa7hNKvBLVTXLyTdTXVVFr9RUOrqQP6lYEy5XuZS7iS5kLz4R6Scir4nIdv/2qoj0i0dwBlRVWzPkuNHFPCALa22bqtpa12JIVJHmj4ikAo8CZwHDgctEZHizt90JvKyqRwGXAo/FOv62oKq6mhwX8ycHqIrz+MW2wk4386eBN4Ej/dtb/ueMOKiqqgLc62IekAdUmbEckYg0f44BylS1QlVrgReB85u9R4FO/p87Y837ZzRTXVXlaoHKBapNgYqInQLVU1WfVtV6/zYT4rIopdGvH9vT04GD83q5JQ+oAOjRw+VIEk6k+dMX2NTk983+55q6B7hSRDYDs4EfxyDeNqXhyCMpr6kh159HbsgB1qtS06uXazEkKjsFapeIXOkf05EqIldi3fQ1nPbcc/x53Dj6pKbS3uVQcrEWNdp3330uR5JwnMyfy4CZqtoPOBt4VkQOy2kRmS4ii0Vk8Y4dO2J06MSw7je/4UogJ8XRSXOCygGuAlb/4heuxZCo7Pxfuxa4BGvtrW3AxYDpOBEn5eXlrozfaC7Qk698/Xo3w0hEkebPFqB/k9/7cfgyHT8AXgZQ1flYtwsPu8RV1RmqOkFVJ/TsmVyNH4ExUDlxnMW8uUDrRyAWw76QBUpVN6jqVFXtqaq9VPUCVd0Yj+CS3q23ctWSJeSJuB0JucBfgA4PPuh2KAklivxZBOSLyGARycDqBPFms/dsBE4DEJFhWAUquS6RQjjigQf4C5DjYi+6HKzcGfTww67FkKha7WYuIrer6oMi8jesm7GHMOvYOK9h6VLyKiv52oUpWprLBfYAmaWlboeSEKLNH1WtF5EfAe9h9VZ+SlULReQ3wGJVfRO4DXhSRH7iP8bVqh643PaQ9qWljMPqneKWjsDRKSl0WrfOxSgSU7BPviL/4+J4BGIcrtrfgy/XA1dQHbEGHAa6vRshRZ0//jFNs5s9d3eTn9cAJ0S6/2RQVV1NFg4v22BDdmpqYz4b9rVaoFT1Lf9YjFGq+rM4xmT4NXYxd/EGb1PZmAJll8kfb6iuqiJLBFy+sMxKS2OfKVBhC/rJp6oNmG9orqny0BUUmAIVLpM/7lJVqqqryXY7ECA7LY2amhrqXOyskYjs3NxYbia7dMeGzEw2ZGYy2SMF6oAIq2trOaq6mqysLLfDSRQmf1zyzTffUNTQQJe0tLivRN1cbadOlOzfT7+NG8nNzQ39BwZgr2m26WSX5/m3c50MyrD8pl8/Hh092u0wGhWnpTEdWGdu9obD5I9LKioquB6oaNfO7VD4auJEKxbT1TwsZrJYDysrK2PixImw2hsLEQfuhZWVlTFs2DCXo0kMJn/cU15eDuDKSrrN5XTuDJgCFS47k8UWiMhHIrLa//toEbnT+dCSW21tLXesW8dPiovdDqXRmPp6noDGJeiN0Ez+uKeiooIngCEemIW/75df8o+UFFOgwmSnie9J4A6gDkBVV2INGjQctGHDBvKB/h5IroBMn4/hqamN30wNW0z+uKSiooKR6emke2CS45Q9exiZmWkKVJjsFKh2qrqw2XNmal6HBYpAdrYX+iBZRITsrCxzBRUekz8uqaioIMtD+ZOdlWUKVJjsFKidIpKLfzS8iFyMNaeY4aBAEfBSgQLIysoyV1DhMfnjkoqKCrI91Ns0Ozub8vJyzGQf9tkpUDcBTwBDRWQLcCvwP3YP4J/BeZmIvB1hjEmpvLyc1JQU0j22xHp2Vhbr16834znsiyp/jMjU1tayadMmT11BZWVns3fvXnbv3u12KAnDToFSVT0daw2boap6os2/C7iFg9O+GDaVlZWxsVs3ZOxYt0M5SISqnBzq6+vZuNHMF2xTtPljRGDjxo34fD6qCgrAA3NZ0rMntf6er6aZzz47ifIqgKpWqup+/3Ov2Nm5f2nrc4B/RBZe8iorK+O1k08GL82AnJHBrp/8BDA9+cIQcf4YkQs0Q++8807o0MHlaICTT6by978HTIEKR7DZzIcCI4DOInJhk5c6YQ0+tONh4HasuUYNmxoaGqioqGDq1Kluh3KYvEGDAMx9qBBilD9GhAJFwEuzNgwePBgwuROOYNe+Q7BGvHfBGv0esB/4Yagdi8i5wHZVXSIik4O8bzowHWDAgAE2Qm77Nm/eTG1tLdd9+ilceaXb4RxUW0uf3/2O7OxscwUVWlT5Y0SnvLycrKwsjrz9dti3z+1w4N136bhzJz179jRXUGEINpv5G8AbInKcf7XOcJ0ATBWRs7G+MXYSkedU9ZBPXFWdAcwAmDBhgunewsHmsx7V1bB5s8vRNKGKbNtGbm6uKVAhxCB/jCiUl5czePBgZMsW8MA4KL79FjZvJjc311xBhcHOPajvikgnEUn3j4jfISIhv9ar6h2q2k9VB2ENTKVV2wwAACAASURBVPy4eXEyWubFMVBN5eXlmQJlX0T5Y0SnoqLCU817AaZAhcdOgTpTVfdhNVesB/KAnzsZVLIrKysjMzOTzMxMt0NpUV5eHhUVFfi88M3U+0z+xJmqUl5e7tkCtWnTJmpra90OJSHYKVDp/sdzgFmqujfcg6jqJ6pqZnC2qaysjJycHLyxyMbhcnNzqampYbOXmh+9K+r8McKzY8cOKisrycnJcTuUw+Tk5KCqrF+/3u1QEoKdAvWWiBQD44GPRKQnYFatc1BZWRl5eXlw3HHW5hUpKTB+PPn5+YDpam6TyZ84a5zFPDfXyh8vjIPq0weOO67xqs4089kTskCp6i+B44EJqlqHteja+U4HlqwCzRN5eXlw//3W5hXp6fCrX1mxYQqUHSZ/4u+QAnX//d4YB3XCCXD//Y0FyvTksyfkVwsRSQeuBE4Sa2XXT4H/cziupLV161YOHDjQeJXiRf379yczM5PS0lK3Q/E8kz/xV1FRgYgwyD9mz0uOOOII2rVrZ66gbLLTxPc4VvPEY/5tnP85wwGBD/28vDy46CJr84qaGrjuOlJSUsjJyTFXUPaY/Imz8vJy+vbtS1ZWlpU/ez1w2+/tt+GiixARcnJyTIGyyU7j7NGqOqbJ7x+LyAqnAkp2gQ/9/Px82LXL5Wha4J/oMj8/3xQoe0z+xFlFRcXBDhK7doEXZg+vrm7M55ycHNPEZ5OdK6gG/3IBAIhIDtDgXEjJrbS0lIyMDPr37+92KEEFxkKZruYhmfyJM692MQ/Izc2loqLCLLthg50rqJ8Dc0SkAhBgIHCNo1ElsUAX89TUVLdDCSo/P5/q6mq2bt1Kv3793A7Hy0z+xFFlZSXbtm1r7MjjRbm5uRw4cICvvvqKPn36uB2Op4UsUKr6kYjkY80tBlCiqjXOhpW8SktLPd1BIiDwAVBaWmoKVBDR5I+ITAEeAVKBf6jqAy285xLgHqwFEVeo6uUxCTxBBe7teL1AgRWrKVDBhWziE5EsrEXX7gF+Ddzgf86IMZ/Pd3AMFMBpp1mbV6SkwIknApixUDZFmj8ikgo8CpwFDAcuE5Hhzd6TD9wBnKCqI7AWQ0xqhxWo006zhke4rX//xlw2wzTss3MP6hmsZQP+Bvzd//OzTgaVrLZt20ZVVdXBK6i77rI2r0hPB/96UP369SMjI8N0NQ8t0vw5BihT1QpVrQVe5PDxUz8EHlXV3QCquj1mUSeowId+4z2ou+6C9u1djMjv2GMbc3ngwIGkpqaaAmWDnXtQI1W16Te3OSKyxqmAktkhXcw9LjU1ldzcXFOgQos0f/oCm5r8vhk4ttl7CgBE5AusZsB7VPXd5jtKpiVtysrK6NGjB507d3Y7lFalp6czaNAgU6BssHMFtVREJgZ+EZFjgcXOhZS8Ah/2jVdQZ51lbV5RUwNXXNH4a35+vilQoTmZP2lAPjAZuAx4UkS6NH+Tqs5Q1QmqOqFnz54xOrQ3Nc7CEnDWWbBnj3sBBbz++iG5bFYEsMdOgRoPzBOR9SKyHpgPHC0iq0RkpaPRJZmysrJDu5hXVVmbl1QfnEauoKDAdDUPLdL82QI0HWvQz/9cU5uBN1W1TlXXAWuxClbSKisrO7SLuVfyp77+kFgCBcp0NQ/OThPfFMejMABYu3YteXl5nu9iHlBQUEBNTQ2bNm1i4MCBbofjVZHmzyIgX0QGYxWmS4HmPfRex7pyelpEemA1+SXtCNCamho2btyYEE3keXl57N27l127dtGjRw+3w/EsO93MN8QjEMMqUAUFBW6HYVsg1rVr15oC1YpI80dV60XkR8B7WPeXnlLVQhH5DbBYVd/0v3am/55WA/BzVfXg9CPxsW7dOlQ1YQoUWE2SpkC1zk4TnxEHDQ0NlJWVJWyBMmJPVWeraoGq5qrq7/3P3e0vTqjlp6o6XFVHqeqL7kbsrkQYAxVguprbYwqUR2zYsIHa2tpDC9S551qbV6SkwOmnN/56xBFH0KFDB1OgDE84rIs5WPmTkeFSRE0MHnxILg8ePBgRMQUqBFOgPCLwIT9kyJCDT/7sZ9bmFenpcMMNjb+KCPn5+aZAGZ5QWlpKp06dDm0y+9nPoF0794IKGD/+kFzOzMykf//+pkCF0Oo9KBHZjzV9SotUtZMjESWpwId8IjXxgRXvokWL3A7Dc0z+xF/gHq5/3S3Py8vLM8M0Qmj1CkpVO/qT6BHgl1gDB/sBvwAejk94yWPt2rV07tyZQ8apTJ5sbV5RU3PY+lQFBQWsX7+emhozPWNTJn/ir8V5LCdP9sY4qFdeOSyXCwoKTIEKwU4T31RVfUxV96vqPlV9HLNkdcwl2re/gIKCAnw+n1nfpnUmf+KgpqaGDRs2JFQLRH5+Pt988w27vLjum0fYKVCVInKFiKSKSIqIXAFUhvojEckSkYUiskJECkXk3ujDbbsSrYt5gOnJF1JE+WOEp7y8HFVNiJUAAkzuhGanQF0OXAJ87d++x+EDBltSA5zqX010LDCl6ZQvxkFVVVVs3Ljx0A4SCSKQZCUlJS5H4lmR5o8RhsOmCUsAgdwxzXytszNQdz0RNEmoNYfHt/5f0/2bmdejBYEpTxIpuQK6dOnCEUccQXFxsduheFKk+WOEJ3AVkkg5NHjwYFJTU80VVBB21oMqEJGPRGS1//fRInKnnZ37mzWWA9uBD1T1yxbeM11EFovI4h07doQbf5sQ+HAfNmzYoS9ccom1eUVqKpx33mFPDx061BSoVkSTP4Z9paWl9OjRg65dux76wiWXQGamO0E1lZ9/WC6np6czePBgU6CCsNPE9yTWomh1AKq6EmtesJBUtUFVx2L1XjpGREa28J6kmWm5NcXFxY1jig5x443W5hVpaXD11Yc9HShQZuLLFkWcP4Z9paWlLd/DvfFGyM6Of0DNjRnTYi6bnnzB2SlQ7VR1YbPn6sM5iKruAeZgJp5tUVFREQMHDqRd8wGFBw5Ym1eothjPkCFD2L17Nzt37nQhKM+LOn+M0NauXdty896BA9Z567a6uhZzJzDQ3Xy5a5mdArVTRHLx3z8SkYuBbaH+SER6BtamEZFs4AzAtAO1oLi4mKFDhx7+wtlnW5tX1NbCVVcd9nQgdtPM16KI8sewr7Kykq1bt7ZcoM4+G/bujX9Qzb3xRou5XFBQwIEDB9i6dasLQXmfnQJ1E/AEMFREtgC3AjcE/xMA+mCtHroSa+mAD1T17YgjbaN8Ph8lJSUtF6gEYQpUUJHmj2FTos7CAqYnXyh2evFVAKeLSHsgRVX329mxv639qCjja/M2b97MgQMHDu8gkUAGDBhAVlaWKVAtiDR/DPsCQxwS8Ute4KqvpKSEyV6aNcYj7PTiaxCRB4ADgeQSkaWOR5YkAh/qiZhcASkpKQwZMsQUqBaY/HFeq52MEkD//v3Jzs424whbYaeJr9D/vvdFpJv/ucSaj8fDioqKgMQuUGDFb5KsRSZ/HFZcXMzgwYPJyspyO5SwmS93wdkpUPWqejvwD+AzERmPGXAbM8XFxXTt2pUWu9hffXWL3bpdk5ra6risoUOHUlFRQVVVVZyD8jyTPw4rKSlpfRaWq68GLxSu4cNbzWUzjrB1dgqUAKjqS8A04Gkgx8mgkklxcTFDhgxpeZJYrxWotDSYNq3Fl4YPH46qmquow5n8cVDITkYJUqDWr19vvty1wE6Bui7wg6quBiYBNzsWUZJZs2YNw4cPb/nFnTutzStUoZWZl0eMGAFY/x7jECZ/HLR582aqqqpav4LauRN8vvgG1ZKqqlZzeejQoaiq6cnXgmALFp6qqh8DA0VkYLOXv23pb4zw7Ny5k+3btzd+uB/m4ovjG1AotbUwfTrMm3fYS/n5+aSmppoC5WfyJz5CdjK6+GLYty+OEbXiv/+FkhL45JPDXmo6TGP06NFxDszbgnUzPxn4GDh88jWrDf0/jkSURAoLCwFaL1AJJCMjg/z8/MZ/k2HyJx7aQi/Y/Px8RMTch2pBqwVKVX/tf7wmfuEkl7ZUoMD6d6xatcrtMDzB5E98lJSU0LlzZ3r16uV2KBFr164dAwcONAWqBcGa+H4a7A9V9c+xDye5rFmzhk6dOtG3b1+3Q4mJ4cOH89prr1FTU0OmF2aQdpHJn/gI2skogZiefC0L1kmiY4jNiFJhYSHDhw9P+OQKGD58OD6fzywfYDH5EweFhYVtogUiMI7Q54UOHR4SrInPLNHusMLCQqZOndr6G27wT9m2sPlk2C5JTYX/9/9afTnwQVFYWMioUaPiFZUnmfxx3q5du/j666+DF6gbboClS63ZxN00ejR8//utvjxs2DAOHDjAhg0bGDx4cBwD8zY7Ux1lichNIvKYiDwV2OIRXFu2Y8cOduzYETy5pk1rddyRK9LS4PzWF4ctKCggJSXFdJRoIpr8EZEpIlIiImUi8ssg77tIRFREJsQucu8LnGetDtMAK3+8MA6qoCBoLjf9cmccZGcc1LPAEcB3gE+xFh80E15GKdAdO2iB2rTJ2rzC54MtW1p9OTMzk4KCAlavXh3HoDwvovwRkVTgUeAsYDhwmYgc9kksIh2BW4DDVqtu62znUENDnCIKYv/+oLlsClTL7BSoPFW9C6hU1X8B5wDHOhtW2xf4EA+aXFdd1eL6S66pq4Obg48xHTVqlOnJd6hI8+cYoExVK1S1FngRaOny9bfAH4DqWAWcKAoLC+nQoQP9+/dv/U1XXWUVB7e9917QXO7SpQt9+/Y1BaoZOwUq0Hi7x79ke2cgcft0esTKlSvp1q0bRx55pNuhxNTo0aMpLy/n22/NWFS/SPOnL9D0K/dm/3ONRGQc0F9V/xtsRyIyXUQWi8jiHTt22I/c49paJ6MRI0aYAtWMnQI1Q0S6AncBbwJrgAcdjSoJrFy5ktGjR7eZ5AoIdI4widbIkfwRkRTgz8Btod6rqjNUdYKqTmhxUuIEtWbNmjbRgy9gxIgRrFmzhgYvNEl6RMgCpar/UNXdqvqpquaoai9V/b94BNdW+Xw+Vq1a1SanNQn8m1auXOlyJN4QRf5sAZq2XfXzPxfQERgJfCIi64GJwJvJ0lHCVg++BDNixAiqq6tZt26d26F4RsgVdUUkE7gIGNT0/ar6G+fCatsqKiqorKxskwVq4MCBdOjQwRQovyjyZxGQLyKDsQrTpcDlTf5+L9CjyXE+AX6mqotjFbuX2erBl2CadpTIy8tzORpvsNPE9wbWzdl6oLLJZkQo8OEdskDddpu1eUVaGlx/fdC3pKSkmI4Sh4oof1S1HvgR8B5QBLysqoUi8hsRCTJ4LjnY6mQEVv5kZ8chohDGjQuZy4Fia5rHDwp5BQX0U9UpjkeSRFauXElKSkro5DqvpXlGXZSaCmeeGfJto0aNYtasWahqm7vHFoGI80dVZwOzmz13dyvvnRzJMRLVypUr6dKlS/AefGDlkBem3crJCZnPnTp1YsCAAWaYRhN2rqDmiUjY0wKISH8RmSMia0SkUERuiSC+NmnlypXk5+fTrl274G8sKbE2r/D5oKws5NtGjx7N7t272RJkzFQSiSh/jOBWrFjBmDFjQn8BKimB+vr4BBXM7t22cnn06NGsWLEiDgElBjsF6kRgiX9E+0oRWSUidm4w1AO3qepwrBu4N7U00DAZrVixwt79p+uvD9mkFld1dfCLX4R825gxYwBYvny50xElgkjzx2hFWJ2Mrr8evDDk4aOPbOXymDFjKCkpobo66Ya1tchOE99ZkexYVbcB2/w/7xeRIqxxHEm9ot2+ffuoqKjgmmva7ioMgW+2y5Yt49xzz3U7HLdFlD9G6wKdjAJfhNqSMWPG0NDQQGFhIePHj3c7HNe1egUlIp38P+5vZbNNRAYBR9HCdCxtdRBhawJXFePGjXM5Eud07NiR/Px8li5d6nYoroll/hiHst3JKAEFiq5p5rMEu4J6ATgXWIK1AmjTxl4FcuwcQEQ6AK8Ct6rqYWsvq+oMYAbAhAkT1F7YiWvJkiVA2y5QYP375rWwNHwSiUn+GIdbsWKFvU5GCSg3N5f27dubAuUXbLmNc/2PEc/9LiLpWMXpeVU1S1wDS5cu5cgjj+SII45wOxRHjRs3jhdffJFdu3bRvXt3t8OJu1jkj9Ey252MElBqaiqjRo0yBcrPzkDdlr7q7wU2+MdqtPZ3AvwTKDKrhx60ZMkS+1dPd95pPQZbMyqe0tLgFnudMQP/xmXLlnH66ac7GZWnRZo/RutWrFjB0Ucfbe/Nd94JF18Me/c6G1QoxxwDN95o661jxozhpZdeMsM0sNeL7zFgAVYz3JP+n2cBJSISbFDMCcBVwKkisty/nR1twImssrKS4uJi+zc/Tz/d2rwiNRVOOsnWW4866iiApL4P5Rdp/hgt2LNnD+vWrbPfQeL00yEjw9mg7BgwwHYujxkzhj179rDJS0vtuMROgdoKHOWfbHI8MBaoAM4gyKSXqvq5qoqqjlbVsf5tdmvvTwYrVqxAVe1fQS1fbm1e4fOBzUGE3bp1Y+DAgSxbtszhoDwvovwxWhb4wjNhgs0pB5cv98Y4qB07bOfy2LFjAfPlDuwVqAJVbZx7Q1XXAENVtcK5sNqmQAcJ21dQt95qbV5RVwe//rXtt48bN47Fi5NiarhgTP7EUEQ55IVxUJ9+ajuXx4wZQ2pqauO/NZnZKVBrRORxETnZvz3mfy6Tg2vdGDYsWbKEXr16tbk1oFpzzDHHUFZWxjfffON2KG4y+RNDS5YsYeDAgW264027du0YMWKE+XKHvQL1faAMuNW/VQBXYyXXKY5F1gZ9+eWXHHPMMUlz4/PYY62FYxcuXOhyJK4y+RNDixcvTooBrBMmTGDx4sWotvmRN0EFLVAikgrMVtWHVPW7/u1PqnpAVX2q6oFr58SwZ88eiouLmThxotuhxM2ECRMQEb788rDx2UnB5E9s7dmzh/Ly8qQpUDt37mTjxo1uh+KqoAVKVRsAn4h0jlM8bVbgKiKZClTHjh0ZPnx40hYokz+xFXYHiQQWKMLJ3sxnZy6+b4FVIvIBTdaxUdWbHYuqDVqwYAEiYn/8BsB991mPXulqnpYGv/xlWH9y7LHH8sYbbyTzmA6TPzESdgcJsHLonHNgzx6HorLp+OPD6vA0evRo0tLSWLx4MRdddJGDgXmbnQL1H/9mRGHBggUMHz6cTp06hX5zwPHHOxdQJFJTIZwCi1WgnnrqKSoqKsjNzXUoME8z+RMjCxcuZNCgQeF1kDj+eEhPdy4ou448Mqx8zsrKYtSoUeYKKtQbVPVf8QikLVNVvvzyS7773e+G94dem8uuoQEWLQprZotjjjkGsDqIJGOBMvkTG6rKvHnzOPnkk8P7w3nzrOERbtu61YoljCJ19NFH89JLL+Hz+UhJsdOfre0J+a8WkXwRecW/8GBFYItHcG1FoKt1oFebbb/6lbV5RX09PPBAWH8ycuRIOnTowBdffOFQUN5m8ic2Nm3axNatWznuuOPC+8Nf/QoqK0O/z2nz5oWdy8cffzx79+6lqKjIoaC8z05Zfhp4HGsBwlOAZ4DnnAyqrQl8OIedXG1AWloaxx13HHPnznU7FLeY/ImB+fPnA8mVQ8f7r7aS9csd2CtQ2ar6ESCqukFV7wHOcTastmXu3Ll069aN4cOTc0HhSZMmsXr16mQdsGvyJwbmz59PdnZ2m1yksDV5eXn06NEjqZetsVOgakQkBSgVkR+JyHeBDg7H1abMnTuXk046KWnbkSdNmgQk7TdBkz8xMG/ePI4++mjSvdDhIU5EhOOPP94UqBBuAdoBNwPjsWYo/76TQbUlW7Zsoby8nJNszgLeFh177LGkp6fz2WefuR2KG0z+RKmqqoply5YlVfNewAknnEBpaSnJsNp4S+z04lvk//Fb4Bpnw2l7AvdeIipQDz9sPXqlu3l6Otx7b9h/lp2dzdFHH52UBcrkT/QWLVpEfX194z2ZsDz8MJx6KuzeHfvAwnHyyfDzn4f9Z4F/87x58zj//PNjHZXntVqgROTNYH+oqh5ZRc/b5s6dS8eOHSNrO/dPu+8ZKSkwcmREfzpp0iQeeughKisrad++fYwD8x6TP7HzySefICKNTcVhGTvWGmDutp49I8rn8ePHk5GRwWeffWYKVDPHAZuAfwNfAkk5DUC0Pv30U0488UTSIkmSDz+MfUDRaGiAuXNh2rSw//SUU07hD3/4A5999hlTpkxxIDjPMfkTI3PmzGHs2LF07do1/D/+8EOorY19UOHauNGKJcxZYbKzs5k4cSJz5sxxKDBvC3YP6gjgV8BI4BGsBdZ2quqnqvppPIJLdFu2bKGoqIhTTolw0urf/c7avKK+Hh55JKI/nTRpEhkZGXzwwQcxDsqzos4fEZkiIiUiUiYih80xJSI/9Y+vWikiH4nIwJj+Czygurqa+fPnR5dDBw7ENqhILFwYcS6feuqpLFu2jN1uN1O6oNUCpaoNqvquqn4fmIi1ZMAnIvKjuEWX4N5//30AvvOd77gcifvatWvHiSeemDQFKtr88c+E/ihwFjAcuExEmo9TWAZMUNXRwCu0wRV6FyxYQE1NTeQFqg049dRTUVU+/TT5rgtCLbeRKSIXYg0svAn4K/BaPAJrC95//32OOOIIRo0a5XYonnDGGWewatUqvvrqK7dDiYso8+cYoExVK1S1FngROOQmhKrOUdXA5cECoF9sIveOOXPmkJKSEtn9pzbi2GOPJTs7m48//tjtUOKu1QIlIs8A84FxwL2qerSq/lZVt8QtugTm8/n44IMPOPPMM5N1Fu/DnHHGGQB89NFHLkfivBjkT1+se1gBm/3PteYHwDutxDJdRBaLyOJE6648Z84cxo0bR+fOybtiSUZGBpMmTTIFqpkrgXyscRzzRGSff9svIvtC7VhEnhKR7SKyOlbBJpJly5axa9cuzjzzTLdD8YyjjjqK7t27NzZ9tnFR5U84RORKYALwx5ZeV9UZqjpBVSf07Nkzlod21N69e5k/fz6ne2W5GRedeuqpFBYWJk3rQ0Cwe1ApqtrRv3VqsnVUVTtrRswEkqK7VkveffddgOiS64knrM0r0tPhD3+I+M9TUlKYMmUKs2fPpqGhIYaBeU8M8mcL0L/J7/38zx1CRE4H/heYqqo1sYneGz788EPq6+s555woZoZ64gno4IGJO047LapcDtzHfuedFi+S2yzH5t5R1blAUk6+BvDGG29wzDHH0Lt378h3MmSItXlFSgrk5UW1i6lTp7Jz587GyT+NVi0C8kVksIhkAJcCh4ytEpGjgCewitN2F2J01OzZs+nSpUt0q1APGeKNcVBdu0aVy2PGjKFv377897//jWFQ3peck8M5bPPmzSxatIgLLrgguh299Za1eUVDA0TZPDdlyhTS09N5882g41iTnqrWAz8C3gOKgJdVtVBEfiMigUG+f8Sa12+WiCwPNTg4kagqs2fP5swzz4xsDGHAW29BjQcuLCsqosplEeGcc87h/fffp9YL47rixPUClcg3cFsT+PCNukA99JC1eUV9fdRNjp06deKUU07hjTfeiFFQbZeqzlbVAlXNVdXf+5+7W1Xf9P98uqr2VtWx/q3NzE6xfPlyvvrqK84+++zodvTQQ1BVFZugorF0adS5fO6557J///6kmjLM9QKVqDdwg3n99dcZMmQIw4YNczsUT5o6dSpr166luLjY7VAMj3rrrbcQkWSZdcSWU089lczMTN5++223Q4kb1wtUW7N7927mzJkT/dVTG3b++ecjIsyaNcvtUAyPmjVrFieeeGJ093DbmPbt23Paaafx+uuvo6puhxMXjhUoEfk31jiQISKyWUR+4NSxvOSVV16hvr6eiy++2O1QPKtfv36cdNJJPP/880mTaIZ9RUVFrF69mu9973tuh+I5l1xyCevXr2fhwoVuhxIXTvbiu0xV+6hquqr2U9V/OnUsL3nuuecYOnQo48ePdzsUT7v88sspKSlh+fLlbodieMysWbMQES666CK3Q/GcCy64gIyMDF588UW3Q4kL08QXQ+vXr2fu3LlceeWVsZk94tlnrc0r0tPhr3+Nya4uuugi0tLSeOGFF2KyP6PtCDTvHXnkkdHv7NlnoWPH6PcTre98Jya53LlzZ6ZMmcKsWbPw+XwxCMzbTIGKocCH7RVXXBGbHfbvb21ekZICfYPNtmNf9+7dmTJlCi+88AL19fUx2aeR+FasWMHq1au55JJLYrPD/v0hNTU2+4pGx44xy+Vp06axZcsWPv/885jsz8tMgYoRn8/H008/zaRJkxg0aFBsdvrSS9bmFfX1EMPu4ddeey1bt25NusGHRuv++c9/kpmZyeWXXx6bHb70ElRXx2Zf0Vi7Nma5PHXqVDp06MBTTz0Vk/15mSlQMfL+++9TVlbGDTfcELudPv64tXlFQwM880zMdnfeeefRt29fHvfSv9FwTVVVFc8++ywXXngh3bp1i81OH3/cGwVq5cqY5XKHDh244ooreOmll9r8GlGmQMXI3//+d3r37m1u7IYhLS2NH/7wh7z33ntUVFS4HY7hstdee409e/Zw3XXXuR2K502fPp3q6mqef/55t0NxlClQMVBRUcHs2bOZPn06GRkZboeTUK677jpSU1N59NFH3Q7FcNnjjz9OTk4OkydPdjsUzxs3bhzjx4/niSeeaNNDNUyBioGHHnqItLQ0rr/+erdDSTh9+/bl0ksv5YknnmDnzp1uh2O4ZP78+Xz++efccsstpKSYjyU7brzxRlavXt2mV6k2Z0KUNm/ezD/+8Q+uueYa+saoh1uyueOOO6isrOSRRx5xOxTDJX/84x/p1q0bP/hBUoznj4krrriCI488kvvvv9/tUBxjClSUHnzwQXw+H3fccUfsd/7KK9bmFRkZMGNGzHc7YsQIvvvd7/K3v/2NPXv2xHz/hrcVFxfzdJ3+OwAADfpJREFU+uuvc9NNN9G+ffvY7vyVV6CTneW3HHbOOTHP5czMTG677TY++eQTFixYENN9e4UpUFFYt24dM2bM4Pvf/37supY31aOHtXmFCHTv7siu7777bvbt28fvf/97R/ZveNf//u//0r59e370ox/Ffuc9eljj99yWne1ILk+fPp1u3bpx9913t8l7UR74P5e4brvtNtLS0rj33nudOcDMmdbmFfX1jo3LGjt2LFdffTWPPPIIZWVljhzD8J4vvviC//znP9x+++306tUr9geYOdMb3czXrHEklzt06MBdd93FBx980LiKd5uiqp7Zxo8fr4ni/fffV0Dvv/9+5w5y8snWlp2tCu5vKSmqxx3n2D9369at2r59ez333HPV5/M5dpxQgMXqgXxwYvNSjjU0NOjEiRO1T58++u233zpzkJNPVk1Pdz93+va1YnFATU2N5ufn67Bhw7Surs6RY8Sa3RwzV1AR2LdvH9OnTycvL4+f/OQnbofTZvTp04d7772Xt99+m2e9NAeh4Yi///3vLFiwgAceeCD2956SSEZGBn/6058oKiriwQcfdDucmDIFKgI333wzGzdu5JlnniEzM9PtcNqUW2+9lUmTJvHjH/+YjRs3uh2O4ZDy8nLuuOMOzjrrLK666iq3w0l4U6dOZdq0adxzzz2sWLHC7XBixhSoMP3zn//kX//6F3feeSfHHXec2+G0OampqcycORNV5cILL+TAgQNuh2TEWFVVFdOmTSMtLY0ZM2bEZuZ/g0cffZRu3bpx+eWXs3//frfDiQlToMLw8ccf8z//8z+ceeaZ3HXXXW6H02bl5OTwwgsvsHTpUq699tqkWFYgWagq119/PUuWLOG5556jX79+bofUZnTv3p3nn3+ekpISLr/8choaGtwOKXp2blTFa/PSDdzmPv30U+3QoYOOGDFC9+zZE5+DVlZam1c6SWRlqZaVxeffrqoPPvigAnrNNddoQ0ND3I6L6SThCJ/Ppz/96U8V0N/85jfxOWhlpWqPHu7nzo03WrHEwd///ncF9Ic//GFc8yYcdnMszdXqmCDeeustpk2bxsCBA3nvvffo3LlzfA7crl18jmOXSFxj+vnPf86BAwe45557qKys5Omnn6ad1/6bGLY0NDTw05/+lL/+9a/cfPPN3HnnnfE5cLt21nnrtvT0uOXOTTfdxNatW7nvvvuoq6tjxowZpKenx+XYMWenisVr89oVVG1trd55550K6Lhx4/Trr7+ObwCPPmptXrmCSk9Xve+++P43UOtKSkR07NixumbNGsePh7mCiqkdO3bomWeeqYD+5Cc/ie8QgkcfVe3Qwf3cmTzZiiVOfD6f3nvvvQroCSecoFu2bInbse2wm2OuJ0zTzUsF6pNPPtGRI0c2NjEdOHAg/kEk2TioYGbPnq3dunXTjIwMveeee3T//v2OHcsUqNhoaGjQ5557Tnv27Knp6en65JNPxu3YjZJgHFQwL7zwgrZr1067dOmiTzzxhNbX18c9hpbYzTHTSaKJuro6Xn/9dU455RQmT57M3r17eeONN3jqqafIzs52O7ykdtZZZ7FmzRouuOAC7rnnHnJycvjtb3/L1q1b3Q7NaKa6upoXXniBo446iiuvvJKcnByWLl1q1nlywWWXXcbSpUsZM2YM119/PSNHjmTmzJkJ0zvW0QIlIlNEpEREykTkl04eKxKqyubNm3nxxRe59tpr6du3L9/97ndZu3Ytf/nLXyguLmbq1Kluh2n49e7dm5deeon58+czbtw47r77bvr378/kyZP54x//yKJFi6itrXU7zJgJlT8ikikiL/lf/1JEBsU/Ssu2bdt4+eWXufrqq+nbty9XXHEF1dXVPPfcc3zxxReMHDnSrdCS3pAhQ5gzZw4vv/wyGRkZXHPNNRxxxBFceeWVvPDCC2zcuBHrosZ7HOskISKpwKPAGcBmYJGIvKmqa2J9LFWlvr6euro66urqqK2tpbq6mgMHDlBZWcnevXvZvXs3O3fuZOvWrWzatIny8nIKCwsb1yDq3LkzZ599Npdeeilnn302aWmm/4hXTZw4kXfffZeysjKeeeYZXnvtNW6//XYA0tPTGTZsGAUFBQwaNIgjjzySXr160aNHDzp37kyHDh1o37492dnZZGZmkpGRQVpaGmlpaaSkpHhmTI7N/PkBsFtV80TkUuAPwLRojuvz+WhoaGjMo5qamsY82rdvX2Mebdu2jY0bN1JWVkZhYWHjlWyXLl0455xzuPrqqzn11FPN2k4eISJ873vf4+KLL+azzz5j5syZvPXWW40r8vbq1YuRI0eSl5fHgAED6NOnDz169KBr16506tSJDh06kJ2dTVZWFhkZGaSnp8clZ8SpyikixwH3qOp3/L/fAaCqrS5ekpqaqoGmtKZxBX5u2jbp8/kOeQwjLnr37k1ubi5Dhw5lzJgxHHfccYwdO9Z7RSmwsujChVBV5WoogNUbauJEmDfP7UgOs23bNj777DOWLVvGqlWrKC0tZePGjVSHOVHoLbfcwiOPPLJEVSc4FKotdvJHRN7zv2e+iKQBXwE9NUhC9O7dW/fv339YW3+gMIWja9eu5OTkMGLECMaMGcOJJ57IUUcd5a0eY5MnW+drXZ27cfTtC3l58Mkn7sbRRENDA8uXL2fevHksX76c1atXs27dOnbs2BHWflJSUhq3QLESkcYtoOnPlZWVtnLMyQJ1MTBFVa/z/34VcKyq/qjZ+6YD0/2/jgRWOxKQPT0AN5d1bfH4QyA/A+Iyp9I3kNoNWvyk8kFDEZT4wMmRs27/PwAYoqod3QzATv6IyGr/ezb7fy/3v2dns315Kcec0up5kwc52eDq+ISv4P+3d2+hVlRxHMe/v8LoZj1UkJh0KiIz0jrYhQizy1OBEVr00EXwxSy6PvTgW0FBQRAUlaRUEhldHk5ESdYRCTETO8c8iaUSZEhKQZkPmvTvYZa1M497H5s9a/ac3wcGZvaZzf+/Zu/Z/zOz115r1x74+RieWofzoRs6OseyXzJExBJgCYCkDTn/cx3v8Q/lsNPHYEPO+GWr0znWLW5Xb+n0HOvmDeIfgSkt2+ekx8ysvU7On7/3Sbf4TufY/ks3q6VuFqgvgQslnSfpBOBOYKCL8cyapJPzZwC4N63PAz472vdPZr2ma7f4IuKgpAeAlcDxwLKIGGnztCXdyqdD4z0+5M8hd3yoQQ6jnT+SnqD4keMAsBRYLmkb8AtFEWsne9u6xO3qLR21q2udJMzMzP4P/0jBzMxqyQXKzMxqqfICVYfhWzrIYZakjZIOpt+jVB3/UUnfSNok6VNJ52bIYaGkryUNSfpc0rQq47fsN1dSSCq1q20H7Z8vaU9q/5CkxgwkJ+l2SSOS/iz7uOZQ9yHVjoWkZZJ2p9+6NYakKZIG0+fbiKSHjvqETkaULWuh+LJ3O3A+cAIwDEw7bJ9FwMtp/U7g7Qw59AHTgTeAeRniXw+cnNbvy3QMTmtZnwN8XGX8tN9EYA2wDphZcfvnAy+UedzrsgAXAxcBq8s8rpna0tF7qdcWYBbQD2zOnUvJ7ZoE9Kf1icC3R3u9qr6CuhLYFhE7IuIAsAK49bB9bgVeT+vvAjeq3MGe2uYQEd9HxCa6M2JCJ/EHI+LQcMPrKH4DU3UOv7VsngKU2Zumk/cBwJMU48uNbbyi8uI3UkRsiYitufMoSSNfy4hYQ9Ezs1EiYldEbEzre4EtwOTR9q+6QE0GfmjZ3sl/k/t7n4g4CPwKnFFxDt001vgLgI9y5CDp/jR8zjPAg1XGl9QPTImID0uM23H8ZG66zfqupClH+Lvll/t8tmOUvr65HPhitH3cSaLGJN0FzASezRE/Il6MiAuAx4GK5ugGSccBzwGPVRXzCD4A+iJiOvAJ/1zV9wRJqyRtPsLS81cX1vsknQq8Bzx82N2af6l6LL6xDN+ys0vDt+Qegqmj+JJuAhYD10XE/hw5tFgBvFRh/IkUg5quTnd3zwYGJM2JiDLGyWvb/ohofc+9SnEV2TMi4qbcOVQk9/lsYyRpAkVxejMi3j/avlVfQdVh+JbcQzC1jS/pcuAVYE5E7M6Uw4Utm7cA31UVPyJ+jYgzI6IvIvoovocrqzi1jQ8gaVLL5hyKe+VWP7nPZxuD1J9gKbAlIp5r+4QMvThupui5sR1YnB57guIDCOBE4B1gG7AeOD9DDldQ3MveR3H1NlJx/FXAT8BQWgYyHIPngZEUfxC4pMr4h+27mpJ7m3XQ/qdT+4dT+6eW/RrkWoDb0vt7f3qfrcydU9mvZa8vwFvALuCP9FotyJ1TSe26lqLD1aaWz7ebR9vfQx2ZmVktuZOEmZnVkguUmZnVkguUmZnVkguUmZnVkguUmZnVkguUmZnVkguUmZnVkgtUQ0iaJ2mdpOE0f9NZuXMyaxJJsyUtz53HeOIC1RyDEXF1RMygGNz0jtwJmTXMDOCr3EmMJy5QzTFf0npJwxSTPpY9h5LZeHcZMDnN9L1D0uzcCTWdC1QDSLqHYuK2G9IV1FaKceTMrDwzgL0RcRWwkGJCTesiF6hmuBRYGxG/S5oLXAN8nTkns8ZIU0ScCTyVHhpK29ZFLlDN8BqwSNJ6ihkqd0TEvrwpmTXKVIqp5Q+k7X6Kke6tizyauZlZG5LuppiOZSowgaIj0iMRsS5rYg1X9Yy6Zma9aAbwPrAWOAl40sWp+3wFZWZmteTvoMzMrJZcoMzMrJZcoMzMrJZcoMzMrJZcoMzMrJZcoMzMrJZcoMzMrJb+AoluM7OJtH1RAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Function which identifies the parameter values enclosing the central 68.27% of the probability (1-sigma region)\n", "# These values divide the total probability into regions containing ~16%, 68%, 16% \n", "def geterr(prob,val):\n", " prob1 = prob/(0.01*np.sum(prob)) # normalise the probability such that it sums to 100\n", " lik = 31.73 # corresponds to 100 minus the 1-sigma confidence probability\n", " likby2 = lik/2. # probability enclosed in each of the 2 tails outside the 68% confidence region\n", " dval = val[1] - val[0]\n", " cumprob = np.cumsum(prob1) # cumulative probability, ranging from 0 to 100\n", " botval = np.interp(likby2,cumprob,val) + 0.5*dval # interpolate to determine the lower limit of the 68% region\n", " topval = np.interp(100.-likby2,cumprob,val) + 0.5*dval # interpolate to determine the upper limit of the 68% region\n", " return botval,topval\n", "\n", "# Determine the lower and upper values of parameter a, enclosing 68% of the posterior probability distribution\n", "abot,atop = geterr(proba,alst)\n", "print('Lower limit of 68% confidence region for a =',abot)\n", "print('Upper limit of 68% confidence region for a =',atop)\n", "# Determine the lower and upper values of parameter b, enclosing 68% of the posterior probability distribution\n", "bbot,btop = geterr(probb,blst)\n", "print('Lower limit of 68% confidence region for b =',bbot)\n", "print('Upper limit of 68% confidence region for b =',btop)\n", "\n", "# Plot the probabilty distributions, shading the 68% confidence regions\n", "fig = plt.figure()\n", "plt.subplot(121)\n", "ymin,ymax = 0.,1.1*np.amax(proba)\n", "i1,i2 = int((abot-amin)/da),int((atop-amin)/da) # array indices corresponding to 68% limits\n", "plt.plot(alst,proba,color='black')\n", "plt.plot([abot,abot],[ymin,ymax],color='red',linestyle='dashed')\n", "plt.plot([atop,atop],[ymin,ymax],color='red',linestyle='dashed')\n", "plt.fill_between(alst[i1:i2],np.zeros(i2-i1),proba[i1:i2],color='red')\n", "plt.xlabel(r'$a$')\n", "plt.ylabel('Marginalized posterior probability')\n", "plt.xlim(amin,amax)\n", "plt.ylim(ymin,ymax)\n", "plt.subplot(122)\n", "ymin,ymax = 0.,1.1*np.amax(probb)\n", "i1,i2 = int((bbot-bmin)/db),int((btop-bmin)/db) # array indices corresponding to 68% limits\n", "plt.plot(blst,probb,color='black')\n", "plt.plot([bbot,bbot],[ymin,ymax],color='red',linestyle='dashed')\n", "plt.plot([btop,btop],[ymin,ymax],color='red',linestyle='dashed')\n", "plt.fill_between(blst[i1:i2],np.zeros(i2-i1),probb[i1:i2],color='red')\n", "plt.xlabel(r'$b$')\n", "plt.ylabel('Marginalized posterior probability')\n", "plt.xlim(bmin,bmax)\n", "plt.ylim(ymin,ymax)\n", "fig.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Activity\n", "\n", "Let's return to the same supernova distance-redshift dataset we were using in Class 3.\n", "\n", "- Convert the $\\chi^2$ values into a **joint 2D probability distribution** in $(\\Omega_m,\\Omega_\\Lambda)$\n", "- Marginalize this probability distribution to obtain the **1D posterior probability distributions** for $\\Omega_m$ and $\\Omega_\\Lambda$\n", "- Determine the **68% confidence regions** for $\\Omega_m$ and $\\Omega_\\Lambda$" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Monte Carlo Markov Chains\n", "\n", "The grid method becomes inefficient as the number of parameters increases. A powerful alternative is to generate a **Monte Carlo Markov Chain** (MCMC) in the parameter space.\n", "\n", "The end result is a \"chain\" (distribution of parameter values) which **samples the underlying probability distribution**.\n", "\n", "Here is a worked example of using python's emcee algorithm to sample the probability distribution of the straight-line fit:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWwAAAGACAYAAACJGVhAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOzdeXhU5fn4//cz2SfJJIRsZCEhCQFZImjYRUBBf6DUWlGpRUttxbp8LK21fvCS6teqtbV1+ViXYlsVra1WsIioSBWQfTEswUBCCAkJSci+TmaSyTy/P5JMERIIgTA5mft1XXM1mXPucx6m8c6T+zyL0lojhBCi7zO5uwFCCCG6RxK2EEIYhCRsIYQwCEnYQghhEJKwhRDCICRhCyGEQUjCFheNUup3SqlNSqm3lVI+3TmulJqulPpCKbVeKXXjxW+1EH2HJGzRK5RSb57y/aVArNZ6KnAImHe240qpAOBBYLbWeobW+sOL0ngh+ihJ2OJimQx83v71Z8CUbhyfBDQBq5VSHyqloi9GQ4XoqyRhewilVFh70mtUShUopW7r5Bw/pdRf24/XK6X2KqVmn3JOwymvVqXUS91owgCgrv3rWiCsG8ejgBRgLvA68Hh3/71C9Efe7m6AuGheBpppS4JjgDVKqX1a629OOscbKASmAceAOcD7SqnRWut8AK11UMfJSqkgoBT4V/v3g4Hl7YeHK6U2tH99DVADWNq/DwGqTmlfZ8drgC1a62al1BfAkp7+44XoD6SH3QcopSKVUh8ppU6092xXK6UsZ4/s9vUDgZuApVrrBq31ZuAj4PaTz9NaN2qtH9da52utnVrrj4GjwOVdXPomoAzY1B5/TGs9XWs9Hfis42utdTOwFZjZHnctsOWUa3V2fBdwiVJK0fZLJq+HH4EQ/YIk7L7BArwEDAYSgHDg7s5OVEp9rJSq6eL1cRfXTwUcWuuck97bB4w8U6OUUlHtsd90ccoPgeW6GyuIaa33AieUUpva77tCKRWtlPp/XR3XWlcAHwIbgd8DT5ztPkL0Z0pW6+t7lFJPAN5a60cu0PWmAv/SWkef9N5dwA/ae8OdxfgAnwJHtNan/fJQSiXQ1uNN0VofvRDtFEKcmfSw+wCl1M1KqS1KqTKlVA3wv0DO2eLOQQP/rQ93sAD1XbTHBLxNW837/i6ueTuwWZK1EBePJGw3U0pdBfwOWAzE0FYOKQP2dnH+p52M1Oh4fdrFbXIAb6XU0JPeu5ROSh3t9eK/0vZw8iatdUsX17wDeKsb/8STry0TZ4Q4D5Kw3e9S2kZm7KNtaNvfgEggq7OTtdaztdZBXbxmdxHTCKwEnlBKBSqlpgA30NaLPtWrwCXAXK11U2fXU0pNBmJpHx3SHTJxRojzJwnb/f4O+NA2jO1j4DCQ1T6y4kK6Fwigrff+D+CejiF97b32R9rr0nfTNiKj9KSe+w9OudYPgZVa605LKl2QiTNCnCcZh+1mWusy2pLVyX7TC/epAr7bxbGTe+aqG9fqdATLWQwAStq/7mrizKnHOybOTKRtyN/jwE97cG8h+gVJ2OKCae8B/7OTQ/ORiTNCnDcpiYgLRmtdetJkmZNfpcjEGSHOm/SwxUWhtd7bPpNzE23T3v/Q3iO/R2v9WGfH23vWHRNnNHCn+/4FQrifTJwRQgiDkJKIEEIYhCRsIYQwCI+tYYeHh+vExER3N6NXFVU2AhA3MBCAI0eOAJCcnNyj67k7vi/5+uuvK7TWEe5uh/AsHpuwExMT2b17t7ub0aseWr4NgGfvmOTmlvQ/SqkCd7dBeB4piQghhEFIwvYgS5YsYcmSns89cXe8EJ7OY0sinmjbtm2GjhfC00kPWwghDEISthBCGIQkbCGEMAipYXuQuLg4Q8cL4ekkYXuQd955x9DxQng6KYkIIYRBSML2IIsXL2bx4sWGjRfC00lJxIPs3dvpRuyGiRfC00nCFp1KTEykoKDz5TKUUiQkJJCfn39xGyWEh5OELTpVUFDAqZtbTJ8+HYANGzbQtmuXEOJi8qgatlJqkVJqt1Jqd3l5ububI4QQ58Sjetha62XAMoD09HSP2xstNTXV0PFCeDqPStiebtmyZYaOF8LTeVRJRAghjEwStgdZtGgRixYtMmy8EJ5OSiIeJCcnx9DxQng66WELIYRBSMIWQgiDkIQteiQhIQGlVKevxMREdzdPiH5JatgeZMyYMRcs/kzT0ruaBXm+9xfC06lTpx97ivT0dL179253N6NXPbS8bdPbZ++YdM6xSqnTpqZfjFijUEp9rbVOd3c7hGeRkogQQhiEJGwPsmDBAhYsWGDYeCE8ndSwPUhRUZGh44XwdNLDFkIIg5CELYQQBiEJWwghDEJq2B5k0qRvD+870zZgCQkJZ40/3/sLIc6NjMPux842Dru3xkvLOGwheoeURIQQwiAkYXuQm266iZtuusmw8UJ4Oqlhe5DKykpDxwvh6Tyqhy27pgshjMyjErbWepnWOl1rnR4REeHu5gghxDnxqIQthBBGJjVsD3L11VdflPiOzQ06ExIScl5tEMKTScL2IEuXLr0o8T3Z3EAIcXZSEhFCCIOQhO1BZs+ezezZs90aL4ToOSmJeJCmpiZDxwvh6aSHLYQQBiEJWwghDEISthBCGITUsD3I9ddf7/b4jRs3ntc1hPBksh52P+au9bDPpL+slS3rYQt3kJKIEEIYhCRsDzJ9+nSmT5/u1nghRM9JwhZCCIOQhC2EEAYhCVsIIQxCErYQQhiEjMP2ILfccovb42UcthA9J+Ow+zEZh917ZBy2cAePKol4+ia8VqsVq9Xq1nghRM95VML29E1458yZw5w5c9waL4ToOY9K2EIIYWSSsIUQwiAkYQshhEHIsL5+bvv27agfTv7Wex07lyckJLijSUKIHpKE3c/ZbHbXMLo333wTgIULF/boWj2NOzlexmEL0XOSsD3IhUi45xv/ox/96LyuIYQnkxq2B6moqKCiosKt8UKInpMetgeZN28eABs2bHBrvBCiZ6SHLYQQBiEJWwghDEISthBCGIQkbCGEMAh56OhB7rnnHrfHyzhsIXpOErYHufXWW90eP3/+/PO6hhCeTEoiHqSwsJDCwkK3xgshek562B7k9ttvB3o+jvpCxQshekZ62EIIYRCSsIUQwiAkYQshhEFIwhZCCIPwqIeOSqlFwCKAwYMHu7k1F9+DDz7o9visrCzXBgqnSkhIID8//7zuIUR/5lEJW2u9DFgGkJ6ert3cnItu7ty5bo8vKyvr8nhXiVwI0UZKIh4kOzub7Oxsw8YL4ek8qoft6e6++26g5+Oo3R0vhKeTHrYQQhiEJGwhhDAISdhCCGEQkrCFEMIg5KGjB3n00UcNHS+Ep5OE7UFmzpxp6HghPJ2URDzI3r172bt3r2HjhfB00sP2IIsXLwZ6Pg7a3fFCeDrpYQshhEFIwhZCCIOQhC2EEAYhCVsIIQxCHjp6kKefftrQ8UJ4OknY/UBiYiIFBQWnvZ/+o9/i7+/n+n7y5MnndR93xwvh6SRh9wMFBQVoffp+DA8t3/at77du3Qr0PHG6O14ITycJ24M88sgjQM/HQbs7XghPJw8dhRDCIDwqYSulFimldiuldpeXl7u7OUIIcU48KmFrrZdprdO11ukRERHubo4QQpwTj0rYQghhZPLQ0YO88MILho4XwtNJwvYgY8aMMXS8EJ5OSiIe5D//+Q//+c9/DBsvhKeTHrYHefLJJ4Ge7/zi7nghPJ30sIUQwiAkYQshhEFIwhZCCIOQhC36jISEBJRSnb4SExPd3Twh3E4eOnqQP//5z306Pj8/v8tjSqnzurcQ/YEkbA8ybNgwQ8cL4emkJOJBVq9ezerVqw0bL4Snkx62B/njH/8IwNy5cw0ZL4Snkx62EEIYhCRsIYQwCEnYQghhEJKwhRDCIOShowd5++23DR0vhKeThO1B4uPjDR0vhKeTkogHee+993jvvfcMGy+Ep/OoHrZSahGwCGDw4MFubs25SUxMpKCgoNNjCQkJ3brGq6++CsCtt97aoza4O14IT+dRCVtrvQxYBpCenq7d3JxzUlBQgNaGarIQ4gKTkogQQhiEJGwhhDAISdhCCGEQHlXD9nQffPCBoeOF8HSSsD1IeHi4YeM7dqM50/EzbYAgRH8gCduDvPnmmwAsXLjQcPFnS8ayI43wBFLD9iBvvvmmK2kaMV4ITyc9bCFEv6CU+h0wGcgH7tRat5x0LAr4EGgBWoEfAM5T39Nal1zkZp8T6WELIQxHKfXmKd9fCsRqracCh4B5p4RUAFdoracBy4Efd/FenyYJW4h+TilVqJS63N3t6GWTgc/bv/4MmHLyQa11q9ba2f5tMPBNZ+9dlJaeB0nYQvRjSqlQIAY42Iv3CFNKfaiUalRKFSilbutGzFCllE0p9U5PjndiAFDX/nUtENbJNccopXYA9wMZXb3Xl0kNu4840+JO0P0Fns7kk08+MXS86JHRwDGttbUX7/Ey0AxEAWOANUqpfVrrM/VYXwZ2nctxpdRg2koXAMOVUhvav74GqAEs7d+HAFWnXlBrvReYoJS6BVgC/LSz987QJreTHnYf0bG4U1evCzHG2Gw2YzabDRsvemQ0cEQp9apSqkopdVgpdcWFurhSKhC4CViqtW7QWm8GPgJuP0PMfNoS7BfnclxrfUxrPV1rPR34rONrrXUzsBWY2X7qtcCWU67pe9K3tYC1s/fO9u91N0nYHuSVV17hlVdeMWy86JHRQDqwBggH3gFe7+xEpdTHSqmaLl4fd3H9VMChtc456b19wMgu7mEBngB+0ZPjXWnvKZ9QSm1qv/cKpVS0Uur/tZ8yRin1lVJqPbAYeLaL9/o0KYl4kPfffx+Ae++915DxZ3KmmZAePgsyDXhOa/0xgFLqL8BjSilvrbXj5BO11tf34PpB/Ld23KGWtod4nfkN8FetdVEX/3+d7XhHWxd28t5Dp7xVCjzWfmwncOUpx0s6ea9Pk4Qt+oUzJWQPnwU5CrjrpO/DgZpTk/V5aOC/teMOFqD+1BOVUmNoK1uM7exCZzsuJGFfVBdi1xghuksplUBb8iw/6e0baSuPdHb+p8DULi63SWs9u5P3cwBvpdRQrfXh9vcupfMhctOBROBY+y/RIMBLKTVCa31ZN46fUQ8mzliBdcAIYKLW+sDZ7uFuUsO+iM70YNGD/2TvdR3lks5eiYmJ7m5ebxoNOIDblFImpdR1tI2CeKKzk7XWs7XWQV28OkvWaK0bgZXAE0qpQKXUFOAG4O1OTl8GJNM2kmQM8Bptvzyu7ebxLvVw4owVuA4wzDKSkrAvsMTExC6Tg/Si3SM/P7/LX5RnGkrZD4ymLTlNAappq+fecFJP+EK5FwgAyoB/APd0DOlTSn2qlHoEQGtt1VqXdrxoK6fYtNbl3Tl+Fj2ZONPSzWv3GcqT9gk8eRNeYBiQfZGbEE7bb3q5p/HvOUxr3dWDNXGRtf9SyNJa/1splQI8obW+7ZRzxgB/BkKBa7TWBe3vvwn8wQglEY+qYZ+8Ca87KKV2a63T5Z7Gv6dSavfFupdoo5SKBv7ZyaGOcdvnPHGml5raazwqYQshjKu9TDK9s2NKqa20jd1eThcTZ9on2IBBJsl0RmrYQgjD6+HEGZRSn9A2tf11pdRCNzT9nEgP++JyRzlG7tk/7ifOogcTZ9Baz7kITbtgPOqhoxBCGJmURIQQwiAkYQshhEF4bA07PDxc9/NZbhRVNgIQNzDQzS3pf77++usKrXVEV8f748/XkSNHAEhOTu7Tsf3h576rny+PTdiJiYns3t2/h9I+tHwbAM/eMcnNLel/lFJnnCLpCT9ffVV/+Lnv6udLSiJCCGEQkrCFEN2yZMkSlixZYqjY/sZjSyJCiHOzbds2w8X2N9LDFkIIg5CELYQQBmH4hK08fP8nIYTnMGwNWynlDzRrrZ1KKS+tdau72yREfxYXF2e42P7GkAlbKXU9cAswQCn1Q631aWvfdhHn2sBg8ODBvdhCIfqfd955x3Cx/Y3hSiJKqRnAU8DrtG1J9Gp3Y7XWy7TW6Vrr9IiILiepCSFEn2S4hA1cDazQWm+ibU3bJqXUw0qpoUopLze3TYh+a/HixSxevNhQsf2NEUsih4CZSqmfAT8H/gXEAS8CDwOZbmzbaZxOJ1arFbPZjMlkxN+PQrTZu3ev4WL7GyMm7C1AIDAU2NaxaLlS6hngIeCOi9WQzpJxQ0MDO3fuZPz48QQFBWG1WqmrqwMgKCjoYjVNCLdKTEzsdEf6jkFdCQkJ5OfnX+RWGZ8hEvbJo0C01keBPyulBgP3KaWitNYngINAwMUcMWK1WqmpqaGhoYHIyEhMJhM7d+50zcy68sorqaurw2w2Yzab0VojoxCFJygoKODkzVGmT58OwIYNGwDkv4Me6tMJWymVqrXO0Vq3dpKIW4DhwC+VUn60bf9z+8Uc3mc2m2loaMDhcGC1WgkKCmL8+PEAjB8/noqKCkpLS4mJiSEkJORiNUsI0U/12YTdPnTvfaXUv7XWt52ctJVSSmtdopR6DEgHEoBbtdbZF7INZ6s/m0wmIiMjXecABAYGMmPGDAD8/f0BCA8Pv5DNEsItUlNTDRfb3/TJhK2UCgTup21348lKqXe01gvak7W31trRfuqh9t2Se0V36s8mk8l1rKqqijVr1nDdddcxYMAAAJqbm8nOziY4OJiYmBi8vfvkRy7EWS1b1vN9h90V29/0yeyhtW5USt0J1AErgddOStoOAKXUGOAKpdRfALvuhd2EO3rNHf97NmvWrOGLL74AYMGCBVRUVLBnzx6Ki4sZNGgQLS0tDBkyREaLCCF6pE8mbACtdXH7lw1KqbuBZR1JWymVBiQD72utbb3VhpN7z91x3XXXuf63urqazz77jDFjxpCSkoLJZMLHx4e6ujpsNhvh4eFn7G3LcEDR1yxatAjoWY/XXbH9TZ9N2CfTWle2J+1nlVLZtE34uVJrXebmpn1LWFgYCxYsANqm027YsAEvLy8WLFiA1tpVYiktLQUgOjq6y2vJcEDR1+Tk5Bgutr8xRMIG0FpXKKX2A7OBWVrrEne3qTMdlZnrrruO5uZmLrnkElpaWvD29sZsNuPv74/JZCI8PJwzVXFOLsfIcEBjkLVqRG8zzN/aSqkBwBzgGq11n5rN6HA4KC0txeFwoJTC6XTS3NzM5MmTKSoqoqCgAKUUSim8vLyIiorCZDLR2NjYZdJWSmE2m7FarWdM7KLvkLVqRG8zTMLWWlcDc7XW+93dllNVVFRQXFxMRUXFt74HiIyMZODAgafFdJQ8rFZrl9ftzjlCCM9hmJIIQG8+YOwJh8NBRUUFoaGhwH/HW1ssFvLy8oiJiSE2NtZVg25tbaWiooKwsDCcTidmsxmn0+kqeXTUuc1ms6uHDd0fpSJEbxozZozhYvsbQyXsvubknnTHA0Sn00lubq7rwWJaWpqr/lxRUUFeXh7btm0jJCSElJQUnE4n+fn5DB8+HKvVSnZ2Ni0tLaSnp+Pn50dgYKB7/nFCnOKFF14wXGx/Iwn7PHT0qMPDw13D8JxOJxaLhcGDBxMYGIivry/Q1rt2OBw0NDSQm5uLt7c3sbGx2O12jh8/TmBgIH5+fmRkZNDQ0IC/v7/rl0B0dLRMuBFCSMI+H97e3kRHR+N0OikrK8PhcGCxWAgPD8dkMnHw4EFMJhNDhw6loqKCsrIyhgwZQmNjI4mJiQwZMgQAPz8/vLy8aG5uZvDgwVRXV7N3716qq6tJTU1lwoQJrvvI2GzhLicPWTVKbH8jCfsCsFqtOBwOvL29CQoKQilFTU0NNpuN+vp6V5ItKSlx1aYtFgs+Pj4AhISEcOTIEQBiYmLQWvPFF19QVFREXl4el1xyiauO3dDQgNPpxGQy4e3tzfHjx4mPj3ddy263k5OTQ2pqKlqDjAYUF0pRUZHhYvsbSdgXwMkPBzvq1dHR0YwbN861tOrOnTspKCjA29ub+vr6b40c6SipdIzRjoqKorq6mrKyMjIzM/nzn//MnXfeSWxsLAB1dXXU19fz6aefYjKZmDVrFiNGjADaJhkcOHDgIn8CQoiLQRL2eTq1TNExZtrLy4tBgwZht9s5cOCAa8Uxq9VKWVkZWVlZJCUlsWfPHmJiYrjkkkvw8fGhurqagwcPMnr0aJRSrFq1igEDBlBXV0dsbCyNjY04nU6+/vpr9u3bh5eXF1OmTAHaJu3ExMTgdDrb7rdnj9s+F9H/dbVJAbRtUCAuPEnY5+lsU8gPHTrEzp07SUpK4oorrqCkpG2CptaalStXsnfvXvz9/ZkzZw5eXl74+fmxevVqrrjiCkJDQ5k7dy4BAQHY7Xaam5sJDAzEbDYzffp0Dh06xIgRIzCbzbS2tmKz2WhubiYlJQU/P7+L+jkIz3PqJgXnIiEhocvZu7IbTdckYZ+nzsZKnzyeOi4ujm+++YaSkhIyMjIoLy8nJCSEQ4cO0draypAhQ6itrWX9+vV4eXm5hgWWlJQwYsQI4uLiqKqqory8nJaWFuLj4/nmm284ceIEXl5eVFdXU1tby9GjR4mIiCAoKOiUtkBNTQ01NTXY7XaGDBniGrkixLmYNGnSBYs9U0I+NZGfz337G0nYZ9CdURmnruintXaNGAHw9fWlsrKS2tpaQkJCaGpqwmKxUFFRwfDhw4mLiyM/Px+LxUJtbS1hYWGUlpZSVVXFtm3biImJwWKxEBQURFlZ21pXzc3NxMfHEx8fT2xsLFu3buXKK6/E6XQSHBwM/Le373S2kp+f7/oFYbfbGTVqlIwyEefst7/9reFi+xtJ2GdwrivmOZ1OSktLXcnZ19eXjz76iG+++YampiauuuoqwsPDOXbsGJWVldTU1ODn50dAQABmsxmHw8H7779PZGQkUVFRFBcXY7FYSExMpLW1lU8++YT6+nomTpzI5ZdfzsCBA/n44485fvw4+fn5XHHFFZw4cYINGzaQnJxMdXW1Kz40NJSqqiosFotrOzMhhLFIwj6Dc50abrVaqa+vx263ExQURFVVFdHR0SQmJrrGandMTY+KisLHx4fa2lpiYmIoLi4mMzOTY8eOYTKZGDFiBHa7ncOHD7umsvv5+VFWVkZtbS0nTpwgIyODAwcOEBERQWhoKEFBQWzatIlVq1YRHBxMddwMWlqaMZvNhISEEB8fT1NTEwEBAa7a48nlG+l1izO56aabAFixYoVhYvsbSdhncK4bGEDbg5jk5GTXcD5oWwshLy+PyspK9u3b5xrel52djVKK4cOHExsbi9lsJjo6muPHj1NVVQVAdXU127dvp66ujoiICOLi4vD29ubo0aN89dVXlJSUMHDgQCorKykrK3OVU7y9vdFOJ06n5sCBA6SkpODv7++aRdmhoaGBEydOEBUVhcViuXAfnuh3KisrDRfb3/SLhK2UMmmtne5ux4EDB8jMzMRutxMVFYXZbCYoKAh/f38SExOJjIwkMzOT3NxcwsLCSEpKori4mPr6enJycjh69CiRkZH4+PiQmZnJ2LFjaWxsxN/fn6amJrTWjBw5ksbGRrKzs4mOjiYqKgovLy+OHTvG1q1bmThxItdeey12u5392s9Vt87LyyMwMJD6+nqgbex3RUWFPIAUwkAMmbCVUjOB8YAZ+G37HpDqbPs69vYC82lpaTidTmJiYti/fz8pKSmu9UOqq6vJz8+nqamJxMRE0tLSaGhooL6+nry8PPz9/SktLXUtEFVcXMzBgwcJDg4mPj6eoKAgfH19sdlsBAYGcuzYMaqrqxk0aBA2m42GhgZKSkpco0u++uortEUTGGh2lWc6Pp7Q0FDXwlXR0dEMGjRIVgQUwgAMl7CVUtcBvwVeBYYBnyulrtJa288Wq7VeBiwDSE9P75VNeydNmsSOHTvYs2cPTqeTlJQUwsLCqK2tZe3atezcuZOrr76a5ORkcnJyOHHiBDt37iQuLg6n08nx48dpbGx0bYRQW1uL2WympaUFm83G4cOH8fLywtvbm7q6Ok6cOEFdXR0BAQGEhISQk5PDnj17yM/PJ3RqKk5nKzZbBMePHyczM5OIiAgaGhqIi4sjMjKSyMhIAMrKys66z6QQwr0M9V+nUmoQcB/wgNZ6A/CqUuotIAX4xp1tO1lKSgqZmZnExcXR3NxMYWEhwcHBaK1dm/EWFRVx6NAhAgMDCQ8Pd20fVlpaitaa8PBwGhoasNvtaK0pKirC6XTS0tLiuo+Pjw+hoaForbHZbNTV1eFwOIiNjSUmJoYmpWhqsrFr1y7XQ8v6+nqOHz9ORUUFkydPxtvbm+LiYvLy8lx/HQjRmauvvtpwsf2NoRI20Ai8rLXeoJTyAjQwELickxK2u2vaubm5NDc3U1xc7FoX5MSJEwwbNoygoCCmTJnCV199xd69exk8eDBTpkwhLCyM3Nxc0tLSKC4uxt/fn8LCwm8l6lN7vy0tLTQ2NlJfX+8a9x0SEkJISAhBQUEc0U6UVmzatAmr1crMmTOZNGkS+/fvdy0F21GCOXToEElJSe74uIRBLF261HCx/Y2hErbWuk4p9UX7t06ttVZK7QVqAZRS/x+wvjvlkd6Ulpbm+l9fX18CAgKAtgkvI0eOpKKigpKSEmJiYrj88ssJCwtj/fr1FBcXc+LECXx9fV29cYvFQlJSEtXV1djtp/+zhg8fTlZWFrW1tQwdOpTDhw+za9eutskxQSZX6aS8vJw1a9aQmprKgAEDOHr0KDt37iQwMNC1BndBQQGDBg3q1r9RNgUW4uIzRMJWSnlprVvhv9uEnfSA0dF+zjzg98DVwFF3tLO9HZjNZiZOnOh6z2w2s3v3burq6ggKCsJut+Pv7+/qRbe0tBAVFeUqgdTU1FBVVcXx48eBtqF9zc3NNDY24nA40Fpjt9tpamrCbreTnZ0NQEBAAElJSWzcuJGqqiqGhk8i2GIhJSmJ8vJy8vLyeP3111m8eDElJSX4+fmRnJxMWFgYhw8fZsiQIRw4cIChQ4fKWiTiNLNnzwbg008/NUxsf9OnE7ZSKlVrnaO1bj05aZ/CAfwRKAeu11q7LVl3xWq14uPjQ3x8PBaLhREjRnDgwOkc+68AACAASURBVAFqamrYs2cPTU1NfO9732PKlClorWlsbCQ3N5eysjJaWlrQWlNTU4PFYnFNe6+srMTb25uoqCgWLlzI7NmzSU5OBmD79u08/vjj2O3NqPoGduzeQVxcHJWVlZSUlLB69WoGDx6Ml5cXFouF8vJyQkND2bNnD7W1tQCMGjXKnR+Z6IOampoMF9vf9NmErZS6HnhfKfVvrfVtZ0jauUALsFBrnX3xW3pmHetcx8TEMHDgQGw2GxkZGZhMJsLCwti4cSOFhYXYbDZCQ0MZMGAAxcXFjBo1CqfTSXV1tas+XVZWRmVlJePGjeO+++5j1qxZ+Pj4uDYv6PD6PieB1y6hYzdI89zHaAQSh9XStPFPZGdnU1ZWhsViYf/+/VRUVBAVFcWQIUPIzs4mIiLCtUmC6L7eHjYqRJ9M2EqpQOB+YDEwWSn1jtZ6QXvS9tZad5RBgoAvgS1a60I3NrlLVquV8vJympqaiIiIwGKxMGbMGEpKSrDb7QwcONA1eaa+vp6SkhKqqqpcsxXNZjNeXl4ArtmOq1atIiQkBIDGxsbT7lljbTntPQBvcwj19fXk5uYSHh5OYGAgwcHBOBwObrzxRkpLS8nLy8PhcDBhwgTXQ86OXXTEmfX2sFEh+mTCbp8IcydQB6wEXjspaXck6zHAlcArHe+di6LKRh5avu2CtrszWrc9bGxtdeCb0YCXl4nm5hbq6wfR2NhI4KQ7GDeuFV9fX+x2GwAWm43WVietrQ5A4efnh3+rg6jmFgICAnhmzWHX9VtbWzGZup9MQ6f+GH+bjYEDB5LV3IxPizdaa0o3l+Hr60NtbTAtheWsyNmByaTw8vLGy8uEn58fLS0OvLxMeHl5y9ZjQrhBn0zYAFrr4vYvG5RSdwPLOpK2UioNSAb+2ZNkfTEp1bbEqtPphclkwul0ApqAgAAaGurx8/UlIMRMYKCZpiYbVquVgAAzdXW1NDU5aBu5+N9rtSXxnktPT+f48SKOHSvE6XTi5dVW9vD19SMoKBCnU9PY2AAKAs2Brkk6Dkcrra1tDzyDg4Nlgo0Huv766w0X29+onu4YcbEppcKBZ4HJgAm4Umtd0tPrpaen6927d1+o5nWL1hqn00lTUxOZmZm899571NbWMmvWLMaMGUNhYSFZWVm0tLRQWFjI2rVrsVqtjBgxgqKiIgoLC/H393eNHoG2ksipNeybnt/YZRs+/t9rACgpKeHXv/41y5cvJyIigqqqKtLS0hgyZAg5OTkMHjyYlJQUoqKiGDBgADNnzqS0tJS6ujrS0tKIi4vrnQ/JIJRSX2ut07s67o6fr4tNKdXjHWd687odfzk/e4dxNz7o6ufLMN0krXWFUmo/MBuYdT7J2l2UUnh5eREUFMSll15KU1MTTqfTtd51x8JMHUP6UlNTOXHiBK2tbc9ZTSYTDQ0N37pmx84z3ZWTk+P6+sEHH+Suu+7ioYceoqysjODgYMrLy0lMTGTo0KEkJyeTmJhIeHg4/v7+hIeHA//dJb7jvicv0aqUknq3EL3EMAlbKTUAmANco7XOdHd7zpfZbOaqq66ioaGBuro6/P39GTlyJC0tLezdu9c129FkMlFXV0dTUxNKKWw2Gw6H44wliWB/L+ptp4+A9MXhSqwdxo4dy7p161i+fDmLFy8mJiaGsLAwgoOD8fLyorW1lZSUFFdMxy+NjpElp+6wExgYeNp9Rf8wffp0ADZs2GCY2P7GMAlba12tlJrbMXHGqKxWK/v37yctLe1ba2Y7nU4CAwNJTU0lPz8fL6+2mndeXh7V1dUArmF2DQ0NhIaGnnbtrKwsnnjiCUpLS6msrGT0bY8DsPuNJa5zvnljCM8++yzDhw93vaeU4oc//CEjRoxg/vz55OTk0NTUxJ49e4iLiyM3N5cf/OAH+Pv7ExwczN69ewkLC3P1rDt+gciKf0L0LkMNtDV6sgbYv38/X3/9Nfv37wf+u0lCYGAgjY2NbNiwgeLiYlJSUkhOTmbMmDGEhYW5FocC2LdvX6fXfumll8jIyGD8+PH84Ac/ID4+nsQhifzpT3/i3Xff5eWXX6ahoYH58+fzxBNPuPaI7DBu3Dg2b97M8OHDOXjwICNGjGD37t1s2bKFVatWUVtby7p16/j888/58MMPKSsrw9/fn9DQUCIjI6UUIkQvM0wPu784eZ0Rh8NBRUUFFouFkpISHA4Hra2t+Pv7o5Ri1KhR7Ny5k8suu4wjR45w7NgxfH19ue+++9i3b99ptetRo0bxxRdfuEobL3xRiNZOpl852nXOihUreOWVV/jggw9YtWoV99xzD7/4xS8YOHAgAIMGDWLt2rX88Y9/5MUXX8Rut3PgwAEsFgvV1dX4+flRU1NDUFAQtbW1eHt7U1tbS3x8/GkPP4UQF5ahetj9Qcc6I2az2bWJwP79+zly5AgAkyZNYsKECQwYMIBvvvmGEydOYLVaGTNmDCaTicbGRg4fPsz69etPu/a1114LwNatW7u8f1hYGI8++iirV69m1qxZvPjii4waNYpnnnnGNQnH39+fRx55hMOHD/OTn/wEh8NBbm4ua9asYcWKFeTn55OXl0dmZiavvvoqW7dupaCgAGh7ANnQ0NA+fFEIcSFJwnaj8PBwYmJiXMPpOurAl1xyCXFxcYwZM4a4uDhiYmIoLCwkKCiI6OhoQkNDWb58+WnXS05OJjw8nO3bt5/13vHx8Tz99NPs2LGDadOm8Zvf/IbJkyeTkZHhOicqKoqXX36Z7du309TUhMPhoKqqivLycnJycli9ejVfffUVu3btwmaz0draSl1dHbm5udTU1Eji7mduueUWbrnlFkPF9jdSEnEjb29voqOjgbZSRFVVFb6+vsTHxxMfH4/D4SAiIoKtW7eilHLtwG632/n3v/9NVVUVDQ0N3yqNpKens3XrVurr610JdOPGnE7v39rayqBBg7jjjju4/PLLeemll5g+fTrf//73mTdvHuPHjwcgNjaW999/n1tuucX1F4K/vz/R0dF4eXkRERFBUVERgYGB+Pr6UldXR35+PoGBgQwaNIjg4GDXPaXObVz33nuv4WL7G+lh9xFms5mwsDAGDRrEoUOHcDgc5OTkkJycTGxsLAEBAcTHx7t2jrHb7axYsYKEhARXgo+Pj+eaa66hoqKifQnXtiVStdadvoqLi7HZbNhsNpKTk3nyyScZO3Ys77zzDo8++ignTpxwtW/48OGsWLECk8nE1q1b2bhxI3v27MFkMmGz2Vzrd5vNZkJDQ/Hx8ZFV1voZq9WK1Wo1VGx/Iwm7j+gYLZKbm8uBAwdYs2YNX3/9Nbt37yYlJYUZM2YQGhrq2n0mLS2NN99887TrTJ06FYBNmzadcxuCgoJ44IEHuPPOO8nJyeHqq6/m888/dx0fOnQoK1aswNfX1zXsLzMzk507d5KVlcXGjRspKCjAYrEwcOBAkpKSCAoK6vFnIvqWOXPmMGfOHEPF9jdSEuljUlNTAYiLiyMnJ4fExEROnDiBxWJxbdKrtSYnJwebzcbGjRtdEwsAV4/7s88+I+0HE7u4S9eUUsyYMYMRI0bwxhtvsHDhQm6++WaWLFlCdHQ0SUlJrFixgvnz51NbW0tTUxP19fV88skn+Pv7k5eXx9ixYxkwYADjxo2TEog4ZwkJCZ3+3CilSEhIID8//+I3qo+QHnYf4+/vT1paGmFhYUycOJHIyEiioqI4cuQIFRUVBAQEEBwcTEpKCiNGjODHP/4xBw8edMV3TIJZv349VVXV3b7vunXr2LZtm+shYWxsLB9//DH/8z//w6pVq5g1axYda2MkJiayYcMGEhMTOXDgAPn5+ZSWlmKz2aitrWXPnj188sknPPXUU2zcuLHTrc2E6Ep+fv63SnfTpk1j2rRpaK1do5E8lSTsPs5kMnHs2DGampoICQkhNDSUlJQUBgwYQGlpKQEBAcyfP5+Skv8urXLPPfeQmprK0aN53RqlkZmZyfLly3nllVf41a9+xYYNG3A4HPj5+bFkyRLWrVtHcHAwN998MytXrgTaekGbN2/mtttuo6mpierqamJiYrjkkku46qqr8PLyYvv27fztb38jKyur1z4fITyJJGwDGDlyJBMnTuTqq68mLi6OpKQkTCaTqzdeW1vL97//ferq6oC25VyfffZZbDY7VVVVZ7y2w+HgnXfeISoqivvvv5+AgAD++te/snjxYv72t7/R1NTE0KFD+fjjjxk7diz3338/f/jDH9BaYzabeeutt1i6dCmlpaVUV1ezY8cOSktLSU1N5bLLLmvbDNhk4vDhw6xdu5bS0lIZ6idEDxm2hq2UUtooa8P2QHNzM4WFhcTHxxMQEEB6ejp2ux1fX182bdpEfHw8NpuNUaNGUV1dzd69e1m4cCHvvvsuvr6+beWUrRWUl5eRnV1MVFTUafc4fvw4a9eupbi4mKlTp1JRUcGECRNISEggJyeHRx99lKeffprLL7+cSy+9lAkTJlBfX89zzz1HdnY2zzzzDP7+/txxxx34+Pjw+OOPM3r0aNfO71OmTGHMmDEcOnSI3NxcbDYbdXV1zJkzp8tFoqTm3XctXLjQcLH9jWETNhAC1Li7Eb2lsLCQ3NxcAJKSkgDw8/MjICDANRY7Pj4eq9VKY2Mj/v7+bNq0if/93//lrbfeQilFSkoKVVWVbNmyhTfeeOO0PRrvuusuMjMzGThwIH5+flRWVgJtPfSYmBiGDx9OVlYWmzZtIiMjgyuvvJKRI0cyYMAA1qxZQ1FREa+88goRERHMnz+f8PBwfv7zn39rqdWamhpqamqIiIhg+PDhXHnllQQEBFzcD1NcEJKw3c+QJRGl1LXAX5RSke5uS2+Jj48nJSWF+Ph413tWqxWLxUJKSgrf+c53sFgspKamcsMNNzB06FCioqL4+9//zj/+8Q8AfHx8iIuPJyMjg1WrVp12j/z8fFpbWxk2bNhpPVulFJGRkUyfPp0ZM2bgcDj44osvKC8vZ/z48bz00ktkZ2dz8803c+jQIQBmzpzJ8uXLqayspLS0lIMHD7J27VqOHj2KyWTiiiuuID8/n9raWvbt20dlZaWUR/q4xMRE1xrnp74SEhK6fZ2KigoqKip61Ibzie1vDJewlVLTgD8Dr2uty852/imxi5RSu5VSu8vLy3ungReIr68vycnJ+Pr6ut4zm82Eh4eTlJREXl7ef5dRHT2aGTNmcNVVV3H55Zfz8MMPu37AwweGc9lll/Hcc899ayLM4cOHKS0tJT4+/qxrWEdGRjJz5kz8/PzYuHEjBw8e5Nprr+Xdd9/F4XDw/e9/n82bNwNw2WWXsXHjRvz8/CgqKqK+vp4BAwZgt9t55513WLNmDUuWLOHTTz/l008/pba2thc+PXGhFBQUdDpaQ2t9TsPr5s2bx7x583rUhvOJ7W8Ml7CBYcDvtNZrlVLRSqnJSqkruxOotV6mtU7XWqdHRET0cjMvPKUUgYGB+Pv7k5yczKRJkxg1ahRpaWmMGDGCKVOm4OvrS2VlJd/97ndxOltRCh577DGam5u59957XTvWdJRHWlo632H9VEFBQVx99dUMHDiQTz/9lNdee42RI0fywQcfEBcXx1133cW7774LwCWXXMLmzZtJTk4mNzeXL7/8ktWrV7Nnzx527drF4cOHXUvEdpRhhBBnZ8SE3QxcrpQaAnwCzAfeVkr90r3Nuniqqqpobm5m0qRJjBw5kurq6radz0tL8fX1JTg4mO3bt3Pw4CG01iQlJfH8889z5MgRnn32WaBtoaj4+HhKSko4218bTU1NWK1W/Pz8mDZtGsOHD+e5555j6dKlhIeH849//IOpU6fy+OOP8+STT9La2kpMTAwbNmxgwoQJlJeXc+zYMQ4fPkxLSwtVVVXExsYSGhrKwIEDsdvt5OXldfuXhxCeqtsJWyl1pVLqK6XUN0qpd5VS43uzYWewC2gCfgC8rbV+gLatw+5WSs12U5suqo5V/iwWCzabjaamJoKCgggKCnLtXDNgwAAqKio4dqwQrTWTJ09m4cKFrFy50rX8akJCAoGBgRw8eLDLZKm1ZsOGDaxZs4aDBw+ilGL27Nncc889vP/++9x9990AvPrqqyxcuJC33nqLG2+8kfr6eoKDg13DAa1WK4WFhWRmZtLc3ExkZCReXl5UV1eTk5NDbm4uhYWFF+0z7A1GKrkJYzqXHvbfgN8A04HlwAtKqYu+5qHW+hvaRod8BxiklApqf+8DwCOGH3h5eREVFYXFYiE2NpbExETS0tLw8fFBKUVwcDAjR47E19eHsrIy3njjDVpbW1m0aBFDhgzh8ccfp66ujpaWFpKSkrDb7WRlZdHQ0OB62Ww29uzZw5YtW6irq8Pb25v9+/fz0UcfkZWVhdVqZdKkSWzevJmZM2fy61//Gq01c+fOZe3atUyYMIEvvviCoqIinn/+eZKSklyLQfn4+HDkyBHXuilDhgxh6NChDBo0yLUka2eLVfV1Ri+5ib7vXIb1VWit17V//ZlSajOwHXj/wjerjVJqGBAG7AacWutWAK31UqVUMzAYeEAp1UBbaeQvvdUWd+pqbLJSCj8/P4KCgjCbzUyePJmysjLy8vLw9vam1M+PkJBQnv/9Ei677DJuvfVW3nrrLaZNm8by5cv55S9/iclk4p133uHDDz9k0aJFjB07FoCbbrrpW/dqbm4GwG63s3btWuLj44mNjWX48OFkZ2fz0UcfMWLECNLS0njppZf41a9+xW233cYLL7zA6NGjefXVV/nJT35CSUmJa0nY+vp6zGYzWVlZ7UMQq/jyyy+ZOXMmgwYN6t0PVZyze+65x3Cx/Y06W89FKbUcyADigDrgaa21QynlA2zTWqf3SsOU+h7wNHC8/bUbeFNrXXfSOVcBycAIYJnW+mBn1+pMenq67lgbw8haW1tpamoiICAAu93O5s2bcTqdHDp0iP+cCKHV2UrjtuXs2LGDNWvWMG3aNH71q1/xwgsv8Nhjj5GWlkZzczMPPfQQNpuN5557jsDAwNMS9skCAgJoamoiMDCQlJQUnE6naz2Ta6+9lqeeeoojR47wwAMPUFFRwf/93/8xYcIEKisruffeeykuLmbChAmMGjUKs9lMSEgIwcHBFBUVUVVVxZw5c7j++utPu29fmlSjlPr6TD/7/eXnSynVp/666U57Hlq+DYBn75h0MZrUK7r6+epOSeSvgJO2nu4NQK5S6j/AIeD0faougPZfBrcCP9ZaXw2sAuKBh5VSIR3naa2/1Fq/Djx0Lsm6PzGZTAQGBmIymSgpKaG1tZX6+noOHjyI1hp/Pz+UUiQlJTFv3jwOHTrE448/TkpKCq+99hpWq9W1T2RVVRVvvPHGWf+DCAsLY9iwYdhsNvbv309dXR2jR4/Gy8uLtWvXsmPHDpKTk3n77bcZPHgwP/vZz9izZw8DBw7k888/JyIigt27d/Pll1+Sm5vL/v37WblyJYcOHXKt8rdt2zbKy8v7VLLwdIWFhT1+zuCu2P7mrAlba71Ra/1/Wus7tdaXAynAYuCxXm6bBRja/vWHwMeAD/B9AKXUOKXUZe3HW3u5LYYQHx9PamoqYWFh+Pr64nS24uPji8ViweFw4OPjw0033YTdbmfZsmWUl5fzzDPPYLfbSU1N5Xvf+x7r16/nlVdeOeN9qqqqOHr0KK2trWitqa6uJiAggNGjRxMcHMwvf/lLioqKCAsL47XXXiMyMpKHH36Y2tpaYmNjWbduHWazmZaWFtdfB1VVVTgcDmJjY8nKymLbtm3s2LFDFq7vQ26//XZuv/12Q8X2N+c8rE9r7dBaH9Bav6O1fqg3GqW1bgGeA76nlJqqtXYCm4G9wJVKqQDgCqC4/XzphtH2MC8pKYlx48YREhKCycuLFkcLMTExJCYmYrFYyM/PZ8GCBUycOJH777+frKwsnnvuORwOB/Pnz+fmm2/myy+/PON97HY7wcHBJCUlMXbsWEaOHAm0TfaZMWMGAA899BB2u52BAwfyu9/9jqqqKn7zm9+gtSYhIYG//OUvHDx4kK+++oqMjAxaWlpwOBzs27ePmpoa12bEZrO51z83IYyiL4/D3gR8DtyulLpSa92qtX4XiAFitNbPa61L3dvEvslqtTJu3Di8TF6YAwJITk5m5syZ3HjjjcTExLBu3Toefvhhpk6dyl133cXu3bt56aWXcDqdzJ8/n0WLFp3x+tHR0QwbNozo6GgCAgK+VVsODg7mySef5ODBg/z+978H2ibS3HffffznP//hrbfeAuCaa67hiSeeoLa2lry8PJqbmzGZTLS2thIaGsr48eOJiIjA4XCQl5eH3W6XTX2Fx+uziz9prW1Kqb8DGliilBoO2IEIoMGtjevjwsLCXCv7OZ2awYmDqays5PDhwyQlJdHY2MhLL72E3W5n1qxZNDY28ve//x2lFPfccw+zZs1i2bJlXV5fKfWtae4ni4qKorS0lEmTJrFixQrsdjujR4/G39/fVc8ePnw4Q4YM4c477yQrK4t//vOfDBkyhEGDBhEeHs6WLVuIjIwkPDycmpoaDh8+TENDAx1D5TqbSt+XHkgK0Vv6bMIG0FpXK6VeB7KAuwEbsEBr3Xm28ECdJarm5mYuu+wyPi06hJ+fH5deGsfWrVupra3F19eXyZMnY7PZeOONN7jlllv429/+xrBhw/j1r3+Nt7c37733Hp999hkHDx5k6dKl+Pr68tRTT5GSkgLA0qVLu1xxr6ysjMjISEaNGkV+fj6fffYZAQEBhIWFMXXqVNcuNh9++CE+Pj4888wzHD9+nG3btmGxWKiqqnIN97vqqqsICAggIiKC2NhY/P39pUQiPFpfLokAoLVu1lqvp21m451a6z3ublNf17FIVHBwEP7+/sTGxnLDDTcwY8YMbrjhBuLi4ti5cyeJiYnMnz+f/Px8Hn74Yf7617+yadMmZsyYQXl5OaNGjeL555/Hy8uLn/3sZ6xatarbozZMJhMzZ87Ex8eHdevW0dzcTFBQEM888wx79uzhxRdfBNrq7suWLSM5OZmMjAzq6uoIDg4mISEBrTVhYWGMGDGCAQMGEBgYKD1pN3rwwQd58MEHDRXb3/T5hN2hvYYtBcxuMJlM2Gw2bDYbNTU1lJSUkJub6yo7JCQkEBQU5HrQ973vfY/6+noWLFjAqlWrOHbsGD//+c/Jy8sjMTGRl19+mUsvvZQ//elPLF26tNtrfpjNZmbOnEltbS1fffUVAN/5zneYN28eL774Irt27QLAYrHw9ttvExAQQHZ2NuXl5ezcuZP169eTk5NDQUEBDoej1z4v0T1z585l7ty5hortbwyTsMW5CQ8Px+FoxW63c+DAAbTWeHt709TUxMSJE5kxYwZTpkwhMjKSQ4cOceONN9LQ0MDMmTNZv75teP3Pf/5zMjIyGDBgAE899RT33XcfGRkZ7Nu3r9vrE3fs83jkyBHXZrxPPvkkgwYN4re//a3rvNjYWJYvX055eTlRUVF89tlnbNmyhby8PI4fP05OTo6MyXaz7OxssrOzDRXb30jC7qe8vb0JDQ3FYglm1qxZREZGUlFRQVZWFuvWrSMuLo5x48Zxxx13kJ6ezpYtW5gxYwZHjhxh9OjRPP/880RGRvLII4+wbt06lFJ897vf5eWXX8bX15fMzEwOHTrUrZ5vQ0MDFosFPz8/oG0kyW233caOHTs4fvy467zRo0fz4IMP8q9//YuqqioKCwvZtm0bgYGBmM1mMjIyePPNN6mu7v5u8OLCufvuu12LfRkltr+RhN2PeXmZCA62EBwcTGJiIjNmzCA6OprGxkYCAwNpaWmhvr6e6dOnM3fuXAoLC5k0aRKrV68mIiKC559/ntGjR/P73/+eN954g+bmZoYMGcKoUaMYPHgwJSUl7Nq164wJ1OFwcPz4cQYPHvyt97/73e+ilOL111//1vsPPvggAwcOZMSIEVitVnbt2sWHH37I4cOHue+++/jXv/7FmjVreuXzEqKvk4TtIUwmE1FRUYSGhlJcXExQUBCXXnop0dHReHt7ExcXx6BBg0hOTmbevHm88cYbBAQE8PTTT7t2l/npT39KZmYmJpOJ5ORkxo4di1KKvXv3cuDAgU5nJRYXF9Pa2npawk5MTOTWW2/lzTff5OjRo673Q0JCWLJkCevWraOsrIzS0lIyMjL4wx/+QGFhIVarlalTp0p5RHgkSdj9nFJ8ax++adOmMWnSJNLS0jCZTJw4cYLjx4/j6+tLQkICRUVFmEwm3nvvPZ555hkuvfRS1qxZw+rVq/Hy8uIXv/gF8fHx/OEPf+CDDz5g165dPPDAA1itVnbv3k1zczPx8fGkpaWRlpaG1WrFx8eHa665hrS0NPLz812v+fPn4+3tzaOPPkp+fj4FBQU0NTXxox/9iCFDhqCUIiwsDB8fH0aOHMnYsWN5+OGHcTgcMmVdeCRJ2B4mJiaGm2++2fUwcOLEiaSnpzNu3DhuuOEGbr75ZlJTUwkNDWX79u2MGzeOLVu2MHv2bPbt28cvf/lLVq5cyaxZs1i1ahX+/v787Gc/48svv2TBggVkZGTw+OOPs3LlShoaGvjmm28YNmwYPj4+p7UlIiKC2267jf+/vTMPb6u88v/nSLIt2/Ie73aUxXHikA0SgrORQoAOS5gZoIG0JVP2lg4ttDAMSzvDrw2dLjC0pfRXpsMyLHlY8wv0AVoCCQNZTEKBLCQ4IXESL4l3K15kWdL5/XEl1XEWstiW5byf59Fj3at7dd57ffW97z3vec9ZvXp1pGQYgNPp5JVXXolkImxra+PgwYOce65VCe7AgQOHVYAPBoNmJqRh2DOkJ84Y+h+bzUZeXh6pqal4PB6SkpJITk5m3759uFwuysvL8fl8rF69mqamJhITE1mwYAE/+clPuOOOO/iP//gPzj33XH70ox/xgx/8gOXLl/PAAw/gdrv5/OQQygAAIABJREFU8Y9/TF5eHqtXr+bdd9/lgw8+oLu7mwULFhy1PYsXL+a1117jt7/97SH+7PHjx/Piiy9y2WWXkZGRwd69ewErnazP5wtNCJpKU1NTpMyYx2Nl3nW5XAN7Ek9T7r///pjbd7hhBPs0xGaz4XK5SEhIoKysjJKSEtrb29m7dy/r16/n888/x+VykZiYSGNjIwsXLuSee+7hnXfe4cUXX6SsrIyXXnqJ559/nl/96ldceuml3HfffVxzzTVkZmayZMkSLrjgAl5//XUqKyuZPHnyUduSmJjI4sWL+c1vfkNlZWWkgALArFmzeOqpp/j6179OTk4OO3bsYOfOnWRnZ+N2uxkxYgQtLS0Eg8FIiTQzE3LguOCCC2Ju3+GGEezTGIfDQWFhYWTZ5/PR0NBAQUEBgUCAsrIy/vKXv7Bnzx4mTZrEqlWrOO+883j00UcpKCjg2muv5cILL+Tuu+/m/vvvZ/Xq1cyaNQuXy0VBQQG33HILqvqlsxP/+te/kpaWxqhRow77bOHChVx66aVs3LiR+vp62tvbSUxMJDc3l0AgQE5ODp2dnaSnp1NcXIzNZsPv99PY2MiIESNwOMwl3l988sknAEybNi1m9h1umKvZECFcGT0zM5OpU6eSnZ3N/v37qaqqwufzMXLkSL744otI/pHS0lLy8vJ48skneeqpp/jlL39JRUUFS5YsYeLEicCXJ2Xas2cPH3zwAddddx1Op/OI2yxatIjXX3+d+Ph4urq66Orq4v3338fv95OQkEBGRgbl5eUUFxcD0NjYSG1tLWBlFjT0D7fffjsAq1evjpl9hxtm0NEAWIN248ePp7y8nGnTplFWVkZnZyc2m42EhAT8fj9FRUWkpKQQCARYtGgR69ZZpZhsNhvXX389r776KklJSfzud7/j5ZdfPq4p7MuWLSM+Pv6YJckuvvhiXC4XgcDf6lR0d3fT2NiI2+0mPT0dsJ4Qmpub8fl8ZGZm4vf7zZT2YYbb7T4k6qn360hPaMMN08M2ANDV1UVnZyfjxo3D6XTi8Xjw+Xw4HA4WLFhAa2srVVVVNDc3U1xczI4dO7juuuv4+c9/zuWXXw7AhAkT+M53vsM777zDqlWr2L59O9deey0FBQURO70n2bS0tPDmm29y4YUXAtDa2kp7+5Ez51588cW8/fbbpKSkICL09PREkkOlp6dTXV3N5s2bOeuss0hKSiIlJQW/34/D4TC97GFEVVXVUT87HRKDxWQPW0SmiUiZiJRFuy2xTO/eSVJSEh6Ph4aGhkhV84aGhkj+j3HjxrFr1y66urpQ1Uia1h/84Ac8//zzFBcXM3r0aG677TZee+01XnjhBfx+P7/4xS949NFHee+992hqaqKhoYHs7GwmTpzI+vXrCQQC3HXXXZxxxhmMGTOGhISEI76uvvpqWltb8Xg8eL1eOjs7qaur49133yUQCPDZZ5+xfv16qqurcbvdFBQUkJeXR1ZW1mCez5tFZKOIbGxoaBg0u4bTh5jrYYvIxcDjwP8DzhORh1T1ySg3K+ax2WxMmDCB5ORk8vPz8Xg8XHXVVTz99NPY7XbWrVvHjh078Pv9iAjl5eVUVlayePFifvzjH/Pqq6/y6KOPUlZm3UMvuugi1qxZw5tvvsnKlSv5+OOPWb58ecReVlYWHR0dXHTRRYwePfpL2zd//nyys7PxeDwkJCQQHx9Pa2srn376KXa7PVI1vqioCIfDgdfrJT09fVAHHVX1caxrkxkzZpipmIZ+J2YEW6znnWTgNuC7qvqaiJQDz4pIgqr+3+P4jpuBm4HDpkobrJqMY8eOBaxwu7q6OiZMmEBlZSVZWVlMmDABh8PB1772NSZMmMCKFSt48803yc7OpqGhgblz5/LDH/6Q++67DxEhOzubJUuWsGTJEgDa2tqoqKjgs88+Y9u2bVRVVfH973//uNrmcDgiU+ZVlYyMDHJycggEAsTHx7Nnzx46Ojqoq6tj3LhxeDwecnJyjvhdJork5HjwwQdjbt/hRsxcraFCu+0ishFIFZE4VV0vItcAL4mIV1Wf+pLvMD2g4yQYDJKUlMSsWbOYNGkSH374IVOnTiU5OZn58+dz8OBBenp6yM/PJy4ujqamJpKSknjooYeor6/n4YcfPkwM09LSmDlzJjNnzjypNgUCAZxOJ3FxcSQkJFBfXx9x5XR1dREIBLDZbHR1dVFbW0tCQgIlJSXY7fZDvsdEkZwcs2fPjrl9hxsxI9i92A8sAF4DelR1o4hcCzwiIu+p6u5j7244HsKDkCNGjKCoqIikpCRUlfb2dmpqasjMzGTmzJmUlZUxduxYnnjiCXbu3EljYyPPPPMMjY2N/PGPf+zXiSwVFRWcffbZfPrppxw4cABVpbu7m56eHjIzM3E6nbS1tbFp0yZUlba2NhoaGiI97c7Ozkg1HiDy13B8rF27Fjg5AY3WvsONmBFsERG1eExEXgB+LyK3Ap2q+oGIbMIq2GvoB8I1GxMTE2loaKCtrY1p06ZRXV2Ny+XC4XAwatQoMjMzyc3NZfTo0TzxxBPs2rWLYDDIW2+9xRVXXMGyZcvIyMg45fa0t7ezefNmnE4nLpeL7u5u4uLiIhEju3fvpqSkhG3btuF2uykqKopU3vF4PDQ1NUXymbhcrsN61sFgMCLoffOUGCzuvfde4OTioaO173BjSAu2iIwHMoGNQBAIAKjq1SKyDHgEWC8iDmA+YIJu+wmbzRapTh7uiaanp1NUVERGRgZer5fc3NzINPeamhr27t1LMBikpKSExMREPv74Y7761a/yr//6r1x66aXEx8cfEkvdl56eHnw+3xE/+/DDDwkGg2RlZZGSkhJZn52djcPhoLa2Fo/Hg9/vp7q6mr1792K320lPT6e+vp6VK1cyZ84cioqKDknNGg4F6+zsPG1zkYwaNYo9e/Yc8TO32z3IrTEciyEr2CJyBfAgUBN6bRSRp1TVA6Cqi0XkeqAAmApcrqrVUWvwMKN3TGvvWOa8vLxIWF8YVaWxsZG8vDz27NlDcXFxZJ3f7+eGG24gLy+PG2+8keuvv56ioqIj2uzu7j7qIOCWLVsAy1UTDjWMj4/H4XDgdDojoX47d+6M+NQ7OzvJyMhg+/btfPTRR7jdbqZMmXLE7w+7bk7HXCR79uwx+cVjhCEp2CISB1wN3KCqa0TkSqAcuFtEfqGqbQCq+kRo+wRV7Y5ei09vOjs7KSoqYv78+fj9fpKTk9m/fz+vvPJKxI0xffp0li5dyoMPPshll13Gt7/9bc4///zjdj9s2LCB0tJSWltbI+t8Ph8+n4/s7GwADh48CFgTeHbu3ElCQgJvv/02Xq+XlpYW2traaGtrIzU19bBJFuEnBYNhKDOUnXWpwLjQ++XAn4A4YDGAiMwUkbNCnx/5OdowKDidTpKSkhg5ciQ7duxgy5Yt9PT0kJiYyJgxY8jNzWXdunUkJiaSnJzM2rVrueSSS5g0aRJLly5l586dx/z+Xbt2sW7dukiK1TBxcXEEAgFqamrYvXs3TU1NiAjp6ekkJCTQ3t5Obm4u+fn5OBwO9uzZQ1VVlSl+YIhZhmQPW1V7RORh4DYR+UJV3xeRD4BC4DIReRqYAywLbW+e56JIZ2cnBw8eRESIj4+no6ODiRMnMmPGDJqbm8nNzaWiogKAzMzMSErUvLw8HnjgAR544AGmT5/OlVdeyaJFiyIZBHt6evj1r3/N0qVLiYuLiySHCveQk5OT6ejowOv14vV6cTqdEfFOSUmhqKgo4uOuqalhzJgxZGZmsmrVKs455xz2799PaWnpUZNOGQ7lkUceibl9hxtDUrBDvA+MB64NRYj8L/B8aPJLgar+Z3SbZ+hLfn4+CxcupKWlhdLSUiZPnszmzZvZvXs3NTU1qCpJSUm0tLSQlJREVVUVeXl5+Hw+gsEg9957L/fddx9z586ltLSUiooKtmzZwt///d+zZ88e6urqyM/Pp7m5Gb/fj9PpjLhUwgOSdrudnTt3YrfbSUlJoaOjg5kzZ5KSksI555zDhg0bWLNmDVVVVZFp60fzaxsO5VTSm0Zr3+HGkBVsVfWKyHNYoXr3iMgEoBvIBo6cIcgQFVwuFzabLZJ0yel0Eh8fT3Nzc0SUx48fz549e/B6vYcMcIXdJC0tLZx55plUVVVx4MAB3n//fYqLi8nLy2Pbtm2kp6eTmpoaKQEWrjxTUFBAY2MjwWCQnp4e0tLSAMufXVdXx4IFC3A4HNTV1WGz2Zg1axbNzc2cd955VFdXU1JSEpVzFousXLkSOLmCAtHad7gxZAUbQFVbROS/gM+AWwAv8E1VPRDdlhnAig4Jxy6HQwCbmpqoq6sDrIgOr9eL3+8nKyuL6upqEhMTKSwspKurC4fDERmg7O7uZvTo0QQCAVSViRMnYrfbEREcDkck98jIkSMpKipi3759eL1eOjo6KCwsxOPxEAgEsNvteDweVBW/389zzz1HYWEhNTU1xMXFMX36dCZMmEB3dzcjR44kGAyaGOzj5Kc//SlwcsIZrX2HG0NasAFU1QesEpH/tRbVVFmNMuEIi46OjkNil4PBIE6nk9zcXEaMGEF6ejpxcXGR91OmTMHr9WK329m+fTt1dXW0tbXx+eefY7PZqKurw+Px0NbWFkkyFR8fj4hEbgIFBQW0trbS3NxMT08PwWAQr9dLd3c3CQkJ2O127HY7EydOpL6+nqqqKmpqapg9ezbJycm0trYSHx9PYWEhPp8vErfd2NjImDFjIrm1DYahyJAX7DCqevQZF4ao0Dd2ubOzk87OTlJTU3E4HDgcDqZMmYLP58PtdjN37lyam5vZsWNHJA/I3r17ycjIwOfzEQgE6OnpweVyYbfbaW9vJxAI4Pf7IzUbm5ubqa6ujhTidTqdZGZmkpycTCAQICUlhVGjRnHVVVdRVVXFihUrGDduHCNHjkRE8Pv9jBo1Cr/fj81m48CBAxw4cIDOzk7y8vKMYBuGNDEj2IahR9/Y5aNNPmlubqatrY3k5GRGjBjBypUr2bBhA7m5ufj9ftxudyTHx9q1a+nu7sblclFdXR3JCRKegl5bW0sgEMDlclFYWEhaWhppaWmkp6dTWVlJS0sL2dnZ1NTUkJGRwQUXXEAwGGT27NmkpqZSXFyM3++PTL3Pzc0lKysLr9d71Ox+BsNQwQi2od/oK+A+n499+/aRn58P/G2KeyAQwOFw4HK5UFXGjRuHy+ViypQpFBQUUFFRQXx8PAkJCbS0tOBwOCLT2sNZ+QoLCykqKiItLY2cnBzS0tIIBoP4/X7sdjsbNmwgPT09kt0vNzeX1NRU3nnnHebNm4fNZiMQCEQSXGVmZkblnBkMJ4IRbMOAsW/fvsikmHCe7f379zN27FhcLhfz5s2jsrISm80W6UVPmjSJ2tpaGhoacDqdFBQUYLPZEBG8Xm8kF0laWholJSV4PB4KCwvJycnB4XAQCATYtm0b27dvp7CwkPz8fOx2O08//TQdHR10dXUBcMkllxwzzWr4ZlNcXHzINPzTmT/84Q8xt+9wwwi2YcAIVzEP/wWrl+33+yMFJMIpUZOTkyktLcXr9TJ27FieffZZamtr8fl8ZGRkUF9fz+TJk3nrrbdoa2uLhPUlJCSgqpxxxhm4XC62bt1KVVUVwWAQt9vNvHnzqKioYMeOHXR0dFBaWkphYSGvv/466enpZGVlHbF3Hb7ZBINB8vPzTQQJMH78+Jjbd7hhBNswYPSuYBPGbrfjcDgiubN37dpFdnY25eXlpKenR0Txzjvv5NVXX0VEeP/998nNzcXhcHDmmWdSUVERiRIZOXIk+/fvZ/PmzeTn53PRRRdRUFDAwYMHGTlyJGvWrKGrq4s5c+aQmZnJlClTWLVqFRUVFWRlZXHRRRfhdrsjqVfr6up48skn+cY3vkFJSQkZGRmnbRa/vrz++usALFy4MGb2HW4YwTYMKiIS6WXHxcXR1dVFRkYGxcXFkcowgUAAr9fL5ZdfzqZNmxARvvjiC9xuN7t378bj8dDd3U17ezuqyhdffEFbWxsjR45kzpw5XH755fj9ftauXRuJ9x4/fjxz5swhEAhQXl5OTk5OJCKkd8/5ySef5K233gKsPMy9Y81Pdx566CHg5IQzWvsON4xgGwYdh8NBeno6wWAQVcVut9Pa2kpubm6kArrX66WoqIji4mI8Hg+TJ09m7969LF68mK1bt3LgwAEqKioYNWoUycnJTJ06FafTSUNDA1988QVNTU08/vjjzJs3j2AwSHp6emQiT21tLePGjaO7u5vt27ezfft2MjMzaW5uZvHixQBcd911gMniZxhaGME2RIVwjzUvL4/m5mYyMzMjE3HCebGDwSDFxcXExcXxxhtvRMqCXXPNNbz99tu43W5aW1vJyckhISEh4sP2+XysWLGCjRs3cuDAAWbNmsXy5cuZMWMGaWlprFu3jhUrVnDTTTcxadIkUlNTWb9+Pc3NzcydOzdS4cRgGGoYwTZElXBxhPb2djweD0lJSZFY6draWpqbmyPhfwBTp04lLS2NuXPnsmzZMkSEtrY25syZQ0tLCxs2bGDixIlkZ2cTDAa55JJLqKmpYevWraxfv56bbropMquyoqKCW265hWAwSHl5Oc3NzZSWlgKHTrs/3QcbYwW3242IMOO6nwEg/zT7kM+qqqqi1LL+wwi2ISr0LcnVe9JNWloagUCA5uZmDh48yGeffYbf76esrAy3201JSQlbtmyJVLgJBAL8+c9/pqamhk8//RS73c4ZZ5zBvHnzmDt3LoFAINJ7b2lp4corryQxMTFSo7K7uxu3231IOayTKRkWyiR5MxCJgjEMHmFBvut/1gGw4YnDS8HFOkawDVGh76zI3jUkgUg+kI0bN5KYmMjEiRMZN25cZKAwIyMDt9tNQkICKSkpTJs2jcbGRlJTUyMujnDekDPPPJMrrriC4uJixo8fT1xcHN/85jcBqzhCfX09XV1dTJo0KfLDPpmSYar6OPA4wIwZM4ZdjvZnnnkm5vYdbhjBNkSF4xnMs9vtnHXWWZEJLD6fD4/Hg8PhID8/n7KyMux2O8nJyRQWFjJ9+nRmz56N2+3G6/WydetW5s+fT319PSUlJYwZMyYy6BgO5QtXYU9NTaWzszNy0xARkpKSjFukF73j6WNl3+GGEWzDkCHcuw3PMiwsLMTv9zN69GhsNlukQG9YQM8++2zq6urYt28flZWVZGdn09LSQmZmJpmZmZFZkx6Ph9bWVgoKCmhoaODAgQOkpKQwduxYMjMzSU9Pjwhz70fnvtkIY5n+qIz+wgsvAHD11VefsP1o7TvciFnBFpF8Va2LdjsM/U94lmFHR0ck/0i4SEJv4YyPj+eSSy5h06ZNlJaWsmbNGjo6OmhqaorMXkxKSjrER36k2Zc2m434+Hh27959yFT0pKQk/H4/Ho8Hp9N51IrusUB/VEb//e9/D5yccEZr3+FGTD7nicjlwF9EZIyIHPevSERuFpGNIrKxoaFhAFtoOBn8fj/79+8nPz+fkpISSktLSU1NPaYfOSkpifLyclpaWlBVUlJSKC4uprOzk/Xr1+NwOCgpKSE1NRX42+zLvvlBwjeJffv2RdbZbDa8Xi/79++nsbFxYA7aYDgBYq7LICIzgKXAHaq660T2He6DQrFO72RMfae0Hwmv10tlZSWlpaUUFxfT0dFBamoqfr+fLVu28NFHHwFQXl4OWDeExsZGRowYcVhvuXfPu/d24R5++K/BEE1iTrABO/CSqq4UETewCKgD9qnqe9FtmuFUOFFxrKysZMuWLYBVSHfSpEkRX3S4sG7vArvHys7XO+/J/v37qa2tJRgMkpqaSk5Ojhl0NAwJYlGw04GLQwV6HwM+BcYAdhFJVtU3oto6w0kTnkRzvIQnuYT/9vZxh10lvel9Q+hbx7Gzs5NNmzYxZcqUyHZOp3PYDDoahgcxJ9iq+mcR+QrwMLBWVR8QkRzgu8DxDXcbhgVOp/OQHnRv+gpy7/zWDoeD5uZmtm/fzoQJE8jMzGTTpk2HuFDy8vIIBoM4HA6T+CnEyy+/HHP7DjeGtGCLiOiRh7ZfAn4ILBKRR1W1XkS6gTIRsWEV6zU+6tMUn8/H9u3bSU1NZcSIEbhcrsOKKTQ1NVFfX092dnYk7Soc6kIxiZ8O5VT8+NHad7gxpAUby1/tDy+IiF1VA6r6VxH5GXAjsEJEXgOuAy4zVdUN+/btY9++fWRlZVFUVAQcHs7ndrux2WyR5SO5UAyH8tRTTwHwrW99K2b2HW4MWcEWkYuBG0TkI6BaVZ9R1YCIOFTVr6pbgNtF5B8BAS5V1cqoNtowKPSODnE6nYd93jtixOv14nK5Dium4HA4IkURDMeHEezoMySvVhGZCfwG+HcgCNwjImWqeq+q+sOiDaCqy6PYVEMU6Bsd0pf4+PhDIkaOxMkkdzIYos2QFGwgHlitqs8BiMhqYG3IpX1fSLTnA2eq6iPRbKhh8OkbHXIkvsz/fDLJnQyxSzj16rE+j4X0q0NVsLuAXBHJVNVmVa0TkVnAGyKySVVfABoBM3x8GnKs6JDjxQwonl58mRjHSvrVITkbQFU/AqqBt3qt2w/8DsgJLW9V1erotNBgMBgGnyHXwxaReFX1qeqtIvInEfkAuCok2COAiSZ0zzAY9I7d7pt7ZKhxrGx8cPwZ+Y7FG2+c/Jy0aO073BhSgi0iNlX1hd7fAzwEXAE8JiI+YCqWeJvQPcOA0zd2eyjTH9n4voxT8fdHa9/hxpAR7JBYB0PvfwHMUdWfAatEpAxIAppUtSqKzTScRhwpFevpzGOPPQbArbfeGjP7DjeGhA+7j1j/CpgMzA9/rqrbVPUjI9aGweRoqVhPV1588UVefPHFmNp3uDEkBLuXWD8ETAQWhkL37NFtmcFgOB0Ih/0d6TVq1KhoNy/CUHKJjATGA5eHxVpVA9Ful8EwVOiPMl+GI3OssL+hFPI3ZARbVfeKyEJVVSPWBsPhDMbAomFoMyRcImHCYXpGrA3DmVGjRh318ftYL9OLjg7Hcpcc6zUQrhQ5Xe/YItIAHD1wdWAYgTVD09iMfZvjVTWl9woRuRm4Ofw58Hk/2TodzufpcIwnYtOtqtl9V562gh0NRGSjqs4wNmPf5mDaM+fT2AwzpFwiBoPBYDg6RrANBoMhRjCCPbg8bmwOG5uDac+cT2MTMD5sg8FgiBlMD9tgMBhiBCPYBoNhSCBDaUrhEMUItsFwkohIQrTbMMxIi3YDBouTvTkZwR4kRCTly7fqd5uTRGR8KD3tYNk8U0Smi8hZg2jzAhH5ymAmCxOR84EbRWTQUvmJSP4g2ZkmImWDfN18FfijiOQMls0+9gdFC0XEGcpOqidzvRrBHgRE5B+Ap0Vk9mA99onIJcAy4IfAEyLyd4Ng8++AZ4GvAc+IyDmDYDMO+BmwFJgpIgOeHyd0nI8Am8IFNwbB5uXAX0RkzEAeo4hcDLwO3Aq8JCLXDZStXjbnA38A/ktV6wfaXsjmBSJyr4j8VESSVTU40L9NEbkMK0pkRahe7Ymn4FBV8xrAFzAO2AGsxhKVckLROQNocwawHTgHEGAJ8OvQe9sA2twKnBtaXgqcDWQN8LEK8FvgPeBPwFfC6wfI3hSgBavyEUAW1nTj0QP8/9wMXDDA59EFvIGVMZPQtboT+PYA/w9vBr4Tep8HzA5fRwNk71JgE/Ad4GlgDZAwwMd4HvApMA/4b+CFk/ke08MeeLqB64B/ABKAq4Fzwo9DA/QYPwJYqqoVal0tu7CEJlIoYgBIBq5X1f8NPbp/F7gNWC4i3xogm4SO7w3g34FXgDtE5HvA7QN0bp3Ai0CeiMwAnsMqZffnATxOO/CSqq4UEbeI3CUi3wz1TPsFtWgHNgKpIhKnquuBa4C7B/J/CPiA6SIyGut/eQ3WE9qd/W2o17X5PVX9var+E9ZNqaS/bfVhAfCKqr4P/BLoEpG7RWTciVynRrAHGFXdC3yiqq1YohLAEu2ZoU1yB8DmW8DKXqs+Abo09AgmInkDYPM9Va0IuSi+DvyLqi4B7gF+LiKT+9tmH25T1SeBGuA/gSQdgKyPqvoh8D9YP/CVwGvADcD1wE9F5Iz+tgmkAxeLyFisR+psYA7wjZDrqz/ZjyUuiQCquhG4FvjnkKAOBBuALuAbwDOq+j3gEuCWkIumP+kAfqeqq0XEHvJdZwHTe280AD7t7cAYEfk+8BbQABRhPflOPN4vMYI9CKhqu4hIqAfzEyzR/qqIPAy8KyIp/eU/C3+Pqtb1Wu0AikIX6LeA/xaRpIHw2alqD/CEqj4eGlxZA6zA6kX1K73a/zawSURmYYnNfwPni8jsfrZnAwgd0wvADar6GBBQ1Q+wfojd/WkzZO/PWC61h4G1qvovwL8BdUC/5Fztdd08hlU/9fcikhbqaX+A5UIYkFl2qroVaAUuB/JFxBVa9zKhG0c/2vIA74QWg6Enzk+ANrDGJ0QkYQCeRNcA64BiYJ2q3qWqt2Gd17uO90uMYPczoaiMWSIS18vtIaqqIQFrU9U7ga8CVwJXq+rB0KP9Kdsk9D/t00PownKL3A18G6v329lfNvu6d1S1JfQ3KCKLgLOA9pO1dTSboXMqqurH8pevAe5U1ZuBV4Hq/rSJ5ecFQFXXYQ3OEWrH1cCZWOf6VOwd7Sb6EtY5XCQiWWoNznUDZSJiO5mb75GuGwBVvTq0/AhwvYh8F6vGqv9EbRzLZm9XgKr+COt8ZgDfC7m1rsHy+54yfWx5Q3/D178/tM1VwGNAwQDY3K2qfwB+A+wVkfCT9Tag5bjdIgPpaD/dXsAVWI8+72A9Nn8PSA19Zuu13WSsR8/Jg2hzHfAZUDYYNrF82rdg/eDOGEiboc9dwPRB/H+G0zo4sAZ1N/XTcTr6LNt7vZ+EJaIfAP+CNZg9fiDOZ2ib64H7sW4WA/5qqXuSAAAFEElEQVQ/DG1zPnATllurP67V0iOdyz7b3I/lx14HTBxIm0A+1hPnL7HE+5MT0YF+u7hP9xcQh/WYPCe0fGXon7IUSOuzbRpQMMg27wUmDLLNxcC4AbaZfoTtTzkS5gSP8x+Bkn6weTGWG+Ae4Npe6/uK+D+GxK/0JO0c97GFPj/lCIqTsOnoB5uXAZ3A873WHSbaWD35bZzkze94bfK3m/w04EYs9+gJ2TQukf4lFSuMD2A5VphZHJZwISJni8hUtdwitYNk8xwRKVXVB1V1+yDaLFPVZaq6Y4BtXhOyOUNEpoHlihlgm+HjnBk6zuWquvNUDInITKwe13KgCrhLRB4EUKsodST2OmTvVVWtPAWTx3Ns4clP/TX+cDy/j7DNUxowFpFk4J+B2wGfiDwLVvnB3udSRFzAu8BFqnpKFYK+zKaG1BrYrqp/VNUfnahNI9j9hFqDbQ8DV4jIvJBofID1yHOuiCQCc4EDg2xzNnAwCjZbB9nmPCw302DanEP/HWc8sFpVn1PVZVhjHItFZGmoPX4RmS8it5+qoRM4ttrQ9qc82HgCv49+samqHVguneeBOwFnLwEN+6ynhbZpVtV9p2LvBG3eKNaMxxMf9D/VRwDzOuRxyIl1h32cXoH/WCP8Y41NY/MYtqZjhQhm9lqXB/wVa2Aa4AygKNaOLZo2e9nIworTfza0PAXLLZMTSzYHfBrv6YSqekXkOazwp3tEZALWSH42/RAlYWwOX5uq+pGIVGOFBs4MrdsvIr8DckLLW/vR3rA+n0ew3SQitwC/FJHPsbwL5+oAToUfCJumgMEAIFZCoDlYURJe4Neq+rGxaWwe7fs1lJNERP6ENVHmqpBg3401seI6QhMS+8tu2DbD7Hx+ie07sMJbL1TVzbFm0wj2ABKKrVQduOngxmaM2wzF5gdD7+8B1mNFgBRiDfZNxRLvfutdH6Udw+J8fom9DKy0Aj9U1U2xaNMItsEQJfqI9S+wQt7mhJbLsGYcNqlqVfRaObwQEaeGJs7Eok0j2AZDFOgj1r/CGlBcqKFoAoPhSJiwPoMhCvQS64ewfNQL1QrdG7QiDIbYwwi2wRAlRGQkMB4r/7Q/lB+l3zMMGoYPxiViMESRUPIqNWJtOB6MYBsMBkOMYFwiBoPBECMYwTYYDIYYwQi2wWAwxAhGsA0GgyFGMIJtMBgMMYIRbIPBYIgRjGAbDAZDjGAE22AwGGIEI9gGg8EQIxjBNhgMhhjBCLbBYDDECEawDQaDIUYwgm0wGAwxghFsg8FgiBGMYBsMBkOMYATbYDAYYgQj2IbDEJGrRGS9iHwqIh+ISHa022QYHojIV0TkmWi3I1Yxgm04EqtUtVxVpwJvA4ui3SDDsGEq8HG0GxGrGME2HIlviciHIvIpcCvgjXaDDMOGaUChiFSIyC4R+Uq0GxRLGME2HIKILAFmAueHetifA1uj2yrDMGIqcFBVzwG+Dfwkyu2JKYxgG/oyGVirqu0iciUwG9gc5TYZhgEiEgeMAB4MrfoktGw4ToxgG/ryFHCriHwInAnsUtWO6DbJMEyYAOxUVV9o+Szg0yi2J+YQVY12GwwGw2mAiFwL/B8s4Y7DGtC+Q1XXR7VhMYQj2g0wGAynDVOBV4G1QCLwEyPWJ4bpYRsMBkOMYHzYBoPBECMYwTYYDIYYwQi2wWAwxAhGsA0GgyFGMIJtMBgMMYIRbIPBYIgRjGAbDAZDjPD/AddRM/4axaK8AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "import emcee\n", "import corner\n", "\n", "# Uniform prior\n", "def lnpriorline(p,amin,amax,bmin,bmax):\n", " a,b = p\n", " if ((amin < a < amax) and (bmin < b < bmax)):\n", " return 0.0\n", " return -np.inf\n", "\n", "# Likelihood\n", "def lnlikeline(p,x,y,yerr):\n", " chisq = chisqstraightline2(p,x,y,yerr)\n", " return -0.5*chisq\n", "\n", "# Posterior\n", "def lnprobline(p,x,y,yerr,amin,amax,bmin,bmax):\n", " lp = lnpriorline(p,amin,amax,bmin,bmax)\n", " if not np.isfinite(lp):\n", " return -np.inf\n", " return lp + lnlikeline(p,x,y,yerr)\n", "\n", "ndim = 2 # Number of parameters\n", "nwalkers = 100 # Number of walkers\n", "a0,b0 = 0.2,1. # Fiducial values\n", "# Priors\n", "amin,amax = 0.,0.5\n", "bmin,bmax = -1.,2.\n", "# Initialize the chain\n", "pmin = np.array([amin,bmin])\n", "pmax = np.array([amax,bmax])\n", "p = [pmin + (pmax-pmin)*np.random.rand(ndim) for i in range(nwalkers)]\n", "# Set up sampler\n", "sampler = emcee.EnsembleSampler(nwalkers,ndim,lnprobline,args=(x,y,yerr,amin,amax,bmin,bmax))\n", "# Burn-in phase (use 300 steps)\n", "p,prob,state = sampler.run_mcmc(p,300)\n", "sampler.reset()\n", "# Perform MCMC (use 700 steps)\n", "p,prob,state = sampler.run_mcmc(p,700)\n", "samples = sampler.flatchain\n", "achain,bchain = samples[:,0],samples[:,1]\n", "# Plot the chains\n", "fig = corner.corner(samples, labels=[\"$a$\", \"$b$\"], range=[[amin,amax], [bmin,bmax]], truths=[a0,b0], quantiles=[0.16, 0.5, 0.84], show_titles=True, labels_args={\"fontsize\": 40})\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Activity\n", "\n", "For the supernova distance-redshift dataset:\n", "\n", "- Run an **MCMC analysis** for parameters $(\\Omega_m,\\Omega_\\Lambda)$\n", "- Determine the **68% confidence regions** of $\\Omega_m$ and $\\Omega_\\Lambda$" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Model selection\n", "\n", "Since Bayesian statistics is related to the probability of models, it allows us to perform **model selection**.\n", "\n", "A typical example is how many model parameters does a dataset justify including in a fit.\n", "\n", "A common approach to this problem is to calculate the **Akaike information criteria** for each model:\n", "\n", "$AIC = \\chi_{\\rm min}^2 + 2p + \\frac{2p(p+1)}{N - p - 1}$\n", "\n", "where $p$ is the number of model parameters, and $N$ is the number of data bins. This formulation penalizes models with more parameters.\n", "\n", "## Activity\n", "\n", "Returning to the same supernova distance-redshift dataset, compute the **Akaike information criteria** for a flat model (where $\\Omega_m + \\Omega_\\Lambda = 1$) and a **curved model** (where $\\Omega_m$ and $\\Omega_\\Lambda$ can take any value). Which model is preferred by this metric?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "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.7.6" } }, "nbformat": 4, "nbformat_minor": 2 }