协方差矩阵是 lmfit Python3-6 中的 NoneType

Question

我想使用 lmfit 拟合 f(x,y) 函数。数据集小，拟合参数多（x轴6个点，y轴11个点，16个无约束拟合参数）。使用 Model.fit 中的所有默认值，我无法获得协方差矩阵，并且在拟合过程中，自由参数的值根本没有改变。

我试图更改参数的初始值。然而，当我在 OriginPro 曲面拟合功能中设置相同类型的问题时，Levenberg-Marquardt 算法设法拟合数据并估计误差（尽管某些参数的值相当大）。这意味着我的代码一定有问题。我找不到问题所在。我不是Python大师

MWE如下。

import numpy as np
from lmfit import Model, Parameters
import numdifftools # not calling this doesn't change anything

x, y = np.array([226.5, 361.05, 404.41, 589, 632.8, 1013.98]), np.linspace(0,100,11)
X, Y = np.meshgrid(x, y)

Z = np.array([[1.3945, 1.34896, 1.34415, 1.33432, 1.33306, 1.32612],\
[1.39422, 1.3487, 1.34389, 1.33408, 1.33282, 1.32591],\
[1.39336, 1.34795, 1.34315, 1.33336, 1.33211, 1.32524],\
[1.39208, 1.34682, 1.34205, 1.3323, 1.33105, 1.32424],\
[1.39046, 1.3454, 1.34065, 1.33095, 1.32972, 1.32296],\
[1.38854, 1.34373, 1.33901, 1.32937, 1.32814, 1.32145],\
[1.38636, 1.34184, 1.33714, 1.32757, 1.32636, 1.31974],\
[1.38395, 1.33974, 1.33508, 1.32559, 1.32438, 1.31784],\
[1.38132, 1.33746, 1.33284, 1.32342, 1.32223, 1.31576],\
[1.37849, 1.33501, 1.33042, 1.32109, 1.31991, 1.31353],\
[1.37547, 1.33239, 1.32784, 1.31861, 1.31744, 1.31114]])

#This has to be defined beforehand (otherwise parameters names are not defined error)
a1,a2,a3,a4 = 1.3208, -1.2325E-5, -1.8674E-6, 5.0233E-9
b1,b2,b3,b4 = 5208.2413, -0.5179, -2.284E-2, 6.9608E-5
c1,c2,c3,c4 = -2.5551E8, -18341.336, -920, 2.7729
d1,d2,d3,d4 = 9.3495, 2E-3, 3.6733E-5, -1.2932E-7

# Function to fit
def model(x, y, *args):
    return a1+a2*y+a3*np.power(y,2)+a4*np.power(y,3)+\
    (b1+b2*y+b3*np.power(y,2)+b4*np.power(y,3))/np.power(x,2)+\
    (c1+c2*y+c3*np.power(y,2)+c4*np.power(y,3))/np.power(x,4)+\
    (d1+d2*y+d3*np.power(y,2)+d4*np.power(y,3))/np.power(x,6)

# This is the callable that is passed to Model.fit. M is a (2,N) array
# where N is the total number of data points in Z, which will be ravelled
# to one dimension.
def _model(M, **args):
    x, y = M
    arr = model(x, y, params)
    return arr

# We need to ravel the meshgrids of X, Y points to a pair of 1-D arrays.
xdata = np.vstack((X.ravel(), Y.ravel()))

# Fitting parameters.
fmodel = Model(_model)
params = Parameters()
params.add_many(('a1',1.3208,True,1,np.inf,None,None),\
                ('a2',-1.2325E-5,True,-np.inf,np.inf,None,None),\
                ('a3',-1.8674E-6,True,-np.inf,np.inf,None,None),\
                ('a4',5.0233E-9,True,-np.inf,np.inf,None,None),\
                ('b1',5208.2413,True,-np.inf,np.inf,None,None),\
                ('b2',-0.5179,True,-np.inf,np.inf,None,None),\
                ('b3',-2.284E-2,True,-np.inf,np.inf,None,None),\
                ('b4',6.9608E-5,True,-np.inf,np.inf,None,None),\
                ('c1',-2.5551E8,True,-np.inf,np.inf,None,None),\
                ('c2',-18341.336,True,-np.inf,np.inf,None,None),\
                ('c3',-920,True,-np.inf,np.inf,None,None),\
                ('c4',2.7729,True,-np.inf,np.inf,None,None),\
                ('d1',9.3495,True,-np.inf,np.inf,None,None),\
                ('d2',2E-3,True,-np.inf,np.inf,None,None),\
                ('d3',3.6733E-5,True,-np.inf,np.inf,None,None),\
                ('d4',-1.2932E-7,True,-np.inf,np.inf,None,None))

result = fmodel.fit(Z.ravel(), params, M=xdata)
fit = model(X, Y, result.params)

print(result.covar)

此代码导致协方差为 NoneType。我希望它最终会被计算出来，因为 Origin 可以以某种方式管理。如果需要，我可以提供 Origin Surface Fitting Parameters 中的所有参数。

绘制 Z 拟合差异时，低 x 值存在相当大的差异（在 Origin 中没有发生）。

Answer 1

您没有以 lmfit 可以合理使用的方式定义您的模型函数。你有：

def _model(M, **args):
    x, y = M
    arr = model(x, y, params)
    return arr


def model(x, y, *args):
    return a1+a2*y+a3*np.power(y,2)+a4*np.power(y,3)+\
    (b1+b2*y+b3*np.power(y,2)+b4*np.power(y,3))/np.power(x,2)+\
    (c1+c2*y+c3*np.power(y,2)+c4*np.power(y,3))/np.power(x,4)+\
    (d1+d2*y+d3*np.power(y,2)+d4*np.power(y,3))/np.power(x,6)


model = Model(_model)

有几个问题：

args 未在 _model 中使用，并且 params 未在函数中定义，因此将是模块级的。
与 model 类似，args 未使用，a1、a2 等将从模块级（编程）变量和（重要的是!!) 这些将不会更新。

简而言之，您的模型函数永远不会看到参数的不同值。

lmfit.Model 采用 命名函数参数 并将其转换为参数名称。它不会将 **kws 或 *position_args 转换为参数名称。所以我认为你想要做的是编写一个这样的模型函数：

def model(x, y, a1, a2, a2, a4, b1, b2, b3 ,b3, c1, c2, c3, c4,
          d1, d2, d3, d4):
    return a1+a2*y+a3*np.power(y,2)+a4*np.power(y,3)+\
    (b1+b2*y+b3*np.power(y,2)+b4*np.power(y,3))/np.power(x,2)+\
    (c1+c2*y+c3*np.power(y,2)+c4*np.power(y,3))/np.power(x,4)+\
    (d1+d2*y+d3*np.power(y,2)+d4*np.power(y,3))/np.power(x,6)

然后创建一个模型：

# Note: don't give a function and Model instance the same name!!
my_model = Model(model, independent_vars=('x', 'y'))

定义该模型后，您可以运行进行拟合，而无需解开数据（lmfit 中的独立数据几乎可以是任何数据类型，数据数组可以是多个维):

result = my_model.fit(Z, params, x=X, y=Y)

就其价值而言，从适合运行完成的意义上来说，进行此类更改对我有用。拟合仍然停留在某些参数未从其初始值更新的情况下，但这与设置和运行拟合的机制是一个独立的问题，可能是由于多项式非常不稳定或初步估计不佳。

顺便说一句：np.power(y,n) 可以拼写为 y**n 并且可读性很重要。此外，有时会通过替换

来提高数值稳定性

a + b*x + c*x**2 + d*x**3

和

a + x*(b + x*(c + x*d))

虽然我不知道这对你的情况是否有帮助。

协方差矩阵是 lmfit Python3-6 中的 NoneType

Covariance matrix is NoneType in lmfit Python3-6

python-3.x

lmfit

scipy-optimize