ENH: Beta negative binomial distribution

[ferrine]

Before

There is no such a distribution that represents Beta Negative Binomial distribution.

The beta negative binomial is non-identifiable which can be seen easily by simply swapping $r$ and
$\beta$ in the above density or characteristic function and noting that it is unchanged. Thus estimation demands that a constraint be placed on $r$,
$\beta$ or both.

After

with pm.Model() as model:
    pm.BetaNegativeBinomial("bnb", alpha=alpha, beta=1, r=r)

Context for the issue:

No response

[ricardoV94]

I would add it to pymc-experimental. Otherwise sounds good

[williambdean]

There seems to be pr on scipy.stats (will link)
But the sampling step might not be the difficult point and doesn't need to be bottleneck
I can give it a try without waiting for scipy