[FR] Alternative parametrisations for numerical stability (e.g. for NegativeBinomial distribution)

Some of the distributions provide multiple different parametrisation options. For example, the MultivariateNormal distribution can be parametrised in terms of covariance, precision, or Cholesky decomposition of the covariance matrix. However, irrespective of the parametrisation used, the distribution is ultimately parametrised using the Cholesky decomposition.

This approach works well in principle but can suffer from numerical issues. For example, consider the NegativeBinomial distribution with concentration a and rate b. In the limit a / b ** 2 -> 0, the NegativeBinomial distribution should reduce to the Poisson distribution with rate a / b. Of course, this limit is problematic in practice because of the limits of floating point arithmetic. Stan has implemented multiple different parametrisations for the negative-binomial distribution (see here for details).

Is there interest in implementing different versions of distributions that could benefit from alternative parametrisations? If yes, should they reuse the interface already used by the MultivariateNormal distribution? Using a similar approach to Stan's neg_binom_2 is also an option but seems clunky.

Note: this may be a discussion for torch rather than pyro to get alignment on interface (if this is at all desirable).