pymc3: model hierarchiczny z wieloma zmiennymi obsesrved

Question

Mam prosty model hierarchiczny z dużą liczbą osób, dla których mam małe próbki z rozkładu normalnego. Środki tych rozkładów również mają rozkład normalny.

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

Chcę użyć PyMC3 do obliczenia parametrów modelu z próbki.

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

Spodziewałem Niedrażniącym model_means wyglądać mojego pierwotnego podziału środków. Wydaje się jednak, że średnia środków jest zbieżna z 30. Jak mogę odzyskać oryginalne odchylenie standardowe środków (12 w moim przykładzie) z modelu pymc3?

Answer 1

To pytanie mnie zmagało się z koncepcjami PyMC3.

Potrzebuję n_individuals obserwowanych zmiennych losowych do modelowania stochastycznych zmiennych losowych y i n_individual w celu modelowania means. Te również wymagają priors hyper_mean i hyper_sigma dla swoich parametrów. sigmas jest poprzednikiem standardowego odchylenia 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()