Tengo un modelo jerárquico simple con muchas personas para las que tengo pequeñas muestras de una distribución normal. Los medios de estas distribuciones también siguen una distribución normal.

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))

Quiero usar PyMC3 para calcular los parámetros del modelo de la muestra.

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()

Esperaba el posterior demodel_means para parecerse a mi distribución original de medios. Pero parece converger a30 La media de los medios. ¿Cómo recupero la desviación estándar original de las medias (12 en mi ejemplo) del modelo pymc3?