diff --git a/stan/math/fwd/prob.hpp b/stan/math/fwd/prob.hpp
index 5d280c06525..c390e191db5 100644
--- a/stan/math/fwd/prob.hpp
+++ b/stan/math/fwd/prob.hpp
@@ -5,5 +5,6 @@
 #include <stan/math/fwd/fun/Eigen_NumTraits.hpp>
 
 #include <stan/math/fwd/prob/std_normal_log_qf.hpp>
+#include <stan/math/fwd/prob/student_t_qf.hpp>
 
 #endif
diff --git a/stan/math/fwd/prob/student_t_qf.hpp b/stan/math/fwd/prob/student_t_qf.hpp
new file mode 100644
index 00000000000..0b87c60171d
--- /dev/null
+++ b/stan/math/fwd/prob/student_t_qf.hpp
@@ -0,0 +1,88 @@
+#ifndef STAN_MATH_FWD_PROB_STUDENT_T_QF_HPP
+#define STAN_MATH_FWD_PROB_STUDENT_T_QF_HPP
+
+#include <stan/math/fwd/meta.hpp>
+#include <stan/math/fwd/fun/sqrt.hpp>
+#include <stan/math/fwd/fun/inv_inc_beta.hpp>
+#include <stan/math/prim/prob/student_t_lpdf.hpp>
+
+namespace stan {
+namespace math {
+
+template <typename T_p, typename T_nu, typename T_mu, typename T_sigma,
+          typename FvarT = return_type_t<T_p, T_mu, T_sigma, T_nu>,
+          require_all_stan_scalar_t<T_p, T_mu, T_sigma, T_nu>* = nullptr,
+          require_fvar_t<FvarT>* = nullptr>
+FvarT student_t_qf(const T_p& p, const T_nu& nu, const T_mu& mu,
+                   const T_sigma& sigma) {
+  static constexpr const char* function = "student_t_qf";
+  using T_partials = partials_type_t<FvarT>;
+
+  const auto& p_val = value_of(p);
+  const auto& nu_val = value_of(nu);
+  const auto& mu_val = value_of(mu);
+  const auto& sigma_val = value_of(sigma);
+
+  check_nonnegative(function, "Degrees of freedom parameter", nu_val);
+  check_positive(function, "Scale parameter", sigma_val);
+  check_bounded(function, "Probability parameter", p_val, 0.0, 1.0);
+
+  if (p_val == 0.0) {
+    return {NEGATIVE_INFTY, 0.0};
+  } else if (p_val == 1.0) {
+    return {INFTY, 0.0};
+  } else if (p_val == 0.5) {
+    return {mu_val, 0.0};
+  }
+
+  const T_partials p_val_flip = p_val < 0.5 ? p_val : 1.0 - p_val;
+  const double p_sign = value_of_rec(p_val) < 0.5 ? -1.0 : 1.0;
+
+  T_partials ibeta_arg = inv_inc_beta(0.5 * nu_val, 0.5, 2 * p_val_flip);
+  T_partials rtn_val = mu_val
+                       + p_sign * sigma_val * sqrt(nu_val)
+                             * math::sqrt(-1.0 + 1.0 / ibeta_arg);
+
+  FvarT rtn(rtn_val, 0.0);
+
+  if constexpr (!is_constant<T_p>::value) {
+    rtn.d_ += p.d_ * exp(-student_t_lpdf(rtn_val, nu_val, mu_val, sigma_val));
+  }
+
+  if constexpr (!is_constant<T_nu>::value) {
+    const T_partials half_nu = nu_val / 2.0;
+    Eigen::Matrix<T_partials, -1, 1> hyper_arg_a(3);
+    hyper_arg_a << 0.5, half_nu, half_nu;
+    Eigen::Matrix<T_partials, -1, 1> hyper_arg_b(2);
+    hyper_arg_b << 1 + half_nu, 1 + half_nu;
+    const T_partials hyper_arg
+        = hypergeometric_pFq(hyper_arg_a, hyper_arg_b, ibeta_arg);
+    const T_partials hyper2f1
+        = hypergeometric_2F1(1, (1 + nu_val) / 2, (2 + nu_val) / 2, ibeta_arg);
+    const T_partials digamma_a1 = digamma(half_nu);
+    const T_partials digamma_a2 = digamma((1 + nu_val) / 2);
+
+    const T_partials arg_1 = (4 * hyper_arg * sqrt(1 - ibeta_arg)) / nu_val;
+    const T_partials arg_2 = -2 * hyper2f1 * (-1 + ibeta_arg)
+                             * (log(ibeta_arg) - digamma_a1 + digamma_a2);
+
+    const T_partials num1 = sigma_val * (-2 + (2 - arg_1 + arg_2) / ibeta_arg);
+    const T_partials den1 = 4 * sqrt(nu_val) * sqrt(-1 + 1 / ibeta_arg);
+    rtn.d_ += nu.d_ * p_sign * num1 / den1;
+  }
+
+  if constexpr (!is_constant<T_mu>::value) {
+    rtn.d_ += mu.d_;
+  }
+
+  if constexpr (!is_constant<T_sigma>::value) {
+    rtn.d_ += sigma.d_ * p_sign * sqrt(nu_val)
+              * math::sqrt(-1.0 + 1.0 / ibeta_arg);
+  }
+
+  return rtn;
+}
+}  // namespace math
+}  // namespace stan
+
+#endif
diff --git a/stan/math/prim/meta.hpp b/stan/math/prim/meta.hpp
index dd6e3053eb6..5a36e92ba69 100644
--- a/stan/math/prim/meta.hpp
+++ b/stan/math/prim/meta.hpp
@@ -71,6 +71,7 @@
 #include <stan/math/prim/meta/ad_promotable.hpp>
 #include <stan/math/prim/meta/append_return_type.hpp>
 #include <stan/math/prim/meta/base_type.hpp>
+#include <stan/math/prim/meta/common_container_type.hpp>
 #include <stan/math/prim/meta/contains_std_vector.hpp>
 #include <stan/math/prim/meta/contains_tuple.hpp>
 #include <stan/math/prim/meta/error_index.hpp>
diff --git a/stan/math/prim/meta/common_container_type.hpp b/stan/math/prim/meta/common_container_type.hpp
new file mode 100644
index 00000000000..df74fddff5a
--- /dev/null
+++ b/stan/math/prim/meta/common_container_type.hpp
@@ -0,0 +1,77 @@
+#ifndef STAN_MATH_PRIM_META_COMMON_CONTAINER_TYPE_HPP
+#define STAN_MATH_PRIM_META_COMMON_CONTAINER_TYPE_HPP
+
+#include <stan/math/prim/meta/is_container.hpp>
+#include <stan/math/prim/meta/is_tuple.hpp>
+#include <stan/math/prim/meta/is_detected.hpp>
+#include <stan/math/prim/meta/is_stan_scalar.hpp>
+#include <stan/math/prim/meta/is_var_matrix.hpp>
+#include <stan/math/prim/meta/plain_type.hpp>
+#include <stan/math/prim/meta/return_type.hpp>
+#include <stan/math/prim/meta/promote_scalar_type.hpp>
+#include <type_traits>
+
+namespace stan {
+namespace internal {
+template <typename T1, typename T2, typename = void, typename = void>
+struct common_container_type_impl;
+
+template <typename T1, typename T2>
+struct common_container_type_impl<T1, T2, require_stan_scalar_t<T1>,
+                                  require_stan_scalar_t<T2>> {
+  using type = return_type_t<T1, T2>;
+};
+
+template <typename T1, typename T2>
+struct common_container_type_impl<T1, T2, require_container_t<T1>,
+                                  require_container_t<T2>> {
+  using return_t = return_type_t<T1, T2>;
+  using container_type_1 = math::promote_scalar_t<return_t, plain_type_t<T1>>;
+  using container_type_2 = math::promote_scalar_t<return_t, plain_type_t<T2>>;
+  using type = std::conditional_t<
+      std::is_same<container_type_1, container_type_2>::value, container_type_1,
+      void>;
+};
+
+template <typename T1, typename T2>
+struct common_container_type_impl<T1, T2, require_stan_scalar_t<T1>,
+                                  require_container_t<T2>> {
+  using type = math::promote_scalar_t<return_type_t<T1, T2>, plain_type_t<T2>>;
+};
+
+template <typename T1, typename T2>
+struct common_container_type_impl<T1, T2, require_container_t<T1>,
+                                  require_stan_scalar_t<T2>> {
+  using type = math::promote_scalar_t<return_type_t<T1, T2>, plain_type_t<T1>>;
+};
+}  // namespace internal
+
+template <typename... Ts>
+struct common_container_type;
+
+template <typename T>
+struct common_container_type<T> {
+  using type = typename internal::common_container_type_impl<
+      T, double>::type;  // Use double for base case
+};
+
+/**
+ * Determine the common container type for a variadic list of types.
+ * If all types are scalars, then the common scalar type is returned.
+ * If all container types the same, but not necessarily the same scalar type,
+ * the common container type with the common scalar type is returned.
+ *
+ * If different container types are present, the result is `void`.
+ */
+template <typename T1, typename... Ts>
+struct common_container_type<T1, Ts...> {
+  using type = typename internal::common_container_type_impl<
+      T1, typename common_container_type<Ts...>::type>::type;
+};
+
+template <typename... Ts>
+using common_container_t = typename common_container_type<Ts...>::type;
+
+}  // namespace stan
+
+#endif  // STAN_MATH_PRIM_META_PLAIN_TYPE_HPP
diff --git a/stan/math/prim/prob.hpp b/stan/math/prim/prob.hpp
index e04c7e87172..4e3207f0c68 100644
--- a/stan/math/prim/prob.hpp
+++ b/stan/math/prim/prob.hpp
@@ -287,6 +287,7 @@
 #include <stan/math/prim/prob/student_t_lccdf.hpp>
 #include <stan/math/prim/prob/student_t_lcdf.hpp>
 #include <stan/math/prim/prob/student_t_lpdf.hpp>
+#include <stan/math/prim/prob/student_t_qf.hpp>
 #include <stan/math/prim/prob/student_t_rng.hpp>
 #include <stan/math/prim/prob/uniform_ccdf_log.hpp>
 #include <stan/math/prim/prob/uniform_cdf.hpp>
diff --git a/stan/math/prim/prob/student_t_qf.hpp b/stan/math/prim/prob/student_t_qf.hpp
new file mode 100644
index 00000000000..8a134498d75
--- /dev/null
+++ b/stan/math/prim/prob/student_t_qf.hpp
@@ -0,0 +1,99 @@
+#ifndef STAN_MATH_PRIM_PROB_STUDENT_T_QF_HPP
+#define STAN_MATH_PRIM_PROB_STUDENT_T_QF_HPP
+
+#include <stan/math/prim/meta.hpp>
+#include <stan/math/prim/err.hpp>
+#include <stan/math/prim/fun/sqrt.hpp>
+#include <stan/math/prim/fun/inv_inc_beta.hpp>
+#include <stan/math/prim/fun/max_size.hpp>
+
+namespace stan {
+namespace math {
+
+/**
+ * The quantile function of the Student's t-distribution.
+ *
+ * @tparam T_p type of the probability parameter
+ * @tparam T_nu type of the degrees of freedom parameter
+ * @tparam T_mu type of the location parameter
+ * @tparam T_sigma type of the scale parameter
+ * @param p Probability in the range [0, 1].
+ * @param nu Degrees of freedom, must be non-negative.
+ * @param mu Location parameter.
+ * @param sigma Scale parameter, must be positive.
+ * @return Quantile function value.
+ * @throw std::domain_error if `nu` is negative or `sigma` is not positive,
+ * or if `p` is not in [0, 1].
+ */
+template <typename T_p, typename T_nu, typename T_mu, typename T_sigma,
+          require_all_stan_scalar_t<T_p, T_nu, T_mu, T_sigma>* = nullptr,
+          require_all_arithmetic_t<T_p, T_nu, T_mu, T_sigma>* = nullptr>
+double student_t_qf(const T_p& p, const T_nu& nu, const T_mu& mu,
+                    const T_sigma& sigma) {
+  static constexpr const char* function = "student_t_qf";
+  check_nonnegative(function, "Degrees of freedom parameter", nu);
+  check_positive(function, "Scale parameter", sigma);
+  check_bounded(function, "Probability parameter", p, 0.0, 1.0);
+
+  if (p == 0.0) {
+    return NEGATIVE_INFTY;
+  } else if (p == 1.0) {
+    return INFTY;
+  } else if (p == 0.5) {
+    return mu;
+  }
+
+  const double p_val_flip = p < 0.5 ? p : 1.0 - p;
+  const double p_sign = p < 0.5 ? -1.0 : 1.0;
+  const auto& ibeta_arg = inv_inc_beta(0.5 * nu, 0.5, 2 * p_val_flip);
+
+  return mu + p_sign * sigma * sqrt(nu) * sqrt(-1.0 + 1.0 / ibeta_arg);
+}
+
+/**
+ * A vectorized version of the Student's t quantile function that accepts
+ * std::vectors, Eigen Matrix/Array objects, or expressions, and containers of
+ * these.
+ *
+ * @tparam T_p type of the probability parameter
+ * @tparam T_nu type of the degrees of freedom parameter
+ * @tparam T_mu type of the location parameter
+ * @tparam T_sigma type of the scale parameter
+ * @tparam T_container type of the container to hold results
+ * @param p Probability in the range [0, 1].
+ * @param nu Degrees of freedom, must be non-negative.
+ * @param mu Location parameter.
+ * @param sigma Scale parameter, must be positive.
+ * @return Container with quantile function values for each input.
+ */
+template <typename T_p, typename T_nu, typename T_mu, typename T_sigma,
+          typename T_container = common_container_t<T_p, T_nu, T_mu, T_sigma>,
+          require_any_vector_t<T_p, T_nu, T_mu, T_sigma>* = nullptr,
+          require_not_t<std::is_void<T_container>>* = nullptr>
+auto student_t_qf(const T_p& p, const T_nu& nu, const T_mu& mu,
+                  const T_sigma& sigma) {
+  static constexpr const char* function = "student_t_qf";
+  const size_t max_size_all = max_size(p, nu, mu, sigma);
+  T_container result(max_size_all);
+
+  ref_type_t<T_p> p_ref = p;
+  ref_type_t<T_nu> nu_ref = nu;
+  ref_type_t<T_mu> mu_ref = mu;
+  ref_type_t<T_sigma> sigma_ref = sigma;
+
+  scalar_seq_view<ref_type_t<T_p>> p_vec(p_ref);
+  scalar_seq_view<ref_type_t<T_nu>> nu_vec(nu_ref);
+  scalar_seq_view<ref_type_t<T_mu>> mu_vec(mu_ref);
+  scalar_seq_view<ref_type_t<T_sigma>> sigma_vec(sigma_ref);
+
+  for (size_t i = 0; i < max_size_all; ++i) {
+    result[i] = student_t_qf(p_vec[i], nu_vec[i], mu_vec[i], sigma_vec[i]);
+  }
+
+  return result;
+}
+
+}  // namespace math
+}  // namespace stan
+
+#endif
diff --git a/stan/math/rev/prob.hpp b/stan/math/rev/prob.hpp
index 98a6fee613f..f1f2427fdd3 100644
--- a/stan/math/rev/prob.hpp
+++ b/stan/math/rev/prob.hpp
@@ -4,5 +4,6 @@
 #include <stan/math/prim/fun/Eigen.hpp>
 
 #include <stan/math/rev/prob/std_normal_log_qf.hpp>
+#include <stan/math/rev/prob/student_t_qf.hpp>
 
 #endif
diff --git a/stan/math/rev/prob/student_t_qf.hpp b/stan/math/rev/prob/student_t_qf.hpp
new file mode 100644
index 00000000000..baaae8377e2
--- /dev/null
+++ b/stan/math/rev/prob/student_t_qf.hpp
@@ -0,0 +1,104 @@
+#ifndef STAN_MATH_REV_PROB_STUDENT_T_QF_HPP
+#define STAN_MATH_REV_PROB_STUDENT_T_QF_HPP
+
+#include <stan/math/rev/meta.hpp>
+#include <stan/math/rev/fun/sqrt.hpp>
+#include <stan/math/rev/fun/inv_inc_beta.hpp>
+#include <stan/math/prim/prob/student_t_lpdf.hpp>
+
+namespace stan {
+namespace math {
+
+template <typename T_p, typename T_nu, typename T_mu, typename T_sigma,
+          require_all_stan_scalar_t<T_p, T_mu, T_sigma, T_nu>* = nullptr,
+          require_any_var_t<T_p, T_mu, T_sigma, T_nu>* = nullptr>
+var student_t_qf(const T_p& p, const T_nu& nu, const T_mu& mu,
+                 const T_sigma& sigma) {
+  static constexpr const char* function = "student_t_qf";
+  const double p_val = value_of(p);
+  const double nu_val = value_of(nu);
+  const double mu_val = value_of(mu);
+  const double sigma_val = value_of(sigma);
+
+  check_nonnegative(function, "Degrees of freedom parameter", nu_val);
+  check_positive(function, "Scale parameter", sigma_val);
+  check_bounded(function, "Probability parameter", p_val, 0.0, 1.0);
+
+  double ret_val(0);
+  double ibeta_arg(0);
+
+  if (p == 0.0) {
+    ret_val = NEGATIVE_INFTY;
+  } else if (p == 1.0) {
+    ret_val = INFTY;
+  } else if (p == 0.5) {
+    ret_val = mu_val;
+  } else {
+    const double sqrt_inv_ibeta_m1 = sqrt(inv(ibeta_arg) - 1.0);
+    const double p_val_flip = p_val < 0.5 ? p_val : 1.0 - p_val;
+    const double p_sign = p_val < 0.5 ? -1.0 : 1.0;
+    ibeta_arg = inv_inc_beta(0.5 * nu_val, 0.5, 2 * p_val_flip);
+
+    ret_val = mu_val
+              + p_sign * sigma_val * math::sqrt(nu_val)
+                    * math::sqrt(-1.0 + 1.0 / ibeta_arg);
+  }
+
+  return make_callback_var(ret_val, [p, mu, sigma, nu,
+                                     ibeta_arg](auto& vi) mutable {
+    const double p_val = value_of(p);
+    const double mu_val = value_of(mu);
+    const double sigma_val = value_of(sigma);
+    const double nu_val = value_of(nu);
+    const double p_val_flip = p_val < 0.5 ? p_val : 1.0 - p_val;
+    const double p_sign = p_val < 0.5 ? -1.0 : 1.0;
+    const double sqrt_inv_ibeta_m1 = sqrt(inv(ibeta_arg) - 1.0);
+
+    if constexpr (!is_constant<T_p>::value) {
+      if (p.val() > 0.0 && p.val() < 1.0) {
+        p.adj() += vi.adj()
+                   * exp(-student_t_lpdf(vi.val(), nu_val, mu_val, sigma_val));
+      }
+    }
+
+    if constexpr (!is_constant<T_nu>::value) {
+      const double half_nu = nu_val / 2.0;
+      const double nu_p1_div_2 = (1 + nu_val) / 2.0;
+      Eigen::VectorXd hyper_arg_a(3);
+      hyper_arg_a << 0.5, half_nu, half_nu;
+      Eigen::VectorXd hyper_arg_b(2);
+      hyper_arg_b << 1 + half_nu, 1 + half_nu;
+      const double hyper_arg
+          = hypergeometric_pFq(hyper_arg_a, hyper_arg_b, ibeta_arg);
+
+      const double hyper2f1 = hypergeometric_2F1(
+          1.0, nu_p1_div_2, (2.0 + nu_val) / 2.0, ibeta_arg);
+      const double digamma_a1 = digamma(half_nu);
+      const double digamma_a2 = digamma(nu_p1_div_2);
+
+      const double arg_1 = (4 * hyper_arg * sqrt(1 - ibeta_arg)) / nu_val;
+      const double arg_2 = -2 * hyper2f1 * (-1 + ibeta_arg)
+                           * (log(ibeta_arg) - digamma_a1 + digamma_a2);
+
+      const double num1 = sigma_val * (-2 + (2 - arg_1 + arg_2) / ibeta_arg);
+      const double den1 = 4 * sqrt(nu_val) * sqrt_inv_ibeta_m1;
+      nu.adj() += vi.adj() * p_sign * num1 / den1;
+    }
+
+    if constexpr (!is_constant<T_mu>::value) {
+      if (p_val > 0.0 && p_val < 1.0) {
+        mu.adj() += vi.adj();
+      }
+    }
+
+    if constexpr (!is_constant<T_sigma>::value) {
+      if (p_val > 0.0 && p_val < 1.0) {
+        sigma.adj() += vi.adj() * p_sign * sqrt(nu_val) * sqrt_inv_ibeta_m1;
+      }
+    }
+  });
+}
+}  // namespace math
+}  // namespace stan
+
+#endif
diff --git a/test/unit/math/mix/prob/student_t_qf_test.cpp b/test/unit/math/mix/prob/student_t_qf_test.cpp
new file mode 100644
index 00000000000..1e5dc213b3c
--- /dev/null
+++ b/test/unit/math/mix/prob/student_t_qf_test.cpp
@@ -0,0 +1,29 @@
+#include <stan/math/rev.hpp>
+#include <stan/math/prim.hpp>
+#include <test/unit/math/rev/util.hpp>
+#include <test/unit/math/test_ad.hpp>
+
+TEST_F(AgradRev, mathMixProb_student_t_qf) {
+  auto f = [](const auto& p, const auto& nu, const auto& sigma) {
+    const auto& p_scaled = stan::math::inv_logit(p);
+    const auto& nu_scaled = stan::math::exp(nu);
+    const auto& sigma_scaled = stan::math::exp(sigma);
+    return stan::math::sum(
+        stan::math::student_t_qf(p_scaled, nu_scaled, 0, sigma_scaled));
+  };
+
+  stan::test::expect_ad(f, 0.3, 0.5, 3);
+  stan::test::expect_ad(f, 0.8, 0.5, 0.1);
+  stan::test::expect_ad(f, 0.1, 3, 3);
+
+  Eigen::VectorXd p(3);
+  p << 0.3, 0.8, 0.1;
+
+  Eigen::VectorXd nu(3);
+  nu << 0.5, 0.5, 3;
+
+  Eigen::VectorXd sigma(3);
+  sigma << 3, 0.1, 3;
+
+  stan::test::expect_ad(f, p, nu, sigma);
+}
diff --git a/test/unit/math/prim/prob/student_t_qf_test.cpp b/test/unit/math/prim/prob/student_t_qf_test.cpp
new file mode 100644
index 00000000000..67009d0c771
--- /dev/null
+++ b/test/unit/math/prim/prob/student_t_qf_test.cpp
@@ -0,0 +1,56 @@
+#include <stan/math/prim.hpp>
+#include <test/unit/util.hpp>
+#include <gtest/gtest.h>
+#include <limits>
+
+TEST(MathFunctions, student_t_qf_vals) {
+  using stan::math::student_t_qf;
+
+  Eigen::VectorXd p(5);
+  p << 0.1, 0.2, 0.5, 0.8, 0.9;
+
+  Eigen::VectorXd res = student_t_qf(p, 4, 2, 3);
+
+  EXPECT_FLOAT_EQ(student_t_qf(0.2, 4, 2, 3), -0.8228937317);
+  EXPECT_FLOAT_EQ(student_t_qf(0.8, 4, 2, 3), 4.822893732);
+  EXPECT_FLOAT_EQ(student_t_qf(0.8, 3, 0, 1), 0.9784723124);
+
+  EXPECT_FLOAT_EQ(student_t_qf(0.5, 3, 0, 1), 0.0);
+}
+
+TEST(MathFunctions, student_t_qf_vec) {
+  using stan::math::student_t_qf;
+
+  Eigen::VectorXd p(5);
+  p << 0.1, 0.2, 0.5, 0.8, 0.9;
+
+  Eigen::VectorXd nu(5);
+  nu << 4, 4, 3, 3, 3;
+
+  Eigen::VectorXd res(5);
+  for (int i = 0; i < 5; ++i) {
+    res[i] = student_t_qf(p[i], nu[i], 2, 3);
+  }
+  EXPECT_MATRIX_FLOAT_EQ(res, student_t_qf(p, nu, 2, 3));
+
+  Eigen::RowVectorXd p_rowvec = p.transpose();
+  Eigen::RowVectorXd nu_rowvec = nu.transpose();
+  Eigen::RowVectorXd res_rowvec = res.transpose();
+
+  EXPECT_MATRIX_FLOAT_EQ(res_rowvec, student_t_qf(p_rowvec, nu_rowvec, 2, 3));
+}
+
+TEST(MathFunctions, student_t_qf_inf) {
+  using stan::math::student_t_qf;
+  long double log_p = std::numeric_limits<long double>::min();
+  const double inf = std::numeric_limits<double>::infinity();
+  EXPECT_EQ(student_t_qf(0, 2, 0, 1), -inf);
+  EXPECT_EQ(student_t_qf(1, 2, 0, 1), inf);
+}
+
+TEST(MathFunctions, student_t_qf_nan) {
+  using stan::math::student_t_qf;
+  double nan = std::numeric_limits<double>::quiet_NaN();
+  EXPECT_THROW(student_t_qf(nan, 5, 0, 1), std::domain_error);
+  EXPECT_THROW(student_t_qf(2.1, 5, 0, 1), std::domain_error);
+}