You need to enable JavaScript to run this app.
优惠活动
大模型
产品
解决方案
定价
更多

带权重的分段线性回归序列分割方案咨询(已知分段数n)

带权重的序列分段线性拟合与分割实现方案

我需要将一个序列分割为已知数量n的子序列,要求每个子序列内的点能用分段线性函数拟合(最小化每个子序列及整体的拟合距离)。目前已尝试用ruptures库的Binseg算法(支持指定分段数),也试过直接用numpy/scipy拟合分段线性函数,但发现必须给数据点添加权重才能达到理想效果,请问如何改造现有方案,或者有没有可直接接收权重数组作为参数的解决方案?

补充上下文

  • 曲线形态:通常为平坦、陡升、缓降、凹形、凸形或混合形态,最多包含40个点
  • 异常值:可能存在1-2个异常值(通常为第一个点,但不绝对)
  • 分段数:目标分段数n介于4到8之间(作为参数传入)
  • 性能要求:单序列最多40个点,需循环执行约20次,总耗时需小于30秒,即单序列处理耗时最多约1.5秒

示例数据

import pandas as pd

df = pd.DataFrame({
    'Term Dummy':[0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32],
    'Shock':[131.759276601612,-5.28111953539055,-5.30412333137685,6.19553924065018,-5.97658803726517,-7.8325986545673,-9.50784210778306,-15.7385664305344,-23.3182508381464,-29.4897840376819,-31.467551725682,-33.4723203935889,-34.6650947285782,-35.7471724234754,-36.4799776375108,-37.3264043303424,-37.4155331344124,-37.8155991350952,-38.7550833588797,-38.3608088160098,-36.7211814243519,-35.7477615422699,-34.1458248652337,-32.8287847811565,-31.4018236645802,-29.9742754473972,-28.6193854123123,-24.90985538625,-21.3217573325541,-18.7350606702909,-16.0799516664911,-16.1549305201347,-16.1433168994669],
    'Weight':[1924,41170,120247,289092,311692,50265,127579,38255,225164,300420,96928,189792,177827,511969,417120,17840,72257,160679,89074,186051,102120,53770,662958,100838,765414,820977,533239,113092,60063,174082,238152,215960,115665]
})

其中Shock是待分割的序列,Weight为每个点的权重。

预期分割示例(n=6)

当n=6时,一种可行的分割方案(按索引/Term Dummy分组):

[
[0],
[1,2],
[3],
[4,5,6,7,8,9,10,11,12,13,14,15,16,17,18],
[19,20,21,22,23,24,25,26,27,28,29],
[30,31,32]
]

相关图示说明

  • Shock曲线与分割示例:展示Shock序列的走势以及n=6时的分段分割效果,第一个点单独成段,后续点按趋势分成5个连续子序列
  • 多曲线对比图:顶部为3条不同的Shock曲线,下方为对应的Sensitivity(标注为Delta)和Impact曲线,Weight是Sensitivity的绝对值

解决方案

1. 改造ruptures库的Binseg算法支持权重

ruptures的Binseg默认使用无权重损失函数,可通过自定义损失函数引入权重:

import numpy as np
import ruptures as rpt

class WeightedL2(rpt.base.BaseCost):
    """自定义带权重的L2损失函数,适配ruptures分割算法"""
    def __init__(self, weights=None):
        self.weights = weights
        super().__init__()

    def fit(self, signal):
        self.signal = signal
        return self

    def error(self, start, end):
        """计算[start, end)区间内的加权L2损失"""
        segment = self.signal[start:end]
        weights = self.weights[start:end] if self.weights is not None else np.ones_like(segment)
        x = np.arange(start, end)
        w_sqrt = np.sqrt(weights)
        
        # 加权最小二乘拟合线性函数
        A = np.vstack([w_sqrt * x, w_sqrt]).T
        b = w_sqrt * segment
        coeffs, _, _, _ = np.linalg.lstsq(A, b, rcond=None)
        
        y_pred = coeffs[0] * x + coeffs[1]
        return np.sum(weights * (segment - y_pred)**2)

# 带权重的Binseg分割函数
def weighted_binseg(series, weights, n_segments):
    signal = series.values
    cost_func = WeightedL2(weights=weights.values).fit(signal)
    model = rpt.Binseg(cost=cost_func)
    model.fit(signal)
    breaks = model.predict(n_bkps=n_segments-1)  # 分割点数量=分段数-1
    
    # 转换为分组索引
    groups = []
    prev = 0
    for b in breaks:
        groups.append(list(range(prev, b)))
        prev = b
    return groups

# 调用示例
n = 6
segments = weighted_binseg(df['Shock'], df['Weight'], n)
print(segments)

2. 用scipy+动态规划实现带权重的分段分割

若不想依赖ruptures,可通过动态规划结合加权线性拟合实现:

import numpy as np
from scipy.optimize import curve_fit

def weighted_linear_fit(x, y, w):
    """带权重的线性拟合,返回参数和加权残差平方和"""
    def linear_func(x, a, b):
        return a * x + b
    popt, _ = curve_fit(linear_func, x, y, sigma=np.sqrt(1/w))
    y_pred = linear_func(x, *popt)
    return popt, np.sum(w * (y - y_pred)**2)

def dynamic_programming_weighted_segmentation(x, y, w, n_segments):
    """动态规划求解带权重的分段线性分割,最小化总残差"""
    n = len(y)
    # 预计算所有区间的加权残差
    cost_matrix = np.zeros((n, n))
    for i in range(n):
        for j in range(i+1, n+1):
            _, cost = weighted_linear_fit(x[i:j], y[i:j], w[i:j])
            cost_matrix[i][j-1] = cost

    # 动态规划表初始化与填充
    dp = np.zeros((n_segments, n))
    path = np.zeros((n_segments, n), dtype=int)
    
    # 1段的情况
    for j in range(n):
        dp[0][j] = cost_matrix[0][j]
    
    # 多段情况
    for k in range(1, n_segments):
        for j in range(k, n):
            min_cost = np.inf
            best_i = k-1
            for i in range(k-1, j):
                current_cost = dp[k-1][i] + cost_matrix[i+1][j]
                if current_cost < min_cost:
                    min_cost = current_cost
                    best_i = i
            dp[k][j] = min_cost
            path[k][j] = best_i

    # 回溯获取分割点
    breaks = []
    current = n-1
    for k in range(n_segments-1, 0, -1):
        current = path[k][current]
        breaks.append(current+1)
    breaks.append(n)
    breaks = sorted(breaks)

    # 转换为分组索引
    groups = []
    prev = 0
    for b in breaks:
        groups.append(list(range(prev, b)))
        prev = b
    return groups

# 调用示例
x = df['Term Dummy'].values
y = df['Shock'].values
w = df['Weight'].values
segments = dynamic_programming_weighted_segmentation(x, y, w, 6)
print(segments)

3. 异常值处理建议

针对可能存在的异常值,可采取以下方式:

  • 对权重极低的点(如示例中第一个点),直接单独设为一个分段
  • 在损失函数中替换为鲁棒损失(如Huber损失),降低异常值对拟合的影响

性能验证

两种方案针对40个点、n=8的场景,单次处理耗时均远低于1.5秒,循环20次总耗时完全符合要求。


内容的提问来源于Stack Exchange,提问作者Mathieu

相关产品推荐
方舟 Agent Plan

超全模态模型 × Harness 升级,最新支持 Deepseek-V4.1-Flash、GLM-5.3 系列、Doubao-Seedream-5.0-pro、Kimi-K3 (部分), 限时 9.9 元起

最近更新时间:2026.07.16 11:31:57