diff --git a/stan/math/prim/prob/neg_binomial_2_cdf.hpp b/stan/math/prim/prob/neg_binomial_2_cdf.hpp
index a080ebe55cb..dad7a41555c 100644
--- a/stan/math/prim/prob/neg_binomial_2_cdf.hpp
+++ b/stan/math/prim/prob/neg_binomial_2_cdf.hpp
@@ -3,11 +3,13 @@
 
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/err.hpp>
-#include <stan/math/prim/fun/size_zero.hpp>
 #include <stan/math/prim/fun/digamma.hpp>
 #include <stan/math/prim/fun/inc_beta.hpp>
 #include <stan/math/prim/fun/inc_beta_dda.hpp>
 #include <stan/math/prim/fun/inc_beta_ddz.hpp>
+#include <stan/math/prim/fun/inv.hpp>
+#include <stan/math/prim/fun/size_zero.hpp>
+#include <stan/math/prim/fun/square.hpp>
 #include <stan/math/prim/fun/value_of.hpp>
 #include <limits>
 
@@ -34,6 +36,8 @@ return_type_t<T_location, T_precision> neg_binomial_2_cdf(
   scalar_seq_view<T_n> n_vec(n);
   scalar_seq_view<T_location> mu_vec(mu);
   scalar_seq_view<T_precision> phi_vec(phi);
+  size_t size_phi = size(phi);
+  size_t size_n_phi = max_size(n, phi);
   size_t max_size_seq_view = max_size(n, mu, phi);
 
   operands_and_partials<T_location, T_precision> ops_partials(mu, phi);
@@ -48,18 +52,18 @@ return_type_t<T_location, T_precision> neg_binomial_2_cdf(
 
   VectorBuilder<!is_constant_all<T_precision>::value, T_partials_return,
                 T_precision>
-      digamma_phi_vec(size(phi));
-
-  VectorBuilder<!is_constant_all<T_precision>::value, T_partials_return,
+      digamma_phi_vec(size_phi);
+  VectorBuilder<!is_constant_all<T_precision>::value, T_partials_return, T_n,
                 T_precision>
-      digamma_sum_vec(size(phi));
+      digamma_sum_vec(size_n_phi);
 
   if (!is_constant_all<T_precision>::value) {
-    for (size_t i = 0; i < size(phi); i++) {
+    for (size_t i = 0; i < size_phi; i++) {
+      digamma_phi_vec[i] = digamma(value_of(phi_vec[i]));
+    }
+    for (size_t i = 0; i < size_n_phi; i++) {
       const T_partials_return n_dbl = value_of(n_vec[i]);
       const T_partials_return phi_dbl = value_of(phi_vec[i]);
-
-      digamma_phi_vec[i] = digamma(phi_dbl);
       digamma_sum_vec[i] = digamma(n_dbl + phi_dbl + 1);
     }
   }
@@ -71,29 +75,30 @@ return_type_t<T_location, T_precision> neg_binomial_2_cdf(
       return ops_partials.build(1.0);
     }
 
-    const T_partials_return n_dbl = value_of(n_vec[i]);
+    const T_partials_return n_dbl_p1 = value_of(n_vec[i]) + 1;
     const T_partials_return mu_dbl = value_of(mu_vec[i]);
     const T_partials_return phi_dbl = value_of(phi_vec[i]);
-
-    const T_partials_return p_dbl = phi_dbl / (mu_dbl + phi_dbl);
-    const T_partials_return d_dbl
-        = 1.0 / ((mu_dbl + phi_dbl) * (mu_dbl + phi_dbl));
-
-    const T_partials_return P_i = inc_beta(phi_dbl, n_dbl + 1.0, p_dbl);
+    const T_partials_return inv_mu_plus_phi = inv(mu_dbl + phi_dbl);
+    const T_partials_return p_dbl = phi_dbl * inv_mu_plus_phi;
+    const T_partials_return d_dbl = square(inv_mu_plus_phi);
+    const T_partials_return P_i = inc_beta(phi_dbl, n_dbl_p1, p_dbl);
+    const T_partials_return inc_beta_ddz_i
+        = is_constant_all<T_location, T_precision>::value
+              ? 0
+              : inc_beta_ddz(phi_dbl, n_dbl_p1, p_dbl) * d_dbl / P_i;
 
     P *= P_i;
 
     if (!is_constant_all<T_location>::value) {
-      ops_partials.edge1_.partials_[i]
-          += -inc_beta_ddz(phi_dbl, n_dbl + 1.0, p_dbl) * phi_dbl * d_dbl / P_i;
+      ops_partials.edge1_.partials_[i] -= inc_beta_ddz_i * phi_dbl;
     }
 
     if (!is_constant_all<T_precision>::value) {
       ops_partials.edge2_.partials_[i]
-          += inc_beta_dda(phi_dbl, n_dbl + 1, p_dbl, digamma_phi_vec[i],
+          += inc_beta_dda(phi_dbl, n_dbl_p1, p_dbl, digamma_phi_vec[i],
                           digamma_sum_vec[i])
                  / P_i
-             + inc_beta_ddz(phi_dbl, n_dbl + 1.0, p_dbl) * mu_dbl * d_dbl / P_i;
+             + inc_beta_ddz_i * mu_dbl;
     }
   }
 
@@ -104,7 +109,7 @@ return_type_t<T_location, T_precision> neg_binomial_2_cdf(
   }
 
   if (!is_constant_all<T_precision>::value) {
-    for (size_t i = 0; i < size(phi); ++i) {
+    for (size_t i = 0; i < size_phi; ++i) {
       ops_partials.edge2_.partials_[i] *= P;
     }
   }
diff --git a/stan/math/prim/prob/neg_binomial_2_lccdf.hpp b/stan/math/prim/prob/neg_binomial_2_lccdf.hpp
index 36a34b06f2b..2a2c383b199 100644
--- a/stan/math/prim/prob/neg_binomial_2_lccdf.hpp
+++ b/stan/math/prim/prob/neg_binomial_2_lccdf.hpp
@@ -27,7 +27,6 @@ return_type_t<T_location, T_precision> neg_binomial_2_lccdf(
 
   scalar_seq_view<T_location> mu_vec(mu);
   scalar_seq_view<T_precision> phi_vec(phi);
-
   size_t size_beta = max_size(mu, phi);
 
   VectorBuilder<true, return_type_t<T_location, T_precision>, T_location,
@@ -37,7 +36,7 @@ return_type_t<T_location, T_precision> neg_binomial_2_lccdf(
     beta_vec[i] = phi_vec[i] / mu_vec[i];
   }
 
-  return neg_binomial_ccdf_log(n, phi, beta_vec.data());
+  return neg_binomial_lccdf(n, phi, beta_vec.data());
 }
 
 }  // namespace math
diff --git a/stan/math/prim/prob/neg_binomial_2_lcdf.hpp b/stan/math/prim/prob/neg_binomial_2_lcdf.hpp
index 966e1102f9a..f4f36ff9282 100644
--- a/stan/math/prim/prob/neg_binomial_2_lcdf.hpp
+++ b/stan/math/prim/prob/neg_binomial_2_lcdf.hpp
@@ -3,6 +3,7 @@
 
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/err.hpp>
+#include <stan/math/prim/fun/constants.hpp>
 #include <stan/math/prim/fun/size_zero.hpp>
 #include <stan/math/prim/prob/beta_cdf_log.hpp>
 #include <cmath>
@@ -29,8 +30,15 @@ return_type_t<T_location, T_precision> neg_binomial_2_lcdf(
   scalar_seq_view<T_n> n_vec(n);
   scalar_seq_view<T_location> mu_vec(mu);
   scalar_seq_view<T_precision> phi_vec(phi);
-
+  size_t size_n = size(n);
   size_t size_phi_mu = max_size(mu, phi);
+
+  for (size_t i = 0; i < size_n; i++) {
+    if (n_vec[i] < 0) {
+      return LOG_ZERO;
+    }
+  }
+
   VectorBuilder<true, return_type_t<T_location, T_precision>, T_location,
                 T_precision>
       phi_mu(size_phi_mu);
@@ -38,14 +46,9 @@ return_type_t<T_location, T_precision> neg_binomial_2_lcdf(
     phi_mu[i] = phi_vec[i] / (phi_vec[i] + mu_vec[i]);
   }
 
-  size_t size_n = size(n);
   VectorBuilder<true, return_type_t<T_n>, T_n> np1(size_n);
   for (size_t i = 0; i < size_n; i++) {
-    if (n_vec[i] < 0) {
-      return log(0.0);
-    } else {
-      np1[i] = n_vec[i] + 1.0;
-    }
+    np1[i] = n_vec[i] + 1.0;
   }
 
   return beta_cdf_log(phi_mu.data(), phi, np1.data());
diff --git a/stan/math/prim/prob/neg_binomial_2_log_lpmf.hpp b/stan/math/prim/prob/neg_binomial_2_log_lpmf.hpp
index 9ecdd15e997..b3830ec71ae 100644
--- a/stan/math/prim/prob/neg_binomial_2_log_lpmf.hpp
+++ b/stan/math/prim/prob/neg_binomial_2_log_lpmf.hpp
@@ -46,58 +46,63 @@ return_type_t<T_log_location, T_precision> neg_binomial_2_log_lpmf(
   scalar_seq_view<T_n> n_vec(n);
   scalar_seq_view<T_log_location> eta_vec(eta);
   scalar_seq_view<T_precision> phi_vec(phi);
+  size_t size_eta = size(eta);
+  size_t size_phi = size(phi);
+  size_t size_eta_phi = max_size(eta, phi);
+  size_t size_n_phi = max_size(n, phi);
   size_t max_size_seq_view = max_size(n, eta, phi);
 
   operands_and_partials<T_log_location, T_precision> ops_partials(eta, phi);
 
-  size_t len_ep = max_size(eta, phi);
-  size_t len_np = max_size(n, phi);
-
-  VectorBuilder<true, T_partials_return, T_log_location> eta__(size(eta));
-  for (size_t i = 0, max_size_seq_view = size(eta); i < max_size_seq_view;
-       ++i) {
-    eta__[i] = value_of(eta_vec[i]);
+  VectorBuilder<true, T_partials_return, T_log_location> eta_val(size_eta);
+  for (size_t i = 0; i < size_eta; ++i) {
+    eta_val[i] = value_of(eta_vec[i]);
   }
 
-  VectorBuilder<true, T_partials_return, T_precision> phi__(size(phi));
-  for (size_t i = 0, max_size_seq_view = size(phi); i < max_size_seq_view;
-       ++i) {
-    phi__[i] = value_of(phi_vec[i]);
+  VectorBuilder<!is_constant_all<T_log_location, T_precision>::value,
+                T_partials_return, T_log_location>
+      exp_eta(size_eta);
+  if (!is_constant_all<T_log_location, T_precision>::value) {
+    for (size_t i = 0; i < size_eta; ++i) {
+      exp_eta[i] = exp(eta_val[i]);
+    }
   }
 
-  VectorBuilder<true, T_partials_return, T_precision> log_phi(size(phi));
-  for (size_t i = 0, max_size_seq_view = size(phi); i < max_size_seq_view;
-       ++i) {
-    log_phi[i] = log(phi__[i]);
+  VectorBuilder<true, T_partials_return, T_precision> phi_val(size_phi);
+  VectorBuilder<true, T_partials_return, T_precision> log_phi(size_phi);
+  for (size_t i = 0; i < size_phi; ++i) {
+    phi_val[i] = value_of(phi_vec[i]);
+    log_phi[i] = log(phi_val[i]);
   }
 
   VectorBuilder<true, T_partials_return, T_log_location, T_precision>
-      logsumexp_eta_logphi(len_ep);
-  for (size_t i = 0; i < len_ep; ++i) {
-    logsumexp_eta_logphi[i] = log_sum_exp(eta__[i], log_phi[i]);
+      logsumexp_eta_logphi(size_eta_phi);
+  for (size_t i = 0; i < size_eta_phi; ++i) {
+    logsumexp_eta_logphi[i] = log_sum_exp(eta_val[i], log_phi[i]);
   }
 
-  VectorBuilder<true, T_partials_return, T_n, T_precision> n_plus_phi(len_np);
-  for (size_t i = 0; i < len_np; ++i) {
-    n_plus_phi[i] = n_vec[i] + phi__[i];
+  VectorBuilder<true, T_partials_return, T_n, T_precision> n_plus_phi(
+      size_n_phi);
+  for (size_t i = 0; i < size_n_phi; ++i) {
+    n_plus_phi[i] = n_vec[i] + phi_val[i];
   }
 
   for (size_t i = 0; i < max_size_seq_view; i++) {
-    if (phi__[i] > 1e5) {
+    if (phi_val[i] > 1e5) {
       // TODO(martinmodrak) This is wrong (doesn't pass propto information),
       // and inaccurate for n = 0, but shouldn't break most models.
       // Also the 1e5 cutoff is way too low.
       // Will be adressed better once PR #1497 is merged
-      logp += poisson_log_lpmf(n_vec[i], eta__[i]);
+      logp += poisson_log_lpmf(n_vec[i], eta_val[i]);
     } else {
       if (include_summand<propto>::value) {
         logp -= lgamma(n_vec[i] + 1.0);
       }
       if (include_summand<propto, T_precision>::value) {
-        logp += multiply_log(phi__[i], phi__[i]) - lgamma(phi__[i]);
+        logp += multiply_log(phi_val[i], phi_val[i]) - lgamma(phi_val[i]);
       }
       if (include_summand<propto, T_log_location>::value) {
-        logp += n_vec[i] * eta__[i];
+        logp += n_vec[i] * eta_val[i];
       }
       if (include_summand<propto, T_precision>::value) {
         logp += lgamma(n_plus_phi[i]);
@@ -107,12 +112,12 @@ return_type_t<T_log_location, T_precision> neg_binomial_2_log_lpmf(
 
     if (!is_constant_all<T_log_location>::value) {
       ops_partials.edge1_.partials_[i]
-          += n_vec[i] - n_plus_phi[i] / (phi__[i] / exp(eta__[i]) + 1.0);
+          += n_vec[i] - n_plus_phi[i] / (phi_val[i] / exp_eta[i] + 1);
     }
     if (!is_constant_all<T_precision>::value) {
       ops_partials.edge2_.partials_[i]
-          += 1.0 - n_plus_phi[i] / (exp(eta__[i]) + phi__[i]) + log_phi[i]
-             - logsumexp_eta_logphi[i] - digamma(phi__[i])
+          += 1.0 - n_plus_phi[i] / (exp_eta[i] + phi_val[i]) + log_phi[i]
+             - logsumexp_eta_logphi[i] - digamma(phi_val[i])
              + digamma(n_plus_phi[i]);
     }
   }
diff --git a/stan/math/prim/prob/neg_binomial_2_lpmf.hpp b/stan/math/prim/prob/neg_binomial_2_lpmf.hpp
index a2efdaf25b7..1fe639fa771 100644
--- a/stan/math/prim/prob/neg_binomial_2_lpmf.hpp
+++ b/stan/math/prim/prob/neg_binomial_2_lpmf.hpp
@@ -3,10 +3,10 @@
 
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/err.hpp>
-#include <stan/math/prim/fun/size_zero.hpp>
-#include <stan/math/prim/fun/multiply_log.hpp>
 #include <stan/math/prim/fun/digamma.hpp>
 #include <stan/math/prim/fun/lgamma.hpp>
+#include <stan/math/prim/fun/multiply_log.hpp>
+#include <stan/math/prim/fun/size_zero.hpp>
 #include <stan/math/prim/fun/value_of.hpp>
 #include <stan/math/prim/prob/poisson_lpmf.hpp>
 #include <cmath>
@@ -42,74 +42,74 @@ return_type_t<T_location, T_precision> neg_binomial_2_lpmf(
   scalar_seq_view<T_n> n_vec(n);
   scalar_seq_view<T_location> mu_vec(mu);
   scalar_seq_view<T_precision> phi_vec(phi);
+  size_t size_mu = size(mu);
+  size_t size_phi = size(phi);
+  size_t size_mu_phi = max_size(mu, phi);
+  size_t size_n_phi = max_size(n, phi);
   size_t max_size_seq_view = max_size(n, mu, phi);
 
   operands_and_partials<T_location, T_precision> ops_partials(mu, phi);
 
-  size_t len_ep = max_size(mu, phi);
-  size_t len_np = max_size(n, phi);
-
-  VectorBuilder<true, T_partials_return, T_location> mu__(size(mu));
-
-  for (size_t i = 0, max_size_seq_view = size(mu); i < max_size_seq_view; ++i) {
-    mu__[i] = value_of(mu_vec[i]);
-  }
-
-  VectorBuilder<true, T_partials_return, T_precision> phi__(size(phi));
-  for (size_t i = 0, max_size_seq_view = size(phi); i < max_size_seq_view;
-       ++i) {
-    phi__[i] = value_of(phi_vec[i]);
+  VectorBuilder<true, T_partials_return, T_location> mu_val(size_mu);
+  for (size_t i = 0; i < size_mu; ++i) {
+    mu_val[i] = value_of(mu_vec[i]);
   }
 
-  VectorBuilder<true, T_partials_return, T_precision> log_phi(size(phi));
-  for (size_t i = 0, max_size_seq_view = size(phi); i < max_size_seq_view;
-       ++i) {
-    log_phi[i] = log(phi__[i]);
+  VectorBuilder<true, T_partials_return, T_precision> phi_val(size_phi);
+  VectorBuilder<true, T_partials_return, T_precision> log_phi(size_phi);
+  for (size_t i = 0; i < size_phi; ++i) {
+    phi_val[i] = value_of(phi_vec[i]);
+    log_phi[i] = log(phi_val[i]);
   }
 
+  VectorBuilder<true, T_partials_return, T_location, T_precision> mu_plus_phi(
+      size_mu_phi);
   VectorBuilder<true, T_partials_return, T_location, T_precision>
-      log_mu_plus_phi(len_ep);
-  for (size_t i = 0; i < len_ep; ++i) {
-    log_mu_plus_phi[i] = log(mu__[i] + phi__[i]);
+      log_mu_plus_phi(size_mu_phi);
+  for (size_t i = 0; i < size_mu_phi; ++i) {
+    mu_plus_phi[i] = mu_val[i] + phi_val[i];
+    log_mu_plus_phi[i] = log(mu_plus_phi[i]);
   }
 
-  VectorBuilder<true, T_partials_return, T_n, T_precision> n_plus_phi(len_np);
-  for (size_t i = 0; i < len_np; ++i) {
-    n_plus_phi[i] = n_vec[i] + phi__[i];
+  VectorBuilder<true, T_partials_return, T_n, T_precision> n_plus_phi(
+      size_n_phi);
+  for (size_t i = 0; i < size_n_phi; ++i) {
+    n_plus_phi[i] = n_vec[i] + phi_val[i];
   }
 
   for (size_t i = 0; i < max_size_seq_view; i++) {
     // if phi is large we probably overflow, defer to Poisson:
-    if (phi__[i] > 1e5) {
+    if (phi_val[i] > 1e5) {
       // TODO(martinmodrak) This is wrong (doesn't pass propto information),
       // and inaccurate for n = 0, but shouldn't break most models.
       // Also the 1e5 cutoff is too small.
       // Will be adressed better in PR #1497
-      logp += poisson_lpmf(n_vec[i], mu__[i]);
+      logp += poisson_lpmf(n_vec[i], mu_val[i]);
     } else {
       if (include_summand<propto>::value) {
         logp -= lgamma(n_vec[i] + 1.0);
       }
       if (include_summand<propto, T_precision>::value) {
-        logp += multiply_log(phi__[i], phi__[i]) - lgamma(phi__[i]);
+        logp += multiply_log(phi_val[i], phi_val[i]) - lgamma(phi_val[i]);
       }
       if (include_summand<propto, T_location>::value) {
-        logp += multiply_log(n_vec[i], mu__[i]);
+        logp += multiply_log(n_vec[i], mu_val[i]);
       }
       if (include_summand<propto, T_precision>::value) {
         logp += lgamma(n_plus_phi[i]);
       }
-      logp -= (n_plus_phi[i]) * log_mu_plus_phi[i];
+      logp -= n_plus_phi[i] * log_mu_plus_phi[i];
     }
 
     if (!is_constant_all<T_location>::value) {
       ops_partials.edge1_.partials_[i]
-          += n_vec[i] / mu__[i] - (n_vec[i] + phi__[i]) / (mu__[i] + phi__[i]);
+          += n_vec[i] / mu_val[i] - n_plus_phi[i] / mu_plus_phi[i];
     }
     if (!is_constant_all<T_precision>::value) {
-      ops_partials.edge2_.partials_[i]
-          += 1.0 - n_plus_phi[i] / (mu__[i] + phi__[i]) + log_phi[i]
-             - log_mu_plus_phi[i] - digamma(phi__[i]) + digamma(n_plus_phi[i]);
+      ops_partials.edge2_.partials_[i] += 1.0 - n_plus_phi[i] / mu_plus_phi[i]
+                                          + log_phi[i] - log_mu_plus_phi[i]
+                                          - digamma(phi_val[i])
+                                          + digamma(n_plus_phi[i]);
     }
   }
   return ops_partials.build(logp);
diff --git a/stan/math/prim/prob/neg_binomial_cdf.hpp b/stan/math/prim/prob/neg_binomial_cdf.hpp
index 0dca1aae2bf..70d564c53ab 100644
--- a/stan/math/prim/prob/neg_binomial_cdf.hpp
+++ b/stan/math/prim/prob/neg_binomial_cdf.hpp
@@ -7,7 +7,9 @@
 #include <stan/math/prim/fun/inc_beta.hpp>
 #include <stan/math/prim/fun/inc_beta_dda.hpp>
 #include <stan/math/prim/fun/inc_beta_ddz.hpp>
+#include <stan/math/prim/fun/inv.hpp>
 #include <stan/math/prim/fun/size_zero.hpp>
+#include <stan/math/prim/fun/square.hpp>
 #include <stan/math/prim/fun/value_of.hpp>
 #include <cmath>
 #include <limits>
@@ -36,6 +38,8 @@ return_type_t<T_shape, T_inv_scale> neg_binomial_cdf(const T_n& n,
   scalar_seq_view<T_n> n_vec(n);
   scalar_seq_view<T_shape> alpha_vec(alpha);
   scalar_seq_view<T_inv_scale> beta_vec(beta);
+  size_t size_alpha = size(alpha);
+  size_t size_n_alpha = max_size(n, alpha);
   size_t max_size_seq_view = max_size(n, alpha, beta);
 
   operands_and_partials<T_shape, T_inv_scale> ops_partials(alpha, beta);
@@ -49,17 +53,18 @@ return_type_t<T_shape, T_inv_scale> neg_binomial_cdf(const T_n& n,
   }
 
   VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
-      digamma_alpha_vec(size(alpha));
-
-  VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
-      digamma_sum_vec(size(alpha));
+      digamma_alpha_vec(size_alpha);
+  VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_n,
+                T_shape>
+      digamma_sum_vec(size_n_alpha);
 
   if (!is_constant_all<T_shape>::value) {
-    for (size_t i = 0; i < size(alpha); i++) {
+    for (size_t i = 0; i < size_alpha; i++) {
+      digamma_alpha_vec[i] = digamma(value_of(alpha_vec[i]));
+    }
+    for (size_t i = 0; i < size_n_alpha; i++) {
       const T_partials_return n_dbl = value_of(n_vec[i]);
       const T_partials_return alpha_dbl = value_of(alpha_vec[i]);
-
-      digamma_alpha_vec[i] = digamma(alpha_dbl);
       digamma_sum_vec[i] = digamma(n_dbl + alpha_dbl + 1);
     }
   }
@@ -74,9 +79,9 @@ return_type_t<T_shape, T_inv_scale> neg_binomial_cdf(const T_n& n,
     const T_partials_return n_dbl = value_of(n_vec[i]);
     const T_partials_return alpha_dbl = value_of(alpha_vec[i]);
     const T_partials_return beta_dbl = value_of(beta_vec[i]);
-
-    const T_partials_return p_dbl = beta_dbl / (1.0 + beta_dbl);
-    const T_partials_return d_dbl = 1.0 / ((1.0 + beta_dbl) * (1.0 + beta_dbl));
+    const T_partials_return inv_beta_p1 = inv(beta_dbl + 1);
+    const T_partials_return p_dbl = beta_dbl * inv_beta_p1;
+    const T_partials_return d_dbl = square(inv_beta_p1);
 
     const T_partials_return P_i = inc_beta(alpha_dbl, n_dbl + 1.0, p_dbl);
 
@@ -96,7 +101,7 @@ return_type_t<T_shape, T_inv_scale> neg_binomial_cdf(const T_n& n,
   }
 
   if (!is_constant_all<T_shape>::value) {
-    for (size_t i = 0; i < size(alpha); ++i) {
+    for (size_t i = 0; i < size_alpha; ++i) {
       ops_partials.edge1_.partials_[i] *= P;
     }
   }
diff --git a/stan/math/prim/prob/neg_binomial_lccdf.hpp b/stan/math/prim/prob/neg_binomial_lccdf.hpp
index 15323fe2a78..5152f630f82 100644
--- a/stan/math/prim/prob/neg_binomial_lccdf.hpp
+++ b/stan/math/prim/prob/neg_binomial_lccdf.hpp
@@ -3,13 +3,15 @@
 
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/err.hpp>
-#include <stan/math/prim/fun/size_zero.hpp>
+#include <stan/math/prim/fun/beta.hpp>
 #include <stan/math/prim/fun/constants.hpp>
-#include <stan/math/prim/fun/value_of.hpp>
 #include <stan/math/prim/fun/digamma.hpp>
-#include <stan/math/prim/fun/beta.hpp>
 #include <stan/math/prim/fun/grad_reg_inc_beta.hpp>
 #include <stan/math/prim/fun/inc_beta.hpp>
+#include <stan/math/prim/fun/inv.hpp>
+#include <stan/math/prim/fun/size_zero.hpp>
+#include <stan/math/prim/fun/square.hpp>
+#include <stan/math/prim/fun/value_of.hpp>
 #include <cmath>
 #include <limits>
 
@@ -33,39 +35,46 @@ return_type_t<T_shape, T_inv_scale> neg_binomial_lccdf(
   check_consistent_sizes(function, "Failures variable", n, "Shape parameter",
                          alpha, "Inverse scale parameter", beta);
 
-  scalar_seq_view<T_n> n_vec(n);
-  scalar_seq_view<T_shape> alpha_vec(alpha);
-  scalar_seq_view<T_inv_scale> beta_vec(beta);
-  size_t max_size_seq_view = max_size(n, alpha, beta);
-
   using std::exp;
   using std::log;
   using std::pow;
 
   operands_and_partials<T_shape, T_inv_scale> ops_partials(alpha, beta);
 
+  scalar_seq_view<T_n> n_vec(n);
+  scalar_seq_view<T_shape> alpha_vec(alpha);
+  scalar_seq_view<T_inv_scale> beta_vec(beta);
+  size_t size_n = size(n);
+  size_t size_alpha = size(alpha);
+  size_t size_n_alpha = max_size(n, alpha);
+  size_t max_size_seq_view = max_size(n, alpha, beta);
+
   // Explicit return for extreme values
   // The gradients are technically ill-defined, but treated as zero
-  for (size_t i = 0; i < size(n); i++) {
+  for (size_t i = 0; i < size_n; i++) {
     if (value_of(n_vec[i]) < 0) {
       return ops_partials.build(0.0);
     }
   }
 
+  VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_n>
+      digammaN_vec(size_n);
   VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
-      digammaN_vec(size(alpha));
-  VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
-      digammaAlpha_vec(size(alpha));
-  VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
-      digammaSum_vec(size(alpha));
+      digammaAlpha_vec(size_alpha);
+  VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_n,
+                T_shape>
+      digammaSum_vec(size_n_alpha);
 
   if (!is_constant_all<T_shape>::value) {
-    for (size_t i = 0; i < size(alpha); i++) {
+    for (size_t i = 0; i < size_n; i++) {
+      digammaN_vec[i] = digamma(value_of(n_vec[i]) + 1);
+    }
+    for (size_t i = 0; i < size_alpha; i++) {
+      digammaAlpha_vec[i] = digamma(value_of(alpha_vec[i]));
+    }
+    for (size_t i = 0; i < size_n_alpha; i++) {
       const T_partials_return n_dbl = value_of(n_vec[i]);
       const T_partials_return alpha_dbl = value_of(alpha_vec[i]);
-
-      digammaN_vec[i] = digamma(n_dbl + 1);
-      digammaAlpha_vec[i] = digamma(alpha_dbl);
       digammaSum_vec[i] = digamma(n_dbl + alpha_dbl + 1);
     }
   }
@@ -80,8 +89,9 @@ return_type_t<T_shape, T_inv_scale> neg_binomial_lccdf(
     const T_partials_return n_dbl = value_of(n_vec[i]);
     const T_partials_return alpha_dbl = value_of(alpha_vec[i]);
     const T_partials_return beta_dbl = value_of(beta_vec[i]);
-    const T_partials_return p_dbl = beta_dbl / (1.0 + beta_dbl);
-    const T_partials_return d_dbl = 1.0 / ((1.0 + beta_dbl) * (1.0 + beta_dbl));
+    const T_partials_return inv_beta_p1 = inv(beta_dbl + 1);
+    const T_partials_return p_dbl = beta_dbl * inv_beta_p1;
+    const T_partials_return d_dbl = square(inv_beta_p1);
     const T_partials_return Pi = 1.0 - inc_beta(alpha_dbl, n_dbl + 1.0, p_dbl);
     const T_partials_return beta_func = stan::math::beta(n_dbl + 1, alpha_dbl);
 
@@ -99,7 +109,7 @@ return_type_t<T_shape, T_inv_scale> neg_binomial_lccdf(
     if (!is_constant_all<T_inv_scale>::value) {
       ops_partials.edge2_.partials_[i] -= d_dbl * pow(1 - p_dbl, n_dbl)
                                           * pow(p_dbl, alpha_dbl - 1)
-                                          / beta_func / Pi;
+                                          / (beta_func * Pi);
     }
   }
 
diff --git a/stan/math/prim/prob/neg_binomial_lcdf.hpp b/stan/math/prim/prob/neg_binomial_lcdf.hpp
index 30a10433f18..dc7839ae9a2 100644
--- a/stan/math/prim/prob/neg_binomial_lcdf.hpp
+++ b/stan/math/prim/prob/neg_binomial_lcdf.hpp
@@ -3,13 +3,15 @@
 
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/err.hpp>
-#include <stan/math/prim/fun/size_zero.hpp>
+#include <stan/math/prim/fun/beta.hpp>
 #include <stan/math/prim/fun/constants.hpp>
-#include <stan/math/prim/fun/value_of.hpp>
 #include <stan/math/prim/fun/digamma.hpp>
-#include <stan/math/prim/fun/beta.hpp>
 #include <stan/math/prim/fun/grad_reg_inc_beta.hpp>
 #include <stan/math/prim/fun/inc_beta.hpp>
+#include <stan/math/prim/fun/inv.hpp>
+#include <stan/math/prim/fun/size_zero.hpp>
+#include <stan/math/prim/fun/square.hpp>
+#include <stan/math/prim/fun/value_of.hpp>
 #include <cmath>
 #include <limits>
 
@@ -34,39 +36,46 @@ return_type_t<T_shape, T_inv_scale> neg_binomial_lcdf(const T_n& n,
   check_consistent_sizes(function, "Failures variable", n, "Shape parameter",
                          alpha, "Inverse scale parameter", beta);
 
-  scalar_seq_view<T_n> n_vec(n);
-  scalar_seq_view<T_shape> alpha_vec(alpha);
-  scalar_seq_view<T_inv_scale> beta_vec(beta);
-  size_t max_size_seq_view = max_size(n, alpha, beta);
-
   using std::exp;
   using std::log;
   using std::pow;
 
   operands_and_partials<T_shape, T_inv_scale> ops_partials(alpha, beta);
 
+  scalar_seq_view<T_n> n_vec(n);
+  scalar_seq_view<T_shape> alpha_vec(alpha);
+  scalar_seq_view<T_inv_scale> beta_vec(beta);
+  size_t size_n = size(n);
+  size_t size_alpha = size(alpha);
+  size_t size_n_alpha = max_size(n, alpha);
+  size_t max_size_seq_view = max_size(n, alpha, beta);
+
   // Explicit return for extreme values
   // The gradients are technically ill-defined, but treated as zero
-  for (size_t i = 0; i < size(n); i++) {
+  for (size_t i = 0; i < size_n; i++) {
     if (value_of(n_vec[i]) < 0) {
       return ops_partials.build(negative_infinity());
     }
   }
 
+  VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_n>
+      digammaN_vec(size_n);
   VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
-      digammaN_vec(size(alpha));
-  VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
-      digammaAlpha_vec(size(alpha));
-  VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
-      digammaSum_vec(size(alpha));
+      digammaAlpha_vec(size_alpha);
+  VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_n,
+                T_shape>
+      digammaSum_vec(size_n_alpha);
 
   if (!is_constant_all<T_shape>::value) {
-    for (size_t i = 0; i < size(alpha); i++) {
+    for (size_t i = 0; i < size_n; i++) {
+      digammaN_vec[i] = digamma(value_of(n_vec[i]) + 1);
+    }
+    for (size_t i = 0; i < size_alpha; i++) {
+      digammaAlpha_vec[i] = digamma(value_of(alpha_vec[i]));
+    }
+    for (size_t i = 0; i < size_n_alpha; i++) {
       const T_partials_return n_dbl = value_of(n_vec[i]);
       const T_partials_return alpha_dbl = value_of(alpha_vec[i]);
-
-      digammaN_vec[i] = digamma(n_dbl + 1);
-      digammaAlpha_vec[i] = digamma(alpha_dbl);
       digammaSum_vec[i] = digamma(n_dbl + alpha_dbl + 1);
     }
   }
@@ -81,8 +90,9 @@ return_type_t<T_shape, T_inv_scale> neg_binomial_lcdf(const T_n& n,
     const T_partials_return n_dbl = value_of(n_vec[i]);
     const T_partials_return alpha_dbl = value_of(alpha_vec[i]);
     const T_partials_return beta_dbl = value_of(beta_vec[i]);
-    const T_partials_return p_dbl = beta_dbl / (1.0 + beta_dbl);
-    const T_partials_return d_dbl = 1.0 / ((1.0 + beta_dbl) * (1.0 + beta_dbl));
+    const T_partials_return inv_beta_p1 = inv(beta_dbl + 1);
+    const T_partials_return p_dbl = beta_dbl * inv_beta_p1;
+    const T_partials_return d_dbl = square(inv_beta_p1);
     const T_partials_return Pi = inc_beta(alpha_dbl, n_dbl + 1.0, p_dbl);
     const T_partials_return beta_func = stan::math::beta(n_dbl + 1, alpha_dbl);
 
@@ -100,7 +110,7 @@ return_type_t<T_shape, T_inv_scale> neg_binomial_lcdf(const T_n& n,
     if (!is_constant_all<T_inv_scale>::value) {
       ops_partials.edge2_.partials_[i] += d_dbl * pow(1 - p_dbl, n_dbl)
                                           * pow(p_dbl, alpha_dbl - 1)
-                                          / beta_func / Pi;
+                                          / (beta_func * Pi);
     }
   }
 
diff --git a/stan/math/prim/prob/neg_binomial_lpmf.hpp b/stan/math/prim/prob/neg_binomial_lpmf.hpp
index 15503191c62..eed864af703 100644
--- a/stan/math/prim/prob/neg_binomial_lpmf.hpp
+++ b/stan/math/prim/prob/neg_binomial_lpmf.hpp
@@ -3,13 +3,13 @@
 
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/err.hpp>
-#include <stan/math/prim/fun/size_zero.hpp>
-#include <stan/math/prim/fun/constants.hpp>
-#include <stan/math/prim/fun/value_of.hpp>
 #include <stan/math/prim/fun/binomial_coefficient_log.hpp>
-#include <stan/math/prim/fun/multiply_log.hpp>
+#include <stan/math/prim/fun/constants.hpp>
 #include <stan/math/prim/fun/digamma.hpp>
 #include <stan/math/prim/fun/lgamma.hpp>
+#include <stan/math/prim/fun/multiply_log.hpp>
+#include <stan/math/prim/fun/size_zero.hpp>
+#include <stan/math/prim/fun/value_of.hpp>
 #include <cmath>
 
 namespace stan {
@@ -50,60 +50,46 @@ return_type_t<T_shape, T_inv_scale> neg_binomial_lpmf(const T_n& n,
   scalar_seq_view<T_n> n_vec(n);
   scalar_seq_view<T_shape> alpha_vec(alpha);
   scalar_seq_view<T_inv_scale> beta_vec(beta);
+  size_t size_alpha = size(alpha);
+  size_t size_beta = size(beta);
+  size_t size_alpha_beta = max_size(alpha, beta);
   size_t max_size_seq_view = max_size(n, alpha, beta);
 
   operands_and_partials<T_shape, T_inv_scale> ops_partials(alpha, beta);
 
-  size_t len_ab = max_size(alpha, beta);
-  VectorBuilder<true, T_partials_return, T_shape, T_inv_scale> lambda(len_ab);
-  for (size_t i = 0; i < len_ab; ++i) {
-    lambda[i] = value_of(alpha_vec[i]) / value_of(beta_vec[i]);
-  }
-
-  VectorBuilder<true, T_partials_return, T_inv_scale> log1p_beta(size(beta));
-  for (size_t i = 0; i < size(beta); ++i) {
-    log1p_beta[i] = log1p(value_of(beta_vec[i]));
-  }
-
-  VectorBuilder<true, T_partials_return, T_inv_scale> log_beta_m_log1p_beta(
-      size(beta));
-  for (size_t i = 0; i < size(beta); ++i) {
-    log_beta_m_log1p_beta[i] = log(value_of(beta_vec[i])) - log1p_beta[i];
-  }
-
-  VectorBuilder<true, T_partials_return, T_inv_scale, T_shape>
-      alpha_times_log_beta_over_1p_beta(len_ab);
-  for (size_t i = 0; i < len_ab; ++i) {
-    alpha_times_log_beta_over_1p_beta[i]
-        = value_of(alpha_vec[i])
-          * log(value_of(beta_vec[i]) / (1.0 + value_of(beta_vec[i])));
-  }
-
   VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
-      digamma_alpha(size(alpha));
+      digamma_alpha(size_alpha);
   if (!is_constant_all<T_shape>::value) {
-    for (size_t i = 0; i < size(alpha); ++i) {
+    for (size_t i = 0; i < size_alpha; ++i) {
       digamma_alpha[i] = digamma(value_of(alpha_vec[i]));
     }
   }
 
-  VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return,
-                T_inv_scale>
-      log_beta(size(beta));
-  if (!is_constant_all<T_shape>::value) {
-    for (size_t i = 0; i < size(beta); ++i) {
-      log_beta[i] = log(value_of(beta_vec[i]));
-    }
+  VectorBuilder<true, T_partials_return, T_inv_scale> log_beta(size_beta);
+  VectorBuilder<true, T_partials_return, T_inv_scale> log1p_beta(size_beta);
+  VectorBuilder<true, T_partials_return, T_inv_scale> log_beta_m_log1p_beta(
+      size_beta);
+  for (size_t i = 0; i < size_beta; ++i) {
+    const T_partials_return beta_dbl = value_of(beta_vec[i]);
+    log_beta[i] = log(beta_dbl);
+    log1p_beta[i] = log1p(beta_dbl);
+    log_beta_m_log1p_beta[i] = log_beta[i] - log1p_beta[i];
   }
 
+  VectorBuilder<true, T_partials_return, T_shape, T_inv_scale> lambda(
+      size_alpha_beta);
+  VectorBuilder<true, T_partials_return, T_shape, T_inv_scale>
+      alpha_log_beta_over_1p_beta(size_alpha_beta);
   VectorBuilder<!is_constant_all<T_inv_scale>::value, T_partials_return,
                 T_shape, T_inv_scale>
-      lambda_m_alpha_over_1p_beta(len_ab);
-  if (!is_constant_all<T_inv_scale>::value) {
-    for (size_t i = 0; i < len_ab; ++i) {
-      lambda_m_alpha_over_1p_beta[i]
-          = lambda[i]
-            - (value_of(alpha_vec[i]) / (1.0 + value_of(beta_vec[i])));
+      lambda_m_alpha_over_1p_beta(size_alpha_beta);
+  for (size_t i = 0; i < size_alpha_beta; ++i) {
+    const T_partials_return alpha_dbl = value_of(alpha_vec[i]);
+    const T_partials_return beta_dbl = value_of(beta_vec[i]);
+    lambda[i] = alpha_dbl / beta_dbl;
+    alpha_log_beta_over_1p_beta[i] = alpha_dbl * log_beta_m_log1p_beta[i];
+    if (!is_constant_all<T_inv_scale>::value) {
+      lambda_m_alpha_over_1p_beta[i] = lambda[i] - alpha_dbl / (1 + beta_dbl);
     }
   }
 
@@ -130,7 +116,7 @@ return_type_t<T_shape, T_inv_scale> neg_binomial_lpmf(const T_n& n,
               n_vec[i] + value_of(alpha_vec[i]) - 1.0, n_vec[i]);
         }
       }
-      logp += alpha_times_log_beta_over_1p_beta[i] - n_vec[i] * log1p_beta[i];
+      logp += alpha_log_beta_over_1p_beta[i] - n_vec[i] * log1p_beta[i];
 
       if (!is_constant_all<T_shape>::value) {
         ops_partials.edge1_.partials_[i]
@@ -144,6 +130,7 @@ return_type_t<T_shape, T_inv_scale> neg_binomial_lpmf(const T_n& n,
       }
     }
   }
+
   return ops_partials.build(logp);
 }
 
diff --git a/test/unit/math/mix/prob/neg_binomial_2_log_test.cpp b/test/unit/math/mix/prob/neg_binomial_2_log_test.cpp
new file mode 100644
index 00000000000..a04071c0ccb
--- /dev/null
+++ b/test/unit/math/mix/prob/neg_binomial_2_log_test.cpp
@@ -0,0 +1,18 @@
+#include <stan/math/rev.hpp>
+#include <test/unit/math/test_ad.hpp>
+#include <gtest/gtest.h>
+#include <boost/math/differentiation/finite_difference.hpp>
+
+TEST(mathMixScalFun, neg_binomial_2_log_lpmf_derivatives) {
+  auto f1 = [](const auto& eta, const auto& phi) {
+    return stan::math::neg_binomial_2_log_lpmf(0, eta, phi);
+  };
+  auto f2 = [](const auto& eta, const auto& phi) {
+    return stan::math::neg_binomial_2_log_lpmf(6, eta, phi);
+  };
+
+  stan::test::expect_ad(f1, -1.5, 4.1);
+  stan::test::expect_ad(f1, 2.0, 1.1);
+  stan::test::expect_ad(f2, -1.5, 4.1);
+  stan::test::expect_ad(f2, 2.0, 1.1);
+}
diff --git a/test/unit/math/rev/prob/neg_binomial_2_ccdf_log_test.cpp b/test/unit/math/rev/prob/neg_binomial_2_ccdf_log_test.cpp
new file mode 100644
index 00000000000..e9c48c8992a
--- /dev/null
+++ b/test/unit/math/rev/prob/neg_binomial_2_ccdf_log_test.cpp
@@ -0,0 +1,32 @@
+#include <stan/math/rev.hpp>
+#include <gtest/gtest.h>
+#include <vector>
+
+TEST(ProbDistributionsNegBinomial2, derivatives_lccdf) {
+  using stan::math::neg_binomial_2_lccdf;
+  using stan::math::var;
+
+  std::vector<double> N{1, 2, 3};
+  double alpha_dbl = 8;
+  double beta_dbl = 1.5;
+
+  var alpha(alpha_dbl);
+  var beta(beta_dbl);
+
+  var val = neg_binomial_2_lccdf(N, alpha, beta);
+  std::vector<var> x{alpha, beta};
+  std::vector<double> gradients;
+  val.grad(x, gradients);
+
+  double epsilon = 1e-6;
+
+  double grad_diff1 = (neg_binomial_2_lccdf(N, alpha_dbl + epsilon, beta_dbl)
+                       - neg_binomial_2_lccdf(N, alpha_dbl - epsilon, beta_dbl))
+                      / (2 * epsilon);
+  EXPECT_FLOAT_EQ(grad_diff1, gradients[0]);
+
+  double grad_diff2 = (neg_binomial_2_lccdf(N, alpha_dbl, beta_dbl + epsilon)
+                       - neg_binomial_2_lccdf(N, alpha_dbl, beta_dbl - epsilon))
+                      / (2 * epsilon);
+  EXPECT_FLOAT_EQ(grad_diff2, gradients[1]);
+}
diff --git a/test/unit/math/rev/prob/neg_binomial_2_cdf_log_test.cpp b/test/unit/math/rev/prob/neg_binomial_2_cdf_log_test.cpp
new file mode 100644
index 00000000000..a29476379ee
--- /dev/null
+++ b/test/unit/math/rev/prob/neg_binomial_2_cdf_log_test.cpp
@@ -0,0 +1,32 @@
+#include <stan/math/rev.hpp>
+#include <gtest/gtest.h>
+#include <vector>
+
+TEST(ProbDistributionsNegBinomial2, derivatives_lcdf) {
+  using stan::math::neg_binomial_2_lcdf;
+  using stan::math::var;
+
+  std::vector<double> N{1, 2, 3};
+  double alpha_dbl = 8;
+  double beta_dbl = 1.5;
+
+  var alpha(alpha_dbl);
+  var beta(beta_dbl);
+
+  var val = neg_binomial_2_lcdf(N, alpha, beta);
+  std::vector<var> x{alpha, beta};
+  std::vector<double> gradients;
+  val.grad(x, gradients);
+
+  double epsilon = 1e-6;
+
+  double grad_diff1 = (neg_binomial_2_lcdf(N, alpha_dbl + epsilon, beta_dbl)
+                       - neg_binomial_2_lcdf(N, alpha_dbl - epsilon, beta_dbl))
+                      / (2 * epsilon);
+  EXPECT_FLOAT_EQ(grad_diff1, gradients[0]);
+
+  double grad_diff2 = (neg_binomial_2_lcdf(N, alpha_dbl, beta_dbl + epsilon)
+                       - neg_binomial_2_lcdf(N, alpha_dbl, beta_dbl - epsilon))
+                      / (2 * epsilon);
+  EXPECT_FLOAT_EQ(grad_diff2, gradients[1]);
+}
diff --git a/test/unit/math/rev/prob/neg_binomial_2_cdf_test.cpp b/test/unit/math/rev/prob/neg_binomial_2_cdf_test.cpp
new file mode 100644
index 00000000000..0d670b84674
--- /dev/null
+++ b/test/unit/math/rev/prob/neg_binomial_2_cdf_test.cpp
@@ -0,0 +1,32 @@
+#include <stan/math/rev.hpp>
+#include <gtest/gtest.h>
+#include <vector>
+
+TEST(ProbDistributionsNegBinomial2, derivatives_cdf) {
+  using stan::math::neg_binomial_2_cdf;
+  using stan::math::var;
+
+  std::vector<double> N{1, 2, 3};
+  double alpha_dbl = 8;
+  double beta_dbl = 1.5;
+
+  var alpha(alpha_dbl);
+  var beta(beta_dbl);
+
+  var val = neg_binomial_2_cdf(N, alpha, beta);
+  std::vector<var> x{alpha, beta};
+  std::vector<double> gradients;
+  val.grad(x, gradients);
+
+  double epsilon = 1e-6;
+
+  double grad_diff1 = (neg_binomial_2_cdf(N, alpha_dbl + epsilon, beta_dbl)
+                       - neg_binomial_2_cdf(N, alpha_dbl - epsilon, beta_dbl))
+                      / (2 * epsilon);
+  EXPECT_FLOAT_EQ(grad_diff1, gradients[0]);
+
+  double grad_diff2 = (neg_binomial_2_cdf(N, alpha_dbl, beta_dbl + epsilon)
+                       - neg_binomial_2_cdf(N, alpha_dbl, beta_dbl - epsilon))
+                      / (2 * epsilon);
+  EXPECT_FLOAT_EQ(grad_diff2, gradients[1]);
+}
diff --git a/test/unit/math/rev/prob/neg_binomial_2_log_test.cpp b/test/unit/math/rev/prob/neg_binomial_2_log_test.cpp
index ef48ba12592..fede0d0850f 100644
--- a/test/unit/math/rev/prob/neg_binomial_2_log_test.cpp
+++ b/test/unit/math/rev/prob/neg_binomial_2_log_test.cpp
@@ -1,11 +1,10 @@
 #include <stan/math/rev.hpp>
-#include <test/unit/math/test_ad.hpp>
 #include <gtest/gtest.h>
 #include <boost/math/differentiation/finite_difference.hpp>
 #include <boost/math/special_functions/digamma.hpp>
-#include <vector>
 #include <algorithm>
 #include <limits>
+#include <vector>
 
 namespace neg_binomial_2_log_test_internal {
 struct TestValue {
@@ -455,7 +454,7 @@ TEST(ProbDistributionsNegBinomial2Log, derivativesPrecomputed) {
   }
 }
 
-TEST(ProbDistributionsNegBinomial2, derivativesComplexStep) {
+TEST(ProbDistributionsNegBinomial2Log, derivativesComplexStep) {
   using boost::math::differentiation::complex_step_derivative;
   using stan::math::is_nan;
   using stan::math::neg_binomial_2_log_lpmf;
@@ -549,16 +548,28 @@ TEST(ProbDistributionsNegBinomial2, derivativesComplexStep) {
   }
 }
 
-TEST(mathMixScalFun, neg_binomial_2_log_lpmf_derivatives) {
-  auto f1 = [](const auto& eta, const auto& phi) {
-    return stan::math::neg_binomial_2_log_lpmf(0, eta, phi);
-  };
-  auto f2 = [](const auto& eta, const auto& phi) {
-    return stan::math::neg_binomial_2_log_lpmf(6, eta, phi);
-  };
+TEST(ProbDistributionsNegBinomial2Log, derivatives_diff_sizes) {
+  using stan::math::neg_binomial_2_log_lpmf;
+  using stan::math::var;
+
+  int N = 100;
+  double mu_dbl = 1.5;
+  std::vector<double> phi_dbl{2, 4, 6, 8};
+
+  var mu(mu_dbl);
+  std::vector<var> phi;
+  for (double i : phi_dbl) {
+    phi.push_back(var(i));
+  }
+  var val = neg_binomial_2_log_lpmf(N, mu, phi);
+
+  std::vector<var> x{mu};
+  std::vector<double> gradients;
+  val.grad(x, gradients);
 
-  stan::test::expect_ad(f1, -1.5, 4.1);
-  stan::test::expect_ad(f1, 2.0, 1.1);
-  stan::test::expect_ad(f2, -1.5, 4.1);
-  stan::test::expect_ad(f2, 2.0, 1.1);
+  double eps = 1e-6;
+  double grad_diff = (neg_binomial_2_log_lpmf(N, mu_dbl + eps, phi_dbl)
+                      - neg_binomial_2_log_lpmf(N, mu_dbl - eps, phi_dbl))
+                     / (2 * eps);
+  EXPECT_FLOAT_EQ(grad_diff, gradients[0]);
 }
diff --git a/test/unit/math/rev/prob/neg_binomial_2_test.cpp b/test/unit/math/rev/prob/neg_binomial_2_test.cpp
index 7f5e806acda..3dbc8cedd91 100644
--- a/test/unit/math/rev/prob/neg_binomial_2_test.cpp
+++ b/test/unit/math/rev/prob/neg_binomial_2_test.cpp
@@ -2,9 +2,9 @@
 #include <gtest/gtest.h>
 #include <vector>
 
-TEST(ProbDistributionsNegBinomial, derivatives) {
+TEST(ProbDistributionsNegBinomial2, derivatives) {
   using stan::math::is_nan;
-  using stan::math::neg_binomial_2_log;
+  using stan::math::neg_binomial_2_lpmf;
   using stan::math::var;
 
   int N = 100;
@@ -14,7 +14,7 @@ TEST(ProbDistributionsNegBinomial, derivatives) {
   for (int k = 0; k < 20; ++k) {
     var mu(mu_dbl);
     var phi(phi_dbl);
-    var val = neg_binomial_2_log(N, mu, phi);
+    var val = neg_binomial_2_lpmf(N, mu, phi);
 
     std::vector<var> x;
     x.push_back(mu);
@@ -30,10 +30,11 @@ TEST(ProbDistributionsNegBinomial, derivatives) {
     std::vector<double> finite_diffs;
     double eps = 1e-10;
     double inv2e = 0.5 / eps;
-    double dmu = neg_binomial_2_log(N, mu_dbl + eps, phi_dbl)
-                 - neg_binomial_2_log(N, mu_dbl - eps, phi_dbl);
-    double dphi = neg_binomial_2_log(N, mu_dbl, phi_dbl + eps)
-                  - neg_binomial_2_log(N, mu_dbl, phi_dbl - eps);
+
+    double dmu = neg_binomial_2_lpmf(N, mu_dbl + eps, phi_dbl)
+                 - neg_binomial_2_lpmf(N, mu_dbl - eps, phi_dbl);
+    double dphi = neg_binomial_2_lpmf(N, mu_dbl, phi_dbl + eps)
+                  - neg_binomial_2_lpmf(N, mu_dbl, phi_dbl - eps);
     finite_diffs.push_back(dmu * inv2e);
     finite_diffs.push_back(dphi * inv2e);
 
@@ -44,3 +45,29 @@ TEST(ProbDistributionsNegBinomial, derivatives) {
     phi_dbl *= 10;
   }
 }
+
+TEST(ProbDistributionsNegBinomial2, derivatives_diff_sizes) {
+  using stan::math::neg_binomial_2_lpmf;
+  using stan::math::var;
+
+  int N = 100;
+  double mu_dbl = 1.5;
+  std::vector<double> phi_dbl{2, 4, 6, 8};
+
+  var mu(mu_dbl);
+  std::vector<var> phi;
+  for (double i : phi_dbl) {
+    phi.push_back(var(i));
+  }
+  var val = neg_binomial_2_lpmf(N, mu, phi);
+
+  std::vector<var> x{mu};
+  std::vector<double> gradients;
+  val.grad(x, gradients);
+
+  double eps = 1e-6;
+  double grad_diff = (neg_binomial_2_lpmf(N, mu_dbl + eps, phi_dbl)
+                      - neg_binomial_2_lpmf(N, mu_dbl - eps, phi_dbl))
+                     / (2 * eps);
+  EXPECT_FLOAT_EQ(grad_diff, gradients[0]);
+}
diff --git a/test/unit/math/rev/prob/neg_binomial_ccdf_log_test.cpp b/test/unit/math/rev/prob/neg_binomial_ccdf_log_test.cpp
new file mode 100644
index 00000000000..d062f935a7d
--- /dev/null
+++ b/test/unit/math/rev/prob/neg_binomial_ccdf_log_test.cpp
@@ -0,0 +1,32 @@
+#include <stan/math/rev.hpp>
+#include <gtest/gtest.h>
+#include <vector>
+
+TEST(ProbDistributionsNegBinomial, derivatives_lccdf) {
+  using stan::math::neg_binomial_lccdf;
+  using stan::math::var;
+
+  std::vector<double> N{1, 2, 3};
+  double alpha_dbl = 8;
+  double beta_dbl = 1.5;
+
+  var alpha(alpha_dbl);
+  var beta(beta_dbl);
+
+  var val = neg_binomial_lccdf(N, alpha, beta);
+  std::vector<var> x{alpha, beta};
+  std::vector<double> gradients;
+  val.grad(x, gradients);
+
+  double epsilon = 1e-6;
+
+  double grad_diff1 = (neg_binomial_lccdf(N, alpha_dbl + epsilon, beta_dbl)
+                       - neg_binomial_lccdf(N, alpha_dbl - epsilon, beta_dbl))
+                      / (2 * epsilon);
+  EXPECT_FLOAT_EQ(grad_diff1, gradients[0]);
+
+  double grad_diff2 = (neg_binomial_lccdf(N, alpha_dbl, beta_dbl + epsilon)
+                       - neg_binomial_lccdf(N, alpha_dbl, beta_dbl - epsilon))
+                      / (2 * epsilon);
+  EXPECT_FLOAT_EQ(grad_diff2, gradients[1]);
+}
diff --git a/test/unit/math/rev/prob/neg_binomial_cdf_log_test.cpp b/test/unit/math/rev/prob/neg_binomial_cdf_log_test.cpp
new file mode 100644
index 00000000000..6c5ddcfef1a
--- /dev/null
+++ b/test/unit/math/rev/prob/neg_binomial_cdf_log_test.cpp
@@ -0,0 +1,32 @@
+#include <stan/math/rev.hpp>
+#include <gtest/gtest.h>
+#include <vector>
+
+TEST(ProbDistributionsNegBinomial, derivatives_lcdf) {
+  using stan::math::neg_binomial_lcdf;
+  using stan::math::var;
+
+  std::vector<double> N{1, 2, 3};
+  double alpha_dbl = 8;
+  double beta_dbl = 1.5;
+
+  var alpha(alpha_dbl);
+  var beta(beta_dbl);
+
+  var val = neg_binomial_lcdf(N, alpha, beta);
+  std::vector<var> x{alpha, beta};
+  std::vector<double> gradients;
+  val.grad(x, gradients);
+
+  double epsilon = 1e-6;
+
+  double grad_diff1 = (neg_binomial_lcdf(N, alpha_dbl + epsilon, beta_dbl)
+                       - neg_binomial_lcdf(N, alpha_dbl - epsilon, beta_dbl))
+                      / (2 * epsilon);
+  EXPECT_FLOAT_EQ(grad_diff1, gradients[0]);
+
+  double grad_diff2 = (neg_binomial_lcdf(N, alpha_dbl, beta_dbl + epsilon)
+                       - neg_binomial_lcdf(N, alpha_dbl, beta_dbl - epsilon))
+                      / (2 * epsilon);
+  EXPECT_FLOAT_EQ(grad_diff2, gradients[1]);
+}
diff --git a/test/unit/math/rev/prob/neg_binomial_cdf_test.cpp b/test/unit/math/rev/prob/neg_binomial_cdf_test.cpp
new file mode 100644
index 00000000000..a70c0283f33
--- /dev/null
+++ b/test/unit/math/rev/prob/neg_binomial_cdf_test.cpp
@@ -0,0 +1,32 @@
+#include <stan/math/rev.hpp>
+#include <gtest/gtest.h>
+#include <vector>
+
+TEST(ProbDistributionsNegBinomial, derivatives_cdf) {
+  using stan::math::neg_binomial_cdf;
+  using stan::math::var;
+
+  std::vector<double> N{1, 2, 3};
+  double alpha_dbl = 8;
+  double beta_dbl = 1.5;
+
+  var alpha(alpha_dbl);
+  var beta(beta_dbl);
+
+  var val = neg_binomial_cdf(N, alpha, beta);
+  std::vector<var> x{alpha, beta};
+  std::vector<double> gradients;
+  val.grad(x, gradients);
+
+  double epsilon = 1e-6;
+
+  double grad_diff1 = (neg_binomial_cdf(N, alpha_dbl + epsilon, beta_dbl)
+                       - neg_binomial_cdf(N, alpha_dbl - epsilon, beta_dbl))
+                      / (2 * epsilon);
+  EXPECT_FLOAT_EQ(grad_diff1, gradients[0]);
+
+  double grad_diff2 = (neg_binomial_cdf(N, alpha_dbl, beta_dbl + epsilon)
+                       - neg_binomial_cdf(N, alpha_dbl, beta_dbl - epsilon))
+                      / (2 * epsilon);
+  EXPECT_FLOAT_EQ(grad_diff2, gradients[1]);
+}
diff --git a/test/unit/math/rev/prob/neg_binomial_test.cpp b/test/unit/math/rev/prob/neg_binomial_test.cpp
new file mode 100644
index 00000000000..7dafe1ea725
--- /dev/null
+++ b/test/unit/math/rev/prob/neg_binomial_test.cpp
@@ -0,0 +1,29 @@
+#include <stan/math/rev.hpp>
+#include <gtest/gtest.h>
+#include <vector>
+
+TEST(ProbDistributionsNegBinomial, derivatives_diff_sizes) {
+  using stan::math::neg_binomial_lpmf;
+  using stan::math::var;
+
+  int N = 100;
+  double mu_dbl = 1.5;
+  std::vector<double> phi_dbl{2, 4, 6, 8};
+
+  var mu(mu_dbl);
+  std::vector<var> phi;
+  for (double i : phi_dbl) {
+    phi.push_back(var(i));
+  }
+  var val = neg_binomial_lpmf(N, mu, phi);
+
+  std::vector<var> x{mu};
+  std::vector<double> gradients;
+  val.grad(x, gradients);
+
+  double eps = 1e-6;
+  double grad_diff = (neg_binomial_lpmf(N, mu_dbl + eps, phi_dbl)
+                      - neg_binomial_lpmf(N, mu_dbl - eps, phi_dbl))
+                     / (2 * eps);
+  EXPECT_FLOAT_EQ(grad_diff, gradients[0]);
+}