python时间序列峰值检测_算法 - 实时时间序列中的峰值信号检测d

算法 - 实时时间序列中的峰值信号检测d

更新：迄今为止表现最佳的算法就是这个算法。

该问题探讨了用于检测实时时间序列数据中的突然峰值的稳健算法。

请考虑以下数据集：

p = [1 1 1.1 1 0.9 1 1 1.1 1 0.9 1 1.1 1 1 0.9 1 1 1.1 1 1 1 1 1.1 0.9 1 1.1 1 1 0.9 1, ...

1.1 1 1 1.1 1 0.8 0.9 1 1.2 0.9 1 1 1.1 1.2 1 1.5 1 3 2 5 3 2 1 1 1 0.9 1 1 3, ...

2.6 4 3 3.2 2 1 1 0.8 4 4 2 2.5 1 1 1];

(Matlab格式，但它不是关于语言，而是关于算法)

你可以清楚地看到有三个大峰和一些小峰。 此数据集是问题所涉及的时间序列数据集类的特定示例。 这类数据集有两个一般特征：

基本噪音具有一般意义

存在明显偏离噪声的大“峰值”或“更高数据点”。

我们还假设如下：

峰的宽度不能预先确定

峰的高度明显偏离其他值

使用的算法必须计算实时(因此每个新数据点都会改变)

对于这种情况，需要构造触发信号的边界值。 但是，边界值不能是静态的，必须基于算法实时确定。

我的问题：什么是实时计算这些阈值的好算法？ 这种情况是否有特定的算法？ 什么是最知名的算法？

强大的算法或有用的见解都受到高度赞赏。 (可以用任何语言回答：它是关于算法的)

22个解决方案

181 votes

平滑的z-score算法(具有稳健阈值的峰值检测)

我构建了一种适用于这些类型数据集的算法。它基于分散原理：如果新的数据点是距离某些移动平均值的给定x个标准偏差，则该算法发信号(也称为z分数)。该算法非常稳健，因为它构造了单独的移动平均值和偏差，使得信号不会破坏阈值。因此，无论先前信号的量如何，未来信号都以大致相同的精度被识别。该算法需要3个输入：filteredY,avgFilter和stdFilter。例如，2的lag将使用最后5个观察来平滑数据。如果数据点与移动平均值相差3.5个标准偏差，则3.5的threshold将发出信号。并且一个0.5的influence给出了正常数据点具有一半影响的信号。同样，influence为0会完全忽略信号以重新计算新阈值。因此，0的影响是最稳健的选择(但假设平稳性);将影响选项置于1是最不稳健的。因此，对于非平稳数据，影响选项应该放在介于0和1之间的位置。

它的工作原理如下：

伪代码

# Let y be a vector of timeseries data of at least length lag+2

# Let mean() be a function that calculates the mean

# Let std() be a function that calculates the standard deviaton

# Let absolute() be the absolute value function

# Settings (the ones below are examples: choose what is best for your data)

set lag to 5; # lag 5 for the smoothing functions

set threshold to 3.5; # 3.5 standard deviations for signal

set influence to 0.5; # between 0 and 1, where 1 is normal influence, 0.5 is half

# Initialise variables

set signals to vector 0,...,0 of length of y; # Initialise signal results

set filteredY to y(1),...,y(lag) # Initialise filtered series

set avgFilter to null; # Initialise average filter

set stdFilter to null; # Initialise std. filter

set avgFilter(lag) to mean(y(1),...,y(lag)); # Initialise first value

set stdFilter(lag) to std(y(1),...,y(lag)); # Initialise first value

for i=lag+1,...,t do

if absolute(y(i) - avgFilter(i-1)) > threshold*stdFilter(i-1) then

if y(i) > avgFilter(i-1) then

set signals(i) to +1; # Positive signal

else

set signals(i) to -1; # Negative signal

end

# Make influence lower

set filteredY(i) to influence*y(i) + (1-influence)*filteredY(i-1);

else

set signals(i) to 0; # No signal

set filteredY(i) to y(i);

end

# Adjust the filters

set avgFilter(i) to mean(filteredY(i-lag),...,filteredY(i));

set stdFilter(i) to std(filteredY(i-lag),...,filteredY(i));

end

为您的数据选择好参数的经验法则可以在附录3(下面)中找到。

演示

该演示的Matlab代码可以在本答案的最后找到。 要使用该演示，只需运行它并通过单击上方图表自行创建时间序列。 在绘制filteredY个观测值后，算法开始工作。

附录1：算法的Matlab和R代码

Matlab代码

function [signals,avgFilter,stdFilter] = ThresholdingAlgo(y,lag,threshold,influence)

% Initialise signal results

signals = zeros(length(y),1);

% Initialise filtered series

filteredY = y(1:lag+1);

% Initialise filters

avgFilter(lag+1,1) = mean(y(1:lag+1));

stdFilter(lag+1,1) = std(y(1:lag+1));

% Loop over all datapoints y(lag+2),...,y(t)

for i=lag+2:length(y)

% If new value is a specified number of deviations away

if abs(y(i)-avgFilter(i-1)) > threshold*stdFilter(i-1)

if y(i) > avgFilter(i-1)

% Positive signal

signals(i) = 1;

else

% Negative signal

signals(i) = -1;

end

% Make influence lower

filteredY(i) = influence*y(i)+(1-influence)*filteredY(i-1);

else

% No signal

signals(i) = 0;

filteredY(i) = y(i);

end

% Adjust the filters

avgFilter(i) = mean(filteredY(i-lag:i));

stdFilter(i) = std(filteredY(i-lag:i));

end

% Done, now return results

end

例：

% Data

y = [1 1 1.1 1 0.9 1 1 1.1 1 0.9 1 1.1 1 1 0.9 1 1 1.1 1 1,...

1 1 1.1 0.9 1 1.1 1 1 0.9 1 1.1 1 1 1.1 1 0.8 0.9 1 1.2 0.9 1,...

1 1.1 1.2 1 1.5 1 3 2 5 3 2 1 1 1 0.9 1,...

1 3 2.6 4 3 3.2 2 1 1 0.8 4 4 2 2.5 1 1 1];

% Settings

lag = 30;

threshold = 5;

influence = 0;

% Get results

[signals,avg,dev] = ThresholdingAlgo(y,lag,threshold,influence);

figure; subplot(2,1,1); hold on;

x = 1:length(y); ix = lag+1:length(y);

area(x(ix),avg(ix)+threshold*dev(ix),'FaceColor',[0.9 0.9 0.9],'EdgeColor','none');

area(x(ix),avg(ix)-threshold*dev(ix),'FaceColor',[1 1 1],'EdgeColor','none');

plot(x(ix),avg(ix),'LineWidth',1,'Color','cyan','LineWidth',1.5);

plot(x(ix),avg(ix)+threshold*dev(ix),'LineWidth',1,'Color','green','LineWidth',1.5);

plot(x(ix),avg(ix)-threshold*dev(ix),'LineWidth',1,'Color','green','LineWidth',1.5);

plot(1:length(y),y,'b');

subplot(2,1,2);

stairs(signals,'r','LineWidth',1.5); ylim([-1.5 1.5]);

R代码

ThresholdingAlgo

signals

filteredY

avgFilter

stdFilter

avgFilter[lag]

stdFilter[lag]

for (i in (lag+1):length(y)){

if (abs(y[i]-avgFilter[i-1]) > threshold*stdFilter[i-1]) {

if (y[i] > avgFilter[i-1]) {

signals[i]

} else {

signals[i]

}

filteredY[i]

} else {

signals[i]

filteredY[i]

}

avgFilter[i]

stdFilter[i]

}

return(list("signals"=signals,"avgFilter"=avgFilter,"stdFilter"=stdFilter))

}

例：

# Data

y

1,1.1,1,1,1.1,1,0.8,0.9,1,1.2,0.9,1,1,1.1,1.2,1,1.5,1,3,2,5,3,2,1,1,1,0.9,1,1,3,

2.6,4,3,3.2,2,1,1,0.8,4,4,2,2.5,1,1,1)

lag

threshold

influence

# Run algo with lag = 30, threshold = 5, influence = 0

result

# Plot result

par(mfrow = c(2,1),oma = c(2,2,0,0) + 0.1,mar = c(0,0,2,1) + 0.2)

plot(1:length(y),y,type="l",ylab="",xlab="")

lines(1:length(y),result$avgFilter,type="l",col="cyan",lwd=2)

lines(1:length(y),result$avgFilter+threshold*result$stdFilter,type="l",col="green",lwd=2)

lines(1:length(y),result$avgFilter-threshold*result$stdFilter,type="l",col="green",lwd=2)

plot(result$signals,type="S",col="red",ylab="",xlab="",ylim=c(-1.5,1.5),lwd=2)

此代码(两种语言)将为原始问题的数据生成以下结果：

其他语言的实现：

Golang(Xeoncross)

Python(R Kiselev)

斯威夫特(我)

Groovy(JoshuaCWebDeveloper)

C ++(布拉德)

C ++(Animesh Pandey)

锈(swizard)

斯卡拉(迈克罗伯茨)

Kotlin(leoderprofi)

Ruby(Kimmo Lehto)

Fortran [用于共振检测](THo)

朱莉娅(马特营)

C＃(海洋空投)

附录2：Matlab演示代码(点击制作数据)

function [] = RobustThresholdingDemo()

%% SPECIFICATIONS

lag = 5; % lag for the smoothing

threshold = 3.5; % number of st.dev. away from the mean to signal

influence = 0.3; % when signal: how much influence for new data? (between 0 and 1)

% 1 is normal influence, 0.5 is half

%% START DEMO

DemoScreen(30,lag,threshold,influence);

end

function [signals,avgFilter,stdFilter] = ThresholdingAlgo(y,lag,threshold,influence)

signals = zeros(length(y),1);

filteredY = y(1:lag+1);

avgFilter(lag+1,1) = mean(y(1:lag+1));

stdFilter(lag+1,1) = std(y(1:lag+1));

for i=lag+2:length(y)

if abs(y(i)-avgFilter(i-1)) > threshold*stdFilter(i-1)

if y(i) > avgFilter(i-1)

signals(i) = 1;

else

signals(i) = -1;

end

filteredY(i) = influence*y(i)+(1-influence)*filteredY(i-1);

else

signals(i) = 0;

filteredY(i) = y(i);

end

avgFilter(i) = mean(filteredY(i-lag:i));

stdFilter(i) = std(filteredY(i-lag:i));

end

end

% Demo screen function

function [] = DemoScreen(n,lag,threshold,influence)

figure('Position',[200 100,1000,500]);

subplot(2,1,1);

title(sprintf(['Draw data points (%.0f max) [settings: lag = %.0f, '...

'threshold = %.2f, influence = %.2f]'],n,lag,threshold,influence));

ylim([0 5]); xlim([0 50]);

H = gca; subplot(2,1,1);

set(H, 'YLimMode', 'manual'); set(H, 'XLimMode', 'manual');

set(H, 'YLim', get(H,'YLim')); set(H, 'XLim', get(H,'XLim'));

xg = []; yg = [];

for i=1:n

try

[xi,yi] = ginput(1);

catch

return;

end

xg = [xg xi]; yg = [yg yi];

if i == 1

subplot(2,1,1); hold on;

plot(H, xg(i),yg(i),'r.');

text(xg(i),yg(i),num2str(i),'FontSize',7);

end

if length(xg) > lag

[signals,avg,dev] = ...

ThresholdingAlgo(yg,lag,threshold,influence);

area(xg(lag+1:end),avg(lag+1:end)+threshold*dev(lag+1:end),...

'FaceColor',[0.9 0.9 0.9],'EdgeColor','none');

area(xg(lag+1:end),avg(lag+1:end)-threshold*dev(lag+1:end),...

'FaceColor',[1 1 1],'EdgeColor','none');

plot(xg(lag+1:end),avg(lag+1:end),'LineWidth',1,'Color','cyan');

plot(xg(lag+1:end),avg(lag+1:end)+threshold*dev(lag+1:end),...

'LineWidth',1,'Color','green');

plot(xg(lag+1:end),avg(lag+1:end)-threshold*dev(lag+1:end),...

'LineWidth',1,'Color','green');

subplot(2,1,2); hold on; title('Signal output');

stairs(xg(lag+1:end),signals(lag+1:end),'LineWidth',2,'Color','blue');

ylim([-2 2]); xlim([0 50]); hold off;

end

subplot(2,1,1); hold on;

for j=2:i

plot(xg([j-1:j