Hierarchical models with several of coefficients

Are you sure n_stores and n_features are int and not float?
Also, try casting mu and sigma to floats: pm.Normal('b_offset', mu=0.0, sigma=1.0, shape=(n_stores, n_features))