diff --git a/stan/math/prim/mat.hpp b/stan/math/prim/mat.hpp
index 387e8a9cd3f..aa83d74f811 100644
--- a/stan/math/prim/mat.hpp
+++ b/stan/math/prim/mat.hpp
@@ -316,6 +316,7 @@
 #include <stan/math/prim/mat/prob/normal_id_glm_log.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_glm_lpmf.hpp>
 #include <stan/math/prim/mat/prob/ordered_logistic_lpmf.hpp>
 #include <stan/math/prim/mat/prob/ordered_logistic_rng.hpp>
 #include <stan/math/prim/mat/prob/ordered_probit_log.hpp>
diff --git a/stan/math/prim/mat/prob/ordered_logistic_glm_lpmf.hpp b/stan/math/prim/mat/prob/ordered_logistic_glm_lpmf.hpp
new file mode 100644
index 00000000000..9a5263c6bdc
--- /dev/null
+++ b/stan/math/prim/mat/prob/ordered_logistic_glm_lpmf.hpp
@@ -0,0 +1,213 @@
+#ifndef STAN_MATH_PRIM_MAT_PROB_ORDERED_LOGISTIC_GLM_LPMF_HPP
+#define STAN_MATH_PRIM_MAT_PROB_ORDERED_LOGISTIC_GLM_LPMF_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/mat/fun/value_of_rec.hpp>
+#include <stan/math/prim/arr/fun/value_of_rec.hpp>
+#include <stan/math/prim/scal/fun/size_zero.hpp>
+#include <stan/math/prim/mat/err/check_ordered.hpp>
+#include <stan/math/prim/arr/err/check_ordered.hpp>
+#include <stan/math/prim/mat/fun/log1m_exp.hpp>
+#include <stan/math/prim/meta.hpp>
+#include <cmath>
+
+namespace stan {
+namespace math {
+
+/**
+ * Returns the log PMF of the ordinal regression Generalized Linear Model (GLM).
+ * This is equivalent to and faster than ordered_logistic_lpmf(y, x * beta,
+ * cuts).
+ *
+ * @tparam T_y type of integer vector of classes. It can be either
+ * `std::vector<int>` or `int`.
+ * @tparam T_x_scalar type of elements in the matrix of independent variables
+ * (features)
+ * @tparam T_x_rows compile-time number of rows of `x`. It can be either
+ * `Eigen::Dynamic` or 1.
+ * @tparam T_beta_scalar type of a scalar in the vector of weights
+ * @tparam T_cuts_scalar type of a scalar in the vector of cutpoints
+ * @param y a scalar or vector of classes. If it is a scalar it will be
+ * broadcast - used for all instances. Values should be between 1 and number of
+ * classes, including endpoints.
+ * @param x design matrix or row vector. If it is a row vector it will be
+ * broadcast - used for all instances.
+ * @param beta weight vector
+ * @param cuts cutpoints vector
+ * @return log probability
+ * @throw std::domain_error If any class is not between 1 and
+ * the number of cutpoints plus 2 or if the cutpoint vector is not sorted in
+ * ascending order or any input is not finite
+ * @throw std::invalid_argument if container sizes mismatch.
+ */
+template <bool propto, typename T_y, typename T_x_scalar, int T_x_rows,
+          typename T_beta_scalar, typename T_cuts_scalar>
+typename stan::return_type_t<T_x_scalar, T_beta_scalar, T_cuts_scalar>
+ordered_logistic_glm_lpmf(
+    const T_y& y, const Eigen::Matrix<T_x_scalar, T_x_rows, Eigen::Dynamic>& x,
+    const Eigen::Matrix<T_beta_scalar, Eigen::Dynamic, 1>& beta,
+    const Eigen::Matrix<T_cuts_scalar, Eigen::Dynamic, 1>& cuts) {
+  using Eigen::Array;
+  using Eigen::Dynamic;
+  using Eigen::Matrix;
+  using Eigen::VectorXd;
+  using std::exp;
+  using std::isfinite;
+
+  typedef typename partials_return_type<T_y, T_x_scalar, T_beta_scalar,
+                                        T_cuts_scalar>::type T_partials_return;
+  typedef typename std::conditional_t<T_x_rows == 1, double,
+                                      Array<double, Dynamic, 1>>
+      T_location;
+
+  static const char* function = "ordered_logistic_glm_lpmf";
+
+  const size_t N_instances = T_x_rows == 1 ? length(y) : x.rows();
+  const size_t N_attributes = x.cols();
+  const size_t N_classes = length(cuts) + 1;
+
+  if (is_vector<T_y>::value && T_x_rows != 1) {
+    check_consistent_size(function, "Vector of dependent variables", y,
+                          N_instances);
+  }
+  check_consistent_size(function, "Weight vector", beta, N_attributes);
+  check_bounded(function, "Vector of dependent variables", y, 1, N_classes);
+  check_ordered(function, "Cut-points", cuts);
+  check_finite(function, "Final cut-point", cuts[N_classes - 2]);
+  check_finite(function, "First cut-point", cuts[0]);
+
+  if (size_zero(y, cuts))
+    return 0;
+
+  if (!include_summand<propto, T_x_scalar, T_beta_scalar, T_cuts_scalar>::value)
+    return 0;
+
+  const auto& x_val = value_of_rec(x);
+  const auto& beta_val = value_of_rec(beta);
+  const auto& cuts_val = value_of_rec(cuts);
+
+  const auto& beta_val_vec = as_column_vector_or_scalar(beta_val);
+  const auto& cuts_val_vec = as_column_vector_or_scalar(cuts_val);
+
+  scalar_seq_view<T_y> y_seq(y);
+  Array<double, Dynamic, 1> cuts_y1(N_instances), cuts_y2(N_instances);
+  for (int i = 0; i < N_instances; i++) {
+    int c = y_seq[i];
+    if (c != N_classes) {
+      cuts_y1[i] = cuts_val_vec[c - 1];
+    } else {
+      cuts_y1[i] = INFINITY;
+    }
+    if (c != 1) {
+      cuts_y2[i] = cuts_val_vec[c - 2];
+    } else {
+      cuts_y2[i] = -INFINITY;
+    }
+  }
+
+  T_location location = x_val * beta_val_vec;
+  if (!isfinite(sum(location))) {
+    check_finite(function, "Weight vector", beta);
+    check_finite(function, "Matrix of independent variables", x);
+  }
+
+  Array<double, Dynamic, 1> cut2 = location - cuts_y2;
+  Array<double, Dynamic, 1> cut1 = location - cuts_y1;
+
+  // Not immediately evaluating next two expressions benefits performance
+  auto m_log_1p_exp_cut1
+      = (cut1 > 0.0).select(-cut1, 0) - (-cut1.abs()).exp().log1p();
+  auto m_log_1p_exp_m_cut2
+      = (cut2 <= 0.0).select(cut2, 0) - (-cut2.abs()).exp().log1p();
+
+  T_partials_return logp(0);
+  if (is_vector<T_y>::value) {
+    Eigen::Map<const Eigen::Matrix<int, Eigen::Dynamic, 1>> y_vec(&y_seq[0],
+                                                                  y_seq.size());
+    logp = y_vec.cwiseEqual(1)
+               .select(m_log_1p_exp_cut1,
+                       y_vec.cwiseEqual(N_classes).select(
+                           m_log_1p_exp_m_cut2,
+                           m_log_1p_exp_m_cut2 + log1m_exp(cut1 - cut2).array()
+                               + m_log_1p_exp_cut1))
+               .sum();
+  } else {
+    if (y_seq[0] == 1) {
+      logp = m_log_1p_exp_cut1.sum();
+    } else if (y_seq[0] == N_classes) {
+      logp = m_log_1p_exp_m_cut2.sum();
+    } else {
+      logp = (m_log_1p_exp_m_cut2 + log1m_exp(cut1 - cut2).array()
+              + m_log_1p_exp_cut1)
+                 .sum();
+    }
+  }
+
+  operands_and_partials<Matrix<T_x_scalar, T_x_rows, Dynamic>,
+                        Eigen::Matrix<T_beta_scalar, Eigen::Dynamic, 1>,
+                        Eigen::Matrix<T_cuts_scalar, Eigen::Dynamic, 1>>
+      ops_partials(x, beta, cuts);
+  if (!is_constant_all<T_x_scalar, T_beta_scalar, T_cuts_scalar>::value) {
+    Array<double, Dynamic, 1> exp_m_cut1 = exp(-cut1);
+    Array<double, Dynamic, 1> exp_m_cut2 = exp(-cut2);
+    Array<double, Dynamic, 1> exp_cuts_diff = exp(cuts_y2 - cuts_y1);
+    Array<double, Dynamic, 1> d1
+        = (cut2 > 0).select(exp_m_cut2 / (1 + exp_m_cut2), 1 / (1 + exp(cut2)))
+          - exp_cuts_diff / (exp_cuts_diff - 1);
+    Array<double, Dynamic, 1> d2
+        = 1 / (1 - exp_cuts_diff)
+          - (cut1 > 0).select(exp_m_cut1 / (1 + exp_m_cut1),
+                              1 / (1 + exp(cut1)));
+    if (!is_constant_all<T_x_scalar, T_beta_scalar>::value) {
+      Matrix<double, 1, Dynamic> location_derivative = d1 - d2;
+      if (!is_constant_all<T_x_scalar>::value) {
+        if (T_x_rows == 1) {
+          ops_partials.edge1_.partials_
+              = beta_val_vec * location_derivative.sum();
+        } else {
+          ops_partials.edge1_.partials_
+              = (beta_val_vec * location_derivative).transpose();
+        }
+      }
+      if (!is_constant_all<T_beta_scalar>::value) {
+        if (T_x_rows == 1) {
+          ops_partials.edge2_.partials_
+              = (location_derivative * x_val.replicate(N_instances, 1))
+                    .transpose();
+        } else {
+          ops_partials.edge2_.partials_
+              = (location_derivative * x_val).transpose();
+        }
+      }
+    }
+    if (!is_constant_all<T_cuts_scalar>::value) {
+      for (int i = 0; i < N_instances; i++) {
+        int c = y_seq[i];
+        if (c != N_classes) {
+          ops_partials.edge3_.partials_[c - 1] += d2[i];
+        }
+        if (c != 1) {
+          ops_partials.edge3_.partials_[c - 2] -= d1[i];
+        }
+      }
+    }
+  }
+  return ops_partials.build(logp);
+}
+
+template <typename T_y, typename T_x_scalar, int T_x_rows,
+          typename T_beta_scalar, typename T_cuts_scalar>
+typename stan::return_type_t<T_x_scalar, T_beta_scalar, T_cuts_scalar>
+ordered_logistic_glm_lpmf(
+    const T_y& y, const Eigen::Matrix<T_x_scalar, T_x_rows, Eigen::Dynamic>& x,
+    const Eigen::Matrix<T_beta_scalar, Eigen::Dynamic, 1>& beta,
+    const Eigen::Matrix<T_cuts_scalar, Eigen::Dynamic, 1>& cuts) {
+  return ordered_logistic_glm_lpmf<false>(y, x, beta, cuts);
+}
+
+}  // namespace math
+}  // namespace stan
+
+#endif
diff --git a/test/unit/math/rev/mat/prob/ordered_logistic_glm_lpmf_test.cpp b/test/unit/math/rev/mat/prob/ordered_logistic_glm_lpmf_test.cpp
new file mode 100644
index 00000000000..dfeedbd49ee
--- /dev/null
+++ b/test/unit/math/rev/mat/prob/ordered_logistic_glm_lpmf_test.cpp
@@ -0,0 +1,447 @@
+#include <stan/math.hpp>
+#include <gtest/gtest.h>
+#include <Eigen/Core>
+#include <vector>
+
+using Eigen::Array;
+using Eigen::Dynamic;
+using Eigen::Matrix;
+using Eigen::MatrixXd;
+using Eigen::RowVectorXd;
+using Eigen::VectorXd;
+using stan::math::ordered_logistic_glm_lpmf;
+using stan::math::ordered_logistic_lpmf;
+using stan::math::var;
+using std::vector;
+
+template <bool propto, typename T_x, typename T_beta, typename T_cuts>
+typename stan::return_type<T_x, T_beta, T_cuts>::type
+ordered_logistic_glm_simple_lpmf(const vector<int>& y,
+                                 const Matrix<T_x, Dynamic, Dynamic>& x,
+                                 const T_beta& beta, const T_cuts& cuts) {
+  typedef typename stan::return_type<T_x, T_beta>::type T_x_beta;
+  using stan::math::as_column_vector_or_scalar;
+
+  auto& beta_col = as_column_vector_or_scalar(beta);
+
+  Eigen::Matrix<T_x_beta, Dynamic, 1> location
+      = x.template cast<T_x_beta>() * beta_col.template cast<T_x_beta>();
+
+  return ordered_logistic_lpmf<propto>(y, location, cuts);
+}
+
+TEST(ProbDistributionsOrderedLogisticGLM,
+     glm_matches_ordered_logistic_doubles) {
+  double eps = 1e-13;
+  int N = 5;
+  int M = 2;
+  int C = 4;
+  vector<int> y{1, 4, 3, 3, 2};
+  VectorXd cuts(C - 1);
+  cuts << 0.9, 1.1, 7;
+  VectorXd beta(M);
+  beta << 1.1, 0.4;
+  MatrixXd x(N, M);
+  x << 1, 2, 3, 4, 5, 6, 7, 8, 9, 0;
+  EXPECT_FLOAT_EQ(ordered_logistic_glm_lpmf(y, x, beta, cuts),
+                  ordered_logistic_glm_simple_lpmf<false>(y, x, beta, cuts));
+}
+
+TEST(ProbDistributionsOrderedLogisticGLM,
+     glm_matches_ordered_logistic_doubles_broadcast_y) {
+  double eps = 1e-13;
+  int N = 5;
+  int M = 2;
+  int C = 4;
+  VectorXd cuts(C - 1);
+  cuts << 0.9, 1.1, 7;
+  VectorXd beta(M);
+  beta << 1.1, 0.4;
+  MatrixXd x(N, M);
+  x << 1, 2, 3, 4, 5, 6, 7, 8, 9, 0;
+  for (int y_scal = 1; y_scal <= C; y_scal++) {
+    vector<int> y(N, y_scal);
+    EXPECT_FLOAT_EQ(ordered_logistic_glm_lpmf(y_scal, x, beta, cuts),
+                    ordered_logistic_glm_lpmf(y, x, beta, cuts));
+  }
+}
+
+TEST(ProbDistributionsOrderedLogisticGLM, glm_matches_ordered_logistic_vars) {
+  double eps = 1e-13;
+  int N = 5;
+  int M = 2;
+  int C = 3;
+  vector<int> y{1, 1, 2, 4, 4};
+  Matrix<var, Dynamic, 1> cuts1(C), cuts2(C);
+  cuts1 << 0.9, 1.1, 7;
+  cuts2 << 0.9, 1.1, 7;
+  Matrix<var, Dynamic, 1> beta1(M), beta2(M);
+  beta1 << 1.1, 0.4;
+  beta2 << 1.1, 0.4;
+  Matrix<var, Dynamic, Dynamic> x1(N, M), x2(N, M);
+  x1 << 1, 2, 3, 4, 5, 6, 7, 8, 9, 0;
+  x2 << 1, 2, 3, 4, 5, 6, 7, 8, 9, 0;
+  var res1 = ordered_logistic_glm_lpmf(y, x1, beta1, cuts1);
+  var res2 = ordered_logistic_glm_simple_lpmf<false>(y, x2, beta2, cuts2);
+  (res1 + res2).grad();
+
+  EXPECT_NEAR(res1.val(), res2.val(), eps);
+  for (int i = 0; i < M; i++) {
+    for (int j = 0; j < N; j++) {
+      EXPECT_NEAR(x1(j, i).adj(), x2(j, i).adj(), eps);
+    }
+  }
+  for (int i = 0; i < M; i++) {
+    EXPECT_NEAR(beta1[i].adj(), beta2[i].adj(), eps);
+  }
+  for (int i = 0; i < C; i++) {
+    EXPECT_NEAR(cuts1[i].adj(), cuts2[i].adj(), eps);
+  }
+}
+
+TEST(ProbDistributionsOrderedLogisticGLM,
+     glm_matches_ordered_logistic_vars_broadcast_y) {
+  double eps = 1e-13;
+  int N = 5;
+  int M = 2;
+  int C = 3;
+  Matrix<var, Dynamic, 1> cuts1(C), cuts2(C);
+  cuts1 << 0.9, 1.1, 7;
+  cuts2 << 0.9, 1.1, 7;
+  Matrix<var, Dynamic, 1> beta1(M), beta2(M);
+  beta1 << 1.1, 0.4;
+  beta2 << 1.1, 0.4;
+  Matrix<var, Dynamic, Dynamic> x1(N, M), x2(N, M);
+  x1 << 1, 2, 3, 4, 5, 6, 7, 8, 9, 0;
+  x2 << 1, 2, 3, 4, 5, 6, 7, 8, 9, 0;
+
+  for (int y_scal = 1; y_scal <= C; y_scal++) {
+    vector<int> y(N, y_scal);
+    var res1 = ordered_logistic_glm_lpmf(y, x1, beta1, cuts1);
+    var res2 = ordered_logistic_glm_lpmf(y_scal, x2, beta2, cuts2);
+    (res1 + res2).grad();
+
+    Matrix<double, Dynamic, Dynamic> x_adj(N, M);
+
+    EXPECT_NEAR(res1.val(), res2.val(), eps);
+    for (int i = 0; i < M; i++) {
+      for (int j = 0; j < N; j++) {
+        x_adj(j, i) = x2(j, i).adj();
+        EXPECT_NEAR(x1(j, i).adj(), x2(j, i).adj(), eps);
+      }
+    }
+    for (int i = 0; i < M; i++) {
+      EXPECT_NEAR(beta1[i].adj(), beta2[i].adj(), eps);
+    }
+    for (int i = 0; i < C; i++) {
+      EXPECT_NEAR(cuts1[i].adj(), cuts2[i].adj(), eps);
+    }
+  }
+}
+
+TEST(ProbDistributionsOrderedLogisticGLM,
+     glm_matches_ordered_logistic_vars_broadcast_x) {
+  double eps = 1e-13;
+  int N = 5;
+  int M = 2;
+  int C = 3;
+  vector<int> y{1, 1, 2, 4, 4};
+  Matrix<var, Dynamic, 1> cuts1(C), cuts2(C);
+  cuts1 << 0.9, 1.1, 7;
+  cuts2 << 0.9, 1.1, 7;
+  Matrix<var, Dynamic, 1> beta1(M), beta2(M);
+  beta1 << 1.1, 0.4;
+  beta2 << 1.1, 0.4;
+  RowVectorXd x_double(M);
+  x_double << 1, 2;
+  Matrix<var, 1, Dynamic> x_row = x_double;
+  Matrix<var, Dynamic, Dynamic> x = x_double.replicate(N, 1);
+
+  var res1 = ordered_logistic_glm_lpmf(y, x_row, beta1, cuts1);
+  var res2 = ordered_logistic_glm_lpmf(y, x, beta2, cuts2);
+  (res1 + res2).grad();
+
+  EXPECT_NEAR(res1.val(), res2.val(), eps);
+  for (int i = 0; i < M; i++) {
+    double x_sum = 0;
+    for (int j = 0; j < N; j++) {
+      x_sum += x(j, i).adj();
+    }
+    EXPECT_NEAR(x_row[i].adj(), x_sum, eps);
+  }
+  for (int i = 0; i < M; i++) {
+    EXPECT_NEAR(beta1[i].adj(), beta2[i].adj(), eps);
+  }
+  for (int i = 0; i < C; i++) {
+    EXPECT_NEAR(cuts1[i].adj(), cuts2[i].adj(), eps);
+  }
+}
+
+TEST(ProbDistributionsOrderedLogisticGLM,
+     glm_matches_ordered_logistic_single_instance) {
+  double eps = 1e-13;
+  int N = 1;
+  int M = 2;
+  int C = 3;
+  vector<int> y{1};
+  Matrix<var, Dynamic, 1> cuts1(C), cuts2(C);
+  cuts1 << 0.9, 1.1, 7;
+  cuts2 << 0.9, 1.1, 7;
+  Matrix<var, Dynamic, 1> beta1(M), beta2(M);
+  beta1 << 1.1, 0.4;
+  beta2 << 1.1, 0.4;
+  Matrix<var, Dynamic, Dynamic> x1(N, M), x2(N, M);
+  x1 << 1, 2;
+  x2 << 1, 2;
+  var res1 = ordered_logistic_glm_lpmf(y, x1, beta1, cuts1);
+  var res2 = ordered_logistic_glm_simple_lpmf<false>(y, x2, beta2, cuts2);
+  (res1 + res2).grad();
+
+  EXPECT_NEAR(res1.val(), res2.val(), eps);
+  for (int i = 0; i < M; i++) {
+    EXPECT_NEAR(x1(0, i).adj(), x2(0, i).adj(), eps);
+  }
+  for (int i = 0; i < M; i++) {
+    EXPECT_NEAR(beta1[i].adj(), beta2[i].adj(), eps);
+  }
+  for (int i = 0; i < C; i++) {
+    EXPECT_NEAR(cuts1[i].adj(), cuts2[i].adj(), eps);
+  }
+}
+
+TEST(ProbDistributionsOrderedLogisticGLM,
+     glm_matches_ordered_logistic_no_attributes) {
+  double eps = 1e-13;
+  int N = 5;
+  int M = 0;
+  int C = 3;
+  vector<int> y{1, 1, 2, 4, 4};
+  Matrix<var, Dynamic, 1> cuts1(C), cuts2(C);
+  cuts1 << 0.9, 1.1, 7;
+  cuts2 << 0.9, 1.1, 7;
+  Matrix<var, Dynamic, 1> beta1(M), beta2(M);
+  Matrix<var, Dynamic, Dynamic> x1(N, M), x2(N, M);
+  var res1 = ordered_logistic_glm_lpmf(y, x1, beta1, cuts1);
+  var res2 = ordered_logistic_glm_simple_lpmf<false>(y, x2, beta2, cuts2);
+  (res1 + res2).grad();
+
+  EXPECT_NEAR(res1.val(), res2.val(), eps);
+  for (int i = 0; i < C; i++) {
+    EXPECT_NEAR(cuts1[i].adj(), cuts2[i].adj(), eps);
+  }
+}
+
+TEST(ProbDistributionsOrderedLogisticGLM,
+     glm_matches_ordered_logistic_single_class) {
+  double eps = 1e-13;
+  int N = 5;
+  int M = 2;
+  int C = 1;
+  vector<int> y{1, 1, 1, 1, 1};
+  Matrix<var, Dynamic, 1> cuts1(C), cuts2(C);
+  cuts1 << 0.9;
+  cuts2 << 0.9;
+  Matrix<var, Dynamic, 1> beta1(M), beta2(M);
+  beta1 << 1.1, 0.4;
+  beta2 << 1.1, 0.4;
+  Matrix<var, Dynamic, Dynamic> x1(N, M), x2(N, M);
+  x1 << 1, 2, 3, 4, 5, 6, 7, 8, 9, 0;
+  x2 << 1, 2, 3, 4, 5, 6, 7, 8, 9, 0;
+  var res1 = ordered_logistic_glm_lpmf(y, x1, beta1, cuts1);
+  var res2 = ordered_logistic_glm_simple_lpmf<false>(y, x2, beta2, cuts2);
+  (res1 + res2).grad();
+
+  EXPECT_NEAR(res1.val(), res2.val(), eps);
+  for (int i = 0; i < M; i++) {
+    for (int j = 0; j < N; j++) {
+      EXPECT_NEAR(x1(j, i).adj(), x2(j, i).adj(), eps);
+    }
+  }
+  for (int i = 0; i < M; i++) {
+    EXPECT_NEAR(beta1[i].adj(), beta2[i].adj(), eps);
+  }
+  EXPECT_NEAR(cuts1[0].adj(), cuts2[0].adj(), eps);
+}
+
+TEST(ProbDistributionsOrderedLogisticGLM,
+     glm_matches_ordered_logistic_vars_big) {
+  double eps = 1e-7;
+  int N = 155;
+  int M = 15;
+  int C = 10;
+  vector<int> y(N);
+  for (int i = 0; i < N; i++) {
+    y[i] = Matrix<unsigned int, Dynamic, 1>::Random(1)[0] % (C + 1) + 1;
+  }
+  VectorXd cuts_double = (VectorXd::Random(C).array() + 1) / 2 / C;
+  for (int i = 1; i < C; i++) {
+    cuts_double[i] += cuts_double[i - 1];
+  }
+  Matrix<var, Dynamic, 1> cuts1 = cuts_double, cuts2 = cuts_double;
+  VectorXd beta_double = VectorXd::Random(M);
+  Matrix<var, Dynamic, 1> beta1 = beta_double, beta2 = beta_double;
+  MatrixXd x_double = MatrixXd::Random(N, M);
+  Matrix<var, Dynamic, Dynamic> x1 = x_double, x2 = x_double;
+  var res1 = ordered_logistic_glm_lpmf(y, x1, beta1, cuts1);
+  var res2 = ordered_logistic_glm_simple_lpmf<false>(y, x2, beta2, cuts2);
+  (res1 + res2).grad();
+
+  Matrix<double, Dynamic, Dynamic> x_adj(N, M);
+
+  EXPECT_NEAR(res1.val(), res2.val(), eps);
+  for (int i = 0; i < M; i++) {
+    for (int j = 0; j < N; j++) {
+      x_adj(j, i) = x2(j, i).adj();
+      EXPECT_NEAR(x1(j, i).adj(), x2(j, i).adj(), eps);
+    }
+  }
+  for (int i = 0; i < M; i++) {
+    EXPECT_NEAR(beta1[i].adj(), beta2[i].adj(), eps);
+  }
+  for (int i = 0; i < C; i++) {
+    EXPECT_NEAR(cuts1[i].adj(), cuts2[i].adj(), eps);
+  }
+}
+
+TEST(ProbDistributionsOrderedLogisticGLM, glm_interfaces) {
+  int N = 5;
+  int M = 2;
+  int C = 3;
+  vector<int> y = {1, 1, 2, 4, 4};
+  int y_scal = 1;
+  VectorXd cuts_double(C);
+  cuts_double << 0.9, 1.1, 7;
+  Matrix<var, Dynamic, 1> cuts_var = cuts_double;
+  VectorXd beta_double(M);
+  beta_double << 1.1, 0.4;
+  Matrix<var, Dynamic, 1> beta_var = beta_double;
+  Matrix<double, Dynamic, Dynamic> x_double(N, M);
+  x_double << 1, 2, 3, 4, 5, 6, 7, 8, 9, 0;
+  Matrix<var, Dynamic, Dynamic> x_var = x_double;
+  Matrix<double, 1, Dynamic> x_row_double(M);
+  x_row_double << 1, 2;
+  Matrix<var, 1, Dynamic> x_row_var = x_row_double;
+
+  EXPECT_NO_THROW(
+      ordered_logistic_glm_lpmf(y, x_double, beta_double, cuts_double));
+  EXPECT_NO_THROW(
+      ordered_logistic_glm_lpmf(y, x_var, beta_double, cuts_double));
+  EXPECT_NO_THROW(
+      ordered_logistic_glm_lpmf(y, x_double, beta_var, cuts_double));
+  EXPECT_NO_THROW(
+      ordered_logistic_glm_lpmf(y, x_double, beta_double, cuts_var));
+  EXPECT_NO_THROW(ordered_logistic_glm_lpmf(y, x_double, beta_var, cuts_var));
+  EXPECT_NO_THROW(ordered_logistic_glm_lpmf(y, x_var, beta_double, cuts_var));
+  EXPECT_NO_THROW(ordered_logistic_glm_lpmf(y, x_var, beta_var, cuts_double));
+
+  EXPECT_NO_THROW(
+      ordered_logistic_glm_lpmf(y_scal, x_double, beta_double, cuts_double));
+  EXPECT_NO_THROW(
+      ordered_logistic_glm_lpmf(y_scal, x_var, beta_double, cuts_double));
+  EXPECT_NO_THROW(
+      ordered_logistic_glm_lpmf(y_scal, x_double, beta_var, cuts_double));
+  EXPECT_NO_THROW(
+      ordered_logistic_glm_lpmf(y_scal, x_double, beta_double, cuts_var));
+  EXPECT_NO_THROW(
+      ordered_logistic_glm_lpmf(y_scal, x_double, beta_var, cuts_var));
+  EXPECT_NO_THROW(
+      ordered_logistic_glm_lpmf(y_scal, x_var, beta_double, cuts_var));
+  EXPECT_NO_THROW(
+      ordered_logistic_glm_lpmf(y_scal, x_var, beta_var, cuts_double));
+
+  EXPECT_NO_THROW(
+      ordered_logistic_glm_lpmf(y, x_row_double, beta_double, cuts_double));
+  EXPECT_NO_THROW(
+      ordered_logistic_glm_lpmf(y, x_row_var, beta_double, cuts_double));
+  EXPECT_NO_THROW(
+      ordered_logistic_glm_lpmf(y, x_row_double, beta_var, cuts_double));
+  EXPECT_NO_THROW(
+      ordered_logistic_glm_lpmf(y, x_row_double, beta_double, cuts_var));
+  EXPECT_NO_THROW(
+      ordered_logistic_glm_lpmf(y, x_row_double, beta_var, cuts_var));
+  EXPECT_NO_THROW(
+      ordered_logistic_glm_lpmf(y, x_row_var, beta_double, cuts_var));
+  EXPECT_NO_THROW(
+      ordered_logistic_glm_lpmf(y, x_row_var, beta_var, cuts_double));
+
+  EXPECT_NO_THROW(ordered_logistic_glm_lpmf(y_scal, x_row_double, beta_double,
+                                            cuts_double));
+  EXPECT_NO_THROW(
+      ordered_logistic_glm_lpmf(y_scal, x_row_var, beta_double, cuts_double));
+  EXPECT_NO_THROW(
+      ordered_logistic_glm_lpmf(y_scal, x_row_double, beta_var, cuts_double));
+  EXPECT_NO_THROW(
+      ordered_logistic_glm_lpmf(y_scal, x_row_double, beta_double, cuts_var));
+  EXPECT_NO_THROW(
+      ordered_logistic_glm_lpmf(y_scal, x_row_double, beta_var, cuts_var));
+  EXPECT_NO_THROW(
+      ordered_logistic_glm_lpmf(y_scal, x_row_var, beta_double, cuts_var));
+  EXPECT_NO_THROW(
+      ordered_logistic_glm_lpmf(y_scal, x_row_var, beta_var, cuts_double));
+}
+
+TEST(ProbDistributionsOrderedLogisticGLM, glm_errors) {
+  int N = 5;
+  int M = 2;
+  int C = 3;
+  vector<int> y{1, 1, 2, 4, 4};
+  vector<int> y_size{1, 1, 2, 4, 4, 2};
+  vector<int> y_val1{1, 1, 2, 4, 5};
+  vector<int> y_val2{1, 1, 2, 4, 0};
+  VectorXd cuts(C);
+  cuts << 0.9, 1.1, 7;
+  VectorXd cuts_val1(C);
+  cuts_val1 << 0.9, 1.1, INFINITY;
+  VectorXd cuts_val2(C);
+  cuts_val2 << -INFINITY, 1.1, 4;
+  VectorXd cuts_val3(C);
+  cuts_val3 << 0, 1.1, 0.5;
+  VectorXd cuts_val4(C);
+  cuts_val4 << 0, NAN, 0.5;
+  VectorXd beta(M);
+  beta << 1.1, 0.4;
+  VectorXd beta_size(M + 1);
+  beta_size << 1.1, 0.4, 0;
+  VectorXd beta_val(M);
+  beta_val << 1.1, INFINITY;
+  Matrix<double, Dynamic, Dynamic> x(N, M);
+  x << 1, 2, 3, 4, 5, 6, 7, 8, 9, 0;
+  Matrix<double, Dynamic, Dynamic> x_size1(N + 1, M);
+  x_size1 << 1, 2, 3, 4, 5, 6, 7, 8, 9, 0, 9, 8;
+  Matrix<double, Dynamic, Dynamic> x_size2(N, M + 1);
+  x_size2 << 1, 2, 3, 4, 5, 6, 7, 8, 9, 0, 4, 5, 6, 7, 9;
+  Matrix<double, 1, Dynamic> x_size3(M + 1);
+  x_size3 << 1, 2, 3;
+  Matrix<double, Dynamic, Dynamic> x_val(N, M);
+  x_val << 1, 2, 3, 4, 5, 6, 7, INFINITY, 9, 0;
+
+  EXPECT_NO_THROW(ordered_logistic_glm_lpmf(y, x, beta, cuts));
+
+  EXPECT_THROW(ordered_logistic_glm_lpmf(y_size, x, beta, cuts),
+               std::invalid_argument);
+  EXPECT_THROW(ordered_logistic_glm_lpmf(y, x_size1, beta, cuts),
+               std::invalid_argument);
+  EXPECT_THROW(ordered_logistic_glm_lpmf(y, x_size2, beta, cuts),
+               std::invalid_argument);
+  EXPECT_THROW(ordered_logistic_glm_lpmf(y, x_size3, beta, cuts),
+               std::invalid_argument);
+  EXPECT_THROW(ordered_logistic_glm_lpmf(y, x, beta_size, cuts),
+               std::invalid_argument);
+
+  EXPECT_THROW(ordered_logistic_glm_lpmf(y_val1, x, beta, cuts),
+               std::domain_error);
+  EXPECT_THROW(ordered_logistic_glm_lpmf(y_val2, x, beta, cuts),
+               std::domain_error);
+  EXPECT_THROW(ordered_logistic_glm_lpmf(y, x_val, beta, cuts),
+               std::domain_error);
+  EXPECT_THROW(ordered_logistic_glm_lpmf(y, x, beta_val, cuts),
+               std::domain_error);
+  EXPECT_THROW(ordered_logistic_glm_lpmf(y, x, beta, cuts_val1),
+               std::domain_error);
+  EXPECT_THROW(ordered_logistic_glm_lpmf(y, x, beta, cuts_val2),
+               std::domain_error);
+  EXPECT_THROW(ordered_logistic_glm_lpmf(y, x, beta, cuts_val3),
+               std::domain_error);
+  EXPECT_THROW(ordered_logistic_glm_lpmf(y, x, beta, cuts_val4),
+               std::domain_error);
+}