Advice with hierarchical model

Thanks, again. Sorry to have taken so long to respond. Your advice about divergence was very helpful, and between your site and the two articles and the lecture video it pointed me to, I was able to solve my divergence problems.

A quick follow-up – do you know of resources like that for dealing with failure to get sufficient effective samples?

Here’s my case – I did the divergence trick of replacing

sigmaRVs = pm.Gamma(self.sigma[i],
                   mu=hyper_sigmaRVs_mu[i],
                   sd=hyper_sigmaRVs_sigma[i],
                   shape=self.num_replicates)

with this:

sigmaRVs = [pm.Deterministic('%s replicate %d' % (self.sigma[i], j),
                            largest(hyper_sigmaRVs_mu[i] +
                                    (sigma_offset[j] * hyper_sigmaRVs_sigma[i]),
                                    0.3))
                    for j in range(self.num_replicates)] 

The good news is that this solved my divergence problem (well, that and raising my variance priors, which were too low).

The bad news is that the sigmaRVs in the new models are getting hardly any effective samples. The sigma_offset RVs get between 75-150 effective samples in chains of 3000, and the sigmaRVs approximately 2!

These sigmas, and corresponding mus, fit as hyperparameters into a normal RV that is observed.

The mu RVs show more effective samples, but the same relative relationship: their offset variables have two orders of magnitude more effective samples. The offset RVs have about 2200 effective samples in 4 chains of 1000 samples, each, and the mu RVs only 20-25.

Can anyone suggest something I might do to this model to improve things, or do I just need more iterations?

Thanks