I don’t think there’s anything wrong with your current method. It won’t be able to compute the gradient of the posterior density, so you can’t use samplers like NUTS, but since you only seem to have 4 parameters, metropolis or smc should probably still work (but make sure to check convergence carefully!). You’ll need many more than 1000 draws though, and also more tuning. Metropolis draws are not nearly as efficient as NUTS draws. You can also help the sampler by rescaling your variables a bit, for instance
lambda_0_raw = pm.Uniform("lambda_0_raw", 0, 1)
lambda_0 = pm.Deterministic("lambda_0", lambda_0_raw * 7000 + 1000)
If you want to use gradient based samplers, I’d currently go for sunode or diffrax, both should work fine. I didn’t have much time to work on sunode for a while now however, and diffrax is for the most part I think better by now. With diffrax, you’ll still have to wrap the solver manually in an op though, (unless someone finishes this PR). For an example, see PYMC Labs | Blog/aaz05zflou9tru7gm475dxwl