zoukankan      html  css  js  c++  java
  • Python 预测[周期性时间序列]

    1、背景

    公司平台上有不同的api,供内部或外部调用,这些api承担着不同的功能,如查询账号、发版、抢红包等等。日志会记录下每分钟某api被访问了多少次,即一个api每天会有1440条记录(1440分钟),将每天的数据连起来观察,有点类似于股票走势的意思。我想通过前N天的历史数据预测出第N+1天的流量访问情况,预测值即作为合理参考,供新一天与真实值做实时对比。当真实流量跟预测值有较大出入,则认为有异常访问,触发报警。

    2、数据探索

    我放了一份样例数据在data文件夹下,
    看一下数据大小和结构

    data = pd.read_csv(filename)
    print('size: ',data.shape)
    print(data.head())
    
     
    data.png

    数据大小:
    共10080条记录,即10080分钟,七天的数据。
    字段含义:
    date:时间,单位分钟
    count:该分钟该api被访问的次数

    画图看一下序列的走势:(一些画图等探索类的方法放在了test_stationarity.py 文件中,包含时间序列图,移动平均图,有兴趣的可以自己尝试下)。

    def draw_ts(timeseries):
        timeseries.plot()
        plt.show()
    
    data = pd.read_csv(path)
    data = data.set_index('date')
    data.index = pd.to_datetime(data.index)
    ts = data['count']
    draw_ts(ts)
    
     
    序列.png

    看这糟心的图,那些骤降为0的点这就是我遇到的第一个坑,我当初一拿到这份数据就开始做了。后来折腾了好久才发现,那些骤降为0的点是由于数据缺失,ETL的同学自动补零造成的,沟通晚了(TДT)。

    把坑填上,用前后值的均值把缺失值补上,再看一眼:


     
    填充好缺失值的序列.png

    发现这份数据有这样几个特点,在模型设计和数据预处理的时候要考虑到:

    1、这是一个周期性的时间序列,数值有规律的以天为周期上下波动,图中这个api,在每天下午和晚上访问较为活跃,在早上和凌晨较为稀少。在建模之前需要做分解。

    2、我的第二个坑:数据本身并不平滑,骤突骤降较多,而这样是不利于预测的,毕竟模型需要学习好正常的序列才能对未知数据给出客观判断,否则会出现频繁的误报,令气氛变得十分尴尬( ´Д`),所以必须进行平滑处理。

    3、这只是一个api的序列图,而不同的api的形态差距是很大的,毕竟承担的功能不同,如何使模型适应不同形态的api也是需要考虑的问题。

    3、预处理

    3.1 划分训练测试集

    前六天的数据做训练,第七天做测试集。

    class ModelDecomp(object):
        def __init__(self, file, test_size=1440):
            self.ts = self.read_data(file)
            self.test_size = test_size
            self.train_size = len(self.ts) - self.test_size
            self.train = self.ts[:len(self.ts)-test_size]
            self.test = self.ts[-test_size:]
    

    3.2 对训练数据进行平滑处理

    消除数据的毛刺,可以用移动平均法,我这里没有采用,因为我试过发现对于我的数据来说,移动平均处理完后并不能使数据平滑,我这里采用的方法很简单,但效果还不错:把每个点与上一点的变化值作为一个新的序列,对这里边的异常值,也就是变化比较离谱的值剃掉,用前后数据的均值填充,注意可能会连续出现变化较大的点:

    def _diff_smooth(self, ts):
        dif = ts.diff().dropna() # 差分序列
        td = dif.describe() # 描述性统计得到:min,25%,50%,75%,max值
        high = td['75%'] + 1.5 * (td['75%'] - td['25%']) # 定义高点阈值,1.5倍四分位距之外
        low = td['25%'] - 1.5 * (td['75%'] - td['25%']) # 定义低点阈值,同上
    
        # 变化幅度超过阈值的点的索引
        forbid_index = dif[(dif > high) | (dif < low)].index 
        i = 0
        while i < len(forbid_index) - 1:
            n = 1 # 发现连续多少个点变化幅度过大,大部分只有单个点
            start = forbid_index[i] # 异常点的起始索引
            while forbid_index[i+n] == start + timedelta(minutes=n):
                n += 1
            i += n - 1
    
            end = forbid_index[i] # 异常点的结束索引
            # 用前后值的中间值均匀填充
            value = np.linspace(ts[start - timedelta(minutes=1)], ts[end + timedelta(minutes=1)], n)
            ts[start: end] = value
            i += 1
    
    self.train = self._diff_smooth(self.train)
    draw_ts(self.train)
    

    平滑后的训练数据:


     
    平滑后的训练序列.png

    3.3 将训练数据进行周期性分解

    采用statsmodels工具包:

    from statsmodels.tsa.seasonal import seasonal_decompose
    
    decomposition = seasonal_decompose(self.ts, freq=freq, two_sided=False)
    # self.ts:时间序列,series类型; 
    # freq:周期,这里为1440分钟,即一天; 
    # two_sided:观察下图2、4行图,左边空了一段,如果设为True,则会出现左右两边都空出来的情况,False保证序列在最后的时间也有数据,方便预测。
    
    self.trend = decomposition.trend
    self.seasonal = decomposition.seasonal
    self.residual = decomposition.resid
    decomposition.plot()
    plt.show()
    
     
    分解图.png

    第一行observed:原始数据;第二行trend:分解出来的趋势部分;第三行seasonal:周期部分;最后residual:残差部分。
    我采用的是seasonal_decompose的加法模型进行的分解,即 observed = trend + seasonal + residual,另还有乘法模型。在建模的时候,只针对trend部分学习和预测,如何将trend的预测结果加工成合理的最终结果?当然是再做加法,后面会详细写。

    4、模型

    4.1 训练

    对分解出来的趋势部分单独用arima模型做训练:

    def trend_model(self, order):
        self.trend.dropna(inplace=True)
        train = self.trend[:len(self.trend)-self.test_size]
        #arima的训练参数order =(p,d,q),具体意义查看官方文档,调参过程略。
        self.trend_model = ARIMA(train, order).fit(disp=-1, method='css')
    

    4.2 预测

    预测出趋势数据后,加上周期数据即作为最终的预测结果,但更重要的是,我们要得到的不是具体的值,而是一个合理区间,当真实数据超过了这个区间,则触发报警,误差高低区间的设定来自刚刚分解出来的残差residual数据:

    d = self.residual.describe()
    delta = d['75%'] - d['25%']
    self.low_error, self.high_error = (d['25%'] - 1 * delta, d['75%'] + 1 * delta)
    

    预测并完成最后的加法处理,得到第七天的预测值即高低置信区间:

        def predict_new(self):
            '''
            预测新数据
            '''
            #续接train,生成长度为n的时间索引,赋给预测序列
            n = self.test_size
            self.pred_time_index= pd.date_range(start=self.train.index[-1], periods=n+1, freq='1min')[1:]
            self.trend_pred= self.trend_model.forecast(n)[0]
            self.add_season()
    
    
        def add_season(self):
            '''
            为预测出的趋势数据添加周期数据和残差数据
            '''
            self.train_season = self.seasonal[:self.train_size]
            values = []
            low_conf_values = []
            high_conf_values = []
    
            for i, t in enumerate(self.pred_time_index):
                trend_part = self.trend_pred[i]
    
                # 相同时间点的周期数据均值
                season_part = self.train_season[
                    self.train_season.index.time == t.time()
                    ].mean()
    
                # 趋势 + 周期 + 误差界限
                predict = trend_part + season_part
                low_bound = trend_part + season_part + self.low_error
                high_bound = trend_part + season_part + self.high_error
    
                values.append(predict)
                low_conf_values.append(low_bound)
                high_conf_values.append(high_bound)
    
            # 得到预测值,误差上界和下界
            self.final_pred = pd.Series(values, index=self.pred_time_index, name='predict')
            self.low_conf = pd.Series(low_conf_values, index=self.pred_time_index, name='low_conf')
            self.high_conf = pd.Series(high_conf_values, index=self.pred_time_index, name='high_conf')
    

    4.3 评估:

    对第七天作出预测,评估的指标为均方根误差rmse,画图对比和真实值的差距:

        md = ModelDecomp(file=filename, test_size=1440)
        md.decomp(freq=1440)
        md.trend_model(order=(1, 1, 3)) # arima模型的参数order
        md.predict_new() 
        pred = md.final_pred
        test = md.test
    
        plt.subplot(211)
        plt.plot(md.ts) # 平滑过的训练数据加未做处理的测试数据
        plt.title(filename.split('.')[0])
    
        plt.subplot(212)
        pred.plot(color='blue', label='Predict') # 预测值
        test.plot(color='red', label='Original') # 真实值
        md.low_conf.plot(color='grey', label='low') # 低置信区间
        md.high_conf.plot(color='grey', label='high') # 高置信区间
    
        plt.legend(loc='best')
        plt.title('RMSE: %.4f' % np.sqrt(sum((pred.values - test.values) ** 2) / test.size))
        plt.tight_layout()
        plt.show()
    
     
    预测结果.png

    可以看到,均方根误差462.8,相对于原始数据几千的量级,还是可以的。测试数据中的两个突变的点,也超过了置信区间,能准确报出来。




    参考链接:https://www.jianshu.com/p/09e5218f58b4
    code链接:  https://github.com/scarlettgin/cyclical_series_predict

  • 相关阅读:
    xls与csv文件的区别
    青音,经典爱情语录
    win7用户账户自动登录方法汇总
    How to using Procedure found Lead Blocker
    FTS(3) BSD 库函数手册 遍历文件夹(二)
    FTS(3) BSD 库函数手册 遍历文件夹(一)
    DisplayMetrics类 获取手机显示屏的基本信息 包括尺寸、密度、字体缩放等信息
    About App Distribution 关于应用发布
    FTS(3) 遍历文件夹实例
    OpenCV 2.1.0 with Visual Studio 2008
  • 原文地址:https://www.cnblogs.com/zhukaijian/p/13203463.html
Copyright © 2011-2022 走看看