Temporal Point Processes

目录

• TPP
• Evolutionary point processes
  • Conditional intensity function [$t$]
  • Conditional intensity function [$t, \kappa$]
• Inference
• Simulation
  • Inverse Method
  • Ogata’s modified thinning algorithm
• 例子

Rasmussen J. G. Lecture notes: Temporal point processes and the conditional intensity function. 2018.

TPP

• 点过程, 里面的内容感觉和 生存模型 很像.

• 在日常生活中, 我们常常会遇到和时间相关的事件序列:

  $$h_0=(t_0, \kappa_0), h_1=(t_1, \kappa_1), \ldots, h_n=(t_n, \kappa_n).$$

  其中 $t$ 表示事件 $\kappa$ 发生的时间.

• 一个惯常的任务是, 基于上述的历史信息, 预测下一可能发生的事件 $\kappa_{n+1}$, 以及它所发生的具体时间 $t_{n+1}$.

• 假设没有两个事件在同一时刻发生, 则必有

  $$t_0 < t_1 < \cdots < t_{n} < t_{n+1}.$$

• 倘若我们用 $f(\cdot)$ 来表示概率密度, 则有

  $$\tag{1} f(h_0, h_1, \ldots, h_n, h_{n+1}) = \prod_{i=0}^{n+1} f(h_i| \mathcal{H}_{i-1}),$$

  其中 $\mathcal{H}_{i-1} := [h_0, h_1, \ldots, h_{i-1}]$.

• 倘若我们能够建模每一个 $f(h_i |h_{<i})$, 则我们就能够很自然地建模联合概率. TPP 就是这样的一个工具.

Evolutionary point processes

• 这部分, 我们仅考虑时间的预测, 即建模:

  $$f(t_0, t_1, \ldots, t_{n+1}) = \prod_{i=1}^{n+1} f(t=t_{i} | \mathcal{T}_{i-1}).$$

  其中 $\mathcal{T}_{i-1} := [t_0, t_1, \ldots, t_{i-1}] = \{t_j \in \mathcal{T}: t_j \le t_{i-1}\}$.

• 例子 [Renewal process]: 我们以 Renewal Process 这种最简单的随机过程为例. Renewal Process 需要满足: $\{t_{i} - t_{i-1}\}_{i=1}^{n+1}$ 为独立同分布序列, 即两个事件间的间隔是互相独立且服从统一分布的. 此时我们有:

  $$f(t=t_i|\mathcal{T}_{i-1}) = g(t_i - t_{i-1}),$$

  这里 $g$ 是某个满足分布条件的密度函数. 比如, 我们可以取 Gamma 分布:

  $$g(\Delta t; \beta, \alpha) = \frac{\beta^{\alpha}}{\Gamma (\alpha)} \Delta t^{\alpha - 1} e^{-\beta \Delta t}.$$

• 下图是设定不同的 $\beta, \alpha$ 后的结果 (注意 $\beta=0.1, \alpha=1$ 实际上是泊松过程).

# %%

import numpy as np
from scipy.stats import gamma
from freeplot import FreePlot

# %%

def sample_from_gamma(beta: float, alpha: float, nums: int = 10) -> float:
    return gamma.rvs(
            a = alpha, scale=1 / beta, size=(nums,)
        )


# %%


nums = 20

fp = FreePlot(
    figsize=(2, 5)
)

# Gamma(0.02, 0.2)
times = np.cumsum(sample_from_gamma(0.02, 0.2, nums))
fp.scatterplot(
    x=times, y=np.ones_like(times) * 1, label='Gamma(0.02, 0.2)', style=['scatter', {'markers.fillstyle': 'none'}]
)

# Gamma(0.1, 1)
times = np.cumsum(sample_from_gamma(0.1, 1, nums))
fp.scatterplot(
    x=times, y=np.ones_like(times) * 0, label='Gamma(0.1, 1)'
)

# Gamma(2, 20)
times = np.cumsum(sample_from_gamma(2, 20, nums))
fp.scatterplot(
    x=times, y=np.ones_like(times) * -1, label='Gamma(2, 20)'
)

fp[0, 0].legend()
fp.show()

Conditional intensity function [$t$]

• 实际中, 往往不似 Renewal process 这般简单, TPP 求助于如下的 conditional intensity function 来确定最终的密度函数 $f$ (下面的 * 用来强调它是基于 $\mathcal{T}_{t^-}$ 的):

  $$\lambda^*(t) := \lim_{\Delta t \rightarrow 0} \frac{\mathbb{P}(t_{n+1} \in [t, t + \Delta t]| \mathcal{T}_{t^-})}{\Delta t},$$

  其中 $\mathcal{T}_{t^-} := \{t_i \in \mathcal{T}: t_i < t\}$. 后面我们令 $\mathcal{T}_{t_n} = \mathcal{T}_n$.

• 令,

  $$F(t|\mathcal{T}_{t_n}) = \int_{t_n}^t f(t|\mathcal{T}_{t_n}) dt = \mathbb{P}(t_{n+1} \in [t_n, t]| \mathcal{T}_{t_n})$$

  为分布函数.

• 则, 我们有如下关系成立:

  $$\begin{array}{ll} \lambda^*(t) &= \lim_{\Delta t \rightarrow 0} \frac{\mathbb{P}(t_{n+1} \in [t, t + \Delta t]| \mathcal{T}_{t^-})}{\Delta t} \\ &= \lim_{\Delta t \rightarrow 0} \frac{\mathbb{P}(t_{n+1} \in [t, t + \Delta t]|t_{n+1} \not \in (t_{n}, t), \mathcal{T}_{t_n})}{\Delta t} \\ &= \lim_{\Delta t \rightarrow 0} \frac{\mathbb{P}(t_{n+1} \in [t, t + \Delta t]|\mathcal{T}_{t_n})}{\Delta t \cdot \mathbb{P}(t_{n+1} \not \in (t_{n}, t)| \mathcal{T}_{t_n})} \\ &= \lim_{\Delta t \rightarrow 0} \frac{\mathbb{P}(t_{n+1} \in [t, t + \Delta t]|\mathcal{T}_{t_n})}{\Delta t \cdot \mathbb{P}(t_{n+1} \not \in (t_{n}, t)| \mathcal{T}_{t_n})} \\ &= \frac{1}{1 - F(t|\mathcal{T}_{t_n})} \cdot \lim_{\Delta t \rightarrow 0} \frac{\mathbb{P}(t_{n+1} \in [t, t + \Delta t]|\mathcal{T}_{t_n})}{\Delta t} \\ &= \frac{1}{1 - F(t|\mathcal{T}_{t_n})} \cdot f(t|\mathcal{T}_{t_n}). \\ \end{array}$$

• 一旦我们确定(或者给定) $\lambda^*(t)$ (且满足一定的条件), 我们就可以唯一确定 $f, F$, 注意到

  $$\lambda^*(t) =\frac{f(t|\mathcal{H}_{t_n})}{1 - F(t| \mathcal{T}_{t_n})} =\frac{\frac{\mathrm{d}F(t|\mathcal{H}_{t_n})}{dt}}{1 - F(t| \mathcal{T}_{t_n})} =-\frac{\mathrm{d}}{\mathrm{d}t} \log (1 - F(t|\mathcal{T}_{t_n})).$$

  两边对 $[t_n, t]$ 积分, 得到:

  $$\begin{array}{ll} &-\log (1 - F(t|\mathcal{T}_{t_n})) = \int_{t_n}^t \lambda^*(s) \mathrm{d}s \\ \Rightarrow & F(t|\mathcal{T}_{t_n}) = 1 - \exp(-\int_{t_n}^t \lambda^*(s) \mathrm{d}s) \\ \Rightarrow & f(t|\mathcal{T}_{t_n}) = \lambda^*(t) \exp(-\int_{t_n}^t \lambda^*(s) \mathrm{d}s). \end{array}$$

  注: $\mathcal{T}_{t_n}$ 建模在 $*$ 中.

• $\lambda^*(t)$ 所满足的条件, 即令所推到得到的 $f$ 满足密度函数应有的性质:

  1. $\lambda^*(t)$ 在 $t > t_n$ 上非负可积;
  2. $\int_{t_n}^t \lambda^* (s) \mathrm{d}s \mathop{\longrightarrow} \limits^{t \rightarrow \infty} \infty$.

Conditional intensity function [$t, \kappa$]

• 这里, 我们同时考虑发生时间和事件, 此时我们类似地定义:

  $$\lambda^*(t, \kappa) := \lim_{\Delta t \rightarrow 0} \frac{\mathbb{P}(t_{n+1} \in [t, t + \Delta t], \kappa_{n+1} = \kappa| \mathcal{H}_{t^-})}{\Delta t},$$

  这里 $\mathcal{H}_{t^-} := \{h_i \in \mathcal{H}: i < t\}$. 注意, 这里我们考虑事件 $\kappa$ 所属的空间是离散的, 实际上, 连续的也可以 (改成 $\kappa + \Delta \kappa$) (但是这里再用区间符号就怪怪的, 但是以离散的眼光审视吧, 结果是一样的).

• 则, 我们有如下关系成立:

  $$\begin{array}{ll} \lambda^*(t, \kappa) &= \lim_{\Delta t \rightarrow 0} \frac{\mathbb{P}(t_{n+1} \in [t, t + \Delta t], \kappa| \mathcal{H}_{t^-})}{\Delta t} \\ &= \lim_{\Delta t \rightarrow 0} \frac{\mathbb{P}(t_{n+1} \in [t, t + \Delta t], \kappa|t_{n+1} \not \in (t_{n}, t), \mathcal{H}_{t_n})}{\Delta t} \\ &= \lim_{\Delta t \rightarrow 0} \frac{\mathbb{P}(t_{n+1} \in [t, t + \Delta t], \kappa|\mathcal{H}_{t_n})}{\Delta t \cdot \mathbb{P}(t_{n+1} \not \in (t_{n}, t)| \mathcal{H}_{t_n})} \\ &= \lim_{\Delta t \rightarrow 0} \frac{\mathbb{P}(t_{n+1} \in [t, t + \Delta t], \kappa|\mathcal{T}_{t_n})}{\Delta t \cdot \mathbb{P}(t_{n+1} \not \in (t_{n}, t)| \mathcal{H}_{t_n})} \\ &= \frac{1}{1 - F(t|\mathcal{H}_{t_n})} \cdot \lim_{\Delta t \rightarrow 0} \frac{\mathbb{P}(t_{n+1} \in [t, t + \Delta t], \kappa|\mathcal{H}_{t_n})}{\Delta t} \\ &= \frac{1}{1 - F(t|\mathcal{H}_{t_n})} \cdot f(t, \kappa|\mathcal{H}_{t_n}). \\ &= \frac{1}{1 - F(t|\mathcal{H}_{t_n})} \cdot f(t|\mathcal{H}_{t_n}) \cdot f(\kappa|t, \mathcal{H}_{t_n}) \\ &= \lambda^*(t) \cdot f^*(\kappa|t). \end{array}$$

  这里我们定义 $f^*(\kappa|t) := f(\kappa|t, \mathcal{H}_{t_n})$.

• 换言之, 我们只要确定 $\lambda^*, f^*$ 即可确定发生时间 + 发生时间的随机过程.

Inference

• 我不清楚为什么不直接确定 $f$ 而需要一个跳板 $\lambda^*$, 想必是它在推理和抽样上具有特别的性质.

• 很多时候, 我们需要通过极大似然来拟合模型, 所以似然函数的计算是一个比较重要的问题. 假设我们在 $[0, T]$ 时间段内观测到 $[(t_0, \kappa_0), \ldots, (t_n, \kappa_n)]$, 下面介绍其上的似然函数的计算.

• 定义:

  $$\Lambda^*(t) = \int_0^t \lambda^*(s) \mathrm{d}s.$$

• 根据 (1) 可得:

  $$L = (1 - F(T|\mathcal{H}_{t_n})) \cdot \prod_{i=0}^{n+1} f(h_i| \mathcal{H}_{i-1})$$

  比较特别的是 $(1 - F(T|\mathcal{H}_{t_n}))$, 即我们可以确定, $t_{n+1}$ 一定发生在 $>T$ 的时间段内.

• 于是, 我们有:

  $$L = \bigg( \prod_{i=0}^n \lambda^*(t_i) \bigg) \bigg( \prod_{i=0}^n f^*(\kappa_i| t_i) \bigg) \exp \bigg(-\int_0^T \lambda^*(s) \mathrm{d}s \bigg).$$

Simulation

Inverse Method

• 假设序列 $[s_0, s_1, \ldots, t_n]$, 我们有序列 (假设 $\Lambda^*$ 可逆)

  $$[s_0= \Lambda^*(t_0), \ldots, s_n = \Lambda^*(t_n)]$$

  服从 unit rate Poisson process. 即

  $$f(s_i|s_{i-1}, \ldots, s_{0}) = f(s_i|s_{i-1}) = e^{-(s_i - s_{i-1})}, \quad s_i > s_{i-1}.$$

  即间隔时间是互相独立的 (且服从 $\lambda=1$ 的泊松分布).

• 我们有

  $$\begin{array}{ll} F(s|\mathcal{H}_{s_n}) &= \mathbb{P}(s_{n+1} \le s| s_{n}, \ldots, s_0) \\ &= \mathbb{P}(\Lambda^*(t_{n+1}) \le s| t_{n}, \ldots, t_0) \\ &= \mathbb{P}(\Lambda^*(t_{n+1}) \le s| \mathcal{H}_{t_n}) \\ &= \mathbb{P}(t_{n+1} \le \Lambda^{*-1}(s)| \mathcal{H}_{t_n}) \\ &= F(\Lambda^{*-1}(s)| \mathcal{H}_{t_n}) \\ &= 1 - \exp(- \int_{t_n}^{\Lambda^{*-1}(s)} \Lambda^*(t) \mathrm{d}t) \\ &= 1 - \exp\bigg(- \big(\Lambda^*(\Lambda^{*-1}(s)) - \Lambda^*(t_n) \big)\bigg) \\ &= 1 - \exp\bigg(- \big(s - s_n \big)\bigg) \\ &= \int_{s_n}^s \exp(-t) \mathrm{d}t. \\ \end{array}$$

  故成立.

• 同理, 如果 $[s_0, s_1, \ldots, s_n]$ 服从 unit rate Poisson process, 则

  $$[t_0 = \Lambda^{*-1}(s_0), \ldots, t_n=\Lambda^{*-1}(s_n)]$$

  服从密度函数为 $\lambda^*$ 的点过程.

• 如此一来, 我们想要从该过程中抽样, 只需要先从 unit rate Poisson process 中抽样然后再做变换即可.

Ogata’s modified thinning algorithm

• 这部分没看明白, 有兴趣的可以回看原文.

例子

• Hawkes process: 它的密度函数定义为:

  $$\lambda^*(t) = \mu + \alpha \sum_{t_i < t} \exp(-(t - t_i)),$$

  我们可以得到:

  $$\begin{array}{ll} \Lambda^*(t) &=\int_0^t \lambda^*(s) \mathrm{d}s \\ &=\int_0^t \mu + \alpha \sum_{t_i < s} \exp(-(s - t_i)) \mathrm{d}s \\ &=\mu t + \alpha \int_0^t \sum_{t_i < s} \exp(-(s - t_i)) \mathrm{d}s \\ &=\mu t + \alpha \sum_{t_i < t} \int_{t_i}^t \exp(-(s - t_i)) \mathrm{d}s \\ &=\mu t + \alpha \sum_{t_i < t} \bigg(1 - \exp(-(t - t_i)) \bigg). \\ \end{array}$$

• 可惜的是, $\Lambda^{*-1}$ 没有显式解, 一个近似策略是:

  1. 采样 $[s_0, s_1, s_2, \ldots]$;
  2. $t$ 从 $0$ 开始逐步增加, 一旦
     $$\Lambda^*(t) \approx s_0 \Rightarrow t_0 \approx t,$$
     然后继续进行逐步得到 $t_1, t_2, \ldots$.

• 上图是在 $\mu=0.5, \alpha=0.9$ 的情况下采样的 20 个点, 以及所对应的 $\lambda^*(t)$ 的可视化. 可以发现, Hawkes process 有这样的性质: 在 $t_i$ 处倘若发生了事件, 则 $t_i + \Delta t$ 一个时间内发生的概率往往会变大, 从而更容易出现聚合的效应.

• 下面是代码:

# %%

import numpy as np
from scipy.stats import gamma
from freeplot import FreePlot

# %%

def sample_from_gamma(beta: float, alpha: float, nums: int = 10) -> float:
    return gamma.rvs(
            a = alpha, scale=1 / beta, size=(nums,)
        )


def hawkes_conditional_intensity_function(t, ts: list, mu: float, alpha: float):
    # \lambda^*(t)
    ts = np.array(ts)
    ts = ts[ts < t]
    return mu + alpha * np.sum(np.exp(ts - t))

def hawkes_integrated_conditional_intensity_function(t, ts: list, mu: float, alpha: float):
    # \Lambda^*(t)
    ts = np.array(ts)
    return mu * t + alpha * np.sum(1 - np.exp(ts - t))

# %%



nums = 20
mu = 0.5
alpha = 0.9
step_size = 0.01

ss = np.cumsum(sample_from_gamma(1, 1, nums))
ts = []
t = 0

for s in ss:
    s_hat = hawkes_integrated_conditional_intensity_function(
        t, ts, mu, alpha
    )
    while s_hat < s:
        t += step_size
        s_hat = hawkes_integrated_conditional_intensity_function(
            t, ts, mu, alpha
        )
    ts.append(t)

points = np.linspace(0, ts[-1] + 1, 100)
densities = [
    hawkes_conditional_intensity_function(
        point, ts, mu, alpha
    ) for point in points
]

fp = FreePlot(
    figsize=(2, 5)
)

fp.scatterplot(
    x=ts, y=np.zeros_like(ts), style=['scatter', {'markers.fillstyle': 'none'}]
)

fp.lineplot(
    points, densities, marker=''
)

fp.show()