爲什麼跟蹤值具有（不需要的）穩定期？

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的可能性很小。