如何平滑数据集的曲线
- 2025-02-10 08:57:00
- admin 原创
- 55
问题描述：
假设我们有一个数据集，其大致内容如下：
import numpy as np
x = np.linspace(0,2*np.pi,100)
y = np.sin(x) + np.random.random(100) * 0.2
因此，我们有 20% 的数据集变化。我的第一个想法是使用UnivariateSplinescipy 函数，但问题是它不能很好地考虑小噪音。如果考虑频率，背景比信号小得多，因此仅截断的样条可能是一个想法，但这将涉及来回傅立叶变换，这可能会导致不良行为。另一种方法是移动平均，但这也需要正确选择延迟。
关于如何解决这个问题有任何提示/书籍或链接吗？
解决方案 1：
我更喜欢Savitzky-Golay 滤波器。它在 scipy 中可用。它使用最小二乘法将数据的一个小窗口回归到多项式上，然后使用多项式估计窗口中心的点。最后将窗口向前移动一个数据点并重复该过程。这一直持续到每个点都相对于其邻居进行了最佳调整。即使对于来自非周期性和非线性源的噪声样本，它也能很好地工作。
这是一个详尽的食谱示例，尽管现在已经过时了。注意：我省略了定义savitzky_golay()函数的代码，因为您可以从我上面链接的食谱示例中复制/粘贴它。
import numpy as np
import matplotlib.pyplot as plt
x = np.linspace(0,2*np.pi,100)
y = np.sin(x) + np.random.random(100) * 0.2
yhat = savitzky_golay(y, 51, 3) # window size 51, polynomial order 3
plt.plot(x,y)
plt.plot(x,yhat, color='red')
plt.show()
更新：我注意到我链接的食谱示例已被删除。幸运的是，正如@dodohjk指出的那样， Savitzky-Golay 过滤器已并入SciPy 库中（感谢@bicarlsen提供更新的链接）。要使用 SciPy 源调整上述代码，请键入：
from scipy.signal import savgol_filter
yhat = savgol_filter(y, 51, 3) # window size 51, polynomial order 3
解决方案 2：
编辑：看看这个答案。使用np.cumsum比np.convolve
基于移动平均框（通过卷积）来平滑我使用的一种快速而粗糙的数据的方法：
x = np.linspace(0,2*np.pi,100)
y = np.sin(x) + np.random.random(100) * 0.8
def smooth(y, box_pts):
box = np.ones(box_pts)/box_pts
y_smooth = np.convolve(y, box, mode='same')
return y_smooth
plot(x, y,'o')
plot(x, smooth(y,3), 'r-', lw=2)
plot(x, smooth(y,19), 'g-', lw=2)
解决方案 3：
如果您对周期性信号的“平滑”版本感兴趣（如您的示例），那么 FFT 是正确的选择。进行傅里叶变换并减去低贡献频率：
import numpy as np
import scipy.fftpack
N = 100
x = np.linspace(0,2*np.pi,N)
y = np.sin(x) + np.random.random(N) * 0.2
w = scipy.fftpack.rfft(y)
f = scipy.fftpack.rfftfreq(N, x[1]-x[0])
spectrum = w**2
cutoff_idx = spectrum < (spectrum.max()/5)
w2 = w.copy()
w2[cutoff_idx] = 0
y2 = scipy.fftpack.irfft(w2)
即使您的信号不是完全周期性的，这也将很好地减去白噪声。有许多类型的滤波器可供使用（高通、低通等...），合适的滤波器取决于您要寻找的内容。
解决方案 4：
这个问题已经得到了彻底的回答，所以我认为对所提出的方法进行运行时分析会很有趣（至少对我来说是这样）。我还将研究这些方法在嘈杂数据集的中心和边缘的行为。
总结
| runtime in s | runtime in s
method | python list | numpy array
--------------------|--------------|------------
kernel regression | 23.93405 | 22.75967
lowess | 0.61351 | 0.61524
naive average | 0.02485 | 0.02326
others* | 0.00150 | 0.00150
fft | 0.00021 | 0.00021
numpy convolve | 0.00017 | 0.00015
*savgol with different fit functions and some numpy methods
核回归的扩展性很差，Lowess 稍快一些，但两者都能产生平滑的曲线。Savgol 的速度处于中间水平，可以产生跳跃和平滑的输出，具体取决于多项式的等级。FFT 速度极快，但仅适用于周期性数据。
使用 numpy 的移动平均方法速度更快，但显然会生成带有步骤的图形。
设置
我生成了 1000 个正弦曲线形状的数据点：
size = 1000
x = np.linspace(0, 4 * np.pi, size)
y = np.sin(x) + np.random.random(size) * 0.2
data = {"x": x, "y": y}
我将它们传递到一个函数中来测量运行时间并绘制最终的拟合结果：
def test_func(f, label): # f: function handle to one of the smoothing methods
start = time()
for i in range(5):
arr = f(data["y"], 20)
print(f"{label:26s} - time: {time() - start:8.5f} ")
plt.plot(data["x"], arr, "-", label=label)
我测试了许多不同的平滑函数。arr是需要平滑的 y 值数组和span平滑参数。越低，拟合效果越接近原始数据，越高，得到的曲线越平滑。
def smooth_data_convolve_my_average(arr, span):
re = np.convolve(arr, np.ones(span * 2 + 1) / (span * 2 + 1), mode="same")
# The "my_average" part: shrinks the averaging window on the side that
# reaches beyond the data, keeps the other side the same size as given
# by "span"
re[0] = np.average(arr[:span])
for i in range(1, span + 1):
re[i] = np.average(arr[:i + span])
re[-i] = np.average(arr[-i - span:])
return re
def smooth_data_np_average(arr, span): # my original, naive approach
return [np.average(arr[val - span:val + span + 1]) for val in range(len(arr))]
def smooth_data_np_convolve(arr, span):
return np.convolve(arr, np.ones(span * 2 + 1) / (span * 2 + 1), mode="same")
def smooth_data_np_cumsum_my_average(arr, span):
cumsum_vec = np.cumsum(arr)
moving_average = (cumsum_vec[2 * span:] - cumsum_vec[:-2 * span]) / (2 * span)
# The "my_average" part again. Slightly different to before, because the
# moving average from cumsum is shorter than the input and needs to be padded
front, back = [np.average(arr[:span])], []
for i in range(1, span):
front.append(np.average(arr[:i + span]))
back.insert(0, np.average(arr[-i - span:]))
back.insert(0, np.average(arr[-2 * span:]))
return np.concatenate((front, moving_average, back))
def smooth_data_lowess(arr, span):
x = np.linspace(0, 1, len(arr))
return sm.nonparametric.lowess(arr, x, frac=(5*span / len(arr)), return_sorted=False)
def smooth_data_kernel_regression(arr, span):
# "span" smoothing parameter is ignored. If you know how to
# incorporate that with kernel regression, please comment below.
kr = KernelReg(arr, np.linspace(0, 1, len(arr)), 'c')
return kr.fit()[0]
def smooth_data_savgol_0(arr, span):
return savgol_filter(arr, span * 2 + 1, 0)
def smooth_data_savgol_1(arr, span):
return savgol_filter(arr, span * 2 + 1, 1)
def smooth_data_savgol_2(arr, span):
return savgol_filter(arr, span * 2 + 1, 2)
def smooth_data_fft(arr, span): # the scaling of "span" is open to suggestions
w = fftpack.rfft(arr)
spectrum = w ** 2
cutoff_idx = spectrum < (spectrum.max() * (1 - np.exp(-span / 2000)))
w[cutoff_idx] = 0
return fftpack.irfft(w)
结果
速度
运行时超过 1000 个元素，在 python 列表以及 numpy 数组上进行测试以保存值。
method | python list | numpy array
--------------------|-------------|------------
kernel regression | 23.93405 s | 22.75967 s
lowess | 0.61351 s | 0.61524 s
numpy average | 0.02485 s | 0.02326 s
savgol 2 | 0.00186 s | 0.00196 s
savgol 1 | 0.00157 s | 0.00161 s
savgol 0 | 0.00155 s | 0.00151 s
numpy convolve + me | 0.00121 s | 0.00115 s
numpy cumsum + me | 0.00114 s | 0.00105 s
fft | 0.00021 s | 0.00021 s
numpy convolve | 0.00017 s | 0.00015 s
尤其kernel regression是计算超过 1k 个元素时非常慢，lowess当数据集变得更大时也会失败。numpy convolve并且fft特别快。我没有研究增加或减少样本大小时的运行时行为（O（n））。
边缘行为
我将把这一部分分成两部分，以便于理解图像。
基于 Numpy 的方法 + savgol 0：
这些方法计算数据的平均值，图形没有经过平滑处理。numpy.cumsum当用于计算平均值的窗口不接触数据的边缘时，它们（除了）都会产生相同的图形。与的差异numpy.cumsum很可能是由于窗口大小的“偏差”错误造成的。
当该方法需要处理较少的数据时，会出现不同的边缘行为：
savgol 0：以常数继续到数据的边缘（分别savgol 1以savgol 2一条线和抛物线结束）numpy average：当窗口到达数据的左侧时停止，并用填充数组中的那些位置，与右侧的方法Nan行为相同my_averagenumpy convolve：非常准确地跟踪数据。我怀疑当窗口的一侧到达数据边缘时，窗口大小会对称减小my_average/me：我自己实现的方法，因为我对其他方法不满意。只需将超出数据范围的窗口部分缩小到数据的边缘，但将另一侧的窗口保持为给定的原始大小span
复杂的方法：
这些方法都以与数据的良好拟合而结束。savgol 1以一条线、savgol 2一条抛物线结束。
曲线行为
展示数据中间不同方法的行为。
不同的savgol和average过滤器产生粗糙的线条，，lowess并fft产生kernel regression平滑的拟合。lowess当数据发生变化时似乎会偷工减料。
动机
我有一个 Raspberry Pi 用来记录数据，但可视化却是个小挑战。除了 RAM 使用情况和以太网流量之外，所有数据点都只以离散步骤记录，并且/或者本身就很嘈杂。例如，温度传感器只输出整度，但连续测量之间的差异最多为两度。从这样的散点图中无法获得有用的信息。因此，为了可视化数据，我需要某种计算成本不太高且能产生移动平均值的方法。我还希望数据边缘有良好的行为，因为这在查看实时数据时尤其会影响最新信息。我决定采用这种numpy convolve方法my_average来改善边缘行为。
解决方案 5：
对数据进行移动平均可以消除噪音，请参阅此答案以了解如何操作。
如果您想使用LOWESS来拟合您的数据（它类似于移动平均线但更为复杂），您可以使用statsmodels库来实现：
import numpy as np
import pylab as plt
import statsmodels.api as sm
x = np.linspace(0,2*np.pi,100)
y = np.sin(x) + np.random.random(100) * 0.2
lowess = sm.nonparametric.lowess(y, x, frac=0.1)
plt.plot(x, y, '+')
plt.plot(lowess[:, 0], lowess[:, 1])
plt.show()
最后，如果您知道信号的函数形式，您可以根据数据拟合曲线，这可能是最好的做法。
解决方案 6：
另一种选择是在statsmodels中使用KernelReg：
from statsmodels.nonparametric.kernel_regression import KernelReg
import numpy as np
import matplotlib.pyplot as plt
x = np.linspace(0,2*np.pi,100)
y = np.sin(x) + np.random.random(100) * 0.2
# The third parameter specifies the type of the variable x;
# 'c' stands for continuous
kr = KernelReg(y,x,'c')
plt.plot(x, y, '+')
y_pred, y_std = kr.fit(x)
plt.plot(x, y_pred)
plt.show()
解决方案 7：
SciPy Cookbook对一维信号平滑的清晰定义向您展示了其工作原理。
捷径：
import numpy
def smooth(x,window_len=11,window='hanning'):
"""smooth the data using a window with requested size.
This method is based on the convolution of a scaled window with the signal.
The signal is prepared by introducing reflected copies of the signal
(with the window size) in both ends so that transient parts are minimized
in the begining and end part of the output signal.
input:
x: the input signal
window_len: the dimension of the smoothing window; should be an odd integer
window: the type of window from 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'
flat window will produce a moving average smoothing.
output:
the smoothed signal
example:
t=linspace(-2,2,0.1)
x=sin(t)+randn(len(t))*0.1
y=smooth(x)
see also:
numpy.hanning, numpy.hamming, numpy.bartlett, numpy.blackman, numpy.convolve
scipy.signal.lfilter
TODO: the window parameter could be the window itself if an array instead of a string
NOTE: length(output) != length(input), to correct this: return y[(window_len/2-1):-(window_len/2)] instead of just y.
"""
if x.ndim != 1:
raise ValueError, "smooth only accepts 1 dimension arrays."
if x.size < window_len:
raise ValueError, "Input vector needs to be bigger than window size."
if window_len<3:
return x
if not window in ['flat', 'hanning', 'hamming', 'bartlett', 'blackman']:
raise ValueError, "Window is on of 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'"
s=numpy.r_[x[window_len-1:0:-1],x,x[-2:-window_len-1:-1]]
#print(len(s))
if window == 'flat': #moving average
w=numpy.ones(window_len,'d')
else:
w=eval('numpy.'+window+'(window_len)')
y=numpy.convolve(w/w.sum(),s,mode='valid')
return y
from numpy import *
from pylab import *
def smooth_demo():
t=linspace(-4,4,100)
x=sin(t)
xn=x+randn(len(t))*0.1
y=smooth(x)
ws=31
subplot(211)
plot(ones(ws))
windows=['flat', 'hanning', 'hamming', 'bartlett', 'blackman']
hold(True)
for w in windows[1:]:
eval('plot('+w+'(ws) )')
axis([0,30,0,1.1])
legend(windows)
title("The smoothing windows")
subplot(212)
plot(x)
plot(xn)
for w in windows:
plot(smooth(xn,10,w))
l=['original signal', 'signal with noise']
l.extend(windows)
legend(l)
title("Smoothing a noisy signal")
show()
if __name__=='__main__':
smooth_demo()
解决方案 8：
对于我的一个项目，我需要创建时间序列建模的间隔，为了使过程更高效，我创建了tsmoothie：一个用于以矢量化方式进行时间序列平滑和异常值检测的 Python 库。
它提供了不同的平滑算法以及计算间隔的可能性。
这里我使用了ConvolutionSmoother但是你也可以测试其他的。
import numpy as np
import matplotlib.pyplot as plt
from tsmoothie.smoother import *
x = np.linspace(0,2*np.pi,100)
y = np.sin(x) + np.random.random(100) * 0.2
# operate smoothing
smoother = ConvolutionSmoother(window_len=5, window_type='ones')
smoother.smooth(y)
# generate intervals
low, up = smoother.get_intervals('sigma_interval', n_sigma=2)
# plot the smoothed timeseries with intervals
plt.figure(figsize=(11,6))
plt.plot(smoother.smooth_data[0], linewidth=3, color='blue')
plt.plot(smoother.data[0], '.k')
plt.fill_between(range(len(smoother.data[0])), low[0], up[0], alpha=0.3)
我还指出，tsmoothie 可以以矢量化方式对多个时间序列进行平滑处理
解决方案 9：
其中有一个简单的功能scipy.ndimage对我来说也很好用：
from scipy.ndimage import uniform_filter1d
y_smooth = uniform_filter1d(y,size=15)
解决方案 10：
使用移动平均线，一种快速的方法（也适用于非双射函数）是
def smoothen(x, winsize=5):
return np.array(pd.Series(x).rolling(winsize).mean())[winsize-1:]
该代码基于https://towardsdatascience.com/data-smoothing-for-data-science-visualization-the-goldilocks-trio-part-1-867765050615。其中还讨论了更高级的解决方案。
解决方案 11：
你也可以使用这个：
def smooth(scalars, weight = 0.8): # Weight between 0 and 1
return [scalars[i] * weight + (1 - weight) * scalars[i+1] for i in range(len(scalars)) if i < len(scalars)-1]
解决方案 12：
如果你正在绘制时间序列图，并且使用 mtplotlib 绘制图形，那么使用中值法来平滑图形
smotDeriv = timeseries.rolling(window=20, min_periods=5, center=True).median()
timeseries您传递的数据集在哪里，您可以对其进行修改windowsize以使其更加平滑。
