Unfortunately it’s not a truncation, I tried that initially. The data likelihood factorizes in the physical process and the rejection separately. Also, the time series are of varying length and I need to marginalize over all possible lengths. You really need to caclulate p = p(accept) + p(accept)*p(reject) + p(accept)*p(reject)^2 + … to get the correct answer.
If I understand pm.CustomDist correctly, this would then just mean running a simple Markov Chain sampler that’s not using the gradients, right? In that case, I already have that implemented myself. I think PyMC might unfortunately not be the right tool for this.