如何将分段阶跃函数拟合到我的数据点下方？

Question

我有一个带有 6 个参数的自定义函数。它创建了一条阶梯“曲线”。
我想优化/拟合最后 5 个参数（m1、m2、FL1、FL2 和 FL3）。

def steps(x, FL1, m1, FL2, m2, FL3):
    if x > m1:
        return FL1
    if x > m2:
        return FL2
    else:
        return FL3

我希望这个函数始终位于低于我的数据点。

示例：

真实数据看起来稍微复杂一点，像这样，但仍然不复杂：

我可以使用什么功能来适应这个功能？

我认为

scipy.optimize.curve_fit

不起作用，因为使用了

least_squares

。这样它将以数据点为中心。

from scipy.optimize import curve_fit  

def steps(m, FL1, m1, FL2, m2, FL3):
    if m > m1:
        return FL1
    if mass > m2:
        return FL2
    else:
        return FL3

x_data = [8000, 7000, 6000, 5000]  
y_data = [300, 350, 400, 450]  

popt, _ = curve_fit(steps, x_data, y_data)  

这失败了：

ValueError：具有多个元素的数组的真值不明确。使用 a.any() 或 a.all()

它似乎不喜欢函数中的

if

语句。

Answer 1

你的优化只有2个参数需要优化；

m1

，

m2

。台阶的高度根据您的插值数据得出。 例如，总是

FL1 == y_data[0]

。

您需要在函数中表达这一点，并且可以在

curve_fit

中使用它。它还希望你对你的函数进行向量化，但我只是不介意。

如果我们表达这个函数

import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
import numpy as np

x_data = np.linspace(8000, 5000, 10)
y_data = np.array([400, 420, 440, 450, 453, 454, 452, 453, 453, 452])


def limit(m):
    return np.interp(m, np.flip(x_data), np.flip(y_data))


def constrained_steps(m, m1, m2):
    results = list()  # Could be vectorized if you cared enough.
    for i in range(len(m)):
        if m[i] > m1:
            results.append(y_data[0])
        elif m[i] > m2:
            results.append(limit(m1))
        else:
            results.append(limit(m2))
    return np.array(results)


# Initialize initial guess at 2/3 and 1/3's through the x_data.
bounds = x_data[[-1,0]]
p0 = ((bounds[1]-bounds[0])*2/3 + bounds[0], (bounds[1]-bounds[0])*1/3 + bounds[0])
popt, _ = curve_fit(constrained_steps, x_data, y_data, p0=p0, bounds=bounds)
opt_m1, opt_m2 = popt

ms = np.linspace(bounds[0], bounds[1], 1000)
plt.scatter(x_data, y_data)
plt.plot(ms, constrained_steps(ms, opt_m1, opt_m2))
plt.gca().invert_xaxis()
plt.show()

Answer 2

我同意Mikael的回答，你应该尽量让你的表达尽可能简洁。我有一个真实数据形状的示例，它结合了线性函数和常数函数。

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

# Define step function
def step_func(x, a, c, f1, f2):
    result = []
    for xi in x:
        if xi >= 0 and xi < f1:
            result.append(a * xi)
        elif xi >= f1 and xi < f2:
            result.append(a * f1)
        elif xi >= f2:
            result.append(c * xi + a * f1 - f2 * c)
        else:
            result.append(0)
    return np.asarray(result)

# Generating some demo data
x_data = np.linspace(0, 10, 100)
y_data = step_func(x_data, 2, 1, 3, 7) + np.random.normal(0, 0.5, x_data.size)  # 添加噪声

# Fitting function
def fitting_func(x, a, c, f1, f2):
    return step_func(x, a, c, f1, f2)

# Fitting for data
popt, pcov = curve_fit(fitting_func, x_data, y_data, p0=[1, -1, 4, 6])

# Obtain fitting parameters
a_fit, c_fit, f1_fit, f2_fit = popt
print(f'Fitting param: a={a_fit}, c={c_fit}, f1={f1_fit}, f2={f2_fit}')

# Generat fitting results
y_fit = fitting_func(x_data, *popt)

# Plot
plt.scatter(x_data, y_data, label='Data', s=10)
plt.plot(x_data, y_fit, color='red', label='Fitted')
plt.vlines([f1_fit, f2_fit], ymin = min(x_data), ymax = max(x_data) )
# plt.xlim(0,16)
# plt.ylim(0,1)
plt.legend()
plt.xlabel('x')
plt.ylabel('y')
plt.title('Step func fitting')
plt.show()

以上返回：

Fitting param: a=1.9532567200579396, c=1.085588494015425, f1=3.084102835100397, f2=7.1628504872999095