总结

本文给出了时空点过程的一阶残差分析,类似于Baddeley及其同事开发的空间点过程的残差分析。给出了一阶和二阶残差的示例。特别是,残差分析可以用作模型改进的强大工具。以地震发生的时空流行病学型余震序列模型为基线模型,二阶残差分析可用于识别基线模型中未隐含的数据的许多特征,为我们提供如何构建更好模型的线索。

1.介绍和动机

时间、空间和时空点过程在流行病学、生物学、环境科学和地球科学等许多领域得到了越来越广泛的应用。在相关的统计推断技术中,如模型说明、参数估计、模型选择、拟合优度测试和模型评估,用于测试拟合优度和模型评估的工具相当落后。这是本文的动机之一。

模型选择程序可用于测试拟合优度和模型评估。给定几个适用于同一数据集的显式模型,我们可以使用一些模型选择标准,例如Akaike的信息标准(例如,请参见Akaike's information criteria(1974))和交叉验证(例如Stone(1977))找出其中最好的模型。为了找到一个比当前最佳模型更好的模型,我们总是可以尝试几个新模型的可能版本,将它们适应相同的数据集,然后再次使用模型选择过程来查看其中一个新模型是否成为性能最佳的模型。然而,这些程序并不总是易于实施。制定模型并将其与数据集相匹配可能需要大量的编程和计算任务。最后,模型选择程序只给了我们一些数量,这些数量表明了每个模型的总体拟合程度。很难从这些数量中推断出一个模型(即使它不是最好的)是否比模型选择过程中排名较高的其他模型具有更好的特性。如果可以简化模型改进过程,这将非常有帮助。本文中开发的残差分析可用于此目的。

为了促进这项工作,我们首先描述了研究中使用的数据集。这里的发展与寻找一系列问题的答案有关,这些问题是建立与地震群有关的现象的模型。虽然地震数据来自地球物理学领域,但在传染病的流行病学建模、生物学、生态学和环境科学中也出现了类似的问题。

地震震中并非均匀分布在地球表面。在全球范围内,地震主要发生在板块边界之间的俯冲带。在局部地区,地震沿着活动断层或在火山区积聚。它们的深度从几公里到700公里不等。尽管一次地震可能大到7级,造成巨大的灾难,但大多数地震都很小,只有灵敏的地震仪才能探测到。

地震活动在空间和时间上都是聚集的。地震群之间以及与背景地震活动的重叠使我们的分析变得复杂。对于长期地震预测,即在大约10年的时间尺度内评估强烈地震的风险,需要很好地估计背景地震活动率。相比之下,对于短期预测(以一小时或一天为单位),有必要充分了解地震集群。

地震目录包括一个列表{(t吨,x个,)}和其他相关信息,其中t吨,x个分别记录发生时间、震中位置和震级第个事件。图。1显示了本分析中使用的日本气象局(JMA)目录中的浅层地震(深度小于100 km)。本目录的时间跨度为1926年1月1日至1999年12月31日。本文选择顶点为(134.0)的多边形中的数据E、 31.9日N) ,(137.9E、 33.0岁N) ,(143.1E、 33.2款N) ,(144.9E、 35.2款N) ,(147.8E、 41.3款N) ,(137.8E、 44.2款N) ,(137.4E、 40.2款N) ,(135.1E、 38.0岁N) 和(130.6E、 35.4N) ●●●●。以1926年1月1日后第10000天至1991年12月31日的时间段为目标范围,采用最大似然法进行参数估计。

日本中部浅层地震的位置和时间(○,目标事件,互补事件;不同大小的圆圈代表4.2至8.1级地震:(a)震中(多边形代表目标区域);(b) 震中纬度与发生时间
图1

日本中部浅层地震的位置和时间(○,目标事件;图解的互补事件;不同大小的圆圈代表4.2至8.1级地震:(a)震中(多边形代表目标区域);(b) 震中纬度与发生时间

从图中很容易看出。1地震是成群的。通常,时空流行病型余震序列(ETAS)模型用于描述地震集群(Kagan,1991; 拉斯本,1993; Musmeci和Vere-Jones,1992; 绪方,1988,1998,2004; 尾形等。,2003; .,2002,2004; 控制台和Murru,2001; 慰问.,2003; Helmstetter和Sornette,2003年a,b条; 赫尔姆斯特特.,2003). 在这个模型中,地震活动被分为两个部分:背景和簇。背景地震活动被建模为一个泊松过程,该过程在时间上是平稳的,但在空间上不是均匀的。一旦事件发生,无论它是一个背景事件还是由另一个前一个事件生成(触发),它都会根据特定规则产生(触发)自己的子事件。这种模型是一个具有移民(背景)的连续型分支过程。该模型可以完全由条件强度函数定义(以给定历史为条件的危险率t吨截至当前时间t吨; 参见附录A或者Daley和Vere-Jones(2003),第7章,了解更多详细信息)

λ(t吨,x个,)=Δt吨0,Δx个0,Δ0[公共关系{N个((t吨,t吨+Δt吨]×(x个,x个+Δx个]×(,+Δ])1t吨}Δt吨Δx个Δ].

本文使用的ETAS模型的条件强度函数采用Ogata给出的形式(1998),即。

λ(t吨,x个,)=γ(){u个(x个)+:t吨<t吨κ()(t吨t吨)(f)(x个x个)}
1

哪里

κ()=A类经验(α),
2
γ()=β经验(β)H(H)(),
(f)(x个)=q个1πC类经验(α){1+x个2C类经验(α)}q个

(t吨)=第页1c(c)(1+t吨c(c))第页H(H)(t吨),第页>1,

H(H)是Heaviside函数。在上面,震级分布γ()基于古腾堡-里希特定律(古腾堡和里希特,1956),预期的儿童数量κ()以山中和岛崎为基地(1990)以及时间密度(t吨)基于修改的Omori公式(Utsu,1969)都是地震活动性研究中的经验法则。背景速率和参数A类,α,β,c(c),第页C类可以通过迭代算法(庄. (2002,2004); 参见附录C).

然而,关于这个模型的公式有很多问题。

  • (a)

    后台进程是否静止?

  • (b)

    背景事件和触发事件在量级分布或触发后代方面是否不同?

  • (c)

    触发事件的大小分布是否取决于其父事件的大小?

  • (d)

    应用相同的指数函数exp合理吗(α)在两者中κ()和(f)(x个|)?

在以前的研究中,残差分析是通过将点过程转换为标准泊松过程(Ogata,1988). 勋伯格(2004)使用稀释残差分析ETAS模型与加利福尼亚地震数据的拟合优度。巴德利. (2005)取得了更加显著和普遍的发展。但它们的残差分析方法都是一阶的,远远不足以解决涉及聚类和抑制等二阶特性的问题。为了回答这些问题,我们有必要将残差分析的概念推广到更高阶。

. (2004)开发了一种随机重建方法,以ETAS模型作为参考模型,测试与地震簇相关的上述假设。他们的方法完全基于直觉,而不是严格的理论基础。正如我们在本文中所示,他们的方法可以通过使用残差分析工具进行验证。本文的另一个动机是为随机重建方法提供理论基础,该方法也可以应用于更广泛的点过程模型。

在接下来的章节中,我们首先回顾了Baddeley开发的一阶残差. (2005)然后提出二阶残差原理。通过一些简单的例子,以及通过解决ETAS模型公式中提出的问题来改进模型,说明了这些残差的用途和功效。

2.一阶残差和示例

巴德利. (2005)给出了空间点过程的一阶残差X(X)={x个1,x个2,…}基于Nguyen和Zessin(1979)公式

E类[x个X(X)D类小时(x个;X(X)\{x个})]=E类[D类小时(x个;X(X))λ(x个;X(X))μ(d日x个)],

对于任何可测集合D类,其中λ是Papangelou条件强度。因此,他们定义创新关于(D类,小时)将成为

V(V)(D类,小时,λ)=x个X(X)D类小时(x个;X(X)\{x个})D类小时(x个;X(X))λ(x个;X(X))μ(d日x个)

残余沉积物对应于V(V)(D类,小时,λ)

R(右)(D类,小时^,λ^)=x个X(X)D类小时^(x个;X(X)\{x个})D类小时^(x个;X(X))λ^(x个;X(X))μ(d日x个),

哪里λ^是拟合的条件强度。如果小时也取决于参数θ在里面λ,小时^通过替换估计参数获得θ^在里面小时.

中定义的时间或时空点过程的条件强度附录A比空间点过程的Papangelou条件强度更简单。在前一种情况下,事件在某个时间或时空位置的发生仅取决于早期发生的事件,而在后一种情况中,事件在特定位置的发生取决于其他地方发生的所有其他事件。原则上,我们可以简单地获得方程中定义的条件强度(1)和附录A通过对Papangelou条件强度的期望,这被定义为基于σ-由当前时间和地点以外的所有事件生成的代数σ-直到当前时间的观测历史的代数,以定义时空点过程的一阶残差。然而,用这种方法将一阶残差推广到二阶残差是困难的,这是本文的主要关注点。我们通过鞅理论的观点利用了这个过程的演化特征。

N个是区间[0,]和ad日-维度区域X(X)⊂ℝd日,承认条件强度λ(t吨,x个)(见Daley和Vere-Jones(1988),第13章,或附录A条件强度的定义)。根据条件强度的鞅性质,对于任何可预测过程(相对于σ-由形式集生成的代数(,t吨B类×E类哪里E类∈ ℱB类∈ 𝒳)小时(t吨,x个)0、几乎无处不在),

E类[D类小时(t吨,x个)N个( d日t吨×d日x个)]=E类[D类小时(t吨,x个)λ(t吨,x个)μ(d日t吨×d日x个)],D类X(X)
4

哪里μ是勒贝格测量×d日(另请参见布瑞莫德(1981)第2章和卡尔(1991),第5章)。

使用巴德利采用的术语. (2005)对于空间点处理,定义一级创新关于可预测函数小时(t吨,x个)0,几乎处处可见,并且是一个可测量的集合D类×X(X)通过

V(V)(D类,小时,λ)=D类{小时(t吨,x个)N个( d日t吨×d日x个)小时(t吨,x个)λ(t吨,x个)μ(d日t吨×d日x个)}.
5

因为方程式(4),E类[V(V)(D类,小时,λ)]=0. 这个一阶残差通过替换方程中的真实模型来定义(5)通过拟合模型,即。

R(右)(D类,小时^,λ^)=D类{小时^(t吨,x个)N个( d日t吨×d日x个)小时^(t吨,x个)λ^(t吨,x个)μ(d日t吨×d日x个)}.

如果拟合模型与真实模型非常接近,则R(右)(D类,小时^,λ^)0.

以下六个示例给出了不同类型的一级创新。由于与这些创新相对应的残差可以很容易地定义,因此我们只给出了创新的形式。

2.1. 示例1:原始创新

小时(t吨,x个)=1; 那么相应的创新是

V(V)(D类,小时,λ)=N个(D类)D类λ(t吨,x个)μ(d日t吨×d日x个),

它被称为原始创新应用中原始创新或残差的一个重要特性是,为了实现过程{(t吨,x个):N个},

{τ=0t吨X(X)λ(t吨,x个)d日( d日x个)( d日t吨)}

是一个标准的泊松过程(Meyer,1971; 帕潘格鲁,1972; 绪方,1988; Vere-Jones和Schoenberg,2004).

2.2. 示例2:互惠λ创新(巴德利.,2005; 勋伯格,2004)

如果小时(t吨,x个)=1/λ(t吨,x个),

V(V)(D类,小时,λ)=D类N个( d日t吨×d日x个)λ(t吨,x个)μ(D类).

这是Stoyan和Grabarnik的类似物(1991)吉布斯点过程的权重。Schoenberg时空ETAS模型的残差分析(2004)基本上是基于互惠λ创新或残差。

2.3. 示例3:培生创新(Baddeley.,2005)

如果小时(t吨,x个)=1/√λ(t吨,x个),

V(V)(D类,小时,λ)=D类N个( d日t吨×d日x个)/λ(t吨,x个)D类λ(t吨,x个)μ(d日t吨×d日x个).

2.4. 示例4:分数创新(Baddeley.,2005)

这个得分创新由定义

V(V)(D类,[日志{λ(t吨,x个)}]θ,λ)=D类[日志{λ(t吨,x个)}]θN个( d日t吨×d日x个)D类λ(t吨,x个)θμ(d日t吨×d日x个),

哪里θ是模型中的任何常规参数。分数残差方程

R(右)(D类,[日志{λ^(t吨,x个)}]θ,λ^)=0

是最大化对数似然函数的条件D类,即。

日志(L(左))=D类日志{λ(t吨,x个)}N个( d日t吨×d日x个)D类λ(t吨,x个)μ(d日t吨×d日x个).

2.5. 示例5:加权得分创新和局部最大似然估计

创新的定义是

V(V)(D类,w个(t吨,x个;t吨0,x个0)[日志{λ(t吨,x个)}]θ,λ)                                                              =D类w个(t吨,x个;t吨0,x个0)[日志{λ(t吨,x个)}]θ{N个( d日t吨×d日x个)λ(t吨,x个)μ(d日t吨×d日x个)}

被称为加权得分创新,其中w个(t吨,x个;t吨0,x个0)是以给定的(t吨0,x个0). 加权分数残差上的等式

R(右)(D类,w个(t吨,x个;t吨0,x个0)[日志{λ^(t吨,x个)}]θ,λ^)=0

是最大化局部对数似然的条件

WLL(WLL)=D类w个(t吨,x个;t吨0,x个0)日志{λ(t吨,x个)}N个( d日t吨×d日x个)D类w个(t吨,x个;t吨0,x个0)λ(t吨,x个)μ(d日t吨×d日x个).

2.6. 例6:检验流行病学型余震序列模型中背景过程的平稳性

定义背景创新ETAS模型

V(V)(D类,φ,λ)=φ(t吨,x个,)D类(t吨,x个,)D类u个(x个)γ()d日t吨 d日x个 d日,
6

哪里ϕ(t吨,x个,)=u个(x个)γ()/λ(t吨,x个,),λ(t吨,x个,)和u个(x个)在方程式中给出(1),D类是以下项的可测量子集×X(X)×M(M)D类是的指数函数D类。要测试后台进程的平稳性,请选择D类(t吨)=(0,t吨B类x个×B类,B类x个X(X)B类S公司.来自E类[V(V)(D类,ϕ,λ)]=0,

V(V)(D类,φ,λ)(t吨)=φ(t吨,x个,)D类(t吨)(t吨,x个,)t吨B类x个u个(x个)μ(d日x个)B类γ()μ(d日),
7

它是关于的线性函数t吨.此类测试和地球物理解释的应用可在壮族地区找到. (2005)Hainzl和Ogata(2005).

3.二阶残差分析

给定二阶ℱ-可预测函数H(H)(t吨,x个;t吨,x个),即可以从表单中的流程生成的流程小时1(t吨,x个)小时2(t吨,x个),其中小时1小时2是一阶可预测过程,通过线性组合和单调极限(详见附录B),的二阶创新关于H(H)(t吨,x个;t吨,x个)的D类∈(𝒯𝒳)为

V(V)2(D类,H(H),λ)=D类\诊断2(D类)H(H)(t吨,x个;t吨,x个)N个( d日t吨×d日x个)N个( d日t吨×d日x个)                                   D类H(H)(t吨,x个;t吨,x个)λ(t吨,x个)λ(t吨,x个)μ(d日t吨×d日x个)μ(d日t吨×d日x个),
8

其中诊断2(D类)={(t吨,x个,t吨,x个) ∈D类}. 如所示附录B,E类[V(V)2(D类,H(H),λ)]=0. 这个二阶余量可以由定义R(右)2(D类,H(H)^,λ^)=V(V)2(D类,H(H)^,λ^).

3.1. 例7:一阶创新的方差

对于一阶可预测过程小时(t吨,x个)以及任何B类∈ 𝒯⊗𝒳,

无功功率,无功功率{B类小时(t吨)N个( d日t吨)}=E类[{B类小时(t吨)N个( d日t吨)}2]E类[B类小时(t吨)N个( d日t吨)]2=E类[B类×B类小时(t吨)小时()N个( d日t吨)N个( d日)]E类[B类小时(t吨)N个( d日t吨)]2=E类[B类×B类小时(t吨)小时()λ(t吨)λ()μ(d日t吨)μ(d日)]+E类[B类小时2(t吨)λ(t吨)μ(d日t吨)]E类[B类小时(t吨)λ(t吨)μ(d日t吨)]2=E类[B类小时2(t吨)λ(t吨)μ(d日t吨)],

无功功率,无功功率{B类小时(t吨)N个( d日t吨)}=E类[B类小时2(t吨)N个( d日t吨)].

该公式可对一阶残差的标准误差进行近似估计。在示例6中,如果ETAS模型与真实模型足够接近,则背景残差的标准误差R(右)(D类,φ^,λ^)(t吨)偏离预期的线性函数t吨近似值(7)可以近似为

{φ^2(t吨,x个,)D类(t吨)(t吨,x个,)}1/2.

4.二阶残差在地震数据中的应用:模型改进

现在我们有了工具来解决第节中列出的与地震集群相关的问题1在本节中,在测试这些假设之前,我们简要回顾了庄提出的随机重建方法的基本思想. (2004).

方程中的条件强度函数(1)经估计,它为我们提供了一种很好的方法来评估事件是否可能是背景事件或由其他事件触发的概率(Kagan和Knopoff,1981; .,2002). 考虑地震发生时背景地震活动率相对于总地震活动率的贡献第次活动,

φ=u个(x个)γ()λ(t吨,x个,).
9

如果我们删除概率为1−的第th个事件ϕ对于流程中的所有事件,我们可以实现具有发生率的流程μ(x个)γ()(另请参阅绪方(1981)和卡尔(1991),第5章,理由)。因此,很自然地ϕ作为概率该事件是一个背景事件。同样,

ρj个=γ(j个)(t吨j个t吨)(f)(x个j个x个)λ(t吨j个,x个j个,j个)
10

可以视为j个事件直接由第个事件,即给定固定,如果我们保留每个事件j个有可能ρij公司,我们可以得到强度的泊松过程γ()(t吨t吨)(f)(x个x个),这是由目录中的第个事件。基于这些想法,庄. (2004)通过用估计概率加权的事件建立经验函数,检验了与地震聚类特征相关的各种假设φ^ρ^j个根据ETAS模型的拟合条件强度。这种重建方法完全基于直觉。在下面的部分中,我们将展示它的工作原理。

为了简化起见,在下面的章节中,我们使用了黎曼积分的符号,毫无疑问,我们假设积分在数据对应的整个范围内。

4.1. 改善个人功能

针对特定数据集,使用条件强度函数建立点过程模型主要是从经验结果和研究人员的经验开始的。一旦给出了模型的显式形式,第一个问题是这种形式是否合适,或者是否可以以某种方式加以改进。在本小节中,我们将使用κ()在方程式中(1)作为一个简单的例子来展示如何改进模型公式中的个体函数。

考虑ETAS模型的条件强度的以下形式:

λ1(t吨,x个,)=γ()[u个(x个)+:t吨<t吨K(K)()(t吨t吨)(f)(x个x个)],

哪里K(K)()假设是比κ()对于重大事件的平均儿童人数.尽管如此K(K)()未知,我们将在下面说明,我们仍然可以构建一个所谓的“无偏比率”(两个无偏估计的比率)估计K(K)()基于二阶新息或残差的性质。

H(H)j个=K(K)()(t吨j个t吨)(f)(x个j个x个)γ(j个)λ1(t吨j个,x个j个,j个).

然后根据方程式(8)

E类[,j个([0δ,0+δ])H(H)j个]=E类[([0δ,0+δ])      ×K(K)()(t吨t吨)(f)(x个x个)γ()λ1(t吨,x个,)d日t吨d日x个d日d日t吨 d日x个 d日]=K(K)()γ()([0δ,0+δ])d日×E类[λ1(t吨,x个)d日t吨 d日x个]2δK(K)(0)γ(0)E类[λ1(t吨,x个)d日t吨 d日x个],
11

哪里

λ1(t吨,x个)=u个(x个)+:t吨<t吨K(K)()(t吨t吨)(f)(x个x个).

使用一阶创新公式,

E类[,j个([0δ,0+δ])]=E类[([0δ,0+δ])λ1(t吨,x个,)d日t吨 d日x个 d日]=γ()([0δ,0+δ])d日×E类[λ1(t吨,x个)d日t吨 d日x个]2δγ(0)E类[λ1(t吨,x个)d日t吨 d日x个].
12

除法表达式(11)按表达式(12)收益率

K(K)^(0),j个H(H)j个([0δ,0+δ])([0δ,0+δ]).
13

尽管K(K)()未知,我们可以设置K(K)^(0)()=κ^()因此H(H)j个(0)=ρ^j个第一步,其中κ^()来自ETAS模型的拟合条件强度ρ^j个是估计的ρij公司在方程式中定义(10). 然后我们可以使用近似值(13)来获得更新的K(K)(),即。

K(K)^(1)(0)=,j个ρ^j个([0δ,0+δ])([0δ,0+δ]).
14

Σj个ρ^j个可以解释为来自第个事件,方程式的右侧(14)可以被视为母亲所生孩子的平均数量0。此过程可以迭代直至收敛。

在这个应用中,二阶创新在保证近似(11)中数量的估计与K(K)()γ(). 将ETAS模型拟合到地震数据的估计值仅用作近似值中的初始值(13)。这种程序在精神上类似于期望最大化算法,期望步长由二阶新息或残差的特性执行,最大化步长由近似中的比率无偏估计的非参数方法取代(13)。

4.1.1. 例8:日本气象厅目录中地震产生的儿童平均数量

我们将上述重建程序应用于JMA目录,重建结果如图所示。2很明显K(K)^(1)()κ^()可以忽略不计,这意味着没有必要改进κ()在这个阶段。我们将讨论图。2在第节示例9中进一步说明4.2.1.

重建所有事件K^(1)(≠)的触发能力(由相同量级事件触发的平均儿童数),背景事件K^(1)(×)和触发事件κ^1(1)
图2

重建所有事件的触发能力(由相同大小的事件触发的平均儿童数)K(K)^(1)(⑪),背景事件K(K)^(1)(×)和触发事件κ^1(1)()(+),如方程式中所示(14), (20)和(21)分别(庄.,2004)(水平轴上的4.2级对应于=0):--,首字母K(K)^(0)()=κ^0(0)()=κ^1(0)()=κ^()通过拟合模型(1)到JMA目录的最大似然估计获得

4.2. 背景事件和触发事件之间的差异

ETAS模型假设背景事件和触发事件之间没有区别。一旦一个事件发生,它的震级来自于共同的独特震级分布,它会以与其他事件相同的方式触发后代。测试背景事件与触发事件是否不同是有趣且重要的。在本节中,我们将测试他们触发儿童的行为是否不同。

考虑一个更复杂的模型,其中包含一个条件强度

λ1(t吨,x个,)=λ1(t吨,x个,,ω)(ω=0)+λ1(t吨,x个,,ω)(ω=1),
15

哪里

λ1(t吨,x个,,ω)={u个(x个)γ0(),如果ω=0,γ1()t吨<t吨ξ(t吨,x个;t吨,x个,,ω),如果ω=1,
16

ξ(t吨,x个;t吨,x个,,ω)={κ0()0(t吨t吨)(f)0(x个x个),如果ω=0,κ1()1(t吨t吨)(f)1(x个x个),如果ω=1,
17

即,事件是背景事件,如果ω=0或触发事件,如果ω=1.在这个混合类型模型中,背景事件(t吨,x个,)有一个量级来自γ0()和触发器平均值κ0()儿童,其数量、时间和位置分布的特点是γ1,0(f)0分别是。对于触发的事件(t吨,x个,),大小来自密度γ1()它会触发κ1()儿童的平均数量、时间和地点密度为γ1,1(f)1分别是。

构造以下各项的比率无偏估计κ0,让

H(H)(t吨,x个,;t吨,x个,)=γ1()ξ(t吨,x个;t吨,x个,,ω)λ1(t吨,x个,,ω)(ω=0)λ1(t吨,x个,)λ1(t吨,x个,)=γ1()κ0()0(t吨t吨)(f)0(x个x个)u个(x个)γ0()λ1(t吨,x个,)λ1(t吨,x个,)

小时(t吨,x个,)=λ1(t吨,x个,,ω)(ω=0)λ1(t吨,x个,)=u个(x个)γ0()λ1(t吨,x个,).

然后

E类[,j个H(H)(t吨,x个,;t吨j个,x个j个,j个)([0δ,0+δ])]                                               =E类[H(H)(t吨,x个,;t吨,x个,)([0δ,0+δ])d日t吨d日x个d日d日t吨 d日x个 d日]                                               =κ0()γ0()([0δ,0+δ])d日u个(x个)d日t吨 d日x个                                               2δκ0(0)γ0(0)u个(x个)d日t吨 d日x个
18

E类[小时(t吨,x个,)([0δ,0+δ])]=E类[小时(t吨,x个,)([0δ,0+δ])d日t吨 d日x个 d日]=γ0()([0δ,0+δ])d日u个(x个)d日t吨 d日x个
2δγ0(0)u个(x个)d日t吨 d日x个.
19

比率无偏估计κ0()由近似值(18)和(19)给出

κ^0(0)=,j个H(H)(t吨,x个,;t吨j个,x个j个,j个)([0δ,0+δ])小时(t吨,x个,)([0δ,0+δ]).

如果在初始阶段,我们让γ0(0)=γ1(0)=γ^,κ0(0)=κ1(0)=κ^,(f)0(0)=(f)1(0)=(f)^0(0)=1(0)=^,其中γ^,κ^,(f)^^从ETAS模型的拟合条件强度,我们可以得到

κ^0(1)(0)=,j个ρ^j个φ^([0δ,0+δ])φ^([0δ,0+δ]),
20

哪里φ^ρ^j个估计值是ϕρij公司在方程式中定义(9)和(10)分别是。

同样,我们可以发现

κ^1(0)=,j个H(H)(t吨,x个,;t吨j个,x个j个,j个)([0δ,0+δ]){1小时(t吨,x个,)}([0δ,0+δ]),

哪里

H(H)(t吨,x个,;t吨,x个,)=γ1()ξ(t吨,x个;t吨,x个,,ω)λ1(t吨,x个,,ω)(ω=1)λ1(t吨,x个,)λ1(t吨,x个,)=γ1()κ1()1(t吨t吨)(f)1(x个x个)λ1(t吨,x个,){1u个(x个)γ0()λ1(t吨,x个,)}

κ^1(1)(0)=,j个ρ^j个(1φ^)([0δ,0+δ])(1φ^)([0δ,0+δ])
21

方程中的其他函数(15)–(17)例如1,2,(f)1(f)2可以用类似的方式重建。

4.2.1. 例9:日本气象局目录中背景事件和触发事件触发的平均儿童人数差异

将上述程序应用于JMA目录,图中绘制了背景事件和每个量级触发事件的触发能力(平均儿童数)。2发现背景事件和触发事件根据不同的正指数规律生成子事件。然而,与同等规模的触发事件相比,背景事件平均触发的子事件更少。这种特征在评估与前震有关的概率方面起着重要作用(每个前震的背景事件都至少有一个较大的后代)(参见庄和绪方的详细讨论(2006)).

4.3. 空间比例因子和触发能力

如第节所述1,关于模型(3)的使用的直接问题包括以下内容。

  • (a)

    是比例因子C类经验(α)必要的,即可以用常量替换吗D类0?

  • (b)

    比例因子应该C类经验(α)具有相同的指数α比如模型(2)的触发能力?

  • (c)

    比例因子是指数定律的形式吗?

假设真实条件强度为

λ1(t吨,x个,)=γ(){u个(x个)+:t吨<t吨κ()(t吨t吨)(f)*(x个x个)},

哪里

(f)*(x个)=(f)*{x个,σ()}                      =π(q个1)σ(){1+x个2σ()}q个

是子代位置的概率分布函数σ()是…的函数.考虑

H(H)(t吨,x个,;t吨,x个,)=γ()κ()(t吨t吨)(f)*{x个x个,σ()}λ1(t吨,x个,)(日志[(f)*{x个x个,σ()}])σ(),

然后,对于任何正数δ,

E类[,j个H(H)(t吨,x个,;t吨j个,x个j个,j个)([0δ,0+δ])]                    =E类[γ()κ()(t吨t吨)(f){x个x个,σ()}σ()                          ×([0δ,0+δ])λ1(t吨,x个,)d日t吨d日x个d日d日t吨 d日x个 d日]                    =κ()([0δ,0+δ])γ()d日×E类[λ1(t吨,x个)]σ()(f)*{x个x个,σ()}d日x个                    =0,

自从(f)*集成到1。因此,对于所有人来说0,如果δ足够小,

,j个H(H)(t吨,x个,;t吨j个,x个j个,j个)([0δ,0+δ])0

此外,该等式给出了

,j个γ(j个)κ()(t吨j个t吨)(f)*{x个j个x个,σ()}λ1(t吨j个,x个j个,j个)([0δ,0+δ])                                                              ×(日志[(f)*{x个j个x个,σ()}])σ()|σ()=σ()0

要更新σ(),我们将此等式转换为

,j个γ(j个)κ()(t吨j个t吨)(f)*{x个j个x个,σ(0)()}λ1(t吨j个,x个j个,j个)([0δ,0+δ])×(日志[(f)*{x个j个x个,σ()}])σ()|σ()=σ(1)()0
22

σ(0)()=C类经验(α)在这个方程式中,即。(f)*{x个,σ(0)()}=(f)^,其中(f)^是估计的(f)在方程式中定义(),近似值(22)给出

,j个ρ^j个([0δ,0+δ])(日志[(f)*{x个j个x个,σ()}]σ()|σ()=σ(1)(0)=0
23

解决方案σ(1)(0)对于每个0在方程式中(23)更新的版本来自σ(0)(见庄. (2004)关于如何迭代求解这个方程)。这种估计可以被视为最大化伪长似然

锁相环{σ(0)}=,j个ρ^j个日志[(f)*{x个j个x个,σ()}]([0δ,0+δ])

哪里δ是一个小正数。

4.3.1. 例10:日本气象局目录中集群位置的比例定律

使用原始模型和配备的模型将上述程序应用于JMA目录

(f)1(x个)=q个1πC类经验(α)经验{1+x个2C类经验(α)}q个,
24

我们得到了如图所示的结果。为了进行比较,同样的程序也适用于根据模型生成的模拟目录。可以看出,后代空间位置的比例律仍然是指数律,但与触发能力的比例律不同。这一结论导致了以方程式形式对时空ETAS模型公式的新修订(24)(详见绪方图和庄中的详细讨论(2006)和庄. (2005)).

根据母地震的相应震级(4.2级对应s=0)重建σ^(1)(s)()结果:(a)JMA目录原始模型的重建结果;(b) 与(a)相同,但目录是对原始模型的模拟;(c) 根据方程(24)(--,c exp(αs);---C exp(α′s));(d) 与(c)相同,但目录是配备有方程(24)的模型的模拟(--,Cexp(αs);--,C exp(α′s))
图3

的重建结果σ^(1)()对应震级(震级4.2对应于=0)母地震:(a)根据JMA目录的原始模型重建结果;(b) 与(a)相同,但目录是对原始模型的模拟;(c) 基于方程模型的重建结果(24) (——,C类经验(α); - - - - - -C类经验(α)); (d) 与(c)相同,但目录是对配备方程式的模型的模拟(24) (——,C类经验(α); - - - - - -,C类经验(α))

5.结论

利用鞅理论的观点,可以推广时空点过程的残差分析思想。对于一阶新息或残差,如示例所示,几乎可以说,与条件强度相关的每个鞅都对应一种特定的残差。将此思想扩展到高阶情况,当关注二阶特性(如聚类或抑制)时,二阶残差分析对于改进模型公式特别有用。本文使用具有非均匀背景率的时空ETAS模型作为描述地震事件聚类特征的初始模型。一阶和二阶残差分析有助于我们找到改进模型公式的适当方向。

致谢

这项工作是在作者得到日本科学促进会资助的博士后项目P04039和美国国家科学基金会拨款EAR 0409890的支持下进行的。作者感谢Ogata Yosihiko教授和David Vere Jones教授向他介绍了点过程这一有趣的研究领域。阿德里安·巴德利教授亲切地寄给了他关于空间点过程残差分析的手稿,这激发了这项研究。事实证明,与Daryl Daley、Estate Khmaladze、Frederic Schoenberg和David Vere-Jones的讨论非常有帮助,联合编辑以及两名裁判的鼓励和评论也很有帮助。

工具书类

Akaike公司
,
小时。
(
1974
)
统计模型识别的新视角
.
IEEE传输。自动。控制
,
19
,
716
723
.

巴德利
,
答:。
,
特纳
,
R。
,
莫勒
,
J。
哈泽尔顿
,
M。
(
2005
)
空间点过程的残差分析(讨论)
.
J.R.统计。Soc.B公司
,
67
,
617
666
.

布瑞莫德
,
第页。
(
1981
)
点过程和队列:鞅动力学
.
纽约
:
施普林格
.

慰问
,
R。
默鲁
,
M。
(
2001
)
一个简单且可测试的地震聚类模型
.
《地球物理学杂志》。物件。
,
106
,
8699
8711
.

慰问
,
R。
,
默鲁
,
M。
伦巴第
,
上午。
(
2003
)
细化地震聚类模型
.
《地球物理学杂志》。物件。
,
108
,
2468
(doi:).

戴利
,
D。
Vere-Jones公司
,
D。
(
2003
)
点过程理论导论
.
纽约
:
施普林格
.

古腾堡
,
B。
里希特
,
成本加运费。
(
1956
)
地震震级、强度、能量和加速度
.
牛市。地震。美国南部。
,
46
,
105
145
.

海因策
,
秒。
尾形
,
年。
(
2005
)
利用统计地震模型检测地震活动数据中的流体信号
.
《地球物理学杂志》。物件。
,
110
,(doi:).

赫尔姆斯特特
,
答:。
梭奈特
,
D。
(
2003年a
)
触发地震活动级联解释的前震
.
《地球物理学杂志》。物件。
,
108
,
2457
(doi:).

赫尔姆斯特特
,
答:。
梭奈特
,
D。
(
2003年3月
)
相互作用触发地震活动的ETAS模型的可预测性
.
《地球物理学杂志》。物件。
,
108
,
2482
(doi:).

赫尔姆斯特特
,
答:。
,
梭奈特
,
D。
格拉索
,
J.-R.公司。
(
2003
)
主震是有条件前震的余震:前震统计特性如何从余震定律中显现出来
.
《地球物理学杂志》。物件。
,
108
,
2046
(doi:).

卡根
,
年。
(
1991
)
地震目录的可能性分析
.
地球物理学。J.国际。
,
106
,
135
148
.

卡根
,
年。
诺波夫
,
L。
(
1981
)
地震目录的随机综合
.
《地球物理学杂志》。物件。
,
86
,
2853
2862
.

卡尔
,
A.F.公司。
(
1991
)
点过程及其统计推断
,第2版。
纽约
:
德克尔
.

迈耶
第页。
(
1971
)
演示simplifie d’un thee orème de Knight
.
勒克特。数学笔记。
,
191
,
191
195
.

马斯梅基
,
F、。
Vere-Jones公司
,
D。
(
1992
)
历史地震的时空聚类模型
.
Ann.Inst.Statist公司。数学。
,
44
,
1
11
.

阮(Nguyen)
,
十、十、。
泽辛
,
小时。
(
1979
)
吉布斯过程的积分和微分特征
.
数学。纳赫特。
,
88
,
105
115
.

尾形
,
年。
(
1981
)
关于点过程的Lewis模拟方法
.
IEEE传输。通知。理论
,
27
,
23
31
.

尾形
,
年。
(
1988
)
地震发生统计模型及点过程残差分析
.
《美国统计杂志》。助理。
,
83
,
9
27
.

尾形
,
年。
(
1998
)
地震发生的时空点过程模型
.
Ann.Inst.统计。数学。
,
50
,
379
402
.

尾形
,
年。
(
2004
)
区域地震活动的时空模型和晶体应力变化检测
.
《地球物理学杂志》。物件。
,
109
,B03308(doi:).

尾形
,
年。
,
桂树
,
英国。
种村
,
M。
(
2003
)
非均匀地震时空发生建模及其残差分析
.
申请。统计师。
,
52
,
499
509
.

尾形
,
年。
,
J。
(
2006
)
时空ETAS模型及其改进扩展
.
构造物理学
,没有。
413
,
13
23
.

帕潘格鲁
,
F、。
(
1972
)
期望增量点过程的可积性和相关的标度随机变化
.
事务处理。梅尔。数学。Soc公司。
,
165
,
483
506
.

Rathbun公司
,
S.L.公司。
(
1993
)
建模标记的时空点模式
.
牛市。国际统计。仪器。
,
55
,第2册,
379
396
.

勋伯格
,
F.P.公司。
(
2004
)
地震发生点过程模型的多维残差分析
.
《美国统计杂志》。助理。
,
98
,
789
795
.

石头
,
M。
(
1977
)
支持和反对交叉验证的渐近
.
生物特征
,
64
,
29
35
.

斯托扬
,
D。
格拉巴尼克
,
第页。
(
1991
)
与Gibbs点过程相关的随机结构的二阶特征
.
数学。纳赫特。
,
151
,
95
100
.

乌苏
,
T。
(
1969
)
余震和地震统计(I):表征余震序列及其相互关系的一些参数
.
J.工厂。科学。北海道第七大学
,
,
129
195
.

Vere-Jones公司
,
D。
勋伯格
,
F.P.公司。
(
2004
)
重新缩放标记点进程
.
澳大利亚。新西兰。J.统计。
,
46
,
133
143
.

山中县
,
年。
岛崎
,
英国。
(
1990
)
余震次数与主震规模的标度关系
.
《物理学杂志》。地球
,
38
,
305
324
.

,
J。
,
,
C.-P.公司。
,
尾形
,
年。
,
Y.-I。
(
2005
)
用点过程模型研究台湾地区地震活动的背景和集群
.
《地球物理学杂志》。物件。
,
110
,B05S18(doi:).

,
J。
尾形
,
年。
(
2006
)
地震群中最大事件的概率分布特征及其对前震的影响
.
物理。版本E
,
73
,046134(doi:).

,
J。
,
尾形
,
年。
Vere-Jones公司
,
D。
(
2002
)
时空地震事件的随机分布
.
《美国统计杂志》。助理。
,
97
,
369
380
.

,
J。
,
尾形
,
年。
Vere-Jones公司
,
D。
(
2004
)
利用随机重建分析地震聚类特征
.
《地球物理学杂志》。物件。
,
109
,B05301(doi:).

附录A:点过程和条件强度

考虑由事件组成的时空点过程{t吨}间隔[0,]和相应位置{x个}在一个地区X(X)R(右)d日假设边际时间点过程为有秩序的,即发生多个事件的概率[t吨,t吨+δB类o个(δ)为所有人t吨0和有界B类X(X)这种空间-时间点过程可以通过随机计数方法定义N个关于[0,∞)×X(X).给,N个(C类×B类)是一个区域内的事件数B类∈𝒳,有时在集合中C类∈\119983;,其中\119987;和𝒯是Borelσ-子集代数X(X)设Φ是局部有界计数测度的集合N个关于[0,∞)×X(X)并且成为σ-Φ上的代数由{N个((第页,t吨B类):B类∈ 𝒳,0第页t吨}. 然后是一个空间-时间点过程N个是概率空间(Ω,ℱ,)在(Φ,𝒜)上。

对于任何可测量的B类∈\119987,有一个ℱ-补偿器A类(t吨,B类)这样的话N个([0,t吨B类)−A类(t吨,B类)是一个ℱ鞅。μμd日表示勒贝格测度R(右)R(右)d日分别是。假设,对于每个x个X(X),有一个可积的、非负的、适应过程λ(t吨,x个)这样,概率为1t吨R(右)+B类∈ 𝒳,

B类0t吨λ(,x个)μ(d日)μd日( d日x个)=A类(t吨,B类).

存在这种过程的条件和理由可以在Vere Jones和Schoenberg中找到(2004). 如果存在,λ(t吨,x个)被称为ℱ-条件强度,满足

B类λ(t吨,x个)μd日( d日x个)=Δt吨0{E类[N个(t吨+Δt吨,B类)N个(t吨,B类)t吨]Δt吨}.

附录B:二阶可预测性和二阶残差

呼叫潜艇-σ-由以下形式的集合生成的代数(,t吨U型×(,t吨V(V)theℱ⊗\8497]-可预测的σ-代数,表示为Ψℱ⊗ℱ,其中U型∈ ℱ,0t吨<∞,V(V)∈ ℱ,0t吨<∞,𝒯是Borelσ-上的代数R(右)+.一个过程H(H)(t吨,ω,t吨,ω)就是可预测的当它是Ψ时ℱ⊗ℱ可测量,即A类∈ 𝒯, {(t吨,ω,t吨,ω):H(H)(t吨,ω,t吨,ω) ∈A类} ∈ Ψℱ⊗ℱ换句话说,所有可预测的过程都可以从形式的过程中生成

(,b条](t吨)(,b条](t吨)U型(ω)V(V)(ω),U型,V(V),

通过线性组合和单调极限。我们称之为过程小时(t吨,t吨,ω)=H(H)(t吨,ω,t吨,ω)二阶可预测如果H(H)(t吨,ω,t吨,ω)是可以预测的。

给定一个可预测的过程H(H)(t吨,ω,t吨,ω),很容易验证该过程可以分解为

H(H)(t吨,ω,t吨,ω)=H(H)0(t吨,ω,t吨,ω)+H(H)(t吨,ω,t吨,ω)+H(H)+(t吨,ω,t吨,ω)
25

哪里H(H)0(t吨,ω,t吨,ω)=H(H)(t吨,ω,t吨,ω){t吨}(t吨)、和H(H)(t吨,ω,t吨,ω)=H(H)(t吨,ω,t吨,ω)(t吨,∞](t吨)在以下方面是可预测的t吨ω什么时候t吨ω是固定的,并且H(H)+(t吨,ω,t吨,ω)=H(H)(t吨,ω,t吨,ω)(t吨′,∞](t吨)在以下方面是可预测的t吨ω什么时候t吨ω都是固定的。

引理1。使用符号H(H)H(H)+在方程式中(25)对于可预测过程H(H)(t吨,ω,t吨,ω),给定一个适应的点过程,N个(t吨,ω)带ℱ-条件强度λ(t吨,ω),定义

F类(t吨,ω,t吨,ω)E类[(t吨,)H(H)(t吨,ω,,ω)N个( d日,ω)t吨]
26

G公司(t吨,ω,t吨,ω)E类[(t吨,)H(H)+(,ω,t吨,ω)N个( d日,ω)t吨].

然后

F类(t吨,ω,t吨,ω)=E类[(t吨,)H(H)(t吨,ω,,ω)λ(,ω)d日t吨]
27

G公司(t吨,ω,t吨,ω)=E类[(t吨,)H(H)+(,ω,t吨,ω)λ(,ω)d日t吨],
28

假设方程右边的积分(27)和(28)存在。此外,F类(t吨,ω,t吨,ω)和G公司(t吨,ω,t吨,ω)都是可预测的过程。

证明。对于F类在方程式中(26),我们只需要在以下情况下考虑t吨>t吨。对于固定t吨以及任何u个>t吨,H(H)(t吨,ω,u个,ω)(t吨′,∞)(u个)在以下方面是可预测的u个.基于与条件强度相关联的鞅性质,

J型(u个)=(0,u个]H(H)(t吨,ω,,ω)(t吨,)(){N个( d日,ω)λ(,ω)μ(d日)}

也是一个零鞅,当t吨t吨是固定的,即。

E类[J型(u个)t吨]=J型(t吨)=0,u个>t吨:
E类[(t吨,u个]H(H)(t吨,ω,,ω)N个( d日,ω)t吨]=E类[(t吨,u个]H(H)(t吨,ω,,ω)λ(,ω)μ(d日)t吨].

由于条件强度测量是扩散的,因此该方程右侧的积分对于t吨因此,如果u个→∞, 用于固定t吨ω,F类(t吨,ω,t吨,ω)保持连续并适应ℱt吨。来自Brémaud(1981),定理T5,第I.3章,F类(t吨,ω,t吨,ω)在以下方面是可预测的(t吨,ω).

此外,如果H(H)(t吨,ω,t吨,ω)采取以下形式(,b条] U型(ω)(′,b条′](t吨)V(V)(ω),U型∈ ℱ,V(V)∈ ℱ,然后F类采取形式(,b条] U型(ω)(f)(t吨,ω)其中(f)(t吨,ω)是可预测的,即。F类是可以预测的。

相关结论G公司也可以得到类似的证明。

提议1。假设二阶ℱ-可预测过程H(H)(t吨,ω,t吨,ω)可以类似地分解为H(H)(t吨,ω,t吨,ω),H(H)+(t吨,ω,t吨,ω)和H(H)0(t吨,ω,t吨,ω)如方程式所示(25). 然后

E类[D类H(H)(t吨,ω,t吨,ω)N个( d日t吨,ω)N个( d日t吨,ω)]=E类[D类H(H)(t吨,ω,t吨,ω,)λ(t吨,ω)λ(t吨,ω)μ(d日t吨)μ(d日t吨)]
29

如果右侧积分存在,其中D类是以下项的可测量子集×.

证明。请注意

E类[D类H(H)(t吨,ω,t吨,ω)N个( d日t吨,ω)N个( d日t吨,ω)]=E类[dom公司(D类)E类[0H(H)(t吨,ω,t吨,ω)D类(t吨)(t吨)N个( d日t吨,ω)t吨]N个( d日t吨,ω)]

其中dom(D类)={t吨:∃t吨这样的话(t吨,t吨) ∈D类}和D类(t吨)={t吨:(t吨,t吨) ∈D类}. 从引理1,

F类(t吨,ω,t吨,ω)E类[0H(H)(t吨,ω,,ω)D类(t吨)()N个( d日,ω)t吨]

是可预测的吗,这意味着

F类(t吨,ω)F类(t吨,ω,t吨,ω)=E类[0H(H)(t吨,ω,t吨,ω)D类(t吨)(t吨)N个( d日t吨,ω)t吨]

是可以预测的。

因此,

=E类[dom公司(D类)F类(t吨,ω)N个( d日t吨,ω)]=E类[dom公司(D类)F类(t吨,ω)λ(t吨,ω)μ(d日t吨)].

然而,从引理1来看,方程的右侧(29)是

E类[D类H(H)(t吨,ω,t吨,ω)λ(t吨,ω)λ(t吨,ω)μ(d日t吨)μ(d日t吨)]                                                            =E类[dom公司(D类)E类[D类(t吨)H(H)(t吨,ω,t吨,ω)λ(t吨,ω)μ(d日t吨)t吨]λ(t吨,ω)μ(d日t吨)]
=E类[dom公司(D类)F类(t吨,ω)λ(t吨,ω)μ(d日t吨)].

类似地,方程式(29)如果我们更换H(H)通过H(H)+,即。

E类[D类H(H)+(t吨,ω,t吨,ω)N个( d日t吨,ω)N个( d日t吨,ω)]=E类[D类H(H)+(t吨,ω,t吨,ω,)λ(t吨,ω)λ(t吨,ω)μ(d日t吨)μ(d日t吨)].
30

什么时候?t吨=t吨,H(H)0(t吨,ω)≡H(H)0(t吨,ω,t吨,ω)是一个可预测的过程N个(d)t吨,ω)N个(d)t吨,ω)=N个(d)t吨,ω); 什么时候t吨t吨,H(H)0(t吨,ω,t吨,ω)=0. 因此

E类[D类H(H)0(t吨,ω,t吨,ω)N个( d日t吨,ω)N个( d日t吨,ω)]=E类[诊断(D类)H(H)0(t吨,ω,t吨,ω)N个( d日t吨,ω)]=E类[诊断(D类)H(H)0(t吨,ω,t吨,ω)λ(d日t吨,ω)μ(d日t吨)].
31

其中诊断(D类)={t吨:(t吨,t吨) ∈D类}.

组合方程式(29)–(31),我们有以下定理。

定理1。对于任何二阶可预测过程H(H)(t吨,ω,t吨,ω),给定积分过程N个(ω)具有条件强度λ(t吨,ω)、等式

E类[D类H(H)(t吨,ω,t吨,ω)N个( d日t吨,ω)N个( d日t吨,ω)]=E类[D类H(H)(t吨,ω,t吨,ω)λ(t吨,ω)λ(t吨,ω)μ(d日t吨)μ(d日t吨)]                                                                                                      +E类[诊断(D类)H(H)0(t吨,ω,t吨,ω)λ(d日t吨,ω)μ(d日t吨)],

E类[D类\诊断2(D类)H(H)(t吨,ω,t吨,ω)N个( d日t吨,ω)N个( d日t吨,ω)]=E类[D类H(H)(t吨,ω,t吨,ω)λ(t吨,ω)λ(t吨,ω)μ(d日t吨)μ(d日t吨)]

如果右侧积分存在,则保持不变,其中diag2(D类)={(t吨,t吨) ∈D类}.

附录C:流行型余震序列模型的估计程序(庄.,2002,2004)

  • (a)

    建立初始背景地震活动率;例如,让u个^(x个)=1.

  • (b)
    设置
    λ(t吨,x个)=γ(){νu个^(x个)+:t吨<t吨κ()(t吨t吨)(f)(x个x个)},
    和估计ν以及通过最大化log-likelihood函数的其他模型参数
    日志(L(左))=(t吨,x个,)D类日志{λ(t吨,x个,)}D类λ(t吨,x个,)d日t吨 d日x个 d日,
    哪里D类是指定的感兴趣的时空区域。
  • (c)
    对于每个事件,套
    φ^=νu个^(x个)γ()/λ(t吨,x个,).
    32
  • (d)
    使用加权核估计获得更好的背景速率估计
    u个^(x个)φ^Z轴(x个x个,小时)
    33
    哪里Z轴表示高斯核函数和带宽小时是到n个第页第个最接近的事件如果此距离小于阈值带宽,则为阈值带宽。
  • (e)

    用这个更好的速率替换背景速率并返回步骤(b)。重复上述步骤,直到结果收敛。

本文根据牛津大学出版社标准期刊出版模式的条款出版和发行(https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model)