Created
January 3, 2020 17:02
-
-
Save megbedell/760edb8472935e966170bbccd0419a28 to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
{ | |
"cells": [ | |
{ | |
"cell_type": "code", | |
"execution_count": 1, | |
"metadata": {}, | |
"outputs": [], | |
"source": [ | |
"import numpy as np\n", | |
"import matplotlib.pyplot as plt\n", | |
"from scipy.optimize import least_squares, curve_fit\n", | |
"\n", | |
"%matplotlib inline" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 2, | |
"metadata": {}, | |
"outputs": [], | |
"source": [ | |
"tc = np.array([953,\n", | |
"1327,\n", | |
"1641,\n", | |
"1302,\n", | |
"1505,\n", | |
"1647,\n", | |
"1573,\n", | |
"1427,\n", | |
"1291,\n", | |
"1150,\n", | |
"1347,\n", | |
"1348,\n", | |
"1033,\n", | |
"1328,\n", | |
"1455,\n", | |
"1647,\n", | |
"1736,\n", | |
"1447,\n", | |
"1570,\n", | |
"1477,\n", | |
"1574,\n", | |
"1594,\n", | |
"1580,\n", | |
"1347,\n", | |
"1647,\n", | |
"1647], dtype=np.float64)\n", | |
"\n", | |
"ab = np.array([-0.005,\n", | |
"0.062,\n", | |
"0.054,\n", | |
"0.063,\n", | |
"0.112,\n", | |
"0.08,\n", | |
"0.097,\n", | |
"0.093,\n", | |
"0.084,\n", | |
"0.042,\n", | |
"0.05,\n", | |
"0.05,\n", | |
"0.003,\n", | |
"0.09,\n", | |
"0.19,\n", | |
"0.17,\n", | |
"0.17,\n", | |
"0.23,\n", | |
"0.19,\n", | |
"0.21,\n", | |
"0.18,\n", | |
"0.22,\n", | |
"0.14,\n", | |
"0.16,\n", | |
"0.14,\n", | |
"0.15])\n", | |
"\n", | |
"err = np.array([0.012, 0.008, 0.002, 0.004, 0.005, 0.017, \n", | |
" 0.013, 0.009, 0.016, 0.006, 0.005, 0.005, \n", | |
" 0.031, 0.004, 0.006, 0.011, 0.01, 0.01,\n", | |
" 0.031, 0.012, 0.014, 0.01, 0.021, 0.014, \n", | |
" 0.006, 0.066])" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"##### least-squares fit without uncertainties:" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 3, | |
"metadata": {}, | |
"outputs": [], | |
"source": [ | |
"def resid(par, x, y):\n", | |
" \"\"\"residuals of linear fit\"\"\" \n", | |
" m, b = par\n", | |
" return y - (m*x + b)" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 4, | |
"metadata": {}, | |
"outputs": [], | |
"source": [ | |
"soln = least_squares(resid, [0.,0.], args=(tc, ab))" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 5, | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"name": "stdout", | |
"output_type": "stream", | |
"text": [ | |
"best-fit slope is: 2.20e-04\n" | |
] | |
} | |
], | |
"source": [ | |
"print('best-fit slope is: {0:.2e}'.format(soln['x'][0]))" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 6, | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"data": { | |
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEOCAYAAAB4nTvgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XucVXW9//HXx+E2agoKXrhMgzKQmhdsQqpfnsoLUCpUdtKs7GRpJtapR5z02M/KslRONxFTUtMuimVGVBrqQe2myCAaqY0M92FQwRFKGYGZ+Zw/1hrYs9kzaw/stddee7+fj8d+zN5rfdfsz8xWPrO+t4+5OyIiIr3ZJ+kARESk9ClZiIhIJCULERGJpGQhIiKRlCxERCSSkoWIiERSshARkUhKFiIiEknJQkREIilZiIhIpH5JB1AoQ4cO9dra2qTDEBFJlSVLlmxy92FR7WJNFmY2GfgBUAXc4u7XZJ3/IvApoB3YCHzS3deE5zqAZWHTte5+Vm/vVVtbS0NDQ4F/AhGR8mZma/JpF1uyMLMqYDZwGtAMLDaz+e7+bEazpUC9u281s4uB64APh+fa3P2EuOITEZH8xTlmMQFocveV7r4dmAtMzWzg7g+7+9bw5ePAyBjjERGRPRRnshgBrMt43Rwe68kFwP0ZrweZWYOZPW5m0+IIUERE8hPnmIXlOJazeIaZfRSoB/4t43CNu7eY2RHAQjNb5u4rsq67ELgQoKampjBRi4jIbuK8s2gGRmW8Hgm0ZDcys1OBK4Cz3H1b13F3bwm/rgQeAcZnX+vuc9y93t3rhw2LHMwXEZE9FGeyWAzUmdloMxsAnAPMz2xgZuOBmwkSxUsZx4eY2cDw+VDgHUDmwLiIiBRRbN1Q7t5uZtOBBQRTZ29z92fM7Cqgwd3nAzOB/YFfmhnsmiJ7FHCzmXUSJLRrsmZRiYhIEVm51OCur693rbMQEekbM1vi7vVR7bTdh4iIRFKyEBGRSEoWIiISSclCREQiKVmIiEgkJQsREYmkZCEiIpHKpviRiPTdvKXrmbmgkZbNbQwfXM2MSeOYNr63/T6lUilZiFSoeUvXc/m9y2jb0QHA+s1tXH5vUG9MCUOyqRtKpELNXNC4M1F0advRwcwFjQlFJKVMyUKkQrVsbuvTcalsShYiFWr44Oo+HZfKpmQhUqFmTBpHdf+qbseq+1cxY9K4hCKSUqYBbpEK1TWIrdlQkg8lC5EKNm38iLJPDpoeXBhKFiJStjQ9uHA0ZiEiZUvTgwtHyUJEypamBxeOkoWIlC1NDy4cJQsRKVuaHlw4GuAWkbKl6cGFo2QhImWtEqYHF4OShYjkResVKpuShYhE0noF0QC3iETSegVRshCRSFqvIEoWIhJJ6xVEyUJEImm9gmiAW0Qiab2CxJoszGwy8AOgCrjF3a/JOv9F4FNAO7AR+KS7rwnPnQ98JWz6TXe/I85YRaR3Wq9Q2WLrhjKzKmA2MAU4GjjXzI7OarYUqHf344B7gOvCaw8CvgqcBEwAvmpmQ+KKVUREehfnmMUEoMndV7r7dmAuMDWzgbs/7O5bw5ePAyPD55OAB9291d1fAR4EJscYq4iI9CLOZDECWJfxujk81pMLgPv7cq2ZXWhmDWbWsHHjxr0MV0REehJnsrAcxzxnQ7OPAvXAzL5c6+5z3L3e3euHDRu2x4GKiEjv4kwWzcCojNcjgZbsRmZ2KnAFcJa7b+vLtSIiUhxxJovFQJ2ZjTazAcA5wPzMBmY2HriZIFG8lHFqAXC6mQ0JB7ZPD4+JiEgCYps66+7tZjad4B/5KuA2d3/GzK4CGtx9PkG30/7AL80MYK27n+XurWb2DYKEA3CVu7fGFauIiPTO3HMOI6ROfX29NzQ0JB2GSGqlcQvyNMZcasxsibvXR7XTCm4RSeUW5GmMOc20N5SIpHIL8jTGnGZKFiKSyi3I0xhzmilZiEgqtyBPY8xppmQhIqncgjyNMaeZBrhFJJVbkKcx5jRTshARIJ1bkGcnjK7B7bT9HGmgZCEiqaXps8WjMQsRSS1NnwVan4SW+6Pb7SXdWYhIalXs9NmO7bDuV/D8DbDpr3DgMXD4ZLBcG3YXhpKFiKTW8MHVrM+RGMp2+mzbC9B0c/Bo2wD7j4ETvw9HfCLWRAFKFiKSYjMmjes2ZgFlOH3WHV5eBI2zYN0voXMHHD4FTroVDp8EVpzRBCULEUmtsp4+27EN1twNz8+C1gbofwDUfRbqLoED6ooejpKFiKRaGqf89mrrelj+Q2iaA9s2wgFHQf1sGP0x6P+GxMJSshAR6YNYtkV3h41/Cgas190L3gkjzoRxl8Khp/Q4HlHMLdqVLERE8lTwdR3tW2H1nUGS2Pw0DBgCb/pC0N20/+jixhJB6yxERPJUsHUdr66Gpf8F80bCE58GOmHCHJjWDONnRiaKgsaSJ91ZiIjkaa/WdbjDiwuDAev1vwUMRk6DsZfCISf3eeprsdeYKFmIiORpj9Z17HgVVv806Gra8iwMHApHfRnqLob9RhU3lr2gbigRkTz1aVv0fzXBki8EXU2LPwv7DIKJP4Zp6+CEb+1VouhzLAWgOwsRkTxFruvwTtjwQNDV1HIfWD+oORvGfg6GTizoKutirzExd4/lGxdbfX29NzQ0JB2GiFSi7Vtg1R1BV9O/lsOgw2DMRVB3EVQfnnR0vTKzJe5eH9VOdxYiMSvmXHgpsi3PBQli1U+g/VUY+jY49usw6oNQNSDp6ApKyUIkRqq3UIY6O6Dl90FX0wsPwT4D4I3nwtjpcHDkH+ippWQhEqPe5sIrWaTM9ldgxa3w/I3w2irYdyQcfzUc+WkYNCzp6GKnZCESo4qtt1BONi8Ldnxd/TPoaAvWRIy/LlgjsU/l/BNaOT+pSAIqrt5Cuehsh+bfBF1NLz0KVdVQe17Q1TTk+KSjS4SShUiMKqLeQjl5fROs+FGw6+vWdbDfG+GE6+DIC2DgQUlHl6hYk4WZTQZ+AFQBt7j7NVnnTwa+DxwHnOPu92Sc6wCWhS/XuvtZccYqEoeyrrdQTlqfDO4iVt8FnduCnV7rZ8HwM2CfqujrK0BsycLMqoDZwGlAM7DYzOa7+7MZzdYCnwC+lONbtLn7CXHFJ1IsZVdvoVxk17Hutx8c+cmgq+nAo3u8rFKnQsd5ZzEBaHL3lQBmNheYCuxMFu6+OjzXGWMcIiK7dNWxXn4TvP5CWMf6e0Ed6wGDe720kqdCx5ksRgDrMl43Ayf14fpBZtYAtAPXuPu8QgYnIhUkZx3rycGOr8Mn513HupKnQseZLHJtgtKXvUVq3L3FzI4AFprZMndf0e0NzC4ELgSoqanZ80hFpDxl17Hu9wYYczGMvQQOGNvnb1fJU6HjTBbNQOa2iiOBlnwvdveW8OtKM3sEGA+syGozB5gDwd5QexmviJSLrc1BN9POOtZvgvobYPTH96qOdSVPhY5zi/LFQJ2ZjTazAcA5wPx8LjSzIWY2MHw+FHgHGWMdIiK7cYeX/gR//nf4TS08861gp9d3PwDveza4m9iLRAHF3xa8lMR2Z+Hu7WY2HVhAMHX2Nnd/xsyuAhrcfb6ZvRX4NTAEONPMvu7uxwBHATeHA9/7EIxZKFmIyO7a22DNncF4RB/rWPdVJU+F1hblIpJOr64OFs+tuAW2t8Lg44Jpr7XnQb99k44uNbRFuYiUH3d48eGwjvV8dtaxHvc55jUfwcx7nqdl88Ox/sWvdRYiIqUqVx3roy+DMZ+B/UYF6x9+Hf/6h0peZ6Ea3CJSurLrWFdVw8TbgzrWx1+9s451b+sfCqlY71OKdGchIqWlWx3r+8GqoOZDwQK6HupYF2v9g9ZZiIgkbcc/YeXt3etYv/nKvOpYF2v9g9ZZiIgkZctzsHg6/HoELPk8DDgY3v5zmLoGjvtaZKKA4q1/0DoLEZFiylnH+pygq2kP6lgXa/2D1lmUAa2zEEmB7DrW1SOg7mIY82kYdEjS0VUkrbMQkdKRs471tWEd6/5JRyd5ULIQkXjsVsd6ENR+tKLrWKeZkoWIFNbrGzPqWDerjnWZULIQkcJoXRJ0Na2ZG9axfg+8ZRaMOFN1rMuAkoWI7LmddaxnwabHdtWxrrsEBh+TdHRSQJHJwsz+lsf32ejupxQgHhFJg72oYy3plM+dRRXw3l7OG3kWNRKRFHOHl5+Axusz6lhPgXGXwuGT8q5jLemUT7K4yN3X9NbAzD5boHhEpNRk17Huf0BQWKjuEjigLunopEgik4W7/7kQbUQkZXarY30U1M+G0R/b6/Kkkj75jFksA3It8zbA3f24gkclIslwh41/Du4i1t0b7AA74sygq+nQU3Lu+CqVIZ9uqDPCrwb8nt7HL0QkjbLrWPcfHFsda0mnfLqhdo5XmNm2qPELEUmR19YE+zTtrGN9LEyYozrWshutsxCpNO7w4sKgbkRmHeuxlwZ7NqmrSXLIZ8zixIyX1WY2nqBLCgB3fzKOwESkwHLVsT7qy8Gur2F50kKZt3R9RW7jXc7yubP4TsbzF4DvZrx24D0FjUhECutfTfD8bFh5W1CNbsiJMPHHQf2IqkEFf7t5S9dz+b3LdtaqXr+5jcvvXQaghJFi+SSL89y9JfZIRKRwvBM2LAgGrDfcD9YPas4O61i/LdauppkLGncmii5tOzqYuaBRySLF8kkWt5rZEOAR4A/An929PdaoRGTPbN8S1LFePjusY31oUMd6zEWw7/CihNCSo0Z1b8clHfKZDTXFzAYB7wLeD/yPma0lSBx/cPe18YYoIpG2PBeMRaz6CbS/CgdPhLd9FWo+BFUDihrK8MHVrM+RGIYPri5qHFJYec2GcvfXCZMDgJmNBqYAN5jZYe4+Ib4QRSSnAtexLpQZk8Z1G7MAqO5fxYxJ4xKLSfZePrOhRrn7usxj7r4KuDFc3b0oruBEJIfsOtb7joTjr4YjP1USday7xiU0G6q85HNn8aiZ3QR8t2uswswOJZglNc7d3xpngCISylnH+rqwjnVpLZmaNn6EkkOZyWdP4bcARwJLzew9ZvZ54AngMeCk3i40s8lm1mhmTWZ2WY7zJ5vZk2bWbmZnZ50738yWh4/z8/+RRMpIZzus/RU89C6477ggUdSeB1OeglMfDWY4lViikPKUzwD3K8BFYZJ4CGgBJrp7c2/XmVkVMBs4DWgGFpvZfHd/NqPZWuATwJeyrj0I+CpQT7CWY0l47Sv5/mAiqfb6pow61utgv1rVsZZE5TNmMRi4luAuYjLBRoL3m9nn3X1hL5dOAJrcfWX4feYCU4GdycLdV4fnOrOunQQ86O6t4fkHw/e+K78fSySlWp8MBqxX3xXWsT4F6mfB8DNUx1oSlc/965PAjcAl4ZjFA2Z2AsEA9xp3P7eH60YAmQPjzUR0W0Vcqw5QKU+dO4Kupudnwaa/7qpjPXY6HHh00tGVPG0tUhz5JIuTs7uc3P0pM3sH8Klersu1RDRXXYw9vtbMLgQuBKipqcnzW4uUiLYXgsJCTTdB24aMOtb/AQMOTDq6VNDWIsWTzwB3zvraHvgRgJnl2kywGcjcnWwkwXhHPvK61t3nuHu9u9cPGzYsz28tkiB32LQI/vpR+E0NLPsqDD4e/u33cGYjvOk/lSj6oLetRaSw8rmzOMrM/tbLeQNy/de9GKgLF/CtB84BPpJnXAuAb4XbjACcDlye57UipafjdVjzi111rPu9oSh1rMu9i0ZbixRPPsniTXm06cg+4O7tZjad4B/+KuA2d3/GzK4CGtx9vpm9Ffg1MAQ408y+7u7HuHurmX2DIOEAXNU12C2SKlubgxlNTT8K61i/qWh1rCuhi0ZbixSPuec7jFDa6uvrvaGhIekwpALt9tf76WOZNnLV7nWsx06Hw04tWnGhd1yzMOc/pCMGV/OXy8qjskB2QoRga5Fvf+DYskmIcTOzJe4euT+MVvOI7IXMf6wG2jbeYQ8w7qnfwT9WwoAhidaxroQuGm0tUjxKFiJ7YeaCRg5iAx877Pd8+KAHGNLvXzzXVss1L3+Ryy7+RqJ1rCuli0ZbixSHkoXIngjrWF954Fc4ddQTOPDAlonc8fKZLHrtzRjGZQkmCtDur1JYShYifZFVx3rC/gfyw5fO5uetU9iwY9f07VL4611dNFJIShYioV6nme6sY/1j2LEFDnoLTLydP73y/5g9b3nJ/vWuLhopFCULEXJPM/3ve5/mkFcf5u0dd0HL/WBVQeW5sZfC0IlgxllA5z6D9Ne7lD0lCxG6rwTef5+tnD3kIT4+9Hccsa4FBh0W1LGuuwiqD9/tWv31LpVAyUKEYDrpkQPX8fGDf8cHhyxk/6o2nnxtHJ9/8Uv8YMbVRa9jLaWr3FfF90TJQipbZwe03MfddV9nQvUStnX243ebT+b2l89kWVsdIwZXK1HITpWwKr4nShZSmba/AituCwatX1vFsW84jO+3fJyfbjydlzsGA6U1UC2lobeNC5UsRMrJ5r8H23Cs+hl0bIVh74Tx11I9chq1T7/EoAWNWIV1L0j+KmFVfE+ULKT8dbZD82+CtREvPQJVg6D2o8FeTUOO39lMA9USpVJWxeeiZCHla7c61m+EE64N61gfnHR0kkKVvCpeyULKz251rN8Db7k+2PlVdaxlL1TyqnglCykPHduD7cCfvx42PaY61hKbSu2uVLKQdGt7AZpuDh5tG2D/I8M61p+AAYOTjk6kbChZSPq4w8uLoHEWrPsldO6AwyfDhFtg+GSwfErLi0hfKFlIenRsgzV3d69jPeZiGHsJHDA26ehEypqShZS+revDOtZzMupY3wCjPx57HWsRCShZSGlyh41/3r2O9bhL4dBTilbHWkQCShZSWtrbYM2dwQK6V55KtI51vhvGVerGclJZlCykNLy2Bp6/EVbcAttbYfCxMGEO1J6XSB3rfDeMq+SN5aSyaNqIJMcdXlgIf3w/zD8C/vEdOPTdcMojMOVpGPPpRBIF9L5h3J60E0k73VlI8WXVsWbgUDjqy1B3Mew3KunogPw3jKvkjeWksihZSPFk17EeciJM/DG88Zxgc78Sku+GcZW8sZxUFnVDSby8E1r+AI+8D347NribGD4FTvsLTG4IVlqXWKKAYMO46v7d95HKtWFcvu1E0k53FhKPHf+ElbcHdxL/eh4GHQpv/v8w5iLYd3jS0UXKd8O4St5YrtJU+qw3c/ekYyiI+vp6b2hoSDoM2fKP4O5h1R3Q/iocfBKM+xyMOlvlSSW1sme9QXAH+e0PHJv6hGFmS9y9PqpdrN1QZjbZzBrNrMnMLstxfqCZ3R2eX2RmteHxWjNrM7OnwsdNccYpe6mzA5p/CwtPh98fFdSQGPUBmLQYJj0OtR9RopBU06y3GLuhzKwKmA2cBjQDi81svrs/m9HsAuAVdx9jZucA1wIfDs+tcPcT4opPCmD7K7Di1mB9xGuroHoEHPfNYMrroEOSjq5kfGXeMu5atI4Od6rMOPekUXxz2rFJhyV9oFlv8Y5ZTACa3H0lgJnNBaYCmcliKvC18Pk9wA1m2seh5G1eFuz4uvpn0NG2s441I6fBPv2Tjq6kfGXeMn72+Nqdrzvcd75WwkgPzXqLtxtqBLAu43VzeCxnG3dvB7YAXfUuR5vZUjN71MzeGWOcko/Odlj7K3joXXDfccE6idqPwJSlcNofoeZDShQ53LVoXZ+OS2nSrLd47yxy3SFkj6b31GYDUOPuL5vZW4B5ZnaMu/+z28VmFwIXAtTU1BQgZNlNzjrW1wVV6FTHOlJHDxNIejoupUmz3uJNFs1A5nLckUBLD22azawfcCDQ6sEUrW0A7r7EzFYAY4Fu053cfQ4wB4LZUHH8EBVrtzrWp0D9LBh+hupY90GVWc7EUKXe1tSp1HKqXeJMFouBOjMbDawHzgE+ktVmPnA+8BhwNrDQ3d3MhhEkjQ4zOwKoA1bGGKtAUHFu7a+CJLHpr6pjXQDnnjSq25hF5nGRNIktWbh7u5lNBxYAVcBt7v6MmV0FNLj7fOBW4Kdm1gS0EiQUgJOBq8ysHegAPuPurXHFWvHaXggKCzXdFFnHutIXJvVV1yC2ZkNJ2mlRXiXbtAgar+9ex3rspT3WsS7nhUkilSrfRXna7qPS7KxjfQO0Lu5THeveFiYpWYiUNyWLSrG1GZbftFd1rLUwSaRyKVmUswLXsdbCJJHKpWRRjtq3wuqwjvXmp6H/4ILUsZ4xaVzOMYtKWpgkUqmULMrJq6th+Y3Bfk0x1LHWwiSRyqVkkXbu8OLDQVfT+vmABXs0jb0UDjm5z11NUSp9YZJIpVKySKsU1LEWkfKhZJE2KapjLSLlQ8kiDbwTNjwQdDW13A9WBTVnB11NQ99W8K4mEZFsShalbLc61ofBm6+Euoug+vCkoxORCqJkUYp2q2M9Ed7+c9WxFpHEKFmUis4OaLkv6Gp64UHYZwDUfBjGfQ4Ojty2RUQkVkoWSVMdaxFJASWLpGxeFnY1/TSoY33IyapjLSIlS8mimDrbg4VzjbPgpUeCqa615wWzmoYcn3R0IiI9UrIoBtWxFpGUU7KIU+uS4C5izVzVsRaRVFOyKLSO7bCuq471Y6pjLSJlQcmiUNpegKabg0dEHWsRkbRRstgb7vDyoqCrKbOO9YRbeqxjLSKSRkoWe2JnHetZ0NrQpzrWIiJppGTRFwWoYy0ikkZKFlFy1rE+I1gbcdip2vFVRCqCkkVP2ttgzZ3BeERXHetx/wljPwv7H5F0dCIiRaVkke21NcE+TStuCepYH/jmgtaxFhFJIyULKHodaxGRtFGyeHU1PPq+sI71wWEd68/AfjVJRyYiUjKULPYdCfsdAUfNUB1rEZEeKFns0w/e9dtE3nre0vXMXNBIy+Y2hg+uZsakcUwbPyKRWEREehPrEmMzm2xmjWbWZGaX5Tg/0MzuDs8vMrPajHOXh8cbzWxSnHEmYd7S9Vx+7zLWb27DgfWb27j83mXMW7o+6dBERHYTW7IwsypgNjAFOBo418yyd9K7AHjF3ccA3wOuDa89GjgHOAaYDNwYfr+yMXNBI207Oroda9vRwcwFjQlFJCLSszjvLCYATe6+0t23A3OBqVltpgJ3hM/vAU4xMwuPz3X3be6+CmgKv1/ZaNnc1qfjIiJJijNZjADWZbxuDo/lbOPu7cAW4OA8r0214YOr+3RcRCRJcSaLXIsTPM82+VyLmV1oZg1m1rBx48Y9CDE5MyaNo7p/95616v5VzJg0LqGIRER6FudsqGZgVMbrkUBLD22azawfcCDQmue1uPscYA5AfX39bskkH0nNSOp6D82GEpE0iDNZLAbqzGw0sJ5gwPojWW3mA+cDjwFnAwvd3c1sPnCnmX0XGA7UAU8UOsCuGUldA81dM5KAoiUMJQcRSYPYuqHCMYjpwALgOeAX7v6MmV1lZmeFzW4FDjazJuCLwGXhtc8AvwCeBf4AXOLuHdnvsbc0I0lEJD+xLspz9/uA+7KOXZnx/HXgQz1cezVwdZzxaUaSiEh+Krrup2YkiYjkp6KThWYkiYjkp6L3htKMJBGR/FR0sgDNSBIRyUdFd0OJiEh+lCxERCSSkoWIiERSshARkUhKFiIiEknJQkREIilZiIhIJCULERGJpGQhIiKRlCxERCSSkoWIiERSshARkUhKFiIiEknJQkREIilZiIhIJHP3pGMoCDPbCKxJOg5gKLAp6SB6oNj2TKnGVqpxgWLbU0nE9kZ3HxbVqGySRakwswZ3r086jlwU254p1dhKNS5QbHuqlGNTN5SIiERSshARkUhKFoU3J+kAeqHY9kypxlaqcYFi21MlG5vGLEREJJLuLEREJJKSRR7M7DYze8nM/p5x7CAze9DMlodfh4THzcyuN7MmM/ubmZ2Ycc35YfvlZnZ+THF9yMyeMbNOM6vPan95GFejmU3KOD45PNZkZpftbVy9xDbTzP4R/l5+bWaDSyi2b4RxPWVmD5jZ8PB40T7PnmLLOPclM3MzG1oqsZnZ18xsffh7e8rM3ptxLtHPNDx+afhez5jZdaUSm5ndnfE7W21mTyURW5+4ux4RD+Bk4ETg7xnHrgMuC59fBlwbPn8vcD9gwERgUXj8IGBl+HVI+HxIDHEdBYwDHgHqM44fDTwNDARGAyuAqvCxAjgCGBC2OTqm39npQL/w+bUZv7NSiO2AjOefA24q9ufZU2zh8VHAAoK1RENLJTbga8CXcrQthc/03cBDwMDw9SGlElvW+e8AVyYRW18eurPIg7v/EWjNOjwVuCN8fgcwLeP4TzzwODDYzA4HJgEPunuru78CPAhMLnRc7v6cuzfmaD4VmOvu29x9FdAETAgfTe6+0t23A3PDtnulh9gecPf28OXjwMgSiu2fGS/3A7oG84r2efYUW+h7wH9lxFVKseWS+GcKXAxc4+7bwjYvlVBsQHB3CPw7cFcSsfWFksWeO9TdNwCEXw8Jj48A1mW0aw6P9XS8WEotrk8S/FVcMrGZ2dVmtg44D7iyVGIzs7OA9e7+dNapxGMLTQ+7wW6zsDu2RGIbC7zTzBaZ2aNm9tYSiq3LO4EX3X15CcbWjZJF4VmOY97L8WIpmbjM7AqgHfh516EeYihqbO5+hbuPCuOaXgqxmdm+wBXsSl7dTvcQQzF/bz8EjgROADYQdKnQSwzFjK0fQTfcRGAG8IvwL/lSiK3Luey6q6CXGJL+90PJYi+8GN7yE37tusVtJuhf7jISaOnleLGURFzhYOsZwHkedtKWSmwZ7gQ+WCKxHUnQd/20ma0O3+dJMzusBGLD3V909w537wR+RNBdQinEFr7XvWE33RNAJ8HeS6UQG2bWD/gAcHdWzInHllMxB0jS/ABq6T54NpPuA9zXhc/fR/dBxyfC4wcBqwj+0hkSPj+o0HFlHH+E7gPcx9B94GwlwaBZv/D5aHYNnB0T0+9sMvAsMCyrXSnEVpfx/FLgniQ+z94+0/DcanYNcCceG3B4xvMvEPS3l8pn+hngqvD5WIJuHCuF2DL+f3g06f8X8v4ZivlmaX0Q3CZuAHYQZPgLgIOB/wWWh18PCtsaMJtg5sIyuv+D/UmCAasm4D9iiuv94fNtwIvAgoyBVFQHAAABTElEQVT2V4RxNQJTMo6/F3g+PHdFjL+zpvB/2KfCx00lFNuvgL8DfwN+C4wo9ufZU2xZ51ezK1kkHhvw0/C9/wbMp3vySPozHQD8LPxcnwTeUyqxhcdvBz6To33RYuvLQyu4RUQkksYsREQkkpKFiIhEUrIQEZFIShYiIhJJyUJERCIpWYiISCQlCxERidQv6QBEyomZdS3WBDgM6AA2hq8neLBjqEjqaFGeSEzM7GvAq+7+P0nHIrK31A0lIiKRlCxERCSSkoWIiETSALdIkYT1C64jKFqzxt2vTzgkkbwpWYgUz8XAb9z90aQDEekrdUOJFM+JwF+SDkJkT2jqrEiRmNlU4CygFfi2u7cmHJJI3pQsREQkkrqhREQkkpKFiIhEUrIQEZFIShYiIhJJyUJERCIpWYiISCQlCxERiaRkISIikZQsREQk0v8B9V/O3aus+NIAAAAASUVORK5CYII=\n", | |
"text/plain": [ | |
"<Figure size 432x288 with 1 Axes>" | |
] | |
}, | |
"metadata": { | |
"needs_background": "light" | |
}, | |
"output_type": "display_data" | |
} | |
], | |
"source": [ | |
"plt.scatter(tc, ab)\n", | |
"xx = np.linspace(min(tc), max(tc), 100)\n", | |
"plt.plot(xx, soln['x'][0]*xx + soln['x'][1], c='orange')\n", | |
"plt.xlabel(r'T$_c$')\n", | |
"plt.ylabel(r'[X/H]');" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"##### least-squares fit with uncertainties:" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 7, | |
"metadata": {}, | |
"outputs": [], | |
"source": [ | |
"def chi(par, x, y, yerr):\n", | |
" \"\"\"error-weighted residuals of linear fit\"\"\" \n", | |
" r = resid(par, x, y)\n", | |
" return r/yerr" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 8, | |
"metadata": {}, | |
"outputs": [], | |
"source": [ | |
"soln = least_squares(chi, [0., 0.], args=(tc, ab, err))" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 9, | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"name": "stdout", | |
"output_type": "stream", | |
"text": [ | |
"best-fit slope is: 4.23e-05\n" | |
] | |
} | |
], | |
"source": [ | |
"print('best-fit slope is: {0:.2e}'.format(soln['x'][0]))" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 10, | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"data": { | |
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEOCAYAAAB4nTvgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAHdJJREFUeJzt3X2QXXWd5/H3J52nFhI6T47QEIkYMsTFJdqCsxkfBh8SnVmScXWB1SrGcQt1h9mtmTJUUjiMi+sYyTxUzS6uYcrUOroDKIsxuwsTGTO6tS6EBIJE4kZCwNAdViMhAhrS6e7v/nHPTU5f7u1zutPn3nP7fl5Vt/qc3znn3m/3rT7fc35PRxGBmZnZWKa1OgAzMys/JwszM8vkZGFmZpmcLMzMLJOThZmZZXKyMDOzTE4WZmaWycnCzMwyFZosJK2WtF/SAUnr62z/Y0n7JD0m6TuSXpvaNizp0eS1rcg4zcxsbCpqBLekLuDHwHuAfmAXcG1E7Evt81vAzoj4laRPAu+MiKuTbS9FxNmFBGdmZuMyvcD3vhw4EBEHASTdCawBTiWLiPjH1P4PAh+Z6IctXLgwLrzwwokebmbWkR5++OGfR8SirP2KTBa9wDOp9X7gijH2/xhwX2p9tqTdwBCwMSK2jvVhF154Ibt3755orGZmHUnST/LsV2SyUJ2yunVekj4C9AHvSBUvjojDkl4H7JC0NyKerDnueuB6gMWLF09O1GZm9gpFNnD3Axek1s8HDtfuJOndwE3AVRFxoloeEYeTnweB7wIrao+NiNsjoi8i+hYtyryLMjOzCSoyWewClkpaImkmcA0wqleTpBXAZiqJ4mep8nmSZiXLC4GVpNo6zMysuQqrhoqIIUk3ANuBLmBLRDwu6RZgd0RsAzYBZwPfkARwKCKuAi4BNksaoZLQNqZ7UZmZWXMV1nW22fr6+sIN3GZm4yPp4Yjoy9rPI7jNzCyTk4WZmWVysjDrcFdvfoCrNz/Q6jCs5JwszMwsk5OFmZllcrIwM7NMThZmZpbJycKsg23dM8CeQ8fY+dRRVm7cwdY9A60OyUrKycKsQ23dM8CGe/YyODwCwMCx42y4Z68ThtXlZGHWoTZt38/xk8Ojyo6fHGbT9v0tisjKzMnCrEMdPnZ8XOXtzGNJzpyThVmHOq+ne1zl1tmcLMw61LpVy+ie0TWqrHtGF+tWLWtRRFZmRT4pz8xKbO2KXgBuvPsxBodH6O3pZt2qZafKzdKcLMw62NoVvdzx0CEA7vr4b7Q4mmJUuwcPDo+wcuMOJ8QJcjWUmU1Z7h48eZwszCy3MvYqGismdw+ePK6GMutwU7X6CTqre3DRfGdhZlOWuwdPHicLM5uy3D148rgaysymrKK6B1fbSKZyFV4tJwszm9I6oXtwM7gaysxy8XTmnc3JwswyebyCOVmYWSaPVzC3WZhZpnYfr+C2ijPnOwszy+TxCuZkYWaZPF7BCk0WklZL2i/pgKT1dbb/saR9kh6T9B1Jr01tu07SE8nruiLjNLOxrV3Ry+c/cCkzuyqnjN6ebj7/gUs9e2sHKazNQlIXcBvwHqAf2CVpW0TsS+22B+iLiF9J+iRwK3C1pPnAnwJ9QAAPJ8c+X1S8ZjY2j1fobEXeWVwOHIiIgxExCNwJrEnvEBH/GBG/SlYfBM5PllcB90fE0SRB3A+sLjBWMzMbQ5HJohd4JrXen5Q18jHgvgkea2ZmBSqy66zqlEXdHaWPUKlyesd4jpV0PXA9wOLFiycWpZmZZSoyWfQDF6TWzwcO1+4k6d3ATcA7IuJE6th31hz73dpjI+J24HaAvr6+uonIzCaP2yo6V5HVULuApZKWSJoJXANsS+8gaQWwGbgqIn6W2rQdeK+keZLmAe9NyszMrAUKu7OIiCFJN1A5yXcBWyLicUm3ALsjYhuwCTgb+IYkgEMRcVVEHJX0WSoJB+CWiDhaVKxmVtGJU29bPoVO9xER9wL31pTdnFp+9xjHbgG2FBedmaVVZ5UdHB5h5cYdk/Lch6KlY774pvu49YNvLH3M7cojuM2sLWeVrY15cHik9DG3MycLM2vLWWXbMeZ25mRhZm05q2w7xtzOnCzMrC1nlW3HmNuZk4WZteWssu0Yczvzw4/M7FQPohvvfozB4RF6e7pL3xuqNuaZXdM8E26BnCzMDGjPWWWrMe979gWWnzvXiaJAThZmdkq7JIm0n794gpdeHmLnU0fbZnxIO3KbhZm1ra17BnjquV+emmW0HcaHtCsnCzNrW5u272ekZgrRZoy12PfsC+x79oVCP6NsnCzMrG15rEVlPq/qnF5FcrIws7blsRbN42RhZm1r3aplTKt5VJrHWhTDycLM2tbaFb0sWXDWqUdr9vZ0e6xFQZwszJqgWfXKnWjhnFmcPXs6VyyZz/fXX1l4oti6Z4CXXh7ixZeHWLlxR0t7XlWnaK92Gy4yFicLM7OcqtOil6GrbrOnlXeyMDPLqUzTojc7FicLM7OcytRVt9mxOFmYmeVUpq66zY7FycLMLKcyTYve7Fg8kaCZtb3l585tyiSI1Z5Wf3TXowS0dCr3Zk8r72RhVrBq98bB4RHPijoFrF3Ry59864cAfH/9lS2PpVnTyrsayqxAze7eaFYUJwuzApWpq6XZmXA1lFmBytTVshNUR8m340OcJqpZv6uThVmBzuvpZqBOYvCsqJOnFYlh+blzm/6ZreZqKLMClamrpdmZ8J2FWYGa3b3RrCiF3llIWi1pv6QDktbX2f52SY9IGpL0wZptw5IeTV7biozTrEhrV/SyYnFP02ZFNStCYclCUhdwG/A+YDlwraTlNbsdAn4P+Ls6b3E8Ii5LXlcVFaeZ2Xg0c1rwMimyGupy4EBEHASQdCewBthX3SEink62jRQYh5nZpGg0bgaY8neMRVZD9QLPpNb7k7K8ZkvaLelBSWsnNzQzs/Hr5HEzRd5ZqE5Z1ClrZHFEHJb0OmCHpL0R8eSoD5CuB64HWLx48cQjNTPLoZPHzRR5Z9EPXJBaPx84nPfgiDic/DwIfBdYUWef2yOiLyL6Fi1adGbRmpllKNMU5c1WZLLYBSyVtETSTOAaIFevJknzJM1KlhcCK0m1dZiZtUInj5sprBoqIoYk3QBsB7qALRHxuKRbgN0RsU3SW4BvAvOAfy7p30fEG4BLgM1Jw/c0YGNEOFlY2+qk6SemsjJNUd5shQ7Ki4h7gXtrym5OLe+iUj1Ve9z/AS4tMjYzm1o8FXyxPN2HmbW9Zk0FX/2cak+dTppy3snCzNpes7q0dnLXWScLM2t7zerS6q6zZmZtrFldWt111sysjTWrS6u7zpqZtbFmTQXvrrNmZm1u7Ype7njoEFDsuJa1K3r5k2/9EIDvr7+ysM8pGycLM7Nx8mNVzczM6nCyMDOzTE4WZmaWycnCzMwyOVmYmVkmJwszM8vkZGFmZpmcLMzMLJOThZmZZcocwS3psRzvcyQi3jUJ8ZiZWQnlme6jC3j/GNsFbJuccMzMJs7POi9OnmTx8Yj4yVg7SPo3kxSPmU0RV29+APAJfKrIbLOIiP89GfuYWefYumeAPYeOsfOpo6zcuKMjnlE91eVps9gLp55PPmoTEBHxxkmPysza1tY9A2y4Zy+DwyMADBw7zoZ79gJ0xHMfpqo81VC/k/wU8D8Zu/3CzDrcpu37OX5yeFTZ8ZPDbNq+38mijWUmi3R7haQTWe0XZtbZDh87Pq5yaw8eZ2Fmk+q8nu5xlU+Gqzc/cKpB3YqRmSwkvan6ArolragpMzM7Zd2qZXTP6BpV1j2ji3WrlrUoIpsMedos/iK1/P+Av0ytB9A5D6E1s0zVdokb736MweERenu6Wbdqmdsr2lyeZPHhiDhceCRmNmWsXdHr5DDF5Gmz+LKkByVtlPROSXkSDACSVkvaL+mApPV1tr9d0iOShiR9sGbbdZKeSF7X5f1MMzObfHl6Q71P0mzgncDvAn8u6RDw98DfR8ShesdJ6gJuA94D9AO7JG2LiH2p3Q4Bvwd8qubY+cCfAn1UqroeTo59fny/npmZTYZcdwkR8TJJcgCQtAR4H/CfJL0mIi6vc9jlwIGIOJgccyewBjiVLCLi6WTbSM2xq4D7I+Josv1+YDVwR+7fzMwAT7thkyNPb6gLassi4qmI+CKwCfjNBof2As+k1vuTsjzO5Fgz6yCeWqQ58rRZfE/Sjem2Ckm/JulrwF9GxGCD41SnrN60IRM+VtL1knZL2n3kyJGcb23WmPvrt5dGU4s4YUy+PMnizcBFwB5JV0r6d8BDwAPAFWMc1w+k70rOB/L2qsp1bETcHhF9EdG3aNGinG9tZlPFWFOL2OTK08D9PPDxJEn8A5WT9lsjoj/j0F3A0qR9YwC4BvhXOePaDvyZpHnJ+nuBDTmPNbNEtYpmcHiElRt3TLnxDq2aWmTfsy8U+v5llKfNokfSZuCjVBqZ7wbukzTmYLyIGAJuoHLi/xHw9Yh4XNItkq5K3vstkvqBDwGbJT2eHHsU+CyVhLMLuKXa2G1m+XRCFU0rphbpVHl6Qz0CfBH4gyQBfFvSZcAXJf0kIq5tdGBE3AvcW1N2c2p5F5UqpnrHbgG25IjPrKW27hko5WjlTpj9dd2qZWy4Z++o33PKTC0SI3DiKJz4OZw4cvrny0dS6z+vrM9ZCr95Z6Hh5EkWb6+tcoqIRyWtBP51MWGZtYcyP7uhE2Z/baupRYZfTk70qZP/qRN/ej3ZPni0kjDqmT4HZi2E2Yug+zUw56LCw8+TLLYBr5gwMCIC+BsASY9EhCcVtI5T5qv383q6GaiTGKZaFc3aFb3c8VBlbHCzxpKIEeZM+yW88OOaK/46y9WfQ79s8GbTYOaCyol/1iI4Z3nl56yFpxPCrEWnf85aAF2zm/J7puVJFpdIemyM7QLOmaR4zNpKma/ep3QVzWQbPnG6Wqfeyb5m/ZGlR5iuEfgfdd6r61WjT/Jzfz052Scn/1k1J/+ZPZWEUXJ5ksWv59hnOHsXs6mnzFfvbVVFM5ki4OQLdU76tSf+VGIYerHBmwlmzT99kp9zMSz8Z2x56EWODc/lxjVvSyWC5G5g+qua+us2y7ielGc2lU2km2nZr95bUUUz6UZOwonnXtm4m0oEn+5+krk6Bt/8VaVs5GT99+qaPfrq/uyLkiv8hamr/dSJf+Z8mNb1irf5j/dtB+DGJauK/M1LJfcMsmZT2UQbqqvbNm3fz+FjxzmvU67eJyoChl7KbtxNr5881vj9Zs6DWQuZqdn8bOQ8XnvexacTQboqqLrc9SpQvQkiLIuThRln1lDd0c9uGBlKrvpru3fWrKdP/iMn6r/XtBmpht1FsKCvwRV/dXkBTKucwm6uTpZ4RfF3T1v3DPDSy0METMmBjo04WZhR7obqMzWu6qehX2Y27o4qH3yehlO+zTjn9Am++3yYt6L+yb+6PH3OhK/6m1XFVr0Drf7GZeoqXTQnCzPK3VA9YacGdeWo6qmWDTdIjpo+uitnzz9tXM9ffXXNbO7v2wRl7ipdNCcLM8rfUA3A0PF8V/vjGdQ1a2FlUFfPpaPr99Mn/9kLYUaP6/qZ2negWZwszGhBN9MYgcFjOa74j5zuBTT8q/rvpWmjr+irg7pqT/4tHtQ1FUzJO9CcnCzMEmfUzfTUoK6cJ/8Tz0E0GJ7U9arUiX0RzL2kfs+eahKYOa8tBnVNBW1xB1oQJwuzV4jkqj+ruie1POagrgWnq3XmLoNZK1/Zsyd98p+ig7qmgo4d6IiThXWC4UEYfC5XVc+XznqWOToGdze66p89+iQ/Z+noOv6cg7qsfU2JgY4T4GRh7SWichXfsKqnzvLJXzR+v5nzT5/g57yeR44s4cXoYe1bLxt9tX9qKoez3NBrHcnJwlorc1BXnUndRho89n3azNFX92cvqdO4m7oDmDn/1KCuqtt/VBnctfaSzrliNMvDycImT0RlUNcrZu8cq3vn843fb8Y5p0/sZy2G+W+uf/KvLk8/21f9ZgVxsrDGRoYrJ/NGT+mqN6Pn8Mv130vTR1frzLusQbfOVGPvtBnN/X3PwNXV6SY6qA7bOouTRSdpOKir0UCv52g4lcP0OadP6t3nQs8bR3f3rE0CM+a2xVW/T/Zm9TlZtKv0oK7aqp7aSdyqSaDhoK6upHtn9Uld/+SVjbujEsEC6JrV3N+3Bcr6bG2zVnCyKIvhE2P06qnTz3+sQV3Tz0r14Hk1zF3eYEBXkgRmnONBXTXGM2X5RJ6DYdZunCyKEFHprlmvqiedDNLbh15q8GZKXfUvTB7RWG9AV3pQ19SfeqBoeSeMm+hzMMzajZNFHsODo3v45JnILYbqv1dX9+ir+rkXN5i5szqoa54HdbVA3gnjOnkW0k7WiW1bThZDv4RD3xh7gFfWoK7q1f3Zr4MFVzS+4p+9qFJFZKWXd8K4Tp6F1DqLk8XQcXjwo5Xl6qCuUyf/Jamr/gW8oqdP6kldNrXknTCuk2chtc7iM92s+XDVwaSu34O6rCLvs7U7eRZS6yxOFppWuYMwq5Hn2dqdPAtpp9m6ZyDz4mEqc7IwO0OdOgtpJ6n2eqveQXZir7dCO9dLWi1pv6QDktbX2T5L0l3J9p2SLkzKL5R0XNKjyetLRcZpVrRLP7OdSz+zvdVh2ASN1eutUxR2ZyGpC7gNeA/QD+yStC0i9qV2+xjwfES8XtI1wBeAq5NtT0bEZUXFZzaZxrqj+PTWvbz4cqUr9UUb7uXaKy7gP6y9tFmh2SRwr7di7ywuBw5ExMGIGATuBNbU7LMG+EqyfDfwLsktzDZ1fHrrXr724KFT68MRfO3BQ3x6694WRmXj1ah3Wyf1eisyWfQCz6TW+5OyuvtExBDwC2BBsm2JpD2SvifpbQXGaVaYO3Y+M65yK6d1q5bRPWP04NhO6/VWZAN3vTuE2ilMG+3zLLA4Ip6T9GZgq6Q3RMQLow6WrgeuB1i8ePEkhGw2uYaj/qy9jcqtnPJ2pZ7KikwW/cAFqfXzgcMN9umXNB04BzgaEQGcAIiIhyU9CVwM7E4fHBG3A7cD9PX1+b+vSfzshvy6pLqJocu1rW0nT1fqqazIaqhdwFJJSyTNBK4BttXssw24Lln+ILAjIkLSoqSBHEmvA5YCBwuM1XKqzrC686mjrNy4g617BlodUqlde8UF4yo3K6vC7iwiYkjSDcB2oAvYEhGPS7oF2B0R24AvA1+VdAA4SiWhALwduEXSEDAMfCIijhYVq+XjGVbHr9rr6Y6dzzAcQZfk3lDWlhRTpO60r68vdu/enb2jTdjKjTvqzoPU29PN99df2YKIzOxMSXo4Ivqy9vMTbyw39zU361xOFpab+5qbdS4nC8vNfc3NOpcnErTcPMOqWedysrBx8QyrZp3JycLGzUnCrPO4zaLFrt78wKkR0WZmZeVkYWZmmZwszMwsk5OFmZllcrIwM7NMThZmZpbJycLMzDI5WZiZWSYnCzMzy+RkYWZmmZwszMwsk5OFmZllcrIwM7NMThZmZpbJycLMzDI5WZiZWSYnCzMzy+RkYWZmmZwszMwsk5OFmZllcrKgdc/B3rpngD2HjrHzqaOs3LiDrXsGmh6DmVkeThYtsnXPABvu2cvg8AgAA8eOs+GevU4YZlZKhSYLSasl7Zd0QNL6OttnSbor2b5T0oWpbRuS8v2SVhUZZyts2r6f4yeHR5UdPznMpu37WxSRmVljhSULSV3AbcD7gOXAtZKW1+z2MeD5iHg98FfAF5JjlwPXAG8AVgNfTN5vyjh87Pi4ys3MWqnIO4vLgQMRcTAiBoE7gTU1+6wBvpIs3w28S5KS8jsj4kREPAUcSN5vyjivp3tc5WZmrVRksugFnkmt9ydldfeJiCHgF8CCnMci6XpJuyXtPnLkyCSGXrx1q5bRPWP0zVL3jC7WrVrWoojMzBorMlmoTlnk3CfPsUTE7RHRFxF9ixYtmkCIreuRtHZFL5//wKXM7Kp8Bb093Xz+A5eydsUrcqKZWctNL/C9+4ELUuvnA4cb7NMvaTpwDnA057FnrFGPJKApJ+21K3q546FDANz18d8o/PPMzCaqyDuLXcBSSUskzaTSYL2tZp9twHXJ8geBHRERSfk1SW+pJcBS4KHJDtA9kszM8insziIihiTdAGwHuoAtEfG4pFuA3RGxDfgy8FVJB6jcUVyTHPu4pK8D+4Ah4A8iYrjuB50B90gyM8unyGooIuJe4N6asptTyy8DH2pw7OeAzxUZ33k93QzUSQzukWRmNlpHj+B2jyQzs3wKvbMou2oj9o13P8bg8Ai9Pd2sW7XMPZLMzGp0dLIA90gyM8ujo6uhzMwsHycLMzPL5GRhZmaZnCzMzCxTxzdwgxu2zcyy+M7CzMwyOVmYmVkmJwszM8vkZGFmZpmcLMzMLJOThZmZZXKyMDOzTB5n0WIe42Fm7cB3FmZmlsnJwszMMjlZmJlZJicLMzPL5GRhZmaZnCzMzCyTk4WZmWVysjAzs0xOFmZmlkkR0eoYJoWkI8BPWh0HsBD4eauDaMCxTUxZYytrXODYJqoVsb02IhZl7TRlkkVZSNodEX2tjqMexzYxZY2trHGBY5uoMsfmaigzM8vkZGFmZpmcLCbf7a0OYAyObWLKGltZ4wLHNlGljc1tFmZmlsl3FmZmlsnJIgdJWyT9TNIPU2XzJd0v6Ynk57ykXJL+WtIBSY9JelPqmOuS/Z+QdF1BcX1I0uOSRiT11ey/IYlrv6RVqfLVSdkBSevPNK4xYtsk6f8mf5dvSuopUWyfTeJ6VNK3JZ2XlDft+2wUW2rbpySFpIVliU3SZyQNJH+3RyW9P7Wtpd9pUv6HyWc9LunWssQm6a7U3+xpSY+2IrZxiQi/Ml7A24E3AT9Mld0KrE+W1wNfSJbfD9wHCHgrsDMpnw8cTH7OS5bnFRDXJcAy4LtAX6p8OfADYBawBHgS6EpeTwKvA2Ym+ywv6G/2XmB6svyF1N+sDLHNTS3/W+BLzf4+G8WWlF8AbKcylmhhWWIDPgN8qs6+ZfhOfwv4B2BWsv7qssRWs/0vgJtbEdt4Xr6zyCEi/hdwtKZ4DfCVZPkrwNpU+d9GxYNAj6RzgVXA/RFxNCKeB+4HVk92XBHxo4jYX2f3NcCdEXEiIp4CDgCXJ68DEXEwIgaBO5N9z0iD2L4dEUPJ6oPA+SWK7YXU6llAtTGvad9no9gSfwXcmIqrTLHV0/LvFPgksDEiTiT7/KxEsQGVu0PgXwJ3tCK28XCymLhfi4hnAZKfr07Ke4FnUvv1J2WNypulbHH9PpWr4tLEJulzkp4BPgzcXJbYJF0FDETED2o2tTy2xA1JNdgWJdWxJYntYuBtknZK+p6kt5Qotqq3AT+NiCdKGNsoThaTT3XKYozyZilNXJJuAoaA/1otahBDU2OLiJsi4oIkrhvKEJukVwE3cTp5jdrcIIZm/t3+M3ARcBnwLJUqFcaIoZmxTadSDfdWYB3w9eRKvgyxVV3L6bsKxoih1ecPJ4sz8NPklp/kZ/UWt59K/XLV+cDhMcqbpRRxJY2tvwN8OJJK2rLElvJ3wL8oSWwXUam7/oGkp5PPeUTSa0oQGxHx04gYjogR4G+oVJdQhtiSz7onqaZ7CBihMvdSGWJD0nTgA8BdNTG3PLa6mtlA0s4v4EJGN55tYnQD963J8m8zutHxoaR8PvAUlSudecny/MmOK1X+XUY3cL+B0Q1nB6k0mk1PlpdwuuHsDQX9zVYD+4BFNfuVIbalqeU/BO5uxfc51neabHua0w3cLY8NODe1/EdU6tvL8p1+ArglWb6YSjWOyhBb6v/he63+X8j9OzTzw9r1ReU28VngJJUM/zFgAfAd4Ink5/xkXwG3Uem5sJfRJ+zfp9JgdQD4aEFx/W6yfAL4KbA9tf9NSVz7gfelyt8P/DjZdlOBf7MDyT/so8nrSyWK7b8BPwQeA/470Nvs77NRbDXbn+Z0smh5bMBXk89+DNjG6OTR6u90JvC15Ht9BLiyLLEl5f8F+ESd/ZsW23heHsFtZmaZ3GZhZmaZnCzMzCyTk4WZmWVysjAzs0xOFmZmlsnJwszMMjlZmJlZpumtDsBsKpFUHawJ8BpgGDiSrF8elRlDzdqOB+WZFUTSZ4CXIuLPWx2L2ZlyNZSZmWVysjAzs0xOFmZmlskN3GZNkjy/4FYqD635SUT8dYtDMsvNycKseT4JfCsivtfqQMzGy9VQZs3zJuD7rQ7CbCLcddasSSStAa4CjgKfj4ijLQ7JLDcnCzMzy+RqKDMzy+RkYWZmmZwszMwsk5OFmZllcrIwM7NMThZmZpbJycLMzDI5WZiZWSYnCzMzy/T/AYao0Zs9GsXdAAAAAElFTkSuQmCC\n", | |
"text/plain": [ | |
"<Figure size 432x288 with 1 Axes>" | |
] | |
}, | |
"metadata": { | |
"needs_background": "light" | |
}, | |
"output_type": "display_data" | |
} | |
], | |
"source": [ | |
"plt.errorbar(tc, ab, err, fmt='o', ls='')\n", | |
"xx = np.linspace(min(tc), max(tc), 100)\n", | |
"plt.plot(xx, soln['x'][0]*xx + soln['x'][1], c='orange')\n", | |
"plt.xlabel(r'T$_c$')\n", | |
"plt.ylabel(r'[X/H]');" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"#### checking with linear algebra:" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 11, | |
"metadata": {}, | |
"outputs": [], | |
"source": [ | |
"Y = np.copy(ab)\n", | |
"A = np.ones((len(ab),2))\n", | |
"A[:,0] = tc\n", | |
"Cinv = np.diag(1./err**2) # sorry Hogg" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 12, | |
"metadata": {}, | |
"outputs": [], | |
"source": [ | |
"ACA = np.dot(A.T, np.dot(Cinv, A))\n", | |
"ACY = np.dot(A.T, np.dot(Cinv, Y))\n", | |
"soln = np.linalg.solve(ACA, ACY)" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 13, | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"name": "stdout", | |
"output_type": "stream", | |
"text": [ | |
"best-fit slope is: 4.23e-05\n" | |
] | |
} | |
], | |
"source": [ | |
"print('best-fit slope is: {0:.2e}'.format(soln[0]))" | |
] | |
} | |
], | |
"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.7.1" | |
} | |
}, | |
"nbformat": 4, | |
"nbformat_minor": 2 | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment