2017-09-06 3 views
2

私はpymc3に適合しようとしているかなり単純なテストデータセットを持っています。 traceplotによって生成されたトレース値に(望ましくない)安定性の期間があるのはなぜですか?

結果は基本的にthis. 再び毛虫に続く750回の反復のための平坦な線が続く100回の反復のための標準的な「毛虫」は、そこにあるように、すべてのパラメータの外観の跡のようになります。

最初の100回の反復は、25,000回のADVI反復と1万回の繰り返し反復の後に発生します。これらの金額を変更すると、私は無作為に安定していない/しないだろう。

誰かが私がこのことが起こることをどうやって止めることができるのか、それを引き起こしていることについてアドバイスがありますか?

ありがとうございました。

完全なコードは以下のとおりです。要約すると、私は対応する値のセットy = a(j)* sin(位相)+ b(j)* sin(位相)を有する一組の「位相」(-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) 
+0

この投稿は実際に私が気付いたよりも微妙な問題を解決しています。私は直接持っていました:http://twiecki.github.io/blog/2017/02/08/bayesian-hierchical-non-centered/ – James

答えて

0

ハミルトニアンモンテカルロの注意事項および文献で説明したように、私は、例えば、herehereを参照してください、divergenciesで見てね。

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

おそらく、target_acceptを増やそうとすることができます。

幸運を祈る!

1

部分的に自分の質問に答えるために:問題は0に行くhyperprior標準偏差の原因である可能性がありますようにしばらくの間、これを再生した後、それが見えます、私はわかりませんなぜアルゴリズムがそこに詰まっているのでしょうか(小さな標準偏差をテストすることは珍しいことではありません...)。いずれの場合においても

、(彼らはそれを完全に削除しませんが)問題を軽減するように見える2の解決策は以下のとおりです。

1)標準偏差の定義にオフセットを追加します。例えば:

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

2)の代わりに前SDためHalfCauchy又はHalfNormalを用い、0はそうであるように設定対数正規分布を使用します。

関連する問題