Hi again!
I tried it out and it sill threw an error regarding broadcasting lambda.
I gave it some thought and, correct me if I am wrong, we only need to index the data on the likelihood. Also I realized that what I actually want is that all articles within a family(the dataset I am using) share a prior of what I call pdiffs, just assume these are coefficients regarding sales.
This way if for a specific article I do not have information on a particular pdiff, it can use the average distribution given the hyperprior. Is that assumption right?
with pm.Model(coords=coords) as allprod_model:
alpha_a = pm.Uniform("alpha_a", lower=0, upper=50)
beta_a = pm.Uniform("beta_a", lower=0, upper=50)
alpha_b = pm.Uniform("alpha_b", lower=0, upper=50)
beta_b = pm.Uniform("beta_b", lower=0, upper=50)
alpha_pdiff = pm.Gamma("alpha_pdiff", alpha=alpha_a, beta=beta_a, dims = 'pdiff')
beta_pdiff = pm.Gamma("beta_pdiff", alpha=alpha_b, beta=beta_b, dims = 'pdiff')
lbda = pm.Gamma('lbda', alpha=alpha_pdiff[None, :], beta=beta_pdiff[None, :], dims=('art', 'pdiff'))
y = pm.Poisson('y', mu=lbda[art_idx,pdiff_idx], observed = data['Uds'])
idata_allarts_hier = pm.sample(draws=4000, tune=4000, nuts_sampler='nutpie')
I compared the results with a model at an article level and in some pdiffs for the same article one can see the shrinkage that should happen when using hierarchical models. Also for each article I have all pdiffs, with can prove useful to estimate sales at a particular pdiff.
I also changed the sampler since it was a bit too slow, it still is unpractical as it takes 90 minutes for 200 articles and I have approximately 30k, but at least is an improvement. I assume that by choosing more informative priors using domain knowledge and maybe changing other things it might improve futher.
Any thoughts/corrections would be highly appreciated.
Kind regards,
Gabriel.