杂项操作¶
pyro.ops
模块实现了大部分独立于 Pyro 其他部分的张量工具。
HMC的实用工具¶
- class DualAveraging(prox_center=0, t0=10, kappa=0.75, gamma=0.05)[source]¶
基础类:
object
双重平均法是一种解决凸优化问题的方案。它属于一类使用次梯度来更新模型参数(在原始空间中)的次梯度方法。在某些条件下,该方案中生成的参数的平均值保证会收敛到一个最优值。然而,传统次梯度方法的一个反直觉方面是“新的次梯度以递减的权重进入模型”(参见\([1]\))。双重平均法通过使用相等的权重更新次梯度(位于对偶空间中)的参数来解决这一现象,因此我们称之为“双重平均”。
该类实现了一种适用于马尔可夫链蒙特卡罗(MCMC)算法的双重平均方案。更准确地说,我们将用MCMC轨迹中计算的一些统计量来替代次梯度。此外,引入一些自由参数如
t0
和kappa
是有帮助的,并且仍然保证方案的收敛性。参考文献
[1] 凸问题的原始对偶次梯度方法, 尤里·涅斯捷罗夫
[2] 无掉头采样器:在哈密顿蒙特卡洛中自适应设置路径长度, Matthew D. Hoffman, Andrew Gelman
- Parameters
- velocity_verlet(z, r, potential_fn, kinetic_grad, step_size, num_steps=1, z_grads=None)[source]¶
使用速度Verlet算法的二阶辛积分器。
- Parameters
potential_fn (callable) – 返回给定z的势能的函数 对于每个采样点。函数的负梯度相对于
z
决定了相应站点的变化率 动量r
。kinetic_grad (callable) – 一个计算动能相对于动量变量的梯度的函数。
step_size (float) – 每次时间步迭代的步长。
num_steps (int) – 进行积分的离散时间步数。
z_grads (torch.Tensor) – 当前
z
处的势能梯度(可选)。
- Return tuple (z_next, r_next, z_grads, potential_energy)
下一个位置和动量, 以及势能及其相对于
z_next
的梯度。
- potential_grad(potential_fn, z)[源代码]¶
potential_fn 关于参数 z 的梯度。
- Parameters
potential_fn – 一个可调用的python函数,它接收一个参数字典并返回势能。
z (dict) – 以站点名称为键的参数值字典。
- Returns
元组 (z_grads, potential_energy),其中 z_grads 是一个字典, 其键与 z 相同,包含梯度,而 potential_energy 是一个 torch 标量。
- register_exception_handler(name: str, handler: Callable[[Exception], bool], warn_on_overwrite: bool = True) None [source]¶
注册一个异常处理器,用于在评估势函数时处理(主要是数值)错误。
- Parameters
name – 处理程序的名称(必须唯一)。
handler – 一个可调用的映射,将异常映射到布尔值。在任何handler中评估为true的异常将在势能的计算中被处理。
warn_on_overwrite – 如果为True,当覆盖已注册的处理程序时会发出警告。
牛顿优化器¶
- newton_step(loss, x, trust_radius=None)[source]¶
执行牛顿更新步骤以最小化一批变量的损失,可选择约束到信任区域 [1]。
这尤其有用,因为牛顿迭代的最终解相对于输入是可微的,即使除了最终的
x
之外的所有内容都被分离,由于该方法的二次收敛性[2]。loss
必须作为x
的函数是二次可微的。如果loss
是2+d
次可微的,那么该函数的返回值是d
次可微的。当
loss
被解释为负对数概率密度时,可以使用此函数的返回值mode,cov
来构建拉普拉斯近似MultivariateNormal(mode,cov)
。警告
在优化循环中使用此函数时,请注意分离其结果。如果在优化过程中忘记分离此函数的结果,那么反向传播将贯穿整个迭代过程,更糟糕的是,每一步都会计算两个额外的导数。
在循环中的示例用法:
x = torch.zeros(1000, 2) # arbitrary initial value for step in range(100): x = x.detach() # block gradients through previous steps x.requires_grad = True # ensure loss is differentiable wrt x loss = my_loss_function(x) x = newton_step(loss, x, trust_radius=1.0) # the final x is still differentiable
- [1] Yuan, Ya-xiang. Iciam. Vol. 99. 2000.
“优化中信任区域算法的回顾。” ftp://ftp.cc.ac.cn/pub/yyx/papers/p995.pdf
- [2] Christianson, Bruce. Optimization Methods and Software 3.4 (1994)
“反向累积和吸引固定点。” http://uhra.herts.ac.uk/bitstream/handle/2299/4338/903839.pdf
- Parameters
loss (torch.Tensor) – 一个关于
x
的标量函数,需要被最小化。x (torch.Tensor) – 一个形状为
(N, D)
的因变量,其中N
是批量大小,D
是一个较小的数字。trust_radius (float) – 一个可选的信任区域 trust_radius。该函数的更新值
mode
将在输入x
的trust_radius
范围内。
- Returns
一对
(mode, cov)
,其中mode
是一个与原始值x
形状相同的更新张量,而cov
是一个估计的协方差 DxD 矩阵,其形状为cov.shape == x.shape[:-1] + (D,D)
。- Return type
- newton_step_1d(loss, x, trust_radius=None)[source]¶
执行牛顿更新步骤以最小化一批一维变量的损失,可选择正则化以约束到信任区域。
详情请参见
newton_step()
。- Parameters
loss (torch.Tensor) – 一个关于
x
的标量函数,需要被最小化。x (torch.Tensor) – 一个最右侧大小为1的因变量。
trust_radius (float) – 一个可选的信任区域 trust_radius。此函数的更新值
mode
将在输入x
的trust_radius
范围内。
- Returns
一对
(mode, cov)
其中mode
是一个与原始值x
形状相同的更新张量,而cov
是一个估计的协方差 1x1 矩阵,其形状为cov.shape == x.shape[:-1] + (1,1)
。- Return type
- newton_step_2d(loss, x, trust_radius=None)[source]¶
执行牛顿更新步骤以最小化一批二维变量的损失,可选择正则化以约束到信任区域。
详情请参见
newton_step()
。- Parameters
loss (torch.Tensor) – 一个关于
x
的标量函数,需要被最小化。x (torch.Tensor) – 一个最右侧大小为2的因变量。
trust_radius (float) – 一个可选的信任区域 trust_radius。此函数的更新值
mode
将在输入x
的trust_radius
范围内。
- Returns
一对
(mode, cov)
其中mode
是一个与原始值x
形状相同的更新张量,而cov
是一个估计的协方差 2x2 矩阵,其形状为cov.shape == x.shape[:-1] + (2,2)
。- Return type
- newton_step_3d(loss, x, trust_radius=None)[源代码]¶
执行牛顿更新步骤以最小化一批三维变量的损失,可选择正则化以约束到信任区域。
详情请参见
newton_step()
。- Parameters
loss (torch.Tensor) – 一个关于
x
的标量函数,需要被最小化。x (torch.Tensor) – 一个最右侧大小为2的因变量。
trust_radius (float) – 一个可选的信任区域 trust_radius。此函数的更新值
mode
将在输入x
的trust_radius
范围内。
- Returns
一对
(mode, cov)
其中mode
是一个与原始值x
形状相同的更新张量,而cov
是一个估计的协方差 3x3 矩阵,其形状为cov.shape == x.shape[:-1] + (3,3)
。- Return type
特殊函数¶
- safe_log(x)[source]¶
类似于
torch.log()
,但通过将 log(0) 的梯度限制在最多1 / finfo.eps
来避免无限梯度。
- log_beta(x, y, tol=0.0)[source]¶
计算对数Beta函数。
当
tol >= 0.02
时,这使用了一个移位的斯特林近似来计算对数贝塔函数。该近似方法调整了斯特林近似来计算对数伽玛函数:lgamma(z) ≈ (z - 1/2) * log(z) - z + log(2 * pi) / 2
近似对数Beta函数:
log_beta(x, y) ≈ ((x-1/2) * log(x) + (y-1/2) * log(y) - (x+y-1/2) * log(x+y) + log(2*pi)/2)
近似值通过使用递归迭代地移动对数伽马近似值,从而在接近零时进一步提高准确性:
lgamma(x) = lgamma(x + 1) - log(x)
如果此递归应用了
n
次,则绝对误差被限制在error < 0.082 / n < tol
,因此我们根据用户提供的tol
选择n
。- Parameters
x (torch.Tensor) – 一个正张量。
y (torch.Tensor) – 一个正张量。
tol (float) – 最大绝对误差的界限。默认为0.1。对于非常小的
tol
,此函数会直接调用log_beta()
。
- Return type
- log_binomial(n, k, tol=0.0)[source]¶
计算对数二项式系数。
当
tol >= 0.02
时,这使用了一个移位的斯特林近似来计算对数 Beta 函数,通过log_beta()
。- Parameters
n (torch.Tensor) – 一个非负整数张量。
k (torch.Tensor) – 一个整数张量,范围在
[0, n]
。
- Return type
- log_I1(orders: int, value: torch.Tensor, terms=250)[source]¶
计算第一类修正贝塞尔函数的前n个对数 .. math
\log(I_v(z)) = v*\log(z/2) + \log(\sum_{k=0}^\inf \exp\left[2*k*\log(z/2) - \sum_kk^k log(kk) - \lgamma(v + k + 1)\right])
- Parameters
orders – 对数修正贝塞尔函数的阶数。
value – 用于计算修正贝塞尔函数的值
terms – 求和的截断
- Returns
0 到 orders 修改的贝塞尔函数
- get_quad_rule(num_quad, prototype_tensor)[source]¶
获取具有指定数量积分点的高斯-埃尔米特积分规则的积分点及相应的对数权重。
示例用法:
quad_points, log_weights = get_quad_rule(32, prototype_tensor) # transform to N(0, 4.0) Normal distribution quad_points *= 4.0 # compute variance integral in log-space using logsumexp and exponentiate variance = torch.logsumexp(quad_points.pow(2.0).log() + log_weights, axis=0).exp() assert (variance - 16.0).abs().item() < 1.0e-6
- Parameters
num_quad (int) – 积分点的数量。
prototype_tensor (torch.Tensor) – 用于确定返回张量的dtype和device。
- Returns
元组形式为torch.Tensor`s的`(quad_points, log_weights)
- sparse_multinomial_likelihood(total_count, nonzero_logits, nonzero_value)[source]¶
以下是等价的:
# Version 1. dense log_prob = Multinomial(logits=logits).log_prob(value).sum() # Version 2. sparse nnz = value.nonzero(as_tuple=True) log_prob = sparse_multinomial_likelihood( value.sum(-1), (logits - logits.logsumexp(-1))[nnz], value[nnz], )
张量工具¶
- as_complex(x)[source]¶
类似于
torch.view_as_complex()
,但在步幅不是两的倍数的情况下会复制数据。
- block_diag_embed(mat)[source]¶
接受一个形状为 (…, B, M, N) 的张量,并返回一个形状为 (…, B x M, B x N) 的块对角张量。
- Parameters
mat (torch.Tensor) – 一个具有3个或更多维度的输入张量
- Returns torch.Tensor
一个具有维度 m.dim() - 1 的块对角张量
- block_diagonal(mat, block_size)[源代码]¶
接受一个形状为 (…, B x M, B x N) 的块对角张量,并返回一个形状为 (…, B, M, N) 的张量。
- Parameters
mat (torch.Tensor) – 一个具有2个或更多维度的输入张量
block_size (int) – 块的数量 B。
- Returns torch.Tensor
一个维度为 mat.dim() + 1 的张量
- periodic_repeat(tensor, size, dim)[source]¶
重复一个
period
大小的张量直到达到给定的size
。例如:>>> x = torch.tensor([[1, 2, 3], [4, 5, 6]]) >>> periodic_repeat(x, 4, 0) tensor([[1, 2, 3], [4, 5, 6], [1, 2, 3], [4, 5, 6]]) >>> periodic_repeat(x, 4, 1) tensor([[1, 2, 3, 1], [4, 5, 6, 4]])
这对于计算时间序列模型中的静态季节性非常有用。
- Parameters
张量 (torch.Tensor) – 一个差异的张量。
size (int) – 沿着维度
dim
的结果的期望大小。dim (int) – 沿着该张量维度进行重复。
- periodic_cumsum(tensor, period, dim)[source]¶
沿给定维度计算周期性累积和。例如,如果 dim=0:
for t in range(period): assert result[t] == tensor[t] for t in range(period, len(tensor)): assert result[t] == tensor[t] + result[t - period]
这对于计算时间序列模型中的漂移季节性非常有用。
- Parameters
tensor (torch.Tensor) – 一个差异的张量。
period (int) – 重复的周期。
dim (int) – 沿此张量维度进行累积。
- periodic_features(duration, max_period=None, min_period=None, **options)[source]¶
从
max_period
到min_period
创建周期性(sin,cos)特征。这在时间序列模型中非常有用,其中长期不均匀的季节性可以通过回归来处理。当仅指定
max_period
时,这会在所有长度尺度上生成周期性特征。当还指定min_period
时,这会在较大的长度尺度上生成周期性特征,但会忽略高频特征。这在将长期季节性的回归与其他技术(如periodic_repeat()
和periodic_cumsum()
)结合用于短时间尺度时非常有用。例如,要将年度季节性回归结合到一周的尺度,可以设置max_period=365.25
和min_period=7
。
- next_fast_len(size)[source]¶
返回下一个最大的数字
n >= size
,其质因数全部为2、3或5。这些大小对于快速傅里叶变换是高效的。等同于scipy.fftpack.next_fast_len()
。- Parameters
size (int) – 一个正数。
- Returns
可能更大的数字。
- Rtype int
- convolve(signal, kernel, mode='full')[source]¶
使用FFT计算信号与核的一维卷积。 两个参数应具有相同的最右侧维度,但除此之外可以是任意可广播的。
- Parameters
signal (torch.Tensor) – 一个要进行卷积的信号。
kernel (torch.Tensor) – 一个卷积核。
mode (str) – 其中之一:'full', 'valid', 'same'。
- Returns
一个具有广播形状的张量。假设
m = signal.size(-1)
且n = kernel.size(-1)
,结果的最右边大小将是:m + n - 1
如果模式是‘full’;max(m, n) - min(m, n) + 1
如果模式是‘valid’;或者max(m, n)
如果模式是‘same’。- Rtype torch.Tensor
- repeated_matmul(M, n)[source]¶
将一批矩阵 M 作为输入,并返回执行 n 次矩阵乘法 \(M\), \(M^2\), …, \(M^n\) 的堆叠结果。 并行成本在 n 中是对数的。
- Parameters
M (torch.Tensor) – 一批形状为 (…, N, N) 的方形张量。
n (int) – 最大乘积的阶数 \(M^n\)
- Returns torch.Tensor
一批形状为 (n, …, N, N) 的方形张量
- dct(x, dim=- 1)[source]¶
类型II的离散余弦变换,缩放为正交归一化。
这是
idct_ii()
的逆操作,并且等同于scipy.fftpack.dct()
使用norm="ortho"
。- Parameters
x (Tensor) – 输入信号。
dim (int) – 计算DCT的维度。
- Return type
张量
- idct(x, dim=- 1)[source]¶
类型II的逆离散余弦变换,缩放为正交归一化。
这是
dct_ii()
的逆函数,并且等同于scipy.fftpack.idct()
使用norm="ortho"
。- Parameters
x (Tensor) – 输入信号。
dim (int) – 计算DCT的维度。
- Return type
张量
- haar_transform(x)[source]¶
离散哈尔变换。
沿最后一个维度执行Haar变换。 这是
inverse_haar_transform()
的逆操作。- Parameters
x (Tensor) – 输入信号。
- Return type
张量
- inverse_haar_transform(x)[source]¶
沿最后一个维度执行逆哈尔变换。 这是
haar_transform()
的逆操作。- Parameters
x (Tensor) – 输入信号。
- Return type
张量
- safe_normalize(x, *, p=2)[source]¶
安全地将向量投影到球体上,相对于
p
-范数。这通过将零映射到向量[1, 0, 0, ..., 0]
来避免在零处的奇点。- Parameters
x (torch.Tensor) – 一个向量
p (float) – 范数指数,默认为2,即欧几里得范数。
- Returns
归一化版本
x / ||x||_p
。- Return type
张量
张量索引¶
- index(tensor, args)[source]¶
使用嵌套元组进行索引。
另请参阅便捷包装器
Index
。这对于编写与多种解释兼容的索引代码非常有用,例如标量评估、向量化评估或重塑。
例如,假设
x
是一个参数,且x.dim() == 2
,我们希望 推广表达式x[..., t]
,其中t
可以是以下任意一种:一个标量
t=1
如x[..., 1]
;一个切片
t=slice(None)
等同于x[..., :]
;或者一个重塑操作
t=(Ellipsis, None)
等同于x.unsqueeze(-1)
。
虽然简单的索引适用于前两个例子,但第三个例子会导致一个嵌套的元组
(Ellipsis, (Ellipsis, None))
。这个辅助函数会将该嵌套元组展平并合并连续的Ellipsis
。- Parameters
张量 (torch.Tensor) – 一个需要被索引的张量。
args (tuple) – 一个索引,作为
__getitem__
的参数。
- Returns
对
tensor[args]
的扁平化解释。- Return type
- class Index(tensor)[source]¶
基础类:
object
围绕
index()
的便捷包装器。以下是等价的:
Index(x)[..., i, j, :] index(x, (Ellipsis, i, j, slice(None)))
- Parameters
tensor (torch.Tensor) – 一个要被索引的张量。
- Returns
一个具有特殊
__getitem__()
方法的对象。
- vindex(tensor, args)[source]¶
具有广播语义的矢量化高级索引。
另请参阅便捷包装器
Vindex
。这对于编写与批处理和枚举兼容的索引代码非常有用,特别是对于使用离散随机变量选择混合组件的情况。
例如,假设
x
是一个参数,且x.dim() == 3
,我们希望将表达式x[i, :, j]
从整数i,j
推广到具有批次维度和枚举维度(但没有事件维度)的张量i,j
。然后我们可以使用Vindex
来编写推广版本。xij = Vindex(x)[i, :, j] batch_shape = broadcast_shape(i.shape, j.shape) event_shape = (x.size(1),) assert xij.shape == batch_shape + event_shape
处理当
x
可能还包含批次维度的情况时(例如,如果x
在使用向量化粒子时在板状上下文中采样),vindex()
使用特殊约定,即Ellipsis
表示批次维度(因此...
只能出现在左侧,不能出现在中间或右侧)。假设x
具有事件维度3。那么我们可以写成:old_batch_shape = x.shape[:-3] old_event_shape = x.shape[-3:] xij = Vindex(x)[..., i, :, j] # The ... denotes unknown batch shape. new_batch_shape = broadcast_shape(old_batch_shape, i.shape, j.shape) new_event_shape = (x.size(1),) assert xij.shape = new_batch_shape + new_event_shape
请注意,这种对
Ellipsis
的特殊处理与NEP [1]不同。正式地,此函数假设:
每个参数要么是
Ellipsis
,要么是slice(None)
,一个整数,或者一个批处理的torch.LongTensor
(即具有空事件形状)。此函数不支持非平凡的切片或torch.BoolTensor
掩码。Ellipsis
只能作为args[0]
出现在左侧。如果
args[0] is not Ellipsis
那么tensor
不是 批处理的,并且其事件维度等于len(args)
。如果
args[0] is Ellipsis
那么tensor
是批处理的,并且 其事件维度等于len(args[1:])
。tensor
的事件维度左侧的维度 被视为批处理维度,并将与张量参数的维度进行广播。
请注意,如果没有任何参数是具有
.dim() > 0
的张量,那么此函数的行为类似于标准索引:if not any(isinstance(a, torch.Tensor) and a.dim() for a in args): assert Vindex(x)[args] == x[args]
参考文献
- [1] https://www.numpy.org/neps/nep-0021-advanced-indexing.html
介绍了
vindex
作为向量化索引的辅助工具。 Pyro的实现与提议的符号x.vindex[]
类似,除了对Ellipsis
的处理略有不同。
- Parameters
tensor (torch.Tensor) – 一个要被索引的张量。
args (tuple) – 一个索引,作为
__getitem__
的参数。
- Returns
对
tensor[args]
的非标准解释。- Return type
- class Vindex(tensor)[source]¶
基础类:
object
围绕
vindex()
的便捷包装器。以下是等价的:
Vindex(x)[..., i, j, :] vindex(x, (Ellipsis, i, j, slice(None)))
- Parameters
tensor (torch.Tensor) – 一个要被索引的张量。
- Returns
一个具有特殊
__getitem__()
方法的对象。
张量收缩¶
- contract(equation, *operands, **kwargs)[source]¶
围绕
opt_einsum.contract()
的包装器,可选择使用Pyro的廉价优化器,并可选择缓存收缩路径。- Parameters
cache_path (bool) – 是否缓存收缩路径。 默认为 True。
- contract_expression(equation, *shapes, **kwargs)[source]¶
围绕
opt_einsum.contract_expression()
的包装器,可选择使用Pyro的廉价优化器,并可选择缓存收缩路径。- Parameters
cache_path (bool) – 是否缓存收缩路径。 默认为 True。
- einsum(equation, *operands, **kwargs)[source]¶
通过张量变量消除的广义板求和积算法。
这以两种方式概括了
contract()
:允许多个输出,并且可以共享中间结果。
输入和输出可以沿着
plates
中给出的符号进行排列;沿着plates
的减少是乘积减少。
理解此函数的最佳方法是尝试下面的示例,这些示例展示了如何将
einsum()
调用实现为对contract()
的多次调用(这通常更昂贵)。为了说明多个输出,请注意以下是等价的:
z1, z2, z3 = einsum('ab,bc->a,b,c', x, y) # multiple outputs z1 = contract('ab,bc->a', x, y) z2 = contract('ab,bc->b', x, y) z3 = contract('ab,bc->c', x, y)
为了说明plated输入,请注意以下是等价的:
assert len(x) == 3 and len(y) == 3 z = einsum('ab,ai,bi->b', w, x, y, plates='i') z = contract('ab,a,a,a,b,b,b->b', w, *x, *y)
当求和维度 a 总是与板维度 i 一起出现时,a 对应于 a 的每个切片的独特符号。因此,以下是等价的:
assert len(x) == 3 and len(y) == 3 z = einsum('ai,ai->', x, y, plates='i') z = contract('a,b,c,a,b,c->', *x, *y)
当这样的求和维度出现在输出中时,它必须伴随其所有的板维度,例如以下是等价的:
assert len(x) == 3 and len(y) == 3 z = einsum('abi,abi->bi', x, y, plates='i') z0 = contract('ab,ac,ad,ab,ac,ad->b', *x, *y) z1 = contract('ab,ac,ad,ab,ac,ad->c', *x, *y) z2 = contract('ab,ac,ad,ab,ac,ad->d', *x, *y) z = torch.stack([z0, z1, z2])
请注意,输出中的每个板块切片在所有输入的所有板块切片中都是多线性的,因此例如批量矩阵乘法将不会使用
plates
来实现,因此以下内容都是等效的:xy = einsum('abc,acd->abd', x, y, plates='') xy = torch.stack([xa.mm(ya) for xa, ya in zip(x, y)]) xy = torch.bmm(x, y)
在所有有效的方程中,有些计算是输入张量大小的多项式,而其他计算是输入张量大小的指数。当计算是指数时,此函数会引发
NotImplementedError
。- Parameters
- Returns
一个请求形状的张量元组,每个输出一个条目。
- Return type
- Raises
ValueError – 如果张量大小不匹配或输出请求了一个没有该维度板的维度。
NotImplementedError – 如果收缩在任何输入张量的大小上会导致指数级成本。
高斯收缩¶
- class Gaussian(log_normalizer: torch.Tensor, info_vec: torch.Tensor, precision: torch.Tensor)[source]¶
基础类:
object
非归一化的高斯分布。
这表示一个任意的半定二次函数,可以解释为一个秩不足的缩放高斯分布。精度矩阵可能具有零特征值,因此可能无法直接使用协方差矩阵。
- Parameters
log_normalizer (torch.Tensor) – 一个归一化常数,主要用于在收缩过程中跟踪归一化项。
info_vec (torch.Tensor) – 信息向量,它是均值的缩放版本
info_vec = precision @ mean
。我们使用这种表示法来使高斯收缩快速且稳定。precision (torch.Tensor) – 该高斯分布的精度矩阵。
- property batch_shape¶
- expand(batch_shape) pyro.ops.gaussian.Gaussian [source]¶
- reshape(batch_shape) pyro.ops.gaussian.Gaussian [source]¶
- __getitem__(index) pyro.ops.gaussian.Gaussian [来源]¶
索引到高斯分布的批次形状中。
- static cat(parts, dim=0) pyro.ops.gaussian.Gaussian [source]¶
沿给定的批次维度连接高斯分布列表。
- event_pad(left=0, right=0) pyro.ops.gaussian.Gaussian [source]¶
沿事件维度填充。
- event_permute(perm) pyro.ops.gaussian.Gaussian [source]¶
沿着事件维度进行置换。
- __add__(other: Union[pyro.ops.gaussian.Gaussian, int, float, torch.Tensor]) pyro.ops.gaussian.Gaussian [source]¶
在log-density空间中添加两个高斯分布。
- log_density(value: torch.Tensor) torch.Tensor [source]¶
评估此高斯分布在某点的对数密度:
-0.5 * value.T @ precision @ value + value.T @ info_vec + log_normalizer
这主要用于测试。
- rsample(sample_shape=torch.Size([]), noise: Optional[torch.Tensor] = None) torch.Tensor [source]¶
重新参数化的采样器。
- condition(value: torch.Tensor) pyro.ops.gaussian.Gaussian [source]¶
在这个高斯分布的尾部子集上进行条件化。 这应该满足:
g.condition(y).dim() == g.dim() - y.size(-1)
请注意,由于这是一个非归一化的高斯分布,我们在结果中包含了
y
的密度。因此,condition()
类似于functools.partial
的参数绑定:left = x[..., :n] right = x[..., n:] g.log_density(x) == g.condition(right).log_density(left)
- left_condition(value: torch.Tensor) pyro.ops.gaussian.Gaussian [source]¶
将这个高斯条件设置在其状态的前导子集上。 这应该满足:
g.condition(y).dim() == g.dim() - y.size(-1)
请注意,由于这是一个非归一化的高斯分布,我们在结果中包含了
y
的密度。因此,condition()
类似于functools.partial
的参数绑定:left = x[..., :n] right = x[..., n:] g.log_density(x) == g.left_condition(left).log_density(right)
- marginalize(left=0, right=0) pyro.ops.gaussian.Gaussian [source]¶
在事件维度的任一侧边缘化变量:
g.marginalize(left=n).event_logsumexp() = g.logsumexp() g.marginalize(right=n).event_logsumexp() = g.logsumexp()
对于数据
x
:- g.condition(x).event_logsumexp()
= g.marginalize(left=g.dim() - x.size(-1)).log_density(x)
- event_logsumexp() torch.Tensor [来源]¶
整合所有潜在状态(即在事件维度上操作)。
- class AffineNormal(matrix, loc, scale)[source]¶
基础类:
object
表示一个随机变量
Y
的条件对角正态分布,其均值是随机变量X
的仿射函数。 因此,X
的似然函数为:AffineNormal(matrix, loc, scale).condition(y).log_density(x)
这相当于:
Normal(x @ matrix + loc, scale).to_event(1).log_prob(y)
- Parameters
matrix (torch.Tensor) – 从
X
到Y
的转换。 应具有最右侧的形状(x_dim, y_dim)
。loc (torch.Tensor) –
Y
的均值的常数偏移量。 应具有最右侧形状(y_dim,)
。scale (torch.Tensor) –
Y
的标准差。 应具有最右侧形状(y_dim,)
。
- property batch_shape¶
- left_condition(value)[source]¶
如果
value.size(-1) == x_dim
,这将返回一个具有event_dim=1
的正态分布。应用此方法后,绘制样本的成本是O(y_dim)
而不是O(y_dim ** 3)
。
- rsample(sample_shape=torch.Size([]), noise: Optional[torch.Tensor] = None) torch.Tensor [source]¶
重新参数化的采样器。
- mvn_to_gaussian(mvn)[source]¶
将多元正态分布转换为高斯分布。
- Parameters
mvn (MultivariateNormal) – 一个多元正态分布。
- Returns
一个等价的高斯对象。
- Return type
- matrix_and_gaussian_to_gaussian(matrix: torch.Tensor, y_gaussian: pyro.ops.gaussian.Gaussian) pyro.ops.gaussian.Gaussian [source]¶
构建一个条件高斯分布
p(y|x)
,其中y - x @ matrix ~ y_gaussian
。- Parameters
matrix (torch.Tensor) – 一个右作用的变换矩阵。
y_gaussian (Gaussian) – 一个关于噪声的分布,噪声为
y - x@matrix
。
- Return type
- matrix_and_mvn_to_gaussian(matrix, mvn)[source]¶
将噪声仿射函数转换为高斯函数。噪声仿射函数定义为:
y = x @ matrix + mvn.sample()
- gaussian_tensordot(x: pyro.ops.gaussian.Gaussian, y: pyro.ops.gaussian.Gaussian, dims: int = 0) pyro.ops.gaussian.Gaussian [source]¶
计算两个高斯函数的积分:
(x @ y)(a,c) = log(integral(exp(x(a,b) + y(b,c)), b)),
其中 x 是变量 (a,b) 上的高斯分布,y 是变量 (b,c) 上的高斯分布,(a,b,c) 可以是零个或多个变量的集合,dims 是 b 的大小。
- Parameters
x – 一个高斯实例
y – 一个高斯实例
dims – 要收缩的变量数量
- sequential_gaussian_tensordot(gaussian: pyro.ops.gaussian.Gaussian) pyro.ops.gaussian.Gaussian [源代码]¶
对高斯
x
进行积分,其最右边的批次维度是时间,计算如下:x[..., 0] @ x[..., 1] @ ... @ x[..., T-1]
- sequential_gaussian_filter_sample(init: pyro.ops.gaussian.Gaussian, trans: pyro.ops.gaussian.Gaussian, sample_shape: Tuple[int, ...] = (), noise: Optional[torch.Tensor] = None) torch.Tensor [source]¶
通过并行扫描前向滤波后向采样,从高斯的马尔可夫积中绘制一个重新参数化的样本。
- Parameters
init (Gaussian) – 表示初始状态的高斯分布。
trans (Gaussian) – 一个表示一系列状态转换的高斯分布,时间作为最右边的批次维度。这个高斯分布的事件维度必须是
init
的两倍:trans.dim() == 2 * init.dim()
。sample_shape (tuple) – 一个可选的额外形状,用于绘制样本。
噪声 (torch.Tensor) – 一个可选的标准白噪声张量,形状为
sample_shape + batch_shape + (duration, state_dim)
,其中duration = 1 + trans.batch_shape[-1]
是要采样的时间点数,state_dim = init.dim()
是状态维度。 这对于计算均值(传递零)、变化温度(传递缩放噪声)和对立采样(传递cat([z,-z])
)非常有用。
- Returns
一个重新参数化的样本,形状为
sample_shape + batch_shape + (duration, state_dim)
。- Return type
统计工具¶
- gelman_rubin(input, chain_dim=0, sample_dim=1)[source]¶
计算样本链的R-hat。要求
input.size(sample_dim) >= 2
和input.size(chain_dim) >= 2
。- Parameters
输入 (torch.Tensor) – 输入张量。
chain_dim (int) – 链的维度。
sample_dim (int) – 样本维度。
- Returns torch.Tensor
input
的 R-hat。
- split_gelman_rubin(input, chain_dim=0, sample_dim=1)[source]¶
计算样本链上的R-hat。要求
input.size(sample_dim) >= 4
。- Parameters
输入 (torch.Tensor) – 输入张量。
chain_dim (int) – 链的维度。
sample_dim (int) – 样本维度。
- Returns torch.Tensor
分割
input
的 R-hat。
- autocorrelation(input, dim=0)[source]¶
计算样本在维度
dim
上的自相关。参考:https://en.wikipedia.org/wiki/Autocorrelation#Efficient_computation
- Parameters
输入 (torch.Tensor) – 输入张量。
dim (int) – 计算自相关的维度。
- Returns torch.Tensor
input
的自相关。
- autocovariance(input, dim=0)[source]¶
计算样本在维度
dim
上的自协方差。- Parameters
输入 (torch.Tensor) – 输入张量。
dim (int) – 计算自相关的维度。
- Returns torch.Tensor
input
的自相关。
- effective_sample_size(input, chain_dim=0, sample_dim=1)[source]¶
计算输入的有效样本大小。
参考:
- [1] Introduction to Markov Chain Monte Carlo,
查尔斯·J·盖尔
- [2] Stan Reference Manual version 2.18,
斯坦开发团队
- Parameters
输入 (torch.Tensor) – 输入张量。
chain_dim (int) – 链的维度。
sample_dim (int) – 样本维度。
- Returns torch.Tensor
input
的有效样本量。
- resample(input, num_samples, dim=0, replacement=False)[source]¶
从
input
的维度dim
中抽取num_samples
个样本。- Parameters
输入 (torch.Tensor) – 输入张量。
num_samples (int) – 从
input
中抽取的样本数量。dim (int) – 从
input
中提取的维度。
- Returns torch.Tensor
从
input
中随机抽取的样本。
- quantile(input, probs, dim=0)[source]¶
计算
input
在probs
处的分位数。如果probs
是标量,输出将在dim
处被压缩。- Parameters
输入 (torch.Tensor) – 输入张量。
probs (list) – 分位数位置。
dim (int) – 从
input
中获取分位数的维度。
- Returns torch.Tensor
input
在probs
的分位数。
- weighed_quantile(input: torch.Tensor, probs: Union[List[float], Tuple[float, ...], torch.Tensor], log_weights: torch.Tensor, dim: int = 0) torch.Tensor [源代码]¶
计算在
probs
处的加权input
样本的分位数。- Parameters
输入 (torch.Tensor) – 输入张量。
probs (list) – 分位数位置。
log_weights (torch.Tensor) – 样本权重张量。
dim (int) – 从
input
中获取分位数的维度。
- Returns torch.Tensor
input
在probs
的分位数。
示例:
>>> from pyro.ops.stats import weighed_quantile >>> import torch >>> input = torch.Tensor([[10, 50, 40], [20, 30, 0]]) >>> probs = torch.Tensor([0.2, 0.8]) >>> log_weights = torch.Tensor([0.4, 0.5, 0.1]).log() >>> result = weighed_quantile(input, probs, log_weights, -1) >>> torch.testing.assert_close(result, torch.Tensor([[40.4, 47.6], [9.0, 26.4]]))
- pi(input, prob, dim=0)[source]¶
计算百分位数区间,该区间为区间的每个尾部分配相等的概率质量。
- Parameters
输入 (torch.Tensor) – 输入张量。
prob (float) – 区间内样本的概率质量。
dim (int) – 从
input
计算百分位区间的维度。
- Returns torch.Tensor
input
在probs
的分位数。
- hpdi(input, prob, dim=0)[source]¶
计算“最高后验密度区间”,这是概率质量
prob
的最窄区间。- Parameters
输入 (torch.Tensor) – 输入张量。
prob (float) – 区间内样本的概率质量。
dim (int) – 从
input
计算百分位区间的维度。
- Returns torch.Tensor
input
在probs
的分位数。
- waic(input, log_weights=None, pointwise=False, dim=0)[source]¶
计算“广泛适用/渡边-赤池信息准则”(WAIC)及其对应的有效参数数量。
参考:
[1] WAIC 和 Stan 中的交叉验证, Aki Vehtari, Andrew Gelman
- Parameters
输入 (torch.Tensor) – 输入张量,它是模型的对数似然。
log_weights (torch.Tensor) – 沿着
dim
维度的样本权重。dim (int) –
input
的样本维度。
- Returns tuple
WAIC 和有效参数数量的元组。
- fit_generalized_pareto(X)[source]¶
给定一个假设从广义帕累托分布中抽取的数据集X,使用参考文献[1]中描述的技术变体来估计分布参数k和sigma,如参考文献[2]中所述。
参考文献 [1] ‘一种新的高效广义帕累托分布估计方法。’ 张杰和斯蒂芬斯,M.A. (2009). [2] ‘帕累托平滑重要性采样。’ 阿基·维塔里,安德鲁·格尔曼,乔纳·加布里
- Parameters
torch.Tensor – 输入数据 X
- Returns tuple
浮点数元组 (k, sigma) 对应于拟合参数
- crps_empirical(pred, truth)[source]¶
计算一组样本
pred
和真实数据truth
之间的负连续排名概率评分 CRPS* [1]。这使用了一个n log(n)
时间算法来计算一个量,该量在样本数量n
上原本具有二次复杂度:CRPS* = E|pred - truth| - 1/2 E|pred - pred'| = (pred - truth).abs().mean(0) - (pred - pred.unsqueeze(1)).abs().mean([0, 1]) / 2
请注意,对于单个样本,这简化为绝对误差。
参考文献
- [1] Tilmann Gneiting, Adrian E. Raftery (2007)
严格适当的评分规则、预测和估计 https://www.stat.washington.edu/raftery/Research/PDF/Gneiting2007jasa.pdf
- Parameters
pred (torch.Tensor) – 一组在最右侧维度上批处理的样本预测。 这应该具有形状
(num_samples,) + truth.shape
。真实值 (torch.Tensor) – 一个包含真实观测值的张量。
- Returns
一个形状为
truth.shape
的张量。- Return type
- energy_score_empirical(pred: torch.Tensor, truth: torch.Tensor) torch.Tensor [source]¶
计算一组多元样本
pred
和真实数据向量truth
之间的负能量评分 ES*(参见[1]中的公式22)。运行时间与样本数量n
的平方成正比。在单变量样本的情况下,输出与 CRPS 一致:ES* = E|pred - truth| - 1/2 E|pred - pred'|
请注意,对于单个样本,这简化为样本
pred
和truth
之间差异的欧几里得范数。这是一个严格的适当分数,因此对于根据分布 \(P\) 分布的
pred
和根据分布 \(Q\) 分布的truth
,我们有 \(ES^{*}(P,Q) \ge ES^{*}(Q,Q)\),当且仅当 \(P=Q\) 时等号成立,即如果 \(P\) 和 \(Q\) 具有相同的多元分布(仅当 \(P\) 和 \(Q\) 具有相同的边缘分布时,等号成立是不够的)。参考文献
- [1] Tilmann Gneiting, Adrian E. Raftery (2007)
严格适当的评分规则、预测和估计 https://www.stat.washington.edu/raftery/Research/PDF/Gneiting2007jasa.pdf
- Parameters
pred (torch.Tensor) – 一组样本预测,批处理在第二个最左边的维度上。 最左边的维度是多变量样本的维度。
真实值 (torch.Tensor) – 一个与
pred
形状相同的真实观测值的张量,除了最左边的第二个维度可以有任何值或被省略。
- Returns
一个形状为
truth.shape
的张量。- Return type
流式统计¶
- class StreamingStats[source]¶
基础类:
abc.ABC
用于张量树的可流式统计的抽象基类。
派生类必须实现
update()
,merge()
, 和get()
.- abstract update(sample) None [source]¶
从单个样本更新状态。
这会改变
self
并且不返回任何内容。更新应该是独立于顺序的,即样本应该是可交换的。- Parameters
sample – 一个样本值,它是一个嵌套的字典,包含
torch.Tensor
叶子节点。它可以有任意的嵌套和 形状,但假设在调用.update()
时形状保持不变。
- abstract merge(other) pyro.ops.streaming.StreamingStats [source]¶
选择两个聚合统计量,例如来自不同的MCMC链。
这是一个纯函数:它返回一个新的
StreamingStats
对象,并且不会修改self
或other
。- Parameters
other – 另一个相同类型的流统计实例。
- class StatsOfDict(types: Dict[Hashable, Callable[[], pyro.ops.streaming.StreamingStats]] = {}, default: Callable[[], pyro.ops.streaming.StreamingStats] = <class 'pyro.ops.streaming.CountStats'>)[source]¶
基础类:
pyro.ops.streaming.StreamingStats
统计样本,这些样本是具有恒定键集的字典。
例如,以下是等价的:
# Version 1. Hand encode statistics. >>> a_stats = CountStats() >>> b_stats = CountMeanStats() >>> a_stats.update(torch.tensor(0.)) >>> b_stats.update(torch.tensor([1., 2.])) >>> summary = {"a": a_stats.get(), "b": b_stats.get()} # Version 2. Collect samples into dictionaries. >>> stats = StatsOfDict({"a": CountStats, "b": CountMeanStats}) >>> stats.update({"a": torch.tensor(0.), "b": torch.tensor([1., 2.])}) >>> summary = stats.get() >>> summary {'a': {'count': 1}, 'b': {'count': 1, 'mean': tensor([1., 2.])}}
- Parameters
default – 字典值的统计默认类型。默认为廉价的
CountStats
。types (dict) – 字典映射键到应记录的统计类型,该统计类型对应于该键的值。
- class StackStats[source]¶
基础类:
pyro.ops.streaming.StreamingStats
统计收集一系列张量到一个单一的堆叠张量。
- update(sample: torch.Tensor) None [source]¶
- class CountStats[source]¶
基础类:
pyro.ops.streaming.StreamingStats
仅统计样本数量的统计跟踪。
例如:
>>> stats = CountStats() >>> stats.update(torch.randn(3, 3)) >>> stats.get() {'count': 1}
- merge(other: pyro.ops.streaming.CountStats) pyro.ops.streaming.CountStats [来源]¶
- class CountMeanStats[source]¶
基础类:
pyro.ops.streaming.StreamingStats
统计跟踪单个
torch.Tensor
的计数和平均值。- update(sample: torch.Tensor) None [来源]¶
- class CountMeanVarianceStats[source]¶
基础类:
pyro.ops.streaming.StreamingStats
统计跟踪单个
torch.Tensor
的计数、均值和(对角线)方差。- update(sample: torch.Tensor) None [source]¶
状态空间模型和GP工具¶
- class MaternKernel(nu=1.5, num_gps=1, length_scale_init=None, kernel_scale_init=None)[源代码]¶
-
提供了将具有Matern核的单变量高斯过程(GPs)表示为状态空间模型的基本构建块。
- Parameters
nu (float) – Matern核的阶数(0.5、1.5或2.5之一)
num_gps (int) – GP的数量
length_scale_init (torch.Tensor) – 可选的num_gps维向量,用于初始化长度尺度
kernel_scale_init (torch.Tensor) – 可选的 num_gps 维向量,用于初始化核尺度
参考文献
- [1] Kalman Filtering and Smoothing Solutions to Temporal Gaussian Process Regression Models,
Jouni Hartikainen 和 Simo Sarkka。
- [2] Stochastic Differential Equation Methods for Spatio-Temporal Gaussian Process Regression,
阿尔诺·索林。
- transition_matrix(dt)[source]¶
计算GP潜在空间的(指数化的)转移矩阵。 生成的矩阵布局为(num_gps, old_state, new_state),即此 矩阵从右侧乘以状态。
详情请参见参考文献[1]中的第5节。
- Parameters
dt (float) – GP潜在空间演化的时间间隔。
- Returns torch.Tensor
一个形状为 (num_gps, state_dim, state_dim) 的转移矩阵的三维张量。
- stationary_covariance()[source]¶
计算稳态协方差。参见参考文献[2]中的公式3.26。
- Returns torch.Tensor
一个形状为 (num_gps, state_dim, state_dim) 的协方差矩阵的三维张量。
- process_covariance(A)[source]¶
给定一个使用transition_matrix计算的转移矩阵A,按照参考文献[2]中的公式3.11计算过程协方差。
- Returns torch.Tensor
一个形状为 (num_gps, state_dim, state_dim) 的批量协方差矩阵