{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "Parabolic fit\n", "==============\n", "\n", "#### Given two list of x and y coordinates, and a list of errors on y, find the parabola that minimizes the sum of distances (in the y direction) between the parabola 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\n",
"from sympy import *"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"xx, yy, aa, bb, cc = symbols('xx yy aa bb cc')"
]
},
{
"cell_type": "code",
"execution_count": 3,
"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$, $b$ and $c$ of the parabola $y = c x^2 + b x + a$ which minimize:\n",
"$ \\chi^2 = \\sum_i \\frac{(y_i - c x_i^2 - b x_i - a)^2}{\\sigma_i^2}$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The minimum condition can be derived expanding $(yy - cc \\cdot xx^2 - bb \\cdot xx - aa)^2$, and equating the derivatives with respect to $aa$, $bb$, $cc$ to zero. "
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"aa**2 + 2*aa*bb*xx + 2*aa*cc*xx**2 - 2*aa*yy + bb**2*xx**2 + 2*bb*cc*xx**3 - 2*bb*xx*yy + cc**2*xx**4 - 2*cc*xx**2*yy + yy**2"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"s = expand((yy - cc*xx**2 - bb*xx - aa)**2)\n",
"s"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Compute the three derivatives:"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2*aa + 2*bb*xx + 2*cc*xx**2 - 2*yy"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"s_aa = diff(s,aa)\n",
"s_aa"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2*aa*xx + 2*bb*xx**2 + 2*cc*xx**3 - 2*xx*yy"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"s_bb = diff(s,bb)\n",
"s_bb"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2*aa*xx**2 + 2*bb*xx**3 + 2*cc*xx**4 - 2*xx**2*yy"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"s_cc = diff(s,cc)\n",
"s_cc"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$a$, $b$ and $c$ 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} & \\sum_i \\frac{x_i^2}{\\sigma_i^2}\\\\\n",
"\\sum_i \\frac{x_i}{\\sigma_i^2} & \\sum_i \\frac{x_i^2}{\\sigma_i^2} & \\sum_i \\frac{x_i^3}{\\sigma_i^2}\\\\\n",
"\\sum_i \\frac{x_i^2}{\\sigma_i^2} & \\sum_i \\frac{x_i^3}{\\sigma_i^2} & \\sum_i \\frac{x_i^4}{\\sigma_i^2}\\\\\n",
"\\end{bmatrix}$\n",
"$\\begin{bmatrix} a\\\\b\\\\c \\end{bmatrix}$ = \n",
"$\\begin{bmatrix} \n",
"\\sum_i \\frac{y_i}{\\sigma_i^2} \\\\ \\sum_i \\frac{x_i y_i}{\\sigma_i^2}\\\\\\sum_i \\frac{x_i^2 y_i}{\\sigma_i^2}\n",
"\\end{bmatrix}$\n",
"\n",
"Compute numerically the entries of the left-hand-side matrix:"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"yerrSq = yerr*yerr\n",
"sum_one_over_yerrSq = (1./yerrSq).sum()\n",
"sum_x_over_yerrSq = (xdata/yerrSq).sum()\n",
"sum_x2_over_yerrSq = (xdata*xdata/yerrSq).sum()\n",
"sum_x3_over_yerrSq = (xdata*xdata*xdata/yerrSq).sum()\n",
"sum_x4_over_yerrSq = (xdata*xdata*xdata*xdata/yerrSq).sum()\n",
"sum_y_over_yerrSq = (ydata/yerrSq).sum()\n",
"sum_xy_over_yerrSq = (xdata*ydata/yerrSq).sum()\n",
"sum_x2y_over_yerrSq = (xdata*xdata*ydata/yerrSq).sum()"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([[ 20., 60., 220.],\n",
" [ 60., 220., 900.],\n",
" [ 220., 900., 3916.]])"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"mat = np.array([[sum_one_over_yerrSq,sum_x_over_yerrSq,sum_x2_over_yerrSq],\n",
" [sum_x_over_yerrSq,sum_x2_over_yerrSq,sum_x3_over_yerrSq],\n",
" [sum_x2_over_yerrSq,sum_x3_over_yerrSq,sum_x4_over_yerrSq]])\n",
"mat"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Invert the matrix:"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([[ 1.15 , -0.825 , 0.125 ],\n",
" [-0.825 , 0.66785714, -0.10714286],\n",
" [ 0.125 , -0.10714286, 0.01785714]])"
]
},
"execution_count": 10,
"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": 11,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([[ 1.00000000e+00, 3.66373598e-15, -9.02056208e-16],\n",
" [-3.10862447e-15, 1.00000000e+00, -2.31759056e-15],\n",
" [ 2.62012634e-14, 7.51620988e-14, 1.00000000e+00]])"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"np.dot(mat,mat_inv)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Compute the entries of the right-hand-side vector:"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([ 52. , 177.6, 685.6])"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"noti = np.array([sum_y_over_yerrSq,sum_xy_over_yerrSq,sum_x2y_over_yerrSq])\n",
"noti"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Compute $a$, $b$ and $c$:"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([-1.02 , 2.25428571, -0.28571429])"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"np.dot(mat_inv,noti)"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"-1.0199999999999534\n",
"2.2542857142857002\n",
"-0.2857142857142829\n"
]
}
],
"source": [
"a, b, c = np.dot(mat_inv,noti)[0],np.dot(mat_inv,noti)[1],np.dot(mat_inv,noti)[2]\n",
"print(a)\n",
"print(b)\n",
"print(c)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Plot, adding extra points for the parabola:"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[