调试pymc概率计算

Question

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

这个想法是，它是三个指数的混合物，具有从分配给scales的Gamma取得的比例参数。然而，在ElemwiseCategoricalStep期间，所有的观测值都被分配到第零指数。你可以看到，观测到的指数成分的分配是通过看initial_assignments开始多样化，你可以看到，所有的意见都分配到零组件的事实都interations是set(tr['exp'].flatten())只包含0

我认为这是因为在ElemwiseCategoricalStep.astep的表达式array([logp(v * self.sh) for v in self.values])中分配给p的所有值都是负无穷大。我想知道这是为什么，以及如何纠正它，但更多的是，我想知道哪些工具可用于调试这种事情。我有什么办法可以逐步计算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()) 

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需要一些工作的一个领域。