Thanks a lot for your answer @ricardoV94
Indeed, your fix works, although only with initval bigger than my maximum samples value. I think the reason is a methodological error of mine, a quite silly one. The model I wrote implies that the value of counts is the same for all the observations in samples, which is not true. There should be a hidden count for each sample that is reduced by the action of the BetaBinomial. So there are two stochastic processes at action. One generating counts through my custom distribution of which I want to determine the three parameters, and another stochastic process that for every draw of the previous, generate another value through the BetaBinomial.
I can restate the model as a DensityDist which is a convolution of the BetaBinomial and my custom probability distributions I wrote for counts. However, the new logpfunction is now terribly slow, a single call takes seconds. Is there a way to approach this kind of problem in a smarter way using PyMC? Thanks!