投稿/爆料
厂商入驻

机器学习开放课程:九、基于Python分析真实手游时序数据

论智|2018-09-02 21:17

【编者按】Zeptolab数据科学家Dmitriy Sergeev介绍了分析和预测时序数据的主要方法。

题图

大家好!

这次的开放机器学习课程的内容是时序数据。

我们将查看如何使用Python处理时序数据,哪些方法和模型可以用来预测;什么是双指数平滑和三指数平滑;如果平稳(stationarity)不是你的菜,该怎么办;如何创建SARIMA并且活下来;如何使用XGBoost做出预测。所有这些都将应用于(严酷的)真实世界例子。

概览

  1. 导言

    • 基本定义
    • 量化指标
  2. 平移、平滑、评估

    • 滚动窗口估计
    • 指数平滑,Holt-Winters模型
    • 时序交叉验证,参数选择
  3. 计量经济学方法

    • 平稳,单位根
    • 摆脱不平稳性
    • SARIMA的直觉和建模
  4. 时序数据上的线性(以及不那么线性)模型

    • 特征提取
    • 线性模型,特征重要性
    • 正则化,特征选取
    • XGBoost
  5. 相关资源

导言

在我的工作中,我几乎每天都会碰到和时序有关的任务。最频繁的问题是——明天/下一周/下个月/等等,我们的指标将是什么样——有多少玩家会安装应用,他们的在线时长会是多少,他们会进行多少次操作,取决于预测所需的质量,预测周期的长度,以及时刻,我们需要选择特征,调整参数,以取得所需结果。

基本定义

时序的简单定义

时序——一系列以时间顺序为索引(或列出、绘出)的数据点。

因此,数据以相对确定的时刻组织。所以,和随机样本相比,可能包含我们将尝试提取的额外信息。

让我们导入一些库。首先我们需要statsmodels库,它包含了一大堆统计学建模函数,包括时序。对不得不迁移到Python的R粉来说,绝对会感到statsmodels很熟悉,因为它支持类似Wage ~ Age + Education这样的模型定义。

import numpy as np                               # 向量和矩阵
import pandas as pd                              # 表格和数据处理
import matplotlib.pyplot as plt                  # 绘图
import seaborn as sns                            # 更多绘图功能

from dateutil.relativedelta import relativedelta # 处理不同格式的时间日期
from scipy.optimize import minimize              # 最小化函数

import statsmodels.formula.api as smf            # 统计学和计量经济学
import statsmodels.tsa.api as smt
import statsmodels.api as sm
import scipy.stats as scs

from itertools import product                    # 一些有用的函数
from tqdm import tqdm_notebook

import warnings                                  # 勿扰模式
warnings.filterwarnings('ignore')

%matplotlib inline

作为例子,让我们使用一些真实的手游数据,玩家每小时观看的广告,以及每天游戏币的花费:

ads = pd.read_csv('../../data/ads.csv', index_col=['Time'], parse_dates=['Time'])
currency = pd.read_csv('../../data/currency.csv', index_col=['Time'], parse_dates=['Time'])

plt.figure(figsize=(15, 7))
plt.plot(ads.Ads)
plt.title('Ads watched (hourly data)')
plt.grid(True)
plt.show()

plt.figure(figsize=(15, 7))
plt.plot(currency.GEMS_GEMS_SPENT)
plt.title('In-game currency spent (daily data)')
plt.grid(True)
plt.show()

广告观看

游戏币花费

预测质量指标

在实际开始预测之前,先让我们理解下如何衡量预测的质量,查看下最常见、使用最广泛的测度:

  • R2,决定系数(在经济学中,可以理解为模型能够解释的方差比例),(-inf, 1] sklearn.metrics.r2_score
  • 平均绝对误差(Mean Absolute Error),这是一个易于解释的测度,因为它的计量单位和初始序列相同,[0, +inf) sklearn.metrics.mean_absolute_error
  • 中位绝对误差(Median Absolute Error),同样是一个易于解释的测度,对离群值的鲁棒性很好,[0, +inf) sklearn.metrics.median_absolute_error
  • 均方误差(Mean Squared Error),最常用的测度,给较大的错误更高的惩罚,[0, +inf) sklearn.metrics.mean_squared_error
  • 均方对数误差(Mean Squared Logarithmic Error),和MSE差不多,只不过先对序列取对数,因此能够照顾到较小的错误,通常用于具有指数趋势的数据,[0, +inf) sklearn.metrics.mean_squared_log_error
  • 平均绝对百分误差(Mean Absolute Percentage Error),类似MAE不过基于百分比——当你需要向管理层解释模型的质量时很方便——[0, +inf),sklearn中没有实现。
# 引入上面提到的所有测度
from sklearn.metrics import r2_score, median_absolute_error, mean_absolute_error
from sklearn.metrics import median_absolute_error, mean_squared_error, mean_squared_log_error

# 自行实现sklearn没有提供的平均绝对百分误差很容易
def mean_absolute_percentage_error(y_true, y_pred): 
    return np.mean(np.abs((y_true - y_pred) / y_true)) * 100

棒极了,现在我们知道如何测量预测的质量了,可以使用哪些测度,以及如何向老板翻译结果。剩下的就是创建模型了。

平移、平滑、评估

让我们从一个朴素的假设开始——“明天会和今天一样”,但是我们并不使用类似y^t=y(t-1)这样的模型(这其实是一个适用于任意时序预测问题的很好的基线,有时任何模型都无法战胜这一模型),相反,我们将假定变量未来的值取决于前n个值的平均,所以我们将使用的是移动平均(moving average)

移动平均公式

def moving_average(series, n):
    """
        计算前n项观测的平均数
    """
    return np.average(series[-n:])

# 根据前24小时的数据预测
moving_average(ads, 24)

结果:116805.0

不幸的是这样我们无法做出长期预测——为了预测下一步的数据我们需要实际观测的之前的数据。不过移动平均还有一种用途——平滑原时序以显示趋势。pandas提供了实现DataFrame.rolling(window).mean()。窗口越宽,趋势就越平滑。遇到噪声很大的数据时(财经数据十分常见),这一过程有助于侦测常见模式。

def plotMovingAverage(series, window, plot_intervals=False, scale=1.96, plot_anomalies=False):

    """
        series - 时序dateframe
        window - 滑窗大小 
        plot_intervals - 显示置信区间
        plot_anomalies - 显示异常值 
    """
    rolling_mean = series.rolling(window=window).mean()

    plt.figure(figsize=(15,5))
    plt.title("Moving average\n window size = {}".format(window))
    plt.plot(rolling_mean, "g", label="Rolling mean trend")

    # 绘制平滑后的数据的置信区间
    if plot_intervals:
        mae = mean_absolute_error(series[window:], rolling_mean[window:])
        deviation = np.std(series[window:] - rolling_mean[window:])
        lower_bond = rolling_mean - (mae + scale * deviation)
        upper_bond = rolling_mean + (mae + scale * deviation)
        plt.plot(upper_bond, "r--", label="Upper Bond / Lower Bond")
        plt.plot(lower_bond, "r--")

        # 得到区间后,找出异常值
        if plot_anomalies:
            anomalies = pd.DataFrame(index=series.index, columns=series.columns)
            anomalies[series<lower_bond] = series[series<lower_bond]
            anomalies[series>upper_bond] = series[series>upper_bond]
            plt.plot(anomalies, "ro", markersize=10)

    plt.plot(series[window:], label="Actual values")
    plt.legend(loc="upper left")
plt.grid(True)

平滑(窗口大小为4小时):

plotMovingAverage(ads, 4)

平滑窗口4小时

平滑(窗口大小为12小时):

plotMovingAverage(ads, 12)

平滑窗口12小时

平滑(窗口大小为24小时):

plotMovingAverage(ads, 24)

平滑窗口24小时

如你所见,在小时数据上按日平滑让我们可以清楚地看到浏览广告的趋势。周末数值较高(周末是娱乐时间),工作日一般数值较低。

我们可以同时绘制平滑值的置信区间:

plotMovingAverage(ads, 4, plot_intervals=True)

平滑值的置信区间

现在让我们在移动平均的帮助下创建一个简单的异常检测系统。不幸的是,在这段时序数据中,一切都比较正常,所以让我们故意弄出点异常