Scipy两类拟合工具及lmfit估算非线性函数参数结果不一致求助

待估计函数 #

log(VA) = gamma - (1/eta)log[alpha*L^(-eta) + beta*K^(-eta)]

尝试用非线性最小二乘法估算上述函数的参数，先后使用Scipy-minimize、Scipy-curve_fit、lmfit-Model三种工具实现拟合，但三种工具返回的参数结果差异较大，无法排查问题原因，需要可行的解决方案或其他适用的求解方法。

测试数据与工具实现 #

基础测试数据 #

import numpy as np
from scipy.optimize import minimize, curve_fit
from lmfit import Model, Parameters

L  = np.array([0.299, 0.295, 0.290, 0.284, 0.279, 0.273, 0.268, 0.262, 0.256, 0.250])
K  = np.array([2.954, 3.056, 3.119, 3.163, 3.215, 3.274, 3.351, 3.410, 3.446, 3.416])
VA = np.array([0.919, 0.727, 0.928, 0.629, 0.656, 0.854, 0.955, 0.981, 0.908, 0.794])

Scipy-minimize实现 #

def f(param):
    gamma = param[0]
    alpha = param[1]
    beta  = param[2]
    eta   = param[3]
    VA_est = gamma - (1/eta)*np.log(alpha*L**-eta + beta*K**-eta)
    
    return np.sum((np.log(VA) - VA_est)**2)

bnds = [(1, np.inf), (0,1),(0,1),(-1, np.inf)]
x0 = (1,0.01,0.98, 1)

result = minimize(f, x0, bounds = bnds)

print(result.fun)
print(result.message)
print(result.x[0],result.x[1],result.x[2],result.x[3])

输出结果

0.30666062040617503
CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL
1.0 0.5587147011643757 0.9371430857380681 5.873041615873815

Scipy-curve_fit实现 #

def f(X, gamma, alpha, beta, eta):
    L,K = X
  
    return gamma - (1/eta) * np.log(alpha*L**-eta + beta*K**-eta)

p0 = 1,0.01,0.98, 1

res, cov = curve_fit(f, (L, K), np.log(VA), p0,  bounds = ((1,0,0,-1),(np.inf,1,1,np.inf)) )
gamma, alpha, beta, eta = res[0],res[1],res[2],res[3] 
gamma, alpha, beta, eta

输出结果

(1.000000000062141,
 0.26366547263939205,
 0.9804436474926481,
 13.449747863921704)

LMFIT-Model实现 #

def f(x, gamma, alpha, beta, eta):
    L = x[0]
    K = x[1]
    
    return gamma - (1/eta)*np.log(alpha*L**-eta + beta*K**-eta)

fmodel = Model(f)
params = Parameters()
params.add('gamma', value = 1,    vary=True, min = 1)
params.add('alpha', value = 0.01, vary=True, max = 1, min = 0)
params.add('beta',  value = 0.98, vary=True, max = 1, min = 0)
params.add('eta',   value = 1,    vary=True, min = -1)

result = fmodel.fit(np.log(VA), params, x=(L,K))
print(result.fit_report())

输出结果

[[Model]]
    Model(f)
[[Fit Statistics]]
    # fitting method   = leastsq
    # function evals   = 103
    # data points      = 10
    # variables        = 4
    chi-square         = 0.31749840
    reduced chi-square = 0.05291640
    Akaike info crit   = -26.4986758
    Bayesian info crit = -25.2883354
##  Warning: uncertainties could not be estimated:
    gamma:  at initial value
    gamma:  at boundary
    alpha:  at boundary
[[Variables]]
    gamma:  1.00000000 (init = 1)
    alpha:  1.3245e-13 (init = 0.01)
    beta:   0.20130064 (init = 0.98)
    eta:    447.960413 (init = 1)

问题原因与解决方案 #

核心问题原因 #

• 拟合目标函数存在参数不可识别性：当eta趋近于无穷大时，CES函数会退化为Leontief生产函数，此时alpha和beta的取值对拟合结果的影响会被大幅压缩，不同参数组合可以得到几乎相同的拟合残差，也就是目标函数存在非常多的局部最优点，三个优化器的收敛逻辑不同，自然会得到完全不同的局部最优解。代码里预设的等式约束既没有定义也没有传入优化函数，等于没有添加额外约束，进一步放大了解的自由度。
• 人为设置的gamma下界为1存在明显不合理性：log(VA)的最大值仅为log(0.981)≈-0.019，这个不符合数据特征的约束直接导致所有结果中gamma都被卡在边界1上，进一步加剧了参数识别的困难。
• 缺少CES生产函数的常规约束：通常CES函数默认假设规模报酬不变，也就是alpha + beta = 1，缺少这个约束会导致参数自由度太高，解不唯一。

可行解决方法 #

1. 修正不合理的约束：如果没有明确的理论依据证明gamma必须大于1，直接去掉gamma≥1的限制。
2. 添加规模报酬约束：增加alpha + beta = 1的等式约束传入优化器，减少参数自由度，大幅提升参数可识别性。
3. 改用全局优化算法：非线性最小二乘的局部优化算法非常依赖初值，可以改用scipy.optimize.dual_annealing这类全局优化算法寻找全局最优解，避免陷入局部最优。
4. 网格搜索校准初值：先对eta在合理区间（比如0到20）做网格搜索，固定eta时其余参数可以用线性方法估算，找到残差最小的eta作为优化初值，再做局部优化，结果稳定性会高很多。

内容的提问来源于stack exchange，提问作者Muhammet Rıdvan İNCE