diff --git a/stan/math/fwd/functor/operands_and_partials.hpp b/stan/math/fwd/functor/operands_and_partials.hpp
index bf8e087d521..e8d53d4a4a6 100644
--- a/stan/math/fwd/functor/operands_and_partials.hpp
+++ b/stan/math/fwd/functor/operands_and_partials.hpp
@@ -27,6 +27,7 @@ class ops_partials_edge<Dx, fvar<Dx>> {
   const Op& operand_;
 
   Dx dx() { return this->partials_[0] * this->operand_.d_; }
+  Dx dx_v() { return this->partials_[0] * this->operand_.d_; }
 };
 }  // namespace internal
 
@@ -106,6 +107,14 @@ class operands_and_partials<Op1, Op2, Op3, Op4, Op5, fvar<Dx>> {
         = edge1_.dx() + edge2_.dx() + edge3_.dx() + edge4_.dx() + edge5_.dx();
     return T_return_type(value, deriv);
   }
+
+  template <typename EigVec, require_eigen_vector_t<EigVec>* = nullptr>
+  auto build(EigVec&& value) {
+    Eigen::Array<fvar<Dx>, -1, 1> ret(value.template cast<fvar<Dx>>());
+    ret.d() = (edge1_.dx_v() + edge2_.dx_v() + edge3_.dx_v() + edge4_.dx_v()
+               + edge5_.dx_v());
+    return ret;
+  }
 };
 
 namespace internal {
@@ -135,6 +144,13 @@ class ops_partials_edge<Dx, std::vector<fvar<Dx>>> {
     }
     return derivative;
   }
+  Eigen::Array<Dx, -1, 1> dx_v() {
+    Eigen::Array<Dx, -1, 1> derivative(this->operands_.size());
+    for (size_t i = 0; i < this->operands_.size(); ++i) {
+      derivative[i] = this->partials_[i] * this->operands_[i].d_;
+    }
+    return derivative;
+  }
 };
 
 template <typename Dx, int R, int C>
@@ -161,6 +177,10 @@ class ops_partials_edge<Dx, Eigen::Matrix<fvar<Dx>, R, C>> {
     }
     return derivative;
   }
+
+  Eigen::Array<Dx, -1, 1> dx_v() {
+    return this->partials_.array() * this->operands_.d().array();
+  }
 };
 
 // Multivariate; vectors of eigen types
@@ -191,6 +211,14 @@ class ops_partials_edge<Dx, std::vector<Eigen::Matrix<fvar<Dx>, R, C>>> {
     }
     return derivative;
   }
+
+  Eigen::Array<Dx, -1, 1> dx_v() {
+    Eigen::Array<Dx, -1, 1> derivative(this->operands_.size());
+    for (size_t i = 0; i < this->operands_.size(); ++i) {
+      derivative[i] = this->partials_vec_[i].dot(this->operands_[i].d());
+    }
+    return derivative;
+  }
 };
 
 template <typename Dx>
@@ -220,6 +248,16 @@ class ops_partials_edge<Dx, std::vector<std::vector<fvar<Dx>>>> {
     }
     return derivative;
   }
+  Eigen::Array<Dx, -1, 1> dx_v() {
+    Eigen::Array<Dx, -1, 1> derivative
+        = Eigen::Array<Dx, -1, 1>::Zero(this->operands_.size());
+    for (size_t i = 0; i < this->operands_.size(); ++i) {
+      for (int j = 0; j < this->operands_[i].size(); ++j) {
+        derivative[i] = this->partials_vec_[i][j] * this->operands_[i][j].d_;
+      }
+    }
+    return derivative;
+  }
 };
 
 }  // namespace internal
diff --git a/stan/math/prim/functor/operands_and_partials.hpp b/stan/math/prim/functor/operands_and_partials.hpp
index abfff95ff9d..e8e3819480b 100644
--- a/stan/math/prim/functor/operands_and_partials.hpp
+++ b/stan/math/prim/functor/operands_and_partials.hpp
@@ -63,6 +63,7 @@ class ops_partials_edge<ViewElt, Op, require_st_arithmetic<Op>> {
    * expression returning zero.
    */
   static constexpr double dx() noexcept { return 0.0; }
+  static constexpr double dx_v() noexcept { return 0.0; }
   /**
    * Return the size of the operand for the edge. For doubles this is a compile
    * time expression returning zero.
@@ -140,7 +141,10 @@ class operands_and_partials {
    * @param value the return value of the function we are compressing
    * @return the value with its derivative
    */
-  inline double build(double value) const noexcept { return value; }
+  template <typename T>
+  inline auto build(T&& value) const noexcept {
+    return std::forward<T>(value);
+  }
 
   // These will always be 0 size base template instantiations (above).
   internal::ops_partials_edge<double, std::decay_t<Op1>> edge1_;
diff --git a/stan/math/prim/functor/prob_reducer.hpp b/stan/math/prim/functor/prob_reducer.hpp
new file mode 100644
index 00000000000..4509b7b85c8
--- /dev/null
+++ b/stan/math/prim/functor/prob_reducer.hpp
@@ -0,0 +1,226 @@
+#ifndef STAN_MATH_PRIM_FUNCTOR_PROB_REDUCER_HPP
+#define STAN_MATH_PRIM_FUNCTOR_PROB_REDUCER_HPP
+
+#include <stan/math/prim/meta.hpp>
+
+namespace stan {
+namespace math {
+
+/**
+ * Used by distributions to decide whether return type shoudl be Scalar or
+ * Vector.
+ */
+enum class ProbReturnType { Scalar, Vector };
+
+/**
+ * For scalars performs summations and is a no-op reducer for eigen vectors.
+ * Used in the probability distributions for scalar or vector return types.
+ */
+template <typename T, typename = void>
+struct prob_reducer;
+
+/**
+ * For scalars performs summations when given eigen types.
+ * @tparam A stan scalar type
+ */
+template <typename T>
+struct prob_reducer<T, require_stan_scalar_t<T>> {
+  T ret_;  // Underlying return type
+
+  /**
+   * Construct from an Eigen type while ignoring size argument passed.
+   * @tparam EigArr A type inheriting from `Eigen::EigenBase`
+   * @tparam Tossed An integral type
+   * @param x will be summed and passed into `ret_`.
+   */
+  template <typename EigArr, typename Tossed,
+            require_eigen_t<EigArr>* = nullptr,
+            require_integral_t<Tossed>* = nullptr>
+  prob_reducer(EigArr&& x, Tossed&& /* */)
+      : ret_(sum(std::forward<EigArr>(x))) {}
+
+  /**
+   * Construct from an Eigen type while ignoring size argument passed.
+   * @tparam Scalar a scalar
+   * @tparam Tossed an integral type
+   * @param x will be summed and inserted into `ret_`.
+   */
+  template <typename Scalar, typename Tossed,
+            require_stan_scalar_t<Scalar>* = nullptr>
+  prob_reducer(Scalar&& x, Tossed&& /* */) : ret_(x) {}
+
+  /**
+   * Perform summation and then assignment
+   * @tparam EigArr A type inheriting from `Eigen::EigenBase`
+   */
+  template <typename EigArr, require_eigen_t<EigArr>* = nullptr>
+  inline auto operator=(EigArr&& x) {
+    ret_ = sum(x);
+    return *this;
+  }
+
+  /**
+   * Assignment
+   * @tparam Scalar A scalar type
+   */
+  template <typename Scalar, require_stan_scalar_t<Scalar>* = nullptr>
+  inline auto operator=(Scalar x) {
+    ret_ = x;
+    return *this;
+  }
+
+  /**
+   * Perform summation and then `+=`
+   * @tparam EigArr A type inheriting from `Eigen::EigenBase`
+   * @param x Eigen object to be summed.
+   */
+  template <typename EigArr, require_eigen_t<EigArr>* = nullptr>
+  inline auto operator+=(EigArr&& x) {
+    ret_ += sum(x);
+    return *this;
+  }
+
+  template <typename Scalar, require_stan_scalar_t<Scalar>* = nullptr>
+  inline auto operator+=(Scalar&& x) {
+    ret_ += x;
+    return *this;
+  }
+
+  /**
+   * Perform summation and then `-=`
+   * @tparam EigArr A type inheriting from `Eigen::EigenBase`
+   * @param x Eigen object to be summed.
+   */
+  template <typename EigArr, require_eigen_t<EigArr>* = nullptr>
+  inline auto operator-=(EigArr&& x) {
+    ret_ -= sum(x);
+    return *this;
+  }
+
+  template <typename Scalar, require_stan_scalar_t<Scalar>* = nullptr>
+  inline auto operator-=(Scalar&& x) {
+    ret_ -= x;
+    return *this;
+  }
+
+  /**
+   * Return the underlying scalar return type.
+   */
+  inline auto ret() noexcept { return ret_; }
+
+  /**
+   * Return a zero value, used when distribution has special cases that
+   *  immedietly return zero.
+   * @tparam Types types to deduce the overall return type of the function.
+   */
+  template <typename... Types>
+  static auto zero(int /* */) {
+    return return_type_t<Types...>(0);
+  }
+};
+
+template <typename T>
+struct prob_reducer<T, require_eigen_t<T>> {
+  T ret_;
+
+  /**
+   * Construct from an Eigen type while ignoring size argument passed.
+   * @tparam EigArr A type inheriting from `Eigen::EigenBase`
+   * @tparam Size An integral type
+   * @param x will be forwarded to `ret_`.
+   */
+  template <typename EigArr, typename Size, require_eigen_t<EigArr>* = nullptr,
+            require_integral_t<Size>* = nullptr>
+  prob_reducer(EigArr&& x, Size /* x */) : ret_(std::forward<EigArr>(x)) {}
+
+  /**
+   * Construct from a scalar type.
+   * @tparam Scalar a scalar
+   * @tparam Size An integral type
+   * @param x passed to `ret_` along with size to fill with a base value.
+   * @param n The size `ret_` should be
+   */
+  template <typename Scalar, typename Size,
+            require_stan_scalar_t<Scalar>* = nullptr,
+            require_integral_t<Size>* = nullptr>
+  prob_reducer(Scalar constant, Size n) : ret_(T::Constant(n, constant)) {}
+
+  /**
+   * assignment
+   * @tparam EigArr A type inheriting from `Eigen::EigenBase`
+   */
+  template <typename EigArr, require_eigen_t<EigArr>* = nullptr>
+  inline auto operator=(EigArr&& x) {
+    ret_ = std::forward<EigArr>(x);
+    return *this;
+  }
+
+  /**
+   * assignm scalar by propogating value over `ret_`
+   * @tparam Scalar a stan scalar
+   * @param x The value to fill `ret_` with.
+   */
+  template <typename Scalar, require_stan_scalar_t<Scalar>* = nullptr>
+  inline auto operator=(Scalar x) {
+    ret_ = Eigen::Array<value_type_t<T>, -1, 1>::Constant(x, ret_.size());
+    return *this;
+  }
+
+  template <typename EigArr, require_eigen_t<EigArr>* = nullptr>
+  inline auto operator+=(EigArr&& x) {
+    ret_ += std::forward<EigArr>(x);
+    return *this;
+  }
+
+  template <typename Scalar, require_stan_scalar_t<Scalar>* = nullptr>
+  inline auto operator+=(Scalar&& x) {
+    ret_ += x;
+    return *this;
+  }
+
+  template <typename EigArr, require_eigen_t<EigArr>* = nullptr>
+  inline auto operator-=(EigArr&& x) {
+    ret_ -= std::forward<EigArr>(x);
+    return *this;
+  }
+
+  template <typename Scalar, require_stan_scalar_t<Scalar>* = nullptr>
+  inline auto operator-=(Scalar&& x) {
+    ret_ -= x;
+    return *this;
+  }
+
+  /**
+   * Return the underlying scalar return type. Passed the underlying by
+   * moving it which can cause `ret_` to be uninitialized after.
+   */
+  inline auto&& ret() noexcept { return std::move(ret_); }
+
+  /**
+   * Return a zero value, used when distribution has special cases that
+   *  immedietly return zero.
+   * @tparam Types types to deduce the overall return type of the function.
+   * @param size The size of the array to return
+   */
+  template <typename... Types>
+  static auto zero(int size) {
+    return Eigen::Array<return_type_t<Types...>, -1, 1>::Constant(0, size)
+        .eval();
+  }
+};
+
+/**
+ * Generate a reducer with correct return type.
+ * @tparam ReturnType Either Scalar or Vector.
+ * @tparam Types A parameter pack of types to deduce the underlying scalar type
+ * from
+ */
+template <ProbReturnType ReturnType, typename... Types>
+using prob_return_t = prob_reducer<std::conditional_t<
+    ReturnType == ProbReturnType::Scalar, return_type_t<Types...>,
+    Eigen::Array<return_type_t<Types...>, -1, 1>>>;
+
+}  // namespace math
+}  // namespace stan
+
+#endif
diff --git a/stan/math/prim/prob/normal_log.hpp b/stan/math/prim/prob/normal_log.hpp
index 962d6657f08..bb16d728ae9 100644
--- a/stan/math/prim/prob/normal_log.hpp
+++ b/stan/math/prim/prob/normal_log.hpp
@@ -29,20 +29,16 @@ namespace math {
  * @tparam T_loc Type of location parameter.
  */
 template <bool propto, typename T_y, typename T_loc, typename T_scale>
-inline return_type_t<T_y, T_loc, T_scale> normal_log(const T_y& y,
-                                                     const T_loc& mu,
-                                                     const T_scale& sigma) {
-  return normal_lpdf<propto, T_y, T_loc, T_scale>(y, mu, sigma);
+inline auto normal_log(const T_y& y, const T_loc& mu, const T_scale& sigma) {
+  return normal_lpdf<propto>(y, mu, sigma);
 }
 
 /** \ingroup prob_dists
  * @deprecated use <code>normal_lpdf</code>
  */
 template <typename T_y, typename T_loc, typename T_scale>
-inline return_type_t<T_y, T_loc, T_scale> normal_log(const T_y& y,
-                                                     const T_loc& mu,
-                                                     const T_scale& sigma) {
-  return normal_lpdf<T_y, T_loc, T_scale>(y, mu, sigma);
+inline auto normal_log(const T_y& y, const T_loc& mu, const T_scale& sigma) {
+  return normal_lpdf<false>(y, mu, sigma);
 }
 
 }  // namespace math
diff --git a/stan/math/prim/prob/normal_lpdf.hpp b/stan/math/prim/prob/normal_lpdf.hpp
index 5d688964549..1a8bc855c1c 100644
--- a/stan/math/prim/prob/normal_lpdf.hpp
+++ b/stan/math/prim/prob/normal_lpdf.hpp
@@ -14,6 +14,7 @@
 #include <stan/math/prim/fun/to_ref.hpp>
 #include <stan/math/prim/fun/value_of.hpp>
 #include <stan/math/prim/functor/operands_and_partials.hpp>
+#include <stan/math/prim/functor/prob_reducer.hpp>
 #include <cmath>
 
 namespace stan {
@@ -38,12 +39,11 @@ namespace math {
  * @return The log of the product of the densities.
  * @throw std::domain_error if the scale is not positive.
  */
-template <bool propto, typename T_y, typename T_loc, typename T_scale,
+template <bool propto, ProbReturnType RetType = ProbReturnType::Scalar,
+          typename T_y, typename T_loc, typename T_scale,
           require_all_not_nonscalar_prim_or_rev_kernel_expression_t<
               T_y, T_loc, T_scale>* = nullptr>
-inline return_type_t<T_y, T_loc, T_scale> normal_lpdf(const T_y& y,
-                                                      const T_loc& mu,
-                                                      const T_scale& sigma) {
+inline auto normal_lpdf(const T_y& y, const T_loc& mu, const T_scale& sigma) {
   using T_partials_return = partials_return_t<T_y, T_loc, T_scale>;
   using T_y_ref = ref_type_if_t<!is_constant<T_y>::value, T_y>;
   using T_mu_ref = ref_type_if_t<!is_constant<T_loc>::value, T_loc>;
@@ -62,12 +62,13 @@ inline return_type_t<T_y, T_loc, T_scale> normal_lpdf(const T_y& y,
   check_not_nan(function, "Random variable", y_val);
   check_finite(function, "Location parameter", mu_val);
   check_positive(function, "Scale parameter", sigma_val);
-
+  using ret_t = prob_return_t<RetType, T_partials_return>;
+  const size_t N = max_size(y, mu, sigma);
   if (size_zero(y, mu, sigma)) {
-    return 0.0;
+    return ret_t::template zero<T_y, T_loc, T_scale>(N);
   }
   if (!include_summand<propto, T_y, T_loc, T_scale>::value) {
-    return 0.0;
+    return ret_t::template zero<T_y, T_loc, T_scale>(N);
   }
 
   operands_and_partials<T_y_ref, T_mu_ref, T_sigma_ref> ops_partials(
@@ -79,13 +80,20 @@ inline return_type_t<T_y, T_loc, T_scale> normal_lpdf(const T_y& y,
   const auto& y_scaled_sq
       = to_ref_if<!is_constant_all<T_scale>::value>(y_scaled * y_scaled);
 
-  size_t N = max_size(y, mu, sigma);
-  T_partials_return logp = -0.5 * sum(y_scaled_sq);
+  prob_return_t<RetType, T_partials_return> logp(-0.5 * y_scaled_sq, N);
   if (include_summand<propto>::value) {
-    logp += NEG_LOG_SQRT_TWO_PI * N;
+    if (RetType == ProbReturnType::Scalar) {
+      logp += NEG_LOG_SQRT_TWO_PI * N;
+    } else {
+      logp += NEG_LOG_SQRT_TWO_PI;
+    }
   }
   if (include_summand<propto, T_scale>::value) {
-    logp -= sum(log(sigma_val)) * N / math::size(sigma);
+    if (RetType == ProbReturnType::Scalar) {
+      logp -= sum(log(sigma_val)) * N / math::size(sigma);
+    } else {
+      logp -= log(sigma_val);
+    }
   }
 
   if (!is_constant_all<T_y, T_scale, T_loc>::value) {
@@ -103,13 +111,11 @@ inline return_type_t<T_y, T_loc, T_scale> normal_lpdf(const T_y& y,
       ops_partials.edge2_.partials_ = std::move(scaled_diff);
     }
   }
-  return ops_partials.build(logp);
+  return ops_partials.build(logp.ret());
 }
 
 template <typename T_y, typename T_loc, typename T_scale>
-inline return_type_t<T_y, T_loc, T_scale> normal_lpdf(const T_y& y,
-                                                      const T_loc& mu,
-                                                      const T_scale& sigma) {
+inline auto normal_lpdf(const T_y& y, const T_loc& mu, const T_scale& sigma) {
   return normal_lpdf<false>(y, mu, sigma);
 }
 
diff --git a/stan/math/rev/functor/operands_and_partials.hpp b/stan/math/rev/functor/operands_and_partials.hpp
index 908d5803174..fec51b1725c 100644
--- a/stan/math/rev/functor/operands_and_partials.hpp
+++ b/stan/math/rev/functor/operands_and_partials.hpp
@@ -48,12 +48,25 @@ inline void update_adjoints(var_value<T1>& x, const T2& y, const vari& z) {
   x.adj() += z.adj() * y;
 }
 
+template <typename T1, typename T2,
+          require_all_kernel_expressions_and_none_scalar_t<T1, T2>* = nullptr>
+inline void update_adjoints(var_value<T1>& x, const T2& y, const var& z) {
+  x.adj() += z.adj() * y;
+}
+
 template <typename Scalar1, typename Scalar2, require_var_t<Scalar1>* = nullptr,
           require_not_var_matrix_t<Scalar1>* = nullptr,
           require_arithmetic_t<Scalar2>* = nullptr>
 inline void update_adjoints(Scalar1 x, Scalar2 y, const vari& z) noexcept {
   x.adj() += z.adj() * y;
 }
+template <typename Scalar1, typename Scalar2, require_var_t<Scalar1>* = nullptr,
+          require_not_var_matrix_t<Scalar1>* = nullptr,
+          require_arithmetic_t<Scalar2>* = nullptr>
+inline void update_adjoints(Scalar1 x, Scalar2 y, const var& z) noexcept {
+  x.adj() += z.adj() * y;
+}
+
 template <typename Matrix1, typename Matrix2,
           require_rev_matrix_t<Matrix1>* = nullptr,
           require_st_arithmetic<Matrix2>* = nullptr>
@@ -61,10 +74,21 @@ inline void update_adjoints(Matrix1& x, const Matrix2& y, const vari& z) {
   x.adj().array() += z.adj() * y.array();
 }
 
+template <typename Matrix1, typename Matrix2,
+          require_rev_matrix_t<Matrix1>* = nullptr,
+          require_st_arithmetic<Matrix2>* = nullptr>
+inline void update_adjoints(Matrix1& x, const Matrix2& y, const var& z) {
+  x.adj().array() += z.adj() * y.array();
+}
+
 template <typename Arith, typename Alt, require_st_arithmetic<Arith>* = nullptr>
 inline constexpr void update_adjoints(Arith&& /* x */, Alt&& /* y */,
                                       const vari& /* z */) noexcept {}
 
+template <typename Arith, typename Alt, require_st_arithmetic<Arith>* = nullptr>
+inline constexpr void update_adjoints(Arith&& /* x */, Alt&& /* y */,
+                                      const var& /* z */) noexcept {}
+
 template <typename StdVec1, typename Vec2,
           require_std_vector_t<StdVec1>* = nullptr,
           require_st_arithmetic<Vec2>* = nullptr>
@@ -74,6 +98,54 @@ inline void update_adjoints(StdVec1& x, const Vec2& y, const vari& z) {
   }
 }
 
+template <typename StdVec1, typename Vec2,
+          require_std_vector_t<StdVec1>* = nullptr,
+          require_st_arithmetic<Vec2>* = nullptr>
+inline void update_adjoints(StdVec1& x, const Vec2& y, const var& z) {
+  for (size_t i = 0; i < x.size(); ++i) {
+    update_adjoints(x[i], y[i], z);
+  }
+}
+
+template <
+    typename T1, typename T2, typename T3,
+    require_all_kernel_expressions_and_none_scalar_t<T1, T2, T3>* = nullptr>
+inline void update_adjoints(var_value<T1>& x, const T2& y, const T3& z) {
+  x.adj() += z.adj() * y;
+}
+
+template <typename Scalar1, typename Scalar2, typename T3,
+          require_var_t<Scalar1>* = nullptr,
+          require_not_var_matrix_t<Scalar1>* = nullptr,
+          require_arithmetic_t<Scalar2>* = nullptr,
+          require_eigen_t<T3>* = nullptr>
+inline void update_adjoints(Scalar1 x, Scalar2 y, const T3& z) noexcept {
+  x.adj() += sum(z.adj() * y);
+}
+template <typename Matrix1, typename Matrix2, typename T3,
+          require_rev_matrix_t<Matrix1>* = nullptr,
+          require_st_arithmetic<Matrix2>* = nullptr,
+          require_eigen_t<T3>* = nullptr>
+inline void update_adjoints(Matrix1& x, const Matrix2& y, const T3& z) {
+  x.adj().array() += z.adj() * y.array();
+}
+
+template <typename Arith, typename Alt, typename T3,
+          require_st_arithmetic<Arith>* = nullptr,
+          require_eigen_t<T3>* = nullptr>
+inline constexpr void update_adjoints(Arith&& /* x */, Alt&& /* y */,
+                                      const T3& /* z */) noexcept {}
+
+template <typename StdVec1, typename Vec2, typename T3,
+          require_std_vector_t<StdVec1>* = nullptr,
+          require_st_arithmetic<Vec2>* = nullptr,
+          require_eigen_t<T3>* = nullptr>
+inline void update_adjoints(StdVec1& x, const Vec2& y, const T3& z) {
+  for (size_t i = 0; i < x.size(); ++i) {
+    update_adjoints(x[i], y[i], z[i]);
+  }
+}
+
 }  // namespace internal
 
 /** \ingroup type_trait
@@ -169,6 +241,33 @@ class operands_and_partials<Op1, Op2, Op3, Op4, Op5, var> {
           }
         });
   }
+  template <typename T, require_eigen_t<T>* = nullptr>
+  auto build(T&& value) {
+    arena_t<promote_scalar_t<var, T>> ret = value;
+    reverse_pass_callback(
+        [ret, operand1 = edge1_.operand(), partial1 = edge1_.partial(),
+         operand2 = edge2_.operand(), partial2 = edge2_.partial(),
+         operand3 = edge3_.operand(), partial3 = edge3_.partial(),
+         operand4 = edge4_.operand(), partial4 = edge4_.partial(),
+         operand5 = edge5_.operand(), partial5 = edge5_.partial()]() mutable {
+          if (!is_constant<Op1>::value) {
+            internal::update_adjoints(operand1, partial1, ret);
+          }
+          if (!is_constant<Op2>::value) {
+            internal::update_adjoints(operand2, partial2, ret);
+          }
+          if (!is_constant<Op3>::value) {
+            internal::update_adjoints(operand3, partial3, ret);
+          }
+          if (!is_constant<Op4>::value) {
+            internal::update_adjoints(operand4, partial4, ret);
+          }
+          if (!is_constant<Op5>::value) {
+            internal::update_adjoints(operand5, partial5, ret);
+          }
+        });
+    return plain_type_t<promote_scalar_t<var, T>>(ret);
+  }
 };
 
 namespace internal {
diff --git a/test/unit/math/fwd/prob/normal_vectorized_test.cpp b/test/unit/math/fwd/prob/normal_vectorized_test.cpp
new file mode 100644
index 00000000000..b10f2fc833b
--- /dev/null
+++ b/test/unit/math/fwd/prob/normal_vectorized_test.cpp
@@ -0,0 +1,15 @@
+#include <stan/math/fwd.hpp>
+#include <gtest/gtest.h>
+
+TEST(ProbNormal, fvar_test_vlpdf) {
+  using stan::math::fvar;
+  Eigen::Matrix<fvar<double>, -1, 1> Y
+      = Eigen::Matrix<double, -1, 1>::Random(5);
+  Eigen::Matrix<fvar<double>, -1, 1> Mu
+      = Eigen::Matrix<double, -1, 1>::Random(5);
+  Eigen::Matrix<fvar<double>, -1, 1> Sigma
+      = stan::math::abs(Eigen::Matrix<double, -1, 1>::Random(5));
+  Eigen::Matrix<fvar<double>, -1, 1> A
+      = stan::math::normal_lpdf<false, stan::math::ProbReturnType::Vector>(
+          Y, Mu, Sigma);
+}
diff --git a/test/unit/math/mix/prob/normal_test.cpp b/test/unit/math/mix/prob/normal_test.cpp
index 81136b202d2..156f498e3ab 100644
--- a/test/unit/math/mix/prob/normal_test.cpp
+++ b/test/unit/math/mix/prob/normal_test.cpp
@@ -1,5 +1,5 @@
 #include <stan/math/mix.hpp>
-#include <test/unit/math/test_ad.hpp>
+#include <test/unit/math/mix/prob/test_distribution_ad.hpp>
 
 TEST(mathMixScalFun, normal_lpdf) {
   auto f = [](const double mu, const double sigma) {
@@ -10,3 +10,14 @@ TEST(mathMixScalFun, normal_lpdf) {
   stan::test::expect_ad(f(0, 1), 0.0);
   stan::test::expect_ad(f(0, 1), 1.7);
 }
+
+TEST(mathMixScalFun, vnormal_lpdf) {
+  auto f = [](const auto& y, const auto& mu, const auto& sigma) {
+    return stan::math::normal_lpdf<true, stan::math::ProbReturnType::Vector>(
+        y, mu, sigma);
+  };
+  Eigen::VectorXd y = Eigen::VectorXd::Random(5);
+  Eigen::VectorXd mu = Eigen::VectorXd::Random(5);
+  Eigen::VectorXd sigma = stan::math::abs(Eigen::VectorXd::Random(5));
+  stan::test::expect_ad_distribution(f, y, mu, sigma);
+}
diff --git a/test/unit/math/mix/prob/test_distribution_ad.hpp b/test/unit/math/mix/prob/test_distribution_ad.hpp
new file mode 100644
index 00000000000..70fc16fc5d0
--- /dev/null
+++ b/test/unit/math/mix/prob/test_distribution_ad.hpp
@@ -0,0 +1,282 @@
+#ifndef TEST_UNIT_MATH_TEST_DISTRIBUTION_AD_HPP
+#define TEST_UNIT_MATH_TEST_DISTRIBUTION_AD_HPP
+
+#include <stan/math/mix.hpp>
+#include <test/unit/math/test_ad.hpp>
+namespace stan {
+namespace test {
+template <typename Cast, typename T, stan::require_eigen_t<Cast>* = nullptr>
+Cast do_cast(T&& x) {
+  return x;
+}
+template <typename Cast, typename T,
+          stan::require_std_vector_t<Cast>* = nullptr>
+Cast do_cast(T&& x) {
+  return Cast(x.data(), x.data() + x.size());
+}
+template <typename Cast, typename T,
+          stan::require_arithmetic_t<Cast>* = nullptr>
+Cast do_cast(T&& x) {
+  return x(0);
+}
+
+using unary_base_set
+    = std::tuple<std::tuple<double>, std::tuple<std::vector<double>>,
+                 std::tuple<Eigen::Matrix<double, Eigen::Dynamic, 1>>>;
+
+using binary_base_set = std::tuple<
+    std::tuple<double, double>, std::tuple<double, std::vector<double>>,
+    std::tuple<double, Eigen::Matrix<double, Eigen::Dynamic, 1>>,
+    std::tuple<std::vector<double>, double>,
+    std::tuple<std::vector<double>, std::vector<double>>,
+    std::tuple<std::vector<double>, Eigen::Matrix<double, Eigen::Dynamic, 1>>,
+    std::tuple<Eigen::Matrix<double, Eigen::Dynamic, 1>, double>,
+    std::tuple<Eigen::Matrix<double, Eigen::Dynamic, 1>, std::vector<double>>,
+    std::tuple<Eigen::Matrix<double, Eigen::Dynamic, 1>,
+               Eigen::Matrix<double, Eigen::Dynamic, 1>>>;
+
+using ternary_base_sets = std::tuple<
+    std::tuple<double, double, double>,
+    std::tuple<double, double, std::vector<double>>,
+    std::tuple<double, double, Eigen::Matrix<double, Eigen::Dynamic, 1>>,
+    std::tuple<double, std::vector<double>, double>,
+    std::tuple<double, std::vector<double>, std::vector<double>>,
+    std::tuple<double, std::vector<double>,
+               Eigen::Matrix<double, Eigen::Dynamic, 1>>,
+    std::tuple<double, Eigen::Matrix<double, Eigen::Dynamic, 1>, double>,
+    std::tuple<double, Eigen::Matrix<double, Eigen::Dynamic, 1>,
+               std::vector<double>>,
+    std::tuple<double, Eigen::Matrix<double, Eigen::Dynamic, 1>,
+               Eigen::Matrix<double, Eigen::Dynamic, 1>>,
+    std::tuple<std::vector<double>, double, double>,
+    std::tuple<std::vector<double>, double, std::vector<double>>,
+    std::tuple<std::vector<double>, double,
+               Eigen::Matrix<double, Eigen::Dynamic, 1>>,
+    std::tuple<std::vector<double>, std::vector<double>, double>,
+    std::tuple<std::vector<double>, std::vector<double>, std::vector<double>>,
+    std::tuple<std::vector<double>, std::vector<double>,
+               Eigen::Matrix<double, Eigen::Dynamic, 1>>,
+    std::tuple<std::vector<double>, Eigen::Matrix<double, Eigen::Dynamic, 1>,
+               double>,
+    std::tuple<std::vector<double>, Eigen::Matrix<double, Eigen::Dynamic, 1>,
+               std::vector<double>>,
+    std::tuple<std::vector<double>, Eigen::Matrix<double, Eigen::Dynamic, 1>,
+               Eigen::Matrix<double, Eigen::Dynamic, 1>>,
+    std::tuple<Eigen::Matrix<double, Eigen::Dynamic, 1>, double, double>,
+    std::tuple<Eigen::Matrix<double, Eigen::Dynamic, 1>, double,
+               std::vector<double>>,
+    std::tuple<Eigen::Matrix<double, Eigen::Dynamic, 1>, double,
+               Eigen::Matrix<double, Eigen::Dynamic, 1>>,
+    std::tuple<Eigen::Matrix<double, Eigen::Dynamic, 1>, std::vector<double>,
+               double>,
+    std::tuple<Eigen::Matrix<double, Eigen::Dynamic, 1>, std::vector<double>,
+               std::vector<double>>,
+    std::tuple<Eigen::Matrix<double, Eigen::Dynamic, 1>, std::vector<double>,
+               Eigen::Matrix<double, Eigen::Dynamic, 1>>,
+    std::tuple<Eigen::Matrix<double, Eigen::Dynamic, 1>,
+               Eigen::Matrix<double, Eigen::Dynamic, 1>, double>,
+    std::tuple<Eigen::Matrix<double, Eigen::Dynamic, 1>,
+               Eigen::Matrix<double, Eigen::Dynamic, 1>, std::vector<double>>,
+    std::tuple<Eigen::Matrix<double, Eigen::Dynamic, 1>,
+               Eigen::Matrix<double, Eigen::Dynamic, 1>,
+               Eigen::Matrix<double, Eigen::Dynamic, 1>>>;
+
+using quad_base_set = std::tuple<
+    std::tuple<double, double, double, double>,
+    std::tuple<double, double, double, std::vector<double>>,
+    std::tuple<double, double, double,
+               Eigen::Matrix<double, Eigen::Dynamic, 1>>,
+    std::tuple<double, double, std::vector<double>, double>,
+    std::tuple<double, double, std::vector<double>, std::vector<double>>,
+    std::tuple<double, double, std::vector<double>,
+               Eigen::Matrix<double, Eigen::Dynamic, 1>>,
+    std::tuple<double, double, Eigen::Matrix<double, Eigen::Dynamic, 1>,
+               double>,
+    std::tuple<double, double, Eigen::Matrix<double, Eigen::Dynamic, 1>,
+               std::vector<double>>,
+    std::tuple<double, double, Eigen::Matrix<double, Eigen::Dynamic, 1>,
+               Eigen::Matrix<double, Eigen::Dynamic, 1>>,
+    std::tuple<double, std::vector<double>, double, double>,
+    std::tuple<double, std::vector<double>, double, std::vector<double>>,
+    std::tuple<double, std::vector<double>, double,
+               Eigen::Matrix<double, Eigen::Dynamic, 1>>,
+    std::tuple<double, std::vector<double>, std::vector<double>, double>,
+    std::tuple<double, std::vector<double>, std::vector<double>,
+               std::vector<double>>,
+    std::tuple<double, std::vector<double>, std::vector<double>,
+               Eigen::Matrix<double, Eigen::Dynamic, 1>>,
+    std::tuple<double, std::vector<double>,
+               Eigen::Matrix<double, Eigen::Dynamic, 1>, double>,
+    std::tuple<double, std::vector<double>,
+               Eigen::Matrix<double, Eigen::Dynamic, 1>, std::vector<double>>,
+    std::tuple<double, std::vector<double>,
+               Eigen::Matrix<double, Eigen::Dynamic, 1>,
+               Eigen::Matrix<double, Eigen::Dynamic, 1>>,
+    std::tuple<double, Eigen::Matrix<double, Eigen::Dynamic, 1>, double,
+               double>,
+    std::tuple<double, Eigen::Matrix<double, Eigen::Dynamic, 1>, double,
+               std::vector<double>>,
+    std::tuple<double, Eigen::Matrix<double, Eigen::Dynamic, 1>, double,
+               Eigen::Matrix<double, Eigen::Dynamic, 1>>,
+    std::tuple<double, Eigen::Matrix<double, Eigen::Dynamic, 1>,
+               std::vector<double>, double>,
+    std::tuple<double, Eigen::Matrix<double, Eigen::Dynamic, 1>,
+               std::vector<double>, std::vector<double>>,
+    std::tuple<double, Eigen::Matrix<double, Eigen::Dynamic, 1>,
+               std::vector<double>, Eigen::Matrix<double, Eigen::Dynamic, 1>>,
+    std::tuple<double, Eigen::Matrix<double, Eigen::Dynamic, 1>,
+               Eigen::Matrix<double, Eigen::Dynamic, 1>, double>,
+    std::tuple<double, Eigen::Matrix<double, Eigen::Dynamic, 1>,
+               Eigen::Matrix<double, Eigen::Dynamic, 1>, std::vector<double>>,
+    std::tuple<double, Eigen::Matrix<double, Eigen::Dynamic, 1>,
+               Eigen::Matrix<double, Eigen::Dynamic, 1>,
+               Eigen::Matrix<double, Eigen::Dynamic, 1>>,
+    std::tuple<std::vector<double>, double, double, double>,
+    std::tuple<std::vector<double>, double, double, std::vector<double>>,
+    std::tuple<std::vector<double>, double, double,
+               Eigen::Matrix<double, Eigen::Dynamic, 1>>,
+    std::tuple<std::vector<double>, double, std::vector<double>, double>,
+    std::tuple<std::vector<double>, double, std::vector<double>,
+               std::vector<double>>,
+    std::tuple<std::vector<double>, double, std::vector<double>,
+               Eigen::Matrix<double, Eigen::Dynamic, 1>>,
+    std::tuple<std::vector<double>, double,
+               Eigen::Matrix<double, Eigen::Dynamic, 1>, double>,
+    std::tuple<std::vector<double>, double,
+               Eigen::Matrix<double, Eigen::Dynamic, 1>, std::vector<double>>,
+    std::tuple<std::vector<double>, double,
+               Eigen::Matrix<double, Eigen::Dynamic, 1>,
+               Eigen::Matrix<double, Eigen::Dynamic, 1>>,
+    std::tuple<std::vector<double>, std::vector<double>, double, double>,
+    std::tuple<std::vector<double>, std::vector<double>, double,
+               std::vector<double>>,
+    std::tuple<std::vector<double>, std::vector<double>, double,
+               Eigen::Matrix<double, Eigen::Dynamic, 1>>,
+    std::tuple<std::vector<double>, std::vector<double>, std::vector<double>,
+               double>,
+    std::tuple<std::vector<double>, std::vector<double>, std::vector<double>,
+               std::vector<double>>,
+    std::tuple<std::vector<double>, std::vector<double>, std::vector<double>,
+               Eigen::Matrix<double, Eigen::Dynamic, 1>>,
+    std::tuple<std::vector<double>, std::vector<double>,
+               Eigen::Matrix<double, Eigen::Dynamic, 1>, double>,
+    std::tuple<std::vector<double>, std::vector<double>,
+               Eigen::Matrix<double, Eigen::Dynamic, 1>, std::vector<double>>,
+    std::tuple<std::vector<double>, std::vector<double>,
+               Eigen::Matrix<double, Eigen::Dynamic, 1>,
+               Eigen::Matrix<double, Eigen::Dynamic, 1>>,
+    std::tuple<std::vector<double>, Eigen::Matrix<double, Eigen::Dynamic, 1>,
+               double, double>,
+    std::tuple<std::vector<double>, Eigen::Matrix<double, Eigen::Dynamic, 1>,
+               double, std::vector<double>>,
+    std::tuple<std::vector<double>, Eigen::Matrix<double, Eigen::Dynamic, 1>,
+               double, Eigen::Matrix<double, Eigen::Dynamic, 1>>,
+    std::tuple<std::vector<double>, Eigen::Matrix<double, Eigen::Dynamic, 1>,
+               std::vector<double>, double>,
+    std::tuple<std::vector<double>, Eigen::Matrix<double, Eigen::Dynamic, 1>,
+               std::vector<double>, std::vector<double>>,
+    std::tuple<std::vector<double>, Eigen::Matrix<double, Eigen::Dynamic, 1>,
+               std::vector<double>, Eigen::Matrix<double, Eigen::Dynamic, 1>>,
+    std::tuple<std::vector<double>, Eigen::Matrix<double, Eigen::Dynamic, 1>,
+               Eigen::Matrix<double, Eigen::Dynamic, 1>, double>,
+    std::tuple<std::vector<double>, Eigen::Matrix<double, Eigen::Dynamic, 1>,
+               Eigen::Matrix<double, Eigen::Dynamic, 1>, std::vector<double>>,
+    std::tuple<std::vector<double>, Eigen::Matrix<double, Eigen::Dynamic, 1>,
+               Eigen::Matrix<double, Eigen::Dynamic, 1>,
+               Eigen::Matrix<double, Eigen::Dynamic, 1>>,
+    std::tuple<Eigen::Matrix<double, Eigen::Dynamic, 1>, double, double,
+               double>,
+    std::tuple<Eigen::Matrix<double, Eigen::Dynamic, 1>, double, double,
+               std::vector<double>>,
+    std::tuple<Eigen::Matrix<double, Eigen::Dynamic, 1>, double, double,
+               Eigen::Matrix<double, Eigen::Dynamic, 1>>,
+    std::tuple<Eigen::Matrix<double, Eigen::Dynamic, 1>, double,
+               std::vector<double>, double>,
+    std::tuple<Eigen::Matrix<double, Eigen::Dynamic, 1>, double,
+               std::vector<double>, std::vector<double>>,
+    std::tuple<Eigen::Matrix<double, Eigen::Dynamic, 1>, double,
+               std::vector<double>, Eigen::Matrix<double, Eigen::Dynamic, 1>>,
+    std::tuple<Eigen::Matrix<double, Eigen::Dynamic, 1>, double,
+               Eigen::Matrix<double, Eigen::Dynamic, 1>, double>,
+    std::tuple<Eigen::Matrix<double, Eigen::Dynamic, 1>, double,
+               Eigen::Matrix<double, Eigen::Dynamic, 1>, std::vector<double>>,
+    std::tuple<Eigen::Matrix<double, Eigen::Dynamic, 1>, double,
+               Eigen::Matrix<double, Eigen::Dynamic, 1>,
+               Eigen::Matrix<double, Eigen::Dynamic, 1>>,
+    std::tuple<Eigen::Matrix<double, Eigen::Dynamic, 1>, std::vector<double>,
+               double, double>,
+    std::tuple<Eigen::Matrix<double, Eigen::Dynamic, 1>, std::vector<double>,
+               double, std::vector<double>>,
+    std::tuple<Eigen::Matrix<double, Eigen::Dynamic, 1>, std::vector<double>,
+               double, Eigen::Matrix<double, Eigen::Dynamic, 1>>,
+    std::tuple<Eigen::Matrix<double, Eigen::Dynamic, 1>, std::vector<double>,
+               std::vector<double>, double>,
+    std::tuple<Eigen::Matrix<double, Eigen::Dynamic, 1>, std::vector<double>,
+               std::vector<double>, std::vector<double>>,
+    std::tuple<Eigen::Matrix<double, Eigen::Dynamic, 1>, std::vector<double>,
+               std::vector<double>, Eigen::Matrix<double, Eigen::Dynamic, 1>>,
+    std::tuple<Eigen::Matrix<double, Eigen::Dynamic, 1>, std::vector<double>,
+               Eigen::Matrix<double, Eigen::Dynamic, 1>, double>,
+    std::tuple<Eigen::Matrix<double, Eigen::Dynamic, 1>, std::vector<double>,
+               Eigen::Matrix<double, Eigen::Dynamic, 1>, std::vector<double>>,
+    std::tuple<Eigen::Matrix<double, Eigen::Dynamic, 1>, std::vector<double>,
+               Eigen::Matrix<double, Eigen::Dynamic, 1>,
+               Eigen::Matrix<double, Eigen::Dynamic, 1>>,
+    std::tuple<Eigen::Matrix<double, Eigen::Dynamic, 1>,
+               Eigen::Matrix<double, Eigen::Dynamic, 1>, double, double>,
+    std::tuple<Eigen::Matrix<double, Eigen::Dynamic, 1>,
+               Eigen::Matrix<double, Eigen::Dynamic, 1>, double,
+               std::vector<double>>,
+    std::tuple<Eigen::Matrix<double, Eigen::Dynamic, 1>,
+               Eigen::Matrix<double, Eigen::Dynamic, 1>, double,
+               Eigen::Matrix<double, Eigen::Dynamic, 1>>,
+    std::tuple<Eigen::Matrix<double, Eigen::Dynamic, 1>,
+               Eigen::Matrix<double, Eigen::Dynamic, 1>, std::vector<double>,
+               double>,
+    std::tuple<Eigen::Matrix<double, Eigen::Dynamic, 1>,
+               Eigen::Matrix<double, Eigen::Dynamic, 1>, std::vector<double>,
+               std::vector<double>>,
+    std::tuple<Eigen::Matrix<double, Eigen::Dynamic, 1>,
+               Eigen::Matrix<double, Eigen::Dynamic, 1>, std::vector<double>,
+               Eigen::Matrix<double, Eigen::Dynamic, 1>>,
+    std::tuple<Eigen::Matrix<double, Eigen::Dynamic, 1>,
+               Eigen::Matrix<double, Eigen::Dynamic, 1>,
+               Eigen::Matrix<double, Eigen::Dynamic, 1>, double>,
+    std::tuple<Eigen::Matrix<double, Eigen::Dynamic, 1>,
+               Eigen::Matrix<double, Eigen::Dynamic, 1>,
+               Eigen::Matrix<double, Eigen::Dynamic, 1>, std::vector<double>>,
+    std::tuple<Eigen::Matrix<double, Eigen::Dynamic, 1>,
+               Eigen::Matrix<double, Eigen::Dynamic, 1>,
+               Eigen::Matrix<double, Eigen::Dynamic, 1>,
+               Eigen::Matrix<double, Eigen::Dynamic, 1>>>;
+
+template <size_t sizeof_args>
+using base_set_t = std::conditional_t<
+    sizeof_args == 1, unary_base_set,
+    std::conditional_t<
+        sizeof_args == 2, binary_base_set,
+        std::conditional_t<sizeof_args == 3, ternary_base_sets,
+                           std::conditional_t<sizeof_args == 4, quad_base_set,
+                                              quad_base_set>>>>;
+
+template <typename F, typename... EigVecs>
+void expect_ad_distribution(F&& f, const EigVecs&... args) {
+  // Need to test all combinations of scalars, Eigen vectors, and std vectors
+  using base_set = base_set_t<sizeof...(args)>;
+  // Would be nice to have version of for_each that just took in a template
+  stan::math::for_each(
+      [&f, &args...](auto&& type_tuple) {
+        stan::math::apply(
+            [&f, &args...](auto&&... types) {
+              stan::test::expect_ad(
+                  ad_tolerances{}, f,
+                  do_cast<std::decay_t<decltype(types)>>(args)...);
+            },
+            type_tuple);
+      },
+      base_set{});
+}
+
+}  // namespace test
+}  // namespace stan
+#endif
diff --git a/test/unit/math/prim/prob/normal_log_test.cpp b/test/unit/math/prim/prob/normal_log_test.cpp
index fc4906fdc7b..e637d1055b1 100644
--- a/test/unit/math/prim/prob/normal_log_test.cpp
+++ b/test/unit/math/prim/prob/normal_log_test.cpp
@@ -12,13 +12,20 @@ TEST(ProbNormal, log_matches_lpdf) {
                   (stan::math::normal_log<true>(y, mu, sigma)));
   EXPECT_FLOAT_EQ((stan::math::normal_lpdf<false>(y, mu, sigma)),
                   (stan::math::normal_log<false>(y, mu, sigma)));
-  EXPECT_FLOAT_EQ(
-      (stan::math::normal_lpdf<true, double, double, double>(y, mu, sigma)),
-      (stan::math::normal_log<true, double, double, double>(y, mu, sigma)));
-  EXPECT_FLOAT_EQ(
-      (stan::math::normal_lpdf<false, double, double, double>(y, mu, sigma)),
-      (stan::math::normal_log<false, double, double, double>(y, mu, sigma)));
-  EXPECT_FLOAT_EQ(
-      (stan::math::normal_lpdf<double, double, double>(y, mu, sigma)),
-      (stan::math::normal_log<double, double, double>(y, mu, sigma)));
+  EXPECT_FLOAT_EQ((stan::math::normal_lpdf<true>(y, mu, sigma)),
+                  (stan::math::normal_log<true>(y, mu, sigma)));
+  EXPECT_FLOAT_EQ((stan::math::normal_lpdf<false>(y, mu, sigma)),
+                  (stan::math::normal_log<false>(y, mu, sigma)));
+  EXPECT_FLOAT_EQ((stan::math::normal_lpdf(y, mu, sigma)),
+                  (stan::math::normal_log(y, mu, sigma)));
+}
+
+TEST(ProbNormal, test_vlpdf) {
+  Eigen::Matrix<double, -1, 1> Y = Eigen::Matrix<double, -1, 1>::Random(5);
+  Eigen::Matrix<double, -1, 1> Mu = Eigen::Matrix<double, -1, 1>::Random(5);
+  Eigen::Matrix<double, -1, 1> Sigma
+      = stan::math::abs(Eigen::Matrix<double, -1, 1>::Random(5));
+  Eigen::Matrix<double, -1, 1> A
+      = stan::math::normal_lpdf<false, stan::math::ProbReturnType::Vector>(
+          Y, Mu, Sigma);
 }
diff --git a/test/unit/math/rev/prob/normal_log_test.cpp b/test/unit/math/rev/prob/normal_log_test.cpp
index da5d17a6f3b..9bad8dbe184 100644
--- a/test/unit/math/rev/prob/normal_log_test.cpp
+++ b/test/unit/math/rev/prob/normal_log_test.cpp
@@ -1,7 +1,7 @@
 #include <stan/math/rev.hpp>
 #include <gtest/gtest.h>
 
-TEST(ProbDistributionsNormal, intVsDouble) {
+TEST(ProbDistributionsNormal, intVsDoublerev) {
   using stan::math::var;
   for (double thetaval = -5.0; thetaval < 6.0; thetaval += 0.5) {
     var theta(thetaval);
@@ -21,3 +21,14 @@ TEST(ProbDistributionsNormal, intVsDouble) {
     EXPECT_FLOAT_EQ(lp1adj, lp2adj);
   }
 }
+
+TEST(ProbNormal, test_vlpdf_rev) {
+  using stan::math::var;
+  Eigen::Matrix<var, -1, 1> Y = Eigen::Matrix<double, -1, 1>::Random(5);
+  Eigen::Matrix<var, -1, 1> Mu = Eigen::Matrix<double, -1, 1>::Random(5);
+  Eigen::Matrix<var, -1, 1> Sigma
+      = stan::math::abs(Eigen::Matrix<double, -1, 1>::Random(5));
+  Eigen::Matrix<var, -1, 1> A
+      = stan::math::normal_lpdf<false, stan::math::ProbReturnType::Vector>(
+          Y, Mu, Sigma);
+}