Practical_Statistics.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": [],
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/eKeiran/astro-codes/blob/main/Practical_Statistics.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Frequentist vs Bayesian: philosophical distinction\n",
        "\n",
        "### <ins>Frequentist statistical inference</ins>\n",
        "*   Probabilities of relative frequency of measurements\n",
        "*   Physical parameters are fixed, constant, and unknown\n",
        "\n",
        "### <ins>Bayesian statistical inference</ins>\n",
        "*   Probability distribution of physical parameters\n",
        "*   Degrees of uncertainty of physical parameters; can be updated!"
      ],
      "metadata": {
        "id": "YFpSU-6a5369"
      }
    },
    {
      "cell_type": "code",
      "source": [],
      "metadata": {
        "id": "STOvcIOctnw6"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "\n",
        "\n",
        "---\n",
        "\n"
      ],
      "metadata": {
        "id": "lcwXAbi_Ui9c"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Practical Bayesian inference example\n",
        "A problem with two state random variable"
      ],
      "metadata": {
        "id": "JzagudMp6QTT"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "### <ins>Problem Statement\n",
        "\n",
        "*   A disease or infection that is prevalent in 3% of the population\n",
        "*   A test for positive or negative identification of the disease<br>\n",
        "    – In sick people it tests positive 90% of the time: P(positive|sick) = 0.9<br>\n",
        "    – In healthy people it tests negative 85% of the time: P(negative|healthy) = 0.85\n",
        "\n",
        "Q. What is the probability that a person testing positive is infected?\n"
      ],
      "metadata": {
        "id": "DZlHFADlKVj_"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "### <ins>The probability matrix\n",
        "\n",
        "|            | **Sick**     | **Healthy** |\n",
        "| :---       |  :----:  |   ---:  |\n",
        "| **Positive**   | 0.9      | 0.15    |\n",
        "| **Negative**   | 0.1      | 0.85    |"
      ],
      "metadata": {
        "id": "LpPjOVoALBRm"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "### <ins>Frequentist vs Bayesian question<br>\n",
        "<mark>Frequentist</mark>: What is the probability of testing positive given one is healthy?<br>\n",
        "<mark>Bayesian</mark>: What is the probability that one is sick after having tested positive?\n"
      ],
      "metadata": {
        "id": "g0Pxu71WMkYk"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "### <ins>Computation for when tested positive</ins>\n",
        "\n",
        "Probability of infered condition given an observed state\n",
        "\n",
        "*   $$P(\\text{Sick}|\\text{positive}) = \\frac{0.9}{(0.9+0.15)}$$\n",
        "\n",
        "*   $$P(\\text{Healthy}|\\text{positive}) = \\frac{0.15}{(0.9+0.15)}$$\n",
        "\n",
        "Incorporating prior information<br>\n",
        "$$P(\\text{Sick}) = 0.03$$\n",
        "$$P(\\text{Healthy}) = 0.97$$\n",
        "$$\\implies$$\n",
        "\n",
        "*   $$P(\\text{Sick}|\\text{positive}) = \\frac{(0.9\\times0.03)}{(0.9\\times0.03+0.15\\times0.97)} = 0.16$$\n",
        "\n",
        "*   $$P(\\text{Healthy}|\\text{positive}) = \\frac{(0.15\\times0.97)}{(0.9\\times0.03+0.15\\times0.97)} = 0.84$$\n",
        "\n",
        "\n"
      ],
      "metadata": {
        "id": "cG-2rcdbNLzu"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "### <ins>General formula\n",
        "$$P(\\text{Sick}|\\text{positive}) = \\frac{P(\\text{positive}|\\text{Sick})\\cdot P(\\text{Sick})}{\\left[P(\\text{positive}|\\text{Sick})+P(\\text{positive}|\\text{Healthy})\\right]}$$<br><br>\n",
        "<ins>Bayes' theorem</ins>\n",
        "Posterior probability of model parameters $\\theta$ for observed data $D$,\n",
        "$$P(\\theta|D) = \\frac{P(D|\\theta)\\cdot P(\\theta)}{P(D)}$$\n",
        "where,\n",
        "$P(D|\\theta)$ is the likelihood,<br>\n",
        "$P(\\theta)$ is the prior probability of model parameters, and<br>\n",
        "$P(D)$ is the probability of the data (normalization)"
      ],
      "metadata": {
        "id": "SjBERudURfnM"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "\n",
        "\n",
        "---\n",
        "\n"
      ],
      "metadata": {
        "id": "B-rcpcrpUe1V"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Generalization to continuous variable"
      ],
      "metadata": {
        "id": "_9j6-KfQ6Zwv"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "<ins>Example:</ins> Measuring flux from a source or spectrum<br>\n",
        "Bayesian: What is the probability distribution of the source flux given a set of observations?\n"
      ],
      "metadata": {
        "id": "fZV649d6QS6F"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "### <ins>Probabilities matrix\n",
        "\n",
        "\n",
        "*   Health $\\Longleftrightarrow$ True flux<br>\n",
        "    Sick, Healthy $\\Longleftrightarrow$ Range of possible source flux values\n",
        "\n",
        "*   Test $\\Longleftrightarrow$ Observed flux<br>\n",
        "    Positive, Negative $\\Longleftrightarrow$ Range of possible observed flux values\n",
        "\n",
        "*   Probability distributions: probability of observing a flux value given a source flux<br>\n",
        "    $P(\\text{positive}|\\text{Sick})\\Longleftrightarrow P(\\text{f}_\\text{obs}|~\\text{F}_\\text{true})$\n",
        "\n",
        "*   $$P(\\text{Sick}~|~\\text{positive}) = \\frac{P(\\text{positive}~|~\\text{Sick})\\cdot P(\\text{Sick})}{P(\\text{positive}~|~\\text{Sick})+P(\\text{positive}~|~\\text{Healthy})}$$\n",
        "$\\Longleftrightarrow$\n",
        "    \n",
        "$$P(\\text{F}_\\text{true}~|~\\text{f}_\\text{obs}) = \\frac{P(\\text{f}_\\text{obs}~|~\\text{F}_\\text{true})\\cdot P(\\text{F}_\\text{true})}{\\left[\\sum_\\text{F}P(\\text{f}_\\text{obs}~|~\\text{F}_\\text{true})\\cdot P(\\text{F}_\\text{true})\\right]}$$\n",
        "\n",
        "What is the probable source flux given an observed flux?"
      ],
      "metadata": {
        "id": "Y4khrVbPQav6"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Python introduction tasks"
      ],
      "metadata": {
        "id": "FsS6e5r56luy"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "\n",
        "1. Colab introduction\n",
        "2. Import libraries\n",
        "3. Load and read data\n",
        "4. Data manipulation\n",
        "5. Plotting data and printing results"
      ],
      "metadata": {
        "id": "5R20pPrt9G5D"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "\n",
        "### Importing libaries\n",
        "<ins>Some basic libraries:</ins><br>\n",
        "Modules for simple math and random number generators\n",
        "```\n",
        "import math\n",
        "import random\n",
        "```\n",
        "<ins>Data handling and array computations libaries:</ins><br>\n",
        "Pandas for data storage and manupulation<br>\n",
        "Numpy for mathematical operations with large arrays\n",
        "```\n",
        "import pandas as pd\n",
        "import numpy as np\n",
        "```\n",
        "<ins>Scientific computing libraries:</ins><br>\n",
        "Scipy for general-purpose scientific computing\n",
        "Astropy for specialized astronomy-themed computations\n",
        "```\n",
        "import scipy as scp\n",
        "import astropy.constants as asc\n",
        "```\n",
        "<ins>Plotting libraries:</ins><br>\n",
        "Matplotlib for plotting\n",
        "```\n",
        "import matplotlib.pyplot as plt\n",
        "```"
      ],
      "metadata": {
        "id": "aI5-iwD6zCNI"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Simple statistical computation in Python"
      ],
      "metadata": {
        "id": "hEqzO4zC6vS-"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "<ins>Initializing task</ins>\n",
        "\n",
        "Import libaries for\n",
        "\n",
        "*   Array manipulation/computations\n",
        "*   Scientific computing\n",
        "*   Plotting\n",
        "\n",
        "<mark>User input required!!</mark>"
      ],
      "metadata": {
        "id": "zdTGaKfe3UwN"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# Library for random number generator\n",
        "import random\n",
        "# Library for array manipulation/computation\n",
        "import numpy as np\n",
        "# Library for scientific computing\n",
        "import scipy as scp\n",
        "# Library for plotting modules\n",
        "import matplotlib.pyplot as plt"
      ],
      "metadata": {
        "id": "Jd1oxmvA31Jm"
      },
      "execution_count": 1,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Flux measurements"
      ],
      "metadata": {
        "id": "hdALZiTy8m5e"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "<ins>Tasks</ins>\n",
        "\n",
        "1.   Generate flux data by deciding on the number of data points, a true flux value, and the sampling distribution\n",
        "2.   Evaluate the uncertainties of the flux values\n",
        "\n",
        "<ins>Algorithm</ins>\n",
        "\n",
        "1.   Set true flux and number of data points\n",
        "2.   Pick the statistic for generating data (e.g., [Poisson](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.poisson.html)) and generate flux values from the distribution (using scipy.stats modules)\n",
        "3.   Compute the uncertainties for the generated data by taking the [square root](https://numpy.org/doc/stable/reference/generated/numpy.sqrt.html) of the data values\n",
        "\n",
        "<mark>User input required!!</mark>"
      ],
      "metadata": {
        "id": "umWs2b_X3KPM"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# True flux (F_true) and number of measuremets (N)\n",
        "F_true = 1000\n",
        "N = 50"
      ],
      "metadata": {
        "id": "fj08Rw_s8vhd"
      },
      "execution_count": 2,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# scipy.stats.poisson.rvs(mu, size=1)\n",
        "flux = scp.stats.poisson.rvs(F_true, size=N)\n",
        "sigma = np.sqrt(flux)"
      ],
      "metadata": {
        "id": "-FEB9Ei85qnE"
      },
      "execution_count": 3,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Visualize data"
      ],
      "metadata": {
        "id": "QXV7t0Yb8rF4"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "<ins>Algorithm</ins>\n",
        "\n",
        "1.   Initialize a figure class that holds plot elements\n",
        "2.   Use a plotting method to plot observed flux values with error bars\n",
        "3.   Add axes titles\n",
        "\n",
        "<mark>User input required!!</mark>"
      ],
      "metadata": {
        "id": "cZUMYdES7RuD"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# Initialize a figure class that holds plot elements\n",
        "fig = plt.figure()\n",
        "ax = plt.axes()\n",
        "# Plot the flux values with error bars\n",
        "ax.errorbar( np.arange(N), flux, yerr=sigma , fmt='ok', ecolor='gray', alpha=0.5)\n",
        "## Add X-axis and Y-axis labels\n",
        "plt.xlabel(\"Observation number\")\n",
        "plt.ylabel(\"Flux (counts)\")\n",
        "txt = plt.title(\"Observed fluxes\")"
      ],
      "metadata": {
        "id": "8Ymnwq-38wAD",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 472
        },
        "outputId": "cfacffa8-311a-48df-f412-4bc288c225c8"
      },
      "execution_count": 4,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 640x480 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkQAAAHHCAYAAABeLEexAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABN80lEQVR4nO3deVhUZf8/8PcMyoAIg4iyKJsbDD5qoEW4pySauWU/11zKNEtzycp8ylyyNC0r1LSeSnNJs1LKFo1cIA1JDdMEF9wYleUhZFUUmfv3R1/mcQRkRmY5M+f9uq65LufcZ87c53Cc+czn3hRCCAEiIiIiGVPaugJEREREtsaAiIiIiGSPARERERHJHgMiIiIikj0GRERERCR7DIiIiIhI9hgQERERkewxICIiIiLZY0BEREREsseAiIhMMn/+fCgUCuTl5dm6KlanUCgwf/78WvfLycnB448/jsaNG0OhUOD999/Hvn37oFAosG/fPovXk4hMx4CIiHDixAk88cQTaNasGVQqFfz9/TF69GicOHHC1lWzSzNnzsSuXbswZ84cbNiwAX379rV1lYioFvVsXQEisq1t27Zh5MiR8PLywoQJExASEoILFy7g008/xddff40tW7ZgyJAhtq6mXdmzZw8GDRqEF198Ub8tOzvbhjUiotowICKSsbNnz2LMmDFo0aIFkpKS0KRJE33Z9OnT0a1bN4wZMwbHjh1DixYtbFjTqnQ6HW7evAkXFxdbV6WK3NxceHp62roaRGQCNpkRydiyZctw7do1fPzxxwbBEAB4e3vjo48+QmlpKZYuXVrltXl5eRg2bBg8PDzQuHFjTJ8+HWVlZQb7JCQkoGvXrvD09ETDhg0RGhqKf//73wb73LhxA/PmzUOrVq2gUqkQEBCAl19+GTdu3DDYT6FQYOrUqdi0aRPatm0LlUqFHTt2wMvLC08++WSV+hUVFcHFxcUgS2Pse924cQMzZ85EkyZN4O7ujoEDB+LSpUu1Xs9169ZBoVBACIFVq1ZBoVBAoVDUuH9wcDDGjx9fZXvPnj3Rs2dP/fNx48bBxcUF6enpBvvFxsaiUaNGuHLlin7bTz/9hG7dusHNzQ3u7u7o379/labP7OxsPPnkk2jevDlUKhX8/PwwaNAgXLhwodZzJHJUzBARydiOHTsQHByMbt26VVvevXt3BAcH44cffqhSNmzYMAQHB2Px4sU4ePAg4uLicPXqVaxfvx7AP/2SHn30UbRv3x4LFy6ESqVCRkYGDhw4oD+GTqfDwIEDsX//fkyaNAkajQbHjx/He++9h9OnTyM+Pt7gPffs2YOtW7di6tSp8Pb2RuvWrTFkyBBs27YNH330EZydnfX7xsfH48aNGxgxYoTJ7/X0009j48aNGDVqFDp37ow9e/agf//+tV7P7t27Y8OGDRgzZgwefvhhjB07ttbXGOODDz7Anj17MG7cOCQnJ8PJyQkfffQRfv75Z2zYsAH+/v4AgA0bNmDcuHGIjY3F22+/jWvXrmH16tXo2rUrUlNTERwcDAAYOnQoTpw4geeffx7BwcHIzc1FQkICMjMz9fsQyY4gIlkqKCgQAMSgQYPuut/AgQMFAFFUVCSEEGLevHkCgBg4cKDBfs8995wAIP78808hhBDvvfeeACD++9//1njsDRs2CKVSKX799VeD7WvWrBEAxIEDB/TbAAilUilOnDhhsO+uXbsEALFjxw6D7Y888oho0aKFye919OhRAUA899xzBvuNGjVKABDz5s2r8Xxur+uUKVMMtu3du1cAEHv37tVvCwoKEuPGjavy+h49eogePXpUe56LFi0S586dEw0bNhSDBw/WlxcXFwtPT08xceJEg9dlZ2cLtVqt33716lUBQCxbtqzW8yCSEzaZEclUcXExAMDd3f2u+1WWFxUVGWyfMmWKwfPnn38eAPDjjz8CgL4PzbfffgudTlftsb/66itoNBqEhYUhLy9P/+jVqxcAYO/evQb79+jRA+Hh4QbbevXqBW9vb3z55Zf6bVevXkVCQgKGDx9u8ntV1n/atGkG7zNjxoxqz8Fa+vTpg2eeeQYLFy7EY489BhcXF3z00Uf68oSEBBQUFGDkyJEG5+fk5ISoqCj9+bm6usLZ2Rn79u3D1atXbXU6RJLDJjMimaoMdCoDo5rUFDi1bt3a4HnLli2hVCr1/VCGDx+OTz75BE8//TReeeUV9O7dG4899hgef/xxKJX//BY7c+YM0tPTq/RfqpSbm2vwPCQkpMo+9erVw9ChQ/HFF1/gxo0bUKlU2LZtG8rLyw0CImPf6+LFi1AqlWjZsqVBeWhoaLWvs6Z33nkH3377LY4ePYovvvgCTZs21ZedOXMGAPQB3p08PDwAACqVCm+//TZmzZoFHx8fPPjgg3j00UcxduxY+Pr6Wv4kiCSKARGRTKnVavj5+eHYsWN33e/YsWNo1qyZ/gu1Jnd2HnZ1dUVSUhL27t2LH374ATt37sSXX36JXr164eeff4aTkxN0Oh3atWuH5cuXV3vMgICAKseszogRI/DRRx/hp59+wuDBg7F161aEhYWhQ4cO+n1MfS9rqKnDdUVFBZycnKpsT01N1Qdux48fx8iRI/VllVm4DRs2VBvY1Kv3v4/7GTNmYMCAAYiPj8euXbswd+5cLF68GHv27EFERESdzonIXjEgIpKxRx99FP/5z3+wf/9+dO3atUr5r7/+igsXLuCZZ56pUnbmzBmDjE1GRgZ0Op1Bp1ylUonevXujd+/eWL58Od566y28+uqr2Lt3L2JiYtCyZUv8+eef6N27911HY9Wme/fu8PPzw5dffomuXbtiz549ePXVVw32Mfa9goKCoNPpcPbsWYOs0KlTp+65fjVp1KgRCgoKqmy/ePFilWkOSktL8eSTTyI8PBydO3fG0qVLMWTIENx///0AoM9oNW3aFDExMbW+d8uWLTFr1izMmjULZ86cwX333Yd3330XGzdurPuJEdkh9iEikrGXXnoJrq6ueOaZZ/D3338blOXn52Py5Mlo0KABXnrppSqvXbVqlcHzFStWAAD69eunf/2d7rvvPgDQD3MfNmwYLl++jP/85z9V9r1+/TpKS0uNOg+lUonHH38cO3bswIYNG3Dr1i2D5jJT3quy/nFxcQb7vP/++0bVxRQtW7bEwYMHcfPmTf2277//Hlqttsq+s2fPRmZmJj7//HMsX74cwcHBGDdunP5axsbGwsPDA2+99RbKy8urvP6///0vAODatWtVpkdo2bIl3N3dq0w/QCQnzBARyVjr1q3x+eefY/To0WjXrl2Vmarz8vKwefPmKv1pAOD8+fMYOHAg+vbti+TkZP0w9cpmqoULFyIpKQn9+/dHUFAQcnNz8eGHH6J58+b6bNSYMWOwdetWTJ48GXv37kWXLl1QUVGBkydPYuvWrdi1axc6depk1LkMHz4cK1aswLx589CuXTtoNBqDcmPf67777sPIkSPx4YcforCwEJ07d8bu3buRkZFRx6td1dNPP42vv/4affv2xbBhw3D27Fls3LixyvXes2cPPvzwQ8ybNw+RkZEAgLVr16Jnz56YO3culi5dCg8PD6xevRpjxoxBZGQkRowYgSZNmiAzMxM//PADunTpgpUrV+L06dPo3bs3hg0bhvDwcNSrVw/bt29HTk6OfooCIlmy9TA3IrK9Y8eOiZEjRwo/Pz9Rv3594evrK0aOHCmOHz9eZd/KYfdpaWni8ccfF+7u7qJRo0Zi6tSp4vr16/r9du/eLQYNGiT8/f2Fs7Oz8Pf3FyNHjhSnT582ON7NmzfF22+/Ldq2bStUKpVo1KiR6Nixo1iwYIEoLCzU74dqhrLfTqfTiYCAAP3Q9OoY+17Xr18X06ZNE40bNxZubm5iwIABQqvVmn3YvRBCvPvuu6JZs2ZCpVKJLl26iMOHDxsMuy8qKhJBQUEiMjJSlJeXG7x25syZQqlUiuTkZIP3iY2NFWq1Wri4uIiWLVuK8ePHi8OHDwshhMjLyxNTpkwRYWFhws3NTajVahEVFSW2bt1a63kROTKFEELYMiAjIiIisjX2ISIiIiLZY0BEREREsseAiIiIiGSPARERERHJHgMiIiIikj0GRERERCR7nJjRSDqdDleuXIG7u3udlhggIiIi6xFCoLi4GP7+/vqFpavDgMhIV65cscnij0RERFR3Wq0WzZs3r7GcAZGR3N3dAfxzQWtb9ZuIiIikoaioCAEBAfrv8ZowIDJSZTOZh4cHAyIiIiI7U1t3F3aqJiIiItljQERERESyx4CIiIiIZI8BEREREckeAyIiIiKSPQZEREREJHsMiIiIiEj2GBARERGR7DEgIiIiItljQERERESyx4CIiIiIZI8BEREREckeAyIiIiKSPQZEREREJHv1bF0BIiIiR1RRUYFff/0VANCtWzc4OTnZuEZ0N8wQERERkewxICIiIiLZY0BEREREsseAiIiIiGSPARERERHJHgMiIiIikj0GRERERCR7DIiIiIhI9hgQERERkewxICIiIiLZY0BEREREsseAiIiIiGSPARERERHJHgMiIiIikj0GRERERCR7DIiIiIhI9hgQERERkewxICIiIiLZY0BEREREsseAiIiIiGSPARERERHJHgMiIiIikj0GRERERCR7DIiIiIhI9hgQERERkewxICIiIiLZY0BEREREsmfTgCgpKQkDBgyAv78/FAoF4uPjDcq3bduGPn36oHHjxlAoFDh69GiVY5SVlWHKlClo3LgxGjZsiKFDhyInJ8dgn8zMTPTv3x8NGjRA06ZN8dJLL+HWrVsWPDMiIiKyJzYNiEpLS9GhQwesWrWqxvKuXbvi7bffrvEYM2fOxI4dO/DVV18hMTERV65cwWOPPaYvr6ioQP/+/XHz5k389ttv+Pzzz7Fu3Tq8/vrrZj8fIiIisk/1bPnm/fr1Q79+/WosHzNmDADgwoUL1ZYXFhbi008/xRdffIFevXoBANauXQuNRoODBw/iwQcfxM8//4y0tDT88ssv8PHxwX333Yc33ngDs2fPxvz58+Hs7Gz28yIiIiL7Ytd9iI4cOYLy8nLExMTot4WFhSEwMBDJyckAgOTkZLRr1w4+Pj76fWJjY1FUVIQTJ05Yvc5EREQkPTbNENVVdnY2nJ2d4enpabDdx8cH2dnZ+n1uD4YqyyvLanLjxg3cuHFD/7yoqMhMtSYiIiKpsesMkSUtXrwYarVa/wgICLB1lYiIiMhC7Dog8vX1xc2bN1FQUGCwPScnB76+vvp97hx1Vvm8cp/qzJkzB4WFhfqHVqs1b+WJiIhIMuw6IOrYsSPq16+P3bt367edOnUKmZmZiI6OBgBER0fj+PHjyM3N1e+TkJAADw8PhIeH13hslUoFDw8PgwcRERE5Jpv2ISopKUFGRob++fnz53H06FF4eXkhMDAQ+fn5yMzMxJUrVwD8E+wA/2R2fH19oVarMWHCBLzwwgvw8vKCh4cHnn/+eURHR+PBBx8EAPTp0wfh4eEYM2YMli5diuzsbLz22muYMmUKVCqV9U+aiIiIJMemGaLDhw8jIiICERERAIAXXngBERER+jmCvvvuO0RERKB///4AgBEjRiAiIgJr1qzRH+O9997Do48+iqFDh6J79+7w9fXFtm3b9OVOTk74/vvv4eTkhOjoaDzxxBMYO3YsFi5caMUzJSIiIilTCCGErSthD4qKiqBWq1FYWMjmMyIiqlVFRQV+/fVXAEC3bt3g5ORk4xrJk7Hf33bdh4iIiIjIHBgQERERkewxICIiIiLZY0BEREREsseAiIiIiGTPrtcyIyIiqsRRXVQXzBARERGR7DEgIiIisgCdTofs7GycO3cOFy5cgE6ns3WV6C7YZEZERGRm6enp+Oabb5CUlISbN29i79690Gg0GDJkCDQaja2rR9VgQERERGRG6enpiIuLQ25uLtRqNVxcXNC4cWOkpqZCq9Vi2rRpDIokiE1mREREZqLT6bB9+3bk5eUhPDwcbm5ucHJygoeHB8LDw5GXl4f4+Hg2n0kQAyIiIiIzyczMxMmTJxEQEACFQmFQplAo0Lx5c6SnpyMzM9NGNaSaMCAiIiIyk+LiYpSVlcHNza3acjc3N5SVlaG4uNjKNaPaMCAiIiIyE3d3d7i4uKC0tLTa8tLSUri4uMDd3d3KNaPaMCAiIiIyk8DAQISFhUGr1UIIYVAmhMClS5eg0WgQGBhooxpSTTjKjIiIyEyUSiWGDBkCrVaLtLQ0KJVKuLi4oLCwEFeuXIG3tzcGDx4MpZL5CKnhX4SIiMiMNBoNpk2bhoiICBQWFkKr1SI/Px+RkZEcci9hzBARERGZmUajwezZs9GqVStcu3YNPXr0QEhICDNDEsaASOK4WCERkX1SKpXw9fUFAAQHBzMYkjgGREREMsEfWEQ1Y7hKREREsseAiIiIiGSPARERERHJHgMiIiIikj12qiYiIpIJdqyvGTNEREREJHsMiIiIiEj22GTmAJgCJSIiqhtmiIiIiEj2GBARERGR7DEgIiIiItljQERERESyx4CIiIiIZI8BEREREckeh92TAQ7ht0/8uxER1Q0zRERERCR7DIiIiIhI9thkRhbBJhwiIrInzBARERGR7DEgIiIiItljQERERESyx4CIiIiIZI8BEREREckeAyIiIiKSPQZEREREJHsMiIiIiEj2GBARERGR7DEgIiIiItljQERERESyx4CIiIiIZI+LuxIRkWxw4WmqCTNEJHkVFRXYt28f9u3bh4qKCltXh4gkSqfTITs7G+fOncOFCxeg0+lsXSWyI8wQERGR3UtPT8c333yDpKQk3Lx5E3v37oVGo8GQIUOg0WhsXT2yAwyIiIjIrqWnpyMuLg65ublQq9VwcXFB48aNkZqaCq1Wi2nTpjEoolqxyYyIiOyWTqfD9u3bkZeXh/DwcLi5ucHJyQkeHh4IDw9HXl4e4uPj2XxGtWJAREREdiszMxMnT55EQEAAFAqFQZlCoUDz5s2Rnp6OzMxMG9WQ7AUDIiIisghrDIgoLi5GWVkZ3Nzcqi13c3NDWVkZiouLLfL+5DgYEBERkd1yd3eHi4sLSktLqy0vLS2Fi4sL3N3drVwzsjcMiIiIyG4FBgYiLCwMWq0WQgiDMiEELl26BI1Gg8DAQBvVkOwFR5kREZHdUiqVGDJkCLRaLdLS0qBUKuHi4oLCwkJcuXIF3t7eGDx4MJRK/v6nu+MdQkREdk2j0WDatGmIiIhAYWEhtFot8vPzERkZySH3ZDRmiIiIyO5pNBrMnj0brVq1wrVr19CjRw+EhIQwM0RGY0BEREQOQalUwtfXFwAQHBzMYMhOSGV9Od4tEse1eYiIiCyPGSIJ49o8RGQsqfzKpv9xcnJCz549bV0Nm7Kn+5IBkURxbR4iMrfKjPO1a9cQEBDAPjZEt2FAJEF3rs2j1WoBAB4eHlCr1UhLS0N8fDxCQ0P5YUZERmHGmejubPptmpSUhAEDBsDf3x8KhQLx8fEG5UIIvP766/Dz84OrqytiYmJw5swZg32Cg4OhUCgMHkuWLDHY59ixY+jWrRtcXFwQEBCApUuXWvrU6kQua/NYY1p/Ivpfxjk1NRVqtRoBAQH6jHNcXBzS09NtXUUim7NpQFRaWooOHTpg1apV1ZYvXboUcXFxWLNmDVJSUuDm5obY2FiUlZUZ7Ldw4UJkZWXpH88//7y+rKioCH369EFQUBCOHDmCZcuWYf78+fj4448tem51wbV5iMhcuBo8kXFs2mTWr18/9OvXr9oyIQTef/99vPbaaxg0aBAAYP369fDx8UF8fDxGjBih39fd3V0/1PJOmzZtws2bN/HZZ5/B2dkZbdu2xdGjR7F8+XJMmjTJ/CdlBrevzdOwYcMq5Vybh4iMZUrGOTg42DaVJJIAyXZAOX/+PLKzsxETE6PfplarERUVheTkZIN9lyxZgsaNGyMiIgLLli3DrVu39GXJycno3r07nJ2d9dtiY2Nx6tQpXL161fIncg+4Ng8RmQszznQ7TuVSM8l2qs7OzgYA+Pj4GGz38fHRlwHAtGnTEBkZCS8vL/z222+YM2cOsrKysHz5cv1xQkJCqhyjsqxRo0bVvv+NGzdw48YN/fOioqK6n5SRuDYPEZkLM85UiR3r706yAZGxXnjhBf2/27dvD2dnZzzzzDNYvHgxVCrVPR938eLFWLBggTmqeE8q1+apvHn/+9//AgAiIyMxePBg3rxEZJTKjHNqairCwsIMyiozzpGRkcw4OzhO5VI7yQZElX2CcnJy4Ofnp9+ek5OD++67r8bXRUVF4datW7hw4QJCQ0Ph6+uLnJwcg30qn9fU7wgA5syZYxBsFRUVISAg4F5O5Z4ZuzYP5xYhKbKnCdkcGTPOxKlcjCPZgCgkJAS+vr7YvXu3PgAqKipCSkoKnn322Rpfd/ToUSiVSjRt2hQAEB0djVdffRXl5eWoX78+ACAhIQGhoaE1NpcBgEqlqlOGyVxqW5uHKVAiqg0zzvLGjvXGsWlAVFJSgoyMDP3z8+fP4+jRo/Dy8kJgYCBmzJiBRYsWoXXr1ggJCcHcuXPh7++PwYMHA/inw3RKSgoeeughuLu7Izk5GTNnzsQTTzyhD3ZGjRqFBQsWYMKECZg9ezb++usvfPDBB3jvvfdsccpmxRQoERmLq8Ebz9Gym8Z0rL98+bLsO9bbNCA6fPgwHnroIf3zyiaqcePGYd26dXj55ZdRWlqKSZMmoaCgAF27dsXOnTvh4uIC4J8szpYtWzB//nzcuHEDISEhmDlzpkFTl1qtxs8//4wpU6agY8eO8Pb2xuuvvy7ZIffGYgqUiEzF1eDliR3rjWPTgKhnz55VhpXfTqFQYOHChVi4cGG15ZGRkTh48GCt79O+fXt9tO8omAIlIiJjsGO9cfjzwE5xbhEiIjJGZcd6b29vpKWlobS0FBUVFSgsLERaWho71v8feZ+9Hbs9BVqde02BctIu+8S/GxHdTWXH+oiICBQWFkKr1SI/Px+RkZHsb/p/JDvKjO7OEilQjlizT/y7EZEx2LH+7hgQ2Slzzy3CEWv2iX83IjIFO9bXjFfCjpkrBWqJ1bDZhGN5XMWciMh8mCGyc+ZIgZp7xBqbcKyDIw1tx9HmqSEiBkQOoa4pUHNO2sUmHOvhZGuGGKQQUV3cU5NZeXk5tFotTp06hfz8fHPXiazMXCPW2IRjXZYaaUhEJEdGB0TFxcVYvXo1evToAQ8PDwQHB0Oj0aBJkyYICgrCxIkTcejQIUvWlSykcsSaVqutMlFm5Yg1jUZT64g1U5pwqO7M9XcjIiIjA6Lly5cjODgYa9euRUxMDOLj43H06FGcPn0aycnJmDdvHm7duoU+ffqgb9++OHPmjKXrTWZkrkm7OFmkdXGyNSIi8zGqD9GhQ4eQlJSEtm3bVlv+wAMP4KmnnsKaNWuwdu1a/Prrr2jdurVZK0qWZY7VsLlejvVxFXMiIvMwKiDavHmzUQdTqVSYPHlynSpEtlPXEWtcL8c2ONka2TN2hiepqPMos6KiIuzZswehoaH8NeoA6jJizdyTRZLxONkaEVHdmPypOWzYMKxcuRIAcP36dXTq1AnDhg1D+/bt8c0335i9gmRfuF4OWUJFRQX27duHffv2oaKiwtbVISIHZHKGKCkpCa+++ioAYPv27RBCoKCgAJ9//jkWLVqEoUOHmr2SZF/YhEPk+NjURY7G5G+owsJCeHl5AQB27tyJoUOHokGDBujfvz9Hl5FeZRNOixYt2IRDRESSZ/K3VEBAAJKTk1FaWoqdO3eiT58+AICrV6/CxcXF7BUkIiIisjSTm8xmzJiB0aNHo2HDhggKCkLPnj0B/NOU1q5dO3PXj4jIbrFZSXoqF56+du0aAgIC2JxPeiYHRM899xyioqKQmZmJhx9+WH8jtWjRAm+++abZK0jEDzAiMgcuPE13Y3JAtHDhQrz44ovo2LGjwfZevXph2bJl6Ny5s9kqR8QPMHJ0zCJZBxeeptqY/DN7wYIFKCkpqbL92rVrWLBggVkqRQT87wMsNTUVarUaAQEB+g+wuLg4pKen27qKRGQHuPA0GcPkgEgIUWXhTgD4888/9aPPiOqKH2BkCZzPSJ648DQZw+gms0aNGkGhUEChUKBNmzYGN1VFRQVKSkq4bAeZjSkfYMHBwbapJNUZm4vIGoxZePry5ctceFrmjA6I3n//fQgh8NRTT2HBggVQq9X6MmdnZwQHByM6OtoilST54QcYkf2TyoAILjxNxjA6IBo3bhwAICQkBJ07d0b9+vUtVikifoAR2TcpDYjgwtNkDJNHmfXo0QM6nQ6nT59Gbm5ulT4c3bt3N1vlSL74AUZkv6Q2oosLT5MxTA6IDh48iFGjRuHixYsQQhiUKRQKdlQks+AHGJF9unNAhFarBQB4eHhArVYjLS0N8fHxCA0Nter/38qFpyuzVv/9738BAJGRkRg8eDCH3JPpAdHkyZPRqVMn/PDDD/Dz86t2xBmROfADjMg26tLZXcoDIrjwNN2NyQHRmTNn8PXXX6NVq1aWqA+RAX6AEdkXqQ+IqFx4GgAXniYDJgdEUVFRyMjIYEBEVsMPMCL7YcsBEU5OTvr1Ncl+SGU0oskB0fPPP49Zs2YhOzsb7dq1qzLarH379marHBER2RcOiCBTSGk0oskB0dChQwEATz31lH6bQqHQz2DNTtVERNYnlV/ZHBBBt7vbfSm10YgmB0Tnz5+3RD2IHA5nYSZrkdKvbIADIugfd7svQ0NDJTca0eSAKCgoyBL1ICKieyC1X9mVOCBC3mq7L4cOHSq50YgmB0Tr16+/a/nYsWPvuTIkL1JJ8RPZK6nO+VOJAyLkyZj7cseOHbh+/bqkRiOaHBBNnz7d4Hl5eTmuXbsGZ2dnNGjQgAERGUVqKX4ieyTlOX9Ivoy5Ly9evAgAklqeyeRw/erVqwaPkpISnDp1Cl27dsXmzZstUUdyMJWp1NTUVKjVagQEBOhTqXFxcUhPT7d1FYnsgjFz/pSVlXERZLIqY+5LJycnBAQEQKvVVln1onI0okajsepoRLPkL1u3bo0lS5ZUyR4R3enOVGrlfwwPDw+Eh4cjLy8P8fHxVdbII6Kqbp/zpzpcBJlswZj70tXVFQMGDIC3tzfS0tJQWlqKiooKFBYWIi0tzSajEc32TvXq1cOVK1fMdThyUKak+Ino7irn/JHSr2yynYqKCuzbtw/79u2z6RQ4xt6XvXr1wrRp0xAREYHCwkJotVrk5+cjMjLSJoMBTO5D9N133xk8F0IgKysLK1euRJcuXcxWMXJMUp/Wn8iecM4fkiJT7kspjUY0OSAaPHiwwXOFQoEmTZqgV69eePfdd81VL3JQtpzWn6yLowitg3P+kBSZcl9KZTSiyQER+3ZYl6OtzWPLaf05UaL1cBShdUnpVzZRJXu7L00OiG5X2TZ4Z18Qopowxe/4pDpRoKOTyq9sR8DspvnY0315TwHR+vXrsWzZMpw5cwYA0KZNG7z00ksYM2aMWStH0lTXTAtT/I5L6hMFyp2jZZwtgdlN+TI5IFq+fDnmzp2LqVOn6jtR79+/H5MnT0ZeXh5mzpxp9kqS47G3VCoZhxMFkj1jdlPeTA6IVqxYgdWrVxvMSD1w4EC0bdsW8+fPZ0BERrOnVCoZh6MIyV4xu0km/1WzsrLQuXPnKts7d+6MrKwss1SKiOwTJwoke8U50sjkgKhVq1bYunVrle1ffvklWrdubZZKEZF94kSBZK/ksgxKZT+ynj17cqTtHUxuMluwYAGGDx+OpKQkfR+iAwcOYPfu3dUGSkQkHxxFSPbK1DnSOI2H4zH5U2no0KFISUmBt7c34uPjER8fD29vb/z+++8YMmSIJepIRHakchShVKbjt4TKYdnnzp3DhQsXOD+bA7BldlMqS27I3T0Nu+/YsSM2btxo7rrQPeJQWpIaRx5FyGHZjonZTTI5IPrxxx/h5OSE2NhYg+27du2CTqdDv379zFY5IrJf1h5FaI3J9Dgs27FxjjR5MzkgeuWVV7BkyZIq24UQeOWVVxgQEZHVWSNrw2HZ8uDI2U26O5MDojNnziA8PLzK9rCwMGRkZJilUkRExrJW1oaTTsoH50iTJ5P/ymq1GufOnauyPSMjo8bhikRElnBn1sbNzQ1OTk7w8PBAeHg48vLyEB8fb5ZOz3IZlk0kVyYHRIMGDcKMGTNw9uxZ/baMjAzMmjULAwcONGvliIjuxpqT6XHSSSLHZnJAtHTpUri5uSEsLAwhISEICQmBRqNB48aN8c4771iijmRFnLTLPsn172bNrA0nnbQMTmFAUmFyHyK1Wo3ffvsNCQkJ+PPPP+Hq6or27duje/fulqgfkVVxsjVpqmkEmamT6dWFLYdlO+p9ySkMSEruaR4ihUKBPn36oE+fPuauDxGRgbt9aYaGhiIsLAypqakICwszeF1l1iYyMtJsWRsOyzYfTmFAUmNUQLRlyxaMGDHCqANqtVpkZmbql/UgeeJkkWQOxnxpWjtrw2HZ/6jLvE+cwoCkyKg7bfXq1dBoNFi6dCnS09OrlBcWFuLHH3/EqFGjEBkZib///tvsFSUieTF2BFloaKjVlwqpHJbdokULWQ7LTk9Px9tvv43PPvsM69evx/z587FkyZJqvx+qw5XlHYOjLTliVIYoMTER3333HVasWIE5c+bAzc0NPj4+cHFxwdWrV5GdnQ1vb2+MHz8ef/31F3x8fCxdbyJycKZ8aTJrYz3maOoypjP85cuXOYUBWZXRfYgGDhyIgQMHIi8vD/v378fFixdx/fp1eHt7IyIiAhEREfzwISKzMfVLU4qT6VljORFrMldTlzU7wxMZy+RO1ZVt8kRElmTvX5rGjqCyp6DJXLN1V05hYK3O8ETGuKdRZkRElmbPX5rGNivZ27BzczV1cWV5kiIGREQmctQ5YaTGEl+a1sjGGNuspNPpsHLlSrsadm7OrB2nMJAuuX7GMSAiIpupLUAx55emtbIxxjQrpaWlIT8/3+6GnZs7a8fO8CQlNr3rkpKSMGDAAPj7+0OhUCA+Pt6gXAiB119/HX5+fnB1dUVMTAzOnDljsE9+fj5Gjx4NDw8PeHp6YsKECSgpKTHY59ixY+jWrRtcXFwQEBCApUuXWvrUjOJoQxaJTGHs0O3KL82nnnoKY8eOxfz58zF79myTg6G4uDikpqZCrVYjICBAn42Ji4szeri4MYxpVsrPz7fLYeeVWTtvb2+kpaWhtLQUFRUVKCwsRFpa2j1l7eQ+hQFJh8l3XllZWY1lWVlZJh2rtLQUHTp0wKpVq6otX7p0KeLi4rBmzRqkpKTAzc0NsbGxBnUYPXo0Tpw4gYSEBHz//fdISkrCpEmT9OVFRUXo06cPgoKCcOTIESxbtgzz58/Hxx9/bFJdiSxBrkGxqQFKXb40jZ3PyFxraBmzCKxSqYROp7PKGmzmVpm1s+a8T0TWYHKTWWRkJL744gvcd999Btu/+eYbTJ48WZ/SNka/fv3Qr1+/asuEEHj//ffx2muvYdCgQQCA9evXw8fHB/Hx8RgxYgTS09Oxc+dOHDp0CJ06dQIArFixAo888gjeeecd+Pv7Y9OmTbh58yY+++wzODs7o23btjh69CiWL19uEDgRkXVYe5Zic42MMpYxzUoajQZZWVl2O4KOTV3kiEy+e3v27IkHH3wQb7/9NoB//uOOHz8eY8aMwb///W+zVez8+fPIzs5GTEyMfptarUZUVBSSk5MBAMnJyfD09NQHQwAQExMDpVKJlJQU/T7du3eHs7Ozfp/Y2FicOnUKV69erfH9b9y4gaKiIoMHEdWdtWcpNqYJy5zZGGOalcaNGweNRgOtVgshhMHrbw+apDiCrhKbusjRmJwh+vDDD9G/f388/fTT+P7775GVlYWGDRvi999/x7/+9S+zVSw7OxsAqsx67ePjoy/Lzs5G06ZNDcrr1asHLy8vg31CQkKqHKOyrFGjRtW+/+LFi7FgwYK6nwgRGbD2LMW2mM/ImM7gSqVSssPO7WluJCJzuadRZv369cNjjz2G1atXo169etixY4dZgyEpmDNnDl544QX986KiIgQEBNiwRkSOwdoBiq3mM6qtWUmqw87tbW4kU3HhaaqJyQHR2bNnMWrUKGRnZ2PXrl1ITEzEwIEDMX36dLz55puoX7++WSpWOQV/Tk4O/Pz89NtzcnL0/Zd8fX2Rm5tr8Lpbt24hPz9f/3pfX1/k5OQY7FP5vHKf6qhUKqhUqjqfB9UdP8Aci7UDFFtOAljbciJS64tjjnXKiOyVyf/r7rvvPoSEhODPP//Eww8/jEWLFmHv3r3Ytm0bHnjgAbNVLCQkBL6+vti9e7d+W1FREVJSUhAdHQ0AiI6ORkFBAY4cOaLfZ8+ePdDpdIiKitLvk5SUhPLycv0+CQkJCA0NrbG5jIgsxxJDt2sj5ZFRUumLY+3ReERSc099iMaMGWOwrXPnzkhNTcWMGTNMOlZJSQkyMjL0z8+fP4+jR4/Cy8sLgYGBmDFjBhYtWoTWrVsjJCQEc+fOhb+/v34tNY1Gg759+2LixIlYs2YNysvLMXXqVIwYMQL+/v4AgFGjRmHBggWYMGECZs+ejb/++gsffPAB3nvvPVNPnYjMxBbNRVLLxkiNtUfjEUmNyQHRncFQJXd3d3z66acmHevw4cN46KGH9M8r++yMGzcO69atw8svv4zS0lJMmjQJBQUF6Nq1K3bu3AkXFxf9azZt2oSpU6eid+/eUCqVGDp0KOLi4vTlarUaP//8M6ZMmYKOHTvC29sbr7/+Oofck8WxY+rd2SJAqa0JS86s3dmdSGpMDojWr19fY5lCoagxYKpOz549qww5vfN4CxcuxMKFC2vcx8vLC1988cVd36d9+/b6dVmIrMHRO6aaCwMU6bDFaDwiKTE5IJo+fbrB8/Lycly7dg3Ozs5o0KCBSQERkSNix1SyR5bo7M4BEeYj1wVXrcnkn2NXr141eJSUlODUqVPo2rUrNm/ebIk6EtkNdky1f3JdTsUWnd2JpMQsd3br1q2xZMmSKtkjIrmx9izMROYk5dF4RJZ2TxMzVnugevVw5coVcx2OyC7JpWMq0/eOi6PxSK5MDoi+++47g+dCCGRlZWHlypXo0qWL2SpGZI+k3jGVgQwZg53d7RNHttaNyQFR5RxAlRQKBZo0aYJevXrh3XffNVe9iOySrZaJICJ548jWujM5IGJnUKKa2XKZCHNhFonIvtj7yFapjEY0Wx8iIqmzVjpZqot2kmOTypcKWdedI1u1Wi0AwMPDA2q1GmlpaYiPj0doaKhNfojZ031pVEB0+6rvtVm+fPk9V4bIUqydTmbHVCKyBi65Yj5GBUSpqalGHezOPwY5Jil23LtbM4+t0snsmEpEliaXka3WYFRAtHfvXkvXg+yEvXXck3o6mYioLqQ+stWeGP0NcO7cubuuO0bSV9cZeCszLampqVCr1QgICNBnWuLi4pCenm7S8Srblnv27GmxjrucKJGIHFnlyFatVlvlO7pyZKtGo+HIViMYHRC1bt1a3zkUAIYPH46cnByLVIqkx16XpDAmnVxWVibpdLJcl5IgotpxyRXzMfoK3Rl5/vjjjygtLTV7hUia7DXTcns6uTpMJxPRvajsS3nu3DlcuHDBpj8GueSKeXDYPRnFXjvucaJEIjI3Kfal5MjWujM6IFIoFNVmBkge7LXjniNMlEhE0mGJUavmGrnLka11Y3RAJITA+PHjoVKpAABlZWWYPHlylYzBtm3bzFtDkgR7zrRwokQiMgdLjFqVYrZJrowOiMaNG2fw/IknnjB7ZUi67D3TwnSy9dnTDLVExjD3JIj2vuSGozE6IFq7dq0l60F2wN4zLUwnE5ExagrmzdmXknOkSQ87VZNJmGkhIrkyZ19KLrkhPfwWI5NVZlpatGjBTAsRyYY5J0F0hDnSHA2/yYiIiIxgzkkQOUea9DAgIiIiMpK5JkHkkhvSwz5ERLcx13wgRObE+1JazNGX0t5H7joiBkRE/4fzgZAU8b6UJnOMWrX3kbuOhgERETgfCEkT70vHx5G70sErTrJ353wgbm5ucHJygoeHB8LDw5GXl4f4+HibLt5I0lI5T03Pnj3h5ORkkffgfSkf9jpyV0oL3JoDM0Qke5wPxDE42szYvC9JyhyxKdc+wlAiC+J8ICRFvC9JqiqbclNTU6FWqxEQEKBvyo2Li0N6erqtq3hPGBCR7HE+EJIi3pckRY7clMuAiGSP84GQFPG+JCkypSnX3jAgItkzdfZZR+tISNJkzlmRiUxxt884R27KZadqG+Jka9Jh7HwgjtiRkKRLDvPUOFpneHtX22ecORe4lRoGRDbCL1bpqW0+EM4JQ7bAeWrIWoz5jAsNDUVYWBhSU1MRFhZm8PrKptzIyEi7bMrl/ygbcNQe+o6gpvlALNWRkM1vxrHGvD9SZq/z1JD9MPYzDoDDNuUyQ2Rld950Wq0WAODh4QG1Wo20tDTEx8cjNDTULm8oR2WJOWGYJSR7xqYux2LKZ5yjNuUyILIyTrZmn4zpSHj58mWjOxKy+Y2IpMTUzzhHbMq135rbKUfuoe/IzDknjCPP40FE9ulePuMcrSnXvmtvhzjZmn0y55wwjjyPBxmSe98nsh+c94oBkdXxprNP5pwThllCIpIaznvFPkRWV3nTabVapKWlQalUwsXFBYWFhbhy5Yosbjp7Za6OhI48jweZjp2TSSoctbO0sRgQ2YDcbzp7Zo6OhJVZQkecx0PqGHwQ3Z0jdpY2FgMiG5HzTWfvKjsSArinjoTMEhKRlNX1M85eMSCyIbnedMQsYSVmbIiMw6WeLI8BEZGNMEtIRMbgJK7WwYCIyIaYJZQe/hInKeEkrtbDgIiITOaoQQN/iUtTRUUFfv31VwBAt27dZDOnE5d6si4GRERkEkcNGvhL/H/Yt0sauNSTdTGkJJKRyszOuXPncOHCBZOXB6kMGlJTU6FWqxEQEKAPGuLi4pCenm6hmlsWl1MhKeIkrtbFDBGRTBib2ampOcyR0/f8JU5SxElcrYsBEZEMGNscdLegydXV1WGDBlNX+iayBk7ial0MiIgcnLGZHZ1Oh5UrV9YYNPXt29dhgwb+Eicp4iSu1sWAiByCFEc9SaVjqjHNQWlpacjPz79r0HTgwAGoVCqHDBr4S5ykipO4Wg8DIrJ7jjrqyVyMaQ46ffo0CgoK0Lp16xqDpqysLPj6+uLChQsOFzTwlzhJGSdxtQ4GRDIixSxKXXGodO2MaQ6q7DRdW3NY165dUVJS4pBBA3+Jk5RxElfLY0AkE46YRXHkUU/mZExzkEajQVZWVq3NYR06dECbNm0cNmjgL3Ei+WJAJAOOmkWRy1Dpumb2jGkOGjduHL799luj+tAolUqHDhr4S5xInhgQOThHzqLIYai0uTJ7xjQHKZVKo/vQMGggIkfDgMjBOXIWxdGHSps7s1dbcxD70BDZN6mMbLVXDIgcnCNnURx5qLSlMnu1ZXbYh4aI5IoBkYNz5CyKIw+VtmVmj81hRCRH/KRzcJVZFK1WCyGEQdntI4zsMYsC/K+ZJyIiAoWFhdBqtcjPz0dkZGSVJqXKdHLPnj3h5ORkw1rXjos6EhFZFzNEDs6RsyiVHLGZx5Eze0REUsSASAbM3VlWih33HK2Zx5H7RxERSREDIplwxCyKI5NDZo+ISEoYEMmIo2VRHB2HwRMRWQ8DIiIJY2aPiMg6JP+pWlxcjBkzZiAoKAiurq7o3LkzDh06pC8fP348FAqFwaNv374Gx8jPz8fo0aPh4eEBT09PTJgwASUlJdY+FaJ7UpnZa9GiBTN7RBJhT6NWyTiSzxA9/fTT+Ouvv7Bhwwb4+/tj48aNiImJQVpaGpo1awYA6Nu3L9auXat/jUqlMjjG6NGjkZWVhYSEBJSXl+PJJ5/EpEmT8MUXX1j1XIiIiEiaJB0QXb9+Hd988w2+/fZbdO/eHQAwf/587NixA6tXr8aiRYsA/BMAVfaNuVN6ejp27tyJQ4cOoVOnTgCAFStW4JFHHsE777wDf39/65wMERFRNaQ4cleOJJ17v3XrFioqKuDi4mKw3dXVFfv379c/37dvH5o2bYrQ0FA8++yz+Pvvv/VlycnJ8PT01AdDABATEwOlUomUlJQa3/vGjRsoKioyeBAREZFjknRA5O7ujujoaLzxxhu4cuUKKioqsHHjRiQnJyMrKwvAP81l69evx+7du/H2228jMTER/fr1Q0VFBQAgOzsbTZs2NThuvXr14OXlhezs7Brfe/HixVCr1fpHQECA5U6UiIiIbErSAREAbNiwAUIINGvWDCqVCnFxcRg5cqS+Y+mIESMwcOBAtGvXDoMHD8b333+PQ4cOYd++fXV63zlz5qCwsFD/qFxck4iIiByP5AOili1bIjExESUlJdBqtfj9999RXl6OFi1aVLt/ixYt4O3tjYyMDACAr68vcnNzDfa5desW8vPza+x3BPzTL8nDw8PgQURERI5J8gFRJTc3N/j5+eHq1avYtWsXBg0aVO1+ly5dwt9//w0/Pz8AQHR0NAoKCnDkyBH9Pnv27IFOp0NUVJRV6k5ERETSJulRZgCwa9cuCCEQGhqKjIwMvPTSSwgLC8OTTz6JkpISLFiwAEOHDoWvry/Onj2Ll19+Ga1atUJsbCyAfya269u3LyZOnIg1a9agvLwcU6dOxYgRIzjCjIiIiADYQYaosLAQU6ZMQVhYGMaOHYuuXbti165dqF+/PpycnHDs2DEMHDgQbdq0wYQJE9CxY0f8+uuvBnMRbdq0CWFhYejduzceeeQRdO3aFR9//LENz4qIiIikRPIZomHDhmHYsGHVlrm6umLXrl21HsPLy4uTMJJROB8ISRHvSyLLk3yGiIiIiMjSGBARERGR7Em+yYyIiIikx9GacpkhIiIiItljQERERESyx4CIiIiIZI99iIhM5Gjt5kREt5PrZxwzRERERCR7zBDZkFyjcCIiIqlhQERkQwyKiYikgU1mREQkaTqdDtnZ2Th37hwuXLgAnU5n6yqRA2KGiIiIJCs9PR3ffPMNkpKScPPmTezduxcajQZDhgyBRqOxdfXIgTAgIiIiSUpPT0dcXBxyc3OhVqvh4uKCxo0bIzU1FVqtFtOmTWNQRGbDJjMiIpIcnU6H7du3Iy8vD+Hh4XBzc4OTkxM8PDwQHh6OvLw8xMfHs/mMzIYBERERSU5mZiZOnjyJgIAAKBQKgzKFQoHmzZsjPT0dmZmZNqohORoGREREJDnFxcUoKyuDm5tbteVubm4oKytDcXGxlWtGjooBERERSY67uztcXFxQWlpabXlpaSlcXFzg7u5u5ZqRo2JAREREkhMYGIiwsDBotVoIIQzKhBC4dOkSNBoNAgMDbVRDcjQcZUZERJKjVCoxZMgQaLVapKWlQalUwsXFBYWFhbhy5Qq8vb0xePBgKJXy+F3PSVwtTx53EhER2R2NRoNp06YhIiIChYWF0Gq1yM/PR2RkJIfck9kxQ0REdBv+EpcWjUaD2bNno1WrVrh27Rp69OiBkJAQ2WSGyHoYEBERkaQplUr4+voCAIKDgxkMkUUwICIikzGLQkSOhgERkYwwkCEiqh7zjkRERCR7zBARkQFmkYhIjpghIiIiItljQERERESyxyYzcghs5iEiorpgQCQjDBqIiIiqxyYzIiIikj1miIgkjpk9IiLLY4aIiIiIZI8BEREREckeAyIiIiKSPQZEREREJHsMiIiIiEj2GBARERGR7DEgIiIiItljQERERESyx4CIiIiIZI8BEREREckeAyIiIiKSPQZEREREJHsMiIiIiEj2GBARERGR7DEgIiIiItljQERERESyx4CIiIiIZI8BEREREckeAyIiIiKSPQZEREREJHsMiIiIiEj2GBARERGR7DEgIiIiItljQERERESyx4CIiIiIZI8BEREREckeAyIiIiKSPQZEREREJHsMiIiIiEj2GBARERGR7DEgIiIiItmrZ+sKEBER3Y2TkxN69uxp62qQg2OGiIiIiGSPARERERHJHgMiIiIikj3JB0TFxcWYMWMGgoKC4Orqis6dO+PQoUP6ciEEXn/9dfj5+cHV1RUxMTE4c+aMwTHy8/MxevRoeHh4wNPTExMmTEBJSYm1T4WIiIgkSvIB0dNPP42EhARs2LABx48fR58+fRATE4PLly8DAJYuXYq4uDisWbMGKSkpcHNzQ2xsLMrKyvTHGD16NE6cOIGEhAR8//33SEpKwqRJk2x1SkRERCQxCiGEsHUlanL9+nW4u7vj22+/Rf/+/fXbO3bsiH79+uGNN96Av78/Zs2ahRdffBEAUFhYCB8fH6xbtw4jRoxAeno6wsPDcejQIXTq1AkAsHPnTjzyyCO4dOkS/P39japLUVER1Go1CgsL4eHhYf6TJSIiIrMz9vtb0hmiW7duoaKiAi4uLgbbXV1dsX//fpw/fx7Z2dmIiYnRl6nVakRFRSE5ORkAkJycDE9PT30wBAAxMTFQKpVISUmp8b1v3LiBoqIigwcRERE5JkkHRO7u7oiOjsYbb7yBK1euoKKiAhs3bkRycjKysrKQnZ0NAPDx8TF4nY+Pj74sOzsbTZs2NSivV68evLy89PtUZ/HixVCr1fpHQECAmc+OiIiIpELSAREAbNiwAUIINGvWDCqVCnFxcRg5ciSUSstWfc6cOSgsLNQ/tFqtRd+PiIiIbEfyAVHLli2RmJiIkpISaLVa/P777ygvL0eLFi3g6+sLAMjJyTF4TU5Ojr7M19cXubm5BuW3bt1Cfn6+fp/qqFQqeHh4GDyIiIjIMUk+IKrk5uYGPz8/XL16Fbt27cKgQYMQEhICX19f7N69W79fUVERUlJSEB0dDQCIjo5GQUEBjhw5ot9nz5490Ol0iIqKsvp5EBERkfRIfi2zXbt2QQiB0NBQZGRk4KWXXkJYWBiefPJJKBQKzJgxA4sWLULr1q0REhKCuXPnwt/fH4MHDwYAaDQa9O3bFxMnTsSaNWtQXl6OqVOnYsSIEUaPMCMiIiLHJvmAqLCwEHPmzMGlS5fg5eWFoUOH4s0330T9+vUBAC+//DJKS0sxadIkFBQUoGvXrti5c6fByLRNmzZh6tSp6N27N5RKJYYOHYq4uDhbnRIRERFJjKTnIZISzkNERERkfxxiHiIiIiIia5B8k5lUVCbSOEEjERGR/aj83q6tQYwBkZGKi4sBgBM0EhER2aHi4mKo1eoay9mHyEg6nQ5XrlyBu7s7FAqF2Y5bVFSEgIAAaLVa9k2yAl5v6+L1ti5eb+vi9baue73eQggUFxfD39//rpM6M0NkJKVSiebNm1vs+Jz80bp4va2L19u6eL2ti9fbuu7let8tM1SJnaqJiIhI9hgQERERkewxILIxlUqFefPmQaVS2boqssDrbV283tbF621dvN7WZenrzU7VREREJHvMEBEREZHsMSAiIiIi2WNARERERLLHgIiIiIhkjwGRja1atQrBwcFwcXFBVFQUfv/9d1tXySEkJSVhwIAB8Pf3h0KhQHx8vEG5EAKvv/46/Pz84OrqipiYGJw5c8Y2lbVzixcvxv333w93d3c0bdoUgwcPxqlTpwz2KSsrw5QpU9C4cWM0bNgQQ4cORU5Ojo1qbN9Wr16N9u3b6yeni46Oxk8//aQv57W2rCVLlkChUGDGjBn6bbzm5jN//nwoFAqDR1hYmL7ckteaAZENffnll3jhhRcwb948/PHHH+jQoQNiY2ORm5tr66rZvdLSUnTo0AGrVq2qtnzp0qWIi4vDmjVrkJKSAjc3N8TGxqKsrMzKNbV/iYmJmDJlCg4ePIiEhASUl5ejT58+KC0t1e8zc+ZM7NixA1999RUSExNx5coVPPbYYzastf1q3rw5lixZgiNHjuDw4cPo1asXBg0ahBMnTgDgtbakQ4cO4aOPPkL79u0NtvOam1fbtm2RlZWlf+zfv19fZtFrLchmHnjgATFlyhT984qKCuHv7y8WL15sw1o5HgBi+/bt+uc6nU74+vqKZcuW6bcVFBQIlUolNm/ebIMaOpbc3FwBQCQmJgoh/rm29evXF1999ZV+n/T0dAFAJCcn26qaDqVRo0bik08+4bW2oOLiYtG6dWuRkJAgevToIaZPny6E4P1tbvPmzRMdOnSotszS15oZIhu5efMmjhw5gpiYGP02pVKJmJgYJCcn27Bmju/8+fPIzs42uPZqtRpRUVG89mZQWFgIAPDy8gIAHDlyBOXl5QbXOywsDIGBgbzedVRRUYEtW7agtLQU0dHRvNYWNGXKFPTv39/g2gK8vy3hzJkz8Pf3R4sWLTB69GhkZmYCsPy15uKuNpKXl4eKigr4+PgYbPfx8cHJkydtVCt5yM7OBoBqr31lGd0bnU6HGTNmoEuXLvjXv/4F4J/r7ezsDE9PT4N9eb3v3fHjxxEdHY2ysjI0bNgQ27dvR3h4OI4ePcprbQFbtmzBH3/8gUOHDlUp4/1tXlFRUVi3bh1CQ0ORlZWFBQsWoFu3bvjrr78sfq0ZEBGR2UyZMgV//fWXQZs/mV9oaCiOHj2KwsJCfP311xg3bhwSExNtXS2HpNVqMX36dCQkJMDFxcXW1XF4/fr10/+7ffv2iIqKQlBQELZu3QpXV1eLvjebzGzE29sbTk5OVXrH5+TkwNfX10a1kofK68trb15Tp07F999/j71796J58+b67b6+vrh58yYKCgoM9uf1vnfOzs5o1aoVOnbsiMWLF6NDhw744IMPeK0t4MiRI8jNzUVkZCTq1auHevXqITExEXFxcahXrx58fHx4zS3I09MTbdq0QUZGhsXvbwZENuLs7IyOHTti9+7d+m06nQ67d+9GdHS0DWvm+EJCQuDr62tw7YuKipCSksJrfw+EEJg6dSq2b9+OPXv2ICQkxKC8Y8eOqF+/vsH1PnXqFDIzM3m9zUSn0+HGjRu81hbQu3dvHD9+HEePHtU/OnXqhNGjR+v/zWtuOSUlJTh79iz8/Pwsf3/XuVs23bMtW7YIlUol1q1bJ9LS0sSkSZOEp6enyM7OtnXV7F5xcbFITU0VqampAoBYvny5SE1NFRcvXhRCCLFkyRLh6ekpvv32W3Hs2DExaNAgERISIq5fv27jmtufZ599VqjVarFv3z6RlZWlf1y7dk2/z+TJk0VgYKDYs2ePOHz4sIiOjhbR0dE2rLX9euWVV0RiYqI4f/68OHbsmHjllVeEQqEQP//8sxCC19oabh9lJgSvuTnNmjVL7Nu3T5w/f14cOHBAxMTECG9vb5GbmyuEsOy1ZkBkYytWrBCBgYHC2dlZPPDAA+LgwYO2rpJD2Lt3rwBQ5TFu3DghxD9D7+fOnSt8fHyESqUSvXv3FqdOnbJtpe1UddcZgFi7dq1+n+vXr4vnnntONGrUSDRo0EAMGTJEZGVl2a7Sduypp54SQUFBwtnZWTRp0kT07t1bHwwJwWttDXcGRLzm5jN8+HDh5+cnnJ2dRbNmzcTw4cNFRkaGvtyS11ohhBB1zzMRERER2S/2ISIiIiLZY0BEREREsseAiIiIiGSPARERERHJHgMiIiIikj0GRERERCR7DIiIiIhI9hgQEVG1goOD8f7779u6Gmazb98+KBSKKusgOZKePXtixowZtq4GkV1iQEQkM1qtFk899RT8/f3h7OyMoKAgTJ8+HX///betq2Y21QUGnTt3RlZWFtRqtW0qRUSSxoCISEbOnTuHTp064cyZM9i8eTMyMjKwZs0a/aLC+fn5NqtbRUUFdDqdxY7v7OwMX19fKBQKi72HI7L034VIKhgQEcnIlClT4OzsjJ9//hk9evRAYGAg+vXrh19++QWXL1/Gq6++arB/cXExRo4cCTc3NzRr1gyrVq3SlwkhMH/+fAQGBkKlUsHf3x/Tpk3Tl9+4cQMvvvgimjVrBjc3N0RFRWHfvn368nXr1sHT0xPfffcdwsPDoVKp8Mknn8DFxaVKs9b06dPRq1cvAMDff/+NkSNHolmzZmjQoAHatWuHzZs36/cdP348EhMT8cEHH0ChUEChUODChQvVNpl98803aNu2LVQqFYKDg/Huu+8avG9wcDDeeustPPXUU3B3d0dgYCA+/vjju17jnj17Ytq0aXj55Zfh5eUFX19fzJ8/X19+4cIFKBQKHD16VL+toKAACoVCf30q67pr1y5ERETA1dUVvXr1Qm5uLn766SdoNBp4eHhg1KhRuHbtmsH737p1C1OnToVarYa3tzfmzp2L21doupe/S2Zm5l3PmcghmGVFNCKSvL///lsoFArx1ltvVVs+ceJE0ahRI6HT6YQQQgQFBQl3d3exePFicerUKREXFyecnJz0C4l+9dVXwsPDQ/z444/i4sWLIiUlRXz88cf64z399NOic+fOIikpSWRkZIhly5YJlUolTp8+LYQQYu3ataJ+/fqic+fO4sCBA+LkyZOipKRE+Pj4iE8++UR/nFu3bhlsu3Tpkli2bJlITU0VZ8+e1dcrJSVFCCFEQUGBiI6OFhMnThRZWVkiKytL3Lp1S7/g79WrV4UQQhw+fFgolUqxcOFCcerUKbF27Vrh6upqsChtUFCQ8PLyEqtWrRJnzpwRixcvFkqlUpw8ebLG69yjRw/h4eEh5s+fL06fPi0+//xzg9Xoz58/LwCI1NRU/WuuXr0qAIi9e/cKIf63OPGDDz4o9u/fL/744w/RqlUr0aNHD9GnTx/xxx9/iKSkJNG4cWOxZMkSg/du2LChmD59ujh58qTYuHGjaNCgQZ3/LqWlpTWeL5GjYEBEJBMHDx4UAMT27durLV++fLkAIHJycoQQ/wQDffv2Ndhn+PDhol+/fkIIId59913Rpk0bcfPmzSrHunjxonBychKXL1822N67d28xZ84cIcQ/X7wAxNGjRw32mT59uujVq5f++a5du4RKpdIHMtXp37+/mDVrlv75nauRCyGqBESjRo0SDz/8sME+L730kggPD9c/DwoKEk888YT+uU6nE02bNhWrV6+usS49evQQXbt2Ndh2//33i9mzZwshTAuIfvnlF/0+ixcvFgDE2bNn9dueeeYZERsba/DeGo1GH9QKIcTs2bOFRqMRQtTt70Lk6NhkRiQz4rbmk9pER0dXeZ6eng4A+H//7//h+vXraNGiBSZOnIjt27fj1q1bAIDjx4+joqICbdq0QcOGDfWPxMREnD17Vn88Z2dntG/f3uA9Ro8ejX379uHKlSsAgE2bNqF///7w9PQE8E+fljfeeAPt2rWDl5cXGjZsiF27dpncrJOeno4uXboYbOvSpQvOnDmDiooK/bbb66dQKODr64vc3Ny7HvvOc/Lz86v1NbUdx8fHBw0aNECLFi0Mtt153AcffNCgn1R0dLT+nOrydyFydPVsXQEiso5WrVpBoVAgPT0dQ4YMqVKenp6ORo0aoUmTJkYdLyAgAKdOncIvv/yChIQEPPfcc1i2bBkSExNRUlICJycnHDlyBE5OTgava9iwof7frq6uVTo533///WjZsiW2bNmCZ599Ftu3b8e6dev05cuWLcMHH3yA999/H+3atYObmxtmzJiBmzdvmnA1jFe/fn2D5wqFotZOxnd7jVL5z+/Q2wPT8vLyWo+jUCjuqS63q8vfhcjRMSAikonGjRvj4YcfxocffoiZM2fC1dVVX5adnY1NmzZh7NixBl+EBw8eNDjGwYMHodFo9M9dXV0xYMAADBgwAFOmTEFYWBiOHz+OiIgIVFRUIDc3F926dTO5rqNHj8amTZvQvHlzKJVK9O/fX1924MABDBo0CE888QQAQKfT4fTp0wgPD9fv4+zsbJDlqY5Go8GBAwcMth04cABt2rSpEiyYU2XAmZWVhYiICAAw6GBdVykpKQbPDx48iNatW8PJyanOfxciR8YmMyIZWblyJW7cuIHY2FgkJSVBq9Vi586dePjhh9GsWTO8+eabBvsfOHAAS5cuxenTp7Fq1Sp89dVXmD59OoB/RiN9+umn+Ouvv3Du3Dls3LgRrq6uCAoKQps2bTB69GiMHTsW27Ztw/nz5/H7779j8eLF+OGHH2qt5+jRo/HHH3/gzTffxOOPPw6VSqUva926NRISEvDbb78hPT0dzzzzDHJycgxeHxwcjJSUFFy4cAF5eXnVZlFmzZqF3bt344033sDp06fx+eefY+XKlXjxxRfv5dIazdXVFQ8++CCWLFmC9PR0JCYm4rXXXjPb8TMzM/HCCy/g1KlT2Lx5M1asWKH/m9X170LkyBgQEclI69atcfjwYbRo0QLDhg1Dy5YtMWnSJDz00ENITk6Gl5eXwf6zZs3C4cOHERERgUWLFmH58uWIjY0FAHh6euI///kPunTpgvbt2+OXX37Bjh070LhxYwDA2rVrMXbsWMyaNQuhoaEYPHgwDh06hMDAwFrr2apVKzzwwAM4duwYRo8ebVD22muvITIyErGxsejZsyd8fX0xePBgg31efPFFODk5ITw8HE2aNKm2f1FkZCS2bt2KLVu24F//+hdef/11LFy4EOPHjzfhit6bzz77DLdu3ULHjh0xY8YMLFq0yGzHHjt2LK5fv44HHngAU6ZMwfTp0zFp0iR9eV3+LkSOTCFM6WFJRERE5ICYISIiIiLZY0BEREREsseAiIiIiGSPARERERHJHgMiIiIikj0GRERERCR7DIiIiIhI9hgQERERkewxICIiIiLZY0BEREREsseAiIiIiGSPARERERHJ3v8H4QjDIinCgL0AAAAASUVORK5CYII=\n"
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Computing the flux using Frequentist (Maximum Likelihood) method"
      ],
      "metadata": {
        "id": "e8D6bttq8t6O"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "1.   Probability of measuring a flux $f_i$ given a true flux $F_\\text{true}$,\n",
        "$$P(f_i~|~F_{\\rm true}) = \\frac{1}{\\sqrt{2\\pi \\sigma_i^2}} \\exp{\\left[-\\frac{(f_i - F_{\\rm true})^2}{2 \\sigma_i^2}\\right]}$$\n",
        "2.   Likelihood function for the observed data points,\n",
        "$$\\mathtt{L}(f~|~F_\\text{true}) = \\prod_i P(f_i~|~F_{\\rm true})$$\n",
        "3. Taking the logarithm of the function to deal with very small values,\n",
        "$$\\log\\mathtt{L} = -\\frac{1}{2}\\sum_i\\left[\\log(2\\pi\\sigma_i^2)+\\frac{(f_i-F_\\text{true})^2}{\\sigma_i^2}\\right]$$\n",
        "4. Maximizing likelihood by taking the derivative w.r.t $F_\\text{true}$ and setting it to zero, the condition for maximum likelihood is, the estimate of the true value,\n",
        "$$F_\\text{est} = \\frac{\\sum w_if_i}{\\sum w_i};\\qquad w_i = \\frac{1}{\\sigma_i^2}$$,\n",
        "and the uncertainty on the estimate,\n",
        "$$\\sigma_\\text{est}=\\sqrt\\frac{1}{\\left(\\sum_i w_i\\right)}$$\n",
        "\n",
        "<ins>Frequentist inference</ins>: In the limit of large number of repeated observations, 68\\% of the observed flux values will lie within the range $F_\\text{est}\\pm\\sigma_\\text{est}$.\n",
        "\n",
        "[Normal distribution](https://en.wikipedia.org/wiki/68%E2%80%9395%E2%80%9399.7_rule#/media/File:Empirical_rule_histogram.svg)\n",
        "\n",
        "\n"
      ],
      "metadata": {
        "id": "Adm_BP3eKhZN"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "<ins>Algorithm</ins>\n",
        "\n",
        "1.   Compute inverse variance (sigma [squared](https://numpy.org/doc/stable/reference/generated/numpy.square.html)) and store in variable $w$\n",
        "2.   Estimate the maximum likelihood flux value and the (1-sigma) uncertainty (Numpy arrays are multiplied element-wise by a*b; Sum over the array using [np.sum](https://numpy.org/doc/stable/reference/generated/numpy.sum.html))\n",
        "3. Print the inferred result\n",
        "\n",
        "<mark>User input required!!</mark>"
      ],
      "metadata": {
        "id": "rg1FJygb8sxx"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# Compute inverse variance square\n",
        "w = 1.0/np.square(sigma)\n",
        "# Estimate the maximum likelihood flux value and the (1-sigma) uncertainty\n",
        "F_est = np.sum(w*flux)/np.sum(w)\n",
        "sigma_est = np.sqrt(1/np.sum(w))"
      ],
      "metadata": {
        "id": "ApGVLUl28jZw"
      },
      "execution_count": 6,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Print the frequentist inferred result\n",
        "print('Based on {0} measurements, 68% of repeated flux measurements will lie in the range {1:.0f}±{2:.0f}'.format(N,F_est,sigma_est))"
      ],
      "metadata": {
        "id": "q8OrM89ZfKtf",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "177ba0c4-8422-4694-ecf6-deabb948c9ca"
      },
      "execution_count": 7,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Based on 50 measurements, 68% of repeated flux measurements will lie in the range 1001±4\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Computing the flux using Bayesian method"
      ],
      "metadata": {
        "id": "U8NFy4DK82gU"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "1.   From Bayes' theorem, the posterior probability of the true flux, $F_\\text{true}$,\n",
        "$$P(\\text{F}_\\text{true}~|~\\text{f}_\\text{obs}) = \\frac{P(\\text{f}_\\text{obs}~|~\\text{F}_\\text{true})\\cdot P(\\text{F}_\\text{true})}{\\left[\\sum_\\text{F}P(\\text{f}_\\text{obs}~|~\\text{F}_\\text{true})\\cdot P(\\text{F}_\\text{true})\\right]}$$\n",
        "\n",
        "2.   $P(\\text{F}_\\text{true})$ is the prior knowledge of the true flux which, in the absence of prior knowledge, may be taken as a uniform distribution over any flux value. $P(\\text{F}_\\text{true})\\propto1$\n",
        "\n",
        "3. $P(\\text{f}_\\text{obs}~|~\\text{F}_\\text{true})$ is the likelihood function\n",
        "$$\\log\\mathtt{L} = -\\frac{1}{2}\\sum_i\\left[\\log(2\\pi\\sigma_i^2)+\\frac{(f_i-F_\\text{true})^2}{\\sigma_i^2}\\right]$$\n",
        "\n",
        "\n",
        "4. Since the denominator term corresponds to a normalization, and is independent of $F_\\text{true}$,\n",
        "$$P(\\text{F}_\\text{true}~|~\\text{f}_\\text{obs}) \\propto P(\\text{f}_\\text{obs}~|~\\text{F}_\\text{true})\\cdot P(\\text{F}_\\text{true})$$,\n",
        "or equivalently,\n",
        "$$\\log P(\\text{F}_\\text{true}~|~\\text{f}_\\text{obs}) \\propto \\log P(\\text{f}_\\text{obs}~|~\\text{F}_\\text{true}) + \\log P(\\text{F}_\\text{true})$$\n"
      ],
      "metadata": {
        "id": "FerTo1FDgsS7"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "<ins>Algorithm</ins>\n",
        "1.   'Define' methods to compute the log of the<br>\n",
        "      (a) prior distribution for given model parameters<br>\n",
        "      Uniform if no prior knowledge<br>\n",
        "      (b) likelihood function for given model parameters and data<br>\n",
        "      Using methods for [squaring](https://numpy.org/doc/stable/reference/generated/numpy.square.html), [log](https://numpy.org/doc/stable/reference/generated/numpy.log.html), and [summing](https://numpy.org/doc/stable/reference/generated/numpy.sum.html)<br>\n",
        "      (c) non-normalized posterior distribution\n",
        "2.   Use Markov Chain Monte Carlo (MCMC) methods to sample the posterior distribution<br>\n",
        "      (a) Set the number of model parameters (`ndim`), the number of sampling chains (walkers – `nwalkers`), and the length of the chain (number of walker steps – `nsteps`)<br>\n",
        "      (b) Use `EnsembleSampler` from the emcee package to initialize the sampler, provide starting guesses for parameters, and `run_mcmc` method to sample the posterior<br>\n",
        "      (c) Combine the returned samples from the sampler after removing the initial burn-in lengths\n",
        "5.   Plot the posterior distribution<br>\n",
        "     A [histrogram](https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.hist.html) works well for visualizing the ditribution of samples.\n",
        "6.   Determine a characteristic value (e.g., [median](https://numpy.org/doc/stable/reference/generated/numpy.median.html)) and credible intervals (10%–90% [limits](https://numpy.org/doc/stable/reference/generated/numpy.percentile.html)) to characterize the distribution"
      ],
      "metadata": {
        "id": "EVPmICLgl-FG"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "<mark>User input required!!</mark>"
      ],
      "metadata": {
        "id": "1mqOHMyzztgQ"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "#log of prior distribution for given model parameters\n",
        "def log_prior(theta):\n",
        "    return\n",
        "\n",
        "# log likelihood function for given model parameters and data\n",
        "def log_likelihood(theta, flux, sigma):\n",
        "    return\n",
        "\n",
        "# log of posterior distribution\n",
        "def log_posterior(theta, flux, sigma):\n",
        "    return"
      ],
      "metadata": {
        "id": "wjWP-duL85ME"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "ndim =         # number of model parameters (F_true)\n",
        "nwalkers =     # number of MCMC walkers (say 30)\n",
        "nsteps =       # number of MCMC steps to take (a few thousands)"
      ],
      "metadata": {
        "id": "zrtTPuOlsFeX"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "<ins>emcee – The Markov chain Monte Carlo (MCMC) Ensemble sampler</ins>\n",
        "\n",
        "The emcee [Ensemble Sampler](https://emcee.readthedocs.io/en/stable/user/sampler/) is initialized by passing the number of walkers, number of parameters (dimensionality), log posterior probability function, a list of arguments passed to the log prob function in addition to the model parameter array, in our case the data.\n",
        "```\n",
        "sampler = emcee.EnsembleSampler(nwalkers, ndim, log_posterior, args=[])\n",
        "```\n",
        "The sampling is undertaken by calling the [run_mcmc](https://emcee.readthedocs.io/en/stable/user/sampler/#emcee.EnsembleSampler.run_mcmc) method of the emcee Ensemble Sampler\n",
        "```\n",
        "sampler.run_mcmc(initial_state, nsteps)\n",
        "```\n",
        "The `initial_state` is a set of starting guesses for the model parameters with some randomness introduced across the samplers/walkers. The array has a dimension of (nwalkers, ndim), and can be generated by using the [numpy.random](https://numpy.org/doc/stable/reference/random/generated/numpy.random.rand.html) methods.\n",
        "\n",
        "The results of the mcmc runs will be stored in the `sampler` variable and the samples can be accessed using the [get_chain](https://emcee.readthedocs.io/en/stable/user/sampler/#emcee.EnsembleSampler.get_chain) method."
      ],
      "metadata": {
        "id": "Pnf9fZigvrta"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import emcee"
      ],
      "metadata": {
        "id": "GAb1-aTBDZqp"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "<mark>User input required!!</mark>"
      ],
      "metadata": {
        "id": "59Ar1aMizvUh"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# Initial guesses for the model parameter – F_true\n",
        "initial_state = 2 * F_true * np.random.rand(nwalkers, ndim)\n",
        "\n",
        "# Initialize the sampler\n",
        "sampler =\n",
        "\n",
        "# Run the sampler\n",
        "sampler.run_mcmc(, progress=True)"
      ],
      "metadata": {
        "id": "MQCElAiatrRQ"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Plot the walker/sampler steps\n",
        "fig, ax = plt.subplots(1, figsize=(10, 3), sharex=True)\n",
        "samples = sampler.get_chain()\n",
        "ax.set_ylabel(\"flux\")\n",
        "ax.set_xlabel(\"step number\")\n",
        "p = ax.plot(samples[:,:,0], \"k\", alpha=0.3)"
      ],
      "metadata": {
        "id": "aMPCL5Vbz6KZ"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Use [`get_chain`](https://emcee.readthedocs.io/en/stable/user/sampler/#emcee.EnsembleSampler.get_chain) method of the sampler to discard non-steady initial values, reduce the total number of samples, and \"flatten\" by combining walker steps/samples for each parameter"
      ],
      "metadata": {
        "id": "ef_d1Vi9eiJT"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# Discard initial length of chains before stabilization\n",
        "flat_samples = sampler.get_chain()"
      ],
      "metadata": {
        "id": "ILpwFj-J08vL"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Plotting the posterior probability distribution"
      ],
      "metadata": {
        "id": "WK-RV6xU851V"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "<mark>User input required!!</mark>"
      ],
      "metadata": {
        "id": "9OPRQmklzw3E"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# plot a histogram of the posterior distribution samples\n",
        "plt.hist( , bins=, histtype=\"stepfilled\", alpha=0.3, normed=True)\n",
        "plt.xlabel(\" \")\n",
        "plt.ylabel(\" \")\n",
        "plt.show()"
      ],
      "metadata": {
        "id": "SJuOBVgt88hl"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Median value of the posterior distribution\n",
        "F_med =\n",
        "# 10% – 90% credible limit (10–90 percentiles)\n",
        "F_10 =\n",
        "F_90 =\n",
        "# Print\n",
        "print('The median flux value is {0:.0f} and the 10–90% credible limit is {1:.0f}–{2:.0f}'.format(F_med,F_10,F_90))"
      ],
      "metadata": {
        "id": "F2ZUisOs8Zxx"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Propagation of uncertainties"
      ],
      "metadata": {
        "id": "EJvH4FWG883X"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "Let us define a magnitude for the observed flux as\n",
        "$$M_\\text{arb} = -2.5\\log(F)$$\n",
        "What is the probable value of the magnitude, for the observed flux above?"
      ],
      "metadata": {
        "id": "tLC1jCshqt-g"
      }
    },
    {
      "cell_type": "code",
      "source": [],
      "metadata": {
        "id": "5DZXXcg78_ut"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# References\n",
        "\n",
        "\n",
        "1.   https://jakevdp.github.io/blog/2014/03/11/frequentism-and-bayesianism-a-practical-intro/\n",
        "2.   https://betterexplained.com/articles/an-intuitive-and-short-explanation-of-bayes-theorem/\n",
        "\n"
      ],
      "metadata": {
        "id": "72KVPLxCgCIo"
      }
    }
  ]
}