Improve numerical stability of fwd gradients for std_normal_log_qf

Author: andrjohns

stan/math/fwd/prob/std_normal_log_qf.hpp

@@ -16,15 +16,9 @@ namespace math {
 template <typename T>
 inline fvar<T> std_normal_log_qf(const fvar<T>& p) {
   const T xv = std_normal_log_qf(p.val_);
-  int p_sign = 1;
-  auto p_d = p.d_;
-  if (p.d_ < 0) {
-    p_sign = -1;
-    p_d *= -1;
-  }
   return fvar<T>(
       xv,
-      p_sign * exp(p.val_ + log(p_d) - NEG_LOG_SQRT_TWO_PI + 0.5 * square(xv)));
+      p.d_ * exp(p.val_ - std_normal_lpdf(xv)));
 }
 }  // namespace math
 }  // namespace stan

test/unit/math/mix/prob/std_normal_log_qf_test.cpp

@@ -4,6 +4,7 @@
 #include <test/unit/math/test_ad.hpp>
 #include <stan/math/fwd/prob/std_normal_log_qf.hpp>
 
+
 TEST_F(AgradRev, mathMixLogFun_stdNormalLogQf) {
   auto f = [](const auto& x1) { return stan::math::std_normal_log_qf(x1); };
   stan::test::expect_ad(f, -100.25);
@@ -46,3 +47,14 @@ TEST_F(AgradRev, mathMixMatFunLog_stdNormalLogQfVarmat) {
   }
   expect_ad_vector_matvar(f, A);
 }
+
+TEST_F(AgradRev, GradientStabilityStdNormalLogQf) {
+  auto f = [](const auto& y) {
+    return stan::math::sum(stan::math::std_normal_log_qf(y));
+  };
+
+  Eigen::VectorXd y1(2);
+  y1 << -10, -2;
+
+  stan::test::expect_ad(f, y1);
+}