diff --git a/stan/math/prim/prob.hpp b/stan/math/prim/prob.hpp
index 04416a966ec..2fce0b34e1d 100644
--- a/stan/math/prim/prob.hpp
+++ b/stan/math/prim/prob.hpp
@@ -175,7 +175,9 @@
 #include <stan/math/prim/prob/logistic_log.hpp>
 #include <stan/math/prim/prob/logistic_lpdf.hpp>
 #include <stan/math/prim/prob/logistic_rng.hpp>
+#include <stan/math/prim/prob/loglogistic_cdf.hpp>
 #include <stan/math/prim/prob/loglogistic_lpdf.hpp>
+#include <stan/math/prim/prob/loglogistic_rng.hpp>
 #include <stan/math/prim/prob/lognormal_ccdf_log.hpp>
 #include <stan/math/prim/prob/lognormal_cdf.hpp>
 #include <stan/math/prim/prob/lognormal_cdf_log.hpp>
diff --git a/stan/math/prim/prob/loglogistic_cdf.hpp b/stan/math/prim/prob/loglogistic_cdf.hpp
new file mode 100644
index 00000000000..33d4ebb161c
--- /dev/null
+++ b/stan/math/prim/prob/loglogistic_cdf.hpp
@@ -0,0 +1,126 @@
+#ifndef STAN_MATH_PRIM_PROB_LOGLOGISTIC_CDF_HPP
+#define STAN_MATH_PRIM_PROB_LOGLOGISTIC_CDF_HPP
+
+#include <stan/math/prim/meta.hpp>
+#include <stan/math/prim/err.hpp>
+#include <stan/math/prim/fun/as_column_vector_or_scalar.hpp>
+#include <stan/math/prim/fun/as_array_or_scalar.hpp>
+#include <stan/math/prim/fun/as_value_column_array_or_scalar.hpp>
+#include <stan/math/prim/fun/exp.hpp>
+#include <stan/math/prim/fun/log.hpp>
+#include <stan/math/prim/fun/max_size.hpp>
+#include <stan/math/prim/fun/size.hpp>
+#include <stan/math/prim/fun/size_zero.hpp>
+#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/fun/promote_scalar.hpp>
+#include <cmath>
+#include <iostream>
+
+namespace stan {
+namespace math {
+
+/** \ingroup prob_dists
+ * The loglogistic cumulative distribution function for the specified
+ * scalar(s) given the specified scales(s) and shape(s). y, alpha, or
+ * beta can each be either a scalar or a vector. Any vector inputs
+ * must be the same length.
+ *
+ *
+ * @tparam T_y type of scalar.
+ * @tparam T_scale type of scale parameter.
+ * @tparam T_shape type of shape parameter.
+ * @param y (Sequence of) scalar(s).
+ * @param alpha (Sequence of) scale parameter(s)
+ * for the loglogistic distribution.
+ * @param beta (Sequence of) shape parameter(s) for the
+ * loglogistic distribution.
+ * @return The loglogistic cdf evaluated at the specified arguments.
+ * @throw std::domain_error if any of the inputs are not positive or
+ * if and of the parameters are not finite.
+ */
+template <typename T_y, typename T_scale, typename T_shape,
+          require_all_not_nonscalar_prim_or_rev_kernel_expression_t<
+              T_y, T_scale, T_shape>* = nullptr>
+return_type_t<T_y, T_scale, T_shape> loglogistic_cdf(const T_y& y,
+                                                     const T_scale& alpha,
+                                                     const T_shape& beta) {
+  using T_partials_return = partials_return_t<T_y, T_scale, T_shape>;
+  using T_y_ref = ref_type_t<T_y>;
+  using T_alpha_ref = ref_type_t<T_scale>;
+  using T_beta_ref = ref_type_t<T_shape>;
+  using std::pow;
+  static const char* function = "loglogistic_cdf";
+  check_consistent_sizes(function, "Random variable", y, "Scale parameter",
+                         alpha, "Shape parameter", beta);
+  T_y_ref y_ref = y;
+  T_alpha_ref alpha_ref = alpha;
+  T_beta_ref beta_ref = beta;
+
+  decltype(auto) y_val = to_ref(as_value_column_array_or_scalar(y_ref));
+  decltype(auto) alpha_val = to_ref(as_value_column_array_or_scalar(alpha_ref));
+  decltype(auto) beta_val = to_ref(as_value_column_array_or_scalar(beta_ref));
+
+  check_nonnegative(function, "Random variable", y_val);
+  check_positive_finite(function, "Scale parameter", alpha_val);
+  check_positive_finite(function, "Shape parameter", beta_val);
+
+  if (size_zero(y, alpha, beta)) {
+    return 1.0;
+  }
+
+  operands_and_partials<T_y_ref, T_alpha_ref, T_beta_ref> ops_partials(
+      y_ref, alpha_ref, beta_ref);
+
+  if (sum(promote_scalar<int>(y_val == 0))) {
+    return ops_partials.build(0.0);
+  }
+
+  const auto& alpha_div_y
+      = to_ref_if<!is_constant_all<T_shape>::value>(alpha_val / y_val);
+  const auto& alpha_div_y_pow_beta
+      = to_ref_if<!is_constant_all<T_y, T_scale, T_shape>::value>(
+          pow(alpha_div_y, beta_val));
+  const auto& prod_all
+      = to_ref_if<!is_constant_all<T_y, T_scale, T_shape>::value>(
+          1 / (1 + alpha_div_y_pow_beta));
+
+  T_partials_return cdf = prod(prod_all);
+
+  if (!is_constant_all<T_y, T_scale, T_shape>::value) {
+    const auto& prod_all_sq = to_ref_if<!is_constant_all<T_y>::value
+                                            + !is_constant_all<T_scale>::value
+                                            + !is_constant_all<T_shape>::value
+                                        >= 2>(square(prod_all));
+    const auto& cdf_div_elt = to_ref_if<!is_constant_all<T_y>::value
+                                            + !is_constant_all<T_scale>::value
+                                            + !is_constant_all<T_shape>::value
+                                        >= 2>(cdf / prod_all);
+    if (!is_constant_all<T_y, T_scale>::value) {
+      const auto& alpha_div_times_beta = to_ref_if<
+          !is_constant_all<T_y>::value + !is_constant_all<T_scale>::value == 2>(
+          alpha_div_y_pow_beta * beta_val);
+      if (!is_constant_all<T_y>::value) {
+        const auto& y_deriv = alpha_div_times_beta / y_val * prod_all_sq;
+        ops_partials.edge1_.partials_ = y_deriv * cdf_div_elt;
+      }
+      if (!is_constant_all<T_scale>::value) {
+        const auto& alpha_deriv
+            = -alpha_div_times_beta / alpha_val * prod_all_sq;
+        ops_partials.edge2_.partials_ = alpha_deriv * cdf_div_elt;
+      }
+    }
+    if (!is_constant_all<T_shape>::value) {
+      const auto& beta_deriv
+          = -multiply_log(alpha_div_y_pow_beta, alpha_div_y) * prod_all_sq;
+      ops_partials.edge3_.partials_ = beta_deriv * cdf_div_elt;
+    }
+  }
+
+  return ops_partials.build(cdf);
+}
+
+}  // namespace math
+}  // namespace stan
+#endif
diff --git a/stan/math/prim/prob/loglogistic_rng.hpp b/stan/math/prim/prob/loglogistic_rng.hpp
new file mode 100644
index 00000000000..06b43e7e3a1
--- /dev/null
+++ b/stan/math/prim/prob/loglogistic_rng.hpp
@@ -0,0 +1,63 @@
+#ifndef STAN_MATH_PRIM_PROB_LOGLOGISTIC_RNG_HPP
+#define STAN_MATH_PRIM_PROB_LOGLOGISTIC_RNG_HPP
+
+#include <stan/math/prim/meta.hpp>
+#include <stan/math/prim/err.hpp>
+#include <stan/math/prim/fun/scalar_seq_view.hpp>
+#include <stan/math/prim/fun/max_size.hpp>
+#include <boost/random/uniform_01.hpp>
+#include <boost/random/variate_generator.hpp>
+
+namespace stan {
+namespace math {
+
+/** \ingroup prob_dists
+ * Return a loglogistic random variate for the given scale and
+ * shape parameters using the specified random number generator.
+ *
+ * alpha and beta can each be a scalar or a one-dimensional container. Any
+ * non-scalar inputs must be the same size.
+ *
+ * @tparam T_scale type of scale parameter
+ * @tparam T_shape type of shape parameter
+ * @tparam RNG type of random number generator
+ * @param alpha (Sequence of) positive scale parameter(s)
+ * @param beta (Sequence of) positive shape parameter(s)
+ * @param rng random number generator
+ * @return (Sequence of) loglogistic random variate(s)
+ * @throw std::domain_error if alpha or beta are nonpositive
+ * @throw std::invalid_argument if non-scalar arguments are of different
+ * sizes
+ */
+template <typename T_scale, typename T_shape, class RNG>
+inline typename VectorBuilder<true, double, T_scale, T_shape>::type
+loglogistic_rng(const T_scale& alpha, const T_shape& beta, RNG& rng) {
+  using boost::uniform_01;
+  using boost::variate_generator;
+  using std::pow;
+  using T_alpha_ref = ref_type_t<T_scale>;
+  using T_beta_ref = ref_type_t<T_shape>;
+  static const char* function = "loglogistic_rng";
+  check_consistent_sizes(function, "Scale parameter", alpha, "Shape Parameter",
+                         beta);
+  T_alpha_ref alpha_ref = alpha;
+  T_beta_ref beta_ref = beta;
+  check_positive_finite(function, "Scale parameter", alpha_ref);
+  check_positive_finite(function, "Shape parameter", beta_ref);
+
+  scalar_seq_view<T_alpha_ref> alpha_vec(alpha_ref);
+  scalar_seq_view<T_beta_ref> beta_vec(beta_ref);
+  size_t N = max_size(alpha, beta);
+  VectorBuilder<true, double, T_scale, T_shape> output(N);
+
+  for (size_t n = 0; n < N; ++n) {
+    variate_generator<RNG&, uniform_01<> > uniform01_rng(rng, uniform_01<>());
+    const double tmp = uniform01_rng();
+    output[n] = alpha_vec[n] * pow(tmp / (1 - tmp), 1 / beta_vec[n]);
+  }
+  return output.data();
+}
+
+}  // namespace math
+}  // namespace stan
+#endif
diff --git a/test/prob/loglogistic/loglogistic_cdf_test.hpp b/test/prob/loglogistic/loglogistic_cdf_test.hpp
new file mode 100644
index 00000000000..4f36f6316eb
--- /dev/null
+++ b/test/prob/loglogistic/loglogistic_cdf_test.hpp
@@ -0,0 +1,86 @@
+// Arguments: Doubles, Doubles, Doubles
+#include <stan/math/prim.hpp>
+
+using stan::math::var;
+using std::numeric_limits;
+using std::pow;
+using std::vector;
+
+class AgradCdfLogistic : public AgradCdfTest {
+ public:
+  void valid_values(vector<vector<double> >& parameters, vector<double>& cdf) {
+    vector<double> param(3);
+
+    param[0] = 2.0;  // y
+    param[1] = 3.0;  // Scale
+    param[2] = 2.0;  // Shape
+    parameters.push_back(param);
+    cdf.push_back(0.30769230769230765388);  // expected cdf
+
+    param[0] = 20.0;  // y
+    param[1] = 11.5;  // Scale
+    param[2] = 15.0;  // Shape
+    parameters.push_back(param);
+    cdf.push_back(0.99975173823523311167);  // expected cdf
+
+    param[0] = 0.0001;  // y
+    param[1] = 1.2;     // Scale
+    param[2] = 1.0;     // Shape
+    parameters.push_back(param);
+    cdf.push_back(0.00008332638946754438);  // expected cdf
+  }
+
+  void invalid_values(vector<size_t>& index, vector<double>& value) {
+    // y
+    index.push_back(0U);
+    value.push_back(-1.0);
+
+    index.push_back(0U);
+    value.push_back(-numeric_limits<double>::infinity());
+
+    // alpha
+    index.push_back(1U);
+    value.push_back(-1);
+
+    index.push_back(1U);
+    value.push_back(0);
+
+    index.push_back(1U);
+    value.push_back(numeric_limits<double>::infinity());
+
+    // beta
+    index.push_back(2U);
+    value.push_back(-1.0);
+
+    index.push_back(2U);
+    value.push_back(0.0);
+
+    index.push_back(2U);
+    value.push_back(numeric_limits<double>::infinity());
+  }
+
+  bool has_lower_bound() { return true; }
+
+  double lower_bound() { return 0.0; }
+
+  bool has_upper_bound() { return false; }
+
+  template <typename T_y, typename T_scale, typename T_shape, typename T3,
+            typename T4, typename T5>
+  stan::return_type_t<T_y, T_scale, T_shape> cdf(const T_y& y,
+                                                 const T_scale& alpha,
+                                                 const T_shape& beta, const T3&,
+                                                 const T4&, const T5&) {
+    return stan::math::loglogistic_cdf(y, alpha, beta);
+  }
+
+  template <typename T_y, typename T_scale, typename T_shape, typename T3,
+            typename T4, typename T5>
+  stan::return_type_t<T_y, T_scale, T_shape> cdf_function(const T_y& y,
+                                                          const T_scale& alpha,
+                                                          const T_shape& beta,
+                                                          const T3&, const T4&,
+                                                          const T5&) {
+    return 1.0 / (1.0 + pow(alpha / y, beta));
+  }
+};
diff --git a/test/unit/math/prim/prob/loglogistic_test.cpp b/test/unit/math/prim/prob/loglogistic_test.cpp
new file mode 100644
index 00000000000..94a38ea4cc7
--- /dev/null
+++ b/test/unit/math/prim/prob/loglogistic_test.cpp
@@ -0,0 +1,77 @@
+#include <stan/math/prim.hpp>
+#include <test/unit/math/prim/prob/vector_rng_test_helper.hpp>
+#include <test/unit/math/prim/prob/util.hpp>
+#include <gtest/gtest.h>
+#include <boost/math/distributions.hpp>
+#include <boost/random/mersenne_twister.hpp>
+#include <limits>
+#include <vector>
+
+class LoglogisticTestRig : public VectorRNGTestRig {
+ public:
+  LoglogisticTestRig()
+      : VectorRNGTestRig(10000, 10, {2.5, 1.7, 0.2, 0.1, 2.0},
+                         {3, 2, 1, 5, 10, 6}, {-2.5, -1.7, -0.2, -0.1, 0.0},
+                         {-3, -2, -1, -4, -10, 0}, {0.1, 1.0, 2.5, 4.0},
+                         {1, 2, 3, 4}, {-2.7, -1.5, -0.5, 0.0},
+                         {-3, -2, -1, 0}) {}
+
+  template <typename T1, typename T2, typename T3, typename T_rng>
+  auto generate_samples(const T1& alpha, const T2& beta, const T3& unused,
+                        T_rng& rng) const {
+    return stan::math::loglogistic_rng(alpha, beta, rng);
+  }
+};
+
+double icdf(double x, double alpha, double beta) {
+  return alpha * pow(x / (1 - x), 1 / beta);
+}
+
+TEST(ProbDistributionsLoglogistic, errorCheck) {
+  check_dist_throws_all_types(LoglogisticTestRig());
+}
+
+TEST(ProbDistributionsLoglogistic, error_check) {
+  boost::random::mt19937 rng;
+  EXPECT_NO_THROW(stan::math::loglogistic_rng(10.0, 2.0, rng));
+
+  EXPECT_THROW(stan::math::loglogistic_rng(2.0, -1.0, rng), std::domain_error);
+  EXPECT_THROW(stan::math::loglogistic_rng(-2.0, 1.0, rng), std::domain_error);
+  EXPECT_THROW(
+      stan::math::loglogistic_rng(10, stan::math::positive_infinity(), rng),
+      std::domain_error);
+  EXPECT_THROW(
+      stan::math::loglogistic_rng(stan::math::positive_infinity(), 2, rng),
+      std::domain_error);
+}
+
+TEST(ProbDistributionsLoglogistic, test_sampling_icdf) {
+  for (double p : {0.0, 0.1, 0.2, 0.5, 0.7, 0.9, 0.99}) {
+    for (double alpha : {1.11, 0.13, 1.2, 4.67}) {
+      for (double beta : {0.11, 1.33, 2.0, 3.2}) {
+        double x = icdf(p, alpha, beta);
+        EXPECT_FLOAT_EQ(stan::math::loglogistic_cdf(x, alpha, beta), p);
+      }
+    }
+  }
+}
+
+TEST(ProbDistributionsLoglogistic, chiSquareGoodnessFitTest) {
+  boost::random::mt19937 rng;
+  int N = 10000;
+  int K = stan::math::round(2 * std::pow(N, 0.4));
+
+  std::vector<double> samples;
+  for (int i = 0; i < N; ++i) {
+    samples.push_back(stan::math::loglogistic_rng(1.2, 2.0, rng));
+  }
+  std::vector<double> quantiles;
+  for (int i = 1; i < K; ++i) {
+    double frac = static_cast<double>(i) / K;
+    quantiles.push_back(icdf(frac, 1.2, 2.0));
+  }
+  quantiles.push_back(std::numeric_limits<double>::max());
+
+  // Assert that they match
+  assert_matches_quantiles(samples, quantiles, 1e-6);
+}