josef-pkt/Mixture_MLE.ipynb

## Mixture_MLE.ipynb
{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fit the Heights of Two Peaks\n",
    "\n",
    "Assume we have two components A and B -- a signal peak over a flat background -- and measure their sum. We know the shape and location of both, but we are interested in the relative contributions from either component:\n",
    "if we see 1000 events, how many of them are generated by process A and B, respectively?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from scipy import stats\n",
    "from statsmodels.base.model import GenericLikelihoodModel"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate some pseudodata."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(42)\n",
    "pdf_a = stats.halfcauchy(loc=0, scale=1)\n",
    "pdf_b = stats.uniform(loc=0, scale=100)\n",
    "\n",
    "n_a = 30\n",
    "n_b = 1000\n",
    "\n",
    "X = np.concatenate([\n",
    "    pdf_a.rvs(size=n_a),\n",
    "    pdf_b.rvs(size=n_b),\n",
    "])[:, np.newaxis]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABFMAAAM9CAYAAABQfTw1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAewgAAHsIBbtB1PgAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3X+slvV9//EXP84B9NjWKgXPgNCSbnoWjY0BEZJGq6uTw+kW2UbhzM1SFzGtP+KUuiyd07WbO2mcZnXoOuncBng2zrHqgdbV6pYGcTRLaiwHF0UtoEeEKQoKPSjn+4df7nHK4Xg+CpxzHx6PhOQ69+dzX/f7Pn8YzzPXdd8jenp6egIAAADAgIwc7AEAAAAAqomYAgAAAFBATAEAAAAoIKYAAAAAFBBTAAAAAAqIKQAAAAAFxBQAAACAAmIKAAAAQAExBQAAAKCAmAIAAABQQEwBAAAAKCCmAAAAABQQUwAAAAAKiCkAAAAABcQUAAAAgAJiCgAAAEABMQUAAACggJgCAAAAUEBMAQAAACggpgAAAAAUGD3YAxwv9u7dm6effjpJMn78+Iwe7VcPAAAAR9o777yT7du3J0nOPPPMjB079oi/hr/oj5Gnn346M2bMGOwxAAAA4Lixfv36TJ8+/Yif120+AAAAAAVcmXKMjB8/vnK8fv36nHbaaYM4DQAAAAxPXV1dlTtDDv5b/EgSU46Rgz8j5bTTTsukSZMGcRoAAAAY/o7W55W6zQcAAACggJgCAAAAUEBMAQAAACggpgAAAAAUEFMAAAAACogpAAAAAAXEFAAAAIACYgoAAABAATEFAAAAoICYAgAAAFBATAEAAAAoIKYAAAAAFBBTAAAAAAqIKQAAAAAFxBQAAACAAmIKAAAAQAExBQAAAKCAmAIAAABQQEwBAAAAKCCmAAAAABQQUwAAAAAKiCkAAAAABcQUAAAAgAJiCgAAAEABMQUAAACggJgCAAAAUEBMAQAAACggpgAAAAAUEFMAAAAACogpAAAAAAXEFAAAAIACYgoAAABAATEFAAAAoMDowR6AoW/qTasHe4Sq8uJtjYM9AgAAAEeRK1MAAAAACogpAAAAAAXEFAAAAIACYgoAAABAATEFAAAAoICYAgAAAFBATAEAAAAoIKYAAAAAFBBTAAAAAAqIKQAAAAAFxBQAAACAAmIKAAAAQAExBQAAAKCAmAIAAABQQEwBAAAAKCCmAAAAABQQUwAAAAAKiCkAAAAABcQUAAAAgAJiCgAAAEABMQUAAACggJgCAAAAUEBMAQAAACggpgAAAAAUEFMAAAAACogpAAAAAAXEFAAAAIACYgoAAABAATEFAAAAoICYAgAAAFBATAEAAAAoIKYAAAAAFBBTAAAAAAqIKQAAAAAFxBQAAACAAmIKAAAAQAExBQAAAKCAmAIAAABQQEwBAAAAKCCmAAAAABQQUwAAAAAKiCkAAAAABcQUAAAAgAJiCgAAAEABMQUAAACggJgCAAAAUEBMAQAAACggpgAAAAAUEFMAAAAACogpAAAAAAXEFAAAAIACYgoAAABAATEFAAAAoICYAgAAAFBATAEAAAAoIKYAAAAAFBBTAAAAAAqIKQAAAAAFxBQAAACAAmIKAAAAQAExBQAAAKCAmAIAAABQQEwBAAAAKCCmAAAAABQQUwAAAAAKiCkAAAAABcQUAAAAgAJiCgAAAEABMQUAAACggJgCAAAAUEBMAQAAACggpgAAAAAUEFMAAAAACogpAAAAAAXEFAAAAIACYgoAAABAgaMaU15++eU89NBD+frXv55LLrkk48ePz4gRIzJixIhcfvnlxedrb29PU1NT6uvrM2bMmNTX16epqSkPPPDAgM/x5ptv5hvf+EbOOeecnHzyyTnxxBPzq7/6q7n66qvz7LPPFs8EAAAAHF9GH82T/8qv/MoROU93d3cWLFiQ9vb2Xo93dXWlo6MjHR0dmTdvXlauXJmamprDnuenP/1pmpqasnXr1l6PP/vss3n22Wdz77335jvf+U6am5uPyNwAAADA8HNUr0w57bTT8oUvfCG33nprvv/97+cnP/nJBzrPVVddVQkpCxcuzIYNG7Jnz55s2LAhCxcuTJK0tbVl8eLFhz3Htm3bMmfOnGzdujUnnnhi7rrrrrz66qt58803s3r16kybNi179uzJ5Zdfnscff/wDzQkAAAAMf0f1ypSXX365188vvvhi8TnWrVuXZcuWJUnmz5+f5cuXV9YaGhqyfPnyvPvuu2ltbc2yZcvyR3/0R5k5c+Yh57n55pvT1dWVJPm3f/u3XHLJJZW1OXPm5Oyzz86ZZ56Z1157LVdffXWeeuqpjBo1qnheAAAAYHgb8h9Ae8cddyRJRo0alZaWlj73tLS0ZOTI997KnXfeecj6G2+8ke9+97tJkosuuqhXSDmgvr4+119/fZJkw4YNefTRR4/I/AAAAMDwMqRjSnd3d9asWZMkmT17dqZMmdLnvilTpmT27NlJko6OjnR3d/daX716deWxA7cF9WXBggWV41/+fBYAAACAZIjHlM7OzuzevTtJMmvWrH73HljfvXt3nnnmmV5r69evP2RfXz71qU9l4sSJSfKBP98FAAAAGN6O6memfFgbN26sHE+bNq3fvQevd3Z25qyzzjrkPCNHjswnP/nJfs/zqU99Kq+88ko2btyYnp6ejBgxYkCz/vI3BP2yA5/XAgAAAFS3IR1Ttm/fXjmeMGFCv3sPXt+xY0ef5/nYxz6W2traAZ1n7969eeutt1JXVzegWSdPnjygfQAAAEB1G9K3+Ry4xSdJxo0b1+/eg9d37drV53ne7xzvdx4AAACAIX1lSjXZsmVLv+tdXV2ZMWPGMZoGAAAAOFqGdEw5+BabPXv29Lv34PWTTjqpz/O83zne7zz9mTRp0oD3AgAAANVrSN/mM378+Mrxtm3b+t178Pqpp57a53l27tx5yNcmH+48Y8eOzYknnlg0LwAAADD8DemYcsYZZ1SON23a1O/eg9cbGhr6PM/+/fvzwgsv9Hue559/vvKcgX6TDwAAAHD8GNIxpaGhoXKLzhNPPNHv3gPrdXV1Of3003utHfxZJf2d5/nnn88rr7ySJJk+ffoHmhkAAAAY3oZ0TKmtrc2cOXOSJGvXrs3mzZv73Ld58+asXbs2STJ37txDvv64sbGx8tiKFSsO+3orV66sHF966aUfanYAAABgeBrSMSVJrrvuuiTJu+++myVLlvS5Z8mSJdm/f3+S5Nprrz1k/aMf/Wi+9KUvJUkeffTR/OAHPzhkz8svv5zbb789SfLrv/7rueiii47I/AAAAMDwMuRjynnnnZdFixYlSVpbW9Pc3JzOzs7s3bs3nZ2daW5uTmtra5Jk0aJFmTlzZp/nueWWW3LaaaclSX7nd34nS5cuzY4dO7J79+6sWbMmn/3sZ/Paa69l9OjR+du//duMGjXq2LxBAAAAoKqM6Onp6TlaJ7/88stz3333DXj/Cy+8kKlTpx7yeHd3dxYsWJD29vbDPnfevHlZuXJlampqDrvnpz/9aZqamrJ169Y+18eNG5fvfOc7aW5uHvDMA7V169ZMnjw5SbJly5aq+irlqTetHuwRqsqLtzUO9ggAAADHrWPx9/eQvzIlee+zU9ra2tLW1pbGxsZMnDgxNTU1mThxYhobG9PW1pZVq1b1G1KS5Oyzz87Pfvaz3HrrrfnMZz6Tj370oxk3blw+/elP5ytf+UqeeuqpoxJSAAAAgOHjqF6Zwv9xZcrxw5UpAAAAg+dY/P09+oifEY5z4tPACU8AAEA1qorbfAAAAACGCjEFAAAAoICYAgAAAFBATAEAAAAoIKYAAAAAFBBTAAAAAAqIKQAAAAAFxBQAAACAAmIKAAAAQAExBQAAAKCAmAIAAABQQEwBAAAAKCCmAAAAABQQUwAAAAAKiCkAAAAABcQUAAAAgAJiCgAAAEABMQUAAACggJgCAAAAUEBMAQAAACggpgAAAAAUEFMAAAAACogpAAAAAAXEFAAAAIACYgoAAABAATEFAAAAoICYAgAAAFBATAEAAAAoIKYAAAAAFBBTAAAAAAqIKQAAAAAFxBQAAACAAmIKAAAAQAExBQAAAKCAmAIAAABQQEwBAAAAKDB6sAcAAGD4mHrT6sEeoaq8eFvjYI8AwAfgyhQAAACAAmIKAAAAQAExBQAAAKCAmAIAAABQQEwBAAAAKCCmAAAAABQQUwAAAAAKiCkAAAAABcQUAAAAgAJiCgAAAEABMQUAAACggJgCAAAAUEBMAQAAACggpgAAAAAUEFMAAAAACogpAAAAAAXEFAAAAIACYgoAAABAATEFAAAAoICYAgAAAFBATAEAAAAoIKYAAAAAFBBTAAAAAAqIKQAAAAAFxBQAAACAAmIKAAAAQAExBQAAAKCAmAIAAABQQEwBAAAAKCCmAAAAABQYPdgDAAAMdVNvWj3YIwAMmP9mDdyLtzUO9ghUKVemAAAAABQQUwAAAAAKiCkAAAAABcQUAAAAgAJiCgAAAEABMQUAAACggJgCAAAAUEBMAQAAACggpgAAAAAUEFMAAAAACogpAAAAAAXEFAAAAIACowd7AAA4kqbetHqwR6gaL97WONgjAMCg8v8NA+f/G3pzZQoAAABAATEFAAAAoICYAgAAAFBATAEAAAAoIKYAAAAAFBBTAAAAAAqIKQAAAAAFxBQAAACAAmIKAAAAQAExBQAAAKCAmAIAAABQYPRgDwAcv6betHqwR6gaL97WONgjAAAA/58rUwAAAAAKiCkAAAAABcQUAAAAgAJiCgAAAEABMQUAAACggJgCAAAAUEBMAQAAACggpgAAAAAUqKqY8uabb+b222/PBRdckE984hOpra3NiSeemGnTpuWLX/xiVq9ePaDztLe3p6mpKfX19RkzZkzq6+vT1NSUBx544Ci/AwAAAKDajR7sAQZq/fr1+a3f+q288sorvR7ft29fnn/++Tz//PNpbW1NU1NT7r///pxwwgmHnKO7uzsLFixIe3t7r8e7urrS0dGRjo6OzJs3LytXrkxNTc1RfT8AAABAdaqKK1P+93//N42NjZWQsnjx4jz11FPZvXt3XnnllTz88MM5++yzkyQPP/xwvvrVr/Z5nquuuqoSUhYuXJgNGzZkz5492bBhQxYuXJgkaWtry+LFi4/BuwIAAACqUVXElOXLl2fHjh1JkmuuuSZLly7NWWedlRNPPDETJkzI3Llz89hjj2Xy5MlJkn/+53/Oa6+91usc69aty7Jly5Ik8+fPz/Lly9PQ0JCxY8emoaEhy5cvz/z585Mky5Yty5NPPnkM3yEAAABQLaoipmzcuLFyvGDBgj73nHzyybnkkkuSJO+8806ee+65Xut33HFHkmTUqFFpaWnp8xwtLS0ZOfK9X8mdd975oecGAAAAhp+qiCkTJkz4UM/p7u7OmjVrkiSzZ8/OlClT+nzOlClTMnv27CRJR0dHuru7P8C0AAAAwHBWFTGlsbExI0aMSJLcf//9fe7ZuXNnvv/97ydJzjrrrF7BpLOzM7t3706SzJo1q9/XOrC+e/fuPPPMMx96dgAAAGB4qYqYMn369Nx8881J3rv95qtf/WqefvrpvPXWW3n11VezZs2aXHjhhdmyZUtOPfXULFu2rBJfkt63CU2bNq3f1zp4vbOz8wi/EwAAAKDaVc1XI998880599xz8+1vfzt/93d/l7vuuqvX+sc//vFcc801ufHGGzNp0qRea9u3b68cv98tQwevH/jQ24HYunVrv+tdXV0DPhcAAAAwdFVNTNm/f382btyY5557Lj09PYesv/HGG3n++eezefPmQ2LKgVt8kmTcuHH9vs7B67t27RrwfAe+SQgAAAAY3qoipuzduze/+7u/m46OjowaNSo33HBDFi1alGnTpmXv3r158sknc8stt6SjoyOPPPJI7r333lx22WWDPTYAAHCETL1p9WCPAFBRFTHlG9/4Rjo6OpIk99xzT7785S9X1mpra/P5z38+F1xwQS688ML8+Mc/zpe//OXMmDEjv/Zrv5Ykqaurq+zfs2dPv6918PpJJ5004Bm3bNnS73pXV1dmzJgx4PMBAAAAQ9OQjyk9PT35+7//+yTvfTjsokWL+txXU1OTv/iLv8j555+fffv25d57701LS0uSZPz48ZV927Zt6/f1Dl4/9dRTBzznL99aBAAAAAxPQ/7bfF599dXKB8iec845vb6l55dNnz69cvyzn/2scnzGGWdUjjdt2tTv6x283tDQUDwvAAAAMLwN+ZgycuQHG/Hg5zU0NFRu9XniiSf6fd6B9bq6upx++ukf6LUBAACA4WvIx5RTTjmlEkL++7//u89v8jlg/fr1leOpU6dWjmtrazNnzpwkydq1a7N58+Y+n7958+asXbs2STJ37tzU1tZ+2PEBAACAYWbIx5SRI0dWQsimTZuybNmyPvft27cvf/Znf1b5ee7cub3Wr7vuuiTJu+++myVLlvR5jiVLlmT//v1JkmuvvfZDzw4AAAAMP0M+piTJzTffnBNOOCFJcuWVV+bGG2/Mxo0bs2/fvrz55pv593//95x//vn58Y9/nCT53Oc+l9/8zd/sdY7zzjuv8uG1ra2taW5uTmdnZ/bu3ZvOzs40NzentbU1SbJo0aLMnDnzGL5DAAAAoFoM+W/zSd77zJMHH3wwCxYsyI4dO/Ktb30r3/rWt/rc+7nPfS6rVq3qc23p0qXZuXNn2tvbs2LFiqxYseKQPfPmzcvdd999ROcHAAAAho+quDIlSS666KI888wz+eu//uucf/75GT9+fGpqajJu3Lh88pOfzO/93u/le9/7Xh599NGcfPLJfZ6jtrY2bW1taWtrS2NjYyZOnJiamppMnDgxjY2NaWtry6pVq1JTU3OM3x0AAABQLariypQDTjnllCxZsuSwn3kyUJdeemkuvfTSIzQVAAAAcDypmitTAAAAAIYCMQUAAACggJgCAAAAUEBMAQAAACggpgAAAAAUEFMAAAAACogpAAAAAAXEFAAAAIACYgoAAABAATEFAAAAoICYAgAAAFBATAEAAAAoIKYAAAAAFBg92AMA8P6m3rR6sEcAAAD+P1emAAAAABQQUwAAAAAKiCkAAAAABcQUAAAAgAJiCgAAAEABMQUAAACggJgCAAAAUEBMAQAAACggpgAAAAAUEFMAAAAACogpAAAAAAXEFAAAAIACYgoAAABAATEFAAAAoICYAgAAAFBATAEAAAAoIKYAAAAAFBBTAAAAAAqMHuwBAADgeDX1ptWDPQIAH4ArUwAAAAAKiCkAAAAABcQUAAAAgAJiCgAAAEABMQUAAACggJgCAAAAUEBMAQAAACggpgAAAAAUEFMAAAAACogpAAAAAAXEFAAAAIACYgoAAABAATEFAAAAoICYAgAAAFBATAEAAAAoIKYAAAAAFBBTAAAAAAqMHuwBAIDBMfWm1YM9AgBAVXJlCgAAAEABMQUAAACggJgCAAAAUEBMAQAAACggpgAAAAAUEFMAAAAACogpAAAAAAXEFAAAAIACYgoAAABAATEFAAAAoICYAgAAAFBATAEAAAAoIKYAAAAAFBBTAAAAAAqIKQAAAAAFxBQAAACAAmIKAAAAQAExBQAAAKCAmAIAAABQQEwBAAAAKCCmAAAAABQQUwAAAAAKiCkAAAAABcQUAAAAgAJiCgAAAEABMQUAAACggJgCAAAAUEBMAQAAACggpgAAAAAUEFMAAAAACogpAAAAAAXEFAAAAIACYgoAAABAATEFAAAAoICYAgAAAFBATAEAAAAoIKYAAAAAFBBTAAAAAAqIKQAAAAAFxBQAAACAAmIKAAAAQAExBQAAAKCAmAIAAABQQEwBAAAAKCCmAAAAABQQUwAAAAAKiCkAAAAABcQUAAAAgAJiCgAAAEABMQUAAACggJgCAAAAUEBMAQAAAChQlTHlf/7nf/K1r30tZ599dk455ZTU1tZm0qRJ+e3f/u0sW7Ysu3fv7vf57e3taWpqSn19fcaMGZP6+vo0NTXlgQceOEbvAAAAAKhWowd7gBI9PT35+te/npaWluzbt6/X2ksvvZSXXnopDz74YDZv3pw///M/P+T53d3dWbBgQdrb23s93tXVlY6OjnR0dGTevHlZuXJlampqjuZbAQAAAKpUVV2ZcsUVV+Sb3/xm9u3blwsuuCAPP/xwtm/fnu7u7nR1deXBBx/MwoULU1dX1+fzr7rqqkpIWbhwYTZs2JA9e/Zkw4YNWbhwYZKkra0tixcvPmbvCQAAAKguI3p6enoGe4iBuPfee3PFFVckSa655prceeedRc9ft25dZs2alSSZP39+7r///kP2fPGLX0xra2tl/8yZMz/k1P9n69atmTx5cpJky5YtmTRp0hE799E29abVgz0CAAAAg+jF2xoHe4QBOxZ/f1fFlSlvv/12brzxxiTJZz7zmdx+++3F57jjjjuSJKNGjUpLS0ufe1paWjJy5Hu/ktJYAwAAABwfqiKmrFixIq+//nqS5IYbbsioUaOKnt/d3Z01a9YkSWbPnp0pU6b0uW/KlCmZPXt2kqSjoyPd3d0fYmoAAABgOKqKmPLwww8nee+qki984QvFz+/s7Kx8w8+BW30O58D67t2788wzzxS/FgAAADC8VcW3+fzXf/1XkqShoSF1dXV56KGHcs899+QnP/lJ3njjjZxyyimZPn16/uAP/iCXXnppRowY0ev5GzdurBxPmzat39c6eL2zszNnnXXWgGbcunVrv+tdXV0DOg8AAAAwtA35mPKLX/wi27ZtS5JMnjw5l19+ee67775ee7q6uvLQQw/loYceysUXX5x//dd/zUc+8pHK+vbt2yvHEyZM6Pf1Dl7fsWPHgOc88OE2AAAAwPA25G/zOfBZKUnyox/9KPfdd1+mTZuWVatW5fXXX8+ePXuydu3aXHjhhUmSRx55JM3Nzb3OceAWnyQZN25cv6938PquXbuOxFsAAAAAhpEhf2XK/v37K8e/+MUvMmHChKxdu7bXFSSzZs3KI488ks9//vN57LHH0tHRkR/+8If5jd/4jWM255YtW/pd7+rqyowZM47RNAAAAMDRMuRjykknndTr5+uvv77PW3VGjRqVv/zLv8zMmTOTJK2trZWYUldXV9m3Z8+efl/v4PVffu3+HI3vrQYAAACGniF/m09dXV3Gjh1b+fn8888/7N7p06fnhBNOSJI89dRTlcfHjx9fOT7w+SuHc/D6qaeeWjouAAAAMMwN+ZgyYsSINDQ0VH7++Mc/fti9I0eOzMc+9rEkyRtvvFF5/Iwzzqgcb9q0qd/XO3j94NcFAAAASKogpiSp3LqTJK+99tph9+3fvz87d+5MkkpUSf7vK5WT5Iknnuj3tQ6s19XV5fTTT//AMwMAAADDU1XElHnz5lWOH3/88cPuW79+fd5+++0kyTnnnFN5vLa2NnPmzEmSrF27Nps3b+7z+Zs3b87atWuTJHPnzk1tbe2Hnh0AAAAYXqoiplxwwQWVOPI3f/M3fX7uybvvvps//dM/rfx82WWX9Vq/7rrrKvuWLFnS5+ssWbKk8u1B11577RGZHQAAABheqiKmjBgxIkuXLs3YsWOzbdu2zJ49O+3t7XnzzTezd+/erFu3LhdffHEee+yxJMnll1+eWbNm9TrHeeedl0WLFiV575t+mpub09nZmb1796azszPNzc1pbW1NkixatKjXrUUAAAAABwz5r0Y+YPr06Wlra0tzc3M2bdrU69afg1122WW55557+lxbunRpdu7cmfb29qxYsSIrVqw4ZM+8efNy9913H9HZAQAAgOGjKq5MOWDOnDnp7OzMn/zJn+TMM8/MRz7ykYwZMyZTpkzJggUL8qMf/Sj/9E//dNjPOqmtrU1bW1va2trS2NiYiRMnpqamJhMnTkxjY2Pa2tqyatWq1NTUHON3BgAAAFSLET09PT2DPcTxYOvWrZk8eXKSZMuWLZk0adIgTzRwU29aPdgjAAAAMIhevK1xsEcYsGPx93dVXZkCAAAAMNjEFAAAAIACYgoAAABAATEFAAAAoICYAgAAAFBATAEAAAAoIKYAAAAAFBBTAAAAAAqIKQAAAAAFxBQAAACAAmIKAAAAQAExBQAAAKCAmAIAAABQQEwBAAAAKCCmAAAAABQQUwAAAAAKiCkAAAAABcQUAAAAgAJiCgAAAEABMQUAAACggJgCAAAAUEBMAQAAACggpgAAAAAUEFMAAAAACogpAAAAAAXEFAAAAIACYgoAAABAATEFAAAAoICYAgAAAFBATAEAAAAoIKYAAAAAFBBTAAAAAAqIKQAAAAAFxBQAAACAAmIKAAAAQAExBQAAAKCAmAIAAABQQEwBAAAAKCCmAAAAABQQUwAAAAAKiCkAAAAABcQUAAAAgAJiCgAAAEABMQUAAACggJgCAAAAUEBMAQAAACggpgAAAAAUEFMAAAAACogpAAAAAAXEFAAAAIACYgoAAABAATEFAAAAoICYAgAAAFBATAEAAAAoIKYAAAAAFBBTAAAAAAqIKQAAAAAFxBQAAACAAmIKAAAAQAExBQAAAKCAmAIAAABQQEwBAAAAKCCmAAAAABQQUwAAAAAKiCkAAAAABcQUAAAAgAJiCgAAAEABMQUAAACggJgCAAAAUEBMAQAAACggpgAAAAAUEFMAAAAACogpAAAAAAXEFAAAAIACYgoAAABAATEFAAAAoICYAgAAAFBATAEAAAAoIKYAAAAAFBBTAAAAAAqIKQAAAAAFxBQAAACAAmIKAAAAQAExBQAAAKCAmAIAAABQQEwBAAAAKCCmAAAAABQQUwAAAAAKiCkAAAAABcQUAAAAgAJiCgAAAEABMQUAAACggJgCAAAAUEBMAQAAACggpgAAAAAUEFMAAAAACogpAAAAAAXEFAAAAIACYgoAAABAATEFAAAAoICYAgAAAFBATAEAAAAoUNUx5cYbb8yIESMq//7jP/7jfZ/T3t6epqam1NfXZ8yYMamvr09TU1MeeOCBoz8wAAAAUPVGD/YAH9QTTzyR22+/fcD7u7u7s2DBgrS3t/d6vKurKx0dHeno6Mi8efMsQrjkAAAgAElEQVSycuXK1NTUHOlxAQAAgGGiKq9Mefvtt/OHf/iH2b9/f+rr6wf0nKuuuqoSUhYuXJgNGzZkz5492bBhQxYuXJgkaWtry+LFi4/a3AAAAED1q8qY8rWvfS3PPfdcZs6cmUWLFr3v/nXr1mXZsmVJkvnz52f58uVpaGjI2LFj09DQkOXLl2f+/PlJkmXLluXJJ588qvMDAAAA1avqYsrjjz+eu+66K2PHjs13v/vdjBo16n2fc8cddyRJRo0alZaWlj73tLS0ZOTI934dd95555EbGAAAABhWqiqm7Nq1K1/60pfS09OTW2+9Naeffvr7Pqe7uztr1qxJksyePTtTpkzpc9+UKVMye/bsJElHR0e6u7uP3OAAAADAsFFVMeX666/Pz3/+85x77rm5/vrrB/Sczs7O7N69O0kya9asfvceWN+9e3eeeeaZDzcsAAAAMCxVTUz5wQ9+kH/4h3/ImDFjBnx7T5Js3Lixcjxt2rR+9x683tnZ+cEGBQAAAIa1qvhq5J07d+aKK65Iktx6660544wzBvzc7du3V44nTJjQ796D13fs2FE049atW/td7+rqKjofAAAAMDRVRUy5+uqr89JLL+Xcc8/NH//xHxc998AtPkkybty4fvcevL5r166i15k8eXLRfgAAAKA6DfnbfL73ve/lX/7lX4pv7wEAAAA4Gob0lSk7duzIlVdemSS55ZZbim7vOaCurq5yvGfPnn73Hrx+0kknFb3Oli1b+l3v6urKjBkzis4JAAAADD1DOqZ885vfzKuvvpoZM2bkhhtu+EDnGD9+fOV427Zt/e49eP3UU08tep1JkyaVDQYAAABUpSEdU15//fUkyfr16zN69PuPesEFF1SOX3jhhUydOrXX1SybNm3q9/kHrzc0NJSOCwAAABwHhvxnpnxYDQ0NlVt9nnjiiX73Hlivq6vL6aefftRnAwAAAKrPkI4p//iP/5ienp5+/918882V/Y8//njl8alTpyZJamtrM2fOnCTJ2rVrs3nz5j5fa/PmzVm7dm2SZO7cuamtrT26bw4AAACoSkM6phwp1113XZLk3XffzZIlS/rcs2TJkuzfvz9Jcu211x6z2QAAAIDqclzElPPOOy+LFi1KkrS2tqa5uTmdnZ3Zu3dvOjs709zcnNbW1iTJokWLMnPmzMEcFwAAABjChvQH0B5JS5cuzc6dO9Pe3p4VK1ZkxYoVh+yZN29e7r777kGYDgAAAKgWx8WVKcl7n53S1taWtra2NDY2ZuLEiampqcnEiRPT2NiYtra2rFq1KjU1NYM9KgAAADCEjejp6ekZ7CGOB1u3bs3kyZOTJFu2bMmkSZMGeaKBm3rT6sEeAQAAgEH04m2Ngz3CgB2Lv7+PmytTAAAAAI4EMQUAAACggJgCAAAAUEBMAQAAACggpgAAAAAUEFMAAAAACogpAAAAAAXEFAAAAIACYgoAAABAATEFAAAAoICYAgAAAFBATAEAAAAoIKYAAAAAFBBTAAAAAAqIKQAAAAAFxBQAAACAAmIKAAAAQAExBQAAAKCAmAIAAABQQEwBAAAAKCCmAAAAABQQUwAAAAAKiCkAAAAABcQUAAAAgAJiCgAAAEABMQUAAACggJgCAAAAUEBMAQAAACggpgAAAAAUEFMAAAAACogpAAAAAAXEFAAAAIACYgoAAABAATEFAAAAoICYAgAAAFBATAEAAAAoIKYAAAAAFBBTAAAAAAqIKQAAAAAFxBQAAACAAmIKAAAAQAExBQAAAKCAmAIAAABQQEwBAAAAKCCmAAAAABQQUwAAAAAKiCkAAAAABcQUAAAAgAJiCgAAAEABMQUAAACggJgCAAAAUEBMAQAAACggpgAAAAAUEFMAAAAACogpAAAAAAXEFAAAAIACYgoAAABAATEFAAAAoICYAgAAAFBATAEAAAAoIKYAAAAAFBBTAAAAAAqIKQAAAAAFxBQAAACAAmIKAAAAQAExBQAAAKCAmAIAAABQQEwBAAAAKCCmAAAAABQQUwAAAAAKiCkAAAAABcQUAAAAgAJiCgAAAEABMQUAAACggJgCAAAAUEBMAQAAACggpgAAAAAUEFMAAAAACogpAAAAAAXEFAAAAIACYgoAAABAATEFAAAAoICYAgAAAFBATAEAAAAoIKYAAAAAFBBTAAAAAAqIKQAAAAAFxBQAAACAAmIKAAAAQAExBQAAAKCAmAIAAABQQEwBAAAAKCCmAAAAABQQUwAAAAAKiCkAAAAABcQUAAAAgAJiCgAAAEABMQUAAACggJgCAAAAUEBMAQAAACggpgAAAAAUEFMAAAAACogpAAAAAAWGfEx5+eWXc9999+XKK6/Mueeem0984hOpqalJXV1dPv3pT+f3f//388Mf/rDonO3t7Wlqakp9fX3GjBmT+vr6NDU15YEHHjhK7wIAAAAYLkYP9gDvZ9asWfn5z39+yOPvvPNOnnvuuTz33HNZvnx5mpqasnz58px00kmHPVd3d3cWLFiQ9vb2Xo93dXWlo6MjHR0dmTdvXlauXJmampoj/l4AAACA6jfkr0w5YPr06fn2t7+dzs7OvPXWW9m1a1ceffTRfPazn02SPPzww7n00kv7PcdVV11VCSkLFy7Mhg0bsmfPnmzYsCELFy5MkrS1tWXx4sVH980AAAAAVWvIx5Rzzz03//mf/5n169fnK1/5Ss4444yccMIJqaury4UXXpjHHnssF198cZLk0UcfzUMPPdTnedatW5dly5YlSebPn5/ly5enoaEhY8eOTUNDQ5YvX5758+cnSZYtW5Ynn3zy2LxBAAAAoKoM+ZjS2tpaufqkL6NGjcpf/dVfVX4+XEy54447KvtbWlr63NPS0pKRI9/7ldx5550fdGQAAABgGBvyMWUgGhoaKsdbt249ZL27+/+1d+8xXpV3/sA/X4GBkeFakMEwrDpqlKZeQlBXU6NoutGp2oKrFbUqbcGgdA2yvcR2qUgbGxMTusVbbF1XIzWVXsSBul7obgKNmNWSlZFGBRW5swoz0GGm03l+f/Cbs0OZ28EZ5jvM65VM8jDP5zzzfMkHZs57zjnfxlixYkVERFx00UUxceLENteZOHFiXHTRRRER8cILL0RjY2MP7BYAAADoy46JMGXr1q3ZeMSIEYfN19TUxL59+yLi4ANtO9Iyv2/fvtiwYUM37hIAAAA4FhT9u/l0xdKlS7NxW7cEvf3229m4srKyw7Vaz9fU1MRZZ53VpT20dUVMa9u2bevSOgAAAEBx6/NhyubNm+P++++PiIjy8vK4+eabD6vZtWtXNh43blyH67We3717d5f3UVFR0eVaAAAAoO/q07f5HDhwIK699tqoq6uLQqEQjz76aAwfPvywupZbfCIiSktLO1yz9XxdXV33bRYAAAA4JvTZK1Oam5vjq1/9aqxduzYiIu677764+uqre20/mzdv7nB+27Ztcd555x2l3QAAAAA9pU+GKSmlmDVrVvzyl7+MiIgf/OAHcc8997RbX1ZWlo3r6+s7XLv1/LBhw7q8pwkTJnS5FgAAAOi7+lyYklKKOXPmxM9+9rOIiHjggQdi/vz5HR4zduzYbLxjx44Oa1vPjxkz5lPsFAAAADgW9bkwZe7cufHII49EoVCIn/70pzFnzpxOjznzzDOz8Xvvvddhbev5SZMmHflGAQAAgGNSn3oA7V133RVLlizJHjbblSAl4mAo0nKrz5o1azqsbZkvKyuLM84449NtGAAAADjm9Jkw5e67747FixdHoVCIxx57LL7xjW90+diSkpK48sorIyJi9erV8eGHH7ZZ9+GHH8bq1asjIuKLX/xilJSUfPqNAwAAAMeUPhGmfOc734kHH3wwCoVCPPzww/H1r3899xp33XVXRET89a9/jW9961tt1nzrW9+K5ubmiIj4p3/6pyPfMAAAAHDMKvow5V/+5V/ixz/+cURELF68OGbPnn1E6/z93/99zJw5MyIinn322bjxxhujpqYmDhw4EDU1NXHjjTfGs88+GxERM2fOjAsuuKB7XgAAAABwTCmklFJvb6IjhUIhV31lZWW8++67bc41NjbGDTfcEL/61a/aPX769OmxdOnSGDRoUK6v25mPPvooKioqIiJi8+bNfeqtlE/6TnVvbwEAAIBe9P79Vb29hS47GuffRX9lSncqKSmJZcuWxbJly6KqqirKy8tj0KBBUV5eHlVVVbFs2bJ47rnnuj1IAQAAAI4dRf/WyD1x4cy0adNi2rRp3b4uAAAAcOzrV1emAAAAAHxawhQAAACAHIQpAAAAADkIUwAAAAByEKYAAAAA5CBMAQAAAMhBmAIAAACQgzAFAAAAIAdhCgAAAEAOwhQAAACAHIQpAAAAADkIUwAAAAByEKYAAAAA5CBMAQAAAMhBmAIAAACQgzAFAAAAIAdhCgAAAEAOwhQAAACAHIQpAAAAADkIUwAAAAByEKYAAAAA5CBMAQAAAMhBmAIAAACQgzAFAAAAIAdhCgAAAEAOwhQAAACAHIQpAAAAADkIUwAAAAByEKYAAAAA5CBMAQAAAMhBmAIAAACQgzAFAAAAIAdhCgAAAEAOwhQAAACAHIQpAAAAADkIUwAAAAByEKYAAAAA5CBMAQAAAMhBmAIAAACQgzAFAAAAIAdhCgAAAEAOwhQAAACAHIQpAAAAADkIUwAAAAByEKYAAAAA5CBMAQAAAMhBmAIAAACQgzAFAAAAIAdhCgAAAEAOwhQAAACAHIQpAAAAADkIUwAAAAByEKYAAAAA5CBMAQAAAMhBmAIAAACQgzAFAAAAIAdhCgAAAEAOwhQAAACAHIQpAAAAADkIUwAAAAByEKYAAAAA5CBMAQAAAMhBmAIAAACQgzAFAAAAIAdhCgAAAEAOwhQAAACAHIQpAAAAADkIUwAAAAByEKYAAAAA5CBMAQAAAMhBmAIAAACQgzAFAAAAIAdhCgAAAEAOwhQAAACAHIQpAAAAADkIUwAAAAByEKYAAAAA5CBMAQAAAMhBmAIAAACQgzAFAAAAIAdhCgAAAEAOwhQAAACAHIQpAAAAADkIUwAAAAByEKYAAAAA5CBMAQAAAMhBmAIAAACQgzAFAAAAIAdhCgAAAEAOwhQAAACAHIQpAAAAADkIUwAAAAByEKYAAAAA5CBMAQAAAMhBmAIAAACQgzAFAAAAIAdhCgAAAEAOwhQAAACAHIQpAAAAADkIUwAAAAByEKYAAAAA5CBMAQAAAMihX4Yp77zzTsydOzdOP/30GDp0aIwaNSomT54cixYtitra2t7eHgAAAFDEBvb2Bo62p59+OmbNmhX19fXZ5/785z/HG2+8EW+88UY89thjsXz58jj77LN7cZcAAABAsepXV6asWrUqbrvttqivr4/Kysqorq6O2tra2LlzZyxZsiSOP/742Lx5c1xxxRWxc+fO3t4uAAAAUIT6TZjy17/+Ne68885oamqK0aNHx3/913/FlVdeGcOGDYuxY8fGnDlz4rnnnouIiG3btsWCBQt6eccAAABAMeo3Ycp//Md/RE1NTUREzJs3L0488cTDaq644oq4/PLLIyLiiSee8PwUAAAA4DD9Jkz59a9/nY1nzJjRbt0NN9wQERENDQ1RXV3d4/sCAAAA+pZ+E6asXbs2IiLKy8vj5JNPbrfuwgsvzMavv/56j+8LAAAA6Fv6RZiSUooNGzZERERlZWWHtaecckoUCoWIiOy2IAAAAIAW/eKtkevq6qKhoSEiIsaNG9dhbUlJSYwcOTI++eST2L17d5e/xkcffdTh/ObNm7Pxtm3burxuMWiq7frfAwAAAMeezs55i0nrc+6mpqYe+Rr9IkzZt29fNi4tLe20vrS0ND755JOoq6vr8teoqKjocu15553X5VoAAADobRUP9/YOjsyuXbvipJNO6vZ1+8VtPgAAAADdpV9cmVJWVpaN6+vrO61vqRk2bFiXv0br23jacuDAgdiwYUOMGzcuxo4dGwMHFv9f/bZt27KraNauXRvjx4/v5R3Bp6OnOZboZ441eppjiX7mWNIX+7mpqSl27doVERGf+9zneuRrFP8ZfTcYNmxYDB48OBoaGmLHjh0d1jY2NsaePXsiImLMmDFd/hoTJkzotObUU0/t8nrFZvz48V16jdBX6GmOJfqZY42e5liinzmW9KV+7olbe1rrF7f5FAqFOOOMMyIi4r333uuwduPGjZFSioiISZMm9fjeAAAAgL6lX4QpEf/30Nft27fHpk2b2q1bs2ZNNp4yZUqP7wsAAADoW/pNmPLlL385Gz/zzDPt1i1dujQiIgYPHhxVVVU9vi8AAACgb+k3YcoXvvCF7LadBx98MLZu3XpYzcqVK+Pll1+OiIjbbrsthg8fflT3CAAAABS/fhOmDBgwIP71X/81Bg4cGB9//HFcfPHFsWLFiti3b1/s3r07Hn744fjHf/zHiDj4UJ177723l3cMAAAAFKN+8W4+LaZOnRpPPPFEzJo1K9577702b+OpqKiI5cuXxwknnNALOwQAAACKXb+5MqXFTTfdFOvWrYs77rgjTjvttCgtLY0RI0bEueeeGwsXLoy33norzj777N7eJgAAAFCkCqnlfYABAAAA6FS/uzIFAAAA4NMQpgAAAADkIEwBAAAAyEGYAgAAAJCDMAUAAAAgB2EKAAAAQA7CFAAAAIAchCkAAAAAOQhTAAAAAHIQptCmd955J+bOnRunn356DB06NEaNGhWTJ0+ORYsWRW1tbW9vj35u69at8eSTT8bs2bPj/PPPjxNOOCEGDRoUZWVlcdppp8VNN90UL730Uq41f/WrX8VVV10VJ554YgwePDhOPPHEuOqqq+LXv/51D70K6Jp//ud/jkKhkH38/ve/7/QY/Uyx+dOf/hTf/va345xzzonPfOYzUVJSEhMmTIgvfelL8fOf/zz27dvX4fF6mmJQW1sbDz74YFx66aVxwgknRElJSQwdOjQqKyvjK1/5SlRXV3dpHf1MT9q6dWs8//zz8f3vfz+uuOKKGDt2bPYzxK233pp7ve7o19ra2li0aFFMnjw5Ro0aFUOHDo3TTz895s6dG++8807uPRWNBH/jqaeeSqWlpSki2vyoqKhIf/zjH3t7m/Rjf/d3f9duf7b+uOqqq1JtbW2HazU0NKRp06Z1uM706dNTY2PjUXp18H9Wr16djjvuuEP6cdWqVe3W62eKTXNzc7rnnnvSoEGDOuzLBQsWtHm8nqZYvPbaa6m8vLxLP3vs37+/zTX0M0dDR/11yy23dHmd7urXN998M02YMKHdNUpLS9PTTz/9KV917xCmcIhXX301DRw4MEVEqqysTNXV1am2tjbt3LkzLVmyJB1//PEpItL48ePTjh07enu79FMtYcqUKVPST3/601RTU5P279+f6urq0ssvv5wuvvji7D/oyy+/vMO1Zs6cmdXOmDEjrV+/PtXX16f169enGTNmZHMzZ848Sq8ODtq/f3869dRTU0SkE088sUthin6m2LTuyUsvvTQtX7487dq1KzU2NqZt27al3/72t2nGjBnpgQce6PR4PU1v2b17dxozZkzWb7fffntat25d2rdvX9q+fXtavnx5Ouecc7L52267rc119DNHw/jx49PVV1+dFi5cmFauXJlef/31IwpTuqNft2/fnsaPH58iIg0dOjQtWbIk7dy5M9XW1qbq6upUWVmZIiINHDgwvfrqq93w6o8uYQqZpqamNGnSpBQRafTo0WnLli2H1axYseKQbyTQG6677rr0n//5n+3ONzU1pX/4h3/IevW3v/1tm3Vr1qzJaq6//vo2a66//vqs5g9/+EO37B+64s4770wRkS644IL0ve99r9MwRT9TbB5//PGs3775zW/mPl5PUywWL17caS9//PHHqaKiIjsx/N///d9D5vUzvWXTpk25w5Tu6tfZs2dnNStWrDhsfsuWLWn06NEpItJnP/vZ1NTU1OXXVQyEKWRaByWLFi1qt+7yyy9PEZEGDx6c9u7dexR3CF33xhtvZP38ta99rc2a6667LkVEGjBgQPrggw/arPnggw+y2yy+8pWv9OSWIfPqq6+mQqGQhgwZkt5+++20YMGCTsMU/Uwx2b9/fxo1alSKiHTuuece0Q/Ieppicfvtt3cp5Jg1a1ZW99prrx0yp5/pLUcSpnRHv+7ZsyeVlJR0eqX4okWLsv397ne/69L+ioUH0JJp/RChGTNmtFt3ww03REREQ0NDlx+0BUfbpEmTsvFHH3102HxjY2OsWLEiIiIuuuiimDhxYpvrTJw4MS666KKIiHjhhReisbGxB3YL/6euri5uu+22SCnFwoUL44wzzuj0GP1MsXnmmWfik08+iYiI+fPnx4ABA3Idr6cpJuPGjftUx+hn+pLu6tfq6ursc105t4w4+LDbvkSYQmbt2rUREVFeXh4nn3xyu3UXXnhhNn799dd7fF9wJLZu3ZqNR4wYcdh8TU1N9u4RrXu6LS3z+/btiw0bNnTjLuFw8+bNiw8++CDOP//8mDdvXpeO0c8Um+XLl0dExIABA+Lqq6/OfbyepphUVVVFoVCIiIhf/OIXbdbs2bMnVq5cGRERZ5111iEnoPqZvqS7+rXl3LKzdU455ZQoLy+PiL53bilMISIiUkrZP4DKysoOa0855ZTsG0pNTU2P7w2OxNKlS7PxxRdffNj822+/nY076/nW83qenvS73/0uHn/88Rg8eHA88cQTXf5tvn6m2Lz22msRcfAqwbKysnj++eejqqoqTjjhhOytNa+55ppYtmxZpJQOO15PU0ymTJkSCxYsiIiIxYsXx5133hn/8z//E/v374+dO3fGihUr4rLLLovNmzfHmDFj4uc//3n2s3KEfqZv6a5+bVnnuOOO6/AX9REHzy9bjmnre0KxEqYQEQcvK29oaIiIzi9lLCkpiZEjR0ZExO7du3t8b5DX5s2b4/7774+Ig1da3XzzzYfV7Nq1Kxt31vOt5/U8PWXPnj3x9a9/PSIiFi5cGGeeeWaXj9XPFJOGhobYsWNHRERUVFTErbfeGtdcc02sWLEidu3aFY2NjbFt27Z4/vnn49prr40rrrgiamtrD1lDT1NsFixYECtXroyqqqp46KGH4qyzzoqysrIYN25cVFVVxfvvvx/f/OY3480334zJkycfcqx+pi/prn5tWWfkyJFRUlLSpXUOHDgQ+/fvz7Xf3iRMISIiu5QrIqK0tLTT+paaurq6HtsTHIkDBw7EtddeG3V1dVEoFOLRRx+N4cOHH1aXp+dbz+t5esrcuXNjy5Ytcf7558fdd9+d61j9TDFpeVZKRMQrr7wSTz75ZFRWVsZzzz0Xn3zySdTX18fq1avjsssui4iIF198MW688cZD1tDTFJvm5uZ4++234913323zN+d79+6NjRs3xocffnjYnH6mL+mufm1ZJ8+5ZVvrFDNhCnDMaG5ujq9+9avZPZr33XffEd2rD0fbb37zm3j66adz394Dxai5uTkbNzQ0xLhx42L16tUxffr0GDlyZAwZMiQuvPDCePHFF2Pq1KkRcfDhhS+99FJvbRk6dODAgbjmmmti3rx58e6778b8+fOjpqYmGhoaYu/evfHiiy/G+eefHy+88EJccskl8dRTT/X2loGjQJhCRESUlZVl4/r6+k7rW2qGDRvWY3uCPFJKMWvWrPjlL38ZERE/+MEP4p577mm3Pk/Pt57X83S33bt3x+zZsyMi4t577811e08L/Uwx+du+mjdvXpuXig8YMCB+9KMfZX9+9tlns7GeppgsWrQoXnjhhYiIePTRR+OBBx6IM888M0pKSmL48OHxhS98IX7/+9/H5z//+fjLX/4SX/va1+JPf/pTdrx+pi/prn5tWSfPuWVb6xQzYQoRcbBpBw8eHBGR3efcnsbGxtizZ09ERIwZM6bH9wadSSnFnDlz4mc/+1lERDzwwAPZg+LaM3bs2GzcWc+3ntfzdLcf/vCHsXPnzjjvvPNi/vz5R7SGfqaYlJWVxZAhQ7I/X3LJJe3WTpkyJY4//viIiFi3bl32eT1NsUgpxWOPPRYRBx+2OXPmzDbrBg0aFPfdd19ERPzlL3/JfiaJ0M/0Ld3Vry3r7Nmzp9O3+W5ZZ8iQITF06NBc++1NwhQiIqJQKMQZZ5wRERHvvfdeh7UbN27M7hWdNGlSj+8NOjN37tx45JFHolAoxJIlS7p0Qtr6t/+d9XzreT1Pd2t5vsTatWtj4MCBUSgUDvu49957s/pLL700+/z7778fEfqZ4lIoFA7prdGjR7dbe9xxx2UPtd+7d2/2eT1Nsdi5c2f2IM3Jkycf8i49f2vKlCnZ+K233srG+pm+pLv6tWWd5ubm2LRpU4frbNy4MTumo39jxUaYQua8886LiIjt27d32PBr1qzJxq2/aUBvuOuuu2LJkiXZw2bnzJnTpeNa3q4z4tCebkvLfFlZWRY6QjHRzxSbCy64IBt//PHH7dY1NzdnV7u2hCoRepricdxxR3a61Po4/Uxf0l392nJu2dk6GzdujO3bt0dE3zu3FKaQ+fKXv5yNn3nmmXbrli5dGhERgwcPjqqqqh7fF7Tn7rvvjsWLF0ehUIjHHnssvvGNb3T52JKSkrjyyisjImL16tVtPn0/IuLDDz+M1atXR0TEF7/4xU7f2g3y+rd/+7dIKXX40fq2tVWrVmWfP+mkkyJCP1N8pk+fno1XrVrVbt3atWvjz3/+c0TEIW8nq6cpFp/5zGeyE8v//u//bvOdfFq0PAA/IrL/nyP0M31Ld/VrVVVV9rmunFtGREybNu1T7f2oS/D/NTU1pUmTJqWISKNHj05btmw5rGbFihUpIlJEpNtvv70XdgkHffvb304RkQqFQnrkkUeOaI01a9Zk/Xz99de3WXP99ddnNX/4wx8+zZbhiC1YsCDrw1WrVrVZo58pJs3NzUFn1xwAAAPHSURBVGny5MkpItK4cePS9u3bD6tpampKU6dOzXpy9erVh8zraYrFddddl/XZ448/3mZNY2Nj+vznP5/VrVy58pB5/Uxv2bRpU9ZXt9xyS5eO6a5+nT17drv/JlJKacuWLWn06NEpItJnP/vZ1NTU1OXXVQyEKRzilVdeSQMHDkwRkSorK1N1dXWqq6tLu3btSg899FAaOnRoiog0fvz4tGPHjt7eLv3U97///ew/5p/85Cefaq2ZM2dma82YMSOtX78+1dfXp/Xr16cZM2ZkczNnzuym3UN+XQlTUtLPFJe1a9emIUOGZD9TLFu2LO3duzfV19enNWvWpMsuuyzryVtvvbXNNfQ0xWD9+vXp+OOPTxGRBgwYkObPn59qampSY2Nj2rt3b3rxxRfThRdemPXj1KlT21xHP9MbjiRMSal7+nX79u1p/PjxKSLS0KFD00MPPZR27dqV6urqUnV1daqsrEwRkQYOHJheffXVbni1R5cwhcM89dRTqbS0NPsH8rcfFRUV6Y9//GNvb5N+rL3ebO+jsrKy3bUaGhrStGnTOjx++vTpqbGx8Si+QjhUV8MU/Uyxqa6uTiNHjuywJ2+++ebU0NDQ5vF6mmLx0ksvpTFjxnT6M8fUqVPTxx9/3OYa+pmj4ZZbbsn1c/KmTZvaXKe7+vXNN99MEyZMaHeN0tLS9PTTT/fA30TP88wUDnPTTTfFunXr4o477ojTTjstSktLY8SIEXHuuefGwoUL46233oqzzz67t7cJ3aKkpCSWLVsWy5Yti6qqqigvL49BgwZFeXl5VFVVxbJly+K5556LQYMG9fZWoVP6mWJz5ZVXRk1NTXz3u9+Nz33uczF8+PAYPHhwTJw4MW644YZ45ZVX4t///d/bfTaEnqZYXH755bFhw4b48Y9/HJdcckmMHTs2Bg0aFKWlpXHyySfHddddF7/5zW/i5ZdfjlGjRrW5hn6mL+mufj3nnHPirbfeioULF8a5554bI0aMiNLS0jjttNPijjvuiHXr1sWNN954lF5V9yqk1MFTlAAAAAA4hCtTAAAAAHIQpgAAAADkIEwBAAAAyEGYAgAAAJCDMAUAAAAgB2EKAAAAQA7CFAAAAIAchCkAAAAAOQhTAAAAAHIQpgAAAADkIEwBAAAAyEGYAgAAAJCDMAUAAAAgB2EKAAAAQA7CFAAAAIAchCkAAAAAOQhTAAAAAHIQpgAAAADkIEwBAAAAyEGYAgAAAJCDMAUAAAAgB2EKAAAAQA7CFAAAAIAchCkAAAAAOQhTAAAAAHL4f65+tFEQVKmOAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 414,
       "width": 553
      }
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(X, bins='auto');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use a custom MLE here, without exogenous variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "class TwoPeakLLH(GenericLikelihoodModel):\n",
    "    \"\"\"Fit height of signal peak over background.\"\"\"\n",
    "    start_params = [10, 1000]\n",
    "    cloneattr = ['start_params', 'signal', 'background']\n",
    "    exog_names = ['n_signal', 'n_background']\n",
    "    endog_names = ['alpha']\n",
    "    \n",
    "    def __init__(self, endog, exog=None, signal=None, background=None,\n",
    "                 *args, **kwargs):\n",
    "        # assume we know the shape + location of the two components,\n",
    "        # so we re-use their PDFs here\n",
    "        self.signal = signal\n",
    "        self.background = background\n",
    "        super(TwoPeakLLH, self).__init__(endog=endog, exog=exog, \n",
    "                                         *args, **kwargs)\n",
    "\n",
    "    def loglike(self, params):        # pylint: disable=E0202\n",
    "        return -self.nloglike(params)\n",
    "\n",
    "    def nloglike(self, params):\n",
    "        endog = self.endog\n",
    "        return self.nlnlike(params, endog)\n",
    "\n",
    "    def nlnlike(self, params, endog):\n",
    "        n_sig = params[0]\n",
    "        n_bkg = params[1]\n",
    "        if (n_sig < 0) or n_bkg < 0:\n",
    "            return np.inf\n",
    "        n_tot = n_bkg + n_sig\n",
    "        alpha = endog\n",
    "        sig = self.signal.pdf(alpha)\n",
    "        bkg = self.background.pdf(alpha)\n",
    "        sumlogl = np.sum(\n",
    "            np.ma.log(\n",
    "                (n_sig * sig) + (n_bkg * bkg)\n",
    "            )\n",
    "        )\n",
    "        sumlogl -= n_tot\n",
    "        return -sumlogl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = TwoPeakLLH(endog=X, \n",
    "                   signal=pdf_a, \n",
    "                   background=pdf_b, \n",
    "                   )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: -1.342181\n",
      "         Iterations: 75\n",
      "         Function evaluations: 143\n"
     ]
    }
   ],
   "source": [
    "res = model.fit()\n",
    "_ = res.bootstrap()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>TwoPeakLLH Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>         <td>['alpha']</td>     <th>  Log-Likelihood:    </th> <td>  1382.4</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                <td>TwoPeakLLH</td>     <th>  AIC:               </th> <td>     nan</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>           <td>Maximum Likelihood</td> <th>  BIC:               </th> <td>     nan</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>              <td>Tue, 22 May 2018</td>  <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                  <td>01:08:48</td>      <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>       <td>  1030</td>       <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>           <td>   NaN</td>       <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>               <td>   NaN</td>       <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "        <td></td>          <th>coef</th>     <th>std err</th>      <th>z</th>      <th>P>|z|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>n_signal</th>     <td>   31.7378</td> <td>    8.912</td> <td>    3.561</td> <td> 0.000</td> <td>   14.270</td> <td>   49.206</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>n_background</th> <td>  998.2622</td> <td>   32.341</td> <td>   30.867</td> <td> 0.000</td> <td>  934.875</td> <td> 1061.650</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                              TwoPeakLLH Results                              \n",
       "==============================================================================\n",
       "Dep. Variable:              ['alpha']   Log-Likelihood:                 1382.4\n",
       "Model:                     TwoPeakLLH   AIC:                               nan\n",
       "Method:            Maximum Likelihood   BIC:                               nan\n",
       "Date:                Tue, 22 May 2018                                         \n",
       "Time:                        01:08:48                                         \n",
       "No. Observations:                1030                                         \n",
       "Df Residuals:                     NaN                                         \n",
       "Df Model:                         NaN                                         \n",
       "================================================================================\n",
       "                   coef    std err          z      P>|z|      [0.025      0.975]\n",
       "--------------------------------------------------------------------------------\n",
       "n_signal        31.7378      8.912      3.561      0.000      14.270      49.206\n",
       "n_background   998.2622     32.341     30.867      0.000     934.875    1061.650\n",
       "================================================================================\n",
       "\"\"\""
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}