{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "73f0db79",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "iminuit version: 2.16.0\n"
     ]
    }
   ],
   "source": [
    "# Example of Bayesian parameter estimation using iminuit to find\n",
    "# posterior mode (MAP) estimators and MCMC for intervals.\n",
    "# pdf is a mixture of Gaussian (signal) and exponential (background),\n",
    "# truncated in [xMin,xMax].\n",
    "# G. Cowan / RHUL Physics / December 2022\n",
    "\n",
    "import numpy as np\n",
    "import scipy.stats as stats\n",
    "from scipy.stats import truncexpon\n",
    "from scipy.stats import truncnorm\n",
    "from scipy.stats import chi2\n",
    "from scipy.special import logsumexp\n",
    "from scipy.signal import correlate\n",
    "import iminuit\n",
    "from iminuit import Minuit\n",
    "import matplotlib.pyplot as plt\n",
    "plt.rcParams[\"font.size\"] = 14\n",
    "print(\"iminuit version:\", iminuit.__version__)  # need 2.x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "6c1df433",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define pdf\n",
    "np.random.seed(seed=1234567)        # fix random seed\n",
    "theta = 0.2                         # fraction of signal\n",
    "mu = 10.                            # mean of Gaussian\n",
    "sigma = 2.                          # std. dev. of Gaussian\n",
    "xi = 5.                             # mean of exponential\n",
    "xMin = 0.\n",
    "xMax = 20.\n",
    "\n",
    "def f(x, par):\n",
    "    theta   = par[0]\n",
    "    mu      = par[1]\n",
    "    sigma   = par[2]\n",
    "    xi      = par[3]\n",
    "    fs = stats.truncnorm.pdf(x, a=(xMin-mu)/sigma, b=(xMax-mu)/sigma, loc=mu, scale=sigma)\n",
    "    fb = stats.truncexpon.pdf(x, b=(xMax-xMin)/xi, loc=xMin, scale=xi)\n",
    "    return theta*fs + (1-theta)*fb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "1428ade9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate data\n",
    "numVal = 400\n",
    "xData = np.empty([numVal])\n",
    "for i in range (numVal):\n",
    "    r = np.random.uniform();\n",
    "    if r < theta:\n",
    "        xData[i] = stats.truncnorm.rvs(a=(xMin-mu)/sigma, b=(xMax-mu)/sigma, loc=mu, scale=sigma)\n",
    "    else:\n",
    "        xData[i] = stats.truncexpon.rvs(b=(xMax-xMin)/xi, loc=xMin, scale=xi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "872cae6f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Specify parameters, initial values, limits, etc.\n",
    "parin   = np.array([theta, mu, sigma, xi]) # initial values (here = true values)\n",
    "parname = [r'theta', r'mu', r'sigma', r'xi']\n",
    "parname_latex = [r'$\\theta$', r'$\\mu$', r'$\\sigma$', r'$\\xi$']\n",
    "parstep = np.array([0.1, 1., 1., 1.])      # initial step sizes for iminuit\n",
    "parfix  = [False, False, False, False]     # change these to fix/free parameters\n",
    "numFreePar = sum(1 for i in parfix if parfix[i] == False)\n",
    "parlim  = [(0.,1.), (None, None), (0., None), (0., None)]    # set limits"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "bbb208cc",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Negative log-likelihood\n",
    "def negLogL(par):\n",
    "    fx = f(xData, par)\n",
    "    return -np.sum(np.log(fx))\n",
    "\n",
    "# Prior pdf\n",
    "def prior(par):\n",
    "    theta    = par[0]\n",
    "    mu       = par[1]\n",
    "    sigma    = par[2]\n",
    "    xi       = par[3]\n",
    "    pi_theta = 1. if theta >= 0. and theta <= 1. else 0.\n",
    "    pi_mu    = 1. if mu >= 0. else 0.\n",
    "    pi_sigma = 1. if sigma > 0. else 0.\n",
    "    pi_xi    = 1. if xi > 0. else 0.\n",
    "    piArr = np.array([pi_theta, pi_mu, pi_sigma, pi_xi])\n",
    "    pi = np.product(piArr[np.array(parfix) == False])   # exclude fixed par\n",
    "    return pi\n",
    "    \n",
    "# Negative log of posterior pdf\n",
    "def negLogPost(par):\n",
    "    return negLogL(par) - np.log(prior(par))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "e53c149d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "par index, name, MAP estimate, standard deviation:\n",
      "   0 theta       =  0.197936  +/-  0.057880\n",
      "   1 mu          =  9.309286  +/-  0.555300\n",
      "   2 sigma       =  2.347174  +/-  0.509533\n",
      "   3 xi          =  5.054644  +/-  0.534132\n"
     ]
    }
   ],
   "source": [
    "# Find maximum of posterior and its covariance\n",
    "m = Minuit(negLogPost, parin, name=parname)\n",
    "m.errors = parstep\n",
    "m.fixed = parfix\n",
    "m.limits = parlim\n",
    "m.errordef = 0.5                          # errors from lnL = lnLmax - 0.5\n",
    "m.migrad()                                # minimize -logPost\n",
    "MAP = m.values                            # posterior mode\n",
    "sigmaMAP = m.errors                       # standard deviations\n",
    "cov = m.covariance                        # covariance matrix\n",
    "rho = m.covariance.correlation()          # correlation coeffs.\n",
    "    \n",
    "print(r\"par index, name, MAP estimate, standard deviation:\")\n",
    "for i in range(m.npar):\n",
    "    if not m.fixed[i]:\n",
    "        print(\"{:4d}\".format(i), \"{:<10s}\".format(m.parameters[i]), \" = \",\n",
    "         \"{:.6f}\".format(MAP[i]), \" +/- \", \"{:.6f}\".format(sigmaMAP[i]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "20a67395",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgYAAAFzCAYAAABFOMFPAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAAsTAAALEwEAmpwYAABHrklEQVR4nO3deXwV5dn/8c9FQthXWQSCLIIsCSFAZFHABVEWFREVUBFFann602Ktrfbpos/TRWvtU2y1C3XBrYqoVFQWRUEoRREQEUQQECRAAdl3Erh+f8zhGLKebOdk+b5fr3klZ+aeme9wEs6VmXvuMXdHREREBKBKrAOIiIhI2aHCQERERMJUGIiIiEiYCgMREREJU2EgIiIiYSoMREREJCw+1gFKS6NGjbx169axjiEiIhIVy5Yt+8bdGxd3OxW2MGjdujVLly6NdQwREZGoMLPNJbEdXUoQERGRMBUGIiIiEqbCQERERMJUGIiIiEiYCgMREREJU2EgIiIiYSoMREREJEyFgYiIiISpMBAREZEwFQYiIiISpsJAREREwlQYiIiISJgKAxEREQlTYSAiIiJhKgxEREQkTIWBiIiIhKkwEBERkTAVBiIiIhKmwkBERETCVBiIiIhImAoDERERCVNhICIiImEqDERERCRMhYGIiIiEqTAQERGRMBUGIiJ5MDPGjBkTfp2ZmUnjxo258sorz2g3bNgw+vTpc8a8Bx98kBYtWpCamkpycjIzZszIsf0pU6bQuHFjUlNTw9Pnn3+eZ57f/OY3Z7y+4IILinJYOcyfP59///vfJbItKf8qbGGQeSrWCUSkvKtVqxarVq3i6NGjALz77ru0aNHijDb79u1j+fLl7Nu3j6+++uqMZT/4wQ9YsWIF06ZNY9y4cZw6lfM/ppEjR7JixYrw1Llz5zzzZC8MSurDvCiFQWZmZonsW8qeqBUGZjbIzNaa2Xozuz+X5R3NbLGZHTeze7Mtq29mr5rZF2a2xsz6ZF8/u73HSjK9iFRWgwcP5u233wbgpZdeYvTo0Wcsf+2117jqqqsYNWoUL7/8cq7b6NSpE/Hx8XzzzTcR7XP79u30798/fLZh4cKF3H///Rw9epTU1FRuuukmAGrXrg0EH+wXXXQRN9xwA+eddx73338/L774Ij179qRLly5s2LABgDfffJNevXrRrVs3LrvsMnbs2MGmTZv461//yh/+8AdSU1NZuHAhmzdvZsCAAaSkpDBgwAC+/vprAG699VbuueceLrnkEu677z4++OCD8JmObt26cfDgwcL/A0vZ4+6lPgFxwAagLZAAfAp0ztamCXA+8Gvg3mzLngXGh75PAOoXtM/arXv4qVMuIlJktWrV8k8//dRHjBjhR48e9a5du/q8efN86NCh4TYDBgzwBQsW+Nq1a71Lly7h+Q888ID/7ne/c3f3Dz/80Js1a+ansv2n9Mwzz3ijRo28a9eu4enIkSP+6KOP+q9+9St3d8/MzPQDBw6E82TP5+4+b948r1evnm/bts2PHTvmzZs391/84hfu7j5p0iSfOHGiu7vv2bMnnOHvf/+733PPPTmyurtfeeWVPmXKFHd3f+qpp3zYsGHu7j527FgfOnSoZ2Zmhtv961//cnf3gwcPekZGRqH/jaXkAEu9BD6z46NUf/QE1rv7RgAzexkYBoQvprn7TmCnmQ3NuqKZ1QX6A7eG2p0AThS0wxMn4eNt0LNFQS1FRPKWkpLCpk2beOmllxgyZMgZy3bs2MH69evp27cvZkZ8fDyrVq0iOTkZgD/84Q+88MIL1KlTh6lTp2JmObY/cuRIHn/88TPmnX/++YwbN46MjAyuueYaUlNTC8x5/vnn06xZMwDOPfdcLr/8cgC6dOnCvHnzAEhPT2fkyJFs376dEydO0KZNm1y3tXjxYl5//XUAxowZw49//OPwsuuvv564uDgALrzwQu655x5uuukmrr32WhITEwvMKWVftC4ltAC2ZHmdHpoXibbALuAZM/vEzJ40s1oFrVTFYOrqwgcVEcnu6quv5t57781xGWHq1Kns3buXNm3a0Lp1azZt2nTG5YTTfQwWLlxIv379It5f//79WbBgAS1atGDMmDE899xzBa5TrVq18PdVqlQJv65SpUq4P8Bdd93FnXfeyWeffcbf/vY3jh2L7Jpr1oKmVq1v//u9//77efLJJzl69Ci9e/fmiy++iGh7UrZFqzDIWSaDR7huPNAd+Iu7dwMOAzn6KACY2R1mttTMllbzY7z1JRw4XrTAIiKnjRs3jl/84hd06dLljPkvvfQSs2fPZtOmTWzatIlly5bl2c+gMDZv3kyTJk34zne+w+23387y5csBqFq1KhkZGUXe7v79+8OdJ5999tnw/Dp16pzRP+CCCy4IH8eLL75I3759c93ehg0b6NKlC/fddx9paWkqDCqIaBUG6UDLLK8TgW2FWDfd3T8KvX6VoFDIwd0nu3uau6c1b1CdY5nw5roiZxYRASAxMZGJEyeeMW/Tpk18/fXX9O7dOzyvTZs21K1bl48++ij7JvI0derUM25X/Pe//838+fPDHfpee+218L7vuOMOUlJSwp0PC+vBBx/k+uuvp1+/fjRq1Cg8/6qrrmL69Onhzod//OMfeeaZZ0hJSeH555/nsccey3V7kyZNIjk5ma5du1KjRg0GDx5cpFxStljQX6GUd2IWD6wDBgBbgY+BG909x8l+M3sQOOTuj2aZt5Cg8+Ha0PJa7v6j/PaZlpbmZ929lIQ4eHN0fi1FRETKPzNb5u5pxd1OVM4YuHsmcCcwB1gDvOLuq81sgplNADCzs80sHbgH+JmZpYc6HgLcBbxoZiuBVOA3OXaSi5FJsHInrNlVwgckIiJSQUXljEEspKWl+bsLl9LzKbipCzx4UawTiYiIlJ5ydcYgVhrUgCvOhdfXwDEN0iUiIlKgCl0YAIxKgv3H4Z0NsU4iIiJS9lX4wuCClpBYF15aFeskIiIiZV+FLwyqWNAJ8d/psGlfrNOIiIiUbRW+MAAY2RniDP6hswYiIiL5qhSFQdPaMLAtTPscjqsTooiISJ4qRWEAwS2Le47CbHVCFBERyVOlKQz6ngPn1IMXP4t1EhERkbKr0hQGVQxuTIaPtsKXe2KdRkREpGyqNIUBwPWdoWoV+IfOGoiIiOSqUhUGjWrC4HbwqkZCFBERyVWlKgwAbuwCB47DW3ocs4iISA6VrjDo3QLObQAv6HKCiIhIDpWuMLBQJ8RP/gOf63HMIiIiZ6h0hQHAdZ2hWhw8vzLWSURERMqWSlkY1K8OwzrA9C+CJy+KiIhIoFIWBgC3dIWjmfDq57FOIiIiUnZU2sKgSxPo0Qye+xROeazTiIiIlA2VtjAAGNsVNu2HBZtjnURERKRsqNSFweB20LgmPPtprJOIiIiUDZW6MEiIC25dnLcJNu+LdRoREZHYq9SFAQSPY46rAs9rwCMREREVBk1rw6Bz4ZXVcDQj1mlERERiq9IXBhB0Qtx/HN5YG+skIlLWnDx5kokTJ5KUlESXLl3YuHFjsbY3e/ZsOnToQLt27Xj44YfzbDdu3DiaNGlCcnLyGfMfe+wxkpOTSUpKYtKkSeH5a9euJTU1NTzVrVs3vDyvdUoqb35tWrduTZcuXUhNTSUtLS3fnMXNKiXE3Svk1KNHD4/UqVPulz/vPuiF4HsRkdN++ctf+qRJk9zdffLkyf7DH/6wyNvKzMz0tm3b+oYNG/z48eOekpLiq1evzrXtBx984MuWLfOkpKTwvM8++8yTkpL88OHDnpGR4QMGDPB169blup+mTZv6pk2bIlpn3rx5Pnbs2CLlLahNq1atfNeuXXn+e5zOWZjjk9wBS70EPj91xoDg+Qm3psLn38CHW2OdRkTKisOHDzN9+nQmTpwIQJs2bVi/fn2Rt7dkyRLatWtH27ZtSUhIYNSoUbzxxhu5tu3fvz8NGzY8Y96aNWvo3bs3NWvWJD4+nosuuojp06fnWPe9997j3HPPpVWrVhGvU9S8hTmm/HIW5vikdKkwCBneERrWgKc/iXUSESkr5s6dy5YtW8KnvceNG5fjwxqgX79+Z5wePz3NnTv3jHZbt26lZcuW4deJiYls3Rr5XyPJycksWLCA3bt3c+TIEWbOnMmWLVtytHv55ZcZPXp0gev06tWL1NRUxo8fz4wZM8K558yZE3HegtqYGZdffjk9evRg8uTJeeYszPFJ6YqPdYCyonp8cIfC40uCWxdb1Y91IhGJtRUrVvC///u/TJgwAYDx48eTkpKSo93ChQsj2l5wtvdMZhZxnk6dOnHfffcxcOBAateuTdeuXYmPP/O/8RMnTjBjxgweeuihAtf56KOPAJg/fz5TpkxhypQphc5bUJtFixbRvHlzdu7cycCBA+nYsSP9+/fPkTPS45PSF7UzBmY2yMzWmtl6M7s/l+UdzWyxmR03s3tzWR5nZp+Y2VullXFMF4ivAs9owCMRAfbu3UvNmjUByMzM5J133uGqq67K0S7SMwaJiYln/AWcnp5O8+bNC5Xp9ttvZ/ny5SxYsICGDRvSvn37M5bPmjWL7t2707Rp04jXyUskeQtqc/r7Jk2aMHz4cJYsWZJnzuJklRJUEh0VCpqAOGAD0BZIAD4FOmdr0wQ4H/g1cG8u27gH+AfwViT7LEznw6wmznbv9IT7/mNFWl1EKpAnnnjC/+u//svd3R955BGfMGFCsbaXkZHhbdq08Y0bN4Y76q1atSrP9l999dUZnQ/d3Xfs2OHu7ps3b/YOHTr4nj17zlg+cuRIf/rppwu1TnHy5tfm0KFDfuDAgfD3ffr08VmzZuWZszhZpeQ6H0arMOgDzMny+ifAT/Jo+2D2wgBIBN4DLi3twmDlf9zPmeT+92VFWl1EKpA9e/Z4r169/Nxzz/Wbb77Zjxw5Uuxtvv32296+fXtv27at/+pXvzpj2eDBg33r1q3u7j5q1Cg/++yzPT4+3lu0aOFPPvmku7v37dvXO3Xq5CkpKT537twz1j98+LA3bNjQ9+3bd8b8vNbp2bOnd+3aNcc0e/bsAvNmzZpXmw0bNnhKSoqnpKR4586dw8vyylnQ8Un+SqowMM/l+lBJM7PrgEHuPj70egzQy93vzKXtg8Ahd380y7xXgYeAOgRFw5V57OcO4A6Ac845p8fmzUV7OtL102DbIVgwNhgVUUREpKwzs2Xunlbc7UTrYy+33jURVSRmdiWw092XFdTW3Se7e5q7pzVu3LiwGcPGdYP0A/BO8cYxERERKXeiVRikAy2zvE4EtkW47oXA1Wa2CXgZuNTMXijZeGe6vC0k1tWtiyIiUvlEqzD4GGhvZm3MLAEYBcyIZEV3/4m7J7p769B677v7zaUXNbh8cFtXWLINPttRmnsSEREpW6JSGLh7JnAnMAdYA7zi7qvNbIKZTQAws7PNLJ3g7oOfmVm6mdWNRr7c3JAEtRNg8vJYJRAREYm+qHQ+jIW0tDRfunRpsbbxq4XB5YQFtwaXFkRERMqq8tb5sFwalxo8R+Ep9TUQEZFKQoVBPprXgavPg5dXw/5jsU4jIiJS+lQYFOCO7nAkA55fGeskIiIipU+FQQE6NYaLWsGUT+FYZqzTiIiIlC4VBhG4ozvsOgLTv4h1EhERkdKlwiACF7aEpMbw9+VwqmLexCEiIgKoMIiIGXy3B2zYC+99Fes0IiIipUeFQYSGtofEOvDX4g2NICIiUqapMIhQfBUY3x2WboclW2OdRkREpHSoMCiEUUlwVg144uNYJxERESkdKgwKoUZVuL0bzN8Mq3bGOo2IiEjJU2FQSGNSoE4C/Fl9DUREpAJSYVBIdavBLSkw88vgLgUREZGKRIVBEYzrBglx8BedNRARkQpGhUERNKoJo5ODkRC3Hoh1GhERkZKjwqCI7ugefJ28PLY5RERESpIKgyJqUReGdwgeyfzNkVinERERKRkqDIphQhocz4QnddZAREQqCBUGxdCuIVx1Hjy7EvYcjXUaERGR4lNhUEx39YSjGTprICIiFYMKg2I676zgAUtTPoW9OmsgIiLlnAqDEvD9nnA4A578JNZJRKQ8eOyxx0hOTiYpKYlJkyblWH7s2DF69uxJ165dSUpK4oEHHjhj+ezZs+nQoQPt2rXj4Ycfzne7BW2rIHntK7vWrVvTpUsXUlNTSUtLK/L84uaVEuDuFXLq0aOHR9OEt9w7/9l979Go7lZEypnPPvvMk5KS/PDhw56RkeEDBgzwdevWndHm1KlTfvDgQXd3P3HihPfs2dMXL17s7u6ZmZnetm1b37Bhgx8/ftxTUlJ89erVeW43v21lNW/ePB87duwZ8/LaV25atWrlu3btKvb8SPNKTsBSL4HPT50xKCETe8GhE/CUzhqIVCgXX3wxa9euBWD37t0kJycXa3tr1qyhd+/e1KxZk/j4eC666CKmT59+Rhszo3bt2gBkZGSQkZGBmQGwZMkS2rVrR9u2bUlISGDUqFG88cYbeW43v20VJK99labi5JWSocKghHRsBIPbwTMrYP+xWKcRkZKyfv162rdvD8DKlSvp0qVLjjb9+vUjNTU1xzR37twcbZOTk1mwYAG7d+/myJEjzJw5ky1btuRod/LkSVJTU2nSpAkDBw6kV69eAGzdupWWLVuG2yUmJrJ169Z8t5vXtgB69epFamoq48ePZ8aMGeHsc+bMyXNfuTEzLr/8cnr06MHkyZOLPL+gvFL64mMdoCKZ2BNmrQ/OGtzTJ9ZpRKS4Nm/eTIsWLahSJfgbauXKlaSkpORot3Dhwoi32alTJ+677z4GDhxI7dq16dq1K/HxOf8rjouLY8WKFezbt4/hw4ezatUqkpOTCc4Yn8nM8t1uXtsC+OijjwCYP38+U6ZMYcqUKeHtTps2Ldd95WbRokU0b96cnTt3MnDgQDp27Ej//v0LPb+gvFL6dMagBHVqDFecC0+vgH06ayBS7q1YseKMQmDZsmW5FgaFOWMAcPvtt7N8+XIWLFhAw4YNw2ckclO/fn0uvvhiZs+eDQR/tWc9w5Cenk7z5s0j2m72bRUkv31ld3p+kyZNGD58OEuWLCnS/OLklRJSEh0VIpmAQcBaYD1wfy7LOwKLgePAvVnmtwTmAWuA1cDESPYX7c6Hp63Z5d5qkvtv/xWT3YtICfqf//kfv/nmm93dfd26dV63bl3fvHlzsbe7Y8cOd3ffvHmzd+jQwffs2XPG8p07d/revXvd3f3IkSPet29ff/PNN93dPSMjw9u0aeMbN24MdwhctWpVntvNb1sFyW9fWR06dMgPHDgQ/r5Pnz4+a9asQs8v6Nglf5RQ58OoXEowszjgCWAgkA58bGYz3P3zLM32AN8Hrsm2eibwQ3dfbmZ1gGVm9m62dcuMjo2C0RCfXgG3pULjWrFOJCJFtWLFCmrUqEHXrl1JSUmhU6dOPPvss/z85z8v1nZHjBjB7t27qVq1Kk888QQNGjQAYMiQITz55JN88803jB07lpMnT3Lq1CluuOEGrrzySgDi4+N5/PHHueKKKzh58iTjxo0jKSkpz+2uXLkyz21B0Mfg+PHjOTL+9re/5YorrshzX6ezNm/enB07djB8+HAAMjMzufHGGxk0aBAbN24s1HyA7du355tXSp95LterSnwnZn2AB939itDrnwC4+0O5tH0QOOTuj+axrTeAx9393fz2mZaW5kuXLi1u9CLZuBcuex5u7Qq/uCgmEUSkBLRr145PPvmEOnXqxDqKSIHMbJm7pxXcMn/R6mPQAsja7TY9NK9QzKw10A34qGRilY62DWBEJ3jhM9h+MNZpRKQoDh48SJUqVVQUSKUTrcIgt26shTpVYWa1gdeAu939QB5t7jCzpWa2dNeuXUWIWXK+3wtOOfwxZ38aESkH6tSpw7p162IdQyTqolUYpBN0IjwtEdgW6cpmVpWgKHjR3V/Pq527T3b3NHdPa9y4cZHDloSWdWF0MrzyOWzeF9MoIiIiEYtWYfAx0N7M2phZAjAKmBHJihbcNPsUsMbd/68UM5a4O8+HOINJZfrCh4iIyLeiUhi4eyZwJzCH4LbDV9x9tZlNMLMJAGZ2tpmlA/cAPzOzdDOrC1wIjAEuNbMVoWlINHIXV9PaMLYrTP8C1u2OdRoREZGCReWuhFiI5V0JWe05Cv2mwIUtYbLuuBERkVJS3u5KqLQa1oDv9oA5G2BpxL0qREREYkOFQRSM7waNa8LDi6CCnqAREZEKQoVBFNSsCnf3go+3wXtfxTqNiIhI3lQYRMnIJGhbPzhrcPJUrNOIiIjkToVBlFSNgx9dAF/ugVfXxDqNiIhI7lQYRNHgdtDtbPi/D+FYZqzTiIiI5KTCIIrM4P4L4T+HYMqKWKcRERHJSYVBlPVOhEtbwxMfB2MciIiIlCUqDGLgJ33hcIaGShYRkbJHhUEMnHcW3NgFXlgZdEYUEREpK1QYxMgPegXjG/xmYayTiIiIfEuFQYycVRPu7Anvb4KFm2OdRkREJKDCIIZu6wot68IvF2rQIxERKRtUGMRQtfigI+La3TB1dazTiIiIqDCIuSHt4Pzm8PvFcPB4rNOIiEhlp8IgxszgF/3hm6Pwp49jnUZERCo7FQZlQEpTuL4TPP0JbNgb6zQiIlKZqTAoI+67EKrHw/98AO6xTiMiIpWVCoMyonEt+EFv+GAzvLsx1mlERKSyUmFQhtySAu0bwv8u0NMXRSqqLVu2cMkll9CpUyeSkpJ47LHHcrQ5duwYPXv2pGvXriQlJfHAAw+csfzkyZN069aNK6+8Mjxv7dq1pKamhqe6desyadKkiPaXn9mzZ9OhQwfatWvHww8/nGub1q1b06VLF1JTU0lLS4to/ccee4zk5GSSkpKYNGlSxNuTKHD3Cjn16NHDy6N/fe1+ziT3SR/GOomIlIZt27b5smXL3N39wIED3r59e1+9evUZbU6dOuUHDx50d/cTJ054z549ffHixeHlv//973306NE+dOjQXPeRmZnpTZs29U2bNkW0P3f3efPm+dixY3Nsp23btr5hwwY/fvy4p6Sk5Lpuq1atfNeuXbnmyG39zz77zJOSkvzw4cOekZHhAwYM8HXr1hW4PckfsNRL4PNTZwzKmAtbBrcw/nkppB+IdRoR+fTTT+nfvz+dO3emSpUqmFmOv+ALo1mzZnTv3h2AOnXq0KlTJ7Zu3XpGGzOjdu3aAGRkZJCRkYGZAZCens7bb7/N+PHj89zHe++9x7nnnkurVq0i2l9elixZQrt27Wjbti0JCQmMGjWKN954I+JjzWv9NWvW0Lt3b2rWrEl8fDwXXXQR06dPj3i7UrpUGJRBP+sXfP2VnqMgElPHjh1j5MiRPProo3z++ef89Kc/5d577+XBBx88o12/fv3OOI1/epo7d26+29+0aROffPIJvXr1yrHs5MmTpKam0qRJEwYOHBhuc/fdd/PII49QpUre/32//PLLjB49OqL99erVi9TUVMaPH8+MGTPC2efMmcPWrVtp2bJluG1iYmKuRYWZcfnll9OjRw8mT54cnp/X+snJySxYsIDdu3dz5MgRZs6cyZYtWwrcnkRHfKwDSE4t6sL/Oz8Y9Gj+Jri4dawTiVROc+fOpXv37vTs2ROAlJQUZs+eHf7r/bSFCwtfxR86dIgRI0YwadIk6tatm2N5XFwcK1asYN++fQwfPpxVq1axadMmmjRpQo8ePZg/f36u2z1x4gQzZszgoYceimh/H30UPP99/vz5TJkyhSlTpoSXTZs2Lcf2sx87wKJFi2jevDk7d+5k4MCBdOzYkf79++O53GJlZnTq1In77ruPgQMHUrt2bbp27Up8fHyB25Po0BmDMuq73eHcBvDz+eqIKBIrq1atokuXLuHXy5cvD5+Wz6qwZwwyMjIYMWIEN910E9dee22+GerXr8/FF1/M7NmzWbRoETNmzKB169aMGjWK999/n5tvvvmM9rNmzaJ79+40bdq0SPvLKjEx8Yy/5NPT02nevHmOdqfnNWnShOHDh7NkyZIC17/99ttZvnw5CxYsoGHDhrRv377A7UmUlERHhbI4ldfOh1ktCnVEfGRRrJOIVE6TJ0/2UaNGubv72rVr/bzzzvNvvvmmWNs8deqUjxkzxidOnJhnm507d/revXvd3f3IkSPet29ff/PNN89oM2/evFw7H44cOdKffvrpQu0vLxkZGd6mTRvfuHFjuPPgqlWrzmhz6NAhP3DgQPj7Pn36+KxZswpcf8eOHe7uvnnzZu/QoYPv2bOnwO1J/iihzoe6lFCGXdASru0If1sG13QMbmUUkegZPXo0M2bMIDk5mUaNGvHSSy9x1llnFWubixYt4vnnnw/fjgfwm9/8hiFDhjBkyBCefPJJvvnmG8aOHcvJkyc5deoUN9xwwxm3JublyJEjvPvuu/ztb3+LaH8Q9DE4fjzng1p++9vfcsUVV/D4449zxRVXcPLkScaNG0dSUhJAOOuxY8cYPnw4AJmZmdx4440MGjQIgPj4+DzXHzFiBLt376Zq1ao88cQTNGjQAIAdO3bkuT2JDvNcrgFVBGlpab506dJYxyi2b47Apc9Bp0bw8ojg2QoiIiLZmdkydy/2wA9R62NgZoPMbK2ZrTez+3NZ3tHMFpvZcTO7tzDrVmSNagaPZv5wK7z+RazTiIhIRReVwsDM4oAngMFAZ2C0mXXO1mwP8H3g0SKsW6GNTILuzYLbF/cejXUaERGpyKJ1xqAnsN7dN7r7CeBlYFjWBu6+090/BjIKu25FV8XgN5fA/mPwa41tICIipShahUELYEuW1+mheaW9boXRqTFMSINpa2Dh5linERGRiipahUFuXeYi7fUY8bpmdoeZLTWzpbt27Yo4XHnx/Z7B2Ab3vw+HT8Q6jYiIVETRKgzSgZZZXicC20p6XXef7O5p7p7WuHHjIgUty6rHw28vg60H4NHFsU4jIiIVUbQKg4+B9mbWxswSgFHAjCisW+Gc3xxu6QrPrIBl22OdRkREKpqoFAbungncCcwB1gCvuPtqM5tgZhMAzOxsM0sH7gF+ZmbpZlY3r3Wjkbus+vEF0LwO/HguHNdwySIiUoI0wFE5NX8TjH0D7uoJ9/aJdRoREYm1cjfAkZSsi1vDiE7wl6Xw2Y5YpxERkYpChUE59kB/OKsG3POunsAoIiIlQ4VBOVavOjxyGazbDf/3YazTiIhIRaDCoJy7uDXcmAyTl8HHkd4AKiIikgcVBhXAT/tBYl344Tsa+EhERIpHhUEFUDsBHh0IX++HhxbFOo2IiJRnKgwqiN6JMK4bPL8SFuhZCiIiUkSFLgzMrFboUchSxvz4AmjfEO55B3YfiXUaEREpjwosDMysipndaGZvm9lO4Atgu5mtNrPfmVn70o8pkageD38aBAeOw4/mQgUdu0pEREpRJGcM5gHnAj8Bznb3lu7eBOgHfAg8bGY3l2JGKYROjeH+C+G9r+C5lbFOIyIi5U18BG0uc/eM7DPdfQ/wGvCamVUt8WRSZLelwgeb4dcLoXcL6NAo1olERKS8KPCMwemiwMwmmZnl10bKBrPgLoU6CXDXbI2KKCIikStM58NDwAwzqwVgZpebmW6OK6Ma14LfXw5rdwdnDkRERCIRyaUEANz9Z2Z2IzDfzI4Dh4H7Sy2ZFNvFreH2VHhqBfRJhCHqJioiIgWI+IyBmQ0AvkNQEDQGvu/u+lu0jLu/L6Q2hR/PhU37Yp1GRETKusJcSvgp8HN3vxi4DphqZpeWSiopMQlx8PgQqGLwXzPV30BERPIXcWHg7pe6+79C338GDAZ+VVrBpOS0rAt/uBw+3wX/+0Gs04iISFkWyQBHed2JsB0YkF8bKTsGtIX/6gEvroLpX8Q6jUj58cEHH9CkSRPi4uJo06YNv//974u1vdatW9OlSxdSU1NJS0vLtc3s2bPp0KED7dq14+GHH454/ZMnT9KtWzeuvPLK8Lx9+/Zx3XXX0bFjRzp16sTixYsjzppfjoIyHTt2jJ49e9K1a1eSkpJ44IEHCswKMG7cOJo0aUJycnLEOaWEuXu+E8EAR3cB52SbnwBcCjwL3FrQdqI99ejRw+VMGSfdr3vFvcPj7mu/iXUakfJh6tSp/r3vfc/37t1bIttr1aqV79q1K8/lmZmZ3rZtW9+wYYMfP37cU1JSfPXq1RGt//vf/95Hjx7tQ4cODc+75ZZb/O9//7u7ux8/fjzX45g3b56PHTu2UDkKOqZTp075wYMH3d39xIkT3rNnT1+8eHG+Wd3dP/jgA1+2bJknJSXlui/JG7DUS+DzM5JLCV8CJ4HpZrbNzD43s42h+aOBP7j7lBKvWKTExVeBPw2GWlXhjrdg//FYJxIp+5599lkuu+wy6tWrF5X9LVmyhHbt2tG2bVsSEhIYNWoUb7zxRoHrpaen8/bbbzN+/PjwvAMHDrBgwQJuv/12ABISEqhfv36p5jjNzKhduzYAGRkZZGRkcPrkcm5ZT+vfvz8NGzaMeD9S8iIpDC5w9z8DBpxDcPmgu7u3cvfvuPuK0gwoJevs2vDnobDlAEycDSdPxTqRSNl21113MWbMGOrXr88//vGPHMv79etHampqjmnu3Lm5bs/MuPzyy+nRoweTJ0/OsXzr1q20bNky/DoxMZGtW7cWuP7dd9/NI488QpUq3/63vnHjRho3bsxtt91Gt27dGD9+PIcPHw4v79WrF6mpqYwfP54ZM2aEs8+ZM6fAHJEc08mTJ0lNTaVJkyYMHDiQXr165ZlVyo5IxjGYY2aLgabALcCnwOpSTSWlqlcLeKA//Hw+/N+H8KMLYp1IpGz64osv+PGPf8ybb77JxRdfTG7dqRYuLNxd24sWLaJ58+bs3LmTgQMH0rFjR/r37x9e7rk8/SzrfnNb/8CBAzRp0oQePXowf/78cNvMzEyWL1/On/70J3r16sXEiRN5+OGH+eUvfwnARx99BMD8+fOZMmUKU6ZMCa87bdq0fHNEckxxcXGsWLGCffv2MXz4cFatWsWmTZtyzSplR4GFgbv/0MzaAvOBNsDVQJKZnQBWufvI0o0opWFMCqzeBY9/DJ0bw1ANfiSSw9/+9jfuueceLrnkkjzb9OvXj4MHD+aY/+ijj3LZZZflmN+8eXMAmjRpwvDhw1myZMkZhUFiYiJbtmwJv05PTw+vk9f6u3fvZsaMGcycOZNjx45x4MABbr75Zh599FESExPDf6lfd911+XYizKqgHIU5pvr163PxxRcze/bsPLO+8MILEeWSKIi0MwJwXrbXtYHeJdHRoTQmdT4s2LEM92Evu3d8wn1N3n2hRCqtCRMm+EMPPVRi2zt06JAfOHAg/H2fPn181qxZZ7TJyMjwNm3a+MaNG8Od/latWhXx+vPmzTujQ1/fvn39iy++cHf3Bx54wO+9996IsuaXI5Jj2rlzZ7ij45EjR7xv377+5ptv5pv1tK+++kqdD4uAKHY+PF1ArMv2+pC7f1iSRYpEV7V4+NtQqJ0A33kLdh+JdSKRsuVHP/oR7777LsnJyQwcOJDt27cXa3s7duygb9++dO3alZ49ezJ06FAGDRoEwJAhQ9i2bRvx8fE8/vjjXHHFFXTq1IkbbriBpKSkAtfPy5/+9CduuukmUlJSWLFiBf/93/8dXna6j0H2ac6cOfnmyJo3r0zbt2/nkksuISUlhfPPP5+BAwfmuDUxN6NHj6ZPnz6sXbuWxMREnnrqqaL8U0sxmOdyPasiSEtL86VLl8Y6RrnwyX9g5KuQ3AT+cS1Uj/gJGiKVx9ixY7nhhhsYOnRorKOI5MrMlrl77oNjFIK6hArdzoY/XAHLtsO978KpilkrihTZW2+9xeHDh3PtMyBS0agwECDofPiTC+HNdfD7yAdGE6kUrrzySl599VWqVasW6ygipS5qhYGZDTKztWa23sxyPK7ZAn8MLV9pZt2zLPuBma02s1Vm9pKZVY9W7srkuz1gdHJwp8K0z2OdRkREYiEqhYGZxQFPEDx4qTMw2sw6Z2s2GGgfmu4A/hJatwXwfSDN3ZOBOGBUNHJXNmbwy4uh3zlw/3vwr69jnUhERKItWmcMegLr3X2ju58AXgaGZWszDHgudNfFh0B9M2sWWhYP1DCzeKAmsC1KuSudqnHw5yFwboNg2OSVO2KdSEREoilahUELYEuW1+mheQW2cfetwKPA18B2YL+7v5PbTszsDjNbamZLd+3aVWLhK5u61eD5a6BBDRj7BmzYG+tEIiISLdEqDHIbRzN73/dc25hZA4KzCW2A5kAtM7s5t524+2R3T3P3tMaNGxcrcGXXtDa8cE3wAzJmOmzPObCbiIhUQNEqDNKBllleJ5LzckBebS4DvnL3Xe6eAbwOaHT/KGjTAJ69Bg4chzH/hL1HY51IRERKW7QKg4+B9mbWxswSCDoPzsjWZgZwS+juhN4Elwy2E1xC6G1mNS14gscAYE2Ucld6yU3gyavg6/1w2ww4qEc1i4hUaFEpDNw9E7gTmEPwof6Ku682swlmNiHUbCawEVgP/B34Xmjdj4BXgeXAZ6HMOZ9VKqWmdyI8Phg+2wm3vgGHTsQ6kYiIlBYNiSwRm/kl3DkLujeDZ4dBrYRYJxIRkdM0JLJE3ZD28KfBsHx7cObgsM4ciIhUOCoMpFCGtofHBsHS7TBuBhzJiHUiEREpSSoMpNCuOg8mXQFLtgW3Mu5Xh0QRkQpDhYEUybAOQYfET3fAqFdh5+FYJxIRkZKgwkCKbGh7ePpq+GofXD8NthyIdSIRESkuFQZSLP1bwYvXwt5jMOIVWLc71olERKQ4VBhIsfVoBtOuC8a4vm4afJge60QiIlJUKgykRHRoBK9dD41qws3TYdrnsU4kIiJFocJASsw59WD6SOjVAu59F367CE5VzPGzREQqLBUGUqLqVYMpw+CmZPjzUvivt+GoxjoQESk3VBhIiasaB7++FH7eD+ZsgOGvwFd7Y51KREQiocJASoUZjO8enD34zyG48mV4+8tYpxIRkYKoMJBSdXFreHs0tG8I35sJ//MBnDgZ61QiIpIXFQZS6lrUhVeug9tS4ekVcMOrsGlfjEOJiEiuVBhIVCTEwYMXwZ+HwIY9MOhFeO5TqKBP/RYRKbdUGEhUDW0P79wMPVvAz+fDmH/CtoOxTiUiIqepMJCoa1YHnh0GD10Ky7bD5S/AS6s05oGISFmgwkBiwgxu7AJzboKkJnD/e3DNVFjxn1gnExGp3FQYSEydUw9evhYeuwK2HwqKg/vmwp6jsU4mIlI5qTCQmDODazrC+2OCsQ9eXQMXPQt/WgKHTsQ6nYhI5aLCQMqMOtXgZ/1g1o3Qszk8uhj6TYHJyzWssohItKgwkDLnvLPgqavhnyMhqTH8emFQIPxtGew/Fut0IiIVmwoDKbO6nQ0vDA8GRzq3IfzmX9DrKfjp+/DlnlinExGpmOJjHUCkIL1awNQRsHoXPLMCpn0OL3wGfVvCiE5w+blQOyHWKUVEKgbzCjr0XFpami9dujTWMaQU7D4SjHvw0mpIPwDV4mBgWxjWAfq3guoqd0WkEjKzZe6eVuztqDCQ8so9GCDpjbXw1pfBLY7V4qB3IvQ/JygS2jcM7noQEanoVBgUQIVB5ZJxEv6dDvM3wQebYcPeYH6TWtC1KXRpEkwpTaFRzdjldIdjmbD3WGg6CgeOB7dlHjoBB0NfT5wMpoyTcDz0NMr4Kt9OVeOgXjVoUB3qV4d61aFpLWhVL7i7Q0Qqn5IqDKJ20tXMBgGPAXHAk+7+cLblFlo+BDgC3Oruy0PL6gNPAsmAA+PcfXG0skvZVzUOLmoVTABbD8DCr+HDrbByB8zdGPzgQPBh2qpeMLjSOfUgsS40rPHth2yD6lArAaqGPoSznnFwh8xTwXQ089sP9MOhD/V9x4M7J/aFpj1HgwJg99GgCNhz9NsP+rxUjw/OfCSEpqqhDBmn4GRo38dPwsHj3x5TVqePr1X94K6O04VRLfXDEJEIROWMgZnFAeuAgUA68DEw2t0/z9JmCHAXQWHQC3jM3XuFlj0LLHT3J80sAajp7vvy26fOGEhWh04EnRdX7oCNe+Hr/bB5f/AAp5MF/ApUrQJxVb4tCCJVNwEa1AiKjtNTg+rB1/pZipB61aB2taB9rYSgGInEyVOhYiR09mH7weC4Th/bxr2w7VDQ1oB2DYPxIS5uDRe2VKEgUtGUtzMGPYH17r4RwMxeBoYBn2dpMwx4zoNK5UMzq29mzYDDQH/gVgB3PwFoPDwplNoJwd0NvVqcOT/jJOw4HHyw7jsKe0Ifsscyg2UZJ4O/1DNPfXsGIT70V3z1+GC7tapCnYTgw71+teADv261oJgoTXFVvi0wWhPc3pndN0eCYmjlDlixA95YBy+uCvKf3xwuaQNXtQ8ebCUiAtErDFoAW7K8Tic4K1BQmxZAJrALeMbMugLLgInufrj04kplUTUuuJSQWDfWSUpHo5pwaZtggqDfwtJtMH9z0B/j1wvhNwuDMwgjOsEV5+pMgkhlF60BjnLrF579BG5ebeKB7sBf3L0bwRmE+3PdidkdZrbUzJbu2rWrOHlFKqSEOLigJfx3X3jnZlgwFu7uBV8fgB+8A2lPwo/nwrrdsU4qIrESrcIgHWiZ5XUisC3CNulAurt/FJr/KkGhkIO7T3b3NHdPa9y4cYkEF6nIWtWHu3sHBcKr18FV5wW3fw58Acb+ExZtCTpcikjlEa3C4GOgvZm1CXUeHAXMyNZmBnCLBXoD+919u7v/B9hiZh1C7QZwZt8EESkmMzi/BTxyGSweB/f0hs92wo2vw9CX4L2NKhBEKouo9DFw90wzuxOYQ3C74tPuvtrMJoSW/xWYSXBHwnqC2xVvy7KJu4AXQ0XFxmzLRKQENawBE3vBd3vA9C/gr0th3JtBZ8X7L4S05rFOKCKlSQMciUi+Mk7C1NUw6SPYdSQYfvq+C4NRJUWk7Cip2xX1dEURyVfVOLg5BRbcCvf2gQ/TYdCL8MgiOJoR63QiUtJUGIhIRGpWhbt6wgdjYXgHeGIpXPYCvP9VrJOJSElSYSAihXJWTXj08uBR2NXj4bYZMOFt2KmRRUQqBBUGIlIkvRNh1o3woz7BWYPLX4DZ62OdSkSKS4WBiBRZQhzc2RPeGg0t6sJ334Z73w0e8CQi5ZMKAxEptvPOguk3wJ3nw2trYNA/4OOtsU4lIkWhwkBESkRCHPzoAph2HcQZjHwtGAOhgt4RLVJhqTAQkRKV1hzeHg2D2sFDi+A7b8H+Y7FOJSKRUmEgIiWuTjV4YjA8eBHM2xQMq/zZzlinEpFIqDAQkVJhBrelwivXQeYpGPEKvKqnnIiUeSoMRKRU9WgGM28Mvv7wXfj1Qjh5KtapRCQvKgxEpNQ1rAHPXQO3pMDk5XD7m7qlUaSsUmEgIlFRNQ5+eQn8+hJY+DVc8wps3hfrVCKSnQoDEYmqm1Pg+WvgmyMwbCos3x7rRCKSlQoDEYm6C1rCGyOhXjUY9RrM2RDrRCJymgoDEYmJ1vXh9Rugc2P47lvw7KexTiQioMJARGLorJrw0rUwsC38Yj785l9wSiMlisSUCgMRiakaVeGvQ2FMCvxtGdzzDmScjHUqkcorPtYBRETiqsAvL4amteDRxcEQyn8eEhQNIhJdOmMgImWCGdzVE35zaTCM8ph/wn6NdSASdSoMRKRMuakLPD4YVvwHRr0KOw/HOpFI5aLCQETKnCvPg2euhk374bppkH4g1olEKg8VBiJSJvVrBS8Oh73H4Ppp8NXeWCcSqRxUGIhImdW9Gbw8Ao6fhOtfhbXfxDqRSMWnwkBEyrSkxsGjm6sY3PAafLYj1olEKjYVBiJS5rVrCK9eD7UTYPTr8PG2WCcSqbhUGIhIuXBOPXj1OmhcE275JyxOj3UikYpJhYGIlBvN6sDU66BFHbj1DVi4OdaJRCqeqBUGZjbIzNaa2Xozuz+X5WZmfwwtX2lm3bMtjzOzT8zsrWhlFpGyp0ktmDoC2tSH29+E97+KdSKRiiUqhYGZxQFPAIOBzsBoM+ucrdlgoH1ougP4S7blE4E1pRxVRMqB0w9fan8W3PGWHtssUpKidcagJ7De3Te6+wngZWBYtjbDgOc88CFQ38yaAZhZIjAUeDJKeUWkjGtQA/5xLSQ3ge/NhLe/jHUikYohWoVBC2BLltfpoXmRtpkE/Bg4ld9OzOwOM1tqZkt37dpVrMAiUvbVqwbPXwOpTeGuWTBjbawTiZR/0SoMLJd52Z+6nmsbM7sS2OnuywraibtPdvc0d09r3LhxUXKKSDlTpxo8dw2c3xwmzoHXv4h1IpHyLVqFQTrQMsvrRCD7nch5tbkQuNrMNhFcgrjUzF4ovagiUt7USoBnhkHvFnDPHHhldawTiZRf0SoMPgbam1kbM0sARgEzsrWZAdwSujuhN7Df3be7+0/cPdHdW4fWe9/db45SbhEpJ2pWDYqDfufAj+bCPz6LdSKR8ikqhYG7ZwJ3AnMI7ix4xd1Xm9kEM5sQajYT2AisB/4OfC8a2USk4qgeD3+/Ci5tDT95H579NNaJRMofc89+qb9iSEtL86VLl8Y6hojEwPFMuHMWvLMRft4PxncveB2R8s7Mlrl7WnG3o5EPRaTCqRYPfx4CQ9rBLxfCX/Q3gkjE4mMdQESkNFSNgz8Nhvg58PAiOHESJvaKdSqRsk+FgYhUWPFVYNIVwdf/+xCOn4Qf9QHL7eZoEQFUGIhIBRdXBX5/eXB54YmP4Vhm0O9AxYFI7lQYiEiFV8XgoUuDuxae+iTonPjLS4L5InImFQYiUimYwQP9oXoc/GVZUBz89rLgjIKIfEuFgYhUGmZw34XBmYM/fARHM+EPV0BCXKyTiZQdKgxEpFIxg7t7Q42q8Jt/waET8NehwWsR0TgGIlJJfbdH0O/gg80w9g04eDzWiUTKBhUGIlJp3dgFHhsEy7bDja/DnqOxTiQSeyoMRKRSG9YBJg+Ftbvh+ldh28FYJxKJLRUGIlLpDWgLz14DOw7BiFfgyz2xTiQSOyoMRESAPokw9TrIOAXXTYPl22OdSCQ2VBiIiIQkNYbXrod61YI+B/M2xTqRSPSpMBARyaJV/aA4aNsAxr8Jr34e60Qi0aXCQEQkm8a1YOoI6NUCfvguTPoQ3GOdSiQ6VBiIiOSiTjWYMgyu6xSMkvijd4NHN4tUdBr5UEQkDwlx8OhASKwLkz6C7YfgL0OhbrVYJxMpPTpjICKSDzP4Qe+gQPhwK4yYBl/vj3UqkdKjwkBEJALXd4Znh8F/DsGwqfDR1lgnEikdKgxERCLU9xx4YyTUrwY3vQ4vr4p1IpGSp8JARKQQ2jaAf44MBkS67z343w8g81SsU4mUHBUGIiKFVK86PDMMbkuFp1bAmOnwzZFYpxIpGSoMRESKIL4KPHgR/O6y4OmMV74En/wn1qlEik+FgYhIMdyQBK/dAHFV4Ppp8MJKDYYk5ZsKAxGRYurSBN4aBRe0hJ/OC0ZLPHwi1qlEikaFgYhICWhQA565Gu7uBa+vgatehtW7Yp1KpPBUGIiIlJC4KsFgSP+4Fg6dgGumwjMrdGlBypeoFQZmNsjM1prZejO7P5flZmZ/DC1faWbdQ/Nbmtk8M1tjZqvNbGK0MouIFMUFLWH2TdDvHHjwg+ApjXuOxjqVSGSiUhiYWRzwBDAY6AyMNrPO2ZoNBtqHpjuAv4TmZwI/dPdOQG/g/+WyrohImdKwBjx1FTzQHxZ8DQOfh3c2xDqVSMGidcagJ7De3Te6+wngZWBYtjbDgOc88CFQ38yauft2d18O4O4HgTVAiyjlFhEpMjMY1w1mjIQmteE7b8E978D+47FOJpK3aBUGLYAtWV6nk/PDvcA2ZtYa6AZ8lNtOzOwOM1tqZkt37VKvHxEpGzo1DoZSntgT/vkFXP4CzN8U61QiuYtWYWC5zMveHSffNmZWG3gNuNvdD+S2E3ef7O5p7p7WuHHjIocVESlpCXFwT59gOOXaCTD2Dfj+bI2YKGVPtAqDdKBllteJwLZI25hZVYKi4EV3f70Uc4qIlKqUpvD26ODswcwv4dLn4KVVcEp3LkgZEa3C4GOgvZm1MbMEYBQwI1ubGcAtobsTegP73X27mRnwFLDG3f8vSnlFREpN9fjg7MHsm6BjI7j/PbjhVfjim1gnE4lSYeDumcCdwByCzoOvuPtqM5tgZhNCzWYCG4H1wN+B74XmXwiMAS41sxWhaUg0couIlKZ2DWHqiOB5C1/ugcH/gJ++D7t1eUFiyLyCjryRlpbmS5cujXUMEZGI7DsGf/gQnl8JtarCxF5wS9egb4JIJMxsmbunFXc7GvlQRKQMqF8d/ufi4PJCt2bwy4XB2Acz1qr/gUSXCgMRkTLkvLPg2WHw9NVQLR7umg1D/wHvfaWhlSU6VBiIiJQxZjCgDcy6ER67Ag5nwLgZMGIaLNysAkFKlwoDEZEyKq4KXNMR3hsDD10K2w7Czf8Mntw480tdYpDSocJARKSMqxoHN3aBD8bCbwfAwePwXzPhsudh6mo4lhnrhFKRqDAQESknqsXDqGR4/xZ4fDBUi4Mfz4ULnobf/Rv+cyjWCaUi0O2KIiLllDv8Ox2mrIB3NwaXHga3g5u6QO8WQV8FqTxK6nbF+JIIIyIi0WcGF7YMpq/3w3Mrg0sLb66DVvXghs5wXWc4u3ask0p5ojMGIiIVyNEMmLUhKBA+TIcqBv3PgWEd4fK2wQOcpGLSGQMREcmhRlW4tmMwbd4Hr3wO07+AH8wJ+iQMaANXd4BLWgfPbBDJTmcMREQquFMOy7YHoyi+/SXsPgo14qF/q+AswoA20KBGrFNKcemMgYiIRKSKwfnNg+mBi2DxFpizMeiwOGcDxBmkNYeLWgVT58bBOlI56YyBiEgldcrhsx3wzkZ4/yv4PPTY50Y1oF8r6NsSeidCYt3Y5pTI6IyBiIgUSxWDrmcH048ugB2HgyGXPwhN078I2iXWhT4toGcL6NEc2tbXrZAVWYU9Y2BmUT+w0/+WZoa7Y1H8zcltf1nzFHUbkSzLvr/89pm9zentZs+afX+5Lc9tP7mtl1uWrPvNb5+5vc7tOLNvL5J/y7yOPZKsBS0rKE9u+84tR2HWPz0/r3+b/OT3fmbNl9/PSiTHnP1Ys8rveHI7voK2lduxZ/95y01uP2OR5ChMvvzynnbKYd1uWJwe3N0w+SrjnEnB8q/vNsb+0+nRDFKbQpemwdMhS0qkuUt6O4VpH8nPUWEU9HMR4TZ0xkBEREpHFYOOjYLptlSYDLx7MyzfDqOBLQdg3qZv27esCylNIbkxdGoMnRpB01o6s1AeqTAQEZGInHdWMI0meLDT/mPw2c5vp5U7grseTmtQPSgsOpwF7RtC+9DXhroDokxTYSAiIkVSrzr0PSeYTtt/HNZ+E3Rk/OIbWLMLpn0ePDr6tIY1gn4KbRp8+7VVvWCqpQGYYk6FgYiIlJh61YJOij1bfDvPHbYfgi/3wJe7Yf1e+Gpv0MFx2udnrt8odDZh4uyg02NiXUisE3xtVkeDMkWD/olFRKRUmUHzOsF0Uaszlx06AV/tC571sHkfbN4Py4Gl24JnPpzM1hevUY2gQGheB5rVDqamtYPnQTSrDU1qQc2qUTqwCkqFgYiIxEztBOjSJJhOewRYNA4yTwWPkk4/EExbDwZnHrYdhI17YdGWoLDIrk5CUCAAfH82NK4JjWuFvtaERjXhrJrQsDpUjYvKYZYrKgxERKRMiq/y7eWEvBw6ERQPp6edh4PxGHYehnnAJ/+BXYfhaGbu6zeoHvR5OCt0CeO/34OGNYPXp5c1CH3foHpwNqKi32mhwkBERMqt2gnQrmEwZfdXYOGtwfeHT8CuI0GR8M1R+OYI7D4SzNt9FPYcDdrN3gB7jwXjOOQmIe7bMRtueDX4vl614Ovp+W+uC+bVqxZ00KxbLTiLUV7OTqgwEBGRCq9WQjC1rp93GwOW3wEnT8GB47DnWFAw7A0VDvuOBUXD3mOwNNR+8z7YdzxYdix0VuLOWXlkqBoqEqpB3YTg+9OvT3t+ZVBE1AkVE3USoHbo+1pVo1NcqDAQERHJIq5K6PJBDTi3Qe5tfgdMve7MeccyocbdwUBQ+48Ft24eOJ7t6zE4cCJ4vfNIcIfGwePfbuNn8/LPVj3+2yKhVkJwxqRW1eBrSVFhICIiUgJO30p53lmFX9cmBF8/Hg8HT8Ch48HXgyeCwuHQCTiUEfoaWnY4I7hEsvNIcGdHSVFhICIiUkY0qfXtHRWFZbeWTIYqJbOZgpnZIDNba2brzez+XJabmf0xtHylmXWPdF0REREpGVEpDMwsDngCGAx0BkabWedszQYD7UPTHcBfCrGuiIiIlIBonTHoCax3943ufgJ4GRiWrc0w4DkPfAjUN7NmEa4rIiIiJSBahUELYEuW1+mheZG0iWRdERERKQHR6nyY2zhR2YePyKtNJOsGGzC7g+AyBMBxYFWkAUuCZRkOy0p3aKxGwDd57buoGfJrH8m2itLm9Ots8xuZWb7Hl9e+Ctsuv/cswqy5Livg37IR8E0kbQu730jzRHL8hV0/y+szfj6L87OTW75IjzW/97OoPz+FPb6iZChEjgLXK2g/ubTN8X9LXtsq5f/jSmz72f9vIZ/jK+x+C/M+FHebEepQnJVPi1ZhkA60zPI6EdgWYZuECNYFwN0nA5MBzGypu6cVL3bZVJGPDXR85Z2Or/yqyMcGleP4SmI70bqU8DHQ3szamFkCMAqYka3NDOCW0N0JvYH97r49wnVFRESkBETljIG7Z5rZncAcIA542t1XmwVDOrj7X4GZwBBgPXAEuC2/daORW0REpLKJ2gBH7j6T4MM/67y/Zvnegf8X6boRmFzYjOVIRT420PGVdzq+8qsiHxvo+CJiweexiIiISBRHPhQREZGyr1wXBsUZZrmsM7OWZjbPzNaY2Wozm5hLm4vNbL+ZrQhNv4hF1qIys01m9lkoe47etOX8/euQ5X1ZYWYHzOzubG3K1ftnZk+b2U4zW5VlXkMze9fMvgx9zfVZdAX9rpYFeRzf78zsi9DP33Qzq5/Huvn+LMdaHsf2oJltzfLzNySPdcvrezc1y7FtMrMVeaxbpt87yPvzoNR+/9y9XE4EHRE3AG0Jbmn8FOicrc0QYBbBWAi9gY9inbsQx9cM6B76vg6wLpfjuxh4K9ZZi3GMm4BG+Swvt+9ftuOIA/4DtCrP7x/QH+gOrMoy7xHg/tD39wO/zeP48/1dLQtTHsd3ORAf+v63uR1faFm+P8uxnvI4tgeBewtYr9y+d9mW/x74RXl870IZc/08KK3fv/J8xqA4wyyXee6+3d2Xh74/CKyh8o34WG7fv2wGABvcfXOsgxSHuy8A9mSbPQx4NvT9s8A1uaxaLoY1z+343P0dd88MvfyQYByVcieP9y4S5fa9O83MDLgBeCmqoUpQPp8HpfL7V54Lg+IMs1yumFlroBvwUS6L+5jZp2Y2y8ySopus2Bx4x8yWWTBqZXYV4v0jGHsjr/+UyvP7B9DUg/FGCH1tkkubivI+jiM4g5Wbgn6Wy6o7Q5dJns7jNHRFeO/6ATvc/cs8lper9y7b50Gp/P6V58KgOMMslxtmVht4Dbjb3Q9kW7yc4PR0V+BPwD+jHK+4LnT37gRPzvx/ZtY/2/KK8P4lAFcD03JZXN7fv0hVhPfxp0Am8GIeTQr6WS6L/gKcC6QC2wlOt2dX7t87YDT5ny0oN+9dAZ8Hea6Wy7x838PyXBgUZ5jlcsHMqhL8ELzo7q9nX+7uB9z9UOj7mUBVC8Y6LxfcfVvo605gOsEpr6zK9fsXMhhY7u47si8o7+9fyI7Tl3dCX3fm0qZcv49mNha4ErjJQxdts4vgZ7nMcfcd7n7S3U8Bfyf3zOX9vYsHrgWm5tWmvLx3eXwelMrvX3kuDIozzHKZF7ou9hSwxt3/L482Z4faYWY9Cd7P3dFLWXRmVsvM6pz+nqCTV/aHXpXb9y+LPP9aKc/vXxYzgLGh78cCb+TSptwOa25mg4D7gKvd/UgebSL5WS5zsvXXGU7umcvtexdyGfCFu6fntrC8vHf5fB6Uzu9frHtbFmci6LW+jqDH5U9D8yYAE0LfG/BEaPlnQFqsMxfi2PoSnO5ZCawITUOyHd+dwGqCXqYfAhfEOnchjq9tKPenoWOoUO9fKH9Ngg/6elnmldv3j6DA2Q5kEPwVcjtwFvAe8GXoa8NQ2+bAzCzr5vhdLWtTHse3nuD67Onfwb9mP768fpbL0pTHsT0f+r1aSfBB0awivXeh+VNO/75laVuu3rtQzrw+D0rl908jH4qIiEhYeb6UICIiIiVMhYGIiIiEqTAQERGRMBUGIiIiEqbCQERERMJUGIiIiEiYCgMREREJU2EgIiUi9Lz4gaHvf2Vmf4x1JhEpvPhYBxCRCuMB4H/NrAnB09+ujnEeESkCjXwoIiXGzD4AagMXe/DceBEpZ3QpQURKhJl1AZoBx1UUiJRfKgxEpNhCT+p7ERgGHDazK2IcSUSKSIWBiBSLmdUEXgd+6O5rgF8CD8Y0lIgUmfoYiIiISJjOGIiIiEiYCgMREREJU2EgIiIiYSoMREREJEyFgYiIiISpMBAREZEwFQYiIiISpsJAREREwv4/wF3pVjPquf8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot fitted pdf (MAP estimators)\n",
    "yMin = 0.\n",
    "yMax = f(0., MAP)*1.1\n",
    "fig = plt.figure(figsize=(8,6))\n",
    "xCurve = np.linspace(xMin, xMax, 100)\n",
    "yCurve = f(xCurve, MAP)\n",
    "plt.plot(xCurve, yCurve, color='dodgerblue')\n",
    "plt.xlabel(r'$x$')\n",
    "plt.ylabel(r'$f(x)$')\n",
    "y_fitval = 0.8\n",
    "delta_y_fitval = 0.08\n",
    "plt.figtext(0.6, y_fitval, 'MAP Estimators')\n",
    "for i in range(len(parin)):\n",
    "    y_fitval -= delta_y_fitval\n",
    "    if not parfix[i]:\n",
    "        plt.figtext(0.6, y_fitval, parname_latex[i] + ' = ' + f'{MAP[i]:.4f}' + r'$\\pm$' + f'{sigmaMAP[i]:.4f}')\n",
    "\n",
    "# Add data as tick marks\n",
    "tick_height = 0.05*(yMax - yMin)\n",
    "xvals = [xData, xData]\n",
    "yvals = [np.zeros_like(xData), tick_height * np.ones_like(xData)]\n",
    "plt.plot(xvals, yvals, color='black', linewidth=1)\n",
    "plt.xlim(xMin, xMax)\n",
    "plt.ylim(yMin, yMax)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "ed1a0725",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "posterior covariance = \n",
      "┌───────┬─────────────────────────────────┐\n",
      "│       │   theta      mu   sigma      xi │\n",
      "├───────┼─────────────────────────────────┤\n",
      "│ theta │ 0.00337 -0.0118  0.0197 -0.0208 │\n",
      "│    mu │ -0.0118   0.308  -0.114  0.0137 │\n",
      "│ sigma │  0.0197  -0.114    0.26  -0.124 │\n",
      "│    xi │ -0.0208  0.0137  -0.124   0.285 │\n",
      "└───────┴─────────────────────────────────┘\n",
      "proposal covariance = \n",
      "┌───────┬─────────────────────────────────┐\n",
      "│       │   theta      mu   sigma      xi │\n",
      "├───────┼─────────────────────────────────┤\n",
      "│ theta │ 0.00478 -0.0167  0.0279 -0.0294 │\n",
      "│    mu │ -0.0167   0.437  -0.162  0.0193 │\n",
      "│ sigma │  0.0279  -0.162   0.369  -0.176 │\n",
      "│    xi │ -0.0294  0.0193  -0.176   0.404 │\n",
      "└───────┴─────────────────────────────────┘\n"
     ]
    }
   ],
   "source": [
    "# Set up MCMC\n",
    "# see J. Rosenthal, http://probability.ca/jeff/ftpdir/galinart.pdf\n",
    "scale_prop = 2.38**2/numFreePar      # or adjust to minimize acf\n",
    "cov_prop = scale_prop*cov\n",
    "print('posterior covariance = ')\n",
    "print(cov)\n",
    "print('proposal covariance = ')\n",
    "print(cov_prop)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "9b2a0a74",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Start MCMC iterations: ...............................................................................................\n",
      " MCMC acceptance fraction =  0.3737\n"
     ]
    }
   ],
   "source": [
    "# Iterate with Metropolis-Hastings algorithm\n",
    "chain = [np.array(MAP)]         # start point is MAP estimate\n",
    "numIterate = 10000\n",
    "numBurn = 100\n",
    "numAccept = 0\n",
    "print(\"Start MCMC iterations: \", end=\"\")\n",
    "while len(chain) < numIterate:\n",
    "    par = chain[-1]\n",
    "    log_post = -negLogL(par) + np.log(prior(par))\n",
    "    par_prop = np.random.multivariate_normal(par, cov_prop)\n",
    "    if prior(par_prop) <= 0:\n",
    "        chain.append(chain[-1])    # never accept if prob<=0.\n",
    "    else:\n",
    "        log_post_prop = -negLogL(par_prop) + np.log(prior(par_prop))\n",
    "        alpha = np.exp(log_post_prop - log_post)\n",
    "        u = np.random.uniform(0, 1)\n",
    "        if u <= alpha:\n",
    "            chain.append(par_prop)\n",
    "            numAccept += 1\n",
    "        else:\n",
    "            chain.append(chain[-1])\n",
    "        if len(chain)%(numIterate/100) == 0:\n",
    "            print(\".\", end=\"\", flush=True)\n",
    "chain = np.array(chain)\n",
    "print('\\n MCMC acceptance fraction = ', numAccept/numIterate)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b66bba58",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Show trace plots\n",
    "fig, axes = plt.subplots(sharex=True, nrows=4, figsize=(8,6))\n",
    "plt.xlim(0, numIterate)\n",
    "for i, (ax, label) in enumerate(zip(axes, parname_latex)):\n",
    "    ax.plot(chain[:, i], color=\"dodgerblue\")\n",
    "    ax.axvline(numBurn, color=\"orange\")\n",
    "    ax.set_ylabel(label)\n",
    "axes[-1].set_xlabel(\"iteration number\")\n",
    "fig.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bc4f8db3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Show marginal distributions\n",
    "fig, axes = plt.subplots(nrows=4, figsize=(8,6))\n",
    "for i, (ax, label, sim) in enumerate(zip(axes, parname_latex, np.array(MAP))):\n",
    "    ax.hist(chain[numBurn:, i], bins=50, density=True, color=\"dodgerblue\")\n",
    "    ax.set(xlabel=label, ylabel=\"pdf\")\n",
    "    ax.axvline(sim, color=\"orange\")\n",
    "fig.tight_layout(pad=0.3)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b79529e7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find autocorrelation and plot\n",
    "var = chain[numBurn:,].var(axis=0)\n",
    "x   = chain[numBurn:,] - chain[numBurn:,].mean(axis=0)\n",
    "max_lag = 200\n",
    "N   = len(x[:,0])\n",
    "acf = np.zeros(shape=[max_lag, len(chain[0])])\n",
    "for i in range(len(chain[0])):\n",
    "    if var[i] > 0:\n",
    "        acf[:,i] = np.correlate(x[:,i], x[:,i], mode='full')[N-1:N-1+max_lag]\n",
    "        acf[:,i] /= (var[i]*N)\n",
    "    \n",
    "fig, axes = plt.subplots(sharex=True, nrows=4, figsize=(8,6))\n",
    "for i, (ax, label) in enumerate(zip(axes, parname_latex)):\n",
    "    ax.plot(acf[:,i], color=\"dodgerblue\")\n",
    "    ax.set_ylabel('ACF[' + label + ']')\n",
    "    ax.set_yticks([0., 0.5, 1.])\n",
    "    ax.grid()\n",
    "axes[-1].set_xlabel(\"lag\")\n",
    "plt.xlim(0, max_lag)\n",
    "plt.subplots_adjust(left=0.15, bottom=0.1, right=0.9, top=0.9, wspace=0.4, hspace=0.4)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1fc33d96",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Functions for computing intervals\n",
    "# input:  values = numpy array with MCMC chain, e.g., chain[:,0]\n",
    "#         CL = credibility level, e.g., 0.683\n",
    "\n",
    "# Central credible interval\n",
    "def cc_interval(values, CL):\n",
    "    alpha = np.array([(1.-CL)/2., (1.+CL)/2.])\n",
    "    return np.quantile(values, alpha)\n",
    "\n",
    "# Highest Probability Density (HPD) interval from MCMC chain\n",
    "def HPD_interval(values, CL):\n",
    "    sorted = np.sort(np.copy(values))\n",
    "    nPoints = len(values)\n",
    "    nCovered = np.floor(0.683 * nPoints).astype(int)\n",
    "    intWidth = sorted[nCovered:] - sorted[:nPoints-nCovered]\n",
    "    minWidth = np.argmin(intWidth)\n",
    "    hpdLo = sorted[minWidth]\n",
    "    hpdHi = sorted[minWidth+nCovered]\n",
    "    return np.array([hpdLo, hpdHi])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fe79498a",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}