I would suggest getting rid of the zeros. If you have many in your larger dataset that reduces the number of parameters considerably. You know their posterior is exactly zero anyway.
Otherwise, you may need to implement a better sampler for the Poisson variable. The Metropolis sampler doesn’t even know Poisson can’t have negative values so it will be proposing useless values many times. This is however something I would only suggest as a last resource.
Otherwise, discrete variables are just trickier to sample, I’m afraid you just have to throw compute at it. You can try to use the JAX backend for the PyMC sampler in case you have a GPU. That could be faster.
You can do that by wrapping your call to pm.sample like this:
with pytensor.config.change_flags(mode="JAX"):
pm.sample()
However I am not sure if that will work with multiple chains in parallel, so you may need to revert to a single chain.