{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "Restrictions via algorithm with convenience function\n", "===========================================" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "------------------- Parameter summary --------------------\n", " A = 0.536203, free: T, restricted: F, related: F\n", " mu = 10.3404, free: T, restricted: F, related: F\n", " sig = 5, free: T, restricted: T, related: F\n", " off = 1.0047, free: T, restricted: F, related: F\n", " lin = 0, free: F, restricted: F, related: F\n", "----------------------------------------------------------\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEICAYAAABRSj9aAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deZgU1bn48e87w4AsgWEZlJ2ICAybKAYMrlejaBLJjTfucY8hmh+a6I1xJ4tXDUajxjXGuCRq1BjFqBE0qDGuoCyCKCAomwKyCCog8P7+ONVOTU8v1TNV3dU17+d5+pnpquqqU9XVb50659Q5oqoYY4xJropSJ8AYY0y0LNAbY0zCWaA3xpiEs0BvjDEJZ4HeGGMSzgK9McYknAX6GBKRJSJySJZ5rUXkcRHZICIPicgJIjKlgHWfIiIvhpfa8iAiA0RkpohsFJEJpU6PqU9EnhKRkxvxud4isklEKqNIV1K0KHUCkkBElgA7A9uBL4CXgPGqujSCzf2Pt63OqrrNm/YXX1oU6K+qCyPYdjn7GTBNVfcodULCFufvXEQmArup6om5llPVwwOubwlwhqo+433uA6BdE5OZeJajD8+3VbUd0A34CLgxou30Ad71BXkTTB9gbjE2VE65SxHJm9kLskwTti8iYnEoaqpqrya+gCXAIb73R+CCcer9N4E3gU+ApcDEtM9/H3gf+Bi4OH19vuV+AWzF3TVsAk4HTgFe9Oa/ACjwqTf/mAzrOAX4D/B7YAMwHzjYN78D8EdgJbAc+DVQ6Zv/A+BtYCMwD9jTm/5zYJFv+n/7PjMR+LPvfV8vnS18aXrP++xi4ATfsqd521sHPA30yfE9HIkL5uuB54BB3vR/4e62NnvHZfcMn30OuBJ4zfueHgM6+eY/BHzoHbMXgMG+eXcBtwBPesf+kFzfuW//T/XmrQPGA3sDs730/z4tfRmPQ7bvHPgWMNNb10vAsLTz9QJvW1tS30Pa9hQ4G1gALA6wzgu882Uj8A5wMDCW+ufrLN+xvgJ3Hn4O7OZNOyPXeQbcC+zwPrMJd5eWOpapc6k7MBlYCywEfpB2Hj4I3OOtdy4wMtc+lDq2hBajSp2AJLzwBWagDXA3cI9v/oHAUNwd1DBcjv873rxa76TdH2gFXAtsI0Og95afSP2geQpeoPfeK+5WOVtaT/HW/xOgCjgGF7w6efP/DtwGtAW64gLfD7153/N+CHsD4v1A+/jmdff28Rhc4OmWJc1f/ji97XwCDPDmdcMLosA478c6yFv2EuClLPu1u7fNb3j79TPvsy29+c/hCyQZPv+ct29DvDT9LS3NpwFf8b6j3wEzffPu8o7hGG//d8rznaf2/1Zv2UNxF6FHvWPeA1gFHBDkOKR/58AI7/OjgErgZNw52sp3vs4EegGtsxwPBaYCnYDWudYJDMBdsLr79q9fpu/ed6w/AAZ7+1Pl/37IfZ4toX6mKnUsU4H+BeBm77juAawG/suXls24jFgl7sL+ijcv6z4k4VXyBCTh5Z18m3A5nS+AFcDQHMv/DrjO+/8y4AHfvLa4XFCUgX4FIL5pr+HuKnbG5fBa++YdhyvbBpeTPCfgMZkJjMuS5i9/nN7+rgeOIi3oAE8Bp/veVwCfkSFXD1wKPJi27HLgQO/9l4EkS3qfA67yva/1vofKDMtWe+nv4L2/C9+FPcB3ntr/Hr75H+O7A8NdaM4NchzSv3Pc3cWv0rb/DnUXjiXAaXnSq3gBMt86cYF4Fe5OpirX+eo71r/MMC0V6LOeZ+QI9LgL13bgK775VwJ3+dLyTNp3/Ln3f9Z9SMLLysbC8x1VrcblJH4MPC8iuwCIyCgRmSYiq0VkA+42vYv3ue64nAQAqvop7kcfpeXqnd2e97109MHlrlaKyHoRWY/L3Xf1luuFK55pQERO8lq1pD43hLp9zMrb32Nwx2SliDwhIgO92X2A633rXIvL4fXIsKru3n6k1rsDd1wzLZuNv/L8fdyx6CIilSJylYgsEpFPcMGGtP2rV/Ge5ztP+cj3/+cZ3qcqGQs5Dqnlz0st732mF+4YZUxvFv5lsq5TXSXwubhAukpEHhCR7g1Xl3Xd6bKeZ3l0B9aq6kbftPepf5w+9P3/GbCTiLRo5D6UDQv0IVPV7ar6CC5nsa83+T5cuWEvVe2Au2UXb95K3IkNgIi0ATpHnMweIiK+971xufyluBx9F1Wt9l7tVXWwt9xSoF/6ykSkD/AH3AWus3fBe4u6ffwUV6SVsov/86r6tKp+A1dsM99bV2p7P/SlpVpVW6vqSxn2aQUuGKXSJLjjujzfwfDp5fu/N+7ubA1wPK745BBcHUbf1Gb8u5G2rlzfeaEKOQ6p5a9IW76Nqt6fI72Z+JfJuU5VvU9V98V9BwpcnWc7ubaf8TwL8LkVQCcR+YpvWm8CngM59qHsWaAPmdeKYBzQEVeZBK5sd62qbhaRr+ECR8rDwLdEZF8RaQn8kqZ9Lx8Bu+ZZpiswQUSqROR7uLLfJ1V1JTAF+K2ItBeRChHpJyIHeJ+7AzhfRPby9nM3L8i3xf0wVnvH4FRcjj5lJrC/1+a5A3BhaoaI7Cwi40SkLe4iswlX4QYuOF4oIoO9ZTt46c3kQeCbInKwiFQB53nryxYMMzlRRGq9i+0vgYdVdTvu+9uCu9NqA/xfgHXl+s4Lle84pH/nfwDGe3cVIiJtReSbaQGwUFnX6T2j8F8i0gpXBv45dd/hR0DfAlvWZDvPMu3rl9Q1Z34JuFJEdhKRYbgGC3/Ot8E8+1D2LNCH53ER2YSrWLwCOFlVU835zgJ+KSIbcWXyD6Y+5C1zNi4HuBLXqmJZE9IxEbjbu70+OssyrwL9cbnVK4D/UdVUcdFJQEtcS4d1uAtRNy+tD3nL34drmfAorhJ3HvBb4GXcD3EorkVFah+nAn/FtfKYAfzDl5YK4Ke43NhaXJnvj7zP/R2Xq3rAKzJ5C8jY3lpV3wFOxDVrXQN8G9fkdWuug5XmXlx5+4e4IrjUg1X34IoAlnvH5ZUA68r6nRcqwHGYiO87V9XpuFYrv8d9hwtxdTONlmedrYCrcMf9Q1xGInUxf8j7+7GIvBFwWxnPM2/2lcAl3r6en+Hjx+HuuFbgGhZcrl6b+zxy7UPZk/pFtcY0TyLyHK7S8I5Sp8WYsFmO3hhjEs4CvTHGJJwV3RhjTMJZjt4YYxIudr1XdunSRfv27VvqZBhjTFmZMWPGGlWtyTQvdoG+b9++TJ8+vdTJMMaYsiIi72ebZ0U3xhiTcBbojTEm4SzQG2NMwlmgN8aYhLNAb4wxCWeB3hhjEs4CvTHGJJwFemPKyXPPwYoVpU6FKTMW6I0pF9u3wzHHQM+ecMAB8PDDpU6RKRMW6I0pBzt2QGWly9FffjmsXAlHHw1LlpQ6ZaYMWKA3Ju5WrYKBA+Hpp2HQIBfon30WRODWW0udOlMGYtfXjTEmzaWXwuLF4O/sr1cv+L//g1GjSpYsUz4s0BsTZ7NmwR13wIQJMGBA/XkXXFCaNJmyY0U3xsTZr34F1dVw2WWZ5y9eDNdeW9w0mbJjgd6YuNq61ZXLH300dOyYeZknnoDzzoPXXitu2kxZsUBvTFxVVsLkyXDWWdmXOekkaNcObrqpeOkyZccCvTFxVVkJBx0EQ4dmX6Z9exfsH3gANm0qXtpMWbFAb0xcXX89BBlt7cgjXTHPK69EnyZTlizQGxNHa9bAT34CTz6Zf9mvfx06dIBly6JPlylL1rzSmDh65hlQhcMOy7/sV74CH3/sinqMycBy9MbE0ZQprqXNyJHBlrcgb3KwQG9M3Ki6ZpWHHBI8gM+eDSNGwMsvR5s2U5Ys0BsTN8uXw6efBiu2SenWDWbOhOefjy5dpmxZoDcmbnr2hDVr2H7ciVxzDXTpAr/9reulOKuaGtfh2QsvFC2ZpnxYoDcmhhYsbsHIMa2YONHVs15+Oey9NyxYkOND++8P//lPniuCaY4s0BsTJ6pwwAHcsNddzJ7tSnDA/Z01C8aMyfHZ/feHTz5x5fXG+FigNyZOPvgAXniBXXf5nB076s/asQOGDMnx2f33d/3iVNjP2tRnZ4QxcfLGGwAMPH5P2rWrP6tdOzjttByf7dkT/vpXGD48uvSZsmSB3pg4eeMNqKzk6+OH0SLtccYWLeDb3w6wjqVLXRGQMR4L9MbEyRtvwKBBdNilNevWuXideq1b53o6yOnee6F3b1i0qCjJNeXBAr0xcbLrrq6TssYaOND9nTUrnPSYRLC+boyJkxtvbNrnhwxxlbGzZ8NRR4WTJlP2LEdvTFxs28b2bRr8IalMWreG3Xe3HL2pJ2+gF5E7RWSViLyVZb6IyA0islBEZovInmnz24vIMhH5fViJNiaJPj7/Sta27sEVl28N/pBUJsOHW6A39QTJ0d8FjM0x/3Cgv/c6E7glbf6vAHsu25g8Xr3lDdZs68D6z1oCwR6S2r6dhncAP/iBG1TcWt4YT95Ar6ovAGtzLDIOuEedV4BqEekGICJ7ATsDU8JIrDFJtpfM4A3q3RBnfEgqFdw7dnQNbBp0k9D7YDjxRBApXuJNrIVRGdsDWOp7vwzoISIfAb8FTgQOybUCETkTdzdA7969Q0iSMWVm9Wp23rKUuS33hK11k9MfklqwwD38umCBy/GvX18378s7gK8rq56aAW3aQG1t8fbBxFaUlbFnAU+qat7xzVT1dlUdqaoja2pqIkySMTH15psAzKmqn6NPf0hqzBjq9YGT7ss7gMMOg+uuiyixptyEEeiXA71873t60/YBfiwiS4BrgJNE5KoQtmdM8nTvDj/5CY8vG5HzIanBg2nQB45fu3Zw2umCDhvOyqdnN771jkmUMAL9ZFwQFxEZDWxQ1ZWqeoKq9lbVvsD5uHL8n4ewPWOSZ8gQuPZaqK7Oudjpp9OgDxy/Fi1cac1f5g6nw9I5rPt4e+Nb75jECNK88n7gZWCA10zydBEZLyLjvUWeBN4DFgJ/wBXZGGMKMW8ebN6cd7Fvf5sGfeBUV7uy+tQdwNix8K81w2nD5+zGwmBdHJtEC9Lq5jhV7aaqVaraU1X/qKq3quqt3nxV1bNVtZ+qDlXV6RnWcZeq/jiKHTCm3G3f9Dk6ZAhXd7o6azFLqqVNv35wySWwbVvu4p031fVgOQzXN33eLo5NotmTscaU0IIFcMJe8xFVpn9em7GYZcECGDkyQzPKtKKY1MVg+nRY2HIwB/MMUzgUCNDFsUk0C/TGlEAqKA8YAC3enQfAPGozFrOkt7TJtIz/YrBpE2za2pJ/cTCf4LL6/tY7GR+yMolmgd6YIvMHZVWoZS5f0IIF9AcaFrNkammTvkymZpcjZQaXtL2uXvFO0LsDkywW6I0psvSgXMs8FtCfL3BdH6QXs2RqaZO+TKaLwUH6LL/69Kewtu7B9iB3ByZ5LNAbU2TpQXkS/8vP+M2X79MfksrU0iZ9mUwXg4U7DXX/vFXXH2GQuwOTPBbojSmy9KD8EmN4vt23uPfezK1oOnQg72hTmS4GC1t63R+8/XbWbYNV1DYHFuiNKTJ/UN6ZDzmCJ+hY+Umw8WCzyHQxmL2ul+vvxhfog9wdmOSxEaaMKbJUUAbg/mlw/PHw79nQYWi4G6qocEMLvvtu5m2bZsMCvTGlNHcuVFa6UaGi8NRT0KlTNOs2ZcMCvTGlNG8e7LYbtGoVzfq7do1mvaasWBm9MaU0b55rChOV996Ds86C+fOj24aJPQv0xhRJgydSP9sCCxdGOzjIli1wyy3w+uvRbcPEnhXdGFME6SNDXX453PfnKh55fDZ9Bufod7ipdtvNNavxtbwxzY8FemOKYMwY1+VA6mGlTz+FmbMr2PvkWlatinDDVVXQv78F+mbOim6MKYJMT6QeuuMpftbpjug3PmiQqwswzZYFemOKINMTqT+s/CNnbvhN5g+EqbbWdWBv3VQ2WxbojSmCTE+kDtR5tN4rwhY3Kb/8JSxa5Nrrm2bJyuiNKYIGT6R+8QW0WQDDvhP9xkWi34aJNcvRG1MK773nilMGDox+Wzt2wLhxcPPN0W/LxJIFemNKYeFC97cYgb6iwnVC/+9/R78tE0tWdGNMKXzzm7B+vetdshgGDbImls2Y5eiNiVDO8Vk7dHDt3Ith0CB45x1redNMWaA3JiI5x2e96CL405+Kl5hBg2DzZnj//eJt08SGBXpjIpJ1fNava/H7nxk+3CVo06bibdPEhgV6YyKSbXzW/QascuXzxaiITdl7b3jxRRg2rHjbNLFhgd6YiGQbn/WHB3hdBhcz0JtmzQK9MSFLVcBOmOCayvu1aAH71ZQo0J9xBowdW9xtmliwQG9MiPwVsOvWuV4HRoxww7aqummtt2+CXXaBnj2Lm7iKCpgxo7jbNLFggd6YEGWtgB3jW+i882DFChd4i2T7dpi2ciCsWcPNv1xjrSybmbxnmojcKSKrROStLPNFRG4QkYUiMltE9vSm9xGRN0RkpojMFZHxYSfemLjJVgE7ZEjagkXsfyZ1l3HjVFdU9MiV79Q18zTNQpAsxV1AroK9w4H+3utM4BZv+kpgH1XdAxgF/FxEujc+qcbEX7YK2NNO8958/jnssw889ljR0pS6y5i1ZQAAfTbPb3iXYRItb6BX1ReAtTkWGQfco84rQLWIdFPVraq6xVumVZBtGVPuMnVH3KKFmw64wvpXXnFjuRZJ6i5jCX25j+NYSq/MdxkmscIIvj2Apb73y7xpiEgvEZntzb9aVVdkWoGInCki00Vk+urVq0NIkjGlkeqOWLXutW6dmw7A/OK3uEndZeygkhO4j6kcWv8uw5OzuwZT1iLNZavqUlUdBuwGnCwiO2dZ7nZVHamqI2tqaqJMkjGlNX++K5/v379om0y/y2jPhvp3GeTprsGUvTAC/XKgl+99T2/al7yc/FvAfiFsz5jyNX8+9O0LrVsXbZP17jJ+fQUbKjuz7sMtdOhQl4sfMMC1DsrZWsiUrTAC/WTgJK/1zWhgg6quFJGeItIaQEQ6AvsC74SwPWPKV48ecNhhpdt+374uui9aVC8Xnypm8rNy/OTI2x+9iNwPHAh0EZFlwOVAFYCq3go8CRwBLAQ+A071PjoI+K2IKCDANao6J+wdMKasXHNNabefqhuYP58x42v5+OOGzUFTMpXjm/KUN9Cr6nF55itwdobpUwHrQcmYOBngmlgyfz6DB8Nzz2VfdNs2OOcc+OgjOPdcG1u8nFmTR2OK5Z//dN0ezCnhjW27di4N8+dnbfM/aRLssYcL7GvXWsVsEligN6ZY5s2D5cuhW7fSpuPii+Goo7K2+b/66gDdOJiyYmPGGlMsb7/tGql36VLadIx3vZF0wLXGSXfQQQ2LdKxitrxZjt6YYnn7bTekX6lt2VI/y54mbzcOpuxYoDemGFTjE+hffNENLfjqqxln5+3GwZQdC/TGFMMXX8Dxx5e2DX1K6mLz9tsZZ+ftxsGUHQv0xoQgbz8xLVvCjTfCd79bkvTV062bi9rz5pU6JaZILNAbE1C2YB6on5hPPnG5+jgQgdpaC/TNiAV6YwLIFcwDjSp16aXQtWvDfgZKxQJ9s2KB3pgAcgXzQKNKvf029OtX1JGlcvrRj+Duu+Nz4TGRskBvTAC5gnmg5ohxaXGTstdeMHZsfC48JlIW6I0JIFcwz9scceNGWLYsXoF++3Z48kl3W2ISzwK9MQHkCuaBR5WKU6CvqIBjjoE77yx1SkwRWBcIxgSQCuaNsssucNVVrvY2LkTchccqZJsFC/TGRK1XL7jgglKnoqHaWpg6tdSpMEVgRTfGNEJBA2nPmgUrVxYtbYHV1sKKFbB+falTYiJmgd6YAhU8kPbRR8OPf1zMJAZTW+v+ZukKwSSHBXpjChToAamUrVth0aJ4VcSm7L+/GwRlr71KnRITMQv0xhQo0ANSKe+848p1UrnnOGnf3iW6ZctSp8REzAK9MQUqqL/21LCBQ4dGnq5GefRRuO22wIsXVDdhYsMCvTEFKqi/9jlz3MzUoNxx89BDcOWVORdJBfeOHaF37wLqJkxsWKA3JodMOdiC+ms/4wx4+OH4Fo/U1sL778OmTRln+yue1693jXRsLNnyY4HemCwKbl2TSb9+MG5cVElsulTdQerp3TTpFc/pbCzZ8mCB3pgsCmpdk8nGjfCnP7lscFylAv3cuRlnZ6p49rOxZMuDBXpjsiiodU0ms2a5KPjmm6GnLTT9+sFOO8HixRlnZ6p49rOxZMuDBXpjsiiodU0mcW9xAy5Sf/SRK5/KIFPFc3W1K6+3sWTLhwV6Y9KkKmAnTIBt2+rPKygHO2eOi4K9eoWexlC1b591Vq6KZ2tqWT4s0Bvj46+AXbcOKithxAh4991G5GDnzHHlPHEf3OP11+Goo+DDDwN/JJSKalM0FuiN8WlyBWyKKrz1VryLbVI+/xweeaSguoTQjpMpiryBXkTuFJFVIvJWlvkiIjeIyEIRmS0ie3rT9xCRl0Vkrjf9mLATb0zYmlwBmyLibgMuvTS0tEVm2DD3t4DRpkI7TqYoguTo7wLG5ph/ONDfe50J3OJN/ww4SVUHe5//nYhUNz6pxkSvyRWwfjU10L17KOmKVHU19OlTUKAP9TiZyOUN9Kr6ArA2xyLjgHvUeQWoFpFuqvquqi7w1rECWAXUhJFoY6JSUPcGuTz2GPz617kbocfJ8OEFBfrQjpMpijDK6HsAS33vl3nTviQiXwNaAotC2J4xkSmoe4NcHnoIbr/djc1aDkaPdjn7gE1nQjtOpigiPwtFpBtwL3CqqmbM3ojImSIyXUSmr169OuokGRO9OXPKoyI25cIL4aWXXDMjkzhhBPrlgL+hcE9vGiLSHngCuNgr1slIVW9X1ZGqOrKmxkp3TJn74gs3alM5BXqTaGEE+snASV7rm9HABlVdKSItgb/jyu8fDmE7xpSH+fNdsC+3QH/YYS5nbxKnRb4FROR+4ECgi4gsAy4HqgBU9VbgSeAIYCGupc2p3kePBvYHOovIKd60U1R1ZojpNyZ+Fi+GVq3cE0XlZNMmV3xjEidvoFfV4/LMV+DsDNP/DPy58UkzpkwdeaTruTK9WUrcDR8O993nalYLfJp3+3a47jq46ip3U3DuuVbcHydl0iTAmDJTVRX/rg/SDR8OGza4gUgKYN0hxJ8FemPCtGUL7LcfTJ5c6pQUbvhw97eA9vRg3SGUAwv0xhBiT4xz5sCLL8LWraGmryiGDnVPPHXsWNDHrDuE+LNAb5q9UIseXn/d/S23iliAtm3dncj++xf0sVzdIVhXxvFggd40e6EWPUyf7qJanz6hprGo1q4tqOuGbN0h1NZa2X1cWKA3zV6oRQ+vv+6iWblVxKY8+CB07lxQNM7WHcLYsVZ2HxcW6E2zF1pPjDt2QP/+cOihoaWt6FJXt1eyPsgemJXdx4cFetPshdYTY0UF/O1vrhF5uRo40A0tGEKgt66M48MCvWn2QuuJMX2A2XJUUQGjRoUS6K0r4/iwQG9MWE45xbWhL3ejR9cvXG8k68o4PsrsGW1jYuy111zBdLk76ijo1ctFZpMIFuiNCcP69a6lyimnlDolTTd8eN1TsiYRrOjGNFuhPszz6qvu79e+FkraSu6992DatFKnwoTEAr1plkLviGvaNNeR2T77hJnM0pk4EY4/3opvEsKKbkyzNGaMC/Cpdt7+h3lWrWrECg84wDVLbNs21HSWzOjRcO+98MEH5f2UrwEsR2+aqdAf5jn8cLjooianKzZGj3Z/Q2hmaUrPAr1plkJ9mGfJEjdGbJKKOYYOhTZtXE+cpuxZoDfNUqgP89x4I+yxB2zeHEraYqGqyvViOXVqqVNiQmCB3jRLoT7MM20afP3r0Lp16Oksqd/9Dl54oeCPWdfE8WOB3pimWLsWZs6Egw4qdUrCN2AAdO1a0EdsWMF4skBvTFM8/7y7HUhioAe45x648srAi9uwgvFkgd6Yppg2zRXZjBpV6pRE4/nn4eqrA3fYZl0Tx5MFemOaYuJEmDIFWrYsdUqicdhhsGFD3RCJeVjXxPFkgd6YpujUCfbdt9SpiM7BB7vRsp5+OtDi1jVxPFmgN4kXWSuQKVPgN79JVrPKdJ07u9rUKVMCLZ6rNZO1xikdC/Qm0dJbgVx2GfTuDR07hhBsbrkFbrghucU2KUccAZWVTRpYxVrjlJZozJ7mGzlypE6fPr3UyTAJ0bVr/T5t/Nq2hd13h7/+1Q31WpCNG6GmBn74Q7j++lDSGluqTR7sPNP3UFHhbhga1beQaUBEZqjqyEzzLEdvEi1TK5CUJjX9+8c/YMsW+N73mpS+spAK8k0oorLWOKVlgd4kWqZWIH6NDjYPPgjdu7snYpuD6693o041Mthba5zSskBvEi1TKxC/RgUbVXeFOPZYV/7QHAweDGvWwOOPN+rj1hqntPKepSJyp4isEpG3sswXEblBRBaKyGwR2dM3758isl5E/hFmoo0Jyt8KZP16qK6uP79RwUYEHnvMNSFpLg46CHr0cH3UN4INFF5aQbIjdwFjc8w/HOjvvc4EbvHNmwR8v7GJMyZMoQWbjRvd3yZWUJaVyko44QR46ilYvbrUqTEFyhvoVfUFYG2ORcYB96jzClAtIt28zz4LbAwlpcbEwaZN0K2b69mxufn+910TywceKHVKTIHCKGDsASz1vV/mTQtMRM4UkekiMn215RZMBrF52ObRR11znREjSpSAEhoyBG6+GY48MtTVxua7TbBY1CSp6u2qOlJVR9bU1JQ6OSZmYvOwjSpcd53rvne//Yq88Zj40Y9CHUM2Nt9twoUxOPhyoJfvfU9vmjGhCH0g78aaNg3eeANuv735tLbJJHVXc8IJTV5VbL7bhAvjbJ0MnOS1vhkNbFDVlSGs1xggRg/bXH+9e8Tz+828fcHdd8P/+3+uV8smis13m3BBmlfeD7wMDBCRZSJyuoiMF5Hx3iJPAu8BC4E/AGf5Pvtv4CHgYO+zh4W+Bybxgj5sE3lZ7x/+4B6U2mmnkFdcZi67zDVXCqHrB3uQqjisrxsTexs2QN++rh18SnU1LFlS1zRywQI4+mj399NPm9iPjcnvu9+Ff/0LFv5iO4sAABBSSURBVC92PcQ1UpDv1gRjfd2Yshak/XukQ9itWAEHHABvvhnCyhJi4kQXpa+9tkmrsQepisMCvUmESMt6L70U/vMfiz5+w4bBhAmw666lTokJwAK9SYTIynqnToU774Tzz7eglu766+HUUxv1UWs7X1xWRm8SIZKy3k2b3C3BTju5YpvWrUNIacKouqeER40K3JOn1adEw8roTeJFUtb7+9/DBx/AH/9oQT6bTz91x+nYY2Ftrp5S6uSrT7Hcfvgs0BuTzXnnuUGxQ6nRTah27VxW/MMPXTFOgBKCXPUpQYd+tItBYSzQG5Nu+nRYuhSqquAb3yh1auJv5EiYNAkmT3ZdROSRqz4lPbf/2Weu0dP69XXdI0ydat0mFMrK6I3xe/11OOQQV+Y8ZUqpU1M+VF3b+ieegIULXTY8i1z1Kd/5Djz3XPbNVFTUDWFr48/WZ2X0JlH8t+2TJsFvfhPSLfzLL8Ohh0KnTnDHHaGlt1kQcd0XT5lSF+SzZCJz1acEGfqxQwfrNqFQFuhNLGUrg00vw73gArjwwoa38AWV4W7e7Mrj993XZS2nTcuZIzVZtGoFBx7o/n/4YfjWt1zuvgBBhn489ljrNqFgqhqr11577aWmeXv3XdU99lBt29bl99q2VR0xwk2vqVGtqPDnBeu/KipUO3bM/vl6tm93f7dtUx09WnX8eNX164u+v4l0223uwFdVqU6YoLp6dcGrWL9etbq6/vdbXa36/vuZpzf3rw6YrlniqpXRm9jp2rV+17VQVwY7eHDuMlxwdajbt2cpw1221XU1/OCDrrXIzJlQUwNbtrgcqQnPypXwi1+4YrCqKrj4YrjkklKnKrFyldGH0R+9MaHKFMxTZbCnneYaxWzalOmTSse2XzCo+wZWL1hHR9bxPn34iF0YuGMu928dDx2mu6KaqioYN86N/1pTY0E+Ct26wa23wjnnwC231D0NtWQJHH447LEHDB3qism6d3fvO3VyX7ZI8xqTN2LJCvRPPw0//WnD6ffdB8OHwyOPuH5L0j32GOy2mxvh/qqrGs5/5hl30t52G9xwQ8P5L73kaoiuvdY9XJNu5kwXWH79a7j//vrzWras6yzrwgtdEzW/jh3hxRfd/xMmwLPP1p/fs6fbb4AzznAVin4DBrj9BjjuONd2zW/ECPjzn93/Rx4JixbVnz9mjBtoA+Dgg117ab9DD61rUjd6NHzySf35//3fcMUV7v9hw9yYo34nnggXXeRy1N7wfI9uUD4U2KFwK+O5gXPo3fZjHp03mraXKft+tgNhO5Vs50ou5GbOph8LeYcBVH66A3zN7H7EzdzKj2jdpoJuNdvgyLNgn31cWXKXLpgiGDSo/u9m82b3KOxLL9Uff/bvf3fNbv75T1e+37KluwC3auV+Pw895J6+ffxx+MlP3G1a6mIg4uoFhgxxd2sTJzZMxxNPwFe/Cn/6k6vFT/fcc+528qab3JCJ6V57zT3GO2kS3HVXw/lz57q/Eye6tPq1bu1yKAD/+7/w5JP153fpAs8/33CdIUlWoG/fHmprG05PPdVYXZ15fio316lT5vmp2qEuXTLPr6x0f3feOfP81Mm4yy4N51dV1f3fvXvD+e3b1/3fq1fD+V271v3fp0/DwSD8w7599asNA23fvnX/9+vXMGfr//zuu7tj5NezZ93/Awe6hs/Z5g8e3LBWtFs391fky2YTrb+AuU8IX3wBH7EzANqiip3225vKVhX03bfC/chbtOCmo3blpsOBdZ3htxdC69Z8XtWec3/RkWWfdWQWwwFY1HIQLae/DNYvWekNHOgyV+BuzVascK9Us5ldd3VFPJs3w9atLhOwbZsrewP3d9SourK5VFF9mzbufadOmZvgpM7tLl0yz0/9zmtqMv+OU6OKZfudp3Tr1nC+fwyDnj0bzq+uzr6+EFgZvTHGJIC1ozfGmDyS3K2CBXoTmcgebDLGp9AAnWn59OczEtetQrZ2l6V6hd2Ofts21UmTVDt3Vr3mGvfeRCd1vKurVbt3r2vLLlLX/j1ru3ZjCpTrmQu/bOdlavlOnRo+n1FR4Z7bKBfkaEdf8sCe/goz0Ac9CeIurItV1Be99OOd65X6EdmF2DRFtgfoROrOp3znZUWFe64r07yDDir1HgbXbAN9ppMg7Kt0sYNnYy9WYawn377me2o1/TVqVP00tWnjclvV1Rb0TTAHHpj9/MqVW09/DRqk2q5d/Wnt2qnee2+p9zC4Zhvos50EYV2l04NnFIEqrItVU9cT5EKR60eX/mrXTvUrX8n+AyzXuy9TXPfe2zBAB82t+8/F224r/24VcgX6RFfGBh1HtLGVhvn6zh450nW61ZQKyLAGvW7qevKNCgT5ex70a9HCPcOWnqaUTOs3Jl2+TtB27HDPQuY6L1u0gGOOafoIZbFutZPtClCqV5g5+mydIvmv0uk51UIqDQvJwTY2h5opx9KYW8p868lXLBPk7ijI8fZvq1071Vatch+3ciojNaWV7RwvRm49DvWBNNeimyCC9IaYrXgj321jIevKJmjwzCYVVDt2VN1pp8zrCXKS5rpQFFJPUUiFbbmVkZrSasxvJaw6tmLUB+bTLAJ9Y7+wILnybLnKTCdWY9cVhaC5jCAnabYf0YwZheVksm2rc+fyLyM18ZIvJoSZCw9aHxhl443EB/qmfGH5cuVBc5VBgn6+dYV9EgTNZTSl0rrQnEzUFeTGqAaLCWHmwoMUsUZdvJP4QN+ULyxfgG5srrLQ28gwW/CkLhgtWgQLqk2pByg0cIdV52BMLkFiQq5zt9BMV5Dfe9TFO4kP9EnIJeaqKyjkyp+vDDxTUG1KPUChgbupdQ7GBBEkJmQ7dydNiibnHXWcalKgB+4EVgFvZZkvwA3AQmA2sKdv3sm43sEXACfn25Y2MtAXI5cYVdlavtx3oVf+fJXLhQbVfPttgdvEUZCYkO3c7dKl6TnvTL+bqONUUwP9/sCeOQL9EcBTXsAfDbzqTe8EvOf97ej93zHf9hoT6KMONlGVrRXSAiXTlT/TyRRmriHXflvXBSbOmhITmvobyva7mTEj2jjV5KIboG+OQH8bcJzv/TtAN+A44LZsy2V7Rdm8srHBKeyytVQ6RNwrSJAPWrGTap8eRq4h234HHnzbmDKUKefdsqWbFiRuBOl/JwpRB/p/APv63j8LjATOBy7xTb8UOD/LOs4EpgPTe/fuHclBaEquPMpccqbXfvs1vmKnkGaKjX1Iqqqq9G2GjYlKrgYaQeJGkP53osgU5Qr0segCQVVvV9WRqjqypqYmkm0EeYQ/m6BdKTQmHenatYMzz8z/OHa2Lg2GDQv2KHeQ/rez7fduu4XTLYMxcdShQ91vqKambgRBaBg3MnV7kKsrkGxxJ+ruE8II9MuBXr73Pb1p2aaXRLbAuG5d9gObOvgTJjQcarVFC9fPRhjpaMx6M51MLVvC668HO1GCXPgy9SPSogWce254Fz5j4ixXH1HZMku1tfn73/Fniooy6Em2rL7/Re6im29SvzL2NW96J2AxriK2o/d/p3zbiqqMPteDUZlup6KqgA2r5j2q28sgxVHW0sY0F7l+r0Hq7oL83sOqA6SJrW7uB1YCXwDLgNOB8cB4b74ANwGLgDnASN9nT8M1u1wInJpvWxphoM/3YFT6gY3q4YYogmRj0moPLhmTX6bfa4cOqr/6VbAHEoP83sOqA2xSoC/2qxidmgU5sOX0EFZj0mq5cmMK15gHErPJ1YtrYzJduQJ9LCpjw5avYiNI5WqYFbBRa0xa/RVOqVeh/W8b09zka0wRtI7NXy6/aRNs2dK49QSVuEAfpGIjWyWj/8AGWSYuyimtxjRWHAb2yNaY4qCDCsssZbpgVFS4Vj5RZLoSF+iDtCYJkpstpxxvOaXVmMYoSsuUAMK60w9r5LigEhfoi30Aiy0OuRpjiq0pz8GEKay752IXDScu0JdT2XqhipWrsYuJiZu4ZODCunsudnFr4gJ9ksuri5GricstsjF+ScvAFbu4NXGBPq7l1WHkkouRq4nLLbIxfknOwBVD4gJ9HIWVSy5GriYut8jG+MU1A1cuLNA3QqG587ByycXI1STtFtkYY4G+YI3JnYeVSy5GrsZukY1JHgv0BWpM7ryccsl2i2xM8ligL1BjcueWSzbGlJIF+gJZvzLGmHJjgb5Aljs3xpSbHOOgmExSuXNjjCkXlqM3xpiEs0BvjDEJZ4HeGGMSzgK9McYknAV6Y4xJOAv0xhiTcOIGD48PEVkNvF/qdDRCF2BNqRNRAs1xv5vjPkPz3O9y2uc+qlqTaUbsAn25EpHpqjqy1Okotua4381xn6F57ndS9tmKbowxJuEs0BtjTMJZoA/P7aVOQIk0x/1ujvsMzXO/E7HPVkZvjDEJZzl6Y4xJOAv0xhiTcBboQyIi54mIikgX772IyA0islBEZovInqVOY5hEZJKIzPf27e8iUu2bd6G33++IyGGlTGfYRGSst18LReTnpU5PFESkl4hME5F5IjJXRM7xpncSkakissD727HUaY2CiFSKyJsi8g/v/VdF5FXvO/+riLQsdRoLZYE+BCLSCzgU+MA3+XCgv/c6E7ilBEmL0lRgiKoOA94FLgQQkVrgWGAwMBa4WUQqS5bKEHn7cRPuu60FjvP2N2m2Aeepai0wGjjb28+fA8+qan/gWe99Ep0DvO17fzVwnaruBqwDTi9JqprAAn04rgN+BvhrtscB96jzClAtIt1KkroIqOoUVd3mvX0F6On9Pw54QFW3qOpiYCHwtVKkMQJfAxaq6nuquhV4ALe/iaKqK1X1De//jbig1wO3r3d7i90NfKc0KYyOiPQEvgnc4b0X4L+Ah71FynK/LdA3kYiMA5ar6qy0WT2Apb73y7xpSXQa8JT3f5L3O8n7lpGI9AVGAK8CO6vqSm/Wh8DOJUpWlH6Hy7Tt8N53Btb7MjVl+Z3bUIIBiMgzwC4ZZl0MXIQrtkmcXPutqo95y1yMu9X/SzHTZqInIu2AvwHnquonLnPrqKqKSKLaZovIt4BVqjpDRA4sdXrCZIE+AFU9JNN0ERkKfBWY5f0IegJviMjXgOVAL9/iPb1pZSPbfqeIyCnAt4CDte6BjLLf7xySvG/1iEgVLsj/RVUf8SZ/JCLdVHWlVwy5qnQpjMQY4EgROQLYCWgPXI8rdm3h5erL8ju3opsmUNU5qtpVVfuqal/cbd2eqvohMBk4yWt9MxrY4LvtLXsiMhZ3i3ukqn7mmzUZOFZEWonIV3GV0a+VIo0ReB3o77XCaImrdJ5c4jSFziuX/iPwtqpe65s1GTjZ+/9k4LFipy1Kqnqhqvb0fsvHAv9S1ROAacD/eIuV5X5bjj46TwJH4CojPwNOLW1yQvd7oBUw1bubeUVVx6vqXBF5EJiHK9I5W1W3lzCdoVHVbSLyY+BpoBK4U1XnljhZURgDfB+YIyIzvWkXAVcBD4rI6biuxI8uUfqK7QLgARH5NfAm7iJYVqwLBGOMSTgrujHGmISzQG+MMQlngd4YYxLOAr0xxiScBXpjjEk4C/TGGJNwFuiNMSbh/j/+mUBwcjAoswAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "from __future__ import print_function, division\n", "from numpy import arange, sqrt, exp, pi, random, ones_like\n", "import matplotlib.pylab as plt\n", "from PyAstronomy import funcFit2 as fuf2\n", "import scipy.optimize as sco\n", "\n", "random.seed(1234)\n", "\n", "# Creating a Gaussian with some noise\n", "# Choose some parameters...\n", "gPar = {\"A\":1.0, \"sig\":10.0, \"mu\":10.0, \"off\":1.0, \"lin\":0.0}\n", "# Calculate profile\n", "x = arange(100) - 50.0\n", "y = gPar[\"off\"] + gPar[\"A\"] / sqrt(2*pi*gPar[\"sig\"]**2) \\\n", " * exp(-(x-gPar[\"mu\"])**2/(2*gPar[\"sig\"]**2))\n", "# Add some noise...\n", "y += random.normal(0.0, 0.002, x.size)\n", "# ...and save the error bars\n", "yerr = ones_like(x)*0.002\n", "# Let us see what we have done...\n", "plt.plot(x, y, 'bp')\n", "\n", "# Create a model object\n", "gf = fuf2.GaussFit()\n", "\n", "# Set guess values for the parameters\n", "gf.assignValues({\"A\":3, \"sig\":3.77, \"off\":0.96, \"mu\":10.5})\n", "\n", "# 'Thaw' those (the order is irrelevant)\n", "gf.thaw([\"mu\", \"sig\", \"off\", \"A\"])\n", "\n", "# Restrict the valid range for sigma\n", "gf.setRestriction({\"sig\":[0.01,5]})\n", "\n", "# The convenience function 'fitfmin_l_bfgs_b1d' automatically channels\n", "# the restrictions from the model to the algorithm.\n", "fuf2.fitfmin_l_bfgs_b(gf, gf.chisqr, x, y, yerr=yerr)\n", "\n", "gf.parameterSummary()\n", "plt.title(\"Bad fit because of parameter restrictions\")\n", "plt.plot(x, gf.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 }