判断峰值-可排除噪声毛刺

峰值判定算法

前言

对于疫情来说，我们有了每个时刻的现存确诊人数数据(Confirmed)，也就是一个一维的数组，记为C（t）。想判断有没二次爆发，那么就需要判断这个序列是不是出现了超过一个峰值。当然啦，这只是这个算法的一个应用而已，这个算法可以用来判断其他一维数组的峰值个数及位置。

算法介绍

由于传染过程会有些波动，现存确诊人数C（t）曲线可能出现毛刺，若直接以C(t)>C(t-1)和C(t)>C(t+1)为条件的话，会出现多判，所以需要对原始的数据进行光滑处理。另外，有些小幅度的感染规模上升，是不能作为峰值的，因此，我们需要增加阈值。
具体操作步骤如下：

• 1， 归一化/计算密度
  这一步并非求峰值的必须，可以先对数据进行归一化，或者将人口数据转化为人口密度数据。
• 2， 光滑处理：
  采用均值滤波对原始进行处理。可根据经验设置滤波器窗口大小window_size，对t时刻的数据取其前后各window_size/2个数据（包含自身）求和，并将平均值作为C_filtered（t）的值，C_filtered（t）是光滑处理后的数据。
• 3， 判断峰值：
  峰值需满足一下两个条件：
  1） C_filtered（t）> C_filtered（t-1）且C_filtered（t）> C_filtered（t+1）
  2） C_filtered(t)>threshold，阈值可以根据情况取，比如说感染人口大于10000。
  另外，需要避免平峰：这里通过峰值宽度peak_width来排除相近的峰值。如果有两个峰值的间隔小于peak_width，则认为他们只有一个有效。
  效果如下图所示。
  窗口大小设置为window size=20，峰值宽度peak width = 10 。（窗口大小在20~40之间均可）

效果

设置参数：

• window_size = 20
• threshold = 0.01，即：对于一百万的人口来说，有一万人意思感染才能作为峰值。
• peak_width = 10

OK，让我们看看效果如何。

上图是未处理的数据，横坐标是时间t，纵坐标是感染人口C（t），右边是对左边的局部放大，可以看到曲线并不平滑，如果直接取峰值的话，可能将毛刺当做峰值。我们可以看出，有两个峰值。
下图是光滑处理后的数据，纵坐标是感染密度。
我用红笔标注了3个位置，我们不希望把3号也当做峰值，因为它只是轻微的起伏，并不能算作爆发。所以需要阈值来排除它，当然啦，仁者见仁智者见智，也许你觉得它是峰值，那么就把阈值调小一些吧。

最后程序是按照上面的说明，程序判断出了两个峰值，t=63, 216这两个位置，也就是图中标注1,2的位置。

源码

Python写的，有一个读取csv文件的代码，不属于本文的主题，主要是读取我存的数据。而且我存的数据格式不太规则，就自己写了。啰嗦一句，我数据保存的格式是，一行一个时间序列，每一行的数据长度不相同。

import numpy as np
import matplotlib.pyplot as plt


def readData(filepath):
    """
    读取csv文件，返回所有的时间序列，该csv文件每一行存储了一个时间序列。
    :param: filepath 文件名字符串
    :return: all_data 大概有几千条时间序列，每个序列是一个一维数组，
    """
    all_data = []

    with open(filepath, encoding='UTF-8') as f:
        sdata = f.readlines()
        all_count = len(sdata)