试图在scipy中获得曲线拟合函数的直觉

问：

我想在scipy函数scipy.optimize.curve_fit（）背后获得一些直觉。下面，我写了一个简单的版本，试图最小化最小二乘函数的总和（权重都是 1），模型是 sin（omega*t）。我采用牛顿-拉夫森方法（中心差分近似）来求导数的零点。下面的代码应显示一些初始值猜测的拟合（根查找）失败（请参阅循环），并显示初始值与（已知）参数值的成功拟合，这些初始值与有根据的猜测足够接近。然后采用最快（迭代次数最少）的收敛结果来绘制图，显示人工数据和最佳拟合曲线。

我一直收到错误：

TypeError: l_sq() missing 1 required positional argument: 'par'

这是我目前所拥有的：

import numpy as np
import matplotlib.pyplot as plt

def my_func(t, freq):
    ''' model function, unit amplitude.'''
    return np.sin(freq * t)

def Newton_Raphson(func, x0, dx, args, eps=1.e-9, Nmax=20):
    ''' Newton-Raphson root finder function. Default precision and maximum number of iterations.'''
    ff = args[0]
    ts = args[1]
    da = args[2]
    for it in range(Nmax + 1):
        F = func(ff, ts, da)
        if ( abs (F) <= eps ): # Converged?
            print ( " Root found , f(root) = ",F, " , eps = ",eps, ", iterations = ",it)
        df = (func(ff, ts, da, x0+dx/2) + func(ff, ts, da, x0-dx/2))/x0 # Central difference
        dx = -F / df
        x0 += dx # New guess
    if it >= Nmax:
        print ( "Newton Failed for Nmax = " , Nmax)
    return x0, it

def l_sq(model, x, data, par):
    ''' Least-Squares objective function to solve for fitting.'''
    return np.sum((data+model(x,par))*np.gradient(model(x,par),x))


# artificial data first
par   = 4.5 # [Hz] YR: correct
times = np.linspace(0,1,1000) # [s] YR: correct
data  = my_func(times, par) # noise-free data first, at set frequency

# small amount of noise added, just 10% error
data += np.random.normal(scale=0.1,size=len(data)) # YR: correct

# find the root of the first derivative of the least_square function - extremum
step = 1.e-4 # YR: correct
roots = []
nits  = []
for val in range(-7,7): # YR: range as intended, correct
    p0 = par+0.01*val*par # YR: correct, get single figure percent changes only
    root, nit = Newton_Raphson(l_sq, p0, step, args=(my_func, times, data))
    print('root = ', root)
    roots.append(root)
    nits.append(nit)
root = roots[np.argmin(nits)] # YR: lowest n iterations, best fit from initial value series, important for test.

# plot data and best fit curve. YR: correct plotting
plt.plot(times,data,'r.') # red
plt.plot(times, my_func(times, root),'b') # blue
plt.xlabel('time [s]')
plt.ylabel('sin($\omega$t)')
plt.show()

以下是更新后的代码：

import numpy as np
import matplotlib.pyplot as plt

def my_func(t, freq):
    ''' Model function, unit amplitude.'''
    return np.sin(freq * t)

def Newton_Raphson(func, x0, dx, args, eps=1.e-9, Nmax=20):
    ''' Newton-Raphson root finder function. Default precision and maximum number of iterations.'''
    ff = args[0]
    ts = args[1]
    da = args[2]
    for it in range(Nmax + 1):
        F = l_sq(ts, da, x0, my_func)  # Calculate F using l_sq function
        if abs(F) <= eps:  # Converged?
            print("Root found , f(root) = ", F, ", eps = ", eps, ", iterations = ", it)
            break  # Exit the loop when converged
        df = (l_sq(ts, da, x0 + dx / 2, my_func) - l_sq(ts, da, x0 - dx / 2, my_func)) / dx  # Central difference
        dx = -F / df
        x0 += dx  # New guess
    if it >= Nmax:
        print("Newton Failed for Nmax = ", Nmax)
    return x0, it

def l_sq(x, data, par, model):
    ''' Least-Squares objective function to solve for fitting.'''
    return np.sum((data - model(x, par))**2)

# Artificial data first
par = 4.5  # [Hz]
times = np.linspace(0, 1, 1000)  # [s]
data = my_func(times, par)  # Noise-free data first, at set frequency

# Add a small amount of noise, just 10% error
data += np.random.normal(scale=0.1, size=len(data))

# Find the root of the first derivative of the least_square function (extremum)
step = 1.e-4
roots = []
nits = []

for val in range(-7, 7):
    p0 = par + 0.01 * val * par
    root, nit = Newton_Raphson(l_sq, p0, step, args=(my_func, times, data))
    print('root = ', root)
    roots.append(root)
    nits.append(nit)

# Find the best fit from the initial value series
root = roots[np.argmin(nits)]

# Plot data and best-fit curve
plt.plot(times, data, 'r.')  # Red dots for data
plt.plot(times, my_func(times, root), 'b')  # Blue curve for best fit
plt.xlabel('time [s]')
plt.ylabel('sin($\omega$t)')
plt.show()