pymc3: modello gerarchico con più variabili ossessiate

Question

Ho un semplice modello gerarchico con un sacco di individui per i quali ho piccoli campioni da una distribuzione normale. Anche i mezzi di queste distribuzioni seguono una distribuzione normale.

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

Desidero utilizzare PyMC3 per calcolare i parametri del modello dal campione.

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

mi aspettavo il posteriore della model_means a guardare come la mia distribuzione originale di mezzi. Ma sembra convergere a 30 il mezzo dei mezzi. Come faccio a recuperare la deviazione standard originale dei mezzi (12 nel mio esempio) dal modello pymc3?

Answer 1

Questa domanda mi ha fatto lottare con i concetti di PyMC3.

Ho bisogno di n_individuals variabili casuali osservate per modellare le variabili casuali stocastiche e n_individual per modellare means. Questi hanno anche bisogno dei priori hyper_mean e hyper_sigma per i loro parametri. sigmas è il precedente per la deviazione standard di .

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