Подгонка трапеции с кривой

Question

Я пытаюсь вписать трапецию в набор временных рядов, используя библиотеку curve_fit из scipy.optimize. Функция, которую я использую для генерации трапеции, следующая:

def trapezoid(x, a, b, c, tau1, tau2):
    y = np.zeros(len(x))
    c = -np.abs(c)
    a = np.abs(a)
    y[:int(tau1)] =  a*x[:int(tau1)] + b
    y[int(tau1):int(tau2)] =  a*tau1 + b 
    y[int(tau2):] = c*(x[int(tau2):]-tau2) + (a*tau1 + b)
    return y

Где a и c - наклоны, а tau1 и tau2 отмечают начало и конец плоской фазы.

А чтобы соответствовать, я просто использую:

popt, pcov = curve_fit(trapezoid, xdata, ydata, method = 'lm')

В большинстве случаев он работает просто отлично, например, в следующем:

Тем не менее, я также получаю некоторые случаи, когда он просто не вписывается в данные, где, похоже, все должно быть в порядке:

Проблема в этих случаях заключается в том, что он устанавливает tau2 (конец плоской фазы) меньше, чем tau1 (начало).

Может ли кто-нибудь предложить способ решения этой проблемы? Будь то путем наложения ограничения или каким-либо другим способом?

Пример массива, для которого не работает подгонка:

массив ([1,2, 1,21, 1,2, 1,19, 1,21, 1,22, 2,47, 2,53, 2,49, 2,39, 2,28, 2,16, 2,07, 1,99, 1,91, 1,83, 1,74, 1,65, 1,57, 1,5, 1,45, 1,41, 1,38, 1,35, 1,33, 1,29, 1,24, 1,19, 1,14, 1,11, 1,07, 1,04, 1, 0,95, 0,91, 0,87, 0,84, 0,8, 0,77, 0,74, 0,72, 0,7, 0,68, 0,66, 0,63, 0,61, 0,59, 0,57, 0,55, 0,52, 0,5, 0,48, 0,45, 0,43, 0,41, 0,39, 0,38, 0,37, 0,37, 0,36, 0,35, 0,34, 0,34, 0,33])

Что дает: тау1: 8,45, тау2: 5,99

Answer 1

Работает только установка минимального и максимального значений t1, t2

def trapezoid(x, a, b, c, tau1, tau2):
    y = np.zeros(len(x))
    c = -np.abs(c)
    a = np.abs(a)
    (tau1,tau2) = (min(tau1,tau2),max(tau1,tau2))

    y[:int(tau1)] =  a*x[:int(tau1)] + b
    y[int(tau1):int(tau2)] =  a*tau1 + b 
    y[int(tau2):] = c*(x[int(tau2):]-tau2) + (a*tau1 + b)



x_data = np.arange(len(A))
popt, pcov = curve_fit(trapezoid, x_data, A, method = 'lm')
print popt
fit = trapezoid(x_data,*popt)

ведет к:

Answer 2

Для этой проблемы может оказаться полезным lmfit (http://lmfit.github.io/lmfit-py/)). Lmfit предоставляет немного более высокий уровень интерфейса для подбора кривой, все еще основанный на оптимизаторах scipy, но с некоторыми лучшими абстракциями и функциями.

В частности, для вашего вопроса lmfit параметры - это объекты Python, которые могут иметь границы, быть фиксированными или записываться как простые алгебраические ограничения в терминах других переменных. Это может поддержать наложение tau2 > tau1. Идея состоит в том, чтобы установить tau2=tau1+taudiff и установить нижнюю границу 0 на taudiff. Хотя вы можете переписать свою функцию, чтобы сделать это в коде, с lmfit вам не нужно это делать, и вы можете вместо этого поместить эту логику в Параметры.

Преобразование вашего скрипта для использования lmfit даст что-то вроде этого:

from lmfit import Model

# use your same model function
def trapezoid(x, a, b, c, tau1, tau2):
    y = np.zeros(len(x))
    c = -np.abs(c)
    a = np.abs(a)
    y[:int(tau1)] =  a*x[:int(tau1)] + b
    y[int(tau1):int(tau2)] =  a*tau1 + b 
    y[int(tau2):] = c*(x[int(tau2):]-tau2) + (a*tau1 + b)
    return y

# turn model function into lmfit Model
tmod = Model(trapezoid)

# create Parameters for this model: they will be *named* according
# to the signature of the model function, and be used as keys in
# an ordered-directory-derived object.  Here you can also give
# initial values
params = tmod.make_params(a=1, b=2, c=0.5, tau1=5, tau2=-1)

# now you can set bounds or constraints.  

# 1st, add a new variable "taudiff"
params.add('taudiff', value=0.1, min=0, vary=True)

# constraint tau2 to be taudiff+tau1 -- this is no longer a "free variable:
params['tau2'].expr = "taudiff + tau1"

# now do fit to data:
result = tmod.fit(ydata, params, x=xdata)

# print report of fit
print(result.fit_report())

# get best fit params:
for parname, param in result.params:
    print(parname, param.value, param.stderr, param.expr)

# get best fit array for plotting
pylab.plot(xdata, ydata)
pylab.plot(xdata, result.best_fit)

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