分段线性拟合

Question

我正在尝试自动拟合功率谱密度。一个典型的情节看起来像这样（loglog scale）：

红色是我尝试安装它的尝试之一。 我想有一种方法可以将最多 3 个线性段拟合到功率谱密度图。 我已经尝试了几种方法，例如使用 curve_fit 查看1：

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))

但这也会产生完全关闭的拟合（如上图中的红色拟合）。 这是因为参数太多了吗？我原以为这个问题很简单，可以通过优化来解决……

1How to apply piecewise linear fit for a line with both positive and negative slopes in Python?

[2]Fit a curve for data made up of two distinct regimes

Answer 1

包含一个显示问题的完整示例总是有帮助的。

当然，您提出问题的方式的部分问题是cutoff 变量不会被优化。 scipy.optimize 或 lmfit 提供的所有求解器都严格处理连续变量，而不是整数或离散变量。求解器通过对变量进行小的更改（如在 1.e-7 级别）并查看它如何改变拟合来工作。在您使用c = int(params['cutoff']) 时，c 的值根本不会改变，并且拟合将决定对cutoff 的微小更改不会改变拟合。

克服这个问题的一种方法是，通常效果很好的方法是不使用离散截止值，而是使用类似阶梯的函数 - 逻辑函数可能是最常见的选择 - 从 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]：Fit a curve for data made up of two distinct regimes