python - 调试pymc概率计算

Question

我尝试通过复制给定 here 的高斯混合示例来对指数混合进行建模。 。代码如下。我知道这里的推论有一些奇怪的方面，但我的问题更多是关于如何调试这样的模型中的计算。

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

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

即使知道如何计算给定模型参数集的 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 属性，并且它们返回一个可以采用点或参数值的函数。

因此你可以说类似 printscales.logp(start) 或 print model.logp(start) 或 printscales.dlogp()(start )。

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

Model、FreeRV 和 ObservedRV 均继承自 Factor它提供了此功能并具有一些其他方法。您可能需要非fast版本，因为它们接受的参数类型更加宽容。

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