调试 pymc 概率计算

Question

我试图通过复制给定here 的混合高斯示例来模拟指数混合。代码如下。我知道这里的推理有一些时髦的方面，但我的问题更多是关于如何在这样的模型中调试计算。

这个想法是它是三个指数的混合，其尺度参数取自 Gamma 分配给scales。但是，在ElemwiseCategoricalStep 期间，所有观测值都被分配给第零个指数。通过查看initial_assignments，您可以看到将观测值分配给指数分量最初是不同的，并且您可以看到所有观测值都分配给所有交互中的第零个分量，因为set(tr['exp'].flatten()) 仅包含 0。

我认为这是因为在ElemwiseCategoricalStep.astep 中的表达式array([logp(v * self.sh) for v in self.values]) 中分配给p 的所有值都是负无穷大。我想知道为什么会这样以及如何纠正它，但更重要的是，我想知道可以使用哪些工具来调试这种事情。有没有办法让我一步步计算logp(v * self.sh)，看看结果是怎么确定的？如果我尝试使用 pdb 执行此操作，我想我会在 theano.compile.function_module.Function.__call__ 中的 outputs = self.fn() 受到阻碍，我想我无法进入，因为它是一个原生函数。

即使知道如何为一组给定的模型参数计算 pdf 也是一个有用的开始。

import numpy as np
import pymc as pm
from pymc import Model, Gamma, Normal, Dirichlet, Exponential
from pymc import Categorical
from pymc import sample, Metropolis, ElemwiseCategoricalStep

durations = np.concatenate(
    [np.random.exponential(1/lam, 10)
     for lam in [1e-3,7e-5,2e-6]])

initial_assignments = np.random.randint(0, 3, len(durations))

print 'initial_assignments', initial_assignments

with Model() as model:
    scales = Gamma('hp', 1, 1, shape=3)
    props = Dirichlet('props', a=np.array([1., 1., 1.]), shape=3)
    category = Categorical('exp', p=props, shape=len(durations))
    points = Exponential('obs', lam=scales[category], observed=durations)
    step1 = pm.Metropolis(vars=[props,scales])
    step2 = ElemwiseCategoricalStep(var=category, values=[0,1,2])
    start = {'exp': initial_assignments,
             'hp':  np.ones(3),
             'props': np.ones(3),}
    tr = sample(3000, step=[step1, step2], start=start)

print set(tr['exp'].flatten())

Answer 1

很好的问题。您可以做的一件事是查看每个组件的 pdf。

模型和每个变量都应该有一个 .logp 和一个 .elemwise_logp 属性，它们返回一个可以接受点或参数值的函数。

因此，您可以说类似print scales.logp(start) 或print model.logp(start) 或print scales.dlogp()(start)。

目前，不幸的是，我认为您必须指定所有参数值（即使是不影响特定变量结果的参数值）。

Model、FreeRV 和 ObservedRV 都继承自 Factor，它提供了此功能并具有一些其他方法。您可能会想要非 fast 版本，因为它们在接受的论点种类上更宽容。

这有帮助吗？如果您对可能有助于调试的事情有其他想法，请告诉我。这是我们知道 pymc3 和 theano 需要做一些工作的领域。