Ковариационная матрица имеет тип NoneType в lmfit Python3-6 - PullRequest
Новая
       35

Ковариационная матрица имеет тип NoneType в lmfit Python3-6

0 голосов
/ 09 октября 2019

Я хочу подогнать функцию af (x, y) с помощью lmfit. Набор данных небольшой, и имеется множество параметров подгонки (6 точек по оси x, 11 точек по оси y и 16 параметров подбора без ограничений). Используя все значения по умолчанию из Model.fit, я не могу получить ковариационную матрицу, и в процессе подбора значения свободных параметров вообще не меняются.

Я пытался изменить начальные значения параметров. Тем не менее, когда я устанавливаю такую ​​же проблему в функциях Подгонки поверхности OriginPro, алгоритм Левенберга-Марквардта удачно подбирает данные и оценивает ошибки (хотя и имеет довольно большие значения для определенных параметров). Это означает, что должна быть некоторая проблема с моим кодом. Я не могу найти, где проблема. Я не владелец 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 может каким-то образом управлять. Если это необходимо, я могу предоставить все параметры из Параметры подгонки поверхности начала координат.

При построении разницы Z-подгонки существует большое расхождение для низких значений x (не происходит в исходной точке).

1 Ответ

0 голосов
/ 12 октября 2019

Вы не определяете функцию своей модели таким образом, чтобы ее можно было разумно использовать 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)

С некоторыми проблемами:

  1. args не используется в _model, а params не определен в функции, поэтому будетбыть на уровне модуля.
  2. Аналогично в 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))

Хотя я не знаю, поможет ли это в вашем случае.

...