我很难找到用 lmfit 拟合振荡函数的良好初始射击

如何解决我很难找到用 lmfit 拟合振荡函数的良好初始射击

我正在解决的问题包括拟合耦合 Van der Pol 振荡器的参数，如引用的论文 link to paper 中所述。我想尝试使用经典的最小二乘法，而不是使用作者提出的伴随方法。为此，我决定用“合成”案例“验证”。让我分享一下这部分的摘录

残差定义

def residual_ode(params,t,data,ID,theta,sol_0):
    """
    Here I define the residual function to minimize.
    """
    A,B,D           = params[ID[0]],params[ID[1]],params[ID[2]]
    c_til,d_len,x_0 = theta
    u0              = ode_solver(A,D,sol_0,t)[0]
    
    #plt.plot(t,u0,'g-')
    
    return u0 - data

以及 ODE 的求解器

def ode_sys(sol,A,D):
    """
    This function returns the order reduction for the Van der Pol's ODEs
    A,D: are scalar parameters to be "determined"
    x_sol,y_sol: are displacements
    u_sol,v_sol: are velocities
    You need to convert the 2-second order ODEs into 4-first order ODEs
    """
    
    x_sol,u_sol,y_sol,v_sol = sol
    
    x_der = u_sol
    u_der = (A - B * (1 + np.power(x_sol,2))) * u_sol + (0.5 *D-1) * x_sol + A * v_sol
    y_der = v_sol
    v_der = (A - B * (1 + np.power(y_sol,2))) * v_sol - (0.5 *D+1) * y_sol + A * u_sol
    
    return [x_der,u_der,y_der,v_der]

def ode_solver(A,t):
    """
    This function integrates the coupled ODEs that represent the Van der Pol oscillators
    A     :: alpha
    B     :: beta
    D     :: delta
    sol_0 :: initial condition [x(0),x'(0),y(0),y'(0)]
    t     :: time array
    u0    :: computed glottal flow
    sol   :: displacements of both vocal folds vs time
    """
    # Here I define constants of the model (physiological props.)
    c_til = 50            # air particle velocity,m/s
    d_len = 0.0175        # length of vocal folds,m
    x_0   = 0.001         # half glottal width at rest position,m
    
    sol   = odeint(ode_sys,args=(A,D))
    u0    = c_til * d_len * (2 * x_0 + sol[:,0] + sol[:,2])
    
    return u0,sol

第一步，定义常量和固定参数：

# Here I define constants of the model
c_til = 50            # air particle velocity,m/s
d_len = 0.0175        # length of vocal folds,m
x_0   = 0.001         # half glottal width at rest position,m
sol_0 = [0,0.1,0.1] # Initial state from papers (reference) (x(0),y'(0))
theta = [c_til,x_0]

if(dataset == 'Synthetic'):
    # Here I define a synthetic dataset just to test the fit
    A     = 0.50 # Reference value
    B     = 0.32 # Reference value
    D     = 0.00 # Reference value

    t     = np.linspace(0,60*pi,1000)
    data,tgt = ode_solver(A,t)
    data  = data + 0.1*np.random.normal(size=t.shape) # Add a noise to see the robustness.
else:
    # Here I upload a given audio file
    t,signal,glot_flow,sr = load_audio('target.WAV')
    t    = t[:900]
    data = glot_flow[:900] * 1.0e-01
else:
    # Here I upload a given audio file
    t,sr = load_audio('target.WAV')
    t    = t[:900]
    data = glot_flow[:900] * 1.0e-01

第 2 步，声明要拟合的参数，加上最小-最大限制和其他道具。

# Define the initial shoot and ranges for the parameters to fit
ID = ['A','B','D']
if(dataset == 'Synthetic'):
    i0 = [0.50,0.50,0.50] # This one works for synthetic
else:
    i0 = [0.50,0.50]

vi,vf = [0,0],[1,1,1] # According to code,they are always between (0,1)

第 3 步，定义并调用 lmfit 最小化工具

# According to lmfit,you must define ranges for parameters and initial guess
params = lmf.Parameters()
params.add(ID[0],value=i0[0],min=vi[0],max=vf[0])  # alpha
params.add(ID[1],value=i0[1],min=vi[1],max=vf[1])  # beta
params.add(ID[2],value=i0[2],min=vi[2],max=vf[2])  # Delta

# lmfit minimizer 
foo = lmf.Minimizer(residual_ode,params,fcn_kws={'t': t,'data':data,'ID':ID,'theta':theta,'sol_0':sol_0})

result = foo.minimize(method='leastsq')
lmf.report_fit(result)

如前所述，我使用合成数据验证拟合脚本。 lmfit 报告显示了非常好的结果，低可变性，

Fitting ODEs running on lmfit v1.0.0
Analyzed dataset Synthetic
[[Fit Statistics]]
    # fitting method   = leastsq
    # function evals   = 101
    # data points      = 1000
    # variables        = 3
    chi-square         = 9.96375601
    reduced chi-square = 0.00999374
    Akaike info crit   = -4602.80117
    Bayesian info crit = -4588.07790
[[Variables]]
    A:  0.50089718 +/- 1.5368e-07 (0.00%) (init = 0.5)
    B:  0.32140460 +/- 2.1184e-05 (0.01%) (init = 0.5)
    D:  0.02961727 +/- 4.6440e-05 (0.16%) (init = 0.5)
[[Correlations]] (unreported correlations are < 0.100)
    C(B,D) = -0.857
    C(A,B) =  0.644
    C(A,D) = -0.373

受这些 GOOD 结果的激励，我决定转向一个“现实”的案例，请看一下健身报告：

Fitting ODEs running on lmfit v1.0.0
Analyzed dataset real
[[Fit Statistics]]
    # fitting method   = leastsq
    # function evals   = 41
    # data points      = 900
    # variables        = 3
    chi-square         = 0.98750335
    reduced chi-square = 0.00110090
    Akaike info crit   = -6127.47314
    Bayesian info crit = -6113.06595
[[Variables]]
    A:  3.6795e-06 +/- 51769552.4 (1406957752181724.50%) (init = 0.5)
    B:  0.99999944 +/- 2.6702e+08 (26701595662.24%) (init = 0.5)
    D:  0.21433002 +/- 2.4401e+11 (113848637625094.62%) (init = 0.5)
[[Correlations]] (unreported correlations are < 0.100)
    C(A,B) =  1.000
    C(B,D) = -0.964

不幸的是，这里有很多变化。让我也分享一些说明问题的数字。

、

。