diff --git a/stan/math/prim/prob.hpp b/stan/math/prim/prob.hpp index 4fae64c1ace..085c9c9810b 100644 --- a/stan/math/prim/prob.hpp +++ b/stan/math/prim/prob.hpp @@ -139,6 +139,11 @@ #include #include #include +#include +#include +#include +#include +#include #include #include #include @@ -277,6 +282,7 @@ #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..025c17d71f1 --- /dev/null +++ b/stan/math/prim/prob/inv_gaussian_cdf.hpp @@ -0,0 +1,139 @@ +#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 + +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. + * + *

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 + * @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); + + // 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; + + 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( + select(is_inf, T_partials_return(0), + log_sum_exp(internal::std_normal_lcdf_value(z1 * INV_SQRT_TWO), + 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& is_underflow = to_ref(lcdf_elt == NEGATIVE_INFTY); + const auto& w_dens + = to_ref(cdf * exp(-0.5 * square(z1) - HALF_LOG_TWO_PI - lcdf_elt)); + const auto& w_upper = to_ref(cdf * exp(log_upper - lcdf_elt)); + if constexpr (is_autodiff_v) { + 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) + = 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) = 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); +} + +} // 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..b0108844552 --- /dev/null +++ b/stan/math/prim/prob/inv_gaussian_lccdf.hpp @@ -0,0 +1,180 @@ +#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 + +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. + * + *

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 + * @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); + + // 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; + + 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(select( + is_inf, T_partials_return(NEGATIVE_INFTY), + internal::log_diff_exp_guarded( + internal::std_normal_lcdf_value(-z1 * INV_SQRT_TWO), 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& is_underflow = to_ref(lccdf_elt == NEGATIVE_INFTY); + const auto& w_dens + = to_ref(exp(-0.5 * square(z1) - HALF_LOG_TWO_PI - lccdf_elt)); + const auto& w_upper = to_ref(exp(log_upper - lccdf_elt)); + if constexpr (is_autodiff_v) { + 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) + = 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) = 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); +} + +} // 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..f6455435e84 --- /dev/null +++ b/stan/math/prim/prob/inv_gaussian_lcdf.hpp @@ -0,0 +1,187 @@ +#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 + +namespace stan { +namespace math { +namespace internal { + +/** + * Return the log of the second inverse Gaussian CDF term, + * \f$\log\left(e^{2\lambda/\mu}\,\Phi(-z_2)\right)\f$. + * + * Both branches come from the shared standard normal log CDF in + * prim/prob/std_normal_lcdf_impl.hpp, evaluated at + * \f$\texttt{scaled} = -z_2/\sqrt{2}\f$. + * + * In the Cody regime the Gaussian exponent is taken separately, because + * \f$z_2^2 - z_1^2 = 4\lambda/\mu\f$ makes \f$2\lambda/\mu\f$ cancel it + * exactly, leaving \f$-z_1^2/2 + C(z_2)\f$ in which no two terms are large and + * opposed. Below that regime \f$z_2 \le 4\sqrt{2}\f$ bounds + * \f$2\lambda/\mu\f$ by 16 and the direct form is accurate. + */ +template * = nullptr> +inline return_type_t log_scaled_upper_term(const T1& z1, const T2& z2) { + const auto scaled = -z2 * INV_SQRT_TWO; + if (value_of_rec(scaled) < -4.0) { + return -0.5 * square(z1) + std_normal_lcdf_cody_correction(scaled); + } + return 0.5 * (square(z2) - square(z1)) + std_normal_lcdf_value(scaled); +} + +/** + * 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 in floating point at \f$y = \mu\f$. + * + *

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 + * @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); + + // 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; + + 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( + select(is_inf, T_partials_return(0), + log_sum_exp(internal::std_normal_lcdf_value(z1 * INV_SQRT_TWO), + 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& is_underflow = to_ref(lcdf_elt == NEGATIVE_INFTY); + const auto& w_dens + = to_ref(exp(-0.5 * square(z1) - HALF_LOG_TWO_PI - lcdf_elt)); + const auto& w_upper = to_ref(exp(log_upper - lcdf_elt)); + if constexpr (is_autodiff_v) { + 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) + = 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) = 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); +} + +} // namespace math +} // namespace stan +#endif 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..a2798ed27fa --- /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$. + * + *

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 + * @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); + + // 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) || (y_val == INFTY)))) { + 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); + 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..1c0a36f0675 --- /dev/null +++ b/stan/math/prim/prob/inv_gaussian_rng.hpp @@ -0,0 +1,94 @@ +#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 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 + * @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); + if (size_zero(mu, lambda)) { + return {}; + } + + 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 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; + } + + return output.data(); +} + +} // namespace math +} // namespace stan +#endif diff --git a/stan/math/prim/prob/std_normal_lcdf.hpp b/stan/math/prim/prob/std_normal_lcdf.hpp index 738fe948765..e3dc56a1a09 100644 --- a/stan/math/prim/prob/std_normal_lcdf.hpp +++ b/stan/math/prim/prob/std_normal_lcdf.hpp @@ -4,19 +4,12 @@ #include #include #include -#include -#include -#include -#include -#include #include #include #include -#include #include #include -#include -#include +#include namespace stan { namespace math { @@ -27,12 +20,15 @@ constexpr char std_normal_lcdf_func[] = "std_normal_lcdf"; /** \ingroup prob_dists * @brief Calculates the log of the cdf of the standard normal distribution * - * The piecewise structure here is the same as `normal_lcdf`, and the - * cutoffs, their provenance and the measurements behind them are documented - * once, on that function. See prim/prob/normal_lcdf.hpp. That covers the - * A&S 7.1.26, Cody (1969) and DLMF 7.12.1 references, the R `pnorm` and - * SciPy `log_ndtr` cross-references, why the erfc/Cody crossover sits at 4, - * the stan-dev/math#1411 origin of the interior Taylor cutoffs, and the two + * The piecewise value and derivative live in `internal::std_normal_lcdf_value` + * and `internal::std_normal_lcdf_grad` in prim/prob/std_normal_lcdf_impl.hpp, + * so that other distributions can share them. Their structure is the same as + * `normal_lcdf`, and the cutoffs, their provenance and the measurements + * behind them are documented once, on that function. See + * prim/prob/normal_lcdf.hpp. That covers the A&S 7.1.26, Cody (1969) and + * DLMF 7.12.1 references, the R `pnorm` and SciPy `log_ndtr` + * cross-references, why the erfc/Cody crossover sits at 4, the + * stan-dev/math#1411 origin of the interior Taylor cutoffs, and the two * cutoff tables. * * Two differences apply when reading it here. The scaled variable is @@ -55,10 +51,6 @@ template < require_all_not_nonscalar_prim_or_rev_kernel_expression_t* = nullptr> inline return_type_t std_normal_lcdf(const T_y& y) { using T_partials_return = partials_return_t; - using std::exp; - using std::fabs; - using std::log; - using std::pow; using T_y_ref = ref_type_t; static constexpr const char* function = func; T_y_ref y_ref = y; @@ -77,160 +69,17 @@ inline return_type_t std_normal_lcdf(const T_y& y) { for (size_t n = 0; n < N; n++) { const T_partials_return y_dbl = y_vec.val(n); const T_partials_return scaled_y = y_dbl * INV_SQRT_TWO; - const T_partials_return x2 = square(scaled_y); - // Rigorous numerical approximations are applied here to deal with values - // of |scaled_y|>>0. This is needed to deal with rare base-rate - // logistic regression problems where it is useful to use an alternative - // link function instead. - // - // use erfc() instead of erf() in order to retain precision - // since for x>0 erfc()->0 - if (scaled_y > 0.0) { - // CDF(x) = 1/2 + 1/2erf(x) = 1 - 1/2erfc(x) - lcdf += log1p(-0.5 * erfc(scaled_y)); - if (!is_not_nan(lcdf)) { - lcdf = 0; - } - } else if (scaled_y > -4.0) { - // CDF(x) = 1/2 - 1/2erf(-x) = 1/2erfc(-x); -4 is R pnorm's M_SQRT_32 - // crossover expressed in scaled_y, sqrt(32)/sqrt(2) - lcdf += log(erfc(-scaled_y)) + LOG_HALF; - } else if (10.0 * log(fabs(scaled_y)) - < log(std::numeric_limits::max())) { - // entering territory where erfc(-x)~0 - // need to use direct numerical approximation of lcdf instead - // the following based on W. J. Cody, Math. Comp. 23(107):631-638 (1969) - // CDF(x) = 1/2erfc(-x) - const T_partials_return x4 = pow(scaled_y, 4); - const T_partials_return x6 = pow(scaled_y, 6); - const T_partials_return x8 = pow(scaled_y, 8); - const T_partials_return x10 = pow(scaled_y, 10); - const T_partials_return temp_p - = 0.000658749161529837803157 + 0.0160837851487422766278 / x2 - + 0.125781726111229246204 / x4 + 0.360344899949804439429 / x6 - + 0.305326634961232344035 / x8 + 0.0163153871373020978498 / x10; - const T_partials_return temp_q - = -0.00233520497626869185443 - 0.0605183413124413191178 / x2 - - 0.527905102951428412248 / x4 - 1.87295284992346047209 / x6 - - 2.56852019228982242072 / x8 - 1.0 / x10; - lcdf += LOG_HALF + log(INV_SQRT_PI + (temp_p / temp_q) / x2) - - log(-scaled_y) - x2; - } else { - // scaled_y^10 term will overflow - lcdf = stan::math::negative_infinity(); + lcdf += internal::std_normal_lcdf_value(scaled_y); + // The upper branch is the only one that can carry a non-finite value + // through from a finite input, so it is guarded exactly where it was. + if (scaled_y > 0.0 && !is_not_nan(lcdf)) { + lcdf = 0; } if constexpr (is_autodiff_v) { - // compute partial derivatives - // based on analytic form given by: - // dln(CDF)/dx = exp(-x^2)/(sqrt(pi)*(1/2+erf(x)/2) - T_partials_return dnlcdf = 0.0; - T_partials_return t = 0.0; - T_partials_return t2 = 0.0; - T_partials_return t4 = 0.0; - - // calculate using piecewise function - // (due to instability / inaccuracy in the various approximations) - if (scaled_y > 2.9) { - // approximation derived from Abramowitz and Stegun (1964) 7.1.26 - t = 1.0 / (1.0 + 0.3275911 * scaled_y); - t2 = square(t); - t4 = pow(t, 4); - // A&S 7.1.26 puts exp(-x2) in the numerator; keep it there so it - // underflows to zero instead of overflowing inside a denominator - const T_partials_return exp_m_x2 = exp(-x2); - dnlcdf = INV_SQRT_PI * exp_m_x2 - / (1.0 - - exp_m_x2 - * (0.254829592 - 0.284496736 * t + 1.421413741 * t2 - - 1.453152027 * t2 * t + 1.061405429 * t4)); - } else if (scaled_y > 2.5) { - // in the trouble area where all of the standard numerical - // approximations are unstable - bridge the gap using Taylor - // expansions of the analytic function - // use Taylor expansion centred around x=2.7 - t = scaled_y - 2.7; - t2 = square(t); - t4 = pow(t, 4); - dnlcdf = 0.0003849882382 - 0.002079084702 * t + 0.005229340880 * t2 - - 0.008029540137 * t2 * t + 0.008232190507 * t4 - - 0.005692364250 * t4 * t + 0.002399496363 * pow(t, 6); - } else if (scaled_y > 2.1) { - // use Taylor expansion centred around x=2.3 - t = scaled_y - 2.3; - t2 = square(t); - t4 = pow(t, 4); - dnlcdf = 0.002846135439 - 0.01310032351 * t + 0.02732189391 * t2 - - 0.03326906904 * t2 * t + 0.02482478940 * t4 - - 0.009883071924 * t4 * t - 0.0002771362254 * pow(t, 6); - } else if (scaled_y > 1.5) { - // use Taylor expansion centred around x=1.85 - t = scaled_y - 1.85; - t2 = square(t); - t4 = pow(t, 4); - dnlcdf = 0.01849212058 - 0.06876280470 * t + 0.1099906382 * t2 - - 0.09274533184 * t2 * t + 0.03543327418 * t4 - + 0.005644855518 * t4 * t - 0.01111434424 * pow(t, 6); - } else if (scaled_y > 0.8) { - // use Taylor expansion centred around x=1.15 - t = scaled_y - 1.15; - t2 = square(t); - t4 = pow(t, 4); - dnlcdf = 0.1585747034 - 0.3898677543 * t + 0.3515963775 * t2 - - 0.09748053605 * t2 * t - 0.04347986191 * t4 - + 0.02182506378 * t4 * t + 0.01074751427 * pow(t, 6); - } else if (scaled_y > 0.1) { - // use Taylor expansion centred around x=0.45 - t = scaled_y - 0.45; - t2 = square(t); - t4 = pow(t, 4); - dnlcdf = 0.6245634904 - 0.9521866949 * t + 0.3986215682 * t2 - + 0.04700850676 * t2 * t - 0.03478651979 * t4 - - 0.01772675404 * t4 * t + 0.0006577254811 * pow(t, 6); - } else if (scaled_y < -29.0) { - // asymptotic Mills ratio, DLMF 7.12.1: dnlcdf grows linearly as - // -2*scaled_y, so no quadratic residual fit can track it - const T_partials_return inv_x2 = 1.0 / x2; - dnlcdf = -2.0 * scaled_y - / (1.0 + inv_x2 * (-0.5 + inv_x2 * (0.75 + inv_x2 * -1.875))); - } else if (10.0 * log(fabs(scaled_y)) - < log(std::numeric_limits::max())) { - // approximation derived from Abramowitz and Stegun (1964) 7.1.26 - // use fact that erf(x)=-erf(-x) - // Abramowitz and Stegun define this for -inf) { - partials<0>(ops_partials)[n] += dnlcdf * INV_SQRT_TWO; - } + partials<0>(ops_partials)[n] + += internal::std_normal_lcdf_grad(scaled_y) * INV_SQRT_TWO; } } diff --git a/stan/math/prim/prob/std_normal_lcdf_impl.hpp b/stan/math/prim/prob/std_normal_lcdf_impl.hpp new file mode 100644 index 00000000000..95719c089cb --- /dev/null +++ b/stan/math/prim/prob/std_normal_lcdf_impl.hpp @@ -0,0 +1,255 @@ +#ifndef STAN_MATH_PRIM_PROB_STD_NORMAL_LCDF_IMPL_HPP +#define STAN_MATH_PRIM_PROB_STD_NORMAL_LCDF_IMPL_HPP + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace stan { +namespace math { +namespace internal { + +/** + * The piecewise value and derivative of the standard normal log CDF. This is + * the body of `std_normal_lcdf`, moved here so that other distributions can + * share it; the inverse Gaussian cdf family is the first. + * + * Both are written in terms of the scaled variable `scaled`, related to the + * argument of \f$\log\Phi\f$ by \f$\texttt{scaled} = z/\sqrt{2}\f$. + * `std_normal_lcdf` passes `y * INV_SQRT_TWO`. `std_normal_lcdf_grad` returns + * the derivative with respect to `scaled`, so callers apply the chain rule for + * their own parameterization. + * + * The piecewise structure is the same as `normal_lcdf`, and the cutoffs, + * their provenance and the measurements behind them are documented once, on + * that function. See prim/prob/normal_lcdf.hpp. + * + * Callers reject NaN before calling. A NaN `scaled` fails every branch test + * in `std_normal_lcdf_value` and returns `-inf`. + */ + +/** + * Return the Cody (1969) lower-tail approximation to + * \f$\log\Phi(\sqrt{2}\,\texttt{scaled})\f$ with the Gaussian exponent + * \f$-\texttt{scaled}^2\f$ left off. + * + * The full branch value is `std_normal_lcdf_cody_correction(scaled) - + * square(scaled)`. The two are exposed separately because a caller that must + * add a large positive quantity to the log CDF -- the inverse Gaussian CDF + * family adds \f$2\lambda/\mu\f$ to \f$\log\Phi(-z_2)\f$ -- can cancel that + * quantity against the exponent analytically rather than forming both and + * subtracting two large opposed terms. + * + * Valid for `scaled < -4`; the caller is responsible for the branch test. + * + * @tparam T a floating point type + * @param scaled the standardized argument divided by sqrt(2) + * @return the Cody branch value excluding `-square(scaled)` + */ +template * = nullptr> +inline T std_normal_lcdf_cody_correction(const T& scaled) { + using std::log; + const T x2 = square(scaled); + const T x4 = pow(scaled, 4); + const T x6 = pow(scaled, 6); + const T x8 = pow(scaled, 8); + const T x10 = pow(scaled, 10); + const T temp_p = 0.000658749161529837803157 + 0.0160837851487422766278 / x2 + + 0.125781726111229246204 / x4 + 0.360344899949804439429 / x6 + + 0.305326634961232344035 / x8 + + 0.0163153871373020978498 / x10; + const T temp_q = -0.00233520497626869185443 - 0.0605183413124413191178 / x2 + - 0.527905102951428412248 / x4 - 1.87295284992346047209 / x6 + - 2.56852019228982242072 / x8 - 1.0 / x10; + return LOG_HALF + log(INV_SQRT_PI + (temp_p / temp_q) / x2) - log(-scaled); +} + +/** + * Return \f$\log\Phi(\sqrt{2}\,\texttt{scaled})\f$. + * + * Rigorous numerical approximations are applied here to deal with values of + * `|scaled| >> 0`. This is needed to deal with rare base-rate logistic + * regression problems where it is useful to use an alternative link function + * instead. `erfc()` is used instead of `erf()` in order to retain precision, + * since for `x > 0` `erfc() -> 0`. + * + * See the branch documentation above. + * + * @tparam T a floating point type + * @param scaled the standardized argument divided by sqrt(2) + * @return the log of the standard normal cdf + */ +template * = nullptr> +inline T std_normal_lcdf_value(const T& scaled) { + using std::fabs; + using std::log; + const T x2 = square(scaled); + if (scaled > 0.0) { + // CDF(x) = 1/2 + 1/2erf(x) = 1 - 1/2erfc(x) + return log1p(-0.5 * erfc(scaled)); + } else if (scaled > -4.0) { + // CDF(x) = 1/2 - 1/2erf(-x) = 1/2erfc(-x); -4 is R pnorm's M_SQRT_32 + // Since we scale by sqrt(2), we use sqrt(32)/sqrt(2) = 4 + return log(erfc(-scaled)) + LOG_HALF; + } else if (10.0 * log(fabs(scaled)) < log(std::numeric_limits::max())) { + // entering territory where erfc(-x)~0 + // need to use direct numerical approximation of the log cdf instead + // the following based on W. J. Cody, Math. Comp. 23(107):631-638 (1969) + // CDF(x) = 1/2erfc(-x) + return std_normal_lcdf_cody_correction(scaled) - x2; + } + // scaled^10 term will overflow + return NEGATIVE_INFTY; +} + +struct std_normal_lcdf_value_fun { + template + static inline auto fun(T&& scaled) { + return std_normal_lcdf_value(std::forward(scaled)); + } +}; + +/** + * A vectorized version of std_normal_lcdf_value(). + * + * @tparam T a container type + * @param scaled the standardized arguments divided by sqrt(2) + * @return elementwise log of the standard normal cdf + */ +template * = nullptr> +inline auto std_normal_lcdf_value(T&& scaled) { + return apply_scalar_unary::apply( + std::forward(scaled)); +} + +/** + * Return the derivative of \f$\log\Phi(\sqrt{2}\,\texttt{scaled})\f$ with + * respect to `scaled`. + * + * Based on the analytic form `dln(CDF)/dx = exp(-x^2)/(sqrt(pi)*(1/2 + + * erf(x)/2))`, calculated using a piecewise function due to instability and + * inaccuracy in the various approximations. See the branch documentation + * above. + * + * @tparam T a floating point type + * @param scaled the standardized argument divided by sqrt(2) + * @return d/d(scaled) of the log of the standard normal cdf + */ +template * = nullptr> +inline T std_normal_lcdf_grad(const T& scaled) { + using std::exp; + using std::fabs; + using std::log; + using std::pow; + const T x2 = square(scaled); + T t = 0.0; + T t2 = 0.0; + T t4 = 0.0; + + if (scaled > 2.9) { + // approximation derived from Abramowitz and Stegun (1964) 7.1.26 + t = 1.0 / (1.0 + 0.3275911 * scaled); + t2 = square(t); + t4 = pow(t, 4); + // A&S 7.1.26 puts exp(-x2) in the numerator; keep it there so it + // underflows to zero instead of overflowing inside a denominator + const T exp_m_x2 = exp(-x2); + return INV_SQRT_PI * exp_m_x2 + / (1.0 + - exp_m_x2 + * (0.254829592 - 0.284496736 * t + 1.421413741 * t2 + - 1.453152027 * t2 * t + 1.061405429 * t4)); + } else if (scaled > 2.5) { + // in the trouble area where all of the standard numerical + // approximations are unstable - bridge the gap using Taylor + // expansions of the analytic function + // use Taylor expansion centred around x=2.7 + t = scaled - 2.7; + t2 = square(t); + t4 = pow(t, 4); + return 0.0003849882382 - 0.002079084702 * t + 0.005229340880 * t2 + - 0.008029540137 * t2 * t + 0.008232190507 * t4 + - 0.005692364250 * t4 * t + 0.002399496363 * pow(t, 6); + } else if (scaled > 2.1) { + // use Taylor expansion centred around x=2.3 + t = scaled - 2.3; + t2 = square(t); + t4 = pow(t, 4); + return 0.002846135439 - 0.01310032351 * t + 0.02732189391 * t2 + - 0.03326906904 * t2 * t + 0.02482478940 * t4 + - 0.009883071924 * t4 * t - 0.0002771362254 * pow(t, 6); + } else if (scaled > 1.5) { + // use Taylor expansion centred around x=1.85 + t = scaled - 1.85; + t2 = square(t); + t4 = pow(t, 4); + return 0.01849212058 - 0.06876280470 * t + 0.1099906382 * t2 + - 0.09274533184 * t2 * t + 0.03543327418 * t4 + + 0.005644855518 * t4 * t - 0.01111434424 * pow(t, 6); + } else if (scaled > 0.8) { + // use Taylor expansion centred around x=1.15 + t = scaled - 1.15; + t2 = square(t); + t4 = pow(t, 4); + return 0.1585747034 - 0.3898677543 * t + 0.3515963775 * t2 + - 0.09748053605 * t2 * t - 0.04347986191 * t4 + + 0.02182506378 * t4 * t + 0.01074751427 * pow(t, 6); + } else if (scaled > 0.1) { + // use Taylor expansion centred around x=0.45 + t = scaled - 0.45; + t2 = square(t); + t4 = pow(t, 4); + return 0.6245634904 - 0.9521866949 * t + 0.3986215682 * t2 + + 0.04700850676 * t2 * t - 0.03478651979 * t4 + - 0.01772675404 * t4 * t + 0.0006577254811 * pow(t, 6); + } else if (scaled < -29.0) { + // asymptotic Mills ratio, DLMF 7.12.1: the derivative grows linearly as + // -2*scaled, so no quadratic residual fit can track it + const T inv_x2 = 1.0 / x2; + return -2.0 * scaled + / (1.0 + inv_x2 * (-0.5 + inv_x2 * (0.75 + inv_x2 * -1.875))); + } else if (10.0 * log(fabs(scaled)) < log(std::numeric_limits::max())) { + // approximation derived from Abramowitz and Stegun (1964) 7.1.26 + // use fact that erf(x)=-erf(-x) + // Abramowitz and Stegun define this for -inf +#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..8a930e55640 --- /dev/null +++ b/test/unit/math/mix/prob/inv_gaussian_test.cpp @@ -0,0 +1,122 @@ +#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 erfc/Cody switch of the shared standard normal log CDF at +// z2 = 4 sqrt(2) +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); +} + +// 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_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); +} 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..a174f71ba3e --- /dev/null +++ b/test/unit/math/prim/prob/inv_gaussian_test.cpp @@ -0,0 +1,276 @@ +#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, 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); +} + +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(-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)); +} + +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. +// reference values from mpmath at 60 digits +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. 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)); +} + +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; + 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, tol); + EXPECT_NEAR(ref_lccdf, b, tol); + } + } + } + // 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) { + 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); + } + } + } +} + +// The cdf family reads log Phi from the shared standard normal log CDF in +// prim/prob/std_normal_lcdf_impl.hpp, at scaled = z / sqrt(2). These +// references are the 60-digit values of log Phi(z) rounded to double; +// tolerances scale with each value's magnitude. +TEST(ProbDistributionsInvGaussian, sharedLogPhi) { + auto log_Phi = [](double z) { + return stan::math::internal::std_normal_lcdf_value( + z * stan::math::INV_SQRT_TWO); + }; + + // erfc branch, up to the switch to Cody at scaled = -4, i.e. z = -4 sqrt(2) + 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(-18.64102922378784691, log_Phi(-5.65), 1e-13); + // Cody branch; the switch sits at z = -4 sqrt(2) = -5.6568542 + EXPECT_NEAR(-18.69925115370277624, log_Phi(-5.66), 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); + EXPECT_NEAR(-454.3212439563431971, log_Phi(-30.0), 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)); + + // branch continuity across the Cody switch: the true slope at + // z = -4 sqrt(2) is about 5.8, so across a 2e-10 interval the honest change + // is about 1.2e-9 + double cut = -4.0 * stan::math::SQRT_TWO; + double eps = 1e-10; + double step = log_Phi(cut + eps) - log_Phi(cut - eps); + EXPECT_LT(std::fabs(step), 1e-7); +} + +TEST(ProbDistributionsInvGaussian, sizeMismatch) { + using stan::math::inv_gaussian_lpdf; + 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..40d5343d114 --- /dev/null +++ b/test/unit/math/rev/prob/inv_gaussian_test.cpp @@ -0,0 +1,287 @@ +#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)); + } +} + +// 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"); + } +} + +// 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(); + + 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"); + // 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"); +} + +// 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(); + + 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 = 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); +}