在曲线拟合函数中实现线性方程

Question

我的教授向我们发送了这段代码，用作拟合曲线或您想要拟合的任何函数的事实。即使写入函数类型是线性的函数，我也看不到没有方程可以在拟合的函数中表达它。有人可以帮我在已经编写的代码中实现这个等式吗？对我来说很好。我已经尝试实现代码，可能在 try except 块之后执行 finally 调用，但没有任何效果。这是教授发给我的代码：

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

此代码是标准线性优化的逐步手工解决方案。一步一步的意思是，如果你把它写成矩阵和向量的形式，很容易看出会发生什么，但是所有的矩阵和向量操作都是在这里手动完成的。通过工作示例和与其他方法的比较，它看起来像这样 （有关更一般的矩阵方法的更多详细信息，请参阅例如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] ) )