使用 scipy 优化 curve_fit 来拟合步长位置变化的步长函数

Question

我正在尝试拟合看起来像

的 x y 数据

x = np.linspace(-2, 2, 1000)
a = 0.5
yl = np.ones_like(x[x < a]) * -0.4 + np.random.normal(0, 0.05, x[x < a].shape[0])
yr = np.ones_like(x[x >= a]) * 0.4 + np.random.normal(0, 0.05, x[x >= a].shape[0])
y = np.concatenate((yl, yr))
plt.scatter(x, y, s=2, color='k')

我正在使用 Heaviside 阶跃函数的变体

def f(x, a, b): return 0.5 * b * (np.sign(x - a))

并适合

popt, pcov = curve_fit(f, x, y, p0=p)

其中 p 是一些初始猜测。对于任何 p curve_fit 只适合 b 而不是 a 例如：

popt, pcov = curve_fit(f, x, y, p0=[-1.0, 0]) 我们得到 popt 是 [-1., 0.20117665]

popt, pcov = curve_fit(f, x, y, p0=[.5, 2]) 我们得到 taht popt 是 [.5, 0.79902]

popt, pcov = curve_fit(f, x, y, p0=[1.5, -2]) 我们得到 taht popt 是 [1.5, 0.40128229]

为什么curve_fit不适合a？

Answer 1

这是使用您的数据和函数的图形 Python 拟合器，scipy 的 differential_evolution 遗传算法模块用于为 curve_fit 提供初始参数估计.该模块使用拉丁超立方体算法来确保彻底搜索参数 space，这需要搜索范围。在此示例中，这些边界取自数据最大值和最小值。

import numpy, scipy, matplotlib
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from scipy.optimize import differential_evolution
import warnings


# generate data for testing
x = numpy.linspace(-2, 2, 1000)
a = 0.5
yl = numpy.ones_like(x[x < a]) * -0.4 + numpy.random.normal(0, 0.05, x[x < a].shape[0])
yr = numpy.ones_like(x[x >= a]) * 0.4 + numpy.random.normal(0, 0.05, x[x >= a].shape[0])
y = numpy.concatenate((yl, yr))

# alias data to match pervious example
xData = x
yData = y


def func(x, a, b): # variation of the Heaviside step function
    return 0.5 * b * (numpy.sign(x - a))


# function for genetic algorithm to minimize (sum of squared error)
def sumOfSquaredError(parameterTuple):
    warnings.filterwarnings("ignore") # do not print warnings by genetic algorithm
    val = func(xData, *parameterTuple)
    return numpy.sum((yData - val) ** 2.0)


def generate_Initial_Parameters():
    # min and max used for bounds
    maxX = max(xData)
    minX = min(xData)

    parameterBounds = []
    parameterBounds.append([minX, maxX]) # search bounds for a
    parameterBounds.append([minX, maxX]) # search bounds for b

    # "seed" the numpy random number generator for repeatable results
    result = differential_evolution(sumOfSquaredError, parameterBounds, seed=3)
    return result.x

# by default, differential_evolution completes by calling curve_fit() using parameter bounds
geneticParameters = generate_Initial_Parameters()

# now call curve_fit without passing bounds from the genetic algorithm,
# just in case the best fit parameters are aoutside those bounds
fittedParameters, pcov = curve_fit(func, xData, yData, geneticParameters)
print('Fitted parameters:', fittedParameters)
print()

modelPredictions = func(xData, *fittedParameters) 

absError = modelPredictions - yData

SE = numpy.square(absError) # squared errors
MSE = numpy.mean(SE) # mean squared errors
RMSE = numpy.sqrt(MSE) # Root Mean Squared Error, RMSE
Rsquared = 1.0 - (numpy.var(absError) / numpy.var(yData))

print()
print('RMSE:', RMSE)
print('R-squared:', Rsquared)

print()


##########################################################
# graphics output section
def ModelAndScatterPlot(graphWidth, graphHeight):
    f = plt.figure(figsize=(graphWidth/100.0, graphHeight/100.0), dpi=100)
    axes = f.add_subplot(111)

    # first the raw data as a scatter plot
    axes.plot(xData, yData,  'D')

    # create data for the fitted equation plot
    xModel = numpy.linspace(min(xData), max(xData))
    yModel = func(xModel, *fittedParameters)

    # now the model as a line plot
    axes.plot(xModel, yModel)

    axes.set_xlabel('X Data') # X axis data label
    axes.set_ylabel('Y Data') # Y axis data label

    plt.show()
    plt.close('all') # clean up after using pyplot

graphWidth = 800
graphHeight = 600
ModelAndScatterPlot(graphWidth, graphHeight)

Answer 2

正如其他人所提到的，curve_fit（以及 scipy.optimize 中的所有其他求解器）可以很好地优化连续而非离散变量。它们都通过对参数值进行小的更改（例如，在 1.e-7 级别）并查看结果中发生了什么（如果有的话）更改，并使用该更改来优化这些值直到最小发现有残留。使用 np.sign:

的模型函数

def f(x, a, b): return 0.5 * b * (np.sign(x - a))

a 值的如此小的变化根本不会改变模型或拟合结果。也就是说，首先拟合将尝试 a=-1.0 或 a=0.5 的起始值，然后将尝试 a=-0.999999995 或 a=0.500000005。对于 np.sign(x-a)，它们都会给出相同的结果。拟合不知道它需要将 a 更改为 1 才能对结果产生任何影响。它不可能知道这一点。 np.sign() 和 np.sin() 相差一个字母，但在这方面的表现却截然不同。

真实数据采取一步但要足够精细地采样以便该步骤不会在一步中完全发生是很常见的。在这种情况下，您将能够使用各种函数形式（线性斜坡、误差函数、反正切、逻辑等）对步骤进行建模。 @JamesPhilipps 的彻底回答给出了一种方法。我可能会使用 lmfit（作为其主要作者之一）并愿意通过查看数据来猜测参数的起始值，也许：

import numpy as np

x = np.linspace(-2, 2, 1000)
a = 0.5
yl = np.ones_like(x[x < a]) * -0.4 + np.random.normal(0, 0.05, x[x < a].shape[0])
yr = np.ones_like(x[x >= a]) * 0.4 + np.random.normal(0, 0.05, x[x >= a].shape[0])
y = np.concatenate((yl, yr))

from lmfit.models import StepModel, ConstantModel

model = StepModel() + ConstantModel()
params = model.make_params(center=0, sigma=1, amplitude=1., c=-0.5)

result = model.fit(y, params, x=x)

print(result.fit_report())

import matplotlib.pyplot as plt
plt.scatter(x, y, label='data')
plt.plot(x, result.best_fit, marker='o', color='r', label='fit')
plt.show()

这将给出一个很好的拟合并打印出

的结果

[[Model]]
    (Model(step, form='linear') + Model(constant))
[[Fit Statistics]]
    # fitting method   = leastsq
    # function evals   = 50
    # data points      = 1000
    # variables        = 4
    chi-square         = 2.32729556
    reduced chi-square = 0.00233664
    Akaike info crit   = -6055.04839
    Bayesian info crit = -6035.41737
##  Warning: uncertainties could not be estimated:
[[Variables]]
    amplitude:  0.80013762 (init = 1)
    center:     0.50083312 (init = 0)
    sigma:      4.6009e-04 (init = 1)
    c:         -0.40006255 (init = -0.5)

请注意，它会找到步长的 center，因为它假设步长有一些有限宽度 (sigma)，但后来发现宽度小于步长在 x。但也请注意，它无法计算参数中的不确定性，因为如上所述，解附近 center（您的 a）的微小变化不会改变结果拟合。 FWIW StepModel 可以使用线性函数、误差函数、反正切函数或逻辑函数作为阶跃函数。

如果您构建的测试数据具有较小的阶梯宽度，例如像

from scipy.special import erf    
y = 0.638  * erf((x-0.574)/0.005)  + np.random.normal(0, 0.05, len(x))

那么拟合就能够找到最佳解决方案并评估不确定性。

我希望这能解释为什么与您的模型函数的拟合无法细化 a 的值，以及对此可以做些什么。

使用 scipy 优化 curve_fit 来拟合步长位置变化的步长函数

fitting step function with variation in the step location with scipy optimize curve_fit

python

numpy

curve-fitting

scipy-optimize