实时时间序列数据中的峰值信号检测


更新:到目前为止表现最好的算法就是这一个


这个问题探讨了用于检测实时时间序列数据中突发峰值的强大算法。

考虑以下数据集:

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格式,但它不是关于语言,而是关于算法)

数据图

你可以清楚地看到有三个大峰和一些小峰。 该数据集是问题所关注的时间序列数据集类的具体示例。 这类数据集有两个一般特征:

  • 有一般意思的基本噪音
  • 有大量的“峰值”或“更高的数据点”明显偏离了噪声。
  • 我们还假设以下几点:

  • 预先不能确定峰的宽度
  • 峰的高度明显偏离其他值
  • 所使用的算法必须计算实时性(因此每个新数据点都要更改)
  • 对于这种情况,需要构建触发信号的边界值。 但是,边界值不能是静态的,必须基于算法实时确定。


    我的问题:什么是实时计算此类阈值的好算法? 有这种情况下的具体算法吗? 什么是最知名的算法?


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


    平滑的z分数算法(非常健壮的阈值算法)

    我已经构建了一种适用于这些类型的数据集的算法。 它基于色散原理:如果一个新的数据点与某个移动均值相距给定的x个标准偏差,算法信号(也称为z-分数)。 该算法非常强大,因为它构造了一个单独的移动均值和偏差,以便信号不会破坏阈值。 因此,无论先前的信号量如何,未来的信号都具有大致相同的精度。 该算法需要3个输入: lag = the lag of the moving window ,门threshold = the z-score at which the algorithm signalsinfluence = the influence (between 0 and 1) of new signals on the mean and standard deviation threshold = the z-score at which the algorithm signals influence = the influence (between 0 and 1) of new signals on the mean and standard deviation 。 例如, lag 5将使用最后5次观察来平滑数据。 如果数据点距移动平均值3.5个标准偏差,则threshold 3.5将发出信号。 0.5的influence给信号一半的正常数据点的影响。 同样,0的influence完全忽略了重新计算新阈值的信号:因此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
    

    演示

    演示鲁棒阈值算法

    此演示的Matlab代码可以在本答案的末尾找到。 要使用演示,只需运行它并通过单击上图创建一个时间序列。 该算法在绘制lag数量的观察值之后开始工作。


    附录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 <- function(y,lag,threshold,influence) {
      signals <- rep(0,length(y))
      filteredY <- y[0:lag]
      avgFilter <- NULL
      stdFilter <- NULL
      avgFilter[lag] <- mean(y[0:lag])
      stdFilter[lag] <- sd(y[0: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] <- 1;
          } else {
            signals[i] <- -1;
          }
          filteredY[i] <- influence*y[i]+(1-influence)*filteredY[i-1]
        } else {
          signals[i] <- 0
          filteredY[i] <- y[i]
        }
        avgFilter[i] <- mean(filteredY[(i-lag):i])
        stdFilter[i] <- sd(filteredY[(i-lag):i])
      }
      return(list("signals"=signals,"avgFilter"=avgFilter,"stdFilter"=stdFilter))
    }
    

    例:

    # Data
    y <- c(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)
    
    lag       <- 30
    threshold <- 5
    influence <- 0
    
    # Run algo with lag = 30, threshold = 5, influence = 0
    result <- ThresholdingAlgo(y,lag,threshold,influence)
    
    # 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)
    

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

    来自Matlab代码的阈值示例


    其他语言的实现:

  • Golang(Xeoncross)

  • Python(R Kiselev)

  • 斯威夫特(我)

  • Groovy(JoshuaCWebDeveloper)

  • C ++(brad)

  • C ++(Animesh Pandey)

  • 锈(swizard)

  • 斯卡拉(Mike Roberts)

  • Kotlin(leoderprofi)

  • Ruby(Kimmo Lehto)


  • 附录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]),yg([j-1:j]),'r'); plot(H,xg(j),yg(j),'r.');
            text(xg(j),yg(j),num2str(j),'FontSize',7);
        end
    end
    end
    

    警告:上面的代码每次运行时都会遍历所有数据点。 在实现这个代码时,确保将信号的计算分解成一个单独的函数(没有循环)。 然后,当新的数据点到达时,更新filteredYavgFilterstdFilter一次。 每次有新的数据点时(例如上面的例子),不要重新计算所有数据的信号,这将非常低效且速度慢!

    修改算法的其他方法(用于可能的改进)是:

  • 使用中位数而不是平均值
  • 使用MAD(中位数绝对偏差)代替标准偏差
  • 使用信令余量,所以信号不会经常切换
  • 改变影响参数的工作方式
  • 以不同的方式处理信号(不对称处理)

  • 如果你在某个地方使用这个功能,请记下我或这个答案。 如果您对此算法有任何疑问,请将它们发布在下面的评论中或通过LinkedIn与我联系。



    一种方法是基于以下观察来检测峰值:

  • 如果(y(t)> y(t-1))&&(y(t)> y(t + 1)),则时间t是峰值
  • 它通过等待上涨趋势结束来避免误报。 从某种意义上来说,它不会完全是“实时”的,它会错过一个dt的峰值。 灵敏度可以通过需要余量进行比较来控制。 噪声检测和检测时间延迟之间存在权衡。 您可以通过添加更多参数来丰富模型:

  • 如果(y(t)-y(t-dt)> m)&&(y(t)-y(t + dt)> m)
  • 其中dtm是用于控制灵敏度与时间延迟的参数

    这就是你用上述算法得到的结果: 在这里输入图像描述

    这里是用python重现剧情的代码:

    import numpy as np
    import matplotlib.pyplot as plt
    input = np.array([ 1. ,  1. ,  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. ,  1. ,  1. ,  1. ])
    signal = (input > np.roll(input,1)) & (input > np.roll(input,-1))
    plt.plot(input)
    plt.plot(signal.nonzero()[0], input[signal], 'ro')
    plt.show()
    

    通过设置m = 0.5 ,您可以得到一个只有一个误报的更清晰的信号: 在这里输入图像描述


    以下是平滑z-score算法的Python / numpy实现(请参阅上面的答案)。 你可以在这里找到要点。

    #!/usr/bin/env python
    # Implementation of algorithm from https://stackoverflow.com/a/22640362/6029703
    import numpy as np
    import pylab
    
    def thresholding_algo(y, lag, threshold, influence):
        signals = np.zeros(len(y))
        filteredY = np.array(y)
        avgFilter = [0]*len(y)
        stdFilter = [0]*len(y)
        avgFilter[lag - 1] = np.mean(y[0:lag])
        stdFilter[lag - 1] = np.std(y[0:lag])
        for i in range(lag, len(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
    
                filteredY[i] = influence * y[i] + (1 - influence) * filteredY[i-1]
                avgFilter[i] = np.mean(filteredY[(i-lag):i])
                stdFilter[i] = np.std(filteredY[(i-lag):i])
            else:
                signals[i] = 0
                filteredY[i] = y[i]
                avgFilter[i] = np.mean(filteredY[(i-lag):i])
                stdFilter[i] = np.std(filteredY[(i-lag):i])
    
        return dict(signals = np.asarray(signals),
                    avgFilter = np.asarray(avgFilter),
                    stdFilter = np.asarray(stdFilter))
    

    以下是在R / Matlab原始答案中产生相同图的相同数据集的测试

    # Data
    y = np.array([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
    lag = 30
    threshold = 5
    influence = 0
    
    # Run algo with settings from above
    result = thresholding_algo(y, lag=lag, threshold=threshold, influence=influence)
    
    # Plot result
    pylab.subplot(211)
    pylab.plot(np.arange(1, len(y)+1), y)
    
    pylab.plot(np.arange(1, len(y)+1),
               result["avgFilter"], color="cyan", lw=2)
    
    pylab.plot(np.arange(1, len(y)+1),
               result["avgFilter"] + threshold * result["stdFilter"], color="green", lw=2)
    
    pylab.plot(np.arange(1, len(y)+1),
               result["avgFilter"] - threshold * result["stdFilter"], color="green", lw=2)
    
    pylab.subplot(212)
    pylab.step(np.arange(1, len(y)+1), result["signals"], color="red", lw=2)
    pylab.ylim(-1.5, 1.5)
    
    链接地址: http://www.djcxy.com/p/26789.html

    上一篇: Peak signal detection in realtime timeseries data

    下一篇: Multiple Inheritance and calling super()