docs/examples/ex_bivariate_likelihood.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "source": "# Example: Likelihood fitting a Bivariate Gaussian\n\nIn this example, we shall perform likelihood fitting to a [bivariate normal\ndistribution](http://mathworld.wolfram.com/BivariateNormalDistribution.html),\nto demonstrate how ``symfit``\u0027s API can easily be used to perform likelihood\nfitting on multivariate problems.\n\nIn this example, we sample from a bivariate normal distribution with a\nsignificant correlation of $\\rho \u003d 0.6$ between $x$ and $y$.\nWe see that this is extracted from the data relatively straightforwardly.\n",
      "metadata": {
        "pycharm": {
          "metadata": false,
          "name": "#%% md\n"
        }
      }
    },
    {
      "cell_type": "code",
      "execution_count": 13,
      "outputs": [],
      "source": "import numpy as np\nfrom symfit import Variable, Parameter, Fit\nfrom symfit.core.objectives import LogLikelihood\nfrom symfit.distributions import BivariateGaussian\nimport matplotlib.pyplot as plt\n",
      "metadata": {
        "pycharm": {
          "metadata": false,
          "name": "#%%\n",
          "is_executing": false
        }
      }
    },
    {
      "cell_type": "markdown",
      "source": "Build a model corresponding to a bivariate normal distribution.",
      "metadata": {
        "pycharm": {
          "metadata": false,
          "name": "#%% md\n"
        }
      }
    },
    {
      "cell_type": "code",
      "execution_count": 14,
      "outputs": [],
      "source": "x \u003d Variable(\u0027x\u0027)\ny \u003d Variable(\u0027y\u0027)\nx0 \u003d Parameter(\u0027x0\u0027, value\u003d0.6, min\u003d0.5, max\u003d0.7)\nsig_x \u003d Parameter(\u0027sig_x\u0027, value\u003d0.1, max\u003d1.0)\ny0 \u003d Parameter(\u0027y0\u0027, value\u003d0.7, min\u003d0.6, max\u003d0.9)\nsig_y \u003d Parameter(\u0027sig_y\u0027, value\u003d0.05, max\u003d1.0)\nrho \u003d Parameter(\u0027rho\u0027, value\u003d0.001, min\u003d-1, max\u003d1)\n\npdf \u003d BivariateGaussian(x\u003dx, mu_x\u003dx0, sig_x\u003dsig_x, y\u003dy, mu_y\u003dy0,\n                       sig_y\u003dsig_y, rho\u003drho)\n",
      "metadata": {
        "pycharm": {
          "metadata": false,
          "name": "#%%\n",
          "is_executing": false
        }
      }
    },
    {
      "cell_type": "markdown",
      "source": "Generate mock data\n",
      "metadata": {
        "pycharm": {
          "metadata": false,
          "name": "#%% md\n"
        }
      }
    },
    {
      "cell_type": "code",
      "execution_count": 15,
      "outputs": [],
      "source": "# Draw 100000 samples from a bivariate distribution\nmean \u003d [0.59, 0.8]\ncorr \u003d 0.6\ncov \u003d np.array([[0.11 ** 2, 0.11 * 0.23 * corr],\n                [0.11 * 0.23 * corr, 0.23 ** 2]])\nnp.random.seed(42)\nxdata, ydata \u003d np.random.multivariate_normal(mean, cov, 100000).T\n",
      "metadata": {
        "pycharm": {
          "metadata": false,
          "name": "#%%\n",
          "is_executing": false
        }
      }
    },
    {
      "cell_type": "code",
      "execution_count": 16,
      "outputs": [
        {
          "data": {
            "text/plain": "\u003cFigure size 432x288 with 1 Axes\u003e",
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAD8CAYAAABw1c+bAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJztvXuQXNd13vut7nk/McAAJAiABCVBD4q2JBtF2nGuRceyTKsUMam4bPLGN1auFFa5rNwq+95U5MotKUWnKk5St1xJRY7EOCxZqbKkxIlsJKFNy7F1pSubDmHJsUTSlEDwgSEI4jGDwbxnunvdP2Yw+1ur+5w509Pzwlm/KhS7++yzzz67ew6/vfZ6iKoiCIIgKA+V3R5AEARBsLPEgz8IgqBkxIM/CIKgZMSDPwiCoGTEgz8IgqBkxIM/CIKgZMSDPwiCoGTEgz8IgqBkdG3UQESeAPBBAJdV9d4Wx/8BgL9N/b0DwGFVnRSRlwHMAKgDqKnq6U4NPAiCIGgP2ShyV0R+GMAsgM+1evC7tn8dwC+o6l9be/8ygNOqenUzg+qRXu3D4GZOCYIgKDUzmLqqqoeLtN1Q8avqV0XkZMFrPwLg8wXbZtKHQdwvP7rVboLg1kectVYbuzOOvUDJ5+IP9LdeKdq2YzZ+ERkA8CCA/0QfK4DfF5E/E5FHNzj/URE5KyJnV7DUqWEFQRAEjg0V/yb46wC+rqqT9NkPqepFETkC4Msi8peq+tVWJ6vq4wAeB4ARORiZ44KgCCVTtbnEXBSmk149D8OZeVT14tp/LwP4EoD7Oni9IAiCoA068uAXkVEA7wXwO/TZoIgM33wN4P0Avt2J6wVBsEmkYv8FpaaIO+fnATwAYFxEJgB8EkA3AKjqp9ea/U0Av6+qc3TqbQC+JCI3r/Obqvp7nRt6EARB0A5FvHoeKdDmswA+6z47D+Bd7Q4sCIIOEvbvgIg1XxAEQcmIB38QBEHJiAd/EARByYgHfxAEQcnoZABXEAT7EXbv3MlN4JKnWNhNQvEHQRCUjFD8QVB2tlNp56n6vOvmBZm1M97dWtXsUULxB0EQlIxQ/EGwk+w3u/ZWlXK7qr6dPvLadWKeb6FVQyj+IAiCkhGKPwj2CFKtrr/Wen2bL5ahXttdkbTTXzt9ezrdR56q3+cqnwnFHwRBUDJC8QfBTpKjGrdd5RehXVWbdV4nPHc6obSL9rFbHk47TCj+IAiCkhEP/iAIgpIRpp4g2G6KugF22l2wqLtkQZNLO5vPfM5qd8pvCp+X2UceGf03jamoiW07XVt3mFD8QRAEJSMUfxDsFQoq4KKqudKd/ry9qjXvC64MtnzOJuB7lIqYY/y+rTEVXTE0n9jeeXuQDWdKRJ4Qkcsi0rJQuog8ICLTIvLna/8+QcceFJEXROSciHy8kwMPgiAI2qOI4v8sgH8N4HM5bb6mqh/kD0SkCuBTAH4MwASAZ0TkjKo+1+ZYg2BrbNVG224gUQeUYju28cZKrWjnhfrLU8pGhberqDPmt2nVQO1y9x2yvu+d3GfZbsx4i5+2oeJX1a8CmGxjSPcBOKeq51V1GcAXADzURj9BEARBB+mUjf8HReR/ArgI4P9S1WcBHANwgdpMALi/Q9cLgs2z19P5Fk0XkLPyYLXNNv7G8nL2ZdvwoPF293Zs7blkeBP5a2lj89eq9PSY92ZuOpG+Yidpc0ydePB/A8BdqjorIh8A8NsATgGQFm0zFyMi8iiARwGgDwMdGFYQBEHQii0/+FX1Br1+UkR+TUTGsarwT1DT41hdEWT18ziAxwFgRA62aSQMgjbYjmIhWf21u0/Adu08ezr10REvnIxViFStas6067vr8irEY5Q3nZdr46e58GPKWuXkrX4Mt3Aq5y2vyUTkdhGRtdf3rfV5DcAzAE6JyN0i0gPgYQBntnq9IAiCYGtsqPhF5PMAHgAwLiITAD4JoBsAVPXTAH4SwM+JSA3AAoCHVVUB1ETkYwCeAlAF8MSa7T8Idp9OKzHXh7VLF/Q8yVCyq90nRc19NNu/2xlv6qMy4Mys9dSOlXKuxxD13WRPp/Oa/PMLehDZVQ19nBczUHBuO+Hxs6Pptdtkwwe/qj6ywfF/jVV3z1bHngTwZHtDC4IgCLaDSNkQBEFQMiJlQ1AetmjeyXUDdJuYeSkHTJ8Zm53tJkHLCuBqMqt09aV2dB+6sGjH5+55/XM3bjteGlOXbdc1PLz+un79ek4fNNaCSdVyvx8+P2dDPJeC7faqeYcJxR8EQVAyQvEH+5+iSr6N+qx5ar3oJh5vQOap17z+qkND6VieOyIp7ErG5wCgS0upHSllcUq+sbCQxjCYNn7Vbe4K3VbeXLDKb1LeNO95bp/S20vjSyuUJoVftA5wlqvrbgbsbTOh+IMgCEpGKP5g/9NptVTQXZCVcmcSmKX+K7195ggrbLNK6OrOvi6Nj5U7YFWzVChYyin5Cq80VlbSgarVjFKhR4nZW3DaksdbWzGH+F5YvXv133D7EHQx+74d232Ge+2m+uv0CnQbCMUfBEFQMkLxB7vLdto58xRbh8sheo8SJTWblzohy2umWVFXMo8ZaqS2yYOm4vcWWL3zvkC/W2mQuq6Qjb9xY9Zet5uUPHsJuRUEK/nqyLA5prQqySsik0lOEJ3pI69d0WCuTYyjo3QocVwo/iAIgpIRij/YXbazKEaH7atN6YvZ+8cp5cY0KdvR0fT5rFXK0pds7azCxa8ElOICeBxODcvoSDpl2drQTbvujL0Btwchw8nGD7a1z81n9m1w88wrCnX7DtkJ4XysAg2XVhC5hd15DDlpJDLHk8dO2uo7dK1Q/EEQBCUjHvxBEAQlI0w9wc6zS4ErWZWmvEkgK2jLu05ysJNPdeA3Lm/C7pEAmkw1Wf2xSUgG+jPbsakmr53QRi0a9B0481Bjer3cBipkRmoiYwO7MjRom1EAV+58svuqc/s0VcbYpdabvdiNNm9/uB03zbzP91De/SxC8QdBEJSMUPzBzrPVpFhtumkWTZyW1c4rT7MB2d9vj0k6rzEzs/66On7IdkGbpBXeSO117qHsPkmq3Ch3d13eEDavAWAxpWxQcgE15wOo3H4kvWHXTn+/Zj4pjYLbwK2Qi6lfhXAqCjOfPl4rw1XWr9yyXEJzUzuYARV0DuhEjeEdZv+NOAiCINgSofiDvUtBJcXpDRpLSR5m2fSB4ikWuA+v6o07olPKbK82tuyGCx5ybqCpnRvf0aS8+Yh4t0pKxQBKxCZHxm27SbK1HxpLn9eca+OMC9S6SY9zBzWBY7RyceebBHM+SK0nYy5yq5txMxcsxqkj8qp7tbOyzKOdPawd3vcKxR8EQVAyQvEH2882ez2wyjeXyUmqJtmLAQjVna2TV0t1wKl6DgTyNnRul6XqAYDVMXvUuP5kIal3dFGxlTHraVOZo7lg9e9hxb5EKtwHdg2SVw6neZDsPRKzanDpoIWO+XQOZjVQz7ah5yVwy8Kmb1jOPGbG025BlU4HG24DGyp+EXlCRC6LyLczjv9tEfmLtX9/LCLvomMvi8i3ROTPReRsJwceBEEQtEeR/11+FqvF1D+XcfwlAO9V1SkR+QkAjwO4n47/iKpe3dIog/2NVzNFvSiK2l7Zbzwr6RmyS/FVvX86qW1WlL5IifG8ccqW7dfSS+q9z6pw5f45EZvbCwAr1l5KqjbnVjusxNmGXnf9jZJ3TReNtebaLdKcUVwArk8ji/rkVBrCkcPmmJLNv2klRPOutbR/wiswAACv8HKKt7CNP+u7B/KT6G0ru1iIZcMHv6p+VURO5hz/Y3r7NIDjWx9WEARBsF102sb/EQC/S+8VwO+LiAL4jKo+nnWiiDwK4FEA6MNAVrNgL7HdhS/MpVp72jSdl+NfnbUaaFLy7PXBnjy+qAircm8bZzs3qfz6ARvJWr00md4M0e9+ynnDHDmYXrNCv27bLb/ltjTcpXrLcwCgspDGV50lv/teex+VGqnh2bn02s8lrTQ4VkHnXQGYA5SwjlYGAFAZO9CyP95nAdz3T6uaxkJOScqs87EJW36nSzTuYoRvxx78IvIjWH3w/1X6+IdU9aKIHAHwZRH5S1X9aqvz1/6n8DgAjMjBvDJFQRAEwRboiEFLRL4XwK8DeEhVr938XFUvrv33MoAvAbivE9cLgiAI2mfLil9E7gTwnwH8b6r6Hfp8EEBFVWfWXr8fwGNbvV6wh2gzl35uMI3ponWtWnFmC06lYIKl/HX9JuFNvAnnQNr4FDZV+A1C3pxccekceOOXTEKVxew+Gv009l4bcFU7QO2603hrd9pkcAOvpvQQvFFbH7SbykIbqXzdyoK7D049QeamypQ1vxg4Qd2Su19O0+A2nBtTFFRWyf79VIaTSYjTYTRVD6MANm1QXWFv2mkn7/4+TNPAbPjgF5HPA3gAwLiITAD4JIBuAFDVTwP4BIBDAH5NVn8kNVU9DeA2AF9a+6wLwG+q6u9twz0EQRAEm6CIV88jGxz/KICPtvj8PIB3NZ8RlA6njuzmWj3jc+eCl1PtSmdI8bN7pAtgsgqQ1JtPU0AJzITTCg+7lMpMjwsyGk2rC63SxvSSVZvLd6Y+G7SS6b1qN0W1izZPF9N4KytWDS/ensbbe42SqjklL6R650+mDdeBV9xGKqnyykxKD6EjdpNarrV272zaSKXvoHLAutE2yNWTvx9OybF6MGOV6L7HrE3bpgpcOa6emeyD1Mt57O/1ShAEQbBpImVDsDGstk1xi2w3OON+6QpucIoFVl++9ikrvbyiJ3bPoEHtXEpgchfUq8mNUlwAlx4il0Nyb6zccCqcAqlkyY19Oqnj2uFk8148YRO99b+e2jW6033MvMmuLvomk2JfOJJWMgMX3Vz0pu9q9s6kygcutU5rAQA915PiXTxqr9v/alLy2pfmQpzrpAnM4lTOzsYvfdnFXDjVM9vufXAcf69FXTH5N+J/ZzadQ5tpGvYZofiDIAhKRij+YGPInmlK2OV4NpgCGb4kXobCaiqOwtdllebtq12cVoFWFz6RGAX7VG4jrxln8105lJRy1w0Kbhqxar0yl+zVy8dGzTG269f60/1WF+3YZ+9K1xqcSOq/e962q/elPlj9L43bfQxW70reP/U+l86gO42vfyKp65VhVzKSbPz1sbRv0eWDpViVc6poV1AGFbL5+5TNnEaaPLB03qWepmtVChbeyUvLkOdZlskuplvoBKH4gyAISkYo/qB9CvrxeyXP6ssUUXHeFZV+SqNMJQaVvG6aLsuJxLxHCfVhCph02XZV8oBZPpSdPqRxJF1rZdD20X+FUgevJEW5PGLbLY+Qfzmp/5VBO2f1nvR+8BJ/7jymhtOKp0GeQFKz/XXNp+9g+p1p76P/ivX+WbozHetmT6OaW8Xx9zOS9gl8oRjl5G55qax5VeiOmZVBnu1+q8XW81KS7EOVz4TiD4IgKBmh+INNUbSEnUmRm1eMwxcw5+64dJ5PU8z9sa+9Ke3nkrmRn7yOJlW6ctAXLE8vKxT9ynZ7AGj0pP6XRq2Gqq4k5b08TNf1ed6oz6WR9LpSt3PL5wlNRXXZ+fEfTHPdM5uOsU1/dUxpNTU4kfYxVkatTb77Rvp+OMJXe613TvVSSrgmPlqXMCsDvxqg784XaTF9UHwGR/hW3V5N3ZelzCCrdOetTCj+IAiCkhEP/iAIgpIRpp5gU5iNWZf0rEFud7xR65fgbMTgoB3NMedwMjNf7dWYDI5S8JBL5rZCG7VdM2l8tQFrEloZIvfLZdqYHXbmnKV0bPiCNW8sHElmkYVDacS9160JZzkjnkka9i6r1P30ya6WnwNA32Saw8UDZAZZVtcufY9Lh5LppO8Na+pg00/PVDpWncp2sUSOe6Q37zAmHz/9Fho3rNsnu5jylXjTd/WDjKRqbmM2M2VDXuW42NwNgiAI9hOh+MvKJqpnZYa0++RmWQFdfpOV4WAc3+4gpVi4sl7mwWwQAi5pG7kI1oZscFNlkRKTHU8bwn7DtbpEQUuUAqHhC3WZzV57sFJL45g9mT5fnrRKfui11O7629LnfVdsuwrHr+VM5/KIXw+twquT1fek+MdSh8tjdnXGid5kJZ3TGHBzezWlVLZBW3awTdXTCHbTNS67VfsFNShlR24KkYzfY24FrjxVv89VPhOKPwiCoGSE4i8rBRU+YBVRdYjSCPskaC5dcuqwtQoFXL1bbxtmd04KCsKcS5Z2JNV4bQwnJdrodcFSo64u7hpdLj3CymBWbVX7dvFguq+FcXuPQiJy/Juk/u+wfVy9n1JbVFK76oId6/yxdGzwFU4HYfur9XP65vT5wmE7F0OvpXscmkjqemHcuXPOpPNWDqeL9b/i0jBzqgz+Hn2BGg7m8wVwODCLk/Ll7P2YAj3ud8vJAfPchstIKP4gCIKSEYo/WIVsm02BWXSMVX6evdbg9gIqnAiMPX66Wyty364xZl1hZDn1vzyWVh31XhdwRcNlD53uOdtu8WA61rWY5mLhsOuP85K9xSrK/lfSvcwf5yM+6C297JpOHXbZRQ26aIzL5PzSO2Xb8V5A9xwdcM40vK9Rp5VRz4y1k9cG0pj6Jij18qwLdGJlzzZ+v9rjxGnL1gsnMxGfs9VzcGBeMkCzGiiYTvxWsuPnEYo/CIKgZBRS/CLyBIAPArisqve2OC4A/iWAD2BVW3xYVb+xduxnAfzfa03/iar+RicGXno67VPM6RZ8abqs8HlfIIPL6pGqN0U1AFRpL4D9umX8oGnHicAag3ROb45nCE0LJykDgLmjnB45fe596WdJofNqYPGwVeuVpXSsa9KlpXh3KmFYfTZdoD5g+5DBNLddb6Q+Fm637Rrd6fupzpN//qK9xxplrxi8RDEIQ9l7EIsH07z0X3GJzshLSqlQjA77zQXy8qH4CVl2tvUeWtU5v3vps55C69fy+wSsVwuq99wkbSWkqOL/LIAHc47/BIBTa/8eBfBvAEBEDmK1OPv9AO4D8EkRGWt3sEEQBMHWKaT4VfWrInIyp8lDAD6nqgrgaRE5ICJHATwA4MuqOgkAIvJlrP4P5PNbGXSAtlR+rg91TqGKpgIpN6nbMXAULtv1q0cO2/PYm4NUflP6XVb5pDaXDlnvIU5a1qA0xT7SdplqpfCaZuZue90q2cPZjl8btWq4i9R6fdlea+lC8kLSu9PyQufsn1y1mu568RgVja+7Oae0yvVRikdw6Zt7r1D5xjvTsb4rtrve6dQHz1PNpZcWShbXRY487NO/Oij6EmZoc8FH03LRHFe8vn45DdL8VnNiTkwhnzz//DxuoYjconTKxn8MwAV6P7H2WdbnQRAEwS7RKa+eVpJQcz5v7kDkUayaidCH7OIXQRAEwdbo1IN/AsAJen8cwMW1zx9wn3+lVQeq+jiAxwFgRA62UQQzAJBfB5fc2/wGrlkWa7YJx+TId5txvAlXYfOOr63Ky3Ny/auP23qvQmkFGlRzttaffY/zt1E7t1fI7pw9M+Tq5ywCnBJh8e2UHmDaupv2HEz3u7BkzU+9J9I9Hx5Jr189f8S0q82nP8Hxo8mWMjNv+1u+Ru9XKG//ip2Len+6r76rpLu85aivtcvqyoAzHd2g38xSul91CfCEzTvd2Ru4Jvna9RvmUIWCA9UFB5ou6DfIv+OiKRtu5VQMRemUqecMgL8jq/wAgGlVfR3AUwDeLyJja5u671/7LAiCINglirpzfh6ryn1cRCaw6qnTDQCq+mkAT2LVlfMcVt05/+7asUkR+WUAz6x19djNjd6gg+So/MJdZCViA9Agl8vqCKnynFQMJmhr0Lr+cV1c3rSt91tF3Rjh9Au0ATlg73futjQODlpitQrYClrXT6XXjWNWXTaWSfKzyh+ym7u3jybFOnLY7p7+xStpK2uhL/XRc9AGPi1dT/d49fW0+1yZdxvxdG2+j65JV+v3anrNwVzewGqC2ygr8eBFm6K463oa79JtSZH3feeSaYdh8iNdzK7Axe67Fa6cBkC5Mhun/x5yrqOURjnT1RjYXiW/zzeEi3r1PLLBcQXw8xnHngDwxOaHFgRBEGwHkbJhL9MBJZ+FLz7BYfBeRbEd1ah1Z7s37pwcpNXrgpsWk614+VgKbqrO2+tOvzmp4QPfTcrTB2YtHKJEYuQhOPVmq8SUgqCE3Cj7+ux169Ru5Lak6ucX7b7I20Yvr7/+vefeaY69467X118/fz5lZqv02NXUfe94af31n72ctskaFSfRlzj4jOzzLhXDErmsDl+gdBMuiRwXcOmZTe04DTUAVGlF1nOFLjbgVPg0BemRHV9dAJdx4XTHGrOu4EoWWQrb/71kqPImt8+swjF5Sn4fqnwmUjYEQRCUjFD8m2UTBUy23H9e31ntcsbHSodT1gKuVOKI9a4xtnyyw/qCKGA1x94/Va/EUn8cILQyYhV1/7U09ht3pxWEL4jC6nWBnWbc9A2MJVt+rZbm4sRBm+ns7SNJyX/7+u3rr3uqVq0/fenk+uv3veMvkUXtZLr/yQU7Zy9Pp0D2Ro0T5bn9ky6+mdRuxX1VA2mhgcUDqY/+q1bVjr6Y5mLxcFpZdd/ITl9cmaYNFB/Ux0F5/BtxK4PGJBVscb8L43XGK1C3EsjbjypCpsJfPbjp/vYjofiDIAhKRij+zbLdiqCd/knlN6dXIJXfS14yC9a7xPj1N6W3pTFR0i7h4tgA0EVKrJ/8q3tcMjdS/Oyh0ui2Y58fp/5IosyctJfllA1Vcht3Lu6Yn6RCImPp/q8vWFU60Z3ua+JaUuT3HLWeLKeGXR4E4sXZ8fXXP33H2fXXX75m9wIW62lurlxOBvqhQ9Z4PzdFeytk/2+4WIUKpXYYuEqrPef8wumWef7qg/a76r6eJrQxlrxwKtesDz6rd6Xfli+iYmI/nFeYkM+/2YPyaZm54A+d06Tks/6W8oqol4Ty3XEQBEHJiQd/EARByQhTz2bZ7s3domReNycDJ+fLdykbuL+mWroUSi+cT73hAr166Rjv9Q3Za1UoFUOdUjHMH7FjX0zWEnTRvqLPaT/4ajIZzHwvmQicS+SbjiXTTF2zNc/sShrv33vn19dfX1i09QL+x9U7119//6EL5thyPd3LFy+eXn+9VLd/cpen0+5spTvNy8J5V2WsK91Lz3XKx+8zItBUL42kdiOv2E3buTvSdzXyEtW6deYS4TQN/P363w+59rJJ0ads0EU3YL4Wbdqya3Ddbe425p0P63oHOTo2z2miJBu6TCj+IAiCkhGKf7MUVQftrgzovKzaokC2S5tvV+m1yb6yxiMD5GboNn7NhhxX3erK/vk0KE2BDwpaGab7os3dvik7Jg7MqtNt1F3qhJnvSSr11F1pA/aVq66iF/G/HH5x/fV/u2A3XB99y/+3/vqnh7+7/vofXnyfaXfvWPKdfPfgq3ZMtTTgmZX0+u7Bq6bdwaNpKfPrz/2V9dcNtxnboE17Vvk+wVyFFjy9N8gd9i7rvjvyMrnvzlFyPZ9ugTZgK1fIFdP9RriWLldcE78yILTmgrvoN21WnXlKnn/Ht3DAVacJxR8EQVAyQvFvlSzbYY7CyA1AYVt7jkubdHGwDwdLORt/xjgqzhVTl9K1KgedmyatPFChtMdHRk2z6mxSgXUK3BE3BK7xusIVpJw3ntB7rp5VvWF/tvWRdP+1RpqX973pBdPuldm0AvhP5969/vo3T/+6aferr78frfind3zZvH/s0gPrr//kxpvNsR8YSSsKtvG/Pm/LfS0sWyV+k5WDLm3GAv1mWIXnpJRWmtrh12x/1SWqsUwqvzFgFXp1koO2qPMD9rtvvJZWP5XDaXOmcdXmZGSV3xREuEQpsPl3vJ1BktvR/z4gFH8QBEHJCMW/VdpQC7kh40yGvR+w6qioNwPb6jk9LuDSL/hUDJx0iwqnVKddOuPB1H+tPym22aN27PMpCwK6Flq/BoClFDuF2p3pfqsXbdRS/3i6l65KGuvz129DFmND6ZxrdZse+PtGkr3+VE/aM/jGklW5Hxz75vrrCyuHzLE/mnr7+uvLs8kr6u2HbNDXN6+n9M0rSxRUVXUppbvT+xoNt+4CuDhQiwPiqotuj4hq5NYOpNVZ17Sz3bMtP88jhwuizKZVQl7q5cJJ2XJqQJu/pbzArKJ7ASUhFH8QBEHJCMW/GxRNvkY0FZzg1QCprbx0y1zOTvptmgLjX+2UnRxK0lsoZUPDFU6RFSqpSEKs5hYQ/SR6Fw9lt1s+kuzB3aTyD7/rsmlXJZV/aiR1fnuPTSvwtavJDv9TJ76x/nqkYlXuB4e/vf76OjnGv63Lfm8/fe5vrr++/+DL5thPHX5m/XUFycZ/oNsua5avpu9h8PakgOdftdnXqotk16eveMiGD6BnLo2x/0p2wrXqXOpEaq1Tcqw2pN8jpwNxq1aTpoHt+G7PyZRN7M5+/OQlX9Oiedky/s6a0jK3kehtvxOKPwiCoGSE4t9r5KVYJpoib29+3p/htw9rhzURuIBNqezs/zpIfXKSthWrlBbvSCpVaqmdTwncoEt3zScVeeMe59fdk+bi9vckWztHuwLAsbHkX36qP60GvnHjTtPuo8eTf/7/e+Nt6697K/a6ry6lZcj3D6ZCKc8s2Gjat4+8sf76zb1vmGNfm0n9P/NaGgePFQD6jyR7+Oy1tOThSF0A6Jmi4ivs4u4UOhep6Z5MDVcOusIpWWmUa9nqt3EjrUgq4zZGokL7R42pdI9Scb9hWpF6pd2W8m6jBGIZFb4nFH8QBEHJKPTgF5EHReQFETknIh9vcfxXReTP1/59R0Su07E6HTvTycEHQRAEm2dDU4+IVAF8CsCPAZgA8IyInFHV5262UdVfoPZ/H8B7qIsFVX03blWyzDE+JUIbVYMy3dZa9L9OzX6ulFbB1MRl045Dj9hlPOfTb1ANVu22997oSe854ZoP4Kr1pfuaP5buq/d1a35avjNtMl++nsw7vT12o/tHjnwXrXjvmA3g+s9Xvn/99S/e8dT66y9Nf79pV6XSXSe7r62/fmrqe0y7fjIR9bjcCWfO3Zve0Fd37qXbTbuu62luKXsFeidtrvrFw6mTQ9+idtN2cquL6T2Oo/52AAAgAElEQVSbdzivPgDIMs0hX3jOtmtcSqazyiiZupwTAZuBeEO47syG/Pfif9NtVdYK18y2KKL47wNwTlXPq+oygC8AeCin/SMAPt+JwQVBEASdp8jm7jEA7DQ2AeD+Vg1F5C4AdwP4Q/q4T0TOAqgB+BVV/e02x7q9tBvGvdUNJV9dKCMxW5PrGyunWk6dVOe2uc6SS8ZFgVneTbN6g1z1RtJGb2XWXrd+OPVR608qsub2m1kcNygwqTFqFWCVNneXp1In77znZdPuTydPrr9+K9XLHajYe/ynJ35n/fXXF+5afz3aZVXpHd1pc/KLU+mnPr1s5/LNB9K1np61KRtGB9LG6vR8Oq93xCrqlfmkcqtzFGznfi59V9J8rlBM1ODr9vdX76M+2E9g2bkDvzSRjo2nzWyfYsEk76ON2saMC74ilc8bus0V4RJe8fPvvbD6b2Nzt2128lrbTJEHf6tvLiv09GEAv6VqPG3vVNWLIvImAH8oIt9S1Rf9iSLyKIBHAaAPA/5wEARB0CGKPPgnAJyg98cBXMxo+zCAn+cPVPXi2n/Pi8hXsGr/b3rwq+rjAB4HgBE5WDCnQQfZI4mgOADLBJrkpECWam/mMVM/N6cPk1a3xwXdjCeJyYE/C8fsPkH3bLpWrS/bxr9C3pjVebIHD9qGjXrSHKNHUzDWuWvjpt3sTFoNzB9Nq5WP3/27pt3XFlKCtIsrKRGdXxm8spT6f2MxDfbHx79t2v2Xy+9CFlcujLU+4GRU92y6/57p9PnAG04N08+p5wa51NZtu/7X0uqlQkVUZM6lYhjnyDly53TulyYwy6di5na86mS1nrc3leOuzOflBlztpPLe5yqfKWLjfwbAKRG5W0R6sPpwb/LOEZG3ARgD8Cf02ZiI9K69HgfwQwCe8+cGQRAEO8eGil9VayLyMQBPYbWu3xOq+qyIPAbgrKre/J/AIwC+oGoiSt4B4DMi0sDq/2R+hb2BbjnybIB5idQy2lWGk9qsT99o1RoA0HUkKVRf2k44/XLN2XkJHUqq2avIrqmkIpdvS2NaGbL3VKeVwoHzSZVNn7SKbYVs+ZxGGGNWeetK6n9hKe0fvOeO10y71waSt8nxwSSb/3TuLabd/YPn1l8/v3DH+ut/eOScafd18oy5tJz6Prdok769MZfm4tKrLqBpiEoWXksrMl7hAEC9P81F34tpOVDvtkuDwcut7dyVFec9xmkzplJBlMYhm2CuMkEBZxT0pz41SAYmSSCACu0FNHhlkOfdlpeskM7TRsG/naAwhSJ3VfVJAE+6zz7h3v/jFuf9MYDv8Z8HQRAEu0ekbOgkpFKa7JJss2T/fO+xQIqfVX51xKYp0JWkKDlEvqnAyvWkgIU8d9iLBwCUC10v2zEtHUuql4uAeBYOpWOLY1RsxYUMVBdSu6WjZIdescquZzitAIYHkor85sVjpt0P35W2jL597ej667MXT5h2N96UvGve3p+2qf7J1bebdn3kn//afJrPJ88fN+1GD6Z0C7LklHwl/WkNXUzHlm3WB4yc46Iq6TdSnbdqeJlWV5x8rbpgPavYrm8+vzZtP+AUC5NTLc8BAGEvs5Xsgj9ZKn9TCdEijfKOEWuoIAiCkhGKP4st+uzm2S/NsZwiE0z9xox5z8nYTESlwxRYIV9r9cVWqJSjzFtbe5Vs93PHk2puVO1Yh8infHmQVi4unxzb+IVSHQ8fsCmLl5bTz7O/OynZO2+3ic441fGJ4XTsB2+zSvbF2bQX0iupv69dsT74y7V03XsPpZKCaNj7nXs+ee6M3mOvdX3C2tRv0mfrsKD/Wrr/roU0L3NHrFIefSl9JxyFW/GRtr0Ug9FLY/juy6ZdhSO3ee/HeX4Zm3/O30FW/EnT30He31Uo+x0jFH8QBEHJiAd/EARByQhTTxZZy868wKy8/Pm0FDZViHz+fFp285I5d3PXJ8Li69Jmb2MsmYSaNgHJvFM/YFMTLB5JY+y5kcZ0/S02tYOS6WeZhqt+yijXfNfraZNxsd+6Eq5MpusuDKbNw2tizVTPv/FOtKLnmO1vdiXZnP7bhXTOjWl7v3ccSeair19IQV9YsOaX2lD67m98xwZsVWkPc/B1DriyYzTvyRP68Flrzmr00VzTfPp8/DKdNpxBFdeQZw7kOg29NhiwMTPjm29I5Lvf+4TiD4IgKBmh+DdLm4FZWSqoMW/VeqWXAql66evx59PKgDd3fQCOHiSVT2pQh63KXTmUnR+JU/02KLBoxZ0ySGkGuEbsnPW+RH2QVjKUmGz5ht0F5lQPFUoPdem6Va/vvD1V57q2mAb19Et3m3Z3HkkplucX07Wq3fY7ff1a2hTVSzRPPbZdz7W0Auh2wphrDs9R3NfoS66PG+n76r2WVjW1UbsS7L5K3x1v4M7MwTakY1RNS2t2ExjkltvglcGC3WDPTJzmE6w1igV+BXuDUPxBEAQlIxT/Vsmy8TftEVCAFKdYLprKwQXCcDoHTp4lB20AFyhsn4N21CXj6ppOffj6rMvD6dp9U2lvYOSC7ePGXVQXloSoy3qM7qukSslDsmvQqsY6FXYxSc/67Ornm99NNW2PH0+qvq/fuqWefylJb65pq3Xrpsmpk+u0FVJxOqk2kr7jgTdcem06j9NS1HvttbpnUsN6f/pdsMIHAKECO/Uh/i252skNUvlLpPJ9IBUlA+TU3X7FyDZ/01+7KUnyuIXSHu91QvEHQRCUjFD8naSNoiyVYWsoN6lv65SoqmZtr2BPDPbY8KlzTaAOqcFuqwCZpVFf9KV1u7nbrG4YuJRU9OJYOsl7stQOJlVZnUnjaFx16aWpSEvXdbIvn7RKvrs7XWDipcPpnEl7H3JbOq/aR+kHXrbfwTKlX6iQoB59wU5ElYbh9ztY5Y99JynllUE770JJ1rqm0nfni+HIbDpWfT4VgCmav1x8Sm4OsnJlFBnl1M57sThK0Bah+IMgCEpGKP4i5NkvMxRNpcclQcsIY9f5bB98Dp+XilOv7P/PIffeD3ssKX5h9eYUG/uJ905bO2+tPym9mePpvg49a5Xi5D3pWDfd1vRbrS6tzHNCOEpFfNDGFvS+SqUch2jOLtk9iMYxsoeT540v7KKUcqF+Pc2Tc+rBoW9QigrySKou2/uo96T++q7bYz0zVPScVD578QBAfSjNe3Uy7cdITjI8o96HbQY8pYR97K0jOemW89S7+d3mqfx24l6K9hF0nFD8QRAEJaPcir+oLbJgdK455tqxn3NWQqvVdmTXPkwFVlwhFl0m9X6IPF5cEezKFL3vSeqyMWBXJJyKeWXEHlsYT4q1/ypF7p6y7aq0DbFIlf0OPG/V69S9dI8k1huTzkOFTqv3p++ge9pF0L5G1cf7ye/ceet0X0vzPvwy+bG77Y4G3dYg1XzxHjmjL6YVj7qEdZWVdI/LI6T4X7XFzHlFppevpgMnXfADedQor/CmbIQvclYKBi6Ortn7PVuOwg0VvycJxR8EQVAy4sEfBEFQMm59U0/e5lLRZWiOSYg3cTn5WsO5yJl2dMzn3+ekbTpLFZ5c/nxTW5dcOHXc1n6VBWpH7pyVGRvCz1W2fM3dwYscZMQblbZdrY9cOOlQbdA0w8AEpX0gT1Te6F19T+dcyDZHcEWvBpvR3K+b+1uhMfW6AlQDV9L3uHSA6gifs3NWG0jH/Jx130hte6ZoY3rQpmIwpjjaqFWuiYtst01O1td0jKu++d96rfVvf9sTrHXa1TNcR9uikOIXkQdF5AUROSciH29x/MMickVE/nzt30fp2M+KyHfX/v1sJwcfBEEQbJ4NFb+IVAF8CsCPAZgA8IyInFHV51zTL6rqx9y5BwF8EsBprIqWP1s7N7vIZ6fZZhXglf1NmmqScog8B191u6+AXThZvS5atckufTpPbnvOnRO0omiQ2tReO76ea8n/sjZiVWm9r7Wy7Zm1c1tdZuWdNEXPjNWr199KdWZp+vqumWZYpMVLD7mHVpzIrbNnK33ee9W265tsrZsHX7ff4dxR+n5IAC8fsJvPPdNpIJVFq5QrlOaak6rJtN18B60S9QCtulxKDdD336CVYN7vm1eTjRx3zsLOC52g03+PofLbosi3fB+Ac6p6XlWXAXwBwEMF+/9xAF9W1cm1h/2XATzY3lCDIAiCTlDExn8MwAV6PwHg/hbt/paI/DCA7wD4BVW9kHHusRbnQkQeBfAoAPQhO0XwjpFlOywYzOUVf6U/qUgOka+7Iiq8GtB56m/A18jlzF+kNhvORkspASo8JhektXIiyeuuGzbtQ6M7Xbtaz04SMPWWNPbh19J158ftnLGyZ1fKwYtuBbFAqwZaXfhUEazKhy+k8S2MF3RtzOHAX2YXIqkNkd+nuxQnVZPZ1IfesIrf7N1cuJheu5WbSfNBK4OGc+fMChTMquW8eg6/yUm+lvd3EMp7X1FE8bf6xfi//v8C4KSqfi+APwDwG5s4d/VD1cdV9bSqnu5Gb6smQRAEQQcoovgnAJyg98cBXOQGqsrW2X8L4J/RuQ+4c7+y2UHuCkXTLdOxvCRWXB6RA7iqg9mrG22QavRePWzn5WvV/fjS/3uVC7v02q++QSmQvf8Me6gs3J7SJfj0Az0kjofPpfEtHLJlI/svU2EXCnyqDVidMPR66r/Wl61RRl9O9784ltp5b52Rc+k7mL8jbQx0TVsb//AypX2gVMlSc547l5N6N95TAMB7N34VRpjvkahN2sEbb6+cNB92pVkw3UIeWUo+FP6+pojifwbAKRG5W0R6ADwM4Aw3EJGj9PZDAJ5fe/0UgPeLyJiIjAF4/9pnQRAEwS6xoeJX1ZqIfAyrD+wqgCdU9VkReQzAWVU9A+D/EJEPYdWpYhLAh9fOnRSRX8bq/zwA4DFVnWy6yG7RTsqGnGNsK22y8bPtvkZuKQ3bt3BhDTrmUzbAlMGjlYEfX19rs5lMWdt1N6VpbvRa75XlsTT23smkjmuDtt3B55N30dS9ySfde+EsHKLVBW0n5GQOMIy8atUrl2jsv0ZpjudcO7K7919OCn3psPVi6r+Q5ro2llZalaXs9MWcUgEAMMeJ47iYuUveN9davVe8txelacj10DGdb6MqDxv/vqZQAJeqPgngSffZJ+j1LwH4pYxznwDwxBbGGARBEHSQSNkQBEFQMm79lA15FHVPK2gSqvQmk4EP7JIsV0+/GVvh6lzJXNJwWTfZpFM5Rlss865SF5/D5h2fxZFv15mfskwfXXPWhrN4mHLcz6WN0OELdi4aXWk+Fw5Tndk5e13OCjpwJc1L1QVLVeeT6WP5AJnUvKzh5Knkltp30dW3XUz31f3KlXTAV7Hiamc+4Irf02+h4Uw7DeqDN3Cb4GuvtM70Crh0IEUrZrVDmHb2NaH4gyAISka5FX9BN81clZ+1aZtzLbMacH1X+w+kdhScU7nrhGmHG6TeOUnbmHWd5FqtIEXZGLbqsjaa3ndftQq40ZfucWWIFPqs3WTsu0obvwPkBumCvpT2NwcupTmrDVgdwrn/uTat769CCl0omX7XvB1fo4crf9HKas7VKebVEKfGuGZdLGWAKoHVcipc8ffjNoF9pbYsGuTCyUq+aaOXfmdG5XdiMzYSot0yhOIPgiAoGeVT/FusuuVD343iytszoIpHbMttuGAcDs6pjFLOYme7V3ItFEq3LOoCo9kldDitJipTVtX3zJDqdfb/pRNpHAMvJJt3fdyuLnifoHeeUznbn1nvFVL5o+QqOmVXTF288qB9gUa/dSOtD1Kt3xtpXqozVskrKeXKDM27d9Ok7wqc9sDZ+BuTtCIbsbVv2ZbPKr8piR6fw3sGOb9NTsWwozb+UPm3DKH4gyAISkb5FP8WVQurrSZyVhNG9bHy7LW2dt4n4Lq6WHb7B6TKOU1vpcfVraVryUxaNTRGbXWUyjSpa2evHngpKVsuJFJZsGPifYKu6aRe62P2Hv0ewk1WDlg1XBlM76vTtBJyRUR0NNnauy5PtxwrAFSmKAiOlHdT2gST6IyS3Lm0GcIrNx+IRXNokqW5uTVKPO+3mbHqbEoLnmXjDwIiFH8QBEHJ2L+KvxMeBgX7yE9pm6GqnI1fukiJsw3Z+fH7VA/rnzslb5K+DZF69yuDroz78v/L50Ri/lpkU6/Mti75CADdC+SvPprUcd/EtGnHJSAbZP/ve8V6zWhG6mlxc951icZE35VcdtlB2LeeC964tBam6A0rarfaM946zqOLv29jh/exAFmeN3k2/jwlH3b4oACh+IMgCErG/lX8nVA2BW2q5hSntlihZ6l1ILv0orfRVkeSpwyr/NpVq16r5PGjbE9esN4/lbEDaEXlVVvMm1W+jjj7/7XWtnEMuZTSvNogO7xP+iaUEK7rOu07DFjlbZQ9ry5yyktimuIbXCRsg1IdV3jOmnzhadVg4jSyffW9tw7v6ZjfTI5NPvfzUPJBBwnFHwRBUDLiwR8EQVAy9q+pp9PkLK1NTVLXjl31tJFtCshK7eADcHRlpeXrpvzsGec01+allAPktqg+idzR8fT69atu8K1TGMBvMrLZaznNha9Oxa6k3K4yZxPRgVMikCml8bo1U7FbZZ65jc1eplatT5THZhW+X2+m4TFU3PdDm7vGBdibbAomA9zWhGtB6QjFHwRBUDJC8d+kYJK2phB52hhkt8+mmrtLnBKBgnG63P97SX0KXaspbIw2GisHx2istiUrW96oNJW+AMh1Uts+PbBPA3ETlyqazxPQZqfbFK1cvJze8KaoDz6bpaAoGoN3vzRVrXJWJJxiwbjUOqSfVhrUR8P1186mrV+RFA3gyg0cLEJsFgdEKP4gCIKSUUjxi8iDAP4lVmvu/rqq/oo7/osAPorVmrtXAPzvqvrK2rE6gG+tNX1VVT/UobEXo6DSaVJirLCyUiq36j+jvywbrQ/8MQm4qq0DjgDrgqiXU+K0yrBNnGbcOSkwqeHTFPBYR13yNbbRc3CTr+fLLpfsIplXYITOYXdLwKZI4GCpJtdJGpNQamej8AGj8oWCuRouVXKD0k1Uh2zyNTO+HLdcO8CM/aLVgWx4Tsv3myUUfkBs+OAXkSqATwH4MQATAJ4RkTOq+hw1+yaA06o6LyI/B+CfA/jptWMLqvruDo87CIIgaJMiiv8+AOdU9TwAiMgXADwEYP3Br6p/RO2fBvAznRzkpilaYIVo21OC1VwjO31zVv/+c07axnZ4k0YA1v7Pq4HGzIxpV6G0BxxUVDkwYtoZ758brswjr15Y5bsSjSYlAo2vcekysuDykuLKF7LKN/fogtTMPNG+g/j9GDqvkZPamFcUfI5fxWUF5fljQbDXKGLjPwbgAr2fWPssi48A+F163yciZ0XkaRH5G22MMQiCIOggRRR/qwxlLV0MRORnAJwG8F76+E5VvSgibwLwhyLyLVV9scW5jwJ4FAD6MOAPb46iRdQ7QJbfuPfCyEr05tMyM3Uqr+jPNx4/XLzF2fg5IZoZnys+Yoq5+Htib6WcvQFTwIWTyA1mf5/q0xlzd2zjp3YNlxyuSiuFOl/X7S2YfRaKufDlC7PUUJ6KL7z3U3AlGATbSRHFPwGAC74eB3DRNxKR9wH4RwA+pKrrNgVVvbj23/MAvgLgPa0uoqqPq+ppVT3djewqRUEQBMHWKKL4nwFwSkTuBvAagIcB/K/cQETeA+AzAB5U1cv0+RiAeVVdEpFxAD+E1Y3fnWObUzYXVWxs/zeXKTi+Jq8jUpg8hlxvHV4lLPiyhNmponk/wEQTs787XLI4o6LdtfhYjj994+q1NHZOc+zs6dyfOebiB5oKk2dgvK5yvl9erZk4DSDzN9Pk1cNEMfNgh9jwwa+qNRH5GICnsOrO+YSqPisijwE4q6pnAPwLAEMA/qOsLvdvum2+A8BnRKSB1dXFrzhvoCAIgmCHKeTHr6pPAnjSffYJev2+jPP+GMD3bGWAQRAEQWeJlA1ZFF1qFw3AIdh9MG9TMLdyE5kxTECTTzjGQ+JrZY271Xm8AUv9qwt8Mue4wDQDm3fYHOPu0VQtM0FQ9cx2fF1zfg6FN1zdnDWZd4JgnxApG4IgCEpGKP4i5KljUqK2GpcL6FnJcB/M6du0cxuTrFKzKwI7ODlcUzroYhuu5pyc1QrPhV/V8EqG1bV4RU3n5dU9znSzzFt1EE2ut8btc4vJ0TZDbOgGO0Qo/iAIgpJRPsXfjstcwdq8xg0wJ4DLqNyc5HBmL8Ap/qwCME3pB8jlsjGb0hl4lWxWEOrGxHNGqwFd8oqa7ov2HaSgy2vTyiArCVrBFNpNLrQFi55kud6GIg9uFULxB0EQlIzyKf5Oq7as/ryKzCjfmGdD5pWBD1rKwnuaVDPOayrlyEnanLdO1oqiKcEcl5fM8fhpZHgXNRW5oXuxdvdiSr7ZWydzSIX6C4JbhfhVB0EQlIzyKf7tpKD3T56iNHZtUrx+ZZBl/2/aM+CCLbwvUHc2/pz9CdNfjr0+S8kXXWXlJUEz13VzluWF0+Tv307B8h1M+BcEO0Uo/iAIgpIRD/4gCIKScWuYejqxBG+3j61mVHTnZKUByAtgsn1kV4kqSlPuf0p9YIKqclxRbR/Z9YcLjy93nlvXRPDf6ZZz3xd0692wbRDsMqH4gyAISsatofg7oa4KBgUVvnaHVyF5G5V5m5ZZFcK8Ss7b0JWsjekspb1Rf10F9UbRmgh8raKumJ3OfR8KP9hHhOIPgiAoGbeG4t8OCob3Fzrf9WHSN3hlzEnfTLvsBGZmfN6unaGGm6pY1TdRT3b9HGczL5iiOtOu3+4qaavfVRCUjFD8QRAEJSMUfxbbaAPOSx2Qaa/PC1rK81bhlAgcHFaw/mzTeQW9cGxSOnewYJqLtuhE4r0guMUJxR8EQVAyCj34ReRBEXlBRM6JyMdbHO8VkS+uHf9TETlJx35p7fMXROTHOzf0DiMV+08b6R9/nndeB9CGrv/LbVevr/9rGnvGfTRWauv/zP212o+gf+Y8blatmn/cnxlfu2Td035gv403KBUb/ipFpArgUwB+AsA9AB4RkXtcs48AmFLVtwD4VQD/bO3cewA8DOCdAB4E8Gtr/QVBEAS7RBE5ch+Ac6p6XlWXAXwBwEOuzUMAfmPt9W8B+FERkbXPv6CqS6r6EoBza/0FmyVPoftj9K/S07P+LxdSqFIR8y/ruqzqm1Yenb7n/cZ+Hntwy1PkL/QYgAv0fmLts5ZtVLUGYBrAoYLnAgBE5FEROSsiZ1eQncs9CIIg2BpFHvytksR4A3RWmyLnrn6o+riqnlbV093obdUkCIIg6ABF3DknAJyg98cBXMxoMyEiXQBGAUwWPHdvUDQtw3Ync+uwaaBwErSC7qa597GNZo1c99AgCDZFEcX/DIBTInK3iPRgdbP2jGtzBsDPrr3+SQB/qKq69vnDa14/dwM4BeB/dGboQRAEQTtsqPhVtSYiHwPwFFazcj2hqs+KyGMAzqrqGQD/DsC/F5FzWFX6D6+d+6yI/AcAzwGoAfh51X2u19pVtftpky9vVVM0NXGH73fLKZWDIFhHVoX53mJEDur98qO7PYzyslu1CYIgaJs/0N/6M1U9XaRtpGzYLvZzYY4yrGqCoMREWGEQBEHJCMW/XeyW+s1baYQpJggChOIPgiAoHaH49zKdKABT9NhuEauQINhxQvEHQRCUjFD8e5nt9JkvuhfQbv9FCZUfBDtOKP4gCIKSEQ/+IAiCkrF3TT03q0cFrdnq3Oy3TeAgCDpGKP4gCIKSsXcf/KE6gyAItoW9++APgiAItoV48AdBEJSMePAHQRCUjL3r1bPT7Oc0ykEQBJsgFH8QBEHJ2JMVuETkCoBXdviy4wCu7vA19zoxJ5aYj2ZiTiy7OR93qerhIg335IN/NxCRs0XLlpWFmBNLzEczMSeW/TIfYeoJgiAoGfHgD4IgKBnx4E88vtsD2IPEnFhiPpqJObHsi/kIG38QBEHJCMUfBEFQMkr34BeRB0XkBRE5JyIfb3H8F0XkORH5CxH57yJy126McyfZaE6o3U+KiIrInvda2ApF5kNEfmrtd/KsiPzmTo9xJynwN3OniPyRiHxz7e/mA7sxzp1CRJ4Qkcsi8u2M4yIi/2ptvv5CRL5vp8e4Iapamn8AqgBeBPAmAD0A/ieAe1ybHwEwsPb65wB8cbfHvdtzstZuGMBXATwN4PRuj3uXfyOnAHwTwNja+yO7Pe5dno/HAfzc2ut7ALy82+Pe5jn5YQDfB+DbGcc/AOB3AQiAHwDwp7s9Zv+vbIr/PgDnVPW8qi4D+AKAh7iBqv6Rqs6vvX0awPEdHuNOs+GcrPHLAP45gMWdHNwuUGQ+/h6AT6nqFACo6uUdHuNOUmQ+FMDI2utRABd3cHw7jqp+FcBkTpOHAHxOV3kawAERObozoytG2R78xwBcoPcTa59l8RGs/p/7VmbDORGR9wA4oar/dScHtksU+Y28FcBbReTrIvK0iDy4Y6PbeYrMxz8G8DMiMgHgSQB/f2eGtmfZ7HNmxylbkjZp8VlLtyYR+RkApwG8d1tHtPvkzomIVAD8KoAP79SAdpkiv5EurJp7HsDqivBrInKvql7f5rHtBkXm4xEAn1XV/0dEfhDAv1+bj7JmOiz8nNktyqb4JwCcoPfH0WJZKiLvA/CPAHxIVZd2aGy7xUZzMgzgXgBfEZGXsWqzPHMLb/AW+Y1MAPgdVV1R1ZcAvIDV/xHcihSZj48A+A8AoKp/AqAPqzlrykqh58xuUrYH/zMATonI3SLSA+BhAGe4wZpZ4zNYfejfyrbbm+TOiapOq+q4qp5U1ZNY3ff4kKqe3Z3hbjsb/kYA/DZWnQAgIuNYNf2c39FR7hxF5uNVAD8KACLyDqw++K/s6Cj3FmcA/J01754fADCtqq/v9qCYUpl6VLUmIh8D8BRWvRWeUNVnReQxAGdV9cMsWWUAAACqSURBVAyAfwFgCMB/FBEAeFVVP7Rrg95mCs5JaSg4H08BeL+IPAegDuAfqOq13Rv19lFwPv5PAP9WRH4BqyaND+uae8utiIh8HqtmvvG1fY1PAugGAFX9NFb3OT4A4ByAeQB/d3dGmk1E7gZBEJSMspl6giAISk88+IMgCEpGPPiDIAhKRjz4gyAISkY8+IMgCEpGPPiDIAhKRjz4gyAISkY8+IMgCErG/w9huOIGImK67AAAAABJRU5ErkJggg\u003d\u003d\n"
          },
          "metadata": {
            "needs_background": "light"
          },
          "output_type": "display_data"
        }
      ],
      "source": "hist \u003d plt.hist2d(xdata, ydata, bins\u003d(100, 100))\n",
      "metadata": {
        "pycharm": {
          "metadata": false,
          "name": "#%%\n",
          "is_executing": false
        }
      }
    },
    {
      "cell_type": "markdown",
      "source": "Finally, we perform the fit to this mock data using the ``LogLikelihood`` objective.",
      "metadata": {
        "pycharm": {
          "metadata": false,
          "name": "#%% md\n"
        }
      }
    },
    {
      "cell_type": "code",
      "execution_count": 17,
      "outputs": [
        {
          "name": "stdout",
          "text": [
            "\nParameter Value        Standard Deviation\nrho       6.026420e-01 2.013810e-03\nsig_x     1.100898e-01 2.461684e-04\nsig_y     2.303400e-01 5.150556e-04\nx0        5.901317e-01 3.481346e-04\ny0        8.014040e-01 7.283990e-04\nStatus message         b\u0027CONVERGENCE: REL_REDUCTION_OF_F_\u003c\u003d_FACTR*EPSMCH\u0027\nNumber of iterations   22\nObjective              \u003csymfit.core.objectives.LogLikelihood object at 0x000001DA8EB8BF28\u003e\nMinimizer              \u003csymfit.core.minimizers.LBFGSB object at 0x000001DA8EB6E7B8\u003e\n\nGoodness of fit qualifiers:\nlikelihood             inf\nlog_likelihood         106241.24669486462\nobjective_value        -106241.24669486462\n"
          ],
          "output_type": "stream"
        }
      ],
      "source": "fit \u003d Fit(pdf, x\u003dxdata, y\u003dydata, objective\u003dLogLikelihood)\nfit_result \u003d fit.execute()\nprint(fit_result)",
      "metadata": {
        "pycharm": {
          "metadata": false,
          "name": "#%%\n",
          "is_executing": false
        }
      }
    },
    {
      "cell_type": "markdown",
      "source": "We see that this result is in agreement with our data.",
      "metadata": {
        "pycharm": {
          "metadata": false,
          "name": "#%% md\n"
        }
      }
    }
  ],
  "metadata": {
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 2
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython2",
      "version": "2.7.6"
    },
    "kernelspec": {
      "name": "python3",
      "language": "python",
      "display_name": "Python 3"
    },
    "stem_cell": {
      "cell_type": "raw",
      "source": "",
      "metadata": {
        "pycharm": {
          "metadata": false
        }
      }
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}