diff --git a/examples/example_combined_fit.ipynb b/examples/example_combined_fit.ipynb new file mode 100644 index 00000000..0f9e0acd --- /dev/null +++ b/examples/example_combined_fit.ipynb @@ -0,0 +1,151 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "id": "ethical-frontier", + "metadata": {}, + "outputs": [], + "source": [ + "import pyerrors as pe\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "incredible-posting", + "metadata": {}, + "outputs": [], + "source": [ + "x_test = {'a':[0,1,2,3,4,5],'b':[0,1,2,3,4,5]}\n", + "y_test = {'a':[pe.Obs([np.random.normal(i, i*1.5, 1000)],['ensemble1']) for i in range(1,7)],\n", + " 'b':[pe.Obs([np.random.normal(val, val*1.5, 1000)],['ensemble1']) for val in [1.0,2.5,4.0,5.5,7.0,8.5]]}\n", + "for key in y_test.keys():\n", + " [item.gamma_method() for item in y_test[key]]" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "subtle-malaysia", + "metadata": {}, + "outputs": [], + "source": [ + "def func_a(a, x):\n", + " return a[1] * x + a[0]\n", + "\n", + "def func_b(a, x):\n", + " return a[2] * x + a[0]\n", + "\n", + "funcs_test = {\"a\": func_a,\"b\": func_b}" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "modern-relay", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Fit with 3 parameters\n", + "Optimization terminated successfully.\n", + "chisquare/d.o.f.: 0.8434205014773611\n", + "fit parameters [1.01510812 0.98190604 1.45453441]\n" + ] + } + ], + "source": [ + "output_test = pe.combined_fits.combined_total_least_squares(x_test,y_test,funcs_test)" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "technological-rolling", + "metadata": {}, + "outputs": [], + "source": [ + "output_test.gamma_method()" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "persistent-mathematics", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Goodness of fit:\n", + "χ²/d.o.f. = 0.843421\n", + "Fit parameters:\n", + "0\t 1.015(32)\n", + "1\t 0.982(32)\n", + "2\t 1.455(41)\n", + "\n" + ] + } + ], + "source": [ + "print(output_test)" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "wooden-potential", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWoAAAD4CAYAAADFAawfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAuMUlEQVR4nO3deVzUdf7A8ddHRcHKymN3MW3R1lxM1BKPykTNY1srNc0y7Vem4p0X2rWrZOXuepWVqHhWm5kHollbHiF2qAWGJ5WKRxogHWIqoMjn98cHDBFlgJn5fmfm/Xw85uEMDDPvqXz34fP9vN9vpbVGCCGEfVWwOgAhhBBXJ4laCCFsThK1EELYnCRqIYSwOUnUQghhc5Vc8aI1a9bUQUFBrnhpIYTwSomJiT9prWsV9z2XJOqgoCASEhJc8dJCCOGVlFJHrvQ92foQQgibk0QthBA2J4laCCFsziV71MU5f/48x44dIzs7211vaWv+/v7UqVMHPz8/q0MRQtic2xL1sWPHuO666wgKCkIp5a63tSWtNT///DPHjh2jXr16VocjhLA5t219ZGdnU6NGDZ9P0gBKKWrUqCG/XQghHOLWPWpJ0r+TfxZCCEfJxUQhhHCCzz77jKlTp7rktSVRCyFEOZw6dYphw4bRtm1b5s2bx5kzZ5z+Hj6VqF9//XWCg4Pp27ev1aEIIbzABx98QKNGjZg3bx5jxoxh165dXHPNNU5/H7ed+rCDqKgoNm7cSJ06dawORQjhwdLT03n66adZvnw5ISEhxMTE0LJlS5e9nzWJevRoSEpy7ms2awavvXbFbw8ZMoSUlBTuu+8+jh49yj//+U8iIiIAaNy4MevWrQPgvvvuo02bNnz55ZfcdNNNrFmzhoCAAA4cOMCQIUPIyMigYsWKrFixgltuueWy9zl9+jTdunXj119/5fz587z88st069bNuZ9VCGEJrTVvv/02Y8aM4cyZM7z88suMHz+eypUru/R9fWbrY+7cudSuXZu4uDjGjBlzxeft37+f4cOHs3fvXm644QZWrVoFQN++fRk+fDg7d+7kyy+/JDAwsNif9/f3Z/Xq1ezYsYO4uDjGjRuHzKUUwvMdOnSILl268OSTT3Lbbbexc+dOXnjhBZcnaXBwRa2UGgMMBDSwG+ivtS77IeCrrHytVq9ePZo1awZA8+bNOXz4ML/99hvHjx+nR48egEnGV6K15vnnn2fLli1UqFCB48ePk56ezp/+9Cd3hC+EcLILFy4wa9Ys/vnPf1KxYkWioqIYPHgwFSq4b51bYqJWSt0EPA000lpnKaWWA48CS1wcm8tUqlSJvLy8i48LF55UqVLl4v2KFSuSlZVVqtd+9913ycjIIDExET8/P4KCgqSwRQgPtWvXLgYOHMjXX3/N/fffT1RUFHXr1nV7HI7+L6ESEKCUqgRUBX50XUiuFxQUxI4dOwDYsWMHhw4duurzr7vuOurUqUNsbCwAOTk5nD17ttjnZmZm8oc//AE/Pz/i4uI4cuSKLWaFEDaVnZ3NP/7xj4u/VS9btoy1a9dakqTBgUSttT4OTAeOAqlAptZ6fdHnKaXClVIJSqmEjIwM50fqRD179uSXX37htttu48033+TWW28t8WfeeecdXn/9dZo0acJdd91FWlpasc/r27cvCQkJhISE8Pbbb/PXv/7V2eELIVzos88+o1mzZrzyyiv07duX5ORkHnnkEUuriVVJF7qUUjcCq4BHgJPACmCl1vq/V/qZ0NBQXXTCS3JyMsHBweWN16vIPxMh7OPUqVM888wzzJ07l6CgIObNm0fnzp3d9v5KqUStdWhx33Nk66MjcEhrnaG1Pg/EAHc5M0AhhLBSQeFKdHQ0Y8aMYc+ePQ4n6cjISJRSl90iIyOdFp8jpz6OAq2VUlWBLOBewOcHIu7evZvHH3/8kq9VqVKF7du3WxSREKK0nFG4EhkZSWRkJO3atQNg8+bNTo+zxESttd6ulFoJ7ABygW+AaKdH4mFCQkJIcnbRjhDCLbTWvPXWW4wdO9athStl5dA5aq31JGCSi2MRQgiXS0lJYfDgwWzcuJE2bdowf/5821/095nKRCGEb7tw4QIzZ84kJCSE7du3ExUVRXx8vO2TNPhYUyYhhG+yS+FKWcmKWgjhtexWuFJWPpWo3dGPevPmzdx///0ue30hhGPcXbiSk5NDUlLSFYvhysOnEnVUVBQbNmzg3XfftToUIYSLnDp1iqFDh9K2bVtycnL45JNPWLJkCTVq1HDp+x45coTMzEwmT57s9Ne2ZI969OjRTj/a1qxZM16zQT9qMP+hdO3alQMHDtC+fXuioqLc2mlLCF/1wQcfMHToUFJTUxkzZgwvvfSSSyauFBYQEHBJ47U5c+YwZ84c/P39S93U7Up8Jnu4qx81wFdffcUbb7zBvn37OHjwIDExMU7/PEKI36Wnp/PII4/w4IMPUr16dbZu3crMmTNdnqTBHPd77LHHLi7GqlatSt++fUts9lYalqyor7bytVp5+1EDtGzZkvr16wPQp08fPv/8c3r16uXSuIXwRXYoXAkMDKRatWrk5eVRoUIFsrOzqVatmlN70PvMirqw0vSjzs3NLfXrF71YYWXXLSHcyR19LwqkpKTQuXNn+vfv7/aJK0Wlp6dTu3Ztbr/9doYMGeL0C4o+mahd2Y8azNbHoUOHyMvL4/3336dNmzZOi10IO4uMjERrTVhYGGFhYWit0Vo7NVHbsXAlJiaGBg0acO211zJ79mynb3f6ZKJ2ZT9qgBYtWjBixAiCg4OpV6/exS0TIUT57Nq1izvvvJNx48bRoUMH9u7dy9ChQ73+Yr1PVSYePnz44v316y+bfQDAnj17Lt4vOBUC0KBBAz799NMS36Ndu3Zs2bKl7EEKIS6TnZ3NSy+9xNSpU7nxxhtZtmwZvXv39pltRZ9K1EIIz/PZZ58xaNAgvvvuO5544glmzJjh8jPRdiOJuoykH7UQrlV04sonn3zi1okrdiKJuoykH7UQV5aTk0NycjJpaWllOqa2du1ahg0b5tbCFTvz7h14IYQlylpOXVC40q1bN7cXrtiZJGohhNMEBASglCI1NRUw5dRKKQICAq76c1prlixZQnBwMLGxsbz88sskJCSUeiyWt7Jnoo6MBKUuv7ng0LwQwnnKUk5tp8KVsigo8omPjyc+Pt4lRT4lJmqlVEOlVFKh2yml1GinRVCcyEjQGsLCzE1rc5NELYStlaacOjc313aFK2VRUORT9ObWKeRa6++AZgBKqYrAcWC10yK4mpwcSE6GtDRwYt28EMJ1CsqpAwMDadWq1cVtkMJ27drFgAEDSEhI8MiJK+5W2q2Pe4GDWusjrgjmMkeOQGYmOKm/a/fu3WnevDm33XYb0dE+P0hdCJe4Wjl1dnY2L7zwAs2bN+fIkSMeO3HF3Up7PO9R4L3ivqGUCgfCAW6++ebyRRUQAIUaJTFnjrn5+0M5+rsuWrSI6tWrk5WVRYsWLejZs6fPHZwXwipSuFJ2Dq+olVKVgQeBFcV9X2sdrbUO1VqH1qpVq3xRpaTAY49BQf1+1arQty+Us7/r66+/TtOmTWndujU//PAD+/fvL1+cQogSZWZmWjJxxZuUZkV9H7BDa53uqmAuCgyEatUgL88k6+xs87gc+9SbN29m48aNbN26lapVq9KuXbtL2psKIZxPClecozSJug9X2PZwifR0qF3bJO1WraCYCxKlkZmZyY033kjVqlX59ttv2bZtm5MCFUIUde7cOQ4cOEC3bt0ICQkhJiZGzkSXg0NbH0qpa4BOgPtmSsXEQIMGcO21MHu2eVwOf/vb38jNzSU4OJhnn32W1q1bOylQIUSBSZMmoZRi69atZGRkAKYvzkcffWRxZJ7NoRW11voM4NEbSlWqVOF///uf1WEI4bVSUlL48ssvAWjTpg3z58/3uDPRdmXvysT4eHOTykQhbCs3N5cZM2bQuHFjjy5csTN7ds+LjJSkLIQHkMIV97DniloIYWtSuOJe9lxRCyFsa8uWLQwaNIjvv/9eClfcRFbUQgiHFBSuhIWFce7cOdavXy+FK25iy0Rd0Daw6M2Z3aiEEI5bu3btxR45Y8eOZc+ePXTq1MnqsHyGbRO11pqwsDDCwsKc0jbw8OHDNG7c2HlBCuEDipu4MmPGDKkudDNbJuoCOTk5JCUlkZaWZnUoQvgUrTWLFy+WiSs2YetEXda5a1eSm5tL3759CQ4OplevXpw9e9YpryuENymYuPLUU0955MQVb2TLRF3WuWsl+e677xg2bBjJyclUq1aNqKgoZ4QrhFeQwhX7smWiLsvcNUfUrVuXu+++G4B+/frx+eeflztWIbzBzp07ufPOO4mIiKBjx47s27ePoUOHXvw7KKxly38LpZm7VhpKqas+FsLXFBSuhIaGXixcWbNmDXXq1LE6NFGIbQteHJm7VlpHjx5l69at3HnnnSxdupQ2bdo4IVIhPFPhwpUnn3yS6dOny5lom7LlihquPnetrBo2bMjs2bMJDg7m119/ZejQoU6IVAjPUlzhyuLFiyVJ25htV9TOFhQUxLfffmt1GEJYas2aNQwbNoy0tDTGjh3L5MmT5Uy0B7DlirqgMjE+Pp74+HipTBSinNLT0+nduzfdu3enRo0aUrjiTAVtmYvenJivlNbaaS9WIDQ0VCckJFzyteTkZIKDg53+Xp5M/pkIV9Nas2TJEsaNG8eZM2eYOHEiEyZMwM/Pz+rQvE+7dubPzZvL9ONKqUStdWhx33Pr1ofWWk5a5HPF/yCFKCwlJYXw8HA2bdokE1c8nKMzE29QSq1USn2rlEpWSt1Z2jfy9/fn559/lgSFSdI///wz/v7+VocivFDhwpWvvvpKCle8gKMr6lnAx1rrXkqpykDV0r5RnTp1OHbs2MWBl77O399fzqoKp9u5cycDBw4kISGBBx54gKioKPnvzAuUmKiVUtcDbYEnAbTW54BzpX0jPz8/6tWrV9ofE0I4IDs7m5deeompU6dSvXp13n//fR5++GHZanSnnBxIToa0NChncV5Rjmx91AMygMVKqW+UUguUUpddKlZKhSulEpRSCbJqFsJ9tmzZQtOmTZkyZQr9+vVj37599O7dW5K0ux05ApmZ4KQmcoU5kqgrAXcAc7TWtwNngGeLPklrHa21DtVah9aqVcvJYQohisrMzGTIkCGEhYVx/vx5KVyxSkCAOY5XUD09Z455XM4mcoU5kqiPAce01tvzH6/EJG4hhEXWrFlDo0aNmD9/PmPHjmX37t0yccXdTp+GN9+E2rUv/XpAAPTtC+VsIldYiYlaa50G/KCUapj/pXuBfU6LQAjhsLS0tIuFKzVr1mTbtm1SuOJuhw7B2LFw000wciTUqgX33mu+V6GC2auuVs2p+9SOViaOBN5VSu0CmgFTnBaBEKJEBRNXGjVqxJo1ay5OXGnRooXVofkGrU0hS48e8Je/wBtvwN//Dlu3wrZtJjHXrg233w5DhpgLik7k0PE8rXUSUGzFjBDCtQoXrtxzzz1ER0fLmWh3yc6GpUth1izYtQtq1IBnn4Vhw8yKukBMzO+VibNnOz0MW/b6EMIbFPSsKXpztGdN0cKVOXPmsHnzZknS7vDjj/CPf0DdujBggFlRL1gAP/wAr7xyaZJ2A7f1+hDCV7XLX2ltLkUPCClcschXX8Frr8GKFXDhAjzwAIwebVbLJR13dGGvD1lRC2Ej2dnZPP/88zRv3pyjR4/y/vvvy8QVVzt/HpYtgzvvhFat4MMPYcQI2L8f1qyB9u1LTtIu5jP9qIWwu6ITV2bMmEH16tWtDst7/fQTREdDVBQcP24uEr7+Ojz5JFx3ndXRXUJW1EK4WE5ODklJSaRd4SRA0cKVDRs2sHjxYknSrrJ7NwwcaPafX3gBGjWCdevgu+/McTubJWmQRC2Eyx05coTMzEwmF1NaXFzhSseOHS2I0stduGC2MTp0gCZNzEmOJ56AvXth/Xro2tWcgS6LgsEB8fHm5smDA4TwNQEBAWRnZ1/2dX9/fw4dOsTIkSNZuXIlTZo0YcGCBXIm2hUyM2HxYnPuOSXFrKJHjDArapv9xmKbwQFC+JKUlBQiIiJYtmwZeXl5VK1ale7du9OiRQsaNWrE2bNneeWVVxg/frxMXHG2/ftNcl682JR63303/PvfpmClkuelPc+LWAgPERgYSLVq1cjLy6NChQpkZWURHx/P0qVLueeee5g/fz4NGzYs+YWEY7SGjRtNccpHH5mE/OijMGoUNG9udXTlIolaCBdKT08nMDCQChUqkJaWRnp6OnPmzCE8PJwKZd0TFZc6exbeecec2Ni3D/7wB5g40ZRyO7kvtFUkUQvhQhMnTuSee+7h9OnTPPjgg8yePVvORDvL0aOmXHv+fPj1V7jjDnjrLXjkEahSxeronEoStRAukJ2dzeTJk5k6dSoVKlQgODiY2NhYaeZfXlrDF1+Y7Y3Vq83jhx4y2xt33215YYqrSKIWwsni4+MJDw+/WLjy/fff4+fnJ0m6PHJyYPlyk6ATE+GGG0yr0eHD4c9/tjo6l5NNMiGcJDMzk8GDB9OuXbtLClfkREc5pKfDiy+aZPx//2f2o+fMgWPHYOpUn0jSICtqIZxizZo1DBs2jLS0NMaNG8eLL74ozfzLY8cOs3petgzOnTO9n0eNgk6dvHZ742pkRS1EOaSlpfHwww9fMnFl+vTpXHPNNRfbnMbHxxMfH1/qNqc+JzcXVq6Ee+4xx+lWrYLwcFPa/eGH0LmzTyZpkMpEIcqkYOLKuHHjyMrKYuLEiVK4Ula//mp6Pb/5pjnJERRkem489ZTZi/YRUpkohBMdPHiQ8PBwPv30UylcKY/kZHP2+e23zd5zu3Zmu+OBB6BiRaujsxWHtj6UUoeVUruVUklKKVkqC5+Um5vL9OnTCQkJ4euvv744cUWSdCEFDYqK3gq2e/LyTNVgly6ma93ixaZ6MCkJ4uKge3dJ0sVwaOtDKXUYCNVa/+TIi8rWh/A2SUlJDBw4kMTERB588EGioqK4yc3jmDxK0Wknp0/DkiWm/8b335tBsMOGmT3oWrUsCtJeZOtDiDLKysripZdeYurUqdSoUYPly5fTq1cvORPtqEOHTHJeuBBOnTITVJYuhZ49oXJlq6PzGI4mag2sV0ppYJ7WOrroE5RS4UA4wM033+y8CIWwSHx8PIMGDWL//v0ycaU0tIaMDLNyvuUWs5XRq5c5Xte6tdXReSRHj+e10VrfAdwHDFdKtS36BK11tNY6VGsdWkt+lREerHDhSm5urkxccVR2NixaBM2ameZIublw++1w+DC8954k6XJwaEWttT6e/+cJpdRqoCWwxZWBCWEFKVwpgx9/NHMH580zcwgL27ED6tQBf3/IyrImPi9Q4opaKXWNUuq6gvtAZ2CPqwMTwp2uVrgirmD7dnjsMVPGPWWKaYq0YgX06fP7WKuqVaFvX7NXLcrMkRX1H4HV+RdPKgFLtdYfuzQqIdykaOGKTFwpwfnzpnpw1iyTqKtVM8UpI0ZA/frmOZs2mWN4FSqY7ZBq1bymL7RVSkzUWusUoKkbYhHCraRwpRR++gmio03/5x9/hAYNzGmOJ564fGp3ero5fhcYaE55pKZaE7MXkeN5wufk5uby2muvMXHiRPz8/Jg7dy6DBg2SiSvF2b3brJ7ffdesjjt1Mgn7vvuuPLU7Jub3c9SzZ7stVG8miVr4FClcccCFC7BunUnQcXEQEGBWzk8/baoJhdvJEkL4hKysLJ577jlCQ0P54YcfWL58ObGxsZKkC8vMhFdfhVtvNaXcBw7Af/5jej/PnStJ2kKyohZer3DhSv/+/Zk+fbqciS5s/37THGnJElPq3aaNSdDdu5tJ3sJysqIWXqu4wpVFixZJkgZTPbh+PXTtalbQ8+ZBjx6QkACffWYqCcuSpAuaMsXHm1vRpkyiTKQftfBKsbGxDB8+nLS0NMaMGSOFKwXOnoV33jEr6H374A9/gKFDYcgQOUJnMWnKJHxGWloaI0eOZOXKlTRp0oTY2FhatGhhdVjWO3rUnMCYP9806r/jDnjrLXjkEahSxeroRAkkUQuvIIUrxdAavvjCnN5Yvdo8fugh0xzp7rt9dqyVJ5JELTyeFK4UkZMDy5ebBJ2YaMZZjR0Lw4f7zNRubyOJWngsKVwpIj3dHKObM8fcDw429x9/HGR/3qNJohYeSQpXCtmxw6yely2Dc+fg73832xudOsn2hpfw0aWH8FRSuJIvN9c0R7rnHmjeHFatMmOtvvsOPvwQOneWJO1FZEUtPIYUrgC//AILFpgTHEePQr16MHMmPPUUXH+91dEJF5FELWwvMzOTCRMmEB0dTb169diwYQMdO3a0Oiz32rfPnH1+5x1zFrp9e/P4/vtlarcPkEQtbC02NpZhw4aRnp7uexNX8vLg44/N/vP69ea8c9++Zv+5SROroxNuJHvUwm0iIyNRSl12iyymvLhg4kqPHj2oVasW27dv952JK7/9Bm++aU5tdO0Ke/bAyy/DDz+Yad6SpH2OlJALt2uX36t48+bNl32vaOHKpEmTiIiI8I3ClZQUk6AXLoRTp0zT/VGjTN8NX/j8Ps4pJeRKqYpAAnBca32/s4ITvicnJ4fk5GTS0tL4U6H+Ej5ZuKI1bN5stjfWrjX7zQ8/bBJ0q1ZWRydsojRbH6OAZFcFInzHkSNHyMzMZPLkyYApXJk2bRohISEkJCQwd+5cNm/e7PlJuqCTXNFbZKSZyL1wITRrBh06mFLv55+Hw4dh6VJJ0uISDm19KKXqAG8BrwBjS1pRy9aHKE5AQADZ2dmXfV0phdbaewtXCsZSbd5s5g1GRZm2oj/9BCEhMHq0mdwdEGBhkMJqV9v6cHRF/RowAchzVlDC96SkpPDYY49dLPGulN/vuHr16t5duJKTY/o89+hhem1MmWKaIn36Kezcac5AS5IWV1HiHrVS6n7ghNY6USnV7irPCwfCAW6++WZnxSe8SGBgINWqVSMvz/z/Pjc3l4YNG/Lll196Z+HK+fOmejAx0dz/8EMYORJGjID69a2OTniQErc+lFL/Ah4HcgF/oBoQo7Xud6Wfka0PUZzMzEyaNWvG4cOHqVy5Ml26dKFSpUrExMRYHZpzZWSYSd3/+Efx3/f3N3vUQhRSrq0PrfVzWus6Wusg4FHg06slaSGKExsbS3BwMEePHqVOnTq0bNmStWvXeleS3rULBgyAunVNkg4LM7eCbn5Vq5qClUOHrI1TeBwpeBEuVVzhyi233EJFbyl7vnAB1qwxJd1Nm5oOdv37w9695uJhcLCpMKxQAbKzoVo1GXklSq1UJeRa683AZpdEIrxK0cKVKVOmcPbs2UvGYqn87m6TJk0qtjrR1jIzYdEieOMNs0KuW9dM7h44EArvt6enQ+3aEBhojtylploXs/BYUpkonM6rC1e+/94k5yVL4PRpaNPGFKd0737lqd2Fj+cJcQUy3Fa4RW5uLq+++iqTJk3yrokrWsOGDaZ68KOPoHJlePRRk6DvuMPq6IQP8PC/QcIukpKSaN26NRMmTKBz587s27ePwYMHe3aSPnPGjLa67Tbo0sUcs4uMNH2g33qr5CRdUJkYH29uhSsThSgFWVGLcsnKymLy5MlMmzaNGjVqsHz5cnr16nVx/9kjHT1qmiMtWAC//momqLz9NvTubVqNOioyUpKycApJ1KLMvGriitam38asWRATY1a+Dz1ktjfuukvGWglLefDvpcIqJ0+eZPDgwbRr144LFy6wYcMGFi1a5JlJOifHrJZDQ838wU2bICLCtBxdvtyUekuSFhaTFbUolcITVyIiInjxxRepWrWq1WGVXlqa2X+eO9ccoQsONvf79QNfGE4gPIokauGQtLQ0Ro4cycqVK2nSpAlr164lNLTYk0T2lphotjeWLTP9N7p2NdsbHTvKylnYliRqcVXFFa543MSV3FxYvdok6C++gGuvhcGDTYOkW2+1OjohSiSJWlxR4cKVtm3bMn/+fG71pMT2yy8wfz7Mnm3mDdarBzNnmrai119vdXRCOEwStbiMxxeu7NsHr79uLhJmZZk+HG+8Afffb0ZdCeFhJFGLSyQlJTFw4EASExPp1q0bs2fP9oxm/nl58L//me2NDRvMeed+/eDpp2Vqt/B4HrJEEq6WlZXFc889R2hoKMeOHWPFihWsXr3a/kn6t9/MavmvfzUr5r174ZVX4NgxU7AiSVp4AVlRi0sKV5566immTZtm/zPRKSkmQS9aBKdOmc50770HPXuCJ13oFMIBkqh92MmTJ3nmmWeIjo6mfv36bNy4kXvvvdfqsK5Ma4iLM9sbH3xg9psfftgcr5Op3cKLSaL2UR5VuJKVBUuXmgS9ezfUrAnPPw9Dh4Ldt2aEcAJJ1D6mcOFK06ZN7V24cvw4REXBvHnw889mv3nhQujTR6Z2C58iidpHaK1ZtGgRERER9i9c2bbNrJ5XrjSjrrp1M9sbYWFSPSh8UomJWinlD2wBquQ/f6XWepKrAxPO4xGFK+fOmcQ8axZ89ZWZLfj00zBihClUEcKHOXI8LwfooLVuCjQD/qaUau3SqIRT5ObmMm3aNEJCQkhISGDevHnExcVZl6QLGukXvXXoAEFBZkL3r7+a0xzHjsGMGZKkhcCBRK2N0/kP/fJvzh+0KJwqKSmJVq1aXTJxJTw83NrqwshIc3KjdWvTb+ORR0xhSlwchITAhx/Ct9+aVfR111kXpxA249DfWqVURaVUEnAC2KC13l7Mc8KVUglKqYSMjAwnhykcVbhw5fjx4/YqXLlwAWJjYccOMxh25Uro39+UfH/yCfz97+ApZepCuFGpppArpW4AVgMjtdZ7rvQ8mUJuDdsWrmRmmtMaERFmRV2Uv785gieED7vaFPJSLV+01ieBOOBvTohLOEnRiSsbN25k4cKF1ifp77832xg33QTjxkHLltCmze+r5qpVzb70oUPWximEzZWYqJVStfJX0iilAoBOwLcujks4KDY2lkaNGrFgwQIiIiLYvXu3tdWFWsP69aYhf8OGps1oz56mYf+2bdC4sWmgVKECZGeb0x1/+pN18QrhARw5Rx0IvKWUqohJ7Mu11utcG5YoSVpaGiNGjGDVqlU0bdqUDz74gObNm1sX0Jkz8M47pr1ocjL88Y/m4uGQIeZ+gfR0qF0bAgNN2XdqqmUhC+EpSkzUWutdwO1uiEU4wHaFK0eOmMb88+fDyZPQvLnpA927tznRUVRMDLRrZ+7Pnu3OSIXwWHKJ3YMcPHiQjh07MnDgQJo0acKuXbt47rnn3J+ktYbPPoNevaB+fTM1pVMn+Pxz+PprePzx4pN0wTnq+HhzKzhHHRnp3viF8DClOvXhKDn14VxFJ65MmzaNgQMHuv9MdE6OGQo7axZ88w3ceCOEh8OwYXDzze6NRQgvc7VTH9Lrw+aSkpIYMGAAO3bssG7iSloazJkDc+fCiRPQqJFplNSvnzm5IYRwKUnUNpWVlcXkyZOZNm0aNWvWZMWKFfTs2RPlzqZEiYlm9bxsGZw/b05yjBoFHTtKcyQh3EgStQ0VLVyZPn06N954o3vePDcXVq82CfqLL0yp95AhMHIkNGjgnhiEEJeQi4kWiYyMRCl12a158+aXFa64JUn/8gv85z/m4mDv3ubY3KuvmuZIr78uSVoIC0mitkhkZCRaawIDAwHo0qULgYGBJCUlubdwZe9eGDwY6tSBZ581CXnNGlNVOHo0XH+962MQQlyVnPqwSEBAANnZ2Zd9vXLlyuTk5Lj2zfPy4KOPzPbGxo2m10a/fqb/c0iIa99bCFEsp/X6EM5z8OBBWhUayOrn50efPn04cuSI6970t99Mr+eGDeGBB0wF4ZQp8MMPpmBFkrQQtiQXEy1w4MABBg8ezPbtplusUooLFy5www038CdX9L1ISTEJetEiOHXK9IN+6SXTg8OOo7iEEJeQRO1GBYUrEydOpHLlyjRr1oz09HRq165Nq1atSHVm3wutTUP+WbPggw+gYkVzkXDUKNPFTgjhMSRRu0nhwpXu3btTv359Zs6cCUBqaiqJiYmAucgYWZ6S6qwsePddc1Jj926oWRNeeAGGDjXNkIQQHkcStYsVLVxZuXIlDz30EEopZsyY4bw3OnYMoqIgOhp+/hmaNjVbHX36mIuFQgiPJYnahTZv3sygQYM4cOCA6wpXtm2D116DVavMqKtu3cz2RliYVA8K4SXk1IcLnDx5kvDwcNq3b09eXp7zC1fOnYOlS00/5zvvhI8/NkfrDh40VYXt2kmSFsKLSKJ2stWrV9OoUSMWLlzI+PHjr1y4UtDys+jtavvTGRnw8ssQFGRGWJ08CW++abY9ZsyAevVc86GEENbSWjv91rx5c+1rfvzxR92zZ08N6KZNm+qEhATHfjAwUGvQeujQKz8nKUnr/v21rlLFPLdLF60/+kjrCxecE7wQwnJAgr5CTpUVdTlprVm4cCGNGjVi3bp1/Otf/+Lrr78ueSxWQIBZQRccyZszxzwOCDCPL1yA2Fho3x6aNYP334f+/WHfPrPVcd99vw+JFUJ4NUeG29ZVSsUppfYppfYqpUa5IzBPcODAAe69914GDhxI06ZN2b17N88++6xjE1dSUuCxxy6fyJ2UZCamNGgAPXqY502darY35syB4GCXfiYhhP04cuojFxintd6hlLoOSFRKbdBa73NxbLZVtHAlOjqaAQMGlG7iSmCgmcBdMJE7Kwt27jQzB8+cgXvugWnTzCmOSnI4Rwhf5shw21QgNf/+b0qpZOAmwCcT9TfffMOAAQP45ptv6N69O7Nnz6Z2WQtJ0tOhenXTA/rUKdPJ7v/+zxyvu13mCQshjFJtciqlgjATybcX871wpVSCUiohIyPDSeHZR1ZWFs8++ywtWrQgNTWVlStXEhMTU7YkfeaMmZayerXpA33qlPm61uZEhyRpIUQhDv9OrZS6FlgFjNZanyr6fa11NBANps2p0yK0gcKFKwMGDGDatGllOxN95Ig5TrdggTla17y5WT337l381G4hhMDBRK2U8sMk6Xe11jGuDck+Tp48yfjx41mwYAH169dn06ZNdOjQoXQvojV8/rlpjrR6tTnZ8dBDJkHfdZcUpgghSlRiolZmmupCIFlrPdP1IdnD6tWrGT58OOnp6YwfP57IyEiqlmbidk6OGQo7axZ88w3ceCOMHw/Dh0Pduq4LXAjhdRxZUd8NPA7sVkol5X/tea31Ry6LykKpqamMGDGCmJgYmjVrxrp167jjjjscf4G0NHOMbu5cOHECGjWCefPMBJXSJHohhMjnyKmPzwGv//1ca82iRYuIiIggOzubf//734wdO9axM9EACQlm9fz+++YUR9euZnvj3ntle0MIUS5yQBdTuBIeHk5cXBxhYWHMnz+fBo5M3c7NhZgYk6C//BKuvdb0fR45Ev7yF9cHLoTwCT6dqHNzc5k5cyaTJk2iSpUqjheu/PKLmTE4e7aZN1i/vmk12r+/KWIRQggn8tlEXbhwpUePHrz55psln4neu9dMTnnnHVNJ2KGDOW7XtasZdSWEEC7gc119srKyeOaZZxwvXMnLg3XroFMnaNwY3n7b9OTYtQs2bYIHH5QkLYRwKZ9aUZeqcOW332DxYjO9+8ABuOkmmDIFBg0ycwiFEMJNfGJFffLkSQYNGkT79u3RWrNp0yYWLFhQfJI+eBBGjzaJedQok5Tfew8OHYLnnpMkLYRwO69fURcUrpw4cYIJEyYwadKkywtXtIZPPzWnN9atM1sZvXubRN2ypTWBCyFEPq9N1A4VrmRlwX//ay4Q7tljVssvvGCO2JW1I54QQjiZ1yXqgokrERER5OTkFF+4cuyYOVoXHW2O2jVtCosWQZ8+4O9vXfBCCFEMr0rUhQtX2rVrR3R09O+FK1rDtm1me2PlSvO4WzezvdG2rVQPCiFsyysSdW5uLjNmzCAyMvLywpVz52DFCpOgv/4arr/eJOcRI2RqtxDCI3h8or5i4cqJE6YZ0pw5ZoDsrbea4pQnnjCl3kII4SE89nhe0cKVVatWmcKVEydMKffNN8PEidCkCXz0ESQnmxajkqSFEB7GIxN1XFwcISEhTJ06lf79+7Nv924eAggLM2Osli+Hp56Cffvg44/hvvt+n/YthBAexqO2PgpPXLnlllvYtGYNHb7/HkJDzZirP//ZTO4eMMA06hdCCC/gMYk6JiaG4cOHk5GRwYSBA4msUIGAxx4zg2LbtoWZM03fjUoe85GEEMIhts9qP/74IyNHjjSFK/Xr82HLltyxYAFUrmzOPY8aJVO7hRBezbaJWmvNggULGD9+PDlnz/LvmjUZm5KC39mz8OKLMHgw/PGPVocphBAuV+IVNqXUIqXUCaXUHpdHk5oKYWEc2rqVDnfdRXh4OLefPs2u8+d5JigIv3feMXvREydKkhZC+AxHVtRLgDeBt10bChyfMIE+W7Zw/1138Q0wXykG9OiBGj0a7rpLqgeFED7JkeG2W5RSQa4M4nzFivjl5fEK8AXQGDgJnAfUihWufGshhLA9px0uVkqFK6USlFIJGRkZpfrZan5+KGAOkJf/pwKqVa7srPCEEMJjOS1Ra62jtdahWuvQWrVqlepnU4YP5zGgoEt0VaAvcGjECGeFJ4QQHssW5XrzrruOn4AszEo6C8gA5kq5txBC2ON4XmRkJA/t2sXQwEDCw8OJjo4mNTWVyMhIq0MTQgjLKa311Z+g1HtAO6AmkA5M0lovvNrPhIaG6oSEBGfFKIQQXk8plai1Di3ue46c+ujj/JCEEEI4yhZ71EIIIa5MErUQQticJGohhLA5SdRCCGFzkqiFEMLmJFELIYTNlXiOukwvqlQGcKSMP14T+MmJ4XgC+czez9c+L8hnLq0/a62L7b/hkkRdHkqphCsd+vZW8pm9n699XpDP7Eyy9SGEEDYniVoIIWzOjok62uoALCCf2fv52ucF+cxOY7s9aiGEEJey44paCCFEIZKohRDC5myTqJVSf1NKfaeUOqCUetbqeNxBKbVIKXVCKbXH6ljcQSlVVykVp5Tap5Taq5QaZXVMrqaU8ldKfaWU2pn/mV+0OiZ3UUpVVEp9o5RaZ3Us7qCUOqyU2q2USlJKObUhvy32qJVSFYHvgU7AMeBroI/Wep+lgbmYUqotcBp4W2vd2Op4XE0pFQgEaq13KKWuAxKB7t7871kppYBrtNanlVJ+wOfAKK31NotDczml1FggFKimtb7f6nhcTSl1GAjVWju9yMcuK+qWwAGtdYrW+hywDOhmcUwup7XeAvxidRzuorVO1VrvyL//G5AM3GRtVK6ljdP5D/3yb9avjlxMKVUH6AossDoWb2CXRH0T8EOhx8fw8r/Avk4pFQTcDmy3OBSXy98CSAJOABu01l7/mYHXgAlAnsVxuJMG1iulEpVS4c58YbskauFDlFLXAquA0VrrU1bH42pa6wta62ZAHaClUsqrt7mUUvcDJ7TWiVbH4mZttNZ3APcBw/O3Np3CLon6OFC30OM6+V8TXiZ/n3YV8K7WOsbqeNxJa30SiAP+ZnEornY38GD+nu0yoINS6r/WhuR6Wuvj+X+eAFZjtnSdwi6J+muggVKqnlKqMvAosNbimIST5V9YWwgka61nWh2POyilaimlbsi/H4C5YP6tpUG5mNb6Oa11Ha11EObv8qda634Wh+VSSqlr8i+Qo5S6BugMOO00ly0StdY6FxgBfIK5wLRca73X2qhcTyn1HrAVaKiUOqaUGmB1TC52N/A4ZoWVlH/7u9VBuVggEKeU2oVZkGzQWvvEcTUf80fgc6XUTuAr4EOt9cfOenFbHM8TQghxZbZYUQshhLgySdRCCGFzkqiFEMLmJFELIYTNSaIWQgibk0QthBA2J4laCCFs7v8BTo9tJdcmdxMAAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "colour= {'a':'red','b':'black'}\n", + "plt.figure()\n", + "for key in funcs_test.keys():\n", + " plt.errorbar(x_test[key],[o.value for o in y_test[key]],ls='none',marker='*',color=colour[key],yerr=[o.dvalue for o in y_test[key]],capsize=3,label=key)\n", + " plt.plot([x_val for x_val in x_test[key]],[funcs_test[key](output_test.fit_parameters,x_val) for x_val in x_test[key]],color=colour[key],label='func_'+key)\n", + "plt.legend()\n", + "plt.show()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.7.3" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/pyerrors/__init__.py b/pyerrors/__init__.py index a6fefe8d..1949e68a 100644 --- a/pyerrors/__init__.py +++ b/pyerrors/__init__.py @@ -457,6 +457,7 @@ from .obs import * from .correlators import * from .fits import * from .misc import * +from . import combined_fits from . import dirac from . import input from . import linalg diff --git a/pyerrors/combined_fits.py b/pyerrors/combined_fits.py new file mode 100644 index 00000000..2f92a9c2 --- /dev/null +++ b/pyerrors/combined_fits.py @@ -0,0 +1,170 @@ +import iminuit +import autograd.numpy as anp +from autograd import jacobian +from pyerrors.fits import Fit_result +import numpy as np +import pyerrors as pe +from autograd import jacobian as auto_jacobian +from autograd import hessian as auto_hessian +from autograd import elementwise_grad as egrad +from numdifftools import Jacobian as num_jacobian +from numdifftools import Hessian as num_hessian +import scipy.optimize +import scipy.stats + +def combined_total_least_squares(x,y,funcs,silent=False,**kwargs): + r'''Performs a combined non-linear fit. + Parameters + ---------- + x : ordered dict + dict of lists. + y : ordered dict + dict of lists of Obs. + funcs : ordered dict + dict of objects + fit functions have to be of the form (here a[0] is the common fit parameter) + ```python + import autograd.numpy as anp + funcs = {"a": func_a, + "b": func_b} + + def func_a(a, x): + return a[1] * anp.exp(-a[0] * x) + + def func_b(a, x): + return a[2] * anp.exp(-a[0] * x) + ``` + It is important that all numpy functions refer to autograd.numpy, otherwise the differentiation + will not work. + silent : bool, optional + If true all output to the console is omitted (default False). + initial_guess : list + can provide an initial guess for the input parameters. Relevant for + non-linear fits with many parameters. + num_grad : bool + Use numerical differentation instead of automatic differentiation to perform the error propagation (default False). + ''' + + output = Fit_result() + output.fit_function = funcs + + if kwargs.get('num_grad') is True: + jacobian = num_jacobian + hessian = num_hessian + else: + jacobian = auto_jacobian + hessian = auto_hessian + + x_all = [] + y_all = [] + for key in x.keys(): + x_all+=x[key] + y_all+=y[key] + + x_all = np.asarray(x_all) + + # number of fit parameters + n_parms_ls = [] + for key in funcs.keys(): + for i in range(42): + try: + funcs[key](np.arange(i), x_all.T[0]) + except TypeError: + continue + except IndexError: + continue + else: + break + else: + raise RuntimeError("Fit function is not valid.") + n_parms = i + n_parms_ls.append(n_parms) + n_parms = max(n_parms_ls) + if not silent: + print('Fit with', n_parms, 'parameter' + 's' * (n_parms > 1)) + + if 'initial_guess' in kwargs: + x0 = kwargs.get('initial_guess') + if len(x0) != n_parms: + raise Exception('Initial guess does not have the correct length: %d vs. %d' % (len(x0), n_parms)) + else: + x0 = [0.1] * n_parms + + def chisqfunc(p): + chisq = 0.0 + for key in funcs.keys(): + x_array = np.asarray(x[key]) + model = anp.array(funcs[key](p,x_array)) + y_obs = y[key] + y_f = [o.value for o in y_obs] + dy_f = [o.dvalue for o in y_obs] + C_inv = np.diag(np.diag(np.ones((len(x_array),len(x_array)))))/dy_f/dy_f + chisq += anp.sum((y_f - model)@ C_inv @(y_f - model)) + return chisq + + if 'tol' in kwargs: + fit_result = iminuit.minimize(chisqfunc, x0,tol=kwargs.get('tol')) + fit_result = iminuit.minimize(chisqfunc, fit_result.x,tol=kwargs.get('tol')) + else: + fit_result = iminuit.minimize(chisqfunc, x0,tol=1e-4) + fit_result = iminuit.minimize(chisqfunc, fit_result.x,tol=1e-4) + + chisquare = fit_result.fun + + output.method = 'migrad' + output.message = fit_result.message + + if x_all.shape[-1] - n_parms > 0: + output.chisquare = chisqfunc(fit_result.x) + output.dof = x_all.shape[-1] - n_parms + output.chisquare_by_dof = output.chisquare/output.dof + else: + output.chisquare_by_dof = float('nan') + + if not silent: + print(fit_result.message) + print('chisquare/d.o.f.:', output.chisquare_by_dof ) + print('fit parameters',fit_result.x) + + # use ordered dicts so the data and fit parameters can be mapped correctly + def chisqfunc_compact(d): + chisq = 0.0 + list_tmp = [] + c1 = 0 + c2 = 0 + for key in funcs.keys(): + x_array = np.asarray(x[key]) + c2+=len(x_array) + model = anp.array(funcs[key](d[:n_parms],x_array)) + y_obs = y[key] + y_f = [o.value for o in y_obs] + dy_f = [o.dvalue for o in y_obs] + C_inv = np.diag(np.diag(np.ones((len(x_array),len(x_array)))))/dy_f/dy_f + list_tmp.append(anp.sum((d[n_parms+c1:n_parms+c2]- model)@ C_inv @(d[n_parms+c1:n_parms+c2]- model))) + c1+=len(x_array) + chisq = anp.sum(list_tmp) + return chisq + + fitp = fit_result.x + y_f = [o.value for o in y_all] # y_f is constructed based on the ordered dictionary if the order is changed then the y values are not allocated to the the correct x and func values in the hessian + dy_f = [o.dvalue for o in y_all] # the same goes for dy_f + try: + hess = hessian(chisqfunc)(fitp) + except TypeError: + raise Exception("It is required to use autograd.numpy instead of numpy within fit functions, see the documentation for details.") from None + + jac_jac_y = hessian(chisqfunc_compact)(np.concatenate((fitp, y_f))) + + # Compute hess^{-1} @ jac_jac_y[:n_parms + m, n_parms + m:] using LAPACK dgesv + try: + deriv_y = -scipy.linalg.solve(hess, jac_jac_y[:n_parms, n_parms:]) + except np.linalg.LinAlgError: + raise Exception("Cannot invert hessian matrix.") + + result = [] + for i in range(n_parms): + result.append(pe.derived_observable(lambda x_all, **kwargs: (x_all[0] + np.finfo(np.float64).eps) / (y_all[0].value + np.finfo(np.float64).eps) * fitp[i], list(y_all), man_grad=list(deriv_y[i]))) + + output.fit_parameters = result + + return output \ No newline at end of file