python - 调试pymc概率计算

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

这个想法是它是三个指数的混合，比例参数取自分配给scales的Gamma。但是，在ElemwiseCategoricalStep期间，所有观察值都分配给第零个指数。通过查看initial_assignments，您可以看到观察值对指数成分的分配最初是不同的，并且由于set(tr['exp'].flatten())仅包含0的事实，您可以看到所有观察值在所有交互中都被分配给第零成分。

我认为这是因为在p的表达式array([logp(v * self.sh) for v in self.values])中分配给ElemwiseCategoricalStep.astep的所有值都是负无穷大。我想知道这是为什么以及如何纠正它，但是，我想知道有什么工具可以用来调试这种事情。我有什么方法可以逐步完成logp(v * self.sh)的计算以查看结果如何确定？如果尝试使用pdb进行操作，我会在outputs = self.fn()的theano.compile.function_module.Function.__call__处受阻，我猜我不能介入，因为它是本机函数。

甚至知道如何为给定的一组模型参数计算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())

最佳答案

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

Model和每个变量都应同时具有.logp和.elemwise_logp属性，并且它们将返回可以采用点或参数值的函数。

因此，您可以说诸如print scales.logp(start)或print model.logp(start)或print scales.dlogp()(start)之类的内容。

就目前而言，我认为您很遗憾必须指定所有参数值（甚至不影响特定变量结果的值）。

Model，FreeRV和ObservedRV均从Factor继承，而提供此功能并具有其他一些方法。您可能希望使用非fast版本，因为它们接受的参数种类更宽容。

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