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

问题描述 投票:0回答:1

我想获得 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()
python math statistics statsmodels numerical-methods
1个回答
0
投票

根据定义,

l_sq()
函数需要 4 个参数。

def l_sq(model, x, data, par)

然而,当它被调用时,只传递了 3 个值。

F = func(ff, ts, da)
© www.soinside.com 2019 - 2024. All rights reserved.