ホークス過程(Hawkes Process)の説明とPythonでのシミュレーション

数学 プログラミング

ホークス過程とは

イベントが起こった直後に、再びそのイベントが起こりやすい現象をモデル化する時に便利なモデル

例えば、大きな地震が起こったときは、直後に余震が起こりやすい

ほかにも、くしゃみをした時は、そのあと連続でくしゃみが続きやすい

ここで、以下の３つの特徴を持つintensity (=強度関数)というものを考える

  1. intensityが大きいほど次のイベントがすぐ起こりやすい
  2. またintensityは時間とともに減っていく
  3. イベントが起こるとintensityが上昇する

このように考えると、イベント発生と、intensityの変化は以下のように図示できる

f:id:babaye:20191017204947j:plain
Rizoiu, et al. (2017)

イベント（地震）が起きるとintensityが上昇し、なにもないときは緩やかに下降しているのが確認できるだろう（一番下のグラフ）

ホークス過程の数式

intensityをλとして表す

指数関数に従ってintensityが減衰すると仮定するのが最もシンプルでよく使われている

\displaystyle λ(t)=ρ +α \sum_{t_{k} \lt t} exp^{- β (t-t_{k} ) }

ただし、ρはbase intensityと呼ばれ、平時にイベントが起こる確率

αはイベントが起こったときにどれくらいintensityが上昇するかを表すパラメータ

βはintensityの減少の早さを示すパラメータ

αとβの設定の仕方次第でintensityが発散してしまうので注意

ここでの \sumでは、過去(つまり t_k \lt t)に起きたすべての各イベントによるintensityの影響を足している

たとえば、計測開始から２秒後と５秒後にイベントが起こったとすると、（６秒の地点のintensity）＝（base intensityのρ） ＋ （１秒前に起こったイベントによる増加分たくさん） + （４秒前に起こったイベントによる増加分すこし）　　　である

コード

import numpy as np
import matplotlib.pyplot as plt

# hyperparameters
mu = 0.1
alpha = 1.5
beta = 30

# 
neuton_iter_num = 200
num_datapoints = 15

# step1
U = np.random.uniform(0, 1)
# step2
t1 = -np.log(U)/mu
t = [t1]


for k in range(num_datapoints):
    # step3
    U = np.random.uniform(0, 1)
    
    # step 4
    if k == 0:
        Sk = 1

    u = t[k] - np.log(U)/mu  # initial value of iteration
    for _ in range(neuton_iter_num):
        fu  = np.log(U) + mu*(u-t[k])  + alpha/beta*Sk*(1-np.exp(-beta*(u-t[k])))
        fu_dash = mu+alpha*Sk*np.exp(-beta*(u-t[k]))
        u -= fu/fu_dash

    # step 5
    t.append(u)
    Sk = np.exp(-beta*(t[k+1]-t[k]))*Sk + 1  # find S(k+1)

t = np.array(t)/max(t)

# plt.figure(figsize=(20, 10), dpi=50)
# plt.scatter(t, np.ones(len(t)), facecolors='none', edgecolors='red')
# plt.show()


# intensity
time = np.linspace(0,1,10000)
increase = list(map(lambda x: alpha*np.exp(-beta * x), time))
intensity = np.zeros(10000*2)
for timepoints in t:
    yuko = np.append(np.append(np.zeros(len(time[time<timepoints])), increase), np.zeros(len(time[time>=timepoints])))    
    intensity += yuko


# plt.figure(figsize=(20, 10), dpi=50)
# plt.plot(time, intensity[:10000])
# plt.show()



fig, ax1 = plt.subplots(figsize=(20, 5))
ax2 = ax1.twinx()
ax1.scatter(t, np.ones(len(t)), facecolors='none', edgecolors='red', marker='*', s=150)
ax2.plot(time, intensity[:10000], c='b', lw=1)
plt.show()

出力結果：

f:id:babaye:20191017214323p:plain
上記コードの出力結果

参考文献