Machine Learning Math Bayesian

BOCD

Bayesian Online Changepoint Detection

Renovamen
2020-08-22
17 min

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

# 问题定义

假设有一个观测序列 x1,x2,,xTRdx_1, x_2, \dots, x_T \in \Reals^d,我们用 xa:bx_{a:b} 来表示 xa,xa+1,xbx_a, x_{a+1}, \dots x_bx1:Tx_{1:T} 可以被划分成一些不重叠的 product partitions,即每个 partition 中的数据点都是 i.i.d(独立同分布)的,都服从参数为 θp\theta_p 的分布 p(xtθp)p(x_t \mid \theta_p),每个 partition 的参数 θpp(ϕ)\theta_p \thicksim p(\phi) 也是 i.i.d 的:

product partitions

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

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

那么一个可能的 rtr_ttt 的变化图如下(rt=0r_t = 0 的时刻(t=5,11t=5, 11)为变点):

run length

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

# 一个递推式

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

# 下一个数据点的分布

p(xt+1x1:t)p(x_{t+1} \mid x_{1:t}) 可以由给定 rtr_txt+1x_{t+1} 的边缘分布算出来:

p(xt+1x1:t)=rtp(xt+1,rtx1:t)(边缘概率公式)=rtp(xt+1rt,x1:t)p(rtx1:t)(链式法则)=r=0tp(xt+1rt=r,x(tr):t)p(rtx1:t)(假设 1)=rtp(xt+1rt,xt(r))p(rtx1:t)\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}

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

# run length 的分布

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

p(rtx1:t)=p(rt,x1:t)p(x1:t)p(r_t \mid x_{1:t}) = \frac{p(r_t, x_{1:t})}{p(x_{1:t})}

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

p(rt,x1:t)=rt1p(rt,rt1,x1:t)(边缘概率公式)=rt1p(rt,rt1,x1:t1,xt)=rt1p(rt,xtrt1,x1:t1)p(rt1,x1:t1)(链式法则)=rt1p(xtrt,rt1,x1:t1)p(rtrt1,x1:t1)p(rt1,x1:t1)(链式法则)=rt1p(xtrt1,x(r))p(rtrt1)p(rt1,x1:t1)(假设 1, 2)\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 的数据来预测 xtx_t 的分布:

    p(xtrt,rt1,x1:t1)=p(xtrt1,xt1(r))(假设 1)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}}

    因为 rtr_t 只依赖于 rt1r_{t-1}(假设 2),所以这里 rtr_t 也可以省掉。

  2. rtr_t 只依赖于 rt1r_{t-1}。即给定 rt1r_{t-1}rtr_t 对其他所有变量都条件独立:

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

# 递推流程

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

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

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

    其中:

    p(rt,x1:t)=rt1p(xtrt1,xt1(r))UPMp(rtrt1)变点先验p(rt1,x1:t1)已求得的信息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(xt+1x1:t)p(x_{t+1} \mid x_{1:t})

    p(xt+1x1:t)=rtp(xt+1rt,xt(r))UPMp(rtx1: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(xtrt1,xt1(r))p(x_t \mid r_{t-1}, x_{t-1}^{(r)})p(xt+1rt,xt(r))p(x_{t+1} \mid r_t, x_t^{(r)}):为了方便,我们按照本文参考的博客的叫法,把它称为 Underlying Probabilistic Model(UPM),它的求法将在这一节解释

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

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

message passing

图片来源:Blog: Bayesian Online Changepoint Detection

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

# 变点先验

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

我们假设:

p(rtrt1)={H(rt1+1)if rt=01H(rt1+1)if rt=rt1+10otherwisep(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(τ)H(\tau) 是 hazard function,描述的是个体在 τ\tau 时刻死亡的概率(...)。当然“死亡”这个词可以换成别的什么特殊事件,比如在这里 H(τ)H(\tau) 描述的就是在 run length =τ= \tau 时(在这之前都没出现变点)出现变点的概率。

# Survival Function

T0T \geq 0 是一个表示当前 run length 长度的随机变量。S(τ)S(\tau) 是 survival function,表示 run length 超过 τ\tau 的概率:

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

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

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

# Hazard Function

那么由 hazard function 的定义:

H(τ)=limΔτ0p(τT<τ+ΔτTτ)Δτ=limΔτ0p(τT<τ+Δτ)p(Tτ)Δτ=f(τ)S(τ)\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(τ)f(\tau) 是我们自己定的一个先验。

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

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

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

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

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

# 递推初始值

我们需要给 p(r0)p(r_0) 赋一个初始值,分为两种情况:

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

    p(r0=0)=1p(r_0 = 0) = 1

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

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

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

    其中 ZZ 是一个归一化常数。

# UPM

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

p(xt+1rt,xt(r))=p(xt+1η)EF modelp(ηt(r)=ηrt,xt(r))EF posteriordη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]

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

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

# 共轭先验

假设 XX 是已经观测到的数据,x^\hat{x} 是要预测的新数据点,η\eta 是模型参数,α\alpha 是超参数。那么:

p(x^X,α)=p(x^η)modelp(ηX,α)posteriordη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(ηX,α)=p(ηα)p(\eta \mid X, \alpha) = p(\eta \mid \alpha')

那么有:

p(x^X,α)=p(x^η)p(ηX,α)dη=p(x^η)p(ηα)dη=p(x^α)\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}

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

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

# 指数族分布

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

p(xη)=h(x)exp(ηU(x)A(η))=1exp(A(η))h(x)exp(ηU(x))=g(η)h(x)exp(ηU(x))\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)h(x) 是 underlying measure(没找到啥合适的翻译);U(x)U(x) 是充分统计量(sufficient statistic),包含样本集合的所有信息,如高斯分布中的均值 μ\mu 和方差 σ\sigmag(η)=1exp(A(η))g(\eta) = \frac{1}{\exp(A(\eta))} 是正则函数(normalizer);A(η)A(\eta) 是对数正则函数(log normalizer),用于保证:

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

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

A(η)=logh(x)exp(ηU(x))dxA(\eta) = \log \int h(x) \exp (\eta^\top U(x)) dx

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

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

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

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

p(ηX,χ,ν)p(Xη)p(ηχ,ν)((i=1Nh(xn))g(η)Nexp(ηn=1Nu(xn)))似然(f(χ,ν)g(η)νexp(ηχ))先验(i=1Nh(xn))f(χ,ν)g(η)N+νexp(ηn=1Nu(xn)+ηχ)g(η)N+νexp(ηn=1Nu(xn)+ηχ)\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}

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

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

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

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

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

# 又一个递推

现在 p(xt+1rt,xt(r))p(x_{t+1} \mid r_t, x_t^{(r)}) 可以被表示成 p(xt+1νt(r),χt(r))p(x_{t+1} \mid \nu_t^{(r)}, \chi_t^{(r)}),相当于我们的目标变成了求 νt(r)\nu_t^{(r)}χt(r)\chi_t^{(r)},而且要对每一个可能的 rt=rtr_t = r \leq t(即 X=xt(r)=x(tr):tX = x_t^{(r)} = x_{(t-r): t})都算一个 νt(r)\nu_t^{(r)}χt(r)\chi_t^{(r)}

νt(r)=νt1(r1)+1\nu_t^{(r)} = \nu_{t-1}^{(r-1)} + 1

χt(r)=χt1(r1)+u(xt)\chi_t^{(r)} = \chi_{t-1}^{(r-1)} + u(x_t)

νt0=νprior\nu_t^0 = \nu_{\text{prior}}

χt0=χprior\chi_t^0 = \chi_{\text{prior}}

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

parameter updates

图片来源:Blog: Bayesian Online Changepoint Detection

# 完整算法

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

# 参考