{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La formula per il fit di un set di dati con una retta usando SymPy\n",
    "===================================================================\n",
    "<P>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La somma degli scarti di dati $(x_i,y_i)$ (pesati con l'errore $\\sigma_i$ sugli $y_i$) rispetto una linea $y = b \\cdot x + a$ è:\n",
    "$$ \\chi^2 = \\sum_i \\frac{(y_i -  b \\cdot x_i - a)^2}{\\sigma_i^2} $$\n",
    "\n",
    "Per trovare il minimo di questa espressione al variare di $b$ e $a$, cerchiamo gli zeri contemporanei delle derivate di $\\chi^2$ rispetto ad $a$ e rispetto a $b$, generalizzando in modo ovvio la formula per gli estremi in una variabile.</br>\n",
    "È sufficiente limitarsi al caso di tre sole misure e poi generalizzare. \n",
    "\n",
    "__Domanda:__ Perchè tre misure e non una oppure due?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "x1,x2,x3,y1,y2,y3,a,b,sq1,sq2,sq3 = symbols('x1 x2 x3 y1 y2 y3 a b sq1 sq2 sq3')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "chi2 = (y1 - b*x1 - a)**2/sq1 + (y2 - b*x2 - a)**2/sq2 + (y3 - b*x3 - a)**2/sq3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Derivata rispetto ad a:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{2 a + 2 b x_{3} - 2 y_{3}}{sq_{3}} + \\frac{2 a + 2 b x_{2} - 2 y_{2}}{sq_{2}} + \\frac{2 a + 2 b x_{1} - 2 y_{1}}{sq_{1}}$"
      ],
      "text/plain": [
       "(2*a + 2*b*x3 - 2*y3)/sq3 + (2*a + 2*b*x2 - 2*y2)/sq2 + (2*a + 2*b*x1 - 2*y1)/sq1"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "chi2_prime_a = diff(chi2,a)\n",
    "chi2_prime_a"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Espandiamo il risultato per poi estrarre i coefficienti di a e b:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{2 a}{sq_{3}} + \\frac{2 a}{sq_{2}} + \\frac{2 a}{sq_{1}} + \\frac{2 b x_{3}}{sq_{3}} + \\frac{2 b x_{2}}{sq_{2}} + \\frac{2 b x_{1}}{sq_{1}} - \\frac{2 y_{3}}{sq_{3}} - \\frac{2 y_{2}}{sq_{2}} - \\frac{2 y_{1}}{sq_{1}}$"
      ],
      "text/plain": [
       "2*a/sq3 + 2*a/sq2 + 2*a/sq1 + 2*b*x3/sq3 + 2*b*x2/sq2 + 2*b*x1/sq1 - 2*y3/sq3 - 2*y2/sq2 - 2*y1/sq1"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "chi2_prime_a_expand = expand(chi2_prime_a)\n",
    "chi2_prime_a_expand"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle a \\left(\\frac{2}{sq_{3}} + \\frac{2}{sq_{2}} + \\frac{2}{sq_{1}}\\right) + b \\left(\\frac{2 x_{3}}{sq_{3}} + \\frac{2 x_{2}}{sq_{2}} + \\frac{2 x_{1}}{sq_{1}}\\right) - \\frac{2 y_{3}}{sq_{3}} - \\frac{2 y_{2}}{sq_{2}} - \\frac{2 y_{1}}{sq_{1}}$"
      ],
      "text/plain": [
       "a*(2/sq3 + 2/sq2 + 2/sq1) + b*(2*x3/sq3 + 2*x2/sq2 + 2*x1/sq1) - 2*y3/sq3 - 2*y2/sq2 - 2*y1/sq1"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "collect(chi2_prime_a_expand,(a,b))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Derivata rispetto a b:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle - \\frac{2 x_{3} \\left(- a - b x_{3} + y_{3}\\right)}{sq_{3}} - \\frac{2 x_{2} \\left(- a - b x_{2} + y_{2}\\right)}{sq_{2}} - \\frac{2 x_{1} \\left(- a - b x_{1} + y_{1}\\right)}{sq_{1}}$"
      ],
      "text/plain": [
       "-2*x3*(-a - b*x3 + y3)/sq3 - 2*x2*(-a - b*x2 + y2)/sq2 - 2*x1*(-a - b*x1 + y1)/sq1"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "chi2_prime_b = diff(chi2,b)\n",
    "chi2_prime_b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Espandiamo e estraiamo i coefficienti di a e b:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{2 a x_{3}}{sq_{3}} + \\frac{2 a x_{2}}{sq_{2}} + \\frac{2 a x_{1}}{sq_{1}} + \\frac{2 b x_{3}^{2}}{sq_{3}} + \\frac{2 b x_{2}^{2}}{sq_{2}} + \\frac{2 b x_{1}^{2}}{sq_{1}} - \\frac{2 x_{3} y_{3}}{sq_{3}} - \\frac{2 x_{2} y_{2}}{sq_{2}} - \\frac{2 x_{1} y_{1}}{sq_{1}}$"
      ],
      "text/plain": [
       "2*a*x3/sq3 + 2*a*x2/sq2 + 2*a*x1/sq1 + 2*b*x3**2/sq3 + 2*b*x2**2/sq2 + 2*b*x1**2/sq1 - 2*x3*y3/sq3 - 2*x2*y2/sq2 - 2*x1*y1/sq1"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "chi2_prime_b_expand = expand(chi2_prime_b)\n",
    "chi2_prime_b_expand"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle a \\left(\\frac{2 x_{3}}{sq_{3}} + \\frac{2 x_{2}}{sq_{2}} + \\frac{2 x_{1}}{sq_{1}}\\right) + b \\left(\\frac{2 x_{3}^{2}}{sq_{3}} + \\frac{2 x_{2}^{2}}{sq_{2}} + \\frac{2 x_{1}^{2}}{sq_{1}}\\right) - \\frac{2 x_{3} y_{3}}{sq_{3}} - \\frac{2 x_{2} y_{2}}{sq_{2}} - \\frac{2 x_{1} y_{1}}{sq_{1}}$"
      ],
      "text/plain": [
       "a*(2*x3/sq3 + 2*x2/sq2 + 2*x1/sq1) + b*(2*x3**2/sq3 + 2*x2**2/sq2 + 2*x1**2/sq1) - 2*x3*y3/sq3 - 2*x2*y2/sq2 - 2*x1*y1/sq1"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "collect(chi2_prime_b_expand,(a,b))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A questo punto è facile scrivere l'equazione matriciale che determina $a$ e $b$:\n",
    "\n",
    "$$\n",
    "\\begin{bmatrix}\n",
    "\\sum \\frac{1}{\\sigma_i^2} & \\sum \\frac{x_i}{\\sigma_i^2} \\\\\n",
    "\\sum \\frac{x_i}{\\sigma_i^2} & \\sum \\frac{x^2_i}{\\sigma_i^2}\n",
    "\\end{bmatrix}\n",
    "\\begin{bmatrix} a \\\\ b\n",
    "\\end{bmatrix}\n",
    "= \n",
    "\\begin{bmatrix} \n",
    "\\sum \\frac{y_i}{\\sigma_i^2}  \\\\ \\sum \\frac{x_i y_i}{\\sigma_i^2} \n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\begin{bmatrix}\n",
    "Sinv & Sx \\\\\n",
    "Sx & Sxsq\n",
    "\\end{bmatrix}\n",
    "\\begin{bmatrix} a \\\\ b\n",
    "\\end{bmatrix}\n",
    "= \n",
    "\\begin{bmatrix} \n",
    "Sy \\\\ Sxy\n",
    "\\end{bmatrix}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{a: (-Sx*Sxy + Sxsq*Sy)/(Sinv*Sxsq - Sx**2),\n",
       " b: (Sinv*Sxy - Sx*Sy)/(Sinv*Sxsq - Sx**2)}"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a, b, Sinv, Sx, Sxsq, Sy, Sxy = symbols('a b Sinv Sx Sxsq Sy Sxy')\n",
    "\n",
    "eq1 = Sinv*a + Sx*b - Sy\n",
    "eq2 = Sx*a  + Sxsq*b -Sxy\n",
    "\n",
    "solve((eq1,eq2),a,b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nel caso particolare in cui tutti i $\\sigma_i$ sono uguali, definendo $\\bar{u} = \\sum u_i /N$ ($u = 1,x,y,xy$)\n",
    "l'equazione che determina $a, b$ diventa:<BR>\n",
    "    \n",
    "$$\n",
    "M \\begin{bmatrix} a \\\\ b\n",
    "\\end{bmatrix}\n",
    "=\n",
    "\\begin{bmatrix}\n",
    "1 & \\overline{x_i} \\\\\n",
    "\\overline{x} & \\overline{x^2}\n",
    "\\end{bmatrix}\n",
    "\\begin{bmatrix} a \\\\ b\n",
    "\\end{bmatrix}\n",
    "= \n",
    "\\begin{bmatrix} \n",
    "\\overline{y}  \\\\ \\overline{xy} \n",
    "\\end{bmatrix},\n",
    "\\qquad\n",
    "M^{-1} = \\frac{1}{\\Delta} \\begin{bmatrix}\n",
    "\\overline{x^2} & -\\overline{x_i} \\\\\n",
    "-\\overline{x} & 1\n",
    "\\end{bmatrix},\\,\\, \\Delta=\\overline{x^2}-\\overline{x}^2\n",
    "$$\n",
    "<BR>\n",
    "con soluzione $a = (\\overline{x^2}\\cdot \\overline{y}-\\overline{x}\\cdot\\overline{xy})/\\Delta$, \n",
    "    $b = (\\overline{xy}-\\overline{x}\\cdot\\overline{y})/\\Delta$\n",
    "<BR>\n",
    "È facile mostrare che il punto $(\\overline{x},\\overline{y})$ appartiene alla retta di fit cioè che\n",
    "    $ \\overline{y} = b \\cdot \\overline{x} +a $.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}