python - 调试 pymc 概率计算

Question

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

这个想法是，它是三个指数的混合，从 Gamma 中获取的尺度参数分配给 尺度。但是，在 ElemwiseCategoricalStep 期间，所有观察值都会分配给第零指数。通过查看 initial_assignments，您可以看到指数分量的观测值分配最初是不同的，并且您可以看到所有观测值都分配给了所有交互的第零分量，因为 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 属性，它们返回一个可以接受点或参数值的函数。

因此你可以这样说 print scales.logp(start) 或 print model.logp(start) 或 print scales.dlogp()(start )。

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

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

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