流行病学

警告

pyro.contrib.epidemiology 中的代码正在开发中。此代码不保证维护向后兼容性。

pyro.contrib.epidemiology 为一类随机离散时间离散计数隔室模型提供建模语言。此模块实现了黑盒推断（包括随机变分推断和哈密顿蒙特卡洛）、潜在变量预测和未来轨迹预测。

用法示例请参阅以下教程

基本隔室模型

class CompartmentalModel(compartments, duration, population, *, approximate=())[source]

Bases: abc.ABC

离散时间离散值随机隔室模型的抽象基类。

派生类必须实现 initialize() 和 transition() 方法。派生类可选实现 global_model()、compute_flows() 和 heuristic()。

用法示例

# First implement a concrete derived class.
class MyModel(CompartmentalModel):
    def __init__(self, ...): ...
    def global_model(self): ...
    def initialize(self, params): ...
    def transition(self, params, state, t): ...

# Run inference to fit the model to data.
model = MyModel(...)
model.fit_svi(num_samples=100)  # or .fit_mcmc(...)
R0 = model.samples["R0"]  # An example parameter.
print("R0 = {:0.3g} ± {:0.3g}".format(R0.mean(), R0.std()))

# Predict latent variables.
samples = model.predict()

# Forecast forward.
samples = model.predict(forecast=30)

# You can assess future interventions (applied after ``duration``) by
# storing them as attributes that are read by your derived methods.
model.my_intervention = False
samples1 = model.predict(forecast=30)
model.my_intervention = True
samples2 = model.predict(forecast=30)
effect = samples2["my_result"].mean() - samples1["my_result"].mean()
print("average effect = {:0.3g}".format(effect))

一个示例工作流程是，在寻找良好的模型结构和先验时使用更廉价的近似推断，一旦模型合理，则转而使用更准确但更昂贵的推断。

  1. 首先使用 .fit_svi(guide_rank=0, num_steps=2000) 进行廉价推断，同时寻找一个好的模型。

  2. 此外，通过 .fit_svi(guide_rank=None, num_steps=5000) 转移到低秩多元正态指南，推断长程相关性。

  3. 可选地，通过转移到更昂贵（但仍通过矩匹配进行近似）的 .fit_mcmc(num_quant_bins=1, num_samples=10000, num_chains=2) 来推断非高斯后验。

  4. 可选地，通过转移到更昂贵的基于枚举的算法 .fit_mcmc(num_quant_bins=4, num_samples=10000, num_chains=2) 来改进小计数附近的拟合（推荐使用 GPU）。

变量

samples (dict) – 后验样本字典。

参数

  • compartments (list) – 隔室名称字符串列表。

  • duration (int) – 此模型中的离散时间步长数量。

  • population (int 或 torch.Tensor) – 单区域模型的总人口或区域模型中各区域人口的张量。

  • approximate (tuple) – 应在 transition() 中为其提供逐点近似的隔室名称，例如，如果指定 approximate=("I")，则 state["I_approx"] 将是 state["I"] 的连续值非枚举点估计。近似有助于降低计算成本。近似是连续值，支持范围为 (-0.5, population + 0.5)。

property time_plate

时间维度上的 pyro.plate。

property region_plate

一个 pyro.plate 或一个微不足道的 ExitStack，具体取决于此模型 .is_regional。

property full_mass

全局随机变量名称的单个元组列表。

property series

每个时间步采样到的采样站点名称的冻结集。

global_model()[源代码]

采样并返回任何全局参数。

返回

任意参数对象（例如 None 或元组）。

abstract initialize(params)[源代码]

返回每个隔室的初始计数。

参数

params – 由 global_model() 返回的全局参数。

返回

将隔室名称映射到初始值的字典。

返回类型

dict

abstract transition(params, state, t)[source]

动力学的正向生成过程。

这将输入当前 state 并在原地随机更新该状态。

请注意，此方法在多种不同解释下被调用，包括批处理和向量化解释。在 generate() 期间，调用此方法生成单个样本。在 heuristic() 期间，调用此方法生成用于 SMC 的一批样本。在 fit_mcmc() 期间，此方法以向量化形式（对时间进行向量化）和顺序形式（针对单个时间步长）调用；这两种形式都对离散潜在变量进行枚举。在 predict() 期间，调用此方法以基于时间间隔 [0:duration] 的后验样本预测一批样本。

参数

  • params – 由 global_model() 返回的全局参数。

  • state (dict) – 将隔室名称映射到当前张量值的字典。此字典应在原地更新。

  • t (int 或 slice) – 时间索引。在推断期间，t 可以是切片（用于向量化推断）或整数时间索引。在预测期间，t 将是整数时间索引。

finalize(params, prev, curr)[源代码]

依赖于整个时间序列的似然的可选方法。

此方法仅应用于耦合跨时间状态的不可分解似然。可分解似然应添加到 transition() 方法中，从而启用其在 heuristic() 初始化中的使用。由于此方法仅在最后一个时间步之后调用，因此未在 heuristic() 初始化中使用。

警告

此方法当前不支持潜在变量。

参数

  • params – 由 global_model() 返回的全局参数。

  • prev (dict) –

  • curr (dict) – 将隔室名称映射到整个时间序列张量的字典。这两个参数偏移 1 步，从而可以轻松计算流量时间序列。对于量化推断，此方法使用近似点估计，因此用户必须在 __init__() 中请求任何需要的时间序列，例如，如果似然依赖于 I 和 E 时间序列，则通过调用 super().__init__(..., approximate=("I", "E")) 来实现。

compute_flows(prev, curr, t)[源代码]

计算隔室之间的流量，给定时间步长 t 前后的隔室人口。

默认实现假定顺序流终止于一个名为“R”的隐式隔室。例如，如果

compartment_names = ("S", "E", "I")

默认实现在时间步长 t = 9 计算

flows["S2E_9"] = prev["S"] - curr["S"]
flows["E2I_9"] = prev["E"] - curr["E"] + flows["S2E_9"]
flows["I2R_9"] = prev["I"] - curr["I"] + flows["E2I_9"]

对于更复杂的流程（非顺序的、分支的、循环的、复制的等），用户可以覆盖此方法。

参数

  • state (dict) – 将隔室名称映射到当前张量值的字典。此字典应在原地更新。

  • t (int 或 slice) – 时间索引。在推断期间，t 可以是切片（用于向量化推断）或整数时间索引。在预测期间，t 将是整数时间索引。

返回

将流名称映射到张量值的字典。

返回类型

dict

generate(fixed={})[源代码]

从先验生成数据。

参数字典 fixed

用于条件的参数字典。这些必须是顶层无父节点，即没有上游随机依赖项。

返回

将采样站点名称映射到采样值的字典。

返回类型

dict

fit_svi(*, num_samples=100, num_steps=2000, num_particles=32, learning_rate=0.1, learning_rate_decay=0.01, betas=(0.8, 0.99), haar=True, init_scale=0.01, guide_rank=0, jit=False, log_every=200, **options)[source]

运行随机变分推断生成后验样本。

此方法运行 SVI，并在完成后设置 .samples 属性。

此近似推断方法有助于快速迭代概率模型。

参数

  • num_samples (int) – 从训练好的指南中要抽取的后验样本数。默认为 100。

  • num_steps (int) – SVI 步骤数。

  • num_particles (int) – 每步 SVI 粒子数。

  • learning_rate (int) – ClippedAdam 优化器的学习率。

  • learning_rate_decay (int) – ClippedAdam 优化器的学习率衰减。注意这是整个计划的衰减，而不是每步衰减。

  • betas (tuple) – ClippedAdam 优化器的动量参数。

  • haar (bool) – 是否使用 Haar 小波重新参数化器。

  • guide_rank (int) – 自动正态指南的秩。如果为零（默认），则使用 AutoNormal 指南。如果是正整数或 None，则使用 AutoLowRankMultivariateNormal 指南。如果字符串为“full”，则使用 AutoMultivariateNormal 指南。后两者需要更多 num_steps 来拟合。

  • init_scale (float) – AutoLowRankMultivariateNormal 指南的初始比例。

  • jit (bool) – 是否使用 jit 编译的 ELBO。

  • log_every (int) – 记录 svi 损失的频率。

  • heuristic_num_particles (int) – 作为 num_particles 传递给 heuristic()。默认为 1024。

返回

SVI 损失时间序列（有助于诊断收敛性）。

返回类型

list

fit_mcmc(**options)[源代码]

运行 NUTS 推断生成后验样本。

此方法使用 NUTS 核运行 MCMC，并在完成后设置 .samples 属性。

当 num_quant_bins > 1 时，此方法使用渐近精确的基于枚举的模型；当 num_quant_bins == 1 时，使用更廉价的矩匹配近似模型。

参数

  • **options – 传递给 MCMC 的选项。其余选项被提取并具有特殊含义。

  • num_samples (int) – 通过 MCMC 抽取的后验样本数。默认为 100。

  • max_tree_depth (int) – （默认为 5）。NUTS 核的最大树深度。

  • full_mass – NUTS 核的质量矩阵规格。默认为全局随机变量的完整质量矩阵。

  • arrowhead_mass (bool) – 是否将 full_mass 视为箭形矩阵的头部，而不仅仅是块。默认为 False。

  • num_quant_bins (int) – 如果大于 1，则通过对这么多量化 bin 进行局部枚举使用渐近精确推断。如果等于 1，则使用连续值松弛近似推断。请注意，计算成本随 num_quant_bins 呈指数级增长。松弛推断默认为 1。

  • haar (bool) – 是否使用 Haar 小波重新参数化器。默认为 True。

  • haar_full_mass (int) – 完整质量矩阵中包含的低频 Haar 分量数量。如果 haar=False，则忽略此参数。默认为 10。

  • heuristic_num_particles (int) – 作为 num_particles 传递给 heuristic()。默认为 1024。

返回

用于诊断的 MCMC 对象，例如 MCMC.summary()。

返回类型

MCMC

predict(forecast=0)[源代码]

预测潜在变量并可选地向前预测。

此方法只能在 fit_mcmc() 之后运行，并抽取与传递给 fit_mcmc() 的 num_samples 相同的数量。

参数

forecast (int) – 向前预测的时间步长数。

返回

将采样站点名称（或隔室名称）映射到张量的字典，其中第一个维度对应于样本批处理。

返回类型

dict

heuristic(num_particles=1024, ess_threshold=0.5, retries=10)[source]

找到所有潜在变量的初始可行猜测，与观测数据一致。这需要是因为并非所有假设都可行，并且 HMC 需要从一个可行的解开始才能进行。

默认实现尝试使用 SMCFilter 并从先验中获取提议来找到一个可行的状态。但是，在 SMC 表现不佳的情况下（例如在高维模型中），此方法可以被覆盖。

参数

  • num_particles (int) – 用于 SMC 的粒子数。

  • ess_threshold (float) – SMC 的有效样本大小阈值。

返回

将采样站点名称映射到张量值的字典。

返回类型

dict

模型示例

简单 SIR

class SimpleSIRModel(population, recovery_time, data)[source]

易感-感染-恢复模型。

要自定义此模型，我们建议 Forking 并编辑此类。

这是一个随机离散时间离散状态模型，包含三个隔室：“S”表示易感者，“I”表示感染者，“R”表示恢复者（恢复者是隐式的：R = population - S - I），转换路径为 S -> I -> R。

参数

  • population (int) – 总人口 population = S + I + R。

  • recovery_time (float) – 平均恢复时间（在状态 I 中的持续时间）。必须大于 1。

  • data (iterable) – 新观测感染的时间序列。每个时间步长在 0 到 S -> I 转换数量之间遵循二项分布。这允许假阴性但不允许假阳性。

简单 SEIR

class SimpleSEIRModel(population, incubation_time, recovery_time, data)[source]

易感-暴露-感染-恢复模型。

要自定义此模型，我们建议 Forking 并编辑此类。

这是一个随机离散时间离散状态模型，包含四个隔室：“S”表示易感者，“E”表示暴露者，“I”表示感染者，“R”表示恢复者（恢复者是隐式的：R = population - S - E - I），转换路径为 S -> E -> I -> R。

参数

  • population (int) – 总人口 population = S + E + I + R。

  • incubation_time (float) – 平均潜伏时间（在状态 E 中的持续时间）。必须大于 1。

  • recovery_time (float) – 平均恢复时间（在状态 I 中的持续时间）。必须大于 1。

  • data (iterable) – 新观测感染的时间序列。每个时间步长在 0 到 S -> E 转换数量之间遵循二项分布。这允许假阴性但不允许假阳性。

简单 SEIRD

class SimpleSEIRDModel(population, incubation_time, recovery_time, mortality_rate, data)[source]

易感-暴露-感染-恢复-死亡模型。

要自定义此模型，我们建议 Forking 并编辑此类。

这是一个随机离散时间离散状态模型，包含四个隔室：“S”表示易感者，“E”表示暴露者，“I”表示感染者，“D”表示死亡个体，“R”表示恢复者（恢复者是隐式的：R = population - S - E - I - D），转换路径为 S -> E -> I -> R 和 I -> D。

由于转换不是简单的线性顺序，此模型实现了自定义的 compute_flows() 方法。

参数

  • population (int) – 总人口 population = S + E + I + R + D。

  • incubation_time (float) – 平均潜伏时间（在状态 E 中的持续时间）。必须大于 1。

  • recovery_time (float) – 平均恢复时间（在状态 I 中的持续时间）。必须大于 1。

  • mortality_rate (float) – 导致死亡的感染比例。必须在开区间 (0, 1) 内。

  • data (iterable) – 新观测感染的时间序列。每个时间步长在 0 到 S -> E 转换数量之间遵循二项分布。这允许假阴性但不允许假阳性。

过度分散的 SIR

class OverdispersedSIRModel(population, recovery_time, data)[source]

通过过度分散分布泛化 SimpleSIRModel。

要自定义此模型，我们建议 Forking 并编辑此类。

此模型添加了一个全局过度分散参数，控制转换和观测分布的过度分散。有关分布详情，请参阅 binomial_dist() 和 beta_binomial_dist()。有关结合过度分散分布的先前工作，请参阅 [1,2,3,4]。

参考文献

[1] D. Champredon, M. Li, B. Bolker. J. Dushoff (2018)

“预测埃博拉合成疫情的两种方法” https://www.sciencedirect.com/science/article/pii/S1755436517300233

[2] Carrie Reed et al. (2015)

“从美国基于人群的监测数据估算流感疾病负担” https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4349859/

[3] A. Leonard, D. Weissman, B. Greenbaum, E. Ghedin, K. Koelle (2017)

“从病原体深度测序数据估算传播瓶颈大小，及其在人甲型流感病毒上的应用” https://jvi.asm.org/content/jvi/91/14/e00171-17.full.pdf

[4] A. Miller, N. Foti, J. Lewnard, N. Jewell, C. Guestrin, E. Fox (2020)

“出行趋势提供了 SARS-CoV-2 传播变化的领先指标” https://www.medrxiv.org/content/medrxiv/early/2020/05/11/2020.05.07.20094441.full.pdf

参数

  • population (int) – 总人口 population = S + I + R。

  • recovery_time (float) – 平均恢复时间（在状态 I 中的持续时间）。必须大于 1。

  • data (iterable) – 新观测感染的时间序列。每个时间步长在 0 到 S -> I 转换数量之间遵循二项分布。这允许假阴性但不允许假阳性。

过度分散的 SEIR

class OverdispersedSEIRModel(population, incubation_time, recovery_time, data)[source]

通过过度分散分布泛化 SimpleSEIRModel。

要自定义此模型，我们建议 Forking 并编辑此类。

此模型添加了一个全局过度分散参数，控制转换和观测分布的过度分散。有关分布详情，请参阅 binomial_dist() 和 beta_binomial_dist()。有关结合过度分散分布的先前工作，请参阅 [1,2,3,4]。

参考文献

[1] D. Champredon, M. Li, B. Bolker. J. Dushoff (2018)

“预测埃博拉合成疫情的两种方法” https://www.sciencedirect.com/science/article/pii/S1755436517300233

[2] Carrie Reed et al. (2015)

“从美国基于人群的监测数据估算流感疾病负担” https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4349859/

[3] A. Leonard, D. Weissman, B. Greenbaum, E. Ghedin, K. Koelle (2017)

“从病原体深度测序数据估算传播瓶颈大小，及其在人甲型流感病毒上的应用” https://jvi.asm.org/content/jvi/91/14/e00171-17.full.pdf

[4] A. Miller, N. Foti, J. Lewnard, N. Jewell, C. Guestrin, E. Fox (2020)

“出行趋势提供了 SARS-CoV-2 传播变化的领先指标” https://www.medrxiv.org/content/medrxiv/early/2020/05/11/2020.05.07.20094441.full.pdf

参数

  • population (int) – 总人口 population = S + E + I + R。

  • incubation_time (float) – 平均潜伏时间（在状态 E 中的持续时间）。必须大于 1。

  • recovery_time (float) – 平均恢复时间（在状态 I 中的持续时间）。必须大于 1。

  • data (iterable) – 新观测感染的时间序列。每个时间步长在 0 到 S -> E 转换数量之间遵循二项分布。这允许假阴性但不允许假阳性。

超级传播 SIR

class SuperspreadingSIRModel(population, recovery_time, data)[source]

通过添加超级传播效应泛化 SimpleSIRModel。

要自定义此模型，我们建议 Forking 并编辑此类。

此模型解释超级传播（个体再生数过度分散），假设每个感染个体感染 BetaBinomial 数量的易感个体，其中 BetaBinomial 分布作为过度分散的二项分布，调整了作为过度分散泊松分布的更标准的负二项分布 [1,2] 到有限人口设定。为保留马尔可夫结构，我们遵循 [2] 的方法，并假设单个个体引起的所有感染都发生在该个体进行 I -> R 转换的单一时间步长。也就是说，SimpleSIRModel 假设感染个体在每个感染时间步长（平均 τ 步）感染 Binomial(S,R/tau) 数量的易感个体，而此模型假设他们感染 BetaBinomial(k,…,S) 数量的易感个体，但仅在恢复前的最后一个时间步长。

参考文献

[1] J. O. Lloyd-Smith, S. J. Schreiber, P. E. Kopp, W. M. Getz (2005)

“超级传播和个体变异对疾病出现的影响” https://www.nature.com/articles/nature04153.pdf

[2] Lucy M. Li, Nicholas C. Grassly, Christophe Fraser (2017)

“使用病原体系统发育和发病时间序列量化传播异质性” https://academic.oup.com/mbe/article/34/11/2982/3952784

参数

  • population (int) – 总人口 population = S + I + R。

  • recovery_time (float) – 平均恢复时间（在状态 I 中的持续时间）。必须大于 1。

  • data (iterable) – 新观测感染的时间序列。每个时间步长在 0 到 S -> I 转换数量之间遵循二项分布。这允许假阴性但不允许假阳性。

超级传播 SEIR

class SuperspreadingSEIRModel(population, incubation_time, recovery_time, data, *, leaf_times=None, coal_times=None)[source]

通过添加超级传播效应泛化 SimpleSEIRModel。

要自定义此模型，我们建议 Forking 并编辑此类。

此模型解释超级传播（个体再生数过度分散），假设每个感染个体感染 BetaBinomial 数量的易感个体，其中 BetaBinomial 分布作为过度分散的二项分布，调整了作为过度分散泊松分布的更标准的负二项分布 [1,2] 到有限人口设定。为保留马尔可夫结构，我们遵循 [2] 的方法，并假设单个个体引起的所有感染都发生在该个体进行 I -> R 转换的单一时间步长。也就是说，SimpleSEIRModel 假设感染个体在每个感染时间步长（平均 τ 步）感染 Binomial(S,R/tau) 数量的易感个体，而此模型假设他们感染 BetaBinomial(k,…,S) 数量的易感个体，但仅在恢复前的最后一个时间步长。

此模型还为合并时间形式的观测系统发育数据添加了一个可选似然。这些数据以一对 (leaf_times, coal_times) 的形式提供，分别表示基因组测序和谱系合并的时间。我们使用 CoalescentRateLikelihood 结合这些数据，其中基础合并率从 S 和 I 人口计算得出。此似然在时间上是独立的，并保留了推断所需的马尔可夫性质。

参考文献

[1] J. O. Lloyd-Smith, S. J. Schreiber, P. E. Kopp, W. M. Getz (2005)

“超级传播和个体变异对疾病出现的影响” https://www.nature.com/articles/nature04153.pdf

[2] Lucy M. Li, Nicholas C. Grassly, Christophe Fraser (2017)

“使用病原体系统发育和发病时间序列量化传播异质性” https://academic.oup.com/mbe/article/34/11/2982/3952784

参数

  • population (int) – 总人口 population = S + E + I + R。

  • incubation_time (float) – 平均潜伏时间（在状态 E 中的持续时间）。必须大于 1。

  • recovery_time (float) – 平均恢复时间（在状态 I 中的持续时间）。必须大于 1。

  • data (iterable) – 新观测感染的时间序列。每个时间步长在 0 到 S -> E 转换数量之间遵循二项分布。这允许假阴性但不允许假阳性。

异构 SIR

class HeterogeneousSIRModel(population, recovery_time, data)[source]

通过允许 Rt 和 rho 随时间变化泛化 SimpleSIRModel。

要自定义此模型，我们建议 Forking 并编辑此类。

在此模型中，响应率 rho 是分段常数，在三个分段上具有未知值。再生数 Rt 是常数 R0 与在对数空间通过布朗运动漂移的因子 beta 的乘积。 rho 和 Rt 均可作为时间序列使用。

参数

  • population (int) – 总人口 population = S + I + R。

  • recovery_time (float) – 平均恢复时间（在状态 I 中的持续时间）。必须大于 1。

  • data (iterable) – 新观测感染的时间序列。每个时间步长在 0 到 S -> I 转换数量之间遵循二项分布。这允许假阴性但不允许假阳性。

稀疏 SIR

class SparseSIRModel(population, recovery_time, data, mask)[source]

泛化 SimpleSIRModel 以允许稀疏观测到的感染。

要自定义此模型，我们建议 Forking 并编辑此类。

此模型允许在不均匀的时间间隔内观测到累积感染。为了保留马尔可夫结构（因此实现可处理的推断），此模型添加了一个辅助隔室 O，表示每个时间点完全观测到的累积观测数。在观测时间（当 mask[t] == True 时），O 必须与提供的数据完全匹配；在观测时间之间，O 随机插补提供的数据。

此模型演示了如何实现自定义的 compute_flows() 方法。在此模型中需要自定义方法，因为 S 隔室的居民可以转换到 I 和 O 隔室，从而允许重复。

参数

  • population (int) – 总人口 population = S + I + R。

  • recovery_time (float) – 平均恢复时间（在状态 I 中的持续时间）。必须大于 1。

  • data (iterable) – 累积观测感染的时间序列。每当 mask[t] == True 时，data[t] 对应于一个观测值；否则 data[t] 可以是任意值，例如 NAN。

  • mask (iterable) – 布尔时间序列，表示是否在每个时间步长进行观测。应满足 len(mask) == len(data)。

未知起始时间 SIR

class UnknownStartSIRModel(population, recovery_time, pre_obs_window, data)[source]

通过允许未知首次感染日期泛化 SimpleSIRModel。

要自定义此模型，我们建议 Forking 并编辑此类。

此模型演示了

  1. 如何结合来自外部来源的自发感染；

  2. 如何通过支持在 transition() 中进行预测来结合随时间变化的分段 rho。

  3. 如何覆盖 predict() 方法来计算额外统计量。

参数

  • population (int) – 总人口 population = S + I + R。

  • recovery_time (float) – 平均恢复时间（在状态 I 中的持续时间）。必须大于 1。

  • pre_obs_window (int) – 数据开始前首次感染可能发生的时间步长数。必须是正数。

  • data (iterable) – 新观测感染的时间序列。每个时间步长在 0 到 S -> I 转换数量之间遵循二项分布。这允许假阴性但不允许假阳性。

区域性 SIR

class RegionalSIRModel(population, coupling, recovery_time, data)[source]

泛化 SimpleSIRModel 以同时对具有弱区域间耦合的多个区域建模。

要自定义此模型，我们建议 Forking 并编辑此类。

区域通过一个 coupling 矩阵耦合，其条目在 [0,1] 范围内。全一矩阵等价于单个区域。单位矩阵等价于一组独立区域。此矩阵不需要对称，但对称矩阵在物理上可能更合理。每个时间步长的新增感染数 S2I 遵循二项分布，均值为

E[S2I] = S (1 - (1 - R0 / (population @ coupling)) ** (I @ coupling))
       ≈ R0 S (I @ coupling) / (population @ coupling)  # for small I

因此，在几乎完全易感的人口中，单个感染个体平均感染大约 R0 个新个体，与 coupling 无关。

此模型演示了

  1. 如何创建具有 population 向量的区域模型。

  2. 如何使用 self.region_plate 对齐次参数（此处为 R0）和具有层次结构的异构参数（此处为 rho）进行建模。

  3. 如何在 transition() 中使用 state["I_approx"] 近似耦合区域。

参数

  • population (torch.Tensor) – 各区域人口张量，定义 population = S + I + R。

  • coupling (torch.Tensor) – 成对耦合矩阵。条目应在 [0,1] 范围内。

  • recovery_time (float) – 平均恢复时间（在状态 I 中的持续时间）。必须大于 1。

  • data (iterable) – 时间 x 区域 大小的新观测感染张量。每个时间步长是二项分布的向量，范围在 0 到 S -> I 转换数量之间。这允许假阴性但不允许假阳性。

异构区域性 SIR

class HeterogeneousRegionalSIRModel(population, coupling, recovery_time, data)[source]

通过允许 Rt 和 rho 随时间变化泛化 RegionalSIRModel。

要自定义此模型，我们建议 Forking 并编辑此类。

在此模型中，响应率 rho 随时间和区域变化，而再生数 Rt 随时间变化但在区域间共享。两个参数都根据变换后的布朗运动漂移，并学习漂移率。

此模型演示了如何对除隔室变量以外的分层潜在时间序列进行建模。

参数

  • population (torch.Tensor) – 各区域人口张量，定义 population = S + I + R。

  • coupling (torch.Tensor) – 成对耦合矩阵。条目应在 [0,1] 范围内。

  • recovery_time (float) – 平均恢复时间（在状态 I 中的持续时间）。必须大于 1。

  • data (iterable) – 时间 x 区域 大小的新观测感染张量。每个时间步长是二项分布的向量，范围在 0 到 S -> I 转换数量之间。这允许假阴性但不允许假阳性。

分布

set_approx_sample_thresh(thresh)[源代码]

实验性上下文管理器/装饰器，用于临时设置 Binomial.approx_sample_thresh 的全局默认值，从而降低从 Binomial、BetaBinomial、ExtendedBinomial、ExtendedBetaBinomial 和 infection_dist() 返回的分布中采样的计算复杂度。

这对于从非常大的 total_count 进行采样很有用。

此方法由 CompartmentalModel 内部使用。

参数

thresh (int 或 float.) – 新临时阈值。

set_approx_log_prob_tol(tol)[源代码]

实验性上下文管理器/装饰器，用于临时设置 Binomial.approx_log_prob_tol 和 BetaBinomial.approx_log_prob_tol 的全局默认值，从而降低对 Binomial 和 BetaBinomial 分布进行评分的计算复杂度。

此方法由 CompartmentalModel 内部使用。

参数

tol (int 或 float.) – 新临时阈值。

binomial_dist(total_count, probs, *, overdispersion=0.0)[source]

返回一个 Beta-Binomial 分布，它是二项分布的过度分散版本，根据参数 overdispersion，通常设置在 0.1 到 0.5 的范围内。

这对于 (1) 拟合相对于二项分布过度分散的真实数据，以及 (2) 放宽大种群模型以改善推断很有用。特别是，overdispersion 参数设定了随机模型中相对不确定性的下界，使得人口增加导致一个具有有限随机性的无标度动力系统，与在大种群极限下收敛到确定性常微分方程的基于二项分布的 SDE 相反。

此参数化满足以下属性

  1. 方差随 overdispersion 单调增加。

  2. overdispersion = 0 导致二项分布。

  3. overdispersion 设定了相对不确定性 std_dev / (total_count * p * q) 的下界，其中 probs = p = 1 - q，并随着 total_count → ∞ 作为相对不确定性的渐近线。这与相对不确定性趋于零的二项分布形成对比。

  4. 如果 X ~ binomial_dist(n, p, overdispersion=σ)，则在大种群极限 n → ∞ 中，缩放随机变量 X / n 在分布上收敛到 LogitNormal(log(p/(1-p)), σ)。

为实现这些属性，我们设置 p = probs，q = 1 - p，以及

concentration = 1 / (p * q * overdispersion**2) - 1
参数

  • total_count (int 或 torch.Tensor) – 伯努利试验次数。

  • probs (float 或 torch.Tensor) – 事件概率。

  • overdispersion (float 或 torch.tensor) – 过度分散量，在半开区间 [0,2) 内。默认为零。

beta_binomial_dist(concentration1, concentration0, total_count, *, overdispersion=0.0)[source]

返回一个 Beta-Binomial 分布，它是通常 Beta-Binomial 分布的过度分散版本，根据额外的参数 overdispersion，通常设置在 0.1 到 0.5 的范围内。

参数

  • concentration1 (float 或 torch.Tensor) – Beta 分布的第一个集中度参数 (alpha)。

  • concentration0 (float 或 torch.Tensor) – Beta 分布的第二个集中度参数 (beta)。

  • total_count (float 或 torch.Tensor) – 伯努利试验次数。

  • overdispersion (float 或 torch.tensor) – 过度分散量，在半开区间 [0,2) 内。默认为零。

infection_dist(*, individual_rate, num_infectious, num_susceptible=inf, population=inf, concentration=inf, overdispersion=0.0)[source]

创建关于离散时间步长新增感染数量的 Distribution。

此函数根据 population 和 concentration 是否有限返回泊松、负二项、二项或 Beta-Binomial 分布。在 Pyro 模型中，人口通常是有限的。在极限 population → ∞ 和 num_susceptible/population → 1 中，二项分布收敛到泊松分布，Beta-Binomial 分布收敛到负二项分布。在极限 concentration → ∞ 中，负二项分布收敛到泊松分布，Beta-Binomial 分布收敛到二项分布。

过度分散的分布（当 concentration < ∞ 时返回的负二项分布和 Beta-二项分布）可用于模拟超级传播者个体 [1,2]。在小群体中以及在截断或审查成本高昂的概率编程系统中，有限支持的分布（二项分布和负二项分布）非常有用 [3]。

参考文献

[1] J. O. Lloyd-Smith, S. J. Schreiber, P. E. Kopp, W. M. Getz (2005)

“超级传播和个体变异对疾病出现的影响” https://www.nature.com/articles/nature04153.pdf

[2] Lucy M. Li, Nicholas C. Grassly, Christophe Fraser (2017)

“使用病原体系统发育和发病时间序列量化传播异质性” https://academic.oup.com/mbe/article/34/11/2982/3952784

[3] Lawrence Murray et al. (2018)

“概率程序的延迟采样和自动 Rao-Blackwellization” https://arxiv.org/pdf/1708.07787.pdf

参数

  • individual_rate – 在大群体极限下，每个感染个体每时间步的平均感染数，等于 R0 / tau，其中 R0 是基本再生数，tau 是平均感染持续时间。

  • num_infectious – 在此时间步的感染个体数，有时是 I，有时是 E+I。

  • num_susceptible – 在此时间步的易感个体数 S。默认为无限大的人口。

  • population – 人口中的总个体数。默认为无限大的人口。

  • concentration – 超级传播者过度分散模型 [1,2] 中的集中度或分散参数 k。默认为最小方差 concentration = ∞。

  • overdispersion (float 或 torch.tensor) – 过度分散量，在半开区间 [0,2) 内。默认为零。

class CoalescentRateLikelihood(leaf_times, coal_times, duration, *, validate_args=None)[source]

基础： object

实验性 这不是一个 Distribution，而是 CoalescentTimesWithRate 的转置版本，使得 rate_grid 的元素独立，从而与 plate 和 poutine.markov 兼容。对于非批处理输入，以下所有似然是等效的

# Version 1.
pyro.sample("coalescent",
            CoalescentTimesWithRate(leaf_times, rate_grid),
            obs=coal_times)

# Version 2. using pyro.plate
likelihood = CoalescentRateLikelihood(leaf_times, coal_times, len(rate_grid))
with pyro.plate("time", len(rate_grid)):
    pyro.factor("coalescent", likelihood(rate_grid))

# Version 3. using pyro.markov
likelihood = CoalescentRateLikelihood(leaf_times, coal_times, len(rate_grid))
for t in pyro.markov(range(len(rate_grid))):
    pyro.factor("coalescent_{}".format(t), likelihood(rate_grid[t], t))

第三种版本适用于例如 SMCFilter，其中 rate_grid 可能会按顺序计算。

参数

  • leaf_times (torch.Tensor) – 采样事件的时间张量，即系统发生树中的叶节点。这些可以是任意实数，顺序任意且可以重复。

  • coal_times (torch.Tensor) – 合并时间张量。这些时间表示尾部维度上大小为 leaf_times.size(-1) - 1 的集合，并且应该沿该维度排序。

  • duration (int) – 速率网格的大小，rate_grid.size(-1)。

__call__(rate_grid, t=slice(None, None, None))[source]

计算 [1] 方程 7-9 在一个或所有时间点的似然。

参考文献

[1] A. Popinga, T. Vaughan, T. Statler, A.J. Drummond (2014)

“使用 Bayesian 合并推断推断流行病学动态：确定性模型和随机模型的优点” https://arxiv.org/pdf/1407.1792.pdf

参数

  • rate_grid (torch.Tensor) – 基本合并速率（成对合并速率）的张量。例如，在简单的 SIR 模型中，这可能是 beta S / I。最右边的维度是时间，该张量表示在时间上分段常量的速率（的批次）。

  • time (int 或 slice) – 可选的时间索引，输入通过该索引进行切片，例如 rate_grid[..., t]。对于顺序模型，这可以是整数；对于向量化模型，可以是 slice(None)。

返回

p(coal_times | leaf_times, rate_grid) 的似然，或者与单个时间步相对应的似然的一部分。

返回类型

torch.Tensor

bio_phylo_to_times(tree, *, get_time=None)[source]

从系统发生树中提取合并汇总统计信息，适合与 CoalescentRateLikelihood 一起使用。

参数

  • tree (Bio.Phylo.BaseTree.Clade) – 系统发生树。

  • get_time (callable) – 可选函数，用于提取每个子 Clade 的时间点。如果缺失，时间将通过累积的 .branch_length 计算。

返回

一对 Tensor s (leaf_times, coal_times)，其中 leaf_times 是采样事件的时间（系统发生树中的叶节点），coal_times 是合并的时间（系统发生二叉树中的叶节点）。

返回类型

元组