{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Class 3: Model-fitting\n", "\n", "In this class we will describe the use of the $\\chi^2$ statistic as a hypothesis test of a model, and in determining best-fitting parameters. After these activities you should be able to:\n", "\n", "- apply the $\\chi^2$ statistic as a hypothesis test\n", "- understand the probability distribution of the $\\chi^2$ statistic and its interpretation as a $p$-value\n", "- apply the $\\chi^2$ statistic in parameter fitting\n", "- determine parameter errors and joint confidence regions using intervals of $\\Delta \\chi^2$\n", "\n", "## Comparing data and models\n", "\n", "A key task in statistics is to build models which describe our data. When comparing data and models, we are typically doing one of two things:\n", "\n", "- **Hypothesis testing:** we have a set of $N$ measurements $x_i \\pm \\sigma_i$, which a theorist says should have values $\\mu_i$. _How probable is it that these measurements would have been obtained, if the theory is correct?_\n", "\n", "- **Parameter estimation:** we have a parameterized model which describes the data, such as $y = ax + b$. _What are the best-fitting parameters and errors in those parameters?_\n", "\n", "The most important statistic to help with these tasks is the **$\\chi^2$ statistic** between the data and the model, which is defined (assuming the data points are statistically independent) by:\n", "\n", "$\\chi^2 = \\sum_{i=1}^N \\left( \\frac{x_i - \\mu_i}{\\sigma_i} \\right)^2$\n", "\n", "$\\chi^2$ is a measure of the **goodness-of-fit** of the data to the model, which the above equation evaluates according to **how many standard deviations** each data point lies from the model. If the data are numbers taken as part of a counting experiment, we could use the Poisson error $\\sigma_i^2 = \\mu_i$.\n", "\n", "## Hypothesis testing using $\\chi^2$\n", "\n", "If we run an experiment many times, sampling $N$ data points $x_i$ from a model using Gaussian statistics and evaluating the $\\chi^2$ statistic, what is the resulting probability distribution of $\\chi^2$? The answer is the **$\\chi^2$ probability distribution**,\n", "\n", "$P(\\chi^2) \\propto \\left( \\chi^2 \\right)^{\\frac{\\nu}{2}-1} e^{-\\chi^2/2}$\n", "\n", "where $\\nu$ is the number of **degrees of freedom**. If the model as no free parameters, then $\\nu = N$. If we are fitting a model with $p$ free parameters (see below), then $\\nu = N - p$. We note that the mean of the $\\chi^2$ distribution, i.e. the average over these experiments, is $\\langle \\chi^2 \\rangle = \\nu$. This makes intuitive sense because each data point should lie about 1-$\\sigma$ from the model, hence contribute $1.0$ to the $\\chi^2$ statistic. The variance of the $\\chi^2$ statistic is ${\\rm Var}(\\chi^2) = 2\\nu$, hence if data is drawn from the model we expect $\\chi^2 \\sim \\nu \\pm \\sqrt{2\\nu}$.\n", "\n", "The $\\chi^2$ distribution looks like this for $\\nu = 10$:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(0.0, 0.15)" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAEMCAYAAADj8ECOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3de5iN9f7/8efbMEjJmXaR2UWFShqTQ0o5REJFDDmENmXrx65UdiodtrY2pYOUNpHDjLOmhNg5Vsok0tQIkUPlkEPkPPP5/bEW32kMZsaada+15vW4rs+11roPc7/uazHvue/PfX9uc84hIiJyrgp4HUBERCKDCoqIiASECoqIiASECoqIiASECoqIiASECoqIiAREQa8D5LUyZcq4ypUrex1DcmPtWt/rFVd4m0MkH/rqq692OefK5mSdiC8olStXJjk52esYkhsNG/peFy3yMoVIvmRmP+V0nYgvKBLGBg70OoGI5IAKioSuxo29TiAiOaBOeQldq1b5moiEBR2hSOjq18/3qj4UkbCgIxQREQkIFRQREQkIFRQREQkIFRQREQkIdcpL6Bo82OsEIpIDnhyhmFkzM1trZuvN7Iks5t9kZivN7LiZtc1ifnEz22pmbwQnsXiiXj1fE5GwEPSCYmZRwAigOVAN6GBm1TItthm4D5h0mh/zPLAkrzJKiPjsM18TkbDgxSmvOGC9c+5HADNLBFoD351YwDm3yT8vPfPKZnY9UB6YC8QGIa945Z//9L3qPhSRsODFKa+LgS0ZPm/1TzsrMysADAMePctyPc0s2cySd+7cmeugIiKSfeF2lVdv4CPn3NYzLeScG+Wci3XOxZYtm6PRl0VEJJe8OOW1DaiY4fMl/mnZURdoYGa9gfOBaDM74Jw7pWNfRESCy4uCsgKoYmYx+ApJPNAxOys65+498d7M7gNiVUxEREJD0AuKc+64mfUB5gFRwBjnXIqZPQckO+eSzKw2MBMoCbQ0s2edc9WDnVU8Nny41wlEJAfMOed1hjwVGxvr9MRGEZGcMbOvnHM5upI23DrlJT9ZsMDXRCQsaOgVCV0vvOB71ZMbRcKCjlBERCQgVFBERCQgVFBERCQgVFBERCQg1Ckvoevtt71OICI5oIIioeuKK7xOICI5oFNeEro++MDXRCQs6AhFQtewYb7Xli29zSEi2aIjFBERCQgVFBERCQgVFBERCQgVFBERCQh1ykvoGj/e6wQikgMqKBK6KlY8+zIiEjJ0yktC1+TJviYiYUFHKBK6Ro70vbZv720OEckWHaGIiEhAqKCIiEhAeFJQzKyZma01s/Vm9kQW828ys5VmdtzM2maYXtPMPjezFDP7xsx0LkREJEQEvaCYWRQwAmgOVAM6mFm1TIttBu4DJmWafhDo4pyrDjQDhptZibxNLCIi2eFFp3wcsN459yOAmSUCrYHvTizgnNvkn5eecUXn3A8Z3v9sZjuAssDevI8tQTdtmtcJRCQHvDjldTGwJcPnrf5pOWJmcUA0sCGLeT3NLNnMknfu3JnroOKxMmV8TUTCQlh2ypvZRcB4oJtzLj3zfOfcKOdcrHMutmzZssEPKIExdqyviUhY8KKgbAMy3gJ9iX9atphZcWA28KRzbnmAs0koUUERCSteFJQVQBUzizGzaCAeSMrOiv7lZwLvOed0gl1EJIQEvaA4544DfYB5wPfAFOdcipk9Z2atAMystpltBe4B3jazFP/q7YCbgPvMbJW/1Qz2PoiIyKk8GXrFOfcR8FGmaU9neL8C36mwzOtNACbkeUAREcmxsOyUFxGR0KPBISV0ffTR2ZcRkZChgiKh67zzvE4gIjmgU14Sut5809dEJCyooEjomjLF10QkLKigiIhIQKigiIhIQKigiIhIQKigiIhIQOiyYQldixZ5nUBEckBHKCIiEhAqKBK6hg71NREJCyooEro+/NDXRCQsqKCIiEhAqKCIiEhAqKCIiEhA6LJhCV1Fi3qdQERyQAVFQtecOV4nEJEc0CkvEREJCBUUCV3PP+9rIhIWPCkoZtbMzNaa2XozeyKL+TeZ2UozO25mbTPN62pm6/yta/BSS9D973++JiJhIegFxcyigBFAc6Aa0MHMqmVabDNwHzAp07qlgGeAG4A44BkzK5nXmUVE5Oy8OEKJA9Y75350zh0FEoHWGRdwzm1yzn0DpGda9zZgvnNut3NuDzAfaBaM0CIicmZeFJSLgS0ZPm/1TwvYumbW08ySzSx5586duQ4qIiLZF5Gd8s65Uc65WOdcbNmyZb2OI7lVurSviUhY8OI+lG1AxQyfL/FPy+66DTOtuyggqSRHDh8+TGpqKjt27GDv3r0UL16cMmXKcOWVV3L++ecHZiPTpwfm54hIUHhRUFYAVcwsBl+BiAc6ZnPdecDgDB3xTYEBgY8oWdm9ezdjx47l/fff54svvuDIkSOnLFOgQAGuvfZamjZtSpcuXahWLfP1FiISqcw5F/yNmt0ODAeigDHOuX+Z2XNAsnMuycxqAzOBksBh4FfnXHX/ut2Bf/p/1L+cc++eaVuxsbEuOTk5r3YlX9i5cydPP/00Y8eO5fDhw1x33XXceuut1KlThwoVKlCiRAl+//13tm/fztdff82yZctYsmQJaWlp1KlTh4EDB3L77bdjZjnb8AD/3wovvhj4nRKRMzKzr5xzsTlax4uCEkwqKLmXnp7OW2+9xZNPPsmBAwfo1q0bffr04Zprrjnrujt27GDixIm8/vrrbNy4kbi4OF577TVuuOGG7Ado2ND3qkcBiwRdbgpKRHbKy7nbt28fd999N3//+9+pVasWq1evZtSoUdkqJgDlypXjH//4B2vXruW///0v27Zto27duvTp04f9+/fncXoR8YIKipzihx9+oHbt2syePZvhw4ezYMGCXPeFFCpUiB49evDdd9/x0EMP8eabb1KrVi1WrlwZ4NQi4jUVFPmTlJQUbrrpJvbu3csnn3xC3759c973kYXixYvz6quvsnjxYg4dOkTdunV56623ApBYREKFCoqctHr1aho2bEiBAgVYsmQJDRo0CPg2GjRowOrVq2ncuDEPPvggffv25fjx41kvfMklviYiYUHPQxEANm/eTLNmzShSpAgLFy7k8ssvz7NtlS5dmqSkJB577DFefvllNmzYwNSpUyma+YFaEybkWQYRCTwdoQj79u2jRYsWHDp0iLlz5+ZpMTkhKiqKYcOGMXLkSD766CNatGjBgQMH8ny7IpJ3VFDyufT0dOLj40lNTWX69OlUr149qNt/4IEHmDBhAkuWLKFp06bs27fv/2b26+drIhIWdMornxsyZAhz585l5MiRNGrUyJMMHTt2pHDhwnTo0IFGjRoxf/58SpYsCatWeZJHRHJHRyj52KeffspTTz1F+/bt6dWrl6dZ2rRpw8yZM1mzZg133HEHf/zxh6d5RCTnVFDyqX379tGhQwcqV67MqFGjAnJp8Llq0aIFkyZNYvny5bRt25b0CB/FQSTSqKDkU4899hjbtm1j0qRJFC9e3Os4J7Vp04a33nqLuXPnkpqaikqKSPhQH0o+tHDhQkaNGkX//v2Ji4vzOs4p/va3v/Hbb7+xbMAA9pUrR12vA4lItugIJZ85ePAg999/P5dffjmDBg3yOs5pPf7443zXty/1vv1Wd9SLhAkdoeQzL774Ij/++CMLFy7kvPPO8zrOaZkZw4YNY926dfTp04fLL7+cxo0bex1LRM5ARyj5yKZNm/jPf/5Dx44daXhiaPgQFvXgg8wsW5arrrqKtm3bkpqa6nUkETkDFZR8pH///kRFRTFkyBCvo2TPDz8QvWkTH3zwAdHR0bRs2ZK9e/d6nUpETkMFJZ9YtGgR06ZN44knnuCSMBtwsXLlysyYMYNNmzbRpUsX0tPTvY4kIllQQckHnHM89thjVKxYkUcffdTrOLly4403MmzYMD744AP+/e9/ex1HRLKgTvl84P3332fFihWMGTPm1BF9w8hDDz3E8uXLeeqpp6hduzZNmjTxOpKIZJCrZ8qbWTHgsHMuLfCRAiu/P1M+LS2Na6+9luPHj/Ptt99SsGAY/Q1xYmDI4cNPTjpw4AB16tTh119/ZeXKlVSqVMmjcCKRLc+eKW9mBcyso5nNNrMdQCrwi5l9Z2b/MbMcjXduZs3MbK2ZrTezJ7KYX9jMJvvnf2Fmlf3TC5nZODNbY2bfm9mAnGw3P5o0aRIpKSk8//zz4VVMwFdIMhQTgPPPP58ZM2Zw9OhR2rZty5EjRzwKJyKZZbcPZSFwGTAAqOCcq+icKwfcCCwHhphZp+z8IDOLAkYAzYFqQAczy/zA8h7AHufc5cArwInLku4BCjvnrgauB3qdKDZyquPHj/Pss89Ss2ZN2rRp43WcgKlatSrjxo1jxYoVYdsnJBKJsvsna2Pn3LHME51zu4HpwHQzK5TNnxUHrHfO/QhgZolAa+C7DMu0Bgb5308D3jDf6IUOKGZmBYGiwFHg92xuN9+ZMmUKGzZsYMaMGRQoEIbXX3Ty/42SxZMb77rrLvr168fw4cNp0qQJrVq1CnI4EcksW79lsiomuVnG72JgS4bPW/3TslzGOXcc2AeUxldc/gB+ATYDQ/1F7U/MrKeZJZtZ8s6dO7MZK7Kkp6czePBgqlevTuvWrb2Okztbt/raafz73//muuuuo1u3bmw9w3IiEhxnLShm1sTM3jGzmv7PPfM+1mnFAWnAX4AY4BEz+2vmhZxzo5xzsc652LJlywY7Y0hISkoiJSWFAQMGhOfRSTYULlyYxMREjhw5QqdOnUhLC/lrREQiWnZ+03QH+gOdzOxWoOY5bnMbUDHD50v807Jcxn9660LgN6AjMNc5d8w5twP4FMjRVQj5gXOOwYMH89e//pX27dt7HSdPVa1alREjRrB48WIGDx7sdRyRfC07BWW/c26vc+5RoClQ+xy3uQKoYmYxZhYNxANJmZZJArr637cFPnG+65s3A7fCyUuX6+C74kwyWLp0KStWrKB///7hd2VXLnTp0oV7772XQYMGsWzZMq/jiORb2Skos0+8cc49Abx3Lhv094n0AeYB3wNTnHMpZvacmZ3oWR0NlDaz9cDDwIlLi0cA55tZCr7C9K5z7ptzyROJhg8fTqlSpejSpYvXUc5N3bq+dhZmxsiRI4mJiaFjx47s3n1Kt5qIBEGObmw0s2jn3NE8zBNw+e3Gxo0bN3LZZZfxxBNP5LtTQMnJydSrV4+WLVsybdq0kHissUi4yrMbGzNYbmZX53AdCaI33niDqKgoevfu7XWUoIuNjWXw4MHMmDGDd9991+s4IvlOTgtKL2CimZ1yN5mZzQ1MJMmt/fv389///pd77rkn7EYUzlKbNr6WAw8//DC33HILffv25ccff8yjYCKSlRwVFOfcCuAGoJaZ/c/MepvZ22a2Jqc/SwJv7Nix/P777/Q7MQZWuPvtN1/LgQIFCjB27FiioqLo3LmzLiUWCaIcFQEzewH4FrgG2A48BRQDmjnnmgY+nmRXeno6r776KnXr1iUuLs7rOJ6qVKkSI0aM4LPPPgufh4mJRICcHlXcB9R2ztVwznXEV1guAIaaWfFAh5Ps++ijj9iwYUPkHJ2co44dO9K+fXueeeYZvvrqK6/jiOQL2R1t+MTlMlUzDnXinNvpnGuNb/DIFXmQT7Lp7bffpkKFCtx1111eRwkJZsabb75J+fLl6dSpE4cOHfI6kkjEy/Zow2b2EFAm40Qzi/bfPX8jMCzQ4SR7tmzZwkcffUSPHj0oVCi7Y3SGgUaNfC2XSpUqxdixY0lNTeXxxx8PYDARyUp2b6Nuhm8IlgT/2Fl7gCJAFPAx8Ipz7uu8iShnM2bMGJxz9OjRw+sogfXUU+f8Ixo3bkzfvn159dVXueOOO2jaVF19Inklx09s9A9TXwY45JzbmyepAijSb2xMS0ujcuXKVK9enblzdeV2Vg4dOkRsbCx79uxhzZo1lC5d2utIIiEvL5/Y2NXMdpnZbuC/wIFwKCb5wdy5c9m6dSs9e3o5CHQead7c185R0aJFmTBhArt27eLBBx8kN4+9FpGzy24fylNAE+BKfAM05q8xPULYic74li1beh0l8A4d8rUAuO6663juueeYOnUqE7J4YJeInLvsFpTfnXNfO+d2OOeewvdcEvHY1q1bmT17Nt27d4+szvg80r9/f2688Ub69OnDTz/95HUckYiT3YJykf8piDeZWVlAv71CwJgxY0hPT+f+++/3OkpYiIqK4r333iM9PZ2uXbuSnp7udSSRiJLdgvIMcDXwPLAWqGFmH5nZi2bWIc/SyWmlp6fz7rvv0rhxY2JiYryOEzZiYmJ47bXXWLx4MS+//LLXcUQiSrYuG3bOjcr42cwuwVdgrgFuBxICH03OZNmyZWzatIkXXnjB6yh554478uTH3nfffSQlJfHkk0/StGlTrrnmmjzZjkh+k+PLhsNNpF423KNHD6ZMmcKvv/5KsWLFvI4Tdnbu3MnVV19NuXLlWLFiBYULF/Y6kkhICcbzUCQEHDx4kKlTp3LPPfeomORS2bJlGT16NGvWrGHgwIFexxGJCCooYWjWrFns378//B/xezYNG/paHmnRogW9evVi2LBhLFq0KM+2I5JfqKCEoffee49LL72Um266yesoYW/YsGFcdtlldO3alX379nkdRySsqaCEmZ9//pn58+fTuXNnChTQ13euihUrxoQJE9i2bRsPPfSQ13FEwponv5HMrJmZrTWz9Wb2RBbzC5vZZP/8L8yscoZ515jZ52aWYmZrzKxIMLN7beLEiaSnp9O5c2evo0SMG264gSeffJLx48czdepUr+OIhK2gFxQziwJGAM2BakAHM6uWabEewB7n3OXAK8AQ/7oFgQnAA8656kBD4FiQonvOOce4ceOoW7cuVatW9TpORBk4cCC1a9fmgQce4Oeff/Y6jkhY8uIIJQ5Y75z70Tl3FEgEWmdapjUwzv9+GtDI/5CvpsA3zrnVAM6535xz+eah4V9//TUpKSmR3xl/Qrt2vhYEhQoVYvz48Rw6dIju3btrAEmRXPCioFwMbMnweat/WpbLOOeOA/uA0kBVwJnZPDNbaWaPZbUB/zAxyWaWvHPnzoDvgFfGjRtHdHQ07YL0S9ZzvXv7WpBcccUVDB06lHnz5vHmm28GbbsikSLcenUL4ns65L3+17vM7JRH+jnnRjnnYp1zsWXLlg12xjxx7NgxJk2aRKtWrShVqpTXcYLj4EFfC6IHH3yQZs2a0b9/f1JTU4O6bZFw50VB2QZUzPD5Ev+0LJfx95tcCPyG72hmiXNul3PuIPARUCvPE4eAjz/+mF27duWvzvjbb/e1IDIzxowZQ9GiRencuTPHjuWbLjqRc+ZFQVkBVDGzGDOLBuKBpEzLJAFd/e/bAp8430ntecDVZnaev9DcDHwXpNyeSkxMpESJEjRr1szrKBHvoosuYtSoUSQnJ/Pcc895HUckbAS9oPj7RPrgKw7fA1Occylm9pyZtfIvNhoobWbrgYeBJ/zr7gFexleUVgErnXOzg70PwXbw4EFmzZpFmzZtiI6O9jpOvtCmTRu6du3K4MGDWbZsmddxRMKCBocMA1OnTqVdu3YsWLCARo1O6TKKXCeGXfFoWJTff/+dWrVqcfToUVatWpV/+q5E0OCQESsxMZHy5cvTMA/HtZJTFS9enISEBH755Rf+9re/6VJikbNQQQlx+/btY/bs2bRr146oqCiv4wTXfff5modq167Niy++yIwZMxg1atTZVxDJx7L1gC3xzqxZszhy5AgdOuTDB2N6XExOePjhh5k/fz79+vWjfv361KhRw+tIIiFJRyghLjExkcqVK1OnTh2vowTfrl2+5rECBQowbtw4ihcvTnx8PIcOHfI6kkhIUkEJYTt37mT+/PnEx8fjG3kmn2nb1tdCQIUKFXjvvfdISUnhkUce8TqOSEhSQQlh06ZNIy0tjfj4eK+jCHDbbbfx6KOPMnLkSGbOnOl1HJGQo4ISwhITE7nqqqu45pprvI4ifv/617+4/vrr6dGjB5s3b/Y6jkhIUUEJUVu3bmXp0qV06NAhf57uClHR0dEkJCRw7NgxOnTooKFZRDJQQQlRkydPxjmn010hqEqVKrzzzjt89tlnDBgwwOs4IiFDlw2HqMTERK6//nqqVKnidRTvPPig1wlOKz4+nmXLljFs2DDq1avH3Xff7XUkEc/pCCUErVu3juTk5Px570lG7dv7WogaNmwYcXFxdOvWjfXr13sdR8RzKighKDExESD/PEjrdLZs8bUQVbhwYaZMmULBggVp27at7k+RfE8FJcQ450hISKBBgwZUrFjx7CtEss6dfS2EXXrppUyYMIHVq1fz0EMPeR1HxFMqKCFmzZo1fP/99zrdFUaaN2/OwIEDGT16NO+++67XcUQ8o4ISYhISEoiKiqJtiNwhLtkzaNAgbr31Vnr37s3XX3/tdRwRT6ighBDnHImJiTRu3JiyZct6HUdyICoqioSEBMqUKcOdd97Jjh07vI4kEnQqKCHkyy+/ZNOmTTrdFabKlSvHrFmz2LFjB23btuXo0aNeRxIJKhWUEJKQkEB0dDR33nmn11FCwyOP+FoYuf766xk9ejRLly6lX79+XscRCSrd2Bgi0tLSmDJlCrfffjsXXnih13FCQ8uWXifIlY4dO7J69Wpeeuklrr32Wnr16uV1JJGg0BFKiFi6dCm//PKLTndltHatr4WhwYMH06xZM/r06cPSpUu9jiMSFJ4UFDNrZmZrzWy9mT2RxfzCZjbZP/8LM6ucaX4lMztgZo8GK3NeS0hIoFixYrRo0cLrKKGjVy9fC0MnOuljYmJo06aNRiaWfCHoBcXMooARQHOgGtDBzKplWqwHsMc5dznwCjAk0/yXgTl5nTVYjh07xrRp02jVqhXFihXzOo4ESIkSJUhKSuLIkSO0bt2aAwcOeB1JJE95cYQSB6x3zv3onDsKJAKtMy3TGhjnfz8NaGT+MdzN7E5gI5ASpLx5bsGCBezevVunuyLQlVdeSWJiImvWrCE+Pp7jx497HUkkz3hRUC4GMg7QtNU/LctlnHPHgX1AaTM7H3gcePZMGzCznmaWbGbJO3fuDFjwvJKYmEiJEiVo2rSp11EkDzRv3pw33niD2bNn07dvX5xzXkcSyRPhdpXXIOAV59yBMz10yjk3ChgFEBsbG9L/ew8dOsTMmTO55557KFy4sNdxJI888MADbNiwgaFDh3LZZZfx8MMPex1JJOC8KCjbgIyjHl7in5bVMlvNrCBwIfAbcAPQ1sxeAkoA6WZ22Dn3Rt7Hzhtz5sxh//79epBWVgYO9DpBQA0ZMoSNGzfy6KOPUrlyZT1DRSKOFwVlBVDFzGLwFY54oGOmZZKArsDnQFvgE+c7T9DgxAJmNgg4EM7FBHynu8qVK8ctt9zidZTQ07ix1wkCqkCBAowfP55t27Zx77338sknn1C3bl2vY4kETND7UPx9In2AecD3wBTnXIqZPWdmrfyLjcbXZ7IeeBg45dLiSLB//34++OAD7rnnHgoWDLezj0GwapWvRZCiRYvy/vvvc8kll9CiRQu+/fZbryOJBIxFegdhbGysS05O9jpGliZOnEinTp1YunQpN954o9dxQk/Dhr7XRYu8TJEnNm7cSP369TEzPv30UypXrux1JJE/MbOvnHOxOVlHd8p7KDExkYoVK1KvXj2vo0iQxcTE8PHHH3Pw4EGaNGmi0YklIqigeGT37t3MmzeP9u3bU6CAvob8qEaNGsyePZtt27bRrFkzfv/9d68jiZwT/SbzyPTp0zl27Jiu7srn6tWrx/Tp01mzZg0tW7bkjz/+8DqSSK6poHhk4sSJVK1alVq1ankdRTzWvHlzxo8fz7Jly2jZsiUHDx70OpJIrujSIg/89NNPLF68mOeff54z3aCZ7w0e7HWCoDkxLEuXLl248847SUpKokiRIl7HEskRFRQPTJw4EYB7773X4yQhLp9drNCpUyfS0tLo1q0bd911F7NmzdLoCRJWdMoryJxzjB8/nhtvvJGYmBiv44S2zz7ztXyka9euvPPOO8ydO5c2bdpw5MgRryOJZJuOUIJs5cqVpKam8vbbb3sdJfT985++1wi8D+VMevTowfHjx3nggQdo1aoVM2fO5LzzzvM6lshZ6QglyCZMmEB0dDT33HOP11EkhPXq1YvRo0ezYMECbrvtNvbt2+d1JJGzUkEJouPHj5OQkMAdd9xByZIlvY4jIa579+4kJCSwfPlyGjVqxK5du7yOJHJGKihBNH/+fLZv306nTp28jiJhol27drz//vukpKRw88038/PPP3sdSeS0VFCCaMKECZQsWZLbb7/d6ygSRm6//XbmzJnD5s2bqV+/PqmpqV5HEsmSCkqQ7N27lxkzZtC+fXtdCppdw4f7mtCwYUMWLlzIwYMHqVevHkuXLvU6ksgpVFCCZNKkSRw+fJj777/f6yjho2ZNXxMAYmNjWb58OeXKlaNx48ZMnjzZ60gif6KCEiSjR4/m2muv1VArObFgga/JSTExMXz22WfExcURHx/Pf/7zHz2jXkKGCkoQrFq1ipUrV9KjRw8NtZITL7zga/InpUqVYv78+bRr147HHnuM7t27c/jwYa9jiaigBMPo0aMpXLiwhlqRgClSpAgJCQk8/fTTjB07loYNG+oKMPGcCkoeO3z4MBMnTuSuu+6iVKlSXseRCFKgQAGeffZZpk+fzrfffktsbCxffPGF17EkH1NByWMzZ85kz5496oyXPHP33Xfz+eefU6RIEW666Sbeeecd9auIJ1RQ8tjo0aOJiYnhlltu8TqKRLCrr76aFStWcPPNN9OzZ086derE/v37vY4l+YwnBcXMmpnZWjNbb2ZPZDG/sJlN9s//wswq+6c3MbOvzGyN//XWYGfPiQ0bNvC///2Pbt266TG/ufH2274m2VK6dGnmzJnD888/T2JiIrGxsaxatcrrWJKPBP23nJlFASOA5kA1oIOZVcu0WA9gj3PucuAVYIh/+i6gpXPuaqArMD44qXNn5MiRFCxYkB49engdJTxdcYWvSbZFRenhgpEAAA80SURBVEUxcOBAPvnkEw4cOECdOnUYOXKkToFJUHjxZ3McsN4596Nz7iiQCLTOtExrYJz//TSgkZmZc+5r59yJS1lSgKJmFpK3nR88eJAxY8Zw11138Ze//MXrOOHpgw98TXLs5ptvZtWqVdx666307t2bFi1a6CowyXNeFJSLgS0ZPm/1T8tyGefccWAfUDrTMm2Alc65U55AZGY9zSzZzJJ37twZsOA5kZCQwJ49e+jTp48n248Iw4b5muRK2bJl+fDDD3n99ddZtGgRNWrUICEhQUcrkmfC8sS+mVXHdxqsV1bznXOjnHOxzrnYsmXLBjecb/u88cYb1KhRgwYNGgR9+yInFChQgD59+rBq1SquuOIKOnbsSPv27fHqDy2JbF4UlG1AxQyfL/FPy3IZMysIXAj85v98CTAT6OKc25DnaXPh888/Z9WqVfTp00d3xktIqFq1KkuXLmXw4MHMmjWLq666infffVdHKxJQXhSUFUAVM4sxs2ggHkjKtEwSvk53gLbAJ845Z2YlgNnAE865T4OWOIeGDRtGiRIldGe8hJSCBQsyYMAAvv76a6688kq6d+/OLbfcouHwJWCCXlD8fSJ9gHnA98AU51yKmT1nZq38i40GSpvZeuBh4MSlxX2Ay4GnzWyVv5UL8i6c0bp165g5cya9e/fm/PPP9zqOyCmqV6/OkiVLGDVqFN988w3XXHMNTz31FH/88YfX0STMWaQf8sbGxrrk5OSgbe+BBx7g3Xff5aeffqJChQpB225E2uK/dqNixTMvJ7m2Y8cOHnnkESZMmMBf/vIXXnzxRTp16qT7pgQz+8o5F5uTdfSvJoB27NjB2LFj6dKli4pJIFSsqGKSx8qVK8f48eNZtmwZF198MV27duWGG25g2bJlXkeTMKSCEkCvv/46R48e5dFHH/U6SmSYPNnXJM/Vr1+f5cuXM378eH755RcaNGjAnXfeyTfffON1NAkjKigBsnv3bl599VXuvvturtDd3YExcqSvSVAUKFCATp06sXbtWp5//nkWLVrEtddeS3x8vDruJVtUUAJk2LBhHDhwgEGDBnkdReScFCtWjIEDB7Jx40aefPJJPvzwQ6pXr07nzp1JSUnxOp6EMBWUANi1axevvfYa7dq1o0aNGl7HEQmIkiVL8sILL7Bx40b+8Y9/MGPGDGrUqEHLli3VxyJZUkEJgKFDh/LHH3/w9NNPex1FJODKli3L0KFD2bx5M4MGDeLzzz+nQYMG1K9fn+nTp3P8+HGvI0qIUEE5Rz/99BOvvvoqHTt2pFq1zIMmi0SO0qVL88wzz/DTTz/x2muv8fPPP9O2bVtiYmJ44YUX2L59u9cRxWO6D+UcxcfHk5SURGpqKpUqVcqz7eRLu3b5XsuU8TaHZCktLY3Zs2czYsQIPv74YwoVKkSbNm3o1q0bjRo1IioqyuuIcg50H0qQLV26lMmTJ9O/f38Vk7xQpoyKSQiLioqiVatWzJs3j9TUVB588EHmzZvHbbfdxqWXXsqAAQN0dVg+oyOUXEpLSyMuLo7t27ezdu1aihUrFvBt5Htjx/pe77vPyxSSA0eOHOGDDz5g7NixzJ07l7S0NG644QY6dOjA3XffTUXdqBo2dIQSRK+88gorV65k6NChKiZ5ZezY/ysqEhYKFy5M27Zt+fDDD9m6dStDhw7l0KFD9OvXj0qVKlGnTh2GDh3Kxo0bvY4qeUBHKLmQmppKzZo1adasGTNnztQQ9XmlYUPf66JFXqaQAPjhhx+YPn0606dP56uvvgLguuuuo3nz5jRv3pw6depQsGBBj1NKRrk5QlFByaG0tDTq16/PunXrSElJ0ZhdeUkFJSJt3LiRGTNm8P777/PZZ5+RlpbGhRdeSJMmTWjevDlNmjTRqbEQkJuCoj8Jcuipp57iiy++YNKkSSomIrkQExPDI488wiOPPMLevXtZsGABc+bMYc6cOUybNu3kMjfffDMNGzbk5ptvpnLlyt6GlmzREUoOTJ06lXbt2tGzZ0/efvvtgPxMOQMdoeQrzjm++eYbFi5cyOLFi1myZAm7d+8GONn/EhcXR1xcHLVq1VLfZR7TKa8sBKqgrF69mnr16lGzZk0++eQTChcuHIB0ckYHD/pezzvP2xziifT0dFJSUli8eDFLly7lyy+/ZNOmTYBvIMvq1asTFxdHbGwsV199NTVq1ODCCy/0NnQEUUHJQiAKyrfffsutt95KdHQ0K1as4KKLLgpQOhHJiR07drBixQq+/PLLk+3EUQz4jmSuvvrqk+2qq66iSpUqenpqLqigZOFcC8qaNWto1KgRhQoVYtGiRVSpUiWA6eSM3nzT99q7t7c5JGQ559i8eTNr1qz5U0tNTf3TGGMXXXQRVapUOdmqVq1KlSpVuPTSS7ngggs83IPQpYKShXMpKNOmTaN79+5ccMEFKiZeUB+K5NLRo0dZu3Yta9euZd26daxbt44ffviBdevWsWPHjj8tW7JkSSpVqkTFihWpVKnSyVaxYkUuuugiypcvny+PcHSVV4Ds3buXgQMHMmLECOrUqcPkyZM1tIpIGImOjj552iuzffv2nSwymzdvZvPmzWzZsoXNmzfz6aefsmfPnlPWKVasGBUqVKB8+fJUqFDh5Pvy5ctTunRpSpUq9adWrFixfHl/micFxcyaAa8CUcB/nXP/zjS/MPAecD3wG9DeObfJP28A0ANIA/6fc25eoHLt27ePMWPG8MILL7Bnzx769evHkCFDiI6ODtQmRMRjF154IbGxscTGZv3H9/79+9myZQtbtmzh119/Zfv27fz6668nW2pqKosWLfpT301mBQsW/FOBKVmyJMWLF+eCCy7g/PPP54ILLjhrO++88yhatCgFCxYMm+IU9IJiZlHACKAJsBVYYWZJzrnvMizWA9jjnLvczOKBIUB7M6sGxAPVgb8AC8ysqnMuLTdZjh49SmpqKl9++SUff/wxSUlJHDlyhCZNmvDSSy9Rs2bNc9lVEQlDF1xwAdWqVTvr4yiOHj3Kzp072b1791nbtm3bSE1NZf/+/ezfv59Dhw5lO0+BAgUoWrQoRYsWpUiRIlm+z/i5SJEiFCpUiEKFChEdHX3yfebPZ5uXG14cocQB651zPwKYWSLQGshYUFoDg/zvpwFvmK9EtwYSnXNHgI1mtt7/8z4/3cbWrl1LvXr1OH78+Ml27Ngxdu3axa4Tw6MD5cqVo2fPnnTq1Im4uLjA7a2IRKTo6GguvvhiLr744hyvm5aWxoEDB04WmKzaoUOHTrbDhw+f9v1vv/32p+mHDx/m2LFjJ1swH4DmRUG5GNiS4fNW4IbTLeOcO25m+4DS/unLM617yrdpZj2Bnv6PRz7//PNvzxZqx44dvP7667z++uvZ3Y9QUQbYddalwlcZzCJ7/yL3+4vkfYPI378rcrpCRHbKO+dGAaMAzCw5p1cqhBPtX3iL5P2L5H2D/LF/OV3Hi+HrtwEZR367xD8ty2XMrCBwIb7O+eysKyIiHvCioKwAqphZjJlF4+tkT8q0TBLQ1f++LfCJ890wkwTEm1lhM4sBqgBfBim3iIicQdBPefn7RPoA8/BdNjzGOZdiZs8Byc65JGA0MN7f6b4bX9HBv9wUfB34x4G/Z+MKr1F5tS8hQvsX3iJ5/yJ530D7d4qIv1NeRESCQ48AFhGRgFBBERGRgIjogmJmzcxsrZmtN7MnvM4TaGa2yczWmNmq3FziF2rMbIyZ7TCzbzNMK2Vm881snf+1pJcZc+s0+zbIzLb5v79VZna7lxnPhZlVNLOFZvadmaWYWV//9Ej5/k63fxHxHZpZETP70sxW+/fvWf/0GDP7wv87dLL/QqrT/5xI7UPxD/HyAxmGeAE6ZBriJayZ2SYg1jkXETdXmdlNwAHgPedcDf+0l4Ddzrl/+/8oKOmce9zLnLlxmn0bBBxwzg31MlsgmNlFwEXOuZVmdgHwFXAncB+R8f2dbv/aEQHfoX8kkmLOuQNmVghYBvQFHgZmOOcSzewtYLVzbuTpfk4kH6GcHOLFOXcUODHEi4Qo59wSfFf1ZdQaGOd/Pw7ff+Kwc5p9ixjOuV+ccyv97/cD3+MbxSJSvr/T7V9EcD4H/B8L+ZsDbsU3/BVk4/uL5IKS1RAvEfMPwM8BH5vZV/7hZiJReefcL/73vwLlvQyTB/qY2Tf+U2JheTooMzOrDFwHfEEEfn+Z9g8i5Ds0sygzWwXsAOYDG4C9zrkTg4Gd9XdoJBeU/OBG51wtoDnwd/9plYjlv7k1ks7RjgQuA2oCvwDDvI1z7szsfGA60M8593vGeZHw/WWxfxHzHTrn0pxzNfGNQBIHXJnTnxHJBSXih2lxzm3zv+4AZuL7RxBptvvPX584j73jLMuHDefcdv9/4nTgHcL8+/Ofe58OTHTOzfBPjpjvL6v9i7TvEMA5txdYCNQFSviHv4Js/A6N5IKSnSFewpaZFfN3DmJmxYCmwFlHVQ5DGYfh6Qq872GWgDrxi9bvLsL4+/N36o4GvnfOvZxhVkR8f6fbv0j5Ds2srJmV8L8viu9ipu/xFZa2/sXO+v1F7FVeAP5L+Ibzf0O8/MvjSAFjZn/Fd1QCviF0JoX7/plZAtAQ37Dg24FngFnAFKAS8BPQzjkXdp3bp9m3hvhOlThgE9ArQ39DWDGzG4GlwBog3T/5n/j6GSLh+zvd/nUgAr5DM7sGX6d7FL4DjSnOuef8v2cSgVLA10An//Oosv45kVxQREQkeCL5lJeIiASRCoqIiASECoqIiASECoqIiASECoqIiASECoqIiARE0B8BLCL/x8zuBFoAxYHRzrmPPY4kkmu6D0UkBPgHFRzqnOvhdRaR3NIpL5HQMBAY4XUIkXOhgiISBP6x17aZ2WD/59r+J/wVNbMhwJwTz9sQCVc65SUSJGZWGkgGquMb46ozcBO+QfdWAKucc295l1Dk3KigiASRmaUA3wBfO+de8jqPSCDplJdIcH0DVADC+hnkIllRQREJEjMrC9wCTPc/kEkkouiUl0iQmFkScAD4wzn3N6/ziASajlBEgsDMegGHgMfxPVpVJOLoCEUkj5lZFeADoK5zbo+ZzQecc66px9FEAkoFRUREAkKnvEREJCBUUEREJCBUUEREJCBUUEREJCBUUEREJCBUUEREJCBUUEREJCBUUEREJCBUUEREJCD+P+WOM7XmTqHOAAAAAElFTkSuQmCC\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", "%matplotlib inline\n", "\n", "# Plot the chi-squared probability distribution for a given number of degrees of freedom\n", "nu = 10 # number of degrees of freedom\n", "xmin,xmax,nx = 0.,30.,301 # range of values for chi-squared\n", "x = np.linspace(xmin,xmax,nx) # variable for chi-squared\n", "y = stats.chi2.pdf(x,nu) # evaluate probability for each value of chi-squared\n", "# Make the plot\n", "ymin,ymax = 0.,0.15\n", "plt.plot(x,y,color='black')\n", "plt.plot([nu,nu],[ymin,ymax],color='red',linestyle='dashed')\n", "plt.xlabel(r'$\\chi^2$')\n", "plt.ylabel(r'$P(\\chi^2)$')\n", "plt.xlim(xmin,xmax)\n", "plt.ylim(ymin,ymax)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To use $\\chi^2$ as a hypothesis test, we compute $\\chi^2$ for a given dataset and model. We then ask the question: **if the model is correct, what is the probability that this value of $\\chi^2$, or a larger one, could arise by chance?** This probability is called the $p$-value and may be calculated from $P(\\chi^2)$ by integrating the area under the curve at values greater than some value of $\\chi^2$.\n", "\n", "For example, suppose that $\\nu = 30$ and consider datasets with two different values, $\\chi^2 = 37.5$ and $\\chi^2 = 52.1$. We would integrate the following areas:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(0.0, 0.06)" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAEQCAYAAACX5IJuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3dd3wUdf7H8dcnoYRQpUqJgoYiCHIQjaIiTQSUoiJFutxhOU7xLKfeySH4886GeIonKCdFRDhEqULARLAASlUBaRaK9A4hhCSf3x+z8ZIYYAO7O7Obz/Px+D52d3Z25p1hyScz35nviKpijDHGXKgotwMYY4yJDFZQjDHGBIQVFGOMMQFhBcUYY0xAWEExxhgTEFZQjDHGBIQrBUVE2onIRhHZIiJP5PN+cRGZ6nt/uYjUzPFeIxFZKiLrRORbEYkJZXZjjDH5C3lBEZFoYDTQHqgP9BSR+nlmGwgcUtV44BXged9niwDvAvepagOgBXA6RNGNMcachRt7KNcAW1T1B1VNB94HOueZpzMwwfd8OtBaRARoC3yjqmsBVPWAqmaGKLcxxpizcKOgVAe253i9wzct33lUNQM4AlQA6gAqIgtEZJWIPB6CvMYYY/xQxO0ABVQEuAG4GkgFPhGRlar6Sc6ZRGQQMAigZMmSTevVqxfyoMYE1MaNzmPduu7m8BjbLMGzcuXK/apaqSCfcaOg7ATicryu4ZuW3zw7fP0mZYEDOHszS1R1P4CIzAOaALkKiqqOBcYCJCQk6IoVK4LwYxgTQi1aOI+ffupmCs+xzRI8IvJzQT/jRkH5GqgtIrVwCkcP4O4888wC+gFLga5AsqqqiCwAHheRWCAduAmn096YyPa3v7mdwJNss3hLyAuKqmaIyGBgARAN/EdV14nIcGCFqs4CxgGTRGQLcBCn6KCqh0RkJE5RUmCeqs4N9c9gTMi1aeN2Ak+yzeItEunD19shLxMR1qxxHhs3djeHx9hmCR5f/3RCQT4Tbp3yxhROQ4Y4j9ZZkIttFm+xoVeMMcYEhBUUY4wxAWEFxRhjTEBYQTHGGBMQ1ilvTDh47jm3E3iSbRZvsYJiTDho1sztBJ5km8Vb7JCXMeHgyy+dZnKxzeIttodiTDh46inn0S64yMU2i7fYHooxxpiAsIJijDEmIKygGGOMCQgrKMYYYwLCOuWNCQejRrmdwJNss3iLFRRjwoGNz54v2yzeYoe8jAkHixY5zeRim8VbbA/FmHDw7LPOo92iMBfbLN5ieyjGGGMCwgqKMcaYgLBDXqZQSE9P58MPPyQpKYktW7aQkZFBfHw8LVu2pGvXrpQqVcrtiMaEPdtDMRFvypQpxMfH06NHDz788ENUlaJFi7JgwQIGDBhA9erVefnll0lPT3c7qjFhzQqKiVhpaWn07t2bu+++m4svvph58+axb98+lixZwqeffsquXbv4/PPPueGGG3j00Ue54YYb2L59u9ux8zdmjNNMLrZZvEVU1e0MQZWQkKArVqxwO4YJsWPHjtG5c2dSUlIYPnw4Tz75JEWKnPkI7/Tp07nnnnuIiYlh4cKFXHXVVSFMa4z3iMhKVU0oyGdsD8VEnNOnT9O1a1eWLFnCu+++y9NPP33WYgLQtWtXvvrqK4oXL07Lli1Zs2ZNiNL6afZsp5lcbLN4ixUUE3H+9Kc/kZSUxFtvvUWvXr38/ly9evVYsmQJpUqVokOHDmzbti2IKQvo5ZedZnKxzeItVlBMRJkyZQpjxozhiSeeYMCAAQX+fK1atZg3bx4nTpygY8eOnDx5MggpjYlMVlBMxPj555+57777aNasGSNGjDjv5Vx55ZVMnTqVb775hiFDhgQwoTGRzZWCIiLtRGSjiGwRkSfyeb+4iEz1vb9cRGr6ptcUkZMissbX3gx1duNdf/rTn8jMzGTy5Mnn7DM5l3bt2vHEE08wduxY3n///QAlNCayhbygiEg0MBpoD9QHeopI/TyzDQQOqWo88ArwfI73tqpqY1+7LyShjefNnDmT2bNnM2zYMGrWrBmQZQ4fPpxmzZpx//33s3v37oAs05hIFvLThkXkOmCYqt7ie/0kgKr+I8c8C3zzLBWRIsBuoBJwKTBHVa/0d3122nDkS01N5YorrqBMmTKsWrWKokWLBmzZmzZtolGjRnTq1Ilp06YFbLkFln19TFycexk8yDZL8ITLacPVgZxXj+3wTct3HlXNAI4AFXzv1RKR1SKyWERuDHZY432vvfYa27ZtY/To0QEtJgB16tRh6NCh/Pe//2XmzJkBXXaBxMXZb8182GbxlnDrlN8FXKKqvwP+DLwnImXyziQig0RkhYis2LdvX8hDmtA5cuQIzz//PB06dKB58+ZBWcdjjz1Gw4YNefDBB90762vqVKeZXGyzeIsbBWUnkPNvihq+afnO4zvkVRY4oKqnVPUAgKquBLYCdfKuQFXHqmqCqiZUqlQpCD+C8YqRI0dy6NChCzqr61yKFi3Kq6++yrZt23jllVeCtp6z+ve/nWZysc3iLW4UlK+B2iJSS0SKAT2AWXnmmQX08z3vCiSrqopIJV+nPiJyGVAb+CFEuY3HHDx4kJEjR3LnnXfSpEmToK6rZcuW3H777Tz33HPs2rUrqOsyJlyFvKD4+kQGAwuADcA0VV0nIsNFpJNvtnFABRHZgnNoK/vU4ubANyKyBpgO3KeqB0P7ExivePPNNzl+/DhDhw4NyfpeeOEF0tPTQ7Y+Y8KNDQ5pwtKpU6eoWbMmV111FfPnzw/Zeh988EHeeOMNvv/+e+Lj40O2Xlq0cB4//TR06wwDtlmCJ1zO8jLmgr333nvs3r2bRx99NKTrffLJJylWrBjDhw8P6XqNCQe2h2LCjqrSsGFDihQpwurVqxGRkK7/scceY+TIkaxbt4569eqFZqX79zuPFSuGZn1hwjZL8NgeiikUPvvsM9atW8eQIUNCXkwAHn/8cUqUKMGwYcNCt9KKFe23Zj5ss3iLFRQTdsaMGUPZsmXp1q2bK+uvVKkSDz30EFOnTuX7778PzUrHj3eaycU2i7dYQTFh5cCBA0yfPp0+ffoQGxvrWo6HHnqImJgYXnzxxdCs0H5z5ss2i7dYQTFhZeLEiaSnpzNo0CBXc1SuXJmBAwcyadIkdu7Me12uMYWTFRQTNlSVsWPHcu2119KwYUO34/DII4+QlZXFqFGj3I5ijCdYQTFh4/PPP+f77793fe8kW61atejWrRtjxozh8OHDbscxxnVWUEzYmDhxIiVLlnStMz4/jz/+OMeOHWPMmDFuRzHGdXYdigkLp06d4uKLL+a2225j0qRJbsfJpVWrVmzdupWtW7de8J0izyg11Xl08UQEL7LNEjx2HYqJWB9//DGHDx+mV69ebkf5jQcffJBt27Yxa1beMU4DKDbWfmvmwzaLt1hBMWFh8uTJVKpUiTZt2rgd5Tc6duzIpZdeyr/+9a/greSNN5xmcrHN4i1WUIznHT16lNmzZ9O9e/fgHVK6ANHR0QwePJjFixezdu3a4Kxk2jSnmVxss3iLFRTjeTNmzODUqVOePNyVbeDAgcTGxvLaa6+5HcUY11hBMZ733nvvcdlll5GYmOh2lDO66KKL6NOnD5MnT2Z/9oiFxhQyVlCMpx08eJDk5GS6devmykCQBTF48GDS0tIYb2OBmELKCorxtFmzZpGZmcmdd97pdpRzuvLKK7n++usZO3YskX46vjH5sYJiPG3GjBnExcXRtGlTt6P4ZdCgQWzevJnFixcHdsGffmq3JcyHbRZvsYJiPOvYsWMkJSVxxx13eP5wV7a77rqLcuXK2ZXzplCygmI8a968eZw6dSosDndlK1GiBH379mXGjBns27cvcAt+6SWnmVxss3iLFRTjWTNmzKBy5co0a9bM7SgFMmjQINLT05k4cWLgFjpnjtNMLrZZvMUKivGktLQ05s6dS5cuXYiOjnY7ToE0aNDAOudNoWQFxXjSJ598wokTJ7j99tvdjnJeBg0axKZNmwLfOW+Mh1lBMZ40d+5cSpYsSYsWLdyOcl6yO+ffeustt6MYEzJWUIznqCpz5szh5ptvJiYmxu0456VEiRLcfffdzJgxIzA33ypRwmkmF9ss3mIFxXjOt99+y/bt27n11lvdjnJB+vfvT1paGtMCMXrhxx87zeRim8VbrKAYz5k7dy4AHTp0cDnJhUlISKBBgwa88847bkcxJiSsoBjPmTNnDk2bNqVatWpuR7kgIsKAAQNYtmwZGzZsuLCFjRjhNJOLbRZvcaWgiEg7EdkoIltE5Il83i8uIlN97y8XkZp53r9ERI6LyKOhymxCY//+/SxdujTsD3dl69WrF9HR0UyYMOHCFvTJJ04zudhm8ZaQFxQRiQZGA+2B+kBPEamfZ7aBwCFVjQdeAZ7P8/5IwI6cRqD58+ejqtx2221uRwmIiy++mA4dOjBx4kQyMjLcjmNMULmxh3INsEVVf1DVdOB9oHOeeToD2X/STQdai28wJxHpAvwIrAtRXhNCc+bMoUqVKmEzGKQ/BgwYwK5du0hKSnI7ijFB5UZBqQ5sz/F6h29avvOoagZwBKggIqWAvwDPnG0FIjJIRFaIyIqAjqdkgiojI4P58+fToUMHoqIip3vv1ltvpWLFinafFBPxwu1/7TDgFVU9fraZVHWsqiaoakKlSpVCk8xcsOXLl3PkyBHat2/vdpSAKlasGL169WLmzJkcPHjw/BZSoYLTTC62WbzFjYKyE4jL8bqGb1q+84hIEaAscABIBF4QkZ+AIcBTIjI42IFNaCQlJREVFUXr1q3djhJwAwYMID09nSlTppzfAj74wGkmF9ss3uJGQfkaqC0itUSkGNADmJVnnllAP9/zrkCyOm5U1ZqqWhMYBTynqq+HKrgJrqSkJK6++mrKly/vdpSAu+qqq2jcuLFdk2IiWsgLiq9PZDCwANgATFPVdSIyXEQ6+WYbh9NnsgX4M/CbU4tNZDl06BBfffUVbdu2dTtK0PTv35+VK1eyfv36gn/4ySedZnKxzeItrvShqOo8Va2jqper6v/5pg1V1Vm+52mqepeqxqvqNar6Qz7LGKaqdmudCJGSkkJWVlZEF5SePXsSHR3NpEmTCv7hpUudZnKxzeIt4dYpbyJUUlISpUuXJjEx0e0oQVO5cmVuueUWJk+eTFZWlttxjAk4KyjGdarKggULaNWqFUWLFnU7TlD16dOH7du3231STESygmJct3XrVn766aeIPtyVrXPnzpQuXfr8DnsZ43FWUIzrsq8gLwwFpUSJEnTt2pXp06eTmprq/wdr1HCaycU2i7dYQTGuS0pKolatWlx++eVuRwmJPn36cOzYMWbNynu2/Fm8+67TTC62WbzFCopx1enTp0lOTqZt27b4hmuLeDfddBNxcXF22MtEHCsoxlXLly/n2LFjheJwV7aoqCh69erFggUL2LNnj38fGjLEaSYX2yzeYgXFuCp7uJVWrVq5HSWk+vTpQ2ZmJu+//75/H1izxmkmF9ss3mIFxbgqKSmJxMREypUr53aUkKpfvz5NmjSxw14molhBMa45cuQIX3/9dUQOBumPPn36sHLlygu/PbAxHmEFxbjms88+Iysrq9Ad7sp2QUOxGONBVlCMa1JSUihevDjXXXed21FcUaVKFdq2bcu777577qFY6tRxmsnFNou3FHE7gCm8kpOTadasGTExMW5HcU2fPn24++67Wbx4MS1btjzzjGPHhi5UGLHN4i22h2JcceDAAdauXXv2X6KFgA3FYiKJFRTjisWLF6Oqhbb/JFtsbKx/Q7EMGuQ0k4ttFm+xgmJckZKSQmxsLFdffbXbUVzn11AsmzY5zeRim8VbrKAYV6SkpHDjjTdSrFgxt6O4zoZiMZHCCooJuT179rBu3bpC33+S7byGYjHGg6ygmJD79NNPAayg5FDgoViM8aDzKigiUlJEogMdxhQOycnJlClThiZNmrgdxTPOORRL48ZOM7nYZvEWv65DEZEooAfQC7gaOAUUF5H9wFxgjKpuCVpKE1FSUlJo3rw5RYrYZVA59enTh4cffpj169dTv3793G+OGuVOKI+zzeIt/u6hpACXA08CF6tqnKpWBm4AlgHPi0jvIGU0EWTHjh1s3ry50J8unJ/soVgmTpzodhRjzou/BaWNqo5Q1W9U9dcxIlT1oKp+oKp3AlODE9FEkpSUFMD6T/JTpUoV2rdvz7vvvktmZmbuN3v3dprJxTaLt/hVUFT1dCDmMSYlJYXy5cvTqFEjt6N4Ut++fdm5cyfJycm539ixw2kmF9ss3nLOgiIiN4vIWyLS2Pfarks15y05OZkWLVoQFWUnGOanY8eOlCtXzg57mbDkz//qe4DHgN4i0gqwcyrMefnxxx/5+eef7XDXWcTExNC9e3dmzJjBsWPH3I5jTIH4U1COqephVX0UaItzlpcxBZbdf2Id8mfXr18/UlNTmT59uttRjCkQfwrK3OwnqvoEcMH74iLSTkQ2isgWEXkin/eLi8hU3/vLRaSmb/o1IrLG19aKyO0XmsWETnJyMlWqVOGKK65wO4qnXXvttdSuXTv3Ya/rrnOaycU2i7eIqvo/s0gxVU2/oBU6F0RuAm4GdgBfAz1VdX2OeR4AGqnqfSLSA7hdVbuLSCyQrqoZIlIVWAtUU9WMM60vISFBV6xYcSGRTQCoKjVq1ODGG2+0q8H98Oyzz/L000/z448/UrNmTbfjmEJIRFaqakJBPlPQntFlItKwgJ/J6xpgi6r+4CtO7wOd88zTGZjgez4daC0ioqqpOYpHDOB/NTSu2rx5M7/88osd7vJTb9+5sO+++67LSYzxX0ELyr3AZBF5NO8bIjLfz2VUB7bneL3DNy3feXwF5AhQwbeeRBFZB3wL3Jff3omIDBKRFSKyYt++fX7GMsGUfRqsdcj7p2bNmrRo0YKJEyeiqnDnnU4zudhm8ZYCFRRV/RpIBJqIyCci8oCIjBGRbwu6rPOlqstVtQHOyQFPishv7h+rqmNVNUFVEypVqhSKWOYcUlJSqF69OvHx8W5HCRt9+/Zl8+bNLFu2DA4ccJrJxTaLtxSoCIjIs8B3QCNgD/A0UBJop6pt/VzMTiAux+savmn5ziMiRYCyQK6vjapuAI4DVxbkZzChp6qkpKTQqlUrRMTtOGGja9eulChRggkTJpx7ZmM8oKB7Ff2Bq1X1SlW9G6ewlAZeEpEyfi7ja6C2iNQSkWI4g07mvVXdLKCf73lXIFlV1feZIgAicilQD/ipgD+DCbF169axb98+O9xVQKVLl+aOO+5g6tSpZGVlnfsDxrjMr4Ii//uzso6qHsyerqr7VLUzzuCRX/uzLF+fx2BgAbABmKaq60RkuIh08s02DqggIluAPwPZpxbfAKwVkTXAh8ADqrrfn/Ua99j1J+evb9++HD58mP12XMeEAX/HD08RkQ+AmcC27Im+PYwbfO1lf1eqqvOAeXmmDc3xPA24K5/PTQLsPqlhJjk5mVq1anHppZe6HSXstG7dmmrVqpGUkUHv1q3djuM5tkm8xd+C0g5nCJYpInIZcAjntN1oIAl4RVVXByeiCWeZmZksXryY22+3a1DPR3R0NL1796b/yy9z86BBVHE7kMc8/bTbCUxO/o42nKaqb6jq9cAlQGugiapeqqp/sGJizmTt2rUcOnTI+k8uQN++fcnMzOS9995zO4oxZ+VvH0o/EdkvIgeBt4Hjqno4uNFMJLD+kwvXoEEDPi9Thqv//ncKMrJFYdC+vdOMN/h7ltfTOEOl1MPpQ3kuaIlMRElJSaFOnTpUq1bN7ShhrWaVKmQcO8bKL790O4qnnDzpNOMN/haUo6q6WlX3qurTOMOnGHNWGRkZLFmyxPZOAqBK5cpEifD2Qw+5HcWYM/K3oFT1DWfSXEQqAUWDGcpEhpUrV3Ls2DHrPwmAIkWKUKlMGaZ88w0nDh1yO44x+fK3oPwdaAiMADYCV4rIPBH5h4j0DFo6E9ay+09atGjhbpAIUbVCBY6ePs30115zO4ox+fL3LK+xqvonVb1JVcsDlwGvAYeBDsEMaMJXcnIyV155JZUrV3Y7Svi77TbKtmpFfGws48aPdzuNZ9x2m9OMN/h7HUouqroDZ5TgjwMbx0SK9PR0vvjiCwYOHOh2lMjw6KNIUhIDV63iyVWr2LhsGXWvvdbtVK579Dfjnhs3hWSEYFP4fPXVV6Smplr/SYD1i48nWoT/vPSS21GM+Q0rKCYokpOTERFuuukmt6NEhhYt4LHHqBoby62XX874efM4nZbmdirXtWjhNOMNVlBMUKSkpNC4cWPKly/vdpSI8/uEBPaePMnct95yO4oxuVhBMQF38uRJli5daoe7gqR97dpUjY3l7TFj3I5iTC5WUEzALV26lFOnTllBCZIiUVH0b9yYj9evZ+f69W7HMeZXVlBMwKWkpBAdHU3z5s3djhKx7mnalCxV3nnhBbejGPMrKygm4FJSUmjatCllyvh7E09zTt26QY4CHV++PK0uuYS3Zswg8/RpF4O5q1s3pxlvsIJiAur48eMsX77cDncF2gMPQMeOuSclJrLt2DHmvf22S6Hc98ADTjPeYAXFBNQXX3xBRkaGDQgZaKmpkOc04U5161I1NpY3Xn/dpVDuS011mvEGKygmoFJSUihatCjXX3+921EiS4cOv7k9YdHoaAY1bcqC9evZumqVS8Hc1aGD04w3WEExAZWcnMw111xDyZIl3Y5SKPwhIYEoEcY8Z7coMu6zgmIC5siRI6xcudIOd4VQ9TJl6BQfz3/mziXt+HG345hCzgqKCZglS5aQlZVlHfIh9sC113IgLY3p//qX21FMIWcFxQTMokWLiImJ4brrrnM7SqHSqlYtapctyxt25bxxmRUUEzALFy6kefPmxMTEuB0l8vTvDzffnO9bUSLcn5jI0m3bWJOUFNpcLuvf32nGG6ygmIDYuXMnGzZs4OYz/NIzF6h/f2jb9sxvN25MiehoXvvHP0KXyQOsoHiLFRQTEIsWLQKwghIs+/fDkSNnfPuiEiXo16gRkz/7jL0//hjCYO7av99pxhusoJiAWLhwIZUqVaJhw4ZuR4lMXbvCs8+edZaHmjXjVGYmbw4fHqJQ7uva1WnGG1wpKCLSTkQ2isgWEXkin/eLi8hU3/vLRaSmb/rNIrJSRL71Pdr5qR6gqixatIg2bdoQFWV/o7ilXsWKtK9VizemTeOUXT5uXBDy//0iEg2MBtoD9YGeIlI/z2wDgUOqGg+8Ajzvm74f6KiqDYF+wKTQpDZn891337Fnzx473OUBQ66/nj2pqUx95RW3o5hCyI0/J68BtqjqD6qaDrwPdM4zT2dggu/5dKC1iIiqrlbVX3zT1wElRKR4SFKbM1q4cCEAbdq0cTmJufmyy6h/0UWMeuMNVNXtOKaQcaOgVAe253i9wzct33lUNQM4AlTIM8+dwCpVPRWknMZPCxcupG7dusTFxbkdpdATEYZcfz2rf/mFzz74wO04ppAJywPeItIA5zDYvWd4f5CIrBCRFfv27QttuELm1KlTLFmyxA53Bdv998Ott/o1a+9GjahQvDijCsEpxPff7zTjDW4UlJ1Azj9la/im5TuPiBQBygIHfK9rAB8CfVV1a34rUNWxqpqgqgmVKlUKcHyT09KlS0lNTbWCEmzdu0OLFn7NWqJoUe5LSOCjVavYunJlcHO5rHt3pxlvcKOgfA3UFpFaIlIM6AHMyjPPLJxOd4CuQLKqqoiUA+YCT6jqFyFLbM5o4cKFREdHc9NNN7kdJbJt3w579/o9+x8TEykaFcVLf/1rEEO5b/t2pxlvCHlB8fWJDAYWABuAaaq6TkSGi0gn32zjgAoisgX4M5B9avFgIB4YKiJrfK1yiH8Ek8OiRYtITEykbNmybkeJbH36wIsv+j171dKl6d+oEe8sWsTuH34IYjB39enjNOMNrvShqOo8Va2jqper6v/5pg1V1Vm+52mqepeqxqvqNar6g2/6s6paUlUb52j+/9lmAurQoUOsWLHCDnd51GM33sjprCxejfC9FOMdYdkpb7xh4cKFZGVl0fYsY0wZ98SXL89ddevyxowZHLGTU0wIWEEx523evHmUL1+exMREt6OYM/hL8+YcTU/n30OHuh3FFAJWUMx5ycrKYv78+bRt25bo6Gi345gz+F3Vqtxy6aWMmjiRk3ZHRxNkVlDMeVm9ejV79uyhQ4cObkcpHB55BO6887w++kSLFuxJTeWdcwwuGY4eecRpxhusoJjz8vHHHwNwyy23uJykkOjYEa699rw+etOll9KsalX+8cYbnDp5MsDB3NWxo9OMN1hBMefl448/JiEhgcqV7aztkNi48bwvuBARhrVuzY5jxxg3YkSAg7lr40anGW+wgmIK7ODBgyxbtswOd4XSvffCv/513h9vc9ll3FCtGs+9/jppETS0/b33Os14gxUUU2BJSUlkZWXRvn17t6MYP4kIz7Ruzc5jx3g7wvZSjHdYQTEFNm/ePCpUqMDVV1/tdhRTAC1r1aJ5tWr84/XXOXnihNtxTASygmIKJCsriwULFnDLLbfY6cJhRkR4pk0bfjl+nLF//7vbcUwEsoJiCmT58uXs3buXW/0cSt14S4tatWhRowb/ePNNjh8+7HYcE2GsoJgCmTlzJkWKFLEO+VD729+gZ8+ALOq5tm3Zc+IEIx97LCDLc9Pf/uY04w1WUEyBfPTRR7Ro0YJy5cq5HaVwadMGmjQJyKKui4vjztq1eXHCBPaE+djvbdo4zXiDFRTjt40bN7Jx40Y6d+7sdpTCZ80a2Jrv/eTOy3O33MLJjAyGDx4csGW6Yc0apxlvsIJi/DZz5kwAKyhuGDIE3nwzYIurU6EC9zZuzNg5c9i0enXAlhtqQ4Y4zXiDFRTjt5kzZ9KkSRPi4uLOPbPxvKGtWhETFcVTDzzgdhQTIaygGL/s2bOHpUuX2t5JBKlSqhSPXncdHyxbxhdz5rgdx0QAKyjGL7Nnz0ZVraBEmEdvuIEaJUsy+P77yczIcDuOCXNWUIxfZs6cSc2aNWnUqJHbUUwAlSxWjJdvuYU1O3YwNgKHtzehZQXFnNPRo/TGH5oAABNbSURBVEdZuHAhXbp0QUTcjlM4Pfcc9O8flEXfdeWVtKxRg7+98AIH9uwJyjqC5bnnnGa8wQqKOadZs2Zx6tQpunXr5naUwqtZM2jQICiLFhFe69iRI2lp/PUPfwjKOoKlWTOnGW+wgmLOadq0acTFxdm949305Zewbl3QFt+gcmX+1KQJY2fPZuWnnwZtPYH25ZdOM95gBcWc1eHDh5k/fz533XUXUVH2dXHNU0/B+PFBXcWwNm2oXKIE9/bvT0aYdNA/9ZTTjDfYbwhzVjNnzuT06dN0797d7SgmyMrGxPBahw6s/PlnXnniCbfjmDBkBcWc1dSpU6lZs6bd+6SQ6NqgAbdffjlDX32Vzd9953YcE2asoJgzOnjwIAsXLqRbt252dlchISKM7tyZ4iL8oXt3srKy3I5kwogVFHNG06dPJyMjw87uKmSqli7NyLZtWbx+PW/ZObmmAFwpKCLSTkQ2isgWEfnNwVoRKS4iU33vLxeRmr7pFUQkRUSOi8jroc5d2EyYMIH69evTJEDDppsLMGoU3HdfyFY3oGlTWteowaMjRvDDhg0hW29BjRrlNOMNIS8oIhINjAbaA/WBniJSP89sA4FDqhoPvAI875ueBjwNPBqiuIXW5s2b+fLLL+nXr58d7vKCxo3h8stDtjoR4T9duxKtSu8uXTx71lfjxk4z3uDGHso1wBZV/UFV04H3gbwDRHUGJvieTwdai4io6glV/RynsJggmjRpElFRUfTq1cvtKAZg0SJYtSqkq7ykbFn+feutLN20if97+OGQrttfixY5zXiDGwWlOpDzNnE7fNPynUdVM4AjQIWQpDNkZWUxceJE2rRpQ/Xqef9pjCuefRamTAn5antedRW969VjxOjRLPvkk5Cv/1yefdZpxhsislNeRAaJyAoRWbFv3z6344SdJUuW8PPPP9OvXz+3oxgPeL1zZ2qULEnPbt04uH+/23GMh7lRUHYCOe/QVMM3Ld95RKQIUBY44O8KVHWsqiaoakKlSpUuMG7hM378eEqXLk2XLl3cjmI8oGxMDFO7dWPn4cP0ve02O5XYnJEbBeVroLaI1BKRYkAPYFaeeWYB2X8edwWSVVVDmLHQOnjwIFOnTqVnz57Exsa6Hcd4RGJcHK+0bs3c5cv5x2OPuR3HeFTIC4qvT2QwsADYAExT1XUiMlxEOvlmGwdUEJEtwJ+BX08tFpGfgJFAfxHZkc8ZYuYCTJgwgbS0NB6w28KaPB647jrurluXp0eOZOFHH7kdx3iQRPof/gkJCbpixQq3Y4SFrKws6tWrR8WKFfnShnD1lo0b4fPPnXbppa7FOJGeTuKbb/JLWhrLvv6aOkEaUt9fGzc6j3XruhojIonISlVNKMhnIrJT3pyf5ORkNm/ebHsnXlS3LsTFnXu+ICtZrBizevcmOiuLW1u35oDLnfR161ox8RIrKOZX//73v6lQoQJdu3Z1O4rJa/ZsWLbM7RQAXFa+PB/16MG2ffvo2qYN6enprmWZPdtpxhusoBgAtm3bxsyZM7nnnnuIiYlxO47J6+WX4YMP3E7xq+tr1mRchw58unYtA++4w7Uzv15+2WnGG6ygGABG+QZEGjx4sMtJTLjo3bQpI66/nnfnzuXPAwcS6f2x5tyKuB3AuO/w4cO89dZb9OjRg0suucTtOCaM/LV1aw6cOMGo8eOpULEiT7/4otuRjIusoBjefPNNjh8/zmN2fYEpIBHh5U6dOJSWxtCXXqJkqVL8+e9/dzuWcYkd8irkTp48yauvvsrNN9/MVVdd5XYcE4aiRHj7rrvoetllPDJsGC8MHep2JOMS20Mp5N588012797N+++/73YUczaTJsHixeDBARoBikRFMaVXL4pMmcJfRowgPT2dv/3zn0Ff76RJQV+FKQDbQynETpw4wT//+U9at27NTTfd5HYcczZxcVC5stspzqpIVBSTevakT506PP388zxy771BP/srLs4Tl+cYHysohdjo0aPZu3cvzzzzjNtRzLlMnQqffup2inMqEhXFO927M7hhQ0aOHUuPjh1JSwve7YumTnWa8QYbeqWQOnToEPHx8Vx99dXMnz/f7TjmXFq0gEOHoEkTV4de8ZeqMjIlhUc/+4wbfvc7Plq4kAoVAn9LoxYtnMcwqLVhx4ZeMX575plnOHz4MC+88ILbUUwEEhEeadWK92+7ja/WruWahg1Zs3q127FMkFlBKYS+//57Ro8eze9//3saNWrkdhwTwbo3bcqnvXqRdvgw1yUmMnHcOLcjmSCyglLIqCoPP/wwsbGxjBgxwu04phC47rLLWHXffVxboQL9fv977uvXj9TUVLdjmSCwglLIvPfee8yfP5/hw4dT2eNnDZnIUaVMGRb+4Q883rgxYyZO5HdXXMFyjwx2aQLHOuULkX379nHFFVdQu3ZtPv/8c6Kjo92OZPy1fz+kpMC8eWHRKX82yRs2MGDuXHakpvLUo4/y9LPPUqxYsfNaVvbo+RUrBjCgAaxT3pyFqvLAAw9w7Ngxxo0bZ8Uk3FSsCGXLup0iIFpdcQXfPPAAfePjefbFF2lUty6LFi48r2VVrGjFxEusoBQSY8aMYfr06YwYMYL69e2uyWFn/HhISnI7RcCUjY3lnbvv5uPbbyfjwAFubtuW7l26sGPHjgItZ/x4pxlvsENehcDatWtJTEykZcuWzJ07l6go+zsi7ITZdSgFkXb6NC8uWMBza9ZAdDR/+uMf+ctf/+rXdSt2HUrw2CEv8xu7d++mc+fOVKhQgYkTJ1oxMZ4TU7QoT992GxsGDaJ7XBwvvfIKl116KSOGDePQoUNuxzMFYL9dIlhqaiqdOnVi3759zJo1i0qVKrkdyZgzqlm5MuN79+bbAQNoWb48Q595hrjq1Xn4wQfZtm2b2/GMH6ygRKiTJ09yxx13sGLFCqZMmULTpk3djmSMXxpccgkf3XMPa/v25Y6qVXl99Gguq1WLHnfeSXJysmu3GzbnZgUlAmXvmSQlJTFu3Dg6derkdiRjCqxRrVpM7NOHHwYN4qErriBpzhxat25NfM2a/N+zz9peiwdZp3yE+eWXX+jSpQsrVqzgnXfeoV+/fm5HMoGQmgqLFsGHH0Zcp7y/0tLS+PCrr3j7++9J3rULgITfNeOOu7px991duLSQbpdgOZ9OeSsoEeTzzz+ne/fuHDlyhMmTJ9O5c2e3I5lASkqCKVMKbUHJ6YcdO5i2ejX//flnVh04AED92vXp0LkD7dq14/rrrycmJsbllOHNCko+CkNBOXnyJMOGDePFF1+kVq1afPjhhzboY6R54w3YsAGOH7eCksMbXyUwe31pfjr8EYdPLWbvqW/J0gxiisXQ7Npm3NjyRpo1a0ZiYiJlI+TC0FA5n4JitwAOY5mZmUyZMoWnnnqK7du3M2jQIF5++WVKlSrldjQTaNOm/e86FPOraeuv5JdjpWhbpwQVYztz+NhJ1u9ey44TX7F62TckL0kBFBGhXnw9Epsl0uiqRjRq1IiGDRvaeHYBZgUlDB04cID33nuPUaNG8cMPP9CkSRMmTZpkt/E1hV650iVoVvpa4Fo0Szlw9ARb929iz8l17N25jmmTZzF+wvhf569UvhL1G9Sndt3aXH755bma7dEUnCsFRUTaAa8C0cDbqvrPPO8XByYCTYEDQHdV/cn33pPAQCATeFBVF4QwuitUlU2bNpGSksKcOXNYsGABGRkZJCYm8vzzz3PHHXfYBYvG5CFRQsVypahYrgng7NmdTlf2HjnCjiM/cfDUVg6d2sL6r7fz1dK1nMw4nOvzpUuWpurFVakeV53qNapTrVo1qlWrRtWqValYsSLly5f/tZUsWRIRceGn9JaQFxQRiQZGAzcDO4CvRWSWqq7PMdtA4JCqxotID+B5oLuI1Ad6AA2AasAiEamjqpmh/SkCKzMzkxMnTnD06FF27drFzp072blzJ5s2beK7777j22+/Zd++fQBccsklPPzww/Ts2ZPGjRvbl9iYAihaTKheqRzVKzUGGgOgWUp6Ohw6kcq+47s5ePIXjmfsJDVjHyd2H+C7nXtZkbmeExkHyNKM/JdbpCjlypbjoosu4qKLLqJ0mdKUKl2KUqVKUbJkSUqVyv08NjaW4sWLU7x4cYoVK/ab53mnFStWjOjo6N80r/3/d2MP5Rpgi6r+ACAi7wOdgZwFpTMwzPd8OvC6OFuuM/C+qp4CfhSRLb7lLT3TyjZu3Ejz5s3JPvlAVUPy/FzzZWZmcvz4cY4fP37Gmw3FxsbSoEEDOnbsyDXXXEOrVq2Ij4/33JfImHAmUULxGLg4piQXV7gcuDzX+9kF51Q6HE07ypGTBzlx+jhpmUc5nXWUU1nHOJV5lFOnj3Jq11G27TxBhu4iU0+SqSc5rWmczjzJ6ay0wGcXITrKKS5RUVH/KzbZ06KjiY6O+vW1iBAVFYWIOC1Kck2LkqhfXwMF/kUT8rO8RKQr0E5Vf+973QdIVNXBOeb5zjfPDt/rrUAiTpFZpqrv+qaPAz5W1el51jEIGOR7eSXwXVB/qOCqCOx3O8QFsPwBUhpKVga/x885AjFlIfC/xULE3/zHKF88lbJFQ5GpIE5zNLooZQpw9ERzPNE8bzivxTeX85v+jL+8Nc/y/F5xzpcn2FsiU9MLtNMRkZ3yqjoWGAsgIisKeuqbl1h+d4VzfhFZsS9Ms0Nk5E/T/WGdv6CfcaMndycQl+N1Dd+0fOcRkSJAWZzOeX8+a4wxxgVuFJSvgdoiUktEiuF0ss/KM88sIHvMkK5AsjrH5mYBPUSkuIjUAmoDX4UotzHGmLMI+SEvVc0QkcHAApzThv+jqutEZDiwQlVnAeOASb5O94M4RQfffNNwOvAzgD/6cYbX2GD9LCFi+d0VzvnDOTtYfrcVOH/ED71ijDEmNOxqOGOMMQFhBcUYY0xARHRBEZF2IrJRRLaIyBNu5zkXEfmPiOz1XYeTPa28iCwUkc2+x4vczHgmIhInIikisl5E1onIQ77p4ZI/RkS+EpG1vvzP+KbXEpHlvu/QVN+JJJ4lItEislpE5vheh01+EflJRL4VkTXZp6yGy/cHQETKich0EfleRDaIyHXhkF9E6vq2eXY7KiJDzid7xBaUHEO8tAfqAz19Q7d42XigXZ5pTwCfqGpt4BPfay/KAB5R1frAtcAffds7XPKfAlqp6lU4Y3K0E5FrcYb9eUVV44FDOMMCedlDwIYcr8Mtf0tVbZzj2p9w+f6AMz7hfFWtB1yF8+/g+fyqutG3zRvjjJ+YCnzI+WTPOVRIJDXgOmBBjtdPAk+6ncuP3DWB73K83ghU9T2vCmx0O6OfP8dMnPHawi4/EAuswhmdYT9QJL/vlNcaznVZnwCtgDk4F1SHU/6fgIp5poXF9wfnWrkf8Z3oFG75c+RtC3xxvtkjdg8FqA5sz/F6h29auKmiqrt8z3cDVdwM4w8RqQn8DlhOGOX3HS5aA+wFFgJbgcOqv44I6PXv0CjgcSDL97oC4ZVfgSQRWekbPgnC5/tTC9gHvOM75Pi2iJQkfPJn6wFM8T0vcPZILigRR50/FTx9nreIlAI+AIao6tGc73k9v6pmqrPbXwNn0NF6Lkfym4jcBuxV1ZVuZ7kAN6hqE5zD1H8UkeY53/T496cIzhj5/1bV3wEnyHOIyOP58fWvdQL+m/c9f7NHckGJlGFa9ohIVQDf416X85yRiBTFKSaTVXWGb3LY5M+mqoeBFJxDROV8w/+At79D1wOdROQn4H2cw16vEj75UdWdvse9OMfwryF8vj87gB2qutz3ejpOgQmX/OAU8lWqusf3usDZI7mg+DPESzjIOQxNP5y+Cc8REcEZ4WCDqo7M8Va45K8kIuV8z0vg9P9swCksXX2zeTa/qj6pqjVUtSbOdz1ZVXsRJvlFpKSIlM5+jnMs/zvC5PujqruB7SJS1zepNc6IHmGR36cn/zvcBeeT3e1OoCB3MHUANuEcC/+r23n8yDsF2AWcxvmLZyDOcfBPgM3AIqC82znPkP0GnF3ib4A1vtYhjPI3Alb78n8HDPVNvwxnvLgtOIcCirud1Y+fpQUwJ5zy+3Ku9bV12f9fw+X748vaGFjh+w59BFwULvmBkjgD8JbNMa3A2W3oFWOMMQERyYe8jDHGhJAVFGOMMQFhBcUYY0xAWEExxhgTEFZQjDHGBIQVFGOMMQER8lsAG2P+R0S6ALcCZYBxqprkciRjzptdh2KMB/juNfGSqnp9eHljzsgOeRnjDX/DuX+PMWHLCooxIeAbq2qniDzne3217+54JUTkeeBjVV3lckxjLogd8jImRESkAs5YTw1w7hXTB2iOM/De18AaVX3TvYTGXBgrKMaEkIiswxk8cLWqvuB2HmMCyQ55GRNa3wAXAy+5HcSYQLOCYkyIiEgloCXwgapmnWt+Y8KNHfIyJkREZBZwHDihqn9wO48xgWZ7KMaEgIjcC5wE/oJza2FjIo7toRgTZCJSG5gNXKeqh0RkIaCq2tblaMYElBUUY4wxAWGHvIwxxgSEFRRjjDEBYQXFGGNMQFhBMcYYExBWUIwxxgSEFRRjjDEBYQXFGGNMQFhBMcYYExBWUIwxxgTE/wMzcO9bQH2TnAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Plot chi-squared probability distribution\n", "# Indicate areas under the curve at values of chi-squared greater than a couple of test values\n", "nu = 30 # number of degrees of freedom\n", "chisq1 = 37.5 # 1st value of chi-squared to test\n", "chisq2 = 52.1 # 2nd value of chi-squared to test\n", "# Range of values for chi-squared\n", "xmin,xmax,nx = 0.,70.,701\n", "x = np.linspace(xmin,xmax,nx)\n", "dx = (xmax-xmin)/nx\n", "ymin,ymax = 0.,0.06\n", "# Determine chi-squared probability distribution\n", "y = stats.chi2.pdf(x,nu)\n", "plt.plot(x,y,color='black')\n", "# Find 1st chi-squared test value in array, and shade area under the curve at greater values\n", "plt.plot([chisq1,chisq1],[ymin,ymax],color='red',linestyle='dashed')\n", "i = int((chisq1-xmin)/dx)\n", "plt.fill_between(x[i:],np.zeros_like(x[i:]),y[i:],color='red',alpha=0.5)\n", "# Find 2nd chi-squared test value in array, and shade area under the curve at greater values\n", "plt.plot([chisq2,chisq2],[ymin,ymax],color='blue',linestyle='dashed')\n", "i = int((chisq2-xmin)/dx)\n", "plt.fill_between(x[i:],np.zeros_like(x[i:]),y[i:],color='blue',alpha=0.5)\n", "# Complete plot\n", "plt.xlabel(r'$\\chi^2$')\n", "plt.ylabel(r'$P(\\chi^2)$')\n", "plt.xlim(xmin,xmax)\n", "plt.ylim(ymin,ymax)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The resulting $p$-values in these cases are:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "p-value for chi2 = 37.5 and nu = 30 is p = 0.16301345461145947\n", "p-value for chi2 = 52.1 and nu = 30 is p = 0.007429708671565914\n" ] } ], "source": [ "nu = 30 # number of degrees of freedom\n", "chisq1 = 37.5 # 1st value of chi-squared to test\n", "chisq2 = 52.1 # 2nd value of chi-squared to test\n", "p1 = stats.chi2.sf(chisq1,nu) # probability of chi-squared exceeding the 1st test value\n", "p2 = stats.chi2.sf(chisq2,nu) # probability of chi-squared exceeding the 2nd test value\n", "print('p-value for chi2 =',chisq1,'and nu =',nu,'is p =',p1)\n", "print('p-value for chi2 =',chisq2,'and nu =',nu,'is p =',p2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The null hypothesis (that the dataset is sampled from the model) has not been rejected for the first value of $\\chi^2$, but has been rejected for the second value of $\\chi^2$. In other words, we say that the first dataset is consistent with being sampled from the model, but the second dataset is not.\n", "\n", "### Caveats when interpreting $\\chi^2$ values\n", "\n", "- We are assuming that the errors in the data are **Gaussian** and **reliable**\n", "- If the errors have been **under-estimated**, then an **improbably high** value of $\\chi^2$ can be obtained\n", "- If the errors have been **over-estimated**, then an **improbably low** value of $\\chi^2$ can be obtained\n", "- Since errors can often be non-Gaussian or approximate, a model is typically only rejected for very low values of $p$ such as $0.001$ or below\n", "- A possible issue in using the $\\chi^2$ statistic is **binning of data**. If the numbers of points in each bin are too small, then the probabilities can become non-Gaussian. As a rule-of-thumb, 80% of bins must have $N > 5$.\n", "\n", "### Reduced $\\chi^2$\n", "\n", "As a way of summarising the model fit, we can quote the **reduced $\\chi^2$ statistic**, $\\chi_r^2 = \\chi^2/\\nu$.\n", "\n", "For a good fit, $\\chi_r^2 \\sim 1$ (because $\\overline{\\chi^2} = \\nu$). But note that the probability of the data being consistent with the model also depends on $\\nu$. **Do not just quote the reduced $\\chi^2$ value**.\n", "\n", "### Modification for correlated data\n", "\n", "As previously mentioned, the above equation for $\\chi^2$ is only valid when the data points are uncorrelated. If you are dealing with correlated data, the equation must be modified:\n", "\n", "$\\chi^2 = \\sum_i \\sum_j (d_i - m_i) \\left( C^{-1} \\right)_{ij} (d_j - m_j) = (d - m)^T C^{-1} (d - m)$\n", "\n", "where the last version is written in vector/matrix notation. Here, $C$ is the covariance matrix of the data:\n", "\n", "$C_{ij} = \\langle x_i x_j \\rangle - \\langle x_i \\rangle \\langle x_j \\rangle$\n", "\n", "Note that the number of degrees of freedom is unchanged for correlated data.\n", "\n", "## Parameter estimation using $\\chi^2$\n", "\n", "A model typically contains free parameters. How do we determine the most likely values of these parameters and their error ranges?\n", "\n", "- Suppose we are fitting a model with 2 free parameters $(a,b)$\n", "- The most likely (\"best-fitting\") values of $(a,b)$ are found by minimizing the $\\chi^2$ statistic\n", "- The joint error distribution of $(a,b)$ can be found by calculating the values of $\\chi^2$ over a grid of $(a,b)$ and enclosing a particular region $\\chi^2 < \\chi_{\\rm min}^2 + \\Delta \\chi^2$ (see below).\n", "\n", "This process is described in the following worked example. Consider the following example dataset containing $N=10$ points:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Text(0, 0.5, '$y$')" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEGCAYAAABlxeIAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAATyElEQVR4nO3db4xdd53f8fcnk7DAeEWqyVSbJnEGlWgrL7vYMMqEIlUpFDnDWskDWDWWym4QK3crINAidSkPyN20D7pVBVsWBIoIu9BGhiqkq3QU7zQq8bBUxZtx7ITEBtWi9iZpKgazJIx3NzTjbx/MsTOezDh2mHvOndz3S7ryvef8fO9X15r5+Pf9nT+pKiRJuqTrAiRJg8FAkCQBBoIkqWEgSJIAA0GS1Li06wJeqSuuuKImJia6LkOSNpWDBw/+qKrG19q3aQNhYmKC+fn5rsuQpE0lyYn19tkykiQBBoIkqWEgSJIAA0GS1DAQJElAC4GQ5LVJ/jzJo0meSPJ7a4y5LclCksPN47f7XZck6VxtHHb6PPDOqlpMchnw7ST7quo7q8Z9vao+3EI9kqQ19D0Qavn62ovNy8uah9fclqQB08oaQpKRJIeBHwIPVtWBNYa9N8ljSe5Ncs0677MnyXyS+YWFhb7WLEmDptfrkeTso9frbej7p80b5CS5HPgvwEeq6vEV28eAxap6Psk/Bf5xVb3zfO81OTlZnqksadjceOONAOzfv/8V/f0kB6tqcq19rR5lVFU/AR4Cblq1/WRVPd+8/BLwtjbrkiS1c5TReDMzIMnrgHcD31s15soVL28Gjva7LknSudqYIVwJPJTkMeBhltcQZpLcmeTmZsztzSGpjwK3A7e1UJckXZB+9+4HRatrCBvJNQRJbfp5e/eDUsfArCFIkgaXgSBJAgwESQNsWHr3g2LT3jFN0qtfr9c72yvvunc/DJwhSJIAA0GS1DAQJEmAgSBJahgIkiTAQJAkNQwESRJgIEiSGgaCJAkwECRJDQNBkgQYCJKkhoEgSQIMBElSo++BkOS1Sf48yaPNfZN/b40xv5Dk60mOJTmQZKLfdUmSztXGDOF54J1V9RZgO3BTkhtWjfkg8JdV9SbgM8Dvt1CXJG0qS0tLnDx5khMnTjAzM8PS0tKGvn/fA6GWLTYvL2setWrYLcBXmuf3Au9Kkn7XJg0i7xKmtSwtLbFz506OHDnC8ePH2b17Nzt37tzQUGjljmlJRoCDwJuAz1fVgVVDrgKeBKiqF5I8C4wBP1r1PnuAPQBbt27td9lSJ7xLmNayb98+Dhw4wOnTpwFYXFzkwIED7Nu3j127dm3IZ7SyqFxVS1W1HbgauD7Jm1/h+9xVVZNVNTk+Pr6xRUrSADt06BCnTp06Z9upU6c4fPjwhn1Gq0cZVdVPgIeAm1btehq4BiDJpcAbgJNt1iZJ6+l37/5C7Nixg9HR0XO2jY6Osn379g37jDaOMhpPcnnz/HXAu4HvrRp2P/BbzfP3Ad+sqtXrDJLUujZ69xdienqaqakpLrlk+df2li1bmJqaYnp6esM+o40ZwpXAQ0keAx4GHqyqmSR3Jrm5GXM3MJbkGPAvgE+0UJckvazz9e7bNDIywuzsLNu2bWNiYoK9e/cyOzvLyMjIhn1G3xeVq+oxYMca2z+14vnfAL/R71okbS5nWjWLi4vMzMwwPT29ob8AL8T5evcbtZh7oUZGRhgbG2NsbKwvn+2ZypIG0qC0atro3Q8KA0HSQBqUVk0bvftBYSBIGkhtHGZ5Idro3Q+KVk5Mk6SLdaZVs7i4eHZbV62afvfuB4UzBEkDaZhaNYPCQJA0kIapVTMobBlJGljD0qoZFM4QJEmAgSBJahgIkiTAQJAkNQwESRJgIEiSGgaCJAkwECRJDQNBkgQYCJKkhoEgSQJaCIQk1yR5KMmRJE8k+egaY25M8mySw83jU2u9lzQMztw28sSJE8zMzLR+hzANrzYubvcC8PGqeiTJLwIHkzxYVUdWjfuzqvLqVRpqK28befr0aXbv3s3U1JRX+VQr+j5DqKpnquqR5vlPgaPAVf3+XGkzGpTbRmo4tbqGkGQC2AEcWGP325M8mmRfkl9psy5pUAzKbSM1nFoLhCRbgG8AH6uq51btfgS4tqreAvwh8CfrvMeeJPNJ5hcWFvpbsNSBM7eNXKmr20Zq8PR6Pebm5pibmyMJvV5vQ9+/lUBIchnLYXBPVd23en9VPVdVi83zB4DLklyxxri7qmqyqibHx8f7XrfUNm8bqfPp9XpU1dnHpguEJAHuBo5W1afXGfNLzTiSXN/UdbLftUmDxttGqkttHGX0DuD9wHeTnGmEfhLYClBVXwTeB/yzJC8Afw3cWlXVQm3SwPG2kepK3wOhqr4N5GXGfA74XL9rkaRX4kzvHiAJd9xxx4a3awZBGzMESdrUer3eqzIAVvPSFZIkwECQJDUMBEkDq9/H3etcriFIGljD0rsfFM4QJEmAgSBJahgIkiTAQJAkNQwESRJgIEiSGgaCJAkwECRJDQNBkgQYCJKkhoEgSQIMBElSw0CQJAEGgiSp0fdASHJNkoeSHEnyRJKPrjEmST6b5FiSx5K8td91SZLO1cb9EF4APl5VjyT5ReBgkger6siKMdPAdc1jCvhC86ckqSV9nyFU1TNV9Ujz/KfAUeCqVcNuAb5ay74DXJ7kyn7XJg0i7xKmrqSq2vuwZAL4FvDmqnpuxfYZ4N9W1beb1/8d+N2qml/19/cAewC2bt36thMnTrRUuSS9OiQ5WFWTa+1rbVE5yRbgG8DHVobBxaiqu6pqsqomx8fHN7ZASRpyrQRCkstYDoN7quq+NYY8DVyz4vXVzTZJUkvaOMoowN3A0ar69DrD7gd+szna6Abg2ap6pt+1SZJe1MZRRu8A3g98N8nhZtsnga0AVfVF4AHgPcAx4K+AD7RQlyRphTaOMvp2VaWqfq2qtjePB6rqi00Y0Bxd9KGq+rtV9aurF5OlNvR6PZKcfXh0j4ZNGzMEaVPo9Xrs378f4Oyf0jDx0hWSJMBA0ACwVSMNBltG6pytGmkwOEOQJAEGgiSpYSBIkgADQdIaXOgfTi4qS3oJF/qHkzMESRJwAYGQ5MEkb2mjGElSdy5khvC7wB8k+SPvYvbqYp9Y0kovu4bQ3P7yHyZ5L/CnSe4D/l1V/XXfq1Nf2SeWtNIFrSE09zT4PvAF4CPA/0ry/n4WJklq14WsIfwPlu9e9hngKuA24Ebg+iR39bM4SVJ7LuSw0z3AkaqqVds/kuRoH2qSJHXgQtYQnjjP7l/fwFokSR36uc5DqKofbFQhkqRu9f3EtCRfTvLDJI+vs//GJM8mOdw8PtXvmiRJL9XGpSv+GPgc8NXzjPmzqtrVQi2SpHX0fYZQVd8Cftzvz5Ek/XwG5VpGb0/yaJJ9SX6lXx/imbmStL5BuNrpI8C1VbWY5D3AnwDXrTUwyR6WD4Nl69atF/1BnpkrSevrfIZQVc9V1WLz/AHgsiRXrDP2rqqarKrJ8fHxVuvUq9/S0hInT57kxIkTzMzMsLS01HVJUqs6D4Qkv9RcGoMk17Nc08luq9KwWVpaYufOnRw5coTjx4+ze/dudu7caShoqLRx2Ole4H8Cv5zkqSQfTPI7SX6nGfI+4PEkjwKfBW5d46xoqa/27dvHgQMHOH36NACLi4scOHCAffv2dVyZ1J6+ryFU1e6X2f85lg9LlTpz6NAhTp06dc62U6dOcfjwYXbt8ohoDYfOW0bSIPTud+zYwejo6DnbRkdH2b59e+u1SF0xENSpQendT09PMzU1xSWXLP9IbNmyhampKaanp1utQ+qSgaBODUrvfmRkhNnZWbZt28bExAR79+5ldnaWkZGRVuuQumQgqFPn6923bWRkhLGxMa699lp27dplGGjoGAjqlL17aXAYCOqUvfvBNAgL/WqfgaBO2bsfPIOy0K/2GQjqnL37wTIoC/1qn4EwxGwLaC2DtNCvdhkIQ8q2gNbjQv/wMhCGlG0BrceF/uFlIAwp2wJajwv9w2sQbpCjDpxpCywuLp7dZltAZ5xZ6B8bG/PifkPEGcKQsi0gaTUDYUjZFpC0mi2jIWZbQNJKzhAkSYCBIElqDFUgeGauJK2v74GQ5MtJfpjk8XX2J8lnkxxL8liSt/ajDs/MlaTza2OG8MfATefZPw1c1zz2AF/oRxGemStJ59f3QKiqbwE/Ps+QW4Cv1rLvAJcnuXKj6/DMXEk6v0FYQ7gKeHLF66eabS+RZE+S+STzCwsLF/UhXrBLks5vEALhglXVXVU1WVWT4+PjF/V3PTNXks5vEALhaeCaFa+vbrZtKM/MlaTzG4RAuB/4zeZooxuAZ6vqmX58kHfmkqT19f3SFUn2AjcCVyR5CrgDuAygqr4IPAC8BzgG/BXwgX7XJEl6qb4HQlXtfpn9BXyo33VIL6fX6zE3NwdAEu644w56vV63RUktGoSW0dDp9XokOfvwl85g6PV6VNXZh/8uGjZe7bQDvV6P/fv3A5z9U5K65gxBnTvTqpmbm3PGJHXIGYI61+v1DAFpADhDkCQBBoIkqWEgSJIAA0HSGlzoH04uKkt6CRf6h5MzBEkSYCAMNdsCklayZTTEbAtIWskZgiQJMBAkSQ0DQZIEGAiSpIaBIEkCDARJUqOVQEhyU5LvJzmW5BNr7L8tyUKSw83jt9uoS5L0or6fh5BkBPg88G7gKeDhJPdX1ZFVQ79eVR/udz2SpLW1MUO4HjhWVT+oqp8BXwNuaeFzJUkXoY1AuAp4csXrp5ptq703yWNJ7k1yzVpvlGRPkvkk8wsLCxddiJdqkKT1Dcqi8n8FJqrq14AHga+sNaiq7qqqyaqaHB8fv+gP6fV6VNXZh4EgSS9qIxCeBlb+j//qZttZVXWyqp5vXn4JeFsLdUmSVmgjEB4GrkvyxiSvAW4F7l85IMmVK17eDBxtoa7OLC0tcfLkSU6cOMHMzAxLS0tdlyRJ/Q+EqnoB+DAwy/Iv+v9cVU8kuTPJzc2w25M8keRR4Hbgtn7X1ZWlpSV27tzJkSNHOH78OLt372bnzp2GgqTOpaq6ruEVmZycrPn5+a7LuGgzMzPs3r2bxcXFs9u2bNnC3r172bVrV4eVSRoGSQ5W1eRa+wZlUXloHDp0iFOnTp2z7dSpUxw+fLijiiRpmYHQsh07djA6OnrOttHRUbZv395RRZK0zEBo2fT0NFNTU1xyyfJXv2XLFqamppienu64MknDzkBo2cjICLOzs2zbto2JiQn27t3L7OwsIyMjXZcmach5T+UOjIyMMDY2xtjYmAvJkgaGMwRJEmAgSJIaBoIkCTAQJEkNA0GSBBgIkqSGgSBJAgwESVLDQJAkAQaCJKlhIEiSAANBktQwECRJgIEgSWq0EghJbkry/STHknxijf2/kOTrzf4DSSbaqEuS9KK+B0KSEeDzwDSwDdidZNuqYR8E/rKq3gR8Bvj9ftclSTpXGzOE64FjVfWDqvoZ8DXgllVjbgG+0jy/F3hXkrRQmySp0UYgXAU8ueL1U822NcdU1QvAs8DY6jdKsifJfJL5hYWFPpUrScNpUy0qV9VdVTVZVZPj4+NdlyNJryptBMLTwDUrXl/dbFtzTJJLgTcAJ1uoTZLUaCMQHgauS/LGJK8BbgXuXzXmfuC3mufvA75ZVdVCbZKkxqX9/oCqeiHJh4FZYAT4clU9keROYL6q7gfuBv5jkmPAj1kODUlSi/oeCABV9QDwwKptn1rx/G+A32ijFknS2jbVovKrRa/XY25ujrm5OZLQ6/W6LkmSyGZt1U9OTtb8/HzXZUjSppLkYFVNrrXPGYIkCTAQJEkNA0GSBBgIkqSGgSBJAgwESVLDQJAkAQaCJKmxaU9MS7IAnOi6jgFwBfCjrosYIH4fL/K7OJffx7Jrq2rN+wds2kDQsiTz6511OIz8Pl7kd3Euv4+XZ8tIkgQYCJKkhoGw+d3VdQEDxu/jRX4X5/L7eBmuIUiSAGcIkqSGgSBJAgyETSvJNUkeSnIkyRNJPtp1TV1LMpLkUJKZrmvpWpLLk9yb5HtJjiZ5e9c1dSnJP29+Th5PsjfJa7uuaRAZCJvXC8DHq2obcAPwoSTbOq6pax8FjnZdxID4D8CfVtXfA97CEH8vSa4Cbgcmq+rNwAhwa7dVDSYDYZOqqmeq6pHm+U9Z/oG/qtuqupPkauDXgS91XUvXkrwB+AfA3QBV9bOq+km3VXXuUuB1SS4FXg/8n47rGUgGwqtAkglgB3Cg20o69QfAvwROd13IAHgjsAD8UdNC+1KS0a6L6kpVPQ38e+AvgGeAZ6vqv3Vb1WAyEDa5JFuAbwAfq6rnuq6nC0l2AT+sqoNd1zIgLgXeCnyhqnYAp4BPdFtSd5L8LeAWloPy7wCjSf5Jt1UNJgNhE0tyGcthcE9V3dd1PR16B3BzkuPA14B3JvlP3ZbUqaeAp6rqzIzxXpYDYlj9I+B/V9VCVf0/4D7g73dc00AyEDapJGG5R3y0qj7ddT1dqqp/VVVXV9UEy4uF36yqof0fYFX9X+DJJL/cbHoXcKTDkrr2F8ANSV7f/Ny8iyFeZD+fS7suQK/YO4D3A99NcrjZ9smqeqDDmjQ4PgLck+Q1wA+AD3RcT2eq6kCSe4FHWD467xBexmJNXrpCkgTYMpIkNQwESRJgIEiSGgaCJAkwECRJDQNBkgQYCJKkhoEgbZDm/hTvbp7/myR/2HVN0sXwTGVp49wB3Jnkb7N89dmbO65HuiieqSxtoCRzwBbgxuY+FdKmYctI2iBJfhW4EviZYaDNyECQNkCSK4F7WL7u/mKSmzouSbpoBoL0c0ryepavsf/xqjoK/GuW1xOkTcU1BEkS4AxBktQwECRJgIEgSWoYCJIkwECQJDUMBEkSYCBIkhr/Hz9kf3Y7EZEXAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Read in and plot 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 = len(x)\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.xlabel(r'$x$')\n", "plt.ylabel(r'$y$')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Could these points be fit by a constant $y = b$, and if so what is the best-fitting value of $b$? We can answer that question by minimizing the $\\chi^2$ statistic:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Best-fitting value of b = 1.7899999897904915\n", "Minimum chi-squared = 31.635200000000005\n", "Number of degrees of freedom = 9\n" ] } ], "source": [ "# Function which evaluates the chi-squared value of a model y=b for the (x,y) dataset\n", "def chisqstraightline1(p,x,y,yerr):\n", " b = p[0]\n", " return np.sum(((y - b)/yerr)**2)\n", "\n", "# Find the best-fitting model y=b to the (x,y) dataset by minimising chi-squared\n", "from scipy.optimize import minimize\n", "p0 = [1.] # initial guess of the value of b (result should not depend on this)\n", "result = minimize(chisqstraightline1,p0,args=(x,y,yerr)) # perform the minimisation\n", "b1 = result['x'][0] # best-fitting value of b\n", "minchisq1 = result['fun'] # chi-squared value corresponding to the best fit\n", "dof1 = n-1 # number of degrees of freedom (data points minus fitted parameters)\n", "print('Best-fitting value of b =',b1)\n", "print('Minimum chi-squared =',minchisq1)\n", "print('Number of degrees of freedom =',dof1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Is this value of $\\chi^2$ likely, given the model? We can answer that by computing the $p$-value which corresponds to the probability of obtaining this value of $\\chi^2$ **or higher**, when a dataset is sampled from a model." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "p-value = 0.00023020030084131723\n" ] } ], "source": [ "# Probability of chi-squared to exceed the value we've found, if the data is sampled from the model\n", "# Low values of this probability indicate that the data is unlikely to be drawn from the model\n", "p1 = stats.chi2.sf(minchisq1,dof1)\n", "print('p-value =',p1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This situation is hence very unlikely, suggesting that the model is not a good fit to the data.\n", "\n", "Instead, could these points be fit by a straight line $y = ax + b$? Minimize $\\chi^2$ as before:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Best-fitting value of a = 0.27042421678017053\n", "Best-fitting value of b = 0.4378789090155025\n", "Minimum chi-squared = 7.502540606060829\n", "Number of degrees of freedom = 8\n", "p-value = 0.48350484463274246\n" ] } ], "source": [ "# 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", "# Find the best fitting model y=a*x+b to the (x,y) dataset by minimising chi-squared\n", "p0 = [0.2,1.] # initial guess of the values of a,b (result should not depend on this)\n", "result = minimize(chisqstraightline2,p0,args=(x,y,yerr)) # perform the minimisation\n", "a2,b2 = result['x'] # best-fitting values of a,b\n", "minchisq2 = result['fun'] # chi-squared value corresponding to the best fit\n", "dof2 = n-2 # number of degrees of freedom (data points minus fitted parameters)\n", "p2 = stats.chi2.sf(minchisq2,dof2) # probability of chi-squared to exceed the value we've found\n", "print('Best-fitting value of a =',a2)\n", "print('Best-fitting value of b =',b2)\n", "print('Minimum chi-squared =',minchisq2)\n", "print('Number of degrees of freedom =',dof2)\n", "print('p-value =',p2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This probability is much higher, meaning that the hypothesis that the data is sampled from this model cannot be rejected. We say that the model is a \"good fit\".\n", "\n", "### Errors in the best-fitting parameters\n", "\n", "We can determine the errors in the parameters by determining the value of $\\chi^2$ across a grid of $(a,b)$, and constructing a region of constant $\\chi^2 = \\chi^2_{\\rm min} + \\Delta \\chi^2$. When constructing a joint confidence region between 2 parameters, the relevant intervals are $\\Delta \\chi^2 = (2.30, 6.17, 11.80)$ for (1,2,3)-$\\sigma$ confidence. \n", "\n", "Let's see how this is done in the above case:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Minimum chi-squared = 7.504212500000001\n", "Best-fitting value of a = 0.2725\n", "Best-fitting value of b = 0.42500000000000016\n" ] }, { "data": { "text/plain": [ "Text(0, 0.5, '$b$')" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEGCAYAAAB2EqL0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3dd3RU5fo24PudSSM9EEJCKkhAupBAJFKlN0FFBCuCFFH6TxAFRAQUFRQRITRBDk16k17EUE2o0iQQepAeSCOTmfv7A5jD8aOMe2YyKc+11qyTmOHZT/Y6a+7svd+iSEIIIYR4Ep2jGxBCCJE/SGAIIYSwiASGEEIIi0hgCCGEsIgEhhBCCIs4OboBe/H392dERISj2xBCiHwlMTHxKsniD/tZgQ2MiIgIJCQkOLoNIYTIV5RSZx71M7klJYQQwiISGEIIISwigSGEEMIiEhhCCCEsIoEhhBDCIhIYQgghLCKBIYQQwiISGEIIISwigSGEEMIiEhhCCCEsIoEhhBDCIhIYQgghLCKBIYQQwiISGEIIISwigSGEEMIiEhhCCCEsIoEhhBDCIhIYQgghLCKBIYQQwiISGEIIISwigSGEEMIiEhhCCCEsIoEhhBDCIhIYQgghLCKBIYQQwiISGEIIISwigSGEEMIiEhhCCCEsIoEhhBDCIhIYQgghLCKBIYQQwiJOjm4AAJRSMwC0AnCZZKWH/Lw+gOUAku/9pyUkR+Reh/nHnTt3kJiYiJ07dyI5ORnnz5/HuXPncOPGDWRnZ8NgMCAnJwdFihSBh4cHPDw8UKxYMQQGBiIwMBAlS5ZEqVKlULp0aZQqVQpeXl6O/pWEEHlEnggMADMB/ADg58e853eSrXKnnfzl6tWrmD17NpYsWYI//vgDd+7cAQD4+voiJCQEoaGhKF++PFxdXeHi4gK9Xo/MzEykpaUhPT0d165dQ3x8PC5duoSsrKz/qR0cHIzy5cujQoUKqFy5MqpXr45KlSrBxcXFEb+qEMKB8kRgkNymlIpwdB/5TWpqKsaNG4dx48YhLS0N1apVw/vvv4/atWsjNjYWJUqU+Ff1SOLGjRtITk7GqVOnkJSUhGPHjuHIkSOYPn060tPTAQAuLi6oXLkyatWqhVq1aiE2Nhbh4eFQStnj1xRC5BGKpKN7AADcC4xVj7kltRjAeQAXAfwfycMPeV83AN0AICwsLOrMmTN27Nixdu7ciVdeeQUXLlxAu3btMGzYMFSuXNluxzOZTDh16hT27t2LxMRE/PHHH9izZ485RIKDg1G3bl3UqVMH9erVQ/ny5SVAhMiHlFKJJKMf+rN8EhjeAEwk05RSLQCMJxn5uHrR0dFMSEiwS6+OtmXLFjRt2hShoaGYO3cuYmJiHNJHTk4O/vzzT2zfvh3x8fHYtm0bLl68CAAoWbIkGjVqhCZNmqBJkyYoXry4Q3oUQvw7+T4wHvLe0wCiSV591HsKamCcOXMG0dHR8Pf3x44dO+Dn5+folsxI4tSpU9iyZQs2bNiATZs24dq1a1BKISoqCs2aNUOLFi1Qs2ZN6PV6R7crhHiIfB8YSqlAAH+TpFKqJoBFAML5mOYLYmCQRK1atXD06FHs2bMH5cqVc3RLj2UymbB3716sXbsWa9euxc6dO2EymVC8eHG0bNkSL7zwApo2bQp3d3dHtyqEuCfPB4ZSah6A+gD8AfwN4FMAzgBAcrJS6gMA7wHIAZAJoD/JHY+rWRADY/v27ahduzamTJmCrl27Orqdf+3GjRtYu3YtVq5ciTVr1uDmzZsoUqQImjZtipdeegmtW7eGr6+vo9sUolDL84FhDwUxMHr37o2pU6fi8uXL+X5+hMFgwLZt27B06VIsXboUFy9ehLOzMxo3boxXXnkFbdu2lfAQwgEkMAqI2NhYuLm5YfPmzY5uxaZMJhP27NmDxYsXY+HChThz5gxcXFzQrFkzdOzYEa1bt4aHh4ej2xSiUHhcYMjSIPnI+fPnERoa6ug2bE6n0+HZZ5/F119/jeTkZOzevRsffPABEhMT0bFjRwQEBOD111/HmjVrkJOT4+h2hSi0JDDykdTU1AJ/m0YphZo1a2Ls2LE4e/YsfvvtN7z55ptYs2YNWrRogZCQEPTv3x/79+93dKtCFDoSGPmIr68vUlNTHd1GrtHpdKhbty4mT56MlJQULF26FM899xx++OEHVKtWDVWrVsW4ceNw+fJlR7cqRKEggZGPFC9eHCkpKY5uwyFcXV3Rtm1bLF68GCkpKZg4cSLc3NwwYMAABAcH4+WXX8aaNWtgNBod3aoQBZYERj5SpUoV7N27FwV1oIKlihUrhp49e2L37t04fPgw+vTpg23btqFFixaIiIjAsGHDcPr0aUe3KUSBI4GRj8TGxuLq1as4ceKEo1vJMypUqIBvvvkGFy5cwMKFC1GxYkWMHDkSpUuXRvPmzbF8+XJ5UC6EjUhg5CMNGzYEAKxYscLBneQ9Li4uaNeuHdauXYvk5GQMHToUhw4dQtu2bVGqVCmMGDGi0N7OE8JWZB5GPhMdHQ2lFP744w9Ht5Ln5eTkYNWqVZg0aRLWr18PJycnvPjii3j//fdRt25dWU1XiIeQeRgFyFtvvYWEhAQkJiY6upU8z8nJCW3btsW6detw4sQJ9O7dGxs2bED9+vVRtWpVTJkyxbw8uxDiySQw8pm3334bHh4emDBhgqNbyVfKlCmDsWPH4sKFC5g2bRr0ej26d++O4OBgDBgwQB6SC2EBCYx8xsfHB507d8acOXNQkDeIshd3d3d06dIFe/fuxfbt29G8eXN8//33eOqpp9CuXTvEx8cX+lFoQjyKBEY+NHDgQOh0OnzxxReObiXfUkohNjYW8+bNQ3JyMgYOHIjNmzejTp06iImJwfz582EwGBzdphB5ijz0zqfef/99TJkyBYcPH0bZsmUd3Y6ZyWTCyZMncfDgQRw4cADHjh1DamoqMjIykJ6eDr1eDy8vL3h5ecHPzw8hISEICwtDWFgYIiMjERER4bDNldLT0/Hzzz/ju+++w19//YXQ0FD06dMHXbt2hbe3t0N6EiK3yWq1BdDff/+NMmXK4Pnnn8fy5csd3Q5SUlIwdepUzJw5E8nJyQDuLu3x1FNPoWjRonB3d4e7uzuMRiNu376N27dv4/r160hJSfmf2dlubm4oW7YsKleujKioKERHR6NatWrw9PTMtd/FZDJh9erVGDt2LH777Td4e3uja9eu6NOnT4Fc/FGIB0lgFFBffvklBg8ejF9//RXNmzd3SA8kMWnSJAwaNAjp6el4/vnn0b59e1SvXh0VK1ZEkSJFHvvvc3JykJKSgrNnz+L48eM4cuQIjh49igMHDuDChQsA7gZP9erVUbduXdStWxd16tRB0aJFc+PXQ0JCAsaOHYuFCxdCKYXXXnsNAwcORMWKFXPl+ELktscFBkgWyFdUVBQLuqysLJYvX55hYWG8ffu2Q3oYNmwYAbBJkyb866+/bFo7JSWFq1at4pAhQ1ivXj26uroSAJVSjIqK4sCBA7lhwwZmZWXZ9LgPc/r0afbp04fu7u4EwNatW3P79u12P64QuQ1AAh/xuerwD3Z7vQpDYJBkfHw8AbB37965fuxffvmFAPjOO+/QaDTa/XhZWVnctm0bhw8fzjp16tDZ2ZkA6OnpyZdeeokzZszg5cuX7drD1atXOXz4cBYrVowAWK9ePa5bt44mk8muxxUit0hgFHC9evUiAG7dujXXjmkymVi9enWWK1eOOTk5uXbcB6WlpXHFihXs3r07g4ODCYA6nY61a9fmN998w1OnTtn12N9++635uFFRUVy8eHGuBKcQ9iSBUcClpaXxqaeeYqlSpXjz5s1cOeb9K5spU6bkyvGexGQyce/evfz0009ZtWpVAiAA1qhRg19//TXPnDljl+PeuXOH06ZNY5kyZQiAFSpU4H/+8x8aDAa7HE8Ie5PAKAS2b99OvV7P9u3b58rtkbi4OALg2bNn7X4sLU6ePMkxY8YwKirKHB61a9fmjz/+yCtXrtj8eAaDgXPnzmWlSpUIgGXKlOH06dOZnZ1t82MJYU8SGIXEF198QQCcNGmS3Y81cuRIAsiVB87WSkpK4siRI1mhQgUCoJOTE1u1asVFixbxzp07Nj2W0Wjk0qVLWb16dQJgeHg4J0+ebPPjCGEvEhiFhNFoZLNmzeji4sKEhAS7HuvLL78kAKamptr1OLZkMpl44MABDhw4kCVLliQAFitWjL169eK+fftsfqxff/2Vzz77LAEwLCyMcXFxEhwiz5PAKESuXLnC0NBQRkRE8Nq1a3Y7ztKlSwmAe/bssdsx7CknJ4dr1qxh+/bt6eLiQgCsVq0aJ0yYwBs3btjsOCaTiWvXrmVMTIz5imPq1Klyq0rkWY8LDFlLqoDx9/fHokWLcPHiRbRv395u6yFVrVoVALBjxw671Lc3vV6PZs2aYcGCBUhJScEPP/wAAOjVqxeCg4Px7rvvwhYTP5VSaNq0KXbu3Ik1a9agRIkS6Nq1K55++mnMmjVLdgMU+cujkiS/vwrrFcZ9M2fOJAD27NnTbseoWLEi69WrZ7f6jpCYmMhu3brRw8PDPFx2+vTpzMjIsEl9k8nEVatWmZ9xlC1blvPmzZPhuCLPgNySKpwGDhxIAPz+++/tUn/o0KHU6XQ8ffq0Xeo70s2bN/nDDz+YH5QXLVqUH374IZOTk21S32QyccmSJeZRVZUqVeLSpUtlAqBwuDwfGABmALgM4M9H/FwB+B5AEoCDAKo/qaYExt379G3atKFOp+Pq1attXv/MmTPU6/UcOHCgzWvnFSaTiVu2bGG7du2o1+up0+nYpk0bbty40SYf7kajkfPmzWPZsmUJgDVr1uSGDRskOITD5IfAqAug+mMCowWANfeC41kAu59UUwLjrrS0NFarVo2enp52GTn12muv0cPDg9evX7d57bzm3Llz/Pjjj+nv708ArFKlCmfOnGmTkU8Gg4HTp09naGgoAbBBgwbcuXOnDboW4t/J84Fxt0dEPCYw4gB0fOD74wCCHldPAuO/Lly4wPDwcBYvXtzmCwTu37+fSil++OGHNq2bl2VmZnLGjBnm20lBQUH84osvbBKaWVlZHD9+PIsXL04AbNu2LQ8fPmyDroWwTEEIjFUAaj/w/SYA0Q95XzcACQASwsLCbHwa87fjx4/T39+fpUqV4oULF2xa++2336arqyuTkpJsWjevuz9ktnHjxgRADw8P9u3b1ybLkNy6dYsjRoygt7c3dTodO3XqlGdn1YuCpdAExoMvucL4/+3Zs4ceHh6sVKmSTedonD9/nl5eXmzUqFGhvfe+f/9+vv7669Tr9XRycuKbb77JQ4cOWV33ypUr7N+/P11cXOjm5saBAwfadJ6IEP9UEAJDbknZyMaNG+ni4sJnn32Wt27dslndiRMnEgBnzJhhs5r50ZkzZ9i3b1/zsNyWLVtyx44dVtc9ffo033zzTSqlWLRoUY4bNy5fLMsi8p+CEBgt//HQe8+T6klgPNrSpUup1+tZv359pqen26Sm0WhknTp16O3tbbeVYfOTq1ev8rPPPjPvm9GgQQNu2rTJ6iuw/fv3m2+BlSpVivPnzy+0V3XCPvJ8YACYByAFgAHAeQBdAPQA0OPezxWAiQBOAjj0pNtRlMB4ojlz5lApxaZNmzIzM9MmNU+ePElPT0/Wr1/fYXtk5DVpaWkcN24cg4KCCIC1atXi2rVrrf6QX7duHatUqUIAjImJYXx8vI06FoVdng8Me7wkMJ5s+vTpBMAWLVrYLDTu1xw6dKhN6lnDZDLx0qVL3L17N1etWsXNmzdzz549PHLkiE1vx1kiMzOTP/74o3nYbExMDNesWWNVcOTk5HDGjBnmhRRffvnlQjfwQNieBIZ4pClTphAAmzdvbpPQMJlMfOeddwiAK1assEGH/47RaOTPP//M2rVr083NzbwXxsNeAQEBjI2NZZcuXRgXF8f9+/fbfeOjO3fuMC4ujuHh4ebgsPaKIy0tjZ999hnd3d3p7OzMAQMGyINxodnjAkPd/XnBEx0dTVssHlcYTJs2DV27dkXTpk2xdOlSFClSxKp6mZmZqF27NpKSkrBr1y6UL1/eRp0+ntFoxLvvvouZM2eibNmyaNmyJUqVKoWIiAgEBAQgKysLaWlpuHXrFs6ePYukpCQkJSXh0KFDuHbtGgDAw8MD9erVQ+PGjdGkSROUL18eSimb95qdnY1Zs2Zh5MiROHv2LGJjYzFixAg0bNhQc82LFy9i6NCh+Omnn1CsWDF89tln6NatG5ycnGzYuSjolFKJJKMf+sNHJUl+f8kVxr8zffp0KqXYsGFDmzwIP3PmDAMCAhgZGZkrs8CNRiPfeOMNAuCnn376r/5iN5lMTEpK4pw5c9izZ0/zMh0AGBoayu7du3PFihU2GyDwoDt37nDy5MkMCQkhANavX5+///67VTX37dvH+vXrEwArVqzI9evX26hbURhAbkkJS8yaNYs6nY516tSxycZIv//+O52dndmgQQO7DwFds2aNOSxs4fTp05w6dSpffPFFenp6EgDd3d3ZoUMHrlixwuYbIWVmZnL8+PEsUaIEAbBZs2ZMTEzUXO/+4oalS5cmALZu3drms/xFwSSBISw2f/58Ojk5sUaNGrx69arV9WbPnk0A7Nixo12X8L6/x/i5c+dsXjsrK4vr169njx49zMNk/fz82K1bN/7+++82/b3S09M5ZswYFi1alAD4yiuv8OjRo1b1PmbMGHp5edHZ2ZkffvhhvtolUeQ+CQzxr6xYsYKurq6sVKmSTZYRub+da79+/ew2Z2DatGkEwNu3b9ul/n3Z2dlcvXo1X3/9dbq7uxMAIyIiOHz4cJsu3XHz5k0OGzaMnp6e1Ol07NKlC8+fP6+5XkpKinkwQokSJfjTTz/JHhzioSQwxL+2ceNGenp6slSpUjxx4oRVtUwmE3v37k0AHDlypI06/F/z588nAO7evdsu9R/m9u3b/Pnnn9m4cWMqpajT6diyZUsuX77cZvNQLl++zL59+9LZ2ZlFihThxx9/zJs3b2qut2fPHvM+4zExMfl2i11hPxIYQpM9e/bQ39+fJUqU4N69e62qZTQa+eabb9ptQ6e///7broH0JKdOneInn3zCwMBAAmBYWBhHjx7NK1eu2Kx+x44dCYD+/v6cMGGC5n3BjUYjZ82axRIlSlApxS5duvDy5cs26VPkfxIYQrOjR48yNDSUXl5e3LRpk1W1srOz2aZNGwLg1KlTbdThf9WpU4chISEOXWMpOzubixcv5vPPP08AdHV1ZadOnbhv3z6b1E9ISGCDBg0IgJGRkVbt0peamsoBAwbQycmJfn5+nDhxoszQFxIYwjrnzp1jxYoV6ezszPnz51tVKysri82aNaNSij///LONOrxr3bp1BMAff/zRpnW1Onz4MHv06GF+1lGvXj0uW7bM6mcH9/cFL1++PAGwbt26Vm2OdfjwYXMIVatWTTZuKuQkMITVrl+/ztq1axMAx40bZ1WtjIwMNmzYkEopm65uazKZWLt2bRYvXtymy7db6/r16/z666/Ns7uffvppzpgxw+qhuQaDgZMnT2ZAQAAB8K233tL8YNxkMnHBggUMDg4mAHbt2tUmo+RE/iOBIWwiMzOTL7/8snnEkzV/Kaenp7NJkyYEwMmTJ9usx/3791Ov17Nbt242q2krBoOB8+bNY9WqVQmAwcHBHDdunNUTAlNTUzlo0CC6uLjQ3d2dI0aMYEZGhqZat27d4oABA6jX61msWDFOnz5dRlMVMhIYwmZycnLYq1cvAuBLL71k1YddZmYmW7VqRQAcO3aszXrs378/AfDXX3+1WU1bMplMXLNmDevWrUsALF68OMeMGWP1kOBTp06xXbt25ofuCxYs0Px84+DBg+Yrytq1a9tkMyiRP0hgCJsymUwcN24clVKsUaMGU1JSNNe6c+eO+UNu+PDhNpmnkZGRwSpVqrBYsWJ2mchnS/Hx8WzatCkBsFixYhw1apTVE+u2bt1qvoqpV68e9+/fr6mO0WjkjBkzWKxYMTo5OXHQoEFMS0uzqjeR90lgCLtYtmwZ3d3dGRYWxoMHD2quYzAY2KlTJwJgnz59bHIL5NixY/T09GRMTIzNlm63p127drFFixYEwKJFi/KLL76w6oojJyeHkydPZrFixajT6dizZ0/Nz3WuXLnCzp07mycp5tUrN2EbEhjCbhISEhgUFERPT0+uWrVKcx2j0ci+ffsSAN944w3NcwwetGTJEgLga6+9lm92pduzZ485OIoXL85x48Zpfh5B3n3g3qtXL+p0Ovr7+3Pq1KmaA3nbtm3mkVmvvvqqVVeWIu+SwBB2de7cOVarVo06nY7jxo3T/OFsMpk4atQo8/4ctljmY+TIkQTAYcOGWV0rN+3cuZONGjUiAJYsWZJxcXFW7dVx4MAB1qlThwBYo0YN/vHHH5rqZGVlccSIEXR1daWPjw/j4uLkoXgBI4Eh7C4tLY0vvfQSAbBLly5WDRmdMmUK9Xo9q1evzosXL1rV14MbOn333XdW1XKELVu2MDY2lgBYrlw5LlmyxKpA/s9//sPAwEAqpfjee+9pXnr++PHj5iXU69ata9UCiSJvkcAQucJoNHLIkCHmkTV///235lqrV6+mh4cHw8PD+eeff1rVl8FgMIfZ9OnTrarlCCaTicuWLTPfDqpVqxZ37NihuV5qair79u1LnU7H4sWLc9asWZpCyGQycfr06fTz86OLiwtHjBhh82XfRe6TwBC5av78+SxSpAjDwsKsWoMqISGBgYGB9PHx4YYNG6zqKSsri02aNKFSirNnz7aqlqMYDAZOnTqVQUFB5ucIycnJmuvt27fPvBBh/fr1NV8lXLp0ie3btycAVq5cWRY0zOckMESuS0xMZEhICIsUKWLVciJnzpxh5cqVqdfrOWXKFKt6ysjIYIMGDajT6Th37lyrajnS7du3OWzYMBYpUoSurq4cPHiw5uc9RqORcXFx9PX1pbOzM4cOHap5VNny5ctZsmRJ6nQ6DhgwwC47FAr7k8AQDnHp0iU+99xzBMCBAwdqXtguNTWVzZo1M88wt2aBvLS0NNarV49KKcbFxWmukxecO3fOvC1tUFAQZ8+erfn5xqVLl/j6668TAMuWLcutW7dqqnPz5k12796dAPjUU09priMcRwJDOMydO3f43nvvEQAbN26seX0ig8HAPn36mEdQWbMnRHp6unno6pdffqm5Tl6xc+dO1qhRw/x8w5qtXdevX2/e1rVLly6aH4pv3rzZXKdnz568deuW5p5E7pLAEA43bdo0uri4MDw83KqVVSdPnkwnJyc+/fTTVu1RnZ2dbd5fom/fvvl+aOj9WdkBAQHmEVBaJ+qlp6dz0KBB1Ov1LFGiBBctWqSpTlpaGvv27UulFMPDw61eHl/kDgkMkSfs3r2boaGhdHV1tWq00pYtW1isWDH6+vpy3bp1musYjUbzVcsrr7ySL2aEP8mNGzfYu3dv80S9n376SfNtqn379rFatWoEwBdffFHzEOf4+HhGRkYSAHv06CFXG3mcBIbIMy5fvsyGDRuab3loncV86tQpVq5cmTqdjmPGjLFqbsLYsWMJgLGxsQVm57n9+/eb52/UrVuXR44c0VTHYDDwyy+/pJubG319fTUHUHp6Ovv370+lFCMiIrh582ZN/Qj7k8AQeUpOTg4//vhjAuAzzzzDpKQkTXVu375tHs758ssvW/WX68KFC+nm5sbSpUtr/nDNa4xGI6dNm8aiRYvS2dmZQ4YM0XwVdfz4cfPqtc2bN+fZs2c11dm+fTvLlClDAOzVq5csZpgHSWCIPGnVqlX08/Ojj48PlyxZoqmGyWTiN998Q71ez/Lly1v1Yb9r1y4GBATQx8eHa9eu1Vwnr7l8+bJ5P/XIyEjNI5eMRiPHjx9Pd3d3ent7c/r06ZqvNnr37k0ALFOmDLdv366pH2EfeT4wADQDcBxAEoCPHvLzTgCuANh/7/Xuk2pKYOQPycnJjI6ONg+Z1TpTePPmzQwICKCHh4dVcyxOnz7NKlWqUKfTcezYsflm0UJLbNiwwTxyqVu3bppHmiUlJbFevXoEwJYtW/LChQua6mzZsoXh4eHU6XQcPHiwQ/diF/+VpwMDgB7ASQClAbgAOACgwj/e0wnAD/+mrgRG/pGVlcUPPviAAPjss8/y9OnTmuqcP3/ePO/j/fff1/wBdPv2bfNSIm+++WaBmoCWlpbGAQMGUKfTMTg4WPNS5Uajkd999x2LFClCPz8/zSGdmppqXjq9SpUqVi2TL2zDJoEBoC6AbQAOA5gLoKal//YJdWsBWPfA94MBDP7HeyQwCoGFCxfS29ubfn5+XLZsmaYa2dnZHDBgAAEwKiqKJ0+e1FTHaDRy+PDhVEqxatWqmp+z5FW7d+9mhQoVCIDvvPOO5quN48ePs1atWualSrTOs1mxYgUDAgLo4uLCb775Jt8Pc87PbBUYSQAaAyh+7xbSDgDtLf33j6nbDsC0B75/85/hcC8wUgAcBLAIQOgjanUDkAAgISwszF7nU9hRUlISo6KiCIC9e/fWfJWwdOlS+vr60tvbW/M8AvLuIoj3n7NoDbG8Kisri4MHD6ZOp2NISIjmIco5OTkcPXo0nZ2dGRQUxDVr1miqc/nyZbZt29a8ttWZM2c01RHWsVVg7PrH954A/rT03z+mriWBUQyA672vuwPY/KS6coWRf2VlZZkfilarVk3zBL1Tp06ZZ0D37NlT8wih5ORkc4j17du3wK3IumfPHvNKuD179tQ8cmnfvn2sWLEiAfCDDz7QNGTaZDJxxowZ9PT0pI+Pj1XrkAltrAoMAD8D6AvgGwDDADjd++/Ojyts6cuSW1L/eL8eQOqT6kpg5H/Lly9n0aJF6enpyZ9//llTjTt37rB///4EwKpVq2pekTUrK4u9evUiAEZHR2u+1ZVXZWRksF+/fua1pLSuOJuZmWneObF8+fKaVytOSkoyr6T71ltvWb3PubCctYFRD0BvADMAJAI4DWDjvQfVXz/p31tQ3wnAKQClHnjoXfEf7wl64OZHP+MAACAASURBVOsX/3m187CXBEbBcPbsWfNOcW+88YbmD47Vq1fT39+f7u7unDZtmubRT4sXL6avry+9vLw4Z84cTTW0ysnJ4cqVKzlixAiuXLnSqkUYH2XTpk0MDQ2lk5MTR48erfkY69evZ1BQEJ2dnTl27FhNzyQMBgM//fRT6nQ6lipVijt37tTUi/h3bHJLiv/7AV8JwBu2CIx7NVsA+OteCH1y77+NAPDCva+/wN2H7QcAbAHw9JNqSmAUHAaDgZ999hl1Oh1Lly7NXbt2aapz4cIF8yzzdu3aaV5YLzk52TyL2poQ+zdycnLYsGFDenp6UilFT09PNmzY0C6hcf36dfOEyHr16mmepHf16lXzM4mmTZtq3gM8Pj6e4eHh1Ov1HDVqlF1+Z/FfNg2M/PKSwCh44uPjGRYWRr1ez5EjR2r64DAajfzyyy/p5OTEkJAQbtmyRVMvBoOBw4cPp06nY0REBH///XdNdSy1cuVKenp6EoD55enpyZUrV9rleCaTiT/99BM9PDxYtGhRLl++XHOdSZMm0c3NjQEBAZonRN64cYOvvvoqAbBBgwaa536IJ5PAEAXGjRs32KFDBwJgnTp1NM/Z+OOPP1i2bFkqpThw4EDND7K3b9/O0qVLU6fT8aOPPrLbA/ERI0ZQKfU/gaGU4ueff26X4933119/mRcg7NOnj+bf788//2SlSpXMe6NkZ2f/6xr3H4i7u7vT39+fq1at0tSLeDwJDFGgmEwmzp49m15eXvT29ta8cVBaWhq7detmXtPq0KFDmvq5desW3333XfOD9QMHDmiq8zi5fYXxoAdHrdWoUUNzSGdkZJg3V4qJidFc5+jRo6xatSoBsH///gVu1JqjSWCIAunUqVPmmd3t27fXvP/D8uXLGRAQQFdXV6smja1YsYIlSpSgs7MzP//8c01/RT9Kbj7DeJQlS5aYJ1auXr1ac50FCxbQy8uLfn5+XLFihaYamZmZ5tUBatasyVOnTmnuR/wvCQxRYN2fNObk5MTg4GCuX79eU52///6bbdq0MT/o1foBdPXqVfMts6ioKJtebdwfJfX555/bbZTUk5w4ccL81/2wYcM0h+uJEyfMt7o+/PBDzeG6ePFi+vj40MfHh4sXL9ZUQ/wvCQxR4CUmJponn73//vuaJp/dv0fu5eVFT09PxsXFaR5+u2jRIgYEBNDJyYlDhw4tUAvrZWRksFOnTgTAtm3bal5WPjMz07x9b506dTQ/yD516hRr1qxpXh1AblFZRwJDFAr3J58ppVimTBnu2LFDU53Tp0+bh982adJE8xIVV69e5RtvvGGexFaQlvE2mUz87rvvqNfrWbFiRasmMs6ZM4fu7u4MCAjQvLHSnTt3zM9Zatasqfn5iJDAEIXMg8tmDxo0SNNf90ajkRMnTqS7uzu9vLw4depUzVcbv/76K0NDQ817bd+4cUNTnbxow4YN9PPzo7+/v1VDiw8fPsynn36aOp2OX331lVVXdrZ4zlKYSWCIQic1NdU8cqlixYpMSEjQVOfkyZOsX78+AbBx48aa/3K9ffs2+/btS51Ox8DAQC5YsKDA7LXx119/MTIyki4uLlbNfr916xZffvllwsodFB98zvLJJ5/IRL9/SQJDFFq//vorS5YsSb1ezyFDhmi+2vjxxx/p6elJDw8PTpgwQfPD3oSEBFavXt08+/nEiROa6uQ1165dM2+qNGrUKKv2WP/666+p0+lYsWJFzcvKZ2RkmPfZaNSoUYHZqz03SGCIQu369et8++23CYCVKlXSfLVx+vRpNm3alAD43HPPaV7I0GAwcPz48fTy8qKrqyuHDRumaWXXvCYrK4uvvfYacW9HP4PBoLnW/VtdRYsW5caNGzXXmTZtGl1dXRkaGsrdu3drrlOYSGAIwbt7iN+/2vj44481LXduMpk4a9Ys+vn50cXFhSNGjNA8KufixYvmD9iIiAguWbIk39+mMplMHDx4sHkEldYl5cm7K9ZWrFiRer2eEyZM0FwnISGB4eHhdHFx4bRp0zTXKSwkMIS45/r16+YhoeXLl9e8AuqlS5fM8y0qVKhg1QioLVu2mJfNaNSoEY8cOaK5Vl4xYcIE85wWrbv5kXefa7Ru3ZoA+N5772mer3H16lU2btyYANi9e/cCNczZ1iQwhPiHNWvWMDQ0lDqdjv369dO8adCqVasYFhZGpRS7d++ueQVcg8HACRMm0NfXl05OTuzdu7fmmet5xdy5c+nk5MTq1atrPi/k3QmLH374oTlQtY4yy8nJ4aBBgwiAtWrV0rx6bkEngSHEQ6SmprJHjx4EwNKlS3PTpk2a6ty+fZv9+vWjTqdjiRIlOHfuXM23li5fvszu3btTp9PRz8+P48ePt+kSI7lt9erVdHFxYXR0tNXDiWfMmEFnZ2eWL1/eqqVAfvnlF7q7uzM4OJh//PGHVT0VRBIYQjzG1q1bGRkZSQDs0qWL5r+G9+7da94StnHjxjx+/Ljmng4ePGi+hVK2bNl8/Xxj5cqVdHZ2ZkxMjOahsvdt2bKFfn5+LF68uOZ9UUhy//79DA8Pp5ubW65vhJXXSWAI8QQZGRn86KOPqNfrWaJECc3zJHJycjhhwgR6e3vTxcWFQ4YM0TwCymQycdWqVXz66acJgLGxsYyPj9dUy9GWLVtGvV7P+vXrWz0i7NixYyxdujSLFCnCJUuWaK5z+fJl81Dgjz76SPNQ6YJGAkMIC+3bt49RUVEEwFatWmleFiQlJYWvv/46AbBUqVJWLUNuMBg4ZcoUBgUFEQBfeOEFuyyhbm9z5syhUoqtW7e2asgteffDPiYmhkopjh8/XnOdO3fumJe4b9OmDW/fvm1VXwWBBIYQ/4LBYODYsWPp7u5ODw8Pfvvtt5pnC2/evNm8KGKrVq00T0Qj7+7fMXLkSPr4+FApxQ4dOlh128sRJk6cSADs2bOn1bfY0tPTzVvAfvzxx1ZNFvz++++p0+lYpUoVzX8kFBQSGEJokJyczBYtWpiXKk9MTNRU586dO/zmm2/o6elJV1dXDhkyRPOoLPLurOrBgwfT3d2dOp2OnTp1smrxv9x2f8TTt99+a3WtnJwc8xXCu+++a9UyIGvXrqW3tzdLlChh1fOR/E4CQwiNTCYTFyxYwMDAQOp0Ovbp00fzg9sLFy6Yb1OFhIRw/vz5Vv2VfenSJfbt25dubm7U6/Xs3LmzVVcwucVoNPKll16iTqfj1q1bra5nMpk4ZMgQAmC7du2sWt788OHDLFWqFN3c3PjLL79Y3Vt+JIEhhJVu3LjB9957j0opBgcHc9GiRZo/7Ldt28ZnnnmGAFi3bl3u27fPqt4uXLjAPn360NXVlXq9nm+//Xaev1V169Ytli1blkFBQfz7779tUnPs2LEEwObNmzM9PV1zncuXLzM2NpYAOHr06Hw7Ok0rCQwhbGTXrl3mlVCbN2+u+S/6nJwcxsXF0d/fn0opdu3a1eoPzosXL7J///4sUqQIdTodO3TowP3791tV054OHDhANzc3vvDCCzb7UI6Li6NSivXr17fqAXZmZiY7duxok3Wx8hsJDCFsyGAwcNy4cfT09KSbmxtHjBiheamJGzdusF+/fnRycqKXlxfHjBlj9bIVly5d4sCBA+np6UkAbNGiBbdt22ZVTXv56quvCIALFiywWc3//Oc/1Ov1jI2NtWpZkgfXxWrVqpVVz53yEwkMIezg/PnzfPXVVwmAkZGRXLNmjeZax44dM6+ZVKpUKZvsl3H9+nV+/vnn9Pf3N+9EN3fu3Dw1c9xgMDAqKooBAQFWT+p70MKFC+nk5MSYmBimpqZaVWvSpEnU6XSsUaNGoVgmXQJDCDtav349y5Yta16hNTk5WXOtDRs2sEqVKub1jrRuM/ug9PR0/vDDD+bZ7MHBwRw1ahSvXLlidW1b2LVrFwHw888/t2ndZcuW0cnJibGxsVaH0fLly+nm5sbIyEirliXJDyQwhLCzrKwsfvHFF3R3d2eRIkX42WefaZ7RnJOTw2nTpjEwMJAA+Morr9hk9JPRaOSqVavYqFEjAqCrqys7deqUJ9ZTatOmDX18fGx6lUHe3bJVr9ezbt26Vj0IJ8n4+Hj6+fkxMDAwTz8bspYEhhC55MyZM2zXrp351tKyZcs031q6ffs2P/30U7q7u9PZ2Zm9e/e22S2Rw4cP87333qOHhwcBMDo6mlOnTnXYffotW7YQgFVLfTzKvHnzqJRiixYtrL4dd/jwYYaEhNDHxyfPPheylgSGELls48aNrFChgnkrVq2785F3Rz91796der2eXl5eHDlypM2WsLh58yYnTJhg3o/D29ub3bt3Z3x8fK4OJ71z5w49PT3ZvXt3u9SPi4sjAHbs2NHqNaPOnDnDcuXK0c3NzaolX/KqPB8YAJoBOA4gCcBHD/m5K4AF936+G0DEk2pKYAhHy87O5nfffUcfHx86OTmxX79+Vo3aOXr0KNu0aUMADAgI4Pfff2+zjYBMJhO3b9/ON998k+7u7gTAp556isOGDcu1OR0xMTFs1KiR3ep/+eWXBMDevXtbHYZXrlxhdHQ0nZycOG/ePBt1mDfk6cAAoAdwEkBpAC4ADgCo8I/39AQw+d7XHQAseFJdCQyRV1y+fJldu3alUorFixdnXFycVUtY7Nixg/Xr1ycAhoeH86effrLpPIFbt25x5syZbNiwIZVS5ltW3377rVUP9J8kKiqKLVu2tFt9k8nEfv36mSfkWSs1NZV169alUopTpkyxQYd5Q14PjFoA1j3w/WAAg//xnnUAat372gnAVQDqcXUlMERek5iYyNq1axMAq1atys2bN2uuZTKZuG7dOvPKuuXKleO8efNsvkT3+fPn+c0337BatWoEQACsVKkSBw8ezPj4eJsF1bVr1+jt7c133nnHJvUexWg0mvdRX7RokdX10tPT2bx5cwKwat/xvCSvB0Y7ANMe+P5NAD/84z1/Agh54PuTAPwfUqsbgAQACWFhYTY/kUJY6/7aVGFhYQTAF1980aoRUCaTiUuXLjU/g6hUqRIXL15sl+cPf/31F8eOHcv69etTr9cTAH18fPjyyy8zLi6Ohw4d0nTllJmZyfbt21MpxUOHDtm873/KyspirVq16O7ubvWyLOTd5y/3V80dN26cDTp0rEITGA++5ApD5GUZGRkcOXIkPTw86OLiwv/7v/+zagtTo9HIuXPnmueDPPPMM1aN0HqSGzducOHChezSpQuDg4PNVx9eXl5s0KABe/fuzYkTJ3LTpk1MSkri9evXzVc/GRkZTE5O5u+//87evXuzaNGiBMAxY8bYpdeHSUlJYUhICMPCwmwyHyU7O9s8Ou6rr76yQYeOk9cDQ25JiULr4sWL7Ny5M5VSLFasGCdMmGDV0E+DwcBZs2bxqaeeIgBWq1aNy5cvt+uIJ5PJxOPHj3PWrFns2bMno6OjzcN1H3wppcwP1O+/XFxc+Oqrr1p1e06rPXv20NXVlQ0bNrTJrTWDwWCe+Z+b4WdreT0wnACcAlDqgYfeFf/xnvf/8dD7lyfVlcAQ+cnevXv5/PPPm59HWPshbzAYOHPmzP8JjiVLluTaNqQmk4nnz5/nxo0bOWvWLH777bccOnQo+/fvz9GjR3P69OlctWoVr169miv9PMpPP/1k3oDJFgwGAzt06EAA/Oabb2xSM7fl6cC42x9aAPjr3q2mT+79txEAXrj3tRuAhbg7rHYPgNJPqimBIfIbk8nEFStWsFy5cualz/fs2WNVzfvBUaZMGfMzjnnz5lk1Squg6dy5M3U6HX///Xeb1DMYDGzfvn2+fRCe5wPDHi8JDJFfZWdn88cff2Tx4sUJgB06dLB6/SKDwcA5c+aYJxNGRkZyypQpNpvHkZ/dunWLpUuXZkREhM2WJsnOzjY/CJ86dapNauYWCQwh8qHU1FR+8sknLFKkCJ2dndm3b1+rb+EYjUYuWrTIPBw3KCiIY8aMsXpF1/wuPj6eSim+//77NquZlZXFZs2aUSmVr3bvk8AQIh87f/48u3TpQp1OR29vb44ePdrqhfRMJhM3btxoXojQ29ubgwYN4oULF2zUdf7Tp08fArDpGlHp6el87rnn6OzszHXr1tmsrj1JYAhRAPz555/mPTOCgoI4adIkm+xtkZCQwPbt21On09HZ2ZmdOnXKlfkQeU1aWhrDw8NZrVo1mw4OuHHjBqtWrUoPDw/u3r3bZnXtRQJDiALk999/53PPPWde72nu3Lk2+YA7efIkP/jgA/PQ16ZNm3L9+vWFak/rOXPmEABnz55t07opKSksXbo0ixUrZtVClLlBAkOIAsZkMnHVqlWsXLkyAbBKlSpcsWKFTT7cr169ypEjR7JEiRIEwAoVKnDSpEmFYotSo9HIatWqMTIy0uZDkJOSkhgQEMBSpUpZvX+7PUlgCFFA3Z/hfX/YbExMDDdu3GiT2pmZmZw1axarV69OAPT19WX//v1tsplTXrZgwQIC4NKlS21ee/fu3SxSpAhr1arFzMxMm9e3BQkMIQq47OxsTpkyhSEhIQTABg0aMD4+3ia17y99/uqrr9LJyYkA2KxZM65cubJAzucwGAyMiIhgw4YN7VJ/0aJF5r058uLtPgkMIQqJzMxMjh8/3nw7qWnTplZP/nvQhQsXOHz4cAYFBZmXVx81ahQvXbpks2PkBUOHDqVOp2NKSopd6o8ePZqwwz7mtiCBIUQhk5aWxq+++orFihUjALZq1YoJCQk2q5+dnc1FixaxYcOGBEAnJye+9NJLXLVqlU335nCUw4cPEwAnTZpkl/omk4lvvPEGAXDZsmV2OYZWEhhCFFK3bt3iyJEj6efnRwB84YUXuHfvXpse49ixYxwwYIB5ZnrJkiU5aNAgHjt2zKbHyU0mk4n+/v7s3Lmz3Y6RmZnJqKgoent788SJE3Y7zr8lgSFEIXfz5k2OGDGCvr6+BMA2bdrY9IqDvHvVsXTpUrZq1cq8X0atWrU4efJkXrt2zabHyg2NGjWivT9HTp8+zaJFi7Jy5crMyMiw67Es9bjA0EEIUeD5+Phg6NChSE5OxvDhw/Hbb78hOjoaLVu2xK5du2xyDGdnZ7Rt2xYrV67E+fPn8dVXXyE1NRU9evRAYGAg2rRpg19++QUZGRk2OZ69+fj42L3X8PBwzJkzB4cOHUKfPn3seiybeFSS5PeXXGEI8Wg3b97kqFGjzM84GjZsyC1btth81I7JZOLevXs5YMAAlixZkgDo4eHBjh07cunSpXnmr+qHadasGaOjo3PlWIMHDyYAzp07N1eO9ziQW1JCiIe5ffs2v/76a/OoqtjYWK5atcouwz1zcnK4efNmdu/e3RxUHh4ebNeuHefMmcObN2/a/JhaZWRk0M/Pj6+99lquHM9gMDA2Npa+vr48d+5crhzzUSQwhBCPlZmZyYkTJ5r3Gq9atSrnz59vt3kW2dnZXL9+PXv06MHAwEDzSKuGDRvyu+++c/hD4FmzZhFAru4EeOLECbq7u7NJkyYOnZ8hgSGEsEh2djZnzpxp3sSpTJkydt83w2g0cvv27Rw4cCDLly9v3r61dOnS7N69OxcuXJirO/Pt2bOHPj4+rFy5cq5/cE+cOJEAOG3atFw97oMeFxjq7s8LnujoaCYkJDi6DSHyJaPRiGXLluGLL75AYmIiAgMD0a9fP3Tv3h0+Pj52PfapU6fw66+/YuPGjdiyZQtu3boFAKhSpQrq16+PevXqoWbNmggODoZSyqbHjo+PR6tWrVC0aFFs3boVYWFhNq3/JCaTCc8//zz279+PI0eOoGTJkrl6fABQSiWSjH7ozyQwhBCPQhKbN2/Gl19+iY0bN8LLywvdunVDnz59EBoaavfj5+TkYM+ePdiyZQu2bt2K7du3IzMzEwBQokQJREVFoWrVqqhYsSIqVaqEcuXKwc3N7V8d49q1a5g/fz5mzpyJhIQElCpVyiFhcd+JEydQpUoVtG3bFvPmzcv140tgCCGslpiYiLFjx+KXX36BUgrt27fHgAEDUL169Vzr4c6dO9i3bx8SEhLMr+PHjyMnJ8f8nqCgIERERCAsLAzFihWDn58ffH19odfrkZ2djezsbFy/fh3Hjx/H0aNHcebMGZDEM888g06dOuHNN99E0aJFc+13epihQ4di5MiR2LlzJ5599tlcPbYEhhDCZs6cOYPx48dj6tSpSEtLQ/369dG/f3+0bNkSOl3uT+3Kzs7GX3/9hT///BNJSUlITk7G6dOncfbsWdy4cQM3b96E0Wj8n3/j7u6OcuXK4emnn0b58uXRunVrPPPMM7ne+6OkpaUhMjISpUuXRnx8vM1vvT2OBIYQwuZSU1MxdepUfP/99zh37hwiIyPRp08fdOrUCR4eHo5uz4wkbt++DZJwcXGBi4sL9Hq9o9t6okmTJqFnz55Yu3YtmjZtmmvHlcAQQtiNwWDAkiVLMG7cOOzZswd+fn7o0qULevbsiVKlSjm6vXwrOzsbkZGRCAoKws6dO3PtKuNxgSFLgwghrOLs7IxXX30Vu3btwvbt29GoUSN8++23eOqpp9CmTRts2LABBfUPU3tycXHBRx99hN27d2PHjh2ObgeABIYQwkaUUoiNjcUvv/yC06dP45NPPsHOnTvRpEkTlC9fHt9//z1u3rzp6Dbzlbfeegt+fn4YN26co1sBIIEhhLCDkJAQfP755zh37hxmz54NX19f9OnTB8HBwejatSv27t3r6BbzBQ8PD3Tu3BkrV67ME2ErgSGEsBtXV1e88cYb2LVrFxISEtCxY0fMmTMHUVFRePbZZ/HTTz8hPT3d0W3mae3bt4fBYMCKFSsc3YoEhhAid0RFRWHatGm4ePEixo8fj9TUVHTu3BklS5bEe++9h8TEREe3mCfVqFEDJUuWxJo1axzdimMDQylVVCm1QSl14t7/+j3ifUal1P57L8fHrBBCM19fX/Tu3RtHjhzBtm3b0LZtW8ycORPR0dGoVq0afvjhB9y4ccPRbeYZSik899xz2Llzp6NbcfgVxkcANpGMBLDp3vcPk0nymXuvF3KvPSGEvSilUKdOHcyaNQspKSn48ccfodPp0KtXLwQFBaFDhw5Ys2bN/8ziLqxiYmJw5swZXLt2zaF9ODow2gCYde/rWQDaOrAXIYSD+Pr6mm9L7d27F127dsWGDRvQokULhIWF4cMPP8TBgwcd3abD3F+3KyUlxaF9ODowSpC8fwYuASjxiPe5KaUSlFK7lFKPDBWlVLd770u4cuWKzZsVQthftWrVMGHCBKSkpGDJkiWoUaMGvvvuO1StWhVVq1bFV199hXPnzjm6zVxVvHhxAMDly5cd2ofdA0MptVEp9edDXm0efN+9ddgfNbsn/N7Mw9cAfKeUeuphbyI5hWQ0yej7J1gIkT+5uLjgxRdfxPLly5GSkoIffvgB7u7uGDRoEMLCwlC3bl1MmjQJheGPw/sjyTw9PR3ah90Dg2QjkpUe8loO4G+lVBAA3Pvfh8YnyQv3/vcUgK0Aqtm7byFE3uHv74/3338fO3fuxIkTJ/D555/j6tWr6NmzJ4KCgtC4cWPExcUV2PD4+++/Adw9D47k6FtSKwC8fe/rtwEs/+cblFJ+SinXe1/7A3gOwJFc61AIkaeUKVMGQ4YMweHDh3HgwAEMGjQIZ86cQY8ePRAYGIgGDRpgwoQJOH/+vKNbtZnNmzfD398f4eHhDu3DoYsPKqWKAfgFQBiAMwDak7yulIoG0IPku0qpWABxAEy4G3DfkZz+pNqy+KAQhQdJHDx4EIsXL8bixYtx5MjdvymjoqLQunVrvPDCC3jmmWdydZlwW8nKykJgYCBefvllTJ/+xI8+q8lqtUKIQuXYsWNYtmwZVqxYgV27doEkSpYsiebNm6NFixZo1KgRvL29Hd2mRYYMGYJRo0Zh69atqFevnt2PJ4EhhCi0Ll++jNWrV2PNmjVYv349UlNTodfrUatWLTRu3BiNGzdGdHQ0nJ2dHd3q/2fZsmV48cUX8c4772DGjBm5ckwJDCGEwN29O3bu3Im1a9diw4YNSExMBEl4eHggNjYW9erVQ926dREdHY0iRYo4rE+TyYTRo0dj2LBhiIqKwrZt23KtHwkMIYR4iGvXrmHLli347bff8Ntvv+HQoUMA7u7xUb16dcTGxiI6OhpRUVGIjIy0+xa0JJGQkIDBgwdj06ZNeP311xEXF5erOxhKYAghhAWuXbuGHTt2YPv27di+fTsSEhKQlZUFAPDy8kLlypVRsWJFVKxYERUqVECZMmUQGhoKJycnzcfMzs7GwYMHsWnTJsyePRuHDx+Gr68vxowZg65du+b6g3oJDCGE0MBgMODIkSPmJUsOHTqEP//8E9evXze/x9nZGREREQgJCUHJkiURFBQEf39/eHt7w8vLC+7u7jCZTDCZTMjJycHVq1dx6dIlXLp0CUePHsW+fftw584dAEBsbCzeeustdOjQAT4+Pg75nSUwhBDCRkji77//xrFjx3Dy5EkkJSXh1KlTuHjxIi5evIgLFy6YA+BRnJycUKJECZQuXRoxMTGIiYnBs88+i5CQkFz6LR7tcYGh/TpKCCEKIaUUAgMDERgYiPr16/9/PyeJzMxM3L59G7du3UJGRgZ0Oh30ej10Oh38/f1RtGhRuz8PsQcJDCGEsCGlFNzd3eHu7o4SJR61nmr+lP8iTgghhENIYAghhLCIBIYQQgiLSGAIIYSwiASGEEIIi0hgCCGEsIgEhhBCCItIYAghhLCIBIYQQgiLSGAIIYSwiASGEEIIi0hgCCGEsIgEhhBCCItIYAghhLCIBIYQQgiLSGAIIYSwiASGEEIIi0hgCCGEsIgEhhBCCItIYAghhLCIBIYQQgiLSGAIIYSwiCLp6B7sQil1BcAZR/eRB/gDuOroJvIQOR//S87Hf8m5uCucZPGH/aDAs5zQeAAAA4tJREFUBoa4SymVQDLa0X3kFXI+/pecj/+Sc/FkcktKCCGERSQwhBBCWEQCo+Cb4ugG8hg5H/9Lzsd/ybl4AnmGIYQQwiJyhSGEEMIiEhhCCCEsIoFRQCilmimljiulkpRSHz3k53WVUnuVUjlKqXaO6DE3WXA++iuljiilDiqlNimlwh3RZ26w4Fz0UEodUkrtV0rFK6UqOKLP3PKk8/HA+15WSlEpJUNt75FnGAWAUkoP4C8AjQGcB/AHgI4kjzzwnggA3gD+D8AKkotyv9PcYeH5aABgN8kMpdR7AOqTfNUhDduRhefCm+Ste1+/AKAnyWaO6NfeLDkf997nBWA1ABcAH5BMyO1e8yK5wigYagJIInmKZDaA+QDaPPgGkqdJHgRgckSDucyS87GFZMa9b3cBCMnlHnOLJefi1gPfegAoyH9FPvF83PM5gDEAsnKzubxOAqNgCAZw7oHvz9/7b4XVvz0fXQCssWtHjmPRuVBKva+UOgngKwC9c6k3R3ji+VBKVQcQSnJ1bjaWH0hgiEJNKfUGgGgAXzu6F0ciOZHkUwAGARji6H4cRSmlAzAOwABH95IXSWAUDBcAhD7wfci9/1ZYWXQ+lFKNAHwC4AWSd3Kpt9z2b/+/MR9AW7t25FhPOh9eACoB2KqUOg3gWQAr5MH3XRIYBcMfACKVUqWUUi4AOgBY4eCeHOmJ50MpVQ1AHO6GxWUH9JhbLDkXkQ982xLAiVzsL7c99nyQTCXpTzKCZATuPt96QR563yWBUQCQzAHwAYB1AI4C+IXkYaXUiHujXqCUqqGUOg/gFQBxSqnDjuvYviw5H7h7C8oTwMJ7w0kLZMBaeC4+UEodVkrtB9AfwNsOatfuLDwf4hFkWK0QQgiLyBWGEEIIi0hgCCGEsIgEhhBCCItIYAghhLCIBIYQQgiLSGAIIYSwiASGEEIIi0hgCJGLlFLtlFK7lFIH7u09UdzRPQlhKZm4J0QuUkoVI3nt3tefArhKcqKD2xLCInKFIUTu6qSU2qOUOgCgJ2S/BZGPODm6ASEKC6XUW7i7gc/zJNOUUtsAFNg1vUTBI1cYQuSeygB23AuLlwHEAjjk4J6EsJg8wxAilyilKgJYAiAVwHoA7UmWdWxXQlhOAkMIIYRF5JaUEEIIi0hgCCGEsIgEhhBCCItIYAghhLCIBIYQQgiLSGAIIYSwiASGEEIIi/w/NK/Ff+WpKwgAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Use the chi-squared method to create joint confidence intervals of model parameters a,b fit to a dataset\n", "# Create grid of values of the model parameters a,b\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", "# (Note that we are not performing a *minimisation* for each set of model parameters, just an evaluation)\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", "minchisq2 = np.amin(chisq2)\n", "print('Minimum chi-squared =',minchisq2)\n", "# Find the values of the model parameters at the point of minimum chi-squared\n", "i,j = np.unravel_index(chisq2.argmin(),chisq2.shape)\n", "a2,b2 = alst[i],blst[j]\n", "print('Best-fitting value of a =',a2)\n", "print('Best-fitting value of b =',b2)\n", "# These values of delta chi-squared correspond to the 1-sigma, 2-sigma and 3-sigma confidence regions\n", "clst = np.array([0.,2.30,6.17,11.8])\n", "# Add the minimum chi-squared to these intervals to obtain the chi-squared contours for the plot\n", "clst += minchisq2\n", "# Make a contour plot of these chi-squared values\n", "alst2,blst2 = np.meshgrid(alst,blst,indexing='ij')\n", "plt.contour(alst2,blst2,chisq2,levels=clst,colors='k')\n", "# Indicate the best-fitting model parameters\n", "plt.plot([a2],[b2],marker='o',markersize=5,color='black',linestyle='None')\n", "plt.xlabel(r'$a$')\n", "plt.ylabel(r'$b$')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can translate this to errors in **individual parameters** by determining for each parameter $a$, the minimum value of $\\chi^2$ varying parameter $b$:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "68% confidence region = 0.18673057259713693 -> 0.2925882056451612\n", "Error in a = 0.05292881652401213\n" ] }, { "data": { "text/plain": [ "Text(0, 0.5, '$\\\\chi^2_{\\\\rm min} \\\\; ({\\\\rm varying} \\\\; b)$')" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEKCAYAAADuEgmxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deZzNdf//8cfLzBi7CkkokUiXfWwtLhV+KntdUYopNRFFKNvMZafoii6yjCg0EemyJLJkScqaXclVrrKracpuxrx+f8wR+ZJzOGfe55x53W+3czv75/PsY5rnfLb3R1QVY4wxxhvZXAcwxhgTOqw0jDHGeM1KwxhjjNesNIwxxnjNSsMYY4zXrDSMMcZ4zUlpiEgOEVkjIptEZJuI9PO8fouIrBaRXSLygYhkd5HPGGPMxbla0zgF3KeqFYFKQAMRqQm8BgxX1VuBX4G2jvIZY4y5CCeloRmOep5GeW4K3Ad86Hl9EtDUQTxjjDGXEOlqxiISAawHbgXeAv4LpKhqmucje4Cil/huHBAHkDt37qply5YNfGATHL79NuO+TBm3OYKRLRtznl27dnH06FHKly9PRETE/3l//fr1P6tqIV+n66w0VPUMUElErgH+A3j9m19VE4FEgJiYGF23bl1gQprgU6dOxv2yZS5TBCdbNsZj7dq1VK9enQEDBhAfH3/Rz4jI/65k2s6PnlLVFGApUAu4RkTOFlkxYK+zYMYYE6ISEhIoWLAgnTp18vu0XR09VcizhoGI5ATqATvIKI9HPB9rA8x2kc8YY0LV559/zqeffkqPHj3Imzev36fvavNUEWCSZ79GNmC6qn4sItuBaSIyEPgamOAonzHGhBxVpVevXhQpUoT27dsHZB5OSkNVNwOVL/L690D1zE9kjDGhb/78+axcuZLRo0eTK1eugMzD+T4NY4wxVy89PZ1evXpRsmRJ2rYN3Cluzo6eMsYY4z/Tp09n06ZNJCUlkT174AbTsDUNY4wJcampqSQkJFChQgVatmwZ0HnZmoYxxoS4d955h127djF37lyyZQvsuoCtaRhjTAg7fvw4/fv358477+Shhx4K+PxsTcMYY0LYyJEj2bt3L1OnTkVEAj4/W9MwxpgQlZyczJAhQ2jUqBH33HNPpszTSsMYY0LUkCFDOHLkCIMHD860eVppGGNMCPrxxx8ZOXIkbdq04W9/+1umzddKwxhjQtA///lPAPr165ep8w2L0jh48KDrCMYYk2m2bNnC5MmTefHFFylevHimzjvkSyM5OZmbbrqJHTt2uI5ijDGZ4pVXXiF//vz06NEj0+cd8qWRL18+cuTIQffu3V1HMcaYgFu8eDELFiwgPj6e6667LtPnH/KlERkZSa9evZg7dy7L7Iplxpgwlp6ezssvv0yJEiXo2LGjkwwhXxrAH9v1unXrRnp6uus4xhgTEO+99x4bN25k8ODBREdHO8kQFqWRM2dOBg0axPr165k2bZrrOMYY43cnTpwgPj6eqlWr0qJFC2c5wqI0AFq1akXlypXp2bMnJ0+edB3HGGP86s033+Snn37i9ddfD/ighH8lbEojW7ZsvP766/z444/8+9//dh3HGGP85tChQwwePJhGjRpRp04dp1nCpjQA7rvvPho2bMigQYM4fPiw6zjGGOMXffr04cSJEwwbNsx1lPAqDYBhw4Zx7Ngx+vbt6zqKMcZctW3btpGYmEj79u0pU6aM6zjhVxply5alffv2jBs3ju3bt7uOY4wxV6Vbt27ky5ePPn36uI4ChGFpQMaqXJ48eXj55ZddRzHGmCv26aefsmDBAhISEihQoIDrOECYlkbBggWJj4/nk08+YeHCha7jGGOMz9LS0ujWrRulSpWiQ4cOruP8ISxLA+CFF17glltuoWvXrqSlpbmOY4wxPhk/fjxbt25l6NChzk7ku5iwLY3o6GiGDRvG1q1bGT9+vOs4xhjjtZSUFBISEqhTpw7NmjVzHedPwrY0AJo3b87f//53EhISSElJcR3HGGO80r9/f5KTkxk+fHimXPfbF2FdGiLCiBEjSE5Opn///q7jGGPMZe3cuZORI0fStm1bKlWq5DrO/+GkNESkuIgsFZHtIrJNRDp5Xu8rIntFZKPn9uDVzqtSpUq0bduWkSNH8u233159eGOMCaCuXbuSM2dOBg4c6DrKRbla00gDuqpqOaAm0EFEynneG66qlTy3T/wxs4EDB5IzZ066du3qj8kZY0xALFy4kI8//pj4+HgKFy7sOs5FOSkNVd2vqhs8j48AO4CigZpf4cKFSUhIYN68eSxYsCBQszHGmCuWmppKp06dKFWqFJ06dXId55Kc79MQkRJAZWC156WOIrJZRCaKyLX+ms+LL75I6dKl6dSpE6dPn/bXZI0xxi9GjRrFN998w/Dhw4PqENsLOS0NEckDzAQ6q+rvwBigFFAJ2A/86xLfixORdSKyztuBCaOjoxkxYgQ7d+60UXCNMUHl4MGD9O3blwYNGtCwYUPXcf6Ss9IQkSgyCiNJVT8CUNWDqnpGVdOB8UD1i31XVRNVNUZVYwoVKuT1PB988EEaNmxIv3792L9/vx/+K4wx5ur17t2b48ePM2LEiKA7xPZCro6eEmACsENV3zjv9SLnfawZsNXf8x4+fDinT5+mR48e/p60Mcb4bO3atUycOJHOnTsHxSi2l+NqTeMu4EngvgsOrx0qIltEZDNwL/CSv2d866230qVLFyZPnsyXX37p78kbY4zX0tPTeeGFF7j++utJSEhwHccrkS5mqqorgYutg/nlENvL6d27N5MnT6Zjx46sWbOGiIiIzJitMcb8yTvvvMPq1auZNGkS+fLlcx3HK86PnnIhT548vPHGG2zYsIHExETXcYwxWVBycjI9evTg7rvv5sknn3Qdx2tZsjQAHn30Ue677z569epll4Y1xmS6+Ph4kpOTGTVqVNDv/D5fli0NEWHUqFEcPXrUdoobYzLV+vXrGTt2LB07dqRixYqu4/gky5YGwO23306XLl2YOHGi7RQ3xmSK9PR0OnbsSKFChejXr5/rOD7L0qUBkJCQQNGiRXn++eftYk3GmICbMGECX331FUOHDuWaa65xHcdnWb408uTJw4gRI9i4cSNvvfWW6zjGmDB2+PBhunfvzt///ndat27tOs4VyfKlAfDwww/ToEED4uPj2bt3r+s4xpgw9fLLL3PkyBFGjx4dUju/z2elwbmd4mlpaXTu3Nl1HGNMGFqxYgWTJk2iW7dulCtX7vJfCFJWGh6lSpUiPj6eDz/8kPnz57uOY4wJI6dPn6Z9+/aUKFEiZM78vhQrjfN069aNsmXL0qFDB44fP+46jjEmTLzxxhts376dkSNHkitXLtdxroqVxnmio6MZM2YMP/zwQ9BeatEYE1r++9//0q9fP5o3bx70w557w0rjAnXq1CE2NpZhw4axefNm13GMMSFMVWnXrh1RUVFhcx0fK42LeP3117nmmmuIi4vjzJkzruMYY0JUUlISixcvZsiQIRQtGrArWmcqK42LKFCgACNGjGD16tWMHTvWdRxjTAj65ZdfeOmll6hRowbt2rVzHcdvrDQu4fHHH6d+/fr07NmTPXv2uI5jjAkx3bp1IyUlhcTExLC6/IKVxiWICGPGjCEtLY0OHTqgqq4jGWNCxKJFi3j33Xfp1q0bFSpUcB3Hr6w0/kLJkiXp378/c+bMYcaMGa7jGGNCwLFjx4iLi+O2226jT58+ruP4nZXGZXTu3JmYmBg6duzIL7/84jqOMSbIJSQksHv3bsaPH0+OHDlcx/E7K43LiIyMZMKECfz666+89JLfL1lujAkja9as4c0336Rdu3bUrl3bdZyAsNLwQoUKFejZsydTpkyxIUaMMRd1+vRp2rZty4033shrr73mOk7AWGl4qXfv3tx+++0899xz/P77767jGGOCzKBBg9i6dStjxowhX758ruMEjJWGl6Kjo5k4cSJ79+7l5Zdfdh3HGBNENm7cyODBg3nyySfDYqiQv2Kl4YOaNWvStWtXEhMTWbRokes4xpggkJqaSmxsLAULFmTEiBGu4wSclYaP+vXrR5kyZXjmmWdsM5UxhiFDhrBp0ybGjRvHdddd5zpOwFlp+Chnzpy888477NmzxzZTGZPFbd68mQEDBtCqVSsaN27sOk6msNK4ArVq1aJLly4kJiaycOFC13GMMQ6cPn2a1q1bU6BAAd58803XcTKNlcYV6t+/P2XLluXpp58mJSXFdRxjTCbr378/mzZtYvz48RQoUMB1nExjpXGFcubMyeTJkzlw4ACdOnVyHccYk4lWr17NkCFDeOqpp2jUqJHrOJnKSWmISHERWSoi20Vkm4h08rx+nYgsEpHvPPfXusjnrWrVqtG7d28mT57MrFmzXMcxxmSCEydO0KZNG4oWLcrw4cNdx8l0rtY00oCuqloOqAl0EJFyQA9giaqWBpZ4nge1+Ph4qlSpQlxcHIcOHXIdxxgTYL169eLbb79l4sSJ5M+f33WcTHdFpSEiuUXkigeIV9X9qrrB8/gIsAMoCjQBJnk+NgloeqXzyCxRUVFMnjyZ33//nbi4OBtC3ZgwtnjxYkaMGEGHDh2oW7eu6zhOeFUaIpJNRB4XkXkicgj4Btjv2bw0TERuvdIAIlICqAysBgqr6n7PWweAwpf4TpyIrBORdYcPH77SWfvNHXfcweDBg5k9ezYTJ050HccYEwDJycnExsZStmxZhg4d6jqOM96uaSwFSgE9gRtUtbiqXg/cDXwFvCYiT/g6cxHJA8wEOqvqn86U04w/2S/6Z7uqJqpqjKrGFCpUyNfZBkTnzp2577776NSpE7t27XIdxxjjR6pK+/btOXjwIO+99x65cuVyHckZb0ujrqoOUNXNqpp+9kVVTVbVmar6MPCBLzMWkSgyCiNJVT/yvHxQRIp43i8ChMxOgmzZsjFp0iSioqJ44oknSEtLcx3JGOMnSUlJTJ8+nX79+lG1alXXcZzyqjRUNdUfnzlLRASYAOxQ1TfOe2sO0MbzuA0w29tpBoNixYoxduxYVq9ezaBBg1zHMcb4we7du+nQoQN33XUX3bt3dx3HOZ92hItIbRFZ4TlM9n0RqX6F870LeBK4T0Q2em4PAq8C9UTkO6Cu53lIadGiBU888QQDBgxg1apVruMYY65CWloaTzzxBKrKlClTiIi44uN/wkakj5+fCLQHNgJVgREiMkJVp/syEVVdCcgl3r7fx0xBZ9SoUXzxxRe0atWKjRs3ZsnD8owJB4MGDeKLL74gKSmJW265xXWcoODrIbc/q+oiVT2sqguA+sA/A5ArpOXPn5+kpCR++ukn2rVrZ4fhGhOCvvjiC/r378+TTz7J448/7jpO0PD2kNvJItIZWCki/xSRs2sop4CTAUsXwmrVqkXfvn2ZNm0aU6ZMcR3HGOODlJQUWrVqRYkSJXjrrbdcxwkq3q5pTADSgevIOAFvl4gsJuN8jaUByhbyevbsSe3atenQoQPfffed6zjGGC+oKs899xx79+5l6tSp5M2b13WkoOLVPg1VXQ4sP/vcs6ZRFqgEVAxMtNAXERHBe++9R8WKFWnZsiWrVq0iOjradSxjzF94++23mT59Oq+++irVq1/psT7h64qGEVHVNFXdqqrvqapdiegvFC9enIkTJ7JhwwZ69Aj6obSMydK2bdvGiy++SL169ewia5dgQ6NngqZNm/LCCy8wYsQI5s6d6zqOMeYijh8/TosWLcifPz9TpkwhWzb79XgxtlQyybBhw6hcuTKxsbHs2bPHdRxjzAU6d+7Mtm3bmDJlCoULX3TYO4OP52mISJeLvPwbsF5VN/onUniKjo5m2rRpVKlShccee4ylS5cSGenraTLGmEBISkpi/Pjx9OjRg3r16rmOE9R8XdOIAdqRMYx5UeA5oAEwXkRe8XO2sHPbbbcxfvx4Vq5cSe/evV3HMcYAO3bs4LnnnuOee+5hwIABruMEPV9LoxhQRVW7qmpXMs4Kvx6oDcT6OVtYeuyxx2jXrh1Dhw5lzpw5ruMYk6UdO3aMf/zjH+TKlYupU6fa2r8XfC2N68k4oe+sVDKugXHigtfNXxg+fDhVqlShTZs27N6923UcY7KsDh06sH37dpKSkihatKjrOCHB19JIAlaLSB8R6QusAt4XkdzAdn+HC1c5cuRgxowZqCr/+Mc/OHXK+taYzPb2228zadIk4uPjbT+GD3wqDVUdAMQBKUAy8Jyq9lfVY6raKhABw1XJkiV59913WbduHS+++KLrOMZkKevWraNDhw7Ur1+fPn36uI4TUnwdGj0auA3IDVwDPCgiNmDhFWratCk9evQgMTHRLhNrTCb5+eefefjhh7nhhhtISkqy4c595Oten9l4DrHF9mH4xYABA1i7di3PP/88lSpVokqVKq4jGRO2zpw5Q6tWrThw4ABffPEFBQsWdB0p5PhaGsVUtUFAkmRRkZGRTJ06lapVq9K8eXPWr19PgQIFXMcyJiz16dOHhQsXkpiYSExMjOs4IcnXHeGrRKR8QJJkYYUKFWLmzJns37+fFi1a2PXFjQmADz/8kEGDBvHMM8/wzDPPuI4TsnwtjbuB9SLyrYhsFpEtIrI5EMGymmrVqjFu3DiWLFnCK6/YeZLG+NOWLVuIjY2lZs2ajBo1CpFLXTjUXI6vm6ceCEgKA0BsbCxff/01w4cPp1KlSrRu3dp1JGNCXnJyMk2bNiVfvnzMnDnTLk9wlXwqDVX9X6CCmAyvv/46W7ZsIS4ujttvv51q1aq5jmRMyEpLS6Nly5bs2bOH5cuXc+ONN7qOFPK8vdzrSs/9ERH5/bzbERH5PbARs5aoqCimT59OkSJFaNq0Kfv27XMdyZiQ1aVLFxYtWsTo0aOpWbOm6zhhwavSUNW7Pfd5VTXfebe8qpovsBGznoIFCzJ79mx+//13mjRpwvHjx11HMibkjBs3jpEjR/LSSy/Rtm1b13HChq8n970gItcEKow5p0KFCrz//vusX7+e2NhY0tPTXUcyJmQsXbqUjh078sADDzBs2DDXccKKr0dPFQbWich0EWkgdghCQDVq1IihQ4cyY8YM+vXr5zqOMSFh165dPPLII5QuXZqpU6faGd9+5uvYU/FAaWACGUOhfycig0WkVACyGaBr16489dRT9O/fn6SkJNdxjAlqv/zyCw8++CDZsmVj7ty55M+f33WksOPz5V5VVYEDnlsacC3woYgM9XM2A4gIY8eO5d577+Xpp59m+fLlriMZE5ROnTpF06ZN+fHHH5k9ezalStnfsoHg6z6NTiKyHhgKfAGUV9X2ZFyM6eEA5DNA9uzZmTlzJiVLlqRZs2Z8++23riMZE1RUlaeffpqVK1cyadIk7rzzTteRwpbXpeHZf1EBaK6q/09VZ6hqKoCqpgMNfZjWRBE5JCJbz3utr4jsFZGNntuDPvx3hL1rr72WTz75hKioKB588EEOHTrkOpIxQSMhIYH333+fwYMH06JFC9dxwprXpeHZLFXjUif4qeoOH+b7LhnXFr/QcFWt5Ll94sP0soRbbrmFuXPnsn//fho2bMjRo0ddRzLGubFjx/4xplSPHj1cxwl7vu7TWC8iV32KsqquIOMiTsZH1atXZ9q0aaxfv55HH32U1NRU15GMcWb27Nl06NCBhx56iDFjxtiYUpnA19KoAXwpIv8N0ICFHT3TnSgi117qQyISJyLrRGTd4cOH/Tj70NC4cWPGjh3L/PnziYuLI2Ml0Jis5csvv6Rly5bExMTwwQcfEBnp61B65kr4upT/X0BSZBgDDADUc/8v4OmLfVBVE4FEgJiYmCz5G/PZZ59l79699OvXjyJFijB48GDXkYzJNNu2baNhw4YUK1aMjz/+mNy5c7uOlGX4PGChZw2gNJDjvLeueiBDVT149rGIjAc+vtpphrs+ffqwf/9+hgwZQoECBejatavrSMYE3O7du6lfvz7R0dF8+umnFCpUyHWkLMWn0hCRZ4BOQDFgI1AT+BK472qDiEgRVd3vedoM2PpXnzcZ53CMHj2aX3/9lW7dunHttdfy9NMXXTkzJiwcPHiQevXqcfz4cVasWEHJkiVdR8pyfN081QmoBnylqveKSFnA5+0iIjIVqAMUFJE9QB+gjohUImPz1G7gOV+nmxVFREQwZcoUfvvtN5599lmuueYamjdv7jqWMX7322+/0aBBA/bt28fixYspX94uIuqCr6VxUlVPiggiEq2q34hIGV9nqqqPXeTlCb5Ox2SIjo7mo48+ol69ejz22GPMnj2bBg3sUu4mfBw9epQHHniAbdu2MWfOHGrVquU6Upbl69FTezyj3M4CFonIbPywP8Ncvdy5czNv3jzKlStHs2bNWLp0qetIxvjFiRMnaNy4MWvWrGHq1Kn2B5Fjvg5Y2ExVU1S1L5BAxtpB00AEM7679tprWbhwISVLlqRRo0asWrXKdSRjrsqpU6do3rw5y5YtY9KkSTz8sI1W5JqvY091EZGiAKq6XFXnqOrpwEQzV6JQoUIsXryYG2+8kQceeIC1a9e6jmTMFTl9+jQtWrRgwYIFjB8/nlatWrmOZPB981ReYKGIfC4iHUWkcCBCmatTpEgRlixZQoECBahXrx7r1q1zHckYn5w+fZpHH32U2bNnM2rUKLvyXhDxdfNUP1W9A+gAFAGWi8jigCQzV6V48eIsXbqU6667jrp161pxmJBxdg3jbGF06NDBdSRzHp+vp+FxiIzrafwCXO+/OMafbr75ZisOE1LOFsasWbMYOXKkFUYQ8nWfxvMisgxYAhQAnlXVCoEIZvzjwuL46quvXEcy5qJOnjxJ8+bN/yiMjh07uo5kLsLXNY3iQGdVvUNV+6rq9kCEMv518803s2zZMgoWLEi9evX4/PPPXUcy5k+OHz9OkyZNmDdvHuPGjbPCCGJelYbnAkyoak9V3fhXnzHB6aabbmLFihUUK1aMBg0asGTJEteRjAEyTtxr2LAhixYtYuLEicTFxbmOZP6Ct2saS0XkBRG56fwXRSS7iNwnIpOANv6PZ/zpxhtvZNmyZZQsWZKHHnqIjz+2MSGNW8nJydSrV4/ly5czZcoUnnrqKdeRzGV4WxoNgDPAVBHZJyLbReR74DvgMWCEqr4boIzGjwoXLszSpUspX748TZs2JSkpyXUkk0UdOHCAOnXqsGHDBj788EM7DyNEeDX2lKqeBEYDo0UkCigInFDVlECGM4FRsGBBlixZQpMmTXjiiSf49ddfbRuyyVS7d++mXr167N+/n3nz5lG3bl3XkYyXfD7kVlVTVXW/FUZoy5cvH/Pnz6dx48a88MIL9O3b164AaDLFpk2buPPOO/n5559ZvHixFUaIudLzNEwYyJEjBzNnziQ2NpZ+/foRFxdHWlqa61gmjC1dupTatWsTERHBypUrqVmzputIxkd2Ud0sLjIykokTJ1KsWDEGDhzI/v37+eCDD+zymcbvpk+fzpNPPsmtt97KggULKF68uOtI5gpc9ZqGiPQSkfdFZKqIvO+PUCZziQgDBgxgzJgxzJ8/n3vvvZcDBw64jmXChKry2muv0aJFC6pXr87KlSutMEKYPzZPiao+rqqPqerjfpiecaRdu3b85z//Ydu2bdSoUYMtW7a4jmRCXGpqKs8++yw9evSgZcuWLFq0iGuvvdZ1LHMV/FEapUTkHyLyoIg86IfpGYcaN27M559/TlpaGnfddRfz5893HcmEqF9//ZUHHniACRMmEB8fT1JSEjly5HAdy1wlf5TGciAXUIiMQ3FNiKtSpQpr1qzh1ltvpWHDhgwfPtyOrDI+2bFjBzVq1GDFihW88847DBgwgGzZ7LibcHBV/4oichdw2HP72XMzYaBo0aKsWLGCZs2a0aVLF9q0acOJEydcxzIhYN68edSsWZOUlBQ+++wzYmNjXUcyfnS11X8dGWsYZ9cybE0jjOTJk4fp06fTv39/pkyZQu3atdmzZ4/rWCZIpaenM2TIEBo1akSpUqVYt24dd999t+tYxs+uqjRUdS6QBGwCdgDf+COUCR7ZsmUjISGB2bNn8+2331K5cmUWL7brbpk/S0lJoVmzZvTq1YsWLVqwcuVKbrrppst/0YQcf2xknE7G2FT3APZnRZhq3Lgxa9asoXDhwtSvX5+BAweSnp7uOpYJAps2bSImJoZPPvmEN998k/fff59cuXK5jmUCxB+lsU1VX1XVf6nqG36YnglSZcuWZfXq1Tz++OMkJCTQsGFDDh8+7DqWcURVSUxMpGbNmpw4cYJly5bx4osvYldJCG+XLQ0RqSci40Wkkuf5hYPdp4rIIs8JfnZyX5jLnTs3U6ZMYfTo0Xz22WdUrFiRzz77zHUsk8lSUlJo0aIFzz33HLVr12bDhg3cddddrmOZTODNmsbTwMvAEyJyH1DpgvdvUNV6nhP87OS+LEBEaN++PatXryZ//vzUrVuXXr16kZqa6jqayQRffvkllStX5qOPPuLVV19l/vz5FC5c2HUsk0m8KY0jqpqiqt2A+kC1C97PJSIt7eS+rKdixYqsW7eOp59+miFDhnDnnXfyzTd2LES4Sk1NJSEh4Y8joj7//HO6d+9u519kMd78a88773FfYPIF7y8Fojl36O1lichEETkkIlvPe+06z2au7zz3NtZACMidOzdvv/02H374IT/88AOVK1fmzTfftJ3kYWbHjh3UqlWLgQMH0qZNGzZt2kStWrVcxzIOXLY0VHX2eU9XAcsueH/S+Tcv5/suGUdcna8HsERVSwNLPM9NiHj44YfZunUr999/P507d+b+++9n165drmOZq5SWlsarr75K5cqV2b17NzNnzmTixInky5fPdTTjiK/rlc8BSSLS7cI3RGSBtxNR1RVA8gUvNwHOls4koKmP2YxjN9xwA3PnzuXtt99mw4YNlC9fntdff92u0RGiNm/eTM2aNenZsyeNGjVi27ZtNG/e3HUs45hPpaGqa4EaQBURWSIiz4vIOBHZ4uu0LqKwqu73PD4AXHLPmojEicg6EVlnh3wGFxGhbdu2bN++nfr16/Pyyy9Tq1Yt1q9f7zqa8dKxY8fo3r07VatW5bTtBMAAAA05SURBVKeffmLGjBnMmDHDdnYbwMdf9CIyENgKVAAOAglAbqCBqtb3VyjNGB3vkiPkqWqiqsaoakyhQl7tRjGZrGjRosyaNYsPPviAPXv2UK1aNTp27EhKil0lOJjNmTOHcuXKMXToUFq3bs327dt55JFHXMcyQcTXtYNYoJqq/s1zeG0FIC/wuohc7UbOgyJSBMBzf+gqp2ccExEeffRRvvnmGzp27MiYMWMoU6YMkyZNsh3lQea7776jYcOGNGnShLx58/L5558zYcIEChQo4DqaCTK+lsZtqvrHvghVPayqTcjYOf7VVWaZA7TxPG4DzP6Lz5oQkj9/fv7973+zdu1abrnlFmJjY6lVqxZr1qxxHS3LO3LkCN27d+eOO+5gxYoVDBs2jK+//toGGjSX5Os+jeOXeH0c0Njb6YjIVOBLoIyI7BGRtsCrQD0R+Q6o63luwkiVKlVYtWoVkyZN4scff6RGjRq0bt2a//3vf66jZTlpaWmMHTuW0qVLM3ToUFq1asXOnTvp1q0bUVFRruOZIOa3s3JU1evjKz2Xhi2iqlGqWkxVJ6jqL6p6v6qWVtW656/RmPCRLVs2Wrduzc6dO+nevTvTp0+nTJkydOvWjeRk+ycPNFVlzpw5lC9fnvbt21O6dGlWr17NO++8ww033OA6ngkBdiqncSJv3ry8+uqr7Ny5k5YtW/LGG29QsmRJ+vXrx2+//eY6XthRVT799FNq1qxJkyZNUFVmzZrFihUrqF69uut4JoRYaRinbrrpJt599102btxInTp16Nu3LyVKlGDgwIFWHn5wtizuueceGjRowMGDB3n77bfZsmULTZo0sRFpjc+sNExQqFChArNmzWL9+vXcc889JCQkULx4cV555RX27dvnOl7IOXPmDNOnT6dq1ao0aNCA3bt3M3r0aHbu3Enbtm1tv4W5YlYaJqhUqVKFOXPmsGHDBh566CH+9a9/UaJECWJjY+0EQS+kpaXxxhtvcNttt9GiRQuOHz/OhAkT+P7772nfvj3Zs2d3HdGEOMk4jy50xeTNq+uqVnUdwwTIiZMn2fPTTxw4cIAz6elUjYggOnt2ImJiiLDRVf9w5OhRIrZs4fTp02wk4zDnYsWKUbBAAdsEZS5Kli9fr6oxvn4vMhBhjPGXnDlyULp0aW655RYOHDyIfv89x0+cYOuqVVxfuDBFbriBPHnzkhV/LaampXHo0CH279/P0aNHqQRkj4qiaoUK5M2Tx3U8E6ZCvzTKlIFly1ynMAEWCRQDtE4dUn79lbHlyzNz5kxO7ttH2bJladmyJS1atKBs2bKuowbUsWPHmDt3LlOnTmX+/PmkpqZSsWJFnnnmGe6aNo2oyEj7/8F45wrXQEN/81RMjK5bt851DJNZ6tTJuF+2jJSUFD744AOmTZvG8uXLUVX+9re/0bhxYxo3bky1atXC4gJBhw8fZu7cucyePZtFixZx4sQJbrzxRlq2bMnjjz9OlSpVMjZBnbdsjLkcEbmizVNWGia0XOIX4759+5g+fTqzZs1i5cqVnDlzhuuvv57777+funXrcv/993PzzTdnetwrceLECb788ksWLVrEokWL2LBhA6pK8eLFadKkCc2bN6d27dpERET8+YtWGsYHVhoma/DiF2NycjLz58/nk08+YcmSJRw8eBCA4sWLU6tWLWrVqkWNGjWoWLEiuXLlCnzmv6Cq/PDDD6xfv56vvvqKL774gg0bNpCamkpkZCS1atWibt26NGrUiEqVKv31Tm0rDeMDKw2TNfj4i1FV2bZtG5999hmrVq1i1apV/PTTT0DGkCa33347lSpVoly5cpQtW5ayZctSokQJv5fJmTNn2LdvH7t27WLHjh1s376dbdu28fXXX/9xEmOOHDmoVq0ad955J/fccw+1a9cmb9683s/ESsP4IOuWhh1ym7Vs3JhxX6nSFU/i1KlTHDlyhCNHj3L0yBGOHj3KqdOn//SZqKgockRHkz06mqioKKIiI4mMiiIiIoKIbNky9pWc91e/pqdzJj2dM2fOcObMGVJTU0k9fZrTqamcOnWKUydP/ukCMREREeTOlYs8efKQJ29e8ubJQ+48ech2NYfH+mHZmKzDDrk1xkvR0dFER0dTsGDBP15LO3OGE8ePc/z4cU56fsmfvT965Aipqamk+/AHVlRkZEbZZM9Ovnz5yFGoEDly5CBHzpzkzpWL7NHRWfIwYRP6Qr807JDbrCVAm2Aiybia2KU2BqkqJ0+e5LinWE6cOEF6evof+xiioqLInTs3uXPnJleuXG6O2rLNU8YXV7hWG/qlYUwmEBFy5sxJzpw57Wp2JksL/YPYjTHGZBorDWOMMV6z0jDGGOM1Kw1jjDFes9IwxhjjNSsNY4wxXrPSMMYY4zUrDWOMMV6z0jDGGOM1Kw1jjDFes9IwxhjjNSsNY4wxXgu6AQtFZDdwBDgDpF3JeO/GGGMCI+hKw+NeVf3ZdQhjjDF/ZpunjDHGeC0YS0OBhSKyXkTiXIcxxhhzTjBunrpbVfeKyPXAIhH5RlVXnP8BT5nEAdx0000uMhpjTJYUdGsaqrrXc38I+A9Q/SKfSVTVGFWNKVSoUGZHNMaYLCuoSkNEcotI3rOPgfrAVrepjDHGnBVsm6cKA/+RjAueRwLvq+oCt5GMMcacFVSloarfAxVd5zDGGHNxQbV5yhhjTHCz0jDGGOM1Kw1jjDFes9IwxhjjNSsNY4wxXrPSMMYY4zUrDWOMMV6z0jDGGOM1Kw1jjDFes9IwxhjjNSsNY4wxXrPSMMYY4zUrDWOMMV6z0jDGGOM1Kw1jjDFes9IwxhjjNSsNY4wxXrPSMMYY4zUrDWOMMV6z0jDGGOM1Kw1jjDFes9IwxhjjNSsNY4wxXrPSMMYY4zUrDWOMMV6z0jDGGOM1Kw1jjDFes9IwxhjjtaArDRFpICLfisguEenhOo8xxphzgqo0RCQCeAt4ACgHPCYi5dymMsYYc1ZQlQZQHdilqt+r6mlgGtDEcSZjjDEeka4DXKAo8NN5z/cANS78kIjEAXGep6dEZGsmZAsFBYGfXYfIFCKX+0TWWRYX+r/LJusui//LlsU5Za7kS8FWGl5R1UQgEUBE1qlqjONIQcGWxTm2LM6xZXGOLYtzRGTdlXwv2DZP7QWKn/e8mOc1Y4wxQSDYSmMtUFpEbhGR7EBLYI7jTMYYYzyCavOUqqaJSEfgUyACmKiq2y7ztcTAJwsZtizOsWVxji2Lc2xZnHNFy0JU1d9BjDHGhKlg2zxljDEmiFlpGGOM8VpIlMblhhYRkWgR+cDz/moRKZH5KTOHF8uitohsEJE0EXnERcbM4sWy6CIi20Vks4gsEZGbXeTMDF4si3YiskVENorIynAeacHboYhE5GERUREJ20Nwvfi5iBWRw56fi40i8sxlJ6qqQX0jY4f4f4GSQHZgE1Dugs88D4z1PG4JfOA6t8NlUQKoAEwGHnGd2fGyuBfI5XncPov/XOQ773FjYIHr3K6WhedzeYEVwFdAjOvcDn8uYoFRvkw3FNY0vBlapAkwyfP4Q+B+kcufMhyCLrssVHW3qm4G0l0EzETeLIulqnrc8/QrMs77CUfeLIvfz3uaGwjXI2C8HYpoAPAacDIzw2WygAzLFAqlcbGhRYpe6jOqmgb8BhTIlHSZy5tlkVX4uizaAvMDmsgdr5aFiHQQkf8CQ4EXMylbZrvsshCRKkBxVZ2XmcEc8Pb/kYc9m3A/FJHiF3n/T0KhNIy5KiLyBBADDHOdxSVVfUtVSwHdgXjXeVwQkWzAG0BX11mCxFyghKpWABZxbovNJYVCaXgztMgfnxGRSCA/8EumpMtcNszKOV4tCxGpC/QGGqvqqUzKltl8/bmYBjQNaCJ3Lrcs8gJ/A5aJyG6gJjAnTHeGX/bnQlV/Oe//i7eBqpebaCiUhjdDi8wB2ngePwJ8pp69PGHGhlk557LLQkQqA+PIKIxDDjJmFm+WRenznj4EfJeJ+TLTXy4LVf1NVQuqaglVLUHGvq7GqnpFg/cFOW9+Loqc97QxsOOyU3W9h9/LowAeBHaScSRAb89r/cn4xwbIAcwAdgFrgJKuMztcFtXI2HZ5jIy1rW2uMztcFouBg8BGz22O68wOl8WbwDbPclgK3OE6s6tlccFnlxGmR095+XMxxPNzscnzc1H2ctO0YUSMMcZ4LRQ2TxljjAkSVhrGGGO8ZqVhjDHGa1YaxhhjvGalYYwxxmtWGsYYY7xmpWGMMcZrVhrG+JmIPCIiX4nIJs+1Kwq5zmSMv9jJfcb4mYgUUNVfPI/7AD+r6luOYxnjF7amYYz/xYrIGhHZRMYFwsL5mg0mi4l0HcCYcCIircm4+M19qnpURFaQMbaPMWHB1jSM8a/ywCpPYTwM3AlscZzJGL+xfRrG+JGI3AF8RMbVIxcCj6rqbW5TGeM/VhrGGGO8ZpunjDHGeM1KwxhjjNesNIwxxnjNSsMYY4zXrDSMMcZ4zUrDGGOM16w0jDHGeO3/AweO8dvytpc9AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Convert the joint chi-squared distribution for 2 parameters into individual errors in 1 parameter\n", "# Minimise the value of chi-squared for each different value of parameter a, varying parameter b\n", "chisq1 = np.amin(chisq2,axis=0)\n", "# Find the value of parameter a corresponding to the minimum value of chi-squared\n", "i = chisq1.argmin()\n", "# The chi-squared interval corresponding to the 68% confidence region varying 1 parameter is, delta chi-squared = 1\n", "# Find the values of a corresponding to chi-squared = minimum chi-squared + 1, by interpolation\n", "abot = np.interp(minchisq2+1.,np.flip(chisq1[:i+1]),np.flip(alst[:i+1]))\n", "atop = np.interp(minchisq2+1.,chisq1[i:],alst[i:])\n", "print('68% confidence region =',abot,'->',atop)\n", "print('Error in a =',0.5*(atop-abot))\n", "# Plot the minimum chi-squared as a function of parameter a, marginalising over parameter b\n", "# Indicate the 68% confidence range of parameter a, corresponding to chi-squared = minimum chi-squared + 1\n", "ymin,ymax = 0.,30.\n", "plt.plot(alst,chisq1,color='black')\n", "plt.plot([amin,amax],[minchisq2,minchisq2],color='red')\n", "plt.plot([amin,amax],[minchisq2+1.,minchisq2+1.],color='red')\n", "plt.plot([abot,abot],[ymin,ymax],color='red')\n", "plt.plot([atop,atop],[ymin,ymax],color='red')\n", "plt.xlim(amin,amax)\n", "plt.ylim(ymin,ymax)\n", "plt.xlabel(r'$a$')\n", "plt.ylabel(r'$\\chi^2_{\\rm min} \\; ({\\rm varying} \\; b)$')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Activity\n", "\n", "In this Activity we will use a recent supernova distance-redshift dataset to determine the cosmological parameters $(\\Omega_m,\\Omega_\\Lambda)$. First let's read in the data and use astropy to compute a fiducial model:" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Text(0, 0.5, '$DM$')" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEGCAYAAABsLkJ6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3dfXxcZZn/8c/VaQgtLSnQlDaUzggEqSAWGxF5mojyE4Efuov4sLhrXZHaVEVRUH/Lmsni7qKsC4tLefBhcXdVdF1AF2SVVRJgESRAeSgstmqCGB6i2Kwttqbp9fvjzJmemcwkkzbzlPN9v17zymTmzJxrpul9nXPf17lvc3dERCR+ZtU6ABERqQ0lABGRmFICEBGJKSUAEZGYUgIQEYmp2bUOYCoWLlzoqVSq1mGIiDSUBx988Nfu3lr4eEMlgFQqRX9/f63DEBFpKGY2WOxxdQGJiMSUEoCISExVLQGYWcLMHjazW7O/m5n9tZn91MyeNLMPVysWERGp7hjABcCTwL7Z31cBBwNHuPtOM1tUxVhERGKvKmcAZrYUOAP4UuThNcBfuftOAHd/oRqxiIhIoFpdQFcCFwM7I48dCrzDzPrN7HYzay/2QjM7P7tN//DwcDViFRGJhYonADM7E3jB3R8seKoZ2ObuHcAXga8Ue727X+/uHe7e0do6roxVRER2UzXGAE4AzjKz04G9gX3N7F+BZ4CbstvcDPxTFWIREZGsip8BuPun3H2pu6eAdwI/cvd3A7cAr89ulgZ+WulYRGRmyGQymFnulslkah1SQ6rllcCXAV8zs48CW4DzahiLiDSQTCZDb28vQO6nTF1VLwRz9153PzN7f7O7n+Hur3T317n7I9WMRUQE4n020VBzAYmITLc4n01oKggRkZhSAhCR2Nu4cSN9fX2sXbu21qFUlRKAiMTe0NAQANddd12NI6kuJQCpS402MNdo8Uq+trY2AFavXl3jSKrM3RvmtnLlSpf4SKfTnk6nax1G2Rot3kbX1tbmgHd1de3xe5X6t+vu7nYgd+vu7p7w8XoF9HuRNrXmjfpUbkoA8dJoDWq9xjuVxqqRGrYwRjPb7ZjL+byl/l0r8u+9dav7Y4+533yz+9/9nfsHPuB+6qnuhxzi/swzu/22pRKAykBFZriplDk2UklkW1sbQ0NDrFmzhg0bNgDjY85kMvT09OR+7+7uzuueCz/v+vXrWbBgAT09PbntC7edNlu2wKZNwW3jxuDnz34W/PzVr/K33W8/OPRQOPZYGB2d9lCUAESkrpVqxNvb22lvb+fqq6+ms7Oz5GsnS2gDAwOMjIwwMjICQDKZZGBgIPf8xo0bGRoaYu3atVx99dVFH29tbc2L8dL/9//Y/8UX+eG119IOtANvXLaM5B/+AM89lx/AokXQ3g5vfGPQ2Le3Bz8PPRT237/cr2m3KAGISEVNdhRezusna8TDxvg1r3kN/f39efuazGj2yLqrq4tbbrmFwcHBvMY+WiEUTQDPDw1xKPD0tddy9ec/z4q2Ng5+6SVWzp8Pf/M3wXuGGx94ICxbFjTu7e1w2GG7Gvr588v+LqabEoCI5Cl1xLu7qtGtFDbSDz/8MOl0Om9fBx100ISfJ3ztunXrco9FG/sjDzyQBc8/z8vHxvicGYcDJyxcyEvAXgA7d8JHP8obZs/ml3PmQDoNhx8Ohx/O+z/3OX41Zw7fu+eeinzuPaUyUBHJ04g18ROVcU72edra2kgAJy5axHv235+LgLva2+GEE+CAA3j8+ee5B/gy8BEz3nrkkbSedBL/dvDBXPbyl8M998DwMEe0tnLk737H2n33hb/8S3jHO9g4fz4vza7f42wlAGlYqr2vjEasiW9vbyedThc9ws/7PFu3wkMPwde+BpdcAmefzT2bN/MScPcLL3DDiy/yOeD4zZthr73gnHP4x0MP5eKjjuJ1ixYxx5216TTcdBNfPOQQ/nPx4iBRLFzI0LPPAvmJpt6vMK7f1CQyiUaqWAntaX94NUQHVxtBJpOhr68PADMjmUxy5NKlcN998MQT/P2sWST335/jbrsNIt08JBL8esECHnnpJb4JPAn8D/AUMPLcc3RnE2DPz36Wt7/rrruO1tbWvH12d3fnqpKiibPU+EG9UAIQqaJaJ61aJqCBgQHMLO+xPd3/3B07yLzpTWSWLoUnnoANGxju7aV1cBBe9zoAziJo1L/+4oss7uzklLVrYflyaG9n4V578VbIxZVIJNixY0fePsL4wrGE1atXk8lkxsXd29s7LnEWSwr1RAlAJEZqlYDCgeW2tja2bt3KihUryt5/JpPh/r4+lgPvMeP8172ORcPD/POmTSwDOP74YMM5c2D5cvoSCfqB5Omns+aqq5iTSrEikWDFBPuYrKHOZDJ5g8Wtra1lJa56P5tSAhCJgYGBAQYHB3NHuslksujZAJDXtRF9bk+O1MPGM/wZ9ouPaxjHxoKLoh59FB57DB5/nA/19vJpdg1Y7vjJT5j9ylfCuefCkUfuuqVSkEhwTng0//3vs+bQQ8uKb7KGutgRf+HzhV1C9da1V1Sxy4Pr9aapIOKlnEvt62n6hXLnppnOmMuduiGdTntLS4un02mfN29e3muSyWTF4wy/m+htf/BO8A+Dfxn8f+bNc58zJ5ihBtxnzXI//HD3s8/2T4P/MfiZ7e3uo6Nl7WsqcwRV6u+oXv4+0VQQIpVViQG/wqP0qPAoc6pdOlu2bMndD2vmi+132o5ox8Z4w0EHcVhLC/v98pccumULRwNLo9u0tsLRR++6vfKV8IpXBN06wF+FR/U//zlMUlY51W6XSh29N8RZQbGsUK83nQHESzlHcvVyhOU+fWcAhUf1yWQyt32x15Z7phSeAUSPxtva2vJeW7hvSpxZlDzz2LrV/b773K+5xn31avfXvjbvqH7UzB+fNcu/v2iR++WX+5sTCV9UxncW7ueoo46acLtyv4+4ocQZgAXPNYaOjg6PXuYtM9tElRmhcA6YeigDLTeWcraLblPqfiisTunq6ip51BtuA7sGPEPpdDrv/cJ9hEevUdGj2DNPOIH2LVu44k//FB5+OLg99VRwZSzAggWwYgWsWMFlt9/OpnnzuGt4mI1PPz3ufSf6N4by/hYK46+Hv4l6YWYPuntH4eO6EEzqViNekDRV5VzMNtnFRBNd6Rq+f7TBHxoaoqmpCdj1HU9mIfCFM88ks9de8La3wSGHcOu993LFo4/CRRfB3XcHc9tccgncfDMMDMCLL8Kdd8IVV/CfixfzX7/+dV7j393dnWvYly9fXnLf0e9kbGys/rpRGljVEoCZJczsYTO7teDxq8xsS6nXSXxNdHVnqBJXWtbbFcbRBr7Y5w0b8WOOOWZc3JlMhnQ6nWvww+3nzp2b6/8vfL99R0fpePFFPgX8OzAADAMfvPVW+Iu/gEcegY4O/nr+fE4FLn7ve+Hpp+E734GeHjLr12OpFDZrVt73l0ql8rofMpkMYQ/Ek08+WfLzZzIZWlpaaGlpyb1OpkmxfqFK3IALga8Dt0Ye6wD+BdhSzntoDCBeyunLJds3nEgkqr7v3X1NYRVOc3Ozt7S05G0THU+I3i/2eaP7LTVG0NLS4i0tLZ5MJvP2vQ/4yeAXmbm//e3uL3vZrioc8I3g3wD/GPgVb3mL+29/m3vfib77aBwTVSo1NTWVNW4Sxl8OjQGMRy1XBCMY8P8hcEqYAIAEcCewRAmg9upxJahy/iNP57KAU9337r4m+j2772oEJ7qFpZrFPu9ECSD675oAP33JEv/uGWf4l8AfBd8Raew9mXQ/+2z/6/nz/RTwBZMMBk/03Zfz7zKVEtbodlMZQK+Hv+N6UOsE8G1gJdAZSQAXAB/N3lcCqAP1duRUKp5iVTLV2veexlKqyqZUQx9W6ZTTwAN5ZxcHgZ8Nfjn4XeBbI439MPht4N3g71m0yP3553PvWW5jO9F3VO6ZWb39zc1UpRJAxa8DMLMzgRfc/UEz68w+1gacQ5AQJnv9+cD5AMuWLatcoNIworXvtdp/tDY/mUySSqXKfu2VV17Jli1bmDdvHitWrMhV27S1teXep729neHhYYaGhnJjANErc8P+/dwVqtu3w0MP8bHjj+c44HXsqrPfBjwEXAf8bvlyPn3rrbztve+F6Lw8ixbl7oZVQuUOEBdT73PgSKAag8AnAGeZ2QBwI0E30AbgMGBT9vG5Zrap2Ivd/Xp373D3jtbW1iqEKxIIJy8rHAwOB1bT6TTuzujoKH19fWzcuLHo+xQOKm/btg2Al156Kdf4JxIJIH9ANlypKpFI4O4kk8nce17b08NfHnkkf2fGvWZs33tvOP54Pk8wsPbAnDlcdeihfPG889iX4D/hhUD3k09ihx7KvT/+ccmY29vbSSaTDA0N0dfXV3QgPLzIqdTz5QzgSx0odlpQqRuRLqCCx9UFVAfq7XS81IAmkS6K8OKm6VbYBTNRbJTRZRLdvrm5edI+/+ita80a9yee8M8dfrjffuCB7occkuvK2W7mj+67r/tFF7nfdJN3LF06rg88vN/U1FRWzNPxd1DOe6i/vnqo5RhAbmdKAHWtERKAe1ARkkgkcs9PJeZyG53JGvViVTphspion7+jo8NbWlryHjMzB3zu3LkO+IJ99vFjCSpvbsn214cN/otNTe5vfav75Ze733uvv/Gkk/LiS6fT4xJM+P7RZFlYiRT9Hvb070ANe/2piwSwpzclgMpqpAQQNi6ljtB3532jJjsDCPdfmIgKyyzDaRyi2xdW/TQTlGJeAn4H+JZIg/8UwURpq8APK9KgFlbbpLMln4lEItfgF54xhdsVJorwfevt70D2nBKATKoa//GncnRYTgIodoQ+0T6KPVfssXDfpcoZw8c7OjqKNvhhI9/W1pY3905XV5c3g6cJKnDuBP99trEfA38Y/B/A3wZ+YOT9woa98LNGE0vh52hubp606qic70canxKATKpaR37l7mdPzgAm2kexi4oKtw9/jzawE20/UUPbBH4i+F8WNPg7wB8gKNM8E3y/Iq+NnokUSwDFElRhwil2ncFEyU1mnlIJQHMBSUNKJBKTVqpMReEUCwMDA3mToR1zzDG5+4UVMGaWKwvt7u5mFvBq4CLgduC3wN1ABpgPXA2cCewPfDyd5iLg1ux2UdE5fKJz+UQVq7aJbhudHygsIw3LOyeaQ0jiQQlAGs62bdsYGxtjcHBw2t6zsDFMpVK5+WfS6TQPPPBArpyzcH7+ZDLJO489Fr/2WjIbNjAMPAh8DjgY+DJwTiLBQuDkuXP5OHAb8L8Un3EzbKjXrFmTa6zb2tpyZaOTzXsUrd9fvXp1rpzU3Umn07S3t+dtp1r9+FICkJoqZ+K1wm3CBq2pqSmvHn9PzgAKJ1Tr6+tjZGSEbdu25a4HCBv+ZDLJ6SecwPmtrawDfjg4yDd+8hP4wAfgvvu4b/FiLj3iCI458ECOIrjk/dtjY/yWoPY/KpFI0NLSkvdYcMYerD0bfsb29vbc5w6TVKla/LCBD4UTwUUnhAu3Kzx7kJgp1i9UrzeNAVRWtcr/SvW3T/Re0UVRKOjLnqgvvtggbzj5WuHjhYPL4W3evHneefLJfv6rX+2fnT/f7wYfzfbj/y/4d8A/CP6KWbPcd+7MqwoKYwv722fPnp333sXKQot97sJ+/YkUViJN9L6q9okHNAgsk5mOQcFyGpVwm92tCCrWoBVOjlY4y2T09dFB4GKlm+Ftf/B3gn8VfHjWrOC/S3bg9jPgJ2UHeKMxlJobKNx/4b66u7vzrmuIKjWwPJHC13R0dJRMuEoA8aEEIJMKG409mVp5KgkgvD/Z1bzllDJGYy7nCDhssAsXKz8a/FPg97BrpswXwP8F/E/AW8m/uCq8eKutrW3Sz1kq2SQSiUkTYbQaqBzF9h19rlgikpmrVALQovCSU6sJvKKDm8X6ozOZDJdeeik7w6UGixgbG8ubLK2URCLB2NgYEPSlN42N8Wbg/xJU5hyc3a4f+Id587jNHV+5kqc2bcpNkNY0PMzo6ChtbW20t7ezfv36cf3uhQuCTySMByi5cHh7eztbt24dt5/JFFuYHIJB7oGBgSm9l8w8GgSWnGoNChaWXBYObhYTbfyjg6ZtbW25Sp10doWriYyNjbEQeA9wx/z5/Br4HvCnwAPAnxMsUPEa4LsrV3L3H/7AnXfdBZArOw3jDUtQR0ZGxg3EhqtdTXYLG+RKKYxDq2lJlBKAlG26lkosLLkMq1MmOvOYNWvXn+rIyAhmllfSCBMvAp4CPgL0As8BNwArRke5ATiNYM3bs4F/yj4PQYlmtLEPy07nzZsHQFdX17gGfarfiRpoqSUlAJlQtNHv6ekhmUzucdllYf152LivW7euZGIJp0tuamqipaUlLyFs27YtdwQe9QrgEoK58H8BXAEsAD5DcKHWAVu38kHg+8D2yOuamppyDXIYa7SxX7lypconZWYo5zS1Xm4aBK6sUgO4xQYUi1W7RH9fsmRJWfsJB4EnGwiGyQd1W/bd118Ffin4E9kB3DHwu8EvBH9ZwQBu4Tw+hbfonEChWsyVszv7nGgwXlNAxA+qApLJTLUev3D7aD37RJVExerjw7nqSzV20flsCqtmVs6e7Z8B32jmnq3e+S/wNeCLJ0ka4fsXVgNBcI1BtRv7UqaaBCZKAOX8G8nMUioBqApIJhQuYQiwefNmOjs7y3pduZVE4XhAuKpWdGqEaEXM3LlzGRkZAYKB3MOBd2Zvy3fsYAy4053PATcDv55kv4lEgh07duR+b29vz8XS1dVVd907uaUfp4GWa5SQxgBknFJLGJYjkUiQTqdpbW3Ne4/wlkql8qYvCAdUC7W0tOQaPTNjZGSEpcDHCObZeQroJhiw/QCwGDgV+CJB49/U1JQbXG5ra8u7H44nREVLItetWzejB2M1BYSEdAYg40QXXe/t7WWvvfZiZGRk0knIJtLc3JxLJAsWLADyzyg2btyYN4vlPvvsQyaT4YqeHt5LUKaZJjhieWDWLD6ycyffAkbmzh03vw7sKi2F/Nkxh4aGMLO8QeRMJpM3sVypWvxGUaz2v5E/j1SOEoBMKlqnf+KJJ5bcLpyls7ArB2Dvvfcmk8nkzaSZyWQYGBgYN6tnAjh6aIgjenp4DpgD/JRgOuWvAz+LXhBWpPGPCo/+o/tw97zrCsLGMYwt/NmojeZ0dhfJzKYuIJlUtE6/8CKuqL333hsIum8KL3Aq1o3U09PD4OBgrktmOfBZ4JcE8+i/kWAq5WOBlwOXAj+bYuxDQ0OkUqlxF165+4QXbqkBlVgoNjJcrzdVAVVWqfljopOnUWJOHfddlTqzZs0qul2xssv54O8Hvy9btvkH8JvB30L+RGvl3sKVs4qtghVWzhRbESxONAlc/KAVwaQc4dz3YTdOJpPhpZdeyo0BhBdGRfvQwwu4wq6iUnP29Pf35+4fR3B0/yxwPbAPcCFwEPBHwHeA0fFvMU7hxV9DQ0OMjIwwOjpKMpnMPZ5IJHJH9dHPIxJnSgCSJ5VKkU6nc/PrZDKZvDGAsDuosJEvXNSkmI+9732sBR4BfgycA3wNeC3wSoIrdYfLjDO6albhUU0YfyqVGjfNRPTzzPRqH5FJFTstqNebuoAqp/BCo7Abp7DbptRUyxNNwbwC/HrwLdlunp+Anwc+b4pdO+EVvOGtubm56GcpFkuxRdPjeiWsuoDih1p3AZlZwsweNrNbs79/zcyeMrPHzewrZtY02XtI5UQHQdPpdO4I/8EHHyzr9YWVPE3Au4D/Bh4GzgW+AawkGNT9ErBlCvENDQ0R/B0Hmpubc4POhcKzmOj20ZlGVQcvEqhmF9AFwJOR378GHEFw9j8HOK+KscgENm7cmLu/cuXKKb12EfBp4GmCks2FBP/wbUB3WxsPTfDa6Bq/hV1Kzc3Neb9v3749NytoqW6c6ONjY2Pq7hEpUJUEYGZLgTMIDvwAcPfvRU5PfgIsrUYsMl7hlb/RC6eiA7cw/kg/dDTwFYKGv4fgat3TCDL8VcAI4xclh6BhT6fTNDc3500HsWDBgrzSze3bd83XGS0xjQ7uFvtc6WlaNF5kRirWLzTdN+DbBGf/ncCtBc81EczYe1KJ155PsEBT/7Jly6a3Y0xySk3Q5r6r9LNwQXPA3wR+R7Zvfwv4F8Bfvhvlm+FtyZIlZU16Vqoff6LJ64p91ripxWymUnvUajZQgpX21mXvF0sAXwSuLOe9NAhcOWGDOtkUyYDPBn83+KPZhv8Z8IvBF0zwmrDRTSQSuQXQo414uF25M1SW04iXauzinAAknkolAPPIQFklmNnfEkzlsgPYG9gXuMnd321m3cAxwB+7e+kFX7M6Ojq8sEtCpkdYVmlmlPqbmAO8D/g4kAQeBy4nGNwtp2Z/1qxZRa8RSCQSHHjggQwNDZU9E2c4h9BEq4BV4rUijcjMHnT3jsLHKz4G4O6fcvel7p4imL33R9nG/zzgTcC7ymn8ZfoU9vlHp3gu1vjPAy4mWFXrCwRTNZxJ0O//zwSNf3iBGIy/OCvsh482/u6etzJYtSpzwonSCtfwFYmjWl4Idi1wIPBjM1tvZp+uYSwzRjnr9kYHR7u7u8dN3BaaD/wFMEgwR8964GTgJOA2gn6VUHTgeM2aNXnTL4fC++HPWpRjas4fkV0q3gU0ndQFVJ7wiL6zszNv9s3otMDz589ny5bilfjzCEo3LwT2B/4D+CuCkfjm5ua8ipzpkkwmSaVSZXfLqBtHpHw16wKS6op2cUQbfwhm3wzPDKKNf3hEPodgwZVfECycfg9B6dZZBI0/MK7xTyQSeReQhSWX0ekkij1feEulUhX5PkSkNCWAGSaTyeS6X7q6umhpack1xGHjW3iR1QtDQ5wPbAL+jqCxPxZ4q9mEF25BbS6wUj++yPRQAmhwxfr8o5OdjYyMMDIywr333ptrMMO1dQHOJqjmuY7gyP9k4M3AAxQfEC5UqgEeGBhgZGSkIo20+vFFpocSQIMrdrVrdKrmUHSJRAimY/5vgiv0Rgm6eU4E7i5zv+FZhLsXnVY5Oh/PZI20juhFaqRYf2y93nQhWHHpdHrC2Tijt2Xg38hewDUE/ufgswou2CqcdRMmnu0zfF3h9tGLrXQFqkjtUKsLwaaTqoCKO+iggxgaGqKpqWnckX7ogDlz+NDvf8/FBC3w5dnb1mnYf7GLt1SlI1I/VAU0g4U1+KUa/z8GHvr97+kGbiFYXzdD8ca/u7s7VxXU1dU1bm3f5ubmcWMDhY2/unREGkSx04J6vcWtC6jcbpN58+YV7ZY5DPz2bHfPevATy+gici89V050Ld3CfapLR6R+UaILaHZ1041MRSaTyXWhTNSVsnLlSu69997cGUAT8AngEmA78GFgHTBW4vVh11FXV1fu6B2CCp/w4rHOzs5c9ZCZ5cpLV6xYoW4ekQalLqAGUljyGd76+vpyjf/xBCtwXQrcTNDd8wXAZo/P9V1dXbg7xx9/fG5KhlIllr29vXmPr1ixojofWkQqRgmggSWTyVw55j4EC6/8N8FUDqcTLMn4XHbbHTt2jHv9unXr6Ozs3K3++krW+YtIdagLqAEMDAyMm2ETdq3O1Ql8GUgB/0AwgVvhAG+xCqFyp14uJpVKTWnuHhGpPzoDaGB7A1cCdxL076eBj1C8uqdYhdCGDRsqGZ6I1DklgDpRahrnjRs3Fl2H99UE62heQHDU/yqCydsKhf384S1a1tnX16euG5E4K1YaVK+3mVoGOlG5JwVlmpZdfvEP4L8Ef0MZpZ3TXaKpq3pFGgslykB1BlBFpY7yM5kM8+bNy23X09OTey66oMoi4D8JFme5mWBFrh9WKfYoTcYmMjNoELiKJqrrL7U4S3t7O+3t7ST6+vga0AK8H/hSkW27u7tLLgAjIlJIZwB1InoGEGXunDs4yA+A3wAdFG/8QUfmIjI1OgOoEwcccMC4s4DLPvlJPnrXXZwFfB04n/EVPmaGu9PV1VWlSEVkplACqAOZTGZcpc+NPT18F3gZcIEZV5WYtdWzj69bt44NGzaoLl9EyqYuoDqQyWRIJpO5398E3AcsAE6Bko1/qLu7G3dX4y8iU6IzgCrbuHEjQ0NDrF27Nu8q3PBq3wuAzwOPAW8Bni54vZlx8sknA5prX0T2TNXOAMwsYWYPm9mt2d9fZmb3m9kmM/umme1VrVhqISwBDefuv+aaa/KeS5hxJcGVvbcAJzC+8QdYs2ZNFaIVkTio2opgZnYhQRHLvu5+ppl9C7jJ3W80s2uBR9z9moneo9FXBAtX7io0B/hXgoVb/h74OMEVVuVQqaeITKamK4KZ2VLgDLIVjBbMbHYKwZrkAF8F3lqNWGqpWOPfAtxB8OEvAD7GxI1/2N+vUk8R2VPV6gK6ErgY2Jn9/QBgs7uHcxQ/AxxU7IVmdr6Z9ZtZ//DwcOUjraDwqt7w5yKgF3gN8HaC6ZxD4Rw+hUsyiohMl4p3AZnZmcDp7t5lZp0EPRyrgPvc/bDsNgcDt7v7URO9VyN3AWUymbyrdA8G/osg6/0RwVlAKJlMMjAwUNX4RGTm2u0uIDN7g5m17sG+TwDOMrMB4EaCrp9/ABaYWViFtBT41R7so+5lMhnS6TQASeAugjOAU8lv/CGY51+LrIhIpZVTBnoH8IKZ7QQeJ6hQfDT7c4O7b5/oxe7+KeBTAOEZgLufa2b/BryNICm8B/jO7n6IRrKMoNtnX4JM+HDkuY6ODh544IFahCUiMVTOGMCHgCGCLurPAP8DrCQoVx8/UX35PgFcaGabCMYEvrwH71XXwhLQn/f10Usw8HsqQePf1tZGS0sL6XRajb+IVNWkCcDdryboxnGCwdxR4AJ3f727L57Kzty9193PzN7/ubsf6+6Hufs5k51JNIpiC7f39PSwbK+9+BGwH0Hj/1B2+6GhIa2tKyI1MaVBYDNrAS4imK3gg+5+f6UCK6ZRBoHDev+2tjba29s57bjjOO2zn6UdeAMQfmnq8hGRaig1CDzpGICZnQwckb0tJxi7/B1Bt40UEdb7Dw0NceQhh/DOG2+kDfi/7Gr8Ac4444xahCciApQ3CNwLrCcYrMgIr8wAABEnSURBVL3K3QcqGdBM0NbWxtDQEAmg6557WAb8CfCDgu3CslB1+4hILZSTANYARxFcyfsxM/sNQQXQY8Dj7n5LBeNrOKlUKncGcBnBFb4fAr5ZZFtN4yAitTRpAnD36yB3sdZ+wIvAKwmWpD2bYO4yIX9e/z8nuOLtC8A/Fmy3ZMmSotNCiIhUUzljACngJuBA4PfAEuBHBPX8n61kcI0qDVxLsID7R4s8r8ZfROpBOdcBfBa4zt0Pyk7d0AL8B/A9MzusotHVucKSz97eXpYB/w5sBN4BjBV5nco9RaQelJMADg+7gQDcfYe7X08wNhDrmcrC6R3S6TTuzhtOOolvEZxWnQX8b8H2yWRSs3iKSN0oJwEUvVDA3X9AUBYaaxs3bsxdxNXymc/wWuC9wM+yz4ezeqbTaVKpVO0CFREpUE4V0GIzex/BPEAb3H1L5LnqrCZTx8L+/LcBHwauAG6OPL9u3TpaW/dkLj0RkcooJwFkgBXAnwFHmdnvCJLB48CUpoKYiebNm8fCLVv4MvBjggmOQtFqn87OzhpEJyJSWjlzAV3v7h9y97S7HwCcCFwN/Bboq3SA9ayzs5OtW7ZwA8Gp0DsBn70rpz777LPq6xeRulXOGUAed3+GYAWv26c/nMaRSqUYHBzkowRln6vILuK+I1jkLJFIsCN7P5PJ0NcX5Eoz0wVgIlIXqrUk5IwRln4ODg6yHPgbgoUMvpp9vqmpCYDVq1fnvUbr+IpIvZnyGYAEZgP/TDAr3vmRx9va2kilUlx99dW1CUxEpEw6A5ii8Oj9Q0AHwcUQL2SfW7JkCYODg5rbX0QaghLAFHR2dmJmtAE9wK0EV/2Gnn32Wbq7u9XVIyINQV1AU9Db24uZ8XmgCbgg+7gWdhGRRqQzgCk6haDc82+Bn2cf6+/v19G+iDQcJYBJRCd828uMfySY5iE6DarKOkWkESkBTCKc8K2pqYnVBJMffRjYDrn+fjX+ItKIlAAmEB799/X10TQ6yiUECyF8r9aBiYhMg4oPApvZ3sBdQHN2f992924zewNwOUES2gKscvdNlY5nKjKZDL29vfT19fFhghVx/mz+fPx/Cyd6FhFpPNU4A9gOnOLuryKYVO40MzsOuAY4191XAF8HLqlCLLtlfzMuBr4L/OB3v1OXj4jMCBVPAB4Ip5Buyt48e9s3+3gLUJfrJN59991c6M5+BBlq9uzZSgAiMiNUZQzAzBJmtp7gotk73P1+4DyCZSWfAf4UuKzEa883s34z6x8eHq5GuHkW7tzJRwhOUR4DduzYoamdRWRGqEoCcPexbFfPUuBYMzuKYL300919KfBPwN+XeO317t7h7h21WFjloqYm5hAsigBB5U9vb2/V4xARmW5VrQJy983AncCbgVdlzwQAvgkcX81YSonW/c81472jo9xCsMg7oMZfRGaMalQBtQKj7r7ZzOYApxJcR9ViZoe7+0+zjz1Z6VjKEVb+AHx83jwOuO02rog8r+4fEZkpqjEX0BLgq2aWIDjj+Ja732pm7wf+3cx2Eqwu9udViGVS4eItBlwL9AP3ZJ/r6urSALCIzBgVTwDu/ihwTJHHbyZ//fS6chpwBHBu5LFwgXclARGZCTQbaEQmk6GnpwcIRqifAf4t+5zm+xGRmUZTQURkMhmampo4imBQ4h+B0exzGvwVkZlGCaDA6Ogo5wHbgOuzjyUSCSUAEZlxlAAiOjs7mQ28C/gPgpFpgEWLFtUuKBGRClECiOjt7eWNwCLgX7OPJZNJhobqcpYKEZE9ogQQ0dnZybuB3wC3Zx8bHBzUAu8iMiOpCigrlUrxm8FBvgd8lV2Dv8lkkoGBgdoFJiJSIToDyFq1ahVvBeayq/unu7tbjb+IzFhKAFmZTIZ3A78A7q11MCIiVaAEQND4LzHjjew6+u/o6FC/v4jMaEoAWecACeBr2d/7+/s18ZuIzGhKAARnAKcDTwBPZR9Lp9O6+EtEZrTYJ4BUKsXeZqSBH9Q6GBGRKop1AshkMgwODnIiMIf8BKDuHxGZ6czdax1D2To6Ory/v39a3zOVStE1OMgFwP7APq2tvPDCC9O6DxGRWjKzB929o/BxnQEMDvImgkVfXgKGh4d19C8isRD7M4AlZjwLfAL4XPaxRvpOREQmozOAAp2dnVi29h929f+bWa1CEhGpqtgmgN7eXsyMNwHPA49kH3d3XQAmIrEQ2wSQyWTAnVOBOwB1+ohI3MQyAYRr/x4NHEh++afW/hWRuIhlAgj9n+zPO7I/29ra1PiLSGzEOgG8Fvgp8Fz296GhIS3+IiKxUfEFYcxsb+AuoDm7v2+7e7cF5TafIZiHbQy4xt2vqnQ8nZ2d9PX1AbACiBaVqvtHROKkGmcA24FT3P1VBG3uaWZ2HLAKOBg4wt2XAzdWIZbcRV77AocC66uxUxGROlTxMwAPrqrakv21KXtzYA3wJ+6+M7tdVedfODr7M0wAuvhLROKmKmMAZpYws/XAC8Ad7n4/wQH4O8ys38xuN7P2Eq89P7tN//Dw8B7HEk7xvCL7e5gA1q5du8fvLSLSSKqSANx9zN1XAEuBY83sKIIxgW3Zy5O/CHylxGuvd/cOd+9obW3dozgymUxe///z7BoAvu666/bovUVEGk1Vq4DcfTNwJ3Aa8AxwU/apm9nVK1MVK8jv/1+9enU1dy8iUnMVTwBm1mpmC7L35wCnAv8D3AK8PrtZmqAis6LCCp/ZwFHkJ4Crr7660rsXEakr1TgDWALcaWaPAg8QjAHcClwGnG1mjwF/C5xXhVgAWE7Q/xRNAKr/F5G4qUYV0KPAMUUe3wycUen9R4UNfOEAcHNzM9u2batmKCIiNRfLK4FXECz+EvY5ffKTn6xhNCIitRGrBHDDDTcAQQJ4DNhZy2BERGosVglg8+bNQJAAHo483tPTo/5/EYmdWCUACOae2J/8AWDNASQicRS7BFA4ANzU1KTGX0RiKVYJYNu2bawg6Pt/LPvY6OioEoCIxFLFy0DrRSqVYvv27SwHfkFQBaTyTxGJs9icAaxatQqANoI5KAC2b9+uo38RiS1rpGmQOzo6vL+/f/INSzAzngIeAt6VfayRPr+IyO4wswezE2/mic0ZwIIFC4BgXornJt5URCQWYpMAVqxYwT7AfODZ7GPNzc01jEhEpLZikQDCdQCWZH8PE8D27dtzS0SKiMRNbBKAmY1LAMlkMrdCmIhI3MQmAbj7uAQQVgaJiMRRLBJAqDABXHbZZbUKRUSk5mKVABYD24EXs79rGmgRibNYJICwnz9aAppMJnURmIjEWiwSQGgJu7p/nnnmmYk2FRGZ8WKXAMIzgNmzYzMNkohIUbFIAGGtf/QMQP3/IhJ3sUgAV155JU3AQnYlABGRuItFAliwYAGLs/fDBKABYBGJu4onADPb28x+YmaPmNkGM+speP4qM9tS6TgKrwFIpVKV3qWISF2rxhnAduAUd38VwYqMp5nZcQBm1gHsV+kAVq1alZcAuru7GRgYqPRuRUTqWsUTgAfCI/ym7M3NLAFcDlxc6RhuuOGGvATQ09OjSeBEJPaqMgZgZgkzWw+8ANzh7vcDHwS+6+4Tjsua2flm1m9m/cPDw7u1/82bN7OEYC3gF3brHUREZp6qJAB3H3P3FcBS4FgzOxk4B/hCGa+93t073L2jtbV1t/a/YsUKlhA0/mNAOp3WLKAiEntVrQJy983AncDrgcOATWY2AMw1s02V2u/AwEDeRWDq/xcRqU4VUKuZLcjenwOcCjzo7ovdPeXuKeAldz+sknEsRtcAiIhEVeMMYAlwp5k9CjxAMAZwaxX2m5NKpfKuAlYJqIgIVHxCHHd/FDhmkm3mVTKG3h/+kLHZs2latgwfHKzkrkREGkYsrgS+/BOfIAHc//TTmJmuAhYRAczdax1D2To6Ory/v3/qL1y/Ho45Br79bTj77OkPTESkjpnZg+7eUfh4LM4AeDbb+79kycTbiYjEiBKAiEhMxSsBLF488XYiIjESnwTQ0gJz5tQ6EhGRuhGfBKDuHxGRPPFYGLejA9rbax2FiEhdiUcC+NSnah2BiEjdiUcXkIiIjKMEICISU0oAIiIxpQQgIhJTSgAiIjGlBCAiElNKACIiMaUEICISUw21HoCZDQNTXdJrIfDrCoQznRTj9Kj3GOs9PlCM06XeYky6e2vhgw2VAHaHmfUXWwihnijG6VHvMdZ7fKAYp0sjxAjqAhIRiS0lABGRmIpDAri+1gGUQTFOj3qPsd7jA8U4XRohxpk/BiAiIsXF4QxARESKUAIQEYmpGZMAzOw0M3vKzDaZ2SeLPN9sZt/MPn+/maXqMMYLzewJM3vUzH5oZsl6ii+y3dlm5mZW9TK3cmI0s7dnv8cNZvb1eovRzJaZ2Z1m9nD23/r0Ksf3FTN7wcweL/G8mdlV2fgfNbNXVzO+MmM8NxvbY2Z2r5m9qt5ijGz3GjPbYWZvq1ZsZXP3hr8BCeBnwCHAXsAjwCsKtukCrs3efyfwzTqM8fXA3Oz9NdWMsZz4stvNB+4C7gM66vA7bAceBvbL/r6oDmO8HliTvf8KYKDKMZ4MvBp4vMTzpwO3AwYcB9xfzfjKjPH4yL/xm+sxxsjfw4+A7wFvq3aMk91myhnAscAmd/+5u/8BuBF4S8E2bwG+mr3/beANZmb1FKO73+nuL2V/vQ9YWk/xZV0KfBbYVsXYQuXE+H7ganf/LYC7v1CHMTqwb/Z+CzBUxfhw97uAFyfY5C3AP3vgPmCBmS2pTnSByWJ093vDf2Oq/38ljGGy7xHgQ8C/A9X+OyzLTEkABwG/jPz+TPaxotu4+w5gBDigKtEV7D+rWIxR7yM4CquWSePLdgUc7O63VTGuqHK+w8OBw83sv83sPjM7rWrRBcqJMQO828yeITgy/FB1QivbVP9Wa63a/1fKYmYHAX8EXFPrWEqJx6LwDcbM3g10AOlaxxIys1nA3wOrahzKZGYTdAN1EhwV3mVmr3T3zTWNKt+7gBvc/fNm9jrgX8zsKHffWevAGo2ZvZ4gAZxY61iKuBL4hLvvrG5nQ/lmSgL4FXBw5Pel2ceKbfOMmc0mOPX+TXXCy9t/qFiMmNkbgb8A0u6+vUqxweTxzQeOAnqzf8yLge+a2Vnu3l8nMUJwtHq/u48CvzCznxIkhAeqE2JZMb4POA3A3X9sZnsTTB5WL90EZf2t1pqZHQ18CXizu1fz/3K5OoAbs/9fFgKnm9kOd7+ltmFF1HoQYjpuBIns58DL2DXwdmTBNmvJHwT+Vh3GeAzBAGJ7PX6HBdv3Uv1B4HK+w9OAr2bvLyToyjigzmK8HViVvb+cYAzAqvxdpig9wHoG+YPAP6n232MZMS4DNgHH1yK2cmIs2O4G6nAQeEacAbj7DjP7IPB9glH3r7j7BjP7K6Df3b8LfJngVHsTwcDNO+swxsuBecC/ZY8annb3s+oovpoqM8bvA//HzJ4AxoCLvIpHh2XG+DHgi2b2UYIB4VWebSWqwcy+QdBFtjA7DtENNGXjv5ZgXOJ0ggb2JeC91YptCjF+mmAMb132/8oOr/Lsm2XEWPc0FYSISEzNlCogERGZIiUAEZGYUgIQEYkpJQARkZhSAhARiSklABGRmFICEBGJKSUAkT1gZj8ys/XZ2zYze3utYxIply4EE5kGZraGYD2Hd7n7WK3jESnHjJgKQqSWzOzPCBYlOVuNvzQSJQCRPWBm5wDnAm/xYAZSkYahBCCym8zsTIKlRs9091qskCayRzQGILKbzOw3BDPLbs0+9AV3/3INQxKZEiUAEZGYUhmoiEhMKQGIiMSUEoCISEwpAYiIxJQSgIhITCkBiIjElBKAiEhM/X/1kGkmeAyxyAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "from astropy.cosmology import LambdaCDM\n", "\n", "# Function which determines the supernova luminosity distance for given cosmological parameters (Omega_m, Omega_L)\n", "def getdmmod(p,z):\n", " om,ol = p[0],p[1]\n", " cosmo = LambdaCDM(H0=70.,Om0=om,Ode0=ol)\n", " r = cosmo.comoving_distance(z).value\n", " dl = r*(1.+z)\n", " return 5.*np.log10(dl) + 25.\n", "\n", "# Read in and plot the supernova dataset\n", "# You will need to change the file path to the location where you've saved the data\n", "datfile = '../datasets/SN_dataset.dat'\n", "data = np.loadtxt(datfile)\n", "z,dm,dmerr = data[:,0],data[:,1],data[:,2]\n", "plt.plot(z,dm,marker='o',markersize=2,color='black',linestyle='None')\n", "plt.errorbar(z,dm,yerr=dmerr,color='black',capsize=2.,linestyle='None')\n", "# Overplot a standard cosmological model as an example\n", "xmin,xmax = 0.,1.5\n", "zmod = np.linspace(0.01,xmax,100)\n", "p = [0.3,0.7] # values of (Omega_m, Omega_L)\n", "dmmod = getdmmod(p,zmod)\n", "plt.plot(zmod,dmmod,color='red')\n", "plt.xlabel(r'$z$')\n", "plt.ylabel(r'$DM$')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now:\n", "\n", "- Minimize $\\chi^2$ and find the best-fitting $(\\Omega_m, \\Omega_\\Lambda)$\n", "- Construct 2D confidence regions in $(\\Omega_m, \\Omega_\\Lambda)$\n", "- Determine the individual errors in $\\Omega_m$ and $\\Omega_\\Lambda$" ] }, { "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 }