Data aggregation in crowdsensing application

Ok I’m trying to model our intuitions, I set up a hierachical model in which each user send me a measure that has a mean that comes from a general distribution:

# Initialize random number generator
np.random.seed(123)

# True parameter values
true_mu = 38.095
true_sigma = 0.005

# Size of dataset
size = 100

# Simulate some data, we suppose each user makes a measurement
y = true_mu + np.random.randn(size)*true_sigma

with pm.Model() as hierachical_model:
    # Latent mean and sigma, I want to estimate these parameters
    mu = pm.Normal('mu', 0., 10.)
    sigma = pm.HalfNormal('sigma', 1.)

    # Each user has a mean and a standard devitation that I called theta and rho
    theta = pm.Normal('theta', mu, sigma, shape=size)
    rho = sigma = pm.HalfNormal('rho', 1., shape=size)

    like = pm.Normal('like', theta, rho, observed=y)

    trace = pm.sample(500)
    print(trace['mu'].mean(axis=0), trace['sigma'].mean(axis=0))

    pm.traceplot(trace)
    plt.show()

Unfortunately, the execution time is very long and I can not get the desired results (the posterior parameters are wrong). What am I doing wrong? How can I model your suggestions?