Eu tenho um modelo hierárquico simples com muitos indivíduos para os quais tenho pequenas amostras de uma distribuição normal. Os meios dessas distribuições também seguem uma distribuição 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))

Eu quero usar o PyMC3 para calcular os parâmetros do modelo a partir da amostra.

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

Eu estava esperando o posterior demodel_means para parecer com minha distribuição original de meios. Mas parece convergir para30 a média dos meios. Como recupero o desvio padrão original das médias (12 no meu exemplo) do modelo pymc3?