预备知识

马尔科夫性质

马尔可夫性质(Markov property)是概率论中的一个概念, 因为俄国数学家安德雷·马尔可夫得名。
当一个随机过程在给定现在状态及所有过去状态情况下, 其未来状态的条件概率分布仅依赖于当前状态；
换句话说, 在给定现在状态时, 它与过去状态(即该过程的历史路径)是条件独立的, 那么此随机过程即具有马尔可夫性质。
具有马尔可夫性质的过程通常称之为马尔可夫过程。
数学上, 如果X(t), t > 0为一个随机过程, 则马尔科夫性质就是指

P(X(t+h) = y | X(s) = x(s), s <= t) = P(X(t+h) = y | X(t) = x(t)), h > 0

马尔科夫链

马尔科夫链是满足马尔科夫性质的随机变量的序列$ X_1 $, $ X_2 $, …, 即给出当前状态，将来状态和过去状态是相互独立的。
从形式上看，如果两边的条件分布有定义(即如果$ P(X_1 = x_1, …, X_n = x_n) > 0 $), 则

$ P(X_{n+1} = x | X_1 = x_1, X_2 = x_2, …, X_n = x_n) = P(X_{n+1} = x | X_n = x_n) $

$ X_i $的可能值构成的可数集S叫做该链的”状态空间”。

m-阶马尔科夫链(记忆为m的马尔科夫链), 其中m有限, 满足

$ P(X_n = x_n | X_{n-1} = x_{n-1}, X_{n-2} = x_{n-2}, …, X_1 = x_1)
= P(X_n = x_n | X_{n-1} = x_{n-1}, X_{n-2} = x_{n-2}, …, X_{n-m} = x_{n-m}), n > m
$

的过程。换句话说, 未来状态取决于其前m个状态。

隐马尔科夫模型

隐马尔可夫模型(Hidden Markov Model, HMM)是统计模型，它用来描述一个含有隐含未知参数的马尔可夫过程。
其难点是从可观察的参数中确定该过程的隐含参数。

x — 隐含状态
y — 可观察的输出
a — 转换概率(transition probabilities)
b — 输出概率(output probabilities)

上图展示了HMM的演化过程，图中箭头方向表示不同状态的关联性，x(t)只由x(t-1)决定，y(t)只由x(t)决定。

HMM有三个典型的问题

1. 已知模型参数，计算某一特定输出序列的概率.通常使用forward算法解决.
2. 已知模型参数，寻找最可能的能产生某一特定输出序列的隐含状态的序列.通常使用Viterbi算法解决.
3. 已知输出序列，寻找最可能的状态转移以及输出概率.通常使用Baum-Welch算法以及Reversed Viterbi算法解决.

维特比算法

维特比算法(Viterbi algorithm)是一种动态规划算法。它用于寻找最有可能产生观测事件序列的维特比路径–隐含状态序列。

假设给定隐式马尔可夫模型状态空间 S，初始状态 i 的概率为 $ \pi_{i} $, 从状态 i 到状态 j 的转移概率为 $ a_{i,j} $.
令观察到的输出为 $ y_{1}, \dots, y_{T}$。 产生观察结果的最有可能的状态序列 $ x_{1}, \dots, x_{T} $ 由递推关系给出：

$ V_{1,k} = P{\big (}y_{1}\ |\ k{\big )} \cdot \pi_{k} $

$ V_{t,k} = P{\big (}y_{t}\ |\ k{\big )} \cdot \max_{x\in S}\left(a_{x,k} \cdot V_{t-1,x}\right) $

此处 $V_{t,k}$ 是前 t 个最终状态为 k 的观测结果最有可能对应的状态序列的概率。
通过保存向后指针记住在第二个等式中用到的状态 x 可以获得维特比路径。
声明一个函数 Ptr(k,t), 当t = 1时返回k, 当t > 1, 返回$V_{t,k}$. 这样:

$ x_{T} = \arg \max_{x\in S}(V_{T,x}) $

$ x_{t} = Ptr(x_{t},t) $

这里我们使用了arg max的标准定义。
算法复杂度为

$ O(T\times \left|{S}\right|^{2}) $

一个例子

想象一个乡村诊所。村民有着非常理想化的特性，要么健康要么发烧。他们只有问诊所的医生的才能知道是否发烧。 聪明的医生通过询问病人的感觉诊断他们是否发烧。村民只回答他们感觉正常、头晕或冷。

假设一个病人每天来到诊所并告诉医生他的感觉。医生相信病人的健康状况如同一个离散马尔可夫链。病人的状态有两种“健康”和“发烧”，但医生不能直接观察到，这意味着状态对他是“隐含”的。每天病人会告诉医生自己有以下几种由他的健康状态决定的感觉的一种：正常、冷或头晕。这些是观察结果。 整个系统为一个隐马尔可夫模型(HMM)。

我们假设观察结果序列为{“Normal”, “Cold”, “Dizzy”}, 那么计算过程如下：

由图可见，最有可能由状态序列{“Healthy”, “Healthy”, “Fever”}产生。

突发事件检测

Introduction

微博空间充斥着海量的短信息,用户无法通过阅读大量微博来获取实时的突发话题,对微博中的海量信息进行挖掘,发现实时发生的突发事件可以有效地帮助用户找到自己感兴趣的话题,改善用户体验; 另一方面,随着近年来的迅猛发展,微博在舆情分析、观点挖掘等领域起到了重要作用,而突发事件检测作为网络舆情分析的重要分支,也逐步受到了国内外学者的关注。因此,研 究在微博空间中的海量信息中挖掘突发事件,具有重要意义。

按照突发特征识别的顺序,突发事件识别可以分为以文本为中心的方法和以突发特征为中心的方法。前者是先进行文本聚类,再在类中抽取出突发特征,进行突发事件的识别; 后者是先抽取出突发特征,再对突发特征进行分组,使用突发特征组进行突发事件的识别。

相关论文研究

Bursty and Hierarchical Structure in Streams, Jon Kleinburg

This work appears in the Proceedings of the 8th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, 2002.

虽然作者提出的方法当时主要是为了发现长文本中的突发事件，但是对我们处理微博这样的短文本有重要的借鉴意义。

或许生成一个消息到达时间序列的最简单随机模型是基于指数分布的：消息以一定的概率模式发出，所以两个相邻消息i和i+1的间隔x有指数分布密度函数

$ f(x) = \alpha e^{-\alpha x}, \alpha > 0 $

所以，

$ F(x) = \int_0^x f(x) = 1 - e^{-\alpha x} $
$ E(x) = \alpha^{-1} $

A two-state model

假设对于事件模型有正常状态q0和突发状态q1，可以看做“低态”和“高态”。当A处在状态q0时，消息以低频率发出，连续消息的
间隔x有密度函数$ f_0(x) = \alpha_0 e^{-\alpha_0 x} $; 当A处在状态q1时，消息以高频发出，连续消息间隔x符合密度函数
$ f_1(x) = \alpha_1 e^{-\alpha_1 x}$, 显然 $\alpha_1 > \alpha_0$。并且，A从q0转换到q1（从q1转换到q0）的概率$p\in(0,1)$.

当给出一个消息集合的时候，我们可以用这个生成模型去寻找一个可能的状态序列。
假设我们知道了n+1个消息，那么可以计算出n个时间间隔$ X = (x_1, x_2, \cdots, x_n)$.
令X对应的状态序列$ q = (q_{i_1}, q_{i_2}, \cdots, q_{i_n}) $, 所以,

$ f_q(x_1, x_2, \cdots, x_n) = \prod_{t=1}^n f_{i_t}(x_t) $

令b为状态序列q中状态的转换次数，所以q的先验概率为，

$ P(q) = (\prod_{i_t \neq i_{t+1}})(\prod_{i_t = i_{t+1}}) = p^b(1-p)^{n-b} $

所以,

$ P(q \mid x) = \frac{P(q)f_q(x)}{\sum_{q^{\prime}}P(q^{\prime})f_{q^{\prime}}(x)} = \frac{1}{Z}(\frac{p}{1-p})^b(1-p)^n \prod_{t=1}^nf_{i_t}(x_t) $

其中，$ Z = \sum_{q^{\prime}}P(q^{\prime})f_{q^{\prime}}(x) $. 对上式中两边同时取对数后取反，得

$ -\ln P(q \mid x) = b\ln \left(\frac{1-p}{p}\right) + \left(\sum_{t=1}^n-\ln f_{i_t}(x_t)\right) - n\ln (1-p) + lnZ $

为求得最可能的状态序列，只需使等式右边最小，而等式右边后两项与状态序列q无关，所以

$ c(q\ \mid\ x) = b\ln \left(\frac{1-p}{p}\right) + \left(\sum_{t=1}^n-\ln f_{i_t}(x_t)\right) $

现在只需使$ c(q \mid x) $最小。从直觉上，使状态变化次数尽量少，同时使状态序列适应时间间隔x的值。具体的解决方法在下面详细介绍。

An infinite-state model

已知n+1个消息的时间序列，总时长为T，令时间间隔的估计量$ \hat{g} = \frac{T}{n} $。相比于上面的两状态，我们在这里把状态扩展到多个，
即状态序列$ q = (q_{i_1}, q_{i_2}, \cdots, q_{i_n})$中的每一个状态都有可能属于多个状态中的一个。并令$ \alpha_0 = \hat{g}^{-1} = \frac{n}{T} $,
对于i > 0, $ \alpha_i = \alpha_0 s^i $.
令代价公式为$ \tau(\cdot, \cdot)$, 表示从低密度状态转换到高密度状态的代价。当q从i转换到j时$(f_i(x) < f_j(x))$，

$\tau(i, j) = (j - i)\gamma \ln n, \gamma > 0 $. 特别地，如果$f_j(x) < f_i(x)$, 代价为0。

$ c(q\mid x) = \left(\sum_{t=0}^{n-1}\tau(i_t, i_{t+1}\right) + \left(\sum_{t=1}^n-\ln f_{i_t}(x_t)\right) $

令$i_0 = 0$, 所以无限状态机$ A_{s,\gamma}^*$ 从$q_0$开始.

Computing a minimum-cost state sequence

已知间隔序列$x = (x_1, x_2, \cdots, x_n)$, 求状态机$ A_{s,\gamma}^* $的状态序列$ q = (q_{i_1}, q_{i_2}, \cdots, q_{i_n}) $，使代价
$c(q\ \mid\ x)$最小。并且$A_{s,\gamma}^k$表示状态机只有k个状态选择。

Theorem

$\delta(x) = min_{i=1}^n x_i$ and $k = \lceil 1 + \log_s T + \log_s \delta(x)^{-1} \rceil $.

如果$q^{*}$是状态机$A_{s,\gamma}^k$的最优状态序列，那么它也是状态机$A_{s,\gamma}^*$的最优状态序列。

证明：

令$ q^{*} = (q_{l_1}, q_{l_2}, \cdots, q_{l_n}) $为状态机$ A_{s,\gamma}^k $的最优状态序列，$ q = (q_{i_1}, q_{i_2}, \cdots, q_{i_n}) $
为状态机$ A_{s,\gamma}^{*} $的最优状态序列，$l_0 = i_0 = 0, l_{n+1} = i_{n+1} = 0 $.所以，我们只需证明$ c(q^{*}\ \mid\ x) \leq c(q\ \mid\ x) $.
令$q^{\prime} = (q_{i^{\prime}_1}, \cdots, q_{i^{\prime}_n}), i^{\prime}_t = min(i_t, k-1)$, 显然

$ \sum_{t=0}^{n-1} \tau(i^{\prime}_t, i^{\prime}_{t+1}) \leq \sum_{t=0}^{n-1} \tau(i_t, i_{t+1}) $

对于代价函数的后一项，$-\ln f_j(x_t) = \alpha_j x_t - \ln \alpha_j, 1 \leq t \leq n$ , 考虑当j为何值时其最小。

而$h(\alpha) = \alpha x_t - \ln\alpha$ 在区间$(0, \infty)$上是凸函数，并且当$\alpha = x_t^{-1}$时取得全局最小值。
我们知道$ k = \lceil 1 + \log_s T + \log_s \delta(x)^{-1} \rceil $ , 所以，

$
\alpha_{k-1} = \hat{g}^{-1}s^{k-1} = \frac{n}{T} \cdot s^{k-1} \geq \frac{1}{T} \cdot s^{\log_s^T + \log_s^{\delta(x)^{-1}}} = \frac{1}{T} \frac{T}{\delta(x)} = \frac{1}{\delta(x)}.
$

因为对于所有的t, 有$\delta(x)^{-1} \geq x_t^{-1}$, 所以当$i_t > k-1$时, 有$\alpha_{i_t} > \alpha_{k-1}$, 所以,

$
-\ln f_{i_t^{\prime}}(x_t) \leq -\ln f_{i_t}(x_t), i_t > i_t^{\prime} = k-1
$


由此，可得，

$
c(q^{\prime}\ \mid\ x) = \left(\sum_{t=0}^{n-1}\tau(i_t^{\prime}, i_{t+1}^{\prime}) \right) + \left(\sum_{t=1}^n-\ln f_{i_t^{\prime}}(x_t)\right)\
\leq \left(\sum_{t=0}^{n-1}\tau(i_t,i_{t+1})\right) + \left(\sum_{t=1}^n-\ln f_{i+t}(x_t)\right) = c(q\ \mid\ x).
$

又因为$q^{\prime}$是自动机$A_{s,\gamma}^k$中的一个状态序列，$q^*$是自动机$A_{s,\gamma}^k$中的最优状态序列,所以有

$c(q^*\ \mid\ x) \leq c(q^{\prime}\ \mid\ x) \leq c(q\ \mid\ x).$

得证。

如何计算得到最优状态序列？

令$C_j(t)$是对于输入$x_1, x_2, \cdots, x_n$以状态$q_j$结尾的最小代价值。
$C_j(t) = -\ln f_j(x_t) + min_l(C_l(t-1) + \tau(l,j))$, 并且$C_0(0)=0, C_j(0)=\infty\ j>0$.

源代码：

#include <iostream>
#include <algorithm>
#include <fstream>
#include <vector>
#include <fstream>
#include <string>
#include <boost/lexical_cast.hpp>
#include <boost/algorithm/string/trim.hpp>

using namespace std;

struct Node {
    int const state;
    double const w;
    Node *parent;
    Node(int s, Node *p, double w): state(s), parent(p), w(w) {}
};

double logf(vector<double> const& alpha, vector<double> const& ln_alpha, int i, double x){
    return ln_alpha[i] - alpha[i] * x;
}

double tou(double gammalogn, int i, int j){
    return i >= j ? 0.0 :(j - i) * gammalogn;
}

vector<int> kleinberg_algorithm(vector<double> timeseries, double alpha0, double const s=2, double const gamma=1.0){
    // The number of events.
    size_t N = timeseries.size();

    // Calculate time intervals between successive events.
    vector<double> intervals(N-1);
    for(size_t i=0; i<intervals.size(); i++){
        intervals[i] = timeseries[i+1] - timeseries[i];
    }

    // The minimum interval.
    double const delta = *min_element(intervals.begin(), intervals.end());

    // The time length of the whole timeseries.
    double const T = timeseries.back() - timeseries.front();

    // The upper limit of burst levels.
    int const K = int(ceil(1 + log(T/delta) / log(s)));
    // int const K = 2;

    // Set alpha and ln_alpha
    vector<double> alpha(K);
    vector<double> ln_alpha(K);
    alpha[0] = N / T;
    // alpha[0] = alpha0;
    ln_alpha[0] = log(alpha[0]);
    for(int i=1; i<K; i++){
        alpha[i] = s * alpha[i-1];
        ln_alpha[i] = log(alpha[i]);
    }

    double const gammalogn = gamma * log(double(N));

    vector<Node*> q(K);
    for(auto &it : q){
        it = NULL;
    }
    vector<double> C(K, numeric_limits<double>::infinity());
    C[0] = 0;

    // Start optimization.
    for(auto it : intervals){
        double interval = it;
        vector<Node*> q_new(K);
        vector<double> C_new(K);
        for(int i=0; i<K; i++){
            vector<double> c(K);
            for(int j=0; j<K; j++){
                c[j] = C[j] + tou(gammalogn, j, i);
            }
            size_t const j_min = min_element(c.begin(), c.end()) - c.begin();
            C_new[i] = -logf(alpha, ln_alpha, i, interval) + c[j_min];
            q_new[i] = new Node(i, q[j_min], C_new[i]);
        }

        q_new.swap(q);
        C_new.swap(C);
    }

    size_t const seq_min = min_element(C.begin(), C.end()) - C.begin();
    vector<int> bursts(N);
    size_t count = 0;
    for(Node* p=q[seq_min]; p!=NULL; p=p->parent){
        bursts[N - ++count] = p->state;
        // cout << p->state << " : " << p->w << endl;
    }

    return bursts;
}
int main(int argc, char** argv){
    if(argc < 4){
        cerr << "Usage: a.out s gamma time-stamp-file" << endl;
        exit(1);
    }

    double const s = boost::lexical_cast<double>(argv[1]);
    double const gamma = boost::lexical_cast<double>(argv[2]);
    double alpha0 = boost::lexical_cast<double>(argv[3]);
    ifstream ifs(argv[4]);

    vector<double> timeseries;

    for(string line; getline(ifs, line);){
        boost::trim(line);
        timeseries.push_back(boost::lexical_cast<double>(line));
    }

    vector<int> bursts = kleinberg_algorithm(timeseries, alpha0, s, gamma);

    for(auto it : bursts){
        cout << it << endl;
    }

}

参考资料

1.维基百科
2.Kleinberg’s Burst Detection Algorithm
3.http://iv.slis.indiana.edu/sw/burst.html
4.https://www.cs.cornell.edu/home/kleinber/bhs.pdf