Я пытаюсь использовать pyMC3 для построения простой байесовской иерархической модели для некоторых экспериментальных данных.У меня есть два набора данных, но для одного из двух сэмплер не сходится, и я не могу найти решение.
Установки следующие:
- Существует два экспериментальных условия (невообразимо называемые A и B) и две группы индивидуумов, испытанных в одном из двух условий (группа A игруппа B).
- Каждый человек выполняет столько испытаний, сколько ему нужно, поэтому не у всех лиц одинаковое количество испытаний
- Каждое испытание имеет двоичный результат (1 или 0).
Данные о производительности каждого субъекта будут представлять собой строку из 1 и 0, и я хотел бы оценить базовую частоту 1 с каждого человека по наблюдаемым данным.
Поскольку по некоторым предметам у меня очень мало испытаний, я решил использовать иерархическую байесовскую модель (см. этот пример ).Модель, которую я решил использовать, основана на модели, показанной здесь [см. Также код ниже].
Теперь модель прекрасно работает для одного из двух наборов данных (B), носэмплер не сходится для другого.Я видел в Интернете, что возможное решение состоит в том, чтобы перейти на нецентрированную модель , но я не знаю, как реализовать это здесь.
Ниже приведен минимальный рабочий пример и результаты.
import numpy as np
import pymc3 as pm
import theano.tensor as tt
import matplotlib.pyplot as plt
def run():
# Define data
datasets_names = ['A', 'B']
number_of_individuals =[22, 17] # per experimental condition
# Number of trials and number of successes (1) of each individual
n_trials_A = [21, 15, 6, 5, 10, 6, 4, 6, 5, 7, 14, 12, 15, 4, 4, 6, 6, 9, 7, 6, 11, 10]
hits_A = [21, 14, 6, 0, 6, 6, 3, 6, 5, 6, 14, 9, 15, 4, 4, 5, 6, 8, 7, 4, 8, 10]
n_trials_B = [5, 5, 33, 4, 13, 18, 24, 8, 8, 9, 9, 7, 14, 8, 15, 9, 11]
hits_B = [2, 5, 26, 3, 7, 7, 13, 6, 1, 5, 4, 2, 7, 5, 9, 4, 1]
datasets = [(number_of_individuals[0], n_trials_A, hits_A), (number_of_individuals[1], n_trials_B, hits_B)]
# Model each dataset separately
for i, (m, n, h) in enumerate(datasets):
print('Modelling dataset: ', datasets_names[i])
# pyMC3 model
with pm.Model() as model:
# The model is from: https://docs.pymc.io/notebooks/hierarchical_partial_pooling.html
# Define hyperpriors
phi = pm.Uniform('phi', lower=0.0, upper=1.0)
kappa_log = pm.Exponential('kappa_log', lam=1.5)
kappa = pm.Deterministic('kappa', tt.exp(kappa_log))
# define second level of hierarchical model
thetas = pm.Beta('thetas', alpha=phi*kappa, beta=(1.0-phi)*kappa, shape=m)
# Likelihood
y = pm.Binomial('y', n=n, p=thetas, observed=h)
# Fit
trace = pm.sample(6000, tune=2000, nuts_kwargs={'target_accept': 0.95})
# Show traceplot
pm.traceplot(trace)
plt.show()
if __name__ == "__main__":
run()
Это то, что выводится на консоль при запуске кода:
Modeeling dataset: A
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [thetas, kappa_log, phi]
Sampling 4 chains: 100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 32000/32000 [00:52<00:00, 610.30draws/s]
There were 928 divergences after tuning. Increase `target_accept` or reparameterize.
There were 818 divergences after tuning. Increase `target_accept` or reparameterize.
There were 885 divergences after tuning. Increase `target_accept` or reparameterize.
There were 842 divergences after tuning. Increase `target_accept` or reparameterize.
The number of effective samples is smaller than 25% for some parameters.
Modeeling dataset: B
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [thetas, kappa_log, phi]
Sampling 4 chains: 100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 32000/32000 [00:35<00:00, 899.07draws/s]
Соответственно, трассировкидля набора данных A покажите, что сходимости не было.
Трассировка для набора данных A
Трассировка для набора данных B
Если кто-то может помочь с советами о том, как провести повторную параметризациюмодель, которая была бы отличной, спасибо!