为什么这种分层泊松模型与生成数据中的真实参数不匹配？

Question

我正在尝试拟合分层泊松回归以估计每个组和全局的 time_delay。我对 pymc 是否自动将 log link 函数应用到 mu 还是我必须明确地这样做感到困惑：

with pm.Model() as model:
    alpha = pm.Gamma('alpha', alpha=1, beta=1)
    beta = pm.Gamma('beta', alpha=1, beta=1)

    a = pm.Gamma('a', alpha=alpha, beta=beta, shape=n_participants)

    mu = a[participants_idx]
    y_est = pm.Poisson('y_est', mu=mu, observed=messages['time_delay'].values)

    start = pm.find_MAP(fmin=scipy.optimize.fmin_powell)
    step = pm.Metropolis(start=start)
    trace = pm.sample(20000, step, start=start, progressbar=True)

下面的跟踪图显示了 a 的估计值。您可以看到 0 到 750 之间的组估计值。

当我使用 alpha 和 beta 的平均值作为参数绘制超参数伽马分布时，我开始感到困惑。下面的分布显示支持大约在 0 到 5 之间。在查看上面 a 的小组估计时，这不符合我的预期。 a代表什么？是 log(a) 还是其他？

感谢您的指点。

按照评论中的要求添加使用假数据的例子：这个例子只有一个组，所以应该更容易看出超参数是否可以合理地产生组的泊松分布。

test_data = []
model = []

for i in np.arange(1):
    # between 1 and 100 messages per conversation
    num_messages = np.random.uniform(1, 100)
    avg_delay = np.random.gamma(15, 1)
    for j in np.arange(num_messages):
        delay = np.random.poisson(avg_delay)

        test_data.append([i, j, delay, i])

    model.append([i, avg_delay])

model_df = pd.DataFrame(model, columns=['conversation_id', 'synthetic_mean_delay'])
test_df = pd.DataFrame(test_data, columns=['conversation_id', 'message_id', 'time_delay', 'participants_str'])
test_df.head()

# Estimate parameters of model using test data
# convert categorical variables to integer
le = preprocessing.LabelEncoder()
test_participants_map = le.fit(test_df['participants_str'])
test_participants_idx = le.fit_transform(test_df['participants_str'])
n_test_participants = len(test_df['participants_str'].unique())

with pm.Model() as model:
    alpha = pm.Gamma('alpha', alpha=1, beta=1)    
    beta = pm.Gamma('beta', alpha=1, beta=1)

    a = pm.Gamma('a', alpha=alpha, beta=beta, shape=n_test_participants)

    mu = a[test_participants_idx]

    y = test_df['time_delay'].values
    y_est = pm.Poisson('y_est', mu=mu, observed=y)

    start = pm.find_MAP(fmin=scipy.optimize.fmin_powell)
    step = pm.Metropolis(start=start)
    trace = pm.sample(20000, step, start=start, progressbar=True)

我不明白下面的超参数如何产生参数在 13 到 17 之间的泊松分布。

Answer 1

答案：pymc 使用与 scipy 不同的参数来表示 Gamma 分布。 scipy 使用 alpha & scale，而 pymc 使用 alpha 和 beta。以下模型按预期工作：

with pm.Model() as model:
    alpha = pm.Gamma('alpha', alpha=1, beta=1)    
    scale = pm.Gamma('scale', alpha=1, beta=1)

    a = pm.Gamma('a', alpha=alpha, beta=1.0/scale, shape=n_test_participants)

    #mu = T.exp(a[test_participants_idx])
    mu = a[test_participants_idx]

    y = test_df['time_delay'].values
    y_est = pm.Poisson('y_est', mu=mu, observed=y)

    start = pm.find_MAP(fmin=scipy.optimize.fmin_powell)
    step = pm.Metropolis(start=start)
    trace = pm.sample(20000, step, start=start, progressbar=True)

为什么这种分层泊松模型与生成数据中的真实参数不匹配？

Why does this hierarchical Poisson model not match true params from generated data?

python

mcmc

pymc