给定一个数据序列,在某个时间点,数据的某个(或某些)参数可能由于系统性因素(而非偶然性因素)而突然发生变化,那么这个时间点被称为变点(changepoint)。变点检测(changepoint detection)就是要估计出变点的位置。

本文主要打算理一下这篇论文:

Bayesian Online Changepoint Detection. Ryan Prescott Adams and David J.C. MacKay. arXiv 2007. [Paper]

翻译过来大概是贝叶斯在线变点检测?“online”这个词指检测变点时只能利用当前已经观测到的数据,不能用未来数据。而很多 offline 的方法是可以拿所有的数据来做检测的。

本文很大程度上(又)参考了一篇讲 BOCD 讲得非常清楚的博客:

Bayesian Online Changepoint Detection (Gregory Gundersen)

问题定义

假设有一个观测序列 ,我们用 来表示 可以被划分成一些不重叠的 product partitions,即每个 partition 中的数据点都是 i.i.d(独立同分布)的,都服从参数为 的分布 ,每个 partition 的参数 也是 i.i.d 的:

product partitions

而 BOCD 会估计当前时刻的 run length(距离上一个变点已经过了多少时刻)的后验分布 时刻的 run length 为 。显然 只可能有两种情况:

那么一个可能的 的变化图如下( 的时刻()为变点):

run length

除了 以外,BOCD 还会估计 ,这样就能既预测变点(run lenght),又预测下一个数据点的值了。不过如果单纯只求变点的话, 是没必要算的。

一个递推式

如果不想看推导过程,可以直接跳到这里看结果。

下一个数据点的分布

可以由给定 的边缘分布算出来:

其中 为 run length 时当前 partion 的观测点集合,因为我们只会用当前 partion 的数据来预测 的分布(假设 1)。因为 可能为 0,所以 可能为空集。

run length 的分布

那么我们就需要计算 run length 的分布

那么我们需要算联合概率

最后一步能成立依赖于以下两个假设:

  1. 前面提到的假设 1,即只用当前 partion 的数据来预测 的分布:

    因为 只依赖于 (假设 2),所以这里 也可以省掉。

  2. 只依赖于 。即给定 对其他所有变量都条件独立:

递推流程

可以看到求 的过程是一个递推的过程。在 时刻, 已经计算出来了,现在我们要求 ,那么计算流程为:

  1. 计算变点的先验

  2. 计算 run length 的后验分布:

    其中:

  3. 计算新数据点的分布

可以看到,我们还需要计算的两个式子是:

  • :为了方便,我们按照本文参考的博客的叫法,把它称为 Underlying Probabilistic Model(UPM),它的求法将在这一节解释

  • 变点先验 :它的求法将在这一节解释

这张图说明了这个递推流程:

message passing

图片来源:Blog: Bayesian Online Changepoint Detection

当然,在开始这个递推流程之前,我们需要手动定义初始值 ,定义方式将在这一节解释。

变点先验

贝叶斯方法的特点就是能利用先验信息,所以我们可以把我们对变点的先验估计通过 hazard function 塞进模型里。

我们假设:

是 hazard function,描述的是个体在 时刻死亡的概率(...)。当然“死亡”这个词可以换成别的什么特殊事件,比如在这里 描述的就是在 run length 时(在这之前都没出现变点)出现变点的概率。

Survival Function

是一个表示当前 run length 长度的随机变量。 是 survival function,表示 run length 超过 的概率:

其中 表示当前 run length 的概率,是一个概率密度函数。可以看到有:

Hazard Function

那么由 hazard function 的定义:

这里 是我们自己定的一个先验。

正好被定为一个指数分布(或几何分布(离散))时:

此时 是一个常数(但论文里面写的是 ,我也不知道我哪里推错了...)。

递推初始值

我们需要给 赋一个初始值,分为两种情况:

  1. 第一个观测到的数据点之前()就是一个变点。论文中给的例子是对游戏数据序列建模,游戏开始的时刻一定是个变点。那么显然有:

    但这个假设并不总是成立。

  2. 假设我们现在有最近一段时间的历史数据,并且我们认为第一个变点会在未来某个时刻出现。论文给的例子是对气象数据建模。这时初始值为归一化后的 survival function:

    其中 是一个归一化常数。

UPM

好的那么现在 就是最后一个要算的东西了,我们把它写成边缘分布:

表示在 时刻 run length 时的超参数,第一项的意义是 服从参数为 指数族分布(Exponential Family),然后对这个指数族分布建模。第二项是 的后验。

的后验并不好计算,而且这里还要对每个 求积分,也不好算。所以我们往往希望 的先验跟它的似然共轭(conjugate),使得先验与后验的形式相同,这样就能简化计算。

共轭先验

假设 是已经观测到的数据, 是要预测的新数据点, 是模型参数, 是超参数。那么:

这里 的先验跟它的似然是共轭的,所以 的后验跟它的先验的形式相同,只是超参数不同:

那么有:

也就是说 的后验跟它的先验的形式也是相同的。如果我们能算出 的后验的超参数 ,就能在不直接算后验 也不算积分的情况下直接算出

而当 的似然是指数族分布时,就一定能写出其共轭先验分布,且后验的超参数 能够很轻松的求出来。

指数族分布

指数族是一类分布,包括高斯分布、伯努利分布、二项分布、泊松分布、Beta 分布、Dirichlet 分布、Gamma 分布等一系列分布。指数族分布可以写为统一的形式:

其中, 是 underlying measure(没找到啥合适的翻译); 是充分统计量(sufficient statistic),包含样本集合的所有信息,如高斯分布中的均值 和方差 是正则函数(normalizer); 是对数正则函数(log normalizer),用于保证:

叫 log normalizer 这个名字是因为由上式可以推出:

跟似然 共轭的 的先验为:

其中 是超参数, 取决于指数族分布的具体形式。

贝叶斯公式,后验正比于似然 先验,所以:

最后一步是因为 是跟 无关的,所以可以去掉。

可以看到后验 跟先验 的形式是一样的,只是超参数上有区别:

所以如果似然是指数族分布,那么在先验的超参数上加一个项就可以得到后验的超参数,是种很简便的方法。

又一个递推

现在 可以被表示成 ,相当于我们的目标变成了求 ,而且要对每一个可能的 (即 )都算一个

这个递推的过程大概是这样:

parameter updates

图片来源:Blog: Bayesian Online Changepoint Detection

完整算法

有空再说,我先摸下🐟...

参考