Reuse intermediate computations in distributions part 2

stan/math/prim/prob/chi_square_cdf.hpp

@@ -52,6 +52,7 @@ return_type_t<T_y, T_dof> chi_square_cdf(const T_y& y, const T_dof& nu) {
 
   scalar_seq_view<T_y> y_vec(y);
   scalar_seq_view<T_dof> nu_vec(nu);
+  size_t size_nu = stan::math::size(nu);
   size_t N = max_size(y, nu);
 
   // Explicit return for extreme values
@@ -62,16 +63,17 @@ return_type_t<T_y, T_dof> chi_square_cdf(const T_y& y, const T_dof& nu) {
     }
   }
 
+  VectorBuilder<!is_constant_all<T_y, T_dof>::value, T_partials_return, T_dof>
+      tgamma_vec(size_nu);
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      gamma_vec(size(nu));
-  VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digamma_vec(size(nu));
-
-  if (!is_constant_all<T_dof>::value) {
-    for (size_t i = 0; i < stan::math::size(nu); i++) {
-      const T_partials_return alpha_dbl = value_of(nu_vec[i]) * 0.5;
-      gamma_vec[i] = tgamma(alpha_dbl);
-      digamma_vec[i] = digamma(alpha_dbl);
+      digamma_vec(size_nu);
+  if (!is_constant_all<T_y, T_dof>::value) {
+    for (size_t i = 0; i < size_nu; i++) {
+      const T_partials_return half_nu_dbl = 0.5 * value_of(nu_vec[i]);
+      tgamma_vec[i] = tgamma(half_nu_dbl);
+      if (!is_constant_all<T_dof>::value) {
+        digamma_vec[i] = digamma(half_nu_dbl);
+      }
     }
   }
 
@@ -82,23 +84,21 @@ return_type_t<T_y, T_dof> chi_square_cdf(const T_y& y, const T_dof& nu) {
       continue;
     }
 
-    const T_partials_return y_dbl = value_of(y_vec[n]);
-    const T_partials_return alpha_dbl = value_of(nu_vec[n]) * 0.5;
-    const T_partials_return beta_dbl = 0.5;
-
-    const T_partials_return Pn = gamma_p(alpha_dbl, beta_dbl * y_dbl);
+    const T_partials_return half_nu_dbl = 0.5 * value_of(nu_vec[n]);
+    const T_partials_return half_y_dbl = 0.5 * value_of(y_vec[n]);
+    const T_partials_return Pn = gamma_p(half_nu_dbl, half_y_dbl);
 
     cdf *= Pn;
 
     if (!is_constant_all<T_y>::value) {
-      ops_partials.edge1_.partials_[n] += beta_dbl * exp(-beta_dbl * y_dbl)
-                                          * pow(beta_dbl * y_dbl, alpha_dbl - 1)
-                                          / tgamma(alpha_dbl) / Pn;
+      ops_partials.edge1_.partials_[n] += 0.5 * exp(-half_y_dbl)
+                                          * pow(half_y_dbl, half_nu_dbl - 1)
+                                          / (tgamma_vec[n] * Pn);
     }
     if (!is_constant_all<T_dof>::value) {
       ops_partials.edge2_.partials_[n]
           -= 0.5
-             * grad_reg_inc_gamma(alpha_dbl, beta_dbl * y_dbl, gamma_vec[n],
+             * grad_reg_inc_gamma(half_nu_dbl, half_y_dbl, tgamma_vec[n],
                                   digamma_vec[n])
              / Pn;
     }
@@ -110,7 +110,7 @@ return_type_t<T_y, T_dof> chi_square_cdf(const T_y& y, const T_dof& nu) {
     }
   }
   if (!is_constant_all<T_dof>::value) {
-    for (size_t n = 0; n < stan::math::size(nu); ++n) {
+    for (size_t n = 0; n < size_nu; ++n) {
       ops_partials.edge2_.partials_[n] *= cdf;
     }
   }

stan/math/prim/prob/chi_square_lccdf.hpp

@@ -54,53 +54,52 @@ return_type_t<T_y, T_dof> chi_square_lccdf(const T_y& y, const T_dof& nu) {
 
   scalar_seq_view<T_y> y_vec(y);
   scalar_seq_view<T_dof> nu_vec(nu);
+  size_t size_nu = stan::math::size(nu);
   size_t N = max_size(y, nu);
 
   // Explicit return for extreme values
   // The gradients are technically ill-defined, but treated as zero
   for (size_t i = 0; i < stan::math::size(y); i++) {
-    if (value_of(y_vec[i]) == 0) {
+    const T_partials_return y_dbl = value_of(y_vec[i]);
+    if (y_dbl == 0) {
       return ops_partials.build(0.0);
     }
+    if (y_dbl == INFTY) {
+      return ops_partials.build(NEGATIVE_INFTY);
+    }
   }
 
+  VectorBuilder<true, T_partials_return, T_dof> half_nu(size_nu);
+  VectorBuilder<!is_constant_all<T_y, T_dof>::value, T_partials_return, T_dof>
+      tgamma_vec(size_nu);
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      gamma_vec(size(nu));
-  VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digamma_vec(size(nu));
-
-  if (!is_constant_all<T_dof>::value) {
-    for (size_t i = 0; i < stan::math::size(nu); i++) {
-      const T_partials_return alpha_dbl = value_of(nu_vec[i]) * 0.5;
-      gamma_vec[i] = tgamma(alpha_dbl);
-      digamma_vec[i] = digamma(alpha_dbl);
+      digamma_vec(size_nu);
+  for (size_t i = 0; i < size_nu; i++) {
+    const T_partials_return half_nu_dbl = 0.5 * value_of(nu_vec[i]);
+    half_nu[i] = half_nu_dbl;
+    if (!is_constant_all<T_y, T_dof>::value) {
+      tgamma_vec[i] = tgamma(half_nu[i]);
+    }
+    if (!is_constant_all<T_dof>::value) {
+      digamma_vec[i] = digamma(half_nu_dbl);
     }
   }
 
   for (size_t n = 0; n < N; n++) {
-    // Explicit results for extreme values
-    // The gradients are technically ill-defined, but treated as zero
-    if (value_of(y_vec[n]) == INFTY) {
-      return ops_partials.build(negative_infinity());
-    }
-
-    const T_partials_return y_dbl = value_of(y_vec[n]);
-    const T_partials_return alpha_dbl = value_of(nu_vec[n]) * 0.5;
-    const T_partials_return beta_dbl = 0.5;
-
-    const T_partials_return Pn = gamma_q(alpha_dbl, beta_dbl * y_dbl);
+    const T_partials_return half_y_dbl = 0.5 * value_of(y_vec[n]);
+    const T_partials_return Pn = gamma_q(half_nu[n], half_y_dbl);
 
     ccdf_log += log(Pn);
 
     if (!is_constant_all<T_y>::value) {
-      ops_partials.edge1_.partials_[n] -= beta_dbl * exp(-beta_dbl * y_dbl)
-                                          * pow(beta_dbl * y_dbl, alpha_dbl - 1)
-                                          / tgamma(alpha_dbl) / Pn;
+      ops_partials.edge1_.partials_[n] -= 0.5 * exp(-half_y_dbl)
+                                          * pow(half_y_dbl, half_nu[n] - 1)
+                                          / (tgamma_vec[n] * Pn);
     }
     if (!is_constant_all<T_dof>::value) {
       ops_partials.edge2_.partials_[n]
           += 0.5
-             * grad_reg_inc_gamma(alpha_dbl, beta_dbl * y_dbl, gamma_vec[n],
+             * grad_reg_inc_gamma(half_nu[n], half_y_dbl, tgamma_vec[n],
                                   digamma_vec[n])
              / Pn;
     }

stan/math/prim/prob/chi_square_lcdf.hpp

@@ -54,53 +54,52 @@ return_type_t<T_y, T_dof> chi_square_lcdf(const T_y& y, const T_dof& nu) {
 
   scalar_seq_view<T_y> y_vec(y);
   scalar_seq_view<T_dof> nu_vec(nu);
+  size_t size_nu = stan::math::size(nu);
   size_t N = max_size(y, nu);
 
   // Explicit return for extreme values
   // The gradients are technically ill-defined, but treated as zero
   for (size_t i = 0; i < stan::math::size(y); i++) {
-    if (value_of(y_vec[i]) == 0) {
-      return ops_partials.build(negative_infinity());
+    const T_partials_return y_dbl = value_of(y_vec[i]);
+    if (y_dbl == 0) {
+      return ops_partials.build(NEGATIVE_INFTY);
+    }
+    if (y_dbl == INFTY) {
+      return ops_partials.build(0.0);
     }
   }
 
+  VectorBuilder<true, T_partials_return, T_dof> half_nu(size_nu);
+  VectorBuilder<!is_constant_all<T_y, T_dof>::value, T_partials_return, T_dof>
+      tgamma_vec(size_nu);
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      gamma_vec(size(nu));
-  VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digamma_vec(size(nu));
-
-  if (!is_constant_all<T_dof>::value) {
-    for (size_t i = 0; i < stan::math::size(nu); i++) {
-      const T_partials_return alpha_dbl = value_of(nu_vec[i]) * 0.5;
-      gamma_vec[i] = tgamma(alpha_dbl);
-      digamma_vec[i] = digamma(alpha_dbl);
+      digamma_vec(size_nu);
+  for (size_t i = 0; i < size_nu; i++) {
+    const T_partials_return half_nu_dbl = 0.5 * value_of(nu_vec[i]);
+    half_nu[i] = half_nu_dbl;
+    if (!is_constant_all<T_y, T_dof>::value) {
+      tgamma_vec[i] = tgamma(half_nu[i]);
+    }
+    if (!is_constant_all<T_dof>::value) {
+      digamma_vec[i] = digamma(half_nu_dbl);
     }
   }
 
   for (size_t n = 0; n < N; n++) {
-    // Explicit results for extreme values
-    // The gradients are technically ill-defined, but treated as zero
-    if (value_of(y_vec[n]) == INFTY) {
-      return ops_partials.build(0.0);
-    }
-
-    const T_partials_return y_dbl = value_of(y_vec[n]);
-    const T_partials_return alpha_dbl = value_of(nu_vec[n]) * 0.5;
-    const T_partials_return beta_dbl = 0.5;
-
-    const T_partials_return Pn = gamma_p(alpha_dbl, beta_dbl * y_dbl);
+    const T_partials_return half_y_dbl = 0.5 * value_of(y_vec[n]);
+    const T_partials_return Pn = gamma_p(half_nu[n], half_y_dbl);
 
     cdf_log += log(Pn);
 
     if (!is_constant_all<T_y>::value) {
-      ops_partials.edge1_.partials_[n] += beta_dbl * exp(-beta_dbl * y_dbl)
-                                          * pow(beta_dbl * y_dbl, alpha_dbl - 1)
-                                          / tgamma(alpha_dbl) / Pn;
+      ops_partials.edge1_.partials_[n] += 0.5 * exp(-half_y_dbl)
+                                          * pow(half_y_dbl, half_nu[n] - 1)
+                                          / (tgamma_vec[n] * Pn);
     }
     if (!is_constant_all<T_dof>::value) {
       ops_partials.edge2_.partials_[n]
           -= 0.5
-             * grad_reg_inc_gamma(alpha_dbl, beta_dbl * y_dbl, gamma_vec[n],
+             * grad_reg_inc_gamma(half_nu[n], half_y_dbl, tgamma_vec[n],
                                   digamma_vec[n])
              / Pn;
     }

stan/math/prim/prob/chi_square_log.hpp

@@ -8,25 +8,7 @@ namespace stan {
 namespace math {
 
 /** \ingroup prob_dists
- * The log of a chi-squared density for y with the specified
- * degrees of freedom parameter.
- * The degrees of freedom parameter must be greater than 0.
- * y must be greater than or equal to 0.
- *
- \f{eqnarray*}{
- y &\sim& \chi^2_\nu \\
- \log (p (y \, |\, \nu)) &=& \log \left( \frac{2^{-\nu / 2}}{\Gamma (\nu / 2)}
- y^{\nu / 2 - 1} \exp^{- y / 2} \right) \\
- &=& - \frac{\nu}{2} \log(2) - \log (\Gamma (\nu / 2)) + (\frac{\nu}{2} - 1)
- \log(y) - \frac{y}{2} \\ & & \mathrm{ where } \; y \ge 0 \f}
- *
  * @deprecated use <code>chi_square_lpdf</code>
- * @param y A scalar variable.
- * @param nu Degrees of freedom.
- * @throw std::domain_error if nu is not greater than or equal to 0
- * @throw std::domain_error if y is not greater than or equal to 0.
- * @tparam T_y Type of scalar.
- * @tparam T_dof Type of degrees of freedom.
  */
 template <bool propto, typename T_y, typename T_dof>
 return_type_t<T_y, T_dof> chi_square_log(const T_y& y, const T_dof& nu) {

stan/math/prim/prob/chi_square_lpdf.hpp

@@ -5,6 +5,7 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/constants.hpp>
 #include <stan/math/prim/fun/digamma.hpp>
+#include <stan/math/prim/fun/inv.hpp>
 #include <stan/math/prim/fun/lgamma.hpp>
 #include <stan/math/prim/fun/log.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
@@ -59,64 +60,59 @@ return_type_t<T_y, T_dof> chi_square_lpdf(const T_y& y, const T_dof& nu) {
 
   scalar_seq_view<T_y> y_vec(y);
   scalar_seq_view<T_dof> nu_vec(nu);
+  size_t size_y = stan::math::size(y);
+  size_t size_nu = stan::math::size(nu);
   size_t N = max_size(y, nu);
 
-  for (size_t n = 0; n < stan::math::size(y); n++) {
+  for (size_t n = 0; n < size_y; n++) {
     if (value_of(y_vec[n]) < 0) {
       return LOG_ZERO;
     }
   }
 
-  VectorBuilder<include_summand<propto, T_y, T_dof>::value, T_partials_return,
-                T_y>
-      log_y(size(y));
-  for (size_t i = 0; i < stan::math::size(y); i++) {
-    log_y[i] = log(value_of(y_vec[i]));
-  }
-
+  VectorBuilder<true, T_partials_return, T_y> log_y(size_y);
   VectorBuilder<include_summand<propto, T_y>::value, T_partials_return, T_y>
-      inv_y(size(y));
-  for (size_t i = 0; i < stan::math::size(y); i++) {
+      inv_y(size_y);
+  for (size_t i = 0; i < size_y; i++) {
+    const T_partials_return y_dbl = value_of(y_vec[i]);
+    log_y[i] = log(y_dbl);
     if (include_summand<propto, T_y>::value) {
-      inv_y[i] = 1.0 / value_of(y_vec[i]);
+      inv_y[i] = inv(y_dbl);
     }
   }
 
   VectorBuilder<include_summand<propto, T_dof>::value, T_partials_return, T_dof>
-      lgamma_half_nu(size(nu));
+      lgamma_half_nu(size_nu);
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digamma_half_nu_over_two(size(nu));
-
-  for (size_t i = 0; i < stan::math::size(nu); i++) {
-    T_partials_return half_nu = 0.5 * value_of(nu_vec[i]);
-    if (include_summand<propto, T_dof>::value) {
+      digamma_half_nu(size_nu);
+  if (include_summand<propto, T_dof>::value) {
+    for (size_t i = 0; i < size_nu; i++) {
+      T_partials_return half_nu = 0.5 * value_of(nu_vec[i]);
       lgamma_half_nu[i] = lgamma(half_nu);
-    }
-    if (!is_constant_all<T_dof>::value) {
-      digamma_half_nu_over_two[i] = digamma(half_nu) * 0.5;
+      if (!is_constant_all<T_dof>::value) {
+        digamma_half_nu[i] = digamma(half_nu);
+      }
     }
   }
 
   for (size_t n = 0; n < N; n++) {
-    const T_partials_return y_dbl = value_of(y_vec[n]);
-    const T_partials_return half_y = 0.5 * y_dbl;
     const T_partials_return nu_dbl = value_of(nu_vec[n]);
-    const T_partials_return half_nu = 0.5 * nu_dbl;
+    const T_partials_return half_nu_m1 = 0.5 * nu_dbl - 1;
+
     if (include_summand<propto, T_dof>::value) {
       logp -= nu_dbl * HALF_LOG_TWO + lgamma_half_nu[n];
     }
-    logp += (half_nu - 1.0) * log_y[n];
-
     if (include_summand<propto, T_y>::value) {
-      logp -= half_y;
+      logp -= 0.5 * value_of(y_vec[n]);
     }
+    logp += half_nu_m1 * log_y[n];
 
     if (!is_constant_all<T_y>::value) {
-      ops_partials.edge1_.partials_[n] += (half_nu - 1.0) * inv_y[n] - 0.5;
+      ops_partials.edge1_.partials_[n] += half_nu_m1 * inv_y[n] - 0.5;
     }
     if (!is_constant_all<T_dof>::value) {
       ops_partials.edge2_.partials_[n]
-          -= HALF_LOG_TWO + digamma_half_nu_over_two[n] - log_y[n] * 0.5;
+          -= 0.5 * (LOG_TWO + digamma_half_nu[n] - log_y[n]);
     }
   }
   return ops_partials.build(logp);