Scipy Curve_fit：我将如何改进这种拟合？

Question

我一直在研究一个标准势，我试图使其适合给定的模型： ax² - bx³ + lx⁴

拟合的 x 和 y 值也是从代码生成的，x 值由

numpy.linspace

生成。我已经限制了 a、b、c 参数，使它们始终为正。我需要拟合来模拟数据，至少局部最大值的高度和全局最小值的位置是准确的。相反，这就是我得到的（蓝色是实际数据，虚线是给定模型的拟合数据）：

这是我的代码的相关部分：

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from scipy.interpolate import *

v =246.
x = np.linspace(0.1, 270.1, 27001)  
Xdata = x/v

def func(x,a,b,l):
    return (a*(x**2)) - (b*(x**3)) + ((l)*(x**4))

temp = np.linspace(80.,80.,1)
b = np.zeros_like(temp)
a = np.zeros_like(temp)
l = np.zeros_like(temp)
Vfit = pot_giver(temp[0])   #External func.
Ydata = (Vfit - Vfit[0])/(pow(v,4))
popt, pcov = curve_fit(func, Xdata, Ydata,bounds=((0.,0.,0.), (np.inf,np.inf,np.inf)))
b[0],a[0],l[0] = popt

plt.plot(Xdata,Ydata)
plt.plot(Xdata, func(Xdata, *popt), 'g--')
plt.show()

我将其作为数组执行，因为我在很大范围内改变

temp

，这个单一元素“temp”是为了故障排除而创建的。

Xdata

和

Ydata

分别位于第一列和第二列： 数据

Answer 1

我不相信还有

curve_fit()

找不到更好的解决方案。

这是我尝试过的事情的列表。

基线

首先，我从您当前的解决方案开始，找到了该解决方案的 MSE（均方误差）。

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

def func(x,a,b,l):
    return (a*(x**2)) - (b*(x**3)) + ((l)*(x**4))

def mse(params, Xdata, Ydata):
    return np.mean((Ydata - func(Xdata, *params)) ** 2)

temp = np.linspace(80.,80.,1)
b = np.zeros_like(temp)
a = np.zeros_like(temp)
l = np.zeros_like(temp)

all_data = np.loadtxt('Pot_data_online.txt')
Xdata = all_data[:, 0]
Ydata = all_data[:, 1]

popt, pcov = curve_fit(func, Xdata, Ydata,bounds=((0.,0.,0.), (np.inf,np.inf,np.inf)))
b[0],a[0],l[0] = popt

plt.plot(Xdata,Ydata, label='data')
plt.plot(Xdata, func(Xdata, *popt), 'g--', label='fit')
print(popt)
plt.legend()
plt.show()

curve_fit_baseline = mse(popt, Xdata, Ydata)

使用不同的

method

我尝试将方法更改为

trf

、

dogbox

和

lm

。这些都没有产生明显更好的配合。

popt, pcov = curve_fit(
    func,
    Xdata,
    Ydata,
    bounds=((0.,0.,0.), (np.inf,np.inf,np.inf)),
    method='dogbox',
    ftol=1e-15,
    gtol=1e-15,
    xtol=1e-15,
)
fun = mse(popt, Xdata, Ydata)
print("curve_fit score", curve_fit_baseline)
print("dogbox score   ", fun)
print("dogbox improvement (higher is better)", (curve_fit_baseline - fun) / curve_fit_baseline)

输出：

curve_fit score 5.214343649304456e-08
dogbox score    5.2143436493044563e-08
dogbox improvement (higher is better) -1.269084921418044e-16

改进是负面的，所以这意味着切换到dogbox比

trf

的默认方法更糟糕。

lm

方法不能容忍边界，但如果删除它们，它也不会做得更好。

trf

从技术上讲可以比默认值好 1.26e-14%，但我不认为你会关心这么小的改进。

最小化器

有时

curve_fit()

可能会被太大或太小的梯度所混淆，而局部最小化器可以做得更好。

from scipy.optimize import minimize, basinhopping, Bounds

result = minimize(
    mse,
    np.zeros(3),
    bounds=Bounds((0.,0.,0.), (np.inf,np.inf,np.inf)),
    args=(Xdata, Ydata),
    options=dict(
        ftol=0,
        gtol=0,
    ),
    method='L-BFGS-B',
)

print("curve_fit score", curve_fit_baseline)
print("minimize score ", result.fun)
print("minimize improvement (higher is better)", (curve_fit_baseline - result.fun) / curve_fit_baseline)

这没有帮助。

我还尝试了 COBYQA 或 SLSQP 的不同方法。这也没有帮助。

全局优化

有时

curve_fit()

可能无法处理局部最小值。我尝试过盆地跳跃来避免这种情况。

示例：

from scipy.optimize import basinhopping, Bounds

result = basinhopping(
    mse,
    np.zeros(3),
    minimizer_kwargs=dict(
        bounds=Bounds((0.,0.,0.), (np.inf,np.inf,np.inf)),
        args=(Xdata, Ydata),
        options=dict(
            ftol=1e-16,
            gtol=1e-7,
        ),
    ),
    niter=1000,
)
print("curve_fit score   ", curve_fit_baseline)
print("basinhopping score", result.fun)
print("basinhopping improvement (higher is better)", (curve_fit_baseline - result.fun) / curve_fit_baseline)

这也没有帮助。

非负最小二乘法

我尝试通过预处理 Xdata 以包含 x 的所有相关幂，将这个问题重新格式化为线性回归。然后，SciPy 有一个专门的求解器来求解非负最小二乘问题。

from sklearn.preprocessing import PolynomialFeatures
from scipy.optimize import nnls

polynomial_transformer = PolynomialFeatures(degree=(2, 4), include_bias=False)
transformed = polynomial_transformer.fit_transform(Xdata.reshape(-1, 1))
# second term has minus sign. Negate feature to make this problem equivalent
transformed[:, 1] *= -1
x, rnorm = nnls(transformed, Ydata)
fun = mse(x, Xdata, Ydata)

print("curve_fit score", curve_fit_baseline)
print("lstsq score    ", fun)
print("lstsq improvement (higher is better)", (curve_fit_baseline - fun) / curve_fit_baseline)

这会将您当前的解决方案与最后一个小数联系起来。它的优点是比您当前的解决方案快 30 倍。

除了盆地跳跃结果之外，这还证明您的问题不包含局部最小值 - nnls 无法陷入局部最小值，因此如果问题可以重新格式化为 nnls 可以解决的问题，那么该问题不包含局部最小值最小值。

结论

您无法比

curve_fit()

为您提供的解决方案做得更好。您可以更改型号，但您已经表明这不是一个选项。