Couple things:
First, do you have a good reason for forcing your background signal to have integer rate? That seems unlikely, particularly if this is just instrument noise. Besides, I believe you’ll in general do better with continuous distributions that discrete ones unless you have no other choice.
I’m a little confused with what you’re doing with p, why is a Binomial coming into play, and why all the 0.5s?
I think you’re also conflating the exponential decay rate with the Poisson rate at each time point. Unless I’m very mistaken, your histogram method is the wrong way to synthesize such data: the distribution of counts over a certain interval is not the same as the number of counts over each of successive intervals (just think, your x-axis units are different).
Further to what @bdyetton said, see this paper, particularly Eq 31 (page 8): https://www.phys.ufl.edu/courses/phy4803L/group_I/gamma_gamma/poisson.pdf
Following that equation, I believe your model should be:
with pm.Model() as model:
Γ = pm.HalfNormal('Γ',1) # Rate of exponential decay for main Poisson rate
λ = pm.HalfNormal('λ',10) # Background Poisson rate
N = pm.HalfNormal('N',1000) # Scale factor / initial Poisson rate
μ = N*Γ*pm.math.exp(-Γ*t)+λ # Observed Poisson (count) rate as a function of time
counts = pm.Poisson('counts', mu=μ, observed = counts_total)
Note that you should use a similar approach to generating the simulated data against which to check your code. Also, that paper includes an ε, the probability of detection, which I don’t include here. If you don’t have very tight a priori estimates or bounds on N and ε, including both in your model would render them un-identifiable and result in very large uncertainties. Putting the 100 multiplier is just so that the Half-Normal can be specified
(P.S., I haven’t actually run that code, so caveat emptor. In particular, those prior means may be too weakly informative, and with means substantially far from 1 you may run into divergences or other sampling issues so you may wish to re-scale them.)