diff --git a/stan/math/prim/mat.hpp b/stan/math/prim/mat.hpp
index 8c9bcfa23ce..f3b26d75972 100644
--- a/stan/math/prim/mat.hpp
+++ b/stan/math/prim/mat.hpp
@@ -289,6 +289,7 @@
 #include <stan/math/prim/mat/prob/multinomial_log.hpp>
 #include <stan/math/prim/mat/prob/multinomial_lpmf.hpp>
 #include <stan/math/prim/mat/prob/multinomial_rng.hpp>
+#include <stan/math/prim/mat/prob/neg_binomial_2_log_glm_lpmf.hpp>
 #include <stan/math/prim/mat/prob/normal_id_glm_lpdf.hpp>
 #include <stan/math/prim/mat/prob/ordered_logistic_log.hpp>
 #include <stan/math/prim/mat/prob/ordered_logistic_lpmf.hpp>
diff --git a/stan/math/prim/mat/prob/bernoulli_logit_glm_lpmf.hpp b/stan/math/prim/mat/prob/bernoulli_logit_glm_lpmf.hpp
index 3f23c3944d3..5c1ab9baceb 100644
--- a/stan/math/prim/mat/prob/bernoulli_logit_glm_lpmf.hpp
+++ b/stan/math/prim/mat/prob/bernoulli_logit_glm_lpmf.hpp
@@ -1,21 +1,19 @@
 #ifndef STAN_MATH_PRIM_MAT_PROB_BERNOULLI_LOGIT_GLM_LPMF_HPP
 #define STAN_MATH_PRIM_MAT_PROB_BERNOULLI_LOGIT_GLM_LPMF_HPP
 
+
 #include <stan/math/prim/scal/meta/is_constant_struct.hpp>
 #include <stan/math/prim/scal/meta/partials_return_type.hpp>
 #include <stan/math/prim/scal/meta/operands_and_partials.hpp>
 #include <stan/math/prim/scal/err/check_consistent_sizes.hpp>
 #include <stan/math/prim/scal/err/check_bounded.hpp>
 #include <stan/math/prim/scal/err/check_finite.hpp>
-#include <stan/math/prim/scal/err/check_not_nan.hpp>
 #include <stan/math/prim/scal/fun/constants.hpp>
-#include <stan/math/prim/scal/fun/log1m.hpp>
 #include <stan/math/prim/scal/fun/value_of.hpp>
 #include <stan/math/prim/mat/fun/value_of.hpp>
 #include <stan/math/prim/scal/meta/include_summand.hpp>
 #include <stan/math/prim/scal/meta/scalar_seq_view.hpp>
 #include <stan/math/prim/scal/fun/size_zero.hpp>
-#include <boost/random/bernoulli_distribution.hpp>
 #include <boost/random/variate_generator.hpp>
 #include <cmath>
 
@@ -32,15 +30,16 @@ namespace stan {
      * should be an Eigen::Matrix type whose number of rows should match the 
      * length of n and whose number of columns should match the length of beta
      * @tparam T_beta type of the weight vector;
-     * this can also be a single double value;
+     * this can also be a single value;
      * @tparam T_alpha type of the intercept;
-     * this can either be a vector of doubles of a single double value;
+     * this has to be single value;
      * @param n binary vector parameter
      * @param x design matrix
      * @param beta weight vector
      * @param alpha intercept (in log odds)
      * @return log probability or log sum of probabilities
-     * @throw std::domain_error if theta is infinite.
+     * @throw std::domain_error if x, beta or alpha is infinite.
+     * @throw std::domain_error if n is not binary.
      * @throw std::invalid_argument if container sizes mismatch.
      */
     template <bool propto, typename T_n, typename T_x, typename T_beta,
@@ -63,9 +62,9 @@ namespace stan {
       T_partials_return logp(0.0);
 
       check_bounded(function, "Vector of dependent variables", n, 0, 1);
-      check_not_nan(function, "Matrix of independent variables", x);
-      check_not_nan(function, "Weight vector", beta);
-      check_not_nan(function, "Intercept", alpha);
+      check_finite(function, "Matrix of independent variables", x);
+      check_finite(function, "Weight vector", beta);
+      check_finite(function, "Intercept", alpha);
       check_consistent_sizes(function,
                              "Rows in matrix of independent variables",
                              x.col(0), "Vector of dependent variables",  n);
@@ -79,7 +78,7 @@ namespace stan {
       const size_t N = x.col(0).size();
       const size_t M = x.row(0).size();
 
-      Matrix<double, Dynamic, 1> signs(N, 1);
+      Matrix<T_partials_return, Dynamic, 1> signs(N, 1);
       {
         scalar_seq_view<T_n> n_vec(n);
         for (size_t n = 0; n < N; ++n) {
diff --git a/stan/math/prim/mat/prob/neg_binomial_2_log_glm_lpmf.hpp b/stan/math/prim/mat/prob/neg_binomial_2_log_glm_lpmf.hpp
new file mode 100644
index 00000000000..7d8ec443bbb
--- /dev/null
+++ b/stan/math/prim/mat/prob/neg_binomial_2_log_glm_lpmf.hpp
@@ -0,0 +1,196 @@
+#ifndef STAN_MATH_PRIM_MAT_PROB_NEG_BINOMIAL_2_LOG_GLM_LPMF_HPP
+#define STAN_MATH_PRIM_MAT_PROB_NEG_BINOMIAL_2_LOG_GLM_LPMF_HPP
+
+#include <stan/math/prim/scal/meta/is_constant_struct.hpp>
+#include <stan/math/prim/scal/meta/partials_return_type.hpp>
+#include <stan/math/prim/scal/meta/operands_and_partials.hpp>
+#include <stan/math/prim/scal/err/check_consistent_sizes.hpp>
+#include <stan/math/prim/scal/err/check_positive_finite.hpp>
+#include <stan/math/prim/scal/err/check_nonnegative.hpp>
+#include <stan/math/prim/scal/err/check_finite.hpp>
+#include <stan/math/prim/scal/fun/constants.hpp>
+#include <stan/math/prim/scal/fun/multiply_log.hpp>
+#include <stan/math/prim/scal/fun/digamma.hpp>
+#include <stan/math/prim/scal/fun/lgamma.hpp>
+#include <stan/math/prim/scal/fun/log_sum_exp.hpp>
+#include <stan/math/prim/scal/fun/value_of.hpp>
+#include <stan/math/prim/mat/fun/value_of.hpp>
+#include <stan/math/prim/scal/meta/include_summand.hpp>
+#include <stan/math/prim/scal/meta/scalar_seq_view.hpp>
+#include <cmath>
+
+namespace stan {
+  namespace math {
+
+    /**
+     * Returns the log PMF of the Generalized Linear Model (GLM)
+     * with Negative-Binomial-2 distribution and log link function.
+     * If containers are supplied, returns the log sum of the probabilities.
+     * @tparam T_n type of positive int vector of variates (labels);
+     * this can also be a single positive integer value;
+     * @tparam T_x type of the matrix of covariates (features); this
+     * should be an Eigen::Matrix type whose number of rows should match the 
+     * length of n and whose number of columns should match the length of beta
+     * @tparam T_beta type of the weight vector;
+     * this can also be a scalar;
+     * @tparam T_alpha type of the intercept;
+     * this should be a scalar;
+     * @tparam T_precision type of the (positive) precision vector phi;
+     * this can also be a scalar;
+     * @param n failures count vector parameter
+     * @param x design matrix
+     * @param beta weight vector
+     * @param alpha intercept (in log odds)
+     * @param phi (vector of) precision parameters
+     * @return log probability or log sum of probabilities
+     * @throw std::invalid_argument if container sizes mismatch.
+     * @throw std::domain_error if x, beta or alpha is infinite.
+     * @throw std::domain_error if phi is infinite or non-positive.
+     * @throw std::domain_error if n is negative.
+     */
+    template <bool propto, typename T_n, typename T_x, typename T_beta,
+              typename T_alpha, typename T_precision>
+    typename return_type<T_x, T_beta, T_alpha, T_precision>::type
+    neg_binomial_2_log_glm_lpmf(const T_n &n, const T_x &x, const T_beta &beta,
+                             const T_alpha &alpha, const T_precision &phi) {
+      static const char* function = "neg_binomial_2_log_glm_lpmf";
+      typedef typename stan::partials_return_type<T_n, T_x, T_beta,
+                                                  T_alpha, T_precision>::type
+        T_partials_return;
+
+      using std::exp;
+      using Eigen::Dynamic;
+      using Eigen::Matrix;
+      using Eigen::Array;
+
+      if (!(stan::length(n) && stan::length(x) && stan::length(beta)
+            && stan::length(phi)))
+        return 0.0;
+
+      T_partials_return logp(0.0);
+
+      check_nonnegative(function, "Failures variables", n);
+      check_finite(function, "Matrix of independent variables", x);
+      check_finite(function, "Weight vector", beta);
+      check_finite(function, "Intercept", alpha);
+      check_positive_finite(function, "Precision parameter", phi);
+      check_consistent_sizes(function,
+                             "Rows in matrix of independent variables",
+                             x.col(0), "Vector of dependent variables",  n);
+      check_consistent_sizes(function,
+                             "Rows in matrix of independent variables",
+                             x.col(0), "Vector of failure counts",  phi);
+      check_consistent_sizes(function,
+                             "Columns in matrix of independent variables",
+                             x.row(0), "Weight vector",  beta);
+
+      if (!include_summand<propto, T_x, T_beta, T_alpha, T_precision>::value)
+        return 0.0;
+
+      const size_t N = x.col(0).size();
+      const size_t M = x.row(0).size();
+
+      Array<T_partials_return, Dynamic, 1> n_arr(N, 1);
+      Array<T_partials_return, Dynamic, 1> phi_arr(N, 1);
+      {
+        scalar_seq_view<T_n> n_vec(n);
+        scalar_seq_view<T_precision> phi_vec(phi);
+        for (size_t n = 0; n < N; ++n) {
+          n_arr[n] = n_vec[n];
+          phi_arr[n] = value_of(phi_vec[n]);
+        }
+      }
+
+      Matrix<T_partials_return, Dynamic, 1> beta_dbl(M, 1);
+      {
+        scalar_seq_view<T_beta> beta_vec(beta);
+        for (size_t m = 0; m < M; ++m) {
+          beta_dbl[m] = value_of(beta_vec[m]);
+        }
+      }
+      Matrix<T_partials_return, Dynamic, Dynamic> x_dbl = value_of(x);
+
+      Array<T_partials_return, Dynamic, 1> theta_dbl = (x_dbl * beta_dbl
+      + Matrix<double, Dynamic, 1>::Ones(N, 1) * value_of(alpha)).array();
+      Array<T_partials_return, Dynamic, 1> log_phi = phi_arr.log();
+      Array<T_partials_return, Dynamic, 1> logsumexp_eta_logphi =
+        theta_dbl.binaryExpr(log_phi,
+                             [](const T_partials_return& xx,
+                                const T_partials_return& yy) {
+                               return log_sum_exp(xx, yy);
+                             });
+      Array<T_partials_return, Dynamic, 1> n_plus_phi = n_arr + phi_arr;
+
+      // Compute the log-density.
+      if (include_summand<propto>::value) {
+        logp -= (n_arr
+          + Array<double, Dynamic, 1>::Ones(N, 1))
+          .unaryExpr([](const T_partials_return& xx) {
+                       return lgamma(xx);
+                     }).sum();
+      }
+      if (include_summand<propto, T_precision>::value) {
+        logp += (phi_arr.binaryExpr(phi_arr, [](const T_partials_return& xx,
+                                                const T_partials_return& yy) {
+                                                  return multiply_log(xx, yy);
+                                              })
+                - phi_arr.unaryExpr([](const T_partials_return& xx) {
+                                      return lgamma(xx);
+                                    })).sum();
+      }
+      if (include_summand<propto, T_x, T_beta, T_alpha, T_precision>::value)
+        logp -= (n_plus_phi * logsumexp_eta_logphi).sum();
+      if (include_summand<propto, T_x, T_beta, T_alpha>::value)
+        logp += (n_arr * theta_dbl).sum();
+      if (include_summand<propto, T_precision>::value) {
+        logp += n_plus_phi.unaryExpr([](const T_partials_return& xx) {
+          return lgamma(xx);
+        }).sum();
+      }
+
+      // Compute the necessary derivatives.
+      operands_and_partials<T_x, T_beta, T_alpha,
+                            T_precision> ops_partials(x, beta, alpha, phi);
+
+      if (!(is_constant_struct<T_x>::value && is_constant_struct<T_beta>::value
+            && is_constant_struct<T_alpha>::value)) {
+        Matrix<T_partials_return, Dynamic, 1> theta_derivative(N, 1);
+        theta_derivative = (n_arr - (n_plus_phi / (phi_arr / (theta_dbl.exp())
+          + Array<double, Dynamic, 1>::Ones(N, 1)))).matrix();
+        if (!is_constant_struct<T_beta>::value) {
+          ops_partials.edge2_.partials_ = x_dbl.transpose() * theta_derivative;
+        }
+        if (!is_constant_struct<T_x>::value) {
+          ops_partials.edge1_.partials_ = theta_derivative
+            * beta_dbl.transpose();
+        }
+        if (!is_constant_struct<T_alpha>::value) {
+          ops_partials.edge3_.partials_[0] = theta_derivative.trace();
+        }
+      }
+      if (!is_constant_struct<T_precision>::value) {
+        ops_partials.edge4_.partials_ = (Array<double, Dynamic, 1>::Ones(N, 1)
+          - n_plus_phi / (theta_dbl.exp() + phi_arr) + log_phi
+          - logsumexp_eta_logphi
+          - phi_arr.unaryExpr([](const T_partials_return& xx) {
+                                return digamma(xx);
+                              })
+          + n_plus_phi.unaryExpr([](const T_partials_return& xx) {
+                                   return digamma(xx);
+                                 })).matrix();
+      }
+      return ops_partials.build(logp);
+    }
+
+    template <typename T_n, typename T_x, typename T_beta, typename T_alpha,
+              typename T_precision>
+    inline
+        typename return_type<T_x, T_beta, T_alpha, T_precision>::type
+        neg_binomial_2_log_glm_lpmf(const T_n &n, const T_x &x,
+                                    const T_beta &beta, const T_alpha &alpha,
+                                    const T_precision& phi) {
+      return neg_binomial_2_log_glm_lpmf<false>(n, x, beta, alpha, phi);
+    }
+  }  // namespace math
+}  // namespace stan
+#endif
diff --git a/test/unit/math/rev/mat/prob/neg_binomial_2_log_glm_lpmf_test.cpp b/test/unit/math/rev/mat/prob/neg_binomial_2_log_glm_lpmf_test.cpp
new file mode 100644
index 00000000000..e5bcf224811
--- /dev/null
+++ b/test/unit/math/rev/mat/prob/neg_binomial_2_log_glm_lpmf_test.cpp
@@ -0,0 +1,166 @@
+#include <stan/math/rev/mat.hpp>
+#include <gtest/gtest.h>
+#include <stan/math/prim/scal/fun/value_of.hpp>
+
+
+using stan::math::var;
+using Eigen::Dynamic;
+using Eigen::Matrix;
+
+//  We check that the values of the new regression match those of one built
+//  from existing primitives.
+TEST(ProbDistributionsNegBinomial2LogGLM,
+     glm_matches_neg_binomial_2_log_doubles) {
+  Matrix<int, Dynamic, 1> n(3, 1);
+  n << 1, 0, 1;
+  Matrix<double, Dynamic, Dynamic> x(3, 2);
+  x << -12, 46, -42,
+      24, 25, 27;
+  Matrix<double, Dynamic, 1> beta(2, 1);
+  beta << 0.3, 2;
+  double alpha = 0.3;
+  Matrix<double, Dynamic, 1> alphavec = alpha * Matrix<double, 3, 1>::Ones();
+  Matrix<double, Dynamic, 1> theta(3, 1);
+  theta = x * beta + alphavec;
+  Matrix<double, Dynamic, 1> phi(3, 1);
+  phi << 2, 1, 0.2;
+
+  EXPECT_FLOAT_EQ((stan::math::neg_binomial_2_log_lpmf(n, theta, phi)),
+                  (stan::math::neg_binomial_2_log_glm_lpmf(n, x, beta, alpha,
+                                                           phi)));
+
+  EXPECT_FLOAT_EQ((stan::math::neg_binomial_2_log_lpmf<true>(n, theta, phi)),
+                  (stan::math::neg_binomial_2_log_glm_lpmf<true>(n, x, beta,
+                                                                 alpha, phi)));
+  EXPECT_FLOAT_EQ((stan::math::neg_binomial_2_log_lpmf<false>(n, theta, phi)),
+                  (stan::math::neg_binomial_2_log_glm_lpmf<false>(n, x, beta,
+                                                                  alpha, phi)));
+  EXPECT_FLOAT_EQ(
+    (stan::math::neg_binomial_2_log_lpmf
+      <true, Matrix<int, Dynamic, 1>>(n, theta, phi)),
+    (stan::math::neg_binomial_2_log_glm_lpmf
+      <true, Matrix<int, Dynamic, 1>>(n, x, beta, alpha, phi)));
+  EXPECT_FLOAT_EQ(
+    (stan::math::neg_binomial_2_log_lpmf
+      <false, Matrix<int, Dynamic, 1>>(n, theta, phi)),
+    (stan::math::neg_binomial_2_log_glm_lpmf
+      <false, Matrix<int, Dynamic, 1>>(n, x, beta, alpha, phi)));
+  EXPECT_FLOAT_EQ(
+    (stan::math::neg_binomial_2_log_lpmf
+      <Matrix<int, Dynamic, 1>>(n, theta, phi)),
+    (stan::math::neg_binomial_2_log_glm_lpmf
+      <Matrix<int, Dynamic, 1>>(n, x, beta, alpha, phi)));
+}
+
+//  We check that the gradients of the new regression match those of one built
+//  from existing primitives.
+TEST(ProbDistributionsNegBinomial2LogGLM, glm_matches_neg_binomial_2_log_vars) {
+  Matrix<int, Dynamic, 1> n(3, 1);
+  n << 1, 0, 1;
+  Matrix<var, Dynamic, Dynamic> x(3, 2);
+  x << -12, 46, -42,
+      24, 25, 27;
+  Matrix<var, Dynamic, 1> beta(2, 1);
+  beta << 0.3, 2;
+  var alpha = 0.3;
+  Matrix<var, Dynamic, 1> alphavec = alpha * Matrix<double, 3, 1>::Ones();
+  Matrix<var, Dynamic, 1> theta(3, 1);
+  theta = x * beta + alphavec;
+  Matrix<var, Dynamic, 1> phi(3, 1);
+  phi << 2, 1, 0.2;
+
+  var lp = stan::math::neg_binomial_2_log_lpmf(n, theta, phi);
+  lp.grad();
+
+  stan::math::recover_memory();
+
+  Matrix<int, Dynamic, 1> n2(3, 1);
+  n2 << 1, 0, 1;
+  Matrix<var, Dynamic, Dynamic> x2(3, 2);
+  x2 << -12, 46, -42,
+      24, 25, 27;
+  Matrix<var, Dynamic, 1> beta2(2, 1);
+  beta2 << 0.3, 2;
+  var alpha2 = 0.3;
+  Matrix<var, Dynamic, 1> phi2(3, 1);
+  phi2 << 2, 1, 0.2;
+
+  var lp2 = stan::math::neg_binomial_2_log_glm_lpmf(n2, x2, beta2, alpha2,
+                                                    phi2);
+  lp2.grad();
+
+  EXPECT_FLOAT_EQ(lp.val(),
+                  lp2.val());
+  for (size_t i = 0; i < 2; i++) {
+    EXPECT_FLOAT_EQ(beta[i].adj(), beta2[i].adj());
+  }
+  EXPECT_FLOAT_EQ(alpha.adj(), alpha2.adj());
+  for (size_t j = 0; j < 3; j++) {
+    EXPECT_FLOAT_EQ(phi[j].adj(), phi2[j].adj());
+    for (size_t i = 0; i < 2; i++) {
+      EXPECT_FLOAT_EQ(x(j, i).adj(), x2(j, i).adj());
+    }
+  }
+}
+
+//  Here, we compare the speed of the new regression to that of one built from
+//  existing primitives.
+
+/*
+#include <chrono>
+typedef std::chrono::high_resolution_clock::time_point TimeVar;
+#define duration(a) \
+  std::chrono::duration_cast<std::chrono::microseconds>(a).count()
+#define timeNow() std::chrono::high_resolution_clock::now()
+
+TEST(ProbDistributionsNegBinomial2LogGLM, glm_matches_neg_binomial_2_log_speed) {
+  for (size_t exponent = 7; exponent < 10; exponent++) {
+    const int R = 10000;
+    const int C = std::pow(3, exponent);
+  
+  Matrix<int,Dynamic,1> n(R, 1);
+  for (size_t i = 0; i < R; i++) {
+    n[i] = rand()%2;
+  }
+  
+  int T1 = 0;
+  int T2 = 0;
+  
+  for (size_t testnumber = 0; testnumber < 30; testnumber++){
+    Matrix<double, Dynamic, Dynamic> xreal = Matrix<double, Dynamic, Dynamic>::Random(R, C);
+    Matrix<double, Dynamic, 1> betareal = Matrix<double, Dynamic, Dynamic>::Random(C, 1);
+    Matrix<double, 1, 1> alphareal = Matrix<double, 1, 1>::Random(1, 1);
+    Matrix<double, Dynamic, 1> alpharealvec = Matrix<double, R, 1>::Ones() * alphareal;
+    Matrix<double, Dynamic, 1> phireal = Matrix<double, Dynamic, Dynamic>::Random(R, 1) + Matrix<double, Dynamic, 1>::Ones(R, 1);  // We want phireal to be positive.
+    
+    Matrix<var, Dynamic, 1> beta = betareal;
+    Matrix<var, Dynamic, 1> phi = phireal;
+    Matrix<var, Dynamic, 1> theta(R, 1);
+
+  
+    TimeVar t1 = timeNow();
+    theta = (xreal * beta) + alpharealvec;                
+    var lp = stan::math::neg_binomial_2_log_lpmf(n, theta, phi);
+
+    lp.grad();
+    TimeVar t2 = timeNow();
+
+    stan::math::recover_memory();
+
+    Matrix<var, Dynamic, 1> beta2 = betareal;
+    Matrix<var, Dynamic, 1> phi2 = phireal;
+    
+    TimeVar t3 = timeNow();
+    var lp2 = stan::math::neg_binomial_2_log_glm_lpmf(n, xreal, beta2, alphareal[0], phi2);
+    lp2.grad();
+    TimeVar t4 = timeNow();
+    stan::math::recover_memory();
+    T1 += duration(t2 - t1);
+    T2 += duration(t4 - t3);
+
+  }
+  
+  std::cout << "Existing Primitives:" << std::endl << T1 << std::endl  << "New Primitives:" << std::endl << T2 << std::endl;    
+}
+}
+*/