图书销量时间序列预测_数学建模_Prophet实现

前言

• 该题来源于数学建模校内选拔赛，给出了简单的数据，恰好我找到prophet这个神器，便用它来完成题目了，一开始只打算来做做第一题预测什么的，结果做着做着居然发现能只靠prophet库把每个小问给圆上。原来时间紧张代码写得东一块西一块，现在有时间了稍微整理了下代码发了上来。

• 说一下对prophet的使用感受。这模型还是非常强的，原理也不难，搞清楚使用后几行代码就能把预测图和时间序列分解图画出来，非常适合需要马上上手预测的任务。虽然给出的可调参数不少，不过需要调整的并不多，像是什么初始影响因子的值的，调了下感觉无论初始设置多少，训练的时候都会很快收敛到相同的解，不过也许是因为我的数据集比较简单。

• 比赛任务一要求对销售数据进行建模预测和模型评估，任务二要求寻求与销量走势相关的其他属性，任务三要求寻分析导致销量量骤升骤降的原因并用于修正原模型。

• 使用prophet解决思路也很简单，不过因为题目给的数据只有2年内逐月的数据，总共24个观察值，采样点太少，可能没法发挥出prophet的真实力量。

• 任务一预测已经是库里面自带的了，评估的话就把数据归一化后计算其预测误差就好；任务二的话，把prophet训练得到的趋势项在变点上的斜率提取出来和周期项的傅里叶展开系数结合起来作为特征，用层次聚类聚集出几个簇，找到类内和类均值比总误差较小的几类就得到走势相同的曲线了，这个任务的完成效果还是不错的；任务三就有点勉强了，大概思路是把对每个样本把骤升骤降的点标记处理，让prophet增加一个残差回归项去学习这些骤变情况（这部分感觉理解不是很到位）。

• 许久后回头看看这份代码建模报告, 感觉也没写啥(捂脸 , 好在校内选拔赛是能过了的

主要参考

• 主要参考的即为github上的发布连接 https://github.com/facebook/prophet

代码

库导入与函数设置

导库

import numpy as np
import pandas as pd
import matplotlib.pylab as plt
import seaborn as sns
from itertools import cycle
from fbprophet import Prophet # 0.7版本
from fbprophet.plot import add_changepoints_to_plot
import copy as cp
import datetime 
import time
import random

pd.set_option('max_columns', 50)
plt.style.use('bmh')
color_pal = plt.rcParams['axes.prop_cycle'].by_key()['color']
color_cycle = cycle(plt.rcParams['axes.prop_cycle'].by_key()['color'])

# 忽略警告
import warnings
warnings.filterwarnings('ignore')

# 绘制更清晰的图片
%config InlineBackend.figure_format = 'svg'

# 设置图片内的中文与负号显示
plt.rcParams['font.sans-serif']=['SimHei']
plt.rcParams['axes.unicode_minus'] = False

import os

INPUT_DIR = os.path.realpath('.')

展示函数

def ShowBiog(data,data_index):
#  展示商品的销售曲线
#  data是数据集，data_index是 索引号
    plt.figure(figsize=(10, 3))
    ax = fig.add_subplot(111)

    example = data.iloc[data_index]
    # 由于是完整的数据集中带有商品品类、出版时间等非销量数据，需去掉
    example[4:-1].plot(
            color=next(color_cycle),
           )
    plt.title('{}'.format(example['品种']))
#     ax.set_xticklabels('')
#     plt.legend('')
    plt.show()

取数据函数

def GetData(data,data_index):
#  该函数主要负责从data中取出索引号为data_index的数据
#  并生成为prophet能识别的形式
    data = data.drop(data.columns[0:4], axis=1)
    if 'SUM' in data.columns:
        data = data.drop(data.columns[-1], axis=1)
    data = data.T.reset_index()

    data0 = pd.DataFrame([data.iloc[:,0].values,data.iloc[:,data_index+1].values]).T
    data0.rename(columns={
   0:'ds', 1:'y'},inplace=True)
    return data0

训练函数

def FlagFunc(ds,train):
#  该函数主要是遍历数据，并数据发生骤升骤降时进行标记，用于残差回归
#  参数train表示作为训练数据的天数，超过该天数的数据标记为0
    tem = list()
    for i in range(ds.shape[0]): 
        if i in [0,1]:  # 前两项默认标识为0
            tem.append(0)
        elif i not in range(train):  #  为非训练数据，记0
            tem.append(0)
        #  骤升骤降的判断标准设定为相对两天前，今天的增/减量超过两天前的40%
        elif (ds['y'][i] - ds['y'][i-2]) > 0.4 * ds['y'][i-2]:  # 骤升，记1
            tem.append(1) 
        elif -(ds['y'][i] - ds['y'][i-2]) > 0.4 * ds['y'][i-2]:  # 骤降，记-1
            tem.append(-1)
        else:  # 正常情况，记0
            tem.append(0)
    return tem

def Fit(data0,
        train=24,periods=6,  #  训练天数以及预测天数
        m=None,holidays=None,
        predict=True,Draw=True,  #  标记是否进行预测和绘图
        length=8,width=4,  # 绘图时代图片尺寸
        Reg=False,  #  是否使用残差回归项
       ):
#     data = pd.read_csv(f'{INPUT_DIR}/Data.csv')
    if holidays==None:   # 若为指定holiday数据，则生成
    #  这里设定为每年元旦、3月开学、6月中、8月中、9月开学、双十一，影响周期为后15天
        holidays = pd.DataFrame({
   
          'holiday': 'Promotion',
          'ds': pd.to_datetime(['2019-01-01','2019-03-01','2019-06-15',
                               '2019-08-15','2019-09-01','2019-11-11',
                               '2020-01-01','2020-03-01','2020-06-15',
                               '2020-08-15','2020-09-01','2020-11-11',
                               '2021-01-01','2021-03-01','2021-06-15',
                               '2021-08-15','2021-09-01','2021-11-11',]),
          'lower_window': 0,
          'upper_window': 15,
        })
    #  未指定训练模型则创建
    if m==None:
        m = Prophet(growth='linear',
                    holidays=holidays,holidays_prior_scale=10,
                    n_changepoints=5,changepoint_range=0.9,changepoint_prior_scale=0.5,
                    daily_seasonality=False,weekly_seasonality=False,yearly_seasonality=False)
        #  趋势项使用分段线性模型，5个变点，均匀分布在前90%数据上，其初始影响强度为0.5
        m.add_seasonality(name='yearly', period=365.25, fourier_order=3,prior_scale=10)
        #  加入年周期项，周期为365.25，傅立叶展开到第3项，初始影响强度为默认值10