为什么跟踪值有（不需要的）稳定期？

Question

我有一个相当简单的测试数据集，我正在尝试使用 pymc3 进行拟合。

traceplot 生成的结果类似于 this. 从本质上讲，所有参数的轨迹看起来像是有一个标准的 'caterpillar' 用于 100 次迭代，然后是一条平坦的线用于 750 次迭代，然后再次是毛毛虫。

最初的 100 次迭代发生在 25,000 次 ADVI 迭代和 10,000 次调整迭代之后。如果我更改这些数量，我随机 will/won 不会有这些不需要的稳定期。

我想知道是否有人对我如何阻止这种情况的发生有任何建议 - 是什么导致了这种情况？

谢谢。

完整代码如下。简而言之，我正在生成一组 'phases' (-pi -> pi) 以及一组相应的值 y = a(j)*sin(phase) + b(j)*sin(phase)。 a 和 b 是为每个主题 j 随机抽取的，但彼此相关。然后我基本上尝试适应这个相同的模型。

编辑：Here is a similar example, running for 25,000 iterations. Something goes wrong around iteration 20,000.

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
import matplotlib.pyplot as plt
import numpy as np
import pymc3 as pm
%matplotlib inline

np.random.seed(0)


n_draw = 2000
n_tune = 10000
n_init = 25000
init_string = 'advi'
target_accept = 0.95

##
# Generate some test data
# Just generates:
# x a vector of phases
# y a vector corresponding to some sinusoidal function of x
# subject_idx a vector corresponding to which subject x is

#9 Subjects
N_j = 9

#Each with 276 measurements
N_i = 276
sigma_y = 1.0
mean = [0.1, 0.1]
cov = [[0.01, 0], [0, 0.01]]  # diagonal covariance

x_sub = np.zeros((N_j,N_i))
y_sub = np.zeros((N_j,N_i))
y_true_sub = np.zeros((N_j,N_i))
ab_sub = np.zeros((N_j,2))
tuning_sub = np.zeros((N_j,1))
sub_ix_sub = np.zeros((N_j,N_i))
for j in range(0,N_j):
    aj,bj = np.random.multivariate_normal(mean, cov)
    #aj = np.abs(aj)
    #bj = np.abs(bj)
    xj = np.random.uniform(-1,1,size = (N_i,1))*np.pi
    xj = np.sort(xj)#for convenience
    yj_true = aj*np.sin(xj) + bj*np.cos(xj)
    yj = yj_true + np.random.normal(scale=sigma_y, size=(N_i,1))

    x_sub[j,:] = xj.ravel()
    y_sub[j,:] = yj.ravel()
    y_true_sub[j,:] = yj_true.ravel()

    ab_sub[j,:] = [aj,bj]
    tuning_sub[j,:] = np.sqrt(((aj**2)+(bj**2)))
    sub_ix_sub[j,:] = [j]*N_i

x = np.ravel(x_sub)
y = np.ravel(y_sub)
subject_idx = np.ravel(sub_ix_sub)
subject_idx = np.asarray(subject_idx, dtype=int)

##
# Fit model
hb1_model = pm.Model()
with hb1_model:
#    Hyperpriors
    hb1_mu_a = pm.Normal('hb1_mu_a', mu=0., sd=100)
    hb1_sigma_a = pm.HalfCauchy('hb1_sigma_a', 4)
    hb1_mu_b = pm.Normal('hb1_mu_b', mu=0., sd=100)
    hb1_sigma_b = pm.HalfCauchy('hb1_sigma_b', 4)

    # We fit a mixture of a sine and cosine with these two coeffieicents
    # allowed to be different for each subject
    hb1_aj = pm.Normal('hb1_aj', mu=hb1_mu_a, sd=hb1_sigma_a, shape = N_j)
    hb1_bj = pm.Normal('hb1_bj', mu=hb1_mu_b, sd=hb1_sigma_b, shape = N_j)

    # Model error
    hb1_eps = pm.HalfCauchy('hb1_eps', 5)

    hb1_linear = hb1_aj[subject_idx]*pm.math.sin(x) + hb1_bj[subject_idx]*pm.math.cos(x)
    hb1_linear_like = pm.Normal('y', mu = hb1_linear, sd=hb1_eps, observed=y)

with hb1_model:
    hb1_trace = pm.sample(draws=n_draw, tune = n_tune,
                      init = init_string, n_init = n_init,
                     target_accept = target_accept)

pm.traceplot(hb1_trace)

Answer 1

部分回答我自己的问题：玩了一段时间后，看起来问题可能是由于超先验标准偏差变为 0。我不确定为什么算法会卡在那里（测试一个小的标准差并不少见...)。

无论如何，有两个似乎可以缓解问题的解决方案（尽管它们并没有完全消除它）是：

1) 在标准偏差的定义中添加一个偏移量。例如：

offset = 1e-2
hb1_sigma_a = offset + pm.HalfCauchy('hb1_sigma_a', 4)

2) 不要先对 SD 使用 HalfCauchy 或 HalfNormal，而是使用 logNormal 分布集，这样 0 就不太可能了。

Answer 2

我会查看分歧，如哈密顿量 Monte Carlo 的注释和文献中所述，请参阅 here and here。

with model:
    np.savetxt('diverging.csv', hb1_trace['diverging'])

作为一个肮脏的解决方案，你可以尝试增加target_accept，也许。

祝你好运！

为什么跟踪值有（不需要的）稳定期？

Why do the trace values have periods of (unwanted) stability?

pymc

pymc3

mcmc