python - pymc3 中具有(大)时间^2 项的分层模型中的 MCMC 收敛

Question

我有一个分层 Logit，可以随着时间的推移进行观察。正在关注Carter 2010 ，我添加了时间、时间^2 和时间^3 术语。在添加时间变量之前，模型会使用 Metropolis 或 NUTS 进行混合。 HamiltonianMC 失败。 NUTS 和 Metropolis 也与时间合作。但是 NUTS 和 Metropolis 在 time^2 和 time^3 上失败了，但它们的失败方式不同且令人费解。然而，与其他因更明显的模型规范原因而失败的模型不同，ADVI 仍然给出了估计(并且 ELBO 不是 inf)。

• NUTS 要么提前停止(上次达到 60 次迭代)，要么进展太快并返回一个空跟踪图，其中包含有关变量名的错误。
• Metropolis 会立即因维度不匹配错误而出错。看起来像 this github error 中的那个，但我使用的是伯努利结果，而不是负二项式。错误的结尾如下所示:ValueError:输入维度不匹配。 (输入[0].shape[0] = 1，输入[4].shape[0] = 18)
• 当我尝试 HamiltonianMC 时，我得到一条空跟踪。它返回没有混合的起始值
• ADVI 给出平均值和标准差。
• 我大幅增加了 ADVI 迭代次数。它给出的起点非常接近，但 NUTS 仍然失败了。
• 我仔细检查了我正在运行的 pymc3 版本中是否已修复 github 问题。这是。

我的直觉是，这与 time^2 和 time^3 变量的大小有关，因为我正在查看一个很大的时间范围。 Time^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())

问题:

• 我对时间变量的大小和分布的直觉是否合理？
• 是否有方法可以更有效地对平方和立方项进行采样？我已经搜索过，但是您能否向我指出有关此问题的资源，以便我可以了解更多信息？
• Metropolis 和尺寸不匹配错误是怎么回事？
• NUTS 的跟踪图为空是怎么回事？通常当它停止时，跟踪直到停止工作。
• 是否有其他方法来处理可能更容易采样的时间？

我还没有发布玩具模型，因为没有数据很难复制。一旦使用模拟数据进行复制，我将添加一个玩具模型。但实际模型如下:

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:大都会错误

奇怪的语法问题。 Theano 似乎对同时具有恒定效应和随机效应的模型感到困惑。我在 data 中创建了一个等于 0 的常量，data['c']=0 并将其用作时间、time^2 和 time^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)

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

更新:已采样!现在将尝试一些使采样器更容易的建议，包括使用规范 suggested here 。但至少它正在发挥作用!

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
       + ...)