docs: tighten inv_gaussian comments

Author: GidonFrischkorn

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 <code>inv_gaussian_lcdf</code>.
  *
- * <p>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.
+ * <p>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

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.
  *
- * <p>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.
+ * <p>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

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$.
  *
- * <p>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.
+ * <p>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

stan/math/prim/prob/inv_gaussian_lpdf.hpp

@@ -98,7 +98,6 @@ inline return_type_t<T_y, T_loc, T_shape> inv_gaussian_lpdf(
   const auto& inv_mu = to_ref_if<any_ad>(inv(mu_val));
   const auto& inv_y = to_ref_if<any_ad>(inv(y_val));
   const auto& y_m_mu = to_ref_if<any_ad>(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);

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.
  *
- * <p>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$.
+ * <p>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

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;