带权重的分段线性回归序列分割方案咨询(已知分段数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
相关产品推荐
相关产品推荐

