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

Question

我有一个相当简单的测试数据集，我试图适合pymc3。

通过traceplot生成的结果看起来像this. 基本上所有的参数看起来像有一个标准的“毛毛虫”为100次迭代，随后是750次迭代的平线，随后再次毛虫的痕迹。

最初的100次迭代发生在25,000次ADVI迭代和10,000次迭代迭代之后。如果我改变这些金额，我随机将/不会有这些不想要的稳定期。

我想知道如果任何人有任何建议我怎么能阻止这种情况发生 - 什么是造成它？

谢谢。

完整代码如下。简而言之，我使用相应的一组值y = a（j）* sin（phase）+ b（j）* sin（phase）生成一组'相位'（-pi-> pi）。 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

我会看分歧，如哈密尔顿蒙特卡洛的笔记和文献中所述，参见例如here和here。

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

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

祝你好运！

Answer 2

要部分地回答我的问题：这个玩了一会儿后，它看起来像这个问题可能是由于hyperprior标准差要0我不知道为什么算法应该卡在那里虽然（测试一个小的标准偏差不能那么罕见...）。

在任何情况下，两种解决方案，似乎缓解这个问题（尽管他们并不完全删除它）是：

1）加一个偏移量的标准偏差的定义。相反之前使用用于SD一个HalfCauchy或HalfNormal的例如为：

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

2）中，使用对数正态分布设定，使得0的可能性很小。