pymc3 中具有（大）time^2 项的分层模型中的 MCMC 收敛

Question

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

NUTS 要么提前停止（上次它进行了 60 次迭代），要么进展太快并且 returns 一个空的跟踪图，其中有关于 varnames 的错误。
Metropolis 会立即出现尺寸不匹配错误。它看起来像 this github error 中的结果，但我使用的是伯努利结果，而不是负二项式。错误的结尾如下所示：ValueError: Input dimension mis-match. (input[0].shape[0] = 1, input[4].shape[0] = 18)
当我尝试 HamiltonianMC 时，我得到一个空的轨迹。它returns没有混合的起始值
ADVI 给出了平均值和标准差。
我大大增加了 ADVI 迭代次数。它给出了非常接近相同的起点，但 NUTS 仍然失败了。
我仔细检查了 github 问题的修复程序是否已在我运行ning 的 pymc3 版本中修复。这是。

我的直觉是，这与 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 完成了上述操作

我也尝试过 Metropolis 和 Hamiltonian。

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)

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

更新：它采样了！现在将尝试一些使采样器更容易的建议，包括使用规范。但至少它在工作！

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

pymc3 中具有（大）time^2 项的分层模型中的 MCMC 收敛

MCMC convergence in hierarchical model with (large) time^2 term in pymc3

python

bayesian

survival-analysis

mcmc

pymc3