From 38edc834a2f74fdf86c8d6942ff4cde4953715f1 Mon Sep 17 00:00:00 2001
From: Brian Ward <bward@flatironinstitute.org>
Date: Thu, 3 Sep 2026 16:09:03 -0400
Subject: [PATCH 1/3] Add size-0 check to vector rng rig

---
 .../math/prim/prob/vector_rng_test_helper.hpp | 22 +++++++++++++++++++
 1 file changed, 22 insertions(+)

diff --git a/test/unit/math/prim/prob/vector_rng_test_helper.hpp b/test/unit/math/prim/prob/vector_rng_test_helper.hpp
index 717ae8a7ec0..27d5d7ae768 100644
--- a/test/unit/math/prim/prob/vector_rng_test_helper.hpp
+++ b/test/unit/math/prim/prob/vector_rng_test_helper.hpp
@@ -421,6 +421,17 @@ struct check_quantiles {
 
       assert_matches_quantiles(samples_to_test, quantiles, 1e-6);
     }
+
+    // check size-0 behavior
+    if constexpr (stan::is_vector<T_param1>::value
+                  || stan::is_vector<T_param2>::value
+                  || stan::is_vector<T_param3>::value) {
+      resize_if_vector(p1, 0);
+      resize_if_vector(p2, 0);
+      resize_if_vector(p3, 0);
+      auto result = rig.generate_samples(p1, p2, p3, rng);
+      EXPECT_EQ(result.size(), 0);
+    }
   }
 };
 
@@ -626,6 +637,17 @@ struct check_counts {
         assert_chi_squared(counts_trimmed, epmf_trimmed, 1e-6);
       }
     }
+
+    // check size-0 behavior
+    if constexpr (stan::is_vector<T_param1>::value
+                  || stan::is_vector<T_param2>::value
+                  || stan::is_vector<T_param3>::value) {
+      resize_if_vector(p1, 0);
+      resize_if_vector(p2, 0);
+      resize_if_vector(p3, 0);
+      auto result = rig.generate_samples(p1, p2, p3, rng);
+      EXPECT_EQ(result.size(), 0);
+    }
   }
 };
 

From 69a0cb4c0b3cb950bc9bcde3368972d6bdac8ade Mon Sep 17 00:00:00 2001
From: Brian Ward <bward@flatironinstitute.org>
Date: Thu, 3 Sep 2026 16:19:10 -0400
Subject: [PATCH 2/3] Fix normal_rng on size 0

---
 stan/math/prim/prob/normal_rng.hpp | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/stan/math/prim/prob/normal_rng.hpp b/stan/math/prim/prob/normal_rng.hpp
index 73c328eb9e4..64993ce5d4a 100644
--- a/stan/math/prim/prob/normal_rng.hpp
+++ b/stan/math/prim/prob/normal_rng.hpp
@@ -44,6 +44,10 @@ inline typename VectorBuilder<true, double, T_loc, T_scale>::type normal_rng(
   check_finite(function, "Location parameter", mu_ref);
   check_positive_finite(function, "Scale parameter", sigma_ref);
 
+  if (size_zero(mu_ref, sigma_ref)) {
+    return {};
+  }
+
   scalar_seq_view<T_mu_ref> mu_vec(mu_ref);
   scalar_seq_view<T_sigma_ref> sigma_vec(sigma_ref);
   size_t N = max_size(mu, sigma);

From 6db3b7395bc9519a0bdc5a60c934eaea803fd149 Mon Sep 17 00:00:00 2001
From: Brian Ward <bward@flatironinstitute.org>
Date: Thu, 3 Sep 2026 16:59:56 -0400
Subject: [PATCH 3/3] Add size_zero check to all multi-arg vectorized rngs

---
 stan/math/prim/prob/beta_binomial_rng.hpp           | 5 +++++
 stan/math/prim/prob/beta_proportion_rng.hpp         | 6 ++++++
 stan/math/prim/prob/beta_rng.hpp                    | 5 +++++
 stan/math/prim/prob/binomial_rng.hpp                | 5 +++++
 stan/math/prim/prob/cauchy_rng.hpp                  | 5 +++++
 stan/math/prim/prob/discrete_range_rng.hpp          | 4 ++++
 stan/math/prim/prob/double_exponential_rng.hpp      | 5 +++++
 stan/math/prim/prob/exp_mod_normal_rng.hpp          | 5 ++++-
 stan/math/prim/prob/frechet_rng.hpp                 | 5 +++++
 stan/math/prim/prob/gamma_rng.hpp                   | 5 +++++
 stan/math/prim/prob/gumbel_rng.hpp                  | 4 ++++
 stan/math/prim/prob/inv_gamma_rng.hpp               | 5 +++++
 stan/math/prim/prob/logistic_rng.hpp                | 5 +++++
 stan/math/prim/prob/loglogistic_rng.hpp             | 4 ++++
 stan/math/prim/prob/lognormal_rng.hpp               | 4 ++++
 stan/math/prim/prob/neg_binomial_2_log_rng.hpp      | 4 ++++
 stan/math/prim/prob/neg_binomial_2_rng.hpp          | 4 ++++
 stan/math/prim/prob/neg_binomial_rng.hpp            | 6 ++++++
 stan/math/prim/prob/normal_rng.hpp                  | 9 +++++----
 stan/math/prim/prob/pareto_rng.hpp                  | 4 ++++
 stan/math/prim/prob/pareto_type_2_rng.hpp           | 4 ++++
 stan/math/prim/prob/scaled_inv_chi_square_rng.hpp   | 4 ++++
 stan/math/prim/prob/skew_double_exponential_rng.hpp | 4 ++++
 stan/math/prim/prob/skew_normal_rng.hpp             | 4 ++++
 stan/math/prim/prob/student_t_rng.hpp               | 4 ++++
 stan/math/prim/prob/uniform_rng.hpp                 | 4 ++++
 stan/math/prim/prob/von_mises_rng.hpp               | 4 ++++
 stan/math/prim/prob/weibull_rng.hpp                 | 4 ++++
 28 files changed, 126 insertions(+), 5 deletions(-)

diff --git a/stan/math/prim/prob/beta_binomial_rng.hpp b/stan/math/prim/prob/beta_binomial_rng.hpp
index a96ecfd5c3c..8deba0cd65a 100644
--- a/stan/math/prim/prob/beta_binomial_rng.hpp
+++ b/stan/math/prim/prob/beta_binomial_rng.hpp
@@ -4,6 +4,7 @@
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/prob/binomial_rng.hpp>
 #include <stan/math/prim/prob/beta_rng.hpp>
+#include <stan/math/prim/fun/size_zero.hpp>
 
 namespace stan {
 namespace math {
@@ -39,6 +40,10 @@ beta_binomial_rng(const T_N &N, const T_shape1 &alpha, const T_shape2 &beta,
   check_consistent_sizes(function, "First prior sample size parameter", alpha,
                          "Second prior sample size parameter", beta);
 
+  if (size_zero(N, alpha, beta)) {
+    return {};
+  }
+
   T_N_ref N_ref = N;
   T_alpha_ref alpha_ref = alpha;
   T_beta_ref beta_ref = beta;
diff --git a/stan/math/prim/prob/beta_proportion_rng.hpp b/stan/math/prim/prob/beta_proportion_rng.hpp
index 61bb114f00f..bd8504152ef 100644
--- a/stan/math/prim/prob/beta_proportion_rng.hpp
+++ b/stan/math/prim/prob/beta_proportion_rng.hpp
@@ -4,6 +4,7 @@
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
+#include <stan/math/prim/fun/size_zero.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <stan/math/prim/fun/to_ref.hpp>
 #include <boost/random/variate_generator.hpp>
@@ -39,6 +40,11 @@ beta_proportion_rng(const T_loc &mu, const T_prec &kappa, RNG &rng) {
   static constexpr const char *function = "beta_proportion_rng";
   check_consistent_sizes(function, "Location parameter", mu,
                          "Precision parameter", kappa);
+
+  if (size_zero(mu, kappa)) {
+    return {};
+  }
+
   T_mu_ref mu_ref = mu;
   T_kappa_ref kappa_ref = kappa;
   check_positive(function, "Location parameter", mu_ref);
diff --git a/stan/math/prim/prob/beta_rng.hpp b/stan/math/prim/prob/beta_rng.hpp
index 9fdc839bf4e..bde99ddbb1f 100644
--- a/stan/math/prim/prob/beta_rng.hpp
+++ b/stan/math/prim/prob/beta_rng.hpp
@@ -5,6 +5,7 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/log_sum_exp.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
+#include <stan/math/prim/fun/size_zero.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <boost/random/gamma_distribution.hpp>
 #include <boost/random/uniform_real_distribution.hpp>
@@ -43,6 +44,10 @@ inline typename VectorBuilder<true, double, T_shape1, T_shape2>::type beta_rng(
   static constexpr const char *function = "beta_rng";
   check_consistent_sizes(function, "First shape parameter", alpha,
                          "Second shape Parameter", beta);
+  if (size_zero(alpha, beta)) {
+    return {};
+  }
+
   T_alpha_ref alpha_ref = alpha;
   T_beta_ref beta_ref = beta;
   check_positive_finite(function, "First shape parameter", alpha_ref);
diff --git a/stan/math/prim/prob/binomial_rng.hpp b/stan/math/prim/prob/binomial_rng.hpp
index 3dab00a3649..8ebbf2922c2 100644
--- a/stan/math/prim/prob/binomial_rng.hpp
+++ b/stan/math/prim/prob/binomial_rng.hpp
@@ -4,6 +4,7 @@
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
+#include <stan/math/prim/fun/size_zero.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <boost/random/binomial_distribution.hpp>
 #include <boost/random/variate_generator.hpp>
@@ -38,6 +39,10 @@ inline typename VectorBuilder<true, int, T_N, T_theta>::type binomial_rng(
   static constexpr const char* function = "binomial_rng";
   check_consistent_sizes(function, "Population size parameter", N,
                          "Probability Parameter", theta);
+  if (size_zero(N, theta)) {
+    return {};
+  }
+
   T_N_ref N_ref = N;
   T_theta_ref theta_ref = theta;
   check_nonnegative(function, "Population size parameter", N_ref);
diff --git a/stan/math/prim/prob/cauchy_rng.hpp b/stan/math/prim/prob/cauchy_rng.hpp
index 21ddc51d505..addc8d36e96 100644
--- a/stan/math/prim/prob/cauchy_rng.hpp
+++ b/stan/math/prim/prob/cauchy_rng.hpp
@@ -4,6 +4,7 @@
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
+#include <stan/math/prim/fun/size_zero.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <boost/random/cauchy_distribution.hpp>
 #include <boost/random/variate_generator.hpp>
@@ -39,6 +40,10 @@ inline typename VectorBuilder<true, double, T_loc, T_scale>::type cauchy_rng(
   using T_sigma_ref = ref_type_t<T_scale>;
   check_consistent_sizes(function, "Location parameter", mu, "Scale Parameter",
                          sigma);
+  if (size_zero(mu, sigma)) {
+    return {};
+  }
+
   T_mu_ref mu_ref = mu;
   T_sigma_ref sigma_ref = sigma;
   check_finite(function, "Location parameter", mu_ref);
diff --git a/stan/math/prim/prob/discrete_range_rng.hpp b/stan/math/prim/prob/discrete_range_rng.hpp
index 2dacf745339..c74f911008e 100644
--- a/stan/math/prim/prob/discrete_range_rng.hpp
+++ b/stan/math/prim/prob/discrete_range_rng.hpp
@@ -4,6 +4,7 @@
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
+#include <stan/math/prim/fun/size_zero.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <boost/random/uniform_int_distribution.hpp>
 #include <boost/random/variate_generator.hpp>
@@ -40,6 +41,9 @@ discrete_range_rng(const T_lower& lower, const T_upper& upper, RNG& rng) {
   check_consistent_sizes(function, "Lower bound parameter", lower,
                          "Upper bound parameter", upper);
   check_greater_or_equal(function, "Upper bound parameter", upper, lower);
+  if (size_zero(upper, lower)) {
+    return {};
+  }
 
   scalar_seq_view<T_lower> lower_vec(lower);
   scalar_seq_view<T_upper> upper_vec(upper);
diff --git a/stan/math/prim/prob/double_exponential_rng.hpp b/stan/math/prim/prob/double_exponential_rng.hpp
index a5ee59a89ac..6913d43a939 100644
--- a/stan/math/prim/prob/double_exponential_rng.hpp
+++ b/stan/math/prim/prob/double_exponential_rng.hpp
@@ -5,6 +5,7 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/log1m.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
+#include <stan/math/prim/fun/size_zero.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <boost/random/uniform_real_distribution.hpp>
 #include <boost/random/variate_generator.hpp>
@@ -41,6 +42,10 @@ double_exponential_rng(const T_loc& mu, const T_scale& sigma, RNG& rng) {
   static constexpr const char* function = "double_exponential_rng";
   check_consistent_sizes(function, "Location parameter", mu, "Scale Parameter",
                          sigma);
+  if (size_zero(mu, sigma)) {
+    return {};
+  }
+
   T_mu_ref mu_ref = mu;
   T_sigma_ref sigma_ref = sigma;
   check_finite(function, "Location parameter", mu_ref);
diff --git a/stan/math/prim/prob/exp_mod_normal_rng.hpp b/stan/math/prim/prob/exp_mod_normal_rng.hpp
index 84df8fd1733..a3d7307b8f6 100644
--- a/stan/math/prim/prob/exp_mod_normal_rng.hpp
+++ b/stan/math/prim/prob/exp_mod_normal_rng.hpp
@@ -5,6 +5,7 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/prob/exponential_rng.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
+#include <stan/math/prim/fun/size_zero.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <stan/math/prim/prob/normal_rng.hpp>
 #include <boost/random/normal_distribution.hpp>
@@ -45,7 +46,9 @@ exp_mod_normal_rng(const T_loc& mu, const T_scale& sigma,
   using T_lambda_ref = ref_type_t<T_inv_scale>;
   check_consistent_sizes(function, "Location parameter", mu, "Scale Parameter",
                          sigma, "Inv_scale Parameter", lambda);
-
+  // if (size_zero(mu, sigma, lambda)) {
+  //   return {};
+  // }
   T_mu_ref mu_ref = mu;
   T_sigma_ref sigma_ref = sigma;
   T_lambda_ref lambda_ref = lambda;
diff --git a/stan/math/prim/prob/frechet_rng.hpp b/stan/math/prim/prob/frechet_rng.hpp
index 0c15df9bc35..821f900c409 100644
--- a/stan/math/prim/prob/frechet_rng.hpp
+++ b/stan/math/prim/prob/frechet_rng.hpp
@@ -4,6 +4,7 @@
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
+#include <stan/math/prim/fun/size_zero.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <boost/random/weibull_distribution.hpp>
 #include <boost/random/variate_generator.hpp>
@@ -38,6 +39,10 @@ inline typename VectorBuilder<true, double, T_shape, T_scale>::type frechet_rng(
   static constexpr const char* function = "frechet_rng";
   check_consistent_sizes(function, "Shape parameter", alpha, "Scale Parameter",
                          sigma);
+  if (size_zero(alpha, sigma)) {
+    return {};
+  }
+
   T_alpha_ref alpha_ref = alpha;
   T_sigma_ref sigma_ref = sigma;
   check_positive_finite(function, "Shape parameter", alpha_ref);
diff --git a/stan/math/prim/prob/gamma_rng.hpp b/stan/math/prim/prob/gamma_rng.hpp
index c609e34214e..0d78382321d 100644
--- a/stan/math/prim/prob/gamma_rng.hpp
+++ b/stan/math/prim/prob/gamma_rng.hpp
@@ -5,6 +5,7 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
+#include <stan/math/prim/fun/size_zero.hpp>
 #include <boost/random/gamma_distribution.hpp>
 #include <boost/random/variate_generator.hpp>
 
@@ -39,6 +40,10 @@ inline typename VectorBuilder<true, double, T_shape, T_inv>::type gamma_rng(
   static constexpr const char* function = "gamma_rng";
   check_consistent_sizes(function, "Shape parameter", alpha,
                          "Inverse scale Parameter", beta);
+  if (size_zero(alpha, beta)) {
+    return {};
+  }
+
   T_alpha_ref alpha_ref = alpha;
   T_beta_ref beta_ref = beta;
   check_positive_finite(function, "Shape parameter", alpha_ref);
diff --git a/stan/math/prim/prob/gumbel_rng.hpp b/stan/math/prim/prob/gumbel_rng.hpp
index b6ece64a4ec..79cbd891a02 100644
--- a/stan/math/prim/prob/gumbel_rng.hpp
+++ b/stan/math/prim/prob/gumbel_rng.hpp
@@ -5,6 +5,7 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
+#include <stan/math/prim/fun/size_zero.hpp>
 #include <boost/random/uniform_01.hpp>
 #include <boost/random/variate_generator.hpp>
 #include <cmath>
@@ -40,6 +41,9 @@ inline typename VectorBuilder<true, double, T_loc, T_scale>::type gumbel_rng(
   static constexpr const char* function = "gumbel_rng";
   check_consistent_sizes(function, "Location parameter", mu, "Scale Parameter",
                          beta);
+  if (size_zero(mu, beta)) {
+    return {};
+  }
   T_mu_ref mu_ref = mu;
   T_beta_ref beta_ref = beta;
   check_finite(function, "Location parameter", mu_ref);
diff --git a/stan/math/prim/prob/inv_gamma_rng.hpp b/stan/math/prim/prob/inv_gamma_rng.hpp
index cd7baeb7bc1..696b6211cb7 100644
--- a/stan/math/prim/prob/inv_gamma_rng.hpp
+++ b/stan/math/prim/prob/inv_gamma_rng.hpp
@@ -5,6 +5,7 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
+#include <stan/math/prim/fun/size_zero.hpp>
 #include <boost/random/gamma_distribution.hpp>
 #include <boost/random/variate_generator.hpp>
 
@@ -39,6 +40,10 @@ inv_gamma_rng(const T_shape& alpha, const T_scale& beta, RNG& rng) {
   static constexpr const char* function = "inv_gamma_rng";
   check_consistent_sizes(function, "Shape parameter", alpha, "Scale Parameter",
                          beta);
+  if (size_zero(alpha, beta)) {
+    return {};
+  }
+
   T_alpha_ref alpha_ref = alpha;
   T_beta_ref beta_ref = beta;
   check_positive_finite(function, "Shape parameter", alpha_ref);
diff --git a/stan/math/prim/prob/logistic_rng.hpp b/stan/math/prim/prob/logistic_rng.hpp
index 9b7e94f1484..f37f3ce9a93 100644
--- a/stan/math/prim/prob/logistic_rng.hpp
+++ b/stan/math/prim/prob/logistic_rng.hpp
@@ -4,6 +4,7 @@
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
+#include <stan/math/prim/fun/size_zero.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <boost/random/exponential_distribution.hpp>
 #include <boost/random/variate_generator.hpp>
@@ -40,6 +41,10 @@ inline typename VectorBuilder<true, double, T_loc, T_scale>::type logistic_rng(
   static constexpr const char* function = "logistic_rng";
   check_consistent_sizes(function, "Location parameter", mu, "Scale Parameter",
                          sigma);
+  if (size_zero(mu, sigma)) {
+    return {};
+  }
+
   T_mu_ref mu_ref = mu;
   T_sigma_ref sigma_ref = sigma;
   check_finite(function, "Location parameter", mu_ref);
diff --git a/stan/math/prim/prob/loglogistic_rng.hpp b/stan/math/prim/prob/loglogistic_rng.hpp
index fa92d93826e..c43aa90c036 100644
--- a/stan/math/prim/prob/loglogistic_rng.hpp
+++ b/stan/math/prim/prob/loglogistic_rng.hpp
@@ -5,6 +5,7 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
+#include <stan/math/prim/fun/size_zero.hpp>
 #include <boost/random/uniform_01.hpp>
 #include <boost/random/variate_generator.hpp>
 
@@ -40,6 +41,9 @@ loglogistic_rng(const T_scale& alpha, const T_shape& beta, RNG& rng) {
   static constexpr const char* function = "loglogistic_rng";
   check_consistent_sizes(function, "Scale parameter", alpha, "Shape Parameter",
                          beta);
+  if (size_zero(alpha, beta)) {
+    return {};
+  }
   T_alpha_ref alpha_ref = alpha;
   T_beta_ref beta_ref = beta;
   check_positive_finite(function, "Scale parameter", alpha_ref);
diff --git a/stan/math/prim/prob/lognormal_rng.hpp b/stan/math/prim/prob/lognormal_rng.hpp
index 64adc0516c5..e09f86b974d 100644
--- a/stan/math/prim/prob/lognormal_rng.hpp
+++ b/stan/math/prim/prob/lognormal_rng.hpp
@@ -4,6 +4,7 @@
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
+#include <stan/math/prim/fun/size_zero.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <boost/random/lognormal_distribution.hpp>
 #include <boost/random/variate_generator.hpp>
@@ -39,6 +40,9 @@ inline typename VectorBuilder<true, double, T_loc, T_scale>::type lognormal_rng(
   static constexpr const char* function = "lognormal_rng";
   check_consistent_sizes(function, "Location parameter", mu, "Scale Parameter",
                          sigma);
+  if (size_zero(mu, sigma)) {
+    return {};
+  }
   T_mu_ref mu_ref = mu;
   T_sigma_ref sigma_ref = sigma;
   check_finite(function, "Location parameter", mu_ref);
diff --git a/stan/math/prim/prob/neg_binomial_2_log_rng.hpp b/stan/math/prim/prob/neg_binomial_2_log_rng.hpp
index aefce01d132..d2a66517e28 100644
--- a/stan/math/prim/prob/neg_binomial_2_log_rng.hpp
+++ b/stan/math/prim/prob/neg_binomial_2_log_rng.hpp
@@ -5,6 +5,7 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/constants.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
+#include <stan/math/prim/fun/size_zero.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <boost/random/gamma_distribution.hpp>
 #include <boost/random/poisson_distribution.hpp>
@@ -43,6 +44,9 @@ neg_binomial_2_log_rng(const T_loc& eta, const T_inv& phi, RNG& rng) {
   static constexpr const char* function = "neg_binomial_2_log_rng";
   check_consistent_sizes(function, "Log-location parameter", eta,
                          "Inverse dispersion parameter", phi);
+  if (size_zero(eta, phi)) {
+    return {};
+  }
   T_eta_ref eta_ref = eta;
   T_phi_ref phi_ref = phi;
   check_finite(function, "Log-location parameter", eta_ref);
diff --git a/stan/math/prim/prob/neg_binomial_2_rng.hpp b/stan/math/prim/prob/neg_binomial_2_rng.hpp
index 67f6b69687d..99492f75cbf 100644
--- a/stan/math/prim/prob/neg_binomial_2_rng.hpp
+++ b/stan/math/prim/prob/neg_binomial_2_rng.hpp
@@ -5,6 +5,7 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/constants.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
+#include <stan/math/prim/fun/size_zero.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <boost/random/gamma_distribution.hpp>
 #include <boost/random/poisson_distribution.hpp>
@@ -42,6 +43,9 @@ neg_binomial_2_rng(const T_loc& mu, const T_prec& phi, RNG& rng) {
   static constexpr const char* function = "neg_binomial_2_rng";
   check_consistent_sizes(function, "Location parameter", mu,
                          "Precision parameter", phi);
+  if (size_zero(mu, phi)) {
+    return {};
+  }
   T_mu_ref mu_ref = mu;
   T_phi_ref phi_ref = phi;
   check_positive_finite(function, "Location parameter", mu_ref);
diff --git a/stan/math/prim/prob/neg_binomial_rng.hpp b/stan/math/prim/prob/neg_binomial_rng.hpp
index 0e01f6693c8..c8537051d2c 100644
--- a/stan/math/prim/prob/neg_binomial_rng.hpp
+++ b/stan/math/prim/prob/neg_binomial_rng.hpp
@@ -5,6 +5,7 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/constants.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
+#include <stan/math/prim/fun/size_zero.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <boost/random/gamma_distribution.hpp>
 #include <boost/random/poisson_distribution.hpp>
@@ -44,6 +45,11 @@ inline typename VectorBuilder<true, int, T_shape, T_inv>::type neg_binomial_rng(
                          "Inverse scale Parameter", beta);
   T_alpha_ref alpha_ref = alpha;
   T_beta_ref beta_ref = beta;
+
+  if (size_zero(alpha, beta)) {
+    return {};
+  }
+
   check_positive_finite(function, "Shape parameter", alpha_ref);
   check_positive_finite(function, "Inverse scale parameter", beta_ref);
 
diff --git a/stan/math/prim/prob/normal_rng.hpp b/stan/math/prim/prob/normal_rng.hpp
index 64993ce5d4a..cf433588201 100644
--- a/stan/math/prim/prob/normal_rng.hpp
+++ b/stan/math/prim/prob/normal_rng.hpp
@@ -5,6 +5,7 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
+#include <stan/math/prim/fun/size_zero.hpp>
 #include <boost/random/normal_distribution.hpp>
 #include <boost/random/variate_generator.hpp>
 
@@ -39,15 +40,15 @@ inline typename VectorBuilder<true, double, T_loc, T_scale>::type normal_rng(
   static constexpr const char* function = "normal_rng";
   check_consistent_sizes(function, "Location parameter", mu, "Scale Parameter",
                          sigma);
+  if (size_zero(mu, sigma)) {
+    return {};
+  }
+
   T_mu_ref mu_ref = mu;
   T_sigma_ref sigma_ref = sigma;
   check_finite(function, "Location parameter", mu_ref);
   check_positive_finite(function, "Scale parameter", sigma_ref);
 
-  if (size_zero(mu_ref, sigma_ref)) {
-    return {};
-  }
-
   scalar_seq_view<T_mu_ref> mu_vec(mu_ref);
   scalar_seq_view<T_sigma_ref> sigma_vec(sigma_ref);
   size_t N = max_size(mu, sigma);
diff --git a/stan/math/prim/prob/pareto_rng.hpp b/stan/math/prim/prob/pareto_rng.hpp
index aac348d7c0e..a3b82731f54 100644
--- a/stan/math/prim/prob/pareto_rng.hpp
+++ b/stan/math/prim/prob/pareto_rng.hpp
@@ -4,6 +4,7 @@
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
+#include <stan/math/prim/fun/size_zero.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <stan/math/prim/fun/to_ref.hpp>
 #include <boost/random/exponential_distribution.hpp>
@@ -39,6 +40,9 @@ inline typename VectorBuilder<true, double, T_shape, T_scale>::type pareto_rng(
   static constexpr const char* function = "pareto_rng";
   check_consistent_sizes(function, "Scale Parameter", y_min, "Shape parameter",
                          alpha);
+  if (size_zero(y_min, alpha)) {
+    return {};
+  }
   const auto& y_min_ref = to_ref(y_min);
   const auto& alpha_ref = to_ref(alpha);
   check_positive_finite(function, "Scale parameter", y_min_ref);
diff --git a/stan/math/prim/prob/pareto_type_2_rng.hpp b/stan/math/prim/prob/pareto_type_2_rng.hpp
index e2bcd5f0db8..e83886039ca 100644
--- a/stan/math/prim/prob/pareto_type_2_rng.hpp
+++ b/stan/math/prim/prob/pareto_type_2_rng.hpp
@@ -4,6 +4,7 @@
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
+#include <stan/math/prim/fun/size_zero.hpp>
 #include <stan/math/prim/fun/to_ref.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <stan/math/prim/prob/exponential_rng.hpp>
@@ -46,6 +47,9 @@ pareto_type_2_rng(const T_loc& mu, const T_scale& lambda, const T_shape& alpha,
   static constexpr const char* function = "pareto_type_2_rng";
   check_consistent_sizes(function, "Location parameter", mu, "Scale Parameter",
                          lambda, "Shape Parameter", alpha);
+  if (size_zero(mu, lambda, alpha)) {
+    return {};
+  }
   const auto& mu_ref = to_ref(mu);
   const auto& lambda_ref = to_ref(lambda);
   const auto& alpha_ref = to_ref(alpha);
diff --git a/stan/math/prim/prob/scaled_inv_chi_square_rng.hpp b/stan/math/prim/prob/scaled_inv_chi_square_rng.hpp
index 819255a94a1..20a8ac3e12c 100644
--- a/stan/math/prim/prob/scaled_inv_chi_square_rng.hpp
+++ b/stan/math/prim/prob/scaled_inv_chi_square_rng.hpp
@@ -4,6 +4,7 @@
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
+#include <stan/math/prim/fun/size_zero.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <stan/math/prim/fun/to_ref.hpp>
 #include <boost/random/chi_squared_distribution.hpp>
@@ -39,6 +40,9 @@ scaled_inv_chi_square_rng(const T_deg& nu, const T_scale& s, RNG& rng) {
   static constexpr const char* function = "scaled_inv_chi_square_rng";
   check_consistent_sizes(function, "Location parameter", nu, "Scale Parameter",
                          s);
+  if (size_zero(nu, s)) {
+    return {};
+  }
   const auto& nu_ref = to_ref(nu);
   const auto& s_ref = to_ref(s);
   check_positive_finite(function, "Degrees of freedom parameter", nu_ref);
diff --git a/stan/math/prim/prob/skew_double_exponential_rng.hpp b/stan/math/prim/prob/skew_double_exponential_rng.hpp
index 7b238d71189..8df1529802c 100644
--- a/stan/math/prim/prob/skew_double_exponential_rng.hpp
+++ b/stan/math/prim/prob/skew_double_exponential_rng.hpp
@@ -5,6 +5,7 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/log1m.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
+#include <stan/math/prim/fun/size_zero.hpp>
 #include <boost/random/uniform_real_distribution.hpp>
 #include <boost/random/variate_generator.hpp>
 #include <cmath>
@@ -45,6 +46,9 @@ skew_double_exponential_rng(const T_loc& mu, const T_scale& sigma,
   static constexpr const char* function = "skew_double_exponential_rng";
   check_consistent_sizes(function, "Location parameter", mu, "Scale Parameter",
                          sigma, "Skewness Parameter", tau);
+  if (size_zero(mu, sigma, tau)) {
+    return {};
+  }
   T_mu_ref mu_ref = mu;
   T_sigma_ref sigma_ref = sigma;
   T_tau_ref tau_ref = tau;
diff --git a/stan/math/prim/prob/skew_normal_rng.hpp b/stan/math/prim/prob/skew_normal_rng.hpp
index d4cc8590de7..b64f7d4e86a 100644
--- a/stan/math/prim/prob/skew_normal_rng.hpp
+++ b/stan/math/prim/prob/skew_normal_rng.hpp
@@ -4,6 +4,7 @@
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
+#include <stan/math/prim/fun/size_zero.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <stan/math/prim/fun/to_ref.hpp>
 #include <stan/math/prim/prob/normal_rng.hpp>
@@ -44,6 +45,9 @@ skew_normal_rng(const T_loc& mu, const T_scale& sigma, const T_shape& alpha,
   static constexpr const char* function = "skew_normal_rng";
   check_consistent_sizes(function, "Location parameter", mu, "Scale Parameter",
                          sigma, "Shape Parameter", alpha);
+  if (size_zero(mu, sigma, alpha)) {
+    return {};
+  }
   const auto& mu_ref = to_ref(mu);
   const auto& sigma_ref = to_ref(sigma);
   const auto& alpha_ref = to_ref(alpha);
diff --git a/stan/math/prim/prob/student_t_rng.hpp b/stan/math/prim/prob/student_t_rng.hpp
index ce7d271ef30..1c080113da7 100644
--- a/stan/math/prim/prob/student_t_rng.hpp
+++ b/stan/math/prim/prob/student_t_rng.hpp
@@ -4,6 +4,7 @@
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
+#include <stan/math/prim/fun/size_zero.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <boost/random/student_t_distribution.hpp>
 #include <boost/random/variate_generator.hpp>
@@ -45,6 +46,9 @@ student_t_rng(const T_deg& nu, const T_loc& mu, const T_scale& sigma,
   static constexpr const char* function = "student_t_rng";
   check_consistent_sizes(function, "Degrees of freedom parameter", nu,
                          "Location parameter", mu, "Scale Parameter", sigma);
+  if (size_zero(nu, mu, sigma)) {
+    return {};
+  }
   T_nu_ref nu_ref = nu;
   T_mu_ref mu_ref = mu;
   T_sigma_ref sigma_ref = sigma;
diff --git a/stan/math/prim/prob/uniform_rng.hpp b/stan/math/prim/prob/uniform_rng.hpp
index ad18cf1add3..f2011c5d189 100644
--- a/stan/math/prim/prob/uniform_rng.hpp
+++ b/stan/math/prim/prob/uniform_rng.hpp
@@ -5,6 +5,7 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/constants.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
+#include <stan/math/prim/fun/size_zero.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <stan/math/prim/fun/value_of.hpp>
 #include <boost/random/uniform_real_distribution.hpp>
@@ -41,6 +42,9 @@ inline typename VectorBuilder<true, double, T_alpha, T_beta>::type uniform_rng(
   static constexpr const char* function = "uniform_rng";
   check_consistent_sizes(function, "Lower bound parameter", alpha,
                          "Upper bound parameter", beta);
+  if (size_zero(alpha, beta)) {
+    return {};
+  }
   T_alpha_ref alpha_ref = alpha;
   T_beta_ref beta_ref = beta;
   check_finite(function, "Lower bound parameter", alpha_ref);
diff --git a/stan/math/prim/prob/von_mises_rng.hpp b/stan/math/prim/prob/von_mises_rng.hpp
index c68a6227f74..2c6eb81a8f6 100644
--- a/stan/math/prim/prob/von_mises_rng.hpp
+++ b/stan/math/prim/prob/von_mises_rng.hpp
@@ -5,6 +5,7 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/constants.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
+#include <stan/math/prim/fun/size_zero.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <boost/random/uniform_real_distribution.hpp>
 #include <boost/random/variate_generator.hpp>
@@ -52,6 +53,9 @@ inline typename VectorBuilder<true, double, T_loc, T_conc>::type von_mises_rng(
   static constexpr const char* function = "von_mises_rng";
   check_consistent_sizes(function, "Location parameter", mu, "Scale parameter",
                          kappa);
+  if (size_zero(mu, kappa)) {
+    return {};
+  }
   T_mu_ref mu_ref = mu;
   T_kappa_ref kappa_ref = kappa;
 
diff --git a/stan/math/prim/prob/weibull_rng.hpp b/stan/math/prim/prob/weibull_rng.hpp
index 0626add9fd8..ce80c2000d9 100644
--- a/stan/math/prim/prob/weibull_rng.hpp
+++ b/stan/math/prim/prob/weibull_rng.hpp
@@ -4,6 +4,7 @@
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
+#include <stan/math/prim/fun/size_zero.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <boost/random/weibull_distribution.hpp>
 #include <boost/random/variate_generator.hpp>
@@ -39,6 +40,9 @@ inline typename VectorBuilder<true, double, T_shape, T_scale>::type weibull_rng(
   static constexpr const char* function = "weibull_rng";
   check_consistent_sizes(function, "Shape parameter", alpha, "Scale Parameter",
                          sigma);
+  if (size_zero(alpha, sigma)) {
+    return {};
+  }
   T_alpha_ref alpha_ref = alpha;
   T_sigma_ref sigma_ref = sigma;
   check_positive_finite(function, "Shape parameter", alpha_ref);