~That looks like the old pymc3 API, I suggest you use PyMC. CustomDist should require much less code: https://www.pymc.io/projects/docs/en/latest/api/distributions/generated/pymc.CustomDist.html~
Nevermind, I got confused. You’re just showing the scipy implementation.
One thing to note is that you won’t be able to use NUTS because you’re using a custom Op. It would be better if you could define the logp using PyTensor operations, in which case you would get gradients for free, but those percentiles may be tricky.