BOCD

给定一个数据序列，在某个时间点，数据的某个（或某些）参数可能由于系统性因素（而非偶然性因素）而突然发生变化，那么这个时间点被称为变点（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)

# 问题定义

假设有一个观测序列 $x_1, x_2, \dots, x_T \in \Reals^d$，我们用 $x_{a:b}$ 来表示 $x_a, x_{a+1}, \dots x_b$。$x_{1:T}$ 可以被划分成一些不重叠的 product partitions，即每个 partition 中的数据点都是 i.i.d（独立同分布）的，都服从参数为 $\theta_p$ 的分布 $p(x_t \mid \theta_p)$，每个 partition 的参数 $\theta_p \thicksim p(\phi)$ 也是 i.i.d 的：

而 BOCD 会估计当前时刻的 run length（距离上一个变点已经过了多少时刻）的后验分布 $p(r_t \mid x_{1:t})$，$t$ 时刻的 run length 为 $r_t$。显然 $r_t$ 只可能有两种情况：

$$x = \begin{cases} 0 &t \text{ 时刻是变点} \\ r_{t-1} + 1 &t \text{ 时刻不是变点} \end{cases}$$

那么一个可能的 $r_t$ 随 $t$ 的变化图如下（$r_t = 0$ 的时刻（$t=5, 11$）为变点）：

除了 $p(r_t \mid x_{1:t})$ 以外，BOCD 还会估计 $p(x_{t+1} \mid x_{1:t})$，这样就能既预测变点（run lenght），又预测下一个数据点的值了。不过如果单纯只求变点的话，$p(x_{t+1} \mid x_{1:t})$ 是没必要算的。

# 一个递推式

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

# 下一个数据点的分布

$p(x_{t+1} \mid x_{1:t})$ 可以由给定 $r_t$ 时 $x_{t+1}$ 的边缘分布算出来：

$$\begin{aligned} p(x_{t+1} \mid x_{1:t})&= \sum_{r_t} p(x_{t+1}, r_t \mid x_{1:t}) & \text{(边缘概率公式)} \\ &= \sum_{r_t} p(x_{t+1} \mid r_t, x_{1:t}) p(r_t \mid x_{1:t}) & \text{(链式法则)}\\ &= \sum_{r = 0}^t p(x_{t+1} \mid r_t = r, x_{(t - r):t}) p(r_t \mid x_{1:t}) & \text{(假设 1)}\\ &= \sum_{r_t} p(x_{t+1} \mid r_t, x_t^{(r)}) p(r_t \mid x_{1:t}) & \end{aligned}$$

其中 $x_t^{(r)} = x_{(t - r):t}$ 为 run length $= r$ 时当前 partion 的观测点集合，因为我们只会用当前 partion 的数据来预测 $x_{t+1}$ 的分布（假设 1）。因为 $r_t$ 可能为 0，所以 $x_t^{(r)}$ 可能为空集。

# run length 的分布

那么我们就需要计算 run length 的分布 $p(r_t \mid x_{1:t})$：

$$p(r_t \mid x_{1:t}) = \frac{p(r_t, x_{1:t})}{p(x_{1:t})}$$

那么我们需要算联合概率 $p(r_t, x_{1:t})$：

$$\begin{aligned} p(r_t, x_{1:t})&= \sum_{r_{t-1}} p(r_t, r_{t-1}, x_{1:t}) & \text{(边缘概率公式)} \\ &= \sum_{r_{t-1}} p(r_t, r_{t-1}, x_{1:t-1}, x_t) & \\ &= \sum_{r_{t-1}} p(r_t, x_t \mid r_{t-1}, x_{1:t-1}) p(r_{t-1}, x_{1:t-1}) & \text{(链式法则)} \\ &= \sum_{r_{t-1}} p(x_t \mid \cancel{r_t}, r_{t-1}, x_{1:t-1}) p(r_t \mid r_{t-1}, \cancel{x_{1:t-1}}) p(r_{t-1}, x_{1:t-1}) & \text{(链式法则)} \\ &= \sum_{r_{t-1}} p(x_t \mid r_{t-1}, x^{(r)}) p(r_t \mid r_{t-1}) p(r_{t-1}, x_{1:t-1}) & \text{(假设 1, 2)} \end{aligned}$$

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

1. 前面提到的假设 1，即只用当前 partion 的数据来预测 $x_t$ 的分布：

   $$p(x_t \mid r_t, r_{t-1}, x_{1:t-1}) = p(x_t \mid r_{t-1}, x_{t-1}^{(r)}) \tag{\text{假设 1}}$$

   因为 $r_t$ 只依赖于 $r_{t-1}$（假设 2），所以这里 $r_t$ 也可以省掉。

2. $r_t$ 只依赖于 $r_{t-1}$。即给定 $r_{t-1}$，$r_t$ 对其他所有变量都条件独立：

   $$p(r_t \mid r_{t-1}, x_{1:t-1}) = p(r_t \mid r_{t-1}) \tag{\text{假设 2}}$$

# 递推流程

可以看到求 $p(x_{t+1} \mid x_{1:t})$ 的过程是一个递推的过程。在 $t$ 时刻，$p(r_{t-1}, x_{1:t-1})$ 已经计算出来了，现在我们要求 $p(x_{t+1} \mid r_t, x_t^{(r)})$，那么计算流程为：

1. 计算变点的先验 $p(r_t \mid r_{t-1})$

2. 计算 run length $r_t$ 的后验分布：$p(r_t \mid x_{1:t}) = \frac{p(r_t, x_{1:t})}{p(x_{1:t})}$

   其中：

   $$p(r_t, x_{1:t}) = \sum_{r_{t-1}} \underbrace{p(x_t \mid r_{t-1}, x_{t-1}^{(r)})}_{\text{UPM}} \underbrace{p(r_t \mid r_{t-1})}_{\text{变点先验}} \underbrace{p(r_{t-1}, x_{1:t-1})}_{\text{已求得的信息}} \\[1pt]$$

3. 计算新数据点的分布 $p(x_{t+1} \mid x_{1:t})$：

   $$p(x_{t+1} \mid x_{1:t}) = \sum_{r_t} \underbrace{p(x_{t+1} \mid r_t, x_t^{(r)})}_{\text{UPM}} \underbrace{p(r_t \mid x_{1:t})}_{\text{已求得的信息}} \\[1pt]$$

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

• $p(x_t \mid r_{t-1}, x_{t-1}^{(r)})$ 和 $p(x_{t+1} \mid r_t, x_t^{(r)})$：为了方便，我们按照本文抄参考的博客的叫法，把它称为 Underlying Probabilistic Model（UPM），它的求法将在这一节解释

• 变点先验 $p(r_t \mid r_{t-1})$：它的求法将在这一节解释

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

图片来源：Blog: Bayesian Online Changepoint Detection

当然，在开始这个递推流程之前，我们需要手动定义初始值 $p(r_0, x = \empty) = p(r_0)$，定义方式将在这一节解释。

# 变点先验

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

我们假设：

$$p(r_t \mid r_{t-1}) = \begin{cases} H(r_{t-1} + 1) &\text{if } r_t = 0 \\ 1 - H(r_{t-1} + 1) &\text{if } r_t = r_{t-1} + 1 \\ 0 &\text{otherwise} \end{cases}$$

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

# Survival Function

设 $T \geq 0$ 是一个表示当前 run length 长度的随机变量。$S(\tau)$ 是 survival function，表示 run length 超过 $\tau$ 的概率：

$$S(\tau) = P(T \geq \tau) = 1 - F(\tau) = \sum_{\tau' = \tau}^\infty f(\tau')$$

其中 $F(\tau) = P(T < \tau)$。$f(\tau)$ 表示当前 run length $= \tau$ 的概率，是一个概率密度函数。可以看到有：

$$f(\tau) = F'(\tau) = -\frac{d S(\tau)}{d \tau}$$

# Hazard Function

那么由 hazard function 的定义：

$$\begin{aligned} H(\tau) &= \lim_{\Delta \tau \rarr 0} \frac{p(\tau \leq T < \tau + \Delta \tau \mid T \geq \tau)}{\Delta \tau} \\ &= \lim_{\Delta \tau \rarr 0} \frac{p(\tau \leq T < \tau + \Delta \tau)}{p(T \geq \tau) \Delta \tau} \\ &= \frac{f(\tau)}{S(\tau)} \end{aligned}$$

这里 $f(\tau)$ 是我们自己定的一个先验。

当 $f(\tau)$ 正好被定为一个指数分布（或几何分布（离散））时：

$$f(\tau) = \lambda e^{- \lambda \tau}$$

$$F(\tau) = 1- e^{- \lambda \tau}$$

$$S(\tau) = 1- F(\tau) = e^{- \lambda \tau}$$

此时 $H(\tau) = \frac{f(\tau)}{S(\tau)} = \lambda$ 是一个常数（但论文里面写的是 $\frac{1}{\lambda}$，我也不知道我哪里推错了...）。

# 递推初始值

我们需要给 $p(r_0)$ 赋一个初始值，分为两种情况：

1. 第一个观测到的数据点之前（$t=0$）就是一个变点。论文中给的例子是对游戏数据序列建模，游戏开始的时刻一定是个变点。那么显然有：

   $$p(r_0 = 0) = 1$$

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

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

   $$p(r_0 = \tau) = \frac{1}{Z} \sum_{\tau' = \tau + 1}^\infty f(\tau')$$

   其中 $Z$ 是一个归一化常数。

# UPM

好的那么现在 $p(x_{t+1} \mid r_t, x_t^{(r)})$ 就是最后一个要算的东西了，我们把它写成边缘分布：

$$p(x_{t+1} \mid r_t, x_t^{(r)}) = \int \underbrace{p(x_{t+1} \mid \eta)}_{\text{EF model}} \underbrace{p(\eta_t^{(r)} = \eta \mid r_t, x_t^{(r)})}_{\text{EF posterior}} d \eta \\[1pt]$$

$\eta_t^{(r)}$ 表示在 $t$ 时刻 run length $= r$ 时的超参数，第一项的意义是 $x$ 服从参数为 $\eta$ 的指数族分布（Exponential Family），然后对这个指数族分布建模。第二项是 $\eta$ 的后验。

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

# 共轭先验

假设 $X$ 是已经观测到的数据，$\hat{x}$ 是要预测的新数据点，$\eta$ 是模型参数，$\alpha$ 是超参数。那么：

$$p(\hat{x} \mid X, \alpha) = \int \underbrace{p(\hat{x} \mid \eta)}_{\text{model}} \underbrace{p(\eta \mid X, \alpha)}_{\text{posterior}} d \eta$$

这里 $\eta$ 的先验跟它的似然是共轭的，所以 $\eta$ 的后验跟它的先验的形式相同，只是超参数不同：

$$p(\eta \mid X, \alpha) = p(\eta \mid \alpha')$$

那么有：

$$\begin{aligned} p(\hat{x} \mid X, \alpha) &= \int p(\hat{x} \mid \eta) p(\eta \mid X, \alpha) d \eta \\ &= \int p(\hat{x} \mid \eta) p(\eta \mid \alpha') d \eta \\ &= p(\hat{x} \mid \alpha') \end{aligned}$$

也就是说 $\hat{x}$ 的后验跟它的先验的形式也是相同的。如果我们能算出 $\eta$ 的后验的超参数 $\alpha'$，就能在不直接算后验 $p(\eta \mid X, \alpha)$ 也不算积分的情况下直接算出 $p(\hat{x} \mid X, \alpha)$。

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

# 指数族分布

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

$$\begin{aligned} p(x \mid \eta) &= h(x) \exp (\eta^\top U(x) - A(\eta)) \\ &= \frac{1}{\exp(A(\eta))} h(x) \exp (\eta^\top U(x)) \\ &= g(\eta) h(x) \exp (\eta^\top U(x)) \end{aligned}$$

其中，$h(x)$ 是 underlying measure（没找到啥合适的翻译）；$U(x)$ 是充分统计量（sufficient statistic），包含样本集合的所有信息，如高斯分布中的均值 $\mu$ 和方差 $\sigma$；$g(\eta) = \frac{1}{\exp(A(\eta))}$ 是正则函数（normalizer）；$A(\eta)$ 是对数正则函数（log normalizer），用于保证：

$$\frac{1}{\exp(A(\eta))} \int h(x) \exp (\eta^\top U(x)) dx = 1$$

$A(\eta)$ 叫 log normalizer 这个名字是因为由上式可以推出：

$$A(\eta) = \log \int h(x) \exp (\eta^\top U(x)) dx$$

跟似然 $p(x \mid \eta)$ 共轭的 $\eta$ 的先验为：

$$p(\eta \mid \chi, \nu) = f(\chi, \nu) g(\eta)^{\nu} \exp ( \eta^{\top} \chi )$$

其中 $\chi, \nu$ 是超参数，$f(\chi, \nu)$ 取决于指数族分布的具体形式。

由贝叶斯公式，后验正比于似然 $\times$ 先验，所以：

$$\begin{aligned} p(\eta \mid X, \chi, \nu) & \propto p(X \mid \eta) p(\eta \mid \chi, \nu) \\ & \propto \underbrace{\vphantom{\Bigg|} \left ( \Big ( \prod_{i=1}^N h(x_n) \Big ) g(\eta)^N \exp (\eta^\top \sum_{n=1}^N u(x_n)) \right )}_{\text{似然}} \underbrace{\vphantom{\Bigg|} \Big ( f(\chi, \nu) g(\eta)^{\nu} \exp ( \eta^{\top} \chi ) \Big )}_{\text{先验}} \\ & \propto \Big ( \prod_{i=1}^N h(x_n) \Big ) f(\chi, \nu) g(\eta)^{N + \nu} \exp \Big (\eta^\top \sum_{n=1}^N u(x_n) + \eta^{\top} \chi \Big ) \\ & \propto g(\eta)^{N + \nu} \exp \Big (\eta^\top \sum_{n=1}^N u(x_n) + \eta^{\top} \chi \Big ) \end{aligned}$$

最后一步是因为 $\Big ( \prod_{i=1}^N h(x_n) \Big ) f(\chi, \nu)$ 是跟 $\eta$ 无关的，所以可以去掉。

可以看到后验 $p(\eta \mid X, \chi, \nu)$ 跟先验 $p(\eta \mid \chi, \nu)$ 的形式是一样的，只是超参数上有区别：

$$\nu' = \nu_{\text{prior}} + N$$

$$\chi' = \chi_{\text{prior}} + \sum_{n=1}^N u(x_n)$$

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

# 又一个递推

现在 $p(x_{t+1} \mid r_t, x_t^{(r)})$ 可以被表示成 $p(x_{t+1} \mid \nu_t^{(r)}, \chi_t^{(r)})$，相当于我们的目标变成了求 $\nu_t^{(r)}$ 和 $\chi_t^{(r)}$，而且要对每一个可能的 $r_t = r \leq t$（即 $X = x_t^{(r)} = x_{(t-r): t}$）都算一个 $\nu_t^{(r)}$ 和 $\chi_t^{(r)}$：

$$\nu_t^{(r)} = \nu_{t-1}^{(r-1)} + 1$$

$$\chi_t^{(r)} = \chi_{t-1}^{(r-1)} + u(x_t)$$

$$\nu_t^0 = \nu_{\text{prior}}$$

$$\chi_t^0 = \chi_{\text{prior}}$$

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

图片来源：Blog: Bayesian Online Changepoint Detection

# 完整算法

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

# 参考

• Bayesian Online Changepoint Detection. Ryan Prescott Adams and David J.C. MacKay. arXiv 2007. (original paper)

• Bayesian Online Changepoint Detection (a blog by Gregory Gundersen)

• Bayesian online change point detection — An intuitive understanding (a blog on Medium)

• Survival Analysis: Techniques for Censored and Truncated Data. John P. Klein and Melvin L. Moeschberger. Springer Science & Business Media, 2006.

• 指数族分布