最优实验设计

在心理学研究中选择下一个要问的问题、设计选举投票策略以及在生物科学中决定合成和测试哪些化合物等任务，本质上都在问同一个问题：我们如何设计一个实验以最大化收集到的信息？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, 和 Isabella Verdinelli. “贝叶斯实验设计：综述.” 统计科学 (1995): 273-304.

[2] Foster, Adam, 等. “变分贝叶斯最优实验设计.” arXiv 预印本 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)[source]

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

Parameters

  • model (function) – Pyro 随机函数，将 design 作为唯一参数。

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

  • observation_labels (list) – 被视为可观测样本点的标签。

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

  • 指南 (函数) – 对应于 模型 的 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 表示是否应尝试使用均值场先验的解析形式。

Returns

EIG估计，可选地包括完整的优化历史

Return type

torch.Tensor

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

自版本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。

Parameters

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

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

  • observation_labels (list) – model中存在的一部分样本站点。这些站点被视为未来的观测值，而其他站点则被视为需要推断后验的潜在变量。

  • target_labels (list) – 用于测量后验熵的样本点的一个子集。

  • vi_parameters (dict) – 变分推断参数，应包括： optim: pyro.Optim 的一个实例，guide: 与 model 兼容的引导函数， num_steps: 要进行的变分推断步骤数， 以及 loss: 用于变分推断的损失函数

  • 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 表示是否应尝试使用均值场先验的解析形式。

Returns

EIG估计，可选地包括完整的优化历史

Return type

torch.Tensor

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

嵌套蒙特卡洛估计期望信息增益（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)\)。 存在随机效应时的估计为

\[\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 时使用后一种形式。

Parameters

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

  • 设计 (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时使用。指示目标变量和干扰变量的先验分布是否独立。在这种情况下，不需要在干扰变量的条件下对目标进行采样。

Returns

EIG估计，可选地包括完整的优化历史

Return type

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)[source]

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 上的损失函数。

Parameters

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

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

  • observation_labels (list) – model中存在的一部分样本站点。这些站点被视为未来的观测值，而其他站点则被视为需要推断后验的潜在变量。

  • target_labels (list) – 用于测量后验熵的样本点的一个子集。

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

  • num_steps (int) – 优化步骤的数量。

  • T (函数 或 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`。

Returns

EIG估计，可选地包括完整的优化历史

Return type

torch.Tensor 或 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, 等. “变分贝叶斯最优实验设计.” arXiv 预印本 arXiv:1903.05480 (2019).

Parameters

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

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

  • observation_labels (list) – model中存在的一部分样本站点。这些站点被视为未来的观测值，而其他站点则被视为需要推断后验的潜在变量。

  • target_labels (list) – 用于测量后验熵的样本点的一个子集。

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

  • num_steps (int) – 优化步骤的数量。

  • 指南 (函数) – 用于（隐式）后验估计的指南系列。 指南的参数被优化以最大化后验目标。

  • 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 表示是否应尝试使用均值场先验的解析形式。

Returns

EIG估计，可选地包括完整的优化历史

Return type

torch.Tensor 或 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\) 上的任意分布。\(q\) 的变分族在 guide 中指定。

警告

此方法在存在随机效应时不能估计正确的数量。

[1] Foster, Adam, 等. “变分贝叶斯最优实验设计.” arXiv 预印本 arXiv:1903.05480 (2019).

Parameters

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

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

  • observation_labels (list) – model中存在的一部分样本站点。这些站点被视为未来的观测值，而其他站点则被视为需要推断后验的潜在变量。

  • target_labels (list) – 用于测量后验熵的样本点的一个子集。

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

  • num_steps (int) – 优化步骤的数量。

  • 指南 (函数) – 用于边际估计的指南系列。 guide 的参数经过优化，以最大化对数似然目标。

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

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

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

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

Returns

EIG估计，可选地包括完整的优化历史

Return type

torch.Tensor 或 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)[source]

使用如[1]所述的基于比率估计的无似然推断方法（LFIRE）来估计EIG。 LFIRE 针对 \(\theta\) 的多个样本分别运行。

[1] Kleinegesse, Steven, 和 Michael Gutmann. “高效贝叶斯实验设计用于隐式模型.” arXiv 预印本 arXiv:1810.09912 (2018).

Parameters

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

  • 设计 (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估计。

Returns

EIG估计，可选地包括完整的优化历史

Return type

torch.Tensor 或 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)[source]

使用变分嵌套蒙特卡洛（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, 等. “变分贝叶斯最优实验设计.” arXiv 预印本 arXiv:1903.05480 (2019).

Parameters

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

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

  • observation_labels (list) – model中存在的一部分样本站点。这些站点被视为未来的观测值，而其他站点则被视为需要推断后验的潜在变量。

  • target_labels (list) – 用于测量后验熵的样本点的一个子集。

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

  • num_steps (int) – 优化步骤的数量。

  • 指南 (函数) – 用于后验估计的指南系列。 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`。

Returns

EIG估计，可选地包括完整的优化历史

Return type

torch.Tensor 或 tuple

广义线性混合模型

警告

该模块最终将被弃用，转而支持 brmp

pyro.contrib.oed.glmm 模块提供了广义线性混合模型（GLMM）的模型和指南。它还包括正态-逆伽马分布族。

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

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()。