Иерархическая биноминальная модель PyMC3 - расхождения после настройки - PullRequest
Новая
       10

Иерархическая биноминальная модель PyMC3 - расхождения после настройки

0 голосов
/ 31 января 2019

Я пытаюсь использовать 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

Если кто-то может помочь с советами о том, как провести повторную параметризациюмодель, которая была бы отличной, спасибо!

...