Implementing linear equation in curve-fitting function

Question

my professor sent us this code to use as a defacto for fitting curves or whatever function you would like to fit. even though is written into the function that the type of function is linear I can't see no equation to express that in the function for the fit. Can someone help me about implementing this equation in the already written code? would be great for me. I've tried implementing the code maybe doing a finally call after the try except block, but nothing is working.Here it is the code that the professor has sent me:

def fit_pesato(x,y,yerr):
  #N.B.: la funizione utilizza come formula della retta y=a+bx !
  if len(x)!=len(y) or len(yerr)!=len(y):
    raise Exception("Le liste di input non sono della stessa lunghezza!")
  try:
    chi_quadro = 0
    y_err2 =[i**2 for i in yerr]
    W,Y,X = 0,0,0
    XX, XY= 0,0
    for i in range(len(x)):
      W += 1/y_err2[i]
      X += x[i]/y_err2[i]
      Y += y[i]/y_err2[i]
      XX += (x[i]**2)/y_err2[i]
      XY += (x[i]*y[i])/y_err2[i]
    delta = W*XX-(X**2)
    A_mc = (XX*Y-X*XY)/delta
    B_mc = (W*XY-X*Y)/delta
    sigmaAmc=np.sqrt(XX/delta)
    sigmaBmc=np.sqrt(W/delta)
    for i in range(len(x)):
      chi_quadro += ((y[i]-A_mc-x[i]*B_mc)/yerr[i])**2
    return ['A_mc:',A_mc,'B_mc:',B_mc,'sigmaA:', sigmaAmc,'sigmaB:', sigmaBmc,'chi^2:', chi_quadro]
  except ZeroDivisionError:
    return ['A_mc:',None,'B_mc:',None, 'sigmaA:', None,'sigmaB', None,'chi^2:',None]
    #questa ultima linea di codice serve a fornire un risultato in casi 

Answer 1

This code is the step by step handmade solution to the standard linear optimization. With step-by-step I mean that it is easy to see what happens if you write it down in matrix and vector form, but all the matrix and vector operations are done manually here. With working example and comparison to other methods, it looks like this (For more details on a more general matrix approach, see e.g. here )

"""
I will assume that x, y and ySig are numpy arrays so,
we can avoid all the loops.

The code follws the main idea that we want to minimize the error in
(y - (c1 + c2 x ) ). Here actually even weighted. So, when expressing
in matrix form this is
( y - A a ).T W ( y - A a ) = min,
where a = (c1, c2) and A.T = ( ( 1, ..., 1 ),( x1, ..., xn ) ),
i.e. A is a n by 2 matrix. A.T is the transposed matrix
setting the derivative to zero results in
A.T W A a = A.T W y
note that W = diag( 1 / sigma_1^2, ..., 1/ sigma_n^2 ), 
i.e. the inverse of the covariance matrix.( here, no correlation is assumed)
Evaluating A.T W A we get

( ( XX, X ),( X, S) ).T with:

    XX = sum x_i^2 / sig_i^2
    X = sum x_i / sig_i^2
    S = sum 1 / sig_i^2

obviously we need the inverse, which is P = ((S, -X),(-X, XX))/ det,
with det = XX * S - X^2
we then have a = P A.T W y 
and A.T W y = ( eta, xi) with
Y = sum y_i / sig_i^2 and
XY  = sum x_i y_i / sig_i^2 )
in general a = J y so if the covariance of y is Vy the covariance of a
is Va = J Vy J.T
here Va = ( A.T W A )^(-1) A.T W Vy W.T A (A.T W.T A )^(-1)
plugging in and noting that W = Vy^(-1)  and W.T = W we have
Va = P

Now putting the variable names from above
"""

import numpy as np


def fit_weighted( x, y, ySig ):
    #N.B.: la funizione utilizza come formula della retta y=a+bx !
    if (
        ( len( x ) != len( y ) ) or 
        ( len( ySig ) != len( y ) )
    ):
        raise ValueError (
            "Input lists of unequal length!"
        )
    try:
        XX = np.dot( x**2, 1 / ySig**2 )
        X = np.dot( x, 1 / ySig**2 )
        S = np.dot( 1 / ySig , 1 / ySig )

        Y = np.dot( 1 / ySig, y / ySig )
        XY = np.dot( x / ySig, y / ySig )

        det = XX * S - X**2  ### problematic for small det

        amc = ( XX * Y - X * XY ) / det
        bmc = ( S * XY - X * Y ) / det
        sigmaAmc = np.sqrt( XX / det )
        sigmaBmc = np.sqrt( S / det )
        ###chi2
        C1 = amc * np.ones( len( x ) )
        C2 = bmc * x
        delta = y - C1 - C2
        chi2 = np.dot( delta / ySig, delta / ySig )
        """
        Note, we already have P and could provide the covariance matrix 
        as output
        """
        return [
            'A_mc:',amc,'B_mc:',bmc,
            'sigmaA:', sigmaAmc,'sigmaB:', sigmaBmc,
            'chi^2:', chi2 
        ]
    except ZeroDivisionError:
        nan = float( 'nan' )
        return [
            'A_mc:', nan,'B_mc:', nan,
            'sigmaA:', nan,'sigmaB', nan,
            'chi^2:', nan
        ]

if __name__ == "__main__":
    """
    showing that we get the same results with curve_fit
    """
    from scipy.optimize import curve_fit
    c1 = 0.34
    c2 = 1.233
    xl = np.linspace( -2,3, 15 )
    yl = c1 + c2 * xl

    s = np.random.normal( size=len( xl ), scale=0.1 )
    yl = yl + s
    sa = np.abs( s )

    result = fit_weighted( xl, yl, sa )
    print( result )
    cfresult = curve_fit(
        lambda x, a, b: a + b * x,
        xl, yl,
        sigma=sa, absolute_sigma=True
    )
    print( cfresult[0] )
    print( np.sqrt( cfresult[1][0,0] ) )
    print( np.sqrt( cfresult[1][1,1] ) )