{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "Custom line with additional jitter term\n", "========================================" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Optimization terminated successfully.\n", " Current function value: 384.634484\n", " Iterations: 2\n", " Function evaluations: 69\n", "-------------------- Parameter summary ----------------------\n", " const = -0.355683, free: T, restricted: F, related: F\n", " slope = 1.10116, free: T, restricted: F, related: F\n", " jitter = 2.97802, free: T, restricted: F, related: F\n", "-------------------------------------------------------------\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deZxcVZn/8c+TdAgQtgRaCCTSAcJmhOC0oIBJB1nC7gIMDCqbRkSQRUaNKHULx1FBBLcBw88A/giIoDMGZItAkiGA0oEYIBATSCKJgTQEWQSy9TN/3Hurb1VX9VZVXUt/369Xv6ruqbpV53VDP3147jnPMXdHRETqy6BKd0BEREpPwV1EpA4puIuI1CEFdxGROqTgLiJShxoq3QGAHXbYwZuamirdDRGRmjJ//vxX3b0x32tVEdybmppobW2tdDdERGqKma0o9JrSMiIidUjBXUSkDnUb3M1supmtMbNnEm23m9mC6Ge5mS2I2pvM7N3Ea9eXs/MiIpJfT3LuNwE/A34VN7j7v8bPzexq4I3E+19w9/Gl6qCIiPRet8Hd3eeaWVO+18zMgFOAw0rbLRERKUaxOfePAa+4+5JE2xgze8rM5pjZxwqdaGZTzKzVzFrb2tqK7IaIiCQVG9xPA25LHK8G3u/uBwCXALea2Tb5TnT3ae7e7O7NjY15p2mKiNS9ICjP5/Y5uJtZA/Ap4Pa4zd3Xuftr0fP5wAvAnsV2UkSkXqXT5fncYkbuhwPPu/vKuMHMGs1scPR8N2As8GJxXRQRqU8tLeX77J5MhbwNeAzYy8xWmtk50Uunkp2SAZgALIymRt4JnOvua0vZYRGRWhcEYAZz5oTHZuFPKVM0Vg07MTU3N7vKD4jIQGQGfQ3DZjbf3ZvzvaYVqiIidUjBXUSkglKp8nyugruISAVV3VRIERHpvXIF81wK7iIi/ahc89pzKbiLiPSDIMie1x4E9H2aTA9UxU5MIiL1Jgg6UjBBkD1iN4uezHuUYNYhZfl+jdxFRMogGcxnzw4ffeMmAN5sGA5A8LEHy/b9Cu4iIiUWp186rURtGAzANhtfD49Tl5d8ZWpMwV1EpERyg3k8eo/nsvvOu5D61NNMnBDm2t3Dn3IEd5UfEBEpgWSOHaKyAr+cDi+8AN/9bnj83joYOrTj9SLDr8oPiIiUWdYUxwULwsdzziG4ZQ/YsCEcvUeBHcq3MjWm4C4iUkBuuiR5nHyemeL4xhtw4YXwL/9Casur+P5eN5L+21kwZEiXn1UOSsuIiBSQmzpJHpuFo+98i5JSzX8geOBgWj45nDlzyjedXWkZEZFeyt1II3mcnA3jzy7Cv3kZEAbx1L+/A8cei40YXtZ67d3RIiYRkYSCC47yHIfP9yW1ecdW0emrtsyaAVOKG6d9oZG7iEhCEHRMUYSO55njdsfvuDN8Thjpg5fOIZUq77Z5vaXgLiIDVp/SJO+8AxdcAMB54x8NP+dnO5BO5982r9yzYgpRcBeRAau7Co2ZwPzOO3DNNaS+tQmGDSP49NMAXLfgo1mfk1mslFic1J959qSebJA93czWmNkzibbAzFaZ2YLo55jEa1PNbKmZLTazo8rVcRGRYvQkhRIEwMyZ8IEPwCWXEBx0b9j+sx3ypm4qFcjz6cnI/SZgcp72a9x9fPRzD4CZ7QucCnwgOue/zGxwqTorIlKsTvVeCs1kWbYMjj8eTjwRhg0Lq38dd1y3n1+pNEyuboO7u88F1vbw804Efu3u69x9GbAUOLCI/omIlFS+G6Zxe5Z/+zd4+GG46ip46imYODHv5+UG82oZvReTcz/fzBZGaZvhUdsuwEuJ96yM2joxsylm1mpmrW1tbUV0Q0Skd/KuLgWYNQteDys2csMN8PzzcOmlMGRIjz6rmvQ1uF8H7A6MB1YDV/f2A9x9mrs3u3tzY2NjH7shItJ76XRHUM5Kzxx5BMEnorow48bBqFEV6V8p9Cm4u/sr7r7J3duBG+hIvawCRifeOipqExGpClmrS9etx7//AwB8iy1JTZpL8MDBnc6p1tF5V/oU3M1sZOLwk0A8k2YmcKqZDTWzMcBY4M/FdVFEpHh5b6QO3YzgG++GDYsWkX54QlblRgj/GPTXptal1JOpkLcBjwF7mdlKMzsHuNLMnjazhcAk4GIAd38W+A2wCLgP+LK7bypb70VEupDcwzR+7n9fHT46+IvLCO5qDleXntlUgR6Wj6pCikjdiYN5XNfFDA6bsJGH5jbg22yLvflG1iyZfCPziRM7RvlJqVT1pGm6qgqp4C4idcesIzjvuiusWNH5PblBuqWFTHnerkr9VhOV/BWRAaOpKXyMR925gT11ueddTTp7dvhYTcW/iqHgLiJ1IU7D5AbzsaPfA8Df/mf4vnRODd/E+VBdxb+KobSMiNSXJ57ADvwwPu6D2DNPZ6VZelrIq1rTMLmUlhGRupQVqNeuJWi+Gw46KDyeOpXU5WGEjkfe1XIjtD9o5C4iNSszwn7mGZg0CXu1Db/oYoLNv0/wvaHdnl9IJUv19oZG7iJSFwrVhAlu34d7Go4PD665pqjAnvs9tUp7qIpIzYjno3fe43QwMD1xXF3z0StBI3cRqQmZmjBjZ+A7jcQtDF/HHPIGUN0bZ1SCgruIVLVONWE+czr28mpaDgiD+j3ztg3breP9ouAuIlUuuLw9HI23h0Nzv/4X+MZNzJ6/NalU9oh9oKdikhTcRaRqZAVmd7jjDth773DLu3ho/sUvwuDBnd+f53ggU3AXkaoR3ygNLngNJk+GU06BLbckuHJLoPBK0VpcQVpuCu4iUhUyUxtTKdI/2x4ef5xg8uN8fNtW0tfvCGSX8E3SiL0zBXcR6Vf5AnPWDdMrwuH7PtusIn3fQWyyzjO2a3HzjP6m4C4i/So3MAdnLAOgeceXstqfX7kV0LmQV1z1Ubqm4C4i/SarnO66dfCd77Buj30BeOKHYRTfddf8506cGD7GVR+TVRulMwV3ESm7vPuXbj4Uu/zbbN4e7mFqn/0MAGeeGb4nOcURwnrr8QKluF2LlQpTcBeRskgG3bgmeiYwX/VDfOye+P0P5A3WubNfNBum97qtCmlm04HjgDXuPi5quwo4HlgPvACc5e7/MLMm4DlgcXT64+5+bnedUFVIkfqS3LIOgPXrsaGb4b+5AzvlZHz9Bmhvh6Fhga9C9dMLVWeslaqN5VZsVcibgMk5bbOAce6+H/BXYGritRfcfXz0021gF5H6Fpy5nOXDx4cHs2aFo/AhQzKBHQqPzAsFcAX27nUb3N19LrA2p+0Bd98YHT4OjCpD30SkBuSmX3Jz6+mbmxjzzqLw+IZppNOap94fSpFzPxu4N3E8xsyeMrM5ZvaxQieZ2RQzazWz1ra2thJ0Q0T6Q24gTk5tDILsm57rLBydH3FoeNM09yaognr5FBXczewyYCMwI2paDbzf3Q8ALgFuNbNt8p3r7tPcvdndmxsbG4vphoj0o2Qwz5raCPDoowSfWJAJ2kN9HQCzHtkC6PoPg5RWnzfrMLMzCW+0ftyju7Luvg5YFz2fb2YvAHsCulsqUgcyNdWDfBtmQIoHSBPgd24EGjJpmnybU3f6wyAl1aeRu5lNBr4GnODu7yTaG81scPR8N2As8GIpOioi/aurXHoc2FPfbgfAh4/AG4aw5+j3whcaGgrm1fPOeddipJLryVTI24AWYAfgFSBFODtmKPBa9LbH3f1cM/s0cAWwAWgHUu5+V3ed0FRIkepTaHpiVvvChdj++5Ha9SbSK87s9N54FkyhwF3oO6RnupoK2W1w7w8K7iLVpdM89QQz8FtmwOmnAxB8YRXBtJ0zuZneBGwF9+IUO89dRAaILlMm7e1w442ktrgyrBHwUljoK7hhl46key9p5Wn5KLiLSEbuVMbM1MVPLCDY9UY4+2yCD82E1lYYPbrTudC7gK08e/kouItIlk4B9623YOJE0ivPgRtvhLlzYf/9C05rVMCuDgruIpIlnQbcSZ2+NBy2b701l+56R/jimWfCoEEd74toWmP1UXAXGeCSUxQzQfqwwwhmjCX4zFLM4OqnjwQ6b5ihaY3VS8FdZICKA3A6HQbrdDoRpGc/jOGw++5ARw4+zqfHG2Zk5rtH7aqxXj0U3EUGqHS6Y6TetGsYvZ1w1ouvacMdZs/NDhEFb7gG/dNn6TkFd5EBJgg60irxSH3O3GiOOmHUbjm5sVO6JT63EE1rrC4K7iIDSFwTJk6rxJJplVSqZ1va5QZzjd6ri4K7yAAye3b46E1jwsfTw31Lc+vI9ISCeXVTcBepU11uorF8Wdi+xy2ZEXihtIrSLbVJtWVE6lSmbsu6dbB+PcHVWxMc/SfsIwfh6zcQfHeIRt81TrVlRAaYzHz1+++HcePg0kvDaYsHHRS2D1Fgr3cK7iI1rsv0y+SjsKVLOPh/vpZ5j9IsA4OCu0iNy93DNPWvz+PDtgIgNWkuAI+tCRcjmZF3g2qpPwruIjUsWdMl+NZGANK37w3HHhu2PTSBVEqLjgYiBXeRGpSvpkv6uw2ZlabcfjupVBj8tQn1wKTgLlKFuhtZBwH4ho34tT8G4Ft8B4AVf+vYDSmdhuXLO85Rrn1gUXAXqULdjraXLIHmZoKLXgfgP/h21ssTJ4aP8UpU5doHnoaevMnMpgPHAWvcfVzUNgK4HWgClgOnuPvrZmbAj4FjgHeAM939ydJ3XaQ+dVkb3T2M1DvuSPDq+QR3DoeFTpC2zLx2s46VqKB9Sgeqno7cbwIm57R9A3jQ3ccCD0bHAEcDY6OfKcB1xXdTpP51WRt90ya4/nqYMAE2bIBttiG96vPw6U8TpLP3L1X6RaCHwd3d5wJrc5pPBG6Ont8MfCLR/isPPQ5sZ2YjS9FZkXoVBPnL6QIEx7USjP4lfOlLMGQIrF2bd3QfB/Xc1IuC/cBUTM59R3dfHT1/Gdgxer4L8FLifSujtixmNsXMWs2sta2trYhuiNS+fPuPTv7YP8MnBx5IevUUuPVWWjY9iO20Y97RfaF8uvLsA1OPcu7dcXc3s15l9dx9GjANwtoypeiHSK0Jguz8ePJG6v2PDAPAvD187+LTmDO3Y0SvXLp0pZiR+ytxuiV6XBO1rwJGJ943KmoTkYS4tnruhhjcdhu+9vVO6ZR8o3uRQooJ7jOBM6LnZwC/T7R/zkIfAd5IpG9EJBKP2DsF8b+eho0YnjnOF+TNOqY7iuTTo+BuZrcBjwF7mdlKMzsH+D5whJktAQ6PjgHuAV4ElgI3AOeVvNciNSLfKDt3Vkxmk2nCJ7729UyJgFSq4zNySwgk0zkiuVTPXaSMcvPiuTc+zcCPP4Fg/vEEvz8A+3Bz3jx6fJ7y7JKkeu4iFZBvumI6Dbz9NsGhf4SlS8PGX/2K4G9nQ3NzwWmL8R8ETWuUnlJwFymxQouRMkW99tmH9LzD4d57w2C93XYweHDm3O4+W6QnFNxFSix3MVI82s4U9VoZLgMJXrtAwVrKRsFdpITyBev4xudU+15WezzrRQFeykHBXaQHehqAs3ZFOm0x0JGe+Z5PBTpG8to4Q8pJwV2kBwqV4E0G5swN1GXL4IQTCH69N37SyVlTGJNTG0XKqSTlB0TqWVcleOOgnwz+ttsYYCapw+cR3Hpg1vs160X6i0buIgV0WYKXjqCfWWQ07Ybw8aST8ZdWEsw6JKziSOdgrtG7lJuCu0gB+UrwJot25daECV46J3xyxx0walSnzxLpTwruIr2Qqbu+bj2pw+dlvZb+zqDMe0QqTcFdpAfitEo6TTi3cfx4gj8eCoC/uCx81OwXqSK6oSrSA0GQuLE6aRI0NcHMmaTmA2PGVKxfIoVo5C6SEI+6k6Pv4PL27Bw7ji1fFhb7it6n2S9SbVQVUiQhrrqYqb746KNw3nnw/PPwwgvYqF0yqRelX6TSVBVSpBvJtEsm/XL22XDIIfDaazBjBuy8c+b9hRY1iVQLBXepWz0dWedud5dJv9w4HcMJPvsCfPrTYEYq1fWiJpFqoeAudasnJQOgo7CXL1sePsbz2p9dFJYL+M/NMufl7nlaqPCXUjZSacq5S11qaQmDcL7/vON8ehysc6VSYXtW7r3AZxSiHZOkPyjnLgNGT0sGQDgLxqffyOtDGgHwCy8i9fX3MnuX9jX9orSNVIM+B3cz28vMFiR+3jSzi8wsMLNVifZjStlhka70pmSADR5EcPYKtvvwnmHjtdcSfH/zHqVf8k197O4Pi0h/KklaxswGA6uAg4CzgLfd/Yc9PV9pGSmH3NSIGfibb2HbbB2233ADwR8+TPC7/QiuGJQ3CPclvdJVSkiklPojLfNx4AV3X1GizxPptdzgnBxdt0yMIu2ee3Y0fuELpH8/HgblD+x9Fd+gFamkUgX3U4HbEsfnm9lCM5tuZsPznWBmU8ys1cxa29raStQNGchyb47GC43MYM7caP/Sl1dnXutJbryvK0+1YlUqrejgbmabAScAd0RN1wG7A+OB1cDV+c5z92nu3uzuzY2NjcV2Qwa4QoE62PpqvGEIPnwEAL5xU2Y2TDmnNCrPLpVWipH70cCT7v4KgLu/4u6b3L0duAE4sMuzRYpQ8Cbm5e1hw267wec+B4vD/UwZPLjgTVcFZKknpQjup5FIyZjZyMRrnwSeKcF3yADWVdDtFKj/ugQ/ajLBZv8ZNnzykwSjfwmNjUqVyIBSVHA3s2HAEcDvEs1XmtnTZrYQmARcXMx3iPRopem774aP48bBY4/BDjt0Oj/fHwkFfKlXRQV3d/+nu2/v7m8k2j7r7h909/3c/QR3X118N2Wg6snm1MFZK+ADHyBFACefHFZwPPfcbs8HpWKkfmmFqlSl3qw0Td+0KwwbRvBwC9xyC4wcqQVFMuApuEtVinPpEyeGx91uTv2phVkRXzdNZaBTcJeqFQSdFwQFBz8AEKZgEtJXmEbmIgkK7lK14px6KgWsXg0nnQRHHQVAcP/BQM9G5rppKgORgrtUpazqjQEEVw3DfnsnRhjN7agjM6/lym3TaF4GIgV3qRpZ5QJyboSyzTakvrkh70g9d2SuLfBEoKHSHRCJxRtk8PbbBC+fi824BW8aA/Pmwc47YzaE4Ludz0uOzFVLXSSkkbtUhTgo/3SPH5O+eiu4IypV9OyzsPPOWUFbtdRFuqdt9qSiCm51d/5rsP32QOGt8AoFbm1xJwOFttmTqhV8OSz3nLpsY1Z7+mfbZ4K65quL9J6Cu1TGpk1w/fW8ufNeAAQf/1+gNEFcUx9FFNyljAoF5mDK3wlG/xL70rlsu3EtAHbYpC7PgZ4HbY3qRZRzlzLKm/t2xwYZvtNI7OXVTJzgzJlrWe+Lp0SKSNeUc5d+lzUlsb0dbr4ZXnuNlklhMZhjd38egNlzrNO5CuwixdM8dymZOCgnZ7eEhb0GMZFdmcP2mfZ75m2beF2jdZFSU3CXksksQgKCS97Ett0GH9wAw4fDlVfCGRNg0KBMukZTFkXKR2kZKYk4DdPSEo3cL7kkbPjCF8L9S886CwbpPzeR/qLfNilK7srQ+DHY+mpSn18F110HI0Zk3gsds140ZVGkfDRbRooWfHM9zJlD+tEjOr0WB/DZs8PAXwX/uYnUjbLOljGz5dGG2AvMrDVqG2Fms8xsSfQ4vNjvkfIo5iZmkHL47W9Jf28zgkfDEry+JlxxmlyEpCqNIv2vVGmZSe4+PvEX5BvAg+4+FngwOpYqVEzgTV9h4QYaAPPmhaP0xsas9zQ1hY8q6CXSv8qVcz8RuDl6fjPwiTJ9jxShN+VxM8H43XdhyZJM0M5snnHIwaTTHfXV41z8ihXZn5NKqTaMSL9w96J+gGXAk8B8YErU9o/E65Y8TrRPAVqB1ve///0u/SeVipMm2T+pVP73T5wYvu4zZ3pq22t6da57dG7iUURKA2j1ArG5FCP3Q939Q8DRwJfNbELOHw8HOt1Gc/dp7t7s7s2NOf8rL+UVBL2rtLjTe8vDJyecAEOG4A89zMSJHefGn9kdzY4R6T9FL2Jy91XR4xoz+2/gQOAVMxvp7qvNbCSwptjvkf6Vvdq0CYhSMK8Cc8LZL/Hq0u6Cdvy6UjEi/aeokbuZDTOzrePnwJHAM8BM4IzobWcAvy/me6R8CgXmOH/uGzfhX/s6QGa0njtfvbugraAu0v+KTcvsCDxiZn8B/gz8wd3vA74PHGFmS4DDo2OpInHAzX1k5UpmN54cPl+7FgYPJtjiB0CeTatFpGppEdMAlVvXxQxSh88j/cdDOr1XC5FEqlNXi5hUOGwAyp0Cedyh/wC2I/jjoaRx/MVl2G5j8gZxjdhFaoOC+wCSu1q0I1BvFx5Hk5qCm8cU/AzNeBGpDUrLDEQbN2JDGvD37QiPPEIwY2zelaqa5SJS3bQT0wDS6QZpbvs5L8GHPxwejB8PDQ2Z1/LNe1ddGJHapOBeZ+JgnBuU02ngi18kPX00vPoqqZMXwX33wZgwBZMv3dKb8gQiUl0U3OtEEGRvmJFpT3nm+NZ7w9w6zz1H8Jt9s+6OJueu59ZoV7EvkdqjnHsdaGrqXKCrJ+JA3hVthSdSvZRzr3NxhcZMznzt6wBsiv55/e4/ZL/eTS0ZEal9Cu41JhmQ86VPAFretwiAwbSH7ccd2+ncntLUR5HapOBeY5I3SvPNcklN/hOzD/waqS++nN2e6qi13hsa3YvUJuXca0hLS/by/8zxRRdj114Ttm/aFA7hB4V/t5UzF6lfyrnXuEKzVzLH116TeR+DB2cCOyitIjJQKbjXgE6bayx6jlTTzThhkt3//ETBG6RKq4gMTAruVaLQitK8liwhvfwMuO668DhecSoiElFwrxJ5V5TG3AlOWQQ//SmpFLT86ISw/dxzlXYRkbxUFbIK5C7zzzpesgQuuID0/ffBo4+SXtXxUrL8rtIvIpKk4F5BhUvwJo/HMtXC/e2C5WcSNHS8plkwIlKI0jIVFM87j/cmjW+axkE7ZVcA8D2fCoANaVCNFxHpEY3cKyydDoO5GbB8Odx1F1xwAQDBktMJdg/fp5G6iPRGn0fuZjbazB42s0Vm9qyZXRi1B2a2yswWRD/HlK679SWTW1+3jtSkubDPPgRffQtefjm8Ubr77gXP1ehdRLpSzMh9I/BVd3/SzLYG5pvZrOi1a9z9h8V3rz7FK0tjtvlQYALLt7iZmzecQrBT5+CtWTEi0ht9Hrm7+2p3fzJ6/hbwHLBLqTpWy7obVcclBI49OKze6GP3xO+7n+UHnlLwszRSF5HeKMkNVTNrAg4A/hQ1nW9mC81supkNL3DOFDNrNbPWtra2UnSjanS1Nd3HJ2wIn7hz97zw0gSnLMImH5V3cwxtcycifVF0cDezrYDfAhe5+5vAdcDuwHhgNXB1vvPcfZq7N7t7c2NjY7HdqBqFtqaL68M89L9DALBBhlk4Uyb4j4bs8gLR89mz+6PHIlKPigruZjaEMLDPcPffAbj7K+6+yd3bgRuAA4vvZvXramu64NK3CZZ+hhQB3hTuWdpVANc2dyJSrD6X/DUzA24G1rr7RYn2ke6+Onp+MXCQu5/a1WfVW8nfrGmL7tggwzcbiq1fh//zHWzYlgWnNQZBdhDXFEgRKaRcJX8PAT4LHJYz7fFKM3vazBYCk4CLi/iOmpI1sv7zn+G992iZFC47PX3/Z8L2LbfscuaLRuciUgp9ngrp7o8Aluele/rendoUj7bTaQjOf5XU+HkEBz1FOpGRuvWJsUDv68FoCqSI9IVWqBYpnrM+5+F2YBDsuSfBW2/B1y4h+PbbsNVWmdRKX1IsGsmLSF8ouBchCMKKAQCz50bb2r2+FoDUFhBsVZl+iYiocFgPFNrhKJ2GFSuy21OXe2ZXpPi8OLWiFIuI9BcF9zySq0JbWvIsJGpvZ/Nf3wSAn31O+BilW4J0R1I9Pk+rTEWkvym45xEH5XyrQ4MvvYINHsTUxWcCYNN/GbYH2SPzQouZRET6g4J7jjgoNzWFj7kLifjFLwDw6TeGjx4G9TgNowVIIlINFNwjuUG5Uy49Bb5yFR/c+dWw4ayzss5NPs9XSkDBXUT6k4J7JF9QBvBJh2UabNQunLTqJ0DHfHUFbRGpRgruhXz96+HjggWkjm0lSHmvR+SaHSMilTIgg3uXo+2FC0lt/SO48kpS438PixcT3N0Mg3p/qTSqF5FKGVDBPQ62ubNgggDYuDE82H13ggkPwSOPEDx1IuQpR6wRuYhUuz5XhSyl/qoKGddPj3dCAuDdd7Ett8D32z8s9jV0aMHzcys2iohUUrmqQlatfAG40NTGr249LWwYNw7eeafLz9WuSCJSK+oyuCeDcDzFMXdq48QRTwPwo00XAmC3zsBGDC84MteiJBGpJXUX3HODcDzFceLE8DhOx8ze7yvh8br1mfZ8s1+0KElEalHdBPfugvDsy2aFT159lVQKWtofCo8326zbz9WiJBGpNXUV3PMG4c+vhFNOgSOPJDX8JwTf2kg6DXPmhquQtBhJROpR3QT3TtwJjnwU9t6b4Hf7wRVXEPx9CsH1O3VK0+SW6C1EUyBFpFbUxVTI3MCcLODlxx2P3X1X3h2QcndG0mbUIlJLKjIV0swmm9liM1tqZt8o1/dAYnbMyy/DGWcQnPp85sZqy1t3FTxPJXpFpF6VJbib2WDg58DRwL7AaWa2bzm+KxOUf/IT2Gsvglv2wPbZO3NjtatZLirRKyL1qlwj9wOBpe7+oruvB34NnFjKL+gUlC/8CvbmG3DeeVk3VvPl1nM/R7NhRKTelCu47wK8lDheGbVlmNkUM2s1s9a2trZef0EmKF/67wD4HXfi7U7w0+2z3jd7dq8/WkSk5lVstoy7T3P3ZndvbsxTnKvH4iH2SSdl7V8Kvd+YWrNhRKRelCu4rwJGJ45HRW2lN2xYJijnrfZIz1MsSsWISL0oV3B/AhhrZmPMbDPgVGBmmb6LINBsFxGRpLIEd3ffCJwP3A88B/zG3Z8tx3dptouISGd1sYgppkVIIjKQDLh67iIiA11dBXfNdhERCdVVcFeeXUQkVFfBXUREQgruIiJ1SMFdRKQOKbiLiNQhBXcRkTqk4C4iUoeqYoWqmbUBK4r4iB2AV/0DgcQAAATISURBVEvUnXKo9v6B+lgq6mNpqI89s6u75y2rWxXBvVhm1lpoCW41qPb+gfpYKupjaaiPxVNaRkSkDim4i4jUoXoJ7tMq3YFuVHv/QH0sFfWxNNTHItVFzl1ERLLVy8hdREQSFNxFROpQTQd3M5tsZovNbKmZfaPS/QEws9Fm9rCZLTKzZ83swqh9hJnNMrMl0ePwCvdzsJk9ZWZ3R8djzOxP0bW8Pdr7tqLMbDszu9PMnjez58zso9V0Hc3s4ujf+Bkzu83MNq+G62hm081sjZk9k2jLe90s9JOovwvN7EMV6t9V0b/zQjP7bzPbLvHa1Kh/i83sqHL3r1AfE6991czczHaIjvv9GvZEzQZ3MxsM/Bw4GtgXOM3M9q1srwDYCHzV3fcFPgJ8OerXN4AH3X0s8GB0XEkXEu5vG/sBcI277wG8DpxTkV5l+zFwn7vvDexP2N+quI5mtgvwFaDZ3ccBgwk3gq+G63gTMDmnrdB1OxoYG/1MAa6rUP9mAePcfT/gr8BUgOh351TgA9E5/xX97leij5jZaOBI4G+J5kpcw+65e03+AB8F7k8cTwWmVrpfefr5e+AIYDEwMmobCSyuYJ9GEf6CHwbcDRjhSruGfNe2Qn3cFlhGdNM/0V4V1xHYBXgJGAE0RNfxqGq5jkAT8Ex31w34BXBavvf1Z/9yXvskMCN6nvV7DdwPfLQS1zBqu5NwoLEc2KGS17C7n5odudPxyxVbGbVVDTNrAg4A/gTs6O6ro5deBnasULcArgW+BrRHx9sD/3D3jdFxNVzLMUAbcGOUPvp/ZjaMKrmO7r4K+CHhCG418AYwn+q7jrFC160af4/OBu6NnldN/8zsRGCVu/8l56Wq6WNSLQf3qmZmWwG/BS5y9zeTr3n4570ic1DN7DhgjbvPr8T390ID8CHgOnc/APgnOSmYCl/H4cCJhH+EdgaGked/46tRJa9bd8zsMsLU5oxK9yXJzLYEvglcXum+9FQtB/dVwOjE8aioreLMbAhhYJ/h7r+Lml8xs5HR6yOBNRXq3iHACWa2HPg1YWrmx8B2ZtYQvacaruVKYKW7/yk6vpMw2FfLdTwcWObube6+Afgd4bWttusYK3Tdqub3yMzOBI4DTo/+AEH19G93wj/kf4l+d0YBT5rZTlRPH7PUcnB/AhgbzU7YjPCmy8wK9wkzM+CXwHPu/qPESzOBM6LnZxDm4vudu09191Hu3kR4zR5y99OBh4GTKt2/mLu/DLxkZntFTR8HFlEl15EwHfMRM9sy+jeP+1dV1zGh0HWbCXwumvHxEeCNRPqm35jZZMJU4Qnu/k7ipZnAqWY21MzGEN60/HN/98/dn3b397l7U/S7sxL4UPTfaVVcw04qnfQv8obHMYR31l8ALqt0f6I+HUr4v7wLgQXRzzGEee0HgSXAH4ERVdDXFuDu6PluhL80S4E7gKFV0L/xQGt0Lf8HGF5N1xFIA88DzwD/HxhaDdcRuI3wPsAGwiB0TqHrRngz/efR79DThLN/KtG/pYR56/h35vrE+y+L+rcYOLpS1zDn9eV03FDt92vYkx+VHxARqUO1nJYREZECFNxFROqQgruISB1ScBcRqUMK7iIidUjBXUSkDim4i4jUof8D7Rwtwlt+5+oAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "from __future__ import print_function, division\n", "import numpy as np\n", "import matplotlib.pylab as plt\n", "from PyAstronomy import funcFit2 as fuf2\n", "import scipy.optimize as sco\n", "\n", "np.random.seed(1234)\n", "\n", "class LinMod(fuf2.MBO):\n", " \"\"\" Linear model with additional jitter \"\"\"\n", "\n", " def __init__(self):\n", " fuf2.MBO.__init__(self, pars=[\"const\", \"slope\", \"jitter\"], rootName=\"LinMod\")\n", "\n", " def evaluate(self, x):\n", " \"\"\" Evaluate model \"\"\"\n", " return self[\"const\"] + x * self[\"slope\"]\n", "\n", " def logL(self, x, y, yerr, **kwargs):\n", " \"\"\" ln(Likelihood) including jitter as additional term \"\"\"\n", " yr = np.sqrt(yerr**2 + self[\"jitter\"]**2)\n", " m = self.evaluate(x)\n", " lnl = -len(x)/2.0*np.log(2.*np.pi) - np.sum(np.log(yr)) - 0.5 * np.sum((m-y)**2/(yr**2))\n", " return lnl\n", "\n", "# Instantiate model\n", "lm = LinMod()\n", "# Starting values\n", "lm[\"slope\"] = 1.1\n", "lm[\"const\"] = -0.5\n", "lm[\"jitter\"] = 1\n", "\n", "# Use -log(L) as SciPy-like objective function and call it logl\n", "lm.addSPLikeObjf(\"-logl\", \"logl\")\n", "\n", "# Get some 'data' and add Gaussian noise with STD 1+3**2=10\n", "x = np.arange(150.)\n", "y = lm.evaluate(x) + np.random.normal(0,np.sqrt(1+3**2),len(x))\n", "# Nominal error has STD 1\n", "yerr = np.ones_like(x)\n", "\n", "lm.thaw([\"slope\", \"const\", \"jitter\"])\n", "\n", "fr = sco.fmin_powell(lm.logl, x0=lm.freeParamVals(), args=(x,y,yerr))\n", "lm.setFreeParamVals(fr)\n", "\n", "lm.parameterSummary()\n", "\n", "plt.errorbar(x, y, yerr=yerr, fmt='b+')\n", "plt.plot(x, lm.evaluate(x), 'r--')\n", "plt.show()" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.10" } }, "nbformat": 4, "nbformat_minor": 4 }