From 284747e302aa779dfe509b1dadf2aef5e3eb812b Mon Sep 17 00:00:00 2001 From: Gidon Frischkorn Date: Sat, 22 Aug 2026 10:14:02 +0200 Subject: [PATCH 01/11] feat: add inv_gaussian_lpdf and inv_gaussian_rng - mean/shape parameterization, vectorized over all three arguments - rng uses the Michael, Schucany & Haas (1976) transform --- stan/math/prim/prob.hpp | 2 + stan/math/prim/prob/inv_gaussian_lpdf.hpp | 138 ++++++++++++++++++++++ stan/math/prim/prob/inv_gaussian_rng.hpp | 92 +++++++++++++++ 3 files changed, 232 insertions(+) create mode 100644 stan/math/prim/prob/inv_gaussian_lpdf.hpp create mode 100644 stan/math/prim/prob/inv_gaussian_rng.hpp diff --git a/stan/math/prim/prob.hpp b/stan/math/prim/prob.hpp index 4fae64c1ace..95399fb3e0f 100644 --- a/stan/math/prim/prob.hpp +++ b/stan/math/prim/prob.hpp @@ -139,6 +139,8 @@ #include #include #include +#include +#include #include #include #include diff --git a/stan/math/prim/prob/inv_gaussian_lpdf.hpp b/stan/math/prim/prob/inv_gaussian_lpdf.hpp new file mode 100644 index 00000000000..10c9bffcf65 --- /dev/null +++ b/stan/math/prim/prob/inv_gaussian_lpdf.hpp @@ -0,0 +1,138 @@ +#ifndef STAN_MATH_PRIM_PROB_INV_GAUSSIAN_LPDF_HPP +#define STAN_MATH_PRIM_PROB_INV_GAUSSIAN_LPDF_HPP + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace stan { +namespace math { + +/** \ingroup prob_dists + * The log of the inverse Gaussian density for the specified scalar(s) given + * the specified mean(s) and shape(s). y, mu, or lambda can each be either a + * scalar or a vector. Any vector inputs must be the same length. + * + *

The result log probability is defined to be the sum of the log + * probabilities for each observation/mean/shape triple. + * + *

The density is parameterized by the mean \f$\mu\f$ and the shape + * \f$\lambda\f$ as \f$\sqrt{\lambda / (2\pi y^3)} + * \exp(-\lambda (y - \mu)^2 / (2 \mu^2 y))\f$, with support + * \f$y \in (0, \infty)\f$ and variance \f$\mu^3 / \lambda\f$. + * + *

\f$y = 0\f$ is accepted and returns \f$-\infty\f$; a negative + * \f$y\f$ throws. + * + * @tparam T_y type of scalar + * @tparam T_loc type of mean parameter + * @tparam T_shape type of shape parameter + * @param y (Sequence of) scalar(s). + * @param mu (Sequence of) mean parameter(s) for the inverse Gaussian + * distribution. + * @param lambda (Sequence of) shape parameter(s) for the inverse Gaussian + * distribution. + * @return The log of the product of the densities. + * @throw std::domain_error if the mean or the shape is not positive and + * finite, if the random variable is negative, or if any argument is NaN. + * @throw std::invalid_argument if container sizes mismatch. + */ +template * = nullptr> +inline return_type_t inv_gaussian_lpdf( + const T_y& y, const T_loc& mu, const T_shape& lambda) { + using T_partials_return = partials_return_t; + using T_y_ref = ref_type_if_not_constant_t; + using T_mu_ref = ref_type_if_not_constant_t; + using T_lambda_ref = ref_type_if_not_constant_t; + static constexpr const char* function = "inv_gaussian_lpdf"; + check_consistent_sizes(function, "Random variable", y, "Mean parameter", mu, + "Shape parameter", lambda); + + T_y_ref y_ref = y; + T_mu_ref mu_ref = mu; + T_lambda_ref lambda_ref = lambda; + + decltype(auto) y_val = to_ref(as_value_column_array_or_scalar(y_ref)); + decltype(auto) mu_val = to_ref(as_value_column_array_or_scalar(mu_ref)); + decltype(auto) lambda_val + = to_ref(as_value_column_array_or_scalar(lambda_ref)); + + check_nonnegative(function, "Random variable", y_val); + check_positive_finite(function, "Mean parameter", mu_val); + check_positive_finite(function, "Shape parameter", lambda_val); + + if (size_zero(y, mu, lambda)) { + return 0; + } + if constexpr (!include_summand::value) { + return 0; + } + + auto ops_partials = make_partials_propagator(y_ref, mu_ref, lambda_ref); + + // y == 0 has zero density and LOG_ZERO absorbs the sum. + // The gradients are technically ill-defined, but treated as zero. + if (sum(promote_scalar(y_val == 0))) { + return ops_partials.build(LOG_ZERO); + } + + constexpr bool any_ad = is_any_autodiff_v; + + const auto& inv_mu = to_ref_if(inv(mu_val)); + const auto& inv_y = to_ref_if(inv(y_val)); + const auto& y_m_mu = to_ref_if(y_val - mu_val); + // (y - mu)^2 / (2 mu^2 y) + const auto& half_sq_scaled = to_ref(0.5 * square(y_m_mu * inv_mu) * inv_y); + + size_t N = max_size(y, mu, lambda); + T_partials_return logp = -sum(lambda_val * half_sq_scaled); + if constexpr (include_summand::value) { + logp += N * NEG_LOG_SQRT_TWO_PI; + } + if constexpr (include_summand::value) { + logp += 0.5 * sum(log(lambda_val)) * N / math::size(lambda); + } + if constexpr (include_summand::value) { + logp -= 1.5 * sum(log(y_val)) * N / math::size(y); + } + + if constexpr (any_ad) { + if constexpr (is_autodiff_v) { + partials<0>(ops_partials) + = -1.5 * inv_y + 0.5 * lambda_val * (square(inv_y) - square(inv_mu)); + } + if constexpr (is_autodiff_v) { + partials<1>(ops_partials) = lambda_val * y_m_mu * inv_mu * square(inv_mu); + } + if constexpr (is_autodiff_v) { + partials<2>(ops_partials) = 0.5 * inv(lambda_val) - half_sq_scaled; + } + } + return ops_partials.build(logp); +} + +template +inline return_type_t inv_gaussian_lpdf( + const T_y& y, const T_loc& mu, const T_shape& lambda) { + return inv_gaussian_lpdf(y, mu, lambda); +} + +} // namespace math +} // namespace stan +#endif diff --git a/stan/math/prim/prob/inv_gaussian_rng.hpp b/stan/math/prim/prob/inv_gaussian_rng.hpp new file mode 100644 index 00000000000..6fefb3aff0a --- /dev/null +++ b/stan/math/prim/prob/inv_gaussian_rng.hpp @@ -0,0 +1,92 @@ +#ifndef STAN_MATH_PRIM_PROB_INV_GAUSSIAN_RNG_HPP +#define STAN_MATH_PRIM_PROB_INV_GAUSSIAN_RNG_HPP + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace stan { +namespace math { + +/** \ingroup prob_dists + * Return an inverse Gaussian random variate for the given mean and shape + * using the specified random number generator. + * + * mu and lambda can each be a scalar or a one-dimensional container. Any + * non-scalar inputs must be the same size. + * + *

The algorithm used in inv_gaussian_rng is the transformation in: + * + * Generating Random Variates Using Transformations with Multiple Roots + * J. R. Michael, W. R. Schucany and R. W. Haas + * The American Statistician, Vol. 30, No. 2 (1976), pp. 88-90 + * + *

A chi-square variate \f$w\f$ on one degree of freedom gives a pair of + * candidate roots and a Bernoulli draw selects between them. The method is + * exact, with no rejection step. + * + *

The smaller root is formed from the discriminant + * \f$\sqrt{u + u^2/4}\f$ with \f$u = \mu w / \lambda\f$, as + * \f$x = \mu(1 + u/2 - \sqrt{u + u^2/4})\f$. + * + * @tparam T_loc type of mean parameter + * @tparam T_shape type of shape parameter + * @tparam RNG type of random number generator + * @param mu (Sequence of) mean parameter(s) + * @param lambda (Sequence of) shape parameter(s) + * @param rng random number generator + * @return (Sequence of) inverse Gaussian random variate(s) + * @throw std::domain_error if mu or lambda is not positive and finite + * @throw std::invalid_argument if non-scalar arguments are of different + * sizes + */ +template +inline typename VectorBuilder::type +inv_gaussian_rng(const T_loc& mu, const T_shape& lambda, RNG& rng) { + using boost::variate_generator; + using boost::random::normal_distribution; + using boost::random::uniform_01; + using T_mu_ref = ref_type_t; + using T_lambda_ref = ref_type_t; + static constexpr const char* function = "inv_gaussian_rng"; + check_consistent_sizes(function, "Mean parameter", mu, "Shape parameter", + lambda); + T_mu_ref mu_ref = mu; + T_lambda_ref lambda_ref = lambda; + check_positive_finite(function, "Mean parameter", mu_ref); + check_positive_finite(function, "Shape Parameter", lambda_ref); + + scalar_seq_view mu_vec(mu_ref); + scalar_seq_view lambda_vec(lambda_ref); + size_t N = max_size(mu, lambda); + VectorBuilder output(N); + + variate_generator > norm_rng( + rng, normal_distribution<>(0, 1)); + variate_generator > uniform01_rng(rng, uniform_01<>()); + + for (size_t n = 0; n < N; ++n) { + const double mu_dbl = mu_vec[n]; + const double lambda_dbl = lambda_vec[n]; + const double nu = norm_rng(); + const double w = nu * nu; + const double u = mu_dbl * w / lambda_dbl; + const double half_ratio = 0.5 * u; + const double disc = std::sqrt(u + half_ratio * half_ratio); + const double x = mu_dbl * (1.0 + half_ratio - disc); + output[n] + = (uniform01_rng() <= mu_dbl / (mu_dbl + x)) ? x : mu_dbl * mu_dbl / x; + } + + return output.data(); +} + +} // namespace math +} // namespace stan +#endif From 162b63e4dea36d24c0d52c20efaf60c2ad5d5c82 Mon Sep 17 00:00:00 2001 From: Gidon Frischkorn Date: Sat, 22 Aug 2026 15:41:19 +0200 Subject: [PATCH 02/11] feat: add inv_gaussian_cdf, inv_gaussian_lcdf and inv_gaussian_lccdf - compute the CDF family in log space; exp(2*lambda/mu) overflows a double above 2*lambda/mu = 710 --- stan/math/prim/prob.hpp | 3 + stan/math/prim/prob/inv_gaussian_cdf.hpp | 130 +++++++++++ stan/math/prim/prob/inv_gaussian_lccdf.hpp | 173 +++++++++++++++ stan/math/prim/prob/inv_gaussian_lcdf.hpp | 240 +++++++++++++++++++++ 4 files changed, 546 insertions(+) create mode 100644 stan/math/prim/prob/inv_gaussian_cdf.hpp create mode 100644 stan/math/prim/prob/inv_gaussian_lccdf.hpp create mode 100644 stan/math/prim/prob/inv_gaussian_lcdf.hpp diff --git a/stan/math/prim/prob.hpp b/stan/math/prim/prob.hpp index 95399fb3e0f..c66b4715187 100644 --- a/stan/math/prim/prob.hpp +++ b/stan/math/prim/prob.hpp @@ -139,6 +139,9 @@ #include #include #include +#include +#include +#include #include #include #include diff --git a/stan/math/prim/prob/inv_gaussian_cdf.hpp b/stan/math/prim/prob/inv_gaussian_cdf.hpp new file mode 100644 index 00000000000..97973700897 --- /dev/null +++ b/stan/math/prim/prob/inv_gaussian_cdf.hpp @@ -0,0 +1,130 @@ +#ifndef STAN_MATH_PRIM_PROB_INV_GAUSSIAN_CDF_HPP +#define STAN_MATH_PRIM_PROB_INV_GAUSSIAN_CDF_HPP + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace stan { +namespace math { + +/** \ingroup prob_dists + * Returns the inverse Gaussian cumulative distribution function for the + * given random variable, mean and shape. Given containers of matching sizes, + * returns the product of probabilities. + * + *

The CDF is computed in log space and exponentiated once at the end, + * since the factor \f$e^{2\lambda/\mu}\f$ in + * \f$\Phi(z_1) + e^{2\lambda/\mu}\Phi(-z_2)\f$ overflows a double above + * \f$2\lambda/\mu = 710\f$. See inv_gaussian_lcdf. + * + *

A boundary observation short-circuits the whole return: + * \f$y = 0\f$ and \f$y = \infty\f$ each take the result to their + * limit. The partials are zero at both. + * + * @tparam T_y type of scalar + * @tparam T_loc type of mean parameter + * @tparam T_shape type of shape parameter + * @param y scalar or container of scalars + * @param mu mean parameter + * @param lambda shape parameter + * @return the product of probabilities + * @throw std::domain_error if the mean or the shape is not positive and + * finite, if the random variable is negative, or if any argument is NaN. + * @throw std::invalid_argument if container sizes mismatch. + */ +template * = nullptr> +inline return_type_t inv_gaussian_cdf( + const T_y& y, const T_loc& mu, const T_shape& lambda) { + using T_partials_return = partials_return_t; + using T_y_ref = ref_type_if_not_constant_t; + using T_mu_ref = ref_type_if_not_constant_t; + using T_lambda_ref = ref_type_if_not_constant_t; + static constexpr const char* function = "inv_gaussian_cdf"; + check_consistent_sizes(function, "Random variable", y, "Mean parameter", mu, + "Shape parameter", lambda); + + T_y_ref y_ref = y; + T_mu_ref mu_ref = mu; + T_lambda_ref lambda_ref = lambda; + + decltype(auto) y_val = to_ref(as_value_column_array_or_scalar(y_ref)); + decltype(auto) mu_val = to_ref(as_value_column_array_or_scalar(mu_ref)); + decltype(auto) lambda_val + = to_ref(as_value_column_array_or_scalar(lambda_ref)); + + check_nonnegative(function, "Random variable", y_val); + check_positive_finite(function, "Mean parameter", mu_val); + check_positive_finite(function, "Shape parameter", lambda_val); + + if (size_zero(y, mu, lambda)) { + return 1.0; + } + + auto ops_partials = make_partials_propagator(y_ref, mu_ref, lambda_ref); + + if (sum(promote_scalar(y_val == 0))) { + return ops_partials.build(0.0); + } + if (sum(promote_scalar(y_val == INFTY))) { + return ops_partials.build(1.0); + } + + constexpr bool any_ad = is_any_autodiff_v; + + const auto& inv_mu = to_ref_if(inv(mu_val)); + const auto& inv_y = to_ref_if(inv(y_val)); + const auto& sqrt_lambda_over_y = to_ref_if(sqrt(lambda_val * inv_y)); + const auto& y_over_mu = to_ref_if(y_val * inv_mu); + const auto& z1 = to_ref(sqrt_lambda_over_y * (y_over_mu - 1.0)); + const auto& z2 = to_ref(sqrt_lambda_over_y * (y_over_mu + 1.0)); + + const auto& log_upper = to_ref(internal::log_scaled_upper_term(z1, z2)); + const auto& lcdf_elt = to_ref(log_sum_exp(internal::log_Phi(z1), log_upper)); + + T_partials_return cdf = exp(sum(lcdf_elt)); + + if constexpr (any_ad) { + // d cdf / d theta_n = cdf * d log cdf_n / d theta_n + // + // Each partial is masked separately: at y == inf the shape partial is + // 0 / 0 on its own. The inner select covers elements whose log CDF has + // saturated to -inf at an interior y. + const auto& w_dens = to_ref(cdf * exp(internal::log_phi(z1) - lcdf_elt)); + const auto& w_upper = to_ref(cdf * exp(log_upper - lcdf_elt)); + if constexpr (is_autodiff_v) { + partials<0>(ops_partials) = w_dens * sqrt_lambda_over_y * inv_y; + } + if constexpr (is_autodiff_v) { + partials<1>(ops_partials) = -2.0 * lambda_val * w_upper * square(inv_mu); + } + if constexpr (is_autodiff_v) { + partials<2>(ops_partials) + = 2.0 * w_upper * inv_mu - w_dens * inv_y / sqrt_lambda_over_y; + } + } + return ops_partials.build(cdf); +} + +} // namespace math +} // namespace stan +#endif diff --git a/stan/math/prim/prob/inv_gaussian_lccdf.hpp b/stan/math/prim/prob/inv_gaussian_lccdf.hpp new file mode 100644 index 00000000000..2f19f1d32aa --- /dev/null +++ b/stan/math/prim/prob/inv_gaussian_lccdf.hpp @@ -0,0 +1,173 @@ +#ifndef STAN_MATH_PRIM_PROB_INV_GAUSSIAN_LCCDF_HPP +#define STAN_MATH_PRIM_PROB_INV_GAUSSIAN_LCCDF_HPP + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace stan { +namespace math { +namespace internal { + +/** + * Return log_diff_exp(a, b), or negative infinity when b >= a. The survivor + * is a difference of two terms that meet in the far upper tail, where + * rounding can invert their order and zero is the correct answer. + */ +template * = nullptr> +inline return_type_t log_diff_exp_guarded(const T1& a, const T2& b) { + if (unlikely(value_of_rec(b) >= value_of_rec(a))) { + return NEGATIVE_INFTY; + } + return log_diff_exp(a, b); +} + +/** + * A vectorized version of log_diff_exp_guarded(). + */ +template * = nullptr> +inline auto log_diff_exp_guarded(T1&& a, T2&& b) { + return apply_scalar_binary( + [](auto&& c, auto&& d) { + return log_diff_exp_guarded(std::forward(c), + std::forward(d)); + }, + std::forward(a), std::forward(b)); +} + +} // namespace internal + +/** \ingroup prob_dists + * Returns the inverse Gaussian log complementary cumulative distribution + * function for the given random variable, mean and shape. Given containers + * of matching sizes, returns the log of the product of complementary + * probabilities. + * + *

The survivor function is \f$S(y \mid \mu, \lambda) = \Phi(-z_1) + * - e^{2\lambda/\mu}\Phi(-z_2)\f$ with + * \f$z_{1,2} = \sqrt{\lambda/y}\,(y/\mu \mp 1)\f$. Both terms are evaluated + * as lower-tail log-CDFs, which stay accurate far into the tail, and the + * difference is taken in log space. See inv_gaussian_lcdf. + * + *

Deep in the upper tail the two terms meet within the spacing of a + * double and the guarded difference returns \f$-\infty\f$. Over + * \f$\mu \in [10^{-3}, 10^3]\f$, \f$\lambda \in [10^{-3}, 10^{16}]\f$ + * that first happens at \f$\log S = -5 \times 10^5\f$, so the survivor is + * zero to any representable precision there. + * + *

A boundary observation short-circuits the whole return: + * \f$y = 0\f$ and \f$y = \infty\f$ each take the result to their + * limit. The partials are zero at both. + * + * @tparam T_y type of scalar + * @tparam T_loc type of mean parameter + * @tparam T_shape type of shape parameter + * @param y scalar or container of scalars + * @param mu mean parameter + * @param lambda shape parameter + * @return the log of the product of complementary probabilities + * @throw std::domain_error if the mean or the shape is not positive and + * finite, if the random variable is negative, or if any argument is NaN. + * @throw std::invalid_argument if container sizes mismatch. + */ +template * = nullptr> +inline return_type_t inv_gaussian_lccdf( + const T_y& y, const T_loc& mu, const T_shape& lambda) { + using T_partials_return = partials_return_t; + using T_y_ref = ref_type_if_not_constant_t; + using T_mu_ref = ref_type_if_not_constant_t; + using T_lambda_ref = ref_type_if_not_constant_t; + static constexpr const char* function = "inv_gaussian_lccdf"; + check_consistent_sizes(function, "Random variable", y, "Mean parameter", mu, + "Shape parameter", lambda); + + T_y_ref y_ref = y; + T_mu_ref mu_ref = mu; + T_lambda_ref lambda_ref = lambda; + + decltype(auto) y_val = to_ref(as_value_column_array_or_scalar(y_ref)); + decltype(auto) mu_val = to_ref(as_value_column_array_or_scalar(mu_ref)); + decltype(auto) lambda_val + = to_ref(as_value_column_array_or_scalar(lambda_ref)); + + check_nonnegative(function, "Random variable", y_val); + check_positive_finite(function, "Mean parameter", mu_val); + check_positive_finite(function, "Shape parameter", lambda_val); + + if (size_zero(y, mu, lambda)) { + return 0; + } + + auto ops_partials = make_partials_propagator(y_ref, mu_ref, lambda_ref); + + if (sum(promote_scalar(y_val == 0))) { + return ops_partials.build(0); + } + if (sum(promote_scalar(y_val == INFTY))) { + return ops_partials.build(NEGATIVE_INFTY); + } + + constexpr bool any_ad = is_any_autodiff_v; + + const auto& inv_mu = to_ref_if(inv(mu_val)); + const auto& inv_y = to_ref_if(inv(y_val)); + const auto& sqrt_lambda_over_y = to_ref_if(sqrt(lambda_val * inv_y)); + const auto& y_over_mu = to_ref_if(y_val * inv_mu); + const auto& z1 = to_ref(sqrt_lambda_over_y * (y_over_mu - 1.0)); + const auto& z2 = to_ref(sqrt_lambda_over_y * (y_over_mu + 1.0)); + + const auto& log_upper = to_ref(internal::log_scaled_upper_term(z1, z2)); + const auto& lccdf_elt = to_ref( + internal::log_diff_exp_guarded(internal::log_Phi(-z1), log_upper)); + + T_partials_return lccdf = sum(lccdf_elt); + + if constexpr (any_ad) { + // Same phi collapse as in the lcdf; the partials are its sign flip, with + // the survivor in place of the CDF. + // + // Each partial is masked separately: at y == inf the shape partial is + // 0 / 0 on its own. The inner select covers elements whose survivor has + // underflowed to -inf. + const auto& w_dens = to_ref(exp(internal::log_phi(z1) - lccdf_elt)); + const auto& w_upper = to_ref(exp(log_upper - lccdf_elt)); + if constexpr (is_autodiff_v) { + partials<0>(ops_partials) = -w_dens * sqrt_lambda_over_y * inv_y; + } + if constexpr (is_autodiff_v) { + partials<1>(ops_partials) = 2.0 * lambda_val * w_upper * square(inv_mu); + } + if constexpr (is_autodiff_v) { + partials<2>(ops_partials) + = w_dens * inv_y / sqrt_lambda_over_y - 2.0 * w_upper * inv_mu; + } + } + return ops_partials.build(lccdf); +} + +} // namespace math +} // namespace stan +#endif diff --git a/stan/math/prim/prob/inv_gaussian_lcdf.hpp b/stan/math/prim/prob/inv_gaussian_lcdf.hpp new file mode 100644 index 00000000000..071ded53012 --- /dev/null +++ b/stan/math/prim/prob/inv_gaussian_lcdf.hpp @@ -0,0 +1,240 @@ +#ifndef STAN_MATH_PRIM_PROB_INV_GAUSSIAN_LCDF_HPP +#define STAN_MATH_PRIM_PROB_INV_GAUSSIAN_LCDF_HPP + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace stan { +namespace math { +namespace internal { + +// erfc underflows to zero below z = -38.5; the two branches agree to a +// relative 3e-15 at this cutoff. +static constexpr double LOG_PHI_ASYMPTOTIC_CUTOFF = -30.0; + +/** + * Return the log of the standard normal cumulative distribution function. At + * or above LOG_PHI_ASYMPTOTIC_CUTOFF this is log(0.5 * erfc(-z / sqrt(2))); + * below it the asymptotic expansion + * \f$-z^2/2 - \log(-z) - \log(2\pi)/2 + * + \log(1 - s + 3s^2 - 15s^3 + 105s^4)\f$, \f$s = z^{-2}\f$, carries the + * lower tail to a relative 1e-15 down to z = -1000. + */ +template * = nullptr> +inline return_type_t log_Phi(const T& z) { + using std::log; + if (value_of_rec(z) >= LOG_PHI_ASYMPTOTIC_CUTOFF) { + return LOG_HALF + log(erfc(-z * INV_SQRT_TWO)); + } + const auto s = inv_square(z); + return -0.5 * square(z) - log(-z) - HALF_LOG_TWO_PI + + log1p(s * (-1.0 + s * (3.0 + s * (-15.0 + s * 105.0)))); +} + +struct log_Phi_fun { + template + static inline auto fun(T&& z) { + return log_Phi(std::forward(z)); + } +}; + +/** + * A vectorized version of log_Phi(). + */ +template * = nullptr> +inline auto log_Phi(T&& z) { + return apply_scalar_unary::apply(std::forward(z)); +} + +/** + * Return the log of the standard normal density at the specified value. + */ +template * = nullptr> +inline return_type_t log_phi(const T& z) { + return -0.5 * square(z) - HALF_LOG_TWO_PI; +} + +struct log_phi_fun { + template + static inline auto fun(T&& z) { + return log_phi(std::forward(z)); + } +}; + +/** + * A vectorized version of log_phi(). + */ +template * = nullptr> +inline auto log_phi(T&& z) { + return apply_scalar_unary::apply(std::forward(z)); +} + +/** + * Return the log of the second inverse Gaussian CDF term, + * \f$\log\left(e^{2\lambda/\mu}\,\Phi(-z_2)\right)\f$. Since + * \f$z_2^2 - z_1^2 = 4\lambda/\mu\f$, the asymptotic expansion of + * \f$\log\Phi(-z_2)\f$ cancels \f$2\lambda/\mu\f$ against \f$-z_2^2/2\f$ and + * leaves \f$\log\phi(z_1) - \log z_2 + \log(1 - s + 3s^2 - 15s^3 + 105s^4)\f$ + * with \f$s = z_2^{-2}\f$, in which no two terms are large and opposed. Below + * the asymptotic regime \f$z_2 \le 30\f$ bounds \f$2\lambda/\mu\f$ by 450 and + * the direct form is accurate. + */ +template * = nullptr> +inline return_type_t log_scaled_upper_term(const T1& z1, const T2& z2) { + using std::log; + if (value_of_rec(z2) > -LOG_PHI_ASYMPTOTIC_CUTOFF) { + const auto s = inv_square(z2); + return log_phi(z1) - log(z2) + + log1p(s * (-1.0 + s * (3.0 + s * (-15.0 + s * 105.0)))); + } + return 0.5 * (square(z2) - square(z1)) + log_Phi(-z2); +} + +/** + * A vectorized version of log_scaled_upper_term(). + */ +template * = nullptr> +inline auto log_scaled_upper_term(T1&& z1, T2&& z2) { + return apply_scalar_binary( + [](auto&& a, auto&& b) { + return log_scaled_upper_term(std::forward(a), + std::forward(b)); + }, + std::forward(z1), std::forward(z2)); +} + +} // namespace internal + +/** \ingroup prob_dists + * Returns the inverse Gaussian log cumulative distribution function for the + * given random variable, mean and shape. Given containers of matching sizes, + * returns the log of the product of probabilities. + * + *

The CDF is \f$F(y \mid \mu, \lambda) = \Phi(z_1) + * + e^{2\lambda/\mu}\Phi(-z_2)\f$ with + * \f$z_{1,2} = \sqrt{\lambda/y}\,(y/\mu \mp 1)\f$. The factor + * \f$e^{2\lambda/\mu}\f$ overflows a double above \f$2\lambda/\mu = 710\f$, + * so the two terms are combined in log space. That grouping of \f$z\f$ is + * exact at \f$y = \mu\f$ and returns \f$z_1 = 0\f$ there for any + * \f$\lambda\f$. + * + *

A boundary observation short-circuits the whole return: + * \f$y = 0\f$ and \f$y = \infty\f$ each take the result to their + * limit. The partials are zero at both. + * + * @tparam T_y type of scalar + * @tparam T_loc type of mean parameter + * @tparam T_shape type of shape parameter + * @param y scalar or container of scalars + * @param mu mean parameter + * @param lambda shape parameter + * @return the log of the product of probabilities + * @throw std::domain_error if the mean or the shape is not positive and + * finite, if the random variable is negative, or if any argument is NaN. + * @throw std::invalid_argument if container sizes mismatch. + */ +template * = nullptr> +inline return_type_t inv_gaussian_lcdf( + const T_y& y, const T_loc& mu, const T_shape& lambda) { + using T_partials_return = partials_return_t; + using T_y_ref = ref_type_if_not_constant_t; + using T_mu_ref = ref_type_if_not_constant_t; + using T_lambda_ref = ref_type_if_not_constant_t; + static constexpr const char* function = "inv_gaussian_lcdf"; + check_consistent_sizes(function, "Random variable", y, "Mean parameter", mu, + "Shape parameter", lambda); + + T_y_ref y_ref = y; + T_mu_ref mu_ref = mu; + T_lambda_ref lambda_ref = lambda; + + decltype(auto) y_val = to_ref(as_value_column_array_or_scalar(y_ref)); + decltype(auto) mu_val = to_ref(as_value_column_array_or_scalar(mu_ref)); + decltype(auto) lambda_val + = to_ref(as_value_column_array_or_scalar(lambda_ref)); + + check_nonnegative(function, "Random variable", y_val); + check_positive_finite(function, "Mean parameter", mu_val); + check_positive_finite(function, "Shape parameter", lambda_val); + + if (size_zero(y, mu, lambda)) { + return 0; + } + + auto ops_partials = make_partials_propagator(y_ref, mu_ref, lambda_ref); + + if (sum(promote_scalar(y_val == 0))) { + return ops_partials.build(NEGATIVE_INFTY); + } + if (sum(promote_scalar(y_val == INFTY))) { + return ops_partials.build(0); + } + + constexpr bool any_ad = is_any_autodiff_v; + + const auto& inv_mu = to_ref_if(inv(mu_val)); + const auto& inv_y = to_ref_if(inv(y_val)); + const auto& sqrt_lambda_over_y = to_ref_if(sqrt(lambda_val * inv_y)); + const auto& y_over_mu = to_ref_if(y_val * inv_mu); + const auto& z1 = to_ref(sqrt_lambda_over_y * (y_over_mu - 1.0)); + const auto& z2 = to_ref(sqrt_lambda_over_y * (y_over_mu + 1.0)); + + const auto& log_upper = to_ref(internal::log_scaled_upper_term(z1, z2)); + const auto& lcdf_elt = to_ref(log_sum_exp(internal::log_Phi(z1), log_upper)); + + T_partials_return lcdf = sum(lcdf_elt); + + if constexpr (any_ad) { + // exp(2 lambda / mu) phi(z2) == phi(z1) exactly, because + // z2^2 - z1^2 == 4 lambda / mu, which collapses both phi terms into one. + // + // Each partial is masked separately: at y == inf the shape partial is + // 0 / 0 on its own. The inner select covers elements whose log CDF has + // saturated to -inf at an interior y. + const auto& w_dens = to_ref(exp(internal::log_phi(z1) - lcdf_elt)); + const auto& w_upper = to_ref(exp(log_upper - lcdf_elt)); + if constexpr (is_autodiff_v) { + partials<0>(ops_partials) = w_dens * sqrt_lambda_over_y * inv_y; + } + if constexpr (is_autodiff_v) { + partials<1>(ops_partials) = -2.0 * lambda_val * w_upper * square(inv_mu); + } + if constexpr (is_autodiff_v) { + partials<2>(ops_partials) + = 2.0 * w_upper * inv_mu - w_dens * inv_y / sqrt_lambda_over_y; + } + } + return ops_partials.build(lcdf); +} + +} // namespace math +} // namespace stan +#endif From 4b633b161202f3f6b35138fd057bdd540580dfcd Mon Sep 17 00:00:00 2001 From: Gidon Frischkorn Date: Sat, 22 Aug 2026 18:57:44 +0200 Subject: [PATCH 03/11] test: add test/prob fixtures and prim, rev and mix unit tests --- .../inv_gaussian_ccdf_log_test.hpp | 100 +++++++ .../inv_gaussian_cdf_log_test.hpp | 99 +++++++ .../inv_gaussian/inv_gaussian_cdf_test.hpp | 97 +++++++ test/prob/inv_gaussian/inv_gaussian_test.hpp | 99 +++++++ test/unit/math/mix/prob/inv_gaussian_test.cpp | 110 ++++++++ .../unit/math/prim/prob/inv_gaussian_test.cpp | 266 ++++++++++++++++++ test/unit/math/rev/prob/inv_gaussian_test.cpp | 161 +++++++++++ 7 files changed, 932 insertions(+) create mode 100644 test/prob/inv_gaussian/inv_gaussian_ccdf_log_test.hpp create mode 100644 test/prob/inv_gaussian/inv_gaussian_cdf_log_test.hpp create mode 100644 test/prob/inv_gaussian/inv_gaussian_cdf_test.hpp create mode 100644 test/prob/inv_gaussian/inv_gaussian_test.hpp create mode 100644 test/unit/math/mix/prob/inv_gaussian_test.cpp create mode 100644 test/unit/math/prim/prob/inv_gaussian_test.cpp create mode 100644 test/unit/math/rev/prob/inv_gaussian_test.cpp diff --git a/test/prob/inv_gaussian/inv_gaussian_ccdf_log_test.hpp b/test/prob/inv_gaussian/inv_gaussian_ccdf_log_test.hpp new file mode 100644 index 00000000000..2505297aac9 --- /dev/null +++ b/test/prob/inv_gaussian/inv_gaussian_ccdf_log_test.hpp @@ -0,0 +1,100 @@ +// Arguments: Doubles, Doubles, Doubles +#include +#include +#include +#include +#include + +using stan::math::var; +using std::numeric_limits; +using std::vector; + +class AgradCcdfLogInvGaussian : public AgradCcdfLogTest { + public: + void valid_values(vector >& parameters, + vector& ccdf_log) { + vector param(3); + + param[0] = 1.2; // y + param[1] = 0.5; // mu + param[2] = 2.0; // lambda + parameters.push_back(param); + ccdf_log.push_back(-3.9949836760335225255); // expected ccdf_log + + param[0] = 0.3; // y + param[1] = 1.0; // mu + param[2] = 5.0; // lambda + parameters.push_back(param); + ccdf_log.push_back(-0.003364845660995599504277); // expected ccdf_log + } + + void invalid_values(vector& index, vector& value) { + // y + index.push_back(0U); + value.push_back(-1.0); + + index.push_back(0U); + value.push_back(-numeric_limits::infinity()); + + // mu + index.push_back(1U); + value.push_back(0.0); + + index.push_back(1U); + value.push_back(-1.0); + + index.push_back(1U); + value.push_back(numeric_limits::infinity()); + + index.push_back(1U); + value.push_back(-numeric_limits::infinity()); + + // lambda + index.push_back(2U); + value.push_back(0.0); + + index.push_back(2U); + value.push_back(-1.0); + + index.push_back(2U); + value.push_back(numeric_limits::infinity()); + + index.push_back(2U); + value.push_back(-numeric_limits::infinity()); + } + + bool has_lower_bound() { return true; } + + double lower_bound() { return 0.0; } + + bool has_upper_bound() { return false; } + + template + stan::return_type_t ccdf_log(const T_y& y, + const T_loc& mu, + const T_shape& lambda, + const T3&, const T4&, + const T5&) { + return stan::math::inv_gaussian_lccdf(y, mu, lambda); + } + + template + stan::return_type_t ccdf_log_function( + const T_y& y, const T_loc& mu, const T_shape& lambda, const T3&, + const T4&, const T5&) { + using stan::math::erfc; + using stan::math::INV_SQRT_TWO; + using std::exp; + using std::log; + using std::sqrt; + + // Linear-space reference built from erfc, independent of the log-space + // implementation under test. + return log(1.0 + - 0.5 * erfc(-sqrt(lambda / y) * (y / mu - 1.0) * INV_SQRT_TWO) + - exp(2.0 * lambda / mu) * 0.5 + * erfc(sqrt(lambda / y) * (y / mu + 1.0) * INV_SQRT_TWO)); + } +}; diff --git a/test/prob/inv_gaussian/inv_gaussian_cdf_log_test.hpp b/test/prob/inv_gaussian/inv_gaussian_cdf_log_test.hpp new file mode 100644 index 00000000000..9d7273e7eec --- /dev/null +++ b/test/prob/inv_gaussian/inv_gaussian_cdf_log_test.hpp @@ -0,0 +1,99 @@ +// Arguments: Doubles, Doubles, Doubles +#include +#include +#include +#include +#include + +using stan::math::var; +using std::numeric_limits; +using std::vector; + +class AgradCdfLogInvGaussian : public AgradCdfLogTest { + public: + void valid_values(vector >& parameters, + vector& cdf_log) { + vector param(3); + + param[0] = 1.2; // y + param[1] = 0.5; // mu + param[2] = 2.0; // lambda + parameters.push_back(param); + cdf_log.push_back(-0.01857927772712011847827); // expected cdf_log + + param[0] = 0.3; // y + param[1] = 1.0; // mu + param[2] = 5.0; // lambda + parameters.push_back(param); + cdf_log.push_back(-5.696055133984662592139); // expected cdf_log + } + + void invalid_values(vector& index, vector& value) { + // y + index.push_back(0U); + value.push_back(-1.0); + + index.push_back(0U); + value.push_back(-numeric_limits::infinity()); + + // mu + index.push_back(1U); + value.push_back(0.0); + + index.push_back(1U); + value.push_back(-1.0); + + index.push_back(1U); + value.push_back(numeric_limits::infinity()); + + index.push_back(1U); + value.push_back(-numeric_limits::infinity()); + + // lambda + index.push_back(2U); + value.push_back(0.0); + + index.push_back(2U); + value.push_back(-1.0); + + index.push_back(2U); + value.push_back(numeric_limits::infinity()); + + index.push_back(2U); + value.push_back(-numeric_limits::infinity()); + } + + bool has_lower_bound() { return true; } + + double lower_bound() { return 0.0; } + + bool has_upper_bound() { return false; } + + template + stan::return_type_t cdf_log(const T_y& y, + const T_loc& mu, + const T_shape& lambda, + const T3&, const T4&, + const T5&) { + return stan::math::inv_gaussian_lcdf(y, mu, lambda); + } + + template + stan::return_type_t cdf_log_function( + const T_y& y, const T_loc& mu, const T_shape& lambda, const T3&, + const T4&, const T5&) { + using stan::math::erfc; + using stan::math::INV_SQRT_TWO; + using std::exp; + using std::log; + using std::sqrt; + + // Linear-space reference built from erfc, independent of the log-space + // implementation under test. + return log(0.5 * erfc(-sqrt(lambda / y) * (y / mu - 1.0) * INV_SQRT_TWO) + + exp(2.0 * lambda / mu) * 0.5 + * erfc(sqrt(lambda / y) * (y / mu + 1.0) * INV_SQRT_TWO)); + } +}; diff --git a/test/prob/inv_gaussian/inv_gaussian_cdf_test.hpp b/test/prob/inv_gaussian/inv_gaussian_cdf_test.hpp new file mode 100644 index 00000000000..0650973b8be --- /dev/null +++ b/test/prob/inv_gaussian/inv_gaussian_cdf_test.hpp @@ -0,0 +1,97 @@ +// Arguments: Doubles, Doubles, Doubles +#include +#include +#include +#include +#include + +using stan::math::var; +using std::numeric_limits; +using std::vector; + +class AgradCdfInvGaussian : public AgradCdfTest { + public: + void valid_values(vector >& parameters, vector& cdf) { + vector param(3); + + param[0] = 1.2; // y + param[1] = 0.5; // mu + param[2] = 2.0; // lambda + parameters.push_back(param); + cdf.push_back(0.9815922531042920910742); // expected cdf + + param[0] = 0.3; // y + param[1] = 1.0; // mu + param[2] = 5.0; // lambda + parameters.push_back(param); + cdf.push_back(0.003359190912064955560317); // expected cdf + } + + void invalid_values(vector& index, vector& value) { + // y + index.push_back(0U); + value.push_back(-1.0); + + index.push_back(0U); + value.push_back(-numeric_limits::infinity()); + + // mu + index.push_back(1U); + value.push_back(0.0); + + index.push_back(1U); + value.push_back(-1.0); + + index.push_back(1U); + value.push_back(numeric_limits::infinity()); + + index.push_back(1U); + value.push_back(-numeric_limits::infinity()); + + // lambda + index.push_back(2U); + value.push_back(0.0); + + index.push_back(2U); + value.push_back(-1.0); + + index.push_back(2U); + value.push_back(numeric_limits::infinity()); + + index.push_back(2U); + value.push_back(-numeric_limits::infinity()); + } + + bool has_lower_bound() { return true; } + + double lower_bound() { return 0.0; } + + bool has_upper_bound() { return false; } + + template + stan::return_type_t cdf(const T_y& y, const T_loc& mu, + const T_shape& lambda, const T3&, + const T4&, const T5&) { + return stan::math::inv_gaussian_cdf(y, mu, lambda); + } + + template + stan::return_type_t cdf_function(const T_y& y, + const T_loc& mu, + const T_shape& lambda, + const T3&, const T4&, + const T5&) { + using stan::math::erfc; + using stan::math::INV_SQRT_TWO; + using std::exp; + using std::sqrt; + + // Linear-space reference built from erfc, independent of the log-space + // implementation under test. + return 0.5 * erfc(-sqrt(lambda / y) * (y / mu - 1.0) * INV_SQRT_TWO) + + exp(2.0 * lambda / mu) * 0.5 + * erfc(sqrt(lambda / y) * (y / mu + 1.0) * INV_SQRT_TWO); + } +}; diff --git a/test/prob/inv_gaussian/inv_gaussian_test.hpp b/test/prob/inv_gaussian/inv_gaussian_test.hpp new file mode 100644 index 00000000000..68f3c5e6d69 --- /dev/null +++ b/test/prob/inv_gaussian/inv_gaussian_test.hpp @@ -0,0 +1,99 @@ +// Arguments: Doubles, Doubles, Doubles +#include +#include +#include +#include +#include + +using stan::math::var; +using std::numeric_limits; +using std::vector; + +class AgradDistributionsInvGaussian : public AgradDistributionTest { + public: + void valid_values(vector >& parameters, + vector& log_prob) { + vector param(3); + + param[0] = 1.2; // y + param[1] = 0.5; // mu + param[2] = 2.0; // lambda + parameters.push_back(param); + log_prob.push_back(-2.479180611448965157415); // expected log_prob + + param[0] = 0.3; // y + param[1] = 1.0; // mu + param[2] = 5.0; // lambda + parameters.push_back(param); + log_prob.push_back(-2.391593703832052124008); // expected log_prob + } + + void invalid_values(vector& index, vector& value) { + // y + index.push_back(0U); + value.push_back(-1.0); + + index.push_back(0U); + value.push_back(-numeric_limits::infinity()); + + // mu + index.push_back(1U); + value.push_back(0.0); + + index.push_back(1U); + value.push_back(-1.0); + + index.push_back(1U); + value.push_back(numeric_limits::infinity()); + + index.push_back(1U); + value.push_back(-numeric_limits::infinity()); + + // lambda + index.push_back(2U); + value.push_back(0.0); + + index.push_back(2U); + value.push_back(-1.0); + + index.push_back(2U); + value.push_back(numeric_limits::infinity()); + + index.push_back(2U); + value.push_back(-numeric_limits::infinity()); + } + + template + stan::return_type_t log_prob(const T_y& y, + const T_loc& mu, + const T_shape& lambda, + const T3&, const T4&, + const T5&) { + return stan::math::inv_gaussian_lpdf(y, mu, lambda); + } + + template + stan::return_type_t log_prob(const T_y& y, + const T_loc& mu, + const T_shape& lambda, + const T3&, const T4&, + const T5&) { + return stan::math::inv_gaussian_lpdf(y, mu, lambda); + } + + template + stan::return_type_t log_prob_function( + const T_y& y, const T_loc& mu, const T_shape& lambda, const T3&, + const T4&, const T5&) { + using stan::math::pi; + using stan::math::square; + using std::log; + using std::sqrt; + + return log(sqrt(lambda / (2.0 * pi() * y * y * y))) + - lambda * square(y - mu) / (2.0 * square(mu) * y); + } +}; diff --git a/test/unit/math/mix/prob/inv_gaussian_test.cpp b/test/unit/math/mix/prob/inv_gaussian_test.cpp new file mode 100644 index 00000000000..c9a6fb49150 --- /dev/null +++ b/test/unit/math/mix/prob/inv_gaussian_test.cpp @@ -0,0 +1,110 @@ +#include +#include +#include + +TEST_F(AgradRev, mathMixScalFun_inv_gaussian_lpdf) { + auto f = [](const auto& y, const auto& mu, const auto& lambda) { + return stan::math::inv_gaussian_lpdf(y, mu, lambda); + }; + + stan::test::expect_ad(f, 0.3, 0.5, 2.0); + stan::test::expect_ad(f, 1.2, 0.5, 2.0); + stan::test::expect_ad(f, 0.3, 1.0, 5.0); + stan::test::expect_ad(f, 5.0, 1.0, 0.5); + stan::test::expect_ad(f, 1e-3, 1.0, 2.0); + // out of support and invalid parameters + stan::test::expect_ad(f, -1.0, 1.0, 2.0); + stan::test::expect_ad(f, 1.0, 0.0, 2.0); + stan::test::expect_ad(f, 1.0, 1.0, 0.0); +} + +TEST_F(AgradRev, mathMixScalFun_inv_gaussian_cdf) { + auto f = [](const auto& y, const auto& mu, const auto& lambda) { + return stan::math::inv_gaussian_cdf(y, mu, lambda); + }; + + stan::test::expect_ad(f, 0.3, 0.5, 2.0); + stan::test::expect_ad(f, 1.2, 0.5, 2.0); + stan::test::expect_ad(f, 0.3, 1.0, 5.0); + stan::test::expect_ad(f, 5.0, 1.0, 0.5); + // out of support and invalid parameters + stan::test::expect_ad(f, -1.0, 1.0, 2.0); + stan::test::expect_ad(f, 1.0, 0.0, 2.0); +} + +TEST_F(AgradRev, mathMixScalFun_inv_gaussian_lcdf) { + auto f = [](const auto& y, const auto& mu, const auto& lambda) { + return stan::math::inv_gaussian_lcdf(y, mu, lambda); + }; + + stan::test::expect_ad(f, 0.3, 0.5, 2.0); + stan::test::expect_ad(f, 1.2, 0.5, 2.0); + stan::test::expect_ad(f, 0.3, 1.0, 5.0); + stan::test::expect_ad(f, 5.0, 1.0, 0.5); + // out of support and invalid parameters + stan::test::expect_ad(f, -1.0, 1.0, 2.0); + stan::test::expect_ad(f, 1.0, 0.0, 2.0); +} + +TEST_F(AgradRev, mathMixScalFun_inv_gaussian_lccdf) { + auto f = [](const auto& y, const auto& mu, const auto& lambda) { + return stan::math::inv_gaussian_lccdf(y, mu, lambda); + }; + + stan::test::expect_ad(f, 0.3, 0.5, 2.0); + stan::test::expect_ad(f, 1.2, 0.5, 2.0); + stan::test::expect_ad(f, 0.3, 1.0, 5.0); + stan::test::expect_ad(f, 5.0, 1.0, 0.5); + // out of support and invalid parameters + stan::test::expect_ad(f, -1.0, 1.0, 2.0); + stan::test::expect_ad(f, 1.0, 0.0, 2.0); +} + +// crosses the internal log_Phi branch at z = -30 +TEST_F(AgradRev, mathMixScalFun_inv_gaussian_lcdf_tails) { + auto f = [](const auto& y, const auto& mu, const auto& lambda) { + return stan::math::inv_gaussian_lcdf(y, mu, lambda); + }; + + stan::test::expect_ad(f, 0.1, 0.1, 100.0); + stan::test::expect_ad(f, 0.05, 1.0, 100.0); + stan::test::expect_ad(f, 0.15, 0.1, 10.0); +} + +TEST_F(AgradRev, mathMixScalFun_inv_gaussian_lccdf_tails) { + auto f = [](const auto& y, const auto& mu, const auto& lambda) { + return stan::math::inv_gaussian_lccdf(y, mu, lambda); + }; + + stan::test::expect_ad(f, 0.1, 0.1, 100.0); + stan::test::expect_ad(f, 0.15, 0.1, 10.0); + stan::test::expect_ad(f, 0.5, 0.1, 50.0); +} + +TEST_F(AgradRev, mathMixScalFun_inv_gaussian_lpdf_vectorized) { + Eigen::VectorXd y(3); + y << 0.3, 1.2, 4.0; + Eigen::VectorXd mu(3); + mu << 0.5, 1.0, 2.0; + Eigen::VectorXd lambda(3); + lambda << 2.0, 5.0, 0.7; + + auto f = [](const auto& y, const auto& mu, const auto& lambda) { + return stan::math::inv_gaussian_lpdf(y, mu, lambda); + }; + stan::test::expect_ad(f, y, mu, lambda); +} + +TEST_F(AgradRev, mathMixScalFun_inv_gaussian_lccdf_vectorized) { + Eigen::VectorXd y(3); + y << 0.3, 1.2, 4.0; + Eigen::VectorXd mu(3); + mu << 0.5, 1.0, 2.0; + Eigen::VectorXd lambda(3); + lambda << 2.0, 5.0, 0.7; + + auto f = [](const auto& y, const auto& mu, const auto& lambda) { + return stan::math::inv_gaussian_lccdf(y, mu, lambda); + }; + stan::test::expect_ad(f, y, mu, lambda); +} diff --git a/test/unit/math/prim/prob/inv_gaussian_test.cpp b/test/unit/math/prim/prob/inv_gaussian_test.cpp new file mode 100644 index 00000000000..8dbaaebba9d --- /dev/null +++ b/test/unit/math/prim/prob/inv_gaussian_test.cpp @@ -0,0 +1,266 @@ +#include +#include +#include +#include +#include +#include +#include +#include + +class InvGaussianTestRig : public VectorRealRNGTestRig { + public: + InvGaussianTestRig() + : VectorRealRNGTestRig(10000, 10, {0.1, 0.5, 1.0, 2.5, 4.0}, {1, 2, 3, 4}, + {0.0, -0.1, -1.0}, {0, -1, -2}, + {0.2, 1.0, 3.0, 7.5}, {1, 2, 3, 8}, + {0.0, -0.1, -1.0}, {0, -1, -2}) {} + + template + auto generate_samples(const T1& mu, const T2& lambda, const T3& unused, + T_rng& rng) const { + return stan::math::inv_gaussian_rng(mu, lambda, rng); + } + + std::vector generate_quantiles(double mu, double lambda, + double unused) const { + std::vector quantiles; + double K = stan::math::round(2 * std::pow(N_, 0.4)); + boost::math::inverse_gaussian_distribution<> dist(mu, lambda); + for (int i = 1; i < K; ++i) { + double frac = i / K; + quantiles.push_back(quantile(dist, frac)); + } + quantiles.push_back(std::numeric_limits::max()); + return quantiles; + } +}; + +TEST(ProbDistributionsInvGaussian, errorCheck) { + check_dist_throws_all_types(InvGaussianTestRig()); +} + +TEST(ProbDistributionsInvGaussian, distributionCheck) { + check_quantiles_real_real(InvGaussianTestRig()); +} + +TEST(ProbDistributionsInvGaussian, error_check) { + boost::random::mt19937 rng; + EXPECT_NO_THROW(stan::math::inv_gaussian_rng(1.0, 2.0, rng)); + EXPECT_THROW(stan::math::inv_gaussian_rng(0.0, 2.0, rng), std::domain_error); + EXPECT_THROW(stan::math::inv_gaussian_rng(-1.0, 2.0, rng), std::domain_error); + EXPECT_THROW(stan::math::inv_gaussian_rng(1.0, 0.0, rng), std::domain_error); + EXPECT_THROW(stan::math::inv_gaussian_rng(1.0, -1.0, rng), std::domain_error); + EXPECT_THROW( + stan::math::inv_gaussian_rng(stan::math::positive_infinity(), 2.0, rng), + std::domain_error); + EXPECT_THROW( + stan::math::inv_gaussian_rng(1.0, stan::math::positive_infinity(), rng), + std::domain_error); +} + +// reference values from mpmath at 60 digits + +TEST(ProbDistributionsInvGaussian, values) { + using stan::math::inv_gaussian_cdf; + using stan::math::inv_gaussian_lccdf; + using stan::math::inv_gaussian_lcdf; + using stan::math::inv_gaussian_lpdf; + + EXPECT_FLOAT_EQ(-2.479180611448965157415, inv_gaussian_lpdf(1.2, 0.5, 2.0)); + EXPECT_FLOAT_EQ(-2.391593703832052124008, inv_gaussian_lpdf(0.3, 1.0, 5.0)); + EXPECT_FLOAT_EQ(0.9815922531042920910742, inv_gaussian_cdf(1.2, 0.5, 2.0)); + EXPECT_FLOAT_EQ(0.003359190912064955560317, inv_gaussian_cdf(0.3, 1.0, 5.0)); + EXPECT_FLOAT_EQ(-0.01857927772712011847827, inv_gaussian_lcdf(1.2, 0.5, 2.0)); + EXPECT_FLOAT_EQ(-5.696055133984662592139, inv_gaussian_lcdf(0.3, 1.0, 5.0)); + EXPECT_FLOAT_EQ(-3.9949836760335225255, inv_gaussian_lccdf(1.2, 0.5, 2.0)); + EXPECT_FLOAT_EQ(-0.003364845660995599504277, + inv_gaussian_lccdf(0.3, 1.0, 5.0)); +} + +TEST(ProbDistributionsInvGaussian, boundaries) { + using stan::math::inv_gaussian_cdf; + using stan::math::inv_gaussian_lccdf; + using stan::math::inv_gaussian_lcdf; + using stan::math::inv_gaussian_lpdf; + double inf = std::numeric_limits::infinity(); + + EXPECT_FLOAT_EQ(-inf, inv_gaussian_lpdf(0.0, 1.0, 2.0)); + EXPECT_FLOAT_EQ(-inf, inv_gaussian_lcdf(0.0, 1.0, 2.0)); + EXPECT_FLOAT_EQ(0.0, inv_gaussian_lccdf(0.0, 1.0, 2.0)); + EXPECT_FLOAT_EQ(0.0, inv_gaussian_cdf(0.0, 1.0, 2.0)); + EXPECT_FLOAT_EQ(0.0, inv_gaussian_lcdf(inf, 1.0, 2.0)); + EXPECT_FLOAT_EQ(-inf, inv_gaussian_lccdf(inf, 1.0, 2.0)); + EXPECT_FLOAT_EQ(1.0, inv_gaussian_cdf(inf, 1.0, 2.0)); +} + +// 2 lambda / mu is past the overflow point of exp for all of these. +TEST(ProbDistributionsInvGaussian, largeExpFactor) { + using stan::math::inv_gaussian_lccdf; + using stan::math::inv_gaussian_lcdf; + + // 2 lambda / mu = 2000 + EXPECT_FLOAT_EQ(-0.68061354470344461635, inv_gaussian_lcdf(0.1, 0.1, 100.0)); + EXPECT_FLOAT_EQ(-0.70583990450544626626, inv_gaussian_lccdf(0.1, 0.1, 100.0)); + EXPECT_FLOAT_EQ(-671.8777811569157, inv_gaussian_lccdf(0.3, 0.1, 100.0)); + EXPECT_FLOAT_EQ(-4057.1238032952494, inv_gaussian_lccdf(1.0, 0.1, 100.0)); + // 2 lambda / mu = 200 + EXPECT_FLOAT_EQ(-805.707747545284162, inv_gaussian_lccdf(0.5, 0.1, 50.0)); +} + +// erfc underflows to zero here; the asymptotic log_Phi branch carries these. +TEST(ProbDistributionsInvGaussian, deepLowerTail) { + using stan::math::inv_gaussian_lcdf; + + EXPECT_FLOAT_EQ(-906.524245083721553, + inv_gaussian_lcdf(0.05, 1.0, 100.0)); // z1 = -42.5 + EXPECT_FLOAT_EQ(-849.06768101617034, + inv_gaussian_lcdf(0.1, 0.75, 225.0)); // z1 = -41.1 + EXPECT_FLOAT_EQ(-5333.88923087778086, + inv_gaussian_lcdf(0.02, 0.75, 225.0)); // z1 = -103.2 + EXPECT_FLOAT_EQ(-4905.33096155861673, + inv_gaussian_lcdf(0.01, 1.0, 100.0)); // z1 = -99.0 +} + +// Once Phi(z1) rounds to exactly one the log CDF lands a few 1e-20 above +// zero. That is bounded by the rounding of Phi and vanishes under exp, so it +// is asserted at that scale. +TEST(ProbDistributionsInvGaussian, probabilityNeverExceedsOne) { + using stan::math::inv_gaussian_lccdf; + using stan::math::inv_gaussian_lcdf; + for (double mu : {1e-3, 1e-2, 0.1, 1.0, 10.0, 1e3}) { + for (double lambda : {1.0, 1e3, 1e8, 1e14, 1e16}) { + for (double rel : {1e-3, 0.1, 0.5, 1.0, 2.0, 10.0, 1e3}) { + double y = mu * rel; + double lf = inv_gaussian_lcdf(y, mu, lambda); + double ls = inv_gaussian_lccdf(y, mu, lambda); + EXPECT_LE(lf, 1e-15); + EXPECT_LE(ls, 1e-15); + EXPECT_LE(std::exp(lf), 1.0); + EXPECT_LE(std::exp(ls), 1.0); + } + } + } + // against mpmath at 60 digits + EXPECT_NEAR(-0.693147179298379049, inv_gaussian_lcdf(0.1, 0.1, 1e16), 1e-12); + EXPECT_FLOAT_EQ(-4.04999999999999999e17, inv_gaussian_lcdf(1e-4, 1e-3, 1e14)); + EXPECT_FLOAT_EQ(-4.99000499999999997e17, inv_gaussian_lcdf(1e-4, 1e-1, 1e14)); +} + +// At y == mu the CDF is a fixed quantity independent of the scale of mu, so +// all of these must agree. +TEST(ProbDistributionsInvGaussian, medianIsScaleInvariant) { + using stan::math::inv_gaussian_lccdf; + using stan::math::inv_gaussian_lcdf; + for (double lambda_over_mu : {1e2, 1e6, 1e11, 1e13}) { + double ref_lcdf = 0; + double ref_lccdf = 0; + bool first = true; + for (double mu : {1e-3, 1e-2, 0.1, 1.0, 10.0}) { + double lambda = lambda_over_mu * mu; + double a = inv_gaussian_lcdf(mu, mu, lambda); + double b = inv_gaussian_lccdf(mu, mu, lambda); + if (first) { + ref_lcdf = a; + ref_lccdf = b; + first = false; + } else { + EXPECT_NEAR(ref_lcdf, a, 1e-13); + EXPECT_NEAR(ref_lccdf, b, 1e-13); + } + } + } + // F(mu) for a very large shape approaches 1/2 from above + EXPECT_NEAR(-0.693147181822, inv_gaussian_lccdf(1e-3, 1e-3, 1e14), 1e-11); + EXPECT_NEAR(-0.693147180686, inv_gaussian_lccdf(1e-3, 1e-3, 1e16), 1e-11); +} + +TEST(ProbDistributionsInvGaussian, cdfCcdfSumToOne) { + using stan::math::inv_gaussian_lccdf; + using stan::math::inv_gaussian_lcdf; + std::vector mus{0.5, 1.0, 0.3, 2.0}; + std::vector lambdas{2.0, 5.0, 10.0, 0.7}; + std::vector ys{0.1, 0.3, 0.8, 1.5, 4.0}; + for (double mu : mus) { + for (double lambda : lambdas) { + for (double y : ys) { + double f = std::exp(inv_gaussian_lcdf(y, mu, lambda)); + double s = std::exp(inv_gaussian_lccdf(y, mu, lambda)); + EXPECT_NEAR(1.0, f + s, 1e-12); + } + } + } +} + +TEST(ProbDistributionsInvGaussian, vectorMatchesScalarSum) { + using stan::math::inv_gaussian_lccdf; + using stan::math::inv_gaussian_lcdf; + using stan::math::inv_gaussian_lpdf; + std::vector y{0.2, 0.5, 1.0, 2.0, 5.0}; + double sum_lpdf = 0; + double sum_lcdf = 0; + double sum_lccdf = 0; + for (double yi : y) { + sum_lpdf += inv_gaussian_lpdf(yi, 1.3, 3.0); + sum_lcdf += inv_gaussian_lcdf(yi, 1.3, 3.0); + sum_lccdf += inv_gaussian_lccdf(yi, 1.3, 3.0); + } + EXPECT_FLOAT_EQ(sum_lpdf, inv_gaussian_lpdf(y, 1.3, 3.0)); + EXPECT_FLOAT_EQ(sum_lcdf, inv_gaussian_lcdf(y, 1.3, 3.0)); + EXPECT_FLOAT_EQ(sum_lccdf, inv_gaussian_lccdf(y, 1.3, 3.0)); +} + +// check helper functions; tolerances scale with each value's magnitude +TEST(ProbDistributionsInvGaussian, internalLogPhi) { + using stan::math::internal::log_Phi; + + // erfc branch, up to the switch at z = -30 + EXPECT_NEAR(-0.6931471805599453094, log_Phi(0.0), 1e-15); + EXPECT_NEAR(-1.841021645009263506, log_Phi(-1.0), 1e-14); + EXPECT_NEAR(-15.06499839398872574, log_Phi(-5.0), 1e-13); + EXPECT_NEAR(-53.23128515051247058, log_Phi(-10.0), 1e-12); + EXPECT_NEAR(-203.9171553710972639, log_Phi(-20.0), 1e-11); + EXPECT_NEAR(-451.3229124585286345, log_Phi(-29.9), 1e-11); + // asymptotic branch + EXPECT_NEAR(-454.3212439563431971, log_Phi(-30.0), 1e-11); + EXPECT_NEAR(-457.3295644163822579, log_Phi(-30.1), 1e-11); + // beyond the point where erfc underflows to zero (z < -37.5) + EXPECT_NEAR(-745.6952702904110813, log_Phi(-38.5), 1e-10); + EXPECT_NEAR(-804.6084420137537882, log_Phi(-40.0), 1e-10); + EXPECT_NEAR(-1805.013560680567139, log_Phi(-60.0), 1e-9); + EXPECT_NEAR(-11255.92961826680818, log_Phi(-150.0), 1e-8); + EXPECT_NEAR(-500007.8266948121843, log_Phi(-1000.0), 1e-6); + // upper tail saturates at log(1) = 0 + EXPECT_NEAR(-2.866516129637635934e-7, log_Phi(5.0), 1e-16); + EXPECT_FLOAT_EQ(0.0, log_Phi(50.0)); + + double inf = std::numeric_limits::infinity(); + EXPECT_FLOAT_EQ(0.0, log_Phi(inf)); + EXPECT_FLOAT_EQ(-inf, log_Phi(-inf)); + EXPECT_TRUE(std::isnan(log_Phi(std::numeric_limits::quiet_NaN()))); +} + +// The true slope at z = -30 is about 30, so across a 2e-10 interval the +// honest change is about 6e-9. +TEST(ProbDistributionsInvGaussian, internalLogPhiBranchContinuity) { + using stan::math::internal::log_Phi; + double eps = 1e-10; + double step = log_Phi(-30.0 + eps) - log_Phi(-30.0 - eps); + EXPECT_LT(std::fabs(step), 1e-7); +} + +TEST(ProbDistributionsInvGaussian, errors) { + using stan::math::inv_gaussian_lpdf; + double inf = std::numeric_limits::infinity(); + + EXPECT_THROW(inv_gaussian_lpdf(-1.0, 1.0, 2.0), std::domain_error); + EXPECT_THROW(inv_gaussian_lpdf(1.0, 0.0, 2.0), std::domain_error); + EXPECT_THROW(inv_gaussian_lpdf(1.0, -1.0, 2.0), std::domain_error); + EXPECT_THROW(inv_gaussian_lpdf(1.0, inf, 2.0), std::domain_error); + EXPECT_THROW(inv_gaussian_lpdf(1.0, 1.0, 0.0), std::domain_error); + EXPECT_THROW(inv_gaussian_lpdf(1.0, 1.0, -1.0), std::domain_error); + EXPECT_THROW(inv_gaussian_lpdf(1.0, 1.0, inf), std::domain_error); + + std::vector y{1.0, 2.0}; + std::vector mu{1.0, 2.0, 3.0}; + EXPECT_THROW(inv_gaussian_lpdf(y, mu, 1.0), std::invalid_argument); +} diff --git a/test/unit/math/rev/prob/inv_gaussian_test.cpp b/test/unit/math/rev/prob/inv_gaussian_test.cpp new file mode 100644 index 00000000000..a6c9493c827 --- /dev/null +++ b/test/unit/math/rev/prob/inv_gaussian_test.cpp @@ -0,0 +1,161 @@ +#include +#include +#include +#include +#include +#include + +namespace { + +struct grad_case { + double y; + double mu; + double lambda; + double d_y; + double d_mu; + double d_lambda; +}; + +template +void check_grads(const F& f, const grad_case& c, double tol) { + stan::math::var y = c.y; + stan::math::var mu = c.mu; + stan::math::var lambda = c.lambda; + stan::math::var out = f(y, mu, lambda); + out.grad(); + EXPECT_NEAR(c.d_y, y.adj(), tol); + EXPECT_NEAR(c.d_mu, mu.adj(), tol); + EXPECT_NEAR(c.d_lambda, lambda.adj(), tol); + stan::math::recover_memory(); +} + +} // namespace + +// reference values from mpmath at 60 digits +TEST_F(AgradRev, inv_gaussian_lpdf_gradients) { + auto f = [](const auto& y, const auto& mu, const auto& lambda) { + return stan::math::inv_gaussian_lpdf(y, mu, lambda); + }; + check_grads(f, + {1.2, 0.5, 2.0, -4.5555555555555556, 11.199999999999999, + -0.56666666666666659}, + 1e-12); + check_grads(f, + {0.3, 1.0, 5.0, 20.27777777777778, -3.5000000000000001, + -0.71666666666666672}, + 1e-12); +} + +TEST_F(AgradRev, inv_gaussian_lcdf_gradients) { + auto f = [](const auto& y, const auto& mu, const auto& lambda) { + return stan::math::inv_gaussian_lcdf(y, mu, lambda); + }; + check_grads(f, + {1.2, 0.5, 2.0, 0.085383591493919688, -0.27616885497947743, + 0.017812058848517546}, + 1e-12); + check_grads(f, + {0.3, 1.0, 5.0, 27.233870297028338, -3.6491778320016701, + -0.90419665142136623}, + 1e-11); +} + +TEST_F(AgradRev, inv_gaussian_lccdf_gradients) { + auto f = [](const auto& y, const auto& mu, const auto& lambda) { + return stan::math::inv_gaussian_lccdf(y, mu, lambda); + }; + check_grads(f, + {1.2, 0.5, 2.0, -4.5530760732154196, 14.72669143770909, + -0.94982721549802085}, + 1e-11); + check_grads(f, + {0.3, 1.0, 5.0, -0.09179211684684445, 0.012299601720089049, + 0.003047606666792857}, + 1e-12); +} + +// A batch of N observations produces one node on the autodiff tape, so tape +// growth beyond the input vars must not depend on N. +TEST_F(AgradRev, inv_gaussian_one_node_per_call) { + using stan::math::var; + using stan::math::vector_v; + + auto tape_growth = [](int n, auto&& f) { + vector_v y = Eigen::VectorXd::LinSpaced(n, 0.4, 2.0); + vector_v mu = Eigen::VectorXd::Constant(n, 1.0); + vector_v lambda = Eigen::VectorXd::Constant(n, 2.0); + std::size_t before + = stan::math::ChainableStack::instance_->var_stack_.size(); + var out = f(y, mu, lambda); + std::size_t after + = stan::math::ChainableStack::instance_->var_stack_.size(); + return after - before; + }; + + auto lpdf = [](const auto& y, const auto& mu, const auto& lambda) { + return stan::math::inv_gaussian_lpdf(y, mu, lambda); + }; + auto lcdf = [](const auto& y, const auto& mu, const auto& lambda) { + return stan::math::inv_gaussian_lcdf(y, mu, lambda); + }; + auto lccdf = [](const auto& y, const auto& mu, const auto& lambda) { + return stan::math::inv_gaussian_lccdf(y, mu, lambda); + }; + auto cdf = [](const auto& y, const auto& mu, const auto& lambda) { + return stan::math::inv_gaussian_cdf(y, mu, lambda); + }; + + for (auto&& f : + {std::function( + lpdf), + std::function( + lcdf), + std::function( + lccdf), + std::function( + cdf)}) { + std::size_t g10 = tape_growth(10, f); + std::size_t g1000 = tape_growth(1000, f); + EXPECT_EQ(g10, g1000); + stan::math::recover_memory(); + } +} + +TEST_F(AgradRev, inv_gaussian_vectorized_matches_scalar) { + using stan::math::var; + using stan::math::vector_v; + Eigen::VectorXd y_d(4); + y_d << 0.3, 0.8, 1.5, 3.0; + Eigen::VectorXd mu_d(4); + mu_d << 0.5, 1.0, 1.0, 2.0; + Eigen::VectorXd lambda_d(4); + lambda_d << 2.0, 5.0, 1.0, 0.7; + + vector_v y = y_d; + vector_v mu = mu_d; + vector_v lambda = lambda_d; + var out = stan::math::inv_gaussian_lccdf(y, mu, lambda); + out.grad(); + Eigen::VectorXd vec_adj_y(4); + for (int i = 0; i < 4; ++i) { + vec_adj_y(i) = y(i).adj(); + } + double vec_val = out.val(); + stan::math::recover_memory(); + + double scalar_sum = 0; + Eigen::VectorXd scalar_adj_y(4); + for (int i = 0; i < 4; ++i) { + var yi = y_d(i); + var out_i = stan::math::inv_gaussian_lccdf(yi, mu_d(i), lambda_d(i)); + out_i.grad(); + scalar_sum += out_i.val(); + scalar_adj_y(i) = yi.adj(); + stan::math::recover_memory(); + } + + EXPECT_FLOAT_EQ(scalar_sum, vec_val); + for (int i = 0; i < 4; ++i) { + EXPECT_FLOAT_EQ(scalar_adj_y(i), vec_adj_y(i)); + } +} From a77ebe0b65efb480e7695eb937e85b893f0d6626 Mon Sep 17 00:00:00 2001 From: Gidon Frischkorn Date: Sun, 23 Aug 2026 09:22:31 +0200 Subject: [PATCH 04/11] fix: stabilize inv_gaussian_rng root against catastrophic cancellation - write the smaller root as a reciprocal; (1 + u/2)^2 - (u + u^2/4) == 1 - pin stability and the first two moments over mu/lambda up to 1e20 --- stan/math/prim/prob/inv_gaussian_rng.hpp | 15 ++++----- .../unit/math/prim/prob/inv_gaussian_test.cpp | 33 +++++++++++++++++++ 2 files changed, 40 insertions(+), 8 deletions(-) diff --git a/stan/math/prim/prob/inv_gaussian_rng.hpp b/stan/math/prim/prob/inv_gaussian_rng.hpp index 6fefb3aff0a..a2a7810369a 100644 --- a/stan/math/prim/prob/inv_gaussian_rng.hpp +++ b/stan/math/prim/prob/inv_gaussian_rng.hpp @@ -5,7 +5,6 @@ #include #include #include -#include #include #include #include @@ -31,9 +30,11 @@ namespace math { * candidate roots and a Bernoulli draw selects between them. The method is * exact, with no rejection step. * - *

The smaller root is formed from the discriminant - * \f$\sqrt{u + u^2/4}\f$ with \f$u = \mu w / \lambda\f$, as - * \f$x = \mu(1 + u/2 - \sqrt{u + u^2/4})\f$. + *

The smaller root is formed as + * \f$x = \mu / (1 + u/2 + \sqrt{u + u^2/4})\f$ with + * \f$u = \mu w / \lambda\f$. Since \f$(1 + u/2)^2 - (u + u^2/4) = 1\f$ + * exactly, this form subtracts nothing and stays accurate for \f$u\f$ up to + * \f$10^{20}\f$. * * @tparam T_loc type of mean parameter * @tparam T_shape type of shape parameter @@ -60,7 +61,7 @@ inv_gaussian_rng(const T_loc& mu, const T_shape& lambda, RNG& rng) { T_mu_ref mu_ref = mu; T_lambda_ref lambda_ref = lambda; check_positive_finite(function, "Mean parameter", mu_ref); - check_positive_finite(function, "Shape Parameter", lambda_ref); + check_positive_finite(function, "Shape parameter", lambda_ref); scalar_seq_view mu_vec(mu_ref); scalar_seq_view lambda_vec(lambda_ref); @@ -77,9 +78,7 @@ inv_gaussian_rng(const T_loc& mu, const T_shape& lambda, RNG& rng) { const double nu = norm_rng(); const double w = nu * nu; const double u = mu_dbl * w / lambda_dbl; - const double half_ratio = 0.5 * u; - const double disc = std::sqrt(u + half_ratio * half_ratio); - const double x = mu_dbl * (1.0 + half_ratio - disc); + const double x = mu_dbl / (1.0 + 0.5 * u + std::sqrt(u + 0.25 * u * u)); output[n] = (uniform01_rng() <= mu_dbl / (mu_dbl + x)) ? x : mu_dbl * mu_dbl / x; } diff --git a/test/unit/math/prim/prob/inv_gaussian_test.cpp b/test/unit/math/prim/prob/inv_gaussian_test.cpp index 8dbaaebba9d..3c0fba1e9bb 100644 --- a/test/unit/math/prim/prob/inv_gaussian_test.cpp +++ b/test/unit/math/prim/prob/inv_gaussian_test.cpp @@ -58,6 +58,39 @@ TEST(ProbDistributionsInvGaussian, error_check) { std::domain_error); } +TEST(ProbDistributionsInvGaussian, rngStableForLargeMuOverLambda) { + boost::random::mt19937 rng(1234); + for (double mu : {1.0, 1e3, 1e6, 1e9, 1e12}) { + for (double lambda : {1e-8, 1e-4, 1.0, 1e4}) { + for (int i = 0; i < 2000; ++i) { + double d = stan::math::inv_gaussian_rng(mu, lambda, rng); + ASSERT_TRUE(std::isfinite(d)) << "mu=" << mu << " lambda=" << lambda; + ASSERT_GT(d, 0.0) << "mu=" << mu << " lambda=" << lambda; + } + } + } +} + +// The sample variance of an inverse Gaussian is heavy tailed, so its +// tolerance is wider than the mean's. +TEST(ProbDistributionsInvGaussian, rngMomentsAtLargeMu) { + boost::random::mt19937 rng(4321); + const double mu = 1e3; + const double lambda = 1e3; + const int N = 200000; + double sum = 0; + double sum_sq = 0; + for (int i = 0; i < N; ++i) { + double d = stan::math::inv_gaussian_rng(mu, lambda, rng); + sum += d; + sum_sq += d * d; + } + double mean = sum / N; + double var = sum_sq / N - mean * mean; + EXPECT_NEAR(mu, mean, 0.05 * mu); + EXPECT_NEAR(mu * mu * mu / lambda, var, 0.15 * mu * mu * mu / lambda); +} + // reference values from mpmath at 60 digits TEST(ProbDistributionsInvGaussian, values) { From 7a60ee538fdaad050bfeb8b3aee02f740a451849 Mon Sep 17 00:00:00 2001 From: Gidon Frischkorn Date: Sun, 23 Aug 2026 11:38:07 +0200 Subject: [PATCH 05/11] fix: return -inf from inv_gaussian_lpdf at y == inf --- stan/math/prim/prob/inv_gaussian_lpdf.hpp | 9 +++++---- test/unit/math/prim/prob/inv_gaussian_test.cpp | 2 ++ 2 files changed, 7 insertions(+), 4 deletions(-) diff --git a/stan/math/prim/prob/inv_gaussian_lpdf.hpp b/stan/math/prim/prob/inv_gaussian_lpdf.hpp index 10c9bffcf65..86eb5f2e477 100644 --- a/stan/math/prim/prob/inv_gaussian_lpdf.hpp +++ b/stan/math/prim/prob/inv_gaussian_lpdf.hpp @@ -35,8 +35,8 @@ namespace math { * \exp(-\lambda (y - \mu)^2 / (2 \mu^2 y))\f$, with support * \f$y \in (0, \infty)\f$ and variance \f$\mu^3 / \lambda\f$. * - *

\f$y = 0\f$ is accepted and returns \f$-\infty\f$; a negative - * \f$y\f$ throws. + *

Both \f$y = 0\f$ and \f$y = \infty\f$ are accepted and return + * \f$-\infty\f$; only a negative \f$y\f$ throws. * * @tparam T_y type of scalar * @tparam T_loc type of mean parameter @@ -86,9 +86,10 @@ inline return_type_t inv_gaussian_lpdf( auto ops_partials = make_partials_propagator(y_ref, mu_ref, lambda_ref); - // y == 0 has zero density and LOG_ZERO absorbs the sum. + // Both boundaries of the support have zero density and LOG_ZERO absorbs the + // sum, so the whole container short-circuits exactly. // The gradients are technically ill-defined, but treated as zero. - if (sum(promote_scalar(y_val == 0))) { + if (sum(promote_scalar((y_val == 0) || (y_val == INFTY)))) { return ops_partials.build(LOG_ZERO); } diff --git a/test/unit/math/prim/prob/inv_gaussian_test.cpp b/test/unit/math/prim/prob/inv_gaussian_test.cpp index 3c0fba1e9bb..47e8c0230eb 100644 --- a/test/unit/math/prim/prob/inv_gaussian_test.cpp +++ b/test/unit/math/prim/prob/inv_gaussian_test.cpp @@ -121,6 +121,8 @@ TEST(ProbDistributionsInvGaussian, boundaries) { EXPECT_FLOAT_EQ(-inf, inv_gaussian_lcdf(0.0, 1.0, 2.0)); EXPECT_FLOAT_EQ(0.0, inv_gaussian_lccdf(0.0, 1.0, 2.0)); EXPECT_FLOAT_EQ(0.0, inv_gaussian_cdf(0.0, 1.0, 2.0)); + // y == inf is in the closure of the support: zero density, saturated CDF + EXPECT_FLOAT_EQ(-inf, inv_gaussian_lpdf(inf, 1.0, 2.0)); EXPECT_FLOAT_EQ(0.0, inv_gaussian_lcdf(inf, 1.0, 2.0)); EXPECT_FLOAT_EQ(-inf, inv_gaussian_lccdf(inf, 1.0, 2.0)); EXPECT_FLOAT_EQ(1.0, inv_gaussian_cdf(inf, 1.0, 2.0)); From 70568684a38e6da0f94c15c95f84f6c080ef4c72 Mon Sep 17 00:00:00 2001 From: Gidon Frischkorn Date: Sun, 23 Aug 2026 14:05:52 +0200 Subject: [PATCH 06/11] fix: handle y boundaries elementwise in inv_gaussian cdf/lcdf/lccdf - a boundary element no longer discards the rest of the container --- stan/math/prim/prob/inv_gaussian_cdf.hpp | 41 +++-- stan/math/prim/prob/inv_gaussian_lccdf.hpp | 42 +++-- stan/math/prim/prob/inv_gaussian_lcdf.hpp | 42 +++-- .../unit/math/prim/prob/inv_gaussian_test.cpp | 25 +++ test/unit/math/rev/prob/inv_gaussian_test.cpp | 165 ++++++++++++++++++ 5 files changed, 266 insertions(+), 49 deletions(-) diff --git a/stan/math/prim/prob/inv_gaussian_cdf.hpp b/stan/math/prim/prob/inv_gaussian_cdf.hpp index 97973700897..fea89062b60 100644 --- a/stan/math/prim/prob/inv_gaussian_cdf.hpp +++ b/stan/math/prim/prob/inv_gaussian_cdf.hpp @@ -10,8 +10,7 @@ #include #include #include -#include -#include +#include #include #include #include @@ -35,9 +34,10 @@ namespace math { * \f$\Phi(z_1) + e^{2\lambda/\mu}\Phi(-z_2)\f$ overflows a double above * \f$2\lambda/\mu = 710\f$. See inv_gaussian_lcdf. * - *

A boundary observation short-circuits the whole return: - * \f$y = 0\f$ and \f$y = \infty\f$ each take the result to their - * limit. The partials are zero at both. + *

Both boundaries are handled elementwise: \f$y = 0\f$ contributes a + * factor of zero and \f$y = \infty\f$ a factor of one, so a container mixing + * a boundary with ordinary observations multiplies to the same value as its + * elements taken one at a time. The partials are zero at both. * * @tparam T_y type of scalar * @tparam T_loc type of mean parameter @@ -82,12 +82,10 @@ inline return_type_t inv_gaussian_cdf( auto ops_partials = make_partials_propagator(y_ref, mu_ref, lambda_ref); - if (sum(promote_scalar(y_val == 0))) { - return ops_partials.build(0.0); - } - if (sum(promote_scalar(y_val == INFTY))) { - return ops_partials.build(1.0); - } + // Boundary masks; see inv_gaussian_lcdf. A y == inf element contributes a + // factor of one, i.e. a log contribution of zero. + const auto& is_inf = to_ref(y_val == INFTY); + const auto& is_bdry = to_ref((y_val == 0) || (y_val == INFTY)); constexpr bool any_ad = is_any_autodiff_v; @@ -99,7 +97,9 @@ inline return_type_t inv_gaussian_cdf( const auto& z2 = to_ref(sqrt_lambda_over_y * (y_over_mu + 1.0)); const auto& log_upper = to_ref(internal::log_scaled_upper_term(z1, z2)); - const auto& lcdf_elt = to_ref(log_sum_exp(internal::log_Phi(z1), log_upper)); + const auto& lcdf_elt + = to_ref(select(is_inf, T_partials_return(0), + log_sum_exp(internal::log_Phi(z1), log_upper))); T_partials_return cdf = exp(sum(lcdf_elt)); @@ -109,17 +109,26 @@ inline return_type_t inv_gaussian_cdf( // Each partial is masked separately: at y == inf the shape partial is // 0 / 0 on its own. The inner select covers elements whose log CDF has // saturated to -inf at an interior y. + const auto& is_underflow = to_ref(lcdf_elt == NEGATIVE_INFTY); const auto& w_dens = to_ref(cdf * exp(internal::log_phi(z1) - lcdf_elt)); const auto& w_upper = to_ref(cdf * exp(log_upper - lcdf_elt)); if constexpr (is_autodiff_v) { - partials<0>(ops_partials) = w_dens * sqrt_lambda_over_y * inv_y; + partials<0>(ops_partials) + = select(is_bdry, T_partials_return(0), + select(is_underflow, T_partials_return(0), + w_dens * sqrt_lambda_over_y * inv_y)); } if constexpr (is_autodiff_v) { - partials<1>(ops_partials) = -2.0 * lambda_val * w_upper * square(inv_mu); + partials<1>(ops_partials) + = select(is_bdry, T_partials_return(0), + select(is_underflow, T_partials_return(0), + -2.0 * lambda_val * w_upper * square(inv_mu))); } if constexpr (is_autodiff_v) { - partials<2>(ops_partials) - = 2.0 * w_upper * inv_mu - w_dens * inv_y / sqrt_lambda_over_y; + partials<2>(ops_partials) = select( + is_bdry, T_partials_return(0), + select(is_underflow, T_partials_return(0), + 2.0 * w_upper * inv_mu - w_dens * inv_y / sqrt_lambda_over_y)); } } return ops_partials.build(cdf); diff --git a/stan/math/prim/prob/inv_gaussian_lccdf.hpp b/stan/math/prim/prob/inv_gaussian_lccdf.hpp index 2f19f1d32aa..ceefe06b45f 100644 --- a/stan/math/prim/prob/inv_gaussian_lccdf.hpp +++ b/stan/math/prim/prob/inv_gaussian_lccdf.hpp @@ -10,8 +10,7 @@ #include #include #include -#include -#include +#include #include #include #include @@ -76,9 +75,11 @@ inline auto log_diff_exp_guarded(T1&& a, T2&& b) { * that first happens at \f$\log S = -5 \times 10^5\f$, so the survivor is * zero to any representable precision there. * - *

A boundary observation short-circuits the whole return: - * \f$y = 0\f$ and \f$y = \infty\f$ each take the result to their - * limit. The partials are zero at both. + *

Both boundaries are handled elementwise: \f$y = 0\f$ contributes + * \f$\log 1 = 0\f$ and \f$y = \infty\f$ contributes \f$-\infty\f$, so a + * container mixing a boundary with ordinary observations sums to the same + * value as its elements taken one at a time. The partials are zero at both, + * and wherever the survivor has underflowed. * * @tparam T_y type of scalar * @tparam T_loc type of mean parameter @@ -123,12 +124,9 @@ inline return_type_t inv_gaussian_lccdf( auto ops_partials = make_partials_propagator(y_ref, mu_ref, lambda_ref); - if (sum(promote_scalar(y_val == 0))) { - return ops_partials.build(0); - } - if (sum(promote_scalar(y_val == INFTY))) { - return ops_partials.build(NEGATIVE_INFTY); - } + // Boundary masks; see inv_gaussian_lcdf. + const auto& is_inf = to_ref(y_val == INFTY); + const auto& is_bdry = to_ref((y_val == 0) || (y_val == INFTY)); constexpr bool any_ad = is_any_autodiff_v; @@ -140,8 +138,9 @@ inline return_type_t inv_gaussian_lccdf( const auto& z2 = to_ref(sqrt_lambda_over_y * (y_over_mu + 1.0)); const auto& log_upper = to_ref(internal::log_scaled_upper_term(z1, z2)); - const auto& lccdf_elt = to_ref( - internal::log_diff_exp_guarded(internal::log_Phi(-z1), log_upper)); + const auto& lccdf_elt = to_ref(select( + is_inf, T_partials_return(NEGATIVE_INFTY), + internal::log_diff_exp_guarded(internal::log_Phi(-z1), log_upper))); T_partials_return lccdf = sum(lccdf_elt); @@ -152,17 +151,26 @@ inline return_type_t inv_gaussian_lccdf( // Each partial is masked separately: at y == inf the shape partial is // 0 / 0 on its own. The inner select covers elements whose survivor has // underflowed to -inf. + const auto& is_underflow = to_ref(lccdf_elt == NEGATIVE_INFTY); const auto& w_dens = to_ref(exp(internal::log_phi(z1) - lccdf_elt)); const auto& w_upper = to_ref(exp(log_upper - lccdf_elt)); if constexpr (is_autodiff_v) { - partials<0>(ops_partials) = -w_dens * sqrt_lambda_over_y * inv_y; + partials<0>(ops_partials) + = select(is_bdry, T_partials_return(0), + select(is_underflow, T_partials_return(0), + -w_dens * sqrt_lambda_over_y * inv_y)); } if constexpr (is_autodiff_v) { - partials<1>(ops_partials) = 2.0 * lambda_val * w_upper * square(inv_mu); + partials<1>(ops_partials) + = select(is_bdry, T_partials_return(0), + select(is_underflow, T_partials_return(0), + 2.0 * lambda_val * w_upper * square(inv_mu))); } if constexpr (is_autodiff_v) { - partials<2>(ops_partials) - = w_dens * inv_y / sqrt_lambda_over_y - 2.0 * w_upper * inv_mu; + partials<2>(ops_partials) = select( + is_bdry, T_partials_return(0), + select(is_underflow, T_partials_return(0), + w_dens * inv_y / sqrt_lambda_over_y - 2.0 * w_upper * inv_mu)); } } return ops_partials.build(lccdf); diff --git a/stan/math/prim/prob/inv_gaussian_lcdf.hpp b/stan/math/prim/prob/inv_gaussian_lcdf.hpp index 071ded53012..ebc4469f4cd 100644 --- a/stan/math/prim/prob/inv_gaussian_lcdf.hpp +++ b/stan/math/prim/prob/inv_gaussian_lcdf.hpp @@ -14,8 +14,7 @@ #include #include #include -#include -#include +#include #include #include #include @@ -144,9 +143,10 @@ inline auto log_scaled_upper_term(T1&& z1, T2&& z2) { * exact at \f$y = \mu\f$ and returns \f$z_1 = 0\f$ there for any * \f$\lambda\f$. * - *

A boundary observation short-circuits the whole return: - * \f$y = 0\f$ and \f$y = \infty\f$ each take the result to their - * limit. The partials are zero at both. + *

Both boundaries are handled elementwise: \f$y = 0\f$ contributes + * \f$-\infty\f$ and \f$y = \infty\f$ contributes \f$\log 1 = 0\f$, so a + * container mixing a boundary with ordinary observations sums to the same + * value as its elements taken one at a time. The partials are zero at both. * * @tparam T_y type of scalar * @tparam T_loc type of mean parameter @@ -191,12 +191,11 @@ inline return_type_t inv_gaussian_lcdf( auto ops_partials = make_partials_propagator(y_ref, mu_ref, lambda_ref); - if (sum(promote_scalar(y_val == 0))) { - return ops_partials.build(NEGATIVE_INFTY); - } - if (sum(promote_scalar(y_val == INFTY))) { - return ops_partials.build(0); - } + // Boundary masks. The standardized arguments are 0 * inf at y == inf, so + // that contribution is substituted; at y == 0 they reach their limits on + // their own. + const auto& is_inf = to_ref(y_val == INFTY); + const auto& is_bdry = to_ref((y_val == 0) || (y_val == INFTY)); constexpr bool any_ad = is_any_autodiff_v; @@ -208,7 +207,9 @@ inline return_type_t inv_gaussian_lcdf( const auto& z2 = to_ref(sqrt_lambda_over_y * (y_over_mu + 1.0)); const auto& log_upper = to_ref(internal::log_scaled_upper_term(z1, z2)); - const auto& lcdf_elt = to_ref(log_sum_exp(internal::log_Phi(z1), log_upper)); + const auto& lcdf_elt + = to_ref(select(is_inf, T_partials_return(0), + log_sum_exp(internal::log_Phi(z1), log_upper))); T_partials_return lcdf = sum(lcdf_elt); @@ -219,17 +220,26 @@ inline return_type_t inv_gaussian_lcdf( // Each partial is masked separately: at y == inf the shape partial is // 0 / 0 on its own. The inner select covers elements whose log CDF has // saturated to -inf at an interior y. + const auto& is_underflow = to_ref(lcdf_elt == NEGATIVE_INFTY); const auto& w_dens = to_ref(exp(internal::log_phi(z1) - lcdf_elt)); const auto& w_upper = to_ref(exp(log_upper - lcdf_elt)); if constexpr (is_autodiff_v) { - partials<0>(ops_partials) = w_dens * sqrt_lambda_over_y * inv_y; + partials<0>(ops_partials) + = select(is_bdry, T_partials_return(0), + select(is_underflow, T_partials_return(0), + w_dens * sqrt_lambda_over_y * inv_y)); } if constexpr (is_autodiff_v) { - partials<1>(ops_partials) = -2.0 * lambda_val * w_upper * square(inv_mu); + partials<1>(ops_partials) + = select(is_bdry, T_partials_return(0), + select(is_underflow, T_partials_return(0), + -2.0 * lambda_val * w_upper * square(inv_mu))); } if constexpr (is_autodiff_v) { - partials<2>(ops_partials) - = 2.0 * w_upper * inv_mu - w_dens * inv_y / sqrt_lambda_over_y; + partials<2>(ops_partials) = select( + is_bdry, T_partials_return(0), + select(is_underflow, T_partials_return(0), + 2.0 * w_upper * inv_mu - w_dens * inv_y / sqrt_lambda_over_y)); } } return ops_partials.build(lcdf); diff --git a/test/unit/math/prim/prob/inv_gaussian_test.cpp b/test/unit/math/prim/prob/inv_gaussian_test.cpp index 47e8c0230eb..fc95dbe2026 100644 --- a/test/unit/math/prim/prob/inv_gaussian_test.cpp +++ b/test/unit/math/prim/prob/inv_gaussian_test.cpp @@ -128,6 +128,31 @@ TEST(ProbDistributionsInvGaussian, boundaries) { EXPECT_FLOAT_EQ(1.0, inv_gaussian_cdf(inf, 1.0, 2.0)); } +TEST(ProbDistributionsInvGaussian, boundariesInContainer) { + using stan::math::inv_gaussian_cdf; + using stan::math::inv_gaussian_lccdf; + using stan::math::inv_gaussian_lcdf; + using stan::math::inv_gaussian_lpdf; + double inf = std::numeric_limits::infinity(); + std::vector y_inf{inf, 0.5}; + std::vector y_zero{0.0, 0.5}; + + // a factor of one / a summand of zero: the finite element carries the value + EXPECT_FLOAT_EQ(inv_gaussian_lcdf(0.5, 1.0, 2.0), + inv_gaussian_lcdf(y_inf, 1.0, 2.0)); + EXPECT_FLOAT_EQ(inv_gaussian_cdf(0.5, 1.0, 2.0), + inv_gaussian_cdf(y_inf, 1.0, 2.0)); + EXPECT_FLOAT_EQ(inv_gaussian_lccdf(0.5, 1.0, 2.0), + inv_gaussian_lccdf(y_zero, 1.0, 2.0)); + + // the boundary absorbs the reduction: the whole return is the boundary + EXPECT_FLOAT_EQ(-inf, inv_gaussian_lpdf(y_inf, 1.0, 2.0)); + EXPECT_FLOAT_EQ(-inf, inv_gaussian_lpdf(y_zero, 1.0, 2.0)); + EXPECT_FLOAT_EQ(-inf, inv_gaussian_lcdf(y_zero, 1.0, 2.0)); + EXPECT_FLOAT_EQ(0.0, inv_gaussian_cdf(y_zero, 1.0, 2.0)); + EXPECT_FLOAT_EQ(-inf, inv_gaussian_lccdf(y_inf, 1.0, 2.0)); +} + // 2 lambda / mu is past the overflow point of exp for all of these. TEST(ProbDistributionsInvGaussian, largeExpFactor) { using stan::math::inv_gaussian_lccdf; diff --git a/test/unit/math/rev/prob/inv_gaussian_test.cpp b/test/unit/math/rev/prob/inv_gaussian_test.cpp index a6c9493c827..629428b7043 100644 --- a/test/unit/math/rev/prob/inv_gaussian_test.cpp +++ b/test/unit/math/rev/prob/inv_gaussian_test.cpp @@ -159,3 +159,168 @@ TEST_F(AgradRev, inv_gaussian_vectorized_matches_scalar) { EXPECT_FLOAT_EQ(scalar_adj_y(i), vec_adj_y(i)); } } + +// Boundary adjoints are asserted directly here; expect_ad's finite-difference +// stencil steps off the support at y == 0 and is NaN at y == inf. +TEST_F(AgradRev, inv_gaussian_boundary_partials_are_zero) { + using stan::math::var; + double inf = std::numeric_limits::infinity(); + + auto check = [](auto f, double yv, const char* what) { + var y = yv; + var mu = 1.0; + var lambda = 2.0; + var out = f(y, mu, lambda); + out.grad(); + EXPECT_FLOAT_EQ(0.0, y.adj()) << what << " d/dy at y=" << yv; + EXPECT_FLOAT_EQ(0.0, mu.adj()) << what << " d/dmu at y=" << yv; + EXPECT_FLOAT_EQ(0.0, lambda.adj()) << what << " d/dlambda at y=" << yv; + stan::math::recover_memory(); + }; + + auto lpdf = [](const auto& y, const auto& mu, const auto& lambda) { + return stan::math::inv_gaussian_lpdf(y, mu, lambda); + }; + auto cdf = [](const auto& y, const auto& mu, const auto& lambda) { + return stan::math::inv_gaussian_cdf(y, mu, lambda); + }; + auto lcdf = [](const auto& y, const auto& mu, const auto& lambda) { + return stan::math::inv_gaussian_lcdf(y, mu, lambda); + }; + auto lccdf = [](const auto& y, const auto& mu, const auto& lambda) { + return stan::math::inv_gaussian_lccdf(y, mu, lambda); + }; + + for (double b : {0.0, inf}) { + check(lpdf, b, "lpdf"); + check(cdf, b, "cdf"); + check(lcdf, b, "lcdf"); + check(lccdf, b, "lccdf"); + } +} + +// A y == 0 element contributes a log survivor of zero and leaves the other +// elements alone, gradients included. +TEST_F(AgradRev, inv_gaussian_lccdf_zero_boundary_is_elementwise) { + using stan::math::var; + using stan::math::vector_v; + + EXPECT_FLOAT_EQ(0.0, stan::math::inv_gaussian_lccdf(0.0, 1.0, 2.0)); + + std::vector y{0.0, 5.0}; + EXPECT_FLOAT_EQ(stan::math::inv_gaussian_lccdf(5.0, 1.0, 2.0), + stan::math::inv_gaussian_lccdf(y, 1.0, 2.0)); + + vector_v y_v(2); + y_v << 0.0, 5.0; + var mu = 1.0; + var lambda = 2.0; + var out = stan::math::inv_gaussian_lccdf(y_v, mu, lambda); + out.grad(); + double adj_boundary = y_v(0).adj(); + double adj_finite = y_v(1).adj(); + double adj_mu = mu.adj(); + double adj_lambda = lambda.adj(); + stan::math::recover_memory(); + + var y_s = 5.0; + var mu_s = 1.0; + var lambda_s = 2.0; + var out_s = stan::math::inv_gaussian_lccdf(y_s, mu_s, lambda_s); + out_s.grad(); + EXPECT_FLOAT_EQ(0.0, adj_boundary); + EXPECT_FLOAT_EQ(y_s.adj(), adj_finite); + EXPECT_FLOAT_EQ(mu_s.adj(), adj_mu); + EXPECT_FLOAT_EQ(lambda_s.adj(), adj_lambda); + stan::math::recover_memory(); +} + +// Deep in the upper tail the survivor's two terms agree to the last bit and +// the guarded difference underflows to -inf at a finite y. The partials are +// zero there, and the guard fires deterministically at these points. +TEST_F(AgradRev, inv_gaussian_lccdf_underflow_partials_are_zero) { + using stan::math::var; + double inf = std::numeric_limits::infinity(); + for (double y_dbl : {1e-2, 1.0, 10.0}) { + var y = y_dbl; + var mu = 1e-3; + var lambda = 1e14; + var out = stan::math::inv_gaussian_lccdf(y, mu, lambda); + EXPECT_FLOAT_EQ(-inf, out.val()) << "y=" << y_dbl; + out.grad(); + EXPECT_FLOAT_EQ(0.0, y.adj()) << "y=" << y_dbl; + EXPECT_FLOAT_EQ(0.0, mu.adj()) << "y=" << y_dbl; + EXPECT_FLOAT_EQ(0.0, lambda.adj()) << "y=" << y_dbl; + stan::math::recover_memory(); + } +} + +// lambda / y overflows for representable arguments, saturating the log CDF to +// -inf at an interior y. The saturation mask zeroes every partial there. +TEST_F(AgradRev, inv_gaussian_saturation_partials_are_zero) { + using stan::math::var; + double inf = std::numeric_limits::infinity(); + + auto check = [](auto f, double y_dbl, double mu_dbl, double lambda_dbl, + double val, const char* what) { + var y = y_dbl; + var mu = mu_dbl; + var lambda = lambda_dbl; + var out = f(y, mu, lambda); + EXPECT_FLOAT_EQ(val, out.val()) << what; + out.grad(); + EXPECT_FLOAT_EQ(0.0, y.adj()) << what; + EXPECT_FLOAT_EQ(0.0, mu.adj()) << what; + EXPECT_FLOAT_EQ(0.0, lambda.adj()) << what; + stan::math::recover_memory(); + }; + + auto lcdf = [](const auto& y, const auto& mu, const auto& lambda) { + return stan::math::inv_gaussian_lcdf(y, mu, lambda); + }; + auto cdf = [](const auto& y, const auto& mu, const auto& lambda) { + return stan::math::inv_gaussian_cdf(y, mu, lambda); + }; + auto lccdf = [](const auto& y, const auto& mu, const auto& lambda) { + return stan::math::inv_gaussian_lccdf(y, mu, lambda); + }; + + // lambda / y = 1e310 overflows; y is far below mu, so the lower tail + // saturates + check(lcdf, 1e-10, 1.0, 1e300, -inf, "lcdf"); + check(cdf, 1e-10, 1.0, 1e300, 0.0, "cdf"); + // y is above mu, so the survivor saturates + check(lccdf, 2e-3, 1e-3, 1e308, -inf, "lccdf"); +} + +// The same for the y == inf boundary of the lcdf and cdf. +TEST_F(AgradRev, inv_gaussian_inf_boundary_is_elementwise) { + using stan::math::var; + using stan::math::vector_v; + double inf = std::numeric_limits::infinity(); + + vector_v y_v(2); + y_v << inf, 0.5; + var mu = 1.0; + var lambda = 2.0; + var out = stan::math::inv_gaussian_lcdf(y_v, mu, lambda); + out.grad(); + double vec_val = out.val(); + double adj_boundary = y_v(0).adj(); + double adj_finite = y_v(1).adj(); + double adj_mu = mu.adj(); + double adj_lambda = lambda.adj(); + stan::math::recover_memory(); + + var y_s = 0.5; + var mu_s = 1.0; + var lambda_s = 2.0; + var out_s = stan::math::inv_gaussian_lcdf(y_s, mu_s, lambda_s); + out_s.grad(); + EXPECT_FLOAT_EQ(out_s.val(), vec_val); + EXPECT_FLOAT_EQ(0.0, adj_boundary); + EXPECT_FLOAT_EQ(y_s.adj(), adj_finite); + EXPECT_FLOAT_EQ(mu_s.adj(), adj_mu); + EXPECT_FLOAT_EQ(lambda_s.adj(), adj_lambda); + stan::math::recover_memory(); +} From 07dd10bc58ee4de4bd2a0ccd33835bb62fa751b0 Mon Sep 17 00:00:00 2001 From: Gidon Frischkorn Date: Sun, 23 Aug 2026 16:47:10 +0200 Subject: [PATCH 07/11] chore: rename internal::log_phi to log_std_normal_density --- stan/math/prim/prob/inv_gaussian_cdf.hpp | 3 ++- stan/math/prim/prob/inv_gaussian_lccdf.hpp | 3 ++- stan/math/prim/prob/inv_gaussian_lcdf.hpp | 18 ++++++++++-------- 3 files changed, 14 insertions(+), 10 deletions(-) diff --git a/stan/math/prim/prob/inv_gaussian_cdf.hpp b/stan/math/prim/prob/inv_gaussian_cdf.hpp index fea89062b60..9e0a1ae79aa 100644 --- a/stan/math/prim/prob/inv_gaussian_cdf.hpp +++ b/stan/math/prim/prob/inv_gaussian_cdf.hpp @@ -110,7 +110,8 @@ inline return_type_t inv_gaussian_cdf( // 0 / 0 on its own. The inner select covers elements whose log CDF has // saturated to -inf at an interior y. const auto& is_underflow = to_ref(lcdf_elt == NEGATIVE_INFTY); - const auto& w_dens = to_ref(cdf * exp(internal::log_phi(z1) - lcdf_elt)); + const auto& w_dens + = to_ref(cdf * exp(internal::log_std_normal_density(z1) - lcdf_elt)); const auto& w_upper = to_ref(cdf * exp(log_upper - lcdf_elt)); if constexpr (is_autodiff_v) { partials<0>(ops_partials) diff --git a/stan/math/prim/prob/inv_gaussian_lccdf.hpp b/stan/math/prim/prob/inv_gaussian_lccdf.hpp index ceefe06b45f..88333efd264 100644 --- a/stan/math/prim/prob/inv_gaussian_lccdf.hpp +++ b/stan/math/prim/prob/inv_gaussian_lccdf.hpp @@ -152,7 +152,8 @@ inline return_type_t inv_gaussian_lccdf( // 0 / 0 on its own. The inner select covers elements whose survivor has // underflowed to -inf. const auto& is_underflow = to_ref(lccdf_elt == NEGATIVE_INFTY); - const auto& w_dens = to_ref(exp(internal::log_phi(z1) - lccdf_elt)); + const auto& w_dens + = to_ref(exp(internal::log_std_normal_density(z1) - lccdf_elt)); const auto& w_upper = to_ref(exp(log_upper - lccdf_elt)); if constexpr (is_autodiff_v) { partials<0>(ops_partials) diff --git a/stan/math/prim/prob/inv_gaussian_lcdf.hpp b/stan/math/prim/prob/inv_gaussian_lcdf.hpp index ebc4469f4cd..85ecc50b348 100644 --- a/stan/math/prim/prob/inv_gaussian_lcdf.hpp +++ b/stan/math/prim/prob/inv_gaussian_lcdf.hpp @@ -74,23 +74,24 @@ inline auto log_Phi(T&& z) { * Return the log of the standard normal density at the specified value. */ template * = nullptr> -inline return_type_t log_phi(const T& z) { +inline return_type_t log_std_normal_density(const T& z) { return -0.5 * square(z) - HALF_LOG_TWO_PI; } -struct log_phi_fun { +struct log_std_normal_density_fun { template static inline auto fun(T&& z) { - return log_phi(std::forward(z)); + return log_std_normal_density(std::forward(z)); } }; /** - * A vectorized version of log_phi(). + * A vectorized version of log_std_normal_density(). */ template * = nullptr> -inline auto log_phi(T&& z) { - return apply_scalar_unary::apply(std::forward(z)); +inline auto log_std_normal_density(T&& z) { + return apply_scalar_unary::apply( + std::forward(z)); } /** @@ -109,7 +110,7 @@ inline return_type_t log_scaled_upper_term(const T1& z1, const T2& z2) { using std::log; if (value_of_rec(z2) > -LOG_PHI_ASYMPTOTIC_CUTOFF) { const auto s = inv_square(z2); - return log_phi(z1) - log(z2) + return log_std_normal_density(z1) - log(z2) + log1p(s * (-1.0 + s * (3.0 + s * (-15.0 + s * 105.0)))); } return 0.5 * (square(z2) - square(z1)) + log_Phi(-z2); @@ -221,7 +222,8 @@ inline return_type_t inv_gaussian_lcdf( // 0 / 0 on its own. The inner select covers elements whose log CDF has // saturated to -inf at an interior y. const auto& is_underflow = to_ref(lcdf_elt == NEGATIVE_INFTY); - const auto& w_dens = to_ref(exp(internal::log_phi(z1) - lcdf_elt)); + const auto& w_dens + = to_ref(exp(internal::log_std_normal_density(z1) - lcdf_elt)); const auto& w_upper = to_ref(exp(log_upper - lcdf_elt)); if constexpr (is_autodiff_v) { partials<0>(ops_partials) From 70c24edf1097e6c918b59d2fb90b2c8dab1c46fa Mon Sep 17 00:00:00 2001 From: Gidon Frischkorn Date: Sun, 23 Aug 2026 18:12:36 +0200 Subject: [PATCH 08/11] test: cover vectorized cdf/lcdf gradients and std::vector inputs --- test/unit/math/mix/prob/inv_gaussian_test.cpp | 53 +++++++++++++++++++ 1 file changed, 53 insertions(+) diff --git a/test/unit/math/mix/prob/inv_gaussian_test.cpp b/test/unit/math/mix/prob/inv_gaussian_test.cpp index c9a6fb49150..3b4aed3d173 100644 --- a/test/unit/math/mix/prob/inv_gaussian_test.cpp +++ b/test/unit/math/mix/prob/inv_gaussian_test.cpp @@ -108,3 +108,56 @@ TEST_F(AgradRev, mathMixScalFun_inv_gaussian_lccdf_vectorized) { }; stan::test::expect_ad(f, y, mu, lambda); } + +// The cdf's partials scale by the whole-container product, which a scalar +// test does not exercise. +TEST_F(AgradRev, mathMixScalFun_inv_gaussian_cdf_vectorized) { + Eigen::VectorXd y(3); + y << 0.3, 1.2, 4.0; + Eigen::VectorXd mu(3); + mu << 0.5, 1.0, 2.0; + Eigen::VectorXd lambda(3); + lambda << 2.0, 5.0, 0.7; + + auto f = [](const auto& y, const auto& mu, const auto& lambda) { + return stan::math::inv_gaussian_cdf(y, mu, lambda); + }; + stan::test::expect_ad(f, y, mu, lambda); +} + +TEST_F(AgradRev, mathMixScalFun_inv_gaussian_lcdf_vectorized) { + Eigen::VectorXd y(3); + y << 0.3, 1.2, 4.0; + Eigen::VectorXd mu(3); + mu << 0.5, 1.0, 2.0; + Eigen::VectorXd lambda(3); + lambda << 2.0, 5.0, 0.7; + + auto f = [](const auto& y, const auto& mu, const auto& lambda) { + return stan::math::inv_gaussian_lcdf(y, mu, lambda); + }; + stan::test::expect_ad(f, y, mu, lambda); +} + +TEST_F(AgradRev, mathMixScalFun_inv_gaussian_std_vector) { + std::vector y{0.3, 1.2, 4.0}; + std::vector mu{0.5, 1.0, 2.0}; + std::vector lambda{2.0, 5.0, 0.7}; + + auto f_lpdf = [](const auto& y, const auto& mu, const auto& lambda) { + return stan::math::inv_gaussian_lpdf(y, mu, lambda); + }; + auto f_cdf = [](const auto& y, const auto& mu, const auto& lambda) { + return stan::math::inv_gaussian_cdf(y, mu, lambda); + }; + auto f_lcdf = [](const auto& y, const auto& mu, const auto& lambda) { + return stan::math::inv_gaussian_lcdf(y, mu, lambda); + }; + auto f_lccdf = [](const auto& y, const auto& mu, const auto& lambda) { + return stan::math::inv_gaussian_lccdf(y, mu, lambda); + }; + stan::test::expect_ad(f_lpdf, y, mu, lambda); + stan::test::expect_ad(f_cdf, y, mu, lambda); + stan::test::expect_ad(f_lcdf, y, mu, lambda); + stan::test::expect_ad(f_lccdf, y, mu, lambda); +} From c082e3986b3f42e20c9db696a5ef3e517415531a Mon Sep 17 00:00:00 2001 From: Gidon Frischkorn Date: Mon, 24 Aug 2026 17:05:44 +0200 Subject: [PATCH 09/11] docs: tighten inv_gaussian comments --- stan/math/prim/prob/inv_gaussian_cdf.hpp | 6 ++---- stan/math/prim/prob/inv_gaussian_lccdf.hpp | 7 ++----- stan/math/prim/prob/inv_gaussian_lcdf.hpp | 9 +++------ stan/math/prim/prob/inv_gaussian_lpdf.hpp | 1 - stan/math/prim/prob/inv_gaussian_rng.hpp | 8 +++----- test/unit/math/prim/prob/inv_gaussian_test.cpp | 3 --- 6 files changed, 10 insertions(+), 24 deletions(-) diff --git a/stan/math/prim/prob/inv_gaussian_cdf.hpp b/stan/math/prim/prob/inv_gaussian_cdf.hpp index 9e0a1ae79aa..e4b1cddd60d 100644 --- a/stan/math/prim/prob/inv_gaussian_cdf.hpp +++ b/stan/math/prim/prob/inv_gaussian_cdf.hpp @@ -34,10 +34,8 @@ namespace math { * \f$\Phi(z_1) + e^{2\lambda/\mu}\Phi(-z_2)\f$ overflows a double above * \f$2\lambda/\mu = 710\f$. See inv_gaussian_lcdf. * - *

Both boundaries are handled elementwise: \f$y = 0\f$ contributes a - * factor of zero and \f$y = \infty\f$ a factor of one, so a container mixing - * a boundary with ordinary observations multiplies to the same value as its - * elements taken one at a time. The partials are zero at both. + *

Both boundaries of the support are handled elementwise; the partials + * are zero at both. * * @tparam T_y type of scalar * @tparam T_loc type of mean parameter diff --git a/stan/math/prim/prob/inv_gaussian_lccdf.hpp b/stan/math/prim/prob/inv_gaussian_lccdf.hpp index 88333efd264..829fd8813d8 100644 --- a/stan/math/prim/prob/inv_gaussian_lccdf.hpp +++ b/stan/math/prim/prob/inv_gaussian_lccdf.hpp @@ -75,11 +75,8 @@ inline auto log_diff_exp_guarded(T1&& a, T2&& b) { * that first happens at \f$\log S = -5 \times 10^5\f$, so the survivor is * zero to any representable precision there. * - *

Both boundaries are handled elementwise: \f$y = 0\f$ contributes - * \f$\log 1 = 0\f$ and \f$y = \infty\f$ contributes \f$-\infty\f$, so a - * container mixing a boundary with ordinary observations sums to the same - * value as its elements taken one at a time. The partials are zero at both, - * and wherever the survivor has underflowed. + *

Both boundaries of the support are handled elementwise; the partials + * are zero at both, and wherever the survivor has underflowed. * * @tparam T_y type of scalar * @tparam T_loc type of mean parameter diff --git a/stan/math/prim/prob/inv_gaussian_lcdf.hpp b/stan/math/prim/prob/inv_gaussian_lcdf.hpp index 85ecc50b348..2ed9bf19bed 100644 --- a/stan/math/prim/prob/inv_gaussian_lcdf.hpp +++ b/stan/math/prim/prob/inv_gaussian_lcdf.hpp @@ -141,13 +141,10 @@ inline auto log_scaled_upper_term(T1&& z1, T2&& z2) { * \f$z_{1,2} = \sqrt{\lambda/y}\,(y/\mu \mp 1)\f$. The factor * \f$e^{2\lambda/\mu}\f$ overflows a double above \f$2\lambda/\mu = 710\f$, * so the two terms are combined in log space. That grouping of \f$z\f$ is - * exact at \f$y = \mu\f$ and returns \f$z_1 = 0\f$ there for any - * \f$\lambda\f$. + * exact in floating point at \f$y = \mu\f$. * - *

Both boundaries are handled elementwise: \f$y = 0\f$ contributes - * \f$-\infty\f$ and \f$y = \infty\f$ contributes \f$\log 1 = 0\f$, so a - * container mixing a boundary with ordinary observations sums to the same - * value as its elements taken one at a time. The partials are zero at both. + *

Both boundaries of the support are handled elementwise; the partials + * are zero at both. * * @tparam T_y type of scalar * @tparam T_loc type of mean parameter diff --git a/stan/math/prim/prob/inv_gaussian_lpdf.hpp b/stan/math/prim/prob/inv_gaussian_lpdf.hpp index 86eb5f2e477..a2798ed27fa 100644 --- a/stan/math/prim/prob/inv_gaussian_lpdf.hpp +++ b/stan/math/prim/prob/inv_gaussian_lpdf.hpp @@ -98,7 +98,6 @@ inline return_type_t inv_gaussian_lpdf( const auto& inv_mu = to_ref_if(inv(mu_val)); const auto& inv_y = to_ref_if(inv(y_val)); const auto& y_m_mu = to_ref_if(y_val - mu_val); - // (y - mu)^2 / (2 mu^2 y) const auto& half_sq_scaled = to_ref(0.5 * square(y_m_mu * inv_mu) * inv_y); size_t N = max_size(y, mu, lambda); diff --git a/stan/math/prim/prob/inv_gaussian_rng.hpp b/stan/math/prim/prob/inv_gaussian_rng.hpp index a2a7810369a..0edfcd17c5d 100644 --- a/stan/math/prim/prob/inv_gaussian_rng.hpp +++ b/stan/math/prim/prob/inv_gaussian_rng.hpp @@ -30,11 +30,9 @@ namespace math { * candidate roots and a Bernoulli draw selects between them. The method is * exact, with no rejection step. * - *

The smaller root is formed as - * \f$x = \mu / (1 + u/2 + \sqrt{u + u^2/4})\f$ with - * \f$u = \mu w / \lambda\f$. Since \f$(1 + u/2)^2 - (u + u^2/4) = 1\f$ - * exactly, this form subtracts nothing and stays accurate for \f$u\f$ up to - * \f$10^{20}\f$. + *

The smaller root is computed in reciprocal form; since + * \f$(1 + u/2)^2 - (u + u^2/4) = 1\f$ exactly, it subtracts nothing and + * stays accurate for \f$u = \mu w / \lambda\f$ up to \f$10^{20}\f$. * * @tparam T_loc type of mean parameter * @tparam T_shape type of shape parameter diff --git a/test/unit/math/prim/prob/inv_gaussian_test.cpp b/test/unit/math/prim/prob/inv_gaussian_test.cpp index fc95dbe2026..ea46c4802ab 100644 --- a/test/unit/math/prim/prob/inv_gaussian_test.cpp +++ b/test/unit/math/prim/prob/inv_gaussian_test.cpp @@ -121,7 +121,6 @@ TEST(ProbDistributionsInvGaussian, boundaries) { EXPECT_FLOAT_EQ(-inf, inv_gaussian_lcdf(0.0, 1.0, 2.0)); EXPECT_FLOAT_EQ(0.0, inv_gaussian_lccdf(0.0, 1.0, 2.0)); EXPECT_FLOAT_EQ(0.0, inv_gaussian_cdf(0.0, 1.0, 2.0)); - // y == inf is in the closure of the support: zero density, saturated CDF EXPECT_FLOAT_EQ(-inf, inv_gaussian_lpdf(inf, 1.0, 2.0)); EXPECT_FLOAT_EQ(0.0, inv_gaussian_lcdf(inf, 1.0, 2.0)); EXPECT_FLOAT_EQ(-inf, inv_gaussian_lccdf(inf, 1.0, 2.0)); @@ -206,8 +205,6 @@ TEST(ProbDistributionsInvGaussian, probabilityNeverExceedsOne) { EXPECT_FLOAT_EQ(-4.99000499999999997e17, inv_gaussian_lcdf(1e-4, 1e-1, 1e14)); } -// At y == mu the CDF is a fixed quantity independent of the scale of mu, so -// all of these must agree. TEST(ProbDistributionsInvGaussian, medianIsScaleInvariant) { using stan::math::inv_gaussian_lccdf; using stan::math::inv_gaussian_lcdf; From 5e6f41b6afcc5d9fe7b25b7e1011b73cf7569133 Mon Sep 17 00:00:00 2001 From: Gidon Frischkorn Date: Tue, 25 Aug 2026 10:31:27 +0200 Subject: [PATCH 10/11] tests: drop coverage duplicated across suites - value pins and error throws already covered by the test/prob fixtures - Eigen-vectorized expect_ad covered by the generated vv/ffv variants --- test/unit/math/mix/prob/inv_gaussian_test.cpp | 42 ------ .../unit/math/prim/prob/inv_gaussian_test.cpp | 72 +--------- test/unit/math/rev/prob/inv_gaussian_test.cpp | 131 ++++++------------ 3 files changed, 50 insertions(+), 195 deletions(-) diff --git a/test/unit/math/mix/prob/inv_gaussian_test.cpp b/test/unit/math/mix/prob/inv_gaussian_test.cpp index 3b4aed3d173..f9f30f3d807 100644 --- a/test/unit/math/mix/prob/inv_gaussian_test.cpp +++ b/test/unit/math/mix/prob/inv_gaussian_test.cpp @@ -81,34 +81,6 @@ TEST_F(AgradRev, mathMixScalFun_inv_gaussian_lccdf_tails) { stan::test::expect_ad(f, 0.5, 0.1, 50.0); } -TEST_F(AgradRev, mathMixScalFun_inv_gaussian_lpdf_vectorized) { - Eigen::VectorXd y(3); - y << 0.3, 1.2, 4.0; - Eigen::VectorXd mu(3); - mu << 0.5, 1.0, 2.0; - Eigen::VectorXd lambda(3); - lambda << 2.0, 5.0, 0.7; - - auto f = [](const auto& y, const auto& mu, const auto& lambda) { - return stan::math::inv_gaussian_lpdf(y, mu, lambda); - }; - stan::test::expect_ad(f, y, mu, lambda); -} - -TEST_F(AgradRev, mathMixScalFun_inv_gaussian_lccdf_vectorized) { - Eigen::VectorXd y(3); - y << 0.3, 1.2, 4.0; - Eigen::VectorXd mu(3); - mu << 0.5, 1.0, 2.0; - Eigen::VectorXd lambda(3); - lambda << 2.0, 5.0, 0.7; - - auto f = [](const auto& y, const auto& mu, const auto& lambda) { - return stan::math::inv_gaussian_lccdf(y, mu, lambda); - }; - stan::test::expect_ad(f, y, mu, lambda); -} - // The cdf's partials scale by the whole-container product, which a scalar // test does not exercise. TEST_F(AgradRev, mathMixScalFun_inv_gaussian_cdf_vectorized) { @@ -125,20 +97,6 @@ TEST_F(AgradRev, mathMixScalFun_inv_gaussian_cdf_vectorized) { stan::test::expect_ad(f, y, mu, lambda); } -TEST_F(AgradRev, mathMixScalFun_inv_gaussian_lcdf_vectorized) { - Eigen::VectorXd y(3); - y << 0.3, 1.2, 4.0; - Eigen::VectorXd mu(3); - mu << 0.5, 1.0, 2.0; - Eigen::VectorXd lambda(3); - lambda << 2.0, 5.0, 0.7; - - auto f = [](const auto& y, const auto& mu, const auto& lambda) { - return stan::math::inv_gaussian_lcdf(y, mu, lambda); - }; - stan::test::expect_ad(f, y, mu, lambda); -} - TEST_F(AgradRev, mathMixScalFun_inv_gaussian_std_vector) { std::vector y{0.3, 1.2, 4.0}; std::vector mu{0.5, 1.0, 2.0}; diff --git a/test/unit/math/prim/prob/inv_gaussian_test.cpp b/test/unit/math/prim/prob/inv_gaussian_test.cpp index ea46c4802ab..cf06170e293 100644 --- a/test/unit/math/prim/prob/inv_gaussian_test.cpp +++ b/test/unit/math/prim/prob/inv_gaussian_test.cpp @@ -43,21 +43,6 @@ TEST(ProbDistributionsInvGaussian, distributionCheck) { check_quantiles_real_real(InvGaussianTestRig()); } -TEST(ProbDistributionsInvGaussian, error_check) { - boost::random::mt19937 rng; - EXPECT_NO_THROW(stan::math::inv_gaussian_rng(1.0, 2.0, rng)); - EXPECT_THROW(stan::math::inv_gaussian_rng(0.0, 2.0, rng), std::domain_error); - EXPECT_THROW(stan::math::inv_gaussian_rng(-1.0, 2.0, rng), std::domain_error); - EXPECT_THROW(stan::math::inv_gaussian_rng(1.0, 0.0, rng), std::domain_error); - EXPECT_THROW(stan::math::inv_gaussian_rng(1.0, -1.0, rng), std::domain_error); - EXPECT_THROW( - stan::math::inv_gaussian_rng(stan::math::positive_infinity(), 2.0, rng), - std::domain_error); - EXPECT_THROW( - stan::math::inv_gaussian_rng(1.0, stan::math::positive_infinity(), rng), - std::domain_error); -} - TEST(ProbDistributionsInvGaussian, rngStableForLargeMuOverLambda) { boost::random::mt19937 rng(1234); for (double mu : {1.0, 1e3, 1e6, 1e9, 1e12}) { @@ -91,25 +76,6 @@ TEST(ProbDistributionsInvGaussian, rngMomentsAtLargeMu) { EXPECT_NEAR(mu * mu * mu / lambda, var, 0.15 * mu * mu * mu / lambda); } -// reference values from mpmath at 60 digits - -TEST(ProbDistributionsInvGaussian, values) { - using stan::math::inv_gaussian_cdf; - using stan::math::inv_gaussian_lccdf; - using stan::math::inv_gaussian_lcdf; - using stan::math::inv_gaussian_lpdf; - - EXPECT_FLOAT_EQ(-2.479180611448965157415, inv_gaussian_lpdf(1.2, 0.5, 2.0)); - EXPECT_FLOAT_EQ(-2.391593703832052124008, inv_gaussian_lpdf(0.3, 1.0, 5.0)); - EXPECT_FLOAT_EQ(0.9815922531042920910742, inv_gaussian_cdf(1.2, 0.5, 2.0)); - EXPECT_FLOAT_EQ(0.003359190912064955560317, inv_gaussian_cdf(0.3, 1.0, 5.0)); - EXPECT_FLOAT_EQ(-0.01857927772712011847827, inv_gaussian_lcdf(1.2, 0.5, 2.0)); - EXPECT_FLOAT_EQ(-5.696055133984662592139, inv_gaussian_lcdf(0.3, 1.0, 5.0)); - EXPECT_FLOAT_EQ(-3.9949836760335225255, inv_gaussian_lccdf(1.2, 0.5, 2.0)); - EXPECT_FLOAT_EQ(-0.003364845660995599504277, - inv_gaussian_lccdf(0.3, 1.0, 5.0)); -} - TEST(ProbDistributionsInvGaussian, boundaries) { using stan::math::inv_gaussian_cdf; using stan::math::inv_gaussian_lccdf; @@ -153,6 +119,7 @@ TEST(ProbDistributionsInvGaussian, boundariesInContainer) { } // 2 lambda / mu is past the overflow point of exp for all of these. +// reference values from mpmath at 60 digits TEST(ProbDistributionsInvGaussian, largeExpFactor) { using stan::math::inv_gaussian_lccdf; using stan::math::inv_gaussian_lcdf; @@ -248,24 +215,6 @@ TEST(ProbDistributionsInvGaussian, cdfCcdfSumToOne) { } } -TEST(ProbDistributionsInvGaussian, vectorMatchesScalarSum) { - using stan::math::inv_gaussian_lccdf; - using stan::math::inv_gaussian_lcdf; - using stan::math::inv_gaussian_lpdf; - std::vector y{0.2, 0.5, 1.0, 2.0, 5.0}; - double sum_lpdf = 0; - double sum_lcdf = 0; - double sum_lccdf = 0; - for (double yi : y) { - sum_lpdf += inv_gaussian_lpdf(yi, 1.3, 3.0); - sum_lcdf += inv_gaussian_lcdf(yi, 1.3, 3.0); - sum_lccdf += inv_gaussian_lccdf(yi, 1.3, 3.0); - } - EXPECT_FLOAT_EQ(sum_lpdf, inv_gaussian_lpdf(y, 1.3, 3.0)); - EXPECT_FLOAT_EQ(sum_lcdf, inv_gaussian_lcdf(y, 1.3, 3.0)); - EXPECT_FLOAT_EQ(sum_lccdf, inv_gaussian_lccdf(y, 1.3, 3.0)); -} - // check helper functions; tolerances scale with each value's magnitude TEST(ProbDistributionsInvGaussian, internalLogPhi) { using stan::math::internal::log_Phi; @@ -294,29 +243,16 @@ TEST(ProbDistributionsInvGaussian, internalLogPhi) { EXPECT_FLOAT_EQ(0.0, log_Phi(inf)); EXPECT_FLOAT_EQ(-inf, log_Phi(-inf)); EXPECT_TRUE(std::isnan(log_Phi(std::numeric_limits::quiet_NaN()))); -} -// The true slope at z = -30 is about 30, so across a 2e-10 interval the -// honest change is about 6e-9. -TEST(ProbDistributionsInvGaussian, internalLogPhiBranchContinuity) { - using stan::math::internal::log_Phi; + // branch continuity: the true slope at z = -30 is about 30, so across a + // 2e-10 interval the honest change is about 6e-9 double eps = 1e-10; double step = log_Phi(-30.0 + eps) - log_Phi(-30.0 - eps); EXPECT_LT(std::fabs(step), 1e-7); } -TEST(ProbDistributionsInvGaussian, errors) { +TEST(ProbDistributionsInvGaussian, sizeMismatch) { using stan::math::inv_gaussian_lpdf; - double inf = std::numeric_limits::infinity(); - - EXPECT_THROW(inv_gaussian_lpdf(-1.0, 1.0, 2.0), std::domain_error); - EXPECT_THROW(inv_gaussian_lpdf(1.0, 0.0, 2.0), std::domain_error); - EXPECT_THROW(inv_gaussian_lpdf(1.0, -1.0, 2.0), std::domain_error); - EXPECT_THROW(inv_gaussian_lpdf(1.0, inf, 2.0), std::domain_error); - EXPECT_THROW(inv_gaussian_lpdf(1.0, 1.0, 0.0), std::domain_error); - EXPECT_THROW(inv_gaussian_lpdf(1.0, 1.0, -1.0), std::domain_error); - EXPECT_THROW(inv_gaussian_lpdf(1.0, 1.0, inf), std::domain_error); - std::vector y{1.0, 2.0}; std::vector mu{1.0, 2.0, 3.0}; EXPECT_THROW(inv_gaussian_lpdf(y, mu, 1.0), std::invalid_argument); diff --git a/test/unit/math/rev/prob/inv_gaussian_test.cpp b/test/unit/math/rev/prob/inv_gaussian_test.cpp index 629428b7043..40d5343d114 100644 --- a/test/unit/math/rev/prob/inv_gaussian_test.cpp +++ b/test/unit/math/rev/prob/inv_gaussian_test.cpp @@ -3,7 +3,6 @@ #include #include #include -#include namespace { @@ -199,65 +198,11 @@ TEST_F(AgradRev, inv_gaussian_boundary_partials_are_zero) { } } -// A y == 0 element contributes a log survivor of zero and leaves the other -// elements alone, gradients included. -TEST_F(AgradRev, inv_gaussian_lccdf_zero_boundary_is_elementwise) { - using stan::math::var; - using stan::math::vector_v; - - EXPECT_FLOAT_EQ(0.0, stan::math::inv_gaussian_lccdf(0.0, 1.0, 2.0)); - - std::vector y{0.0, 5.0}; - EXPECT_FLOAT_EQ(stan::math::inv_gaussian_lccdf(5.0, 1.0, 2.0), - stan::math::inv_gaussian_lccdf(y, 1.0, 2.0)); - - vector_v y_v(2); - y_v << 0.0, 5.0; - var mu = 1.0; - var lambda = 2.0; - var out = stan::math::inv_gaussian_lccdf(y_v, mu, lambda); - out.grad(); - double adj_boundary = y_v(0).adj(); - double adj_finite = y_v(1).adj(); - double adj_mu = mu.adj(); - double adj_lambda = lambda.adj(); - stan::math::recover_memory(); - - var y_s = 5.0; - var mu_s = 1.0; - var lambda_s = 2.0; - var out_s = stan::math::inv_gaussian_lccdf(y_s, mu_s, lambda_s); - out_s.grad(); - EXPECT_FLOAT_EQ(0.0, adj_boundary); - EXPECT_FLOAT_EQ(y_s.adj(), adj_finite); - EXPECT_FLOAT_EQ(mu_s.adj(), adj_mu); - EXPECT_FLOAT_EQ(lambda_s.adj(), adj_lambda); - stan::math::recover_memory(); -} - -// Deep in the upper tail the survivor's two terms agree to the last bit and -// the guarded difference underflows to -inf at a finite y. The partials are -// zero there, and the guard fires deterministically at these points. -TEST_F(AgradRev, inv_gaussian_lccdf_underflow_partials_are_zero) { - using stan::math::var; - double inf = std::numeric_limits::infinity(); - for (double y_dbl : {1e-2, 1.0, 10.0}) { - var y = y_dbl; - var mu = 1e-3; - var lambda = 1e14; - var out = stan::math::inv_gaussian_lccdf(y, mu, lambda); - EXPECT_FLOAT_EQ(-inf, out.val()) << "y=" << y_dbl; - out.grad(); - EXPECT_FLOAT_EQ(0.0, y.adj()) << "y=" << y_dbl; - EXPECT_FLOAT_EQ(0.0, mu.adj()) << "y=" << y_dbl; - EXPECT_FLOAT_EQ(0.0, lambda.adj()) << "y=" << y_dbl; - stan::math::recover_memory(); - } -} - -// lambda / y overflows for representable arguments, saturating the log CDF to -// -inf at an interior y. The saturation mask zeroes every partial there. -TEST_F(AgradRev, inv_gaussian_saturation_partials_are_zero) { +// The log probability saturates to -inf at an interior y in two ways: the +// guarded survivor difference underflows deep in the upper tail, and +// lambda / y overflows for representable arguments. Both trigger the same +// mask, which zeroes every partial. +TEST_F(AgradRev, inv_gaussian_saturated_partials_are_zero) { using stan::math::var; double inf = std::numeric_limits::infinity(); @@ -291,36 +236,52 @@ TEST_F(AgradRev, inv_gaussian_saturation_partials_are_zero) { check(cdf, 1e-10, 1.0, 1e300, 0.0, "cdf"); // y is above mu, so the survivor saturates check(lccdf, 2e-3, 1e-3, 1e308, -inf, "lccdf"); + // survivor underflow: the two lccdf terms agree to the last bit + check(lccdf, 1e-2, 1e-3, 1e14, -inf, "lccdf underflow"); + check(lccdf, 1.0, 1e-3, 1e14, -inf, "lccdf underflow"); + check(lccdf, 10.0, 1e-3, 1e14, -inf, "lccdf underflow"); } -// The same for the y == inf boundary of the lcdf and cdf. -TEST_F(AgradRev, inv_gaussian_inf_boundary_is_elementwise) { +// A boundary element leaves the other elements' adjoints alone; the +// container values themselves are pinned in the prim suite. +TEST_F(AgradRev, inv_gaussian_boundary_adjoints_are_elementwise) { using stan::math::var; using stan::math::vector_v; double inf = std::numeric_limits::infinity(); - vector_v y_v(2); - y_v << inf, 0.5; - var mu = 1.0; - var lambda = 2.0; - var out = stan::math::inv_gaussian_lcdf(y_v, mu, lambda); - out.grad(); - double vec_val = out.val(); - double adj_boundary = y_v(0).adj(); - double adj_finite = y_v(1).adj(); - double adj_mu = mu.adj(); - double adj_lambda = lambda.adj(); - stan::math::recover_memory(); + auto check = [](auto f, double boundary, double finite) { + vector_v y_v(2); + y_v << boundary, finite; + var mu = 1.0; + var lambda = 2.0; + var out = f(y_v, mu, lambda); + out.grad(); + double adj_boundary = y_v(0).adj(); + double adj_finite = y_v(1).adj(); + double adj_mu = mu.adj(); + double adj_lambda = lambda.adj(); + stan::math::recover_memory(); - var y_s = 0.5; - var mu_s = 1.0; - var lambda_s = 2.0; - var out_s = stan::math::inv_gaussian_lcdf(y_s, mu_s, lambda_s); - out_s.grad(); - EXPECT_FLOAT_EQ(out_s.val(), vec_val); - EXPECT_FLOAT_EQ(0.0, adj_boundary); - EXPECT_FLOAT_EQ(y_s.adj(), adj_finite); - EXPECT_FLOAT_EQ(mu_s.adj(), adj_mu); - EXPECT_FLOAT_EQ(lambda_s.adj(), adj_lambda); - stan::math::recover_memory(); + var y_s = finite; + var mu_s = 1.0; + var lambda_s = 2.0; + var out_s = f(y_s, mu_s, lambda_s); + out_s.grad(); + EXPECT_FLOAT_EQ(0.0, adj_boundary); + EXPECT_FLOAT_EQ(y_s.adj(), adj_finite); + EXPECT_FLOAT_EQ(mu_s.adj(), adj_mu); + EXPECT_FLOAT_EQ(lambda_s.adj(), adj_lambda); + stan::math::recover_memory(); + }; + + check( + [](const auto& y, const auto& mu, const auto& lambda) { + return stan::math::inv_gaussian_lccdf(y, mu, lambda); + }, + 0.0, 5.0); + check( + [](const auto& y, const auto& mu, const auto& lambda) { + return stan::math::inv_gaussian_lcdf(y, mu, lambda); + }, + inf, 0.5); } From 45f645fd1349bf4b4ee9e48ac1fb451a39b3cd9c Mon Sep 17 00:00:00 2001 From: Gidon Frischkorn Date: Tue, 1 Sep 2026 20:15:23 +0200 Subject: [PATCH 11/11] Loosen platform-sensitive tolerances at extreme lambda/mu The tests at y == mu probe the rounding of y * (1 / mu) - 1, which z1 and z2 amplify by sqrt(lambda / y) -- up to ~3e9 for lambda / mu = 1e19. That rounding is zero on macOS/arm64 but ~1 ulp with clang 7 on Linux x86, so the tolerances now carry the sqrt(lambda / mu) * eps scale that bounds the error. --- .../unit/math/prim/prob/inv_gaussian_test.cpp | 23 +++++++++++++------ 1 file changed, 16 insertions(+), 7 deletions(-) diff --git a/test/unit/math/prim/prob/inv_gaussian_test.cpp b/test/unit/math/prim/prob/inv_gaussian_test.cpp index cf06170e293..a90d3a8c73d 100644 --- a/test/unit/math/prim/prob/inv_gaussian_test.cpp +++ b/test/unit/math/prim/prob/inv_gaussian_test.cpp @@ -166,8 +166,11 @@ TEST(ProbDistributionsInvGaussian, probabilityNeverExceedsOne) { } } } - // against mpmath at 60 digits - EXPECT_NEAR(-0.693147179298379049, inv_gaussian_lcdf(0.1, 0.1, 1e16), 1e-12); + // against mpmath at 60 digits. At y == mu the rounding of y * (1 / mu) - 1 + // is amplified by sqrt(lambda / y), here to about sqrt(1e17) * 1e-16 ~ 3e-8, + // and the rounding differs across platforms, so the tolerance carries that + // scale. + EXPECT_NEAR(-0.693147179298379049, inv_gaussian_lcdf(0.1, 0.1, 1e16), 1e-6); EXPECT_FLOAT_EQ(-4.04999999999999999e17, inv_gaussian_lcdf(1e-4, 1e-3, 1e14)); EXPECT_FLOAT_EQ(-4.99000499999999997e17, inv_gaussian_lcdf(1e-4, 1e-1, 1e14)); } @@ -176,6 +179,10 @@ TEST(ProbDistributionsInvGaussian, medianIsScaleInvariant) { using stan::math::inv_gaussian_lccdf; using stan::math::inv_gaussian_lcdf; for (double lambda_over_mu : {1e2, 1e6, 1e11, 1e13}) { + // y == mu probes the rounding of y * (1 / mu) - 1, which z1 and z2 + // amplify by sqrt(lambda / mu); the rounding differs across platforms, + // so the tolerance carries that scale. + double tol = 1e-13 + std::sqrt(lambda_over_mu) * 1e-15; double ref_lcdf = 0; double ref_lccdf = 0; bool first = true; @@ -188,14 +195,16 @@ TEST(ProbDistributionsInvGaussian, medianIsScaleInvariant) { ref_lccdf = b; first = false; } else { - EXPECT_NEAR(ref_lcdf, a, 1e-13); - EXPECT_NEAR(ref_lccdf, b, 1e-13); + EXPECT_NEAR(ref_lcdf, a, tol); + EXPECT_NEAR(ref_lccdf, b, tol); } } } - // F(mu) for a very large shape approaches 1/2 from above - EXPECT_NEAR(-0.693147181822, inv_gaussian_lccdf(1e-3, 1e-3, 1e14), 1e-11); - EXPECT_NEAR(-0.693147180686, inv_gaussian_lccdf(1e-3, 1e-3, 1e16), 1e-11); + // F(mu) for a very large shape approaches 1/2 from above. The same + // sqrt(lambda / mu) amplification bounds the error here, at about 3e-7 + // for lambda / mu = 1e17 and 3e-6 for 1e19. + EXPECT_NEAR(-0.693147181822, inv_gaussian_lccdf(1e-3, 1e-3, 1e14), 1e-6); + EXPECT_NEAR(-0.693147180686, inv_gaussian_lccdf(1e-3, 1e-3, 1e16), 1e-5); } TEST(ProbDistributionsInvGaussian, cdfCcdfSumToOne) {