diff --git a/examples/08_combined_fit_example.ipynb b/examples/08_combined_fit_example.ipynb new file mode 100644 index 00000000..a492c674 --- /dev/null +++ b/examples/08_combined_fit_example.ipynb @@ -0,0 +1,134 @@ +{ + "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": "45f67973", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Fit with 3 parameters\n", + "Method: migrad\n", + "Optimization terminated successfully.\n", + "chisquare/d.o.f.: 0.8085703524653507\n", + "fit parameters [0.97737577 1.01063624 1.47900852]\n", + "chisquare/expected_chisquare: 0.8121288230401409\n" + ] + } + ], + "source": [ + "output_test = pe.fits.least_squares(x_test,y_test,funcs_test,method='migrad',expected_chisquare=True)" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "technological-rolling", + "metadata": {}, + "outputs": [], + "source": [ + "output_test.gamma_method()" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "wooden-potential", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWoAAAD4CAYAAADFAawfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA0/ElEQVR4nO3de3zO9f/H8ccbwxajRI2lUQ5DIguhEB2IEJWQIiFttmQ5xUZCyaE5bU79hJRYKYeQr0MizPmwOZ82w3JYTttse//+eA+Fsdl1XZ9r2+t+u123bLuuz/Wabp577/15v19vpbVGCCGE88pjdQFCCCHuTIJaCCGcnAS1EEI4OQlqIYRwchLUQgjh5PLZ46IPPvig9vLysselhRAiR9q8efPfWuvit/uaXYLay8uLiIgIe1xaCCFyJKXU0fS+JlMfQgjh5CSohRDCyUlQCyGEk7PLHPXtXL16lejoaBISEhz1lk6tYMGCeHp64uLiYnUpQggn57Cgjo6OpnDhwnh5eaGUctTbOiWtNWfOnCE6OpoyZcpYXY4QwsllaOpDKeWvlNqllNqtlAq4lzdKSEigWLFiuT6kAZRSFCtWTH67EEJkyF2DWilVBXgfqAk8CTRTSj1+L28mIX2D/F0IITIqIyNqb2CD1vqy1joZWA28Zt+yhBBCXJORoN4FPKuUKqaUcgOaAo/c/CSlVFelVIRSKiIuLs7WdQohhFMKDg5GKXXLIzg42Gbvcdeg1lpHAl8Ay4DfgG1Aym2eN1lr7aO19ile/La7IC0XEhKCt7c37du3t7oUIUQOERwcjNaa+vXrU79+fbTWaK1tGtQZWvWhtZ4GTANQSg0Dom1WgQNNnDiR33//HU9PT6tLEUKIDMtQUCulSmitTyulSmPmp2tn6V0DAmDbtixd4hbVqsHYsel+uXv37hw6dIgmTZpw7NgxBg4cSO/evQGoUqUKCxcuBKBJkybUq1ePdevWUapUKRYsWICrqysHDhyge/fuxMXFkTdvXn788Ucee+yxW97n4sWLtGjRgnPnznH16lWGDh1KixYtbPu9CiFylYzuTJyvlNoD/Ap8qLU+b7+S7CM0NJSSJUuycuVKPvroo3Sft3//fj788EN2795N0aJFmT9/PgDt27fnww8/ZPv27axbtw4PD4/bvr5gwYL89NNPbNmyhZUrV/Lxxx8j51IKIbIio1Mfz9r0Xe8w8rVamTJlqFatGgA1atTgyJEjXLhwgZiYGFq1agWYME6P1pr+/fuzZs0a8uTJQ0xMDKdOneLhhx92RPlCiBzIYTsTnUm+fPlITU29/vG/N54UKFDg+p/z5s3LlStXMnXt2bNnExcXx+bNm3FxccHLy0s2tgghsiRXNmXy8vJiy5YtAGzZsoXDhw/f8fmFCxfG09OTn3/+GYDExEQuX7582+fGx8dTokQJXFxcWLlyJUePpttiVgiRgyQmJrJt2zZOnjxp82vnyqBu3bo1Z8+epXLlyowfP57y5cvf9TUzZ84kJCSEqlWrUqdOnXT/Z7Rv356IiAieeOIJvv32WypWrGjr8oUQTujo0aPEx8czZMgQm19b2eNGl4+Pj775hJfIyEi8vb1t/l7ZmfydCJH9ubq63nZ6s2DBgpmaOlVKbdZa+9zua7lyRC2EELZy6NAh2rVrR548Jk7d3Nxo3779XadUMyNX3ky0hZ07d/L222//53MFChRgw4YNFlUkhLCCh4cH7u7upKamkidPHhISEnB3d7fpSi8J6nv0xBNPsM3Wm3aEENnSqVOnKFmyJB4eHtSqVYvY2FibXl+CWgghsig8PJwGDRoAMGHCBJtfX4JaCCFsQGtNUlKSXa4tQS2EEFm0bds2tm3bRlJSEleuXMHV1dWm15dVH0IIcY/Onj1Ljx49qFGjBleuXKF06dL/2d1sK7kqqB3Rj3rVqlU0a9bMbtcXQlgvJSWFsLAwypcvT1hYGB9++CFPP/00Hh4e15fp2VKuCuqJEyeyfPlyZs+ebXUpQohsat26ddSsWZPu3btTuXJltm7dSkhICC4uLnZ7T0vmqAMCAmy+tK1atWqMdYJ+1AD//PMPr7zyCgcOHKBhw4ZMnDjRLj9lhRCOc/LkSfr06cO3335LqVKlmDNnDm+++aZDDqrONenhqH7UABs3bmTcuHHs2bOHgwcPEh4ebvPvRwjhGFevXmX06NGUL1+eOXPm0LdvX6Kiomjbtq1DQhosGlHfaeRrtaz2owaoWbMmZcuWBeCtt95i7dq1tGnTxq51CyFs7/fff6dnz55ERkbSpEkTvv76a8qVK/ef5wQHBzN48ODrH18L76CgIJudm5ihEbVS6iOl1G6l1C6l1Byl1J2Tysllph91cnJypq9/809ZR/3UFULYxtGjR2nTpg0vvPACiYmJ/PLLLyxatOiWkIYbh9ve/HDoKeRKqVJAT8BHa10FyAu0tVkFFrBnP2owUx+HDx8mNTWVH374gXr16tmsdiGE/Vy5coUhQ4bg7e3N4sWLGTp0KLt376Z58+aWDrgyOkedD3BVSuUD3IAT9ivJ/uzZjxrg6aefxtfXF29vb8qUKXN9ykQI4Zy01ixYsIDKlSsTFBREs2bNiIqKYsCAAXed6nSEDPWjVkr5A58DV4BlWutbFiIrpboCXQFKly5d4+aTTaT38q3k70QI6+3duxd/f3+WLl1KpUqVGDduHM8//7zD68hSP2ql1P1AC6AMUBK4TynV4ebnaa0na619tNY+xYsXz2rNQghhVxcuXKBPnz488cQTrF+/njFjxrBt2zZLQvpuMrLqozFwWGsdB6CUCgfqALPsWZizk37UQmRPWmvmzJlDYGAgJ06c4N1332XEiBE89NBDVpeWrowE9TGgtlLKDTP10QiIuPNLcj7pRy1E9rN9+3b8/Pz4448/qFGjBvPnz6d27dpWl3VXd5360FpvAOYBW4Cdaa+ZbOe6hBDCZs6ePYuvry9PPfUUe/bsYfLkyWzYsCFbhDRkcMOL1joICLJzLUIIYVMpKSlMnz6dfv36ce7cOT744AOGDBnCAw88YHVpmeKcW8iDg0GpWx82XEAuhMjZ/vrrL2rVqkXXrl2pVKkSW7ZsYfz48dkupMGZg1prqF/fPLQ2DwlqIcRdnDp1ik6dOvHMM88QGxvL7NmzWb16NU8++aTVpd0z5wzqaxITYds2uMPmEiGEANM8aezYsZQvX57Zs2fTp08foqKiaNeuXbZv4+DcQX30KMTHw5AhNrlcy5YtqVGjBpUrV2byZLkfKkRO8b///Y/q1avz0UcfUadOHXbu3MmIESMoXLiw1aXZhHMGtaurmZO+duT6pEnm4yyeQzZ9+nQ2b95MREQEISEhnDlzxgbFCiGscuzYMd544w0aNWrE5cuXWbBgAYsXL6ZChQpWl2ZTzhnUhw5Bu3Zwrdm+mxu0bw93aZ50NyEhITz55JPUrl2b48ePs3//fhsUK4RwtISEBD7//HMqVqzIr7/+ypAhQ9i9ezevvvpqtp/muB3nPIXcwwPc3SE11YR1QoL5+OGH7/mSq1at4vfff2f9+vW4ubnRoEGD/7Q3FUJkDwsXLsTf359Dhw7RunVrRo0axaOPPmp1WXblnEENcOoUlCxpQrtWrRvTIPcoPj6e+++/Hzc3N6Kiovjrr79sVKgQwhH2799PQEAAixcvpmLFiixfvpzGjRtbXZZDOG9Qh4dDgwbmzxMmZPlyL7/8MqGhoXh7e1OhQoVssyNJiNzu4sWLDBs2jFGjRlGgQAG++uor/Pz8yJ8/v9WlOYzzBrWNFShQgCVLllhdhhAig7TW/PDDD/Tu3ZuYmBg6duzIiBEj7nheaU7lnDcTr+1MXL3aPGRnohC5ys6dO2nYsCFvvfUWJUqU4M8//2TGjBm5MqTBmYP62m7Efz8kqIXI0c6fP0/Pnj2pXr06O3fuJDQ0lE2bNlGnTh2rS7NUrpn6EEI4r9TUVL755hv69evHmTNn6NatG5999hnFihWzujSn4JwjaiFErrFx40Zq165Nly5dKF++PBEREUycOFFC+l8kqIUQljh9+jTvvfcetWrV4vjx48ycOZM//viD6tWrW12a08nImYkVlFLb/vX4RykVYM+igoODUUrd8giWOWohsr3k5GRCQkIoX7483377Lb1792bv3r106NAhR+4qtIWMnPCyV2tdTWtdDagBXAZ+smdRwcHBaK2pX78+9evXR2uN1jpLQX3kyBGqVKliuyKFEJm2atUqqlevjr+/PzVr1mTnzp2MHDkSd3d3q0tzapmd+mgEHNRaH7VHMTdLTExk27ZtnJQ2p0JkC+n9NtyrVy/atm1Lw4YNuXDhAuHh4SxdupSKFStaXXK2kNmgbgvMud0XlFJdlVIRSqmIuLi4rFcGHD16lPj4eIbYqM1pcnIy7du3x9vbmzZt2nD58mWbXFcIYdz823BCQgLDhg0jLCyMBQsWEBQURGRkJK1atZJpjkzIcFArpfIDrwI/3u7rWuvJWmsfrbVP8eLFs1SUq6srSili0/p7TJo0CaUUrllsc7p371569OhBZGQk7u7uTJw4MUvXE0Kk78yZM1SpUoX+/fvz0ksvERkZSXBwcJb/HedGmRlRNwG2aK1P2auYaw4dOkS7du3Ik9bm1M3Njfbt23M4i21OH3nkEerWrQtAhw4dWLt2bZZrFUL814EDB9i5cye7du0ib968LF26lPDwcLy8vKwuLdvKzIaXt0hn2sPWPDw8cHd3JzU1lTx58pCQkIC7uzsPZ6HNKXDLr1ryq5cQtnPp0iWGDRvGV199RUpKCmXLlmXHjh25qnmSvWRoRK2Uug94AQi3bzk3nDp1ipIlS1K9enW6d+9ukxuKx44dY/369QB899131KtXL8vXFCK301rz448/UrFiRYYNG8Ybb7zBk08+yZkzZzh79qzV5eUIGQpqrfUlrXUxrXW8vQu6Jjw8nHLlylGoUCEmTJhAeHjWf0ZUqFCBCRMm4O3tzblz5/jggw9sUKkQudfu3btp1KgRb7zxBg8++CB//PEHM2fOJDY21qYLAXK7XNPrw8vLi6ioKKvLECJHiI+PJzg4mHHjxl2/Md+1a1cKFSr0n5OTJk2axKRJkyhYsCBXrlyxsOLszSm3kF9bi7l69WpWr14tOxOFcBLXmieVL1+er7/+mi5durBv3z4++OAD8ubNa7eFALmdU46og4ODJZSFcDIRERH4+vqyYcMGnnnmGZYsWcJTTz31n+fYayFAbufQoNZay0qLNFprq0sQIkPi4uLo378/06ZNo0SJEsyYMYMOHTpcHzXf7NpCAA8PD2rVqnV9P4S4dw4L6oIFC3LmzBmKFSuW68Naa82ZM2coWLCg1aUIka7k5GRCQ0MZOHAgFy9e5KOPPiIoKOiufTnCw8NpkHbe6QQbnHcqHBjUnp6eREdHY6vt5dldwYIF8fT0tLoMIW5rzZo1+Pn5sWPHDho3bkxISAje3t5Wl5VrOSyoXVxcKFOmjKPeTghxD2JiYggMDGTOnDmULl2aefPm8dprr+X634Kt5pSrPoQQjpWYmMgXX3xBhQoVCA8PZ+DAgURGRtK6dWsJaSfglKs+hBCO89tvv9GzZ0/279/Pq6++ypgxYyhbtqzVZYl/kRG1ELnUoUOHaNGiBU2aNAFgyZIlLFiwIEshLXsg7EPZY5mYj4+PjoiIsPl1hRBZd/nyZUaMGMGXX35Jvnz5GDhwIAEBARQoUMDq0nI1pdRmrbXP7b4mUx9C5BJaa8LDw+nVqxfHjh2jXbt2fPnll5QqVcrq0sRdyNSHELnAnj17eOGFF2jTpg1FixZl9erVzJ49W0I6m5CgFiIH++eff/j444958skn2bx5M+PHj2fz5s0899xzVpcmMkGmPoTIgVJTU5k1axaffPIJp0+fpkuXLnz++edk9Zg8YY2MHhxQVCk1TykVpZSKVEo9Y+/ChBD3ZsuWLdSrV4933nkHLy8vNm7cyOTJkyWks7GMTn18Dfymta4IPAlE2q8kIcS9OHPmDN27d8fHx4eDBw/yzTffsG7dOnx8bruQQNia1nDihF0ufdegVkoVAZ4DppladJLW+rxdqhFCZFpKSgqTJk2iXLlyTJ06FX9/f/bt28e7776bboc7YUPnzsHXX0PlylC3LqSk2PwtMvJ/sQwQB3yjlNqqlJqadobifyiluiqlIpRSEdJ4SQjHWLt2LT4+PvTo0YNq1aqxfft2xowZQ5EiRawuLWfTGtavh3ffhZIlISAACheGQYMgNdXmb5eRoM4HPAVM0lpXBy4BfW9+ktZ6stbaR2vtI3NhQthXbGwsb7/9Ns8++yxnzpxh7ty5rFixgsqVK1tdWs72zz8waRJUqwZ16sD8+Sast26FDRugUydwcbH522Zk1Uc0EK213pD28TxuE9RCCPtLSkoiJCSEwYMHk5SUxIABA+jXrx/33XfLL7nCljZvhtBQmDMHLl2C6tUhLAzeesuMpO3sriNqrfVJ4LhSqkLapxoBe+xalRDiFsuWLaNq1aoEBgbSsGFD9uzZw9ChQyWk7eXiRZg6FXx8zGP2bHjzTdi40QR3164mpIODQalbHzbsb5KhXh9KqWrAVCA/cAjopLU+l97zpdeHELZz5MgRevXqxU8//cTjjz/O119/TdOmTa0uK+fascOMlmfOhAsXoEoV6NYNOnSAokXTf13aqTasWnVPb5vlXh9a622ArPERwoGuXLnCl19+yYgRI8iTJw/Dhw/no48+kuZJ9nDlCsydawJ6/XooUADeeAO6d4dnnjEjZAvJzkQhnIzWmp9//plevXpx5MgR2rZty8iRI+XoNnuIjDThPGMGnD8PFSrA6NHQsSMUK2Z1dddJUAvhRKKiovD392fZsmVUqVKFlStXXj8oVthIYqJZrREWBmvWmFUarVub6Y369e999JyYaIL/5El4+GGbliyr4YVwAhcuXCAwMJAnnniCDRs2EBISwtatWyWkbWn/fggMBE9PaN8eYmLgiy8gOtqs5mjQIGtTHEePQnw8DBlis5KvkRG1EBbSWjN79mw++eQTTp48SefOnRk2bBglSpSwurSc4epVWLDALK1bsQLy5oWWLc3ouVEjsMXOTVdXSEi48fGkSeZRsKCZ+7YBGVELYZFt27bx7LPP8vbbb+Pp6clff/3F1KlTJaRt4fBh6N8fHnkEXn/djKaHDoXjx2HePHjhBduENMChQ9Cu3Y3rubmZEfvhw7a5PjKiFsLhzp49y6effkpYWBjFihVj6tSpdOrUSfpyZFVyMixaZEbPS5eaaYxXXjErN156yYym7cHDA9zdzdbxPHnM6Nrd3abz1BLUQjhISkoKU6dOZcCAAZw/fx5fX18GDx5M0TutzRV3Fx1tNqZMnWrmnUuWhIEDoUsXM6J2hFOnzPt6eECtWhAba9PLS1AL4QDr16/H19eXLVu2UL9+fcaNG8cTTzxhdVnZV0qKGTWHhcHChaZJ0ksvwYQJZhSdz8HRFh5+Y8PLhAk2v7wEtRB2EhwczODBg2/5fP369SWk71VsLEyfDlOmmFUWDz0EffrA++9DmTJWV2c3EtRC2MHVq1dxd3encOHCXLx4EU9PT/bs2UOhQoWsLi37SU01KzbCwswKjuRks2Jj5Eho0QLy57e6QruToBbCxlasWIGfnx+RkZE0bdqUuLg43NzcJKQzKy4OvvkGJk+GgwfNTsGAANMMqVw5q6u7ITgY/v2b07W12EFBNmvMJLeZhbCRo0eP0qZNGxo3bkxSUhK//vorixYtws3NzerSsg+tYfVq0z60VCkzrVGqlOlcFx1tRtHOFNJgwljrWx827J4nI2ohsighIYGRI0cyfPhwAIYOHcrHH39MwYIFAUhMTCQyMpKTJ0/ysI23FucYZ8+afhuTJ0NUlOlS16OHGT1XqmR1dZaTEbUQ90hrzS+//EKlSpUYNGgQzZo1IyoqigEDBlwPaTAj7fj4eIbYYWtxtqY1rFtnGiCVLAm9epmA/r//M8vsxo6VkE4jQS3EPdi3bx9NmzalRYsWuLm5sWLFCubOnUvp0qWvP8fV1RWlFLFpa2onTZqEUgpXV1eryra/jDTRj483S9iqVjWHwf78M3TuDNu2mRaj77xjdveJG7TWNn/UqFFDC5ETXbhwQffp00e7uLhod3d3PWbMGJ2UlHTb5544cUK3a9dO58mTRwPazc1Nt2/fXsfGxjq4agvUr28e16Smar1xo9bvvae1m5uZxa1RQ+spU7S+cMGqKp0KEKHTydQMzVErpY4AF4AUIFmncwqBEDmV1po5c+YQGBjIiRMn6NSpE8OHD+ehhx5K9zUeHh64u7uTmppKnjx5SEhIwN3dPXfNU1+4YDrThYaaA2Dd3ExfjG7dzPFWIkMyczOxodb6b7tVIoST2rFjB35+fqxZs4YaNWowf/58ateunaHXnjp1ipIlS+Lh4UGtWrWuT4PkeOfOmd7MDz8Mly+baY4JE0yzoiJFrK4u25FVH0Kk49y5cwwaNIiJEydy//33M3nyZDp37kzeTDT3CQ8Pv95TeoIdthY7lcuX4YcfzOh5xw7zubJlzc3BWrUsP84qO8toUGtgmVJKA2Fa68k3P0Ep1RXoCvznhooQ2U1KSgrTp0+nf//+nD17lg8++IAhQ4bwwAMPWF2ac9q92+wa/PZbc6Pw3/buNWcO2rA3c26U0VUf9bTWTwFNgA+VUs/d/ASt9WSttY/W2qd48eI2LVIIR9mwYQO1a9ema9eueHt7s2XLFsaPHy8hfbOEBJg1C5591pzSHRZmmiGFh5vNKnbszZwbZSiotdYxaf89DfwE1LRnUUI42qlTp+jUqRO1a9fmxIkTzJ49m9WrV/Pkk09aXZpz2bsXPv7Y7BZ8+21zPuDIkWbd8+zZ0KqVmYO2Y2/m3OiuQa2Uuk8pVfjan4EXgV32LkwIR7h69Spjx46lfPnyzJ49mz59+hAVFUW7du1QWZxTDQ4ORinF6tWrWb16NUoplFIE23BrsUMkJZm55+efh4oVISTENEX6/XcT3L17w4MP3nj+td7M1aubpv0nT1pXew6hzPK9OzxBqbKYUTSYOe3vtNaf3+k1Pj4+OiIiwjYVCmEn//vf/+jZsye7d+/m5ZdfZuzYsVSoUMHqspzHoUNmS/c338Dp0+DlZbZ0d+p09xHytd7Mq1bZucicQym1Ob2lz3e9mai1PgTI738ixzh27Bi9e/fmxx9/pEyZMixYsIDmzZtneQSdI1y9Cr/+auacly0zx1c1b27WPb/4ou3OGRSZIsvzRK6RkJDAqFGj+Pzzz9FaM2TIEHr37p2zt3Rn1LFjphn/tGmmOb+np2nd+d57Zj5aWEqCWuQKCxcuxN/fn0OHDtG6dWtGjRrFo48+anVZ1kpJgSVLzLrnJUtMk6QmTcxoukmTezvOygG9mXOju85R3wuZoxbOYv/+/QQEBLB48WK8vb0JCQmhcePGVpdlrZgYM3KeOhWOHzfzzV26mEdu/+FloSzNUQuRHV28eJFhw4YxatQoChQowKhRo/Dz88PFxcXq0qyRmgrLl5vR86+/mtH0Cy+YVqLNm0Nu/XvJJiSoRY6iteaHH36gd+/exMTE0LFjR0aMGIGHh4fVpVnj1Kkbx1kdPgzFi5vldO+/D489ZnV1IoMkqEWOsXPnTvz8/Fi9ejXVq1dn7ty51KlTx+qyHE9rWLnSjJ5//tms5GjQAIYPh5YtoUABiwsUmSVBLbK98+fPX2+eVKRIEUJDQ+nSpUummiflCH//bY6zCguD/fvhgQfA19esfa5Y0erqRBbIokiRbaWmpjJt2jTKly/PhAkT6Nq1K/v27aNbt27OEdIZOe0kq7SGP/4w/TRKlTLTGiVKwMyZ5qbh6NES0jmArPoQ2dLGjRvx9fVl06ZN1K1bl3HjxlG9enWry7q9Z54xvZmjomzX8+LcORPGYWGwZ4/pr/H222ZjSpUqtnkP4VB3WvUhI2qRrZw+fZr33nuPWrVqcfz4cWbOnMkff/zhvCENcPSoaf+Z1cNttYa//oJ33zW9NPz9oVAhs9QuJgbGjZOQzqFkRC2yheTkZCZOnMigQYO4dOkSAQEBDBw4EHd3d6tLS5+rq+ked7PM9mb+5x/Tme5aQ/5ChcxUR7dupvGRyBFkRC2ytVWrVlG9enX8/f2pWbMmO3fuZOTIkc4d0mCaGrVrd++9mTdvNjcCS5aEHj1M343QUDhxwvxXQjrXkKAWDnOt7efNj/TafkZHR9O2bVsaNmzIhQsXCA8PZ+nSpVTMLjfHPDxML+bM9Ga+eNHsGPTxMY9Zs+DNN2HDBhPc3bpB4cKO+x6EU5CpD+Fw184QXJVOC8zExERGjx7N0KFDSU1NpU+fPvTp0yd7Nk967TUTsh4e5tzA2FhzCsrNduwwNwZnzjQnd1epYkK5QwcoWtThZQvHky3kwqkkJiYSGRnJyZMnefim0eWiRYsICAjgwIEDtGrVitGjR+Pl5WVNobYQHn6jN/PNh9teuQJz55qAXr/ebER54w0T0HXqyGGw4roMT30opfIqpbYqpRbasyCR8x09epT4+HiG/GsVxIEDB2jevDnNmjUjT548/Pbbb4SHh2fvkE5PZCQEBJi553ffhTNnzHrnmBhzQGzduhLS4j8yPPWhlOoF+ADuWutmd3quTH2I23F1dSXhNqsg8ubNS968ecmfPz+DBg3C39+f/PnzW1Chjd3c8vPfXFygdWszeq5fX4JZZH3Vh1LKE3gFmGrLwkTucujQIdq1a0eetFUQ+fPnx9XVlZSUFF5//XX27t1LYGBgzghpMEG9b99/zxR87DH44guIjoY5c8y0iIS0uIuMzlGPBT4B0r3drJTqCnQFKF26dJYLEzmPh4cH7u7upKamApCUlESxYsVYtmwZ9erVs7g6G0pKggULzNzzihVmWV2LFuag10aN5DgrkWl3DWqlVDPgtNZ6s1KqQXrP01pPBiaDmfqwVYEi5zh//jzLly8HzHRHnTp1KFasWM4J6cOHzXFW06eb9qKlS8PQodC5s1n1IcQ9ysiIui7wqlKqKVAQcFdKzdJad7BvaSKnSE1NZcaMGfTt25e4uDg8PDwoU6YMa9assbq0rEtOhoULzeh56VIzjdGsmZl7fuklM5oWIovu+juY1rqf1tpTa+0FtAX+JyEtMmrTpk3UqVOHzp0789hjj7Fp0ybKly+f/U9aOX7cnAPo5QWtWpl10IMGwZEjZtqjaVMJaWEzso5a2EVcXBz9+/dn2rRplChRghkzZnDw4EF8fG7c1FZpN9GCgoLS3Z3oVFJSzKg5NBQWLTJNkl56yayPfuWVezsMVogMkJ2JwqaSk5MJDQ1l4MCBXLx4kZ49ezJo0CCKFClidWn3LjbWdKibMgWOHYOHHoL33jOHwZYpY3V1IoeQnYnCIdasWYOvry87d+6kUaNGhISEUKlSJavLujepqWbFRmgo/PKLmYtu1AhGjYJXX4WcsoRQZAsS1CLLYmJiCAwMZM6cOZQuXZp58+bx2muvXZ/ayFZOn4b/+z9zGOzBg1CsmNlF2LUrlCtndXUil5KgFvcsMTGRsWPH8tlnn5GcnMzAgQPp27cvbm5uVpeWOVrD6tVm5cb8+eYw2OeeM43+W7eWw2CF5SSoxT357bff6NmzJ/v37+fVV19lzJgxlC1b1uqyMufs2RuHwe7da7rU9ehhRs/ZdcpG5EgS1CJTDh06xEcffcQvv/xCuXLlWLJkCS+//LLVZWWc1rBunQnnuXMhMdGcafh//2c612XHVqoix5OgFhly+fJlRowYwZdffkm+fPkYMWIEAQEBFMgu0wLx8TcOg921yzTff+89szGlalWrqxPijiSoxR1prZk/fz4ff/wxx44do127dnz55ZeUKlXK6tLuTmvYtMmE8/ffw+XL5tSUKVOgbVtz9qAQ2YAEtUjXnj176NmzJytWrKBq1arMnDmT5557zuqy7u7CBfjuOxPQW7fCfffdOAy2Rg2rqxMi0ySoxS3i4+MZPHgw48aNo1ChQowfP55u3bqRz9l33m3dasJ59mxz9mDVqjBxoglpZz8IV4g7cPJ/ecKRUlNTmTlzJn369OH06dN06dKFzz//nOLFi1tdWvouXYIffjABvXEjFCxopjW6dTNnFGbHtdxC3ESCWgCwZcsWfH19Wb9+PbVq1WLhwoX/6cvhdHbtunEYbHw8eHvD11/D22/D/fdbXZ0QNiVBncv9/fffDBgwgClTplC8eHG++eYbOnbseP0UFqeSkADz5plt3X/+abZxt2ljGvLXqyejZ5FjSVDnUikpKYSFhfHpp5/yzz//4O/vT3BwsHM2T9q714yeZ8wwm1TKlYOvvoJ33rlxxJUQOZgEdS60du1afH192b59Ow0bNmTcuHFUrlzZ6rL+KykJfvrJjJ5XrTItRFu1MqPnhg1l9CxyFSf8/VbYy4kTJ+jQoQPPPvssZ8+eZe7cuaxYscJxIR0cbAL25se/e1EfPAh9+4Knp7kpeOQIDBtmGvXPnQvPPy8hLXIfrfUdH5jjtzYC24HdwOC7vaZGjRpaOI/ExET95Zdf6kKFCun8+fPrAQMG6IsXL1pXUO3aWhcponVsrPk4KUnr+fO1fuEFrUHrvHm1btlS699+0zolxbo6hXAgIEKnk6kZmfpIBJ7XWl9USrkAa5VSS7TWf9nrh4ewnaVLl+Lv78/evXtp3rw5Y8aM4bHHHrO2qKNHzUqNwEBzlNW0aaY5v6cnDB5stnZnh52PQjjIXYM6Lekvpn3okvaQU8ad3OHDh+nVqxc///wzjz/+OIsWLaJp06bWFuXqalZuXDNrlvlvnjymOX+TJnKclRC3kaE5aqVUXqXUNuA0sFxrveE2z+mqlIpQSkXExcXZuEyRUZcvXyYoKIhKlSqxbNkyhg8fzq5du6wP6ZgY8PX9b3e6fPmgRQvztebNJaSFSEeGglprnaK1rgZ4AjWVUlVu85zJWmsfrbWPU+9ky6G01oSHh1OpUiWGDBlCy5Yt2bt3L3379rWuw11qKvz2m1mt8eijZkldsWLma0qZr5csCQ8/bE19QmQTmVr1obU+D6wEslED4pwvKiqKl156idatW1O4cGFWrlzJnDlz8PT0tKagU6dg+HB47DEznfHnn9C7Nxw4AE8/bcL5qafMUruTJ62pUYjsJL27jPrGqo/iQNG0P7sCfwDN7vQaWfXhGPHx8frjjz/W+fLl00WKFNEhISH66tWr1hSTkqL1779r/frrWufLZ1ZvNGyo9fffa52QYJ4TFGQ+f/MjKMiamoVwItxh1YcyX0+fUqoqMAPIixmBz9VaD7nTa3x8fHRERESWf4iI29NaM2vWLD755BNOnTpF586dGTZsGCVKlHB8MX//feMw2P374YEH4N13zXFWFSo4vh4hsiml1Gat9W0b7GRk1ccOoLrNqxL3ZOvWrfj6+rJu3TqefvppFixYQM2aNR1bhNawdq3ZNThvntlFWLcuDBpkem8ULOjYeoTI4eQ2ezZx5swZPv30U8LCwnjwwQeZOnUqnTp1cmzzpHPn4NtvTd+NyEgoUsS0E+3aFarccn9ZCGEjEtROLiUlhSlTpjBgwADi4+Px8/Nj8ODBFC1a1DEFaA1//WXC+YcfzDromjVh+nR4801wc3NMHULkYhLUTmzdunX4+vqydetW6tevz7hx43jiiScc8+b//GNOSgkNhR07zPmC77xjRtDVZSZMCEeSpkxOKDY2lo4dO1K3bl1Onz7N999/z8qVKx0T0ps3w/vvmyV0PXpA3rxmNH3ihAltCWkhHE5G1BYJDg5m8ODBt3z+hRde4K+//iIxMZF+/frRv39/Ctn7tOyLF2HOHBPImzeb3YNvvWVGz08/Ld3qhLDYXZfn3QtZnpdxJUuWJDY2lmbNmnHgwAGioqJo2rQpY8eOpVy5cvZ98x07zCh51ixzcneVKiacO3QAR82BCyGAOy/Pk6C2iKurKwn/blCUxsXFhaSkJPu98ZUrpq9zaKi5SVigALzxhgnoOnVk9CyERe4U1DJHbZHdu3f/Z87ZxcWFtm3bcuzYMfu8YWQkBASYued33zVHWo0ebRoiffutWQctIS2EU5I5agfTWvPLL78QEBDAkSNHAFBKkZKSwv3338/DtmxQlJgI8+ebuec1a8DFBVq3NqPn+vUlmIXIJmRE7UB79+6lSZMmtGzZkkuXLl3/vNaa1NRUJk2aRPC/j6W6V/v3m6b8pUpB+/Zm1PzFFxAdbW4aNmggIS1ENiIjage4cOECQ4cOZcyYMbi6ujJmzBg+/PBDXFxcbPcmSUmwYIEZPa9YYZbVtWxpRs+NGpnm/EKIbEmC2o601nz33XcEBgYSGxtLp06dGD58OA899JDt3uTwYZgyxewUPHUKSpeGoUOhc2fw8LDd+wghLCNBbSfbt2/Hz8+PP/74Ax8fH8LDw6ldu7ZtLp6cDAsXmtHz0qVmGuOVV0x/55deMqNpIUSOIUFtY2fPnmXgwIGEhobywAMPMGXKFDp37myb5knHj8PUqeYw2JgYs4Jj4EDo0gUeeSTr1xdCOCUJahtJSUlh2rRp9O/fn3PnztGjRw+GDBnC/fffn9ULm+OswsJg0SLTJOmll2DCBDOKlnMGhcjx5A6TDaxfv55atWrRrVs3KlWqxJYtWxg3blzGQjo21iyVu/lIqthYM9dctiw0awYbN0KfPnDwICxZYg6FlZAWIle46790pdQjwLfAQ4AGJmutv7Z3YdnByZMn6du3LzNmzKBkyZJ89913tG3bFpWRpW/BwfDvXh/Xbvx16ACXL8Mvv5i56EaNzKGwLVpA/vx2+T6EEM4tI0dxeQAeWustSqnCwGagpdZ6T3qvyelbyK9evcr48eMJDg7mypUr9OrVi08//TRzzZNcXU1v59spVgw6dTIN+e3d70MI4RSytIVcax2rtd6S9ucLQCRQyrYlZh8rVqygWrVq9OrVizp16rBr1y5GjBiR+Q53Bw+a0fK/lShh5p6jo2HkSAlpIQSQyTlqpZQX5vzEDbf5WlelVIRSKiIuLs5G5TmPY8eO8frrr9O4cWOuXLnCggULWLx4MeXLl8/chc6ehTFj4PnnzcaUa5Qy27t79JAzB4UQ/5Hhu1FKqULAfCBAa/3PzV/XWk8GJoOZ+rBZhRZLSEjgq6++YtiwYQAMGTKEwMBACmYmTLWGdetMx7offzQ9OJ55Bp56yjTkL1UKatUyNxCFEOJmWuu7PgAXYCnQKyPPr1Gjhs7uUlNT9YIFC3TZsmU1oFu3bq2PHDmSuYucO6f1uHFaV6miNWhduLDWPXpovX271kFB5nM3P4KC7PDdCCGcHRCh08nUjKz6UMA0IFJrPdqePzScxb59+wgICGDJkiV4e3uzfPlyGjdunLEXaw2bNpl1z3PmmP7PNWqYbd5t25qzBwGqVjUrP4QQ4i4yMvVRF3gb2KmU2pb2uf5a68V2q8oiFy9eZOjQoYwePZqCBQsyatQo/Pz8MtY86cIF+O47E9Bbt8J995mldt26maAWQoh7dNeg1lqvBXJ0T0ytNd9//z2BgYHExMTQsWNHvvjii4z1ht661YTz7Nnm7MGqVWHiRNNe1N3d/sULIXK8XL+1bceOHfj5+bFmzRqeeuop5s6dS506de78okuX4IcfTEBv3GhWabz5pmmKVKuW9HoWQthUrg3qc+fOERQUxIQJEyhatCihoaF06dKFvHfqPLdrlwnnmTMhPh68vWHsWOjYEbLa00MIIdKR64I6NTWV6dOn069fP86ePUv37t357LPPeOCBB27/goQEmDfPLK3780+zjbtNGzP3/OyzMnoWQthdrgrqjRs34uvry6ZNm6hbty7jx4+nWrVqt3/y3r1m9DxjhtmkUq6c2S347rvw4IOOLFsIkcvliqA+ffo0/fr1Y/r06Xh4eDBr1izatWt3a/OkpCT46Sczel61ynSna9XKzD03aCDHWQkhLJGjgzo5OZkJEyYQFBTEpUuXCAwMZODAgRQuXPi/Tzx48MZxVnFx4OUFw4aZxki2PBVcCCHuQY4N6lWrVuHn58euXbt48cUX+frrr6lYseKNJ1y9Cr/+akbPy5eb46uaNzdzzy++KKNnIYTTyHFBffz4cXr37s3cuXPx8vLip59+okWLFjemOY4evXGcVWwseHqavtDvvWd6bgghhJPJMUGdmJjIqFGj+Pzzz0lNTSU4OJhPPvkEV1dXc5zV4sXm5uDitA2VTZua0XOTJnJSihDCqeWIhFq0aBH+/v4cPHiQVq1aMXr0aLy8vMwBsNOmmRH08eNmvnnAAHMY7KOPWl22EEJkSLaeiD1w4ADNmjWjWbNmuLi4sGzZMsLnzcMrKsqs1nj0UQgKgooVYf58OHYMPvtMQloIka1ky6C+dOkSAwYMoHLlyqxevZqRI0eyfdkyXoiIgMcfN9MZf/4JvXvDgQOwbBm89hpkpLmSEEI4mWw19aG1Zu7cufTu3Zvo6Gg6tG/PF82bU3L+fOjXzxwG26ABDB8OLVtCgQJWlyyEEFmWbYJ6165d+Pn5sWrVKqpVqcKcV1+l3vLlpmvdAw+An585DPbfS/CEECIHcPqgPn/+PMHBwYwfPx53NzcmPv00XbdtI++uXVC3LgwaZHpvyDmDQogc6q5z1Eqp6Uqp00qpXXavJjYW6teHkyevN08qX64cISEhdClcmH0XLvDBvn3k7dYNdu6EtWtNc34JaSFEDpaREfX/AeOBb+1bCmZFxtq1nPrwQ16NimLjnj08oxS/ac1T5cubdc9vvmlOTxFCiFwiIye8rFFKedm1ChcXSE4mFmgL/BAezgYgGcjz/vvk6d4dqle3awlCCOGsbLY8TynVVSkVoZSKiIuLy9yLjx0jsUULBgFrgSAgqUYN8u3fT56wMAlpIUSupswp5Xd5khlRL9RaV8nIRX18fHRERESGi3DJm5fk1NRbPp8vTx6upqRk+DpCCJFdKaU2a619bvc1p9jwciw6mnaenuRP+9gtXz7ae3pyPCbG0rqEEMIZOEVQh4WF8V10NElpH19OTmZ2dDShoaGW1iWEEM7grjcTlVJzgAbAg0qpaCBIaz3NlkUEBwezY8cOPDw86Nq1K5MnTyY2Npbg4GBbvo0QQmRLGZqjzqzMzlELIURu5/Rz1EIIIdInQS2EEE5OgloIIZycBLUQQjg5CWohhHByEtRCCOHkJKiFEMLJ2WUdtVIqDjh6jy9/EPjbhuVkB/I953y57fsF+Z4z61GtdfHbfcEuQZ0VSqmI9BZ951TyPed8ue37BfmebUmmPoQQwslJUAshhJNzxqCebHUBFpDvOefLbd8vyPdsM043Ry2EEOK/nHFELYQQ4l8kqIUQwsk5TVArpV5WSu1VSh1QSvW1uh5HUEpNV0qdVkrtsroWR1BKPaKUWqmU2qOU2q2U8re6JntTShVUSm1USm1P+54HW12Toyil8iqltiqlFlpdiyMopY4opXYqpbYppWzakN8p5qiVUnmBfcALQDSwCXhLa73H0sLsTCn1HHAR+DajBwdnZ0opD8BDa71FKVUY2Ay0zMn/n5VSCrhPa31RKeUCrAX8tdZ/WVya3SmlegE+gLvWupnV9dibUuoI4KO1tvkmH2cZUdcEDmitD2mtk4DvgRYW12R3Wus1wFmr63AUrXWs1npL2p8vAJFAKWursi9tXEz70CXtYf3oyM6UUp7AK8BUq2vJCZwlqEsBx//1cTQ5/B9wbqeU8gKqAxssLsXu0qYAtgGngeVa6xz/PQNjgU+AVIvrcCQNLFNKbVZKdbXlhZ0lqEUuopQqBMwHArTW/1hdj71prVO01tUAT6CmUipHT3MppZoBp7XWm62uxcHqaa2fApoAH6ZNbdqEswR1DPDIvz72TPucyGHS5mnnA7O11uFW1+NIWuvzwErgZYtLsbe6wKtpc7bfA88rpWZZW5L9aa1j0v57GvgJM6VrE84S1JuAckqpMkqp/EBb4BeLaxI2lnZjbRoQqbUebXU9jqCUKq6UKpr2Z1fMDfMoS4uyM611P621p9baC/Nv+X9a6w4Wl2VXSqn70m6Qo5S6D3gRsNlqLqcIaq11MuALLMXcYJqrtd5tbVX2p5SaA6wHKiilopVS71ldk53VBd7GjLC2pT2aWl2UnXkAK5VSOzADkuVa61yxXC2XeQhYq5TaDmwEFmmtf7PVxZ1ieZ4QQoj0OcWIWgghRPokqIUQwslJUAshhJOToBZCCCcnQS2EEE5OgloIIZycBLUQQji5/weO8BTT/ODOrwAAAABJRU5ErkJggg==\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" + }, + "vscode": { + "interpreter": { + "hash": "31f2aee4e71d21fbe5cf8b01ff0e069b9275f58929596ceb00d14d90e3e16cd6" + } + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/pyerrors/fits.py b/pyerrors/fits.py index d74f5feb..1e73466f 100644 --- a/pyerrors/fits.py +++ b/pyerrors/fits.py @@ -72,9 +72,12 @@ def __repr__(self): def least_squares(x, y, func, priors=None, silent=False, **kwargs): r'''Performs a non-linear fit to y = func(x). + ``` Parameters ---------- + For an uncombined fit: + x : list list of floats. y : list @@ -96,9 +99,32 @@ def func(a, x): (x1, x2) = x return a[0] * x1 ** 2 + a[1] * x2 ``` + It is important that all numpy functions refer to autograd.numpy, otherwise the differentiation + will not work. + + OR For a combined fit: + + x : dict + dict of lists. + y : dict + dict of lists of Obs. + funcs : 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. + priors : list, optional priors has to be a list with an entry for every parameter in the fit. The entries can either be Obs (e.g. results from a previous fit) or strings containing a value and an error formatted like @@ -114,6 +140,11 @@ def func(a, x): The possible methods are the ones which can be used for scipy.optimize.minimize and migrad of iminuit. If no method is specified, Levenberg-Marquard is used. Reliable alternatives are migrad, Powell and Nelder-Mead. + tol: float, optional + can be used (only for combined fits and methods other than Levenberg-Marquard) to set the tolerance for convergence + to a different value to either speed up convergence at the cost of a larger error on the fitted parameters (and possibly + invalid estimates for parameter uncertainties) or smaller values to get more accurate parameter values + The stopping criterion depends on the method, e.g. migrad: edm_max = 0.002 * tol * errordef (EDM criterion: edm < edm_max) correlated_fit : bool If True, use the full inverse covariance matrix in the definition of the chisquare cost function. For details about how the covariance matrix is estimated see `pyerrors.obs.covariance`. @@ -137,6 +168,11 @@ def func(a, x): ''' if priors is not None: return _prior_fit(x, y, func, priors, silent=silent, **kwargs) + + elif (type(x) == dict and type(y) == dict and type(func) == dict): + return _combined_fit(x, y, func, silent=silent, **kwargs) + elif (type(x) == dict or type(y) == dict or type(func) == dict): + raise TypeError("All arguments have to be dictionaries in order to perform a combined fit.") else: return _standard_fit(x, y, func, silent=silent, **kwargs) @@ -474,7 +510,6 @@ def chisqfunc_compact(d): def _standard_fit(x, y, func, silent=False, **kwargs): - output = Fit_result() output.fit_function = func @@ -668,6 +703,180 @@ def chisqfunc_compact(d): return output +def _combined_fit(x, y, func, silent=False, **kwargs): + + if kwargs.get('correlated_fit') is True: + raise Exception("Correlated fit has not been implemented yet") + + output = Fit_result() + output.fit_function = func + + if kwargs.get('num_grad') is True: + jacobian = num_jacobian + hessian = num_hessian + else: + jacobian = auto_jacobian + hessian = auto_hessian + + key_ls = sorted(list(x.keys())) + + if sorted(list(y.keys())) != key_ls: + raise Exception('x and y dictionaries do not contain the same keys.') + + if sorted(list(func.keys())) != key_ls: + raise Exception('x and func dictionaries do not contain the same keys.') + + x_all = np.concatenate([np.array(x[key]) for key in key_ls]) + y_all = np.concatenate([np.array(y[key]) for key in key_ls]) + + y_f = [o.value for o in y_all] + dy_f = [o.dvalue for o in y_all] + + if len(x_all.shape) > 2: + raise Exception('Unknown format for x values') + + # number of fit parameters + n_parms_ls = [] + for key in key_ls: + if not callable(func[key]): + raise TypeError('func (key=' + key + ') is not a function.') + if len(x[key]) != len(y[key]): + raise Exception('x and y input (key=' + key + ') do not have the same length') + for i in range(100): + try: + func[key](np.arange(i), x_all.T[0]) + except TypeError: + continue + except IndexError: + continue + else: + break + else: + raise RuntimeError("Fit function (key=" + key + ") 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): + func_list = np.concatenate([[func[k]] * len(x[k]) for k in key_ls]) + model = anp.array([func_list[i](p, x_all[i]) for i in range(len(x_all))]) + chisq = anp.sum(((y_f - model) / dy_f) ** 2) + return chisq + + output.method = kwargs.get('method', 'Levenberg-Marquardt') + if not silent: + print('Method:', output.method) + + if output.method != 'Levenberg-Marquardt': + if output.method == 'migrad': + tolerance = 1e-4 # default value of 1e-1 set by iminuit can be problematic + if 'tol' in kwargs: + tolerance = kwargs.get('tol') + fit_result = iminuit.minimize(chisqfunc, x0, tol=tolerance) # Stopping criterion 0.002 * tol * errordef + output.iterations = fit_result.nfev + else: + tolerance = 1e-12 + if 'tol' in kwargs: + tolerance = kwargs.get('tol') + fit_result = scipy.optimize.minimize(chisqfunc, x0, method=kwargs.get('method'), tol=tolerance) + output.iterations = fit_result.nit + + chisquare = fit_result.fun + + else: + def chisqfunc_residuals(p): + model = np.concatenate([np.array(func[key](p, np.asarray(x[key]))) for key in key_ls]) + chisq = ((y_f - model) / dy_f) + return chisq + if 'tol' in kwargs: + print('tol cannot be set for Levenberg-Marquardt') + fit_result = scipy.optimize.least_squares(chisqfunc_residuals, x0, method='lm', ftol=1e-15, gtol=1e-15, xtol=1e-15) + + chisquare = np.sum(fit_result.fun ** 2) + assert np.isclose(chisquare, chisqfunc(fit_result.x), atol=1e-14) + output.iterations = fit_result.nfev + + output.message = fit_result.message + + if not fit_result.success: + raise Exception('The minimization procedure did not converge.') + + if x_all.shape[-1] - n_parms > 0: + output.chisquare = chisquare + output.dof = x_all.shape[-1] - n_parms + output.chisquare_by_dof = output.chisquare / output.dof + output.p_value = 1 - scipy.stats.chi2.cdf(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) + + def chisqfunc_compact(d): + func_list = np.concatenate([[func[k]] * len(x[k]) for k in key_ls]) + model = anp.array([func_list[i](d[:n_parms], x_all[i]) for i in range(len(x_all))]) + chisq = anp.sum(((d[n_parms:] - model) / dy_f) ** 2) + return chisq + + def prepare_hat_matrix(): + hat_vector = [] + for key in key_ls: + x_array = np.asarray(x[key]) + if (len(x_array) != 0): + hat_vector.append(jacobian(func[key])(fit_result.x, x_array)) + hat_vector = [item for sublist in hat_vector for item in sublist] + return hat_vector + + fitp = fit_result.x + + if np.any(np.asarray(dy_f) <= 0.0): + raise Exception('No y errors available, run the gamma method first.') + + 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.") + + if kwargs.get('expected_chisquare') is True: + if kwargs.get('correlated_fit') is not True: + W = np.diag(1 / np.asarray(dy_f)) + cov = covariance(y_all) + hat_vector = prepare_hat_matrix() + A = W @ hat_vector # hat_vector = 'jacobian(func)(fit_result.x, x)' + P_phi = A @ np.linalg.pinv(A.T @ A) @ A.T + expected_chisquare = np.trace((np.identity(x_all.shape[-1]) - P_phi) @ W @ cov @ W) + output.chisquare_by_expected_chisquare = output.chisquare / expected_chisquare + if not silent: + print('chisquare/expected_chisquare:', output.chisquare_by_expected_chisquare) + + result = [] + for i in range(n_parms): + result.append(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 + + def fit_lin(x, y, **kwargs): """Performs a linear fit to y = n + m * x and returns two Obs n, m. diff --git a/tests/fits_test.py b/tests/fits_test.py index 828c0cbe..8a9a4f33 100644 --- a/tests/fits_test.py +++ b/tests/fits_test.py @@ -108,24 +108,6 @@ def test_prior_fit_num_grad(): auto = pe.fits.least_squares(x, y, lambda a, x: anp.exp(a[0] * x) + a[1], num_grad=False, piors=y[:2]) -def test_least_squares_num_grad(): - x = [] - y = [] - for i in range(2, 5): - x.append(i * 0.01) - y.append(pe.pseudo_Obs(i * 0.01, 0.0001, "ens")) - - num = pe.fits.least_squares(x, y, lambda a, x: np.exp(a[0] * x) + a[1], num_grad=True) - auto = pe.fits.least_squares(x, y, lambda a, x: anp.exp(a[0] * x) + a[1], num_grad=False) - - assert(num[0] == auto[0]) - assert(num[1] == auto[1]) - - - assert(num[0] == auto[0]) - assert(num[1] == auto[1]) - - def test_total_least_squares_num_grad(): x = [] y = [] @@ -608,6 +590,249 @@ def f(a, x): pe.fits.ks_test(fit_res) +def test_combined_fit_list_v_array(): + res = [] + for y_test in [{'a': [pe.Obs([np.random.normal(i, 0.5, 1000)], ['ensemble1']) for i in range(1, 7)]}, + {'a': np.array([pe.Obs([np.random.normal(i, 0.5, 1000)], ['ensemble1']) for i in range(1, 7)])}]: + for x_test in [{'a': [0, 1, 2, 3, 4, 5]}, {'a': np.arange(6)}]: + for key in y_test.keys(): + [item.gamma_method() for item in y_test[key]] + def func_a(a, x): + return a[1] * x + a[0] + + funcs_test = {"a": func_a} + res.append(pe.fits.least_squares(x_test, y_test, funcs_test)) + + assert (res[0][0] - res[1][0]).is_zero(atol=1e-8) + assert (res[0][1] - res[1][1]).is_zero(atol=1e-8) + + +def test_combined_fit_vs_standard_fit(): + + x_const = {'a':[0, 1, 2, 3, 4, 5, 6, 7, 8, 9], 'b':np.arange(10, 20)} + y_const = {'a':[pe.Obs([np.random.normal(1, val, 1000)], ['ensemble1']) + for val in [0.25, 0.3, 0.01, 0.2, 0.5, 1.3, 0.26, 0.4, 0.1, 1.0]], + 'b':[pe.Obs([np.random.normal(1, val, 1000)], ['ensemble1']) + for val in [0.5, 1.12, 0.26, 0.25, 0.3, 0.01, 0.2, 1.0, 0.38, 0.1]]} + for key in y_const.keys(): + [item.gamma_method() for item in y_const[key]] + y_const_ls = np.concatenate([np.array(o) for o in y_const.values()]) + x_const_ls = np.arange(0, 20) + + def func_const(a,x): + return 0 * x + a[0] + + funcs_const = {"a": func_const,"b": func_const} + for method_kw in ['Levenberg-Marquardt', 'migrad', 'Powell', 'Nelder-Mead']: + res = [] + res.append(pe.fits.least_squares(x_const, y_const, funcs_const, method = method_kw, expected_chisquare=True)) + res.append(pe.fits.least_squares(x_const_ls, y_const_ls, func_const, method = method_kw, expected_chisquare=True)) + [item.gamma_method for item in res] + assert np.isclose(0.0, (res[0].chisquare_by_dof - res[1].chisquare_by_dof), 1e-14, 1e-8) + assert np.isclose(0.0, (res[0].chisquare_by_expected_chisquare - res[1].chisquare_by_expected_chisquare), 1e-14, 1e-8) + assert np.isclose(0.0, (res[0].p_value - res[1].p_value), 1e-14, 1e-8) + assert (res[0][0] - res[1][0]).is_zero(atol=1e-8) + +def test_combined_fit_no_autograd(): + + def func_exp1(x): + return 0.3*np.exp(0.5*x) + + def func_exp2(x): + return 0.3*np.exp(0.8*x) + + xvals_b = np.arange(0,6) + xvals_a = np.arange(0,8) + + def func_a(a,x): + return a[0]*np.exp(a[1]*x) + + def func_b(a,x): + return a[0]*np.exp(a[2]*x) + + funcs = {'a':func_a, 'b':func_b} + xs = {'a':xvals_a, 'b':xvals_b} + ys = {'a':[pe.Obs([np.random.normal(item, item*1.5, 1000)],['ensemble1']) for item in func_exp1(xvals_a)], + 'b':[pe.Obs([np.random.normal(item, item*1.4, 1000)],['ensemble1']) for item in func_exp2(xvals_b)]} + + for key in funcs.keys(): + [item.gamma_method() for item in ys[key]] + + with pytest.raises(Exception): + pe.least_squares(xs, ys, funcs) + + pe.least_squares(xs, ys, funcs, num_grad=True) + +def test_combined_fit_invalid_fit_functions(): + def func1(a, x): + return a[0] + a[1] * x + a[2] * anp.sinh(x) + a[199] + + def func2(a, x, y): + return a[0] + a[1] * x + + def func3(x): + return x + + def func_valid(a,x): + return a[0] + a[1] * x + + xvals =[] + yvals =[] + err = 0.1 + + for x in range(1, 8, 2): + xvals.append(x) + yvals.append(pe.pseudo_Obs(x + np.random.normal(0.0, err), err, 'test1') + pe.pseudo_Obs(0, err / 100, 'test2', samples=87)) + [o.gamma_method() for o in yvals] + for func in [func1, func2, func3]: + with pytest.raises(Exception): + pe.least_squares({'a':xvals}, {'a':yvals}, {'a':func}) + with pytest.raises(Exception): + pe.least_squares({'a':xvals, 'b':xvals}, {'a':yvals, 'b':yvals}, {'a':func, 'b':func_valid}) + with pytest.raises(Exception): + pe.least_squares({'a':xvals, 'b':xvals}, {'a':yvals, 'b':yvals}, {'a':func_valid, 'b':func}) + +def test_combined_fit_invalid_input(): + xvals =[] + yvals =[] + err = 0.1 + def func_valid(a,x): + return a[0] + a[1] * x + for x in range(1, 8, 2): + xvals.append(x) + yvals.append(pe.pseudo_Obs(x + np.random.normal(0.0, err), err, 'test1') + pe.pseudo_Obs(0, err / 100, 'test2', samples=87)) + with pytest.raises(Exception): + pe.least_squares({'a':xvals}, {'b':yvals}, {'a':func_valid}) + +def test_combined_fit_no_autograd(): + + def func_exp1(x): + return 0.3*np.exp(0.5*x) + + def func_exp2(x): + return 0.3*np.exp(0.8*x) + + xvals_b = np.arange(0,6) + xvals_a = np.arange(0,8) + + def func_a(a,x): + return a[0]*np.exp(a[1]*x) + + def func_b(a,x): + return a[0]*np.exp(a[2]*x) + + funcs = {'a':func_a, 'b':func_b} + xs = {'a':xvals_a, 'b':xvals_b} + ys = {'a':[pe.Obs([np.random.normal(item, item*1.5, 1000)],['ensemble1']) for item in func_exp1(xvals_a)], + 'b':[pe.Obs([np.random.normal(item, item*1.4, 1000)],['ensemble1']) for item in func_exp2(xvals_b)]} + + for key in funcs.keys(): + [item.gamma_method() for item in ys[key]] + + with pytest.raises(Exception): + pe.least_squares(xs, ys, funcs) + + pe.least_squares(xs, ys, funcs, num_grad=True) + + +def test_combined_fit_num_grad(): + def func_exp1(x): + return 0.3*np.exp(0.5*x) + + def func_exp2(x): + return 0.3*np.exp(0.8*x) + + xvals_b = np.arange(0,6) + xvals_a = np.arange(0,8) + + def func_num_a(a,x): + return a[0]*np.exp(a[1]*x) + + def func_num_b(a,x): + return a[0]*np.exp(a[2]*x) + + def func_auto_a(a,x): + return a[0]*anp.exp(a[1]*x) + + def func_auto_b(a,x): + return a[0]*anp.exp(a[2]*x) + + funcs_num = {'a':func_num_a, 'b':func_num_b} + funcs_auto = {'a':func_auto_a, 'b':func_auto_b} + xs = {'a':xvals_a, 'b':xvals_b} + ys = {'a':[pe.Obs([np.random.normal(item, item*1.5, 1000)],['ensemble1']) for item in func_exp1(xvals_a)], + 'b':[pe.Obs([np.random.normal(item, item*1.4, 1000)],['ensemble1']) for item in func_exp2(xvals_b)]} + + for key in funcs_num.keys(): + [item.gamma_method() for item in ys[key]] + + num = pe.fits.least_squares(xs, ys, funcs_num, num_grad=True) + auto = pe.fits.least_squares(xs, ys, funcs_auto, num_grad=False) + + assert(num[0] == auto[0]) + assert(num[1] == auto[1]) + +def test_combined_fit_dictkeys_no_order(): + def func_exp1(x): + return 0.3*np.exp(0.5*x) + + def func_exp2(x): + return 0.3*np.exp(0.8*x) + + xvals_b = np.arange(0,6) + xvals_a = np.arange(0,8) + + def func_num_a(a,x): + return a[0]*np.exp(a[1]*x) + + def func_num_b(a,x): + return a[0]*np.exp(a[2]*x) + + def func_auto_a(a,x): + return a[0]*anp.exp(a[1]*x) + + def func_auto_b(a,x): + return a[0]*anp.exp(a[2]*x) + + funcs = {'a':func_auto_a, 'b':func_auto_b} + funcs_no_order = {'b':func_auto_b, 'a':func_auto_a} + xs = {'a':xvals_a, 'b':xvals_b} + xs_no_order = {'b':xvals_b, 'a':xvals_a} + yobs_a = [pe.Obs([np.random.normal(item, item*1.5, 1000)],['ensemble1']) for item in func_exp1(xvals_a)] + yobs_b = [pe.Obs([np.random.normal(item, item*1.4, 1000)],['ensemble1']) for item in func_exp2(xvals_b)] + ys = {'a': yobs_a, 'b': yobs_b} + ys_no_order = {'b': yobs_b, 'a': yobs_a} + + for key in funcs.keys(): + [item.gamma_method() for item in ys[key]] + [item.gamma_method() for item in ys_no_order[key]] + for method_kw in ['Levenberg-Marquardt', 'migrad', 'Powell', 'Nelder-Mead']: + order = pe.fits.least_squares(xs, ys, funcs,method = method_kw) + no_order_func = pe.fits.least_squares(xs, ys, funcs_no_order,method = method_kw) + no_order_x = pe.fits.least_squares(xs_no_order, ys, funcs,method = method_kw) + no_order_y = pe.fits.least_squares(xs, ys_no_order, funcs,method = method_kw) + no_order_func_x = pe.fits.least_squares(xs_no_order, ys, funcs_no_order,method = method_kw) + no_order_func_y = pe.fits.least_squares(xs, ys_no_order, funcs_no_order,method = method_kw) + no_order_x_y = pe.fits.least_squares(xs_no_order, ys_no_order, funcs,method = method_kw) + + assert(no_order_func[0] == order[0]) + assert(no_order_func[1] == order[1]) + + assert(no_order_x[0] == order[0]) + assert(no_order_x[1] == order[1]) + + assert(no_order_y[0] == order[0]) + assert(no_order_y[1] == order[1]) + + assert(no_order_func_x[0] == order[0]) + assert(no_order_func_x[1] == order[1]) + + assert(no_order_func_y[0] == order[0]) + assert(no_order_func_y[1] == order[1]) + + assert(no_order_x_y[0] == order[0]) + assert(no_order_x_y[1] == order[1]) + def fit_general(x, y, func, silent=False, **kwargs): """Performs a non-linear fit to y = func(x) and returns a list of Obs corresponding to the fit parameters.