Это, по сути, пример "Несколько монет от нескольких монетных дворов / бейсболистов" из Анализ байесовских данных, второе издание (DBDA2).Я считаю, что у меня есть код PyMC3, который функционально эквивалентен, но один работает, а другой нет.Это с PyMC версии 3.5.Более подробно,
Допустим, у меня есть следующие данные.Каждая строка является наблюдением:
observations_dict = {
'mint': [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
'coin': [0, 0, 0, 1, 1, 1, 2, 2, 2, 2, 3, 3, 4, 4, 4, 5, 5, 6, 6, 7],
'outcome': [1, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1]
}
observations = pd.DataFrame(observations_dict)
observations
Один монетный двор, несколько монет
Ниже, который реализует DBDA2, рисунок 9.7, работает просто отлично:
num_coins = observations['coin'].nunique()
coin_idx = observations['coin']
with pm.Model() as hierarchical_model:
# mint is characterized by omega and kappa
omega = pm.Beta('omega', 1., 1.)
kappa_minus2 = pm.Gamma('kappa_minus2', 0.01, 0.01)
kappa = pm.Deterministic('kappa', kappa_minus2 + 2)
# each coin is described by a theta
theta = pm.Beta('theta', alpha=omega*(kappa-2)+1, beta=(1-omega)*(kappa-2)+1, shape=num_coins)
# define the likelihood
y = pm.Bernoulli('y', theta[coin_idx], observed=observations['outcome'])
Много монетных дворов, много монет
Однако, как только это превращается в иерархическую модель (как видно из DBDA2, рисунок 9.13):
num_mints = observations['mint'].nunique()
mint_idx = observations['mint']
num_coins = observations['coin'].nunique()
coin_idx = observations['coin']
with pm.Model() as hierarchical_model2:
# Hyper parameters
omega = pm.Beta('omega', 1, 1)
kappa_minus2 = pm.Gamma('kappa_minus2', 0.01, 0.01)
kappa = pm.Deterministic('kappa', kappa_minus2 + 2)
# Parameters for mints
omega_c = pm.Beta('omega_c',
omega*(kappa-2)+1, (1-omega)*(kappa-2)+1,
shape = num_mints)
kappa_c_minus2 = pm.Gamma('kappa_c_minus2',
0.01, 0.01,
shape = num_mints)
kappa_c = pm.Deterministic('kappa_c', kappa_c_minus2 + 2)
# Parameters for coins
theta = pm.Beta('theta',
omega_c[mint_idx]*(kappa_c[mint_idx]-2)+1,
(1-omega_c[mint_idx])*(kappa_c[mint_idx]-2)+1,
shape = num_coins)
y2 = pm.Bernoulli('y2', p=theta[coin_idx], observed=observations['outcome'])
Ошибка:
ValueError: operands could not be broadcast together with shapes (8,) (20,)
, поскольку модель имеет 8 тэта на 8 монет, но видит 20 строк данных.
Однако, если данные сгруппированы таким образом, что каждая строка представляет окончательную статистику отдельной монеты, как в следующем
grouped = observations.groupby(['mint', 'coin']).agg({'outcome': [np.sum, np.size]}).reset_index()
grouped.columns = ['mint', 'coin', 'heads', 'total']
И конечная переменная вероятности изменяется на Бином, как показано ниже
num_mints = grouped['mint'].nunique()
mint_idx = grouped['mint']
num_coins = grouped['coin'].nunique()
coin_idx = grouped['coin']
with pm.Model() as hierarchical_model2:
# Hyper parameters
omega = pm.Beta('omega', 1, 1)
kappa_minus2 = pm.Gamma('kappa_minus2', 0.01, 0.01)
kappa = pm.Deterministic('kappa', kappa_minus2 + 2)
# Parameters for mints
omega_c = pm.Beta('omega_c',
omega*(kappa-2)+1, (1-omega)*(kappa-2)+1,
shape = num_mints)
kappa_c_minus2 = pm.Gamma('kappa_c_minus2',
0.01, 0.01,
shape = num_mints)
kappa_c = pm.Deterministic('kappa_c', kappa_c_minus2 + 2)
# Parameter for coins
theta = pm.Beta('theta',
omega_c[mint_idx]*(kappa_c[mint_idx]-2)+1,
(1-omega_c[mint_idx])*(kappa_c[mint_idx]-2)+1,
shape = num_coins)
y2 = pm.Binomial('y2', n=grouped['total'], p=theta, observed=grouped['heads'])
Все работает.Теперь последняя форма более эффективна и, как правило, предпочтительнее, но я считаю, что первая также должна работать.Поэтому я считаю, что это в первую очередь проблема PyMC3 (или, что более вероятно, ошибка пользователя).
Цитируя DBDA Edition 1,
"Модель BUGS использует биномиальное распределение правдоподобия для полной коррекции вместо использования распределения Бернулли для отдельных испытаний. Это использование биномиальногопросто удобство для сокращения программы. Если данные были указаны как результаты испытания по делу, а не как полностью правильные, тогда модель могла бы включать цикл испытания по делу и использовать функцию вероятности Бернулли "
Что меня беспокоит, так это то, что в самом первом примере (Один монетный двор, несколько монет) похоже, что PyMC3 может обрабатывать отдельные наблюдения вместо агрегированных наблюдений просто отлично.Поэтому я считаю, что первая форма должна работать, но не работает.
Код
http://nbviewer.jupyter.org/github/JWarmenhoven/DBDA-python/blob/master/Notebooks/Chapter%209.ipynb
Ссылки
PyMC3 - Различия в способах передачи наблюдений в модель -> разница в результатах?
https://discourse.pymc.io/t/pymc3-differences-in-ways-observations-are-passed-to-model-difference-in-results/501
http://www.databozo.com/deep-in-the-weeds-complex-hierarchical-models-in-pymc3
https://stats.stackexchange.com/questions/157521/is-this-correct-hierarchical-bernoulli-model