{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "In this notebook, we investigate linear least squares fitting" ] }, { "cell_type": "code", "execution_count": 78, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib\n", "%matplotlib inline\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's define a linear function that is a polynomial, where linear means linear in the parameters, not the independent variable. Write a function to evaluate a general polynomial, where the order of the polynomial is determined by the length of the parameter array that you pass it.\n" ] }, { "cell_type": "code", "execution_count": 94, "metadata": {}, "outputs": [], "source": [ "def func(x,par) :\n", " \"\"\" Function to return a polynomial evaluates at x\n", " y = par[0] + par[1]*x + par[2]*x**2 + par3*x**3 + ....\n", " where the number of terms (order of polynomial) is determined by the lenght\n", " of the parameter array (e.g. if you give it a 2-parameter array, that would\n", " be a straight line)\n", " \"\"\"\n", " n=len(x)\n", " m=len(par)\n", " y=np.zeros(n) + par[0]\n", " for i in range(1,m) :\n", " y+= par[i]*x**i\n", " return y" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's generate a data set. Start with a straight line" ] }, { "cell_type": "code", "execution_count": 100, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 100, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAMn0lEQVR4nO3dXYhd1RnG8edpoilTFS8yIJhkxtJSDFawDGIR6qC2+BH0phfao0gV5kaLgmK1c50rwVpQKhNtKTggxQ8sovWDmkAvIk40ChoVERPjB46UonQgEnx7sU/MnDgxc2avmb3fc/4/GJi15rD2yyY8rKyz116OCAEA8vpe0wUAAOohyAEgOYIcAJIjyAEgOYIcAJJb38RFN27cGOPj401cGgDS2rNnz+cRMXpsfyNBPj4+rrm5uSYuDQBp2d6/VD9LKwCQHEEOAMkR5ACQHEEOAMkR5ACQHEEOAMkR5ACQHEEOAMkR5ACwViYnq5/CCHIASI4gB4DkCHIASI4gB4DkCHIASI4gB4DkCHIASI4gB4DkCHIASI4gB4DkCHIASI4gB4DkCHIASI4gB4DkCHIASI4gB4DkCHIASK5IkNs+3fZjtt+2vc/2z0uMCwADY3ZW2r1b2rVLGh+v2oWsLzTOnyT9MyJ+bftkSSOFxgWA/GZnpakp6dChqr1/f9WWpE6n9vC1Z+S2T5P0C0kPS1JEfBUR/607LgAMjOlpaWGht29hoeovoMTSyg8lzUv6q+3XbD9k+wfHfsj2lO0523Pz8/MFLgsASRw40F9/n0oE+XpJP5P054g4T9L/JN117IciYiYiJiJiYnR0tMBlASCJLVv66+9TiSA/KOlgRLzcbT+mKtgBAJK0fbs0csxXhyMjVX8BtYM8Ij6V9KHtn3S7LpH0Vt1xAWBgdDrSzIy0YUPVHhur2gW+6JTKPbXyO0mz3SdW3pf020LjAsBg6HSkHTuq33fuLDp0kSCPiL2SJkqMBQDoDzs7ASA5ghwAkiPIASA5ghwAkiPIASA5ghwAkiPIASA5ghwAkiPIASA5ghwAkiPIASA5ghwAkiPIASA5ghwAkiPIASC5UgdLAABOpPCBEkcwIweA5AhyAEiOIAeA5AhyAEiOIAeA5AhyAKtncrL6waoiyAEgOYIcAJIjyAEgOYIcwOAb8LV6ghwAkiPIASA5ghwAkiPIASA5ghwAkiPIASA5ghwYRAP+uB16FQty2+tsv2b76VJjAgBOrOSM/FZJ+wqOBwBYhiJBbnuTpCslPVRiPCAtljTQgFIz8vsk3Snp6+N9wPaU7Tnbc/Pz84UuCwCoHeS2t0n6LCL2fNfnImImIiYiYmJ0dLTuZQEAXSVm5BdKusr2B5IelXSx7UcKjAsAWIbaQR4Rd0fEpogYl3SNpH9FxHW1KwOQ2+ystHu3tGuXND5etbEqeI4cQHmzs9LUlHToUNXev79qE+aromiQR8TOiNhWckwACU1PSwsLvX0LC1U/imNGDqC8Awf660ctBDmA8rZs6a9/NQ3BWj1BDqC87dulkZHevpGRqn8tDclaPUEOoLxOR5qZkTZsqNpjY1W701nbOoZkrX590wUAGFCdjrRjR/X7zp3N1DAka/XMyAEMrjat1a8ighzA4GrLWv0qI8gBDK62rNWvMtbIgUFz5HG7Q4eqx+22bx+44OpLG9bqVxkzcmCQDMnjduhFkAODZEget0MvghwYJEPyuB16EeQYDG04Yq0NW8GH5HE79CLIgRLasjY9JI/boRdBDpTQlrXpIXncDr14/BAooU1r00PwuB16MSMHSmBtGg0iyIESWJtGgwhyoATWptEg1siBUlibRkOYkQNAcszIAawe/meyJpiRA0ByBDkAJEeQA0ByrJGjviMvq2I9FG014P82mZEDQHIEOQAkR5ADQHIEOQAkR5AjvzaczAM0iCBHbm05mQdoEEGO3NpyMk/b7Nw58I/c4ajaQW57s+2XbO+z/abtW0sUBixLm07mARpSYkZ+WNLtEXG2pAsk3Wx7a4FxgRPjZB6gfpBHxCcR8Wr39y8l7ZN0Zt1xgWXhZB6g7BZ92+OSzpP0cslxgeM6cgLPTTdVX3iOjVUh3tTJPKxLowHFgtz2KZIel3RbRHyxxN+nJE1J0hb+24uSOJkHQ67IUyu2T1IV4rMR8cRSn4mImYiYiIiJ0dHREpcFAKjMUyuW9LCkfRFxb/2SAAD9KDEjv1DS9ZIutr23+3NFgXEBAMtQe408Iv4tyQVqAQCsADs7M5ucPHqoA4ChRZCjHl5YBTSOIMfK8cIqoBUIcqwcL6wCWoEgx8rxwiqgFQhyrBwvrAJagSDHyvHCKqAVCHKsXKcjzcxIGzZU7bGxqt3UC6uAIVX07YcYQrywCmgcM3IASI4gB4DkCPKVYGs8gBZhjRyDgfV5DDFm5ACQHEEOAMkR5ACQHEEOAMkR5ACQHEGeFQc6AOgiyDPiQAcAixDkGXGgA4BFCPKMONABwCIEeUYc6ABgEYI8Iw50ALAIQZ4RBzoAWISXZmXVpgMdmr4+MOSYkQNAcgQ5ACRHkANAcgR5v9gaD6BlCPJ+sDUeQAsR5P1gazyAFiLI+8HWeAAtRJD3g63xAFqoSJDbvsz2O7bfs31XiTFbia3xAFqodpDbXifpAUmXS9oq6VrbW+uO20psjQfQQiW26J8v6b2IeF+SbD8q6WpJbxUYu33atDUeAFRmaeVMSR8uah/s9vWwPWV7zvbc/Px8gcsCAKQyQe4l+uJbHREzETEREROjo6MFLgsAkMoE+UFJmxe1N0n6uMC4AIBlKBHkr0j6se2zbJ8s6RpJ/ygwLgBgGWp/2RkRh23fIuk5Sesk/SUi3qxdGQBgWYocLBERz0h6psRY6ANPzQAQOzsBID2CHACSyxXkk5PVDwDgG7mCHADwLQQ5ACRHkANAcgQ5ACRHkANAcgQ5ACRHkANAckW26A8dtsYDaBFm5ACQHEEOAMkR5ACQHEEOAMkR5ACQHEEOAMkR5ACQHEEOAMkR5ACQXJ4gn52Vdu+Wdu2SxserNgAgSZDPzkpTU9KhQ1V7//6qTZgDQJIgn56WFhZ6+xYWqn4AGHI5gvzAgf76AWCI5AjyLVv66weAIZIjyLdvl0ZGevtGRqp+ABhyOYK805FmZqQNG6r22FjV7nSarQsAWiDPwRKdjrRjR/U7BzsAwDdyzMgBAMdFkANAcgQ5ACRHkANAcgQ5ACRXK8ht32P7bdtv2H7S9umlCgMALE/dGfkLks6JiHMlvSvp7volAQD6USvII+L5iDjcbe6WtKl+SQCAfpRcI79R0rPH+6PtKdtztufm5+cLXhYAhtsJd3baflHSGUv8aToinup+ZlrSYUnHfUF4RMxImpGkiYmJWFG1AIBvOWGQR8Sl3/V32zdI2ibpkoggoAFgjdV614rtyyT9XtJFEbFwos8DAMqru0Z+v6RTJb1ge6/tBwvUBADoQ60ZeUT8qFQhAICVYWcnACRHkANAcnkOlpA4UAIAlsCMHACSI8gBIDmCHACSI8gBIDmCHACSI8gBIDmCHACSI8gBIDmCHACScxOvELc9L2n/ml+4rI2SPm+6iBbhfhzFvejF/ehV536MRcTosZ2NBPkgsD0XERNN19EW3I+juBe9uB+9VuN+sLQCAMkR5ACQHEG+cjNNF9Ay3I+juBe9uB+9it8P1sgBIDlm5ACQHEEOAMkR5AXYvsN22N7YdC1NsX2P7bdtv2H7SdunN11TE2xfZvsd2+/Zvqvpeppie7Ptl2zvs/2m7VubrqkNbK+z/Zrtp0uOS5DXZHuzpF9KOtB0LQ17QdI5EXGupHcl3d1wPWvO9jpJD0i6XNJWSdfa3tpsVY05LOn2iDhb0gWSbh7ie7HYrZL2lR6UIK/vj5LulDTU3xpHxPMRcbjb3C1pU5P1NOR8Se9FxPsR8ZWkRyVd3XBNjYiITyLi1e7vX6oKrzObrapZtjdJulLSQ6XHJshrsH2VpI8i4vWma2mZGyU923QRDThT0oeL2gc15OElSbbHJZ0n6eVmK2ncfaomfV+XHnh96QEHje0XJZ2xxJ+mJf1B0q/WtqLmfNe9iIinup+ZVvXf6tm1rK0lvETfUP9PzfYpkh6XdFtEfNF0PU2xvU3SZxGxx/Zk6fEJ8hOIiEuX6rf9U0lnSXrdtlQtJbxq+/yI+HQNS1wzx7sXR9i+QdI2SZfEcG5QOChp86L2JkkfN1RL42yfpCrEZyPiiabradiFkq6yfYWk70s6zfYjEXFdicHZEFSI7Q8kTUTEUL7lzfZlku6VdFFEzDddTxNsr1f1Re8lkj6S9Iqk30TEm40W1gBXs5u/SfpPRNzWdD1t0p2R3xER20qNyRo5Srlf0qmSXrC91/aDTRe01rpf9t4i6TlVX+79fRhDvOtCSddLurj772FvdzaKVcCMHACSY0YOAMkR5ACQHEEOAMkR5ACQHEEOAMkR5ACQHEEOAMn9H+4HWM1kJSoLAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# set true parameters\n", "npar=2\n", "truepar=[2.6,0.7] # intercept and slope\n", "\n", "# generate a data set with true underlying relation plus gaussian noise\n", "n=10 # number of points to generate\n", "x=np.arange(0,10,1)-5 # independent variable\n", "sig=1*np.ones(n) # uncertainties on the dependent variable\n", "\n", "# generate the data set using your function + noise\n", "y=func(x,truepar)+np.random.normal(0.,sig)\n", "plt.errorbar(x,y,sig,fmt='ro')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Linear least squares using relations specific to a straight line function, as derived in class:\n", "$$a = {S_{xx} S_y - S_xS_{xy} \\over S S_{xx} - S_x^2}$$\n", "$$b = {S S_{xy} - S_xS_y \\over S S_{xx} - S_x^2}$$\n", "where\n", "$$S = \\sum {1\\over \\sigma_i^2}$$\n", "$$S_x = \\sum {x_i\\over \\sigma_i^2}$$\n", "$$S_{xy} = \\sum {x_i y_i\\over \\sigma_i^2}$$\n", "$$S_{xx} = \\sum {x_i x_i\\over \\sigma_i^2}$$\n", "$$S_y = \\sum {y_i\\over \\sigma_i^2}$$\n", "\n" ] }, { "cell_type": "code", "execution_count": 110, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "par [2.6167978 0.66856591]\n" ] }, { "data": { "text/plain": [ "[]" ] }, "execution_count": 110, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAVRUlEQVR4nO3dfXBV9Z3H8c+XB9EIaC0ICiRRwQewnQVTJHV2pVUUH+pTO101pE86TGdqV1ttq9LZ1rZMu7WPs7rtBLudhmRH0bZTtPZBK9J2J2ENQrA+1FrLTQK0IE8qMQkk3/3jJJiEG5Kbe5Jzf/e+XzMOOSfXc39zhvl4/H3O7xxzdwEAwjUm6QEAALJDkANA4AhyAAgcQQ4AgSPIASBw45L40ilTpnhpaWkSXw0Awdq4ceNr7j61//5Egry0tFQNDQ1JfDUABMvMUun2M7UCAIEjyAEgcAQ5AASOIAeAwBHkABA4ghwAAkeQA0DgCHIACBxBDgCjZfHi6J+YEeQAEDiCHAACR5ADQOAIcgAIHEEOAIEjyAEgcAQ5AASOIAeAwBHkABA4ghwAAkeQA0DgCHIACBxBDgCBI8gBIHAEOQAEjiAHgMAR5AAQuFiC3MxONLNHzOwlM3vRzMrjOC4A5I3aWqm+Xlq/XiotjbZjMi6m43xf0q/d/UNmdoykopiOCwDhq62Vli+X2tuj7VQq2pakioqsD5/1FbmZTZb0L5J+JEnu3uHu+7I9LgDkjRUrpNbWvvtaW6P9MYhjauV0Sbsk/djMNpnZA2Z2fP8PmdlyM2sws4Zdu3bF8LUAEIimpsz2ZyiOIB8naYGkH7j7fEkHJN3Z/0PuXuXuZe5eNnXq1Bi+FgACUVyc2f4MxRHkLZJa3H1D9/YjioIdACBJK1dKRf2qw6KiaH8Msg5yd/+7pGYzO6t710WSXsj2uACQNyoqpKoqacKEaLukJNqOoeiU4rtr5dOSarvvWHlV0sdjOi4A5IeKCmnVqujnp5+O9dCxBLm7b5ZUFsexAACZYWUnAASOIAeAwBHkABA4ghwAAkeQA0DgCHIACBxBDgCBI8gBIHAEOQCMMHfXxtReffaMy7Vz/BEPh81aXEv0AQD9vNXRqbWN21Rdl9Lz21/XxHfM1lWvvaCTY/4eghwAYva31w6opj6lhxua9XrbIZ01bZK+ds25uubumzWx62Ds30eQA0AMOrtcT720U9V1W/WHv7ymcWNMS8+drspFJVp42kkyM2kEQlwiyAEgK7vfbNeDzzTrfzY0adu+tzRt8gR95uIzdcPCWTp58rGjMgaCHAAy5O7a1LxPq+tS+uWWHero7FL56e/UF684RxfPnabxY0f3PhKCHACG6IjycsI43bBwlpYtKtGcaZMSGxdBDgCD6F9enjltor56zbm6dv4MTZyQQYzG/EKJHgQ5AKTRU16urk/p9y/v0rgxpkvPna6P9C4vcwRBDgC97H6zXQ81NKu2PrnyMlMEOYCCN1B5ueKKc7QkgfIyUwQ5gIKVrry8fuEsVSZcXmaKIAdQcLb2lJcbW7T/rYPDLy9zRHgjBhCOxYujP0fobo1MdHa51r20U9UBlJeZIsgB5LV05eVtF8/RDQuLNS1Hy8tMEeQA8k668nLR6ScFU15miiAHkDfe6ujUo43bVV2/VX/aFm55mSmCHEDwBi0vc2iufiQQ5ACClLa8nDddleUlOj/w8jJTBDmAoBRCeZkpghxAzuspL2vqUnqsAMrLTBHkAHJW//Ly+GPG6vrux8aemcflZaYIcgA5p395Oefkifrq1fN07YKZQa68HGmcEQA5gfJy+AhyIB8FdLtd//Ly5EmUl5mKLcjNbKykBknb3P3KuI4LIP+4uzZ3r7zsXV7effk5umQe5WWm4rwiv1XSi5Imx3hMAHmE8nJkxBLkZjZT0hWSVkr6bBzHBIIU0JTGaKK8HFlxncHvSfq8pAH/k2pmyyUtl6Ti4uKYvhZArurscq078Qytnv5PWv+tpykvR1DWQW5mV0ra6e4bzWzxQJ9z9ypJVZJUVlbm2X4vgNy0+812rWloUe2GlFrOvk4nd7xJeTnC4rgiv0DSVWZ2uaRjJU02sxp3XxbDsQEEoE95+dwOdRyKysu7NjyoS/a+ovHfeSrpIea1rIPc3e+SdJckdV+R30GIA4UhXXn5r2WzVFleojOfXCs9Xi21t0ulpdLKlVJFRdJDzku0DAAyNmh5WVsrLV8ehbgkpVLRtkSYj4BYg9zdn5b0dJzHBJAbOrtcT/95p6rrUlo/2MrLFSuk1ta+B2htjfYT5LHjihzAUe050KGHnmmOysu9Q1x52dSU2X5khSAHcIR05eX5p52kuy4b4srL4uJoOiXd/tFWWyvV1+f1XD1BDuCwtoOdWrt5u1bXp/Tctv19y8tMVl6uXBnNifeeXikqivaPpgKZqyfIASi1Oyov1zT0LS+vmT9Dk44dn/kBe0LyppuiEC0pSeZKuEDm6glyoED1Ly/HjjEtnTddyxaVaNHpMay8rKiQVq2Kfk7qkQUFMldPkAMFJl15eetFc3Tj+Xm48jKX5upHEEEOFIDD5WV992NjMy0vQ5Urc/UjjCAH8ljbwU6tbdyu1XV9y8tli0p01vQCeGxsrszVjzCCHMg3tbU6tOH/NLbtLe2Zeqr++M+VanvfB/SVq+fp2uGWlyHLhbn6EUaQA3mis8v14rd/oDlfvF0TOtokSafu36nvPvlfGnP9fFn5hQmPECMlTyfGgMKx50CHfrj+r7rw3nU68WtfPhziPca2vSVbsSKZwWFUcEUOBMjd1diyX9V1W/uUlzPeeC39v5Bnt9uhL4Ic+SEXXrE2CkvB05WXHy6bqcpFpVF5+fXCuN0OfRHkQBxGeCl4/5WXs0+emL68LJDb7dAXQQ7EYQSWgvesvFxdH628HGOmS+dNU+Wi0oFXXhbI7XboiyAH4hDjUvA9Bzq0pqFZNfVvr7z8t/dHj42dfsIQVl4WwO126IsgB+IQw1Lwzc37jigv77zsbF06b3r+rrxELAhyIA7DnJvuKS9r6lPa0pKmvASGgCAH4pDh3PSQy0tgCAhyIC6DzE13drnWv/z2Y2OHVF4CQ0CQAyOsp7ys3ZBS855hlJfAIAhyYISkKy+/sLTAykvumhkVBDkQozYbp7VTzlbNfX+kvMSoIciBGKR2H1DthiatWfBJ7Rt/nGZ3dFJeYtQQ5MAwpS0vX2/Ssn9sUvnXqykvMWoIcmQvFx5YNYr6l5dTJ03Qp98/RzcuLNb0q78ZfYgQzy15/neTIAeGqH95ubAQy0vkJIIcOIq2g516tHG7VrPyEjmMIAfSaNrdqpoNKa1paNa+VlZeIrcR5EC3dOXlJXOnqbK8ROWnv5PyEjmLIEf4snwzz94DHXpooPKSlZcIAEGOsGXxZp7G5n2qrkvp0S3bD5eXn780Ki+PGUd5iXAQ5Ahbhm/m6V9eFnWXl8sWlejs6ZNHadCjIM9vt0NfWQe5mc2SVC1puqQuSVXu/v1sjwsMyRDfzNO0u1W1G1J6qFd5ec9V83TdAspLhC+OK/JDkm5392fNbJKkjWb2hLu/EMOxgaM7ypt5Ortcv395l6rrtuppykvksayD3N13SNrR/fMbZvaipBmSCHKMvDRv5vHjivS7Zbfqnm+tU/OetzRl4gR9+n2zdcP5xTrlhOMSHCwwMmKdIzezUknzJW2I87jAgHq9mcfb27V3yin62gXL9LNDZ2rhCceNfnnJ3DQSEFuQm9lEST+VdJu7v57m98slLZek4gxeSAscTdvBTj161r+o5ub71DjxFBUdM1bXzp+hX5fnWXkJHEUsQW5m4xWFeK27/yzdZ9y9SlKVJJWVlXkc34vC1b+8PGPsMfry357UdT/5piZTXqLAxHHXikn6kaQX3f072Q8JSK+ry7U+XXm5qETlN39IJkmEOApQHFfkF0iqlPScmW3u3ne3uz8ew7EB7e1+bGxN98pLykugrzjuWvmjJO7jQuwam/dpdX1KjzZuV/uhLi0sZeUlkA4rO0OWhy906Fl5WVOfUmP3yssPnTdTlZSXwIAIcmQnywdW9Wje06qa+l7l5dTjdc9V83TtghmUl8AgCHIMXxYPrJIGKS/PYOUlMFQEOYYvwwdW9dh7+J2XTWra00p5CWSJIMfwDfGBVT3SlZefu/QsyksgSwQ5hu8oD6zq0XawU49t2aHVdVspL4ERQpBj+NI8sEpFRdLKlYfLyzUNzdpLeQmMKIIcw9frgVVqb5cXF+uFW+7Utw+eqXX3rqO8BEYJQY7sVFRo749X6+Gp71LN/MvVtLtVU9r3U14Co4ggx7Btael+5+WCT6p9zHgtnHws5SWQAIIcGTlcXtan1Ni8Lyovdz2vyr9v0tnfeCTp4QEFiSAfjjxcGj+Yo5aXS7+R9PCAgkaQY0BdXa71f9ml1XUprfvzztwuLwvoP6pAfwQ5jrCvtfuxsfWsvARCQJDjsMPlJSsvgaAQ5AUubXnJyksgKAR5gWre06qaDSmteYaVl0DoCPICElR5CWDICPJQZfBCB8pLIL8R5CEa4gsdtrTs0+q6lNZSXgJ5jSAP0VFe6ND24ev1yy07VE15CRQMgjxEA7y4wZuaVP7131FeAgWGIA/RAC902D5pihad/k7KS6DAMFEaoANfukcHJxzbZ9/BCceq6Nv/oR8sO0/vnT2FEAcKCEEekOda9utzDzdqwV9P1u1LPqVdE0+SS/LiYo3/0QN6x80fT3qIABLA1EqOazvYmb68vOXfNfX6huhDST8wKunvBwocQZ6j0q28/PIH5uq682ZSXgLogyDPIT0rL2vqUnqKlZcAhoggzwH7Wjv0cEOLajaklNrNyksAmSHIM5XB0vjBPNeyX9V1W/usvLzjElZeAsgMQZ6JIS6NP5p05eUHz5upykUlOucUVl4CyBxBnomjLI0fLMib97SqdkOT1jQ0a8+BDspLALEhyDMxwNL4gfZ3dbl+3/3Y2J7ycsk50/SRcspLAPEhyDMxwNJ4FRf32aS8BDCaYglyM1sq6fuSxkp6wN2/Ecdxc87KldGceO/plaKiaL/Sl5e3X3KWllJeAhhBWQe5mY2VdL+kJZJaJD1jZmvd/YVsj51zeubBb7opKjxLStTxla/qsXMuVPX9/6vNlJcAEhDHFflCSa+4+6uSZGYPSrpaUv4FuRSF+apVap4wWbW3fTMqL19opLwEkJg4gnyGpOZe2y2Szu//ITNbLmm5JBX3m1MOxeHy8qxr9dSJZ2jMH16lvASQuDiCPF16+RE73KskVUlSWVnZEb/PZUeUl8efolu21enG+79IeQkgcXEEeYukWb22Z0raHsNxE/dcy36trt+qX2yOysv3lL4jKi9vW6ZjvEsixAHkgDiC/BlJc8zsNEnbJF0v6cYYjpuItoOdevy5HaquSw1cXnpXsoMEgF6yDnJ3P2Rmt0j6jaLbD//b3Z/PemSjrP/Ky9OnHq8vfWCuPkh5CSDHxXIfubs/LunxOI41mnrKy5r6lH730tsrLyvLS/TeEMpLXugAQAW6snNfa4ce2diimvqUtnavvLzlfbN1IysvAQSooIL8T9uilZe9y8vPsvISQODCCvLFi6M/M5hS6F9eHjc+Ki+XnV+iuaey8hJA+MIK8gxQXgIoFHkV5P3LS5O0ZO40faS8NIzyEgCGIS+CfKDy8oaFxTr1RMpLAPkt6CDvKS/XNm5X28GovPzMkjN12bmnUF4CKBjBBXmbjdXjz7ZodX1Km5qi8vLa+dHKS8pLAIUoqCD/8fQF+s8Z5dqzpvFweXndgpk64TjKSwCFK6ggN3e9540WVX6uUhfMprwEACmwIP/YPzbpY//YJM35TLIDYWk8gBxCIwgAgSPIASBwBDkABI4gB4DAEeQAEDiCHAACR5ADQOAIcgAIHEEOAIELJ8hra6X6emn9eqm0NNoGAAQS5LW10vLlUnt7tJ1KRduEOQAEEuQrVkitrX33tbZG+wGgwIUR5E1Nme0HgAISRpAXF2e2HwAKSBhBvnKlVFTUd19RUbQfAApcGEFeUSFVVUkTJkTbJSXRdkVFsuMCgBwQzoslKiqkVauin3mxAwAcFsYVOQBgQAQ5AASOIAeAwBHkABA4ghwAApdVkJvZvWb2kpltMbOfm9mJcQ0MADA02V6RPyHpXHd/t6SXJd2V/ZAAAJnIKsjd/bfufqh7s17SzOyHBADIRJxz5J+Q9KuBfmlmy82swcwadu3aFePXAkBhG3Rlp5k9KWl6ml+tcPdfdH9mhaRDkgZ8QLi7V0mqkqSysjIf1mgBAEcYNMjd/eKj/d7MPirpSkkXuTsBDQCjLKtnrZjZUklfkHShu7cO9nkAQPyynSO/T9IkSU+Y2WYz+2EMYwIAZCCrK3J3nx3XQAAAw8PKTgAIHEEOAIEL58USEi+UAIA0uCIHgMAR5AAQOIIcAAJHkANA4AhyAAgcQQ4AgSPIASBwBDkABI4gB4DAWRKPEDezXZJSo/7F8Zoi6bWkB5FDOB9v41z0xfnoK5vzUeLuU/vvTCTI84GZNbh7WdLjyBWcj7dxLvrifPQ1EueDqRUACBxBDgCBI8iHryrpAeQYzsfbOBd9cT76iv18MEcOAIHjihwAAkeQA0DgCPIYmNkdZuZmNiXpsSTFzO41s5fMbIuZ/dzMTkx6TEkws6Vm9mcze8XM7kx6PEkxs1lmts7MXjSz583s1qTHlAvMbKyZbTKzx+I8LkGeJTObJWmJpKakx5KwJySd6+7vlvSypLsSHs+oM7Oxku6XdJmkuZJuMLO5yY4qMYck3e7u50haJOlTBXwuertV0otxH5Qgz953JX1eUkG3xu7+W3c/1L1ZL2lmkuNJyEJJr7j7q+7eIelBSVcnPKZEuPsOd3+2++c3FIXXjGRHlSwzmynpCkkPxH1sgjwLZnaVpG3u3pj0WHLMJyT9KulBJGCGpOZe2y0q8PCSJDMrlTRf0oZkR5K47ym66OuK+8Dj4j5gvjGzJyVNT/OrFZLulnTJ6I4oOUc7F+7+i+7PrFD0v9W1ozm2HGFp9hX0/6mZ2URJP5V0m7u/nvR4kmJmV0ra6e4bzWxx3McnyAfh7hen229m75J0mqRGM5OiqYRnzWyhu/99FIc4agY6Fz3M7KOSrpR0kRfmAoUWSbN6bc+UtD2hsSTOzMYrCvFad/9Z0uNJ2AWSrjKzyyUdK2mymdW4+7I4Ds6CoJiY2VZJZe5ekE95M7Olkr4j6UJ335X0eJJgZuMUFb0XSdom6RlJN7r784kOLAEWXd38RNIed78t6fHkku4r8jvc/cq4jskcOeJyn6RJkp4ws81m9sOkBzTausveWyT9RlG5t6YQQ7zbBZIqJb2/++/D5u6rUYwArsgBIHBckQNA4AhyAAgcQQ4AgSPIASBwBDkABI4gB4DAEeQAELj/B+0kbGjBU2n0AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "#sums\n", "s=(1./sig**2).sum()\n", "sx=(x/sig**2).sum()\n", "sxy=(x*y/sig**2).sum()\n", "sxx=(x*x/sig**2).sum()\n", "sy=(y/sig**2).sum()\n", "\n", "# derived parameters\n", "par=np.zeros(npar)\n", "par[0]=(sxx*sy-sx*sxy)/(s*sxx-sx**2)\n", "par[1]=(s*sxy-sx*sy)/(s*sxx-sx**2)\n", "print('par',par)\n", "\n", "#plot them\n", "plt.errorbar(x,y,sig,fmt='ro')\n", "plt.plot(x,func(x,par))\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "Uncertainties in the parameters:\n", "$$\\sigma_a^2 = {S_{xx}\\over S S_{xx} - S_x^2}$$\n", "$$\\sigma_b^2 = {S \\over S S_{xx} - S_x^2}$$\n", "$$ \\sigma_{ab} = -\\sigma_b^2 \\bar{x} $$\n" ] }, { "cell_type": "code", "execution_count": 111, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "sigpar [0.32098334 0.11009638]\n", "covariance -0.00606060606060606\n" ] } ], "source": [ "#uncertainties\n", "sigpar=np.zeros(npar)\n", "sigpar[0]=np.sqrt(sxx/(s*sxx-sx**2))\n", "sigpar[1]=np.sqrt(s/(s*sxx-sx**2))\n", "print('sigpar',sigpar)\n", "print('covariance',sigpar[1]**2*x.mean())\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "How do the derived parameters compare with the true parameters? What happens as you change the size of the uncertainties and/or the number of data points?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's use the more general formalism for linear least squares. Start by writing a function that calculates the derivatives of the function with respect to each parameter, which, for a polynomial, is just 1, $x$, $x^2$, ... Note that this special case of derivatives with respect to the parameters of a polynomial is called a Vandermonde matrix, and numpy has a function that calculates it if you are lazy, \n", " numpy.vander() \n" ] }, { "cell_type": "code", "execution_count": 120, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[ 1. -5.]\n", " [ 1. -4.]\n", " [ 1. -3.]\n", " [ 1. -2.]\n", " [ 1. -1.]\n", " [ 1. 0.]\n", " [ 1. 1.]\n", " [ 1. 2.]\n", " [ 1. 3.]\n", " [ 1. 4.]]\n" ] } ], "source": [ "def deriv(x,npar) :\n", " \"\"\" Return derivatives with respect to each parameter for a polynomial with npar parameters\n", " returns N rows by npar cols\n", " \"\"\"\n", " n=len(x)\n", " d=np.zeros([n,npar])\n", " d[:,0] = 1\n", " for i in range(1,npar) :\n", " d[:,i] = x**i\n", " return d #np.array([np.ones(len(x)),x])\n", "\n", "print(deriv(x,2))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To do the linear least squares, you can caculate the $\\alpha$ matrix and the $\\beta$ vector:\n", "\n", "$$\\alpha_{jk} = \\sum_{i=0}^{N-1} {f_j(x_i) f_k(x_i)\\over \\sigma_i^2}$$\n", "$$\\beta_j = \\sum_{i=0}^{N-1} {f_j(x_i) y_i \\over \\sigma_i^2}$$\n", "Your deriv() routine supplies the $f_j$" ] }, { "cell_type": "code", "execution_count": 121, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "alpha:\n", "[[10. -5.]\n", " [-5. 85.]]\n", "beta\n", "[22.82514847 43.7441136 ]\n" ] } ], "source": [ "# sum over the data points. Note that this can be done better with \n", "# array arithmetic, but I write the for loops to be clear\n", "beta=np.zeros(npar)\n", "alpha=np.zeros([npar,npar])\n", "for k in np.arange(npar) :\n", " beta[k] = (deriv(x,2)[:,k]*y/sig**2).sum()\n", " for j in np.arange(npar) :\n", " alpha[j,k] = (deriv(x,2)[:,j]*deriv(x,2)[:,k]/sig**2).sum()\n", "print('alpha:')\n", "print(alpha)\n", "print('beta')\n", "print(beta)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can now solve for the parameters, $p$, given the set of equations\n", "$$ \\alpha_{jk} p_k = \\beta_{k} $$\n", "by inverting the $\\alpha$ matrix and multiplying by $\\beta$.\n", "Matrix inversion can be done using numpy.linalg.inv() and matrix multiplication by numpy.dot" ] }, { "cell_type": "code", "execution_count": 122, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "pars: [2.6167978 0.66856591]\n" ] } ], "source": [ "# now invert the matrix alpha, and multiply by beta to get the parameters\n", "c=np.linalg.inv(alpha)\n", "par=np.dot(c,beta)\n", "print('pars: ',par)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "How do the parameters compare with those derived specifically for the straight line fit? (They should be the same!)\n", "

ANSWER HERE :" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can calculate the uncertainties on the parameters from the error matrix (the inverse of the $\\alpha$ matrix: diagonal elements give the variances of the parameterws." ] }, { "cell_type": "code", "execution_count": 86, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "inverse matrix: \n", "[[0.1030303 0.00606061]\n", " [0.00606061 0.01212121]]\n", "uncertainties on individual parameters from sqrt of diagonal\n", "0.3209833376209784 0.11009637651263605\n" ] }, { "data": { "text/plain": [ "[]" ] }, "execution_count": 86, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWoAAAD4CAYAAADFAawfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAdgElEQVR4nO3deXTU9b3/8eeH7AmBAAk7SQgJ+05EEJdI0Crua9W4tGJTWxV6q7231ntPf20v3RcDKAraVtt0sdbueqsEAqKABEVRQCYJ2VnCloSErPP5/TFJBQoywEzmO5nX4xwP+X5n/M77zMl5nm++y4yx1iIiIs7VK9ADiIjIp1OoRUQcTqEWEXE4hVpExOEUahERhwv3x0YTExNtamqqPzYtItIjbdmy5YC1NulUj3kVamPMfwAPABbYBnzeWtt8uuenpqZSVFR0LrOKiIQkY0z56R4746EPY8wwYCGQaa2dCIQBd/huPBER+TTeHqMOB2KMMeFALFDjv5FEROR4Zwy1tbYa+DFQAewB6qy1r5/8PGNMrjGmyBhTVFtb6/tJRURClDeHPvoBNwAjgaFAnDHm7pOfZ61dYa3NtNZmJiWd8ni4iIicA28OfcwDdltra621bcArwEX+HUtERLp4E+oKYJYxJtYYY4BsYId/xxIRkS7eHKPeBLwMvIvn0rxewAo/zyUiIp28uo7aWvtN4Jt+nkVERE5Bt5CLiPhCVpbnPz9QqEVEHE6hFhFxOIVaRMThFGoREYdTqEVEHE6hFhFxOIVaRMThFGoREYdTqEVEHE6hFhFxOIVaRMThFGoREYdTqEVEHE6hFhFxOIVaRMThFGoREYdTqEVEHE6hFhFxOIVaROR85efDxo2wdi2kpnqWfUihFhE5H/n5kJsLLS2e5fJyz7IPY61Qi4icjyeegKamE9c1NXnW+4hCLSJyPioqzm79OVCoRUTOR3Ly2a0/Bwq1iMj5WLwYYmNPXBcb61nvIwq1iMj5yMmBFSsgKsqznJLiWc7J8dlLhPtsSyIioSonB1au9PxcWOjzzWuPWkTE4RRqERGHU6hFRBxOoRYRcTiFWkTE4RRqERGHU6hFRBxOoRYRcTivQm2MSTDGvGyM2WmM2WGMme3vwURExMPbOxPzgP+z1t5qjIkEYs/0P4iIiG+cMdTGmD7ApcDnAKy1rUCrf8cSEZEu3hz6SANqgV8YY94zxjxnjIk7+UnGmFxjTJExpqi2ttbng4qIhCpvQh0OTAeWW2unAY3A109+krV2hbU201qbmZSU5OMxRURClzehrgKqrLWbOpdfxhNuERE5TlOvCL9s94zHqK21e40xlcaYMdbaj4FsYLtfphERCUIfVtfx5OgbqYxO4DW3pVcv49Pte3vVxyNAfucVH6XA5306hYhIENpWVUdewS5W7dhPnz4jWLCniDa3m6heYT59Ha9Cba3dCmT69JVFRILU+5VHyCtwsXrnfvrGRPDoFaO577sP0aejFcJ9G2nQN7yIiHjtvYrD5BW4KPy4loTYCB67cjT3XZRKfHQEZL/ut9dVqEVEzmBLuSfQ63bV0i82gq99Zgz3XZRK76juSahCLSJyGkVlh8grcPGm6wD94yL5r6vGcs/slG4LdBeFWkTkJJvLDpG3ysX64gMMiIvk8avHcvesFOK6OdBdFGoRkU6bSg+SV+Di7ZKDJPaO5In548iZlUxsZGBTqVCLSHDLyvL8W1h4zpvYUHKQvIJdbCw9RGLvKP77mnHkXJhCTKTvr+A4Fwq1iIQkay0bSg/y5CoX7+w+RFJ8FP9z7XjumpnsmEB3UahFJKRYa3m75CB5q1y8U3aIgfFRfPO68dw5M5noCGcFuotCLSIhwVrL+uID5K1yUVR+mMF9ovnW9RP47AUjHBvoLgq1iPRo1lrWuQ6Qt2oX71YcYUjfaL5zwwRuy3R+oLso1CLSI1lrKdxVS94qF1srjzC0bzTfuXEit2cOJ8oPt3n7k0ItIj2KtZbCj2t5ssDF+5VHGJYQw+KbJnLrjOALdBeFWkR6BGstq3fuJ6/AxQdVdQzvF8P3bp7ELdOHExnu1fd4O5ZCLSJBzQKr+o1iybK32FZdx4j+MfzglkncPH04EWHBHeguCrVIMPLBTR7BzlrL69v3sWTSvXwUN4jkY2388NbJ3DRtWI8JdBeFWkSCittteX37XvIKitmxp56UsEh+VPIqNy5e2uMC3UWhFpGg4HZb/vnRXvIKXOzc28DIxDh+ctsUbnj0HsKx0EMjDQq1iDic22157cO9LClw8fG+BtIS4/jZZ6dw3eShhIf1wnOUumdTqEXEkTrclle37WHpahe79h1lVFIceXdM5drJQwnz8ZfHOp1CLSKO0uG2/P2DGpauLqZ4/1HSB/ZmyZ3TuGbSkJALdBeFWuRs6GoLv+kK9JICFyW1jYwe1Juld05jfggHuotCLSIB1d7h5m+de9CltY2MGRTPU3dN5+qJg+kV4oHuolCLSEC0d7j5y9Yalq0pZveBRsYOjmd5znQ+M0GBPplCLSLdqr3DzZ/eq+apNcWUHWxi3JA+PHP3DK4cP0iBPg2FWkS6RVtnoJetLqbiUBMThvbh2XtmcMW48wh0fj5s3AgtLZCaCosXQ06OT+d2AoVaRPyqrcPNK+9WsWxNMZWHjjFxWB9W3pvJvHEDMeY89qDz8yE31xNpgPJyzzL0uFgr1CLiF63tbv74bhVPrSmm6vAxJg/vy/+7bgJzx55noLs88QQ0NZ24rqnJs16hFhE5vdZ2N3/YUsnTa0qoPnKMKcP78p0bJpI1Jsk3ge5SUXF264OYQi0iPtHS3sEfiqp4ek0xNXXNTB2RwP/eNJGs0T4OdJfkZM/hjlOt72EUahE5Ly3tHby0uZKnC0vYU9fM9OQEvnfLZC7NSPRPoLssXuw5Jn384Y/YWM/6HkahFpFz0tzWwe83V7K8sIS99c1kpvTjh7dO5uJ0Pwe6S9dx6AULPCcUU1J01YeICHgC/dt3KnhmbQn76lu4ILUfP7l9CheNGtA9gT5eTg6sXOn5uQff1q9Qi4hXmts6yN9UwbNrS9jf0MLMkf352cYXmb2xAvNgYaDH69EUapFg0803eRxr7SB/UznPriultqGFWWn9ybtjGrNHDYDfPu6315VPeB1qY0wYUARUW2uv9d9IInJa3XiTR1NrO/kbK3h2XSkHjrYwO20AS++cxqy0AT59HTmzs9mjXgTsAPr4aRYROZNuuMmjqbWdX20oZ8W6Ug42tjInfQBPZ09n5sj+Ptm+nD2vQm2MGQ5cAywGvurXiUTk9Px4k0djSzsvbihn5ZulHGps5ZKMRBZlZ5CZqkAHmrd71E8C/wnEn+4JxphcIBcguQdecC7iCH64yeNoSzsvbihj5bpSDje1cenoJBZlZzAjpd+5zyk+dcZQG2OuBfZba7cYY7JO9zxr7QpgBUBmZmbP/7ZJkUDw4U0eDc1t/9qDPtLURtaYJBZmZzA9WYF2Gm/2qOcA1xtj5gPRQB9jzK+ttXf7dzQRh3HCR2r64CaP+uY2XnirjOfW76buWBuXj0li0bzRTB2R4Keh5XydMdTW2seBxwE696gfU6Ql5DjpIzXP8SaP+uY2frG+jOfXl1Lf3E722IEszM5gigLteLqOWsQbQfyRmnXH2vj5+t38/K3dNDS3M2/cIBZlZzBpeN9AjyZeOqtQW2sLgUK/TCLiZEH4kZp1TW08/9ZuftEZ6CvHD2JhdgYThynQwUZ71CLeCKKP1DzS1Mrz63fzy7fKaGhp56oJg3kkO50JQxXoYKVQi3gjCD5S83BjK8+tL+WFt8s52tLO1RMHszA7g3FDdI9asFOoRbzh4I/UPNTYynNvlvLC22U0tXUwf+IQHslOZ+xgBbqnUKhFvOWwj9Q8GB7Dytd28uKGMo61dXDNpCEszM5g9KDT3pcmQUqhFgkyB462sDL5Ml4cNJXmdSVcN3koj8xNJ0OB7rEUapEgUdvQwop1Jfx6YwUtQzK57uBOHvnug6QPDPFAO+CvG39TqEUcbn9DM8+uLSV/Uzmt7W5umDqMh3/xLUY1H4aBXwv0eNINFGoRh9pf38wznYFu63Bz47RhPHx5OmlJvWH54UCP54xb6kOEQi3iMPvqm1leWMJv36mg3W25qTPQqYlxgR7tE066pT4EKNQiDrG3rpnlhcX8dnMlHW7LzdOG8fDcdFIGOCjQXYL4lvpgpFCLBFjNkWMsLyzh95srcVvLLdOH89Dl6SQPiA30aKcXhLfUBzOFWiRAqo8cY3lhMS9trsJtLbdlDufLWemM6O/gQHcJolvqewKFWqSbVR1u4unCEv5QVAnAbZkj+HLWKIb3C4JAdwmCW+p7EoVapJtUHmri6cJiXt5SBcBnLxjBl7LSGZYQE+DJzoGDb6nviRRqET+rPNTEU2s8ge5lDHdckMyXskYx9HwC7YSbPBx2S31PplCL+EnFwSaWrXHxyrvV9OplyLkwmQezRjGkbxDuQUtAKdQiPlZ2oJFla4r503vVhPUy3D0rhQcvG8XgvtGBHk2ClEIt4iO7DzSydLWLv2ytIbyX4d7ZnkAP6qNAy/lRqEXOU2ntUZatLubPW6uJDO/F5y5K5YuXpjFQgRYfUahFzlHx/qMsW+3ir+/XEBneiwUXjyT30lEkxUcFejTpYRRqkbNUHN2fJb99j799UEN0eBhfuCSNL1yaRmJvBVr8Q6EW8dKufQ0syf0p/9i2h5gd+8i9NI3cS9IYoECLnynUImfw8d4Glqx28eq2PcRGhPHgZaP4wiVp9I+LDPRoEiIUapHT2Lm3niUFLl7dtpe4yDC+nDWKBy5Oo58CLd1MoRY5yfYaT6D/76O99I4K55G56Sy4eCQJsQq0BIZCLdLpw+o6lhS4eH37PuKjwlk4N537FWhxAIVagkNWludfP3ymxIfVdTy5ysWqHfuIjw5nUXYG988ZSd/YCJ+/lsi5UKglZG2rqiOvYBerduynT3Q4/zFvNJ+bk0rfGAVanEWhlpDzfuUR8gpcrN65n74xETx6xWjum5NKn2gFWpxJoZaQ8V7FYfIKXBR+XEtCbASPXTma+y5KJV6BFodTqKXH21LuCfS6XbX0i43ga58Zw30XpdI7Sr/+Ehz0myo9VlHZIfIKXLzpOkD/uEj+66qx3DM7RYGWoKPfWOlxNpcdIm+Vi/XFBxgQF8njV4/l7lkpxCnQvqdvdukW+s2VHmNT6UHyCly8XXKQxN6RfGO+J9Cxkfo1l+Cm32AJehtKDpJXsIuNpYdI7B3Ff18zjpwLU4iJDAv0aCI+ccZQG2NGAC8CgwE3sMJam+fvwUQ+jbWWDaUHyVvlYtPuQyTFR/E/147nrpnJCrT0ON7sUbcDj1pr3zXGxANbjDFvWGu3+3k2kX9jreXtEk+g3yk7xMD4KL553XjunJlMdIQCLT3TGUNtrd0D7On8ucEYswMYBijU0m0ssL5vCnnPbKCo/DCD+0Tzresn8NkLRijQ0uOd1TFqY0wqMA3YdIrHcoFcgOTkZB+MJuLZg17nOsCSCXexJX4YQ44c49s3TOD2TAVaQofXoTbG9Ab+CHzFWlt/8uPW2hXACoDMzEzrswklJFlrWburlidXudhaeYShkfF8Z/cb3J7/E6LCFWgJLV6F2hgTgSfS+dbaV/w7koQyay2FH9fyZIGL9yuPMCwhhsU3TeTWr99PlO0ARVpCkDdXfRjgeWCHtfan/h9JQpG1ltU795NX4OKDqjqGJcTwvZsnccv04USG9wLbEegRRQLGmz3qOcA9wDZjzNbOdd+w1r7qv7EkVFhrWbVjP0sKXGyrrmN4vxi+f/Mkbu4KtIh4ddXHesB0wywSQqy1vLF9H3kFLj6qqSe5fyw/vGUyN00fRkSYAi1yPN2ZKN3K7ba8vn0fSwpcbN9TT8qAWH5062RunKZAi5yOQi3dwu22/POjveQVuNi5t4HUAbH85LYp3DB1KOFnCnR+PmzcCC0tkJoKixdDTk63zC3iBAq1+JXbbXntw70sXe0JdFpiHD+9fQrXT/Ei0OCJdG6uJ9IA5eWeZVCsJWQYa31/yXNmZqYtKiry+XYleLjdllc/3MOSAhe79h0lLSmOhXMzuG7KUMJ6ncUpj9RUT5xPlpICZWW+Glck4IwxW6y1mad6THvU4lMdbss/tu1haYEL1/6jpA/sTd4dU7l28lkGuktFxdmtF+mBFGrxiQ635e8f1LCkwEVJbSMZA3uz9M5pzJ805NwC3SU5+dR71PqYAgkhCrWcl/YON3/7oIalq4sprW1k9KDePHXXdK6eOJhe5xPoLosXe45JNzV9si421rNeJEQo1HJO2jvc/PV9T6B3H2hk7OB4ns6ZzlUTfBToLl0nDBcs8JxQTEnRVR8SchRqOSvtHW7+vLWGZatdlB1sYtyQPjxz93SuHO/jQB8vJwdWrvT8rO/okxCkUItX2jrc/Om9ap5aU0z5wSbGD+nDs/fM4Ipxg/wXaBEBFGo5g7bL5/KnxPEsu+AWKg41MWFoH1bem8m8cQPxfF6XiPibQi2n1Nru5pV3q1g2ZQFV0QlMionguXszyVagRbqdQi0naG138/KWKp5aU0z1kWNMaT/Gt3cWcPn3nlegRQJEoRYAWto7+ENRFcsLS6g+coypIxL435smkvXF2z0fnahIiwSMQh3iWto7eGlzJcsLS6ipa2ZacgLfvXkSl2Ykag9axCEU6hDV3NbBS0WVPL2mhL31zcxI6cf3b5nMJQq0iOMo1CGmua2D371TwfK1Jeyrb+GC1H78+LYpzEkfoECLOJRCHSKa2zr4zaYKnllbwv6GFmaO7M/Pbp/K7FEKtIjTKdROlZXl+fc878RrbusgvzPQtQ0tXDiyP3l3TGP2qAHnPaKIdA+Fuoc61tpB/qZynllbyoGjLcxOG8DSO6cxK02BFgk2CnUP09Tazq83lrNiXSkHjrYyJ30AT2dPZ+bI/oEeTUTOkUJ9Kj467NCdGlva+dXGclauK+VgYysXpyeyaF4GF6Qq0CLBTqEOco0t7by4oZyVb5ZyqLGVSzIS+cq8DGakKNAiPYWzQh2Ee7KBcrSlnRfeLuO5N0s53NTGZaOTWJidwYyUfoEeTUR8zFmhljNqaG7zBHr9bo40tXH5GE+gpyUr0CI9lUIdJOqb2/jlW2U8v343dcfamDt2IAuzM5g6IiHQo4mInynUDld3rCvQpdQ3tzNvnCfQk4eHWKB1OExCmELtUHVhUfx8yAx+/oPVNDS3c8X4QSzKzmDisL6BHk1EuplC7TB1TW08v76UX0z7Ig3hUXxm1AAWZmcwYagCLRKqFGqHONLUyvPrd/PLt8poaGnn6royHqnewPjv/yFwQ+Xnw8aNnm//Tk3Vt3+LBIhCHWCHG1t5bn0pL7xdztGWduZPGswjczMYd+ePAjtYfj7k5noiDVBe7lkGxVqkmynUAXKosZWVb5by4ttlNLV1MH/SEBbOzWDM4PhAj+bxxBPQ1HTiuqYmz3qFWqRbKdTd7ODRFla8WcqvNpRzrK2DaycP5ZG56Ywe5JBAd6moOLv1IuI3CnU3OXC0hZXrSnlxQznN7R1cP8UT6PSBDgt0l+Rkz+GOU60XkW6lUPtZbUMLK9aV8OuNFbR0BvrhuRmkD+wd6NE+3eLFnmPSxx/+iI31rBeRbuVVqI0xVwF5QBjwnLX2+36dqgfY39DMs2tLyd9UTmu7mxunDuPhuemkJTk80F26jkMvWOA5oZiSoqs+RALkjKE2xoQBTwFXAFXAZmPMX6212/09XDDaX9/M8rUl/GZTBe1u+69Aj0yMC/RoZy8nB1au9PysOwNFAsabPeqZQLG1thTAGPM74AZAoT7O3rpmnllbwm/eqaDDbbl52jAeujyd1HMJtK5fFpHjeBPqYUDlcctVwIUnP8kYkwvkAiSH0AmnPXXHWF5Ywu82V+J2W26ZPpyHLk8neUDsuW1Q1y+LyEm8CfWpvqLa/tsKa1cAKwAyMzP/7fGg4eXebM0RT6B/v7kSt7XcOsMT6BH9zzHQXXT9soicxJtQVwEjjlseDtT4Z5wA82JvtvrIMZ5eU8xLRZ4/Mm6dMYIvZ406/0B30fXLInISb0K9GcgwxowEqoE7gLv8OlWgfMrebNX8m3hqTQkvb/EE+vbMEXwpaxTD+/ko0F10/bKInOSMobbWthtjHgb+iefyvJ9baz/y+SROOIF2mr1WW1FB1o8K6WUMd1yQzJeyRjE0IcY/M+j6ZRE5iVfXUVtrXwVe9dsUTjmBdpq92Zr4JO660BPoIX39FOguun5ZRE7SK9ADAJ9+Aq07LV6MO+bEELdGRRP3kx/w7Rsm+j/SXXJyYNYsuOwyKCtTpEVCnDNC7YATaLsPNPJoxES+esVDVPdJwgIdI0YQ+fxzJDzwuW6bQ0TkZM74rI8AnkArrT3KstXF/HlrNZHhvcj5wueJaNqOaWskTHfjiYgDOCPUATiBVrz/KMtWu/jr+zVEhvfi/jkjyb0sjYHx0fDjRr+9rojI2XJGqLvxBFrx/gaWri7mr+/XEB0exgOXpPGFS9JIio/y+WuJiPiCM0INfv8AINe+BpasLubvH9QQExFG7qWeQCf2VqBFxNmcE2o/+XhvA0tWu3h12x5iI8J48LJRPHDxSAYo0CISJHpsqHfurWdpQTH/2LaHuMgwvnTZKB64JI3+cZGBHk1E5Kz0uFDv2FPPkgIXr324l95R4Tx8eToLLh5JPwVaRIJUjwn1RzV1LClw8c+P9hEfFc7Cuencf/FIEmIVaBEJbkEf6g+r68grcPHG9n3ER4ezKDuD++eMpG9sRKBH6xl0LblIwAVtqLdVeQK9asc++kSH85V5GXx+zkj6xijQItKzBF2oP6g6Qt4qFwU799M3JoKvXjGaz81JpU+0Ai0iPVPQhHpr5RHyVu1izce1JMRG8NiVo7nvolTiFWgR6eEcH+r3Kg6TV+CisDPQX/vMGO6dnaJAi0jIcGyot5R7Ar1uVy39YiP4z6vGcO/sVHpHOXZkERG/cFz1inoPI+/5TbzpOkD/uEi+fvVY7pmVQpwCLSIhyjH1q29u48tjb2N9QiqJe+r5xvyx3D0rhdhIx4zYvXRZnIh0ckwF46PCiXG38UT5GnJe+H7oBlpE5CSOqaExhpW7/uxZUKRFRP5FRTwVHXYQEQdxxncmiojIaSnUIiIO56xDHzrkICLyb7RHLSLicAq1iIjDKdQiIg6nUIuIOJxCLSLicAq1iIjDKdQiIg6nUIuIOJxCLSLicMZa6/uNGlMLlPt8w90rETgQ6CEcQu/FifR+nEjvxyfO571IsdYmneoBv4S6JzDGFFlrMwM9hxPovTiR3o8T6f34hL/eCx36EBFxOIVaRMThFOrTWxHoARxE78WJ9H6cSO/HJ/zyXugYtYiIw2mPWkTE4RRqERGHU6jPwBjzmDHGGmMSAz1LIBljfmSM2WmM+cAY8ydjTEKgZ+puxpirjDEfG2OKjTFfD/Q8gWSMGWGMWWOM2WGM+cgYsyjQMwWaMSbMGPOeMebvvt62Qv0pjDEjgCuAikDP4gBvABOttZOBXcDjAZ6nWxljwoCngKuB8cCdxpjxgZ0qoNqBR62144BZwEMh/n4ALAJ2+GPDCvWn+xnwn0DIn3G11r5urW3vXNwIDA/kPAEwEyi21pZaa1uB3wE3BHimgLHW7rHWvtv5cwOeQA0L7FSBY4wZDlwDPOeP7SvUp2GMuR6otta+H+hZHOh+4LVAD9HNhgGVxy1XEcJhOp4xJhWYBmwK7CQB9SSenTq3PzburG8h72bGmFXA4FM89ATwDeDK7p0osD7t/bDW/qXzOU/g+bM3vztncwBzinUh/5eWMaY38EfgK9ba+kDPEwjGmGuB/dbaLcaYLH+8RkiH2lo771TrjTGTgJHA+8YY8PyZ/64xZqa1dm83jtitTvd+dDHG3AdcC2Tb0LsAvwoYcdzycKAmQLM4gjEmAk+k8621rwR6ngCaA1xvjJkPRAN9jDG/ttbe7asX0A0vXjDGlAGZ1tqQ/YQwY8xVwE+By6y1tYGep7sZY8LxnETNBqqBzcBd1tqPAjpYgBjPHswLwCFr7VcCPY9TdO5RP2atvdaX29UxavHWMiAeeMMYs9UY80ygB+pOnSdSHwb+iefE2UuhGulOc4B7gLmdvw9bO/coxQ+0Ry0i4nDaoxYRcTiFWkTE4RRqERGHU6hFRBxOoRYRcTiFWkTE4RRqERGH+/+cBiP3deOAbwAAAABJRU5ErkJggg==\n", "text/plain": [ "

" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# uncertainties from inverse matrix\n", "print('inverse matrix: ')\n", "print(c)\n", "print('uncertainties on individual parameters from sqrt of diagonal')\n", "print(np.sqrt(c[0,0]),np.sqrt(c[1,1]))\n", "plt.errorbar(x,y,sig,fmt='ro')\n", "plt.plot(x,func(x,par))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can simplify the formulation even further, using the design matrix formalism and matrix arithmetic. Here, the design matrix is the matrix of derivatives for each point, i.e. an N (number of points) by M (number of parameters) matrix, divided by the uncertainty on each point You can use your deriv() routine to supply the design matrix:\n", "\n", "$$A_{ij} = {f_j(x_i) \\over \\sigma_i}$$\n", "with $N$ rows and $M$ columns. \n", "The data vector is the observed data points divided by the uncertainty.\n", "$$b_i = {y_i \\over \\sigma_i}$$\n", "\n" ] }, { "cell_type": "code", "execution_count": 89, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(10, 2)\n", "design\n", "[[ 1. -5.]\n", " [ 1. -4.]\n", " [ 1. -3.]\n", " [ 1. -2.]\n", " [ 1. -1.]\n", " [ 1. 0.]\n", " [ 1. 1.]\n", " [ 1. 2.]\n", " [ 1. 3.]\n", " [ 1. 4.]]\n", "b:\n", "[-0.34247188 -0.2221414 -0.08567665 0.26244425 1.54440844 3.78254897\n", " 4.08243545 3.27141421 4.49359491 6.72487196]\n" ] } ], "source": [ "# if your deriv routine returns a [npar,npoints] array, transpose the matrix (.T) to get [npoints,npar]\n", "design=(deriv(x,2)/sig).T\n", "print(design.shape)\n", "\n", "# design array is the derivatives for all points, divided by the uncertainty\n", "#design=np.array([deriv(x,2)[0]/sig,deriv(x,2)[1]/sig]).T\n", "print('design')\n", "print(design)\n", "\n", "# b array is data points / sig\n", "b=np.array(y/sig)\n", "print('b:')\n", "print(b)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we set up the matrix arithmetic. What we previously called the alpha array is just $A^T A$ (where the multiplication is a matrix dot product), and the beta array is $A^T b$, where $A$ is the design matrix, and $b$ is the vector from above. Confirm that you can reproduce the $\\alpha$ matrix and $\\beta$ array with the design matrix formalism." ] }, { "cell_type": "code", "execution_count": 90, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[10. -5.]\n", " [-5. 85.]]\n", "[[10. -5.]\n", " [-5. 85.]]\n", "beta\n", "[23.51142827 51.79419443]\n", "[23.51142827 51.79419443]\n" ] } ], "source": [ "ATA=np.dot(design.T,design)\n", "print(ATA)\n", "print(alpha)\n", "print('beta')\n", "print(np.dot(design.T,b))\n", "print(beta)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can solve this set of (2) equations using [np.linalg.solve()](https://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.solve.html), which is a solver that does not invert the matrix. Confirm that you get the same answers as above." ] }, { "cell_type": "code", "execution_count": 91, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[2.73629379 0.77030192]\n" ] } ], "source": [ "w = np.linalg.solve(ATA,np.dot(design.T,b))\n", "print(w)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here we invert the matrix using [numpy.linalg.inv()](https://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.inv.html) to get the uncertainties. If we're going to do this anyway, you would just use the inverted matrix to get the solution by multiplying by the beta array. Again, same answers!" ] }, { "cell_type": "code", "execution_count": 92, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[0.1030303 0.00606061]\n", " [0.00606061 0.01212121]]\n", "0.3209833376209784 0.11009637651263605\n", "[2.73629379 0.77030192]\n" ] } ], "source": [ "c = np.linalg.inv(ATA)\n", "print(c)\n", "print(np.sqrt(c[0,0]),np.sqrt(c[1,1]))\n", "w = np.dot(c,np.dot(design.T,b))\n", "print(w)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here we get the solution using [scipy.linalg.lstsq()](https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.lstsq.html), just providing the design matrix and the data array. Very simple!" ] }, { "cell_type": "code", "execution_count": 93, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([2.73629379, 0.77030192])" ] }, "execution_count": 93, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import scipy.linalg\n", "scipy.linalg.lstsq(design,b)[0]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "OK, let's play with this for general linear functions. Write routines to simulate the function\n", "$$ y= a + bx + c sin(x)$$\n", "and the routines to solve for the least squares parameters" ] }, { "cell_type": "code", "execution_count": 74, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "True parameters: [13, 13, 13]\n", "Derived parameters: [14.01412538 11.02617948 12.73189632]\n", "Derived parameters: [14.01412538 11.02617948 12.73189632]\n" ] }, { "data": { "text/plain": [ "[]" ] }, "execution_count": 74, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXxU1f3/8deZ7AmBEBKykJV9XwMEEAUVRVQUv7hgULRYrLWtXbWt/bbW35cutmpb29piXVDjCkXAXREE2cMSFtlDQvaEQCAh62TO748ZFTEhk2Rm7tyZz/PxmEdmbm5y35ckH86ce+45SmuNEEII87EYHUAIIUTnSAEXQgiTkgIuhBAmJQVcCCFMSgq4EEKYVKAnDxYTE6PT0tI8eUghhDC9HTt2nNRax1643aMFPC0tjZycHE8eUgghTE8pVdDadulCEUIIk5ICLoQQJiUFXAghTEoKuBBCmJQUcCGEMCkp4EIIYVJSwIUQwqSkgAshhElJARdCGGPaNPtDdJoUcCGEMCmP3kovhPBP1hYbB0pr2Hq8ivKzDQQHWgjqM5n4phoura4nMSrM6Iim1G4BV0qFAuuBEMf+y7TWv1FKvQBcBpxx7HqX1nq3u4IKIcznaEUtT31yhDUHKqhttAIQGmTB2qKxJk+x7/SHTxgcH8msEQl865J0uoVIu9JZzvxLNQKXa61rlVJBwGdKqfccn/uZ1nqZ++IJIcyo8FQdT358mLd2FRMaFMANoxPJ7NuLiem9iO8RCtnZ6IULobGR2rhEnp/1bZ4oG8+Lm/P54ZUDuW18MoEB0sPbnnYLuLavelzreBnkeMhKyEKIVq3OLeGh5XtosWkWXpLOdy7rR69uIV/tkJ0NixahGhsBiCwv4Qev/4kb//AXfhY8nF+9tY8XN+fzz6xx9O/dzZiTMAnlzKr0SqkAYAfQH/iH1vohRxfKJOwt9DXAz7XWja187SJgEUBKSsq4goJWZ0UUQphcc4uNP7x3kGc/O8641J48NW9M633baWnQWh1ITUUfP85Hn5fzi//upclq4y+3jeaKIXFuz+7tlFI7tNYZ39juTAE/75tEASuA7wNVQBkQDCwBjmmtH73Y12dkZGiZD1wI31PXZGXhCzlszqvirslp/HLWEIID2+gCsVigtbqjFNhsABRX13PvSznsLznLT68axHen9UMp5cYz8G5tFfAOdTJprauBdcBMrXWptmsEngcmuCSpEMJUGppbuGdpDluPV/H4zaN4ZPawtos3QEpKu9v7RIXx5r2TmT0qkT99cIjHPzzs4tS+od0CrpSKdbS8UUqFAVcCB5VSCY5tCrgR2OfOoEII79NktXHfyzvYnFfF47eM4n/GJbX/RYsXQ3j417eFh9u3nycsOIC/3Dqa28Yn8/e1R/nnuqMuTO4bnBmFkgAsdfSDW4A3tNZvK6U+UUrFAgrYDXzHjTmFEF7GZtP88PVdrD1Uye/mjGDOGCeKN0BWlv2jYxQKqan24v3F9vMopVg8ZwR1TS089v4hIoIDWTA5zXUnYXLOjELZA4xpZfvlbkkkhDCFv689yrt7y3h41hBun9hGt0hbsrLgmWfsz9etu+iuARbF47eMor65hd+s2k9iVBgzhsqFTZBb6YUQnbDhSCVPfnyYG0cncs/UdLcfLyjAwlPzxjCiTw9+/MZu8k+ec/sxzUAKuBCiQ0qq6/nBq7sY0Lsbv7tpROdHh6xb127r+3yhQQE8PX8sARbFd17eQX1TS+eO60OkgAshnNbcYuP+V3bS3KJ5ev44woM9e9t7Us9w/nrbGA6V1/DLFXvpyDBoXyQFXAjhtCXr89h1opo//M8I+sUac5fkZQNj+fGVA1mxq5jlO4sNyeAtpIALIZxytKKWv645wqwR8Vw3MtHQLPdP78+EtGh+u3o/ZWcaDM1iJCngQoh22Wyah5bvISwogN/OHm50HCwWxWNzR9LcYuMX/93jt10pUsCFEO16aUsBOwpO8+vrhhIbGdL+F3hAWkwED80czNpDlSzbUWR0HENIARdCXFTpmXr++P5BLhsYy01j+xgd52sWTEpjQno0j779uV92pUgBF0Jc1GPvH8Jq0/zfjcO9bkIpi0Xxp7kjabLaWPzuAaPjeJwUcCFEm3YXVrNiVzH3XJJOcnR4+19ggNReEdx7WT9W55awJa/K6DgeJQVcCNEqrTX/9/bnxHQL4bvT+xsd56Luu6wffaLCeGTVfqwtNqPjeIwUcCFEq97dW0ZOwWl+etVAr1+nMiw4gF9dO4SDZTW8su2E0XE8Rgq4EOIbGppb+P17BxiS0J2bM5KNjuOUmcPjmdK/F3/+4BBVtd9YHMwnSQEXQnzDq9tOUHS6nl9dO4QAi3dduGyLUopHrh/GuaYWnvrEP+YOlwIuhPiahuYW/rnuGJl9o5nSP8boOB0yIC6SWzKSyN5aQOGpOqPjuJ0UcCHE17y8pYDKmkZ+dOVAo6N0ygNXDMSiFE985PvLsEkBF0J8qb6phX99msfkfr2Y2LeX0XE6Jb5HKHdPSeet3cUcKD1rdBy3kgIuhPjSy1sKOFnbyI9mmLP1/YX7LutHZEggf/rgUOs7TJtmf5icFHAhBAB1TVb+9ekxLukfw/i0aKPjdEmP8CDum9afTw5WsO34KaPjuI0UcCEEAK9uK6TqXBM/vHKA0VFc4q7JacRGhvCkD/eFSwEXQtDcYuPZDXlMSIsmw+St7y+EBQdw76V92ZxXxfZ832yFSwEXQvDOnlJKzjRw72V9jY7iUlkTU4npFszf1hz5amN2NmzZAp9+Cmlp9tcmJQVcCD+nteZfnx6jf+9uTB/U2+g4LhUWHMC3p/Zlw5GT7Dxx2l6sFy2CRsedmgUF9tcmLeJSwIXwZ9OmseHGuzlYVsOiS/tiMcldlx0xPzOVnuFBPLXmCDz8MNRdcINPXZ19uwm1W8CVUqFKqW1KqVyl1H6l1G8d29OVUluVUkeUUq8rpYLdH1cI4Wr/TpxA78gQbhht7DqX7hIREsg9U/uy9lAl+kQbE121td3LOdMCbwQu11qPAkYDM5VSmcAfgSe11gOA08BC98UUQrjDvvDebOyRyt1T0gkJDDA6jtvcOSmV7qGBnIqOa32HlBTPBnKRdgu4tqt1vAxyPDRwObDMsX0pcKNbEgoh3CM7m8T1H5H3x+u5987LTdsP7IzI0CDunJTGo5OzsIWFff2T4eGweLExwbrIqT5wpVSAUmo3UAF8BBwDqrXWVscuRUCri+UppRYppXKUUjmVlZWuyCyE6KrsbPS3FxF9rhoLGkvhCVNfzHPGgslpvDfyCt5Y9L8Q4liYOTUVliyBrCxjw3WSUwVca92itR4NJAETgCGt7dbG1y7RWmdorTNiY2M7n1QI4ToPP4yq952Lec6IjQxh7rgkfh0xmoqpV8Bll0F+vmmLN3RwFIrWuhpYB2QCUUqpL5bpSAJKXBtNCOEuvnYxz1mLpval2WbjhfixRkdxCWdGocQqpaIcz8OAK4EDwFpgrmO3BcBKd4UUQrhWQ3wbI05MejHPWWkxEVwzPJ6X+l1C7QcfGx2ny5xpgScAa5VSe4DtwEda67eBh4AfK6WOAr2AZ90XUwjhSv+55h7qg0K+vtHEF/M64t5L+1HTYOXVreZ/t+HMKJQ9WusxWuuRWuvhWutHHdvztNYTtNb9tdY3a639YxE6IUzuYNlZHo8dz8YHf+8zF/M6YlRyFBPTo3lhU77pV7CXOzGF8DPZW04QHGhh3C/uh8xMn7iY11HfuiSd4up6Pvy83OgoXSIFXAg/Utdk5a1dxVw3IoGeEf578/SVQ+JIiQ7nuc+OGx2lSwLb30W4wsnaRvYWn6Gx2YbVZiNAKUYmR9EnKqztL/pixZB16zwRUfiBt/eUUtNoZd5E375Y2Z4Ai+KuyWk8+vbn5BZWMyo5yuhInSIF3I1Kz9Tz+vZC1h6qZE9RNbqVkfLJ0WFM7hvDgslpDE3s7vmQwq+8uu0E/Xt3IyO1p32DHzcObs5I4omPDvPcxuP89bYxRsfpFCngblDbaOXfnx7jmQ15NFptjE6O4kdXDiSzby8iQgIIDrBQ39xCTv5ptuRV8c7eUl7PKWTWiHh+eOVABsZFGn0KwgcdLDvLrhPV/O91Q1HK92Yd7KjI0CBuyUjmxc35/OKaIcT3CDU6UodJAXex9/eV8au39nGytpHZoxL52dWDSI4Ob3XfkUlRfOuSdM7UNfPsZ3k8tzGf9/aVcf+0/vxoxkB8d2ohYYTXthUSHGDhpjGtznrhl+6anMbzm47z4uZ8Hpw52Og4HSYXMV3EZtP85ePDfOflHST0CGXFdyfzt3lj2ize5+sRHsSPrxrEhgenc/O4JP6+9ihPL3qUlq1bfWLVEGG8+qYW/ruziGtGxPv1xcsLpfQKZ8aQOF7bXkhDc4vRcTpMCrgL1DVZ+d6rO/nLx0e4aWwf3vzOJMak9Ozw9+kZEcxjc0fxasQxFi79PQENDfZPmHzVEGG8d/eWcrbByrwJ/n3xsjULJqdx6lwTb+8pNTpKh0kB76KG5hbuen477+0r45ezBvP4zaMIDepa58ekZ58gzHrBfVE+PtGQcK/XcwpJj4lgYrpvLFjsSpP79aJ/724s3ZSPbm2kgReTAt4F1hYb33tlJ9vzT/GXW0ez6NJ+rrk45KcTDQn3KKg6x7bjp5g7LkkuXrZCKcWCSansLT7D7sJqo+N0iBTwTrLZNA8t38vHByp49Ibh3DDahReG2phQqLlPkuuOIfzG8h1FWBTcNFYuXrZlztgkuoUE8uLmAqOjdIgU8E7604eHWL6ziB/PGMgdmamu/eaLF9snFjpPfVAIi6fM5/S5JtceS/g0m02zfGcxlwyIJaHHRW4a83PdQgKZOy6Jt/eUUFljnmmdpIB3wpoD5Ty97hjzJqTw/cv7u/4AWVn2iYXOm2io/M9/45V+U/neqztNPwGP8JxNx6oorq7n5nHy7q09d0xKpblF89o283RVSgHvoOLqen7yZi5DE7rzm+vdeENEVtbXJhpK+8Ei/m/OcDYereKP7x90zzGFz3lzRyHdQwOZMbSNxXzFl/rFdmPqgBiyt54wTSNJCngHNLfY+P4rO7G2aP6RNbbLo0066paMZBZMSuWZDcdZlSsLIImLO9vQzPv7ypg9OtHjv6tmdUdmKmVnG1hzsMLoKE6RAt4BT3x0mJ0nqvn9TSNIj4lw/wHXrfvGXBW/um4o41J78vCKvZSeqXd/BmFab+eW0mi1cfO4ZKOjmMblg3uT0COUl7eY42KmFHAn7S06w78/PcatGclcP6qN5ag8ICjAwhO3jMLaonlw2R7TjVsVnrN8ZxEDendjZFIPo6OYRmCAhdsnpLDhyEmOnzxndJx2SQF3grXFxs//u4eYbiH88tohRschtVcEv7x2CBuOnDRNS0F4VkHVOXYUnGbO2D4y9ruDbp2QTKBFkW2Cvy0p4E54buNx9pec5bezh9EjLMjoOADMn5jC1AEx/O7dg+SboKUgPGvFrmKUghtdeX+Cn+gdGcrVw+N5c0cR9U3ePT+KFPB2nKiq44mPDnPlkDhmDo83Os6XlFI8NnckgQGKh5ZLV4r4itaaFbuKyUzvReLFFgwRbZo/MZUz9c2s3uPdgwWkgLfjf1fuI0ApHr1hmNe9FU3oEcZDMwez9fgpGZUivrSrsJqCqjrmyJ2XnZbZN5r+vbt5fTeKFPCL+PRwJZ8eruRHMwZ6bUtm3oQURvTpweJ3DlDT0Gx0HOEFVuwsJiTQwjVe9I7RbJRSzJ+YQm7RGfYWnTE6TpukgLehxab53TsHSIkO545JLr5V3oUCLIr/d+NwKmsb+evHR4yOIwzWZLWxek8JVw2LJzLUO67XmNWcsUmEBll4ZZv3tsLbLeBKqWSl1Fql1AGl1H6l1AOO7Y8opYqVUrsdj1nuj+s5y3YUcqi8hodmDiYk0LtvghidHMVt45N5flM+h8pqjI4jDLTuUAXVdc2y6o4L9AgLYvaoRFbuLuGsl767daYFbgV+orUeAmQC9yulhjo+96TWerTj8a7bUnpYXZOVxz88zJiUKGaNMMfb0J9dPZjI0EB+u3q/XND0Yyt2FdMrIphLBsQYHcUnZE1Mpa6phZW7io2O0qp2C7jWulRrvdPxvAY4APj0f+/PrD9ORU0jv7p2CGr6dJg2zehI7YqOCOaBKwaw6VgV64+cNDqOMMDZhmbWHKzg+lGJBAVI76grjEqOYkSfHmRvPeGVDaMO/ZSVUmnAGGCrY9P3lFJ7lFLPKaVaXUNMKbVIKZWjlMqprKzsUlhPOH2uiSXrj3HN8HjGpZpr9ZKsiakkR4fxh/cOYrN53y+bcK8P9pXRZLUxe7Rxdwr7oqyJKRwsq2HnidNGR/kGpwu4UqobsBz4odb6LPA00A8YDZQCj7f2dVrrJVrrDK11RmxsrAsiu9d/PsujrrmFH88YaHSUDgsOtPDTqwZxoPQsK3O98y2fcJ9VuSWkRIczJjnK6Cg+5fpRiUSGBJK9xfummXWqgCulgrAX72yt9X8BtNblWusWrbUNeAaY4L6YnlFd18TSTQXMGpHAgLhI+yLCW7aYamX460cmMqJPD/78wWFTrrItOqeipoGNR09yw+hEr7tfwewiQgKZM7YPb+8t9boFVZwZhaKAZ4EDWusnztuecN5uc4B9ro/nWc9+dpzaRis/uHyAvVgvWgSNjtU5TLIyvMWi+Pk1gymurpd5UvzI27ml2DTcIN0nbnH7xBSarDaW7ywyOsrXONMCnwLcAVx+wZDBx5RSe5VSe4DpwI/cGdTdztQ188LGfGaNiGdQfKR9Bfi6uq/vZJKV4af0j2HqgBj+ue4Y5xqtRscRHrAyt4ShCd3p3zvS6Cg+aXB8d8al9vS6i5nOjEL5TGuttNYjzx8yqLW+Q2s9wrF9tta61BOB3eXZjcepabTygysG2DeYfGX4H80YyKlzTbwkrXCfl3/yHLmF1dL6drOsiSkcP3mOzceqjI7yJRlrBNQ0NPP8xuNcPSyOwfHd7RvbWBm+ze1eZmxKTy4dGMuS9XnUNUkr3Jet3F2CUsjoEzebNSKBqPAgsrd6TyNOCjjw+vZCahqs3D/9vAWKW1kZnvBw+3aTeOCKAfZW+GZphfsqrTWrcosZnxYtq867WWhQAHPHJvHB/jIqahqMjgNIAae5xcbzG/OZmB7NyKTzhl+1sjI8S5bYt5vEuNSeTB0QI61wH3agtIZjleeYbeAqUf5k3sQUrDbNmznecTHT7wv4u3tLKa6uZ9Glfb/5yQtWhjdT8f7CD68cQNW5JhmR4qNW5ZYQYFHMGpHQ/s6iy/rFdmNS3168svUELV5ws5xfF3CtNc9syKNfbATTB/U2Oo5bjEuN/rIVLuPCfYvWmtW5JVzSP4boiGCj4/iNrMwUiqvrWX/Y+DvL/bqAb86rYl/xWb49tS8Wi+/e/PDdaf05WdvEsh3e8bZPuMauwmqKq+ul+8TDrhoaT0y3EK94V+vXBfyZ9XnEdAvmxotNvblunf1hYpl9oxmVHMWS9XlYW2xGxxEusmp3CcGBFmYMizM6il8JDrRw2/hkPjlUQdHpuva/wI38toAfrahl7aFK7pyURmiQd8/33VVKKe67rC8nTtXx3r4yo+MIF2ixad7ZW8r0QbF0l4UbPG7exBQU8Nq2QkNz+G0Bf3lLAcEBFm6faI5x3V01Y2g8fWMi+Nenx7zqTjLROVuPV1FZ08jsUT49s7PX6hMVxvRBvXlteyHNBr6r9csCXttoZfmOImaNsPdl+YMAi+Ley/qyv+Qsnx2V+cLNbnVuCRHBAVw+2DcvvpvB/MxUTtY28uH+csMy+GUBX7GrmJpGK3dMSjM6ikfdOKYPcd1DeHrdMaOjiC5obrHx3r4yrhwaR1iwb3f/ebNLB8aS1DPM0IuZflfAtda8tDmf4X26MzbFv+ZNDgkM4O4p6Ww6VsX+koustD1tmilWIfJXnx09SXVdM9ePlNEnRgqwKG6fmMLmvCqOVtQaksHvCvjW46c4XF7LnZlpfjlv8rzxKYQHB/DsZ8eNjiI6aXVuCZGhgUwdKOteGu3WjGSCAyyGtcL9roC/uDmfHmFBXO+nY2d7hAdx87gkVueWUHHWO+ZzEM5raG7ho/3lzBwWT0igdJ8YrVe3EK4dmcDyHUWGTN3sVwW8/GwDH+wv59bxyX7dd3j3lHSsNt36VLMmXIXIn6w/XElNo5Xr/LQB4o3mZ6ZS02hlhQEr1/tVAX8zp5AWm+b2Cf4xdLAtaTERXDE4juytJ75+e71JVyHyJ6v3lBIdEczkfr2MjiIcxqZEMSyxOy9vKfD4EF2/KeA2m+b1nEIm9e1FWkyE0XEMt/CSdE6da/p6q8HEqxD5g/qmFtYcKGfm8HiCAvzmT9frKaW4c1IqB8tq2J7v2ZXr/ea3YOOxkxSeque2CclGR/EKmX2jGZrQnec+O/5Vq8HkqxD5uk8OVlDX1MJ1I2XmQW8ze1QfuocG8uLmfI8e128K+GvbCokKD+LqYfFGR/EKSikWXpLOkYpaNh51LBFl8lWIfN3q3BJiI0OYmC7dJ94mLDiAmzOSeX9fmUcHB/hFAa+qbeTDz8u4aUySz8970hHXjkygV0QwSzfn2zf4wCpEvqq20craQxVcOyKBAB+eOdPM7shMpUVrXvbgkmt+UcCX7yyiuUUzT7pPviY0KIDbJiSz5kA5hafqfGIVIl/18eflNFpt0n3ixdJi7OsKvLK1gEarZ+be9/kCrrXmte2FjEvtyYC4SKPjeJ35makopb66EcEHViHyRatzS0jsEcrYlJ5GRxEXcdfkNE7WNvHOnlKPHM/nC3hOwWnyKs9x23hpfbcmoUcYVw+L47XthdQ3yYo93uhMXTPrj1Ry7cgEn154xBdMHRBDv9gInt+Y75EhhT5fwJflFBERHMC18tazTQsmpXGmvpmVuz1/I4Jo3wefl9HcorlO5j7xekop7pqcxt7iM+w8Ue3247VbwJVSyUqptUqpA0qp/UqpBxzbo5VSHymljjg+et17u7omK+/sLWXWiATCgwONjuO1JqRHMzg+khc2OVoNPrAKkS9ZnVtCSnQ4I5N6GB1FOOGmsUlEhgbywqZ8tx/LmRa4FfiJ1noIkAncr5QaCvwcWKO1HgCscbz2Kh/sL6O20crccUlGR/FqX7QaDpbVsKPAszciiIurqm1k07EqrhuZ4JeTr5lRREggt2Yk897eUkrP1Lv1WO0WcK11qdZ6p+N5DXAA6APcACx17LYUuNFdITtr2Y4iUqLDGZ8WbXQUrzd7dCKRoYGtz48iDPP+/jJabNJ9YjYLJqdh05qlm9z799ShPnClVBowBtgKxGmtS8Fe5IFWlwZRSi1SSuUopXIqKyu7lrYDik7XselYFf8zNkku/DghPDiQueOSeHdvKSdrG42OIxzezi2lb2wEQxJkBJWZJEeHM3N4PK9sLXDrLIVOF3ClVDdgOfBDrfVZZ79Oa71Ea52htc6IjY3tTMZOWbGzGK3hprGyZqCz5mem0tyieX27sQu1CruKsw1sOV7FdSMTpfvEhO6Z2pezDVbeyHHf35NTBVwpFYS9eGdrrf/r2FyulEpwfD4BqHBPxI7TWrNsZxGT+vYiOTq8/S8QAPSL7caU/r14ZesJWmyy8LHR3tlbitYwe5SMoDKjsSk9GZfak+c2Hnfb35Mzo1AU8CxwQGv9xHmfWgUscDxfAKx0fbzOySk4TUFVnVy87IQ7MlMprq5n7UGv+f/Yb63OLWFIQnf695buE7P69tR0Ck/V88H+Mrd8f2da4FOAO4DLlVK7HY9ZwB+AGUqpI8AMx2uv8N+dxYQFBTBzuExc1VFXDokjrnuIXMw0WOGpOnaeqOZ6aX2b2oyh8aT2CueZDXlu+f7OjEL5TGuttNYjtdajHY93tdZVWusrtNYDHB9PuSVhBzVaW3hnTwkzh8cTESJjvzsqMMDC7RNS+fRwJfknzxkdx2+9s9d+K7YsXGxuARbFt6aks+tENTsKXF8ife5OzLUHKzjbYOXGMXLxsrNum5BMgEXx6jaZB9woq3NLGJ0cJddwfMDNGUlkTUyhV0SIy7+3zxXwFbuKiY0MYYosOdVpcd1DmTEkjjdyCj02q5r4yrHKWvaXnJWZB31EeHAgi+eMcMtKYD5VwKvrmvjkYAWzRyUSKEtOdUlWZgqn65p5f597Lr6Itr2dW4pSyM07ol0+VeXe2VtKc4tmjnSfdNmUfjGk9gone4t0o3iS1prVe0oYnxZNfI9Qo+MIL+dTBfytXcUM6N2NYYndjY5iehaLImtiCtvyT3GorMboOH7jYFkNRytquX6UtL5F+3ymgBeeqmN7/mluHNNH7lpzkbnjkgkOsPDKVhlS6Ckrd5cQYFHMkiGwwgk+U8C/mMv6htHScnGV6IhgZo2I5787i6lrct98DsLOZtOszi1h6oAYenVz/YgF4Xt8ooBrrVm5u4QJadEk9ZRhV66UlZlKTaOV1bklRkfxeTtPnKa4up7Z0n0inOQTBfxgWQ1HKmqZLa1vl8tI7cnAuG5ke3ClbX+1KreEkEALVw2T7hPhHJ8o4Ct3lxBoUcwaIeNmXU0pRdbEVPYUnWFv0Rmj4/gs67TpvLP+c64cEkc3uYNYOMn0Bfz8fsPoiGCj4/ikOWP7EBYUwCvb5GKmu2zskUpVUISMPhEdYvoCvsPRb3jDaBn77S7dQ4OYPSqRlbtLONvQbHQcn7QqZgiR1gamDfLcnPnC/ExfwFftLiE0yMKMoXFGR/FpWZkp1DW1sHKXrFzvak0vvsRP3vwzex6fS+iAfpCdbXQkYRKmLuDNLTbe2VvKjKEy86C7jUyKYkSfHmRvPWFfuV64RnY2lnvvJbHmJAqgoAAWLZIiLpxi6gL+2dGTnDrXJMOuPCRrYgoHy2rYeUJWrneZhx8msOGClcvr6uDhh43JI0zF1AV89e4SuocGcunAGKOj+IXrRyUSGRIo86O4kD7Rxr9lW9uFOI9pC3hDcwsffl7ONcMTCAkMMDqOX4gICWTO2D68vbeU0+eajI7jE87FtfHuMSXFs0GEKZm2gK87VEFto1WGXXlY1sRUmpUJlsgAABHKSURBVKw2lu0oMjqKT3hm5j00BF1w23x4OCxebEwgYSqmLeCrc0uJ6RZMZt9oo6P4lUHxkYxP60n21gJssnJ9l5yoquOvceNZ/+DvIMRRxFNTYckSyMoyNpwwBVMW8NpGKx8fKGfWiARZuMEAWRNTya+qY9OxKqOjmNpbjgnYhv30PsjMhMsug/x8Kd7CaeaoftOm2R8OH39eTqPVJt0nBpk5PJ6e4UG8LCvXd5rWmrd2FTMxPZo+UWFGxxEmZY4CfoHVuSUk9ghlXEpPo6P4pdCgAG7JSOajA+WUn20wOo4p5RadIe/kOVl8W3SJ6Qp4dV0T649Uct2oRCwWWbjBKPMmpNBi07JyfSct31FESKCFa79YuHjdOvtDiA5ot4ArpZ5TSlUopfadt+0RpVSxUmq34zHLvTG/8v6+MppbNNfLgq+GSouJ4NKBsby67QTNLTaj45hKo7WFVbklXDUsnu6hQUbHESbmTAv8BWBmK9uf1FqPdjzedW2s82Rnw5Yt8OmnkJZG1ZLnSe0VzvA+su6l0e7MTKX8bCMff15udBRT+eRABWfqm5k7LsnoKMLk2i3gWuv1wCkPZPmm7Gz7vBCNjfbXBQXcvfT3/OzkDln30gtMH9ybPlFhvCQXMztk2Y4i4rqHcEl/uYNYdE1X+sC/p5Ta4+hicc/VxIcfts8LcZ5wayNXv/I3txxOdEyARZGVmcKmY1UcrZCV651RWdPIusOV3DimDwFyDUd0UWcL+NNAP2A0UAo83taOSqlFSqkcpVROZWVlx47SxnwQgcVyF6C3uCXDvnL9yzI/ilNW7i6mxaaZO1a6T0TXdaqAa63LtdYtWmsb8Aww4SL7LtFaZ2itM2JjOzhZfRvzQSiZJ8JrxHQLYdaIeJbvKOJco6xc355lO4oYldSDAXGRRkcRPqBTBVwpdf7ik3OAfW3t2yWLF9vnhTiPLSxM5onwMndMSqOm0frlnYWidftLznCwrIb/kYuXwkWcGUb4KrAZGKSUKlJKLQQeU0rtVUrtAaYDP3JLuqws+7wQISHYgPKecVieeUZuNfYyY1OiGJbYnaWb8mWxh4t4Y3shwYEWmb9euEy7y9horee1svlZN2RpXVYWhS++ztQx9/LgzEF8d1p/jx1aOEcpxV2T0/jZsj1szqticj8ZXXGhhuYWVuwq5prh8USFy+LbwjVMcSfmO9GDAbhuhLRcvNX1oxKJjghm6aZ8o6N4pff3lXG2wcqt45ONjiJ8iCkKePD99zFjaBwpvcLb31kYIjQogHkTkvno83KKTte1/wV+5rXtJ0iJDiczvZfRUYQPMUUB/9Yl6TxzZ4bRMUQ75memopSSG3sucPzkObbkneLW8ckyf49wKVMUcGEOCT3CuHpYHK9tK6S+qcXoOF7jjZxCLAq5dV64nBRw4VJ3TU7nTH2zDCl0aG6x8WZOEZcP7k1c91Cj4wgfIwVcuNT4tJ4MS+zOc58dlyGFwJoDFZysbeTW8XLzmXA9KeDCpZRS3DM1nSMVtaw/ctLoOIbL3lpAfPdQpg/q4F3IQjhBCrhwuWtHJBLXPYT/bMgzOoqh8ipr2XDkJLdPTJG1W4VbyG+VcLngQAsLJqex4chJDpadNTqOYbK3niDQorhtgoz9Fu4hBVy4xe0TUggLCuDZDceNjmKI+qYW3swpZObweHpHysVL4R5SwIVbRIUHc3NGEit3l1BR438LH6/OLeFsg5U7MlONjiJ8mBRw4TZ3T0mn2Wbjpc3+dWOP1poXt+QzMK4bE9KjjY4jfJgUcOE26TERXDU0jhc3F1DrR3OF7y6sZl/xWe5w3JkqhLtIARdudd+0/pypb+bVrf6zYs/zG/PpFhLIHFl1R7iZFHDhVqOTo5jcrxf/+SyPRqvv315fXF3PO3tLuW18Mt1C2p2tWYgukQIu3O6+af0oP9vIip2+f3v9F9Pp3jUlzdAcwj9IARdud0n/GEb06cG/1+fRYvPd2+trG628uvUE1wyPJ6mnTH0s3E8KuHA7pRT3TevH8ZPn+GB/mdFx3OaN7YXUNFq5Z2pfo6MIPyEFXHjE1cPi6RsTwVOfHMXmg61wa4uN5zYeZ3xaT0YnRxkdR/gJKeDCIwIsiu9f0Z8DpWf58HPfa4V/sL+cotP1LLxEWt/Cc6SAC4+ZPaoPfWMj+MvHR3yqFa615p/rjpIeE8GMoXFGxxF+RAq48JgAi+KBKwZwsKyG932oL/yTgxXsLznLd6f1I0CWTBMeJAVceNR1IxPp37sbf/n4sE+0wrXWPPXJUZJ6hnHjmD5GxxF+Rgq48KgAi+IHVwzgcHkt7+wtNTpOl3129CS7C6u5b1o/gmTOb+Fh7f7GKaWeU0pVKKX2nbctWin1kVLqiONjT/fGFL7k2hEJDIzrxpMfHaa5xWZ0nC55as1REnqEyoLFwhDONBleAGZesO3nwBqt9QBgjeO1EE4JsCgevHoweSfP8dr2QqPjdNrWvCq25Z/i3kv7EhIYYHQc4YfaLeBa6/XAqQs23wAsdTxfCtzo4lzCx10xpDcT06P568eHTTlTodaaP394iNjIEG6bIAsWC2N0ttMuTmtdCuD42LutHZVSi5RSOUqpnMrKyk4eTvgapRS/nDWEk7VN/PvTY0bH6bCPD1SwPf80D1wxgNAgaX0LY7j9qovWeonWOkNrnREbKytzi6+MSo7i+lGJPLMhj7Iz5lm1x9pi44/vH6RvTAS3jpf1LoVxOlvAy5VSCQCOjxWuiyT8yYNXD8Jmgz9/eMjoKE5btqOIoxW1PDhzsIw8EYbq7G/fKmCB4/kCYKVr4gh/kxwdzt1T0li2o4gdBaeNjtOu+qYWnvz4MONSe3L1MLnrUhjLmWGErwKbgUFKqSKl1ELgD8AMpdQRYIbjtRCd8oMrBpDQI5RfvbUPq5cPK3z2szzKzzbyi2sGy3JpwnDOjEKZp7VO0FoHaa2TtNbPaq2rtNZXaK0HOD5eOEpFCKdFhATym+uHcqD0LC84FkTwRoWn6vj72qPMHBZPRposViyMJx14witcPSye6YNiefKjw5SeqTc6zjdorfn1yn0EKMVvZg81Oo4QgBRw4SWUUvx29nCsDQ08+tOnjY7zDR/sL2PtoUp+NGMgCT3CjI4jBCAFXHiRlF7h/KBoM+/1GsTK3d6zfmZto5VHVn3OkITu3DU5zeg4QnxJCrjwKveWbGNcTTG/emsfxdXe0ZXy5w8OUV7TwO/mDCdQhg0KLyK/jcKrBKJ58ug72GyaH7++2/BFkNcequCFTfksmJTGmBSZs014FyngwntkZ8OWLaR8uIptz9xD3NvL+c+GPMPiVNQ08NM3chkcH8nPrxlsWA4h2iIFXHiH7GxYtAgaGwGIKCvmTx/+g4NP/IsteVUej2OzaX7yRi7nmqw8NW+MzHcivJIUcOEdHn4Y6uq+timkqYGH1r/EfS/voPBUXRtf6B7/Xp/HhiMn+fV1wxgQF+nRYwvhLCngwjucONHq5rgzFdg0LFy6nZqGZo9EeX9fGY99cJBZI+KZN0EmqxLeSwq48A4prc+prVJS+GfWWI5VnuOB13a7fQWfHQWneOC1XYxOjuLxm0fL7fLCq0kBF95h8WIID//6tvBwWLyYKf1j+O3sYXxysIL7s3fSaG1xS4RjM+ew8G+fkBgVxrMLxhMWLP3ewrtJARfeISsLliyBkBD769RU++usLADmZ6byyPVD+fDzcr7z0g4aml1bxA+WnWX+kFsI1DaW3j2B6Ihgl35/IdxBCrjwHllZkJkJl10G+flfFu8v3DUlnd/NGcG6w5UsXLqd6romlxx2w5FK5j69GZtSvHjwTVJ6hbf/RUJ4ASngwlRun5jCn+aOYmveKWb9dQPbjnd+IkytNa9uO8Hdz29n/rHP2Pj0txj63jJIS7MPaxTCyymtPXenW0ZGhs7JyfHY8YTvyi2s5oHXdnHiVB3fm96f+6b171Cf9dGKGn69cj+bjlXxYFUO973yGOr8YYzh4V/rwhHCSEqpHVrrjG9slwIuzKq20cpvVu5n+c4ioiOCuWtyGndOSiUqvPX+a601h8treSOnkBc35xMWFMDPZg5m/q2Xolobxpiaau/KEcJgUsCFz9p2/BT/+vQYnxysIDTIwsikKEb26cHQxO602DRn6puprGlkzcEKjlbUYlFw09gkfn7NYGK6hYDFAq39HSgFNu9eIUj4h7YKeKARYYRwpQnp0UxIj+Zg2Vle21bI7sJqXtxSQJP1q+JrUfb9Fkwezsxh8cRGhnz1DVJSoKDgm9+4jbHpQngLKeDCZwyO784js4cB0Nxio6DqHCGBAXQPCyIyJBCLpY2bchYvts/DcmEf+OLFHkgtROdJARc+KSjAQv/eTs5h8sWFyoUL7ZNppabai7dcwBReTgq4EGAv1s88Y3++bp2hUYRwlowDF0IIk5ICLoQQJiVdKEJ8QbpOhMl0qYArpfKBGqAFsLY2TlEIIYR7uKIFPl1rfdIF30cIIUQHSB+4EEKYVFcLuAY+VErtUEotam0HpdQipVSOUiqnsrKyi4cTQgjxha4W8Cla67HANcD9SqlLL9xBa71Ea52htc6IjY3t4uGEEEJ8oUsFXGtd4vhYAawAJrgilBBCiPZ1uoArpSKUUpFfPAeuAva5KpgQQoiL68oolDhghWPV7kDgFa31+y5JJYQQol2dLuBa6zxglAuzCCGE6ACPLuiglKoEWpl42SkxgD+ON/fH8/bHcwb/PG9/PGfo+Hmnaq2/MQrEowW8K5RSOf54p6c/nrc/njP453n74zmD685bbuQRQgiTkgIuhBAmZaYCvsToAAbxx/P2x3MG/zxvfzxncNF5m6YPXAghxNeZqQUuhBDiPFLAhRDCpLyugCulZiqlDimljiqlft7K50OUUq87Pr9VKZXm+ZSu5cQ5/1gp9blSao9Sao1SKtWInK7W3nmft99cpZRWSpl+uJkz56yUusXx896vlHrF0xndwYnf8RSl1Fql1C7H7/ksI3K6klLqOaVUhVKq1SlGlN3fHP8me5RSYzt8EK211zyAAOAY0BcIBnKBoRfs813gX47ntwGvG53bA+c8HQh3PL/P7Ofs7Hk79osE1gNbgAyjc3vgZz0A2AX0dLzubXRuD533EuA+x/OhQL7RuV1w3pcCY4F9bXx+FvAeoIBMYGtHj+FtLfAJwFGtdZ7Wugl4Dbjhgn1uAJY6ni8DrlCOCVlMqt1z1lqv1VrXOV5uAZI8nNEdnPlZA/w/4DGgwZPh3MSZc/428A+t9Wn4cqZPs3PmvDXQ3fG8B1DiwXxuobVeD5y6yC43AC9quy1AlFIqoSPH8LYC3gcoPO91kWNbq/tora3AGaCXR9K5hzPnfL6F2P/XNrt2z1spNQZI1lq/7clgbuTMz3ogMFAptVEptUUpNdNj6dzHmfN+BJivlCoC3gW+75lohuro3/43eNuq9K21pC8c5+jMPmbi9PkopeYDGcBlbk3kGRc9b6WUBXgSuMtTgTzAmZ91IPZulGnY32ltUEoN11pXuzmbOzlz3vOAF7TWjyulJgEvOc7b5v54hulyLfO2FngRkHze6yS++Vbqy32UUoHY325d7G2Kt3PmnFFKXQk8DMzWWjd6KJs7tXfekcBwYJ1SKh97H+Eqk1/IdPb3e6XWullrfRw4hL2gm5kz570QeANAa70ZCMU+4ZMvc+pv/2K8rYBvBwYopdKVUsHYL1KuumCfVcACx/O5wCfacUXApNo9Z0dXwr+xF29f6BOFds5ba31Gax2jtU7TWqdh7/ufrbXOMSauSzjz+/0W9ovWKKVisHep5Hk0pes5c94ngCsAlFJDsBdwX19EdxVwp2M0SiZwRmtd2qHvYPSV2jauzB7GftX6Yce2R7H/8YL9B/smcBTYBvQ1OrMHzvljoBzY7XisMjqzJ877gn3XYfJRKE7+rBXwBPA5sBe4zejMHjrvocBG7CNUdgNXGZ3ZBef8KlAKNGNvbS8EvgN857yf9T8c/yZ7O/P7LbfSCyGESXlbF4oQQggnSQEXQgiTkgIuhBAmJQVcCCFMSgq4EEKYlBRwIYQwKSngQghhUv8fm5FZ1QOvfhUAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "def func(x,par):\n", " \"\"\" Another linear function: a + b x + c sin(x) (linear in the parameters!)\n", " \"\"\"\n", " return par[0]+par[1]*x+par[2]*np.sin(10*x)\n", "\n", "def deriv(x):\n", " \"\"\" Derivatives with respect to each parameter at input value(s)\n", " \"\"\"\n", " return [np.ones(len(x)),x,np.sin(10*x)]\n", "\n", "\n", "x=np.arange(0,1,0.1) # desired independent variables\n", "truepar=[13,13,13] # input parameters\n", "sig=1 # uncertainties on data points\n", "# simulate a data set with noise:\n", "y=func(x,truepar)+np.random.normal(0.,sig,size=len(x))\n", "plt.errorbar(x,y,sig,fmt='ro')\n", "\n", "# design matrix\n", "design=np.array(deriv(x)).T/sig\n", "#print(design)\n", "\n", "data=np.array(y)/sig\n", "ATA=np.dot(design.T,design)\n", "par=np.linalg.solve(ATA,np.dot(design.T,data))\n", "print('True parameters: ',truepar)\n", "print('Derived parameters: ',par)\n", "\n", "par=scipy.linalg.lstsq(design,data)[0]\n", "print('Derived parameters: ',par)\n", "\n", "x=np.arange(0,1,0.01)\n", "plt.plot(x,func(x,par))" ] } ], "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.11" } }, "nbformat": 4, "nbformat_minor": 1 }