机器算法验证 - pymc3：调优后的接受概率和分歧 - 吾爱随笔录

我在pymc3中编写了两个模型，我认为这很简单。

逻辑回归

第一个是实验中的逻辑回归，该实验对 2x2 因子设计中特定任务的正确和错误答案进行建模。这是它的样子：

with pm.Model() as m1:
    b_spe = pm.Normal('b_spe', 0, 5)
    b_noi = pm.Normal('b_noi', 0, 5)
    b_sn = pm.Normal('b_sn', 0, 5)
    sigma_sub = pm.HalfCauchy('sigma_sub', 1)
    sigma_run = pm.HalfCauchy('sigma_run', 1)

    a = pm.Normal('a', 0, 5)

    a_sub = pm.Normal('a_sub', 0, sigma_sub, shape=len(d['sub'].unique()))
    a_run = pm.Normal('a_run', 0, sigma_run, shape=len(d['run'].unique()))    

    p = pm.invlogit(a + a_sub[d["sub"].values] + a_run[d["run"].values] + b_spe*d.speech 
                    + b_noi * d.noise + b_sn * d.speech * d.noise)
    hits = pm.Binomial('hits', p=p, n=d.total, observed=d.hit)

    trace_m1 = pm.sample(5000, tune=5000, njobs=4, target_accept=.99)

采样后，我收到以下警告：

采样 4 条链：100%|██████████| 40000/40000 [03:38<00:00, 37.86draws/s] 接受概率与目标不匹配。它是 0.8942185942514685，但应该接近 0.8。尝试增加调整步骤的数量。某些参数的有效样本数小于 25%。

以下是跟踪图：

罗特(a) = 1.0002822268704348

至少在我看来，traceplots 看起来不错，而且 r-hat 非常接近 1。模型有问题吗？我能做些什么来修复它？

线性回归

这与上面的实验有关，但我使用我用学生 T 分布建模的 fMRI 数据用于较长的尾巴。它仍然是一个 2x2 因子设计，我还使用 LKJ 协方差矩阵对每个受试者的截距和斜率之间的协方差进行建模。此外，该模型是偏离中心的，灵感来自Statistical Rethinking第 13.23 章的模型。

with pm.Model() as m2:
    sd_dist = pm.HalfCauchy.dist(beta=2)
    pchol_sub = pm.LKJCholeskyCov('pchol_sub', eta=4, n=n_dim, sd_dist=sd_dist)
    chol = pm.expand_packed_triangular(n_dim, pchol_sub, lower=True)    

    a = pm.Normal('a', 0, 1)

    sigma_sub = pm.HalfCauchy('sigma_sub', 1)
    a_sub = pm.Normal('a_sub', 0, sigma_sub, shape=n_subs)  

    # includes beta for spe, noi, sn
    b_spe = pm.StudentT('b_spe', nu=3, mu=0, lam=1)
    b_noi = pm.StudentT('b_noi', nu=3, mu=0, lam=1)
    b_sn = pm.StudentT('b_sn', nu=3, mu=0, lam=1)

    # per subject
    b_s_sub = tt.dot(chol, pm.StudentT('b_s_sub', nu=3, mu=0, lam=1, shape=(1, n_subs)))
    b_n_sub = tt.dot(chol, pm.StudentT('b_n_sub', nu=3, mu=0, lam=1, shape=(1, n_subs)))
    b_sn_sub = tt.dot(chol, pm.StudentT('b_sn_sub', nu=3, mu=0, lam=1, shape=(1, n_subs)))

    A = a + a_sub[d['sub'].values]
    B_s = b_spe + b_s_sub[0, d['sub'].values]
    B_n = b_noi + b_n_sub[0, d['sub'].values]
    B_sn = b_sn + b_sn_sub[0, d['sub'].values]

    sigma = pm.HalfCauchy('sigma', 2)
    lam = sigma ** -2
    nu = pm.Exponential('ν_minus_one', 1/29.) + 1
    mu = pm.Deterministic('mu', A + B_s * d.speech + B_n * d.noise + B_sn * d.speech * d.noise)
    psc = pm.StudentT('psc', nu=nu, mu=mu, lam=lam, observed=d.psc_c)

    trace_m2 = pm.sample(5000, tune=5000, njobs=4, target_accept=0.99)

以下是我收到的警告/错误：

采样 4 条链：100%|██████████| 40000/40000 [04:55<00:00, 135.45draws/s] 调优后有2个分歧。增加target_accept或重新参数化。错误：pymc3：调整后有 2 个分歧。增加target_accept或重新参数化。接受概率与目标不匹配。它是 0.7183184493954338，但应该接近 0.8。尝试增加调整步骤的数量。

还有其他重复的警告几乎说了同样的话。对于每条链，我都会收到警告。

下面是跟踪图。

如果我在 2000:3000 之间放大，我会得到看起来很时髦的痕迹：

sigma_sub这是和之间的示例配对图a_sub_7。红点是分歧点。我受到Thomas Wiecki 的帖子的启发。

罗特(a) = 1.0014072110125143

再一次，这些痕迹在我看来很好（除了放大的部分）并且 r-hat 非常接近 1。我做错了什么？

这个模型的数据非常嘈杂，样本量是 17，不是很大。

真的很感激这里的一些帮助:-)谢谢！

编辑： 使用pm.sample(..., nuts_kwargs=dict(target_accept=0.95))摆脱了分歧。