Add inverse Gaussian distribution support

[alizma]

Description

Analogous to other distributions, add support for the inverse Gaussian distribution with a sampling statement and Stan functions for working with its log density, CDF, and CCDF.

This can be a wrapper around the Boost implementation of the distribution.

Example

The motivation is from writing a GLM with an inverse Gaussian. Similar to rstanarm, I implemented a custom helper function to perform the likelihood computation and it seems like other people in the past have re-implemented the same set of algorithms. For example, see here and here. It would be nice to have this as a library function.

Current Version:

v4.4.0

[bob-carpenter]

[Edit: this is wrong; see next message]

To code this in Stan, you need to follow the approach taken for lognormal or inverse gamma. We don't directly wrap Boost implementations because we want efficient vectorization and derivatives.

As a workaround, if you want to give y an inverse normal distribution in Stan, you can do this:

inv(y) ~ normal(mu, sigma);
target += -2 * log(y);

The Jacobian adjustment is the log of the absolute value of the derivative of the inverse transform, which is just inverse: inv(y) = 1/y. Because log and inv are vectorized and target += accepts containers, the above code even works in cases where y is a vector or 1D array.

[GidonFrischkorn]

I wanted to ask, if the native support for the inverse Gaussian as suggested by @alizma is still in scope. I agree that it would be great to have the inverse Gaussian natively implemented in Stan.

Unfortunately, I do not have experience with coding in C or implementing custom distributions in Stan. But still if I can be of help to get the inverse Gaussian implemented in Stan, I am happy to contribute.

[bob-carpenter]

This would be easy to implement for specific cases, at least without vectorization of all the arguments.

functions {
  real inverse_normal_lpdf(real y, real mu, real lambda) {
    return 0.5 * log(lambda) + -1.5 * log(y) - lambda * 0.5 * ((y - mu) / mu)^2 /  y);
  }
}

I took the definition of the pdf from the following site (no idea why it's not on the regular wikipedia), so I'd double check the pdf and my conversion to the log scale:

https://www.wikiwand.com/en/Inverse_Gaussian_distribution

For a given application, any of y, mu, or lambda could be turned into a vector or array with a more efficient implementation than just looping, but just the above shouldn't be too bad

[GidonFrischkorn]

Thanks for the quick response and for sharing the inverse_normal_lpdf function. There actually are some implementations similar to the one you suggested that are already used in R packages interfacing with Stan, like here: https://github.com/paul-buerkner/brms/blob/master/inst/chunks/fun_inv_gaussian.stan

For now, I have used this approach and defined the function myself as you suggested. But the question is if it would be possible to integrate the Inverse Gaussian as a native distribution in Stan, so that you do not have to define it in the functions block and would potentially also have optimized sampling that is faster than for the user defined function.

[bob-carpenter]

the question is if it would be possible to integrate the Inverse Gaussian as a native distribution in Stan

Yes. It's just a matter of finding a C++ developer willing to add it.

defined the function myself as you suggested.

I'd make sure to check it twice against another implementation. It's easy to get signs or numerator/denominator or log/exp wrong.

optimized sampling that is faster than for the user defined function.

You can also write the function in terms of vectors when there are vectors for y, mu, or lambda, which will make it a lot faster. It' be nice if you could reduce it to one of our existing distributions with a little algebra, but I don't see a way to do that.