diff --git a/doc/notebooks/correlated_data.ipynb b/doc/notebooks/correlated_data.ipynb new file mode 100644 index 00000000..51d7025c --- /dev/null +++ b/doc/notebooks/correlated_data.ipynb @@ -0,0 +1,269 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Fitting data with correlated uncertainties\n", + "\n", + "We sometimes want to combine results from different studies. If these results have independent uncertainties and can be expected to have the same mean, then the optimal combination is given by a weighted mean, where the weight is inversely proportional to the variance (uncertainty squared) of each individual input value. [This is a well-known result](https://en.wikipedia.org/wiki/Inverse-variance_weighting).\n", + "\n", + "If the uncertainties of the results are correlated, it is more complicated to compute an optimally weighted mean. Instead of deriving analytical formulas, we use a fit here to obtain the mixing weight, which is equivalent. It serves to demonstrate how fits to correlated data values can be carried out.\n", + "\n", + "We consider a toy example where two measurements should be combined which have strongly correlated systematic uncertainties." + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [], + "source": [ + "from iminuit import Minuit\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Here are the results we want to comine. The statistical uncertainties are assumed to be uncorrelated, the systematic uncertainties are assumed to be perfectly correlated (represented by thick bars in the plot)." + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGdCAYAAADAAnMpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAYI0lEQVR4nO3df5BWdd3/8dcuxIrhrmGprBBLYSWj0VpjQ9w1gAqs3tuUUzbaHSrlZOWUUpZMMxBNyX13SzkZZvYDbFIbdZByMIipjJFxMrrZ7oicInAgWag09gfWkuzefzju98vgIhex+9ldH4+Z88c51zl73hfM5T495yxb1dPT0xMAgEKqSw8AALy0iREAoCgxAgAUJUYAgKLECABQlBgBAIoSIwBAUWIEAChqZOkBjkZ3d3d2796dk046KVVVVaXHAQCOQk9PTzo6OlJfX5/q6r6vfwyJGNm9e3cmTJhQegwA4Bjs2rUr48eP7/P1IREjJ510UpLn3kxtbW3haQCAo9He3p4JEyb0fh/vy5CIkedvzdTW1ooRABhiXuwRCw+wAgBFiREAoCgxAgAUJUYAgKLECABQlBgBAIoSIwBAUWIEAChKjAAARVUcIxs2bEhzc3Pq6+tTVVWV1atXv+gxd911V6ZOnZoTTzwx48aNy/z58/PUU08dy7wAwDBTcYzs378/U6dOzfLly49q/40bN2bevHn54Ac/mN/+9re577778thjj+Xqq6+ueFgAYPip+HfTNDU1pamp6aj3f/TRR9PQ0JCPf/zjSZJJkyblwx/+cP7rv/6r0lMDAMNQvz8zMm3atOzatSsPPfRQenp6snfv3tx///256KKL+jymq6sr7e3thywAwPDU7zEyffr03HXXXXnf+96XUaNG5fTTT09dXd0Rb/MsXbo0dXV1vcuECRP6e0wA+tkzB55Nw41r0nDjmjxz4NnS4zCI9HuMbN26NZ/4xCeyaNGi/OpXv8ratWvzxBNP5JprrunzmIULF6atra132bVrV3+PCQAUUvEzI5VaunRppk+fnhtuuCFJ8sY3vjEvf/nL8/a3vz1f+MIXMm7cuMOOqampSU1NTX+PBgAMAv1+ZeSZZ55JdfWhpxkxYkSSpKenp79PDwAMchXHSGdnZ1paWtLS0pIk2bFjR1paWrJz584kz91imTdvXu/+zc3NWbVqVb7+9a9n+/bt2bhxYz7+8Y/nvPPOS319/fF5FwDAkFXxbZpNmzZl5syZvesLFixIklxxxRVZuXJlWltbe8MkSa688sp0dHTka1/7Wj75yU/m5JNPzqxZs/xoLwCQJKnqGQL3Strb21NXV5e2trbU1taWHgeAY/DMgWczZdG6JMnWz8/JiaP6/bFFCjva799+Nw0AUJQYAQCKEiMAQFFiBAAoSowAAEWJEQCgKDECABQlRgCAosQIAFCUGAEAihIjAEBRYgQAKEqMAABFiREAoCgxAgAUJUYAgKLECABQlBgBAIoSIwBAUWIEAChKjAAARYkRAKAoMQIAFCVGAICixAgAUJQYAQCKEiMAQFFiBAAoSowAAEWJEQCgKDECABQlRgDod091dmXKonW961MWrctTnV0FJ2IwESMAQFFiBAAoSowAAEWJEQCgKDECABQlRgCAosQIAFCUGAEAihIjAEBRYgQAKEqMAABFiREAoCgxAgAUJUYAgKLECABQlBgBAIoSIwBAUWIEAChKjAAARYkRAKAoMQIAFCVGAICixAgAUJQYAQCKEiMAQFFiBAAoSowAAEWJEQCgKDECABQlRgCAosQIAFCUGAEAiqo4RjZs2JDm5ubU19enqqoqq1evftFjurq68tnPfjYTJ05MTU1NGhoa8p3vfOdY5gUAhpmRlR6wf//+TJ06NfPnz88ll1xyVMdceuml2bt3b7797W9n8uTJaW1tTXd3d8XDAgDDT8Ux0tTUlKampqPef+3atfn5z3+e7du3Z+zYsUmShoaGSk8LAAxT/f7MyA9/+MO85S1vyZe+9KWcccYZed3rXpdPfepT+fvf/97nMV1dXWlvbz9kAQCGp4qvjFRq+/bteeSRR3LCCSfkgQceyF//+td89KMfzVNPPZUVK1a84DFLly7NkiVL+ns0AGAQ6PcrI93d3amqqspdd92V8847LxdddFG+/OUv58477+zz6sjChQvT1tbWu+zatau/xwQACun3KyPjxo3LGWeckbq6ut5tZ511Vnp6evKnP/0pZ5555mHH1NTUpKampr9HAwAGgX6/MjJ9+vTs3r07nZ2dvdt+//vfp7q6OuPHj+/v0wMAg1zFMdLZ2ZmWlpa0tLQkSXbs2JGWlpbs3LkzyXO3WObNm9e7/+WXX55TTjklV111VbZu3ZoNGzbkhhtuyPz58zN69Ojj8y4AgCGr4hjZtGlTGhsb09jYmCRZsGBBGhsbs2jRoiRJa2trb5gkyZgxY7J+/frs27cvb3nLW/L+978/zc3N+epXv3qc3gIAMJRV/MzIjBkz0tPT0+frK1euPGzbG97whqxfv77SUwEALwF+Nw0AUJQYAQCKEiMAQFFiBAAoSowAAEWJEQCgKDECABQlRgCAosQIAFCUGAEAihIjAEBRYgQAKEqMAABFiREAoCgxAgAUJUYAgKLECABQlBgBAIoSIwBAUWIEAChKjAAARYkRAKAoMQIAFCVGAICixAgAUJQYAaCIg909pUdgkBAjAPS7nz7+58O2/futj2TtltYC0zDYiBEA+tXaLa254f7/PWz7nzu68pHv/Y8gQYwA0H8OdvdkyYNbj7jPkge3umXzEidGAOg3j+14Oq1t/+jz9Z4krW3/yGM7nh64oRh0xAiDyjMHnk3DjWvScOOaPHPg2dLjAP+iP3f0HSLHsh/DkxgBoN+cetIJx3U/hicxAkC/OW/S2Iyr6zs0qpKMqzsh500aO3BDMeiIEQD6zYjqqixunnLEfRY3T8mI6qoBmojBSIwA0K/mnj0u//2eNx62/bTamnz9P87N3LPHFZiKwWRk6QEAGP5mveHUw7Y9eO2/5dRaz4rgyggAhbg1w/PECABQlBgBAIoSIwBAUWIEAChKjAAARYkRAKAoMQIAFCVGAICixAgAUJQYAQCKEiMAQFFiBAAoSowAAEWJEQCgKDECABQlRgCAosQIAFCUGAEAihIjAEBRYgQAKEqMAABFiREAoCgxAgAUJUYAgKLECABQlBgBAIoSIwBAURXHyIYNG9Lc3Jz6+vpUVVVl9erVR33sxo0bM3LkyLzpTW+q9LQAwDBVcYzs378/U6dOzfLlyys6bt++fZk3b17OP//8Sk8JAAxjIys9oKmpKU1NTRWf6Jprrsnll1+eESNGVHQ1BQAY3gbkmZEVK1Zk+/btWbx48VHt39XVlfb29kMWAGB46vcY+cMf/pAbb7wx3/ve9zJy5NFdiFm6dGnq6up6lwkTJvTzlABAKf0aIwcPHszll1+eJUuW5HWve91RH7dw4cK0tbX1Lrt27erHKQGAkip+ZqQSHR0d2bRpUzZv3pxrr702SdLd3Z2enp6MHDkyP/7xjzNr1qzDjqupqUlNTU1/jgYADBL9GiO1tbX5zW9+c8i22267LT/96U9z//33Z9KkSf15egBgCKg4Rjo7O7Nt27be9R07dqSlpSVjx47Nq1/96ixcuDBPPvlkvvvd76a6ujpnn332IcefeuqpOeGEEw7bDgC8NFUcI5s2bcrMmTN71xcsWJAkueKKK7Jy5cq0trZm586dx29CAGBYqzhGZsyYkZ6enj5fX7ly5RGP/9znPpfPfe5zlZ4WABim/G4aAKAoMQIAFCVGAICixAgAUJQYAQCKEiMAQFFiBAAoSowAAEWJEQCgKDECABQlRgCAosQIAFCUGAEAihIjAEBRYgQAKEqMAABFiREAoCgxAgAUJUYAgKLECABQlBgBAIoSIwBAUWIEAChKjAAARYkRAKAoMQIAFCVGAICixAgAUJQYAQCKEiMAQFFihEHjqc6uTFm0rnd9yqJ1eaqzq+BEwPFyypiabP38nN71rZ+fk1PG1BSciMFEjAAARYkRAKAoMQIAFCVGAICixAgAUJQYAQCKEiMAQFFiBAAoSowAAEWJEQCgKDECABQlRgCAosQIAFCUGAEAihIjAEBRYgQAKEqMAABFiREAoCgxAgAUJUYAgKLECABQlBgBAIoSIwBAUWIEAChKjAAARYkRAKAoMQIAFCVGAICixAgAUJQYAQCKEiMAQFFihEHtYHdP6REA6GcVx8iGDRvS3Nyc+vr6VFVVZfXq1Ufcf9WqVbnwwgvzqle9KrW1tZk2bVrWrVt3rPMyjP308T8ftu3fb30ka7e0FpgGgIFScYzs378/U6dOzfLly49q/w0bNuTCCy/MQw89lF/96leZOXNmmpubs3nz5oqHZfhau6U1N9z/v4dt/3NHVz7yvf8RJADD2MhKD2hqakpTU9NR73/LLbccsn7TTTflBz/4QR588ME0NjZWenqGoYPdPVny4NYj7rPkwa25cMrpGVFdNUBTATBQBvyZke7u7nR0dGTs2LF97tPV1ZX29vZDFoavx3Y8nda2f/T5ek+S1rZ/5LEdTw/cUAAMmIqvjPyrbr755nR2dubSSy/tc5+lS5dmyZIlAzgVJf25o+8QOZb9gMHpxFEj88R/Xlx6DAahAb0ycvfdd2fJkiW59957c+qpp/a538KFC9PW1ta77Nq1awCnZKCdetIJx3U/AIaWAbsy8v3vfz8f+tCHct999+WCCy444r41NTWpqakZoMko7bxJYzOu7oQ+b9VUJTm97oScN6nvW3sADF0DcmXknnvuyVVXXZV77rknF1/sEh2HGlFdlcXNU464z+LmKR5eBRimKo6Rzs7OtLS0pKWlJUmyY8eOtLS0ZOfOnUmeu8Uyb9683v3vvvvuzJs3L8uWLctb3/rW7NmzJ3v27ElbW9vxeQcMC3PPHpf/fs8bD9t+Wm1Nvv4f52bu2eMKTAXAQKg4RjZt2pTGxsbeH8tdsGBBGhsbs2jRoiRJa2trb5gkyR133JFnn302H/vYxzJu3Lje5ROf+MRxegsMF7PecPhzRA9e+29CBGCYq/iZkRkzZqSnp+9/onvlypWHrD/88MOVngJ6uTUDMPz53TQAQFFiBAAoSowAAEWJEQCgKDECABQlRgCAosQIAFCUGAEAihIjAEBRYgQAKEqMAABFiREAoCgxAgAUJUYAgKLECABQlBgBAIoSIwBAUWIEAChKjAAARYkRAKAoMQIAFCVGAICixAgAUJQYAQCKEiMAQFFiBAAoSowAAEWJEQCgKDECABQlRgCAosQIAFCUGGHQOGVMTbZ+fk7v+tbPz8kpY2oKTgTAQBAjAEBRYgQAKEqMAABFiREAoCgxAgAUJUYAgKLECABQlBgBAIoSIwBAUWIEAChKjAAARYkRAKAoMQIAFCVGAICixAgAUJQYAQCKEiMAQFFiBAAoSowAAEWJEQCgKDECABQlRgCAosQIAFCUGAEAihIjAEBRYgQAKEqMAABFiREAoCgxAgAUJUYAgKLECABQlBgBAIqqOEY2bNiQ5ubm1NfXp6qqKqtXr37RYx5++OGce+65qampyeTJk7Ny5cpjGBUAGI4qjpH9+/dn6tSpWb58+VHtv2PHjlx88cWZOXNmWlpact111+VDH/pQ1q1bV/GwAMDwM7LSA5qamtLU1HTU+99+++2ZNGlSli1bliQ566yz8sgjj+QrX/lK5syZU+npAYBhpt+fGXn00UdzwQUXHLJtzpw5efTRR/s8pqurK+3t7YcsAMDw1O8xsmfPnpx22mmHbDvttNPS3t6ev//97y94zNKlS1NXV9e7TJgwob/HBAAKGZQ/TbNw4cK0tbX1Lrt27So9EgDQTyp+ZqRSp59+evbu3XvItr1796a2tjajR49+wWNqampSU1PT36MBAINAv18ZmTZtWn7yk58csm39+vWZNm1af58aABgCKo6Rzs7OtLS0pKWlJclzP7rb0tKSnTt3JnnuFsu8efN697/mmmuyffv2fPrTn87jjz+e2267Lffee2+uv/764/MOAIAhreIY2bRpUxobG9PY2JgkWbBgQRobG7No0aIkSWtra2+YJMmkSZOyZs2arF+/PlOnTs2yZcvyrW99y4/1AgBJjuGZkRkzZqSnp6fP11/oX1edMWNGNm/eXOmpAICXgEH50zQAwEuHGAEAihIjAEBRYgQAKEqMAABFiREAoCgxAgAUJUYAgKLECABQlBgBAIoSIwBAUWIEAChKjAAARYkRAKAoMQIAFCVGAICixAgAUJQYAQCKEiMAQFFiBAAoSowAAEWJEQCgKDECABQlRgCAosQIAFCUGAEAihIjAEBRYgQAKEqMAABFiREAoKiRpQeA/9+Jo0bmif+8uPQYAAwgV0YAgKLECABQlBgBAIoSIwBAUWIEAChKjAAARYkRAKAoMQIAFCVGAICixAgAUJQYAQCKEiMAQFFiBAAoSowAAEWJEQCgqJGlBzgaPT09SZL29vbCkwAAR+v579vPfx/vy5CIkY6OjiTJhAkTCk8CAFSqo6MjdXV1fb5e1fNiuTIIdHd3Z/fu3TnppJNSVVVVehz6WXt7eyZMmJBdu3altra29DjAceTz/dLS09OTjo6O1NfXp7q67ydDhsSVkerq6owfP770GAyw2tpa/7GCYcrn+6XjSFdEnucBVgCgKDECABQlRhh0ampqsnjx4tTU1JQeBTjOfL55IUPiAVYAYPhyZQQAKEqMAABFiREAoCgxwpDR0NCQW265pfQYQD/w+X5pEyMMWVVVVVm9evWL7vfFL34xb3vb23LiiSfm5JNP7ve5gH/d0Xy+n3jiiXzwgx/MpEmTMnr06Lz2ta/N4sWLc+DAgYEZkuNmSPwLrAxuBw4cyKhRo0qP0acDBw7kve99b6ZNm5Zvf/vbpceBIWUwf74ff/zxdHd35xvf+EYmT56cLVu25Oqrr87+/ftz8803lx6PCrgyQsVmzJiRa6+9Ntddd11e+cpXZs6cOUmSLVu2pKmpKWPGjMlpp52WD3zgA/nrX//ae9z999+fc845J6NHj84pp5ySCy64IPv37+/9mtddd90h53nXu96VK6+88gVnaGhoSJK8+93vTlVVVe/6C1myZEmuv/76nHPOOcf8nuGlYih9vufOnZsVK1Zk9uzZec1rXpN3vvOd+dSnPpVVq1b9S38GDDwxwjG58847M2rUqGzcuDG333579u3bl1mzZqWxsTGbNm3K2rVrs3fv3lx66aVJktbW1lx22WWZP39+fve73+Xhhx/OJZdc8qK/Vrovv/zlL5MkK1asSGtra+868K8byp/vtra2jB079pjOSzlu03BMzjzzzHzpS1/qXf/CF76QxsbG3HTTTb3bvvOd72TChAn5/e9/n87Ozjz77LO55JJLMnHixCT5l65UvOpVr0qSnHzyyTn99NOP+esAhxuqn+9t27bl1ltvdYtmCBIjHJM3v/nNh6z/+te/zs9+9rOMGTPmsH3/+Mc/Zvbs2Tn//PNzzjnnZM6cOZk9e3be85735BWveMVAjQwcpaH4+X7yySczd+7cvPe9783VV189YOfl+HCbhmPy8pe//JD1zs7ONDc3p6Wl5ZDlD3/4Q97xjndkxIgRWb9+fX70ox9lypQpufXWW/P6178+O3bsSJJUV1cfdkn3n//854C9H+D/GWqf7927d2fmzJl529veljvuuOO4fV0GjhjhuDj33HPz29/+Ng0NDZk8efIhy/P/Yauqqsr06dOzZMmSbN68OaNGjcoDDzyQ5LnLsq2trb1f7+DBg9myZcsRz/myl70sBw8e7L83BSQZ3J/vJ598MjNmzMib3/zmrFixItXVvq0NRf7WOC4+9rGP5emnn85ll12WX/7yl/njH/+YdevW5aqrrsrBgwfzi1/8IjfddFM2bdqUnTt3ZtWqVfnLX/6Ss846K0kya9asrFmzJmvWrMnjjz+ej3zkI9m3b98Rz9nQ0JCf/OQn2bNnT/72t7/1ud/OnTvT0tKSnTt35uDBg73/V9fZ2Xk8/whg2Bqsn+/nQ+TVr351br755vzlL3/Jnj17smfPnuP9R0A/88wIx0V9fX02btyYz3zmM5k9e3a6uroyceLEzJ07N9XV1amtrc2GDRtyyy23pL29PRMnTsyyZcvS1NSUJJk/f35+/etfZ968eRk5cmSuv/76zJw584jnXLZsWRYsWJBvfvObOeOMM/LEE0+84H6LFi3KnXfe2bve2NiYJPnZz36WGTNmHJf3D8PZYP18r1+/Ptu2bcu2bdsyfvz4Q17zC+mHlqoef2MAQEFu0wAARYkRAKAoMQIAFCVGAICixAgAUJQYAQCKEiMAQFFiBAAoSowAAEWJEQCgKDECABQlRgCAov4PBCmUSEx9yeUAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "value = [1.2, 1.5]\n", + "error_sta = [0.3, 0.3]\n", + "error_sys = [0.1, 0.2]\n", + "correlation_sys = 1.0\n", + "\n", + "plt.errorbar((\"result 1\", \"result 2\"), value, error_sta, fmt=\"o\")\n", + "plt.errorbar((\"result 1\", \"result 2\"), value, error_sys, lw=3, fmt=\"none\")\n", + "plt.xlim(-0.5, 1.5);" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "For the combination, we formulate this as a fitting problem. Our objective function is derived from the log-probability of a multivariate normal distribution (in the derivation we dropped constants and scaled the result). We predict the constant mean of this distribution, which is matched to the two observed values while taking their covariance into account. The covariance matrix is the way to include correlated observations.\n", + "\n", + "The simpler special case for uncorrelated observations is handled in `iminuit.cost.LeastSquares`, but for the general case there is no ready-made cost function yet, so we write it here." + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
Migrad
FCN = 0.4737 (χ²/ndof = 0.5) Nfcn = 13
EDM = 1.46e-07 (Goal: 0.0002)
Valid Minimum Below EDM threshold (goal x 10)
No parameters at limit Below call limit
Hesse ok Covariance accurate
\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
Name Value Hesse Error Minos Error- Minos Error+ Limit- Limit+ Fixed
0 z 0.4 0.7 0 1
\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
z
z 0.737
" + ], + "text/plain": [ + "┌─────────────────────────────────────────────────────────────────────────┐\n", + "│ Migrad │\n", + "├──────────────────────────────────┬──────────────────────────────────────┤\n", + "│ FCN = 0.4737 (χ²/ndof = 0.5) │ Nfcn = 13 │\n", + "│ EDM = 1.46e-07 (Goal: 0.0002) │ │\n", + "├──────────────────────────────────┼──────────────────────────────────────┤\n", + "│ Valid Minimum │ Below EDM threshold (goal x 10) │\n", + "├──────────────────────────────────┼──────────────────────────────────────┤\n", + "│ No parameters at limit │ Below call limit │\n", + "├──────────────────────────────────┼──────────────────────────────────────┤\n", + "│ Hesse ok │ Covariance accurate │\n", + "└──────────────────────────────────┴──────────────────────────────────────┘\n", + "┌───┬──────┬───────────┬───────────┬────────────┬────────────┬─────────┬─────────┬───────┐\n", + "│ │ Name │ Value │ Hesse Err │ Minos Err- │ Minos Err+ │ Limit- │ Limit+ │ Fixed │\n", + "├───┼──────┼───────────┼───────────┼────────────┼────────────┼─────────┼─────────┼───────┤\n", + "│ 0 │ z │ 0.4 │ 0.7 │ │ │ 0 │ 1 │ │\n", + "└───┴──────┴───────────┴───────────┴────────────┴────────────┴─────────┴─────────┴───────┘\n", + "┌───┬───────┐\n", + "│ │ z │\n", + "├───┼───────┤\n", + "│ z │ 0.737 │\n", + "└───┴───────┘" + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# construct covariance matrices\n", + "cov_sta = np.diag(np.square(error_sta))\n", + "cov_sys = np.diag(np.square(error_sys))\n", + "cov_sys[0, 1] = error_sys[0] * error_sys[1] * correlation_sys\n", + "cov_sys[1, 0] = cov_sys[0, 1]\n", + "\n", + "# total covariance is sum of individual contributions\n", + "cov = cov_sta + cov_sys\n", + "inv_cov = np.linalg.inv(cov)\n", + "\n", + "def model(x, z):\n", + " \"\"\"Return combined value based on input x and mixing parameter z.\"\"\"\n", + " return z * x[1] + (1 - z) * x[0]\n", + "\n", + "def cost(z):\n", + " \"\"\"Chi-square distributed cost function.\"\"\"\n", + " xp = model(value, z)\n", + " delta = value - xp\n", + " return np.einsum(\"i,j,ij\", delta, delta, inv_cov)\n", + "\n", + "# with this extra information, iminuit will also print the chi2/ndof gof statistic\n", + "cost.errordef = Minuit.LEAST_SQUARES\n", + "cost.ndata = len(value)\n", + "\n", + "m = Minuit(cost, 0.5)\n", + "m.limits[\"z\"] = (0, 1)\n", + "m.migrad()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Our fit gives more weight to the more accurate measurement, as expected. In order to find the statistical and systematic uncertainty of the combined result, we do error propagation. We compute the trivial Jacobian for our model analytically." + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "total = 1.33 +/- 0.21(sta) + 0.16(sys)\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGdCAYAAADAAnMpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAeDklEQVR4nO3dfVBVBf7H8c8F4iIiKJYIScJmD7Iae8u1Qbf5oZVIxW7b0062YVlutrlZlhXThGGbaGk1tfZgtVCbtlOOujk2klsak1sZLrcytFJxoLhoWfJkYcD5/eF058cPMa5x75eH92vmznTOPeee7/VmvDvnAC7HcRwBAAAYCbMeAAAA9G/ECAAAMEWMAAAAU8QIAAAwRYwAAABTxAgAADBFjAAAAFPECAAAMBVhPUBXtLW1qaamRoMGDZLL5bIeBwAAdIHjOGpoaFBSUpLCwjo//9ErYqSmpkbJycnWYwAAgONQXV2tESNGdPp8r4iRQYMGSTryZmJjY42nAQAAXVFfX6/k5GT/1/HO9IoY+fHSTGxsLDECAEAv81O3WHADKwAAMEWMAAAAU8QIAAAwRYwAAABTxAgAADBFjAAAAFPECAAAMEWMAAAAU8QIAAAwFXCMlJaWKicnR0lJSXK5XFq7du1P7rNixQqlp6crOjpaiYmJmjFjhg4cOHA88wIAgD4m4BhpampSenq6li1b1qXtt2zZotzcXN1www365JNP9Oqrr2rr1q2aOXNmwMMCAIC+J+DfTZOdna3s7Owub//uu+8qJSVFt956qyQpNTVVN910kxYvXhzooQEAQB8U9HtGMjIyVF1drddff12O42jfvn1atWqVLrrook73aW5uVn19fbsHAADom4IeIxMnTtSKFSv0hz/8QZGRkRo+fLji4uKOeZmnsLBQcXFx/kdycnKwxwRg4NDhFqXcs14p96zXocMt1uMAMBL0GKmoqNCcOXOUn5+vbdu2acOGDdq7d69mzZrV6T55eXmqq6vzP6qrq4M9JgAAMBLwPSOBKiws1MSJEzVv3jxJ0llnnaWBAwfqvPPO01//+lclJiZ22Mftdsvtdgd7NAAA0AME/czIoUOHFBbW/jDh4eGSJMdxgn14AADQwwUcI42NjfJ6vfJ6vZKkyspKeb1eVVVVSTpyiSU3N9e/fU5OjlavXq2nnnpKe/bs0ZYtW3Trrbdq/PjxSkpK6p53AQAAeq2AL9OUlZVp0qRJ/uW5c+dKkqZPn67i4mL5fD5/mEjSddddp4aGBv3tb3/THXfcocGDB2vy5Ml8ay8AAJAkuZxecK2kvr5ecXFxqqurU2xsrPU4ALrJocMtSssvkSRVLMhSdGTQb2MDEEJd/frN76YBAACmiBEAAGCKGAEAAKaIEQAAYIoYAQAApogRAABgihgBAACmiBEAAGCKGAEAAKaIEQAAYIoYAQAApogRAABgihgBAACmiBEAAGCKGAEAAKaIEQAAYIoYAQAApogRAABgihgBAACmiBEAAGCKGAEAAKaIEQAAYIoYAQAApogRAABgihgBAACmiBEAAGCKGAEAAKaIEQAAYIoYAQAApogRAABgihgBAACmiBEAJg40Nistv8S/nJZfogONzYYTAbBCjAAAAFPECAAAMEWMAAAAU8QIAAAwRYwAAABTxAgAADBFjAAAAFPECAAAMEWMAAAAU8QIAAAwRYwAAABTxAgAADBFjAAAAFPECAAAMEWMAAAAU8QIAAAwRYwAAABTxAgAADBFjAAAAFPECAAAMEWMAAAAU8QIAAAwRYwAAABTxAgAADBFjAAAAFPECAAAMEWMAAAAU8QIAAAwRYwAAABTxAgAADBFjAAAAFPECAAAMBVwjJSWlionJ0dJSUlyuVxau3btT+7T3Nyse++9VyNHjpTb7VZKSor+/ve/H8+8AACgj4kIdIempialp6drxowZuuyyy7q0z1VXXaV9+/bp+eef16hRo+Tz+dTW1hbwsAAAoO8JOEays7OVnZ3d5e03bNigt99+W3v27FF8fLwkKSUlJdDDAgCAPiro94y89tprGjdunB566CGdfPLJOv3003XnnXfqu+++63Sf5uZm1dfXt3sAAIC+KeAzI4Has2eP3nnnHUVFRWnNmjX6+uuv9ec//1kHDhxQUVHRUfcpLCxUQUFBsEcDAAA9QNDPjLS1tcnlcmnFihUaP368LrroIj3yyCN64YUXOj07kpeXp7q6Ov+juro62GMCAAAjQT8zkpiYqJNPPllxcXH+daNHj5bjOPriiy902mmnddjH7XbL7XYHezQAANADBP3MyMSJE1VTU6PGxkb/us8++0xhYWEaMWJEsA8PAAB6uIBjpLGxUV6vV16vV5JUWVkpr9erqqoqSUcuseTm5vq3nzZtmoYOHarrr79eFRUVKi0t1bx58zRjxgwNGDCge94FAADotQKOkbKyMnk8Hnk8HknS3Llz5fF4lJ+fL0ny+Xz+MJGkmJgYbdy4UQcPHtS4ceN0zTXXKCcnR48//ng3vQUAANCbBXzPSGZmphzH6fT54uLiDuvOPPNMbdy4MdBDAQCAfoDfTQMAAEwRIwAAwBQxAgAATBEjAADAFDECAABMESMAAMAUMQIAAEwRIwAAwBQxAgAATBEjAADAFDECAABMESMAAMAUMQIAAEwRIwAAwBQxAgAATBEjAICQOHS4RSn3rFfKPet16HCL9TjoQYgRAABgihgBAACmiBEAAGCKGAEAAKaIEQAAYIoYAQAApogRAABgihgBAACmiBEAAGCKGAHQY7S2OdYjADBAjAAw8dbO/R3WXfLEO9qw3WcwDQBLxAiAkNuw3ad5qz7qsH5/Q7Nufum/BAnQzxAjAEKqtc1RwbqKY25TsK6CSzZAP0KMAAiprZXfyFf3fafPO5J8dd9ra+U3oRsKgCliBD0Kv2K879vf0HmIHM92AHo/YgRASA0bFNWt2wHo/YgRACE1PjVeiXGdh4ZLUmJclManxoduKACmiBEAIRUe5tL8nLRjbjM/J03hYa4QTQTAGjECIOSmjknUw1ec1WF9QqxbT/3xbE0dk2gwFQArEdYDAOifJp85rMO6dbN/o2Gx3CsC9DecGQHQY3Bppu860NistPwS/3JafokONDYbToSehBgBAACmiBEAAGCKGAEAAKaIEQAAYIoYAQAApogRAABgihgBAACmiBEAAGCKGAEAAKaIEQAAYIoYAQAApogRAABgihgBAACmiBEAAGCKGAEAAKaIEQAAYIoYAQAApogRAABgihgBAACmiBEAAGCKGAEAAKaIEQAAYIoYAQAApogRAABgihgBAACmiBEAAGCKGAEAAKYCjpHS0lLl5OQoKSlJLpdLa9eu7fK+W7ZsUUREhH71q18FelgAANBHBRwjTU1NSk9P17JlywLa7+DBg8rNzdX5558f6CEBAEAfFhHoDtnZ2crOzg74QLNmzdK0adMUHh4e0NkUAADQt4XknpGioiLt2bNH8+fP79L2zc3Nqq+vb/cAAAB9U9Bj5PPPP9c999yjl156SRERXTsRU1hYqLi4OP8jOTk5yFMCAAArQY2R1tZWTZs2TQUFBTr99NO7vF9eXp7q6ur8j+rq6iBOCQAALAV8z0ggGhoaVFZWpvLycs2ePVuS1NbWJsdxFBERoTfeeEOTJ0/usJ/b7Zbb7Q7maAAAoIcIaozExsbq448/brfuySef1FtvvaVVq1YpNTU1mIcHAAC9QMAx0tjYqF27dvmXKysr5fV6FR8fr1NOOUV5eXn68ssv9eKLLyosLExjxoxpt/+wYcMUFRXVYT0AAOifAo6RsrIyTZo0yb88d+5cSdL06dNVXFwsn8+nqqqq7psQAAD0aQHHSGZmphzH6fT54uLiY+5///336/777w/0sAAAoI/id9MAAABTxAgAADBFjAAAAFPECAAAMEWMAAAAU8QIAAAwRYwAAABTxAgAADBFjAAAAFPECAAAMEWMAAAAU8QIAAAwRYwAAABTxAgAADBFjAAAAFPECADARGubYz0CeghiBAAQdG/t3N9h3SVPvKMN230G06CnIUYAAEG1YbtP81Z91GH9/oZm3fzSfwkSECMAgOBpbXNUsK7imNsUrKvgkk0/R4wAAIJma+U38tV93+nzjiRf3ffaWvlN6IZCj0OMAACCZn9D5yFyPNuhbyJGAABBM2xQVLduh76JGAEABM341HglxnUeGi5JiXFRGp8aH7qh0OMQIwCAoAkPc2l+Ttoxt5mfk6bwMFeIJkJPRIwAAIJq6phEPXzFWR3WJ8S69dQfz9bUMYkGU6EnibAeAADQ900+c1iHdetm/0bDYrlXBJwZAQAY4dIMfkSMAAAAU8QIAAAwRYwAAABTxAh6jAONzUrLL/Evp+WX6EBjs+FECKahMW5VLMjyL1csyNLQGLfhRACsECMAAMAUMQIAAEwRIwAAwBQxAgAATBEjAADAFDECAABMESMAAMAUMQIAAEwRIwAAwBQxAgAATBEjAADAFDECAABMESMAAMAUMQIAAEwRIwAAwBQxAgAATBEjAADAFDECAABMESMAAMAUMQIAAEwRIwAAwBQxAgAATBEjAADAFDECAABMESMAAMAUMQIAAEwRIwAAwBQxAgAATBEjAADAFDECAABMESMAAMAUMYIerbXNsR4BABBkAcdIaWmpcnJylJSUJJfLpbVr1x5z+9WrV+vCCy/USSedpNjYWGVkZKikpOR450Uf9tbO/R3WXfLEO9qw3WcwDQAgVAKOkaamJqWnp2vZsmVd2r60tFQXXnihXn/9dW3btk2TJk1STk6OysvLAx4WfdeG7T7NW/VRh/X7G5p180v/JUgAoA+LCHSH7OxsZWdnd3n7xx57rN3ywoUL9a9//Uvr1q2Tx+MJ9PDog1rbHBWsqzjmNgXrKnRh2nCFh7lCNBUAIFRCfs9IW1ubGhoaFB8f3+k2zc3Nqq+vb/dA37W18hv56r7v9HlHkq/ue22t/CZ0QwEAQibgMyM/15IlS9TY2Kirrrqq020KCwtVUFAQwqlgaX9D5yFyPNuh94iOjNDeRRdbjwHAWEjPjKxcuVIFBQV65ZVXNGzYsE63y8vLU11dnf9RXV0dwikRasMGRXXrdgCA3iVkZ0b++c9/6sYbb9Srr76qCy644Jjbut1uud3uEE0Ga+NT45UYF9XppRqXpOFxURqf2vmlPQBA7xWSMyMvv/yyrr/+er388su6+GJOyaK98DCX5uekHXOb+Tlp3LwKAH1UwDHS2Ngor9crr9crSaqsrJTX61VVVZWkI5dYcnNz/duvXLlSubm5Wrp0qc4991zV1taqtrZWdXV13fMO0CdMHZOoh684q8P6hFi3nvrj2Zo6JtFgKgBAKAQcI2VlZfJ4PP5vy507d648Ho/y8/MlST6fzx8mkrR8+XK1tLTolltuUWJiov8xZ86cbnoL6Csmn9nxPqJ1s39DiABAHxfwPSOZmZlynM5/RHdxcXG75c2bNwd6CMCPSzMA0Pfxu2kAAIApYgQAAJgiRgAAgCliBAAAmCJGAACAKWIEABB0Q2PcqliQ5V+uWJCloTH8pG0cQYwAAABTxAgAADBFjAAAAFPECAAAMEWMAAAAU8QIAAAwRYwAAABTxAgAADBFjAAAAFPECAAAMEWMAAAAU8QIAAAwRYwAAABTxAgAADBFjAAAAFPECAAAMEWMAAAAU8QIAAAwRYwAAABTxAgAADAVYT0AAKB/iI6M0N5FF1uPgR6IMyMAAMAUMQIAAEwRIwAAwBQxAgAATBEj6DGGxrhVsSDLv1yxIEtDY9yGEwEAQoEYAQAApogRAABgihgBAACmiBEAAGCKGAEAAKaIEQAAYIoYAQAApogRAABgihgBAACmiBEAAGCKGAEAAKaIEQAAYIoYAQAApogRAABgihgBAACmiBEAAGCKGAEAAKaIEQAAYIoYAQAApogRAABgihgBAACmiBEAAGCKGAEAAKaIEQAAYIoYAQAApogRAABgihgBAACmiBEAAGCKGAEAAKaIEQAAYIoYAQAApogRAABgKuAYKS0tVU5OjpKSkuRyubR27dqf3Gfz5s06++yz5Xa7NWrUKBUXFx/HqAAAoC8KOEaampqUnp6uZcuWdWn7yspKXXzxxZo0aZK8Xq9uu+023XjjjSopKQl4WAAA0PdEBLpDdna2srOzu7z9008/rdTUVC1dulSSNHr0aL3zzjt69NFHlZWVFejhAQBAHxP0e0beffddXXDBBe3WZWVl6d133+10n+bmZtXX17d7AACAvinoMVJbW6uEhIR26xISElRfX6/vvvvuqPsUFhYqLi7O/0hOTg72mAAAwEiP/G6avLw81dXV+R/V1dXWIwEAgCAJ+J6RQA0fPlz79u1rt27fvn2KjY3VgAEDjrqP2+2W2+0O9mgAAKAHCPqZkYyMDL355pvt1m3cuFEZGRnBPjQAAOgFAo6RxsZGeb1eeb1eSUe+ddfr9aqqqkrSkUssubm5/u1nzZqlPXv26K677tLOnTv15JNP6pVXXtHtt9/ePe8AAAD0agHHSFlZmTwejzwejyRp7ty58ng8ys/PlyT5fD5/mEhSamqq1q9fr40bNyo9PV1Lly7Vc889x7f1AgAAScdxz0hmZqYcx+n0+aP9dNXMzEyVl5cHeigAANAP9MjvpgEAAP0HMQIAAEwRIwAAwBQxAgAATBEjAADAFDECAABMESMAAMAUMQIAAEwRIwAAwBQxAgAATBEjAADAFDECAABMESMAAMAUMQIAAEwRIwAAwBQxAgAATBEjAADAFDECAABMESMAAMAUMQIAAEwRIwAAwBQxAgAATBEjAADAFDECAABMESMAAMAUMQIAAEwRIwAAwBQxAgAATBEjAADAFDECAABMRVgPAPxf0ZER2rvoYusxAAAhxJkRAABgihgBAACmiBEAAGCKGAEAAKaIEQAAYIoYAQAApogRAABgihgBAACmiBEAAGCKGAEAAKaIEQAAYIoYAQAApogRAABgihgBAACmiBEAAGAqwnqArnAcR5JUX19vPAkAAOiqH79u//h1vDO9IkYaGhokScnJycaTAACAQDU0NCguLq7T513OT+VKD9DW1qaamhoNGjRILpfLehwEWX19vZKTk1VdXa3Y2FjrcRBkfN79C593/+I4jhoaGpSUlKSwsM7vDOkVZ0bCwsI0YsQI6zEQYrGxsfzHqh/h8+5f+Lz7j2OdEfkRN7ACAABTxAgAADBFjKDHcbvdmj9/vtxut/UoCAE+7/6FzxtH0ytuYAUAAH0XZ0YAAIApYgQAAJgiRgAAgCliBL1GSkqKHnvsMesxECJ83n3D3r175XK55PV6O91m8+bNcrlcOnjwoPkssEGMoNdyuVxau3btT2734IMPasKECYqOjtbgwYODPheCoyuf9969e3XDDTcoNTVVAwYM0Kmnnqr58+fr8OHDoRkSx2XChAny+Xxd+uFY6Jt6xU9gRc92+PBhRUZGWo/RqcOHD+vKK69URkaGnn/+eetxer2e/Hnv3LlTbW1teuaZZzRq1Cht375dM2fOVFNTk5YsWWI9HjoRGRmp4cOHW48BQ5wZQcAyMzM1e/Zs3XbbbTrxxBOVlZUlSdq+fbuys7MVExOjhIQEXXvttfr666/9+61atUpjx47VgAEDNHToUF1wwQVqamryv+Ztt93W7jiXXnqprrvuuqPOkJKSIkn6/e9/L5fL5V8+moKCAt1+++0aO3bscb/n/qw3fd5Tp05VUVGRpkyZol/84hf67W9/qzvvvFOrV6/+WX8GfUlbW5seeughjRo1Sm63W6eccooefPBBSdLHH3+syZMn+z+zP/3pT2psbPTve9111+nSSy/VwoULlZCQoMGDB2vBggVqaWnRvHnzFB8frxEjRqioqKjDcXfu3KkJEyYoKipKY8aM0dtvv+1/7v9fpikuLtbgwYNVUlKi0aNHKyYmRlOnTpXP52v3ms8995xGjx6tqKgonXnmmXryySfbPb9161Z5PB5FRUVp3LhxKi8v764/RnQzYgTH5YUXXlBkZKS2bNmip59+WgcPHtTkyZPl8XhUVlamDRs2aN++fbrqqqskST6fT1dffbVmzJihHTt2aPPmzbrssst+8tdKd+aDDz6QJBUVFcnn8/mXERy9+fOuq6tTfHz8cR23L8rLy9OiRYt03333qaKiQitXrlRCQoKampqUlZWlIUOG6IMPPtCrr76qf//735o9e3a7/d966y3V1NSotLRUjzzyiObPn69LLrlEQ4YM0fvvv69Zs2bppptu0hdffNFuv3nz5umOO+5QeXm5MjIylJOTowMHDnQ656FDh7RkyRL94x//UGlpqaqqqnTnnXf6n1+xYoXy8/P14IMPaseOHVq4cKHuu+8+vfDCC5KkxsZGXXLJJUpLS9O2bdt0//33t9sfPYwDBOh//ud/HI/H027dAw884EyZMqXduurqakeS8+mnnzrbtm1zJDl79+7t9DXnzJnTbt3vfvc7Z/r06f7lkSNHOo8++qh/WZKzZs2aLs9dVFTkxMXFdXl7HNFbP2/HcZzPP//ciY2NdZYvXx7Qfn1VfX2943a7nWeffbbDc8uXL3eGDBniNDY2+tetX7/eCQsLc2prax3HcZzp06c7I0eOdFpbW/3bnHHGGc55553nX25paXEGDhzovPzyy47jOE5lZaUjyVm0aJF/mx9++MEZMWKEs3jxYsdxHGfTpk2OJOfbb791HOfI31VJzq5du/z7LFu2zElISPAvn3rqqc7KlSvbvYcHHnjAycjIcBzHcZ555hln6NChznfffed//qmnnnIkOeXl5V37A0PIcM8Ijss555zTbvnDDz/Upk2bFBMT02Hb3bt3a8qUKTr//PM1duxYZWVlacqUKbriiis0ZMiQUI2Mn6E3ft5ffvmlpk6dqiuvvFIzZ84M2XF7sh07dqi5uVnnn3/+UZ9LT0/XwIED/esmTpyotrY2ffrpp0pISJAk/fKXv2z3q+ATEhI0ZswY/3J4eLiGDh2q/fv3t3v9jIwM/z9HRERo3Lhx2rFjR6ezRkdH69RTT/UvJyYm+l+zqalJu3fv1g033NDus21pafHfBLtjxw6dddZZioqKOuoM6FmIERyX//sfLOnIKdGcnBwtXry4w7aJiYkKDw/Xxo0b9Z///EdvvPGGnnjiCd177716//33lZqaqrCwsA6n8H/44Yegvgd0XW/7vGtqajRp0iRNmDBBy5cv77bX7e0GDBjws1/jhBNOaLfscrmOuq6tra3bj/PjvzM/3sfy7LPP6txzz223XXh4+M86Lmxwzwi6xdlnn61PPvlEKSkpGjVqVLvHj1/IXC6XJk6cqIKCApWXlysyMlJr1qyRJJ100kntbk5rbW3V9u3bj3nME044Qa2trcF7U+hUT/68v/zyS2VmZuqcc85RUVFRu/+L7+9OO+00DRgwQG+++WaH50aPHq0PP/zQf5OxJG3ZskVhYWE644wzfvax33vvPf8/t7S0aNu2bRo9evRxvVZCQoKSkpK0Z8+eDv/+paamSjryfj766CN9//33R50BPQt/S9EtbrnlFn3zzTe6+uqr9cEHH2j37t0qKSnR9ddfr9bWVr3//vtauHChysrKVFVVpdWrV+urr77y/8do8uTJWr9+vdavX6+dO3fq5ptv/skfgJSSkqI333xTtbW1+vbbbzvdrqqqSl6vV1VVVWptbZXX65XX6233XQIITE/9vH8MkVNOOUVLlizRV199pdraWtXW1nb3H0GvFBUVpbvvvlt33XWXXnzxRe3evVvvvfeenn/+eV1zzTWKiorS9OnTtX37dm3atEl/+ctfdO211/ov0fwcy5Yt05o1a7Rz507dcsst+vbbbzVjxozjfr2CggIVFhbq8ccf12effaaPP/5YRUVFeuSRRyRJ06ZNk8vl0syZM1VRUaHXX3+db+/uwbhMg26RlJSkLVu26O6779aUKVPU3NyskSNHaurUqQoLC1NsbKxKS0v12GOPqb6+XiNHjtTSpUuVnZ0tSZoxY4Y+/PBD5ebmKiIiQrfffrsmTZp0zGMuXbpUc+fO1bPPPquTTz5Ze/fuPep2+fn5/jvsJcnj8UiSNm3apMzMzG55//1NT/28N27cqF27dmnXrl0aMWJEu+f+/2Wh/uq+++5TRESE8vPzVVNTo8TERM2aNUvR0dEqKSnRnDlz9Otf/1rR0dG6/PLL/V/cf65FixZp0aJF8nq9GjVqlF577TWdeOKJx/16N954o6Kjo/Xwww9r3rx5GjhwoMaOHev/lvGYmBitW7dOs2bNksfjUVpamhYvXqzLL7+8W94PupfL4W8oAAAwxGUaAABgihgBAACmiBEAAGCKGAEAAKaIEQAAYIoYAQAApogRAABgihgBAACmiBEAAGCKGAEAAKaIEQAAYIoYAQAApv4XGVr+5nicYCkAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "z = m.values[0]\n", + "jac = np.array([z, (1-z)])\n", + "total = model(value, z)\n", + "total_err_sta = np.einsum(\"i,j,ij\", jac, jac, cov_sta) ** 0.5\n", + "total_err_sys = np.einsum(\"i,j,ij\", jac, jac, cov_sys) ** 0.5\n", + "\n", + "print(f\"total = {total:.2f} +/- {total_err_sta:.2f}(sta) + {total_err_sys:.2f}(sys)\")\n", + "\n", + "plt.errorbar((\"result 1\", \"result 2\", \"combined\"), value + [total], error_sta + [total_err_sta], fmt=\"o\")\n", + "plt.errorbar((\"result 1\", \"result 2\", \"combined\"), value + [total], error_sys + [total_err_sys], lw=3, fmt=\"none\")\n", + "plt.xlim(-1, 3);" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Note how the systematic uncertainty gets barely reduced by the combination, a consequence of the strong correlation. Try running this example with zero correlation to see how the uncertainty becomes smaller. " + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "py311", + "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.11.5" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/doc/tutorials.rst b/doc/tutorials.rst index 6dbf041e..c5356327 100644 --- a/doc/tutorials.rst +++ b/doc/tutorials.rst @@ -13,6 +13,7 @@ Important for most users are only the first two entries. notebooks/basic notebooks/cost_functions + notebooks/correlated_data notebooks/error_bands notebooks/interactive notebooks/simultaneous_fits