双类型事件标记下带幂律强度函数的二元非齐次泊松过程(Bivariate Non-Homogeneous Poisson Process, BNHPP)建模与极大似然估计

《MethodsX》：A notes bivariate power law processes: conditional intensity and parameter estimation techniques,

　　

摘要：本研究提出一种基于双类型事件标记的二元非齐次泊松过程(Bivariate Non-Homogeneous Poisson Process, BNHPP)模型，用于描述由同一母非齐次泊松过程(Non-Homogeneous Poisson Process, NHPP)经独立稀释(thinning)产生的两类并发点过程。母过程的强度函数(conditional intensity function)设为幂律形式λ(t)=c·tb（c为基线强度参数，b为时间尺度指数），各事件以概率p被标记为类型1、概率q被标记为类型2（p+q≤1）。研究人员推导了BNHPP的联合概率质量函数及对数似然函数，并给出参数p、q、c、b的极大似然估计(Maximum Likelihood Estimation, MLE)解析式（其中b需数值求解或代入c的估计后显式求解）。通过时间重标定定理(Time Rescaling Theorem)进行拟合优度检验，并将模型应用于印度尼西亚苏拉威西岛2018—2020年地震目录（将M≥5.0记为类型1，其余记为类型2）验证实用性。模拟实验表明所提估计方法随观测窗口增大能有效收敛至真值，时间重标定后的残差服从单位指数分布Exp(1)，Kolmogorov-Smirnov检验无法拒绝模型正确性假设(p=0.2540)。

【论文解读】

研究背景方面，传统齐次泊松过程难以刻画强度随时间变化的事件序列，而非齐次泊松过程(NHPP)可通过时变强度函数描述此类现象。当同一观测区间内的事件可划分为多种互斥类型（如高低幅值地震、故障与正常记录），现有文献多分别建模各子过程而忽略其共享母过程结构。将NHPP经独立稀释得到多类型NHPP已有理论支撑，但针对幂律(power-law)强度λ(t)=c·t^b下二元情况的联合似然推导、参数可识别性及拟合诊断尚缺乏系统阐述，且地震学中常需同时分析全量事件与特定阈值事件的关系。为此，研究人员开展了BNHPP模型构建、MLE推导、模拟验证与真实地震数据应用的研究。结论表明该模型能以闭合形式给出联合分布，MLE估计具有一致性趋势，时间重标定定理可有效检验拟合，应用于地震数据可量化年度震情演化特征。该文发表于《MethodsX》。

关键技术方法：①基于母NHPP λ(t)=c·t^b独立稀释构造类型1（概率p，强度p·c·t^b）与类型2（概率q，强度q·c·t^b）两个独立NHPP，形成BNHPP；②推导观测区间(0,T]内恰发生n₁次类型1与n₂次类型2事件的联合概率及对数似然?(θ)=n₁ln p+n₂ln q+n ln c+b·Σ_i=1ⁿln t_i-(p+q)c·T^b+1/(b+1)；③令偏导??/?p=??/?q=??/?c=??/?b=0解得p? =n₁/(n₁+n₂)，? =n(b+1)/((p+q)T^b+1)，b?需解Σln t_i+cT^b+1/(b+1)²-cT^b+1ln T/(b+1)=0（代入?后可显式得b?=n/(n ln T-(p+q)Σln t_i)-1）；④蒙特卡洛模拟生成不同T下BNHPP数据评估估计偏差；⑤采用时间重标定定理——计算Λ(t_i)=∫_ti-1^t_ic?·s^b?ds=c?(t_i^b?+1-t_i-1^b?+1)/(b?+1)，残差τ_i应~Exp(1)，用Q-Q图与Kolmogorov-Smirnov检验评估；⑥实证数据来自印尼BMKG 2018—2020年苏拉威西岛地震（M<5.0为类型2，M≥5.0为类型1），以日为单位时间尺度分年度拟合。

研究结果：

Bivariate point process model：研究人员证明母NHPP N(t)~NHPP(λ(t)=c·t^b)经独立稀释产生N₁(t)~NHPP(p·c·t^b)与N₂(t)~NHPP(q·c·t^b)，其联合分布为BNHPP，P(n₁,n₂)=[ (Λ₁(T))^n₁/n₁! · (Λ₂(T))^n₂/n₂! ]·exp{-(Λ₁(T)+Λ₂(T))}，其中Λ_k(T)=∫₀^Tp_k·c·s^bds=p_k·c·T^b+1/(b+1)。推论1说明若N₁是强度p·c·t^b的NHPP则补过程N₂是强度(1-p)·c·t^b的NHPP；推论2说明两独立此类NHPP构成BNHPP。

Likelihood function：由无穷小区间乘积形式导出似然L*(θ)=(pc)^n₁(qc)^n₂(∏_i=1ⁿt_i)^b·exp{-(p+q)cT^b+1/(b+1)}，取对数得上述?(θ)。

Parameter estimation technique：对?(θ)求偏导得p?=n₁/n，q?=n₂/n（n=n₁+n₂），?=n(b+1)/((p+q)T^b+1)，b? = n / [n ln T - (p+q)Σ_i=1ⁿln t_i] - 1（代入?后显式）。若p,q未知可用p?,q?迭代或直接联立求解。

Method validation — Simulation of experiment two Poisson process：设定真值c=3,b=3,p=0.4,q=0.6，T=3,5,10,15,20分别仿真。随T增大，估计值p?→0.4，b?→3.0，?→3.0，吻合真值，说明MLE具一致性。T=3时p?=0.375,b?=1.868,?=6.875接近真值。

Goodness of fit assesment using the time rescaling theorem：对T=5仿真数据做时间重标定，τ_i~Exp(1)，Q-Q图贴近对角线，K-S统计量D=0.0119(p=0.2540>0.05)，不拒绝模型正确假设。

Case study：苏拉威西岛2018—2020地震(类型1:M≥5.0;类型2:M<5.0)。年度MLE结果为——2018:c?=0.3366,b?=0.7123,p?=0.0147，λ₁(t)=0.0049·t^0.7123，b>0说明震频后期增强；2019:c?=9.9814,b?=-0.1984,p?=0.0128，b<0说明震频递减；2020:c?=3.3218,b?=0.1609,p?=0.0144，微升趋势。p?稳定在约1.3%—1.5%，反映大震比例稳定。

讨论与结论：作者指出BNHPP由独立稀释构造，故类型间无交叉激励(cross-excitation)，无法刻画余震自激发行为，未来可扩展含互激成分的自激励多元点过程(self-exciting multivariate point process, e.g. Hawkes)。局限性为假设条件独立给定强度函数。研究表明所提BNHPP框架结合MLE与时间重标定诊断可有效建模共享幂律强度母过程的双类型点过程，参数估计简便，适用于地震学及其他具分类标记的非齐次事件流分析。