最优实验设计

选择心理学研究中要问的下一个问题、设计选举民意调查策略以及决定在生物科学中合成和测试哪些化合物等任务，从根本上都在问同一个问题：我们如何设计实验以最大化收集的信息？Pyro 旨在支持自动化最优实验设计：指定模型和指南就足以获得许多不同类型实验场景的最优设计。请查看我们的实验设计教程，这些教程使用 Pyro 设计自适应心理学研究，利用过去的数据选择下一个问题，并设计选举民意调查策略，旨在对最终选举获胜者给出最强的预测。

贝叶斯最优实验设计 (BOED) 是一种处理实验设计问题的强大方法，也是 Pyro 采用的框架。在 BOED 框架中，我们首先建立一个贝叶斯模型，该模型具有似然 \(p(y|\theta,d)\) 和目标潜在变量的先验 \(p(\theta)\)。在 Pyro 中，任何全贝叶斯模型都可以用于 BOED 框架。对应于实验结果的采样位点是观测位点，对应于感兴趣的潜在变量的采样位点是目标位点。设计 \(d\) 是模型的参数，而不是随机变量。

在 BOED 框架中，我们选择最大化运行实验后从目标 \(\theta\) 中获得的预期信息增益 (EIG) 的设计

\(\text{EIG}(d) = \mathbf{E}_{p(y|d)} [H[p(\theta)] − H[p(\theta|y, d)]]\) ,

其中 \(H[·]\) 代表熵，\(p(\theta|y, d) \propto p(\theta)p(y|\theta, d)\) 是使用设计 \(d\) 并观测到 \(y\) 后得到的后验。换句话说，最优设计是在预期未来观测值的情况下，对目标潜在变量的后验熵减少最多的设计。如果预测模型是正确的，这将形成一种从信息论角度看是（一步）最优的设计策略。更多细节请参见 [1, 2]。

pyro.contrib.oed 模块提供了为 Pyro 模型创建最优实验设计的工具。特别是，它提供了预期信息增益 (EIG) 的估计器。

要估计特定设计的 EIG，我们首先建立我们的 Pyro 模型。例如

def model(design):

    # This line allows batching of designs, treating all batch dimensions as independent
    with pyro.plate_stack("plate_stack", design.shape):

        # We use a Normal prior for theta
        theta = pyro.sample("theta", dist.Normal(torch.tensor(0.0), torch.tensor(1.0)))

        # We use a simple logistic regression model for the likelihood
        logit_p = theta - design
        y = pyro.sample("y", dist.Bernoulli(logits=logit_p))

        return y

然后我们选择合适的 EIG 估计器，例如

eig = nmc_eig(model, design, observation_labels=["y"], target_labels=["theta"], N=2500, M=50)

可以在设计网格中估计 EIG

designs = torch.stack([design1, design2], dim=0)

从多个选项中找到最佳设计。

[1] Chaloner, Kathryn, and Isabella Verdinelli. “Bayesian experimental design: A review.” Statistical Science (1995): 273-304.

[2] Foster, Adam, et al. “Variational Bayesian Optimal Experimental Design.” arXiv preprint arXiv:1903.05480 (2019).

预期信息增益

laplace_eig(model, design, observation_labels, target_labels, guide, loss, optim, num_steps, final_num_samples, y_dist=None, eig=True, **prior_entropy_kwargs)[源代码]

通过对后验进行重复的拉普拉斯近似来估计预期信息增益 (EIG)。

参数

  • model (function) – Pyro 随机函数，只接受 design 作为参数。

  • design (torch.Tensor) – 可能设计的张量。

  • observation_labels (list) – 被视为可观测变量的采样位点标签。

  • target_labels (list) – 被视为感兴趣的潜在变量的采样位点标签，即我们希望获取信息的位点。

  • guide (function) – 对应于 model 的 Pyro 随机函数。

  • loss – Pyro 损失函数，例如 pyro.infer.Trace_ELBO().differentiable_loss。

  • optim – 损失优化器

  • num_steps (int) – 每个采样伪观测值执行的梯度步数。

  • final_num_samples (int) – 采样的 y 样本（伪观测值）数量。

  • y_dist – 采样 y 的分布 - 如果为 None，则使用贝叶斯边际分布。

  • eig (bool) – 是否计算 EIG 或平均后验熵 (APE)。EIG 由 EIG = 先验熵 - APE 给出。如果为 True，先验熵将根据 model 进行解析估计或蒙特卡洛估计。如果为 False，则返回 APE。

  • prior_entropy_kwargs (dict) – 估计先验熵的参数：num_prior_samples 表示 MC 估计先验熵的样本数量，以及 mean_field 表示是否尝试均场先验的解析形式。

返回

EIG 估计值，可选包含完整的优化历史记录

返回类型

torch.Tensor

vi_eig(model, design, observation_labels, target_labels, vi_parameters, is_parameters, y_dist=None, eig=True, **prior_entropy_kwargs)[源代码]

自版本 0.4.1 起已弃用: 请改用 posterior_eig。

使用变分推断 (VI) 估计预期信息增益 (EIG)。

APE 定义为

\(APE(d)=E_{Y\sim p(y|\theta, d)}[H(p(\theta|Y, d))]\)

其中 \(H[p(x)]\) 是微分熵。APE 与预期信息增益 (EIG) 的关系由公式给出

\(EIG(d)=H[p(\theta)]-APE(d)\)

特别是，最小化 APE 等效于最大化 EIG。

参数

  • model (function) – 接受 design 作为唯一参数的 pyro 模型。

  • design (torch.Tensor) – 设计的张量表示

  • observation_labels (list) – model 中存在的采样位点的子集。这些位点被视为未来观测，其他位点被视为要推断后验的潜在变量。

  • target_labels (list) – 要衡量后验熵的采样位点的子集。

  • vi_parameters (dict) – 变分推断参数，应包括：optim：pyro.Optim 的实例，guide：与 model 兼容的指南函数，num_steps：执行的 VI 步数，以及 loss：用于 VI 的损失函数

  • is_parameters (dict) – 响应变量 \(Y\) 的边际分布的 重要性采样参数。可能包括 num_samples：从边际分布中抽取的样本数量。

  • y_dist (pyro.distributions.Distribution) – (可选) 响应变量 \(Y\) 假定的分布

  • eig (bool) – 是否计算 EIG 或平均后验熵 (APE)。EIG 由 EIG = 先验熵 - APE 给出。如果为 True，先验熵将根据 model 进行解析估计或蒙特卡洛估计。如果为 False，则返回 APE。

  • prior_entropy_kwargs (dict) – 估计先验熵的参数：num_prior_samples 表示 MC 估计先验熵的样本数量，以及 mean_field 表示是否尝试均场先验的解析形式。

返回

EIG 估计值，可选包含完整的优化历史记录

返回类型

torch.Tensor

nmc_eig(model, design, observation_labels, target_labels=None, N=100, M=10, M_prime=None, independent_priors=False)[源代码]

嵌套蒙特卡洛估计预期信息增益 (EIG)。当不存在随机效应时，估计值为

\[\frac{1}{N}\sum_{n=1}^N \log p(y_n | \theta_n, d) - \frac{1}{N}\sum_{n=1}^N \log \left(\frac{1}{M}\sum_{m=1}^M p(y_n | \theta_m, d)\right)\]

其中 \(\theta_n, y_n \sim p(\theta, y | d)\) 且 \(\theta_m \sim p(\theta)\)。当 M_prime != None 时，使用后一种形式。

\[\frac{1}{N}\sum_{n=1}^N \log \left(\frac{1}{M'}\sum_{m=1}^{M'} p(y_n | \theta_n, \widetilde{\theta}_{nm}, d)\right)- \frac{1}{N}\sum_{n=1}^N \log \left(\frac{1}{M}\sum_{m=1}^{M} p(y_n | \theta_m, \widetilde{\theta}_{m}, d)\right)\]

其中 \(\widetilde{\theta}\) 是随机效应，且 \(\widetilde{\theta}_{nm} \sim p(\widetilde{\theta}|\theta=\theta_n)\) 和 \(\theta_m,\widetilde{\theta}_m \sim p(\theta,\widetilde{\theta})\)。当 M_prime != None 时使用后一种形式。

参数

  • model (function) – 接受 design 作为唯一参数的 pyro 模型。

  • design (torch.Tensor) – 设计的张量表示

  • observation_labels (list) – model 中存在的采样位点的子集。这些位点被视为未来观测，其他位点被视为要推断后验的潜在变量。

  • target_labels (list) – 要衡量后验熵的采样位点的子集。

  • N (int) – 外部期望样本数量。

  • M (int) – 用于 p(y|d) 的内部期望样本数量。

  • M_prime (int) – 如果需要，用于 p(y | theta, d) 的样本数量。

  • independent_priors (bool) – 仅当 M_prime 不为 None 时使用。指示目标变量和无关变量的先验分布是否独立。在这种情况下，不需要以无关变量为条件采样目标变量。

返回

EIG 估计值，可选包含完整的优化历史记录

返回类型

torch.Tensor

donsker_varadhan_eig(model, design, observation_labels, target_labels, num_samples, num_steps, T, optim, return_history=False, final_design=None, final_num_samples=None)[源代码]

Donsker-Varadhan 估计预期信息增益 (EIG)。

EIG 的 Donsker-Varadhan 表示为

\[\sup_T E_{p(y, \theta | d)}[T(y, \theta)] - \log E_{p(y|d)p(\theta)}[\exp(T(\bar{y}, \bar{\theta}))]\]

其中 \(T\) 是任何（可测）函数。

此方法在预定义函数类 T 上优化损失函数。

参数

  • model (function) – 接受 design 作为唯一参数的 pyro 模型。

  • design (torch.Tensor) – 设计的张量表示

  • observation_labels (list) – model 中存在的采样位点的子集。这些位点被视为未来观测，其他位点被视为要推断后验的潜在变量。

  • target_labels (list) – 要衡量后验熵的采样位点的子集。

  • num_samples (int) – 每次迭代的样本数量。

  • num_steps (int) – 优化步数。

  • T (function or torch.nn.Module) – Donsker-Varadhan 损失函数中使用的可优化函数 T。

  • optim (pyro.optim.Optim) – 使用的优化器。

  • return_history (bool) – 如果为 True，还返回一个张量，给出优化过程中每一步的损失函数。

  • final_design (torch.Tensor) – 最终评估的设计张量。如果为 None，则使用 design。

  • final_num_samples (int) – 最终评估时使用的样本数量。如果为 None，则使用 num_samples。

返回

EIG 估计值，可选包含完整的优化历史记录

返回类型

torch.Tensor or tuple

posterior_eig(model, design, observation_labels, target_labels, num_samples, num_steps, guide, optim, return_history=False, final_design=None, final_num_samples=None, eig=True, prior_entropy_kwargs={}, *args, **kwargs)[源代码]

根据平均后验熵 (APE) 计算的预期信息增益 (EIG) 的后验估计，使用公式 \(EIG(d) = H[p(\theta)] - APE(d)\)。详见 [1]。

APE 的后验表示为

\(\sup_{q}\ E_{p(y, \theta | d)}[\log q(\theta | y, d)]\)

其中 \(q\) 是 \(\theta\) 上的任何分布。

此方法在给定的代表 \(q\) 的 guide 族上优化损失。

[1] Foster, Adam, et al. “Variational Bayesian Optimal Experimental Design.” arXiv preprint arXiv:1903.05480 (2019).

参数

  • model (function) – 接受 design 作为唯一参数的 pyro 模型。

  • design (torch.Tensor) – 设计的张量表示

  • observation_labels (list) – model 中存在的采样位点的子集。这些位点被视为未来观测，其他位点被视为要推断后验的潜在变量。

  • target_labels (list) – 要衡量后验熵的采样位点的子集。

  • num_samples (int) – 每次迭代的样本数量。

  • num_steps (int) – 优化步数。

  • guide (function) – 用于（隐式）后验估计的指南族。优化 guide 的参数以最大化后验目标。

  • optim (pyro.optim.Optim) – 使用的优化器。

  • return_history (bool) – 如果为 True，还返回一个张量，给出优化过程中每一步的损失函数。

  • final_design (torch.Tensor) – 最终评估的设计张量。如果为 None，则使用 design。

  • final_num_samples (int) – 最终评估时使用的样本数量。如果为 None，则使用 num_samples。

  • eig (bool) – 是否计算 EIG 或平均后验熵 (APE)。EIG 由 EIG = 先验熵 - APE 给出。如果为 True，先验熵将根据 model 进行解析估计或蒙特卡洛估计。如果为 False，则返回 APE。

  • prior_entropy_kwargs (dict) – 估计先验熵的参数：num_prior_samples 表示 MC 估计先验熵的样本数量，以及 mean_field 表示是否尝试均场先验的解析形式。

返回

EIG 估计值，可选包含完整的优化历史记录

返回类型

torch.Tensor or tuple

marginal_eig(model, design, observation_labels, target_labels, num_samples, num_steps, guide, optim, return_history=False, final_design=None, final_num_samples=None)[源代码]

通过估计边际熵 \(p(y|d)\) 来估计 EIG。详见 [1]。

EIG 的边际表示为

\(\inf_{q}\ E_{p(y, \theta | d)}\left[\log \frac{p(y | \theta, d)}{q(y | d)} \right]\)

其中 \(q\) 是 \(y\) 上的任何分布。guide 中指定了一个用于 \(q\) 的变分族。

警告

在存在随机效应的情况下，此方法不会估计正确的量。

[1] Foster, Adam, et al. “Variational Bayesian Optimal Experimental Design.” arXiv preprint arXiv:1903.05480 (2019).

参数

  • model (function) – 接受 design 作为唯一参数的 pyro 模型。

  • design (torch.Tensor) – 设计的张量表示

  • observation_labels (list) – model 中存在的采样位点的子集。这些位点被视为未来观测，其他位点被视为要推断后验的潜在变量。

  • target_labels (list) – 要衡量后验熵的采样位点的子集。

  • num_samples (int) – 每次迭代的样本数量。

  • num_steps (int) – 优化步数。

  • guide (function) – 用于边际估计的指南族。优化 guide 的参数以最大化对数似然目标。

  • optim (pyro.optim.Optim) – 使用的优化器。

  • return_history (bool) – 如果为 True，还返回一个张量，给出优化过程中每一步的损失函数。

  • final_design (torch.Tensor) – 最终评估的设计张量。如果为 None，则使用 design。

  • final_num_samples (int) – 最终评估时使用的样本数量。如果为 None，则使用 num_samples。

返回

EIG 估计值，可选包含完整的优化历史记录

返回类型

torch.Tensor or tuple

lfire_eig(model, design, observation_labels, target_labels, num_y_samples, num_theta_samples, num_steps, classifier, optim, return_history=False, final_design=None, final_num_samples=None)[源代码]

使用似然无关推断通过比率估计 (LFIRE) 的方法估计 EIG，如 [1] 中所述。LFIRE 对 \(\theta\) 的几个样本分别运行。

[1] Kleinegesse, Steven, and Michael Gutmann. “Efficient Bayesian Experimental Design for Implicit Models.” arXiv preprint arXiv:1810.09912 (2018).

参数

  • model (function) – 接受 design 作为唯一参数的 pyro 模型。

  • design (torch.Tensor) – 设计的张量表示

  • observation_labels (list) – model 中存在的采样位点的子集。这些位点被视为未来观测，其他位点被视为要推断后验的潜在变量。

  • target_labels (list) – 要衡量后验熵的采样位点的子集。

  • num_y_samples (int) – 对于每个 \(\theta\)，在 \(y\) 中抽取的样本数量。

  • num_steps (int) – 优化步数。

  • classifier (function) – 用于区分 \(p(y|d)\) 下的 \(y\) 样本和在某个 \(\theta\) 下 \(p(y|\theta,d)\) 样本的 Pytorch 或 Pyro 分类器。

  • optim (pyro.optim.Optim) – 使用的优化器。

  • return_history (bool) – 如果为 True，还返回一个张量，给出优化过程中每一步的损失函数。

  • final_design (torch.Tensor) – 最终评估的设计张量。如果为 None，则使用 design。

  • final_num_samples (int) – 最终评估时使用的样本数量。如果为 None，则使用 num_samples。

Param

int num_theta_samples: 在 \(\theta\) 中抽取的初始样本数量。通过 LFIRE 估计每个样本的似然比。

返回

EIG 估计值，可选包含完整的优化历史记录

返回类型

torch.Tensor or tuple

vnmc_eig(model, design, observation_labels, target_labels, num_samples, num_steps, guide, optim, return_history=False, final_design=None, final_num_samples=None)[源代码]

使用变分嵌套蒙特卡洛 (VNMC) 估计 EIG。VNMC 估计值 [1] 为

\[\frac{1}{N}\sum_{n=1}^N \left[ \log p(y_n | \theta_n, d) - \log \left(\frac{1}{M}\sum_{m=1}^M \frac{p(\theta_{mn})p(y_n | \theta_{mn}, d)} {q(\theta_{mn} | y_n)} \right) \right]\]

其中 \(q(\theta | y)\) 是学习到的变分后验近似，且 \(\theta_n, y_n \sim p(\theta, y | d)\)，\(\theta_{mn} \sim q(\theta|y=y_n)\)。

当 \(N \to \infty\) 时，这是 EIG 的一个上限。我们通过随机梯度下降来最小化此上限。

警告

此方法不能在存在随机效应的情况下使用。

[1] Foster, Adam, et al. “Variational Bayesian Optimal Experimental Design.” arXiv preprint arXiv:1903.05480 (2019).

参数

  • model (function) – 接受 design 作为唯一参数的 pyro 模型。

  • design (torch.Tensor) – 设计的张量表示

  • observation_labels (list) – model 中存在的采样位点的子集。这些位点被视为未来观测，其他位点被视为要推断后验的潜在变量。

  • target_labels (list) – 要衡量后验熵的采样位点的子集。

  • num_samples (tuple) – 每次迭代的 (\(N, M\)) 样本数量。

  • num_steps (int) – 优化步数。

  • guide (function) – 用于后验估计的指南族。优化 guide 的参数以最小化 VNMC 上限。

  • optim (pyro.optim.Optim) – 使用的优化器。

  • return_history (bool) – 如果为 True，还返回一个张量，给出优化过程中每一步的损失函数。

  • final_design (torch.Tensor) – 最终评估的设计张量。如果为 None，则使用 design。

  • final_num_samples (tuple) – 最终评估时使用的 (\(N, M\)) 样本数量。如果为 None，则使用 num_samples。

返回

EIG 估计值，可选包含完整的优化历史记录

返回类型

torch.Tensor or tuple

广义线性混合模型

警告

此模块最终将被弃用，转而使用 brmp

pyro.contrib.oed.glmm 模块提供了用于广义线性混合模型 (GLMM) 的模型和指南。它还包括 Normal-inverse-gamma 族。

要创建经典的贝叶斯线性模型，请使用

from pyro.contrib.oed.glmm import known_covariance_linear_model

# Note: coef is a p-vector, observation_sd is a scalar
# Here, p=1 (one feature)
model = known_covariance_linear_model(coef_mean=torch.tensor([0.]),
                                      coef_sd=torch.tensor([10.]),
                                      observation_sd=torch.tensor(2.))

# An n x p design tensor
# Here, n=2 (two observations)
design = torch.tensor(torch.tensor([[1.], [-1.]]))

model(design)

可以引入非线性链接函数，例如

from pyro.contrib.oed.glmm import logistic_regression_model

# No observation_sd is needed for logistic models
model = logistic_regression_model(coef_mean=torch.tensor([0.]),
                                  coef_sd=torch.tensor([10.]))

随机效应可以作为常规贝叶斯回归系数合并。对于具有共享协方差矩阵的随机效应，请参见 pyro.contrib.oed.glmm.lmer_model()。