python数据分析-时间序列分析

本文已参与「新人创造礼」活动，一同开启创造之路

时刻序列剖析是依据随机过程理论和数理核算学办法

• 每日的均匀气温
• 每天的销售额
• 每月的降水量

时刻序列剖析首要经过statsmodel库的tsa模块完结

• 依据时刻序列的散点图，自相关函数和偏自相关函数图辨认序列是否平稳的非随机序列，假如对错随机序列，调查其平稳性
• 对非平稳的时刻序列数据采用差分进行滑润处理
• 依据辨认出来的特征树立相应的时刻序列模型
• 参数估计，查验是否具有核算含义
• 假定查验，判别模型的残差序列是否为白噪声序列
• 利用已经过查验的模型进行猜测

6.1 时刻序列的相关查验

白噪声查验

假如为白噪声数据(即独立分布的随机数据)，阐明其没有任何有用的信息

## 输出高清图画
%config InlineBackend.figure_format = 'retina'
%matplotlib inline
## 图画显现中文的问题
import matplotlib
matplotlib.rcParams['axes.unicode_minus']=False
import seaborn as sns 
sns.set(font= "Kaiti",style="ticks",font_scale=1.4)

## 导入会运用到的相关库
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.stattools import *
import statsmodels.api as sm
import statsmodels.formula.api as smf
from statsmodels.tsa.api import SimpleExpSmoothing,Holt,ExponentialSmoothing,AR,ARIMA,ARMA
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
import pmdarima as pm
from sklearn.metrics import mean_absolute_error
import pyflux as pf
from fbprophet import Prophet
## 疏忽提醒
import warnings
warnings.filterwarnings("ignore")

## 读取时刻序列数据,该数据包括:X1为飞机乘客数据,X2为一组随机数据
df = pd.read_csv("data/chap6/timeserise.csv")
## 查看数据的改变趋势
df.plot(kind = "line",figsize = (10,6))
plt.grid()
plt.title("时序数据")
plt.show()

## 白噪声查验Ljung-Box查验
## 该查验用来查看序列是否为随机序列，假如是随机序列，那它们的值之间没有任何关系
## 运用LB查验来查验序列是否为白噪声，原假定为在推迟期数内序列之间相互独立。
lags = [4,8,16,32]
LB = sm.stats.diagnostic.acorr_ljungbox(df["X1"],lags = lags,return_df = True)
print("序列X1的查验成果:\n",LB)
LB = sm.stats.diagnostic.acorr_ljungbox(df["X2"],lags = lags,return_df = True)
print("序列X2的查验成果:\n",LB)
## 假如P值小于0.05，阐明序列之间不独立，不是白噪声
'''
序列X1的查验成果:
         lb_stat      lb_pvalue
4    427.738684   2.817731e-91
8    709.484498  6.496271e-148
16  1289.037076  1.137910e-264
32  1792.523003   0.000000e+00
序列X2的查验成果:
       lb_stat  lb_pvalue
4    1.822771   0.768314
8    8.452830   0.390531
16  15.508599   0.487750
32  28.717743   0.633459
'''

在推迟阶数为[4,6,16,32]的状况下，序列X1的LB查验P值均小于0.05，即该数据不是随机的。有规则可循，有剖析价值，而序列X2的LB查验P值均大于0.05，该数据为白噪声，没有剖析价值

平稳性查验

时刻序列是否平稳，对挑选猜测的数学模型非常要害

假如数据是平稳的，能够运用自回归均匀移动模型（ARMA）

假如数据是不平稳的，能够运用差分移动自回归均匀移动模型（ARIMA）

\

## 序列的单位根查验，即查验序列的平稳性
dftest = adfuller(df["X2"],autolag='BIC')
dfoutput = pd.Series(dftest[0:4], index=['adf','p-value','usedlag','Number of Observations Used'])
print("X2单位根查验成果:\n",dfoutput)
dftest = adfuller(df["X1"],autolag='BIC')
dfoutput = pd.Series(dftest[0:4], index=['adf','p-value','usedlag','Number of Observations Used'])
print("X1单位根查验成果:\n",dfoutput)
## 对X1进行一阶差分后的序列进行查验
X1diff = df["X1"].diff().dropna()
dftest = adfuller(X1diff,autolag='BIC')
dfoutput = pd.Series(dftest[0:4], index=['adf','p-value','usedlag','Number of Observations Used'])
print("X1一阶差分单位根查验成果:\n",dfoutput)
## 一阶差分后 P值大于0.05, 小于0.1，能够以为其是平稳的
'''
X2单位根查验成果:
 adf                           -1.124298e+01
p-value                        1.788000e-20
usedlag                        0.000000e+00
Number of Observations Used    1.430000e+02
dtype: float64
X1单位根查验成果:
 adf                              0.815369
p-value                          0.991880
usedlag                         13.000000
Number of Observations Used    130.000000
dtype: float64
X1一阶差分单位根查验成果:
 adf                             -2.829267
p-value                          0.054213
usedlag                         12.000000
Number of Observations Used    130.000000
dtype: float64
'''

序列X2的查验P值小于0.05，阐明X2是一个平稳时刻序列（该序列是白噪声，白噪声序列是平稳序列）

序列X1的查验P值远大于0.05，阐明不平稳，而其一阶差分后的成果，P值大于0.05，但小于0.1，能够以为平稳

\

\

针对数据的平稳性查验，还能够运用KPSS查验，其原假定该序列是平稳的，该查验能够用kpss()函数完结

## KPSS查验的原假定为：序列x是平稳的。
## 对序列X2运用KPSS查验平稳性
dfkpss = kpss(df["X2"])
dfoutput = pd.Series(dfkpss[0:3], index=["kpss_stat"," p-value"," usedlag"])
print("X2 KPSS查验成果:\n",dfoutput)
## 承受序列平稳的原假定
## 对序列X1运用KPSS查验平稳性
dfkpss = kpss(df["X1"])
dfoutput = pd.Series(dfkpss[0:3], index=["kpss_stat"," p-value"," usedlag"])
print("X1 KPSS查验成果:\n",dfoutput)
## 回绝序列平稳的原假定
## 对序列X1运用KPSS查验平稳性
dfkpss = kpss(X1diff)
dfoutput = pd.Series(dfkpss[0:3], index=["kpss_stat"," p-value"," usedlag"])
print("X1一阶差分KPSS查验成果:\n",dfoutput)
## 承受序列平稳的原假定
'''
X2 KPSS查验成果:
 kpss_stat     0.087559
 p-value      0.100000
 usedlag     14.000000
dtype: float64
X1 KPSS查验成果:
 kpss_stat     1.052175
 p-value      0.010000
 usedlag     14.000000
dtype: float64
X1一阶差分KPSS查验成果:
 kpss_stat     0.05301
 p-value      0.10000
 usedlag     14.00000
dtype: float64
'''

ARIMA(p,d,q)模型

## 查验ARIMA模型的参数d
X1d = pm.arima.ndiffs(df["X1"], alpha=0.05, test="kpss", max_d=3)
print("运用KPSS办法对序列X1的参数d取值进行猜测,d = ",X1d)
X1diffd = pm.arima.ndiffs(X1diff, alpha=0.05, test="kpss", max_d=3)
print("运用KPSS办法对序列X1一阶差分后的参数d取值进行猜测,d = ",X1diffd)
X2d = pm.arima.ndiffs(df["X2"], alpha=0.05, test="kpss", max_d=3)
print("运用KPSS办法对序列X2的参数d取值进行猜测,d = ",X2d)
'''
运用KPSS办法对序列X1的参数d取值进行猜测,d =  1
运用KPSS办法对序列X1一阶差分后的参数d取值进行猜测,d =  0
运用KPSS办法对序列X1的参数d取值进行猜测,d =  0
'''

针对SARIMA模型，还有一个时节性平稳性参数D

## 查验SARIMA模型的参数时节阶数D
X1d = pm.arima.nsdiffs(df["X1"], 12, max_D=2)
print("对序列X1的时节阶数D取值进行猜测,D = ",X1d)
X1diffd = pm.arima.nsdiffs(X1diff, 12, max_D=2)
print("序列X1一阶差分后的时节阶数D取值进行猜测,D = ",X1diffd)
'''
对序列X1的时节阶数D取值进行猜测,D =  1
序列X1一阶差分后的时节阶数D取值进行猜测,D =  1
'''

自相关和偏相关剖析

## 对随机序列X2进行自相关和偏相关剖析可视化
fig = plt.figure(figsize=(16,5))
plt.subplot(1,3,1)
plt.plot(df["X2"],"r-")
plt.grid()
plt.title("X2序列动摇")
ax = fig.add_subplot(1,3,2)
plot_acf(df["X2"], lags=60,ax = ax)
plt.grid()
ax = fig.add_subplot(1,3,3)
plot_pacf(df["X2"], lags=60,ax = ax)
plt.grid()
plt.tight_layout()
plt.show()

在图画中滞后0表示自己和自己的相关性，恒等于1。不用于确认p和q。

## 对非随机序列X1进行自相关和偏相关剖析可视化
fig = plt.figure(figsize=(16,5))
plt.subplot(1,3,1)
plt.plot(df["X1"],"r-")
plt.grid()
plt.title("X1序列动摇")
ax = fig.add_subplot(1,3,2)
plot_acf(df["X1"], lags=60,ax = ax)
plt.grid()
ax = fig.add_subplot(1,3,3)
plot_pacf(df["X1"], lags=60,ax = ax)
plt.ylim([-1,1])
plt.grid()
plt.tight_layout()
plt.show()

## 对非随机序列X1一阶差分后的序列进行自相关和偏相关剖析可视化
fig = plt.figure(figsize=(16,5))
plt.subplot(1,3,1)
plt.plot(X1diff,"r-")
plt.grid()
plt.title("X1序列一阶差分后动摇")
ax = fig.add_subplot(1,3,2)
plot_acf(X1diff, lags=60,ax = ax)
plt.grid()
ax = fig.add_subplot(1,3,3)
plot_pacf(X1diff, lags=60,ax = ax)
plt.grid()
plt.tight_layout()
plt.show()

ARMA(p,q)中，自相关系数的滞后，对应着参数q；偏相关系数的滞后对应着参数p。

## 时刻序列的分化
## 经过调查序列X1，能够发现其既有上升的趋势，也有周期性的趋势，所以能够将该序列进行分化
## 运用乘法模型分化成果(一般适用于有增长趋势的序列)
X1decomp = pm.arima.decompose(df["X1"].values,"multiplicative", m=12)
## 可视化出分化的成果
ax = pm.utils.decomposed_plot(X1decomp,figure_kwargs = {"figsize": (10, 6)},
                              show=False)
ax[0].set_title("乘法模型分化成果")
plt.show()

## 运用加法模型分化成果(一般适用于平稳趋势的序列)
X1decomp = pm.arima.decompose(X1diff.values,"additive", m=12)
## 可视化出分化的成果
ax = pm.utils.decomposed_plot(X1decomp,figure_kwargs = {"figsize": (10, 6)},
                              show=False)
ax[0].set_title("加法模型分化成果")
plt.show()

6.2 移动均匀算法

## 数据预备
## 对序列X1进行切分，后边的24个数据用于测验集
train = pd.DataFrame(df["X1"][0:120])
test = pd.DataFrame(df["X1"][120:])
## 可视化切分后的数据
train["X1"].plot(figsize=(14,7), title= "乘客数量数据",label = "X1 train")
test["X1"].plot(label = "X1 test")
plt.legend()
plt.grid()
plt.show()

print(train.shape)
print(test.shape)
df["X1"].shape
'''
(120, 1)
(24, 1)
(144,)
'''

6.2.1 简略移动均匀法

## 简略移动均匀进行猜测
y_hat_avg = test.copy(deep = False)
y_hat_avg["moving_avg_forecast"] =  train["X1"].rolling(24).mean().iloc[-1]
## 可视化出猜测成果
plt.figure(figsize=(14,7))
train["X1"].plot(figsize=(14,7),label = "X1 train")
test["X1"].plot(label = "X1 test")
y_hat_avg["moving_avg_forecast"].plot(style="g--o", lw=2,
                                      label="移动均匀猜测")
plt.legend()
plt.grid()
plt.title("简略移动均匀猜测")
plt.show()

## 核算猜测成果和实在值的差错
print("猜测绝对值差错:",mean_absolute_error(test["X1"],y_hat_avg["moving_avg_forecast"]))
'''
猜测绝对值差错: 82.55208333333336
'''

6.2.2 简略指数滑润法

## 数据预备
y_hat_avg = test.copy(deep = False)
## 模型构建
model1 = SimpleExpSmoothing(train["X1"].values).fit(smoothing_level=0.15)
y_hat_avg["exp_smooth_forecast1"] = model1.forecast(len(test))
model2 = SimpleExpSmoothing(train["X1"].values).fit(smoothing_level=0.5)
y_hat_avg["exp_smooth_forecast2"] = model2.forecast(len(test))
## 可视化出猜测成果
plt.figure(figsize=(14,7))
train["X1"].plot(figsize=(14,7),label = "X1 train")
test["X1"].plot(label = "X1 test")
y_hat_avg["exp_smooth_forecast1"].plot(style="g--o", lw=2,
                                      label="smoothing_level=0.15")
y_hat_avg["exp_smooth_forecast2"].plot(style="g--s", lw=2,
                                      label="smoothing_level=0.5")
plt.legend()
plt.grid()
plt.title("简略指数滑润猜测")
plt.show()
## 核算猜测成果和实在值的差错
print("smoothing_level=0.15,猜测绝对值差错:",
      mean_absolute_error(test["X1"],y_hat_avg["exp_smooth_forecast1"]))
print("smoothing_level=0.5,猜测绝对值差错:",
      mean_absolute_error(test["X1"],y_hat_avg["exp_smooth_forecast2"]))

smoothing_level=0.15,猜测绝对值差错: 81.10115706423566

smoothing_level=0.5,猜测绝对值差错: 106.813228720506

6.2.3 霍尔特(Holt)线性趋势法

## 数据预备
y_hat_avg = test.copy(deep = False)
## 模型构建
model1 = Holt(train["X1"].values).fit(smoothing_level=0.1,
                                      smoothing_slope = 0.05)
y_hat_avg["holt_forecast1"] = model1.forecast(len(test))
model2 = Holt(train["X1"].values).fit(smoothing_level=0.1,
                                      smoothing_slope = 0.25)
y_hat_avg["holt_forecast2"] = model2.forecast(len(test))
## 可视化出猜测成果
plt.figure(figsize=(14,7))
train["X1"].plot(figsize=(14,7),label = "X1 train")
test["X1"].plot(label = "X1 test")
y_hat_avg["holt_forecast1"].plot(style="g--o", lw=2,
                                 label="Holt线性趋势法(1)")
y_hat_avg["holt_forecast2"].plot(style="g--s", lw=2,
                                 label="Holt线性趋势法(2)")
plt.legend()
plt.grid()
plt.title("Holt线性趋势法猜测")
plt.show()
## 核算猜测成果和实在值的差错
print("smoothing_slope = 0.05,猜测绝对值差错:",
      mean_absolute_error(test["X1"],y_hat_avg["holt_forecast1"]))
print("smoothing_slope = 0.25,猜测绝对值差错:",
      mean_absolute_error(test["X1"],y_hat_avg["holt_forecast2"]))

smoothing_slope = 0.05,猜测绝对值差错: 54.727467142360275

smoothing_slope = 0.25,猜测绝对值差错: 69.79052992788556

6.2.4 Holt-Winters时节性猜测模型

## 数据预备
y_hat_avg = test.copy(deep = False)
## 模型构建
model1 = ExponentialSmoothing(train["X1"].values,
                              seasonal_periods=12, # 周期性为12  
                              trend="add", seasonal="add").fit()
y_hat_avg["holt_winter_forecast1"] = model1.forecast(len(test))
## 可视化出猜测成果
plt.figure(figsize=(14,7))
train["X1"].plot(figsize=(14,7),label = "X1 train")
test["X1"].plot(label = "X1 test")
y_hat_avg["holt_winter_forecast1"].plot(style="g--o", lw=2,
                                 label="Holt-Winters")
plt.legend()
plt.grid()
plt.title("Holt-Winters时节性猜测模型")
plt.show()
## 核算猜测成果和实在值的差错
print("Holt-Winters时节性猜测模型,猜测绝对值差错:",
      mean_absolute_error(test["X1"],y_hat_avg["holt_winter_forecast1"]))

Holt-Winters时节性猜测模型,猜测绝对值差错: 30.06821059070873

6.3 ARIMA模型

留意针对乘客数据X1，运用AR模型或许ARMA模型进行猜测，并不对错常的合适，

这里是运用AR和ARMA模型进行猜测的目的首要是为了和更好的模型猜测成果进行对比

## 运用AR模型对乘客数据进行猜测 
## 经过前面序列的偏相关系数的可视化成果，运用AR(2)模型可对序列进行建模
## 数据预备
y_hat = test.copy(deep = False)
## 模型构建
ar_model = ARMA(train["X1"].values,order = (2,0)).fit()
## 输出拟合模型的成果
print(ar_model.summary())
## AIC=1141.989；BIC= 1153.138；两个系数是明显的

## 查看模型的拟合残差分布
fig = plt.figure(figsize=(12,5))
ax = fig.add_subplot(1,2,1)
plt.plot(ar_model.resid)
plt.title("AR(2)残差曲线")
## 查看残差是否符合正太分布
ax = fig.add_subplot(1,2,2)
sm.qqplot(ar_model.resid, line='q', ax=ax)
plt.title("AR(2)残差Q-Q图")
plt.tight_layout()
plt.show()

## 猜测未来24个数据，并输出95%置信区间
pre, se, conf = ar_model.forecast(24, alpha=0.05)  
## 收拾数据
y_hat["ar2_pre"] = pre
y_hat["ar2_pre_lower"] = conf[:,0]
y_hat["ar2_pre_upper"] = conf[:,1]
## 可视化出猜测成果
plt.figure(figsize=(14,7))
train["X1"].plot(figsize=(14,7),label = "X1 train")
test["X1"].plot(label = "X1 test")
y_hat["ar2_pre"].plot(style="g--o", lw=2,label="AR(2)")
## 可视化出置信区间
plt.fill_between(y_hat.index, y_hat["ar2_pre_lower"], 
                 y_hat["ar2_pre_upper"],color='k',alpha=.15,
                 label = "95%置信区间")
plt.legend()
plt.grid()
plt.title("AR(2)模型")
plt.show()
# 核算猜测成果和实在值的差错
print("AR(2)模型猜测的绝对值差错:",
      mean_absolute_error(test["X1"],y_hat["ar2_pre"]))

AR(2)模型猜测的绝对值差错: 165.79608244918572

能够发现运用AR(2)的猜测作用并不好

ARMA模型

## 测验运用ARMA模型进行猜测
## 依据前面的自相关系数和偏相关系数,为了降低模型的杂乱读,能够运用ARMA(2,1)
## 数据预备
y_hat = test.copy(deep = False)
## 模型构建
arma_model = ARMA(train["X1"].values,order = (2,1)).fit()
## 输出拟合模型的成果
print(arma_model.summary())
## AIC=1141.989；BIC= 1153.138；两个系数是明显的

## 查看模型的拟合残差分布
fig = plt.figure(figsize=(12,5))
ax = fig.add_subplot(1,2,1)
plt.plot(arma_model.resid)
plt.title("ARMA(2,1)残差曲线")
## 查看残差是否符合正太分布
ax = fig.add_subplot(1,2,2)
sm.qqplot(arma_model.resid, line='q', ax=ax)
plt.title("ARMA(2,1)残差Q-Q图")
plt.tight_layout()
plt.show()

## 猜测未来24个数据，并输出95%置信区间
pre, se, conf = arma_model.forecast(24, alpha=0.05)  
## 收拾数据
y_hat["arma_pre"] = pre
y_hat["arma_pre_lower"] = conf[:,0]
y_hat["arma_pre_upper"] = conf[:,1]
## 可视化出猜测成果
plt.figure(figsize=(14,7))
train["X1"].plot(figsize=(14,7),label = "X1 train")
test["X1"].plot(label = "X1 test")
y_hat["arma_pre"].plot(style="g--o", lw=2,label="ARMA(2,1)")
## 可视化出置信区间
plt.fill_between(y_hat.index, y_hat["arma_pre_lower"], 
                 y_hat["arma_pre_upper"],color='k',alpha=.15,
                 label = "95%置信区间")
plt.legend()
plt.grid()
plt.title("ARMA(2,1)模型")
plt.show()
# 核算猜测成果和实在值的差错
print("ARMA模型猜测的绝对值差错:",
      mean_absolute_error(test["X1"],y_hat["arma_pre"]))

ARMA模型猜测的绝对值差错: 147.26531763335154

针对ARMA模型主动挑选合适的参数

## 主动搜索合适的参数
model = pm.auto_arima(train["X1"].values,
                      start_p=1, start_q=1, # p,q的开始值
                      max_p=12, max_q=12, # 最大的p和q
                      d = 0,            # 寻觅ARMA模型参数
                      m=1,              # 序列的周期
                      seasonal=False,   # 没有时节性趋势
                      trace=True,error_action='ignore',  
                      suppress_warnings=True, stepwise=True)
print(model.summary())

## 运用ARMA(3,2)对测验集进行猜测
pre, conf = model.predict(n_periods=24, alpha=0.05,
                          return_conf_int=True)
## 可视化ARMA(3,2)的猜测成果，收拾数据
y_hat["arma_pre"] = pre
y_hat["arma_pre_lower"] = conf[:,0]
y_hat["arma_pre_upper"] = conf[:,1]
## 可视化出猜测成果
plt.figure(figsize=(14,7))
train["X1"].plot(figsize=(14,7),label = "X1 train")
test["X1"].plot(label = "X1 test")
y_hat["arma_pre"].plot(style="g--o", lw=2,label="ARMA(3,2)")
## 可视化出置信区间
plt.fill_between(y_hat.index, y_hat["arma_pre_lower"], 
                 y_hat["arma_pre_upper"],color='k',alpha=.15,
                 label = "95%置信区间")
plt.legend()
plt.grid()
plt.title("ARMA(3,2)模型")
plt.show()
# 核算猜测成果和实在值的差错
print("ARMA模型猜测的绝对值差错:",
      mean_absolute_error(test["X1"],y_hat["arma_pre"]))

ARMA模型猜测的绝对值差错: 158.11464180972925

能够发现运用主动ARMA(3,2)模型的作用并没有ARMA(2,1)的猜测作用好

\

\

\

6.7 时序数据的异常值检测

能够将忽然增大或忽然减小的数据无规则看作异常值

## 运用prophet检测时刻序列是否有异常值
## 从1991年2月到2005年5月，每周供给美国成品汽车汽油产品的时刻序列（每天数千桶）
## 数据预备
data = pm.datasets.load_gasoline()
datadf = pd.DataFrame({"y":data})
datadf["ds"] = pd.date_range(start="1991-2",periods=len(data),freq="W")
## 可视化时刻序列的改变状况
datadf.plot(x = "ds",y = "y",style = "b-o",figsize=(14,7))
plt.grid()
plt.title("时刻序列数据的动摇状况")
plt.show()

\

## 对该数据树立一个时刻序列模型
np.random.seed(1234)  ## 设置随机数种子
model = Prophet(growth="linear",daily_seasonality = False,
                weekly_seasonality=False,
                seasonality_mode = 'multiplicative', 
                interval_width = 0.95,   ## 获取95%的置信区间
                )
model = model.fit(datadf)     # 运用数据拟合模型
forecast = model.predict(datadf)  # 运用模型对数据进行猜测
forecast["y"] = datadf["y"].reset_index(drop = True)
forecast[["ds","y","yhat","yhat_lower","yhat_upper"]].head()

ds	y	yhat	yhat_lower	yhat_upper
0	1991-02-03	6621.0	6767.051491	6294.125979	7303.352309
1	1991-02-10	6433.0	6794.736479	6299.430616	7305.414252
2	1991-02-17	6582.0	6855.096282	6352.579489	7379.717614
3	1991-02-24	7224.0	6936.976642	6415.157617	7445.523000
4	1991-03-03	6875.0	6990.511503	6489.781400	7488.240435

## 依据模型猜测值的置信区间"yhat_lower"和"yhat_upper"判别样本是否为异常值
def outlier_detection(forecast):
    index = np.where((forecast["y"] <= forecast["yhat_lower"])|
                     (forecast["y"] >= forecast["yhat_upper"]),True,False)
    return index
outlier_index = outlier_detection(forecast)
outlier_df = datadf[outlier_index]
print("异常值的数量为:",np.sum(outlier_index))
'''
异常值的数量为: 38
'''

\

## 可视化异常值的成果
fig, ax = plt.subplots()
## 可视化猜测值
forecast.plot(x = "ds",y = "yhat",style = "b-",figsize=(14,7),
              label = "猜测值",ax=ax)
## 可视化出置信区间
ax.fill_between(forecast["ds"].values, forecast["yhat_lower"], 
                forecast["yhat_upper"],color='b',alpha=.2,
                label = "95%置信区间")
forecast.plot(kind = "scatter",x = "ds",y = "y",c = "k",
              s = 20,label = "原始数据",ax = ax)
## 可视化出异常值的点
outlier_df.plot(x = "ds",y = "y",style = "rs",ax = ax,
                label = "异常值")
plt.legend(loc = 2)
plt.grid()
plt.title("时刻序列异常值检测成果")
plt.show()

异常值大部分都在置信区间外