Как указать ошибку Y-переменной при подгонке с lmfit?

Question

Я почти новичок в Python и пытаюсь собрать данные из колледжа, используя lmfit. Переменная Y имеет переменную ошибку 3%. Как добавить эту ошибку в процесс подбора? Я изменяю от подгонки кривой scipy, и в scipy это было действительно легко сделать, просто создав массив со значениями ошибок и указав ошибку при подгонке, добавив текст "sigma = [yourarray]". Это мой текущий код:

from lmfit import Minimizer, Parameters, report_fit
import matplotlib.pyplot as plt
w1, V1, phi1, scal1 = np.loadtxt("./FiltroPasaBajo_1.txt", delimiter = "\t", unpack = True)

t = w1
eV= V1*0.03 + 0.01

def funcion(parametros, x, y):
    R = parametros['R'].value
    C = parametros['C'].value

    modelo = 4/((1+(x**2)*(R**2)*(C**2))**1/2)
    return modelo - y
parametros = Parameters()
parametros.add('R', value = 1000, min= 900, max = 1100)
parametros.add('C', value = 1E-6, min = 1E-7, max = 1E-5)

fit = Minimizer(funcion, parametros, fcn_args=(t,V1))
resultado = fit.minimize()


final = V1 + resultado.residual
report_fit(resultado)

try:
    plt.plot(t, V1, 'k+')
    plt.plot(t, final, 'r')
    plt.show()
except ImportError:
    pass

V1 - это значения, которые я измерил, а eV - массив ошибок. t является координатой х. Спасибо за ваше время

Answer 1

minimize() Функция минимизирует массив в смысле наименьших квадратов, настраивая параметры переменных для минимизации (resid**2).sum() для resid массив, возвращаемый вашей целевой функцией. Он действительно ничего не знает о неопределенности в ваших данных или даже о ваших данных. Чтобы использовать неопределенности в вашей форме, вам нужно передать в ваш массив eV так же, как вы проходите в t а также V1 и затем используйте это в своем вычислении массива, который будет минимизирован.

Обычно хочется свести к минимуму Sum[ (data-model)^2/epsilon^2 ], где epsilon неопределенность в данных (ваш eV), поэтому остаточный массив должен быть изменен с data-model в (data-model)/epsilon, Для вашей подгонки, вы хотели бы

def funcion(parametros, x, y, eps):
    R = parametros['R'].value
    C = parametros['C'].value

    modelo = 4/((1+(x**2)*(R**2)*(C**2))**1/2)
    return (modelo - y)/eps

а затем использовать это с

fit = Minimizer(funcion, parametros, fcn_args=(t, V1, eV))
resultado = fit.minimize()
...

Если вы используете lmfit.Model интерфейс (разработан для подгонки кривой), то вы можете перейти в weights массив, который умножается data -modelи так будет 1.0 / eV представлять вес для неопределенности (как указано выше с minimize). С использованием lmfit.Model Интерфейс и обеспечение неопределенности будут выглядеть следующим образом:

from lmfit import Model
# model function, to model the data
def func(t, r, c):
    return 4/((1+(t**2)*(r**2)*(c**2))**1/2)

model  = Model(func)
parametros = model.make_params(r=1000, c=1.e-6)

parametros['r'].set(min=900, max=1100)
parametros['c'].set(min=1.e-7, max=1.e-5)

resultado = model.fit(V1, parametros, t=t, weights=1.0/eV)

print(resultado.fit_report())

plt.errorbar(t, V1, eV, 'k+', label='data')
plt.plot(t, resultado.best_fit, 'r', label='fit')
plt.legend()
plt.show()

надеюсь, это поможет....

Answer 2

Я думаю, что вы не можете предоставить сигму в fit.minimize() напрямую.

Однако я вижу, что fit.minimize() использует метод leastsq scipy (по умолчанию), который является тем же методом, что и scipy's curve_fit.

Если вы посмотрите на источник Scipy's Curve_fit, он будет следовать сигма (для 1-го случая).

 transform = 1.0 / sigma
 jac = _wrap_jac(jac, xdata, transform)
 res = leastsq(func, p0, Dfun=jac, full_output=1, **kwargs)

Поскольку fit.minimize() позволяет вам передавать kwargs (Dfun) для leastsq, вы можете передать jac так, как это делается в scipy curve_fit.