Facebook 时间序列预测算法 Prophet 的研究

Prophet 简介

Facebook 去年开源了一个时间序列预测的算法,叫做 fbprophet,它的官方网址与基本介绍来自于以下几个网站:

  1. Github:https://github.com/facebook/prophet
  2. 官方网址:https://facebook.github.io/prophet/
  3. 论文名字与网址:Forecasting at scale,https://peerj.com/preprints/3190/

从官网的介绍来看,Facebook 所提供的 prophet 算法不仅可以处理时间序列存在一些异常值的情况,也可以处理部分缺失值的情形,还能够几乎全自动地预测时间序列未来的走势。从论文上的描述来看,这个 prophet 算法是基于时间序列分解和机器学习的拟合来做的,其中在拟合模型的时候使用了 pyStan 这个开源工具,因此能够在较快的时间内得到需要预测的结果。除此之外,为了方便统计学家,机器学习从业者等人群的使用,prophet 同时提供了 R 语言和 Python 语言的接口。从整体的介绍来看,如果是一般的商业分析或者数据分析的需求,都可以尝试使用这个开源算法来预测未来时间序列的走势。

Prophet 的算法原理

Prophet 数据的输入和输出

prophetexample1

首先让我们来看一个常见的时间序列场景,黑色表示原始的时间序列离散点,深蓝色的线表示使用时间序列来拟合所得到的取值,而浅蓝色的线表示时间序列的一个置信区间,也就是所谓的合理的上界和下界。prophet 所做的事情就是:

  1. 输入已知的时间序列的时间戳和相应的值;
  2. 输入需要预测的时间序列的长度;
  3. 输出未来的时间序列走势。
  4. 输出结果可以提供必要的统计指标,包括拟合曲线,上界和下界等。

就一般情况而言,时间序列的离线存储格式为时间戳和值这种格式,更多的话可以提供时间序列的 ID,标签等内容。因此,离线存储的时间序列通常都是以下的形式。其中 date 指的是具体的时间戳,category 指的是某条特定的时间序列 id,value 指的是在 date 下这个 category 时间序列的取值,label 指的是人工标记的标签(’0′ 表示异常,’1‘ 表示正常,’unknown’ 表示没有标记或者人工判断不清)。

prophetexample2.png

而 fbprophet 所需要的时间序列也是这种格式的,根据官网的描述,只要用 csv 文件存储两列即可,第一列的名字是 ‘ds’, 第二列的名称是 ‘y’。第一列表示时间序列的时间戳,第二列表示时间序列的取值。通过 prophet 的计算,可以计算出 yhat,yhat_lower,yhat_upper,分别表示时间序列的预测值,预测值的下界,预测值的上界。两份表格如下面的两幅图表示。

prophetexample3

prophetexample4

Prophet 的算法实现

在时间序列分析领域,有一种常见的分析方法叫做时间序列的分解(Decomposition of Time Series),它把时间序列 y_{t} 分成几个部分,分别是季节项 S_{t},趋势项 T_{t},剩余项 R_{t}。也就是说对所有的 t\geq 0,都有

y_{t} = S_{t} + T_{t} + R_{t}.

除了加法的形式,还有乘法的形式,也就是:

y_{t} = S_{t} \times T_{t} \times R_{t}.

以上式子等价于 \ln y_{t} = \ln S_{t} + \ln T_{t} + \ln R_{t}。所以,有的时候在预测模型的时候,会先取对数,然后再进行时间序列的分解,就能得到乘法的形式。在 fbprophet 算法中,作者们基于这种方法进行了必要的改进和优化。

一般来说,在实际生活和生产环节中,除了季节项,趋势项,剩余项之外,通常还有节假日的效应。所以,在 prophet 算法里面,作者同时考虑了以上四项,也就是:

y(t) = g(t) + s(t) + h(t) + \epsilon_{t}.

其中 g(t) 表示趋势项,它表示时间序列在非周期上面的变化趋势;s(t) 表示周期项,或者称为季节项,一般来说是以周或者年为单位;h(t) 表示节假日项,表示在当天是否存在节假日;\epsilon_{t} 表示误差项或者称为剩余项。Prophet 算法就是通过拟合这几项,然后最后把它们累加起来就得到了时间序列的预测值。

趋势项模型 g(t)

在 Prophet 算法里面,趋势项有两个重要的函数,一个是基于逻辑回归函数(logistic function)的,另一个是基于分段线性函数(piecewise linear function)的。

首先,我们来介绍一下基于逻辑回归的趋势项是怎么做的。

如果回顾逻辑回归函数的话,一般都会想起这样的形式:\sigma(x) = 1/(1+e^{-x}), 它的导数是 \sigma'(x) = \sigma(x) \cdot(1-\sigma(x)), 并且 \lim_{x\rightarrow +\infty} \sigma(x) = 1, \lim_{x\rightarrow -\infty} \sigma(x) = 0. 如果增加一些参数的话,那么逻辑回归就可以改写成:

f(x) = C / (1 + e^{-k(x-m)}),

这里的 C 称为曲线的最大渐近值,k 表示曲线的增长率,m 表示曲线的中点。当 C=1, k = 1, m =0 时,恰好就是大家常见的 sigmoid 函数的形式。从 sigmoid 的函数表达式来看,它满足以下的微分方程:y'=y(1-y)

那么,如果使用分离变量法来求解微分方程 y'=y(1-y) 就可以得到:

\frac{y'}{y} + \frac{y'}{1-y} = 1 \Rightarrow \ln\frac{y}{1-y} = 1 \Rightarrow y = 1/(1+K e^{-x}).

但是在现实环境中,函数 f(x) = C / (1+e^{-k(x-m)}) 的三个参数 C, k, m 不可能都是常数,而很有可能是随着时间的迁移而变化的,因此,在 Prophet 里面,作者考虑把这三个参数全部换成了随着时间而变化的函数,也就是 C = C(t), k = k(t), m = m(t)

除此之外,在现实的时间序列中,曲线的走势肯定不会一直保持不变,在某些特定的时候或者有着某种潜在的周期曲线会发生变化,这种时候,就有学者会去研究变点检测,也就是所谓 change point detection。例如下面的这幅图的 t_{1}^{*}, t_{2}^{*} 就是时间序列的两个变点。

prophetchangepoint1

在 Prophet 里面,是需要设置变点的位置的,而每一段的趋势和走势也是会根据变点的情况而改变的。在程序里面有两种方法,一种是通过人工指定的方式指定变点的位置;另外一种是通过算法来自动选择。在默认的函数里面,Prophet 会选择 n_changepoints = 25 个变点,然后设置变点的范围是前 80%(changepoint_range),也就是在时间序列的前 80% 的区间内会设置变点。通过 forecaster.py 里面的 set_changepoints 函数可以知道,首先要看一些边界条件是否合理,例如时间序列的点数是否少于 n_changepoints 等内容;其次如果边界条件符合,那变点的位置就是均匀分布的,这一点可以通过 np.linspace 这个函数看出来。

下面假设已经放置了 S 个变点了,并且变点的位置是在时间戳 s_{j}, 1\leq j\leq S 上,那么在这些时间戳上,我们就需要给出增长率的变化,也就是在时间戳 s_{j} 上发生的 change in rate。可以假设有这样一个向量:\boldsymbol{\delta}\in\mathbb{R}^{S}, 其中 \delta_{j} 表示在时间戳 s_{j} 上的增长率的变化量。如果一开始的增长率我们使用 k 来代替的话,那么在时间戳 t 上的增长率就是 k + \sum_{j:t>s_{j}} \delta_{j},通过一个指示函数 \mathbf{a}(t)\in \{0,1\}^{S} 就是

a_{j}(t) = \begin{cases} 1, \text{ if } t\geq s_{j},\\ 0, \text{ otherwise.} \end{cases}

那么在时间戳 t 上面的增长率就是 k + \mathbf{a}^{T}\boldsymbol{\delta}. 一旦变化量 k 确定了,另外一个参数 m 也要随之确定。在这里需要把线段的边界处理好,因此通过数学计算可以得到:

\gamma_{j} = \bigg(s_{j} - m - \sum_{\ell <j} \gamma_{\ell} \bigg) \cdot \bigg( 1- \frac{k + \sum_{\ell < j} \delta_{\ell}}{k + \sum_{\ell\leq j}\delta_{\ell}} \bigg).

所以,分段的逻辑回归增长模型就是:

g(t) = \frac{C(t)}{1+exp(-(k+\boldsymbol{a}(t)^{t}\boldsymbol{\delta}) \cdot (t - (m+\boldsymbol{a}(t)^{T}\boldsymbol{\gamma})},

其中,

\boldsymbol{a}(t) = (a_{1}(t),\cdots,a_{S}(t))^{T},  \boldsymbol{\delta} = (\delta_{1},\cdots,\delta_{S})^{T}, \boldsymbol{\gamma} = (\gamma_{1},\cdots,\gamma_{S})^{T}.

在逻辑回归函数里面,有一个参数是需要提前设置的,那就是 Capacity,也就是所谓的 C(t) ,在使用 Prophet 的 growth = ‘logistic’ 的时候,需要提前设置好 C(t) 的取值才行。

再次,我们来介绍一下基于分段线性函数的趋势项是怎么做的。众所周知,线性函数指的是 y=kx+b, 而分段线性函数指的是在每一个子区间上,函数都是线性函数,但是在整段区间上,函数并不完全是线性的。正如下图所示,分段线性函数就是一个折线的形状。

prophetpiecewiselinear1

因此,基于分段线性函数的模型形如:

g(t)=(k+\boldsymbol{a}(t)\boldsymbol{\delta})\cdot t+(m+\boldsymbol{a}(t)^{T}\boldsymbol{\gamma}),

其中 k 表示增长率(growth rate),\boldsymbol{\delta} 表示增长率的变化量,m 表示 offset parameter。而这两种方法(分段线性函数与逻辑回归函数)最大的区别就是 \boldsymbol{\gamma} 的设置不一样,在分段线性函数中,\boldsymbol{\gamma}=(\gamma_{1},\cdots,\gamma_{S})^{T}, \gamma_{j}=-s_{j}\delta_{j}. 注意:这与之前逻辑回归函数中的设置是不一样的。

在 prophet 的源代码中,forecast.py 这个函数里面包含了最关键的步骤,其中 piecewise_logistic 函数表示了前面所说的基于逻辑回归的增长函数,它的输入包含了 cap 这个指标,因此需要用户事先指定 capacity。而在 piecewise_linear 这个函数中,是不需要 capacity 这个指标的,因此 m = Prophet() 这个函数默认的使用 growth = ‘linear’ 这个增长函数,也可以写作 m = Prophet(growth = ‘linear’);如果想用 growth = ‘logistic’,就要这样写:

m = Prophet(growth='logistic')
df['cap'] = 6
m.fit(df)
future = m.make_future_dataframe(periods=prediction_length, freq='min')
future['cap'] = 6

变点的选择(Changepoint Selection)

在介绍变点之前,先要介绍一下 Laplace 分布,它的概率密度函数为:

f(x|\mu, b) = exp\bigg(-|x-\mu|/b\bigg)/2b,

其中 \mu 表示位置参数,b>0 表示尺度参数。Laplace 分布与正态分布有一定的差异。

在 Prophet 算法中,是需要给出变点的位置,个数,以及增长的变化率的。因此,有三个比较重要的指标,那就是

  1. changepoint_range,
  2. n_changepoint,
  3. changepoint_prior_scale。

changepoint_range 指的是百分比,需要在前 changepoint_range 那么长的时间序列中设置变点,在默认的函数中是 changepoint_range = 0.8。n_changepoint 表示变点的个数,在默认的函数中是 n_changepoint = 25。changepoint_prior_scale 表示变点增长率的分布情况,在论文中, \delta_{j} \sim Laplace(0,\tau),这里的 \tau 就是 change_point_scale。

在整个开源框架里面,在默认的场景下,变点的选择是基于时间序列的前 80% 的历史数据,然后通过等分的方法找到 25 个变点(change points),而变点的增长率是满足 Laplace 分布 \delta_{j} \sim Laplace (0,0.05) 的。因此,当 \tau 趋近于零的时候,\delta_{j} 也是趋向于零的,此时的增长函数将变成全段的逻辑回归函数或者线性函数。这一点从 g(t) 的定义可以轻易地看出。

对未来的预估(Trend Forecast Uncertainty)

从历史上长度为 T 的数据中,我们可以选择出 S 个变点,它们所对应的增长率的变化量是 \delta_{j} \sim Laplace(0,\tau)。此时我们需要预测未来,因此也需要设置相应的变点的位置,从代码中看,在 forecaster.py 的 sample_predictive_trend 函数中,通过 Poisson 分布等概率分布方法找到新增的 changepoint_ts_new 的位置,然后与 changepoint_t 拼接在一起就得到了整段序列的 changepoint_ts。

changepoint_ts_new = 1 + np.random.rand(n_changes) * (T - 1)
changepoint_ts = np.concatenate((self.changepoints_t, changepoint_ts_new))

第一行代码的 1 保证了 changepoint_ts_new 里面的元素都大于 change_ts 里面的元素。除了变点的位置之外,也需要考虑 \delta 的情况。这里令 \lambda = \sum_{j=1}^{S}|\delta_{j}|/S,于是新的增长率的变化量就是按照下面的规则来选择的:当 j>T 时,

\delta_{j}=\begin{cases} 0 \text{, with probability } (T-S)/T \\ \sim Laplace(0,\lambda) \text{, with probability } S/T \end{cases}.

季节性趋势

几乎所有的时间序列预测模型都会考虑这个因素,因为时间序列通常会随着天,周,月,年等季节性的变化而呈现季节性的变化,也称为周期性的变化。对于周期函数而言,大家能够马上联想到的就是正弦余弦函数。而在数学分析中,区间内的周期性函数是可以通过正弦和余弦的函数来表示的:假设 f(x) 是以 2\pi 为周期的函数,那么它的傅立叶级数就是 a_{0} + \sum_{n=1}^{\infty}(a_{n}\cos(nx) + b_{n}\sin(nx))

在论文中,作者使用傅立叶级数来模拟时间序列的周期性。假设 P 表示时间序列的周期,P = 365.25 表示以年为周期,P = 7 表示以周为周期。它的傅立叶级数的形式都是:

s(t) = \sum_{n=1}^{N}\bigg( a_{n}\cos\bigg(\frac{2\pi n t}{P}\bigg) + b_{n}\sin\bigg(\frac{2\pi n t}{P}\bigg)\bigg).

就作者的经验而言,对于以年为周期的序列(P = 365.25)而言,N = 10;对于以周为周期的序列(P = 7 )而言,N = 3。这里的参数可以形成列向量:

\boldsymbol{\beta} = (a_{1},b_{1},\cdots,a_{N},b_{N})^{T}

N = 10 时,

X(t) = \bigg[\cos(\frac{2\pi(1)t}{365.25}),\cdots,\sin(\frac{2\pi(10)t}{365.25})\bigg]

N = 3 时,

X(t) = \bigg[\cos(\frac{2\pi(1)t}{7}),\cdots,\sin(\frac{2\pi(3)t}{7})\bigg]

因此,时间序列的季节项就是:s(t) = X(t) \boldsymbol{\beta},\boldsymbol{\beta} 的初始化是 \boldsymbol{\beta} \sim Normal(0,\sigma^{2})。这里的 \sigma 是通过 seasonality_prior_scale 来控制的,也就是说 \sigma= seasonality_prior_scale。这个值越大,表示季节的效应越明显;这个值越小,表示季节的效应越不明显。同时,在代码里面,seasonality_mode 也对应着两种模式,分别是加法和乘法,默认是加法的形式。在开源代码中,X(t) 函数是通过 fourier_series 来构建的。

节假日效应(holidays and events)

在现实环境中,除了周末,同样有很多节假日,而且不同的国家有着不同的假期。在 Prophet 里面,通过维基百科里面对各个国家的节假日的描述,hdays.py 收集了各个国家的特殊节假日。除了节假日之外,用户还可以根据自身的情况来设置必要的假期,例如 The Super Bowl,双十一等。

prophetholiday1.png

由于每个节假日对时间序列的影响程度不一样,例如春节,国庆节则是七天的假期,对于劳动节等假期来说则假日较短。因此,不同的节假日可以看成相互独立的模型,并且可以为不同的节假日设置不同的前后窗口值,表示该节假日会影响前后一段时间的时间序列。用数学语言来说,对与第 i 个节假日来说, D_{i} 表示该节假日的前后一段时间。为了表示节假日效应,我们需要一个相应的指示函数(indicator function),同时需要一个参数 \kappa_{i} 来表示节假日的影响范围。假设我们有 L 个节假日,那么

h(t)=Z(t) \boldsymbol{\kappa}=\sum_{i=1}^{L} \kappa_{i}\cdot 1_{\{t\in D_{i}\}},

其中 Z(t)=(1_{\{t\in D_{1}\}},\cdots,1_{\{t\in D_{L}\}})\boldsymbol{\kappa}=(\kappa_{1},\cdots,\kappa_{L})^{T}.

其中 \boldsymbol{\kappa}\sim Normal(0,v^{2}) 并且该正态分布是受到 v = holidays_prior_scale 这个指标影响的。默认值是 10,当值越大时,表示节假日对模型的影响越大;当值越小时,表示节假日对模型的效果越小。用户可以根据自己的情况自行调整。

模型拟合(Model Fitting)

按照以上的解释,我们的时间序列已经可以通过增长项,季节项,节假日项来构建了,i.e.

y(t)=g(t)+s(t)+h(t)+\epsilon

下一步我们只需要拟合函数就可以了,在 Prophet 里面,作者使用了 pyStan 这个开源工具中的 L-BFGS 方法来进行函数的拟合。具体可以参考 forecast.py 里面的 stan_init 函数。

Prophet 中可以设置的参数

在 Prophet 中,用户一般可以设置以下四种参数:

  1. Capacity:在增量函数是逻辑回归函数的时候,需要设置的容量值。
  2. Change Points:可以通过 n_changepoints 和 changepoint_range 来进行等距的变点设置,也可以通过人工设置的方式来指定时间序列的变点。
  3. 季节性和节假日:可以根据实际的业务需求来指定相应的节假日。
  4. 光滑参数:\tau= changepoint_prior_scale 可以用来控制趋势的灵活度,\sigma= seasonality_prior_scale 用来控制季节项的灵活度,v= holidays prior scale 用来控制节假日的灵活度。

如果不想设置的话,使用 Prophet 默认的参数即可。

 

Prophet 的实际使用

Prophet 的简单使用

因为 Prophet 所需要的两列名称是 ‘ds’ 和 ‘y’,其中,’ds’ 表示时间戳,’y’ 表示时间序列的值,因此通常来说都需要修改 pd.dataframe 的列名字。如果原来的两列名字是 ‘timestamp’ 和 ‘value’ 的话,只需要这样写:

df = df.rename(columns={'timestamp':'ds', 'value':'y'})

如果 ‘timestamp’ 是使用 unixtime 来记录的,需要修改成 YYYY-MM-DD hh:mm:ss 的形式:

df['ds'] = pd.to_datetime(df['ds'],unit='s')

在一般情况下,时间序列需要进行归一化的操作,而 pd.dataframe 的归一化操作也十分简单:

df['y'] = (df['y'] - df['y'].mean()) / (df['y'].std())

然后就可以初始化模型,然后拟合模型,并且进行时间序列的预测了。

初始化模型:m = Prophet()
拟合模型:m.fit(df)
计算预测值:periods 表示需要预测的点数,freq 表示时间序列的频率。
future = m.make_future_dataframe(periods=30, freq='min')
future.tail()
forecast = m.predict(future)

而 freq 指的是 pd.dataframe 里面的一个指标,’min’ 表示按分钟来收集的时间序列。具体参见文档:http://pandas.pydata.org/pandas-docs/stable/timeseries.html#offset-aliases

prophetdataframefrequency

在进行了预测操作之后,通常都希望把时间序列的预测趋势画出来:

画出预测图:
m.plot(forecast)
画出时间序列的分量:
m.plot_components(forecast)

prophetexample5prophetexample6

如果要画出更详细的指标,例如中间线,上下界,那么可以这样写:

x1 = forecast['ds']
y1 = forecast['yhat']
y2 = forecast['yhat_lower']
y3 = forecast['yhat_upper']
plt.plot(x1,y1)
plt.plot(x1,y2)
plt.plot(x1,y3)
plt.show()

prophetexample7

其实 Prophet 预测的结果都放在了变量 forecast 里面,打印结果的话可以这样写:第一行是打印所有时间戳的预测结果,第二行是打印最后五个时间戳的预测结果。

print(forecast[['ds', 'yhat', 'yhat_lower', 'yhat_upper']])
print(forecast[['ds', 'yhat', 'yhat_lower', 'yhat_upper']].tail())

 

Prophet 的参数设置

Prophet 的默认参数可以在 forecaster.py 中看到:

def __init__(
    self,
    growth='linear',
    changepoints=None,
    n_changepoints=25, 
    changepoint_range=0.8,
    yearly_seasonality='auto',
    weekly_seasonality='auto',
    daily_seasonality='auto',
    holidays=None,
    seasonality_mode='additive',
    seasonality_prior_scale=10.0,
    holidays_prior_scale=10.0,
    changepoint_prior_scale=0.05,
    mcmc_samples=0,
    interval_width=0.80,
    uncertainty_samples=1000,
):

增长函数的设置

在 Prophet 里面,有两个增长函数,分别是分段线性函数(linear)和逻辑回归函数(logistic)。而 m = Prophet() 默认使用的是分段线性函数(linear),并且如果要是用逻辑回归函数的时候,需要设置 capacity 的值,i.e. df[‘cap’] = 100,否则会出错。

m = Prophet()
m = Prophet(growth='linear')
m = Prophet(growth='logistic')

变点的设置

在 Prophet 里面,变点默认的选择方法是前 80% 的点中等距选择 25 个点作为变点,也可以通过以下方法来自行设置变点,甚至可以人为设置某些点。

m = Prophet(n_changepoints=25)
m = Prophet(changepoint_range=0.8)
m = Prophet(changepoint_prior_scale=0.05)
m = Prophet(changepoints=['2014-01-01'])

而变点的作图可以使用:

from fbprophet.plot import add_changepoints_to_plot
fig = m.plot(forecast)
a = add_changepoints_to_plot(fig.gca(), m, forecast)

prophetexample8

周期性的设置

通常来说,可以在 Prophet 里面设置周期性,无论是按月还是周其实都是可以设置的,例如:

m = Prophet(weekly_seasonality=False)
m.add_seasonality(name='monthly', period=30.5, fourier_order=5)
m = Prophet(weekly_seasonality=True)
m.add_seasonality(name='weekly', period=7, fourier_order=3, prior_scale=0.1)

prophetexample9

节假日的设置

有的时候,由于双十一或者一些特殊节假日,我们可以设置某些天数是节假日,并且设置它的前后影响范围,也就是 lower_window 和 upper_window。

playoffs = pd.DataFrame({
  'holiday': 'playoff',
  'ds': pd.to_datetime(['2008-01-13', '2009-01-03', '2010-01-16',
                        '2010-01-24', '2010-02-07', '2011-01-08',
                        '2013-01-12', '2014-01-12', '2014-01-19',
                        '2014-02-02', '2015-01-11', '2016-01-17',
                        '2016-01-24', '2016-02-07']),
  'lower_window': 0,
  'upper_window': 1,
})
superbowls = pd.DataFrame({
  'holiday': 'superbowl',
  'ds': pd.to_datetime(['2010-02-07', '2014-02-02', '2016-02-07']),
  'lower_window': 0,
  'upper_window': 1,
})
holidays = pd.concat((playoffs, superbowls))

m = Prophet(holidays=holidays, holidays_prior_scale=10.0)

 

结束语

对于商业分析等领域的时间序列,Prophet 可以进行很好的拟合和预测,但是对于一些周期性或者趋势性不是很强的时间序列,用 Prophet 可能就不合适了。但是,Prophet 提供了一种时序预测的方法,在用户不是很懂时间序列的前提下都可以使用这个工具得到一个能接受的结果。具体是否用 Prophet 则需要根据具体的时间序列来确定。

参考文献:

  1. https://otexts.org/fpp2/components.html
  2. https://en.wikipedia.org/wiki/Decomposition_of_time_series
  3. A review of change point detection methods, CTruong, L. Oudre, N.Vayatis
  4. https://github.com/facebook/prophet
  5. https://facebook.github.io/prophet/

 

Advertisements

一篇关于时间序列异常检测的论文

近期阅读了一篇论文《Rapid Deployment of Anomaly Detection Models for Large Number of Emerging KPI Streams》,这篇文章基于之前的 ROCKA 系统做了一些额外的工作。ROCKA 系统是用来做时间序列的实时聚类的,而这篇文章是在 ROCKA 系统的基础上增加了时间序列异常检测的功能。通常来说,时间序列异常检测可以使用有监督的方法来解决,参考 Opperentice 系统。而本篇文章使用了半监督学习的思路来解决异常检测的问题,下面来详细分析一下这篇文章的细节,本文的作者把这个系统称为 ADS(Anomaly Detection through Self-training)。

数据集的情况

在论文中,作者使用了两份数据集,分别是已经历史上的 70 条时间序列,另外还有新来的 81 条时间序列。在 ADS 系统中,历史上的 70 条时间序列被划分成 5 类,并且已经可以找出每个类的质心位置,并且每条历史上的时间序列通常来说会大于三个星期(3 weeks)。本篇论文的评价指标是 F-Score,也属于机器学习领域里面比较常用的衡量模型效果的指标。整体来看,这篇文章的数据集大约是 200 条时间序列,时间序列的时间间隔通常来说是五分钟(不过其余的运维场景会有一分钟的数据采集粒度),而一般来说都拥有大半年甚至一年的时间跨度。那么时间点的个数预估是 200 * (1440 / 5) * 365。假设异常的数据:正常的数据 = 1:10000(也就是说平均每条时间序列每周至少发生一次异常),于是这批时间序列数据的异常数据量大约是 200 * (1440 / 5) * 365 / 10000 = 2102,也就是说异常的样本大约是 2102 个左右,剩下的都是正常的样本。PS:当然如果异常的数据:正常的数据的比例大于 1:10000 的话,异常的样本还会更多一些。整体来看,时间序列异常检测是一个样本极其不均衡的场景。

ADS 的系统架构

按照作者之前论文的经验,时间序列异常检测通常都是先做聚类,然后再根据每一个类的特点来做一个异常检测模型,之前的技术架构就是 ROCKA + Opperentice。因为 ROCKA 可以根据时间序列的走势和趋势来进行时间序列的实时分类/聚类,然后 Opperentice 就是做时间序列异常检测的模型。在本文的场景下,作者把 70 条时间序列分成了5 类,因此只需要维护 5 个时间序列的异常检测模型就可以了。当然把时间序列切分成更多的类也是可以的,只是需要维护的时间序列异常检测就变多了,人工成本会加大。

ADS系统架构
ADS 的系统架构
ROCKA系统架构
ROCKA 的系统架构

如果看到上面两幅图,有心的读者一定会发现其实 ADS 就是基于 ROCKA 所做的工作。ADS 先对时间序列进行了分类,然后进行了特征提取的工作,再通过半监督学习模型,最后进行异常检测。也就是说,ADS 会走下面四个步骤:

  1. ADS 先把历史上的时间序列进行聚类;
  2. 通过算法获得每一个类的质心,并且标记出质心曲线的异常点和正常点;
  3. 对新来的时间序列进行实时聚类,划分到合适的类别;
  4. 基于新来的时间序列(没有标记)和历史上的时间序列(有标记)使用无监督算法来重新训练一个新的模型,进行该类别的时间序列异常检测。

ADS 的细节

对于时间序列的聚类框架 ROCKA,之前的一篇 BLOG 里面已经详细介绍过,这里将不会再赘述。而 ADS 的另一个模块就是半监督学习算法 Contrastive Pessimistic Likelihood Estimation(CPLE),详细的论文细节可以参考论文《Contrastive Pessimistic Likelihood Estimation for Semi-Supervised Classification》。CPLE 有几个好处:

  1. CPLE 是半监督学习算法中比较健壮的,因为它并没有过多的假设条件,并且也符合这篇文章的业务场景,同时拥有质心曲线(有标记)和新来的曲线(无标记),使用半监督学习也是符合常理的。除了 CPLE,其实在实战过程中也可以多尝试其他的半监督模型,具体可以参考周志华的《机器学习》。
  2. CPLE 的复杂度比较低,计算快。
  3. CPLE 支持增量学习,因此,当越来越多新的时间序列进入 ADS,这个模型也会随之而调整并提高准确率。

整体来看,ADS = ROCKA + CPLE,而在论文中,它的对比模型就是 ROCKA + Opperentice。而且在 CPLE 中,也使用了与 Opperentice 系统类似的特征,如下图所示。

ADS特征
特征工程

其实,从本质上来看,就是半监督学习与有监督学习在这份数据集合上面的比较。从这篇论文里面所展示的数据来看,CPLE 有一定的优势。

ADS效果对比1
Average best F-scores of ADS, iForest, Donut, Opperentice, ROCKA + Opperentice
ADS效果对比2
The New KPI Streams

整体来看,本篇文章介绍了时间序列异常检测的一种方案,也就是把时间序列先进行聚类的操作,然后根据不同的类来进行异常检测。在异常检测的方法中,不仅可以使用 Random Forest,GBDT,XGBoost 等有监督学习方法,也可以使用 CPLE 等半监督算法。具体在业务中如何使用,其实只能够根据具体的数据来进行合理地选择了。

 

两篇关于时间序列的论文

这次整理的就是清华大学裴丹教授所著的两篇与时间序列相关的论文。一篇是关于时间序列聚类的,《Robust and Rapid Clustering of KPIs for Large-Scale Anomaly Detection》;另外一篇文章是关于时间序列异常检测的,重点检测时间序列上下平移的,《Robust and Rapid Adaption for Concept Drift in Software System Anomaly Detection》。本文将会整理一下这两篇文章的关键技术点。

Robust and Rapid Clustering of KPIs for Large-Scale Anomaly Detection

在互联网公司中,通常会拥有海量的的时间序列,而海量的时间序列就有着各种各样的形状和走势。因此,就有学者提出可以先对时间序列进行分类,然后根据不同的类使用不同的检测模型来进行异常检测。如果要做时间序列的分类,就先需要做聚类的操作,无论从 KMeans,DBSCAN,还是层次聚类来说,都会消耗一定的运算时间。所以,如何在较短的时间内进行聚类或者分类的操作则是这个系统的关键之处。于是,这篇文章提出了一个将时间序列快速聚类的方法。

时间序列 -> 时间序列分类

-> 根据每一类时间序列使用不同的异常检测模型

而在做时间序列聚类的时候,也有着不少的挑战。通常挑战来自于以下几点:

  1. 形状:通常来说,时间序列随着业务的变化,节假日效应,变更的发布,将会随着时间的迁移而造成形状的变化。
  2. 噪声:无论是从数据采集的角度,还是系统处理的角度,甚至服务器的角度,都有可能给时间序列带来一定的噪声数据,而噪声是需要处理掉的。
  3. 平移:在定时任务中,有可能由于系统或者人为的原因,时间序列的走势可能会出现一定程度的左右偏移,有可能每天 5:00 起的定时任务由于前序任务的原因而推迟了。
  4. 振幅:通常时间序列都存在一条基线,而不同的时间序列有着不同的振幅,振幅决定了这条时间序列的振荡程度,而振幅或者基线其实也是会随着时间的迁移而变化的。

从整篇论文来看,ROCKA 系统是为了做实时的时间序列分类判断的。要想做成实时的分类判断,就需要有离线和在线两个模块。其中离线是为了做模型训练或者聚类的,在线是为了使用离线处理好的模块来做曲线分类的。

ROCKA系统架构
ROCKA系统架构

从整个系统来看,离线模块需要做以下几件事情:首先需要收集一批时间序列数据,也就是所谓的 Raw Time Series Data(Raw),通过预处理模块,实施基线提取,再进行聚类的操作,获得相应的聚类结果和质心。在线模块同样也要做类似的事情:首先对于每一条新来的时间序列数据,也就是所谓的 New Time Series Data(Raw),通过预处理模块,实施基线提取,然后使用已经聚类好的离线模块来进行实时的分类。

下面,我们来逐一分析每个模块的作用。

预处理模块(Preprocessing)

通常预处理模块包含几个关键点:

  1. 缺失值:缺失值指的是在该上报数据的时间戳上并没有相应的数据上报,数据处于缺失的状态。通常的办法就是把数据补齐,而数据补齐的方法有很多种,最简单的就是使用线性插值(Linear Interpolation)的方式来补齐。
  2. 标准化:对于一个时间序列而言,有可能它的均值是10万,有可能只有10,但是它们的走势有可能都是一样的。所以在这个时候需要进行归一化的操作。最常见的有两种归一化方法:一种是标准化,另外一种是最大最小值归一化。如果 [x_{1},x_{2}, \cdots, x_{n}] 表示原始的时间序列的话,标准化指的是 \hat{x}_{i} = (x_{i} -\mu) /\sigma,其中 \mu\sigma 分别表示均值和标准差。最大最小值归一化指的是 \hat{x}_{i} = (x_{i} - \min) / (\max - \min),其中 \max, \min 分别表示这段时间内的最大值与最小值。

基线提取模块(Baseline Extraction)

基线提取指的是把时间序列分成基线和剩余项两个部分,假设时间序列是 [x_{1},x_{2},\cdots,x_{n}],基线提取就是:

x_{i} = baseline_{i} + residual_{i},

其中 baseline_{i}, residual_{i} 分别指的是 x_{i} 的基线和剩余项。

在处理基线的时候,有几件事情需要注意:

  1. 异常值的处理:通常需要移除一些明显异常的值,然后使用线性插值等方法来把这些移除的值补上。
  2. 使用简单的移动平均算法加上一个窗口值 w 来提取基线。假设时间序列是 [x_{1},\cdots,x_{n}]SMA_{i} = \sum_{j=i-w+1}^{i} x_{j} / wR_{i} = x_{i} - SMA_{i}。也就是说 x_{i} = SMA_{i} + R_{i}
  3. 提取基线的方式其实还有很多,使用带权重的移动平均算法(Weighted Moving Average),指数移动平均算法(Exponentially Weighted Moving Average)都可以提取基线,甚至使用深度学习中的 Autoencoder 或者 VAE 算法都能够提取基线。

基于密度的聚类算法(Clustering)

使用预处理和基线提取技术之后,可以得到原始时间序列的基线值,然后根据这些基线值来进行时间序列的聚类操作。一般来说,时间序列的聚类有两种方法:

  1. 通过特征工程的方法,从时间序列中提取出必要的时间序列特征,然后使用 KMeans 等算法来进行聚类。
  2. 通过相似度计算工具,对比两条时间序列之间的相似度,相似的聚成一类,不相似的就分成两类。

对于第一种方法而言,最重要的就是特征工程;对于第二种方法而言,最重要的是相似度函数的选择。在这篇文章中,作者选择了第二种方法来进行时间序列的聚类。对于两条时间序列 X = [x_{1},\cdots,x_{m}]Y = [y_{1},\cdots,y_{m}] 而言,为了解决左右平移的问题,需要考虑一个偏移量 |s|,然后计算它们之间的内积。

CC_{s}(X,Y) = \begin{cases}  \sum_{i=1}^{m-s}x_{s+i}\cdot y_{i}, \text{ where } s\geq 0, \\ \sum_{i=1}^{m+s}x_{i} \cdot y_{i-s} \text{ where } s<0.\end{cases}

通过这个偏移量 |s| 就可以计算出最大的相似度,然后计算出两条时间序列之间的距离。i.e.

NCC(X,Y) = \max_{s}\bigg(\frac{CC_{s}(X,Y)}{||X||_{2}\cdot||Y||_{2}}\bigg),

SBD(X,Y) = 1 - NCC(X,Y).

其中 NCC \in [-1,1] 指的是 Normalized version of Cross-Correlation,SBD \in [0,2] 指的是 Shape-based distance。

而聚类的另一个重要指标就是质心的选择,在这里,每个类的质心可以设置为:

\text{Centroid} = argmin_{X \in \text{cluster}_{i}} \sum_{Y \in \text{cluster}_{i}} SBD(X,Y)^{2}.

实时分类(Assignment)

对于一条新的时间序列(Raw Data),同样需要经过预处理,基线提取等步骤,然后计算它与之前每一个质心的距离。然后进行距离的从小到大排序,最小的那一类可能就是所需要的。当然,当最小距离大于某个 threshold ,就说明这条新的时间序列曲线很可能不属于之前的任何一类。通过人工查看之后,可以考虑新增一类,并且更新之前所做的聚类模型。

聚类与时间序列异常检测

如果要做海量的时间序列异常检测的话,通常有以下两种做法。

  1. 先对时间序列进行聚类或者分类,针对不同的时间序列类型来使用不同的模型;
  2. 在时间序列异常检测中加入分类特征。

对于第一种方法而言,需要针对不同形状的时间序列维护不同的模型,而且如果第一层的聚类/分类错误了,那么使用的模型也会出现错误。对于第二种方法而言,关键在于样本的积累和分类特征的构建。

Robust and Rapid Adaption for Concept Drift in Software System Anomaly Detection

这篇文章是关于时间序列异常检测的,而清华大学的 Netman 实验室做时间序列异常检测的相关工作比较多。从 2015 年开始的 Opperentice(Random Forest 模型),Funnel(SST 模型,变更相关),到后来的 Donut(VAE 模型),都是时间序列异常检测的相关文章。而这篇问题提到的 StepWise 系统针对的场景是关于指标的迁移的,所谓概念漂移(Concept Drift)指的就是时间序列会随着变更,发布,调度或者其他事件的发生而导致上下漂移或者左右漂移。一个比较典型的例子就是关于网络流量的漂移,通常来说在某个特点的时间点,时间序列会出现猛跌或者猛涨的情况,但是下跌或者上涨之后的走势和历史数据是及其相似的,只是绝对值上有所变化。

StepWiseExample.png
概念漂移的例子

正如上图所示,在时间戳 08-02 附近,时间序列出现了一个下跌的情况。但是根据历史数据所计算出来的上界(Upper Bound)和下界(Lower Bound)却会把未来一段时间的序列都判断为异常。但是前后的曲线走势却是一样的,其实只有下跌的那一小段时间应该被判断为异常,其他时间都是正常的。基于以上的业务场景,这篇文章的作者就提出了概念漂移(Concept Drift)的一些方法。

在整个 StepWise 系统背后,有两个比较关键的地方。其中第一个关键之处就是如何判断异常点。在这个业务场景下,需要检测的是上图发生暴跌的点。而暴跌或者暴涨只是业务运维用来描述时间序列的一个词语。在统计学领域,这种点通常称为变点(Change Point)。所以概念漂移的检测可以转化为一个变点检测的问题,正是论文里面写的。

StepWise系统架构
StepWise 的系统架构

Insight 1:Concept drift detection can be converted into a spike detection problem in the change score stream.

如果是进行变点检测的话,其实可以参考 2016 年的论文 Funnel:Assessing Software Changes in Web-based Services。里面使用了 Singular Spectrum Transform 来进行检测,并且使用 Difference in Difference 来判断时间序列是否真正出现了异常。关于 SST(Singular Spectrum Transform)的具体细节可以参考 Yasser Mohammad, Toyoaki Nishida 所撰写的 Robust Singular Spectrum Transform。SST 算法是基于 SVD 算法的,有一定的时间复杂度。而基于 SST 又有学者提出了 Robust SST,具体可以详见论文。而 StepWise 的其中两步就是 Detection 和 Classification,其中的 Detection 使用了 Improved SST,Classification 使用了Difference in Difference。而 Funnel 系统的其中两步也是 Improved SST 和 Difference in Difference,这与 Funnel 有异曲同工之妙,Funnel 系统详情请见下图。

Funnel系统架构
FUNNEL系统架构

而突变前后的时间序列可以分别叫做 Old Concept 和 New Concept,由于前后的走势几乎一致,所以本文的第二个关键点就是相似度的判断。

Insight 2:The relationship between new concept and old concept can be almost perfectly fitted by linear regression.

从数学的角度来讲,Insight 2 的陈述是想判断两条时间序列是否相似。假设 X=[x_{1},\cdots, x_{n}]Y = [y_{1},\cdots, y_{n}],那么以下陈述是等价的。

  1. 存在一个线性关系 y = kx 使得二维点集 \{(x_{i},y_{i}), 1\leq i\leq n\} 能够被很好的拟合好,也就是说此刻的方差较小。
  2. [x_{1},\cdots,x_{n}][y_{1},\cdots,y_{n}] 的 Pearson 系数很高;
  3. [x_{1},\cdots,x_{n}][y_{1},\cdots,y_{n}] 标准化(z-score)之后的曲线几乎一致。
  4. 分别存在两个值 \mu_{1},\mu_{2} 使得 [x_{1}/\mu_{1},\cdots,x_{n}/\mu_{1}][y_{1}/\mu_{2},\cdots,y_{n}/\mu_{2}] 几乎一致。

所以,在 StepWise 中,作者通过线性回归的算法来判断 old concept 与 new concept 是否存在线性关系,也就是说这两者是否只是平移的关系。其实也可以尝试上面所说的其余方法。

整体来看,StepWise 大致分成两个关键步骤。首先使用类似 Funnel 系统的思想,先进行异常检测(SST),再使用判断一下是不是因为软件的变更引起的(DID),最后使用一层过滤逻辑(线性拟合)来判断是否出现了概念漂移的情况。在实际使用过程中,无论是异常检测还是过滤逻辑,都需要根据具体的业务来做,很难找到一个固定的模型来解决所有的难题。

 

Kelly Criterion

典型案例

无论是在数学课上还是日常生活中,相比大家对抛硬币这个游戏并不陌生。通常来说,硬币有正反两面,并且抛一次硬币出现正面和反面的概率是一样的,也就是二分之一。如果用数学符号来表示的话,那就是 p = q = 1/2, 这里的 p,q 分别表示硬币出现正面和反面的概率值。 从直觉上说,在这种情况下,如果我们每次投入的资金都是一样的并且资金链条不会断裂,那么随着投币次数的增多,最终我们应该是既不会赢钱也不会输钱,总资产和一开始的时候一样保持不变。

数学描述

用数学语言来描述这个过程的话,假设我们的初始资产是 X_{0}X_{n} 表示 n 次下注之后我们当前的资产,每次下注的钱是 B_{k},用 T_{k} =1 表示第 k 次下注获胜,T_{k} = -1 表示第 k 次下注失败。那么从数学公式上来说就是:对于所有的 k\geq 1,都有

X_{k} = X_{k-1} + T_{k}B_{k},

进一步可以得到:

X_{n} = X_{0} + \sum_{k=1}^{n} T_{k}B_{k},

所以,

E(X_{n}) = X_{0} + \sum_{k=1}^{n}E(T_{k}B_{k}) = X_{0} + \sum_{k=1}^{n}(p-q)E(B_{k}).

对于上面特殊的例子,当 p = q = 1/2 时,就可以得到 E(X_{n}) = X_{0}.

p>q 的时候,表示赌赢的概率大于赌输的概率。这种时候,如果需要让 E(X_{n}) 最大,那么就需要最大化每次的 E(B_{k})。所以,在这种情况下,我们需要在每次下注的时候,把手上所有的资产全部下注,然后总资产随着次数的增加就会出现几何级数的增长。

相反的,如果 p<q,表示下注获胜的概率小于下注失败的概率。这种时候,如果需要让 E(X_{n}) 最大,那么就需要最小化每次的 E(B_{k}),因为 p-q<0. 所以,在这种情况下,我们在每次下注的时候,其实就不需要投注,每次的 B_{k}=0,在这种情况下,我们的总资产其实就会保持原样不变。i.e. E(X_{n}) = X_{0}.

按比例投注

在每次下注或者投资的时候,我们通常都会想着按照一定的比例来投注自己当前的资产。也就是说,存在一个参数 0\leq f\leq 1,使得 B_{i} = f X_{i-1},也就是说,每次投注的时候,基于上一轮的总资产来投注相应的比例即可。如果当前的资产是 X,如果下注赢了,那么总资产就变成 X(1+f);如果下注输了,那么总资产就变成 X(1-f)。意思就是说,如果赢了,那就在原来资产的基础上乘以 (1+f);如果输了,那就在原来资产的基础上乘以 (1-f)。在这种情况下,如果我们进行了 n 次下注,那么此刻的资产就变成了

X_{n} = X_{0} \cdot(1+f)^{S}\cdot(1-f)^{F}

其中,S, F 分别表示成功和失败的次数,并且 S+F=n

f=0 时,表示永远不下注,此时对于所有的 n\geq 0,都有 X_{n}=X_{0};当 f=1 时,表示每次下注的时候都是全部下注,那么只要 F\geq 1,就会出现 X_{n}=0 的情况。意思就是说,如果有一次失败了,那就全盘皆输,没有任何资产可以继续运营。但是如果运气足够好的话,那就是 S=n, F=0,总资产就是 X_{n} = X_{0}\cdot 2^{n}.

0<f<1 时,从公式中可以得到:

\ln X_{n} = \ln X_{0} + S\cdot \ln(1+f) + F\cdot \ln(1-f)

\Rightarrow G_{n}(f)=\ln(X_{n}/X_{0})^{1/n}

=(\ln X_{n}-\ln X_{0})/n = (S/n) \ln(1+f) + (F/n)\ln(1-f)

进一步可以得到:

g(f) = E(\ln(X_{n}/X_{0})^{1/n})

= E((S/n) \ln(1+f) + (F/n)\ln(1-f))

= p\ln(1+f)+q\ln(1-f).

因此,可以通过研究 g(f) 来确定 E(X_{n}) 的最大值。为了计算 g(f) 的最大值和最小值,则需要通过导数来研究。可以简单的看出来:当 p,q\in(0,1) 时,g(0) = 0\lim_{f\rightarrow 1^{-}} g(f) = -\infty。计算导数可以得到:

Dg(f) = \frac{p}{1+f} - \frac{q}{1-f} = \frac{p-q-f}{1-f^{2}},

D^{2}g(f) = \frac{-(1-f)^{2}-4fq}{(1-f^{2})^{2}}<0.

所以,导数 Dg(f) 是一个严格递减函数,当 f<p-q 时,Dg(f)>0;当 f>p-q 时,Dg(f)<0。因此,函数 g(f) 会在 f^{*} = p-q 的时候达到最大值。所以,每次最佳的投注比例应该是 f^{*} =p-q

从这个结果上来看,如果 p<q,那么意思就是最好不要做投注,因为输的概率比较大。如果 p>q,那么意思就是投注的比例是 f^{*}=p-q。在这种情况下,可以得到 g(f^{*}) = p\ln p + q\ln q + \ln 2. 从以上信息分析,也可以得到存在一个点 f_{c} 使得 g(f_{c}) = 0

kellycriterion1

一般形式(一)

在一般情况下,在现实生活中都会有赔率的概念。也就是说:赢钱的时候需要考虑赔率,押 1b,也就是所谓的赢钱率,资产从 1 增加到 1+b。举个简单的例子来描述那就是:当硬币有偏差,同时有赔率(赢钱率)的时候该怎么办?

在这种情况下,之前的 g(f) 可以写成如下的形式:

g(f)=p\ln(1+bf)+q\ln(1-f).

这里的 p, q 分别表示赢钱和输钱的概率,并且 p+q=1f 表示每次应该押上的当前总资产的比例,并且 0\leq f\leq 1。通过计算可以直接得到:

Dg(f) = \frac{pb-q-bf}{(1+bf)(1-f)}.

并且 g(0) = 0\lim_{f\rightarrow 1^{-}}g(f) = -\infty。从导数上可以看出,当 f>(pb-q)/b 时,g(f) 是增函数;当 f<(pb-q)/b 时,g(f) 是递减函数。因此,当 f = (pb -q)/b 时,g(f) 达到最大值,也就是

g(f) = p\ln p + q\ln q + \ln(1+b) - q\ln b.

在这种情况下,直接带入一些简单的数字可以得到:

  1. p = q = 0.5, b = 1 时,f^{*} = (pb-q)/b = 0
  2. p = 2/3, q= 1/3, b = 1 时,f^{*} = (pb-q)/b = 1/3;
  3. p = q = 0.5, b = 2 时,f^{*} = (pb-q)/b = 1/4;
  4. p = q = 0.5, b = 1/2 时,f^{*} = (pb-q)/b = - 1/2; 此时不该下注。

从公式 f^{*} = p - \frac{q}{b} 可以看出:

  1. p 增加的时候,q 会减少,f^{*} 也会增加,意思就是随着赢钱的机会加大,就该增加投注的比例;反之,当 p 减少的时候,f^{*} 也会随之减少,也就是说当赢钱的机会变小的时候,就该减少投注的比例。
  2. b 增加的时候,f^{*} 会增加,意思就是说当赢钱率增加的时候,应该加大投注的比例;反之,当赢钱率减少的时候,应该减少投注的比例。
  3. f^{*}\leq 0 时,就表示这个游戏不值得投注,因为总资产的期望是负数。

一般形式(二)

除了赢钱率之外,在一般情况下,还有一个损失率的概念。也就是说,当投掷的硬币有偏差,并且同时有赢钱率和损失率的时候该怎么办?

假设

  1. 赌赢的概率是 p
  2. 赌输的概率是 q
  3. 赢钱率是 b,下注的资产会按比例增加,从 1 变成 1+b
  4. 损失率是 c,下注的资产会按比例减少,从 1 变成 1-c

一般情况下,0\leq c\leq 1。在这种情况下,函数 g(f) 的形式又会发生微小的变化,

g(f) = p\ln(1+bf) + q\ln(1-cf).

直接计算可以得到:Dg(f) = \frac{pb-cq-cbf}{(1+bf)(1-cf)}。因此,临界点是 f^{*} = \frac{p}{c} - \frac{q}{b}。从公式上看,

  1. p 增加的时候,q 会减少,f^{*} 也会增加,意思就是随着赢钱的机会加大,就该增加投注的比例;反之,当 p 减少的时候,f^{*} 也会随之减少,也就是说当赢钱的机会变小的时候,就该减少投注的比例。
  2. b 增加的时候,f^{*} 会增加,意思就是说当赢钱率增加的时候,应该加大投注的比例;反之,当赢钱率减少的时候,应该减少投注的比例。
  3. c 增加时,f^{*} 会减少,意思是说当损失率增加的时候,应该减少投注的比例;反之,当输钱率减少的时候,应该增加投注的比例。
  4. f^{*}\leq 0 时,就表示这个游戏不值得投注,因为总资产的期望是负数。

结论

本文从简单的概率论常识出发,介绍了 Kelly 公式的初步概念。对此有兴趣的读者,建议读其它更有深度的文章。

 

浅谈机器学习的一些小众方向

随着 DeepMind 的 AlphaGo 在 2016 年战胜了李世石,“人工智能”这个词开始进入大众的视野。从那时起,不管是大型互联网公司还是初创企业都开始大规模招聘机器学习的相关从业者,无论社招的求职者还是校招的应聘学生都出现了大规模的增长。由于机器学习的人才短缺并且大量应届生涌入,以至于现在某些公司的校园招聘出现了算法工程师简历太多,并且移动端岗位,web 开发岗位的简历略有不足的情况,导致这些互联网公司甚至通过邮件的方式来呼吁应届生尽量修改投递职位。

字节跳动

就这几年的人工智能发展情况和笔者的个人经验而言,人工智能可以大致分成以下几个方向:

  1. 计算机视觉方向
  2. 自然语言处理方向
  3. 语音识别方向
  4. 机器学习方向

CV,NLP & Speech Recognition

计算机视觉方向(Computer Vision)无论是在学校还是在公司,都有着大量的从业者,并且 ImageNet 项目可以提供上千万的标注图片供大家使用。既然 ImageNet 是开源的数据集,那么无论是学校的教授还是学生,不管是大型互联网公司还是初创企业,都可以轻易地获取到这些数据集,不仅可以进行 CV 算法的研究工作,还可以进行相关的工程实践。由于计算机视觉方向的历史悠久,不管是计算机系,工程系,甚至数学系,都有着大量的老师和相应的学生从事该方向的研究工作,因此,学校或者研究所能够对工业界输出的计算机视觉人才数量也是可观的。

与计算机视觉方向相比,自然语言处理方向(Natural Language Processing)在学校里面也有不少的教授从事相关研究。不过要想让计算机理解人类的语言可不是一件容易的事情。尤其是中文还拥有多音字,语义双关等情形,而且理解中文很可能还要基于上下文来前后推敲。如果和聊天机器人聊过就会发现,其实聊天机器人和人类的聊天给用户的感觉是完全不一样的。语音方向笔者不是很了解,也只是道听途说而已,在这里就不在赘述了。

ImageNet

机器学习

除了以上三个方向,人工智能的另外一个研究方向自然就是机器学习了。在周志华老师的教材《机器学习》中,无监督学习,有监督学习,半监督学习,强化学习等方向都已经在该教材中进行了详细的解释。貌似几年前强化学习这个方向也是不温不火,但是在 AlphaGo 崛起之后,深度学习和强化学习就已经开始进入了大多数人的视野。随着围棋被攻克之后,德州扑克AI,或者其他的游戏 AI 也被很多学者和大型游戏公司所关注。DeepMind 也在 2017 年开放了星际争霸的研究平台,今年无论是在 Dota2 还是星际争霸上,游戏 AI 相比之前都有了巨大的突破。

starcraft2

除了强化学习之下的游戏 AI 之外,其实机器学习一直在一个领域发挥着巨大的用处,那就是推荐系统。无论是广告推荐,YouTube 视频推荐,甚至今年非常火的抖音 APP,推荐系统在其中的作用都不容忽视。关于推荐系统的书其实有很多,笔者也没有一一读过,不过就近些年的发展状况来看,无论是在学术界还是工业界,从零到一搭建一套推荐系统已经不是壁垒,如何搭建一套结合业务场景的优秀推荐系统才是难题。而推荐系统中常用的各种模型,例如逻辑回归(logistic regression),SVD,ItemCF & UserCF,甚至深度神经网络,在各种开源框架之下(Spark,Tensorflow等),只要拥有足够的计算资源,训练出一个可以使用的模型已经没有太大的难度。难度在于算法工程师如何贴近业务并且理解业务,在此基础上如何使用机器学习算法将内容库里面的优质内容推荐给用户,而不引起用户的反感,点击率如何在合理的范围内进一步提升。搭建一套推荐系统已经不是难题,如何结合多种多样的推荐场景才是关键,怎么结合业务来使用推荐系统则是算法工程师需要思考的问题。

Tensorflow

机器学习+安全业务

就笔者的个人经验来看,推荐系统或者游戏 AI 其实只是机器学习的一个应用领域。既然机器学习能够应用在推荐系统或者游戏 AI 上,那么为何不能够应用在别的领域上呢?

对于一些大型互联网公司而言,推荐系统能够给用户们带来足够优质的体验,游戏 AI 能够帮助玩家提升自己的技艺。但是在给用户带来优质体验的时候,总有一些黑产用户在伺机而动,通过 APP 的各种 bug 来寻找赚钱的机会,给正常用户带来各种各样的骚扰。在游戏中,有一些人使用了外挂等技术,破坏了游戏中的平衡。在金融行业中,一直都有黑产用户正在进行各种各样违法犯罪的事情,例如信用卡欺诈等,给正常用户带来了不少的损失。在社交网络中,有一些用户通过社交网络传播着各种各样的不良信息,无论是谣言,虚假广告还是各种各样的假冒伪劣产品宣传,都给正常用户带来了不好的体验。因此,安全业务一直是互联网公司和金融公司的重点业务,安全业务一直是保护着互联网公司能够正常运行的基石。各种各样的安全实验室在大型互联网公司里面并不罕见,也是必须要配备的力量。对于业务安全上,无论是盗号,刷帖,传播虚假消息等都是需要关注的对象。在黑产力量日益壮大的情况下,打击黑产的人力也越来越多。随着人力的增多,如何使用机器学习算法来进行人类经验的传承,或者说随着黑产技术的升级如何才能够尽快的提升互联网公司的黑产对抗能力,这些都是值得做的工作。除了互联网公司之外,银行等金融机构也需要进行信用卡的风控评级,打击信用卡盗刷,黑色产业的资金链条挖掘等。因此,银行等金融机构对于业务安全上面的要求有的时候可能比互联网公司还要严格。

黑客图片1

能够用在安全领域上的机器学习算法有很多,最容易想到的当然就是异常检测。无论是高维异常检测,还是图(Graph)上的异常检测,都在业务安全领域有着巨大的应用场景。异常检测算法可以从众多的数据中发现数据中的异常点,然后通过人工审核等方式进行数据的标注,并且可以使用有监督学习模型进行训练和上线预测。整体来说,就是使用无监督算法,有监督算法,图挖掘算法等机器学习常见技术来进行恶意黑产的打击工作。对于从事业务安全+机器学习方向的算法工程师来说有一些潜在的优势,那就是业务安全方向是工业界的刚需。但是学术界并不完全有能力培养相关的人才,因为互联网或者金融公司的数据都具有保密性,很难把数据像 ImageNet 一样开放给全世界,共同享受数据带来的巨大优势。如果没有基础的数据,那么学校的教授或者学生就无法接触到这个领域,也就无法在学校提升相关的技术。虽然异常检测等其他机器学习算法会在学术中有所突破,但是安全的业务经验只有在做过相关业务之后,真正地打击过黑产用户之后才能够有更深层次的体会和理解。一个没有接触过安全业务的人,即使他的学术造诣再高,在短时间内也是很难提出一些靠谱想法或者技术方案的。

机器学习+运维业务

在这里做一个不恰当的比喻来方便大家理解。

如果把 APP 比喻成一栋楼房的话,那么后台开发就是搭建钢筋水泥的人,前台开发就是负责刷墙贴砖的人,设计师是负责把这栋楼设计得更加美观的人,安全人员就好比楼房的保卫人员,那么运维人员就是这栋大楼的检修人员。

在一些互联网公司,运维人员也被称为技术运营人员,整体来说就是保障APP或者业务稳定运营的。例如:网络抖动了该怎么办,交换机何时宕机,大量用户无法登陆APP了该怎么办,APP的某个页面无法打开了该怎么办等诸如此类的问题。为了保障业务的稳定运营,就需要有一定数量的技术运营同事来维护整个业务的正常运行。正所谓“天有不测风云,人有旦夕祸福”,公司拥有安全人员和运维人员好比买保险,在没有黑客攻击或者业务正常运行的时候,通常存在感略低。但是一旦业务出了问题,第一个要召集的人肯定就是安全和运维人员。因此,无论是安全工作还是运维工作,都是大型互联网公司和金融机构必不可少的力量。

随着机器学习的发展,智能运维(Artificial Intelligence Operations),也就是所谓的 AIOps,也开始被众多技术公司所关注。提到技术运营工作,根据 2018 年的《企业级AIOps实施建议白皮书V0.6》 的观点,可以大致分成以下三个方向:

  1. 质量保障;
  2. 效率提升;
  3. 成本管理。

其中质量保障就是为了保证业务的正常,高效,稳定地运转。在质量保障的过程中,无法避免的就需要进行异常检测。在运维领域,异常检测的范围非常广,不仅包括大家耳熟能详的时间序列异常检测,还包括多维数据下钻分析,甚至还包括日志模板提取和异常挖掘。除了质量保障之外,效率提升也是一个方面,无论是自动化运维领域还是使用 NLP 的技术来构建智能聊天机器人,甚至使用机器学习等技术来进行智能扩缩容,机器学习技术在运维领域都有着巨大的发挥空间。

AIOps场景

在智能运维领域,最重要的任务之一就是时间序列异常检测,这里的时间序列不仅包括服务器的各种各样的指标(CPU,进程,PKG等),还有网络出入流量等交换机数据,甚至包括各种各样的业务指标(在线用户数,失败数,请求量等)。各种各样的时间序列组合在一起就形成了一个时间序列数据库,而且这些时间序列通常来说都是按照分钟量级来收集数据,因此,时间序列项目完全符合机器学习项目的各种条件。在时间序列异常检测或者趋势预测中,时间序列和机器学习,甚至深度学习结合的各种技术都可以在这里有着一定的用武之地。

timeseries

除了时间序列之外,服务器的异常挖掘,多维度数据分析都是智能运维中非常有挑战的项目。除了质量保障之外,效率提升中的智能聊天机器人将有希望把运维人员从繁重的客服任务中解放出来,智能扩缩容技术将有机会取代原来很多“拍脑袋”所做出来的容量估计。对于一家正常经营的公司而言,质量保障和效率提升只是其中的两个方面,如何有效地进行成本的管理则是非常重要的项目。如果成本预算过少,那么明年的项目发展将会受到限制;如果成本预算过多,那么明年的资源势必造成各种浪费。因此,无论是质量保障,效率提升,还是成本管理,都是技术运营领域的核心问题。

成本

机器学习+其他领域

除了以上笔者接触过或者略微了解过的领域之外,其实机器学习在其他的领域应该都是有着自己的用武之地。在量化分析方向,据说有的团队已经开始用机器学习的方法进行股票交易。在化学或者生物学领域,也有学者使用机器学习的方法来挖掘数据之间的信息。总之,除了人工智能在那几个经典领域的应用之外,机器学习的方法应该有希望应用到各行各业中,改变原来的工作方式,提升原有学科的效率。机器学习本身并不是一个新的东西,只要运用得当,机器学习在各行各业都有着强大的创造力和生命力。

 

25 岁做什么,可在 5 年后受益匪浅?

25 岁做什么,可在 5 年后受益匪浅?

很久之前在知乎上看到一个问题:“25 岁做什么,可在 5 年后受益匪浅?

25 岁的时候

在写这篇文章之前先回忆一下自己在 25 岁的时候在干什么。

笔者在 25 岁的时候应该是 2013 年,正好是在 NUS 读博士的时候。当时笔者的科研进展缓慢,几乎处于无法自我推进的状态。而在笔者迷茫了大半年之后,碰巧在网上看到一本书叫做《战胜拖拉》,于是笔者花了几天功夫读完了整本书之后就将知识用于实战,目标是战胜自己长期拖延不科研的状态。不过花费了一段时间之后,效果比较明显,没花费多少时间就把当时论文里面的第一步 Real Bound Theorem 搞出来了。事后想起这件事情的时候,最感谢的就是《战胜拖拉》那本书的作者,作者在书中确实提供了不少有用的建议。在解决了毕业论文里面的重大难题之后,笔者写了两篇关于拖延症的文章,希望对大家有帮助。

PHD 身边的时间陷阱

战胜拖延-让PHD达成每天必要的工作时间

在看完拖延症的书籍并做完论文的第一步之后,当然要完成的就是论文的全部。在做科研的过程中整体来说还是比较辛苦,以至于读完博士之后还写了一篇文章来介绍科研整件事情,那就是”科研这条路“。

选择比努力更重要

所有的博士生在拿到博士学位之后自然就是面临就业的问题。笔者 2005 – 2015 年一直在数学系攻读学位,回顾读书这十年给笔者最大的感受就是,学校招聘老师的标准是越来越高。刚开始的时候,如果博士的学校较好,也许没有太多的论文都能够进一个还不错的大学。到了 2010 年之后,即使有了不少的论文也不能够保证一定能够进好大学。同时,博士生的数量也是越来越多,在学校的教师岗位根本无法容纳下那么多博士生的情况下,很多无法找到教职的博士生就要去企业工作或者继续从事博士后的岗位。不过有的专业找工作确实相对容易一些,有的专业找工作相对难一些。假设有一个博士生在 2010 年选择机器学习专业,那么在他毕业的时候,正好就是机器学习刚刚火起来的时候,那么肯定会非常容易就业。不过在 2010 年的时候,机器学习也不算什么热门方向,企业所提供的岗位也相对较少。因此,有的时候要想在未来获得更大的收益,选择当前热门的领域并不是一件很好的事情,选择未来有可能火的领域才是关键。不过要想判断未来哪个领域会火实在是太难了。因此,选择比努力重要的多,无论是整个大的行业,还是研究方向,甚至到每天手上所做的那件事情。

自我的成长

在选择了一个合适自己的方向之后,那就要去努力做这件事情,无论是在上学的时候,还是工作的时候,都要时刻注意自己的产出。在学校的时候,学生需要注意的就是这半年需要做什么,能够达到的目标是什么。然后从后往前反推此刻应该做什么,需要做什么样的事情才能够达到相应的结果。比如,如果在开学的时候就想要在期末考试的时候获得一个不错的成绩,那么在开学的时候不仅要下定决心学好这门课,还要根据课程的实际情况制定出相应的学习计划,最后才能够在期末考试中获得相应的成绩。又比如,如果目标是半年后写一篇论文,那么就需要准备开题报告,论文规划和预期效果,以及相应的时间节点。但是做论文的时候最大的风险点在于不确定性,所以很多时候需要根据论文的实际情况来进行调整。如果是在公司里面的话,通常来说都是季度考核,半年考核,全年考核等几个大的阶段。而且在项目的不同时期,考核的方式也是截然不同的,在项目的初期,可能也只是需要有一些调研方向和小的产出;在项目的中期,可能就需要有阶段性的成绩;在项目的后期,可能就需要把整个项目进行回顾,然后总结成功点和失败点,再让所有项目组成员来学习成功点,总结失败点,避免在未来的项目中走同样的弯路。其实无论是在学校里面还是在公司里面,“项目管理”这类知识还是挺有用的,在学校的时候可以用来管理自己的学业,在公司的时候可以用来管理项目的进展。

工作的意义不仅仅在于获得一份工资,有可能的话最好从工作中获得一定的自我认同感,更重要的是获得个人的持续成长。刚开始进入公司的时候,由于是新人,所以可以学的东西其实非常多。无论是专业技能,业务沟通,还是项目管理,每一个领域都够一个新人学一段时间。因此,在刚进入公司的时候,需要给自己一定的压力,前一年半其实是个人成长最快的时候,在这个时候最好需要充分利用上下班的时间,甚至周末最好也花一定的时间进行自我充电。而到了两年之后,成长的速度就会明显放缓很多,感觉每隔半年才会成长一点。

专业的问题

在学校的时候,通常学生都只会关注本专业的知识,只有遇到了想转行的时候才会去主动学习别的学科。但是在工作的时候,如果只想做本专业的知识,那么在其实就是在限制自己的发展空间,千万不要让自己的专业限制了自己的发展。无论是在学校还是工作中所学到的经验是一种财富,但是这些经验在有的时候也会形成自己的绊脚石。一般来说,在工作中做一些相对成熟的事情是容易出成绩的,做未知的事情是需要承担风险的。在做成熟事情的时候大家都会采用已有的方案继续做下去,在做未知事情的时候就需要有人去做很多的调研工作,看看这件事情是否值得做下去,是否能够达到预期的收益效果。如果调研了之后确实能够达到预期的效果,那就值得做下去;如果有一定的风险,那就要把预期降低,在一个合理的预期之内做适当的事情。

明确不想做的事情

可以在学校期间明确一下自己不想做的事情。通常来说,你问别人,你想做什么,他能够七七八八的说出一堆,但是绝大部分都不会去做。而且随着时间的迁移,每个人的想法都会产生变化,今年想做的事情明年不一定想做。但是就一般情况而言,一个人不想做的事情是不会发生太大变化的。如果一个人不喜欢学习物理,那么十几年过去之后可能还是这样;如果一个人不喜欢写作,那他肯定也很难提起自己的笔去认认真真地写一篇文章。所以,可以在学校或者工作的时候确定自己不喜欢做什么事情,然后在选择专业或者就业的时候避开这些专业或者岗位。因为在学校里面会与这些专业相伴几年,而工作的时候在岗位的时间有可能比在家的时间还要多。

 

不要为已经打翻的牛奶而哭泣

当年,曹操与刘备在汉中决战,两军久久僵持不下。曹操见久攻不下,心中烦闷,此时士兵来询问夜间的口令,曹操顺口说一句:“鸡肋。”而主簿杨修听到这句话,便开始收拾行装,并告诉周边的人一起收拾行装。众人不解,反问之,杨修解释道:“鸡肋鸡肋,食之无肉,弃之有味,今丞相进不能胜,恐人耻笑,明日必令退兵。”而杨修因为这句话而引来了杀身之祸。

三国1

熟读三国演义的人都知道,曹操杀了杨修之后,便令众军前进。其实到了最后,曹操也未能战胜刘备,获得汉中攻略战的胜利,只能退回许昌。虽然现在我们很难猜测曹操当时的想法是什么,不过“夫鸡肋,弃之如可惜,食之无所得”这句话却由于这个故事而流传下来。无论是当年的曹孟德,还是现在的很多平凡的普通人,身边总有一些事物与鸡肋一样,食之无味而弃之可惜。

在经济学中可以用沉没成本来描述鸡肋这个概念。沉没成本指的是已经付出了,但是不能收回的成本。例如,有一家电影院不允许顾客退票,有位顾客买了一张电影票,但是他看了半小时的电影之后觉得电影十分难看,这种时候他就有两种选择:

1. 继续看下去;

2. 中途离场,去做自己想做的事情。

其实绝大多数的人都会忍着看完这部难看的电影。而从经济学上的观点来看,如果人是足够理性的,当这个人在做决策的时候,是应该把沉没成本放在一边而不去考虑的,因为沉没成本是无法被改变的。按照上面的例子,无论他选择是否继续看下去,这个电影票的钱都已经无法退回,而此时需要做决定的事情是是否继续看完电影。通常这种情况下,按照经济学的理论来说,经济学家会建议这个人选择中途离场,去做任何自己想做的事情。因为这样的话他只是浪费了一张电影票的钱,但是省下来的时间却可以做其他更有趣的事情。

COSTS

虽然这只是书本上的一个简单例子,但是这样的例子在生活中比比皆是。无论是对学生还是职场人士,无论是对普通人还是位高权重的决策者,都面临着沉没成本是否放弃的难题。对于学生而言,最常见的情况就是这个学生在大学期间学了四年根本不感兴趣的专业,但是在面临存在转专业的机会时,是否要放弃原有的技能而重新学习一个新的专业便成为了一个难题。对于在职人士而言,在面临一个有挑战但是并不熟悉的行业或者领域的时候,是否愿意放弃原有的一些经验,是否存在勇气进入一个新的领域也是一个难以抉择的问题。

沉没成本其实很影响一个人的决策。对于学生而言,如果在某个方向上花费了巨大的精力和时间,是很难下定决心转一个全新的方向的;对于工作后的人士而言,主动放弃已经拥有的一些经验,放弃已经掌握的一些人脉和资源,也是十分困难的。但是,在人生的十字路口,其实又必须要下定决心做一些事情。众所周知,学校里面的不少专业其实就是“鸡肋”,完全符合“食之无味,弃之可惜”的条件,无论是学生的就业率和成材率都处于所有专业的底部。如果这些专业的学生不放弃自己的专业,将会在这些专业里面越陷越深,最终无法自拔。其实,这些专业的学生继续从事该专业的学习都不能称之为“坚持”,而是在“死扛”,用自己的大好前途来耗费在一些没有任何用处的专业技能上。有的技能虽然看上去比较高大上,但实在是“屠龙之技”,离开了已有的圈子,学校之后就再无任何用武之地。这些专业还能继续招生的原因大概就是学校招聘了不少这些专业的教师。因此,对于这些专业的学生而言,不破不立,只有勇敢地走出自己所在的圈子,才能够体会到其他专业的精妙之处。

CHOICE

而对于职场人士而言,基本上都会想靠一些好项目来升职加薪,但是在整个社会的大环境下,有的方向确实是在走下坡路,行业越来越饱和,竞争越来越激烈,所做的技术难度越来越低。随着技术的发展,原有的一些技术和框架都会逐渐被淘汰,掌握的技能价值也会越来越低,甚至可能出现找一些应届生培训几个月之后就能够达到老员工的水平。在这种情况下,随着行业的整体下滑,如果还抱着原有的技术栈不松手,那只能变得越来越没有竞争力。在这种情况下,就要主动寻求突破,寻找自己所拥有的技能和其他专业的共同点,将自己的技能主动地迁移到更有潜力的方向上。在职场上,一定不要死抱着一个东西不放手,一定不要抱着我就是来做这个方向的想法,而其他的方向都不去了解和尝试。在工作中,应该审时度势,创造或者寻找优质的项目和资源,因为一个有潜力的项目和优质的资源所能够带来的好处有的时候会远远大于自己当年所做的方向,那个自己当年不舍得放弃的方向。

上升通道

整体来说,在一个人面临着决策的时候,沉没成本最会影响一个人的决策。无论是从经济学的原理上来说,还是从日常经验上来讲,其实都应该放下包袱,轻装上阵,寻找一个更有前景和前途的方向去发挥自己的特长。

 

zr9558's Blog