如何基于最优点组自动计算含噪数组的平滑导数？

我有如下一组数据点：

import numpy as np
import matplotlib.pyplot as plt

x = np.array([1,1,1,1,2,3,4,4,4,4,4.5,5,5.5,6,6.5,7,7.5,8,8.1,8.2,8.3,8.4,8.5,8.5,8.5,8.5])
y = np.linspace(10,20,len(x))

plt.plot(x*0.1,y,'o', label='x')
plt.xlabel('x * 0.1')
plt.ylabel('y')

图1是该数据的可视化结果。我希望计算导数时，能自动找到具有共同增长趋势的最优点组，最终得到类似图2的结果——对无噪数据用numpy.diff做点到点求导可以实现这个效果，但在含随机噪声的真实数据中，numpy.diff会得到像图3那样的糟糕结果。请问能否针对含噪数据得到图2那样的平滑导数结果？也就是不用点到点求导，而是通过算法自动识别最优的点组来计算平滑导数。

解决方案 #

1. 分段线性拟合 + 断点检测 #

先检测数据中的趋势变化断点，再对每个分段做线性拟合，用拟合直线的斜率作为该段的导数，能精准匹配分段趋势：

import numpy as np
import matplotlib.pyplot as plt
from sklearn.metrics import mean_squared_error

# 模拟含噪数据
x_noisy = x + np.random.normal(0, 0.2, size=len(x))
y = np.linspace(10,20,len(x))

# 滑动窗口计算残差，识别趋势断点
window_size = 3
residuals = []
for i in range(len(x_noisy)-window_size):
    seg_x = x_noisy[i:i+window_size]
    seg_y = y[i:i+window_size]
    coeffs = np.polyfit(seg_y, seg_x, 1)
    pred_x = np.polyval(coeffs, seg_y)
    residuals.append(mean_squared_error(seg_x, pred_x))

# 基于残差阈值筛选断点
threshold = np.mean(residuals) + 2*np.std(residuals)
breakpoints = [i+window_size//2 for i, res in enumerate(residuals) if res > threshold]
breakpoints = [0] + breakpoints + [len(x_noisy)-1]

# 分段拟合并计算导数
derivatives = []
for i in range(len(breakpoints)-1):
    start, end = breakpoints[i], breakpoints[i+1]+1
    seg_y, seg_x = y[start:end], x_noisy[start:end]
    coeffs = np.polyfit(seg_y, seg_x, 1)
    slope = coeffs[0] * 0.1  # 对应x*0.1的导数
    derivatives.extend([slope]*(end-start))

# 可视化
plt.figure(figsize=(12,6))
plt.subplot(121)
plt.plot(x_noisy*0.1, y, 'o', label='含噪数据')
plt.xlabel('x * 0.1')
plt.ylabel('y')
plt.legend()

plt.subplot(122)
plt.plot(y, derivatives, '-', label='分段平滑导数')
plt.xlabel('y')
plt.ylabel('导数 (dx/dy * 0.1)')
plt.legend()
plt.tight_layout()
plt.show()

2. 滑动窗口加权平均导数 #

用滑动窗口对局部范围内的点计算导数并加权平均，快速抑制噪声，适合趋势变化平缓的数据：

def smooth_derivative(x, y, window_size=5):
    derivatives = np.zeros_like(x)
    half_win = window_size//2
    for i in range(len(x)):
        start = max(0, i-half_win)
        end = min(len(x), i+half_win+1)
        seg_x, seg_y = x[start:end], y[start:end]
        coeffs = np.polyfit(seg_y, seg_x, 1)
        derivatives[i] = coeffs[0] * 0.1
    return derivatives

# 应用到含噪数据
deriv_smooth = smooth_derivative(x_noisy, y, window_size=5)

# 可视化
plt.plot(y, deriv_smooth, '-', label='滑动窗口平滑导数')
plt.xlabel('y')
plt.ylabel('导数')
plt.legend()
plt.show()

3. 贝叶斯分段线性模型 #

通过贝叶斯方法自动推断分段数和断点位置，适合趋势复杂、需要精准自动化识别的场景：

import pymc3 as pm
import arviz as az

with pm.Model() as model:
    # 先验定义：断点数量、位置、分段斜率与截距
    k = pm.DiscreteUniform('k', lower=1, upper=5)
    changepoints = pm.Uniform('changepoints', lower=y.min(), upper=y.max(), shape=k)
    slopes = pm.Normal('slopes', mu=0, sd=1, shape=k+1)
    intercepts = pm.Normal('intercepts', mu=0, sd=10, shape=k+1)
    
    # 定义分段线性函数
    def linear_model(y_val):
        idx = pm.math.sum(changepoints < y_val)
        return slopes[idx] * y_val + intercepts[idx]
    
    # 似然函数
    pm.Normal('x_obs', mu=linear_model(y), sd=0.2, observed=x_noisy)
    
    # MCMC采样
    trace = pm.sample(2000, tune=1000, cores=2)

# 提取后验均值计算导数
post_slopes = az.summary(trace)['mean'].filter(regex='slopes')
changepoints_mean = np.sort(trace['changepoints'].mean(axis=0))
deriv_bayes = np.zeros_like(y)

for i in range(len(changepoints_mean)+1):
    if i == 0:
        mask = y <= changepoints_mean[i]
    elif i == len(changepoints_mean):
        mask = y > changepoints_mean[-1]
    else:
        mask = (y > changepoints_mean[i-1]) & (y <= changepoints_mean[i])
    deriv_bayes[mask] = post_slopes[f'slopes[{i}]'] * 0.1

# 可视化
plt.plot(y, deriv_bayes, '-', label='贝叶斯分段导数')
plt.xlabel('y')
plt.ylabel('导数')
plt.legend()
plt.show()

内容的提问来源于stack exchange，提问作者Lucas Oliveira