pymc3中(大)时间^2项层次模型的MCMC收敛性

Question

我有一个等级逻辑，它有随时间变化的观察结果。在卡特，2010年之后，我包含了时间、时间^2和时间^3项。在我添加时间变量之前，这个模型混合使用大都会或坚果。HamiltonianMC失败。坚果和大都会也随着时间的推移而起作用。但坚果和大都会在时间^2和时间^3中失败，但它们失败的方式不同，而且令人费解。然而，与其他由于更明显的模型规范原因而失败的模型不同，ADVI仍然给出了一个估计，( ELBO不是inf)。

• 坚果要么提前停止(上次达到60次迭代)，要么进展太快，返回一个空的跟踪图，其中包含一个关于变量名称的错误。
• 大都市误差立即与维数不匹配误差。它看起来像这个github错误中的一个，但是我使用的是Bernoulli结果，而不是一个负二项式。错误的结尾类似于：ValueError: Input dimension mis-match. (input[0].shape[0] = 1, input[4].shape[0] = 18)
• 当我尝试HamiltonianMC时，会得到一个空的跟踪。它返回没有混合的起始值。
• 平日给了我一个平均值和一个标准差。
• 我的高级迭代次数增加了很多。它给了相当接近相同的起点和坚果仍然失败。
• 我再次检查了github问题的修复在我正在运行的pymc3版本中是否已经到位。它是。

我的直觉是，这与时间^2和时间^3变量有多大有关，因为我正在查看一个很大的时间框架。时间^3从0开始，到64,000。

这是我到目前为止尝试过的抽样方法。请注意，我在测试时有很小的样本大小，因为运行时间太长了(如果它完成了的话)，我只是想让它完全被采样。一旦我找到一个起作用的，我会增加迭代次数

with my_model:
    mu,sds,elbo = pm.variational.advi(n=500000,learning_rate=1e-1)
    print(mu['mu_b'])
    step = pm.NUTS(scaling=my_model.dict_to_array(sds)**2,
                   is_cov=True)
    my_trace = pm.sample(500, 
                         step=step, 
                         start=mu,
                         tune=100)

我还用tune=1000做了上面的工作

我也试过大都会和哈密顿。

with my_model:
    my_trace = pm.sample(5000,step=pm.Metropolis())

with my_model:
    my_trace = pm.sample(5000,step=pm.HamiltonianMC())

问题：

• 关于时间变量的大小和分布，我的直觉是否合理？
• 是否有办法更有效地取样平方和立方条件？我已经搜索过了，但是你能告诉我关于这个的资源吗，这样我就可以了解更多的信息了吗？
• “大都会”和维数错配误差是怎么回事？
• 空的坚果跟踪图怎么回事？通常当它停下来的时候，跟踪直到货摊工作为止。
• 是否有其他方法可以更容易地处理时间？

我没有贴出玩具模型，因为没有数据很难复制。一旦我用模拟数据进行复制，我将添加一个玩具模型。但实际模式如下：

with pm.Model() as my_model:
    mu_b = pm.Flat('mu_b')
    sig_b = pm.HalfCauchy('sig_b',beta=2.5)
    b_raw = pm.Normal('b_raw',mu=0,sd=1,shape=n_groups)
    b = pm.Deterministic('b',mu_b + sig_b*b_raw)

    t1 = pm.Normal('t1',mu=0,sd=100**2,shape=1)
    t2 = pm.Normal('t2',mu=0,sd=100**2,shape=1)
    t3 = pm.Normal('t3',mu=0,sd=100**2,shape=1)

    est =(b[data.group.values]* data.x.values) +\
        (t1*data.t.values)+\
        (t2*data.t2.values)+\
        (t3*data.t3.values)

    y = pm.Bernoulli('y', p=tt.nnet.sigmoid(est), observed = data.y)

突破1:大都会错误

奇怪的语法问题。西亚诺似乎对一个既有恒定效应又有随机效应的模型感到困惑。我在数据中创建了一个常数，等于0，data['c']=0，并将其用作时间、时间^2和时间^3效果的索引，如下所示：

est =(b[data.group.values]* data.x.values) +\
            (t1[data.c.values]*data.t.values)+\
            (t2[data.c.values]*data.t2.values)+\
            (t3[data.c.values]*data.t3.values)

我不认为这是整个问题，但这是朝着正确方向迈出的一步。我敢打赌这就是为什么我的不对称规范不起作用的原因，如果是这样的话，我怀疑它可能会更好的样本。

更新:它取样了！现在将尝试一些简化采样器操作的建议，包括使用规范建议在这里。但至少它起作用了！

Answer 1

如果没有数据集可玩，很难给出一个明确的答案，但以下是我的最佳猜测：

对我来说，听到里面的三次多项式有点出乎意料。我还没看过那份报纸，所以我不能对它发表评论，但我认为这可能是你出现问题的原因。即使t3的值很小，也会对预测器产生巨大影响。为了保持这一点的合理性，我尝试更改参数化:首先，确保您的预测器是居中的(类似于data['t'] = data['t'] - data['t'].mean()，然后定义data.t2和data.t3)。然后尝试在t2和t3上设置一个更合理的先验。它们应该很小，所以也许可以尝试一下

t1 = pm.Normal('t1',mu=0,sd=1,shape=1)
t2 = pm.Normal('t2',mu=0,sd=1,shape=1)
t2 = t2 / 100
t3 = pm.Normal('t3',mu=0,sd=1,shape=1)
t3 = t3 / 1000

如果要查看其他模型，可以尝试将预测器建模为GaussianRandomWalk或高斯过程。

将pymc3更新到最新版本的候选版本也会有所帮助，采样器被改进了一点。

更新--我刚刚注意到您的模型中没有拦截项。除非有一个很好的理由你可能想添加

intercept = pm.Flat('intercept')
est = (intercept
       + b[..] * data.x
       + ...)