0

我有以下模型，其中将杰弗里的几何分布先验用于其中一种分布。我有一个从其他分布派生的分布，并在 Mixture 分布中使用它。我正进入（状态

IndexError: axis 1 is out of bounds [-1, 1) error.

请在下面找到我的代码。

actual_visit_mod_90_days = df_visit_period['mod_90_days'].values

def jeffreys_logp(theta):
    #log(((1-theta)**-1) * theta**-0.5 )
    return -(T.log(1-theta) + 0.5 * T.log(theta))


with Model() as pattern_study:

    dist1 = Binomial('dist1',n=14,p=0.5)
    dist2 = DiscreteUniform('dist2',lower = 0, upper = 89)
    pp = pm.DensityDist('pp',jeffreys_logp, testval = 0.05)
    dist3 = Geometric('dist3',pp)

    someDervDist =  pm.DensityDist('someDervDist', lambda x: dist2-dist1+dist3 -1)
    likelihoodDist1 = DiscreteUniform.dist(lower=0, upper=89)

    likelihoodDist2 = pm.DensityDist.dist(lambda x: pm.logsumexp(someDervDist))

    p=0.05

    p_vd = Mixture('p_vd',[p,1-p],[likelihoodDist1,likelihoodDist2],observed = actual_visit_mod_90_days)

错误消息是

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-19-1c9d266976ce> in <module>()
     28     p=0.05
     29 
---> 30     p_vd = Mixture('p_vd',[p,1-p],[likelihoodDist1,likelihoodDist2],observed = actual_visit_mod_90_days)
     31 
     32     #print pm.distributions.generate_samples(postponeProbability)

C:\Users\mdevananda\Continuum\Anaconda2\lib\site-packages\pymc3\distributions\distribution.pyc in __new__(cls, name, *args, **kwargs)
     35             total_size = kwargs.pop('total_size', None)
     36             dist = cls.dist(*args, **kwargs)
---> 37             return model.Var(name, dist, data, total_size)
     38         else:
     39             raise TypeError("Name needs to be a string but got: {}".format(name))

C:\Users\mdevananda\Continuum\Anaconda2\lib\site-packages\pymc3\model.pyc in Var(self, name, dist, data, total_size)
    780                 var = ObservedRV(name=name, data=data,
    781                                  distribution=dist,
--> 782                                  total_size=total_size, model=self)
    783             self.observed_RVs.append(var)
    784             if var.missing_values:

C:\Users\mdevananda\Continuum\Anaconda2\lib\site-packages\pymc3\model.pyc in __init__(self, type, owner, index, name, data, distribution, total_size, model)
   1228 
   1229             self.missing_values = data.missing_values
-> 1230             self.logp_elemwiset = distribution.logp(data)
   1231             # The logp might need scaling in minibatches.
   1232             # This is done in `Factor`.

C:\Users\mdevananda\Continuum\Anaconda2\lib\site-packages\pymc3\distributions\mixture.pyc in logp(self, value)
    111         w = self.w
    112 
--> 113         return bound(logsumexp(tt.log(w) + self._comp_logp(value), axis=-1).sum(),
    114                      w >= 0, w <= 1, tt.allclose(w.sum(axis=-1), 1),
    115                      broadcast_conditions=False)

C:\Users\mdevananda\Continuum\Anaconda2\lib\site-packages\pymc3\distributions\mixture.pyc in _comp_logp(self, value)
     83         except AttributeError:
     84             return tt.stack([comp_dist.logp(value) for comp_dist in comp_dists],
---> 85                             axis=1)
     86 
     87     def _comp_means(self):

C:\Users\mdevananda\Continuum\Anaconda2\lib\site-packages\theano\tensor\basic.pyc in stack(*tensors, **kwargs)
   4583         dtype = scal.upcast(*[i.dtype for i in tensors])
   4584         return theano.tensor.opt.MakeVector(dtype)(*tensors)
-> 4585     return join(axis, *[shape_padaxis(t, axis) for t in tensors])
   4586 
   4587 

C:\Users\mdevananda\Continuum\Anaconda2\lib\site-packages\theano\tensor\basic.pyc in shape_padaxis(t, axis)
   4475     if not -ndim <= axis < ndim:
   4476         msg = 'axis {0} is out of bounds [-{1}, {1})'.format(axis, ndim)
-> 4477         raise IndexError(msg)
   4478     if axis < 0:
   4479         axis += ndim

IndexError: axis 1 is out of bounds [-1, 1)

对我来说，问题似乎出在 DensityDist.dist 函数上。我也不确定我是否在这里遗漏了什么。

编辑更新：鉴于 event1 的历史模式（actual_visit_mod_90_days），我正在尝试了解 event1 从开始日期起多少天后发生。最初，我将 someDervDist 定义为一个方程：

someDervDist = dist2-dist1+dist3-1

这给出了关于决定 event1 的事件的各种概率之间的关系。

likelihoodDist2 = pm.DensityDist.dist(lambda x: pm.logsumexp(someDervDist))

将给出 logsumexp(someDervDist.logp()) 我没有写为 .logp 因为我认为它会给出对数概率（嗯，我错了吗？）但是，我意识到这需要对数概率并将 someDervDist eqn 更改为DensityDist 分布（我也尝试过 Deterministic）。然后，我假设 p=0.05 的概率很小，这偏离了 event1 的周期性发生模型。基本上我正在寻找类似 Mixture([p, 1-p], [ DiscreteUniform(lower = 0, upper = 89), someDervDist]) 的东西。

4

1 回答 1

0

更新的答案，Constant毕竟使用发行版不是一个好主意，很快就会被弃用。所以也许一个小技巧会有所帮助，只需使用具有均值someDervDist和小方差的分布。

actual_visit_mod_90_days = [29, 76, 66, 66, 10, 63, 69, 75, 66, 87, 58]

with pm.Model() as pattern_study:

    dist1 = pm.Binomial('dist1', n=14, p=0.5)
    dist2 = pm.DiscreteUniform('dist2', lower = 0, upper = 89)
    pp = pm.DensityDist('pp', jeffreys_logp, testval = 0.5)
    dist3 = pm.Geometric('dist3', pp)

    someDervDist = pm.Deterministic('someDervDist', dist2-dist1+dist3-1)
    likelihoodDist2 =  pm.DiscreteUniform.dist(someDervDist-1, someDervDist + 1)
    likelihoodDist1 = pm.DiscreteUniform.dist(lower=0,upper=89)

    p= 0.05

    p_vd = pm.Mixture('p_vd', [p, 1-p],
                      [likelihoodDist1, likelihoodDist2],
                     observed = actual_visit_mod_90_days)
trace = pm.sample(1000, step=pm.Metropolis())

作为替代方案，您可以尝试使用模型的连续版本，如下所示：

with pm.Model() as pattern_study:

    dist1 = pm.Normal('dist1', 7, 1.5)
    dist2 = pm.Uniform('dist2', lower = 0, upper = 89)
    dist3 = pm.Exponential('dist3', 1)

    someDervDist = pm.Deterministic('someDervDist', dist2-dist1+dist3-1)
    likelihoodDist2 =  pm.Normal.dist(someDervDist, 0.1)
    likelihoodDist1 = pm.Uniform.dist(lower=0,upper=89)

    p = 0.05

    p_vd = pm.Mixture('p_vd', [p, 1-p],
                      [likelihoodDist1, likelihoodDist2],
                      observed = actual_visit_mod_90_days)
   trace = pm.sample(1000)
于 2018-01-20T06:41:21.977 回答