Почему значения трасс имеют периоды (нежелательной) стабильности?
У меня довольно простой набор тестовых данных, который я пытаюсь согласовать с pymc3.
Результат, генерируемый traceplot, выглядит примерно так . По сути, трассировка всех параметров выглядит так, как будто есть стандартная "гусеница" на 100 итераций, за которой следует плоская линия на 750 итераций, а затем снова гусеница.
Первые 100 итераций происходят после 25 000 итераций ADVI и 10 000 итераций настройки. Если я изменю эти суммы, у меня случайно будут / не будут эти периоды нежелательной стабильности.
Мне интересно, есть ли у кого-нибудь совет о том, как я могу предотвратить это - и чем это вызвано?
Благодарю.
Полный код ниже. Вкратце, я генерирую набор "фаз" (-pi -> pi) с соответствующим набором значений y = a(j)*sin(фаза) + b(j)*sin(фаза). a и b нарисованы для каждого субъекта j в произвольном порядке, но связаны друг с другом. Я тогда по существу пытаюсь соответствовать этой же модели.
Изменить: Вот аналогичный пример, работает для 25 000 итераций. Что-то идет не так с итерацией 20000.
#!/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)
2 ответа
Чтобы частично ответить на мой собственный вопрос: после того, как вы поиграете с этим некоторое время, похоже, что проблема может быть связана с тем, что стандартное отклонение гиперприора становится равным 0. Однако я не уверен, почему алгоритм должен застрять там (тестируя небольшое стандартное отклонение не может быть таким необычным...).
В любом случае, два решения, которые, кажется, облегчают проблему (хотя они не устраняют ее полностью):
1) Добавьте смещение к определениям стандартного отклонения. например:
offset = 1e-2
hb1_sigma_a = offset + pm.HalfCauchy('hb1_sigma_a', 4)
2) Вместо того, чтобы использовать HalfCauchy или HalfNormal для SD ранее, используйте набор распределения logNormal, так что 0 маловероятно.