Pymc3: hierarchical model with several observational variables

Question

I have a simple hierarchical model with a lot of people for whom I have small samples from a regular distribution. The means of these distributions also follow the normal distribution.

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

I want to use PyMC3 to calculate model parameters from a sample.

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

I was expecting the back of model_means to look like my original distribution of funds. But he seems to be converging to an average of 30 . How to restore the original standard deviation of funds (12 in my example) from the pymc3 model?

Answer 1

This question was that I was struggling with PyMC3 concepts.

I need n_individuals observable random variables for modeling stochastic random variables y and n_individual for modeling means . They also require hyper_mean and hyper_sigma . sigmas is the forerunner of y standard deviation.

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