{ "cells": [ { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "K7l4AAzxVr-K" }, "source": [ "\n", "\n", "# Parameter estimation for the Cauchy (Lorentzian) distribution" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "Nx6NrJGKVr-Q" }, "source": [ "## Introduction\n", "**Cauchy (Lorentzian) distribution** is a symmetric distribution described by location parameter $\\mu$ and a scale parameter $\\gamma$ as \n", "\n", "$$p(x|\\mu, \\gamma) = \\frac{1}{\\pi \\gamma} (\\frac{\\gamma^2}{\\gamma^2+(x-\\mu)^2})$$ \n", "\n", "Because its tails decrease as slowly as $x^{-2}$ for large $|x|$, the mean, variance, standard deviation, and higher moments do not exist. Therefore we cannot estimate location and scale parameter from the mean, standard deviation. Instead, we can estimate it from median value and interquartile range for {$x_i$} using a Bayesian approach. \n", " \n", "Suppose we want to determine the location of a lighthouse signal along the coastline. Let coastline located at y = 0, and the lighthouse is at distance $\\gamma$ away from the coastline. Define the angle between the line of light and coastline as $\\theta$, then the position is expressed as $x = \\mu + \\gamma tan(\\theta)$. From $-\\pi/2 \\leq \\theta \\leq \\pi/2$, and angle $\\theta$ is distributed uniformly, we can find data likelihood using $p(x) = (\\pi dx/d\\theta)^{-1}p(\\theta)$. \n", "The datalikelihood for a set of data {$x_i$} with Cauchy distribution is\n", "\n", "$$p({x_i}|\\mu,\\gamma,I) = \\prod_{i=1}^N \\frac{1}{\\pi}(\\frac{\\gamma}{\\gamma^2+(x_i-\\mu)^2})$$ \n", " \n", "In this notebook, we will explore how the changes of $\\mu$ and $\\gamma$ affect the estimated probability density function (pdf)." ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "DDikty8OVr-T" }, "source": [ "## Import Functions\n", "The core function used for the Cauchy function is cauchy in scipy.stats." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": {}, "colab_type": "code", "id": "76toF_JMVr-W" }, "outputs": [], "source": [ "import numpy as np\n", "from matplotlib import pyplot as plt\n", "from scipy.stats import cauchy\n", "from astroML.plotting.mcmc import convert_to_stdev\n", "from astroML.stats import median_sigmaG\n", "from astroML.resample import bootstrap" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "z2L8rf7rVr-l" }, "source": [ "## Log-likelihood for Cauchy Distribution\n", "Given a data set of measured positions {$x_i$}, we need to estimate $\\mu$ and $\\gamma$. Analogously to the Gaussian case discussed in the prevoius chapter, we shall adopt a uniform prior distribution for the location parameter $\\mu$, and a uniform prior distribution for ln$\\gamma$, for $\\mu_{min} \\leq \\mu \\leq \\mu_{max}$ and $\\gamma_{min} \\leq \\gamma \\leq \\gamma_{max}$. \n", " \n", "The logarithm of the posterior pdf is\n", "\n", "$$L_p \\equiv ln[p(\\mu,\\gamma|{x_i},I)] = constant + (N-1)ln\\gamma - \\sum^N_{i=1}ln[\\gamma^2 + (x_i-\\mu)^2]$$" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "sIJ7xN0uVr-o" }, "source": [ "### 1. Define function\n", "We define the log-likelihood function as \n", "\n", "$$L_p \\equiv ln[p(\\mu,\\gamma|{x_i},I)] = constant + (N-1)ln\\gamma - \\sum^N_{i=1}ln[\\gamma^2 + (x_i-\\mu)^2]$$\n", "\n", "In this example, we will use N = 10 values of $x_i$, $\\mu$ = 0, and $\\gamma$ = 2." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": {}, "colab_type": "code", "id": "37Q0eG7BVr-q" }, "outputs": [], "source": [ "def cauchy_logL(xi, gamma, mu):\n", " \"\"\"Equation 5.74: cauchy likelihood\"\"\"\n", " xi = np.asarray(xi)\n", " n = xi.size\n", " shape = np.broadcast(gamma, mu).shape\n", "\n", " xi = xi.reshape(xi.shape + tuple([1 for s in shape]))\n", "\n", " return ((n - 1) * np.log(gamma)\n", " - np.sum(np.log(gamma ** 2 + (xi - mu) ** 2), 0))\n", "\n", "# Define the grid and compute logL\n", "gamma = np.linspace(0.1, 5, 70)\n", "mu = np.linspace(-5, 5, 70)\n", "\n", "np.random.seed(44)\n", "mu0 = 0\n", "gamma0 = 2\n", "xi = cauchy(mu0, gamma0).rvs(10)\n", "\n", "logL = cauchy_logL(xi, gamma[:, np.newaxis], mu)\n", "logL -= logL.max()" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "XjiRQRinVr-0" }, "source": [ "### 2. Find mu and gamma at max data likelihood\n", "We can find $\\mu$ and $\\gamma$ that maximize L by setting the derivative of L to 0, as $\\frac{dlnL(\\mu)}{d\\mu}\\Big|_{\\mu^0} \\equiv 0$. Here in this example, the result is $\\mu = -0.36$, and $\\gamma = 0.81$. \n", " \n", "Median and interquartile range imply $\\mu = -0.26$ and $\\gamma = 1.11$. The interquartile range for the Cauchy distribution is equal to $2\\gamma$ and thus\n", "$$\\sigma_G = 1.483\\gamma$$\n", "\n", "When using this shortcut, the bootstrap method can be used to estimate parameter uncertainties as shown in the next section." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": {}, "colab_type": "code", "id": "2gh0hD7iVr-2", "outputId": "550753cd-33ea-4b7a-e205-07388840b9c6" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "mu from likelihood: [-0.36231884]\n", "gamma from likelihood: [0.81014493]\n", "\n", "mu from median -0.2625695623133895\n", "gamma from quartiles: 1.10950042396172\n", "\n" ] } ], "source": [ "i, j = np.where(logL >= np.max(logL))\n", "\n", "print(\"mu from likelihood:\", mu[j])\n", "print(\"gamma from likelihood:\", gamma[i])\n", "print()\n", "\n", "med, sigG = median_sigmaG(xi)\n", "print(\"mu from median\", med)\n", "print(\"gamma from quartiles:\", sigG / 1.483)\n", "print()" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "CWlWJI5EVr_D" }, "source": [ "### 3. Plot the results\n", "Here we show the result of the logarithm of posterior probability distribution for $\\mu$ and $\\gamma$ for $L_p(\\mu,\\gamma)$ given above. \n", "This example uses N = 10, $\\mu$ = 0 and $\\gamma$ = 2. The maximum of L is renormalized to 0, and color coded as shown in the legend. The contours enclose the regions that contain 0.683, 0.955 and 0.997 of the cumulative (integrated) posterior probability." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": {}, "colab_type": "code", "id": "AZ1gtXlkVr_G", "outputId": "73b6142d-956b-47bd-d383-d751ccfb8094" }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUoAAAD/CAYAAACXWZ93AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO2deViUVfvHPwcQUUBEQVxQMRdAUVxRNHctW9TU3LLF8s02e3vbfPu9pWVpZqXlUpltVppauee+pyaK+5a44IY7igjIPuf3BwwNA8MMw+ycz3U9l7M8zzn3jPN8ue9z7nMfIaVEoVAoFIZxs7cBCoVC4egooVQoFAojKKFUKBQKIyihVCgUCiMooVQoFAojKKFUKBQKI9hEKIUQ54QQR4QQB4UQe23Rp0KhcE6EEH2EEHFCiNNCiLfsbQ+AsEUepRDiHNBWSplo9c4UCoXTIoRwB04CvYEEIBYYLqU8bk+7VOitUCgciSjgtJQyXkqZBSwE+tvZJpsJpQTWCyH2CSFG26hPhULhfNQBLuo8T8h/za542KifTlLKy0KIGsAGIcQJKeWfuifkC+hoAG9v7zZhYWE2Mk2hUGg5d+4ciYmJorTXCSFMHcM7BmToPJ8jpZyj21Qx19h9nbVNhFJKeTn/3+tCiKXkudd/6p0zB5gD0LZtW7l3r5rzUShsTdu2bc2+Vgjj+iqlzJBSltRJAlBX53kwcNlsoyyE1UNvIYS3EMJX+xi4Dzhq7X4VCoVtEUIYPUwgFmgshGgghPAEhgErrGq4CdjCowwCluZ/SR7AL1LKtTbo12yioqJISEiwtxmKckZwcDB79uyxtxlmY6IQloiUMkcIMQZYB7gD30spj5W54TJidaGUUsYDkdbux5IkJCRw+bLdvX1FOaN27dr2NsFshBC4uRkPUHNzc42eI6VcDay2gFkWw1aTOQqFwsWxhEfpqCihVCgUFkEJpUKhUBhBCaXCJencuTMpKSl4eHig0rEUZcHUMUpnxXU/mRX4+uuveemllwq9lp6eTteuXU0apDaFrKwsunTpQk5OjsFzZsyYQXh4OCNGjChTX9u3b+fgwYNlEsm1a9cSGhpKo0aN+Oijj8pkjzlcvHiR7t27Ex4eTrNmzZg+fbrNbQB45plnqFGjBhEREYVet/f3Y0sslB7kkCihLAWHDx+mefPmhV77/vvvGThwIO7u7hbpw9PTk549e7Jo0SKD53z55ZesXr2a+fPnG21PSolGo7GIbfrk5uby0ksvsWbNGo4fP86CBQs4fty2tQs8PDyYOnUqf//9NzExMXzxxRc2twFg5MiRrF1bOOvNEb4fW6KEUgHAkSNHigjl/Pnz6d//nzX73bp1Iy4uDoCbN28W8TC0nDp1ipCQEE6fPg1AdnY2kZGRJCQk8MgjjxgUweeff574+Hj69evHZ599xrRp04iIiCAiIoLPP/8cyFuGFh4ezosvvkjr1q25ePFikXa6d+/Ohg0bAHjnnXf497//XcpvA/bs2UOjRo2455578PT0ZNiwYSxfvtzodSV99tJSq1YtWrduDYCvry/h4eFcunTJpGsHDBjAO++8Q+fOnalZsyYbN24sdf9aunTpQrVq1Qq9Zu7346y4slCqMcpScPTo0ULCl5WVRXx8PCEhIQWvnT59msaNGwPFe6BaGjduzOjRo1m3bh2NGjVi1qxZ9O/fn+DgYGrVqkVsbGyx182ePZu1a9eyZcsWzp8/z8iRI9m9ezdSStq3b0/Xrl3x9/cnLi6OH374gS+//LLYdiZMmMD48eO5fv06Bw4cYMWKfxY/aMcu9fn000/p1atXwfNLly5Rt+4/q82Cg4PZvXt3sf2Z+tlLa4Mu586d48CBA7Rv396oDZD3/9mpUye2b9/OkiVLmD9/fpG2zbFDi7nfjzPi6mOUSihN5OLFi/j6+uLn51fwWmJiIlWrVi14fv78eerUqVPwgzl8+DAtWrQw2GZERAQbN27k1q1bfPfddwU3kbu7O56enqSkpODr62vw+h07djBgwAC8vb0BGDhwINu3b6dfv37Ur1+fDh06GLy2S5cuSCmZNm0aW7duLTR0sH37diPfRh7F1TI11Wsw9NlLa4OW1NRUBg0axOeff06VKlWMnn/37l2Sk5N59dVXAcjJySn0f2muHbqU5ftxRlz5symhNJHivMNKlSqRkfFPIZSDBw8WEsZ9+/YxdOhQg202adKEL774gvfee4833nijQPAAMjMz8fLyKtGmkoou67ZVHEeOHOHKlSsEBAQUEWNTvajg4OBCYX1CQoLJq0tK+uylsQHyQvdBgwYxYsQIBg4caFL/x44do02bNgV/IA4fPlzsMElZPMqyfD/OiCsLJVJKhzvatGkj7UmtWrWKvDZ58mT51ltvFXk9ODhYpqenSymlnDBhgnz88cellFKePHlSVqlSRZ4/f15KKWWPHj1kQkJCoWuzsrJkQECAbN++vczNzS14PTExUYaFhRm0r379+vLGjRty3759snnz5jItLU2mpqbKZs2ayf3798uzZ8/KZs2aGbz+8uXLsnnz5vL48eOyV69ecu3atSV8G4bJzs6WDRo0kPHx8TIzM1O2aNFCHj161ODn1cXQZy8tGo1GPvHEE/KVV14p8l5JNnz//feF/j/79u0r9+zZY7YdUsoi33tJ309xFPe7szX5916p71l3d3dZrVo1owew15z27X247qCChTly5Ahz5swhJCSEkJAQoqOjAbjvvvvYsWMHkOdRajQaIiMjef/99wkPD+fHH39Eo9Fw+vTpIoP9FSpUoEqVKnz00UeFxne2bNnCgw8+aNSm1q1bM3LkSKKiomjfvj3/+te/aNWqVYnX3L17l4EDBzJ16lTCw8MZN24c7733Xim/jTw8PDyYNWsW999/P+Hh4QwZMoRmzZoZ/Ly6GPrspWXnzp38/PPPbN68mZYtW9KyZUtWr15t1IYjR47QsmXLguf648+lZfjw4URHRxMXF0dwcDDfffedwe/HFdGOURo7nBZ7K3VxhyN6lIbYv39/gRfZsGFDeefOnSLnHDlyRL766qvFXl+3bl2p0WgKvTZgwAB54sSJUljsWJT0eXUp7rPb2gZHwpk9Sg8PDxkYGGj0wEk9SjVGWUZatWpF9+7dSU5Oxs3NrdjJl4iICKZNm1bk9XPnzlG/fv1CYztZWVk88sgjhIaGWtVua2Lo8+pS3Ge3tQ0Ky+LKY5Q22YWxtNi7wnnt2rVVmTWFzXGE313btm3Zu3dvqRWvQoUKsnr16kbPu3bt2j5ZcoVzh0R5lAqFosyoPEqFQqEwAVcOvZVQKhQKi6CEUqFQKIyghFKhUChKQI1RlkOCg4NdeqmZwjHRLQrijCiPspzhzFuGKhT2QgmlQlECiYmJxMbGcuHCBS5dusTly5dxd3fH398ff39/goKCaNCgASEhIdSpUwcPD/Wzc0WUUCoUOty8eZNNmzaxYcMGduzYwYkTJwrec3NzIygoiNzcXJKSksjOzi50rY+PDw8++CADBgzg/vvvx9/f39bmK6yAGqNUKICzZ8/y/fffs2LFCg4fPgxA1apV6dSpE0899RTR0dE0bNiQmjVrFniMUkrS09O5cuUKZ8+e5dy5c8TGxrJs2TJ+/fVXhBA0b96cTp06ce+999K3b98S628qHBtX9ijVEkaFQaSULF++nC+//JINGzbg5uZG165d6dmzJ927dycqKsqsMDo3N5eYmBg2btzIzp072bVrF6mpqQQGBjJ+/HhGjx6Np6enFT6RwhjmLmH08vKS9erVM3reqVOn1BJGhesQExPD66+/zl9//UW9evWYMGECTz/9dKGtDczF3d2dTp060alTJyCvuviuXbsYP348L7/8Mp9//jkTJ05k8ODBFtu0TWF9XNmjdN1BBYVZ3Lhxg8cee4zo6Gji4+P55ptviI+PZ/z48RYRyeLw8PCgc+fObN68mVWrVlGpUiWGDx9OWFgY33zzDZmZmVbpV2E5XL0epfNarrA4a9asoXnz5ixZsoRx48Zx6tQp/vWvf9nMqxNC8OCDD3Lo0CF+//13/Pz8GD16NKGhoS67KZcr4cq7MCqhVJCens6///1vHnzwQQIDA4mNjeX999/Hx8fHLva4ubkxaNAgYmNjWb9+PUIIOnfuzIwZM3DEMXVFHkooFS7LyZMn6dChAzNnzuQ///kPsbGxBrfYtTVCCHr37s3+/fvp06cPr7zyCsOHDyctLc3epimKQYXeCpdk4cKFtGnThkuXLrFq1So+++wzozs/2gN/f3+WLVvG5MmT+e233+jWrRtXr161t1kKHUzxJpVHqXAqpJS88cYbDB8+nBYtWnDgwAGTNjOzJ25ubrz11lssW7aM48eP0759e44ePWpvsxQ6WFsohRCDhRDHhBAaIYRNU4xsJpRCCHchxAEhxB+26lNRlNzcXEaPHs3UqVMZM2YMW7dutdpstjXo27cvf/75J9nZ2XTq1Ilt27bZ2yRFPjbwKI8CA4E/y25t6bClR/kK8LcN+1PokZOTw5NPPsm3337LO++8w4wZM6hQoYK9zSo1bdq0ISYmhjp16tCnTx/WrFljb5MUWH+MUkr5t5QyzkLmlgqbCKUQIhh4CPjWFv0pipKTk8OwYcP45ZdfmDx5Mh988IFTjxnVq1ePbdu2ER4eTv/+/fn999/tbVK5Ro1RWobPgbGAxkb9KXSQUjJmzBgWL17MtGnTeOutt+xtkkUIDAxk8+bNtGvXjqFDh7Jw4UJ7m1SuMVEoA4QQe3WO0XptbBRCHC3m6G+fT5WH1ZcwCiEeBq5LKfcJIbqVcN5oYDTkeQsKyzFp0iS+/vpr3nrrLV599VV7m2NRqlatyrp163j44Yd5/PHHqVSpEv372/WeKreY6DEmlrTWW0rZy3IWWQ5beJSdgH5CiHPAQqCHEGKe/klSyjlSyrZSyraBgYE2MKt88MMPPzBu3DieeOIJPvzwQ3ubYxV8fHxYuXIlbdq0YciQIWzatMneJpVLVB5lGZBS/p+UMlhKGQIMAzZLKR+3dr8KWL16Nc8++yy9e/fm22+/deoxImP4+vqyZs0aQkNDGTBgAEeOHLG3SeUKW4xRCiEGCCESgGhglRBinUWMNwHnlXhFiezevZvBgwcTGRnJ4sWLy0XZsmrVqrF69Wp8fX15+OGHSUpKsrdJ5QprC6WUcmm+01VRShkkpbzfQqYbxaZCKaXcKqV82JZ9lkdOnTrFQw89RM2aNQuEo7wQHBzM0qVLuXz5Mv/617/U2nAboma9FU5DRkYGgwcPBmDdunUEBQXZ2SLbExUVxYcffsiSJUv48ccf7W1OuUGNUSqchokTJ3Lo0CHmzp1Lo0aN7G2O3Xj99dfp2rUr//73vzl37py9zXF5VB6lwmk4duwYU6ZM4cknn+Thh8v3CIebmxtz584FYOTIkeTm5trXoHKAEkqFw6NNKq9SpQpTp061tzkOQUhICDNmzGDbtm1MmTLF3ua4PEooFQ7P77//ztatW5k0aRIBAQH2NsdheOqppxg+fDjjx49n586d9jbHpVFjlAqH5u7du7z++utERkby7LPP2tsch0IIwezZs6lfvz7Dhw/n9u3b9jbJJVFjlAqH54MPPuDixYvMnDlT7VpYDFWqVGHBggUkJCTw3nvv2dscl0UJpcJhOXr0KJ9++ilPP/00nTt3trc5DktUVBTPPvsss2bN4tixY/Y2xyVRQqlwSKSUvPjii/j5+fHxxx/b2xyHZ+LEifj6+vLKK6+oRHQroMYoFQ7JH3/8wfbt2/nwww/VBI4JBAYGMnHiRDZt2qTqV1oYNUapcEg0Gg3jxo2jYcOGPP300/Y2x2l4/vnnadmyJa+99hqpqan2NselUEKpcDiWLl3KoUOHeO+995xyOwd74e7uzhdffEFCQgKTJk2ytzkuhRJKhUOh0WiYNGkSjRs3Zvjw4fY2x+no2LEjI0aMYPr06dy4ccPe5rgMaoxS4VAsXryYAwcOMG7cOJUOZCbvvPMOGRkZahWThVBjlAqHIjc3l3HjxtGsWTMee+wxe5vjtISFhTFkyBC++OILbt68aW9zXAIllAqHYfHixcTFxfHee+8pb7KMvPPOO6SmpjJz5kx7m+ISKKFUOARSSj766CNCQ0MZOHCgvc1xeiIiIujfvz8zZswgJSXF3uY4PWqMUuEQbNy4kQMHDjB27Fin/tE5Ev/3f/9HUlISc+bMsbcpTo0ao1Q4DLNnzyYwMJARI0bY2xSXoX379nTu3JmvvvoKjUZtO18WlFAq7M61a9dYuXIlTzzxBBUrVrS3OS7Fc889x5kzZ9iyZYu9TXFqlFAq7M7cuXPJzs5WZdSswKBBg/D391fhdxkQQqgxSoV9kVLyww8/cO+99xIWFmZvc1wOLy8vnn76aZYsWUJCQoK9zXFalEepsCu7d+8mLi5Orem2Ii+//DIajUalCpUBJZQKuzJ37lwqV65csA2toyOlLPZwZEJCQnj00Uf5+uuvVbEMM1FCqbAb6enpLFy4kIEDB+Lr62tvc1yaV155heTkZBYuXGhvU5wONUapsCsrVqwgOTmZkSNH2tsUlyc6OpqIiAhmz55tb1OcElf2KD3sbYCiZObNm0fdunXp3r27vU0pkmeYk5NT7HuGwmz9G0V3CaaHh4fB82yFEILnn3+eMWPGsH//flq3bm0XO5wVZxZCYyiP0oHJyclhy5Yt9O3b16nDFmdi6NChCCFYsWKFvU1xOlzZo1R3nwOzb98+0tLS6Nq1q71NKTcEBATQoUMHVq1aZW9TnApXH6NUobcDs3btWoQQdg27s7Ozi30MpQ+99W8U3ee6obduxXbd123FQw89xDvvvENCQgLBwcE2799ZcWaP0RjOK/HlgMWLF9OpUycCAwPtbUq5YsiQIQDMnz/fzpY4Fyr0VticuLg4jhw54jS5k65E48aN6dixIz/++KPD5386EtYOvYUQnwghTgghDgshlgohqlrIdKNYPa4RQngBfwIV8/v7XUr5rrX7dXaWLFkC5K1DtjWZmZkFj7Oysgoe64feubm5xT42NfTWnfXWvd5QGG/LTdSeeuopnnvuOfbv30+bNm1s1q+zYiOPcQPwf1LKHCHEFOD/gP9au1OwjUeZCfSQUkYCLYE+QogONujXqdm0aRMtWrSgTp069jalXDJo0CDc3d1ZunSpvU1xGqwdeksp10sptQPjMYDNBpCtLpQyD+2asAr5h4pnSiArK4u//vpLzXbbkerVq3PvvfeyfPlye5viNJgolAFCiL06x2gzu3sGWGM560vGJlOKQgh3YB/QCPhCSrnbFv06K3v37iU9Pd2mQqkbYhsKvXUfg2VDb93ZbUNt6Xsk1p4Rf+SRR3j11VeJj4/nnnvusWpfroCJY5CJUsq2ht4UQmwEahbz1ttSyuX557wN5AA2m22zyWSOlDJXStmSPFc5SggRoX+OEGK09q9Med9redOmTQghlEdpZx544AEA1q9fb2dLHB9TvElTQm8pZS8pZUQxh1YknwIeBkZIG8602XTWW0p5G9gK9CnmvTlSyrZSyrblPR1m06ZNtGzZkoCAAHubUq5p0qQJdevWZcOGDfY2xSmw9hilEKIPeZM3/aSUdy1itInYYtY7EMiWUt4WQlQCegFTrN2vs5KVlcWuXbt4+eWXrd6XbsJ4RkZGsY9LmvXWvd6SobehvWv0bzTd9qyx6kMIQe/evVmyZAkajcapV5bYAhvMes8iL3tmQ35fMVLK563dKdjGo6wFbBFCHAZigQ1Syj9s0K9TcuHCBbKysmjevLm9TVEAHTt25Pbt25w+fdrepjg81s6jlFI2klLWlVK2zD9sIpJgA49SSnkYaGXtflyF+Ph4AKtMHuh7erqTNqY8tqZHqZsjacpySP32KlWqZPC8shAVFQVATEwMTZo0sUofroCzr7wxhoolHIxTp04B1hFKRelp1qwZVatWZdu2bfY2xeFx5SWMqiiGg7Fz505q1apF7dq17W2KgjyvtVu3bmzevNnepjg8ji6EQghvIENKmWv0ZD2UUDoQUkq2bt1K9+7drfKj08+DNBRi607m2Cr0NrSEURdTC/9aeqljz549WbZsGWfPnqVBgwYWbduVcLTJLiGEGzAMGAG0I2+VYEUhxA1gNTBHSnnKlLaUUDoQZ8+e5cqVK3Tp0sXeppQaKSWXL1/m+PHjnDhxAiEE/v7++Pv7U69ePUJDQ/H09LS3mWbRs2dPADZv3syoUaPsbI1j4qCh9RZgI3lrwo9KKTUAQohqQHfgIyHEUinlPGMNKaF0II4dOwZAZGSknS0xjeTkZGJiYti5cye7d+8mKSkJoEAQdT3YypUr06JFC9q1a0evXr0ICgqyi83mEBYWRrVq1di1a5cSyhJwQKHsJaXM1n9RSnkLWAwsFkKYFH4ooXQgjh8/DkB4eLjF2tQNg3XDaP3npoTh2dnZZGdns2PHDv744w92796NRqPBz8+PqKgoIiIiCA8PJyQkBHd3d+7evUtSUhJnzpzh4MGDHDx4kJkzZzJr1ixatGhBnz59uO+++/D09HTo0FsIQXR0NH/99ZdF23U1HE0oixNJLUKInVLKTiWdo4sSSgfixIkT1KpVCz8/P3ubUgSNRsO8efOYP38+ycnJBAYGMnz4cDp16kRYWFihm0Qrzt7e3nh7exMcHEzXrl1xc3MjISGBTZs2sXHjRqZMmcL8+fN5+eWXHX65ZnR0NKtWrSI5Odkh/38cAUcbozRCqWZLlVA6ECdOnCAsLMzeZhQhLS2N9957j61bt9KxY0cGDRpEVFRUIc/PkBeoT3BwME899RRPPvkk+/btY8aMGfz3v/8lOjqaN954w2HLymlrUh48eNDhRd0eOOIYpRBiJnAk/zgqpUzRebtU68SVUDoIUkri4uIKtiGwFCVV/9ENsdPT0wse64bbJ0+e5H//+x8JCQm88MILPProowghCsJwLYZEU/fm0fc4IiMjmT17NsuXL2fevHk8/vjjvPrqq/Tt27fYz1JS9SDdcFvXLkuF4a1a5a2Z2L9/vxJKAziaUJInkC3Im/WOEELc4R/h9C1NQ0ooHYTExESSkpIIDQ21tykF7Nu3jzfeeAN3d3emTp1KixYtLN5HhQoVePTRR+nevTuTJk3iww8/ZP/+/bzxxht4e3tbvD9zCQoKolatWhw8eNDepjgsjiaUUso5us+FEMHkCWdzYF1p2nKqQQVXJi4uDsBhhHLjxo28/PLLVK9enW+++cbq2yHUrFmTzz//nJEjR7J+/Xqef/55bt26ZdU+S0vz5s05evSovc1wWBxtu1qhp9xSygQp5Wop5RQp5ePFnWMIox6lEKK2lPKyeaYqTMVaQmmoIK/+c93z5s2bx4wZM4iIiGDChAlUqVKFzMzMItebknBeUuhdXML5sGHDCA8PZ/z48bz44ovMnDmzYBdK/esNhd66n8WSM+ARERF8+eWX5ObmFrJd4ZhjlOQV41kMLJdSXtC+KITwBO4FniIv13KusYZMkfjVQoj38kukKaxEfHw87u7uhISE2NWOn376ienTp9O5c2c+++wzqlSpYnMboqKi+Oijj7h+/TpjxozBUQo5N23alIyMDC5evGhvUxwSB1zr3QfIBRYIIS4LIY4LIc4Cp4DhwGdSyrmmNGSKULYFkoHdQognzTRYYYQbN24QEBBgV09lwYIFTJ8+nR49ejBhwgQqVqxoN1tatmzJlClTSExM5OWXXyY1NdX4RVZGOyN/+bIKsIrD0YRSSpkhpfxSStkJqA/0BFpJKetLKZ+VUpo84Gw09M7f9ewzIcRc4D0hxPPAf6WU2820X1EMN2/epHr16hZpSzcMNnXWe+HChXz66ad06dKFsWPHkpubS25ubqHZcP3ry7rW21CxXu31jRs3ZtKkSbzxxhuMHz+eadOmFWrDUJk2XYH38vIq9nxzqFWrFgBXrlwpUzuuiiPnUeYnlpv9H2fKGOU9wP1AaP7RCPghf+nPOSmlypWwAJYUytKyf/9+pk6dSnR0NO+//34hAbQ3rVq14qWXXmLGjBl8/fXXvPDCC3azRVvR6dKlS3azwVFx0DFKAIQQrxXzcjKwz1Sv0pQ/AZsAv/x/XwFq51carg+oUNxCXL9+nRo1ati83+zsbD744AOCgoJ4//33Lb78zxIMGDCAhx56iLlz57Jmjc12KC1CQEAA3t7eBcWVFYVxtNBbh7bA80Cd/GM00A34Rggx1pQGTMmj7C2lLLYOvpTyvGl2Koxx/fp1unXrZpG2DIXbxc16//jjj8THxzNp0iTc3d3JzMwsFG4bWgMOhRO7zZn1NhZ6617//PPPc+nSJSZNmkTt2rUJCwszONNtaJ+fsobeQghCQ0M5ceJEmdpxVRzVowSqA62llKkAQoh3gd+BLuRto/2xsQaMepSGRFJhOTQaDbdu3bJ56H3p0iW+/fZbunTpQnR0tE37Li0VKlTg/fffp2rVqkyePNluwwOhoaGcPHnSLn07Oo6WR6lDPUB3gD0bqC+lTCevRqVRHHf0tRyRnp6OlBJf31KtqioTUkomT56Mh4cHb775ps36LQtVq1bl1Vdf5ezZs/z22292saF+/fpcvHixkAetsNy+3lbiFyBGCPGuEOI9YCd5KUPewHFTGlBLGB0AbeqLj4+PRdrTDTcNhaTLly9nz549vPbaa1StWrVQ+o2hNeD6s96GQm/dMLqkLWVLE3pDXujcrl07oqOj+eGHH7jvvvuoWbMmYDjhXNdG3Rlwc6lfvz45OTlcvXrVYQt42AtHDb2llB8IIVaTl2QugOellHvz3x5hShvKo3QA0tLSAGy2tjk9Pb1g5U3//v1t0qcleemll5BSMmvWLJv3Xa9ePQDOn1fD8/o4sEcJkANo8v81qQalLkooHQCt12atLVf1WbhwIYmJibzwwgsOnftmiKCgIIYOHcqmTZtsXqRCK5RqdU5RHHWMUgjxCjAfCABqAPOEEC+Xpg0VejsAlhZKQ6F3ZmYmqampfP/993Ts2JHQ0NCCMFu3tJqhkmu6j6H0a731Z50Nhd666N5cuuf379+ftWvXMnXqVObMmVMoyVz38xsqBafftqlohfLChZQUFDIAACAASURBVAtGzixfOIDHWBKjgPZSyjQAIcQUYBcw09QGnM+dcEFs6VH+9ttv3Llzh+eee87qfVkTLy8vnn/+eU6fPs2WLVts1m+VKlXw9fUlISHBZn06Cw4cegvy1nxryc1/zWSUUDoAthLKjIwMFixYQMeOHS26L4+96N69OyEhIcybN8/kCuuWoEaNGg5TqMORcGCh/IG8WhXvCSEmALuB70vTgAq9HQBLCKVuGGwo9Fy8eDG3b99mxIgRRcqmlXaPb/22Tclr1A2d9Z+XNvSGvJnuYcOG8dFHH7Fhwwa6d+8OGJ711rfR3O1za9SowfXr18261pVx1PFuKeU0IcRWoBN5nuRTpSmIAcqjdAi0N7M1lw9mZmYyf/58IiMjrVKp3F507dqVWrVqsWDBApv1Wb16dW7evGmz/pwBR8yjFEKkCCHu5G8BsRX4EJgEbM9/zWSUUDoQ1vwhLVu2jMTEREaOHGm1PuyBu7s7ffv25fDhw5w6dcomffr6+hakdCn+wdGEUkrpK6WsonP46hylKrSqQm8XwVC4nZWVRXZ2Nt988w1NmzalRYsWBe+bEmKXNOtd2oRz/VlvXQ/alJny4kJvgG7duvHjjz+ycOFC3nrrLYPfhaVCbx8fH4eoj+loOPCsd5lRHqUDoBUDa61f3r9/P1euXGHIkCEu+WP29fWlT58+bNiwgaSkJJv0d+dOqSK3coGj5lFaAuVROgBaz0Z/iWBp0PXo9D2qzZs3U7FiRSIjI0s9gWOqR6kr8oa8Q32PsLS5l/pjuLr99+/fn+XLl7Ns2TJeeuklozaWBT8/P9LS0sjJySnymcortgithRAfAP3JW2FzHRhpq/28nFfiXQhtsrT+rLIlkFKyfft2oqKiLLLW2VGpX78+bdu2Zfny5YXE0Rpo9xFSXmVhbDBG+YmUsoWUsiXwBzC+7FabhtWFUghRVwixRQjxtxDiWP5yIoUO1hTKM2fOcOXKFTp37mzxth2Nhx9+mJs3b3LkyBGr9qNdk3/37l2r9uNsWDv0llLq/mXyBorfc8QK2CJuyAFel1LuF0L4AvuEEBuklCaVNyoPWCL01g0rdR9v27YNgHbt2pW4Z44p4bbu0kYwLazVvTn0zzGUO2koXC+pelF2djbNmzfH3d2dHTt20Lx58yJ9Wir01g4BWNtzdTZM9BgDhBB7dZ7PkVLOKUUfk8jbWSEZ6F46C83H6h6llPKKlHJ//uMU4G/yyrEr8tF6lGURSkP89ddfNGnShICAAIu37Wh4e3sTHh7Orl27rNqPEsqilCKPMlFK2VbnmKPXzkYhxNFijv4AUsq3pZR1yStyMcZWn8+mY5RCiBCgFXlLiBT5aFfkWDo378aNGxw9epROnTpZtF1Hpm3btpw8edKqCeHaCSZVvLcwlhijlFL2klJGFHMs1zv1F2CQVT5IMdhsyk4I4QMsBv6jN9agfX80eZv+FFRoKS9oJwdSUlLMbqO4LWrXrFmDlJLOnTuTk5NTYuhtKAzXHYcrqXqQobC2pFlvU5Yt6tqln/dY3FLFVq1aMXfuXHbs2MGDDz5osMJRWdDaZ2hb3vKKtdN/hBCNpZTaVQX9AJttXmQTj1LkbW27GJgvpVxS3DlSyjladzwwMNAWZjkM2i0gLD2Lun79eurXr0+DBg0s2q4j06BBA6pXr05MTIzV+tCKvxLKwthg1vuj/DD8MHAfebvC2gSre5Qi79v5DvhbSjnN2v05I97e3ri5uVlUKO/evcuBAwcYPHiwxdp0BoQQtG/fnj///NNqCfxa71XlUP6DLfIopZQ2C7X1scX/dCfgCeCIEEJbseN/UsrVNujbKRBCULVqVW7fvm12G7phZW5uLvv37ycnJ4fmzZsX3Nj6obehrWzNSTg3FNbqhmP655iyra2hvXD0+9d93L59e1avXs3hw4fp0aOHwf51n5dmK1tbV6R3Flxx1ZcWqwullHIHpSySWR6pWrWqRZff7du3D3d3d5o2bWqxNp2Fdu3a4eHhwa5duwoJpaVQQlk8zrxE0RgqdnAQAgICSExMtFh7x48fp2HDhlSuXNlibToLlStXJjw83GqJ51qhLI/frSHsXJjX6iihdBBq1KhRpu0F9EPvuLg4OnToYHD/HCh9uG1O6K178+iv1dad9Ta0vttQEV795/oz8GFhYSxbtoz09PSCfvVt1O2/NKG3NhNAeZSFcWWhdF1f2ckICgri2rVrFmkrKSmJmzdv0qRJE4u054yEhYWRnZ3NyZMnLd52eno6np6eLh1qmoOj1aO0JOp/2kGoXbs2165ds8hM7dmzZwG45557ytyWs9K4cWMAqwjlzZs3qV69usXbdXZUmTWF1alXrx4ajYbLly+blXCvG0ZqhTI4ONjgzDBYbtZbSklycjKZmZlkZ2eTnZ2Nj48P/v7+hUJa/QRz3VlvQwV+dZPM9YuGlLQ3TkBAABUrVuTUqVMFIXdJ/ZeGq1evUrNmTbOudVWc3WM0hhJKB6Fu3bpA3n7RZV2ZdP78eTw9PQkKCrLqJljp6ekkJCRw4cKFIgUzIE/k6tatS0hIiM1LvLm5uVGvXj3OnDlj8bavXr1KUFCQxdt1dpRQKqyOVigvXrxY5rbi4+OpX79+qSYoSkNGRgbHjx/nwoULSCkJCAigUaNGVK5cmQoVKuDh4cHt27e5cuUKZ86c4ezZs9StW5fQ0FCbzhTXrl2b+Ph4i7d7/fp1mjVrZvF2nR0llAqrU1ah1A0r4+PjiYyMJDc3t9CYp37oWtpq52lpaZw9e5b4+Hg0Gg116tQhODiYypUrFwlj/f398ff3p0GDBly4cKHgqFu3LmFhYUXGqwxVMi9pZl1/pl//cbVq1YiJiSE3NxchhEVCbykl169fp7wtszUFZx6DNIYSSgfB19cXPz8/Lly4UKZ2UlJSuHbtmsUncrKysoiJieHOnTsEBQURGhpq0hK+ypUrExYWRv369bl48SIXLlwgIyODli1bWtS+4ggICCAzM5Pk5GSqVq1qkTbT0tLIyMhQQqmHq49Ruu6fACckLCyMY8eOlamNGzduAHl5mZYiMzOTHTt2kJKSQuvWrWndunVBlW9TqVSpEhEREYSHh3P9+nX+/vtvqxeV0NbgtOQ4rdbjDw4OtlibroIrpwcpj9KBaNmyJYsWLUJKWeoflTasvHXrFpDnoeqH3vqpR6bMeicnJxMbG0taWhrNmzfHz8+v4H3d6w0lnOt6nVJKatasSVpaGhcuXKBSpUo0bNgQKBy26c50lzRrb6iEmva70KbwXLlyhUaNGhURZnOE+ty5c0DeHj2KwjizEBpDeZQORGRkJLdv3y5T+K0trKEt3VYWNBoNe/fuJS0tjTZt2lgsd/Cee+6hRo0anDx5ssADtgZaj9JSifyghLIkXDmP0nktd0Fat24N5O3DbS5Xr14FsMjWD3Fxcdy5c4fIyEiLbiUhhCA8PBxvb2+OHz9utUrhfn5+AGWqyqTP33//jY+PD7Vr17ZYm65AKbaCcEpU6O1AtGjRAnd3d/bu3cuAAQPMauPKlSu4ubnh7++PRqMpMXQtaaY7LS2N48ePExAQQNWqVcnIyCiSK6kbehtat21o726ARo0acejQIU6fPk1ERESx7ZZUodzQrLfWFjc3N3x9fbl9+zZSSouMiR49epRmzZo5tXdkLZxZCI2h/rcdCO2ER2xsrNltXL58merVq5epqKyUkkOHDiGEKFgKaA2qV69OQEAA586ds9rWr35+fiQnJ1usvaNHjxYSdcU/uLJHqYTSwejYsSO7du0yOxy9fft2mVNhzp49S2JiIk2bNrX6ihqtEB89etQq7fv6+pZpLyJd0tLSuHHjRsEElKIwrjxGqUJvB6N169Z89dVXXLhwway9bpKTk/Hx8SkIP3VD4pLKlGlD79TUVI4fP05gYCA1a9YsJDKmllkraX23Lu7u7lSoUIH69etz9uxZLl26RI0aNQzO1OvP2hcXbus/9vHxITU1tdiwu7QezuXLlwGoU0fttqyPs3uMxnBeiXdRQkNDAThxwrwN5u7cuYOPj4/Z/WvzOJs2bWqzH36dOnWoXLmyVZYbent7W2wbYG29UCWUxaNCb4XNCA8PB/JmV80hPT3d7HA5MzOTCxcuULt2bZsWsXBzc6NmzZrcvn3b4jPgbm5uJXq1pSEuLg7AquO2zowrC6UKvR2MgIAAatSoYfYKnYyMDDw9PQtCTV3hMZZwfvr0aXJzc6lZs2ZBKF5SmTVDM9K6obfu6/o3iu55vr6+SCm5efNmoVU/hma29Z8bmtHWTd7X77+0Y2Z///033t7eBevyFYVx5jFIY7juJ3NimjVrZvZ+L1qhLC1SSuLj46lWrVqZQndz0eY8WnKDNcCsVU6GOH78OGFhYU7tGVkLV8+jVELpgLRr146DBw8W8eBMITs7u8jeNKZw5coV0tLSaNSoUamvtQSenp74+vpy8+ZNi7abkZFhkWEE7SqlNm3aWMAq10QJpcKmREdHk52dzb59+0y+Rpt+kZOTU0goNRpNwZGTk1PoyMrKKjhOnjxJpUqVqF69OhkZGQWH7jmZmZkGD21l8+zs7ELX6B6652RnZ5Obm1voqFatGsnJyWRkZBSyW3tok8aLSx439Hp6ejqVK1dGCFGmtJXTp09z+/ZtoqKiTP4/KW8ooVTYlOjoaAB27txZquuklOTk5JS6YO+tW7dISkqiQYMGdh1n0q4lt+T677S0NIsUC46JiQGgffv2ZW7LVVF5lAqbEhQURHh4OJs3b2bs2LEmXaP9EXp5eRUK2XVnfPUnc7TPtaXDqlevXuA5atF9rL/drX5epr4t+hjqH/ImZrSeX1JSErVq1QJMr/Cj663oPk5OTiY8PBw3N7cif0BK8wdl27ZtVKtWjaZNm5p8TXnC2T1GYzivxLs4PXv2ZPv27UXEyRje3t7F7l9TEpcvX8bb29vm+9ro4+bmhre3t8WWHEopSUpKwt/fv8xtbd26lS5duji1V2RtVOitsDm9evXi7t27bNu2rVTX+fj4lCrBOicnh+vXrzvM9qs+Pj4WXXKYk5NTZqGMj48nPj6ebt26WcQuV8WVhVKF3g7K/fffj5+fHz/99BO9e/c2er62CIafnx93794t8Hx0Q1f9xOvs7GyuXr2KRqOhSpUqBd6roYK++nvuGNqDXLcgh6GcSn17tHZWqFCBrKwsNBoNQgizQm/tZ9eWVwsKCsLd3d3szdZ+/fVXAB555BGzri8vuLK37bqfzMnx8vJi+PDhLF68mDt37ph8nZ+fH6mpqSaff/XqVdzc3Cy2p0xZ8fDwQEppkRU6iYmJQNm3xVi0aBEdOnRQxXpLQOVRKuzGiBEjSE9PZ+PGjSZfU6VKlVKFrtqw21pb25YWrR2GJopKg3b2vCx7cJ86dYqDBw8yZMiQMtvj6thKKIUQbwghpBDCctWkjaBCbwemXbt2eHp6EhMTw8CBA0s8V5s7Wa1aNVJSUooNg/S9tKysLG7fvk2jRo0Mzm6XNOtd0vJELbphuH7/JYXVbm5uBbmPhvowFOJrX79+/ToeHh7Url27oFJRadm0aRMA/fr1K/W15Q1beIxCiLpAb6Bs25WWEuVROjAVK1akdevWBTl8plC1alWTQ+87d+6g0WgcJuyGf8TTEjfdpUuXCA4OLpO3/Oeff1K7dm2Lb//ritgoj/IzYCxg3S089VBC6eBER0cTGxtrcspPtWrVyMrKMmnmW5uG40hCqbuNQ1m5dOlSmQpY5ObmsnnzZrp06eLU42u2wBZjlEKIfsAlKeUhy1htOlYXSiHE90KI60II65SwdnHuv/9+MjIy2Lp1a4nnVaxYkYoVKxaUALt69SoeHh6F/prrLwlMTk5GCIGXl5fBpY0lHbrLEXWvL275YXGlzoq7kdLT0/H09KRixYpFvBEPDw+DR4UKFQoODw8PhBBcvnyZxo0b4+npWdCm7mGMHTt2cO3aNbP3LypvmCiUAUKIvTrHaL02NgohjhZz9AfeBsbb47PZYoxyLjAL+MkGfbkcXbt2pVKlSqxevZoHHnjA6Pk1a9YETKvCk5qaire3t0Olddy9e9ci1YuuXLlCdnZ2mYp8/P7773h5efHggw+W2Z7ygIm/o0QpZVtDb0opexX3uhCiOdAAOJQvuMHAfiFElJTyqhnmlgqr3yFSyj+BW9bux1Xx8vKid+/eLFmyxKSZ4MDAQKB0QukoSCm5e/euRfYk1xY+btKkiVnXZ2RksGjRIh566CG7lJ1zRqwZekspj0gpa0gpQ6SUIUAC0NoWIgkONOud74KPBqhXr56drXEsRo0axYoVK1i5cqXB2e9KlSoBFOyzk5SUhKenZ6G/8rqzzFJKUlNTCQgIIDc3t5AI685u676un2CuG07rzkDr9mNoDbb+c3d3d1JTU8nJySEgIKBgAkZ3plp/Z0nd93RrcFasWJGdO3dSo0YNWrZsWdBPaZZozp8/nxs3bvDSSy+ZfE15xtnzJI3hMDGXlHKOlLKtlLKt1itS5PHQQw9Rt25dvvrqK6Pnent74+Pjw/Xr10s8LyUlBY1G41DeknYlTUBA2dLjUlNTiY2NpWvXrmbdvFJKPvvsMyIjI9WyxVJgy4TzfM8y0WINGsFhhFJhGHd3d5577jk2btzI4cOHjZ5fo0YNo6XKtKt9LFGCzFIkJSVRsWLFMtu0ceNGsrOz6d69u1nXr127lmPHjvHaa6+5tJdkaVy5zJrzWl7OePHFF/Hy8uL7778v9n3dH2NQUBA3b96kYsWKhWaGdf+ya0uxeXh4FFvU15TD0Iy2IW9Cf6ZaOxvt6emJh4cHSUlJBAYG4uXlVWSGumLFilSqVKnQ4eXlVeyxZs0awsLCaNOmTaHz9fs3xNy5c6lRowbDhw+33H9gOUAtYSwDQogFwC4gVAiRIIQYZe0+XRF/f386d+5csFKkJAIDA7l27VqJ52iF0pzVKtZAuwNjWcPuixcvcuLECR544AGzbsy0tDT++OMPBg0a5DDfjTOg1nqXESnlcCllLSllBSllsJTyO2v36ar07t2bo0ePcvbs2RLPCwoKIjEx0WB1H/hnwsZRxODixYu4u7uXudzbli1bcHd3N6niUnH88ssv3L17V3mTZqCEUuEQDB48GPin7Jch7rnnHnJzc7lz506h8FZbaszd3b1ARLWJ6Lp71xja16YkdG8G3WEA3f71D21InZOTw9WrV2nYsCE+Pj6Fwmjd0Lly5cqFDm9v70JHxYoVCwrs1q9fv8j7xtBoNEydOpXWrVtz7733mvi/otCixigVDkFISAgdOnRgwYIFJZ6nTa/SbvFQHNnZ2QWFJ+zNhQsX8PDwoGHDhmVqZ8+ePdy6dYtBgwaZdf3KlSuJi4vjzTffdIjvxdlQHqXCYXjiiSc4dOhQiWOVISEhQF6JMEPk5uY6xF/45ORkEhMTadiwoVn7kWvJzMzk22+/JSgoiM6dO5f6+uzsbP73v/9xzz338Oijj5ptR3nF1ccoHSbhXGEazzzzDB9//DFjx44lNja2WLFr3LgxNWvW5NixY7Rs2bLgdd2ZXjc3N4PVz0uqiq6Lbnu6Y52666h1H2uT4rV4eXlx4MABvLy8aNGiRUF7uitzdB9XqVKl0PW6z3/++WcuXrzIjBkzCm39oH+NIWbOnMnx48dZsWJFiTPiCsM4sxAaw/4uhaJUeHl5MXHiRPbv328wBBdC0LZtW/bu3WuwHUf4UZ87d47U1FTCw8PLJE7Hjx/n119/5dFHHzVr3+2srCw++eQTevfuTd++fc22o7yjxigVDsVjjz1G8+bNmTZtmsFzoqKiOHfuHLduFb/MXgjT96OxBsnJyZw+fZoaNWqUaauG27dv8+mnn1KrVi1efPFFs9r45ptvuHr1Km+++abZdihce4xSxRhOiJubG48//jj//e9/uXjxYpGai/7+/vTp04fx48dz+vRpevToARRdD63dwEtfNA0JqL5HoOsF6q6j1g2xdWebtcslU1JSOHPmDFWqVKFjx45UqFChUE1MQ4/1d1P09fVl/PjxpKam8t133xEcHAzk1eQ0lbS0ND744AO6dOlCr17FFq5RmICzC6ExlEfppPTv3x+A5cuXF/t+ZGQk/v7+HDhwoNj3tSJna68yPT2dw4cP4+HhQefOnc3O45RSMmPGDA4dOsTYsWNp2rSpWe1Mnz6da9euMXnyZJe+0W2BK3uUSiidlNDQUCIjI/nyyy+LzXHUJl3/9ddfRfa6gX+8y5KS0i1NWloa+/fvR6PREBkZWWRyx1Q0Gg3fffcdy5YtY+jQoWYnlx86dIj333+fAQMG0LFjR7PaUPyDK49RqtDbiXnnnXcYPHgwv/32G8OGDSt4XVuibPTo0fz6668cPHiQ++67r5AwacPYnJwcKleubLAcmq7Hp/9D121PtwpRceFySkoKBw8exM3Nja5du+Ln51doFY7u0kXd6lHaQsSQt+IoJyeHKVOmsG7dOp555hnGjh2LEKLQOKcpyeVZWVkMHz6catWq8fXXXxs9X2EcZ/YYjeG8Eq9g4MCBREREMGHChGKL+nbv3p1KlSoVW3FIK2y6uyxai+TkZA4ePIi7uzvdunXDz8/PrHbu3r3Lu+++y7p16xg1alSBSJrDp59+yt9//813332HKutXdlw9j1IJpRPj5ubGxIkTOXHiBBMnTiz2/UaNGhWbeK7NT9QWx7AWV65c4cCBA3h4eNCtWzez61+eOXOGUaNGsXPnTl5++WWeeuops2+8jRs3Mn78eIYMGWLS9hoK03BloVSht5PTv39/HnvsMSZPnsx//vOfIjPDLVq0YMWKFQQEBBRK3tYmgmdlZRVs5KVF97Fuwrh+hXBD4Xa1atXQaDTExcVx6dIlgoKC6NixY8GstBbdcFs3dA4KCir0eOHChXz11VcEBgayYMEC2rdvX+Q8U6uXnz59miFDhhAWFsa3335r0jUK03DmMUhjuO4nK0e8+uqrZGdns3jx4iLvNWvWjJSUlCLrvrVFKe7evWtxe9LT09m3bx+XLl0iPDycrl27mrTroT7Jycm8+eabzJw5k44dO7J69eoCkTSHpKQk+vXrhxCCFStWWGRvHsU/KI9S4dC0adOGiIgIPv74Y5566qlCEzCtWrUCYNeuXUXCXh8fH1JSUixmh0aj4fr16wW1MJs3b06zZs3Mamv//v3MmjWLtLQ0XnvtNR599NEi3nJpuH37Nvfddx9nzpxh7dq13HPPPWa3pSiKswuhMZRQugBCCCZPnkzfvn355ptvCq1Q6d27N3Xr1mXnzp089thjha6rXr0658+fL6gWrsVQGKs/CaM7a52bm8vRo0dJTU2lXr16tG/fHh8fH4MhNRSe3dael5KSwqJFi1i0aBENGzZkxYoVBWJbmkRyXe7cuUOfPn04dOgQS5cuNXuLCEXJuLJQqtDbRXjooYfo0qULEyZMKLRsUQjBww8/zMaNG4vMcFetWpXMzEzS09PN7jcjI4O4uDhiYmLQaDRERUXRs2fPUk/a5Obmsnr1akaNGsVvv/3GiBEjWLBggdkeqZarV6/ywAMPsG/fPn777TceeuihMrWnMIzKo1Q4PEIIPvvsM6Kjoxk6dChr167F3d0dNzc3HnvsMb766itiY2MZOnRowTUhISHs2bOHtLS0QrmHuhV3dL1LXY8wKyuLK1eucOrUKYQQdOjQgaioKDw8PAp5kbVq1Sr2MeR5kVlZWWzYsIEpU6Zw9OhRoqKi+PrrrwtVPTKH3NxcZs+ezdtvv01GRgYLFy4sWM2ksA6u7FEqoXQhWrduzcyZM3nuuedYtWoV/fr1A+Dee++lR48e/PLLL/Tv379A/IKCgvD09OTChQsmlyPTaDQkJCRw+vRpcnNzadSoEZGRkaXei/3atWv89ttvLF68mFu3btGwYUNmzZrFgAEDqF27duk+uA6ZmZls2rSJd999l71799KrVy+++OILmjRpYnabCuOoMUqFU/H000/z9ttv88MPPxQIJcCECRPo3Lkzy5YtK1jF4+7uTnBwMOfPnyciIsLoDz0pKYkTJ06QkpJCtWrVuPfeewulBRlDSsnBgwdZvHgxsbGxAHTp0oVhw4YVzEabQ0ZGBlu2bGH+/PmsWLGClJQUgoKC+OWXXxg2bJhL38COhCt/z0ooXYwKFSrw4osv8v7777Ns2TIeeeQRIM+rbNOmDTt27ODdd98F8pYHtmvXjkWLFpGVlVXgyemOJQUEBJCVlcWxY8c4d+4cvr6+9O3blyZNmhRUUofCeZC6kzR169YlPT2dFStW8Omnn3LkyBGCgoJ46623eO6556hfv75Zn/PSpUusW7eOP/74g/Xr15OWlkbVqlUZMmQIAwYMoGfPnibnViosgzOPQRpDCaUL8vbbb7Ny5UqeffZZOnToULBeukePHnz66afcvn27wBOMjIxk8+bNnDhxgho1ahQqnSal5OLFixw5coTMzEzatm1Lx44dTdqyQUrJ6dOn+fHHH1m6dCl37tyhRYsWfPXVVwwfPrxg725TycrKYvv27axevZp169Zx7NgxAIKDg3niiSfo27cvPXr0UOJoJ1TorXA6PD09mTdvHm3atOHJJ59k5cqVVKxYkc6dO/PJJ5+wffv2gkrebm5uPPDAA/z000/8+eefhIaG4uXlxc2bNzl//jzp6en4+fnRsWNH2rRpY7Tv1NRU9uzZw9KlS4mPj6dChQr07t2bxx9/nAEDBpT6Zjpx4gRTpkxh8eLFpKSk4OnpSZcuXRg5ciT33XcfzZs3d+kb1Jlw5f8HJZQuStOmTfniiy8YNWoU06dPZ+zYsfTp04c6deowb948nnnmmYIJjiZNmnDnzh22bdvGvn37Ctrw9/fnkUceoVmzZri5uRXZZqFJkyZIKTl//jyHDh1i+fLlbN26lZycnALvcejQoWYliu/bt4/Jkyezj3h7PwAAB91JREFUZMkSvLy8GDFiBP369aNHjx4mVQdS2B5XFkphz+0ADNG2bVtZ0n4vCtPp0qUL165d48SJEwghmDZtGq+//joxMTFcvny54Lxt27YhpSQ+Pp7ExEQqVKhAw4YNC80Wa4UyKSmJ2NhY4uLi2L17N1evXgXy9hMfMGAAQ4cOpW3btibfOBqNhps3b7Jnzx62b9/Otm3biImJwc/PjzFjxvDKK6+oCj82In+vpVIrXuvWreXOnTuNnle5cuV9Usq2ZhlnR5RH6eI8/fTTPPPMM6xYsYL+/fvz7LPPMnHiRMaPH8/o0aMLDcALIWjYsKHBPWyOHTvG3LlzOXjwIFJKAgMDiY6Opn379vTp04ewsDCEECUW5M3NzWXPnj2sXLmSDRs2cP78eW7evFmw26OHhwdt27Zl8uTJvPDCC2aXZFPYFjVGqXBqhg4dyvTp0xk2bBirVq2iR48efPDBB4wZM4bw8HDGjRuHEKKQYKalpRU87ty5M7du3WLy5MksWrSIWrVq8e6779KvXz9atmxp0s2RmprK+vXrWbFiBatWrSIxMRF3d3c6derEwIEDCQwMJCAggBYtWtC+fXsqV65sle9CYV2UUCqclsqVK7Nx40a6d+9O3759Wb9+PS+++CKHDh1i+vTp3Llzh48++qjYa7Ozs5k3bx4ff/wxqamp/Pe//2XcuHEGxwizs7NJTU3lzp07HD58mJiYGGJiYti5cyeZmZlUrVqVBx98kL59+9KnT59S5WAqHB8llAqnJiAggI0bN9KhQweeeeYZYmNjmT17Nl5eXsycOZNTp04xYsSIQksMY2Nj+emnn7h16xbt2rXjww8/LNilUKPRsGrVKn755Rf27t1LYmIiaWlpRaqsu7u7ExkZyQsvvED//v3p1KmT2ZuJKRwfa+dRCiHeA54FbuS/9D8p5WqrdpqPEspyQlBQELNnz+bhhx+mb9++rFmzhhkzZtCqVSteeeUVduzYQZ06dWjQoAHu7u5s27aNyMhIFixYQO/evRFCkJWVxS+//MInn3zC8ePHqVmzZkFB3kqVKuHj41NwhIWF0bp1axVGlxNsOEb5mZTyU1t0pIsSynLE/fffz88//8xjjz3GoEGDWLZsGU8//TS9e/dm3bp1bN++nbNnz3LlyhUmT57M66+/TlJSEj///DMrV65k3bp1pKSkEBkZyfz58xk8eLDyEBUFqNBb4TIMGzaMu3fvMmbMGA4dOkRUVBTBwcGMGjWKUaNGFTo3JyeHRo0akZKSQq1atRg2bBiDBw+mV69eLn1TKMzDRr+JMUKIJ4G9wOtSyiRbdOqQeZRCiBvAeSs0HQAkWqFda+JsNjubveB8NlvT3vpSylInrQoh1pJnlzG8AN0d7eZIKefotLMRqFnkKngbiCHvc0vgA6CWlPKZ0tpqDg4plNZCCLHX2ZJdnc1mZ7MXnM9mZ7PXGgghQoA/pJQRtujPdct9KBQKl0IIoVv5eQBw1FZ9qzFKhULhLHwshGhJXuh9DnjOVh2XN6GcY/wUh8PZbHY2e8H5bHY2ey2ClPIJe/VdrsYoFQqFwhzUGKVCoVAYodwKpRDiDSGEFEKYktJgN4QQnwghTgghDgshlgohHHaBtBCijxAiTghxWgjxlr3tKQkhRF0hxBYhxN9CiGNCiFfsbZOpCCHchRAHhBB/2NuW8kK5FEohRF2gN3DB3raYwAYgQkrZAjgJ/J+d7SkWIYQ78AXwANAUGC6EaGpfq0okh7yE5XCgA/CSg9uryyvA3/Y2ojxRLoUS+AwYS97smUMjpVwvpczJfxoDBNvTnhKIAk5LKeOllFnAQsBhN9KWUl6RUu7Pf5xCnvDUsa9VxhFCBAMPAd/a25byRLkTSiFEP+CSlPKQvW0xg2eANfY2wgB1gIs6zxNwAuGBguTlVsBu+1piEp+T90deY29DyhMumR5kZBnU/4D7bGtRyZRkr5Ryef45b5MXLs63pW2loLiFvg7vsQshfIDFwH+klHfsbU9JCCEeBq5LKfcJIbrZ257yhEsKpZSyV3GvCyGaAw2AQ/kL+IOB/UKIKCnlVRuaWAhD9moRQjwFPAz0lI6bz5UA1NV5HgxcNnCuQyCEqECeSM6XUi6xtz0m0AnoJ4R4kLw101WEEPOklI/b2S6Xp1znUQohzgFtpZQOWxBBCNEHmAZ0lVLeMHa+vRBCeJA32dQTuATEAo9JKY/Z1TADiLy/lD8Ct6SU/7G3PaUl36N8Q0r5sL1tKQ+UuzFKJ2QW4AtsEEIcFELMtrdBxZE/4TQGWEfexMivjiqS+XQCngB65H+vB/M9NYWiCOXao1QoFApTUB6lQqFQGEEJpUKhUBhBCaVCoVAYQQmlQqFQGEEJpUKhUBhBCaVCoVAYQQmlQqFQGEEJpcIiCCG2CiFC8x9XF0LYbOMnhcLaKKFUWIpGwKn8xy2AI3a0RaGwKEooFWVGCFGfvNJ12tJfLYDDdjRJobAoSigVlqAlhYWxDUooFS6EEkqFJYgkr+wXQojG5FU2V6G3wmVQQqmwBC0BNyHEIWA8edWDnrKvSQqF5VDVgxRlRghxGmiVv/eMQuFyKI9SUSaEEL6ARomkwpVRHqVCoVAYQXmUCoVCYQQllAqFQmEEJZQKhUJhBCWUCoVCYQQllAqFQmEEJZQKhUJhBCWUCoVCYQQllAqFQmGE/wdItp1S2l8d7AAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light", "tags": [] }, "output_type": "display_data" } ], "source": [ "fig = plt.figure(figsize=(5, 3.75))\n", "plt.imshow(logL, origin='lower', cmap=plt.cm.binary,\n", " extent=(mu[0], mu[-1], gamma[0], gamma[-1]),\n", " aspect='auto')\n", "plt.colorbar().set_label(r'$\\log(L)$')\n", "plt.clim(-5, 0)\n", "\n", "plt.contour(mu, gamma, convert_to_stdev(logL),\n", " levels=(0.683, 0.955, 0.997),\n", " colors='k')\n", "\n", "plt.text(0.5, 0.93,\n", " r'$L(\\mu,\\gamma)\\ \\mathrm{for}\\ \\bar{x}=0,\\ \\gamma=2,\\ n=10$',\n", " bbox=dict(ec='k', fc='w', alpha=0.9),\n", " ha='center', va='center', transform=plt.gca().transAxes)\n", "\n", "plt.xlabel(r'$\\mu$')\n", "plt.ylabel(r'$\\gamma$')\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "iIX6uze_Vr_R" }, "source": [ "## Posterior for Cauchy Distribution\n", "In the prevoius secsion, we use the interquartile range shortcut expressed by $\\sigma_G = 1.483\\gamma$ Here we will use a bootstrap method to estimate parameter uncertainties. " ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "HmNQV-wbVr_T" }, "source": [ "### 1. define function" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": {}, "colab_type": "code", "id": "C174nTJKVr_V" }, "outputs": [], "source": [ "def estimate_mu_gamma(x, axis=None):\n", " \"\"\"Equation 3.54: Cauchy point estimates\"\"\"\n", " q25, q50, q75 = np.percentile(x, [25, 50, 75], axis=axis)\n", " return q50, 0.5 * (q75 - q25)" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "PGy5bOH0Vr_g" }, "source": [ "### 2. Generate sample and use Bootstrap estimation\n", "Draw a random sample from the cauchy distribution, and compute marginalized posteriors of mu and gamma. Here we have a Cauchy distributed sample with N = 10, $\\mu = 0$, $\\gamma = 2$. " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": {}, "colab_type": "code", "id": "y-mTPq4QVr_i" }, "outputs": [], "source": [ "np.random.seed(44)\n", "\n", "n = 10\n", "mu_0 = 0\n", "gamma_0 = 2\n", "xi = cauchy(mu_0, gamma_0).rvs(n)\n", "\n", "gamma = np.linspace(0.01, 5, 70)\n", "dgamma = gamma[1] - gamma[0]\n", "\n", "mu = np.linspace(-3, 3, 70)\n", "dmu = mu[1] - mu[0]\n", "\n", "likelihood = np.exp(cauchy_logL(xi, gamma[:, np.newaxis], mu))\n", "\n", "pmu = likelihood.sum(0)\n", "pmu /= pmu.sum() * dmu\n", "\n", "pgamma = likelihood.sum(1)\n", "pgamma /= pgamma.sum() * dgamma\n", "\n", "# bootstrap estimate\n", "mu_bins = np.linspace(-3, 3, 21)\n", "gamma_bins = np.linspace(0, 5, 17)\n", "\n", "mu_bootstrap, gamma_bootstrap = bootstrap(xi, 20000, estimate_mu_gamma,\n", " kwargs=dict(axis=1), random_state=0)" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "dawHHBVKVr_r" }, "source": [ "### 3. Plot results\n", "We will draw a bootstrap estimation (blue dotted lines) of this Cauchy distributed sample along with the posterior pdf (black solid lines). The left two plots show how the result changes with the change of posterior $\\mu$, and the right two plots show how the result changes with the posterior $\\gamma$. \n", "The bottom two plots show the cumulative probability distribution." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": {}, "colab_type": "code", "id": "pOOoLbRkVr_t", "outputId": "e294624e-4852-4ffc-e525-914a7c9a4dfe" }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdoAAAHPCAYAAAD0/xuHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdd3xUZfb48c8hIUCABEgAKUKQoqK4kCCgLIsNQSEgP0RBxbYScRcLoKtrCxZW/OqurooIqNgpIgJBBIFVQAGVDlIkNAlFmkCoAXJ+f0wGQkiZTGbmzkzO+/WaF5k7d+49CfPMuc9znyKqijHGGGP8o4zTARhjjDHhzBKtMcYY40eWaI0xxhg/skRrjDHG+JElWmOMMcaPLNEaY4wxfhTpdAAlFR8frwkJCU6HYYzXFi9evEdVqzsdR6BZ2TWhztOyG/KJNiEhgUWLFjkdhjFeE5EtTsfgBCu7JtR5Wnat6dgYY4zxI0u0xhhjjB9ZojXGGGP8yBKtMcYY40eWaI0xxhg/skRrjDEFEIHFi10PkTOPwYNdr9eufWZbUpJrW0oK2Kglk1vID+8xxhh/cifQ/FYU3b793G0jR7oSrzFuVqM1xhhj/MgSrTHGFCA11bv3bdvm2zhMaLNEa4wxBXDfiy2uxYt9GoYJcZZojTGmALVre/e+rl19G4cJbQFLtCLyvojsEpFVBbwuIvKGiKSLyAoRSQxUbMaYgpXmsrtjh9MRmHAQyBrtB0CnQl6/AWic80gBhgcgJmNM0T7Ayq4xXgtYolXVucC+QnbpBnykLguBKiJSKzDRGWMKUprLbqKXdfMRI3wbhwltwXSPtg6wNdfzjJxt5xCRFBFZJCKLdu/eHZDgjDEFCtuy622nppQU38ZhQlswJdr8hnjnM0QcVHWkqrZU1ZbVq5e69bKNCTZhW3a9TZg2YYXJLZgSbQZwfq7ndYF85l0xxgSZsC27o0Y5HYEJB8GUaKcAd+b0YGwDHFBV6/NnTPCzsmtMIQI217GIjAGuAuJFJANIBcoCqOo7wDTgRiAdOALcE6jYjDEFs7JbfF26OB2BCSYBS7Sq2ruI1xX4e4DCMcZ4qDSXXW+nUkxL820cJrQFU9OxMcYEFW97HScn+zYOE9os0RpjTAG8nUpx6lTfxmFCmyVaY4wxxo8s0RpjjDF+ZInWGGMK4O1UiprvdB2mtLJEa4wxBfB2ZqiRI30bhwltlmiNMaYA3k6leP/9vo3DhDZLtMYYY4wfWaI1xhhj/MgSrTHGFMDbqRSnTPFtHCa0WaI1xpgCeDuVYlKSb+Mwoc0SrTHGFMDbqRTr5LvsvSmtLNEaY0wBbCpF4wuWaI0xxhg/skRrjDEF2sWyZcs4efJksd7Vt6+fwjEhyRKtMcbkY+/evVxwwRW0aNGCqlWrct111/HOO+949F6bGcrkZonWGGPyOHHiBD179uS33zJ4/fXXufvuu9m5cycPPPAA48ePL/L91uvY5BbpdADGGBNsBg4cyLfffgt8wMMP3wW4km+7du1ISUmhTZs21KtXr8D3L1kSoEBNSLAarTHG5DJq1CjeeustBg0aBNx1envZsmX59NNPOXXqFH369OHUqVPOBWlCiiVaY4zJkZWVxeOPP87VV1/Nyy+/fM7rDRs2ZNiwYcydOzff191q1fJnlCbUWKI1xpgcM2bM4I8//mDQoEFERETkO5Vinz596NWrF6mpqWzcuDHf42zf7udATUixRGuMMTnGjBlDXFwc119/PZB/pyYR4dVXXwXgzTffzPc4gwf7K0ITiizRGmMMcOjQISZPnszNN99M2bJlgYKnUqxTpw633HIL7733HgcPHjzn9eee82ekJtRYojXGGGDKlCkcOXKE2267zaP9H3nkETIzM3n//ff9HJkJdQFNtCLSSUTWiUi6iDyRz+v1RORbEVkqIitE5MZAxmeMyV9pKLtjxoyhbt26/PnPf/Zo/8svv5y2bdvyxhtvWA9kU6iAJVoRiQCGATcATYHeItI0z25PA+NVtQXQC3g7UPEZY/JXGsru3r17mT59Or1796ZMmTNfi0VNpThgwAA2bdrElDy9phYt8keUJlQFskbbCkhX1Y2qmgWMBbrl2UeBmJyfYwHru2eM88K+7H7xxRecPHmS3r17n7W9qKkUu3XrRv369Xnttdf8GJ0JdYFMtHWArbmeZ+Rsy20wcIeIZADTgAfzO5CIpIjIIhFZtHv3bn/Eaow5I+zL7meffcZFF11E8+bNz9pe1FSKkZGRPPTQQ8ybN49ly5ad3t6ypT+iNKEqkIlW8tmmeZ73Bj5Q1brAjcDHInJOjKo6UlVbqmrL6tWr+yHU4JeQAJLzFx050vWz+5GW5hrHJ/n9xY0pvrAuuwcPHmTu3Ln07NkTyVNoPJlK8a677iIyMpJPP/3UTxGaUBfIRJsBnJ/reV3ObV76KzAeQFUXAOWB+IBEF2K2bAHN+apLSXH97H4kJ0Pt2s7GZ8JKWJfdn376CVWlXbt2Xr0/Li6OTp06MXbsWLKzs30cnQkHgUy0PwONRaSBiETh6jCRd96V34BrAUTkYlyFNXjal4JIly5OR2BKkbAuuwsXLkREaNWq1TmveTqVYu/evcnIyOD7778HIDXVlxGaUBewRKuqJ4H+wAxgDa4eir+IyPMi0jVnt0FAXxFZDowB7lbVvE1UBlfzcFFs8WnjC+FedhcsWEDTpk2JjY095zVPp1Ls2rUr0dHRfPbZZ4DNDGXOJiFSFgrUsmVLXVQK+9InJ3uWbE3wE5HFqlrqus8EQ9lVVeLj4+nevTvvvvvuOa8PHux50rztttuYMWMGO3bsICEhyuY7LgU8Lbs2M1SImjq16H1s8WljCrd+/Xr27dvHFVdcke/rxZlKsXfv3uzbt4+ZM2eyY4ePAjRhwRJtGLPFp40p3IIFCwBo06ZNiY/VsWNHqlatypgxY0p8LBNeLNEaY0qthQsXEhMTw8UXX1ziY0VFRXHzzTczadIkmjc/4oPoTLiwRBuiPLm1botPG1O4BQsW0Lp167OmXcytuLeQb7vtNg4fPswTT1gHCnOGJdoQVdTUcGCLTxtTmEOHDrFy5coC7896o127dlSvXp2nnprks2Oa0GeJNkTdf3/R+9gQA2MK9vPPP5OdnV3o/dniTqUYERFB165d2bDhK44fP17CCE24sEQbxmzxaWMKtnDhQgBat27t0+PedNNNQCbfffedT49rQpcl2jCXmZnJypUrOXHihNOhGOOo3HOAi8CTTy4gNvZCqlWrRlLSme3u6UsHD4b69Yt/nmuvvRaoyKRJ1nxsXCzRhqgpeSfAy5GRkcHTTz9N165dgQbExMRw2WWXkZCQQGpqKhkZGQGN05hg0bWrK4mqQna2Ur36Qrp3d92fXbz4zFzh7r4NgwfD5s3FP0+FChW48cZOTJ482eY+NoAl2pCV32QUy5cvp1WrVgwdOpQNGzbQoUNrXnzxRd5//30uu+wyXnjhBRISEujevTvLly8PfNDGBIktW7awe/dunzcbu110UXd27NjBzz//7Jfjm9AS6XQAxjt16pw9xGfWrFn8v//3/4iJiWHp0qU0a9aMxYvPJOR77rmHjRs3MmLECN577z3atm3LxIkTuf766535BYxx0KpVqwC47LLL/HL8//znRiIjI5k0aZLfkrkJHVajDQOffPIJN9xwA/Xr12fhwoU0a9YMOLfH5AUXXMDLL7/MypUradiwIZ07dz49Cbox4W7EiDM/r169GoCmTZv66WxVueqqq+w+rQEs0Ya8YcOG0adPH9q1a8e8efOoW7duke+pVasWc+fOpW3bttx+++289tprAYjUGGelpJz5+ZdffqF27dpUqVLFb+e76aabWLt2LWvXrvXbOUxosEQbovr2ha1bt/Loo49y44038vXXXxfrSyM2Npbp06fTo0cPBg4cyNNPP+3HaI1xnsiZn3/55RcuueQSv51rxAjo1q0bgNVqjSXaUDVyJDz55JOoKm+//TblypU7Z5+iFp8uX74848aN469//StDhgxh7NixforWmOCRnZ3NmjVr/JpoU1Kgbt26JCYmMtWTpbZMWLNEG6IuuuhnPvnkEwYOHEj9Agb7eTIzVEREBMOHD6dt27bcd9991sxlwt6WLVs4cuSIH+/Pnqk9d+nShQULFrBnzx6/ncsEP0u0IUhVWbduIDVq1OCJJ54ocD/3wPuilC1blrFjx1KhQgV69uzJkSO28ogJP126uP795ZdfAPxao3VLTk4mOzubadOm+f1cJnhZog1BEydOBL7nhRdeICYmpsD9irP4dN26dfn000/55Zdf6N+/f8mDNCbIpOUsqONOtP6s0bolJiZSq1Yt0tJsNZ/SzBJtiDl+/Dj/+Mc/iIy8lHvvvdenx77++ut5+umnGT16NKNHj/bpsY1xWnKy699A9Dh2157LlClDly5dmDFjBllZWX47nwlulmhDzLBhw9i4cSNfffVvIiMLn28kMbH4x09NTeWaa67hb3/7G+vXr/cySmOCj7tP0urVq/3ebJy7AtulSxcyMzOZO3euX89pgpcl2hBy9OhRXnjhBTp16sT8+UXP6LR4cfHPERERwSeffEJUVBT9+/dHPVlh3pgQEYgex3Cm9gxw3XXXUb58eWs+LsWKnWhFpKKIRPgjGFO4r7/+mv379zNo0CCPlsDLPUC/OGrVqsWLL77IN998w4QJE7w7iDFBaPPmzX7vcQxnas8A0dHRXHvttaSlpdmFaylVZKIVkTIicpuIfCUiu4C1wA4R+UVEXhGRxv4P0wCMHTuWGjVqcNVVV3m0/6hR3p/rgQceoHnz5gwYMIDMzEzvD2SCRmm/SFY9M/ViIHoc55acnMymTZtYs2ZNQM9rgoMnNdpvgYbAP4HzVPV8Va0BtAMWAkNF5A5PTiYinURknYiki0i+41JE5BYRWZ2TyG0i3hyZmZlMnTqVnj17Fnlv1hciIyMZPnw427Zt4zlbQT4k+fIiORzK7siRge1xnFvnzp0BrPm4tFLVQh9AWR/tEwFsAC4AooDlQNM8+zQGlgJVc57XKOq4SUlJWhp8+umnCui8efNUVXXRoqLfAyU/73333acRERG6YsWKkh/M5AtYpEV8zr15AHOAZ4DLgDK5tlcDegBfAHd4cJywKLug2qdPH61Tp05Az+vWokULbdu2rSPnNv7hadn1pEb7XxF5QETaiki+gzZV9YQHx2kFpKvqRlXNAsYC3fLs0xcYpqp/5Bx3lwfHLRXGjBlD3bp1ufLKKz1+z7ZtJT/v0KFDqVKlCn/729/s/lLouU5VX1DVFap6egVyVd2nql+oag9gnAfHCZuy+8svvwSkNjty5LnbkpOTbZaoUsqTRLsMuAT4F7BJRDaLSJqIDBGRXsU4Vx1ga67nGTnbcmsCNBGRH0RkoYh0Ksbxw9a+ffuYMWMGt956K2XKuP7L8i6Blx9veh3nFRcXx8svv8z333/PJ598UvIDmkDy1UVymJTdwPQ4Brj//nO32SxRpVeRiVZVR6pqf1Vtr6pxwJ+B4cBBoHMxziX5bMtbRYrE1QR1FdAbeFdEzhlVLiIpIrJIRBbt3r27GCGEpi+//JITJ07Qq1dxrmuga1ffnP+ee+4hMTGRZ555huPHj/vmoCYQfHWRHBZld+TIzRw9ejTgHaHcbJao0qvYw3tUNUNVp6nqy6rapxhvzQDOz/W8LrA9n30mq+oJVd0ErMNVePPGMFJVW6pqy+rVqxf3Vwg5Y8eOpVGjRiQlJTly/jJlyjB06FC2bNnC8OHDHYnBFJ8PL5LDouyWK+fq8XvxxRcH9LxuNktU6RXICSt+BhqLSAMRiQJ6AVPy7DMJuBpAROJxNUdtDGCMQef333/nf//7H7169UJyLahZ1BJ4vtahQweuvfZahgwZwsGDBwN7cuMTJbhIDouye9ddrpnOmjRp4vdzTcn718mRnJxMZmYmc+bM8XsMJngELNGq6kmgPzADWAOMV9VfROR5EXE3cs4A9orIalzDih5T1b2BijEYTZgwgezs7HOajT1ZAm/ECN/GMnToUPbs2cOrr77q2wOboBY+ZTedmJgY4uPj/X6mghqfrr32WpslqhQSb3uSikgtYJ+qOnrTrmXLlrpo0SInQ/Crdu3asX//flauXHnW9tq1YXvexrsAuOWWW5g2bRobNmygZs2agQ8gDInIYlX1oHubz85XKsuuSEeSkvYSiHOKuCbIyE9ycjKrVq1i48aNZ7VSmdDjadktSY32Y2CtiFj1xk+2bt3K999/n28nKE+WwPNHGX7xxRc5duwYL7zwgu8PbgKlVJbdypXX07ix8xPZJScns3nz5tOTZ5jw53GiFZHXJdfll6peh2sAu62n5iczZswAoHv37g5HckaTJk247777GDFiBBs3BtUtOOOh0lh2s7KyOHx4C40aNXI6lNOzRE3NPSGyCWvFqdEeAqaISEUAEekAfK+qdlnmJ9999x01a9bMt5ekN0vg+cqzzz5L2bJlefbZZ50Lwngs70UynJ4zrNSU3U2bNpGdnR2wGm3fvgW/VqdOHZKSkphSUI8pE3Y8TrSq+jQwBvhORL4HHgXynfPUlJyqMmfOHNq3b5/vfRxPJqNwLz7ta7Vr1+bBBx/ks88+s0nSQ0Pei+TrReQHh2MKqPT0dICAJdr8ZobKrVu3bixcuJCdO3cGJB7jrOI0HV+La5q1w0B14CFVneevwEq7jRs3kpGRQfv27fN93ZMl8PzZsfGxxx4jOjra7tWGgHwukgdRyi6S1693De0JVNNxUUPeb7rpJlTVeh+XEsVpOn4KeEZVrwJuBsaJyDV+icrw3XffARS4JJ4nS+DlXnza1+Lj43nwwQcZO3bs6aXHTHCyi2RXohWJDcjQHoAlSwp//dJLL+WCCy5g0qRJAYnHOKs4TcfXqOr3OT+vBG4AXvRXYKXdnDlzqF69eolmsfF3X4tBgwZRsWJFnn/+ef+eyJTUU8CzpfkiOT09naSkxkEznEZEuOmmm5g1a5at91wKeD28R1V3ANf6MBaTQ1X57rvvCrw/Gyzctdrx48fbUIUglnORPC/n51J5kbx+/XqOHQtcj+NatYre56abbiIrK4vp06f7PyDjqCITrYjUK+gBVM/1PN/VQUzxbd68ma1btxbYbAy+WQLPFwYNGkSlSpVscfggVEi5LQv8tbSU3aysLLZs2cKqVYEbQ+vJZDJXXnkl8fHx1nxcCkR6sM+HuFbqKKxqpcAHwEc+iKnUK+r+LLh6HdeuXfhxArF8bFxcHA899BBDhgxh1apVXHrppf4/qfFUYWXXvT3sy657aE8+axz4zeDBRU+TGhERQXJyMhMnTiQrK4uoqKhAhGac4Mnq8MH8SEpK0nBz5513anx8vGZnZxe4z+mRkIUYMcKHQRVi7969WrlyZb355psDc8IwAyxSP5YRoLY/j+/tI1BlNy0tTQGF+QE5n6pn5VNVdfLkyQroN99849+AjF94WnaLM7xnpYh8KiKPi8gNIlJXRJ7yT/ov3QobP1sc+S0+7Q/VqlXj4YcfZsKECefMyWyCwjQRSRWR8k4H4gT30J5Zs5yffjGvDh06EB0dbc3HYa44naHaA6OAo7iWyVpF8da0NB7YvHkzW7ZsKbTZOBgNGDCAypUr27ja4NQS1xq0P4vInU4HE2jp6elUqVKF2Ng4p0M5R4UKFejYsSOTJk3Kad424ag4w3v2qep3qvqGqt4FXA6s919opZP7/mxBE1W4ebIEXv36rkkrtm93LTDgfrgnu0hKgoSEksXrVq1aNR566CE+//xzVq1a5ZuDGp9Q1ZOq+hrwFyBJROaLSDun4wqU9evX06hRIy6/PHA9+IuzQFCPHj3Yvn07Cxcu9F9AxlHFaTo+q91FVdcDl/k8olJuzpw5xMXFcckllxS6nyczQ23e7Jq0onZtV8co98M9PdzixbBli2dxjRx5drLOncDdydpqtcFJRC4QkQeAVOBCoBEwWkS2iEjYr0Cenp4eFKv2FCQ5OZly5crx+eefOx2K8ZPiNB2PFJHfRGSBiIwQkQ+BVSIS7a/gSiP3+NkyZQr/r/HV8FpPxvu5a8G5k3XuBO5O1nFxcTz44IN8/vnnNltUcJkNxOb8+zCuzlGNVLU+ENZNye6hPYFOtC2LsbpwTEwMHTt25PPPP7fm4zBVnKbjq1W1HnAr8BWQDkQDy0VkrZ/iK1W2bNnC5s2bi2w29iVfLB6fe6WSgQMHUrFiRavVBpcOqjpUVdNUdZ2qnnS/oKoetmmEpo0bN5KdnU2jRo1ITXU6moL17NmTbdu2WfNxmPJkwoq8y2v9pqpTVPUFVe2hqo2BVn6LsBSZM8fVihfIjlBFjfXzRO6VSty12nHjxlmt1mHusquq6UXtE65yr9rji8+6vyQnJxMVFWXNx2HKkxrttyLyYM6MMqeJSJSIXJPThPz//BNe6TJv3jyqVq3q0aQPvloCz5MJnYo6V96VSgYOHEh0dDQvvliqZvkLRp6U3bscii0g3Im2UaNGRU7w4kvFrT3HxsbSsWNHJkyYYM3HYciTRNsJOAWMEZEdIrJaRDbh6nHcG3hNVT/wY4ylxvz582nTpk2R92fBv0vgFfdceVcqiY+Pp3///owdO9bWq3VWqS+76enpxMbGEhcXx44dgTuvN7Xnnj17kpGRwY8//ujzeIyzivxGV9Vjqvq2qrYF6uFaSKCFqtZX1b6quszvUZYC+/fvZ/Xq1VxxxRUe7e/PJfB8ca5HH33UVvZxmJVd2LBhA40aNQr44hze1J67du1qzcdhqjjDe24A5gHf4eqB3MZfQZVG7qtYTxOtr5bA82S8X1Hnyq/ncnx8PA899BDjxo2zcbUOK81lNz09nYYNGwKQmBi483pTe7bm4/BVnOE9bwODgDbASOBVEentl6hKoQULFiAitGoVev3KCuq5bCv7BI2wKrsJCWfGc0NhY7xPkp6+mfHjG5KS4ho3HuxuueUWtm7dyvz5850OxfhQcRLt76r6g6r+oaqzgI64FpQ2PrBgwQIuvfRSYmICu2JZccb7FaSg+1HVqlVjwIABTJgwgeXLl5f8RMZbYVV2mzU7M54bCh7jvWHDb8BJ3nuv0Vk94wPB29rzTTfdRHR0NJ988olvAzKOKk6i3SwiL4qIey2nE0BmcU4mIp1EZJ2IpIvIE4Xsd7OIqIj4IA0Ev+zsbH788UePm40hMEvgeXquwiqsAwYMIDY2lsHBPLYi/IVV2fW0I+CGDRsATjcdB5K3tedKlSrRvXt3xo8fz/Hjx30blHFMcRKt4hrGs1VEvsc1YcV3eadmLIiIRADDgBuApkBvEWmaz36VgYeAUtP1bs2aNRw4cKBYiTaQV+glOVeVKlUYNGgQkyZNYnEotN2Fp7Aqu552zss9tCfQPJkitSB33HEHf/zxB9OmTfNdQMZRxZkZqreqNgXqA48AzwEVgXdFZKsHh2gFpKvqRlXNAsYC3fLZ7wXg/4BjnsYW6hYsWADAlVde6fF7fLUEnifj/Up6rocffpiqVauSGsxT84SxcCu7nnYE3LBhA+XLl6eWJ/OM+tg335y5Z7x9u6sWnvs+svviNfc29wXEddddR82aNa35OIxEFvcNqnoMWJTzKI46QO5CnQG0zr2DiLQAzlfVqSLyaHFjC1ULFiwgLi7OkYnPfdGiW1TP5ZiYGB577DGefPJJfvzxR1q3bl34G4xflLay6+5x7Mm4dF/bvPns5+55wfPKuy05GdLSIunduzdvv/02f/zxB1WrVvVbnCYwAvkJzG8g2+mPmYiUAV7D1Tuy8AOJpIjIIhFZtHv3bh+G6IwFCxbQpk2bgI/1A+/G+3mjf//+xMfH8/TTTwfmhMaXQrLsbtiwwZH7syXhrq3fcccdZGVlMWHCBGcDMj4RyESbAZyf63ldIPfAkMrApbjuHW3GNRRhSn6dKlR1pKq2VNWW1atX92PI/vfHH3+wZs2aYt2fBZgyxTfnj4pyNVu5p1FMSTm7OWv79qLP5UnP5cqVK/PUU08xa9YsZs2aVfLATSAFVdn1pCOgqp6erCIUJSYmctFFF/Hxxx87HYrxgUAm2p+BxiLSIKf3Yy/g9Fe4qh5Q1XhVTVDVBGAh0FVVi9vMFVKKO1GFW975hb21ebPri8vdT2nkyLOHStSu7btZqPr160e9evX45z//iQay27QpqaAqu550ztuxYwdHjx4NuRqtm4jQp08f5s2bx6ZNm5wOx5RQwBJtztJc/YEZwBpgvKr+IiLPi0jXQMURbBYsWECZMmWKPVFFnTp+CsiPypcvz3PPPceiRYuYOHGi0+EYDwVb2fWkc56TQ3tKIvf15x133IGIMHr0aOcCMj4R0F4CqjpNVZuoakNVHZKz7VlVPadxUlWvCvfaLLgSbbNmzahUqZLToXitOJ2J+/TpQ9OmTXnqqac4efJk0W8wQSHUyq6TQ3tKIndtvV69enTq1In33nvPykqIC3x3PHOaNxNVBKPBg11NzO77ugXd7wWIiIhgyJAhrFu3jg8//NCxmE1427BhAxEREdSrV6/onYNI3tp6SkoK27dvtzG1Ic4SrYNWr17NwYMHvUq0ffv6IaAS2L79zH3d/O73btt2Zt9u3brRunVrBg8ezNGjR50J2IQsTzoCpqenk5CQQNmyZf0fkB917tyZ8847j1GjRjkdiikBS7QOck8c7k2iDfTcrSWVe1IoEWHo0KFkZGQwbNgw54IyIcmTjoChOLQnP2XLluXee+9l2rRpZGRkOB2O8ZIlWgctWLCA+Ph4r+4j+arXcaB0zdNl5qqrrqJjx47861//Yt++fc4EZUKSJx0B09PTQ+7+LORfW7/vvvvIzs7m/fffD3xAxics0Tpo/vz5Xk9UsWSJHwIKsFdeeYUDBw7Y4vDGp/bt28f+/ftDskab3wV0gwYN6NChA++++y6nTp0KfFCmxCzROmTHjh38+uuvtG/f3ulQHNOsWTPuu+8+hg0bxrp165wOx4SJUO1xDAXX1lNSUti6dStff/11YAMyPmGJ1iFz584F8DrROjBPeomMGJH/9ueff54KFSrw2GOPBTYgE7KK6ggYqmNoC9O1a1dq167NfxJsLNEAACAASURBVP/7X6dDMV6wROuQOXPmULlyZVq0aOHV+7dvL3qfYFLQsmE1a9bkqaeeIi0tjdmzZwc2KBOSiuoI6K7RXnDBBQGIJjCioqLo378/s2bNYuXKlU6HY4rJEq1D5syZQ9u2bYmMLPYCSoBvVt0JpMJuQz/88MMkJCQwcOBAuwdlilRUR8Bff/2VunXrUqFChcAE5EOF1dbvv/9+KlSowOuvvx64gIxPWKJ1wO7du1m9enWJ7s8+95wPA3JY+fLlefnll1mxYoVNN2eKVFRHwHXr1nHhhRcGJhgfK6y2Xq1aNe666y4+/fRTdu3aFbigTIlZonVASe/PhqOePXty5ZVX8vTTT3PgwAGnwzEhSlVDOtEWVVt/5JFHOH78OMOHDw9MQMYnLNE6YM6cOURHR9PSk/XlwkSXLoW/LiL897//ZdeuXaQWZ/JkU+oU1hHw999/5+DBgyGbaIuqrV944YV07tyZt99+m2PHjgUmKFNilmgdMGfOHK688soSTQ+3KMSWW0hLK3qfli1b0q9fP958802WLl3q/6BMSCqsI6B7mFioJlpPDBgwgF27dvHpp586HYrxkCXaANu3bx8rV64sdc3GycnnLjAwcuTZ29LSYMiQIcTHx/PAAw+QnZ3tbNAmKBXWETDUE60nw/auueYaEhMTeemll2xVnxBhiTbA5s2bh6qWONGGWqtzWtrZC8qDa8hP7m1JSVC1alVeffVVfvzxR959911ngzZBqbCOgOvWraN8+fIht2qPmyfD9kSEZ599lg0bNlitNkRYog2wOXPmUL58+WIv9F4auGfFueOOO2jfvj1PPPEEu3fvdjYoE1LWrVtH48aNKVMmNL/aPB2217VrV5o3b86LL75otdoQEJqfxhA2Z84c2rRpQ7ly5ZwOJWiJCG+//TaZmZk8/vjjTodjQkgo9zgGz4ftiQipqamkp6czZswY/wZlSswSbQAdOHCAZcuW+eT+bLh3zG3atCmDBg1i9OjRp4dDGQMFdwTMyspi06ZNIZ1oi6Nbt240b96cF154wWq1Qc4SbQB9//33ZGdn+yTRhtrMUJ7IOyvOM888Q4MGDbj33ns5fPiwM0GZkLFhwwZOnTpVahKt+17t+vXrrVYb5CzRBtB3331HVFQUbdq0KfGxatf2QUBBJu+sOBUrVmT06NFs2LCBf/7zn84EZYJOQR0BQ73HMRR/2J67Vvvss8/auNogZok2gObMmUOrVq18Mgfrjh0+CCjI5DcrTvv27XnwwQd58803mTNnTuCDMiEjHBJtcZUpU4ZXX32VzZs388YbbzgdjimAJdoAOXjwIEuWLCl142eLo6BZcV566SUaNmzIPffcw6FDhwIblAkZ69ato2bNmsTGxjodite8GbZ37bXXkpyczIsvvmhzIAcpS7QB8tVXX3Hq1Ck6derkk+MlJvrkMCHB3YS8efNm64VsCuwIGOo9jkvilVde4ejRozz77LNOh2LyYYk2QL744gvOO+88rrzySp8cb/FinxwmqBQ2K067du14+OGHefvtt/nf//4XuKBM0CmoI2BpTrQXXnghDzzwAKNGjWLVqlVOh2PyCGiiFZFOIrJORNJF5Il8Xh8oIqtFZIWIzBaR+oGMz1+OHDnC119/Tffu3X02kL6ghdRDWVGz4gwZMoQmTZrQp08fm8giwIKp7ObXEXDv3r3s3bs35BNtSYbtpaamEhMTwyOPPIK6p18zQSFgiVZEIoBhwA1AU6C3iDTNs9tSoKWqXgZMAP4vUPH504wZMzhy5Ag9evTw2TFHjfLZoYJGUUOWoqOjGTt2LHv27OHOO++0uZADJNjKbn4dAcOlI1RJhu3FxcUxZMgQZs+ezUcffeSzmEzJBbJG2wpIV9WNqpoFjAW65d5BVb9V1SM5TxcCdQMYn9988cUXVKtWjb/85S9OhxLUPJkVp0WLFrz++utMnz6dV155xf9BGQiBshsuibakw/b69etH27ZtGTBgAL///rtvgjIlFshEWwfYmut5Rs62gvwV+NqvEQVAVlYWaWlpdOvWrUTL4pkz+vXrR8+ePXnqqaf44YcfnA6nNAiqsptfR8B169ZRtmxZGjRo4K/TBkRJh+2VKVOGd999l8OHD/PQQw/5JihTYoFMtJLPtnxvJIjIHUBLIN8qi4ikiMgiEVkU7PfqZs+ezcGDB33abAywbZtPDxdSRIRRo0ZRv359evXqxd69e50OKdwFVdnNryPg2rVradiwIZGRkV4dM5xcdNFFPPPMM4wfP54pU6Y4HY4hsIk2Azg/1/O6wDndX0TkOuApoKuqHs/vQKo6UlVbqmrL6tWr+yVYX/niiy+oXLky1113nU+PG469joszK05sbCzjx49n165d3HnnnZw6dcp/gZmgKrv5dQRcunQpzZs39+p4wcRXw/b+8Y9/0KxZMx544AH27Nnjm4MarwUy0f4MNBaRBiISBfQCzrrcEpEWwAhcBTXkR16fPHmSSZMm0aVLF5+v1tO1q08PF5KSkpL473//y7Rp02x8rX8FVdnN2xFwz549/PbbbySGweByX11AR0VF8dFHH7Fnzx7uvvtu64XssIAlWlU9CfQHZgBrgPGq+ouIPC8i7rTxClAJ+FxElolISLd7zJs3j7179/q82ThceTMrTr9+/XjwwQf597//zahw7IodBIK97C5duhQgLBKtL4ftNW/enH//+9989dVXvPbaa747sCk+VQ3pR1JSkgarv//971qhQgU9dOiQz48NPj+k47z9nU6cOKE33HCDRkZG6qxZs3wbVAAAizQIylKgH96W3byfk6FDhyqge/fu9ep4wcTX5To7O1u7d++uZcuW1Z9++sm3Bzcel12bGcpPsrOzmThxIp06daJixYo+P/6IET4/pOPq1z8zjrB2bRBxPdyLDaSkQELCue+LjIxk7NixXHTRRfTo0YO1a9cGKmRTAitXuv5/k5Ndz5OTz/yfS073q5Ejz96WlnZuR8AlS5aQkJBAtWrVAvsLhAAR4b333qN27drceuut1nHQIaIh3nbfsmVLXVTctaUCYP78+bRt25ZPPvmE22+/3elwwoYIFPSR3bJlC61ataJSpUrMnz+fmjVrBjY4L4nIYlX1ouE8tIm0VNWSl90mTZrQrFkzvvjiCx9E5azCPt8l8eOPP9K+fXtat27NzJkziYqK8v1JSiFPy67VaP3k9ddfJyYmhmT35Xox5b6KL+yK35xRv359pkyZws6dO7n22mttmsYgl9+yiMV18OBB1q9fHxb3Z8F/w/Zat27N6NGjmTt3LikpKdY5KsAs0frBunXrmDBhAn//+9+JiYkp9vtHjnRd1bofaWmu7WlpZ28vjYr6ImrdujVTp05lw4YNdOjQgX379gUmMFNsvrgOWrZsGRAeHaHAv8P2evfuzXPPPceHH37I0KFD/Xcicw5LtH4wdOhQypcvzyOPPOLV+++/38cBhRFPvoiuvvpqJk+ezJo1a7j++uvZv3+//wMzxfbbbyU/xpKcRYzDJdH6e9jeM888w+23386TTz7J+++/79+TmdMs0frYli1b+OSTT+jbty81atRwOpyw4+kX0fXXX8/EiRNZsWIFnTp14uDBg/4NzDhiyZIl1K5dO2TuxzvN3TmqY8eO3Hfffbb4QIBYovWx//u//0NEeOyxx5wOpVRLSYEuXTpz4sR4fvxxMW3bXsWHH+44fX87v97LJvQsWbKEFi1aOB1GSClXrhxffvkl11xzDffccw+fffaZ0yGFPUu0PrRz507ee+897rrrLurW9X7xEpuetGSSknLf576JadOmsGnTr6SmXsHq1WtQhSefdDpK06hRyd5/5MgR1qxZEzbNxhC4YXsVKlRgypQptGvXjj59+vDBBx8E5sSllCVaH/rPf/7DiRMnSjwdoC96Y4YrT76Icm7bnXbDDTcwZ84cjh07Rtu2bZk3b55PZ+Ax3omOLtn7V6xYQXZ2dlgl2kB+LqOjo5k6derpmu3zzz9vvZH9xBKtj+zbt4/hw4dz66230qiEl+p1CluArJTz9osoKSmJBQsWUKNGDTp06IDIGN8GZoptxYqSvT/cOkJB4IftVapUia+++oo777yT1NRU+vbty4kTJwIbRClgidZHXn31VQ4dOsST1ibpV558EdWqlf/2Bg0a8MMPP3D55ZcDt9G/f3+OH893kRkTApYsWUJcXBznn39+0TuHuaSkM+Pr3YvHDx589rj7xYtdj9zbBg92LUDwwQcf8PTTT/Pee+9x3XXXsX37OYszmZLwZJ7GYH4Ew1zH3333nZYpU0b79Onjk+OF4zzGvuKLv83x48cVBimgiYmJmp6eXvKDlgCldK5jKFnZbdasmXbo0KFExwg2Tpf9jz/+WKOjo7VGjRohOW94oHladq1GW0K7d+/mtttuo1GjRgwbNswnx+zb1yeHKbXc8yUXJCoqii5dXmXy5Mls3LiRxMREPv/884DEZs6Ij/f+vVu3bmXlypU+X+fZaV26ePc+b2+puGu/bnfccQc//fQTcXFxdOjQgaefftpafXzBk2wczA8na7SnTp3Sjh07arly5XTZsmWOxVGadOlS9D7FqRVs2rRJW7VqpYDefPPNun37du+D8xKltEZbkrL71ltvKaBr1671+hjBqksX1dxzwKmqjhhx9rYpU1S3bTvzvH59785VUFnJzMzUu+++WwFt2rSpLly40LsThDlPy67jha2kDycT7UsvvaSADh8+3KfHTUz06eFKHU8Sbe6EnZWVpUOGDNFy5cppbGysjhgxQk+dOuW/APMorYk2Otr7snv99ddrkyZNvH6/cSmqrHz11Vdat25dFRF96KGHwmIpQl+yROtn8+bN04iICL3llls0Ozvbp8d2+j5NMPNVjTa/fX799Ve9+uqrFdB27drp4sWLix+gF0provX2Hu2BAwe0bNmy+uijj3r1fnOGJxf1Bw4c0H79+mmZMmW0atWq+vrrr+vx48f9H1wI8LTs2j1aL/z888/06NGDhIQERo0ahdhSOgEzdWrR+3i7amLjxo2ZPXs277//PqtXryYpKYlevXqxfv167w5o/GL69OmcOHGCbt26OR1KyPNk7vCYmBiGDx/O0qVLSUpK4pFHHqFp06a8++67ZGVl+T/IMGCJtpgmTpxI+/btTw/29mZ1nqIUNDzF+J+IcM8997Bhwwaefvpp0tLSaNq0KQ888ABbt251OrywUrasd++bMmUKcXFxXHHFFb4NqBQqTieqyy67jG+++YavvvqKqlWr0rdvXxo1asQbb7xhc4kXxZNqbzA/AtV0nJ2dra+88oqKiLZu3Vp37twZkPOas3nbLOytHTt26N/+9jeNjIzUiIgI7dWrl/7000++O4GW3qZjb8puVlaWVqlSRe+6665iv9ecy9uykp2drdOnT9d27dopoJUqVdJ+/frp8uXLC3xP/fqqtWq5fk5N1bM6dy1a5Hrk3paa6trX/Z5g5GnZdbywlfQRiER76NAhTUlJUUB79uypR44c8ev53B8w4x1PvjxGjCjeMTdt2qQDBw7UmJgYBbRt27Y6btw4PXr0qHdB5lJaE22tWsUvu//73/8U0IkTJxb7veZcJbkoPdMT+ieFu7Vs2XIKKLRQeFUhQ/v2de2bmOj7ntHBwBKtD5w8eVLff/99rV27tgL6xBNPBKQ3ajB/sJzmSYL0Z633wIED+tprr2lCQoICGhsbqykpKTpv3jyvO8WV1kTrTWeoRx55RMuVK6eZmZnFfq85l7flIL/37dmzR19//XW9/PLLFVAR0SuuuEKHDh2qq1evLkH58C7GQLBEW0IzZ87UP/3pTwpo69at9fvvv/fLefITzB8sp3nyt/GkRaCkf+OTJ0/qzJkztU+fPhodHa2AJiQk6EMPPaQzZ84sVq9MS7SeycrK0vr16+uNN95YrPeZgm3b5hqTm7vJ1n0xm3ubu7e/e4xvUbXTdevW6XPPPaeJiYk5tVy0Xr16es899+gnn3yiGRkZHscYzMMdPS274to3dLVs2VIXedvNNI+dO3fy2Wef8fHHH7Ns2TIaNGjASy+9xC233BLQnsUiro+3OZev/ja+/BsfOnSIiRMnMm7cOGbPns3x48epXLkyHTt25JprruEvf/kLF198MWXK5N/3UEQWq2pL30QTOkRaqqrnZffFF1/kmWeeYfLkyXTt2tWPkRlfysjIYOrUqcycOZNvv/2WP/74A4B69epx5ZVX0rp1axITE2nevLlfOpf6k6dlt1QnWlUlPT2dOXPm8MUXX/DNN9+QnZ3N5Zdfzt13381f//pXypUr5+OIi7Z4sS2VVxBPEmTt2lDUnOhpaZCc7Lu43A4fPszs2bOZOnUq06ZNY9u2bQBUq1aNdu3a0aZNG1q0aEGLFi2oUaMGUHoTbdOmLXX1as/K7sqVK0lKSqJHjx6MGWMrL4WqU6dOsXTpUn744Qfmz5/P/PnzycjIOP36BRdcwCWXXELTpk1p2rQpjRs3Zvjwxnz4YVxQDqMMykQrIp2A/wIRwLuqOjTP6+WAj4AkYC9wq6puLuyYniZaVWXXrl2sXr2aVatW8cMPPzB37lx27NgBuK6u7rjjDvr06cNFF13kza/nM5ZoC5aW5vrb5F5KsG9f10LvSUmutWjr14fNmws/zvbt587z6muqysaNG5k3bx5z585l3rx5pKenn369du3a/OlPf+Lrr78O+kTrj7LraaI9ceIEbdq0YevWraxevZr4kkySbILOzp07Wbp0KUuWLGH58uWsXr2aX3/9Nc9yfVWABDp2rE9UVAJpaXWBukAdXnihFo88ch6VK1c6vXeXLq7vCn8LukQrIhHAr0AHIAP4Geitqqtz7fM34DJV7ScivYDuqnprYcd1J9ojR46we/dudu/ezc6dO9m6devpx+bNm1mzZg179+49/b46derQvn17/vKXv9C+fXsuvPDCoLlisqZj/8v9X503UYNrLLM/Vgrbv38/y5YtY8mSJSxdupSVK1eyfPnyoE60/iq7njYd/+tf/+Kpp55iwoQJ9OjRoyS/igkRJ06cYMOGDaSnp7N+/XrS09PZvHkzW7ZsYfPmzRw+fPic90RHR1OjRg3i4+OpXr06K1bE0aNHNapWrUrVqlWJjY2lSpUqxMbGUrlyZSpXrkxMTAwVK1akYsWKREREFDvOYEy0VwCDVbVjzvN/AqjqS7n2mZGzzwIRiQR2AtW1kCCjoqI0MjKSo0ePnvNaZGQkderUoV69elx88cVnNUnUqlUraBJrXpZoS5dgbzr2V9ktLNHu27ePL7/88vR975tvvplx48b58LcyoUpVOXjwIBkZGWzbto2dO3eefuzevZs9e/awe/duli/fx8mT+4ADHh65POXKRVO9ejTR0dFUqFCB8uXLn/63fPnylCtX7qzHW2+95VHZjSzZr1wsdYDcU+tkAK0L2kdVT4rIASAO2FPQQWNiYrj77rtPX8XEx8dTs2ZNzj//fM477zyvrlKMGTy46OX2PLF9u+tWQO6+OyNGuGbkCdLrvPz4pezmNnr0aPr160d2djanTp3CnZ8bNmzIP//5Tx577LGS/g4mTIgIsbGxxMbGcskllxS5/6lTp9i/fz8HDhw4/cjMzDz9OHz4MIcOHeLQoUMcPXqUMWMOk5FxBDhKVNQxWrc+yooVe9m+/ThwFMgiLu44J054vnxgIBNtfl8rea92PdkHEUkBUgBiY+vx73+/evo19+3aunXP7J+a6vrSrF0bcm7Jkpjo+gJMSYFRo87su22bZ1+M7nsAyclnz7+r6mqGvP/+M9umTCn6viKcaa5MTc3nr2AC6oMP4LnnXD+7P1Mtc123Fucz5R4kkZd7WwgkXL+VXbdLL72UgQMHUqZMGSIiIqhQoQLXX389iYmJQdvyZEJDREQEcXFxxMXFebT/m296fmxPP5sh33Tsy+E9xjihtDYdW9k1oc7TshvIRQV+BhqLSAMRiQJ6AVPy7DMFuCvn55uB/xVWUI0xAWFl15gSCFjTcc59m/7ADFxDBN5X1V9E5Hlcs2tMAd4DPhaRdGAfrgJtjHGQlV1jSiaQ92hR1WnAtDzbns318zGgZyBjMsYUzcquMd6z9WiNMcYYP7JEa4wxxviRJVpjjDHGj0J+UQER2Q1sceDU8Xg4GD8M2O/qX/VVtXqAz+k4EckE1jkdRzGEWjkItXgh9GK+UFUrF7VTQDtD+YNTX1AisiiYxz76kv2uxk/WhdLfOtQ+G6EWL4RezCLi0UBwazo2xhhj/MgSrTHGGONHlmi9N9LpAALIflfjD6H2t7Z4/S/UYvYo3pDvDGWMMcYEM6vRGmOMMX5kidZLIvKKiKwVkRUi8qWIVHE6Jl8TkU4isk5E0kXkCafj8RcROV9EvhWRNSLyi4g87HRM4SzUPlci8r6I7BKRVU7H4olQ+zyLSHkR+UlElufE+5zTMXlCRCJEZKmITC1yX2s69o6IXI9rhZKTIvIygKo+7nBYPiMiEcCvQAdcC33/DPRW1dWOBuYHIlILqKWqS0SkMrAYuCkcf1enheLnSkT+AhwCPlLVS52Opyih9nkW16KuFVX1kIiUBb4HHlbVhQ6HVigRGQi0BGJUtUth+1qN1kuq+o2qnsx5uhCoW9j+IagVkK6qG1U1CxgLdHM4Jr9Q1R2quiTn50xgDVDH2ajCVsh9rlR1Lq4ViUJCqH2e1eVQztOyOY+grgGKSF2gM/CuJ/tbovWNe4GvnQ7Cx+oAW3M9zyCIC6uviEgC0AL40dlIwlap/Fw5JVQ+zznNsMuAXcBMVQ3qeIHXgX8A2Z7sbIm2ECIyS0RW5fPolmufp4CTwKfOReoXks+2oL7KLCkRqQR8ATyiqgedjidMlbrPlVNC6fOsqqdUtTmulsFWIhK0TfQi0gXYpaqLPX1PyE/B6E+qel1hr4vIXUAX4FoNv5vdGcD5uZ7XBbY7FIvf5dwb+gL4VFUnOh1PGCtVnyunhOrnWVX3i8h3QCcgWDuftQW6isiNQHkgRkQ+UdU7CnqD1Wi9JCKdgMeBrqp6xOl4/OBnoLGINBCRKKAXMMXhmPwipzPGe8AaVf2P0/GEuVLzuXJKqH2eRaS6e9SGiFQArgPWOhtVwVT1n6paV1UTcH1+/1dYkgVLtCXxFlAZmCkiy0TkHacD8qWcjl79gRm4OlOMV9VfnI3Kb9oCfYBrcv4vl+VcrRofC8XPlYiMARYAF4pIhoj81emYihBqn+dawLcisgLXhdhMVS1yyEwoseE9xhhjjB9ZjdYYY4zxI0u0xhhjjB9ZojXGGGP8yBKtMcYY40eWaI0xxhg/skRrjDHG+JElWmOMMcaPLNGaAonIdyJyYc7PcaGyHqcxpZWINBaRzSLSKOd52Zx1XsNtdbGQYonWFKYRsD7n58uAlQ7GYowpgqquB0YCHXM29Qcmq2qGc1EZW1TA5EtE6gPbVNW9DNRlwAoHQzLGeGYVcJ2IVAP+CrR2OJ5Sz2q0piDNOTuxJmGJ1phQ8CtwITAYeFVVDzsbjrFEawryJ1xLQCEijYFuWNOxMaFgA5AItAI+cjgWgyVaU7DmQBkRWQ48i2ullbucDckYUxRVPQEcBJ7IdevHOMhW7zH5EpF0oIWqZjodizGmeETkN6C+2hd8ULAarTmHiFQGsi3JGhN6RCQB2GJJNnhYjdYYY4zxI6vRGmOMMX5kidYYY4zxI0u0xhhjjB9ZojXGGGP8yBKtMcYY40eWaI0xxhg/CvlFBeLj4zUhIcHpMIzx2uLFi/eoanWn4wg0K7sm1HladkM+0SYkJLBo0SKnwzDGayKyxekYnGBl14Q6T8uuNR0bY4wxfmSJ1hhjjPEjS7TGGGOMH1miNcYYY/zIEq0xxhjjR5ZojTHGGD+yRGuMMcb4kSVaY4wxxo8s0RpjjDF+ZInWGGOM8aOAJVoReV9EdonIqgJeFxF5Q0TSRWSFiCQGKjZjTMGs7BpTMoGs0X4AdCrk9RuAxjmPFGB4AGIyxhTtA6zsGuO1gCVaVZ0L7Ctkl27AR+qyEKgiIrUCE50xpiBWdo0pmWBavacOsDXX84ycbTucCccY4yEruyZkqSonT57kxIkTZGVlceLECU6cOHF628mTJ08/3M9PnTrFqVOnPD5HMCVayWeb5rujSAquJirq1avnz5hMmNq+HerUOfO8b18YORKSkmDJEte2WrVc+w0eDM89d2Zf98puLVue2Zaa6tqvdm3YseMUsIeLLtrJW2/t4aWX9jF79l5gL3CA228/wG+/HWTevANApj9/zUCxsmv85uTJk+zfv5/9+/dz4MABDhw4wMGDB8nMzCQzM5NDhw5x6NAhDh8+zJEjR856HD16lGPHjnH06FGOHz/OsWPHOH78+OlHVlYWWVlZfv8dginRZgDn53peF9ie346qOhIYCdCyZct8C7QpvZKTYerUM89VXUn0/vvPbJsyxbU9r8WLz902eLDr4ZaZmcmvv/7K+PEb2bJlC1u2bGHJki20aLEV1Z2UKbOL7Oxs1q6F6647+1gVKlRg1qwYYmNjadkyhsqVK/PttyX5bYOClV3jMVXljz/+YNu2bWzbto0dO3bw+++/8/vvv7Nr1y727NnDnj172Lt3L/v27SMz07OL0ejoaLKyKnHyZDQQTdmyFWjdugK//16RjRvjgApAObp1K0dUVDk+/7wcEAVEcfXVUXTsWJYhQ6KIjCzLSy+VpWzZskRGRhIZGXnWz5GRkURERBAREcF1eQt4AYIp0U4B+ovIWKA1cEBVrenJFFta2rnbUlJcj+I4duwYK1euZNmyZSxbtozVq1fz66+/sn372TkkJiaG+vXrc/7559OyZUvOO+88zjvvPGrWrEn16tWpVq0acXFxxMXFUa5cuXPOI5JfhTCkWNk1Zzl+/Djr168nPT399OPDDzdz7NgW4DfgyDnviYioRL161fnjFXPu6AAAIABJREFUj+rs338ecAlQjeeeq0bVqlWpUqUKsbGxxMbGEhPjukitVKkSlStXpkKFCpQpU/IuR48/XuJD5CtgiVZExgBXAfEikgGkAmUBVPUdYBpwI5CO63/hnkDFZsJLcnL+ybYov/32G/PmzWPevHnMnz+f1atXn74PU7lyZS655BI6dOhAkyZNuPDCC2nUqBH169enSpUqPv4NgouVXVMQVWXz5s2nL0aXL1/O6tWr2bBhA9nZ2af3i4uLo27dBvzpT5dQv/6NnH/++dSuXZs6depQq1YtatasScWKFR38TfwrYIlWVXsX8boCfw9QOCaM5W42LkxmZiYzZ87kq6++YubMmWzd6urPExMTwxVXXEHXrl1p3rw5LVq0oEGDBj65Yg5FVnaN27Fjx1iwYAHff/89Cxcu5Mcff2Tv3r0AlClThiZNmvCnP/2JXr16cfHFF9OkSRMaNmxIlSpV2L7d1YehNAqmpmNj/G7Xrl2MGzeOyZMnM3fuXE6cOEFsbCwdOnTgscceo127djRr1oyIiAinQzXGcarK8uXLT1+MLly4kOPHj+Pq/3Yx0I3OnVuRmppI166XsnZtBdauhcREeP551+2aUaNyH8+hX8RhlmhN2Dt8+DCTJk3i008/5ZtvvuHUqVM0bdqUAQMG0LlzZ6688koiI60oGANw6tQp5s6dy/jx40lLS2Pbtm0AtGjRgr///e9cffXVVKr0Z6666uxbJjvyuSs/cqTrUdrZt4sJO+6r5vXr1/Pmm28yevRoDh06RP369fnHP/7B7bffziWXXOJskMYEmRUrVvDee+8xfvx4du7cSXR0NJ06daJLly7ccMMNnHfeeaf3za93vimYJVoTUvI2RW3b5ir0Xbu6tygPPzyL9PT/Mm3aNCIjI+nVqxd9+/albdu2pfY+qzH52bjxGA0bfoZrxNWPlCkTxU03dWHp0lvZtKkzEydWZMECuOeec8eTl9ZmYK+oakg/kpKS1IQHV9F1PVJTXdtq1TqzLTGx8PfPnj1bW7durYDWqFFDU1NTdceOHX6Pu6SARRoEZSnQDyu7ztm3b58OGTJEq1SpqYBefPHF+p///Ed3797tdGghxdOyazVaEzQWLXLNzJTb9nynPTjbjz/+yFNPPcXs2bOpW7cuI0eO5M4778x3zKox4Sr3cOwuXVxD3PJO3rJ//wFuu+1Vpk17DThM+fKdmD37Ma6++upwGM8dtCzRmpC1Y8cOBgwYwLhx44iPj+e1116jX79+lC9f3unQjAmokSPzb8p1jyc/fvw4b731Fhdc8C/27dvHrbfeypNPPslll10W2EBLKbthZYJG7rmDC5Odnc0777zDxRdfzKRJk0hNTWXjxo088sgjlmRNqZR7etG8Zs+ezWWXXcajjz7K5ZdfzuLFixk7dqwl2QCyGq0JKatWrSIlJYUFCxZwzTXX8M4779C4cWOnwzIm6OzatYsBAwbw2Wef0bBhQ6ZPn07Hjh2dDqtUshqtCQmqyrBhw0hKSmL9+vV89NFHzJo1y5KsMfmYNm0azZo1Y8KECTz77LOsWrXKkqyDrEZrgkZqav7b9+/fz3333ccXX3zBjTfeyAcffED16tUDG5wxQWzKFNe/x44d4/HHH+eNN96gWbNmzJ49m0svvdTZ4IzVaE3wyL0UndtPP/1EixYtmDx5Mq+++ippaWmWZE3YEjnzcJeH2rXPbHP3yk9JOXvfpCTIyMigbdu2vPHGGzz88MP89NNPlmSDhCVaExDbt5/9xeBesi4p6cy2vBOOf/TRR/z5z39GVZk3bx6DBg2yCSdMWFu06Mxocnei3b79zDb3jEzuXsbux2+/LeTyyy9n/fr1TJkyhddff906BgYRazo2AbF4secLrasqL7zwAqmpqVx77bV8/vnnVK1a1f9BGhOCPvvsM+69917q1KnD7Nmzadq0qdMhmTysemAC4swUiYU7ceIE9957L6mpqdx1111MmzbNkqwpNTwd4uY2bNgwbr/9dtq0acNPP/1kSTZIWaI1QePgwYN07tyZDz74gNTUVEaPHk1UVJTTYRkTlF555RX69+9P165dmT59OnFxcU6HZApgidYEhYMHD9KxY0e+/fZbRo8ezeDBg21KOBOyEhJcszIVp2/C4MFQv75nxx88eDD/+Mc/6NWrFxMmTLD7sUHO7tGagBgxouDXMjMzufHGG1m0aBHjx4+ne/fugQvMGD/YssU1zzB43jdh8OD8e97n9e9//5vnnnuOe+65h1GjRhEREVGSUE0AWKI1AeG+ks/r8OHDdO7cmYULFzJ27FhLssYU4uOPP+bRRx+lZ8+elmRDiCVaExAi517ZHzlyhC5duvDDDz/w2WefcfPNNzsTnDE+VlgLjre+/vpr7r33Xq655ho+/vhjS7IhxO7RGkecPHmSHj16MHfuXD7++GNuvfVWp0MyxmcKasHx1vLly7n55ptp1qwZX375pS0BGWKsRmsCTlV58MEHmT59OiNHjuS2225zOiRjfCq/Fhxv7du3j+7du1O1alWmTZtGTEyMbw5sAsYSrQmILl3O/Pzaa6/xzjvv8Pjjj9O3b1/ngjImyJ06dYrbbruNbdu2MXfuXM477zynQzJesERrAsK9APWkSZN49NFH/3979x5nU73/cfz1IXfjFClmCKfEURGmqy6S0wVD9aCSTnRBnehyuiDn10ydR0mnRxHCiFxSMmXGJo6Q6kTExKmQyEkzQy6NHLfBjM/vjz3DGHPZe89ee+295/N8POZh77XXXvs9tb7z3eu7vhd69uzJyy+/7G4oY8JcYmIiixYtYuLEiVxxxRVuxzEBCuk9WhG5RUQ2icgWERlazOvnicgyEVkrIt+KSJdQ5jPOSUiA9PR0+vTpw2WXXcb06dNt3uIIYmXXP4VbcAK1cOFCXnrpJR588EFr+YlwosG6kVDWB4lUBn4E/gxkAquB3qq6odA+ycBaVR0vIq2ABaratLTjxsfH65o1a5wLboJCZCexse2oUqUKq1at4txzz3U7UtgQkXRV9XPyvdCxsht6e/bs4ZJLLuHss89m9erVNiFFmPK17IbykuJyYIuqblXVo8AsoEeRfRQouNP/B2B7CPMZh+Tl5QH3kp2djcfjsUo28ljZ9VPBZBWBUFUeeeQRfvvtN959912rZKNAKO/RxgEZhZ5nAkVvOiQBn4jIYKAW0Dk00Ux5NG3qnQknK8s7403hBQQmToRdu0YASxg79m1at27tVkwTOCu7fpo/P/D3vvfee3z44YeMGDGCNm3aBC+UcU0or2iLm7i2aLt1b2CqqjYCugAzROS0jCIyQETWiMia3bt3OxDV+GPbNu9QhthY7zf5wutkXnjhZyQmJnLvvffywAMPuB3VBMbKbohkZGTw6KOP0qFDB5555hm345ggCWVFmwk0LvS8Eac3Lz0IzAZQ1a+A6sDZRQ+kqsmqGq+q8fXr13corvFVSf00du7cSe/evWnevDnjx4+3RQIil5XdEBk8eDDHjh1j2rRpNvNTFAllRbsaaC4izUSkKnA34Cmyzy/AjQAi8ie8hdW+9oa55OTTt+Xl5XHvvfeyb98+UlJSqF27duiDmWCxsuunQPqYfvzxx8ydO5fnn3+e888/P/ihjGtCVtGqai4wCFgEbARmq+p6EXlRRAru6j0F9BeR/wDvA/00VN2iTcDatz992+jRo1myZAljxozhkksuCX0oEzRWdv1X3JfP0hw+fJjBgwfzpz/9iSeffNKZUMY1IRve4xQbIuC+otPN/fjjj7Rp04abb76Z1NRUazIuQ7gP73FKNJddf6dgfP755/nHP/7Bp59+yg033OBcMBNU4Ti8x1QAeXl53H///dSoUcPuyxrjg82bNzNy5Ejuueceq2SjlE3BaMqtYcOTj998801WrFjBjBkzaFj4BWNMsZ5++mmqVavGa6+95nYU4xC7ojXltj2//+nmzZt57rnnSEhIoE+fPu6GMsZFnqJdxUqwfPlyPB4PQ4YMsS+mUcwqWlNuSUknm4yrV6/OhAkTrMnYVGjFdRAsSlUZMmQIDRo04IknnnA+lHGNVbSm3F54Ad566y2WL1/O6NGjiY2NdTuSMUGRnu7t2FTwk5Tk3R4be3JbQaU6YMDJbXFxZR97/vz5LF++nMTERGrVquXY72DcZ72OTbmJ/EqdOi248sor+de//mVXs36yXsfhKz3dt6tTf+Xl5dGmTRuOHj3K+vXrqVKlSvA/xDjO17JrnaFMEAzl8OHDjBkzxipZE1Xi4wObfKIs7777LuvXr+eDDz6wSrYCsKZjUy4rVqwApvHUU09x4YUXuh3HmLCXm5vLiy++SLt27ejZs6fbcUwI2BWtCVheXh6PPvoo557biL///e9uxzEmIsyePZutW7eSmppKpUp2rVMRWEVrAjZx4kTWrVsHfGCdOUxUSkwM7vGOHz/OiBEjaNWqFd0Lrydpopp9nTIB2b17N8OHD6dTp05AL7fjGOOIgl7GwTJ//ny+//57hg0bZlezFYj9nzYBGT58OAcOHGDMmDEUv1ypMZEvmCPVVJWXXnqJpk2bcvfddwfvwCbsWUVrytS0qXeMIHiHOoisZ9KkyVSrNohWrVoFvXnNmHCxY0fwjrVs2TK+/vprhgwZwhln2F27isQqWlOmbdtOLvuVng4JCcOoU6c227Z5O0AFu3nNmGg0YsQIGjRoQL9+/dyOYkLMKlrjl3//+9/MmzePoUOHUq9ePbfjGOOodu2Cc5zvv/+eJUuW8Pjjj1O9evXgHNREDKtoTZmysrz/FszNGhsby+OPP+5uKGNCID09OMcZO3Ys1atXp3///sE5oIkoVtGaMhX8sZk7dy5fffUVL7zwAjVr1nQ3lDEhUNA3oTz27t3LjBkz6NOnj7UCVVB+V7QiUktEKjsRxoSn7t29s9kMGzaMli1b2j2mCGVl13+TJpX/GFOmTOHQoUMMHjy4/AczEanMrm8iUgm4G+gDXAYcAaqJyG5gAZCsqpsdTWlc98477/DDDz+QmppqPSYjhJVd9+Xl5TF27Fiuu+462rRp43Yc4xJfrmiXAecDw4AGqtpYVc8BrgVWAq+IyL0OZjSuO0xiYiJXX301PXr0cDuM8Z2VXZd9/PHH/Pzzzzz22GNuRzEu8uXSpLOqHiu6UVWzgY+Aj0TElp+IYr16TSQlZQfvv/++rc4TWazsllNBR8BAvfnmmzRu3Ni+oFZwvlzRjhaRR0Skg4jUKW6H4gqziQ6HDx/m3/8eSceOHbn++uvdjmP8Y2W3nMrT63jz5s0sXbqUhx9+2G63VHC+/N9fB7TGe6/nYhHZD3wHfAt8p6qzHMxnXDZp0iR+/fVX3n//fbejGP9Z2S2n7t0DX492ypQpVK5cmfvvvz+4oUzEKbOiVdXkws9FpBHewnsJ0BWwwhqlcnJyGDlyJHAdHTt2dDuO8ZOVXffk5uYydepUunbtSsOGDd2OY1zm9/AeVc1U1QWqOlJV/+LPe0XkFhHZJCJbRGRoCfvcKSIbRGS9iLznbz4TPG+//Tbbt28HktyOYoLAym7oLFiwgF9//ZUHH3zQ7SgmDITsxkH++L1xwJ+BTGC1iHhUdUOhfZrj7SHZQVX3isg5ocpnTpWTk8OIESO49tprqVOno9txjIsqctmdODGw902ePJkGDRrQpUuX4AYyESngmaFEpKGIVPPjLZcDW1R1q6oexdtsVbQrXn9gnKruBVDVXYHmM+UzefJktm/fTmJiIvPnW0/jaGJl13eBzAy1Y8cOPv74Y/r27WudoAxQvikYZwA/iMhrPu4fB2QUep6Zv62wC4ELRWS5iKwUkVvKkc8E6MiRI7zyyit06NCBTp06kZDgdiITZFZ2fRTIaLZp06aRl5dnzcbmhIC/bqlqZ/EOqmzl41uKO2WL9uc7A2gOdAQaAf8WkYtV9fdTDiQyABgAcN555/kT2/hg6tSpZGZmMmXKFESE+fPdTmSCycquc1SVyZMnc91119G8eXO345gw4fMVrYiMkiKzFajXeh8PkQk0LvS8EbC9mH3mquoxVf0vsAlv4T2FqiararyqxtevX9/XX8H4IC8vj9dee434+Hg6d+7sdhwTBFZ2Q2f58uVs2bKFBx54wO0oJoz403R8APCISC0AEblJRJb78f7VQHMRaSYiVfGO7fMU2ScNuCH/+GfjbY7a6sdnmHJKTU1ly5YtDBkyxGaBih5WdgPUrZt/+8+cOZMaNWpwxx13OBPIRCSfm45V9e8icg/wmYgcAQ4CxXbzL+H9uSIyCFgEVAamqOp6EXkRWKOqnvzXbhKRDUAe8Iyq/ubH72PKQVUZOXIkF1xwAbfffnuh7S6GMuVmZTdw8+b5vu/Ro0eZPXs2PXr0ICYmxrlQJuL4XNGKyI14exYeBBoCD6rqJn8+TFUX4F01pPC25ws9VuBv+T8mxJYtW8aaNWuYMGEClSufXE0tOTk463Iad1jZDVxCgu+V7aJFi8jOzqZPnz7OhjIRx5+m4+HA86raEegJfCAinRxJZVzx6quvcu6559K3b99Ttg8c6FIgEyxWdgPkT0fAmTNnUq9ePW6++WbnApmI5E/TcadCj78TkVvxrgBytRPBTGitW7eORYsW8fLLL1O9enW345ggsrLrvP379+PxeOjXrx9VqtiCSOZUviz8Xlof/AcLvf67qv4vOLFMqL366qvUrl2bRx55xO0oJkis7IZOamoqhw8ftmZjUyxfrmin4R0zV1oXVAWmAtODkMmEUNOmsG3bf4HZPPXUEzz77JlMmnTy9aws8BTtX2oiRWllt2C7ld1S+NoRcObMmTRt2pSrr7ZGAnM6X1bvuSEUQYw7jh6FwYPfYMKESjz55JPExXk7PxUWG+tONlM+BWVXRGJVtei4V+MDXzoC7ty5kyVLljB06FAbEmeK5c+EFd+JyEwRGSIit4pIIxEZ7mQ447wNG35nypQp9O7dm7i4orPqmSixQEQSRcRuvvvJl46AH374IcePH+eee+5xPpCJSP70Or4emAQcxjtg/Xu8a1qaCHbnnW9z8OBBnnzySbejGOfEA//Du+rOfW6HiTYpKSm0atWKiy66yO0oJkz5XNGqaraqfqaqb6pqX+AyYLNz0YzTcnNzWbx4DB07duTSSy91O45xiKrmquobwHVAexFZISLXup0rGvz666988cUX9OrVy+0oJoz5M2FFc1U9UbGq6mYRae1MLBMKqampwC888cSbbkcxDhKRPwI3Ay3yfy4A3hGRKsDPqnq9m/nCWVkdAefMmYOqWkVrSuXP6j3JInI+kAV8C1QHvheRmqp6yJF0xlFvvPEGcD7d/J3Q1USapcDE/H/HAz+pai6AiDRxM1i4a9++9NdTUlJo2bIlrVr5uhCSqYj8mbCioAfjecClQJv8f/8jInmq2tKZiMYJq1at4quvvuLpp988ZbpFE5X+rKpbintBVbeFOkwkiYsreYjPzp07+eKLLxg+fLj1Njal8ns9WlX9BfiFQqt3iEjtYIYyzhs1ahR16tQhIaGf21GMw0qqZE35zJkzh+PHj1uzsSmTP72OS6SqB4JxHBMaGRkZpKSk0L9/f66/3lYZMSYQKSkptGjRgosvvtjtKCbMBaWiNZFl7NixqCqDBw92O4oxYa1//+K379q1i88//5xevXpZs7Epk08VrYjUFJE2RbadJyI2w0GEOXToEJMmTeKOO+6gSRPrBxPtrOyWT9FZ0gqkpqZy/PhxevbsGdpAJiL5ekV7DJgjIrUKbXsb79qWJoLMnDmTvXv38thjjwGQmOhyIOM0K7vlUFKv47S0NM4//3xat7YRjqZsPlW0qnoMSAXughM9j+ur6hoHs5kgU1XGjBlDmzZtuOaaawBISnI3k3GWld3y+eab07ft27ePpUuXctttt1mzsfGJP/do3wbuz398H/BO8OMYJ33xxRd89913DB48+MQfCFswoEKwshtECxcu5NixY9x+++1uRzERwp9xtD+ICCJyIdAbuMa5WMYJY8aMoW7duqdMfr5jh4uBTEhY2Q1cw2Ia2FNTUzn33HO58sorQx/IRCR/ex1Pxvvt+FtV3etAHuOQjIwM0tLSeOihh6hRo4bbcUzoWdkNwPYiiwseOXKEBQsW0L17d5voxfjM34p2Nt4ZoSY7kMU4aPz48agqf/3rX0/Z3q6dS4FMqFnZDUDRPgxLly7lwIED1mxs/OJXRauqh1T1D6q6xKlAJvhycnKYNGkS3bt3P21IT3q6S6FMSFnZDcwLL5z6PC0tjZiYGDp16uROIBORbMKKCmDWrFns2bOn2AkqBgxwIZAxESgvL4+5c+fSpUsXqlWr5nYcE0HKrGjFh/7rvuxj3DFxonL//WOAVtx44w3Mm+e97yTi/fnkE7cTGqdY2Q2ulStXsmvXLm677Ta3o5gI48sV7TIRGZw//u4EEakqIp1EZBrQ15cPE5FbRGSTiGwRkaGl7NdTRFRE4n05rimeCLRuvRL4hrfeGoSqkJDgHdKj6v35+We3UxoHWdktpzWFRhunpaVRpUoVunTp4l4gE5F8qWhvAfKA90Vkh4hsEJGtwGa8QwXeUNWpZR1ERCoD44BbgVZAbxE5bRFHEYkBHgNW+fxbmBKNGzeOOnXq8Je//MXtKCb0rOwW0rSp98tnwe2S9u1PtuwUjCdPSjq5rei1vsfjoVOnTtSpUyeUsU0UKHMcrarmAG8Bb4lIFeBs4LCq/u7nZ10ObFHVrQAiMgvoAWwost8/gFeBp/08vjnNTmbPns0jjzxC7dq2kmFFY2X3VNu2nbq2bHEdAZOSip8tbdOmTfz44488/vjjTsUzUcyXe7R9RWSPiGTjHYd3IICCChAHZBR6npm/rfBntQUaq+r8AI5vimjR4m2OHTt22pAeUzFY2T1VVlbg7/V4vMtvJyQkBCmNqUh8aTr+P+DPQEu8C76/HOBnFdfp4sT3SxGpBLwBPFXmgUQGiMgaEVmze/fuAONEt9zcXA4enEDnzp1p0aKF23GMO6zsFlKeoWwej4e2bdvSuHHj4AUyFYYvFe3/VHWtqu5S1f/D24wUiEyg8FnaCCg870oMcDHwmYj8DFwJeIrrVKGqyaoar6rx9evXDzBOdPN4PGRmZjJo0CC3oxj3WNktpHv3wN63e/duVqxYQfdAD2AqPF/mOm4oIgOAjcAPQJUAP2s10FxEmgFZwN3AiUl3VXUf3ntIAIjIZ8DTtspIYMaNGwecR7du3dyOYtxjZTcIFixYwPHjx62iNQHzpaJNBFoDfYBLgNoisgD4D955U9/35YNUNVdEBgGLgMrAFFVdLyIvAmtU1RPQb2BOs3HjRj799FPgZZuPtWKzshsEHo+HuLg42rZt63YUE6F86XWcXPi5iDTCW3gvAboAPhXW/GMtABYU2fZ8Cft29PW45lTjxo2jatWqHD36kNtRjIus7J5q4kT/35OTk8OiRYu47777bO1ZEzCfl8kroKqZeO/ZLChrXxN6+/fvZ/r06dx1111Mn273r81JFb3sBjLd6LJlyzh48KA1G5tysbmOo8yMGTPYv38/jz76KMnJZe9vTEURyAWpx+Ohdu3a3HDDDcEPZCoMq2ijiKoyduxY2rdvz+WXX87AgW4nMiZyqSrz5s3jpptuskUETLn43XRswteyZcvYuHEjU6dOtftJxpTT2rVrycrKsmZjU252RRtFxo4dS7169bjrrrvcjmJM2PF3pJvH46FSpUq2iIApN6too8Qvv/zC3Llz6d+/P9WrVwfAE9WDLozxz7x5/u3v8Xi46qqrsElxTHlZRRslJkyYAMDDDz98Ylv79m6lMSb8+DNNcWZmJmvXrrVmYxMUVtFGgZycHCZNmkT37t1p0qTJie1xcaW8yZgKZr4fyx3My7/8tUUETDBYRRsFZs+ezZ49e2xeY2OCZN68eVxwwQW0bNnS7SgmClhFGwXGjh1Ly5Yt6dSpk9tRjIl4Bw4cYOnSpSQkJFjvfRMUVtFGuK+//prVq1czaNCg0/4o9O/vUihjwlDhRd9Ls3jxYo4ePWr3Z03QWEUb4caMGUNMTAz33Xffaa/ZzFDGnORrefB4PJx55pl06NDB2UCmwrCKNoLt2LGDDz74gH79+hETE3Pa69br2JiTfJkpLS8vj/nz59OlSxeqVAl0VUFjTmUVbQSbMGECubm5DB48uNjXv/kmxIGMiXArVqxgz5493HbbbW5HMVHEKtoIJZLDmDHjueaablx4YXNEvJOmJyV5X4+NhUIjfYwxPkhLS6Nq1arccsstbkcxUcTmOo5Y77N3724SEx/nxhtPf3X79tAnMiaclTVTmqoyd+5cbrzxxmJvxRgTKLuijUCqCozm4osvtiE9xviorD4L69ev56effrJmYxN0dkUbgT7//HPgPzzxxNs2zs8YH8XFlT7EJy0tDRGxYT0m6OyKNgKNGjWKevXqcc8997gdxZioMXfuXK644goaNGjgdhQTZayijTBbt27F4/Fw5MjD1KhRw+04xkSFjIwM1qxZY83GxhFW0UaYMWPGULlyZQ4c+KvbUYyJKKXNlObJ7yllFa1xglW0EWTfvn1MnjyZO++8E4h1O44xEaW0maHS0tJo2bIlLVq0CF0gU2FYRRtBxo8fz/79+3n66adp187tNMZElpJ6HWdnZ/PZZ5/Ro0eP0AYyFYZVtBEiJyeHUaNGcdNNN9G2bVvS091OZExkKWmmtLlz55Kbm0uvXr1CG8hUGCGtaEXkFhHZJCJbRGRoMa//TUQ2iMi3IrJURGxuo3zTpk1j586dDB3q/c82YIDLgUyFEs1ld/bs2TRr1ox21kxkHBKyilZEKgPjgFuBVkBvEWlVZLe1QLyqtgY+BF4NVb5wlpeXxz//+U8uu+wyOnbsCMCkSe5mMhVHtJTdhg1P35adnc2SJUvo1auXjUk3jgnlFe3lwBZV3aqqR4H7h2W6AAAPAElEQVRZwCk3RVR1maoeyn+6EmgUwnxh66OPPuKnn35i6NCh9sfAuCEqym5x05Jas7EJhVBWtHFARqHnmfnbSvIgsNDRRBFAVXnllVdo0aKFDT0wbomKsluw4EZhKSkpNG3alPa2pqRxUCgr2uIuxYqdEE1E7gXigX+W8PoAEVkjImt2794dxIjhZ/Hixaxdu5Znn32WSpVO/u/KynIxlKlooqLsvvDCqc+zs7NZvHgxd955p7UUGUeFsqLNBBoXet4IOK0xR0Q6A8OB7qp6pLgDqWqyqsaranz9+vUdCRsuXnnlFWJjY+nTp88p263XsQmhsCq7BR0B27fnxPKQsfnDypOSTm4T8ZaT9HTv46LLRlqzsQmVUFa0q4HmItJMRKoCdwOnLFwlIm2BiXgL6q4QZgtLK1euZNmyZfztb3+jWrVqp7xm856bEAqrslvQETA93btIgOrJ+69JSSe3qXor4/btvY9//vnU41izsQmVkFW0qpoLDAIWARuB2aq6XkReFJGCauOfQG0gRUTWiUgZK0hGL1Xlueeeo379+gwcONDtOKYCi8ayW9BsbL2NTSiEdJk8VV0ALCiy7flCjzuHMk84W7JkCcuWLWP06NHUrl3b7Timgou2sjtr1ixyc3Pp3bu321FMBWAzQ4UhVWXYsGE0adKkxKvZiRNDHMqYMBGMjoBTp06ldevWtG3btvwHM6YMVtGGoY8++oj09HRefPHF0+7NFrCZoUxFVd6OgBs2bGD16tX069cvKHmMKYtVtGEmNzeX4cOHc9FFF53W07gwu61kKqrydgScNm0aZ5xxRqnly5hgCuk9WlO2qVOn8uOPP5KWlkblypXdjmNMVMnNzWXGjBnceuutnHPOOW7HMRWEXdGGkcOHD/PCCy9w5ZVX0t3G7xgTdEuWLGHHjh3WbGxCyq5ow8ibb75JZmYmM2bMKHPIQbduIQplTJgpT0fAqVOnUrduXbp27Rq8QMaUwa5ow8TmzZtJSkritttuO7FCT2nmzXM+kzHhKNCOgL///jtpaWncc889JXYyNMYJVtGGgePHj9O/f3+qVavGuHHjfHpPQoLDoYwJU4F2BHznnXc4cuQI999/f3ADGVMGazoOA8nJyXz++edMnjyZ2IJJW8swf77DoYyJIrm5uYwePZprr73WFng3IWcVrcsyMjJ49tln6dy5s33TNsYhqampbNu2jVGjRrkdxVRA1nTsIlVl4MCB5OXlkZycbHOuGuODQDoCvv7665x//vkk2D0X4wK7onXRzJkzWbhwIaNHj6ZZs2Z+vVeLXQ3UmOjnb0fAr776ipUrVzJ27Fgbm25cYVe0Llm3bh2PPPIIV199NY8++qjf709OdiCUMRHA34vS119/nbPOOsvGzhrXWEXrgoyMDLp27cqZZ57J7NmzA/qWbSvnmYrKn46A//3vf5kzZw4DBw6kVq1azoUyphTWdBxi+/bto2vXrhw4cIAvv/ySuLg4tyMZE7WSkpI444wzGDRokNtRTAVmFW0IHTt2jJ49e7Jx40YWLlzIJZdc4nYkY6LWqlWrmD59OkOHDrUvtMZVVtGGSF5eHgMGDGDJkiW88847dO5cvnWyPZ4gBTMmwvjSEVBVeeKJJ2jQoAHPPfec86GMKYXdow2BvXv30q1bN6ZOnUpiYmJQOmW0b1/+XMZEIl86Ar733nusXLmSESNGEBMT43woY0phFa3Dvv32W+Lj41m6dCkTJkwgMTExKMe1ljBTUZXVEfDgwYMMGTKE+Ph47rvvvtCEMqYUVtE6aNasWVx11VXk5OTw+eefM3DgQJuUwph827d7/42N9c5fLHKypWbAgJPbRLz7zpvnfdykSenHHT58OFlZWYwaNYpKlexPnHGfnYUOWL9+PXfccQe9e/emXbt2pKenc9VVV/n8/qZNISnJ+7ikP0Jl/bExJtzt2OH9d/t2731XVUhP925LTj65TdVbDhISvI9//rnkY06bNo3Ro0czePBgOnTo4PjvYIwvRCN8iqH4+Hhds2aN2zEA75i9xMRE3n33XWrXrs0zzzzDkCFDqFq1ql/HEbGZnyoSEUlX1Xi3c4SaSLyqBq/srlixghtuuIFrr72Wf/3rX5xxhvX1NM7ytezamVhOR44c4ZNPPmHWrFmkpKRQuXJlnn76aYYMGUK9evXcjmdMhZCRkcEdd9xB48aN+eCDD6ySNWHFzsYA7N27lxUrVjBnzhzmzJnD77//Tt26dRk4cCDDhg3zeam7koTJBboxjvrTn4JznFWrVnH33Xdz6NAhli5dal9wTdgJaUUrIrcAo4HKwNuq+kqR16sB04H2wG/AXar6cygzFrVnzx42bdrExo0bWblyJStWrGDjxo0AxMTEcPvtt3PXXXfRuXNnv5uIjYkU4Vh2jx8/zhtvvHFiQorFixdz0UUXOfmRxgQkZBWtiFQGxgF/BjKB1SLiUdUNhXZ7ENirqheIyN3ASOCuYGdRVQ4fPsxvv/1GdnY22dnZ7Ny5k6ysLLZv305WVha//PILmzZtIjs7+8T76taty1VXXcW9997L1VdfzZVXXkn16tWDHY/4eLtHa8KHU2U3//uq3w4ePMicOXNITk7myy+/5Pbbb2fy5MmcddZZgR3QGIeF8or2cmCLqm4FEJFZQA+gcGHtASTlP/4QGCsioqX02Nq+fTvDhw/nyJEjHD16lCNHjpCTk8Phw4dP/Bw6dIj9+/dz4MAB9u/fz/79+zl69Gixx6tZsyZxcXE0atSIXr160aJFCy688EJatGjBH//4RxsuYCoiR8puYXv37iUrK+vE8+PHj3Ps2DGOHDnC4cOHycjI4KeffmLTpk0sXLiQAwcO0KxZM8aPH2/D5kzYC2VFGwdkFHqeCVxR0j6qmisi+4B6wJ6SDrpjxw5GjhxJtWrVqFq1KtWqVaNatWrUrFmTGjVqULNmTWrVqkWDBg2IiYmhdu3axMTEcNZZZ1GvXj3q1q1L3bp1Oeecc4iLi6NOnTpWaI05lSNlt7A5c+bw0EMPlbpPpUqVOO+887jzzjvp27cv11xzjX3xNREhlBVtcbVX0W+7vuyDiAwABgD84Q/nsW/fNg4dgkOHTnYkii/U4Tox0TsuNTb25Ni9du28Y/YGDIBJk07um5Xl3d69+8ltEyeeHEBfoFs37wD6hIRTl+1S9Y4BLDx7jcfjHQNbeDan/v29+7VvD998493WsKF3TGGQJo8yJlgcK7sFOnXqREpKyin7Fv7y3KhRI5o0aUKVKlX8zW6M60I2jlZErgKSVPXm/OfDAFR1RKF9FuXv85WInAH8CtQvrfkpnMbRGhOIcB9Ha2XXmOL5WnZD2e6yGmguIs1EpCpwN1B0DRoP0Df/cU/gU1/v8RhjHGNl15hyCFnTcf59m0HAIrxDBKao6noReRFYo6oeYDIwQ0S2ANl4C7QxxkVWdo0pn5COo1XVBcCCItueL/Q4B+gVykzGmLJZ2TUmcNZlzxhjjHGQVbTGGGOMg6yiNcYYYxwU8cvkichuYJsLH302Pg7GjwL2uzqriarWD/Fnuk5E9gOb3M7hh0grB5GWFyIvcwtVjSlrp4hfvcetP1Aisiacxz4Gk/2uxiGbIum/daSdG5GWFyIvs4j4NBDcmo6NMcYYB1lFa4wxxjjIKtrAJbsdIITsdzVOiLT/1pbXeZGW2ae8Ed8ZyhhjjAlndkVrjDHGOMgq2gCJyD9F5AcR+VZEUkXkTLczBZuI3CIim0Rki4gMdTuPU0SksYgsE5GNIrJeRB53O1M0i7TzSkSmiMguEfne7Sy+iLTzWUSqi8jXIvKf/LwvuJ3JFyJSWUTWisj8Mve1puPAiMhNeFcoyRWRkQCqOsTlWEEjIpWBH4E/413oezXQW1U3uBrMASLSEGioqt+ISAyQDtwWjb+r2yLxvBKR64ADwHRVvdjtPGWJtPNZRASopaoHRKQK8CXwuKqudDlaqUTkb0A8UEdVu5W2r13RBkhVP1HV3PynK4FGbuZxwOXAFlXdqqpHgVlAD5czOUJVd6jqN/mP9wMbgTh3U0WtiDuvVPULvCsSRYRIO5/V60D+0yr5P2F9BSgijYCuwNu+7G8VbXA8ACx0O0SQxQEZhZ5nEsaFNVhEpCnQFljlbpKoVSHPK7dEyvmc3wy7DtgFLFbVsM4LjAKeBY77srNVtKUQkSUi8n0xPz0K7TMcyAVmupfUEVLMtrD+llleIlIb+Ah4QlX/53aeKFXhziu3RNL5rKp5qnop3pbBy0UkbJvoRaQbsEtV0319T8RPwegkVe1c2usi0hfoBtyo0XezOxNoXOh5I2C7S1kcl39v6CNgpqrOcTtPFKtQ55VbIvV8VtXfReQz4BYgXDufdQC6i0gXoDpQR0TeVdV7S3qDXdEGSERuAYYA3VX1kNt5HLAaaC4izUSkKnA34HE5kyPyO2NMBjaq6utu54lyFea8ckuknc8iUr9g1IaI1AA6Az+4m6pkqjpMVRupalO85++npVWyYBVteYwFYoDFIrJORCa4HSiY8jt6DQIW4e1MMVtV17ubyjEdgL8AnfL/X67L/7ZqgiwSzysReR/4CmghIpki8qDbmcoQaedzQ2CZiHyL94vYYlUtc8hMJLHhPcYYY4yD7IrWGGOMcZBVtMYYY4yDrKI1xhhjHGQVrTHGGOMgq2iNMcYYB1lFa4wxxjjIKlpjjDHGQVbRmhKJyGci0iL/cb1IWY/TmIpKRJqLyM8ickH+8yr567xG2+piEcUqWlOaC4DN+Y9bA9+5mMUYUwZV3QwkAzfnbxoEzFXVTPdSGVtUwBRLRJoAWapasAxUa+BbFyMZY3zzPdBZROoCDwJXuJynwrMrWlOSSzm1Ym2PVbTGRIIfgRZAEvCaqh50N46xitaUpA3eJaAQkeZAD6zp2JhI8BPQDrgcmO5yFoNVtKZklwKVROQ/wPN4V1rp624kY0xZVPUY8D9gaKFbP8ZFtnqPKZaIbAHaqup+t7MYY/wjIr8ATdT+wIcFu6I1pxGRGOC4VbLGRB4RaQpss0o2fNgVrTHGGOMgu6I1xhhjHGQVrTHGGOMgq2iNMcYYB1lFa4wxxjjIKlpjjDHGQVbRGmOMMQ6yitYYY4xxkFW0xhhjjIP+H1PqVLZd7TzfAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light", "tags": [] }, "output_type": "display_data" } ], "source": [ "fig = plt.figure(figsize=(7, 7))\n", "fig.subplots_adjust(wspace=0.35, right=0.95,\n", " hspace=0.2, top=0.95)\n", "\n", "# first axes: mu posterior\n", "ax1 = fig.add_subplot(221)\n", "ax1.plot(mu, pmu, '-k')\n", "ax1.hist(mu_bootstrap, mu_bins, density=True,\n", " histtype='step', color='b', linestyle='dashed')\n", "ax1.set_xlabel(r'$\\mu$')\n", "ax1.set_ylabel(r'$p(\\mu|x,I)$')\n", "\n", "# second axes: mu cumulative posterior\n", "ax2 = fig.add_subplot(223, sharex=ax1)\n", "ax2.plot(mu, pmu.cumsum() * dmu, '-k')\n", "ax2.hist(mu_bootstrap, mu_bins, density=True, cumulative=True,\n", " histtype='step', color='b', linestyle='dashed')\n", "ax2.set_xlabel(r'$\\mu$')\n", "ax2.set_ylabel(r'$P(<\\mu|x,I)$')\n", "ax2.set_xlim(-3, 3)\n", "\n", "# third axes: gamma posterior\n", "ax3 = fig.add_subplot(222, sharey=ax1)\n", "ax3.plot(gamma, pgamma, '-k')\n", "ax3.hist(gamma_bootstrap, gamma_bins, density=True,\n", " histtype='step', color='b', linestyle='dashed')\n", "ax3.set_xlabel(r'$\\gamma$')\n", "ax3.set_ylabel(r'$p(\\gamma|x,I)$')\n", "ax3.set_ylim(-0.05, 1.1)\n", "\n", "# fourth axes: gamma cumulative posterior\n", "ax4 = fig.add_subplot(224, sharex=ax3, sharey=ax2)\n", "ax4.plot(gamma, pgamma.cumsum() * dgamma, '-k')\n", "ax4.hist(gamma_bootstrap, gamma_bins, density=True, cumulative=True,\n", " histtype='step', color='b', linestyle='dashed')\n", "ax4.set_xlabel(r'$\\gamma$')\n", "ax4.set_ylabel(r'$P(<\\gamma|x,I)$')\n", "ax4.set_ylim(-0.05, 1.1)\n", "ax4.set_xlim(0, 4)\n", "\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": {}, "colab_type": "code", "id": "fWo5_d-UVr_5" }, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": {}, "colab_type": "code", "id": "Zn_DSiZ3VsAE" }, "outputs": [], "source": [] } ], "metadata": { "colab": { "name": "Figure5-10,11.ipynb", "provenance": [] }, "jupytext": { "formats": "ipynb,md:myst" }, "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.1" }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": false, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": false, "toc_position": {}, "toc_section_display": true, "toc_window_display": false } }, "nbformat": 4, "nbformat_minor": 1 }