{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "Linear fit\n", "===========\n", "\n", "#### Given two list of x and y coordinates, and a list of errors on y, find the straight line that minimizes the sum of distances (along the y direction) between the line and the(x,y) points.\n", "
"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\n",
"import numpy as np"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"xdata = np.array([1, 2, 3, 4, 5])\n",
"ydata = np.array([1., 2.3, 3., 3.7, 3.])\n",
"yerr = np.array([0.5, 0.5, 0.5, 0.5, 0.5])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We want to find the parameters $a$ and $b$ of the line $y = b x + a$ which minimize:\n",
"$ \\chi^2 = \\sum_i \\frac{(y_i - b x_i - a)^2}{\\sigma_i^2}$\n",
" \n",
"$a$ and $b$ are the solutions of the linear system:\n",
"$\\begin{bmatrix}\n",
"\\sum_i \\frac{1}{\\sigma_i^2} & \\sum_i \\frac{x_i}{\\sigma_i^2}\\\\\n",
"\\sum_i \\frac{x_i}{\\sigma_i^2} & \\sum_i \\frac{x_i^2}{\\sigma_i^2}\n",
"\\end{bmatrix}$\n",
"$\\begin{bmatrix} a\\\\b \\end{bmatrix}$ = \n",
"$\\begin{bmatrix} \\sum_i \\frac{y_i}{\\sigma_i^2} \\\\ \\sum_i \\frac{x_i y_i}{\\sigma_i^2}\\end{bmatrix}$\n",
"\n",
"Compute the entries which depend on the data:"
]
},
{
"cell_type": "code",
"execution_count": 38,
"metadata": {},
"outputs": [],
"source": [
"yerrSq = yerr*yerr\n",
"sum_one_over_yerrSq = (1./yerrSq).sum()\n",
"sum_x_over_yerrSq = (xdata/yerrSq).sum()\n",
"sum_xSq_over_yerrSq = (xdata*xdata/yerrSq).sum()\n",
"sum_y_over_yerrSq = (ydata/yerrSq).sum()\n",
"sum_xy_over_yerrSq = (xdata*ydata/yerrSq).sum()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Construct the left-hand-side matrix"
]
},
{
"cell_type": "code",
"execution_count": 39,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([[ 20., 60.],\n",
" [ 60., 220.]])"
]
},
"execution_count": 39,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"mat = np.array([[sum_one_over_yerrSq,sum_x_over_yerrSq],[sum_x_over_yerrSq,sum_xSq_over_yerrSq]])\n",
"mat"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Compute the inverse matrix:"
]
},
{
"cell_type": "code",
"execution_count": 40,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([[ 0.275, -0.075],\n",
" [-0.075, 0.025]])"
]
},
"execution_count": 40,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"mat_inv = np.linalg.inv(mat)\n",
"mat_inv"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Check"
]
},
{
"cell_type": "code",
"execution_count": 41,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([[1., 0.],\n",
" [0., 1.]])"
]
},
"execution_count": 41,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"np.dot(mat,mat_inv)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Compute the right-hand-side vector"
]
},
{
"cell_type": "code",
"execution_count": 42,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([ 52. , 177.6])"
]
},
"execution_count": 42,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"noti = np.array([sum_y_over_yerrSq,sum_xy_over_yerrSq])\n",
"noti"
]
},
{
"cell_type": "code",
"execution_count": 43,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([0.98, 0.54])"
]
},
"execution_count": 43,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"np.dot(mat_inv,noti)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Extract a and b multiplying both sides by the inverse:"
]
},
{
"cell_type": "code",
"execution_count": 44,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0.980000000000004\n",
"0.54\n"
]
}
],
"source": [
"a, b = np.dot(mat_inv,noti)[0],np.dot(mat_inv,noti)[1]\n",
"print(a)\n",
"print(b)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Plot the result:"
]
},
{
"cell_type": "code",
"execution_count": 45,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[