Кусочно-линейная посадка

Question

Я пытаюсь автоматизировать подгонку спектральных плотностей мощности. Типичный график выглядит следующим образом (масштаб журнала):

в красном - одна из моих попыток подобрать его. Я хотел бы иметь возможность разместить до 3 линейных сегментов на графиках плотности мощности. Я уже попробовал несколько подходов, например, глядя на 1 с использованием curve_fit:

def logfit(x, a1, a2, b, cutoff):
    cutoff = int(params[3])
    out = np.empty_like(x)
    out[:cutoff] = x[:cutoff]*a1 + b
    out[cutoff:] = x[cutoff]*a1 + b + (x[cutoff:] - x[cutoff])*a2 
    return out   
#the fit is only linear on loglog scale, need to apply np.log10 to x and y.
popt, pcov = curve_fit(logfit, np.log10(x), np.log10(y), bounds = ([-2,-2,-np.inf,99],[0,0,np.inf,100]))

, который обычно обеспечивает плоскую посадку. Я также посмотрел на [2], используя пакет lmfit:

def logfit(params, x, data):
    a1, a2, b = params['a1'],  params['a2'], params['b']
    c = int(params['cutoff'])
    out = np.empty_like(x)
    out[:c] = x[:c]*a1 + b
    out[c:] = x[c]*a1 + b + (x[c:] - x[c])*a2 
    return data - out    

    params = Parameters()
    params.add('a1', value = -0.3, min = -2, max = 0)
    params.add('a2', value = -2, min = -2, max = -0.5)
    params.add('b', min = 0, max = 5)
    params.add('cutoff', min = 150, max = 500)
    params = minimize(logfit, params, args=(f, pxx))

Но это также приводит к полностью отключенным подгонкам (как красная подгонка на графике выше). Это потому что слишком много параметров? Я бы подумал, что эту проблему достаточно просто решить с помощью оптимизации ...

1 Как применить кусочно-линейное соответствие для линии с положительным и отрицательным наклоном в Python?

[2] Подогнать кривую для данных, составленных из двух различных режимов

Answer 1

Всегда полезно включать полный пример, показывающий проблему.

Конечно, отчасти проблема в том, как вы сформулировали проблему, заключается в том, что переменная cutoff не будет оптимизирована в соответствии. Все решатели, доступные с scipy.optimize или lmfit, работают строго с непрерывными переменными, а не с целыми или дискретными переменными. Решатели работают, внося небольшие изменения (например, на уровне 1.e-7) в переменные и наблюдая, как это меняет соответствие. Если вы берете c = int(params['cutoff']), значение c не изменится вообще, и подгонка решит, что небольшие изменения в cutoff не изменят подгонку.

Один из подходов к преодолению этого заключается в том, что обычно работает довольно хорошо, это не использовать дискретную отсечку, а использовать пошаговую функцию - функция logisti c, вероятно, является наиболее распространенным выбором - которая занимает шаг от 0 до 1 в течение небольшого интервала, скажем, 1 или 2 x баллов. С такой функцией, которая непрерывно распространяется на несколько x точек данных, небольшое изменение в cutoff внесет небольшое изменение в расчетную модель, так что подбор сможет найти хорошее решение.

Вот пример этого. FWIW, lmfit включает в себя функцию logistic, но это всего лишь одна строка, поэтому я включу ее использование и выполнение трудным путем:

import numpy as np
import matplotlib.pyplot as plt

from lmfit import Model
from lmfit.lineshapes import logistic

# build some fake data comprised of two line segments
x = np.linspace(0, 50, 101)
y = 37.0  - 0.62*x
y[62:] = y[62] - 5.41*(x[62:] - x[62])
y = y + np.random.normal(size=len(y), scale=0.5)

def two_lines(x, offset1, slope1, offset2, slope2, cutoff, width=1.0):
    """two line segments, joined at a cutoff point by a logistic step"""
    # step = 1 - 1./(1. + np.exp((x - cutoff)/max(1.-5, width)))
    step = logistic(x, center=cutoff, sigma=width)
    return (offset1 + slope1*x)*(1-step)  + (offset2 + slope2*x)*(step)

lmodel = Model(two_lines)
params = lmodel.make_params(offset1=1.8, slope1=0, offset2=10, slope2=-4,
                            cutoff=25.0, width=1.0)
# the width could be varied, or just set to ~1 `x` interval or so.
# here, we fix it at 'half a step'
params['width'].value = (x[1] - x[0]) /2.0
params['width'].vary = False

# do the fit, print fit report   
result = lmodel.fit(y, params, x=x)  
print(result.fit_report())

# plot the results
plt.plot(x, y, label='data')
plt.plot(x, result.best_fit, label='fit')
plt.legend()
plt.show()

Это распечатает отчет как

[[Model]]
    Model(two_lines)
[[Fit Statistics]]
    # fitting method   = leastsq
    # function evals   = 115
    # data points      = 101
    # variables        = 5
    chi-square         = 27.2552501
    reduced chi-square = 0.28390885
    Akaike info crit   = -122.297309
    Bayesian info crit = -109.221707
[[Variables]]
    offset1:  36.8909381 +/- 0.13294218 (0.36%) (init = 1.8)
    slope1:  -0.62144415 +/- 0.00742987 (1.20%) (init = 0)
    offset2:  185.617291 +/- 0.71277881 (0.38%) (init = 10)
    slope2:  -5.41427487 +/- 0.01713805 (0.32%) (init = -4)
    cutoff:   31.3853560 +/- 0.22330118 (0.71%) (init = 25)
    width:    0.25 (fixed)
[[Correlations]] (unreported correlations are < 0.100)
    C(offset2, slope2) = -0.992
    C(offset1, slope1) = -0.862
    C(offset2, cutoff) = -0.338
    C(slope2, cutoff)  =  0.318

и сюжет вроде

Надеюсь, этого достаточно, чтобы вы пошли.

Answer 2

Благодаря @M Newville и используя ответы от [1] я придумал следующее рабочее решение:

def continuousFit(x, y, weight = True):
    """
        fits data using two lines, forcing continuity between the fits.
        if `weights` = true, applies a weights to datapoints, to cope with unevenly distributed data.
    """
    lmodel = Model(continuousLines)
    params = Parameters()
    #Typical values and boundaries for power-spectral densities :
    params.add('a1', value = -1, min = -2, max = 0)
    params.add('a2', value = -1, min = -2, max = -0.5)
    params.add('b', min = 1, max = 10)
    params.add('cutoff', min = x[10], max = x[-10])
    if weight:
        #weights are based on the density of datapoints on the x-axis
        #since there are fewer points on the low-frequencies they get heavier weights.
        weights = np.abs(x[1:] - x[:-1])
        weights = weights/np.sum(weights)
        weights = np.append(weights, weights[-1])
        result = lmodel.fit(y, params, weights = weights, x=x)  
        return result.best_fit
    else:
        result = lmodel.fit(y, params, x=x)  
        return result.best_fit

def continuousLines(x, a1, a2, b, cutoff):
    """two line segments, joined at a cutoff point in a continuous manner"""
    return a1*x + b  + a2*np.maximum(0, x - cutoff)

[1]: Подогнать кривую для данных, составленных из двух отдельные режимы