title: pymc3: hierarchical model with multiple obsesrved variables
body: <p>I have a simple hierarchical model with lots of individuals for which I have small samples from a normal distribution. The means of these distributions also follow a normal distribution.</p>
  
  <pre class="lang-py prettyprint-override"><code>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))
  </code></pre>
  
  <p>I want to use PyMC3 to compute the model parameters from the sample.</p>
  
  <pre class="lang-py prettyprint-override"><code>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()
  </code></pre>
  
  <p><a href="https://i.stack.imgur.com/9SiKS.png" rel="noreferrer"><img src="https://i.stack.imgur.com/9SiKS.png" alt="mcmc samples"></a></p>
  
  <p>I was expecting the posterior of <code>model_means</code> to look like my original distribution of means. But it seems to converge to <code>30</code> the mean of the means. How do I recover the original standard deviation of the means (12 in my example) from the pymc3 model?</p>