diff --git a/stan/math/prim/prob.hpp b/stan/math/prim/prob.hpp
index 308b084a165..33c39a183b8 100644
--- a/stan/math/prim/prob.hpp
+++ b/stan/math/prim/prob.hpp
@@ -113,6 +113,11 @@
 #include <stan/math/prim/prob/gamma_rng.hpp>
 #include <stan/math/prim/prob/gaussian_dlm_obs_lpdf.hpp>
 #include <stan/math/prim/prob/gaussian_dlm_obs_rng.hpp>
+#include <stan/math/prim/prob/geometric_cdf.hpp>
+#include <stan/math/prim/prob/geometric_lccdf.hpp>
+#include <stan/math/prim/prob/geometric_lcdf.hpp>
+#include <stan/math/prim/prob/geometric_lpmf.hpp>
+#include <stan/math/prim/prob/geometric_rng.hpp>
 #include <stan/math/prim/prob/gumbel_ccdf_log.hpp>
 #include <stan/math/prim/prob/gumbel_cdf.hpp>
 #include <stan/math/prim/prob/gumbel_cdf_log.hpp>
diff --git a/stan/math/prim/prob/geometric_cdf.hpp b/stan/math/prim/prob/geometric_cdf.hpp
new file mode 100644
index 00000000000..5e0c2474f67
--- /dev/null
+++ b/stan/math/prim/prob/geometric_cdf.hpp
@@ -0,0 +1,90 @@
+#ifndef STAN_MATH_PRIM_PROB_GEOMETRIC_CDF_HPP
+#define STAN_MATH_PRIM_PROB_GEOMETRIC_CDF_HPP
+
+#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/scalar_seq_view.hpp>
+#include <stan/math/prim/fun/size.hpp>
+#include <stan/math/prim/fun/size_zero.hpp>
+#include <stan/math/prim/fun/value_of.hpp>
+#include <vector>
+#include <stan/math/prim/prob/neg_binomial_cdf.hpp>
+#include <stan/math/prim/fun/elt_divide.hpp>
+#include <stan/math/prim/fun/subtract.hpp>
+
+namespace stan {
+namespace math {
+
+/** \ingroup prob_dists
+ * Returns the CDF of the geometric distribution. Given containers of
+ * matching sizes, returns the product of probabilities.
+ *
+ * Delegates to the negative binomial CDF with alpha = 1 and
+ * beta = theta / (1 - theta).
+ *
+ * @tparam T_n type of outcome variable
+ * @tparam T_prob type of success probability parameter
+ *
+ * @param n outcome variable (number of failures before first success)
+ * @param theta success probability parameter
+ * @return probability or product of probabilities
+ * @throw std::domain_error if theta is not in (0, 1]
+ * @throw std::invalid_argument if container sizes mismatch
+ */
+template <typename T_n, typename T_prob>
+inline return_type_t<T_prob> geometric_cdf(const T_n& n, const T_prob& theta) {
+  using T_n_ref = ref_type_t<T_n>;
+  using T_prob_ref = ref_type_t<T_prob>;
+  static constexpr const char* function = "geometric_cdf";
+
+  check_consistent_sizes(function, "Outcome variable", n,
+                         "Success probability parameter", theta);
+  if (size_zero(n, theta)) {
+    return 1.0;
+  }
+
+  T_n_ref n_ref = n;
+  T_prob_ref theta_ref = theta;
+  check_bounded(function, "Success probability parameter", value_of(theta_ref),
+                0.0, 1.0);
+
+  scalar_seq_view<T_n_ref> n_vec(n_ref);
+  for (int i = 0; i < stan::math::size(n); i++) {
+    if (n_vec.val(i) < 0) {
+      return 0.0;
+    }
+  }
+
+  // theta = 1 => CDF is always 1 for n >= 0
+  scalar_seq_view<T_prob_ref> theta_vec(theta_ref);
+  bool all_theta_one = true;
+  for (size_t i = 0; i < stan::math::size(theta); i++) {
+    if (value_of(theta_vec[i]) != 1.0) {
+      all_theta_one = false;
+      break;
+    }
+  }
+  if (all_theta_one) {
+    return 1.0;
+  }
+
+  if constexpr (is_stan_scalar_v<T_prob>) {
+    const auto beta = theta_ref / (1.0 - theta_ref);
+    return neg_binomial_cdf(n_ref, 1, beta);
+  } else if constexpr (is_std_vector_v<T_prob>) {
+    std::vector<value_type_t<T_prob>> beta;
+    beta.reserve(stan::math::size(theta));
+    for (size_t i = 0; i < stan::math::size(theta); i++) {
+      beta.push_back(theta_vec[i] / (1.0 - theta_vec[i]));
+    }
+    return neg_binomial_cdf(n_ref, 1, beta);
+  } else {
+    const auto beta = elt_divide(theta_ref, subtract(1.0, theta_ref));
+    return neg_binomial_cdf(n_ref, 1, beta);
+  }
+}
+
+}  // namespace math
+}  // namespace stan
+#endif
diff --git a/stan/math/prim/prob/geometric_lccdf.hpp b/stan/math/prim/prob/geometric_lccdf.hpp
new file mode 100644
index 00000000000..9927caaf801
--- /dev/null
+++ b/stan/math/prim/prob/geometric_lccdf.hpp
@@ -0,0 +1,87 @@
+#ifndef STAN_MATH_PRIM_PROB_GEOMETRIC_LCCDF_HPP
+#define STAN_MATH_PRIM_PROB_GEOMETRIC_LCCDF_HPP
+
+#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/scalar_seq_view.hpp>
+#include <stan/math/prim/fun/size.hpp>
+#include <stan/math/prim/fun/size_zero.hpp>
+#include <stan/math/prim/fun/value_of.hpp>
+#include <vector>
+#include <stan/math/prim/prob/neg_binomial_lccdf.hpp>
+#include <stan/math/prim/fun/elt_divide.hpp>
+#include <stan/math/prim/fun/subtract.hpp>
+
+namespace stan {
+namespace math {
+
+/** \ingroup prob_dists
+ * Returns the log CCDF of the geometric distribution. Given containers of
+ * matching sizes, returns the log sum of probabilities.
+ *
+ * Delegates to the negative binomial log CCDF with alpha = 1 and
+ * beta = theta / (1 - theta).
+ *
+ * @tparam T_n type of outcome variable
+ * @tparam T_prob type of success probability parameter
+ *
+ * @param n outcome variable (number of failures before first success)
+ * @param theta success probability parameter
+ * @return log complementary probability or log sum
+ * @throw std::domain_error if theta is not in (0, 1]
+ * @throw std::invalid_argument if container sizes mismatch
+ */
+template <typename T_n, typename T_prob>
+inline return_type_t<T_prob> geometric_lccdf(const T_n& n,
+                                             const T_prob& theta) {
+  using T_n_ref = ref_type_t<T_n>;
+  using T_prob_ref = ref_type_t<T_prob>;
+  static constexpr const char* function = "geometric_lccdf";
+
+  check_consistent_sizes(function, "Outcome variable", n,
+                         "Success probability parameter", theta);
+  if (size_zero(n, theta)) {
+    return 0.0;
+  }
+
+  T_n_ref n_ref = n;
+  T_prob_ref theta_ref = theta;
+  check_bounded(function, "Success probability parameter", value_of(theta_ref),
+                0.0, 1.0);
+
+  scalar_seq_view<T_n_ref> n_vec(n_ref);
+  for (int i = 0; i < stan::math::size(n); i++) {
+    if (n_vec.val(i) < 0) {
+      return 0.0;
+    }
+  }
+
+  // theta = 1 => CCDF = 0 for n >= 0, log CCDF = -inf
+  scalar_seq_view<T_prob_ref> theta_vec(theta_ref);
+  const size_t max_sz = max_size(n_ref, theta_ref);
+  for (size_t i = 0; i < max_sz; i++) {
+    if (value_of(theta_vec[i]) == 1.0 && n_vec.val(i) >= 0) {
+      return negative_infinity();
+    }
+  }
+
+  if constexpr (is_stan_scalar_v<T_prob>) {
+    const auto beta = theta_ref / (1.0 - theta_ref);
+    return neg_binomial_lccdf(n_ref, 1, beta);
+  } else if constexpr (is_std_vector_v<T_prob>) {
+    std::vector<value_type_t<T_prob>> beta;
+    beta.reserve(stan::math::size(theta));
+    for (size_t i = 0; i < stan::math::size(theta); i++) {
+      beta.push_back(theta_vec[i] / (1.0 - theta_vec[i]));
+    }
+    return neg_binomial_lccdf(n_ref, 1, beta);
+  } else {
+    const auto beta = elt_divide(theta_ref, subtract(1.0, theta_ref));
+    return neg_binomial_lccdf(n_ref, 1, beta);
+  }
+}
+
+}  // namespace math
+}  // namespace stan
+#endif
diff --git a/stan/math/prim/prob/geometric_lcdf.hpp b/stan/math/prim/prob/geometric_lcdf.hpp
new file mode 100644
index 00000000000..0032d0c7ec4
--- /dev/null
+++ b/stan/math/prim/prob/geometric_lcdf.hpp
@@ -0,0 +1,90 @@
+#ifndef STAN_MATH_PRIM_PROB_GEOMETRIC_LCDF_HPP
+#define STAN_MATH_PRIM_PROB_GEOMETRIC_LCDF_HPP
+
+#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/scalar_seq_view.hpp>
+#include <stan/math/prim/fun/size.hpp>
+#include <stan/math/prim/fun/size_zero.hpp>
+#include <stan/math/prim/fun/value_of.hpp>
+#include <vector>
+#include <stan/math/prim/prob/neg_binomial_lcdf.hpp>
+#include <stan/math/prim/fun/elt_divide.hpp>
+#include <stan/math/prim/fun/subtract.hpp>
+
+namespace stan {
+namespace math {
+
+/** \ingroup prob_dists
+ * Returns the log CDF of the geometric distribution. Given containers of
+ * matching sizes, returns the log sum of probabilities.
+ *
+ * Delegates to the negative binomial log CDF with alpha = 1 and
+ * beta = theta / (1 - theta).
+ *
+ * @tparam T_n type of outcome variable
+ * @tparam T_prob type of success probability parameter
+ *
+ * @param n outcome variable (number of failures before first success)
+ * @param theta success probability parameter
+ * @return log probability or log sum of probabilities
+ * @throw std::domain_error if theta is not in (0, 1]
+ * @throw std::invalid_argument if container sizes mismatch
+ */
+template <typename T_n, typename T_prob>
+inline return_type_t<T_prob> geometric_lcdf(const T_n& n, const T_prob& theta) {
+  using T_n_ref = ref_type_t<T_n>;
+  using T_prob_ref = ref_type_t<T_prob>;
+  static constexpr const char* function = "geometric_lcdf";
+
+  check_consistent_sizes(function, "Outcome variable", n,
+                         "Success probability parameter", theta);
+  if (size_zero(n, theta)) {
+    return 0.0;
+  }
+
+  T_n_ref n_ref = n;
+  T_prob_ref theta_ref = theta;
+  check_bounded(function, "Success probability parameter", value_of(theta_ref),
+                0.0, 1.0);
+
+  scalar_seq_view<T_n_ref> n_vec(n_ref);
+  for (int i = 0; i < stan::math::size(n); i++) {
+    if (n_vec.val(i) < 0) {
+      return negative_infinity();
+    }
+  }
+
+  // theta = 1 => log CDF is 0 for n >= 0
+  scalar_seq_view<T_prob_ref> theta_vec(theta_ref);
+  bool all_theta_one = true;
+  for (size_t i = 0; i < stan::math::size(theta); i++) {
+    if (value_of(theta_vec[i]) != 1.0) {
+      all_theta_one = false;
+      break;
+    }
+  }
+  if (all_theta_one) {
+    return 0.0;
+  }
+
+  if constexpr (is_stan_scalar_v<T_prob>) {
+    const auto beta = theta_ref / (1.0 - theta_ref);
+    return neg_binomial_lcdf(n_ref, 1, beta);
+  } else if constexpr (is_std_vector_v<T_prob>) {
+    std::vector<value_type_t<T_prob>> beta;
+    beta.reserve(stan::math::size(theta));
+    for (size_t i = 0; i < stan::math::size(theta); i++) {
+      beta.push_back(theta_vec[i] / (1.0 - theta_vec[i]));
+    }
+    return neg_binomial_lcdf(n_ref, 1, beta);
+  } else {
+    const auto beta = elt_divide(theta_ref, subtract(1.0, theta_ref));
+    return neg_binomial_lcdf(n_ref, 1, beta);
+  }
+}
+
+}  // namespace math
+}  // namespace stan
+#endif
diff --git a/stan/math/prim/prob/geometric_lpmf.hpp b/stan/math/prim/prob/geometric_lpmf.hpp
new file mode 100644
index 00000000000..df0a8f2a03e
--- /dev/null
+++ b/stan/math/prim/prob/geometric_lpmf.hpp
@@ -0,0 +1,98 @@
+#ifndef STAN_MATH_PRIM_PROB_GEOMETRIC_LPMF_HPP
+#define STAN_MATH_PRIM_PROB_GEOMETRIC_LPMF_HPP
+
+#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/scalar_seq_view.hpp>
+#include <stan/math/prim/fun/size.hpp>
+#include <stan/math/prim/fun/size_zero.hpp>
+#include <stan/math/prim/fun/value_of.hpp>
+#include <vector>
+#include <stan/math/prim/prob/neg_binomial_lpmf.hpp>
+#include <stan/math/prim/fun/elt_divide.hpp>
+#include <stan/math/prim/fun/subtract.hpp>
+
+namespace stan {
+namespace math {
+
+/** \ingroup prob_dists
+ * Returns the log PMF of the geometric distribution. If containers
+ * of matching sizes are supplied, returns the log sum of probabilities.
+ *
+ * The geometric distribution is a special case of the negative
+ * binomial distribution with alpha = 1 and beta = theta / (1 - theta).
+ *
+ * @tparam T_n type of outcome variable
+ * @tparam T_prob type of success probability parameter
+ *
+ * @param n outcome variable (number of failures before first success)
+ * @param theta success probability parameter
+ * @return log probability or log sum of probabilities
+ * @throw std::domain_error if theta is not in (0, 1]
+ * @throw std::domain_error if n is negative
+ * @throw std::invalid_argument if container sizes mismatch
+ */
+template <bool propto, typename T_n, typename T_prob,
+          require_all_not_nonscalar_prim_or_rev_kernel_expression_t<
+              T_n, T_prob>* = nullptr>
+inline return_type_t<T_prob> geometric_lpmf(const T_n& n, const T_prob& theta) {
+  using T_n_ref = ref_type_t<T_n>;
+  using T_prob_ref = ref_type_t<T_prob>;
+  static constexpr const char* function = "geometric_lpmf";
+  check_consistent_sizes(function, "Outcome variable", n,
+                         "Success probability parameter", theta);
+  if (size_zero(n, theta)) {
+    return 0.0;
+  }
+
+  T_n_ref n_ref = n;
+  T_prob_ref theta_ref = theta;
+  check_nonnegative(function, "Outcome variable", n_ref);
+  check_bounded(function, "Success probability parameter", value_of(theta_ref),
+                0.0, 1.0);
+
+  // theta = 1 => deterministic: P(0) = 1, P(n>0) = 0
+  // Cannot delegate since beta = theta / (1 - theta) diverges
+  scalar_seq_view<T_n_ref> n_vec(n_ref);
+  scalar_seq_view<T_prob_ref> theta_vec(theta_ref);
+  const size_t max_sz = max_size(n_ref, theta_ref);
+  bool all_theta_one = true;
+  for (size_t i = 0; i < max_sz; i++) {
+    if (value_of(theta_vec[i]) == 1.0) {
+      if (n_vec[i] > 0) {
+        return negative_infinity();
+      }
+    } else {
+      all_theta_one = false;
+    }
+  }
+  if (all_theta_one) {
+    return 0.0;
+  }
+
+  // geometric(theta) = neg_binomial(1, theta / (1 - theta))
+  if constexpr (is_stan_scalar_v<T_prob>) {
+    const auto beta = theta_ref / (1.0 - theta_ref);
+    return neg_binomial_lpmf<propto>(n_ref, 1, beta);
+  } else if constexpr (is_std_vector_v<T_prob>) {
+    std::vector<value_type_t<T_prob>> beta;
+    beta.reserve(stan::math::size(theta));
+    for (size_t i = 0; i < stan::math::size(theta); i++) {
+      beta.push_back(theta_vec[i] / (1.0 - theta_vec[i]));
+    }
+    return neg_binomial_lpmf<propto>(n_ref, 1, beta);
+  } else {
+    const auto beta = elt_divide(theta_ref, subtract(1.0, theta_ref));
+    return neg_binomial_lpmf<propto>(n_ref, 1, beta);
+  }
+}
+
+template <typename T_n, typename T_prob>
+inline return_type_t<T_prob> geometric_lpmf(const T_n& n, const T_prob& theta) {
+  return geometric_lpmf<false>(n, theta);
+}
+
+}  // namespace math
+}  // namespace stan
+#endif
diff --git a/stan/math/prim/prob/geometric_rng.hpp b/stan/math/prim/prob/geometric_rng.hpp
new file mode 100644
index 00000000000..fa7619a86c6
--- /dev/null
+++ b/stan/math/prim/prob/geometric_rng.hpp
@@ -0,0 +1,63 @@
+#ifndef STAN_MATH_PRIM_PROB_GEOMETRIC_RNG_HPP
+#define STAN_MATH_PRIM_PROB_GEOMETRIC_RNG_HPP
+
+#include <stan/math/prim/meta.hpp>
+#include <stan/math/prim/err.hpp>
+#include <stan/math/prim/fun/as_array_or_scalar.hpp>
+#include <stan/math/prim/fun/to_ref.hpp>
+#include <stan/math/prim/fun/value_of.hpp>
+#include <stan/math/prim/prob/neg_binomial_rng.hpp>
+#include <utility>
+#include <vector>
+
+namespace stan {
+namespace math {
+
+/** \ingroup prob_dists
+ * Return a geometric random variate with the given success probability
+ * parameter, using the given random number generator.
+ *
+ * Delegates to neg_binomial_rng with alpha = 1 and
+ * beta = theta / (1 - theta), with a special case for theta = 1.
+ *
+ * @tparam T_prob type of success probability parameter
+ * @tparam RNG type of random number generator
+ *
+ * @param theta (Sequence of) success probability parameter(s)
+ * @param rng random number generator
+ * @return (Sequence of) geometric random variate(s)
+ * @throw std::domain_error if theta is not in (0, 1]
+ */
+template <typename T_prob, typename RNG>
+inline auto geometric_rng(T_prob&& theta, RNG& rng) {
+  static constexpr const char* function = "geometric_rng";
+  decltype(auto) theta_ref = to_ref(std::forward<T_prob>(theta));
+  check_positive_finite(function, "Success probability parameter",
+                        value_of(theta_ref));
+  check_less_or_equal(function, "Success probability parameter",
+                      value_of(theta_ref), 1.0);
+
+  if constexpr (is_stan_scalar_v<T_prob>) {
+    if (theta_ref == 1.0) {
+      return 0;
+    }
+    double beta = theta_ref / (1.0 - theta_ref);
+    return neg_binomial_rng(1, beta, rng);
+  } else {
+    auto theta_arr = as_array_or_scalar(theta_ref);
+    std::vector<int> result(theta_arr.size());
+    for (size_t i = 0; i < theta_arr.size(); i++) {
+      if (theta_arr[i] == 1.0) {
+        result[i] = 0;
+      } else {
+        double beta_i = theta_arr[i] / (1.0 - theta_arr[i]);
+        result[i] = neg_binomial_rng(1, beta_i, rng);
+      }
+    }
+    return result;
+  }
+}
+
+}  // namespace math
+}  // namespace stan
+#endif
diff --git a/test/prob/geometric/geometric_ccdf_log_test.hpp b/test/prob/geometric/geometric_ccdf_log_test.hpp
new file mode 100644
index 00000000000..6a5595c47d7
--- /dev/null
+++ b/test/prob/geometric/geometric_ccdf_log_test.hpp
@@ -0,0 +1,63 @@
+// Arguments: Ints, Doubles
+#include <stan/math/prim/prob/geometric_lccdf.hpp>
+#include <stan/math/prim/fun/log1m.hpp>
+
+using std::numeric_limits;
+using std::vector;
+
+class AgradCcdfLogGeometric : public AgradCcdfLogTest {
+ public:
+  void valid_values(vector<vector<double>>& parameters,
+                    vector<double>& ccdf_log) {
+    vector<double> param(2);
+
+    // lccdf(n=0|theta=0.5) = log(0.5) = -0.693147...
+    param[0] = 0;    // n
+    param[1] = 0.5;  // theta
+    parameters.push_back(param);
+    ccdf_log.push_back(-0.69314718055994528623);
+
+    // lccdf(n=2|theta=0.5) = log(0.125) = -2.079441...
+    param[0] = 2;    // n
+    param[1] = 0.5;  // theta
+    parameters.push_back(param);
+    ccdf_log.push_back(-2.07944154167983574766);
+
+    // lccdf(n=4|theta=0.3) = 5*log(0.7) = -1.783374...
+    param[0] = 4;    // n
+    param[1] = 0.3;  // theta
+    parameters.push_back(param);
+    ccdf_log.push_back(-1.78337471969366179181);
+  }
+
+  void invalid_values(vector<size_t>& index, vector<double>& value) {
+    // theta
+    index.push_back(1U);
+    value.push_back(-0.1);
+
+    index.push_back(1U);
+    value.push_back(1.1);
+  }
+
+  bool has_lower_bound() { return false; }
+
+  bool has_upper_bound() { return false; }
+
+  template <class T_n, class T_prob, typename T2, typename T3, typename T4,
+            typename T5>
+  stan::return_type_t<T_n, T_prob> ccdf_log(const T_n& n, const T_prob& theta,
+                                            const T2&, const T3&, const T4&,
+                                            const T5&) {
+    return stan::math::geometric_lccdf(n, theta);
+  }
+
+  template <class T_n, class T_prob, typename T2, typename T3, typename T4,
+            typename T5>
+  stan::return_type_t<T_n, T_prob> ccdf_log_function(const T_n& n,
+                                                     const T_prob& theta,
+                                                     const T2&, const T3&,
+                                                     const T4&, const T5&) {
+    using stan::math::log1m;
+    return (n + 1.0) * log1m(theta);
+  }
+};
diff --git a/test/prob/geometric/geometric_cdf_log_test.hpp b/test/prob/geometric/geometric_cdf_log_test.hpp
new file mode 100644
index 00000000000..9843a30c3c0
--- /dev/null
+++ b/test/prob/geometric/geometric_cdf_log_test.hpp
@@ -0,0 +1,64 @@
+// Arguments: Ints, Doubles
+#include <stan/math/prim/prob/geometric_lcdf.hpp>
+#include <stan/math/prim/fun/log1m.hpp>
+
+using std::numeric_limits;
+using std::vector;
+
+class AgradCdfLogGeometric : public AgradCdfLogTest {
+ public:
+  void valid_values(vector<vector<double>>& parameters,
+                    vector<double>& cdf_log) {
+    vector<double> param(2);
+
+    // lcdf(n=0|theta=0.5) = log(0.5) = -0.693147...
+    param[0] = 0;    // n
+    param[1] = 0.5;  // theta
+    parameters.push_back(param);
+    cdf_log.push_back(-0.69314718055994528623);
+
+    // lcdf(n=2|theta=0.5) = log(0.875) = -0.133531...
+    param[0] = 2;    // n
+    param[1] = 0.5;  // theta
+    parameters.push_back(param);
+    cdf_log.push_back(-0.13353139262452262681);
+
+    // lcdf(n=4|theta=0.3) = log(0.83193) = -0.184006...
+    param[0] = 4;    // n
+    param[1] = 0.3;  // theta
+    parameters.push_back(param);
+    cdf_log.push_back(-0.18400697631582843550);
+  }
+
+  void invalid_values(vector<size_t>& index, vector<double>& value) {
+    // theta
+    index.push_back(1U);
+    value.push_back(-0.1);
+
+    index.push_back(1U);
+    value.push_back(1.1);
+  }
+
+  bool has_lower_bound() { return false; }
+
+  bool has_upper_bound() { return false; }
+
+  template <class T_n, class T_prob, typename T2, typename T3, typename T4,
+            typename T5>
+  stan::return_type_t<T_n, T_prob> cdf_log(const T_n& n, const T_prob& theta,
+                                           const T2&, const T3&, const T4&,
+                                           const T5&) {
+    return stan::math::geometric_lcdf(n, theta);
+  }
+
+  template <class T_n, class T_prob, typename T2, typename T3, typename T4,
+            typename T5>
+  stan::return_type_t<T_n, T_prob> cdf_log_function(const T_n& n,
+                                                    const T_prob& theta,
+                                                    const T2&, const T3&,
+                                                    const T4&, const T5&) {
+    using stan::math::log1m;
+    using stan::math::log1m_exp;
+    return log1m_exp((n + 1.0) * log1m(theta));
+  }
+};
diff --git a/test/prob/geometric/geometric_cdf_test.hpp b/test/prob/geometric/geometric_cdf_test.hpp
new file mode 100644
index 00000000000..cfd36f66c4d
--- /dev/null
+++ b/test/prob/geometric/geometric_cdf_test.hpp
@@ -0,0 +1,62 @@
+// Arguments: Ints, Doubles
+#include <stan/math/prim/prob/geometric_cdf.hpp>
+#include <stan/math/prim/fun/log1m.hpp>
+
+using std::numeric_limits;
+using std::vector;
+
+class AgradCdfGeometric : public AgradCdfTest {
+ public:
+  void valid_values(vector<vector<double>>& parameters, vector<double>& cdf) {
+    vector<double> param(2);
+
+    // CDF(n=0|theta=0.5) = 1 - (0.5)^1 = 0.5
+    param[0] = 0;    // n
+    param[1] = 0.5;  // theta
+    parameters.push_back(param);
+    cdf.push_back(0.5);
+
+    // CDF(n=2|theta=0.5) = 1 - (0.5)^3 = 0.875
+    param[0] = 2;    // n
+    param[1] = 0.5;  // theta
+    parameters.push_back(param);
+    cdf.push_back(0.875);
+
+    // CDF(n=4|theta=0.3) = 1 - 0.7^5 = 0.83193
+    param[0] = 4;    // n
+    param[1] = 0.3;  // theta
+    parameters.push_back(param);
+    cdf.push_back(0.83193);
+  }
+
+  void invalid_values(vector<size_t>& index, vector<double>& value) {
+    // theta
+    index.push_back(1U);
+    value.push_back(-0.1);
+
+    index.push_back(1U);
+    value.push_back(1.1);
+  }
+
+  bool has_lower_bound() { return false; }
+
+  bool has_upper_bound() { return false; }
+
+  template <class T_n, class T_prob, typename T2, typename T3, typename T4,
+            typename T5>
+  stan::return_type_t<T_n, T_prob> cdf(const T_n& n, const T_prob& theta,
+                                       const T2&, const T3&, const T4&,
+                                       const T5&) {
+    return stan::math::geometric_cdf(n, theta);
+  }
+
+  template <class T_n, class T_prob, typename T2, typename T3, typename T4,
+            typename T5>
+  stan::return_type_t<T_n, T_prob> cdf_function(const T_n& n,
+                                                const T_prob& theta, const T2&,
+                                                const T3&, const T4&,
+                                                const T5&) {
+    using stan::math::log1m;
+    return 1.0 - stan::math::exp((n + 1.0) * log1m(theta));
+  }
+};
diff --git a/test/prob/geometric/geometric_test.hpp b/test/prob/geometric/geometric_test.hpp
new file mode 100644
index 00000000000..8acf08fd22f
--- /dev/null
+++ b/test/prob/geometric/geometric_test.hpp
@@ -0,0 +1,73 @@
+// Arguments: Ints, Doubles
+#include <stan/math/prim/prob/geometric_lpmf.hpp>
+#include <stan/math/prim/fun/log1m.hpp>
+
+using stan::math::var;
+using std::numeric_limits;
+using std::vector;
+
+class AgradDistributionsGeometric : public AgradDistributionTest {
+ public:
+  void valid_values(vector<vector<double>>& parameters,
+                    vector<double>& log_prob) {
+    vector<double> param(2);
+
+    // lpmf(n=0|theta=0.5) = log(0.5) = -0.693147...
+    param[0] = 0;    // n
+    param[1] = 0.5;  // theta
+    parameters.push_back(param);
+    log_prob.push_back(-0.69314718055994528623);
+
+    // lpmf(n=2|theta=0.5) = log(0.5) + 2*log(0.5) = -2.079441...
+    param[0] = 2;    // n
+    param[1] = 0.5;  // theta
+    parameters.push_back(param);
+    log_prob.push_back(-2.07944154167983574766);
+
+    // lpmf(n=4|theta=0.3) = log(0.3) + 4*log(0.7) = -2.630672...
+    param[0] = 4;    // n
+    param[1] = 0.3;  // theta
+    parameters.push_back(param);
+    log_prob.push_back(-2.63067258008086568566);
+  }
+
+  void invalid_values(vector<size_t>& index, vector<double>& value) {
+    // n
+    index.push_back(0U);
+    value.push_back(-1);
+
+    // theta
+    index.push_back(1U);
+    value.push_back(-0.1);
+
+    index.push_back(1U);
+    value.push_back(1.1);
+  }
+
+  template <class T_n, class T_prob, typename T2, typename T3, typename T4,
+            typename T5>
+  stan::return_type_t<T_prob> log_prob(const T_n& n, const T_prob& theta,
+                                       const T2&, const T3&, const T4&,
+                                       const T5&) {
+    return stan::math::geometric_lpmf(n, theta);
+  }
+
+  template <bool propto, class T_n, class T_prob, typename T2, typename T3,
+            typename T4, typename T5>
+  stan::return_type_t<T_prob> log_prob(const T_n& n, const T_prob& theta,
+                                       const T2&, const T3&, const T4&,
+                                       const T5&) {
+    return stan::math::geometric_lpmf<propto>(n, theta);
+  }
+
+  template <class T_n, class T_prob, typename T2, typename T3, typename T4,
+            typename T5>
+  stan::return_type_t<T_prob> log_prob_function(const T_n& n,
+                                                const T_prob& theta, const T2&,
+                                                const T3&, const T4&,
+                                                const T5&) {
+    using stan::math::log1m;
+    using std::log;
+    return log(theta) + n * log1m(theta);
+  }
+};
diff --git a/test/unit/math/prim/prob/geometric_test.cpp b/test/unit/math/prim/prob/geometric_test.cpp
new file mode 100644
index 00000000000..e564f78cfbf
--- /dev/null
+++ b/test/unit/math/prim/prob/geometric_test.cpp
@@ -0,0 +1,103 @@
+#include <stan/math/prim.hpp>
+#include <stan/math/prim/prob/geometric_rng.hpp>
+#include <test/unit/math/prim/prob/vector_rng_test_helper.hpp>
+#include <test/unit/math/prim/prob/VectorIntRNGTestRig.hpp>
+#include <gtest/gtest.h>
+#include <boost/random/mersenne_twister.hpp>
+#include <limits>
+#include <vector>
+
+class GeometricTestRig : public VectorIntRNGTestRig {
+ public:
+  GeometricTestRig()
+      : VectorIntRNGTestRig(10000, 10, {0, 1, 2, 3, 4, 5, 6, 7, 8, 9},
+                            {0.1, 0.3, 0.5, 0.7, 0.9}, {1}, {-0.1, 1.1},
+                            {-1, 0, 2}) {}
+
+  template <typename T1, typename T2, typename T3, typename T_rng>
+  auto generate_samples(const T1& theta, const T2&, const T3&,
+                        T_rng& rng) const {
+    return stan::math::geometric_rng(theta, rng);
+  }
+
+  template <typename T1>
+  double pmf(int y, T1 theta, double, double) const {
+    if (y < 0) {
+      return 0.0;
+    }
+    return std::exp(stan::math::geometric_lpmf(y, theta));
+  }
+};
+
+TEST(ProbDistributionsGeometric, errorCheck) {
+  check_dist_throws_all_types(GeometricTestRig());
+}
+
+TEST(ProbDistributionsGeometric, distributionCheck) {
+  check_counts_real(GeometricTestRig());
+}
+
+TEST(ProbDistributionsGeometric, error_check) {
+  boost::random::mt19937 rng;
+
+  EXPECT_NO_THROW(stan::math::geometric_rng(0.5, rng));
+  EXPECT_NO_THROW(stan::math::geometric_rng(0.1, rng));
+  EXPECT_NO_THROW(stan::math::geometric_rng(0.9, rng));
+  EXPECT_NO_THROW(stan::math::geometric_rng(1.0, rng));
+
+  EXPECT_THROW(stan::math::geometric_rng(0.0, rng), std::domain_error);
+  EXPECT_THROW(stan::math::geometric_rng(-0.5, rng), std::domain_error);
+
+  EXPECT_THROW(stan::math::geometric_rng(stan::math::positive_infinity(), rng),
+               std::domain_error);
+  EXPECT_THROW(stan::math::geometric_rng(stan::math::not_a_number(), rng),
+               std::domain_error);
+}
+
+TEST(ProbDistributionsGeometric, lpmf_values) {
+  // P(n=0|theta=0.5) = 0.5
+  EXPECT_NEAR(stan::math::geometric_lpmf(0, 0.5), std::log(0.5), 1e-10);
+  // P(n=2|theta=0.5) = 0.5 * 0.5^2 = 0.125
+  EXPECT_NEAR(stan::math::geometric_lpmf(2, 0.5), std::log(0.125), 1e-10);
+  // P(n=0|theta=1.0) = 1.0
+  EXPECT_NEAR(stan::math::geometric_lpmf(0, 1.0), 0.0, 1e-10);
+  // P(n=1|theta=1.0) = 0
+  EXPECT_EQ(stan::math::geometric_lpmf(1, 1.0),
+            stan::math::negative_infinity());
+  // P(n=0|theta=0.3) = 0.3
+  EXPECT_NEAR(stan::math::geometric_lpmf(0, 0.3), std::log(0.3), 1e-10);
+  // P(n=3|theta=0.3) = 0.3 * 0.7^3 = 0.1029
+  EXPECT_NEAR(stan::math::geometric_lpmf(3, 0.3), std::log(0.1029), 1e-10);
+}
+
+TEST(ProbDistributionsGeometric, cdf_values) {
+  // CDF(n=0|theta=0.5) = 1 - (1-0.5)^1 = 0.5
+  EXPECT_NEAR(stan::math::geometric_cdf(0, 0.5), 0.5, 1e-10);
+  // CDF(n=2|theta=0.5) = 1 - (0.5)^3 = 0.875
+  EXPECT_NEAR(stan::math::geometric_cdf(2, 0.5), 0.875, 1e-10);
+  // CDF(n=0|theta=1.0) = 1.0
+  EXPECT_NEAR(stan::math::geometric_cdf(0, 1.0), 1.0, 1e-10);
+  // CDF(n=4|theta=0.3) = 1 - 0.7^5 = 0.83193
+  EXPECT_NEAR(stan::math::geometric_cdf(4, 0.3), 0.83193, 1e-5);
+}
+
+TEST(ProbDistributionsGeometric, cdf_below_support) {
+  EXPECT_DOUBLE_EQ(stan::math::geometric_cdf(-1, 0.5), 0.0);
+  EXPECT_DOUBLE_EQ(stan::math::geometric_lcdf(-1, 0.5),
+                   stan::math::negative_infinity());
+}
+
+TEST(ProbDistributionsGeometric, lccdf_values) {
+  // CCDF(n=0|theta=0.5) = (1-0.5)^1 = 0.5
+  EXPECT_NEAR(stan::math::geometric_lccdf(0, 0.5), std::log(0.5), 1e-10);
+  // CCDF(n=2|theta=0.5) = (0.5)^3 = 0.125
+  EXPECT_NEAR(stan::math::geometric_lccdf(2, 0.5), std::log(0.125), 1e-10);
+}
+
+TEST(ProbDistributionsGeometric, rng_theta_one) {
+  boost::random::mt19937 rng;
+  // theta=1.0 should always return 0 (always succeed on first trial)
+  for (int i = 0; i < 100; i++) {
+    EXPECT_EQ(stan::math::geometric_rng(1.0, rng), 0);
+  }
+}