Sorry, spoke to soon. Having defined random, I still have the following error when working with incomplete data. Say for example, I adjust the script as follows:
ll = pm.DensityDist('ll', logp, random=gamma.random, observed = {'d': observations})
When working with incomplete data, I still get the error:
ValueError: Mass matrix contains zeros on the diagonal.
The derivative of RV User.ravel()[2] is zero.