【发布时间】：2018-10-08 03:04:04
【问题描述】：
这本质上是Doing Bayesian Data Analysis, Second Edition (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 个 theta，但看到 20 行数据。
但是，如果对数据进行分组，每行代表单个硬币的最终统计数据，如下所示
grouped = observations.groupby(['mint', 'coin']).agg({'outcome': [np.sum, np.size]}).reset_index()
grouped.columns = ['mint', 'coin', 'heads', 'total']
并且最终的似然变量变成了Binomial，如下
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 第 1 版，
"BUGS 模型使用二项似然分布 正确，而不是对个体使用伯努利分布 试验。二项式的这种使用只是为了方便缩短 该程序。如果数据被指定为逐项试验结果 而不是完全正确，那么模型可以包括 试验循环并使用伯努利似然函数"
困扰我的是，在第一个示例中（一个铸币厂，几个硬币），看起来 PyMC3 可以处理单个观察而不是聚合观察。所以我相信第一种形式应该有效，但没有。
代码
http://nbviewer.jupyter.org/github/JWarmenhoven/DBDA-python/blob/master/Notebooks/Chapter%209.ipynb
参考文献
PyMC3 - Differences in ways observations are passed to model -> difference in results?
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
【问题讨论】：
标签： pymc3 hierarchical-bayesian