Add inverse Gaussian distribution (closes #2788)

[GidonFrischkorn]

Summary

Closes #2788.

This adds the inverse Gaussian (Wald) distribution in the mean/shape parameterization inv_gaussian(y | mu, lambda): inv_gaussian_lpdf, inv_gaussian_cdf, inv_gaussian_lcdf, inv_gaussian_lccdf, and inv_gaussian_rng. All are vectorized over the three arguments with analytic partials via make_partials_propagator, following lognormal_lpdf and normal_lcdf.

The issue suggested wrapping the Boost density. I implemented the distribution natively instead, because the textbook CDF form Phi(z1) + exp(2*lambda/mu) * Phi(-z2) overflows at shape-to-mean ratios that are routine in response-time modeling. The CDF family is computed entirely in log space, with the scaled upper term collapsed analytically so that no opposing large quantities remain, and an internal elementwise log_Phi carries the lower tail past the point where erfc underflows. (std_normal_lcdf cannot be reused because it reduces over its container argument; the helper placement follows von_mises_cdf.hpp and gamma_lccdf.hpp.) The rng uses the Michael, Schucany & Haas (1976) transformation with the smaller root in reciprocal form, since the usual subtractive form loses precision and returns invalid variates at large mu*w/lambda. Both y = 0 and y = inf are accepted as the closure of the support and handled elementwise, with partials defined to be zero there and wherever the log probability has saturated.

Test reference values are computed with mpmath at 60 digits. The implementation is additionally cross-validated against
statmod, with both compared against a high-precision python implementation rather than against each other.

Two questions where I would value feedback:

• Underflow convention. At survivor underflow the lccdf returns -inf with zero partials; normal_lccdf instead returns a signed infinity for the gradient there. I am happy to switch; the shape partial would then need an explicitly assigned sign.
• Placement of log_Phi. The elementwise log-CDF helper could be promoted to stan/math/prim/fun/ as a standalone function. I kept it internal to keep this PR self-contained, but I am happy to split it out here or in a follow-up.

This is my first contribution to Stan Math; I work on response-time models in the R package bmm, where the inverse Gaussian currently lives as a user-defined function.

Tests

• test/prob/inv_gaussian/: the generated distribution fixtures for lpdf / cdf / cdf_log / ccdf_log, covering ordinary-point values and domain errors across every AD and container variant. The fixture reference implementations are written in linear space from erfc, so they are an independent code path from the log-space implementation under test.
• test/unit/math/mix/prob/: expect_ad for all four functions over scalars and containers, plus the vectorized cdf, whose partials couple through the whole-container product; invalid arguments; tail points that cross the internal asymptotic branch and the region where exp(2*lambda/mu) would overflow.
• test/unit/math/rev/prob/: hand-derived partials pinned against mpmath; a check that tape growth does not depend on N; vectorized values and adjoints equal to the scalar sums; zero partials at the support boundaries and at saturation, asserted directly because finite differences step off the support there.
• test/unit/math/prim/prob/: RNG rig with quantile agreement against boost::math::inverse_gaussian, RNG stability at extreme parameters, value pins in the overflow region and the deep lower tail, F + S = 1 consistency, and direct tests of the internal log_Phi.

Side Effects

I am not aware of any; only new functions are added. stan/math/prim/prob.hpp gains the five includes.

Release notes

Added the inverse Gaussian (Wald) distribution: inv_gaussian_lpdf, inv_gaussian_cdf, inv_gaussian_lcdf, inv_gaussian_lccdf, inv_gaussian_rng.

AI use disclosure

I used AI assistance (Claude Code) for this contribution. I am at the start of learning C++ and used it to express established mathematics (the density and CDF identities and the Michael, Schucany & Haas (1976) sampler) in the Stan Math template and autodiff idiom. I verified the mathematics and all test reference values myself, have reviewed and understood the code, and can explain and defend the changes during review.

Checklist

• Copyright holder: Gidon T. Frischkorn

  The copyright holder is typically you or your assignee, such as a university or company. By submitting this pull request, the copyright holder is agreeing to the license the submitted work under the following licenses:
  - Code: BSD 3-clause (https://opensource.org/licenses/BSD-3-Clause)
  - Documentation: CC-BY 4.0 (https://creativecommons.org/licenses/by/4.0/)

• the basic tests are passing

  • unit tests pass (to run, use: ./runTests.py test/unit)
  • header checks pass, (make test-headers)
  • dependencies checks pass, (make test-math-dependencies)
  • docs build, (make doxygen)
  • code passes the built in C++ standards checks (make cpplint)

• the code is written in idiomatic C++ and changes are documented in the doxygen

• the new changes are tested

[WardBrian]

Thanks for your submission @GidonFrischkorn!

We have a release feature freeze starting at the end of this week, so don't be too discouraged if we don't look at this for a week or two.

The tests you currently have are failing on jenkins, seemingly due to minor numeric issues. These are probably to be expected if you were setting the tolerances based off e.g. a Mac, and we are testing on Linux, and assuming the values look fine to you it's acceptable to simply loosen the tolerances from 1e-13 to 1e-11 or similar

[WardBrian]

A few of these tolerances seem a bit too tight on our test platform (linux/clang 7). Errors look like:

test/unit/math/prim/prob/inv_gaussian_test.cpp:191: Failure

The difference between ref_lcdf and a is 8.7538865045644343e-12, which exceeds 1e-13, where
ref_lcdf evaluates to -0.69314591898922784,
a evaluates to -0.69314591898047395

[GidonFrischkorn]

Thanks @WardBrian for this info. The tests were indeed run on a Mac.

I checked through the tests and adapted the test tolerances to be scale with lambda over mu. Details are in the comments for the tests

[GidonFrischkorn]

Thanks for your submission @GidonFrischkorn!

We have a release feature freeze starting at the end of this week, so don't be too discouraged if we don't look at this for a week or two.

The tests you currently have are failing on jenkins, seemingly due to minor numeric issues. These are probably to be expected if you were setting the tolerances based off e.g. a Mac, and we are testing on Linux, and assuming the values look fine to you it's acceptable to simply loosen the tolerances from 1e-13 to 1e-11 or similar

No worries about the feature freeze. The feature is not urgent and especially since this is my first PR, I would rather not rush it and make sure everything is in place and well reviewed prior to merging.