1

我想使用 Hawkes 过程对一些数据进行建模。我找不到 PyMC 是否支持霍克斯进程。更具体地说，我想要一个带有霍克斯过程的观察变量，并学习其参数的后验。

如果它不存在，那么我可以在 PyMC 中以某种方式定义它，例如 @deterministic 等吗？

4

1 回答 1

3

自从您提出问题以来已经有很长时间了，但是我今天已经在 PyMC 上解决了它，所以我想我会为可能遇到相同问题的其他人分享我的实施要点。我们将推断霍克斯过程的参数 λ 和 α。我不打算介绍时间尺度参数 β，我将把它作为练习留给读者。

首先让我们生成一些数据：

def hawkes_intensity(mu, alpha, points, t):
    p = np.array(points)
    p = p[p <= t]
    p = np.exp(p - t)
    return mu + alpha * np.sum(p)


def simulate_hawkes(mu, alpha, window):
    t = 0
    points = []
    lambdas = []
    while t < window:
        m = hawkes_intensity(mu, alpha, points, t)
        s = np.random.exponential(scale=1/m)
        ratio = hawkes_intensity(mu, alpha, points, t + s)

        t = t + s
        if t < window:
            points.append(t)
            lambdas.append(ratio)
        else:
            break
    points = np.sort(np.array(points, dtype=np.float32))
    lambdas = np.array(lambdas, dtype=np.float32)
    return points, lambdas


# parameters
window = 1000
mu = 8
alpha = 0.25
points, lambdas = simulate_hawkes(mu, alpha, window)
num_points = len(points)

我们刚刚使用我从那里改编的一些函数生成了一些时间点：https ://nbviewer.jupyter.org/github/MatthewDaws/PointProcesses/blob/master/Temporal%20points%20processes.ipynb

现在，诀窍是创建一个大小为 (num_points, num_points) 的矩阵，其中包含第 i 个点与所有其他点的时间距离。所以矩阵的 (i, j) 点是第 i 个点到第 j 个点的时间间隔。该矩阵将用于计算霍克斯过程的指数之和，即。自励的部分。创建这个矩阵以及指数总和的方法有点棘手。我建议你自己检查每一行，这样你就可以看到他们做了什么。

tile = np.tile(points, num_points).reshape(num_points, num_points)
tile = np.clip(points[:, None] - tile, 0, np.inf)
tile = np.tril(np.exp(-tile), k=-1)
Σ = np.sum(tile, axis=1)[:-1]  # this is our self-exciting sum term

我们有点，我们有一个包含激励项总和的矩阵。霍克斯过程的两个连续事件之间的持续时间遵循参数 λ = λ0 + ∑ 激发的指数分布。这就是我们要建模的内容，但首先我们必须计算生成数据的两个连续点之间的持续时间。

interval = points[1:] - points[:-1]

我们现在可以进行推理了：

with pm.Model() as model:
    λ = pm.Exponential("λ", 1)
    α = pm.Uniform("α", 0, 1)

    lam = pm.Deterministic("lam", λ + α * Σ)
    interarrival = pm.Exponential(
        "interarrival", lam, observed=interval)


    trace = pm.sample(2000, tune=4000)
    pm.plot_posterior(trace, var_names=["λ", "α"])
    plt.show()

    print(np.mean(trace["λ"]))
    print(np.mean(trace["α"]))

7.829 0.284

注意：tile如果您有很多数据点，矩阵可能会变得非常大。

于 2019-07-29T15:22:31.877 回答