酷应用

实例：使用LSTMs和Prophet进行时间序列预测你的电子邮箱负载（附代码）

百家作者：数据分析 2019-01-10 01:37:03

作者：Maximilian Strauß；翻译：笪洁琼；校对：丁楠雅

本文共3400字，建议阅读10分钟。
本文通过基线模型、LSTMs和Facebook的Prophet模型来预测每天的电子邮箱负荷，并详细解析了生成训练数据集的过程以及相应代码。

时间序列预测为数据科学算法提供了一个极好的训练场。

毕竟，如果一个人能够预测未来，那该有多牛逼啊！通常用来演示预测算法的典型数据集是股票图表、销售和气象数据。在这里，我们将尝试一些与每个用户更相关的东西，即你将收到的电子邮件的数量。根据一份电子邮件统计报告，2014年上半年，上班族平均每天收到85封电子邮件(http:/www.Radati.com/wp/wp-content/upploads/2014/01/Email-Statistics-Report-2014-2018-Executive-Summary.pdf）。

我们想尝试根据历史收件数据做出准确的预测。为此，我们将探索运用LSTMs和Facebook的Prophet。这里的目标是如何为不同的算法准备数据，并提供定性的概述，而不是精细化。预测结果将根据历史收到的电子邮件数量以及所准备的训练数据而有很大的不同。

收集数据

在IBM imapclient sql包的帮助下，我们从使用自己的收件箱创建数据集开始。

关于Automatetheboringstuff.com/Chapter16/，可以在Automatetheboringstuff.com一章地址：https：/Automatetheboringstuff.com/中找到很好的介绍。

我们将加载过去三年(从2016年1月1日开始)的所有电子邮件，并获取主题和日期。我们会使用Pandas将其转换为一个数据文件(Dataframe)。

Code:
import imapclient
import pandas as pd
import getpass

youremail = input()
yourpassword = getpass.getpass()

# Replace 'imap.gmail.com' with provider of choice
imapObj = imapclient.IMAPClient("imap.gmail.com", ssl=True)
imapObj.login(youremail, yourpassword)
imapObj.select_folder("INBOX", readonly=True)

UIDs = imapObj.search('(SINCE "01-Jan-2016")')

mails = []
for msgid, data in imapObj.fetch(UIDs, ["ENVELOPE"]).items():
    envelope = data[b"ENVELOPE"]
    date = envelope.date
    if envelope.subject is not None:
        subject = envelope.subject.decode()
    else:
        subject = None
    mails.append((subject, date))

mail_df = pd.DataFrame(mails)
mail_df.columns = ["Subject", "Date"]
mail_df["Date"] = pd.to_datetime(mail_df["Date"])
mail_df = mail_df.set_index("Date")

print("A total of {} e-mails loaded.".format(len(mail_df)))

邮箱负载：总计12738封邮件

我们现在得到了12738封用于训练的电子邮件。请注意，上面的代码是针对一个装满电子邮件的收件箱。如果你已将电子邮件放在不同的文件夹中，请相应地调整代码。包括我在内的一些人会立即删除那些不重要的电子邮件。然后，模型输出的是重要邮件的数量，而不是实际收到的电子邮件数量。还要注意的是，一些电子邮件提供商(如Google)会阻止这个连接，因为他们不允许“不太安全”的应用程序连接到他们的服务。你可以在其设置中启用此功能。原则上，你还可以检验（签出）本地邮箱。检验（签出）统一邮箱程序包可能是个好的开始。

数据探索

现在，让我们首先做一些可视化的数据探索和图表，每小时和每天的电子邮件数量。对于这一点，我们将使用pandas及重新采样的groubpy 函数（聚合分组运算）。通过使用sum()和count()参数，我们可以对每个时间间隔进行标准化。另外，我们还使用了seaborn函数来实现可视化。

Code:
import calendar
import seaborn as sns
import matplotlib.pyplot as plt

weekdays = [calendar.day_name[i] for i in range(7)]

# E-Mails per Hour
per_hour = pd.DataFrame(mail_df["Subject"].resample("h").count())
per_hour_day = (
    per_hour.groupby([per_hour.index.hour]).sum()
    / per_hour.groupby([per_hour.index.hour]).count()
)
per_hour_day.reset_index(inplace=True)
per_hour_day.columns = ["Hour", "Count"]

# E-Mails per day
per_day = pd.DataFrame(mail_df["Subject"].resample("d").count())
per_day_week = (
    per_day.groupby([per_day.index.weekday]).sum()
    / per_day.groupby([per_day.index.weekday]).count()
)
per_day_week.reset_index(inplace=True)
per_day_week.columns = ["Weekday", "Count"]
per_day_week["Weekday"] = weekdays

def return_cmap(data):
    # Function to create a colormap
    v = data["Count"].values
    colors = plt.cm.RdBu_r((v - v.min()) / (v.max() - v.min()))
    return colors

plt.figure(figsize=(12, 10), dpi=600)
plt.subplot(2, 1, 1)
cmap = return_cmap(per_hour_day)
sns.barplot(x="Hour", y="Count", data=per_hour_day, palette=cmap)
plt.title("Emails per hour")

plt.subplot(2, 1, 2)
cmap = return_cmap(per_day_week)
sns.barplot(x="Weekday", y="Count", data=per_day_week, palette=cmap)
plt.title("Emails per weekday")

plt.show()

print(
    "Average number of emails per day: {:.2f}".format(
        per_hour_day.sum()["Count"]
    )
)

图1：探索数据（图）

从数据中我们可以看到一个明显的规律，即分发遵循一个典型模式，从星期一至星期五早9-晚5的时间表，其中上午10点收到的邮件最多。由于每小时电子邮件的最大值仅略高于1，因此对每小时进行预测是没有意义的。我们将尝试预测每天的电子邮件工作量。

初步考虑：预测不可预测的

在深入研究我们的模型之前，有必要考虑一下这个问题的相关假设。接收电子邮件的时间在某种程度上是随机的，因此无法预测。在大多数情况下，发送方互相不认识，因此可以假设它们在统计上是独立的。另外，我们可以将收件近似为泊松过程（Poissonian），它的标准偏差等于平均值。由于分布是随机的，12封电子邮件的期望的RMSE值至少为3.5(sqrt(11.95)~3.461)。

如果你想进一步阅读，我建议阅读有关DB2Erlang发行版

(https：/en.wikipea.org/wiki/Erlang_Distributions)的内容，研究该发行版是为了检验在某一时期内的电话呼出的数量。(https：/en.wikipea.org/wiki/erlang_Distributions)

计算基准线

为了创建基线模型，我们可以使用历史数据中的查询表。对于任意的某一天，计算在一天之前收到的电子邮件数量。为了对模型进行基准测试，我使用一个移动窗口来创建查询表，并计算到第二天的预测差异值。

Code:
from sklearn.metrics import mean_squared_error
import numpy as np

data = per_day.copy()
test_split = int(len(data) * 0.8)

pred_base = []
for i in range(len(data) - test_split):
    train_data = data[i : test_split + i]
    test_data = data.iloc[test_split + i]
    train_data_week = (
        train_data.groupby([train_data.index.weekday]).sum()
        / train_data.groupby([train_data.index.weekday]).count()
    )
    baseline_prediction = train_data_week.loc[test_data.name.weekday()]
    pred_base.append(baseline_prediction.values)

test_data = data[test_split:]

mse_baseline = mean_squared_error(test_data.values, pred_base)

print("RMSE for BASELINE {:.2f}".format(np.sqrt(mse_baseline)))

RMSE 基线 7.39

基线模型得到的RMSE为7.39 ，相对于3.5的预期值，这个结果还不错。

使用LSTM进行预测

作为进阶模型，我们将使用长短时记忆(LSTM)神经网络。在这里可以找到对LSTM的很好的介绍(https：/machinelearningmaster ery.com/time-Series-prediction-lstm-rrurn-neuro-network-python-keras/)。

在这里，我们将窗口设置为6天，并让模型预测第7天。

Code:
import numpy as np
from keras.models import Sequential
from keras.layers import Dense, Dropout, Activation
from keras.layers import LSTM
from sklearn.preprocessing import MinMaxScaler

def get_window_data(data, window):
    # Get window data and scale
    scaler = MinMaxScaler(feature_range=(0, 1))
    data = scaler.fit_transform(data.reshape(-1, 1))

    X = []
    y = []

    for i in range(len(data) - window - 1):
        X.append(data[i : i + window])
        y.append(data[i + window + 1])

    X = np.asarray(X)
    y = np.asarray(y)
    return X, y, scaler

window_size = 6
X, y, scaler = get_window_data(per_day["Subject"].values, window_size)

X_train = X[:test_split]
X_test = X[test_split:]

y_train = y[:test_split]
y_test = y[test_split:]

model = Sequential()
model.add(LSTM(50, input_shape=(window_size, 1)))
model.add(Dropout(0.2))
model.add(Dense(1))
model.add(Activation("linear"))
model.compile(loss="mse", optimizer="adam")

history = model.fit(
    X_train,
    y_train,
    epochs=20,
    batch_size=1,
    validation_data=(X_test, y_test),
    verbose=2,
    shuffle=False,
)

# plot history

plt.figure(figsize=(6, 5), dpi=600)
plt.plot(history.history["loss"], 'darkred', label="Train")
plt.plot(history.history["val_loss"], 'darkblue', label="Test")
plt.title("Loss over epoch")
plt.xlabel('Epoch')
plt.ylabel('Loss')
plt.legend()
plt.show()

mse_lstm = mean_squared_error(
    scaler.inverse_transform(y_test),
    scaler.inverse_transform(model.predict(X_test)),
)
print("RMSE for LSTM {:.2f}".format(np.sqrt(mse_lstm)))

LSTM为7.26RMSE

图2：训练LSTM

通过观察损失函数，我们可以看到LSTM网络通过迭代学习更好地预测了未来。它的RMSE比基线模型要好一些，但也并没有好多少。这一结果表明，LSTM网络能够学习周末和工作日的结构，周末收到的电子邮件更少。另外值得一提的是，我们不能直接将基线的RMSE与LSTM进行比较。因为LSTM实际上只用到了80%的训练数据。我们也可以使用全部数据来训练LSTM，然而从计算上来说，这样做的代价要高得多。

使用Prophet进行预测

接下来，我们将使用facebook的Prophet库(https://github.com/facebook/prophet)。它是一个加性模型，我们可以用年、周和日的季节性来拟合非线性趋势。同样地，我们将把数据分成训练集和测试集，并计算RMSE值。

Code:
from fbprophet import Prophet
from tqdm import tqdm

prophet_data = data.reset_index()
prophet_data["ds"] = prophet_data["Date"]
prophet_data["y"] = prophet_data["Subject"]

pred = []
for i in tqdm(range(len(data) - test_split)):

    data_to_fit = prophet_data[: (test_split + i)]
    prophet_model = Prophet(interval_width=0.95)

    prophet_model.fit(data_to_fit)

    prophet_forecast = prophet_model.make_future_dataframe(periods=1, freq="d")
    prophet_forecast = prophet_model.predict(prophet_forecast)

    pred.append(prophet_forecast["yhat"].iloc[-1])

mse_prophet = mean_squared_error(test_data.values, pred)

print("RMSE for PROPHET {:.2f}".format(np.sqrt(mse_prophet)))

PROPHET为6.96RMSE

从RMSE来看，Prophet模型实现了最佳性能。现在让我们研究一下这是为什么。为此，我们绘制模型的各个构成部分，以便更好地理解模型所做的工作。这只需要使用性能指标函数就可以输出结果。

（输出）

from fbprophet.diagnostics import performance_metricsprophet_model.plot_components(prophet_forecast)

图3：Prophet模型的构成

通过检查Prophet模型的构成，我们可以看到它识别出了数据中的关键趋势。总趋势表明邮件数量在整体上不断增加。每周的季节性准确地描绘了工作日/周末的时间波动。每年的季节性显示主要节日，即新年电子邮件很少，但在圣诞节前有所增加，而低点在9月份。

总结：预测邮件数量。

最后，我们使用Prophet模型作为我们的预测工具。为此，我们再次登录到IMAP服务器，并使用自2016年1月1日以来的所有历史数据来训练我们的预测模型。训练完毕后，我们绘制了前一周的历史数据和下一周的预测情况。

Code:
from datetime import datetime, timedelta
import matplotlib.pyplot as plt
from fbprophet import Prophet
import imapclient
import pandas as pd
import getpass

youremail = input()
yourpassword = getpass.getpass()

imapObj = imapclient.IMAPClient("imap.gmail.com", ssl=True)
imapObj.login(youremail, yourpassword)
imapObj.select_folder("INBOX", readonly=True)

UIDs = imapObj.search('(SINCE "01-Jan-2016")')

mails = []
for msgid, data in imapObj.fetch(UIDs, ["ENVELOPE"]).items():
envelope = data[b"ENVELOPE"]
date = envelope.date

    if envelope.subject is not None:
        subject = envelope.subject.decode()
    else:
        subject = None
    mails.append((subject, date))

mail_df = pd.DataFrame(mails)
mail_df.columns = ["Subject", "Date"]
mail_df["Date"] = pd.to_datetime(mail_df["Date"])
mail_df = mail_df.set_index("Date")

data = pd.DataFrame(mail_df["Subject"].resample("d").count())

prophet_model = Prophet(interval_width=0.95)

prophet_data = data.reset_index()
prophet_data["ds"] = prophet_data["Date"]
prophet_data["y"] = prophet_data["Subject"]

prophet_model.fit(prophet_data)

prophet_forecast = prophet_model.make_future_dataframe(periods=7, freq="d")
prophet_forecast = prophet_model.predict(prophet_forecast)

fig1 = prophet_model.plot(prophet_forecast)

datenow = datetime.now()
dateend = datenow + timedelta(days=7)
datestart = dateend - timedelta(days=14)

plt.xlim([datestart, dateend])
plt.title("Email forecast", fontsize=20)
plt.xlabel("Day", fontsize=20)
plt.ylabel("Emails expected", fontsize=20)
plt.axvline(datenow, color="k", linestyle=":")
plt.show()