点過程の理論 — ポアソン過程の一般化

コールセンターに電話がかかってくるタイミング、地震が発生する時刻、ニューロンが発火する瞬間 — これらはすべて「時間軸上にランダムに散らばった点(イベント)」として捉えることができます。このような「ランダムな点の集合」を数学的に扱う枠組みが点過程(point process)です。

最も基本的な点過程であるポアソン過程では、イベントは一定の率で互いに独立に発生します。しかし、現実のイベントの多くはこの単純なモデルでは捉えきれません。地震は余震が連鎖的に発生し(時間的クラスタリング)、コールセンターへの電話は時間帯によって頻度が変わります(非定常性)。このような現象を扱うために、ポアソン過程を一般化した点過程の理論が必要になります。

点過程の理論は以下の分野で広く応用されています。

  • 地震学: ETASモデル(Epidemic-Type Aftershock Sequence)は、余震の連鎖を自己励起型点過程でモデル化します
  • 神経科学: ニューロンのスパイク列(発火パターン)を点過程として解析し、外部刺激に対する応答特性や神経回路の結合関係を推定します
  • 待ち行列理論: 顧客やジョブの到着過程を非定常ポアソン過程やリニューアル過程として記述し、待ち時間やキュー長を解析します
  • 金融工学: 高頻度取引データにおける注文の到着をホークス過程でモデル化し、流動性やボラティリティのダイナミクスを分析します
  • 疫学: 感染症の発症イベントを空間的・時間的な点過程として分析し、感染の拡散パターンやクラスターを検出します

本記事では、点過程の基本概念をポアソン過程から出発して体系的に構築し、非定常ポアソン過程、リニューアル過程、自己励起型点過程(ホークス過程)へと段階的に一般化します。条件付き強度関数という統一的な記述法を中心に据え、尤度関数に基づくパラメータ推定やモデル診断の方法も解説します。点過程の理論は確率過程論の重要な一分野であると同時に、データ駆動科学の多くの分野で実用的なモデリングツールとして活用されています。

本記事の内容

  • 点過程の定義と計数過程
  • ポアソン過程の復習と特徴づけ
  • 非定常(非斉次)ポアソン過程と強度関数
  • リニューアル過程と到着間隔
  • 条件付き強度関数と自己励起型点過程
  • 補償子とマルチンゲール表現
  • Pythonによるシミュレーションと統計解析

前提知識

この記事を読む前に、以下の記事を読んでおくと理解が深まります。

点過程とは — ランダムなイベントの数学

計数過程としての点過程

点過程を最も自然に捉える方法は、計数過程(counting process)としての見方です。時刻 $t$ までに発生したイベントの数を $N(t)$ と書くと、$N(t)$ は以下の性質を持つ確率過程です。

  1. $N(0) = 0$(時刻0ではまだイベントは発生していない)
  2. $N(t) \in \{0, 1, 2, \dots\}$(非負整数値)
  3. $s < t$ ならば $N(s) \leq N(t)$(単調非減少)
  4. $N(t) – N(s)$ は区間 $(s, t]$ に発生したイベントの数

$N(t)$ の経路は右連続な階段関数(cadlag関数)であり、各イベント発生時刻 $T_1 < T_2 < T_3 < \dots$ でちょうど1だけジャンプします。同時に2つ以上のイベントが発生しないと仮定する場合、これを単純点過程と呼びます。

等価的に、点過程はイベント発生時刻の列 $\{T_n\}_{n \geq 1}$ として定義することもできます。計数過程と発生時刻は一対一に対応しています。

$$ N(t) = \max\{n : T_n \leq t\}, \quad T_n = \inf\{t : N(t) \geq n\} $$

到着間隔

連続するイベント間の時間を到着間隔(inter-arrival time)と呼び、$\tau_n = T_n – T_{n-1}$($T_0 = 0$)で定義します。到着間隔の分布が点過程の性質を大きく左右します。

これらの基本概念を使って、まずポアソン過程を特徴づけましょう。

ポアソン過程の特徴づけ

定義と基本性質

率 $\lambda > 0$ の(斉次)ポアソン過程 $N(t)$ は、以下の条件で特徴づけられます。

  1. $N(0) = 0$
  2. 独立増分: 重なりのない区間の増分は独立
  3. 定常増分: $N(t+s) – N(s) \sim \text{Poisson}(\lambda t)$

この3条件から、ポアソン過程の多くの性質が導かれます。

到着間隔: $\tau_n \sim \text{Exp}(\lambda)$(指数分布)であり、$\{\tau_n\}$ は互いに独立です。指数分布の無記憶性 $P(\tau > t + s | \tau > t) = P(\tau > s)$ は、「過去にどれだけ待っても、次のイベントまでの残り時間の分布は変わらない」ことを意味します。

超位置原理: 率 $\lambda_1$ と $\lambda_2$ の独立なポアソン過程の和は、率 $\lambda_1 + \lambda_2$ のポアソン過程になります。

間引き定理: 率 $\lambda$ のポアソン過程の各イベントを確率 $p$ で独立に保持すると、保持されたイベントは率 $\lambda p$ のポアソン過程になります。破棄されたイベントも率 $\lambda(1-p)$ のポアソン過程であり、保持イベントと破棄イベントの過程は互いに独立です。

順序統計量性: 区間 $[0, T]$ で $N(T) = n$ 個のイベントが発生したという条件のもとで、発生時刻 $T_1 < T_2 < \dots < T_n$ は $[0, T]$ 上の一様分布からの順序統計量と同じ分布を持ちます。この性質はポアソン過程の「完全なランダムさ」を最も端的に表しています。

ポアソン過程の限界

ポアソン過程は定常性と独立増分性を仮定しているため、以下のような現象は捉えられません。

  • 非定常性: 朝の通勤ラッシュ時の電話の方が多い
  • クラスタリング: 地震の後に余震が集中する
  • 到着間隔の依存性: 顧客の到着間隔が前回の間隔に依存する

これらを扱うために、ポアソン過程の一般化が必要です。まず非定常ポアソン過程から見ていきましょう。

非定常ポアソン過程と強度関数

定義

非定常(非斉次)ポアソン過程(non-homogeneous Poisson process)は、イベントの発生率が時間とともに変化するポアソン過程です。強度関数(intensity function)$\lambda(t) \geq 0$ を使って、次のように定義されます。

  1. $N(0) = 0$
  2. 独立増分: 重なりのない区間の増分は独立
  3. $N(t) – N(s) \sim \text{Poisson}\left(\int_s^t \lambda(u)\,du\right)$

累積強度関数(integrated intensity)を

$$ \begin{equation} \Lambda(t) = \int_0^t \lambda(u)\,du \end{equation} $$

と定義すると、$N(t) \sim \text{Poisson}(\Lambda(t))$ です。

直感的な意味

微小時間 $dt$ の間にイベントが発生する確率は $\lambda(t)\,dt + o(dt)$ です。$\lambda(t)$ が大きい時間帯ではイベントが密に発生し、小さい時間帯では疎になります。$\lambda(t) = \lambda$(定数)のとき斉次ポアソン過程に帰着します。

シミュレーション法: 間引き法(thinning algorithm)

非定常ポアソン過程をシミュレーションするための効率的な方法として、間引き法(thinning algorithm, Lewis-Shedler法)があります。

  1. $\lambda_{\max} = \max_{0 \leq t \leq T} \lambda(t)$ を求める
  2. 率 $\lambda_{\max}$ の斉次ポアソン過程を生成する
  3. 各イベント発生時刻 $t_i$ を確率 $\lambda(t_i)/\lambda_{\max}$ で採択する

この方法は、間引き定理の逆操作であり、一様に高い率で生成したイベントを選択的に「間引く」ことで非定常なパターンを作り出します。

非定常ポアソン過程は独立増分性を保持しているため、「過去のイベントが未来に影響する」というクラスタリング現象は扱えません。クラスタリングを扱うには、条件付き強度関数という概念が必要です。

条件付き強度関数と自己励起型点過程

条件付き強度関数

点過程の最も一般的な記述は、条件付き強度関数(conditional intensity function)によるものです。

$$ \begin{equation} \lambda^*(t) = \lim_{\Delta t \to 0} \frac{P(N(t + \Delta t) – N(t) = 1 | \mathcal{H}_t)}{\Delta t} \end{equation} $$

ここで $\mathcal{H}_t = \{T_1, T_2, \dots, T_{N(t)}\}$ は時刻 $t$ までのイベント履歴(内部履歴)です。$\lambda^*$ の上付き asterisk は条件付き強度であることを示す慣例的な記法です。

条件付き強度関数は、「過去のイベント履歴が与えられたとき、次のイベントがいつ発生するか」の確率的な率を記述します。ポアソン過程では $\lambda^*(t) = \lambda(t)$(履歴に依存しない)ですが、一般の点過程では $\lambda^*(t)$ は過去のイベントに依存します。

ホークス過程(自己励起型点過程)

ホークス過程(Hawkes process)は、過去のイベントが将来のイベント発生率を一時的に増加させる自己励起型の点過程です。

$$ \begin{equation} \lambda^*(t) = \mu + \sum_{T_i < t} g(t - T_i) \end{equation} $$

ここで $\mu > 0$ はベースライン強度、$g(s) \geq 0$($s > 0$)は励起関数(excitation function)です。各イベントが「波及効果」$g(t – T_i)$ を生み出し、後続のイベントの発生率を一時的に高めます。

最もよく使われる励起関数は指数型です。

$$ g(s) = \alpha \beta e^{-\beta s}, \quad \alpha > 0, \beta > 0 $$

$\alpha$ は各イベントの励起強度、$\beta$ は減衰率です。$\alpha < 1$(各イベントが平均して1未満の子イベントを生む)のとき過程は定常であり、$\alpha \geq 1$ のとき過程は爆発する可能性があります。この構造は分岐過程と深く関係しています。

ホークス過程は以下のように理解できます。「背景ノイズ」(率 $\mu$ の基本的なイベント)の上に、各イベントが「余震」(励起関数 $g$ による追加イベント)を生み出す。余震がさらに余震を生む連鎖反応が起こり得ますが、$\alpha < 1$ ならば連鎖は有限で止まります。

補償子とマルチンゲール

点過程の理論で重要な概念が補償子(compensator)です。計数過程 $N(t)$ の補償子は、

$$ \begin{equation} \tilde{N}(t) = N(t) – \int_0^t \lambda^*(s)\,ds \end{equation} $$

で定義されます。$\tilde{N}(t)$ はマルチンゲールになります。つまり、

$$ E[\tilde{N}(t) | \mathcal{H}_s] = \tilde{N}(s), \quad s < t $$

が成り立ちます。これは「実際のイベント数から予測される累積発生率を引くと、公平なゲーム(マルチンゲール)になる」ことを意味しています。

補償マルチンゲール $\tilde{N}(t)$ は、点過程の統計的推論やモデル診断において中心的な役割を果たします。

以上の理論をPythonで実装し、シミュレーションで確認しましょう。

Pythonによるシミュレーションと検証

ポアソン過程と非定常ポアソン過程

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(42)

T = 10.0  # 観測時間

# --- 斉次ポアソン過程 ---
lam = 5.0  # 率
n_events_hom = np.random.poisson(lam * T)
event_times_hom = np.sort(np.random.uniform(0, T, n_events_hom))

# --- 非斉次ポアソン過程 ---
# 強度関数: λ(t) = 3 + 4*sin(πt/5)^2
def intensity(t):
    return 3 + 4 * np.sin(np.pi * t / 5)**2

lam_max = 7.0  # max of intensity function

# 間引き法
n_candidates = np.random.poisson(lam_max * T)
candidate_times = np.sort(np.random.uniform(0, T, n_candidates))
accept_probs = intensity(candidate_times) / lam_max
accepted = np.random.random(len(candidate_times)) < accept_probs
event_times_nhom = candidate_times[accepted]

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 左上: 斉次ポアソン過程
ax = axes[0, 0]
N_hom = np.arange(1, len(event_times_hom) + 1)
ax.step(np.concatenate([[0], event_times_hom]), np.concatenate([[0], N_hom]),
        where='post', linewidth=1.5, color='steelblue')
ax.plot(event_times_hom, np.zeros_like(event_times_hom), '|', color='red',
        markersize=15, markeredgewidth=1.5)
ax.plot([0, T], [0, lam * T], 'k--', linewidth=1, alpha=0.5, label=rf'$\lambda t$ ($\lambda={lam}$)')
ax.set_xlabel('Time $t$', fontsize=12)
ax.set_ylabel('$N(t)$', fontsize=12)
ax.set_title('Homogeneous Poisson process', fontsize=13)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)

# 右上: 非斉次ポアソン過程
ax = axes[0, 1]
N_nhom = np.arange(1, len(event_times_nhom) + 1)
ax.step(np.concatenate([[0], event_times_nhom]), np.concatenate([[0], N_nhom]),
        where='post', linewidth=1.5, color='steelblue')
ax.plot(event_times_nhom, np.zeros_like(event_times_nhom), '|', color='red',
        markersize=15, markeredgewidth=1.5)
# 累積強度
t_fine = np.linspace(0, T, 500)
Lambda_t = np.array([np.trapz(intensity(np.linspace(0, t, 200)), np.linspace(0, t, 200))
                      for t in t_fine])
ax.plot(t_fine, Lambda_t, 'k--', linewidth=1, alpha=0.5, label=r'$\Lambda(t)$')
ax.set_xlabel('Time $t$', fontsize=12)
ax.set_ylabel('$N(t)$', fontsize=12)
ax.set_title('Non-homogeneous Poisson process', fontsize=13)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)

# 左下: 強度関数
ax = axes[1, 0]
ax.plot(t_fine, intensity(t_fine), 'r-', linewidth=2)
ax.fill_between(t_fine, 0, intensity(t_fine), alpha=0.2, color='red')
# イベント発生位置
ax.plot(event_times_nhom, intensity(event_times_nhom), 'k|', markersize=10)
ax.set_xlabel('Time $t$', fontsize=12)
ax.set_ylabel(r'$\lambda(t)$', fontsize=12)
ax.set_title('Intensity function and events', fontsize=13)
ax.grid(True, alpha=0.3)

# 右下: 到着間隔のヒストグラム比較
ax = axes[1, 1]
iat_hom = np.diff(np.concatenate([[0], event_times_hom]))
ax.hist(iat_hom, bins=30, density=True, alpha=0.6, color='steelblue',
        edgecolor='white', label='Homogeneous IAT')
t_exp = np.linspace(0, max(iat_hom), 100)
ax.plot(t_exp, lam * np.exp(-lam * t_exp), 'b-', linewidth=2,
        label=rf'Exp($\lambda={lam}$)')
ax.set_xlabel('Inter-arrival time', fontsize=12)
ax.set_ylabel('Density', fontsize=12)
ax.set_title('Inter-arrival time distribution', fontsize=13)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('point_process_basic.png', dpi=150, bbox_inches='tight')
plt.show()

上のグラフから、ポアソン過程の基本的な性質が確認できます。

  1. 左上: 斉次ポアソン過程の計数過程は $\lambda t$ に沿って線形に成長しています。ジャンプの間隔はおおむね等しく、定常性を反映しています

  2. 右上: 非斉次ポアソン過程では、強度の高い時間帯でイベントが密集しています。計数過程は累積強度関数 $\Lambda(t)$ に沿って成長しており、非定常な到着パターンが見えます

  3. 左下: 強度関数の高い部分にイベントが集中しています。間引き法により、$\lambda(t)$ に比例した密度でイベントが配置されていることが分かります

  4. 右下: 斉次ポアソン過程の到着間隔は指数分布に従うことが確認できます。ヒストグラムと理論的な指数分布の密度関数が一致しています

ホークス過程のシミュレーション

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(123)

def simulate_hawkes(mu, alpha, beta, T):
    """ホークス過程のシミュレーション(間引き法)"""
    events = []
    t = 0
    lam_max = mu / (1 - alpha)  # 定常状態での平均強度の上界
    lam_max *= 3  # 安全マージン

    while t < T:
        dt = np.random.exponential(1 / lam_max)
        t += dt
        if t >= T:
            break
        # 条件付き強度の計算
        lam_t = mu + alpha * beta * sum(np.exp(-beta * (t - ti)) for ti in events)
        if np.random.random() < lam_t / lam_max:
            events.append(t)
    return np.array(events)

# --- パラメータ ---
mu = 1.0       # ベースライン強度
alpha = 0.5    # 励起強度 (< 1 for stability)
beta = 2.0     # 減衰率
T = 50.0

events = simulate_hawkes(mu, alpha, beta, T)

# 条件付き強度の計算
t_fine = np.linspace(0, T, 2000)
lam_star = np.zeros_like(t_fine)
for i, t in enumerate(t_fine):
    lam_star[i] = mu + alpha * beta * sum(np.exp(-beta * (t - ti))
                                           for ti in events if ti < t)

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 左上: イベントと条件付き強度
ax = axes[0, 0]
ax.plot(t_fine, lam_star, 'r-', linewidth=1, alpha=0.8, label=r'$\lambda^*(t)$')
ax.axhline(mu, color='gray', linestyle='--', linewidth=1, alpha=0.5,
           label=rf'Baseline $\mu={mu}$')
ax.plot(events, np.zeros_like(events), '|', color='steelblue',
        markersize=12, markeredgewidth=1.5)
ax.set_xlabel('Time $t$', fontsize=12)
ax.set_ylabel(r'$\lambda^*(t)$', fontsize=12)
ax.set_title(rf'Hawkes process ($\mu={mu}$, $\alpha={alpha}$, $\beta={beta}$)', fontsize=13)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
ax.set_xlim(0, 30)

# 右上: 計数過程
ax = axes[0, 1]
N_hawkes = np.arange(1, len(events) + 1)
ax.step(np.concatenate([[0], events]), np.concatenate([[0], N_hawkes]),
        where='post', linewidth=1.5, color='steelblue')
# 定常平均強度
lam_mean = mu / (1 - alpha)
ax.plot([0, T], [0, lam_mean * T], 'k--', linewidth=1, alpha=0.5,
        label=rf'$\lambda_{{mean}} t$ ($\lambda_{{mean}}={lam_mean:.1f}$)')
ax.set_xlabel('Time $t$', fontsize=12)
ax.set_ylabel('$N(t)$', fontsize=12)
ax.set_title('Counting process', fontsize=13)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)

# 左下: 到着間隔の分布
ax = axes[1, 0]
iat = np.diff(np.concatenate([[0], events]))
ax.hist(iat, bins=40, density=True, alpha=0.7, color='steelblue',
        edgecolor='white', linewidth=0.5, label='Hawkes IAT')
# ポアソン過程(同じ平均率)との比較
iat_exp = np.linspace(0, max(iat), 100)
ax.plot(iat_exp, lam_mean * np.exp(-lam_mean * iat_exp), 'r--', linewidth=2,
        label=f'Exp($\\lambda={lam_mean:.1f}$)')
ax.set_xlabel('Inter-arrival time', fontsize=12)
ax.set_ylabel('Density', fontsize=12)
ax.set_title('IAT distribution (Hawkes vs Poisson)', fontsize=13)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)

# 右下: 補償マルチンゲール
ax = axes[1, 1]
# Λ(t) = ∫_0^t λ*(s) ds の数値計算
Lambda = np.cumsum(lam_star) * (t_fine[1] - t_fine[0])
# N(t) - Λ(t)
N_t = np.array([np.sum(events <= t) for t in t_fine])
compensated = N_t - Lambda
ax.plot(t_fine, compensated, 'b-', linewidth=1, alpha=0.8)
ax.axhline(0, color='r', linestyle='--', linewidth=1)
ax.set_xlabel('Time $t$', fontsize=12)
ax.set_ylabel(r'$\tilde{N}(t) = N(t) - \Lambda(t)$', fontsize=12)
ax.set_title('Compensated martingale', fontsize=13)
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('hawkes_process.png', dpi=150, bbox_inches='tight')
plt.show()

print(f"イベント数: {len(events)}")
print(f"経験的平均強度: {len(events)/T:.2f}")
print(f"理論的定常強度: {lam_mean:.2f}")

ホークス過程のシミュレーション結果から、以下の重要な特徴が確認できます。

  1. 左上: 条件付き強度がイベント発生のたびにジャンプし、指数的に減衰している。各イベントが $\alpha\beta$ だけ強度を上げ、率 $\beta$ で減衰する様子が見えます。イベントが連続して発生する「クラスター」では強度が高く保たれています

  2. 右上: 計数過程は定常平均強度 $\mu/(1-\alpha) = 2.0$ に沿って成長していますが、ポアソン過程よりも不規則な振る舞いを示しています。クラスターがある時間帯では急激に増加し、イベントが少ない時間帯では停滞します

  3. 左下: ホークス過程の到着間隔分布は指数分布からずれている。短い到着間隔が過剰に多い(クラスタリング効果)一方、長い間隔も多くなっています。指数分布のフィット(赤破線)とのずれはクラスタリングの証拠です

  4. 右下: 補償マルチンゲール $\tilde{N}(t) = N(t) – \Lambda(t)$ がゼロの周りを揺らいでいる。これは条件付き強度の推定が正しいことを示しています。マルチンゲール性の数値的な確認であり、モデル診断の一手法として使われます

点過程の尤度関数

条件付き強度による尤度の表現

観測区間 $[0, T]$ でイベント発生時刻 $\{t_1, t_2, \dots, t_n\}$ が観測されたとき、点過程の対数尤度は条件付き強度関数を用いて次のように表されます。

$$ \begin{equation} \log L = \sum_{i=1}^n \log \lambda^*(t_i) – \int_0^T \lambda^*(t)\,dt \end{equation} $$

この式は直感的に理解できます。第一項 $\sum_i \log \lambda^*(t_i)$ は、各イベントが発生した時刻での強度が高いほどもっともらしい(尤度が高い)ことを表しています。第二項 $-\int_0^T \lambda^*(t)\,dt$ は、イベントが発生しなかった時間帯での強度が低いほどもっともらしいことを表しています。

ポアソン過程($\lambda^*(t) = \lambda$)の場合、対数尤度は

$$ \log L = n \log\lambda – \lambda T $$

となり、最尤推定量は $\hat{\lambda} = n/T$(観測された率)です。

ホークス過程のパラメータ推定も、上の尤度関数を $(\mu, \alpha, \beta)$ で最大化することで行えます。$\Lambda(T) = \int_0^T \lambda^*(t)\,dt$ はホークス過程の場合に解析的に計算できるため、数値最適化が比較的容易です。

モデル診断: 残差分析

推定されたモデルが正しいかどうかの診断には、残差過程(residual process)が使われます。時間変換定理(random time change theorem)により、正しいモデルのもとで変換後のイベント間隔

$$ \xi_i = \int_{t_{i-1}}^{t_i} \lambda^*(t)\,dt $$

は標準指数分布 $\text{Exp}(1)$ に従います。したがって、$\{\xi_i\}$ のQ-Qプロットを指数分布と比較することで、モデルの妥当性を視覚的に診断できます。

この残差分析はKSテスト(コルモゴロフ・スミルノフ検定)と組み合わせることで定量的な検定にもなります。

マーク付き点過程

定義と動機

これまで扱ってきた点過程は、各イベントの「発生時刻」のみを記録していました。しかし、実際の応用では各イベントに付随する情報(マーク)が重要であることが多いです。

  • 地震: 発生時刻に加えて、マグニチュード(震度)がマーク
  • 金融: 取引時刻に加えて、取引量や価格変動がマーク
  • 疫学: 発症時刻に加えて、発症場所(空間座標)がマーク

マーク付き点過程(marked point process)は、イベント発生時刻 $T_n$ とそれに付随するマーク $K_n \in \mathcal{K}$(マーク空間)の組 $\{(T_n, K_n)\}_{n \geq 1}$ として定義されます。

条件付き強度は時刻とマークの同時的な記述に拡張されます。

$$ \lambda^*(t, k) = \lambda_g^*(t) \cdot f^*(k|t) $$

ここで $\lambda_g^*(t) = \int_\mathcal{K} \lambda^*(t, k)\,dk$ はグラウンド強度(マークを周辺化した発生率)、$f^*(k|t)$ は時刻 $t$ でのマークの条件付き分布です。

ETASモデル — 地震学への応用

地震学で広く使われるETASモデル(Epidemic-Type Aftershock Sequence model)は、マーク付きホークス過程の代表例です。

$$ \lambda^*(t, m) = \left[\mu + \sum_{T_i < t} \frac{K e^{\alpha(M_i - M_0)}}{(t - T_i + c)^p}\right] \cdot \frac{\beta e^{-\beta(m - M_0)}}{\text{(GR法則)}} $$

ここで $M_i$ は第 $i$ 番目の地震のマグニチュード(マーク)、$K, \alpha, c, p$ はパラメータ、$\mu$ は背景地震活動の強度です。大きな地震ほど多くの余震を誘発し($e^{\alpha(M_i – M_0)}$)、余震の率はべき乗則で減衰します(大森・宇津の法則: $(t – T_i + c)^{-p}$)。

リニューアル過程

定義

リニューアル過程(renewal process)は、到着間隔 $\tau_1, \tau_2, \dots$ が互いに独立で同一分布(i.i.d.)に従う点過程です。ポアソン過程は到着間隔が指数分布に従う特殊なリニューアル過程です。

到着間隔の分布を $F$ とすると、$n$ 番目のイベント発生時刻は $T_n = \sum_{k=1}^n \tau_k$ であり、計数過程は $N(t) = \max\{n : T_n \leq t\}$ です。

リニューアル過程は独立増分を持ちませんが(到着間隔が指数分布でない場合)、到着間隔が独立同一分布であるため、大数の法則が適用できます。

リニューアル方程式

$m(t) = E[N(t)]$ をリニューアル関数と呼びます。$m(t)$ は次のリニューアル方程式を満たします。

$$ \begin{equation} m(t) = F(t) + \int_0^t m(t – s)\,dF(s) \end{equation} $$

この方程式は「第1のイベントが時刻 $s$ に発生した(確率 $dF(s)$)後、残りの時間 $t – s$ でさらに $m(t – s)$ 個のイベントが発生する」という再帰的な構造を表しています。リニューアル方程式はラプラス変換やフーリエ変換で解けますが、一般には解析解が得られないため数値的に解かれます。

基本リニューアル定理

到着間隔の期待値を $\mu_\tau = E[\tau]$ とすると、基本リニューアル定理(Elementary Renewal Theorem)により、

$$ \begin{equation} \frac{N(t)}{t} \xrightarrow{t \to \infty} \frac{1}{\mu_\tau} \quad \text{a.s.} \end{equation} $$

が成り立ちます。長期的にはイベントの発生率は到着間隔の期待値の逆数に収束します。この結果は大数の強法則の直接的な帰結です。

さらに、$m(t)/t \to 1/\mu_\tau$($t \to \infty$)も成り立ちます。これはリニューアル過程の漸近的な線形成長を保証し、待ち行列理論のトラフィック強度の計算などで使われます。

まとめ

本記事では、点過程の理論をポアソン過程の一般化として体系的に解説しました。

  • 計数過程と到着間隔の等価な記述により、ランダムなイベントの数学的な枠組みを定義しました。計数過程 $N(t)$ はcadlagな単調非減少整数値過程であり、到着間隔 $\tau_n = T_n – T_{n-1}$ と一対一に対応します
  • ポアソン過程は独立定常増分を持つ最も基本的な点過程であり、到着間隔は指数分布に従います。超位置原理と間引き定理がポアソン過程の重要な演算規則です
  • 非定常ポアソン過程は時変強度関数 $\lambda(t)$ により非定常な到着パターンを記述し、Lewis-Shedlerの間引き法で効率的にシミュレーションできます
  • 条件付き強度関数 $\lambda^*(t)$ は過去のイベント履歴 $\mathcal{H}_t$ に依存する最も一般的な点過程の記述方法です。条件付き強度を指定することで点過程の有限次元分布が完全に決まります
  • ホークス過程は自己励起型の点過程であり、各イベントが後続のイベントの発生率を一時的に高めます。安定性条件は分岐過程の絶滅条件と同じ $\alpha < 1$ であり、地震の余震、金融取引のクラスタリング、ソーシャルメディアの情報拡散などのモデル化に使われます
  • 対数尤度は条件付き強度を使って $\sum_i \log\lambda^*(t_i) – \int_0^T \lambda^*(t)\,dt$ と表され、パラメータの最尤推定とモデル選択の基盤となります
  • マーク付き点過程は各イベントに付随する情報(マグニチュード、取引量など)を組み込んだ拡張であり、ETASモデルなどの実用的なモデルの基礎です
  • 補償マルチンゲール $N(t) – \Lambda(t)$ と時間変換定理に基づく残差分析は、モデル診断の強力な道具です
  • リニューアル過程は独立同一分布の到着間隔を持つ点過程であり、基本リニューアル定理により長期的な発生率が $1/\mu_\tau$ に収束します

次のステップとして、以下の記事も参考にしてください。