pymc3：具有多个观察变量的层次模型答案

【问题标题】：pymc3: hierarchical model with multiple obsesrved variablespymc3：具有多个观察变量的层次模型
【发布时间】：2016-02-13 03:57:14
【问题描述】：

我有一个简单的层次模型，其中包含许多个人，我有来自正态分布的小样本。这些分布的均值也服从正态分布。

import numpy as np

n_individuals = 200
points_per_individual = 10
means = np.random.normal(30, 12, n_individuals)
y = np.random.normal(means, 1, (points_per_individual, n_individuals))

我想使用 PyMC3 从样本中计算模型参数。

import pymc3 as pm
import matplotlib.pyplot as plt

model = pm.Model()
with model:
    model_means = pm.Normal('model_means', mu=35, sd=15)

    y_obs = pm.Normal('y_obs', mu=model_means, sd=1, shape=n_individuals, observed=y)

    trace = pm.sample(1000)

pm.traceplot(trace[100:], vars=['model_means'])
plt.show()

我期待model_means 的后验看起来像我原来的均值分布。但它似乎收敛到均值的30。如何从 pymc3 模型中恢复均值的原始标准差（在我的示例中为 12）？

【问题讨论】：

标签： bayesian mcmc pymc3

【解决方案1】：

这个问题让我在 PyMC3 的概念上苦苦挣扎。

我需要n_individuals 观察随机变量来建模y 和n_individual 随机随机变量来建模means。这些还需要先验 hyper_mean 和 hyper_sigma 的参数。 sigmas 是 y 标准差的先验。

import matplotlib.pyplot as plt

model = pm.Model()
with model:
    hyper_mean = pm.Normal('hyper_mean', mu=0, sd=100)
    hyper_sigma = pm.HalfNormal('hyper_sigma', sd=3)

    means = pm.Normal('means', mu=hyper_mean, sd=hyper_sigma, shape=n_individuals)
    sigmas = pm.HalfNormal('sigmas', sd=100)

    y = pm.Normal('y', mu=means, sd=sigmas, observed=y)

    trace = pm.sample(10000)

pm.traceplot(trace[100:], vars=['hyper_mean', 'hyper_sigma', 'means', 'sigmas'])
plt.show()

【讨论】：