pymc3: hierarchisches Modell mit mehreren Variablen

Question

Ich habe ein einfaches hierarchisches Modell mit vielen Individuen, für die ich kleine Stichproben aus einer Normalverteilung habe. Die Mittel dieser Verteilungen folgen ebenfalls einer Normalverteilung.

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

Ich möchte PyMC3 verwenden, um die Modellparameter aus der Probe zu berechnen.

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

Ich erwarte die hintere von model_means wie meine ursprünglichen Verteilung der Mittel zu suchen. Aber es scheint zu konvergieren 30 das Mittel der Mittel. Wie stelle ich die ursprüngliche Standardabweichung der Mittel (12 in meinem Beispiel) vom pymc3 Modell wieder her?

Answer 1

Diese Frage beschäftigte mich mit den Konzepten von PyMC3.

Ich brauche n_individuals Zufallsvariablen beobachtet die y und n_individual stochastischen Zufallsvariablen zu modellieren die means zu modellieren. Diese benötigen auch die Priors hyper_mean und hyper_sigma für ihre Parameter. sigmas ist der Prior für die Standardabweichung von 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()