diff --git a/test/unit/math/fwd/prob/categorical_test.cpp b/test/unit/math/fwd/prob/categorical_test.cpp
index 99562998d15..e5e0f6cbfa5 100644
--- a/test/unit/math/fwd/prob/categorical_test.cpp
+++ b/test/unit/math/fwd/prob/categorical_test.cpp
@@ -1,6 +1,6 @@
 #include <stan/math/fwd.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <boost/math/distributions.hpp>
 #include <vector>
 
diff --git a/test/unit/math/fwd/prob/dirichlet_test.cpp b/test/unit/math/fwd/prob/dirichlet_test.cpp
index 8edb32581a3..c535728cde2 100644
--- a/test/unit/math/fwd/prob/dirichlet_test.cpp
+++ b/test/unit/math/fwd/prob/dirichlet_test.cpp
@@ -1,6 +1,6 @@
 #include <stan/math/fwd.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <boost/math/distributions.hpp>
 
 TEST(ProbDistributions, fvar_double) {
diff --git a/test/unit/math/fwd/prob/inv_wishart_test.cpp b/test/unit/math/fwd/prob/inv_wishart_test.cpp
index d4b7889d76e..3b32d246785 100644
--- a/test/unit/math/fwd/prob/inv_wishart_test.cpp
+++ b/test/unit/math/fwd/prob/inv_wishart_test.cpp
@@ -1,7 +1,7 @@
 #include <stan/math/fwd.hpp>
 #include <gtest/gtest.h>
 
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <boost/math/distributions.hpp>
 #include <boost/math/special_functions/digamma.hpp>
 
diff --git a/test/unit/math/fwd/prob/lkj_corr_test.cpp b/test/unit/math/fwd/prob/lkj_corr_test.cpp
index 6063bedc0c6..d0fff6b1bbf 100644
--- a/test/unit/math/fwd/prob/lkj_corr_test.cpp
+++ b/test/unit/math/fwd/prob/lkj_corr_test.cpp
@@ -1,11 +1,11 @@
 #include <stan/math/fwd.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <boost/math/distributions.hpp>
 
 TEST(ProbDistributionsLkjCorr, fvar_double) {
   using stan::math::fvar;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int K = 4;
   Eigen::Matrix<fvar<double>, Eigen::Dynamic, Eigen::Dynamic> Sigma(K, K);
   Sigma.setZero();
@@ -24,7 +24,7 @@ TEST(ProbDistributionsLkjCorr, fvar_double) {
 
 TEST(ProbDistributionsLkjCorrCholesky, fvar_double) {
   using stan::math::fvar;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int K = 4;
   Eigen::Matrix<fvar<double>, Eigen::Dynamic, Eigen::Dynamic> Sigma(K, K);
   Sigma.setZero();
@@ -43,7 +43,7 @@ TEST(ProbDistributionsLkjCorrCholesky, fvar_double) {
 
 TEST(ProbDistributionsLkjCorr, fvar_fvar_double) {
   using stan::math::fvar;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int K = 4;
   Eigen::Matrix<fvar<fvar<double> >, Eigen::Dynamic, Eigen::Dynamic> Sigma(K,
                                                                            K);
@@ -63,7 +63,7 @@ TEST(ProbDistributionsLkjCorr, fvar_fvar_double) {
 
 TEST(ProbDistributionsLkjCorrCholesky, fvar_fvar_double) {
   using stan::math::fvar;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int K = 4;
   Eigen::Matrix<fvar<fvar<double> >, Eigen::Dynamic, Eigen::Dynamic> Sigma(K,
                                                                            K);
diff --git a/test/unit/math/fwd/prob/multi_normal_cholesky_test.cpp b/test/unit/math/fwd/prob/multi_normal_cholesky_test.cpp
index af2767268c8..b3f54c2177f 100644
--- a/test/unit/math/fwd/prob/multi_normal_cholesky_test.cpp
+++ b/test/unit/math/fwd/prob/multi_normal_cholesky_test.cpp
@@ -1,6 +1,6 @@
 #include <stan/math/fwd.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <boost/math/distributions.hpp>
 
 TEST(ProbDistributionsMultiNormalCholesky, fvar_double) {
diff --git a/test/unit/math/fwd/prob/multi_normal_test.cpp b/test/unit/math/fwd/prob/multi_normal_test.cpp
index 15ce1d2a565..d5f8c169a69 100644
--- a/test/unit/math/fwd/prob/multi_normal_test.cpp
+++ b/test/unit/math/fwd/prob/multi_normal_test.cpp
@@ -1,6 +1,6 @@
 #include <stan/math/fwd.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <boost/math/distributions.hpp>
 
 TEST(ProbDistributionsMultiNormal, fvar_double) {
diff --git a/test/unit/math/fwd/prob/multi_student_t_test.cpp b/test/unit/math/fwd/prob/multi_student_t_test.cpp
index e0360b2e251..c6d975f1e85 100644
--- a/test/unit/math/fwd/prob/multi_student_t_test.cpp
+++ b/test/unit/math/fwd/prob/multi_student_t_test.cpp
@@ -1,6 +1,6 @@
 #include <stan/math/fwd.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <boost/math/distributions.hpp>
 
 TEST(ProbDistributionsMultiStudentT, fvar_double) {
diff --git a/test/unit/math/fwd/prob/multinomial_test.cpp b/test/unit/math/fwd/prob/multinomial_test.cpp
index bb5a975f576..ef3cd316722 100644
--- a/test/unit/math/fwd/prob/multinomial_test.cpp
+++ b/test/unit/math/fwd/prob/multinomial_test.cpp
@@ -1,6 +1,6 @@
 #include <stan/math/fwd.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <boost/math/distributions.hpp>
 #include <vector>
 
diff --git a/test/unit/math/fwd/prob/wishart_test.cpp b/test/unit/math/fwd/prob/wishart_test.cpp
index b0dc20e74ce..ccea981fbd5 100644
--- a/test/unit/math/fwd/prob/wishart_test.cpp
+++ b/test/unit/math/fwd/prob/wishart_test.cpp
@@ -1,6 +1,6 @@
 #include <stan/math/fwd.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <boost/math/special_functions/digamma.hpp>
 #include <boost/math/distributions.hpp>
 
diff --git a/test/unit/math/mix/prob/categorical_test.cpp b/test/unit/math/mix/prob/categorical_test.cpp
index 01d24102956..920e9683c02 100644
--- a/test/unit/math/mix/prob/categorical_test.cpp
+++ b/test/unit/math/mix/prob/categorical_test.cpp
@@ -1,6 +1,6 @@
 #include <stan/math/mix.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <boost/math/distributions.hpp>
 #include <vector>
 
diff --git a/test/unit/math/mix/prob/lkj_corr_test.cpp b/test/unit/math/mix/prob/lkj_corr_test.cpp
index 64c65d419ea..69e0ff12cd3 100644
--- a/test/unit/math/mix/prob/lkj_corr_test.cpp
+++ b/test/unit/math/mix/prob/lkj_corr_test.cpp
@@ -1,7 +1,7 @@
 #include <stan/math/mix.hpp>
 #include <gtest/gtest.h>
 #include <test/unit/math/rev/prob/lkj_corr_cholesky_test_functors.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <boost/math/distributions.hpp>
 #include <test/unit/math/mix/prob/higher_order_utils.hpp>
 #include <vector>
@@ -9,7 +9,7 @@
 TEST(ProbDistributionsLkjCorr, fvar_var) {
   using stan::math::fvar;
   using stan::math::var;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int K = 4;
   Eigen::Matrix<fvar<var>, Eigen::Dynamic, Eigen::Dynamic> Sigma(K, K);
   Sigma.setZero();
@@ -31,7 +31,7 @@ TEST(ProbDistributionsLkjCorr, fvar_var) {
 TEST(ProbDistributionsLkjCorrCholesky, fvar_var) {
   using stan::math::fvar;
   using stan::math::var;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int K = 4;
   Eigen::Matrix<fvar<var>, Eigen::Dynamic, Eigen::Dynamic> Sigma(K, K);
   Sigma.setZero();
@@ -54,7 +54,7 @@ TEST(ProbDistributionsLkjCorrCholesky, fvar_var) {
 TEST(ProbDistributionsLkjCorr, fvar_fvar_var) {
   using stan::math::fvar;
   using stan::math::var;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int K = 4;
   Eigen::Matrix<fvar<fvar<var> >, Eigen::Dynamic, Eigen::Dynamic> Sigma(K, K);
   Sigma.setZero();
@@ -78,7 +78,7 @@ TEST(ProbDistributionsLkjCorr, fvar_fvar_var) {
 TEST(ProbDistributionsLkjCorrCholesky, fvar_fvar_var) {
   using stan::math::fvar;
   using stan::math::var;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int K = 4;
   Eigen::Matrix<fvar<fvar<var> >, Eigen::Dynamic, Eigen::Dynamic> Sigma(K, K);
   Sigma.setZero();
diff --git a/test/unit/math/mix/prob/multinomial_test.cpp b/test/unit/math/mix/prob/multinomial_test.cpp
index 9289d5d26f6..23b58a03108 100644
--- a/test/unit/math/mix/prob/multinomial_test.cpp
+++ b/test/unit/math/mix/prob/multinomial_test.cpp
@@ -1,6 +1,6 @@
 #include <stan/math/mix.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <boost/math/distributions.hpp>
 #include <vector>
 
diff --git a/test/unit/math/opencl/multiply_transpose_test.cpp b/test/unit/math/opencl/multiply_transpose_test.cpp
index 5f193aaacc1..b36e50bfa54 100644
--- a/test/unit/math/opencl/multiply_transpose_test.cpp
+++ b/test/unit/math/opencl/multiply_transpose_test.cpp
@@ -4,10 +4,10 @@
 #include <stan/math/opencl/multiply_transpose.hpp>
 #include <stan/math/opencl/copy.hpp>
 #include <test/unit/util.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <gtest/gtest.h>
 #include <algorithm>
-boost::random::mt19937 rng;
+boost::random::mixmax rng;
 
 TEST(MathMatrixOpenCL, multiply_transpose_exception_fail_zero) {
   stan::math::row_vector_d rv(0);
diff --git a/test/unit/math/opencl/prim/mdivide_left_tri_low_test.cpp b/test/unit/math/opencl/prim/mdivide_left_tri_low_test.cpp
index 89e66e4f1c0..310dff21f82 100644
--- a/test/unit/math/opencl/prim/mdivide_left_tri_low_test.cpp
+++ b/test/unit/math/opencl/prim/mdivide_left_tri_low_test.cpp
@@ -3,7 +3,7 @@
 #include <stan/math/opencl/prim.hpp>
 #include <test/unit/math/expect_near_rel.hpp>
 #include <test/unit/util.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <gtest/gtest.h>
 #include <algorithm>
 
@@ -29,7 +29,7 @@ TEST(MathMatrixCL, mdivide_left_tri_low_cl_exception) {
 }
 
 void mdivide_left_tri_low_Ab_test(int size) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   auto m1 = stan::math::matrix_d(size, size);
   for (int i = 0; i < size; i++) {
     for (int j = 0; j < i; j++) {
@@ -57,7 +57,7 @@ void mdivide_left_tri_low_Ab_test(int size) {
 }
 
 void mdivide_left_tri_low_A_test(int size) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   auto m1 = stan::math::matrix_d(size, size);
   for (int i = 0; i < size; i++) {
     for (int j = 0; j < i; j++) {
diff --git a/test/unit/math/opencl/prim/mdivide_right_tri_low_test.cpp b/test/unit/math/opencl/prim/mdivide_right_tri_low_test.cpp
index 9c4cc3ceb91..91958ef6879 100644
--- a/test/unit/math/opencl/prim/mdivide_right_tri_low_test.cpp
+++ b/test/unit/math/opencl/prim/mdivide_right_tri_low_test.cpp
@@ -3,7 +3,7 @@
 #include <stan/math/opencl/prim.hpp>
 #include <test/unit/math/expect_near_rel.hpp>
 #include <test/unit/util.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <gtest/gtest.h>
 #include <algorithm>
 
@@ -29,7 +29,7 @@ TEST(MathMatrixCL, mdivide_right_tri_low_cl_exception) {
 }
 
 void mdivide_right_tri_low_Ab_test(int size) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   auto m1 = stan::math::matrix_d(size, size);
   for (int i = 0; i < size; i++) {
     for (int j = 0; j < i; j++) {
diff --git a/test/unit/math/opencl/tri_inverse_test.cpp b/test/unit/math/opencl/tri_inverse_test.cpp
index 9eda4495e94..45f2d2d9810 100644
--- a/test/unit/math/opencl/tri_inverse_test.cpp
+++ b/test/unit/math/opencl/tri_inverse_test.cpp
@@ -2,7 +2,7 @@
 #include <stan/math/prim.hpp>
 #include <stan/math/opencl/prim.hpp>
 #include <test/unit/util.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <gtest/gtest.h>
 #include <algorithm>
 
@@ -25,7 +25,7 @@ TEST(MathMatrixCL, inverse_cl_exception) {
 }
 
 void lower_inverse_test(int size) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   auto m1 = stan::math::matrix_d(size, size);
   for (int i = 0; i < size; i++) {
     for (int j = 0; j < i; j++) {
@@ -56,7 +56,7 @@ void lower_inverse_test(int size) {
 }
 
 void upper_inverse_test(int size) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   auto m1 = stan::math::matrix_d(size, size);
   for (int i = 0; i < size; i++) {
     for (int j = 0; j < i; j++) {
diff --git a/test/unit/math/prim/fun/chol2inv_test.cpp b/test/unit/math/prim/fun/chol2inv_test.cpp
index 9942e366ac1..0a0ba7c48cd 100644
--- a/test/unit/math/prim/fun/chol2inv_test.cpp
+++ b/test/unit/math/prim/fun/chol2inv_test.cpp
@@ -1,5 +1,5 @@
 #include <stan/math/prim.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <gtest/gtest.h>
 
 TEST(MathMatrixPrimMat, chol2inv_exception) {
@@ -25,7 +25,7 @@ TEST(MathMatrixPrimMat, chol2inv) {
   using stan::math::matrix_d;
   using stan::math::wishart_rng;
 
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   matrix_d I(3, 3);
   I.setZero();
   I.diagonal().setOnes();
diff --git a/test/unit/math/prim/prob/bernoulli_logit_glm_rng_test.cpp b/test/unit/math/prim/prob/bernoulli_logit_glm_rng_test.cpp
index 35e800b3a63..9e1b302b099 100644
--- a/test/unit/math/prim/prob/bernoulli_logit_glm_rng_test.cpp
+++ b/test/unit/math/prim/prob/bernoulli_logit_glm_rng_test.cpp
@@ -1,7 +1,7 @@
 #include <gtest/gtest.h>
 #include <stan/math/prim.hpp>
 #include <boost/math/distributions.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <test/unit/math/prim/prob/util.hpp>
 #include <limits>
 #include <vector>
@@ -9,7 +9,7 @@
 TEST(ProbDistributionsBernoulliLogitGlm, vectorized) {
   using stan::math::bernoulli_logit_glm_rng;
   //  Test scalar/vector combinations.
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
 
   Eigen::MatrixXd x(2, 3);
   x << 3.5, -1.5, 0.0, 2.0, 1.0, 3.0;
@@ -70,7 +70,7 @@ TEST(ProbDistributionsBernoulliLogitGlm, vectorized) {
 TEST(ProbDistributionsBernoulliLogitGlm, errorCheck) {
   using stan::math::bernoulli_logit_glm_rng;
   // Check errors for nonfinite and wrong sizes.
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
 
   int N = 3;
   int M = 2;
@@ -107,7 +107,7 @@ TEST(ProbDistributionsBernoulliLogitGlm, errorCheck) {
 TEST(ProbDistributionsBernoulliLogitGlm, marginalChiSquareGoodnessFitTest) {
   using stan::math::bernoulli_logit_glm_rng;
   // Check distribution of result.
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   Eigen::MatrixXd x(2, 2);
   x << 3.5, -1.5, 2.0, -1.2;
   std::vector<double> alpha{2.0, 1.0};
diff --git a/test/unit/math/prim/prob/bernoulli_logit_test.cpp b/test/unit/math/prim/prob/bernoulli_logit_test.cpp
index e608decd807..ebe497cb002 100644
--- a/test/unit/math/prim/prob/bernoulli_logit_test.cpp
+++ b/test/unit/math/prim/prob/bernoulli_logit_test.cpp
@@ -3,7 +3,7 @@
 #include <test/unit/math/prim/prob/VectorIntRNGTestRig.hpp>
 #include <test/unit/math/prim/prob/util.hpp>
 #include <boost/math/distributions.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <gtest/gtest.h>
 #include <limits>
 #include <vector>
@@ -35,7 +35,7 @@ TEST(ProbDistributionsBernoulliLogit, distributionCheck) {
 }
 
 TEST(ProbDistributionsBernoulliLogit, error_check) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
 
   EXPECT_NO_THROW(stan::math::bernoulli_logit_rng(-3.5, rng));
   EXPECT_THROW(
@@ -44,7 +44,7 @@ TEST(ProbDistributionsBernoulliLogit, error_check) {
 }
 
 TEST(ProbDistributionsBernoulliLogit, logitChiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   // number of samples
   int N = 10000;
 
diff --git a/test/unit/math/prim/prob/bernoulli_test.cpp b/test/unit/math/prim/prob/bernoulli_test.cpp
index a1e3b722d7a..cd94a3d34b8 100644
--- a/test/unit/math/prim/prob/bernoulli_test.cpp
+++ b/test/unit/math/prim/prob/bernoulli_test.cpp
@@ -3,7 +3,7 @@
 #include <test/unit/math/prim/prob/VectorIntRNGTestRig.hpp>
 #include <test/unit/math/prim/prob/util.hpp>
 #include <boost/math/distributions.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <gtest/gtest.h>
 #include <limits>
 #include <vector>
@@ -35,7 +35,7 @@ TEST(ProbDistributionsBernoulli, distributionCheck) {
 }
 
 TEST(ProbDistributionsBernoulli, error_check) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   EXPECT_NO_THROW(stan::math::bernoulli_rng(0.6, rng));
 
   EXPECT_THROW(stan::math::bernoulli_rng(1.6, rng), std::domain_error);
@@ -45,7 +45,7 @@ TEST(ProbDistributionsBernoulli, error_check) {
 }
 
 TEST(ProbDistributionsBernoulli, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int N = 10000;
 
   std::vector<double> expected;
diff --git a/test/unit/math/prim/prob/beta_binomial_test.cpp b/test/unit/math/prim/prob/beta_binomial_test.cpp
index 9adc12fe794..289d4b1db83 100644
--- a/test/unit/math/prim/prob/beta_binomial_test.cpp
+++ b/test/unit/math/prim/prob/beta_binomial_test.cpp
@@ -1,7 +1,7 @@
 #include <stan/math/prim.hpp>
 #include <test/unit/math/prim/prob/vector_rng_test_helper.hpp>
 #include <test/unit/math/prim/prob/VectorIntRNGTestRig.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <boost/math/distributions.hpp>
 #include <gtest/gtest.h>
 #include <limits>
@@ -40,7 +40,7 @@ TEST(ProbDistributionsBetaBinomial, distributionCheck) {
 }
 
 TEST(ProbDistributionBetaBinomial, error_check) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   EXPECT_NO_THROW(stan::math::beta_binomial_rng(4, 0.6, 2.0, rng));
 
   EXPECT_THROW(stan::math::beta_binomial_rng(-4, 0.6, 2, rng),
diff --git a/test/unit/math/prim/prob/beta_proportion_test.cpp b/test/unit/math/prim/prob/beta_proportion_test.cpp
index f38f7e71e16..7d5ada3ea60 100644
--- a/test/unit/math/prim/prob/beta_proportion_test.cpp
+++ b/test/unit/math/prim/prob/beta_proportion_test.cpp
@@ -1,7 +1,7 @@
 #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 <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <boost/math/distributions.hpp>
 #include <gtest/gtest.h>
 #include <limits>
@@ -50,7 +50,7 @@ TEST(ProbDistributionsBetaProportion, distributionTest) {
 }
 
 TEST(ProbDistributionsBetaProportion, error_check) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   EXPECT_NO_THROW(stan::math::beta_proportion_rng(0.5, 3.0, rng));
   EXPECT_NO_THROW(stan::math::beta_proportion_rng(1e-10, 1e-10, rng));
   EXPECT_THROW(stan::math::beta_proportion_rng(2.0, -1.0, rng),
@@ -68,7 +68,7 @@ TEST(ProbDistributionsBetaProportion, error_check) {
 }
 
 TEST(ProbDistributionsBetaProportion, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int N = 10000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
 
@@ -99,7 +99,7 @@ TEST(ProbDistributionsBetaProportion, chiSquareGoodnessFitTest) {
 }
 
 TEST(ProbDistributionsBetaProportion, chiSquareGoodnessFitTest2) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int N = 10000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
 
diff --git a/test/unit/math/prim/prob/beta_test.cpp b/test/unit/math/prim/prob/beta_test.cpp
index 2fa451a1288..665797a98b1 100644
--- a/test/unit/math/prim/prob/beta_test.cpp
+++ b/test/unit/math/prim/prob/beta_test.cpp
@@ -1,7 +1,7 @@
 #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 <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <boost/math/distributions.hpp>
 #include <gtest/gtest.h>
 #include <limits>
@@ -46,7 +46,7 @@ TEST(ProbDistributionsBeta, distributionTest) {
 }
 
 TEST(ProbDistributionsBeta, error_check) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   EXPECT_NO_THROW(stan::math::beta_rng(2.0, 1.0, rng));
   EXPECT_NO_THROW(stan::math::beta_rng(1e-10, 1e-10, rng));
 
@@ -60,7 +60,7 @@ TEST(ProbDistributionsBeta, error_check) {
 }
 
 TEST(ProbDistributionsBeta, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int N = 10000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
 
@@ -83,7 +83,7 @@ TEST(ProbDistributionsBeta, chiSquareGoodnessFitTest) {
 }
 
 TEST(ProbDistributionsBeta, chiSquareGoodnessFitTestSmallParameters) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int N = 10000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
 
diff --git a/test/unit/math/prim/prob/binomial_test.cpp b/test/unit/math/prim/prob/binomial_test.cpp
index c12507db534..45622d90af8 100644
--- a/test/unit/math/prim/prob/binomial_test.cpp
+++ b/test/unit/math/prim/prob/binomial_test.cpp
@@ -1,7 +1,7 @@
 #include <stan/math/prim.hpp>
 #include <test/unit/math/prim/prob/vector_rng_test_helper.hpp>
 #include <test/unit/math/prim/prob/VectorIntRNGTestRig.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <boost/math/distributions.hpp>
 #include <gtest/gtest.h>
 #include <limits>
@@ -39,7 +39,7 @@ TEST(ProbDistributionsBinomial, distributionCheck) {
 }
 
 TEST(ProbDistributionBinomiali, error_check) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   EXPECT_NO_THROW(stan::math::binomial_rng(4, 0.6, rng));
   EXPECT_THROW(stan::math::binomial_rng(-4, 0.6, rng), std::domain_error);
   EXPECT_THROW(stan::math::binomial_rng(4, -0.6, rng), std::domain_error);
@@ -50,7 +50,7 @@ TEST(ProbDistributionBinomiali, error_check) {
 }
 
 TEST(ProbDistributionsBinomial, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int N = 10000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
   boost::math::binomial_distribution<> dist(100, 0.6);
diff --git a/test/unit/math/prim/prob/categorical_logit_rng_test.cpp b/test/unit/math/prim/prob/categorical_logit_rng_test.cpp
index b60603f28e7..94c5421d560 100644
--- a/test/unit/math/prim/prob/categorical_logit_rng_test.cpp
+++ b/test/unit/math/prim/prob/categorical_logit_rng_test.cpp
@@ -1,13 +1,13 @@
 #include <stan/math/prim.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <boost/math/distributions.hpp>
 #include <limits>
 
 TEST(ProbDistributionsCategoricalLogit, error_check) {
   using Eigen::VectorXd;
   using stan::math::categorical_logit_rng;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
 
   VectorXd beta(3);
 
@@ -27,7 +27,7 @@ TEST(ProbDistributionsCategoricalLogit, error_check) {
 TEST(ProbDistributionsCategoricalLogit, chiSquareGoodnessFitTest) {
   using Eigen::VectorXd;
   using stan::math::softmax;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int N = 10000;
   int K = 3;
   VectorXd beta(K);
diff --git a/test/unit/math/prim/prob/categorical_test.cpp b/test/unit/math/prim/prob/categorical_test.cpp
index 2cbeeb1f82d..3edf96dab7a 100644
--- a/test/unit/math/prim/prob/categorical_test.cpp
+++ b/test/unit/math/prim/prob/categorical_test.cpp
@@ -1,6 +1,6 @@
 #include <stan/math/prim.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <boost/math/distributions.hpp>
 #include <vector>
 #include <limits>
@@ -88,7 +88,7 @@ TEST(ProbDistributionsCategorical, error_check) {
   using Eigen::Dynamic;
   using Eigen::Matrix;
   using stan::math::categorical_log;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
 
   Matrix<double, Dynamic, Dynamic> theta(3, 1);
   theta << 0.15, 0.45, 0.50;
@@ -100,7 +100,7 @@ TEST(ProbDistributionsCategorical, chiSquareGoodnessFitTest) {
   using Eigen::Dynamic;
   using Eigen::Matrix;
   using stan::math::categorical_log;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
 
   int N = 10000;
   Matrix<double, Dynamic, Dynamic> theta(3, 1);
diff --git a/test/unit/math/prim/prob/cauchy_test.cpp b/test/unit/math/prim/prob/cauchy_test.cpp
index 7409f7bcc73..bcd8ba39b6a 100644
--- a/test/unit/math/prim/prob/cauchy_test.cpp
+++ b/test/unit/math/prim/prob/cauchy_test.cpp
@@ -1,7 +1,7 @@
 #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 <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <boost/math/distributions.hpp>
 #include <gtest/gtest.h>
 #include <limits>
@@ -46,7 +46,7 @@ TEST(ProbDistributionsCauchy, distributionTest) {
 }
 
 TEST(ProbDistributionsCauchy, error_check) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   EXPECT_NO_THROW(stan::math::cauchy_rng(2.0, 1.0, rng));
 
   EXPECT_THROW(stan::math::cauchy_rng(2.0, -1.0, rng), std::domain_error);
@@ -59,7 +59,7 @@ TEST(ProbDistributionsCauchy, error_check) {
 }
 
 TEST(ProbDistributionsCauchy, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int N = 10000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
 
diff --git a/test/unit/math/prim/prob/chi_square_test.cpp b/test/unit/math/prim/prob/chi_square_test.cpp
index 1e19844a18f..f105f304e88 100644
--- a/test/unit/math/prim/prob/chi_square_test.cpp
+++ b/test/unit/math/prim/prob/chi_square_test.cpp
@@ -2,7 +2,7 @@
 #include <test/unit/math/prim/prob/vector_rng_test_helper.hpp>
 #include <test/unit/math/prim/prob/util.hpp>
 #include <boost/math/distributions.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <gtest/gtest.h>
 #include <limits>
 #include <vector>
@@ -42,7 +42,7 @@ TEST(ProbDistributionsChiSquare, distributionTest) {
 }
 
 TEST(ProbDistributionsChiSquare, error_check) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   EXPECT_NO_THROW(stan::math::chi_square_rng(2.0, rng));
 
   EXPECT_THROW(stan::math::chi_square_rng(-2.0, rng), std::domain_error);
@@ -51,7 +51,7 @@ TEST(ProbDistributionsChiSquare, error_check) {
 }
 
 TEST(ProbDistributionsChiSquare, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int N = 10000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
 
diff --git a/test/unit/math/prim/prob/dirichlet_multinomial_test.cpp b/test/unit/math/prim/prob/dirichlet_multinomial_test.cpp
index 74de88fe44b..a88b8c35d9a 100644
--- a/test/unit/math/prim/prob/dirichlet_multinomial_test.cpp
+++ b/test/unit/math/prim/prob/dirichlet_multinomial_test.cpp
@@ -1,6 +1,6 @@
 #include <stan/math/prim.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <boost/math/distributions.hpp>
 #include <limits>
 #include <vector>
@@ -90,7 +90,7 @@ TEST(ProbDistributionsDirichletMultinomial, zeros) {
 TEST(ProbDistributionsDirichletMultinomial, rng) {
   using Eigen::Dynamic;
   using Eigen::Matrix;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
 
   Matrix<double, Dynamic, 1> theta(3);
   theta << 0.5, 5.5, 12.0;
@@ -113,7 +113,7 @@ TEST(ProbDistributionsDirichletMultinomial, rng) {
 TEST(ProbDistributionsDirichletMultinomial, rngError) {
   using Eigen::Dynamic;
   using Eigen::Matrix;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
 
   Matrix<double, Dynamic, 1> theta(3);
   theta << 0.5, 1.5, 4.0;
@@ -135,7 +135,7 @@ TEST(ProbDistributionsDirichletMultinomial, rngError) {
 TEST(ProbDistributionsDirichletMultinomial, chiSquareGoodnessFitTest) {
   using Eigen::Dynamic;
   using Eigen::Matrix;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int M = 100;
   int trials = 100;
   int N = M * trials;
@@ -171,7 +171,7 @@ TEST(ProbDistributionsDirichletMultinomial, equivBetaBinomial) {
   using Eigen::Matrix;
   using stan::math::beta_binomial_lpmf;
 
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   std::vector<int> ns;
   Matrix<double, Dynamic, 1> theta(2, 1);
   double alpha = 4.0;
diff --git a/test/unit/math/prim/prob/dirichlet_test.cpp b/test/unit/math/prim/prob/dirichlet_test.cpp
index 9cf7344cef4..9a531eaa35e 100644
--- a/test/unit/math/prim/prob/dirichlet_test.cpp
+++ b/test/unit/math/prim/prob/dirichlet_test.cpp
@@ -1,6 +1,6 @@
 #include <stan/math/prim.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <boost/math/distributions.hpp>
 #include <vector>
 
@@ -139,7 +139,7 @@ double chi_square(std::vector<int> bin, std::vector<double> expect) {
 }
 
 void test_dirichlet3_1(Eigen::VectorXd alpha) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int N = 10000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
 
@@ -168,7 +168,7 @@ void test_dirichlet3_2(Eigen::VectorXd alpha) {
   using Eigen::Dynamic;
   using Eigen::Matrix;
   using Eigen::VectorXd;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int N = 10000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
   boost::math::beta_distribution<> dist(alpha(1), alpha(0) + alpha(2));
@@ -213,7 +213,7 @@ TEST(ProbDistributionsDirichlet, random) {
   using Eigen::Dynamic;
   using Eigen::Matrix;
   using Eigen::VectorXd;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   VectorXd alpha(3);
   alpha << 2.0, 3.0, 11.0;
   EXPECT_NO_THROW(stan::math::dirichlet_rng(alpha, rng));
diff --git a/test/unit/math/prim/prob/discrete_range_test.cpp b/test/unit/math/prim/prob/discrete_range_test.cpp
index d6535d384f3..6b64e70a060 100644
--- a/test/unit/math/prim/prob/discrete_range_test.cpp
+++ b/test/unit/math/prim/prob/discrete_range_test.cpp
@@ -1,12 +1,12 @@
 #include <stan/math/prim.hpp>
 #include <boost/math/distributions.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <gtest/gtest.h>
 #include <vector>
 
 TEST(ProbDistributionsDiscreteRange, error_check) {
   using stan::math::discrete_range_rng;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
 
   std::vector<int> lower{5, 11, -15};
   std::vector<int> upper{7, 15, -10};
@@ -22,7 +22,7 @@ TEST(ProbDistributionsDiscreteRange, error_check) {
 
 TEST(ProbDistributionsDiscreteRange, boundary_values) {
   using stan::math::discrete_range_rng;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
 
   std::vector<int> lower{-5, 11, 17};
   EXPECT_EQ(lower, discrete_range_rng(lower, lower, rng));
@@ -37,7 +37,7 @@ TEST(ProbDistributionsDiscreteRange, boundary_values) {
 }
 
 TEST(ProbDistributionsDiscreteRange, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
 
   int N = 10000;
   int lower = -3;
diff --git a/test/unit/math/prim/prob/double_exponential_test.cpp b/test/unit/math/prim/prob/double_exponential_test.cpp
index 8e8d677e226..38462da9860 100644
--- a/test/unit/math/prim/prob/double_exponential_test.cpp
+++ b/test/unit/math/prim/prob/double_exponential_test.cpp
@@ -2,7 +2,7 @@
 #include <test/unit/math/prim/prob/vector_rng_test_helper.hpp>
 #include <test/unit/math/prim/prob/util.hpp>
 #include <boost/math/distributions.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <gtest/gtest.h>
 #include <limits>
 #include <vector>
@@ -46,7 +46,7 @@ TEST(ProbDistributionsDoubleExponential, distributionTest) {
 }
 
 TEST(ProbDistributionsDoubleExponential, error_check) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   EXPECT_NO_THROW(stan::math::double_exponential_rng(2.0, 1.0, rng));
 
   EXPECT_THROW(stan::math::double_exponential_rng(2.0, -1.0, rng),
@@ -60,7 +60,7 @@ TEST(ProbDistributionsDoubleExponential, error_check) {
 }
 
 TEST(ProbDistributionsDoubleExponential, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int N = 10000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
 
diff --git a/test/unit/math/prim/prob/exp_mod_normal_test.cpp b/test/unit/math/prim/prob/exp_mod_normal_test.cpp
index d817dd10342..b735237e0a1 100644
--- a/test/unit/math/prim/prob/exp_mod_normal_test.cpp
+++ b/test/unit/math/prim/prob/exp_mod_normal_test.cpp
@@ -1,7 +1,7 @@
 #include <stan/math/prim.hpp>
 #include <test/unit/math/prim/prob/vector_rng_test_helper.hpp>
 #include <boost/math/distributions.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <gtest/gtest.h>
 #include <vector>
 
@@ -30,7 +30,7 @@ TEST(ProbDistributionsExpModNormal, errorCheck) {
  * the distributions manually
  */
 TEST(ProbDistributionsExpModNormal, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int N = 10000;
 
   double mu = 2.0;
@@ -93,7 +93,7 @@ TEST(ProbDistributionsExpModNormal, chiSquareGoodnessFitTest) {
 }
 
 TEST(ProbDistributionsExpModNormal, error_check) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   EXPECT_NO_THROW(stan::math::exp_mod_normal_rng(10.0, 2.0, 1.0, rng));
 
   EXPECT_THROW(stan::math::exp_mod_normal_rng(10.0, 2.0, -1.0, rng),
diff --git a/test/unit/math/prim/prob/exponential_test.cpp b/test/unit/math/prim/prob/exponential_test.cpp
index c581dd64508..73f35e95c0c 100644
--- a/test/unit/math/prim/prob/exponential_test.cpp
+++ b/test/unit/math/prim/prob/exponential_test.cpp
@@ -1,7 +1,7 @@
 #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 <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <boost/math/distributions.hpp>
 #include <gtest/gtest.h>
 #include <limits>
@@ -43,7 +43,7 @@ TEST(ProbDistributionsExponential, distributionTest) {
 }
 
 TEST(ProbDistributionsExponential, error_check) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   EXPECT_NO_THROW(stan::math::exponential_rng(2.0, rng));
 
   EXPECT_THROW(stan::math::exponential_rng(-2.0, rng), std::domain_error);
@@ -53,7 +53,7 @@ TEST(ProbDistributionsExponential, error_check) {
 }
 
 TEST(ProbDistributionsExponential, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int N = 10000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
 
diff --git a/test/unit/math/prim/prob/frechet_test.cpp b/test/unit/math/prim/prob/frechet_test.cpp
index 8dd29a48bfb..298add661b7 100644
--- a/test/unit/math/prim/prob/frechet_test.cpp
+++ b/test/unit/math/prim/prob/frechet_test.cpp
@@ -1,7 +1,7 @@
 #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 <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <boost/math/distributions.hpp>
 #include <gtest/gtest.h>
 #include <limits>
@@ -51,7 +51,7 @@ TEST(ProbDistributionsFrechet, distributionTest) {
 }
 
 TEST(ProbDistributionsFrechet, error_check) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   EXPECT_NO_THROW(stan::math::frechet_rng(2.0, 3.0, rng));
 
   EXPECT_THROW(stan::math::frechet_rng(-2.0, 3.0, rng), std::domain_error);
@@ -62,7 +62,7 @@ TEST(ProbDistributionsFrechet, error_check) {
 }
 
 TEST(ProbDistributionsFrechet, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int N = 10000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
 
diff --git a/test/unit/math/prim/prob/gamma_test.cpp b/test/unit/math/prim/prob/gamma_test.cpp
index 63f7d396ee5..2adf6396167 100644
--- a/test/unit/math/prim/prob/gamma_test.cpp
+++ b/test/unit/math/prim/prob/gamma_test.cpp
@@ -1,7 +1,7 @@
 #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 <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <boost/math/distributions.hpp>
 #include <gtest/gtest.h>
 #include <limits>
@@ -46,7 +46,7 @@ TEST(ProbDistributionsGamma, distributionTest) {
 }
 
 TEST(ProbDistributionGamma, error_check) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   EXPECT_NO_THROW(stan::math::gamma_rng(2.0, 3.0, rng));
 
   EXPECT_THROW(stan::math::gamma_rng(-2.0, 3.0, rng), std::domain_error);
@@ -58,7 +58,7 @@ TEST(ProbDistributionGamma, error_check) {
 }
 
 TEST(ProbDistributionGamma, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int N = 10000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
 
diff --git a/test/unit/math/prim/prob/gaussian_dlm_obs_rng_test.cpp b/test/unit/math/prim/prob/gaussian_dlm_obs_rng_test.cpp
index d49ae9a4d47..9167aa533ed 100644
--- a/test/unit/math/prim/prob/gaussian_dlm_obs_rng_test.cpp
+++ b/test/unit/math/prim/prob/gaussian_dlm_obs_rng_test.cpp
@@ -1,7 +1,7 @@
 #include <gtest/gtest.h>
 #include <stan/math/prim.hpp>
 #include <boost/math/distributions.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <test/unit/math/prim/prob/util.hpp>
 #include <exception>
 #include <limits>
@@ -36,7 +36,7 @@ TEST_F(ProbDistributionsGaussianDLMInputsRng, PoliciesF) {
   using Eigen::MatrixXd;
   using Eigen::VectorXd;
   using stan::math::gaussian_dlm_obs_rng;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   MatrixXd FF_sz1 = MatrixXd::Random(4, 3);
   EXPECT_THROW(gaussian_dlm_obs_rng(FF_sz1, GG, V, W, m0, C0, T, rng),
                std::invalid_argument);
@@ -66,7 +66,7 @@ TEST_F(ProbDistributionsGaussianDLMInputsRng, PoliciesG) {
   using Eigen::MatrixXd;
   using Eigen::VectorXd;
   using stan::math::gaussian_dlm_obs_rng;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   // size
   MatrixXd GG_sz1 = MatrixXd::Random(3, 3);
   EXPECT_THROW(gaussian_dlm_obs_rng(FF, GG_sz1, V, W, m0, C0, T, rng),
@@ -97,7 +97,7 @@ TEST_F(ProbDistributionsGaussianDLMInputsRng, PoliciesW) {
   using Eigen::MatrixXd;
   using Eigen::VectorXd;
   using stan::math::gaussian_dlm_obs_rng;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   // Not symmetric
   MatrixXd W_asym = W;
   W_asym(0, 1) = 1;
@@ -155,7 +155,7 @@ TEST_F(ProbDistributionsGaussianDLMInputsRng, PoliciesVMatrix) {
   using Eigen::MatrixXd;
   using Eigen::VectorXd;
   using stan::math::gaussian_dlm_obs_rng;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   // Not symmetric
   MatrixXd V_asym = V;
   V_asym(0, 2) = 1;
@@ -198,7 +198,7 @@ TEST_F(ProbDistributionsGaussianDLMInputsRng, PoliciesVVector) {
   using Eigen::MatrixXd;
   using Eigen::VectorXd;
   using stan::math::gaussian_dlm_obs_rng;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   // negative
   MatrixXd V_neg = V_vec;
   V_neg(0) = -1;
@@ -230,7 +230,7 @@ TEST_F(ProbDistributionsGaussianDLMInputsRng, Policiesm0) {
   using Eigen::MatrixXd;
   using Eigen::VectorXd;
   using stan::math::gaussian_dlm_obs_rng;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   // size
   VectorXd m0_sz = VectorXd::Zero(4, 1);
   EXPECT_THROW(gaussian_dlm_obs_rng(FF, GG, V, W, m0_sz, C0, T, rng),
@@ -256,7 +256,7 @@ TEST_F(ProbDistributionsGaussianDLMInputsRng, PoliciesC0) {
   using Eigen::MatrixXd;
   using Eigen::VectorXd;
   using stan::math::gaussian_dlm_obs_rng;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   // size
   MatrixXd C0_sz = MatrixXd::Identity(3, 3);
   EXPECT_THROW(gaussian_dlm_obs_rng(FF, GG, V, W, m0, C0_sz, T, rng),
@@ -293,7 +293,7 @@ TEST_F(ProbDistributionsGaussianDLMInputsRng, PoliciesT) {
   using Eigen::MatrixXd;
   using Eigen::VectorXd;
   using stan::math::gaussian_dlm_obs_rng;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   // Must be positive.
   EXPECT_THROW(gaussian_dlm_obs_rng(FF, GG, V, W, m0, C0, 0, rng),
                std::domain_error);
@@ -305,7 +305,7 @@ TEST_F(ProbDistributionsGaussianDLMInputsRng, chiSquaredGoodnessOfFit) {
   using Eigen::MatrixXd;
   using Eigen::VectorXd;
   using stan::math::gaussian_dlm_obs_rng;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
 
   // With identity state transition and initial mean 0 and identity
   // covariance matrices:
diff --git a/test/unit/math/prim/prob/gumbel_test.cpp b/test/unit/math/prim/prob/gumbel_test.cpp
index 9de4ec2c765..b8b26aeb9bb 100644
--- a/test/unit/math/prim/prob/gumbel_test.cpp
+++ b/test/unit/math/prim/prob/gumbel_test.cpp
@@ -2,7 +2,7 @@
 #include <test/unit/math/prim/prob/vector_rng_test_helper.hpp>
 #include <test/unit/math/prim/prob/util.hpp>
 #include <boost/math/distributions.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <gtest/gtest.h>
 #include <limits>
 #include <vector>
@@ -46,7 +46,7 @@ TEST(ProbDistributionsGumbel, distributionTest) {
 }
 
 TEST(ProbDistributionsGumbel, error_check) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   EXPECT_NO_THROW(stan::math::gumbel_rng(10.0, 2.0, rng));
 
   EXPECT_THROW(
@@ -56,7 +56,7 @@ TEST(ProbDistributionsGumbel, error_check) {
 }
 
 TEST(ProbDistributionsGumbel, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int N = 10000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
 
diff --git a/test/unit/math/prim/prob/hmm_latent_rng_test.cpp b/test/unit/math/prim/prob/hmm_latent_rng_test.cpp
index 0a5606e9abd..33d08340fe2 100644
--- a/test/unit/math/prim/prob/hmm_latent_rng_test.cpp
+++ b/test/unit/math/prim/prob/hmm_latent_rng_test.cpp
@@ -26,7 +26,7 @@ TEST(hmm_rng_test, chiSquareGoodnessFitTest) {
   Eigen::MatrixXd log_omegas
       = Eigen::MatrixXd::Ones(n_states, n_transitions + 1);
 
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int N = 10000;
 
   std::vector<double> expected;
@@ -68,7 +68,7 @@ TEST(hmm_rng_test, chiSquareGoodnessFitTest_symmetric) {
   Eigen::MatrixXd log_omegas
       = Eigen::MatrixXd::Ones(n_states, n_transitions + 1);
 
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int N = 10000;
 
   std::vector<double> expected_0;
diff --git a/test/unit/math/prim/prob/hypergeometric_test.cpp b/test/unit/math/prim/prob/hypergeometric_test.cpp
index d84fc39c781..be4d0f6dffa 100644
--- a/test/unit/math/prim/prob/hypergeometric_test.cpp
+++ b/test/unit/math/prim/prob/hypergeometric_test.cpp
@@ -1,11 +1,11 @@
 #include <stan/math/prim.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <boost/math/distributions.hpp>
 #include <vector>
 
 TEST(ProbDistributionsHypergeometric, error_check) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   EXPECT_NO_THROW(stan::math::hypergeometric_rng(10, 10, 15, rng));
 
   EXPECT_THROW(stan::math::hypergeometric_rng(30, 10, 15, rng),
@@ -19,7 +19,7 @@ TEST(ProbDistributionsHypergeometric, error_check) {
 }
 
 TEST(ProbDistributionsHypergeometric, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int N = 10000;
   int num_draws = 10;
   int K = num_draws;
diff --git a/test/unit/math/prim/prob/inv_chi_square_test.cpp b/test/unit/math/prim/prob/inv_chi_square_test.cpp
index f9fd797e1c6..bc95a4b93df 100644
--- a/test/unit/math/prim/prob/inv_chi_square_test.cpp
+++ b/test/unit/math/prim/prob/inv_chi_square_test.cpp
@@ -1,7 +1,7 @@
 #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 <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <boost/math/distributions.hpp>
 #include <gtest/gtest.h>
 #include <limits>
@@ -42,7 +42,7 @@ TEST(ProbDistributionsInvChiSquare, distributionTest) {
 }
 
 TEST(ProbDistributionsInvChiSquare, error_check) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   EXPECT_NO_THROW(stan::math::inv_chi_square_rng(4.0, rng));
 
   EXPECT_THROW(stan::math::inv_chi_square_rng(-4.0, rng), std::domain_error);
@@ -52,7 +52,7 @@ TEST(ProbDistributionsInvChiSquare, error_check) {
 }
 
 TEST(ProbDistributionsInvChiSquare, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int N = 10000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
 
diff --git a/test/unit/math/prim/prob/inv_gamma_test.cpp b/test/unit/math/prim/prob/inv_gamma_test.cpp
index 4687a03df55..8f979dc689d 100644
--- a/test/unit/math/prim/prob/inv_gamma_test.cpp
+++ b/test/unit/math/prim/prob/inv_gamma_test.cpp
@@ -1,7 +1,7 @@
 #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 <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <boost/math/distributions.hpp>
 #include <gtest/gtest.h>
 #include <vector>
@@ -46,7 +46,7 @@ TEST(ProbDistributionsInvGamma, distributionTest) {
 }
 
 TEST(ProbDistributionsInvGamma, error_check) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   EXPECT_NO_THROW(stan::math::inv_gamma_rng(4.0, 3.0, rng));
 
   EXPECT_THROW(stan::math::inv_gamma_rng(-4.0, 3.0, rng), std::domain_error);
@@ -60,7 +60,7 @@ TEST(ProbDistributionsInvGamma, error_check) {
 }
 
 TEST(ProbDistributionsInvGamma, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int N = 10000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
 
diff --git a/test/unit/math/prim/prob/inv_wishart_cholesky_rng_test.cpp b/test/unit/math/prim/prob/inv_wishart_cholesky_rng_test.cpp
index 95369fffa58..42b3a1ddf20 100644
--- a/test/unit/math/prim/prob/inv_wishart_cholesky_rng_test.cpp
+++ b/test/unit/math/prim/prob/inv_wishart_cholesky_rng_test.cpp
@@ -1,6 +1,6 @@
 #include <stan/math/prim.hpp>
 #include <test/unit/math/prim/util.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <boost/math/distributions/chi_squared.hpp>
 #include <gtest/gtest.h>
 #include <stdexcept>
@@ -13,7 +13,7 @@ TEST(ProbDistributionsInvWishartCholesky, rng) {
 
   using stan::math::inv_wishart_cholesky_rng;
 
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
 
   MatrixXd omega(3, 4);
   EXPECT_THROW(inv_wishart_cholesky_rng(3.0, omega, rng),
@@ -35,7 +35,7 @@ TEST(ProbDistributionsInvWishartCholesky, rng_pos_def) {
   using Eigen::MatrixXd;
   using stan::math::inv_wishart_cholesky_rng;
 
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
 
   MatrixXd Sigma(2, 2);
   MatrixXd Sigma_non_pos_def(2, 2);
@@ -58,7 +58,7 @@ TEST(ProbDistributionsInvWishartCholesky, marginalTwoChiSquareGoodnessFitTest) {
   using stan::math::inv_wishart_cholesky_rng;
   using std::log;
 
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   MatrixXd sigma(3, 3);
   sigma << 9.0, -3.0, 2.0, -3.0, 4.0, 0.0, 2.0, 0.0, 3.0;
 
@@ -91,7 +91,7 @@ TEST(ProbDistributionsInvWishartCholesky, SpecialRNGTest) {
   using stan::math::inv_wishart_cholesky_rng;
   using stan::math::multiply_lower_tri_self_transpose;
 
-  boost::random::mt19937 rng(92343U);
+  boost::random::mixmax rng(9234U);
   int N = 1e5;
   double tol = 0.1;
   for (int k = 1; k < 5; k++) {
@@ -123,7 +123,7 @@ TEST(ProbDistributionsInvWishartCholesky, compareToInvWishart) {
   using stan::math::multiply_lower_tri_self_transpose;
   using stan::math::qr_thin_Q;
 
-  boost::random::mt19937 rng(92343U);
+  boost::random::mixmax rng(9234U);
   int N = 1e4;
   double tol = 0.05;
   for (int k = 1; k < 4; k++) {
diff --git a/test/unit/math/prim/prob/inv_wishart_cholesky_test.cpp b/test/unit/math/prim/prob/inv_wishart_cholesky_test.cpp
index 062de6d554c..d625c47b14f 100644
--- a/test/unit/math/prim/prob/inv_wishart_cholesky_test.cpp
+++ b/test/unit/math/prim/prob/inv_wishart_cholesky_test.cpp
@@ -1,6 +1,6 @@
 #include <stan/math/prim.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <boost/math/distributions.hpp>
 #include <boost/math/special_functions/digamma.hpp>
 
diff --git a/test/unit/math/prim/prob/inv_wishart_rng_test.cpp b/test/unit/math/prim/prob/inv_wishart_rng_test.cpp
index 77c953e771a..a2104441ac4 100644
--- a/test/unit/math/prim/prob/inv_wishart_rng_test.cpp
+++ b/test/unit/math/prim/prob/inv_wishart_rng_test.cpp
@@ -1,14 +1,14 @@
 #include <stan/math/prim.hpp>
 #include <test/unit/math/prim/util.hpp>
 #include <boost/math/distributions/chi_squared.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <gtest/gtest.h>
 #include <stdexcept>
 
 TEST(ProbDistributionsInvWishart, rng) {
   using Eigen::MatrixXd;
   using stan::math::inv_wishart_rng;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
 
   MatrixXd omega(3, 4);
   EXPECT_THROW(inv_wishart_rng(3.0, omega, rng), std::invalid_argument);
@@ -26,7 +26,7 @@ TEST(ProbDistributionsInvWishart, rng_pos_def) {
   using Eigen::MatrixXd;
   using stan::math::inv_wishart_rng;
 
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
 
   MatrixXd Sigma(2, 2);
   MatrixXd Sigma_non_pos_def(2, 2);
@@ -46,7 +46,7 @@ TEST(ProbDistributionsInvWishart, rng_symmetry) {
   using stan::test::unit::expect_symmetric;
   using stan::test::unit::spd_rng;
 
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   for (int k = 1; k < 20; ++k)
     for (double nu = k - 0.5; nu < k + 20; ++nu)
       for (int n = 0; n < 10; ++n)
@@ -61,7 +61,7 @@ TEST(ProbDistributionsInvWishart, chiSquareGoodnessFitTest) {
   using stan::math::inv_wishart_rng;
   using std::log;
 
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   MatrixXd sigma(3, 3);
   sigma << 9.0, -3.0, 0.0, -3.0, 4.0, 1.0, 0.0, 1.0, 3.0;
   int N = 10000;
@@ -92,7 +92,7 @@ TEST(ProbDistributionsInvWishart, SpecialRNGTest) {
   using Eigen::MatrixXd;
   using stan::math::inv_wishart_rng;
 
-  boost::random::mt19937 rng(1234U);
+  boost::random::mixmax rng(1234U);
   int N = 1e5;
   double tol = 0.1;
   for (int k = 1; k < 5; k++) {
diff --git a/test/unit/math/prim/prob/inv_wishart_test.cpp b/test/unit/math/prim/prob/inv_wishart_test.cpp
index e7871b6c8ef..76f922294e9 100644
--- a/test/unit/math/prim/prob/inv_wishart_test.cpp
+++ b/test/unit/math/prim/prob/inv_wishart_test.cpp
@@ -1,6 +1,6 @@
 #include <stan/math/prim.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <boost/math/distributions.hpp>
 #include <boost/math/special_functions/digamma.hpp>
 
diff --git a/test/unit/math/prim/prob/lkj_corr_test.cpp b/test/unit/math/prim/prob/lkj_corr_test.cpp
index 49d0400aa18..8b362fffa94 100644
--- a/test/unit/math/prim/prob/lkj_corr_test.cpp
+++ b/test/unit/math/prim/prob/lkj_corr_test.cpp
@@ -1,10 +1,10 @@
 #include <stan/math/prim.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <boost/math/distributions.hpp>
 
 TEST(ProbDistributionsLkjCorr, testIdentity) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   unsigned int K = 4;
   Eigen::MatrixXd Sigma(K, K);
   Sigma.setZero();
@@ -18,7 +18,7 @@ TEST(ProbDistributionsLkjCorr, testIdentity) {
 }
 
 TEST(ProbDistributionsLkjCorr, testHalf) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   unsigned int K = 4;
   Eigen::MatrixXd Sigma(K, K);
   Sigma.setConstant(0.5);
@@ -33,7 +33,7 @@ TEST(ProbDistributionsLkjCorr, testHalf) {
 }
 
 TEST(ProbDistributionsLkjCorr, Sigma) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   unsigned int K = 4;
   Eigen::MatrixXd Sigma(K, K);
   Sigma.setZero();
@@ -55,7 +55,7 @@ TEST(ProbDistributionsLkjCorr, Sigma) {
 }
 
 TEST(ProbDistributionsLKJCorr, error_check) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   EXPECT_NO_THROW(stan::math::lkj_corr_cholesky_rng(5, 1.0, rng));
   EXPECT_NO_THROW(stan::math::lkj_corr_rng(5, 1.0, rng));
 
@@ -65,7 +65,7 @@ TEST(ProbDistributionsLKJCorr, error_check) {
 }
 
 TEST(ProbDistributionsLKJCorr, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int N = 10000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
   boost::math::beta_distribution<> dist(2.5, 2.5);
@@ -101,7 +101,7 @@ TEST(ProbDistributionsLKJCorr, chiSquareGoodnessFitTest) {
 }
 
 TEST(ProbDistributionsLkjCorrCholesky, testIdentity) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   unsigned int K = 4;
   Eigen::MatrixXd Sigma(K, K);
   Sigma.setZero();
@@ -115,7 +115,7 @@ TEST(ProbDistributionsLkjCorrCholesky, testIdentity) {
 }
 
 TEST(ProbDistributionsLkjCorrCholesky, testHalf) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   unsigned int K = 4;
   Eigen::MatrixXd Sigma(K, K);
   Sigma.setConstant(0.5);
diff --git a/test/unit/math/prim/prob/lkj_cov_test.cpp b/test/unit/math/prim/prob/lkj_cov_test.cpp
index 5e10b5d9ce5..ba23dcf53c3 100644
--- a/test/unit/math/prim/prob/lkj_cov_test.cpp
+++ b/test/unit/math/prim/prob/lkj_cov_test.cpp
@@ -1,10 +1,10 @@
 #include <stan/math/prim.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <boost/math/distributions.hpp>
 
 TEST(ProbDistributionsLkjCov, testIdentity) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   unsigned int K = 4;
   Eigen::MatrixXd Sigma(K, K);
   Sigma.setZero();
@@ -26,7 +26,7 @@ TEST(ProbDistributionsLkjCov, testIdentity) {
 }
 
 TEST(ProbDistributionsLkjCov, testHalf) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   unsigned int K = 4;
   Eigen::MatrixXd Sigma(K, K);
   Sigma.setConstant(0.5);
@@ -49,7 +49,7 @@ TEST(ProbDistributionsLkjCov, testHalf) {
 }
 
 TEST(ProbDistributionsLkjCov, ErrorChecks) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   unsigned int K = 4;
   Eigen::MatrixXd Sigma(K, K);
   Sigma.setZero();
diff --git a/test/unit/math/prim/prob/logistic_test.cpp b/test/unit/math/prim/prob/logistic_test.cpp
index 4af72d4c321..5e69483d319 100644
--- a/test/unit/math/prim/prob/logistic_test.cpp
+++ b/test/unit/math/prim/prob/logistic_test.cpp
@@ -1,7 +1,7 @@
 #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 <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <boost/math/distributions.hpp>
 #include <gtest/gtest.h>
 #include <limits>
@@ -46,7 +46,7 @@ TEST(ProbDistributionsLogistic, distributionTest) {
 }
 
 TEST(ProbDistributionsLogistic, error_check) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   EXPECT_NO_THROW(stan::math::logistic_rng(4.0, 3.0, rng));
 
   EXPECT_THROW(stan::math::logistic_rng(4.0, -3.0, rng), std::domain_error);
@@ -59,7 +59,7 @@ TEST(ProbDistributionsLogistic, error_check) {
 }
 
 TEST(ProbDistributionsLogistic, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int N = 10000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
 
diff --git a/test/unit/math/prim/prob/loglogistic_test.cpp b/test/unit/math/prim/prob/loglogistic_test.cpp
index 94a38ea4cc7..a9cfad213ab 100644
--- a/test/unit/math/prim/prob/loglogistic_test.cpp
+++ b/test/unit/math/prim/prob/loglogistic_test.cpp
@@ -3,7 +3,7 @@
 #include <test/unit/math/prim/prob/util.hpp>
 #include <gtest/gtest.h>
 #include <boost/math/distributions.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <limits>
 #include <vector>
 
@@ -32,7 +32,7 @@ TEST(ProbDistributionsLoglogistic, errorCheck) {
 }
 
 TEST(ProbDistributionsLoglogistic, error_check) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax 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);
@@ -57,7 +57,7 @@ TEST(ProbDistributionsLoglogistic, test_sampling_icdf) {
 }
 
 TEST(ProbDistributionsLoglogistic, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int N = 10000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
 
diff --git a/test/unit/math/prim/prob/lognormal_test.cpp b/test/unit/math/prim/prob/lognormal_test.cpp
index bfc16aaa881..a4d8dd37aaf 100644
--- a/test/unit/math/prim/prob/lognormal_test.cpp
+++ b/test/unit/math/prim/prob/lognormal_test.cpp
@@ -2,7 +2,7 @@
 #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/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <boost/math/distributions.hpp>
 #include <limits>
 #include <vector>
@@ -46,7 +46,7 @@ TEST(ProbDistributionsLogNormalMat, chiSquareGoodnessFitTest) {
 }
 
 TEST(ProbDistributionsLogNormalPrim, error_check) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   EXPECT_NO_THROW(stan::math::lognormal_rng(2.0, 1.0, rng));
 
   EXPECT_THROW(stan::math::lognormal_rng(2.0, -1.0, rng), std::domain_error);
@@ -59,7 +59,7 @@ TEST(ProbDistributionsLogNormalPrim, error_check) {
 }
 
 TEST(ProbDistributionsLogNormalPrim, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int N = 10000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
 
diff --git a/test/unit/math/prim/prob/matrix_normal_prec_rng_test.cpp b/test/unit/math/prim/prob/matrix_normal_prec_rng_test.cpp
index 32a550ef72b..e271049bf0d 100644
--- a/test/unit/math/prim/prob/matrix_normal_prec_rng_test.cpp
+++ b/test/unit/math/prim/prob/matrix_normal_prec_rng_test.cpp
@@ -1,6 +1,6 @@
 #include <gtest/gtest.h>
 #include <boost/math/distributions.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <stan/math/prim.hpp>
 #include <test/unit/math/prim/prob/util.hpp>
 #include <exception>
@@ -10,7 +10,7 @@
 TEST(ProbDistributionsMatrixNormalPrecRng, ErrorSigma) {
   using Eigen::MatrixXd;
   using stan::math::matrix_normal_prec_rng;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   double nan = std::numeric_limits<double>::quiet_NaN();
   double inf = std::numeric_limits<double>::infinity();
   double ninf = -std::numeric_limits<double>::infinity();
@@ -58,7 +58,7 @@ TEST(ProbDistributionsMatrixNormalPrecRng, ErrorSigma) {
 TEST(ProbDistributionsMatrixNormalPrecRng, ErrorD) {
   using Eigen::MatrixXd;
   using stan::math::matrix_normal_prec_rng;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   double nan = std::numeric_limits<double>::quiet_NaN();
   double inf = std::numeric_limits<double>::infinity();
   double ninf = -std::numeric_limits<double>::infinity();
@@ -106,7 +106,7 @@ TEST(ProbDistributionsMatrixNormalPrecRng, ErrorD) {
 TEST(ProbDistributionsMatrixNormalPrecRng, ErrorSize) {
   using Eigen::MatrixXd;
   using stan::math::matrix_normal_prec_rng;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
 
   MatrixXd Mu = MatrixXd::Zero(3, 5);
 
@@ -170,7 +170,7 @@ std::vector<double> extract_sum_of_entries(
 TEST(ProbDistributionsMatrixNormalPrecRng, marginalChiSquareGoodnessFitTest) {
   using Eigen::MatrixXd;
   using stan::math::matrix_normal_prec_rng;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
 
   MatrixXd Mu(2, 3);
   Mu << 1, 2, 3, 4, 5, 6;
diff --git a/test/unit/math/prim/prob/multi_normal_cholesky_test.cpp b/test/unit/math/prim/prob/multi_normal_cholesky_test.cpp
index 798741f3b0d..7c96b34ff9f 100644
--- a/test/unit/math/prim/prob/multi_normal_cholesky_test.cpp
+++ b/test/unit/math/prim/prob/multi_normal_cholesky_test.cpp
@@ -1,6 +1,6 @@
 #include <stan/math/prim.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <boost/math/distributions.hpp>
 #include <vector>
 
@@ -8,7 +8,7 @@ TEST(ProbDistributionsMultiNormalCholesky, NotVectorized) {
   using Eigen::Dynamic;
   using Eigen::Matrix;
   using std::vector;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   Matrix<double, Dynamic, 1> y(3, 1);
   y << 2.0, -2.0, 11.0;
   Matrix<double, Dynamic, 1> mu(3, 1);
@@ -23,7 +23,7 @@ TEST(ProbDistributionsMultiNormalCholesky, Vectorized) {
   using Eigen::Dynamic;
   using Eigen::Matrix;
   using std::vector;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   vector<Matrix<double, Dynamic, 1> > vec_y(2);
   vector<Matrix<double, 1, Dynamic> > vec_y_t(2);
   Matrix<double, Dynamic, 1> y(3);
@@ -89,7 +89,7 @@ TEST(ProbDistributionsMultiNormalCholesky, MultiNormalOneRow) {
   using Eigen::Dynamic;
   using Eigen::Matrix;
   using std::vector;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   Matrix<double, 1, Dynamic> y(3);
   y << 2.0, -2.0, 11.0;
   Matrix<double, Dynamic, 1> mu(3);
@@ -105,7 +105,7 @@ TEST(ProbDistributionsMultiNormalCholesky, error_check) {
   using Eigen::Dynamic;
   using Eigen::Matrix;
   using std::vector;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   Matrix<double, Dynamic, 1> mu(3, 1);
   mu << 2.0, -2.0, 11.0;
 
@@ -124,7 +124,7 @@ TEST(ProbDistributionsMultiNormalCholesky,
   using Eigen::Dynamic;
   using Eigen::Matrix;
   using std::vector;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   Matrix<double, Dynamic, Dynamic> sigma(3, 3);
   sigma << 9.0, -3.0, 0.0, -3.0, 4.0, 1.0, 0.0, 1.0, 3.0;
   Matrix<double, Dynamic, Dynamic> L = sigma.llt().matrixL();
@@ -173,7 +173,7 @@ TEST(ProbDistributionsMultiNormalCholesky,
   using Eigen::Dynamic;
   using Eigen::Matrix;
   using std::vector;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   Matrix<double, Dynamic, Dynamic> sigma(3, 3);
   sigma << 9.0, -3.0, 0.0, -3.0, 4.0, 1.0, 0.0, 1.0, 3.0;
   Matrix<double, Dynamic, Dynamic> L = sigma.llt().matrixL();
@@ -222,7 +222,7 @@ TEST(ProbDistributionsMultiNormalCholesky,
   using Eigen::Dynamic;
   using Eigen::Matrix;
   using std::vector;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   Matrix<double, Dynamic, Dynamic> sigma(3, 3);
   sigma << 9.0, -3.0, 0.0, -3.0, 4.0, 1.0, 0.0, 1.0, 16.0;
   Matrix<double, Dynamic, Dynamic> L = sigma.llt().matrixL();
diff --git a/test/unit/math/prim/prob/multi_normal_prec_rng_test.cpp b/test/unit/math/prim/prob/multi_normal_prec_rng_test.cpp
index 0611fb25650..4a34ae1ca1b 100644
--- a/test/unit/math/prim/prob/multi_normal_prec_rng_test.cpp
+++ b/test/unit/math/prim/prob/multi_normal_prec_rng_test.cpp
@@ -1,6 +1,6 @@
 #include <gtest/gtest.h>
 #include <boost/math/distributions.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <stan/math/prim.hpp>
 #include <test/unit/math/prim/prob/util.hpp>
 #include <limits>
@@ -8,7 +8,7 @@
 
 TEST(ProbDistributionsMultiNormalPrec, vectorized) {
   // Test scalar/vector combinations.
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
 
   Eigen::VectorXd mu(3);
   Eigen::RowVectorXd mu_t(3);
@@ -32,7 +32,7 @@ TEST(ProbDistributionsMultiNormalPrec, vectorized) {
 }
 
 TEST(ProbDistributionsMultiNormalPrec, policiesSigma) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
 
   double nan = std::numeric_limits<double>::quiet_NaN();
   double inf = std::numeric_limits<double>::infinity();
@@ -79,7 +79,7 @@ TEST(ProbDistributionsMultiNormalPrec, policiesSigma) {
 }
 
 TEST(ProbDistributionsMultiNormalPrec, policiesMu) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
 
   double nan = std::numeric_limits<double>::quiet_NaN();
   double inf = std::numeric_limits<double>::infinity();
@@ -108,7 +108,7 @@ TEST(ProbDistributionsMultiNormalPrec, policiesMu) {
 }
 
 TEST(ProbDistributionsMultiNormalPrec, SizeMismatch) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   Eigen::VectorXd mu(2);
   mu << 1.0, -1.0;
   Eigen::MatrixXd Sigma(3, 3);
@@ -118,7 +118,7 @@ TEST(ProbDistributionsMultiNormalPrec, SizeMismatch) {
 }
 
 TEST(ProbDistributionsMultiNormalPrec, marginalOneChiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
 
   Eigen::MatrixXd sigma(3, 3);
   sigma << 9.0, -3.0, 0.0, -3.0, 4.0, 1.0, 0.0, 1.0, 3.0;
@@ -152,7 +152,7 @@ TEST(ProbDistributionsMultiNormalPrec, marginalOneChiSquareGoodnessFitTest) {
 }
 
 TEST(ProbDistributionsMultiNormalPrec, marginalTwoChiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
 
   Eigen::MatrixXd sigma(3, 3);
   sigma << 9.0, -3.0, 0.0, -3.0, 4.0, 1.0, 0.0, 1.0, 3.0;
@@ -186,7 +186,7 @@ TEST(ProbDistributionsMultiNormalPrec, marginalTwoChiSquareGoodnessFitTest) {
 }
 
 TEST(ProbDistributionsMultiNormalPrec, marginalThreeChiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
 
   Eigen::MatrixXd sigma(3, 3);
   sigma << 9.0, -3.0, 0.0, -3.0, 4.0, 1.0, 0.0, 1.0, 16.0;
diff --git a/test/unit/math/prim/prob/multi_normal_test.cpp b/test/unit/math/prim/prob/multi_normal_test.cpp
index 4cf0f4bcf3d..54e7e305204 100644
--- a/test/unit/math/prim/prob/multi_normal_test.cpp
+++ b/test/unit/math/prim/prob/multi_normal_test.cpp
@@ -1,6 +1,6 @@
 #include <stan/math/prim.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <boost/math/distributions.hpp>
 #include <limits>
 #include <vector>
@@ -10,7 +10,7 @@ TEST(ProbDistributionsMultiNormal, NotVectorized) {
   using Eigen::Matrix;
   using std::vector;
 
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   Matrix<double, Dynamic, 1> y(3, 1);
   y << 2.0, -2.0, 11.0;
   Matrix<double, Dynamic, 1> mu(3, 1);
@@ -25,7 +25,7 @@ TEST(ProbDistributionsMultiNormal, Vectorized) {
   using Eigen::Dynamic;
   using Eigen::Matrix;
   using std::vector;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   vector<Matrix<double, Dynamic, 1> > vec_y(2);
   vector<Matrix<double, 1, Dynamic> > vec_y_t(2);
   Matrix<double, Dynamic, 1> y(3);
@@ -89,7 +89,7 @@ TEST(ProbDistributionsMultiNormal, Sigma) {
   using Eigen::Dynamic;
   using Eigen::Matrix;
   using std::vector;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   Matrix<double, Dynamic, 1> y(2, 1);
   y << 2.0, -2.0;
   Matrix<double, Dynamic, 1> mu(2, 1);
@@ -108,7 +108,7 @@ TEST(ProbDistributionsMultiNormal, Mu) {
   using Eigen::Dynamic;
   using Eigen::Matrix;
   using std::vector;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   Matrix<double, Dynamic, 1> y(3, 1);
   y << 2.0, -2.0, 11.0;
   Matrix<double, Dynamic, 1> mu(3, 1);
@@ -132,7 +132,7 @@ TEST(ProbDistributionsMultiNormal, MultiNormalOneRow) {
   using Eigen::Dynamic;
   using Eigen::Matrix;
   using std::vector;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   Matrix<double, 1, Dynamic> y(1, 3);
   y << 2.0, -2.0, 11.0;
   Matrix<double, Dynamic, 1> mu(3, 1);
@@ -147,7 +147,7 @@ TEST(ProbDistributionsMultiNormal, SigmaMultiRow) {
   using Eigen::Dynamic;
   using Eigen::Matrix;
   using std::vector;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   Matrix<double, 1, Dynamic> y(1, 2);
   y << 2.0, -2.0;
   Matrix<double, Dynamic, 1> mu(2, 1);
@@ -172,7 +172,7 @@ TEST(ProbDistributionsMultiNormal, MuMultiRow) {
   using Eigen::Dynamic;
   using Eigen::Matrix;
   using std::vector;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   Matrix<double, 1, Dynamic> y(1, 3);
   y << 2.0, -2.0, 11.0;
   Matrix<double, Dynamic, 1> mu(3, 1);
@@ -196,7 +196,7 @@ TEST(ProbDistributionsMultiNormal, SizeMismatch) {
   using Eigen::Dynamic;
   using Eigen::Matrix;
   using std::vector;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   Matrix<double, 1, Dynamic> y(1, 3);
   y << 2.0, -2.0, 11.0;
   Matrix<double, Dynamic, 1> mu(2, 1);
@@ -213,7 +213,7 @@ TEST(ProbDistributionsMultiNormal, marginalOneChiSquareGoodnessFitTest) {
   using Eigen::Dynamic;
   using Eigen::Matrix;
   using std::vector;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   Matrix<double, Dynamic, Dynamic> sigma(3, 3);
   sigma << 9.0, -3.0, 0.0, -3.0, 4.0, 1.0, 0.0, 1.0, 3.0;
   std::vector<Matrix<double, Dynamic, 1> > mu(3);
@@ -260,7 +260,7 @@ TEST(ProbDistributionsMultiNormal, marginalTwoChiSquareGoodnessFitTest) {
   using Eigen::Dynamic;
   using Eigen::Matrix;
   using std::vector;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   Matrix<double, Dynamic, Dynamic> sigma(3, 3);
   sigma << 9.0, -3.0, 0.0, -3.0, 4.0, 1.0, 0.0, 1.0, 3.0;
   std::vector<Matrix<double, 1, Dynamic> > mu(3);
@@ -307,7 +307,7 @@ TEST(ProbDistributionsMultiNormal, marginalThreeChiSquareGoodnessFitTest) {
   using Eigen::Dynamic;
   using Eigen::Matrix;
   using std::vector;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   Matrix<double, Dynamic, Dynamic> sigma(3, 3);
   sigma << 9.0, -3.0, 0.0, -3.0, 4.0, 1.0, 0.0, 1.0, 16.0;
   Matrix<double, Dynamic, 1> mu(3, 1);
@@ -354,7 +354,7 @@ TEST(multiNormalRng, nonPosDefErrorTest) {
   S << 0, 1, 1, 0;  // not pos definite
   Eigen::VectorXd mu(2);
   mu << 1, 2;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   EXPECT_THROW(multi_normal_rng(mu, S, rng), std::domain_error);
 }
 
diff --git a/test/unit/math/prim/prob/multi_student_t_cholesky_test.cpp b/test/unit/math/prim/prob/multi_student_t_cholesky_test.cpp
index 5111af3312b..81579413432 100644
--- a/test/unit/math/prim/prob/multi_student_t_cholesky_test.cpp
+++ b/test/unit/math/prim/prob/multi_student_t_cholesky_test.cpp
@@ -1,6 +1,6 @@
 #include <stan/math/prim.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <boost/math/distributions.hpp>
 #include <vector>
 #include <limits>
@@ -11,7 +11,7 @@ TEST(ProbDistributionsMultiStudentTCholesky, NotVectorized) {
   using stan::math::multi_student_t_cholesky_lpdf;
   using stan::math::multi_student_t_cholesky_rng;
   using std::vector;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   Matrix<double, Dynamic, 1> y(3, 1);
   y << 2.0, -2.0, 11.0;
   Matrix<double, Dynamic, 1> mu(3, 1);
@@ -34,7 +34,7 @@ TEST(ProbDistributionsMultiStudentTCholesky, Vectorized) {
   using stan::math::multi_student_t_cholesky_lpdf;
   using stan::math::multi_student_t_cholesky_rng;
   using std::vector;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   vector<Matrix<double, Dynamic, 1> > vec_y(2);
   vector<Matrix<double, 1, Dynamic> > vec_y_t(2);
   Matrix<double, Dynamic, 1> y(3);
@@ -115,7 +115,7 @@ TEST(ProbDistributionsMultiStudentTCholesky, Sigma) {
   using stan::math::multi_student_t_cholesky_lpdf;
   using stan::math::multi_student_t_cholesky_rng;
   using std::vector;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   Matrix<double, Dynamic, 1> y(3, 1);
   y << 2.0, -2.0, 11.0;
   Matrix<double, Dynamic, 1> mu(3, 1);
@@ -139,7 +139,7 @@ TEST(ProbDistributionsMultiStudentTCholesky, Mu) {
   using stan::math::multi_student_t_cholesky_lpdf;
   using stan::math::multi_student_t_cholesky_rng;
   using std::vector;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   Matrix<double, Dynamic, 1> y(3, 1);
   y << 2.0, -2.0, 11.0;
   Matrix<double, Dynamic, 1> mu(3, 1);
@@ -195,7 +195,7 @@ TEST(ProbDistributionsMultiStudentTCholesky, Nu) {
   using stan::math::multi_student_t_cholesky_lpdf;
   using stan::math::multi_student_t_cholesky_rng;
   using std::vector;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   Matrix<double, Dynamic, 1> y(3, 1);
   y << 2.0, -2.0, 11.0;
   Matrix<double, Dynamic, 1> mu(3, 1);
@@ -253,7 +253,7 @@ TEST(ProbDistributionsMultiStudentTCholesky, ErrorSize0) {
   using stan::math::multi_student_t_cholesky_lpdf;
   using stan::math::multi_student_t_cholesky_rng;
   using std::vector;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   Matrix<double, Dynamic, 1> y_empty(0, 1);
   Matrix<double, Dynamic, 1> y(3, 1);
   y << 2.0, -2.0, 11.0;
@@ -293,7 +293,7 @@ TEST(ProbDistributionsMultiStudentTCholesky, ErrorSize2) {
   using stan::math::multi_student_t_cholesky_lpdf;
   using stan::math::multi_student_t_cholesky_rng;
   using std::vector;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   Matrix<double, Dynamic, 1> y(3, 1);
   y << 2.0, -2.0, 11.0;
   Matrix<double, Dynamic, 1> mu(3, 1);
@@ -315,7 +315,7 @@ TEST(ProbDistributionsMultiStudentTCholesky, ErrorSize3) {
   using stan::math::multi_student_t_cholesky_lpdf;
   using stan::math::multi_student_t_cholesky_rng;
   using std::vector;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   Matrix<double, Dynamic, 1> y(3, 1);
   y << 2.0, -2.0, 11.0;
   Matrix<double, Dynamic, 1> mu(3, 1);
@@ -356,7 +356,7 @@ TEST(ProbDistributionsMultiStudentTCholesky,
   using Eigen::Matrix;
   using stan::math::multi_student_t_cholesky_rng;
   using std::vector;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   Matrix<double, Dynamic, 1> mu(3, 1);
   mu << 2.0, 3.0, 11.0;
 
@@ -401,7 +401,7 @@ TEST(ProbDistributionsMultiStudentTCholesky,
   using stan::math::multi_student_t_cholesky_lpdf;
   using stan::math::multi_student_t_cholesky_rng;
   using std::vector;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   Matrix<double, Dynamic, 1> mu(3, 1);
   mu << 2.0, 3.0, 11.0;
 
diff --git a/test/unit/math/prim/prob/multi_student_t_test.cpp b/test/unit/math/prim/prob/multi_student_t_test.cpp
index 5d0795e9e5e..d49f93cedf5 100644
--- a/test/unit/math/prim/prob/multi_student_t_test.cpp
+++ b/test/unit/math/prim/prob/multi_student_t_test.cpp
@@ -1,6 +1,6 @@
 #include <stan/math/prim.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <boost/math/distributions.hpp>
 #include <vector>
 #include <limits>
@@ -11,7 +11,7 @@ TEST(ProbDistributionsMultiStudentT, NotVectorized) {
   using stan::math::multi_student_t_log;
   using stan::math::multi_student_t_rng;
   using std::vector;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   Matrix<double, Dynamic, 1> y(3, 1);
   y << 2.0, -2.0, 11.0;
   Matrix<double, Dynamic, 1> mu(3, 1);
@@ -30,7 +30,7 @@ TEST(ProbDistributionsMultiStudentT, Vectorized) {
   using stan::math::multi_student_t_log;
   using stan::math::multi_student_t_rng;
   using std::vector;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   vector<Matrix<double, Dynamic, 1> > vec_y(2);
   vector<Matrix<double, 1, Dynamic> > vec_y_t(2);
   Matrix<double, Dynamic, 1> y(3);
@@ -100,7 +100,7 @@ TEST(ProbDistributionsMultiStudentT, Sigma) {
   using stan::math::multi_student_t_log;
   using stan::math::multi_student_t_rng;
   using std::vector;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   Matrix<double, Dynamic, 1> y(3, 1);
   y << 2.0, -2.0, 11.0;
   Matrix<double, Dynamic, 1> mu(3, 1);
@@ -123,7 +123,7 @@ TEST(ProbDistributionsMultiStudentT, Mu) {
   using stan::math::multi_student_t_log;
   using stan::math::multi_student_t_rng;
   using std::vector;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   Matrix<double, Dynamic, 1> y(3, 1);
   y << 2.0, -2.0, 11.0;
   Matrix<double, Dynamic, 1> mu(3, 1);
@@ -178,7 +178,7 @@ TEST(ProbDistributionsMultiStudentT, Nu) {
   using stan::math::multi_student_t_log;
   using stan::math::multi_student_t_rng;
   using std::vector;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   Matrix<double, Dynamic, 1> y(3, 1);
   y << 2.0, -2.0, 11.0;
   Matrix<double, Dynamic, 1> mu(3, 1);
@@ -234,7 +234,7 @@ TEST(ProbDistributionsMultiStudentT, ErrorSize2) {
   using stan::math::multi_student_t_log;
   using stan::math::multi_student_t_rng;
   using std::vector;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   Matrix<double, Dynamic, 1> y(3, 1);
   y << 2.0, -2.0, 11.0;
   Matrix<double, Dynamic, 1> mu(3, 1);
@@ -253,7 +253,7 @@ TEST(ProbDistributionsMultiStudentT, ErrorSize3) {
   using stan::math::multi_student_t_log;
   using stan::math::multi_student_t_rng;
   using std::vector;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   Matrix<double, Dynamic, 1> y(3, 1);
   y << 2.0, -2.0, 11.0;
   Matrix<double, Dynamic, 1> mu(3, 1);
@@ -272,7 +272,7 @@ TEST(ProbDistributionsMultiStudentT, ErrorSizeSigma) {
   using stan::math::multi_student_t_lpdf;
   using stan::math::multi_student_t_rng;
   using std::vector;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   Matrix<double, Dynamic, 1> y(3, 1);
   y << 2.0, -2.0, 11.0;
   Matrix<double, Dynamic, 1> mu(3, 1);
@@ -308,7 +308,7 @@ TEST(ProbDistributionsMultiStudentT, marginalOneChiSquareGoodnessFitTest) {
   using stan::math::multi_student_t_log;
   using stan::math::multi_student_t_rng;
   using std::vector;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   Matrix<double, Dynamic, 1> mu(3, 1);
   mu << 2.0, 3.0, 11.0;
 
@@ -351,7 +351,7 @@ TEST(ProbDistributionsMultiStudentT, marginalTwoChiSquareGoodnessFitTest) {
   using stan::math::multi_student_t_log;
   using stan::math::multi_student_t_rng;
   using std::vector;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   Matrix<double, Dynamic, 1> mu(3, 1);
   mu << 2.0, 3.0, 11.0;
 
diff --git a/test/unit/math/prim/prob/multinomial_logit_test.cpp b/test/unit/math/prim/prob/multinomial_logit_test.cpp
index ea587254e9a..dfefdd5ca3f 100644
--- a/test/unit/math/prim/prob/multinomial_logit_test.cpp
+++ b/test/unit/math/prim/prob/multinomial_logit_test.cpp
@@ -1,6 +1,6 @@
 #include <stan/math/prim.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <boost/math/distributions.hpp>
 #include <limits>
 #include <vector>
@@ -9,7 +9,7 @@ using Eigen::Dynamic;
 using Eigen::Matrix;
 
 TEST(ProbDistributionsMultinomialLogit, RNGZero) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   Matrix<double, Dynamic, 1> beta(3);
   beta << 1.3, 0.1, -2.6;
   // bug in 4.8.1: RNG does not allow a zero total count
@@ -22,7 +22,7 @@ TEST(ProbDistributionsMultinomialLogit, RNGZero) {
 }
 
 TEST(ProbDistributionsMultinomialLogit, RNGSize) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   Matrix<double, Dynamic, 1> beta(5);
   beta << log(0.3), log(0.1), log(0.2), log(0.2), log(0.2);
   std::vector<int> sample = stan::math::multinomial_logit_rng(beta, 10, rng);
@@ -106,7 +106,7 @@ TEST(ProbDistributionsMultinomialLogit, zeros) {
 }
 
 TEST(ProbDistributionsMultinomialLogit, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int M = 10;
   int trials = 1000;
   int N = M * trials;
diff --git a/test/unit/math/prim/prob/multinomial_test.cpp b/test/unit/math/prim/prob/multinomial_test.cpp
index 62e1dbd5a58..3003a557e86 100644
--- a/test/unit/math/prim/prob/multinomial_test.cpp
+++ b/test/unit/math/prim/prob/multinomial_test.cpp
@@ -1,6 +1,6 @@
 #include <stan/math/prim.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <boost/math/distributions.hpp>
 #include <limits>
 #include <vector>
@@ -8,7 +8,7 @@
 TEST(ProbDistributionsMultinomial, RNGZero) {
   using Eigen::Dynamic;
   using Eigen::Matrix;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   Matrix<double, Dynamic, 1> theta(3);
   theta << 0.3, 0.1, 0.6;
   // bug in 4.8.1: RNG does not allow a zero total count
@@ -23,7 +23,7 @@ TEST(ProbDistributionsMultinomial, RNGZero) {
 TEST(ProbDistributionsMultinomial, RNGSize) {
   using Eigen::Dynamic;
   using Eigen::Matrix;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   Matrix<double, Dynamic, 1> theta(5);
   // error in 2.1.0 due to overflow in binomial call due to division
   theta << 0.3, 0.1, 0.2, 0.2, 0.2;
@@ -124,7 +124,7 @@ TEST(ProbDistributionsMultinomial, zeros) {
 TEST(ProbDistributionsMultinomial, error_check) {
   using Eigen::Dynamic;
   using Eigen::Matrix;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
 
   Matrix<double, Dynamic, 1> theta(3);
   theta << 0.15, 0.45, 0.40;
@@ -138,7 +138,7 @@ TEST(ProbDistributionsMultinomial, error_check) {
 TEST(ProbDistributionsMultinomial, chiSquareGoodnessFitTest) {
   using Eigen::Dynamic;
   using Eigen::Matrix;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int M = 10;
   int trials = 1000;
   int N = M * trials;
diff --git a/test/unit/math/prim/prob/neg_binomial_2_ccdf_log_test.cpp b/test/unit/math/prim/prob/neg_binomial_2_ccdf_log_test.cpp
index 309a32c1bdc..3e2cc26f550 100644
--- a/test/unit/math/prim/prob/neg_binomial_2_ccdf_log_test.cpp
+++ b/test/unit/math/prim/prob/neg_binomial_2_ccdf_log_test.cpp
@@ -1,6 +1,6 @@
 #include <stan/math/prim.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <boost/math/distributions.hpp>
 
 TEST(ProbNegBinomial2, ccdf_log_matches_lccdf) {
diff --git a/test/unit/math/prim/prob/neg_binomial_2_log_test.cpp b/test/unit/math/prim/prob/neg_binomial_2_log_test.cpp
index 335ee820f8e..c6af68bff51 100644
--- a/test/unit/math/prim/prob/neg_binomial_2_log_test.cpp
+++ b/test/unit/math/prim/prob/neg_binomial_2_log_test.cpp
@@ -3,7 +3,7 @@
 #include <test/unit/math/prim/prob/NegativeBinomial2LogTestRig.hpp>
 #include <test/unit/math/expect_near_rel.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <boost/math/distributions.hpp>
 #include <algorithm>
 #include <limits>
@@ -17,7 +17,7 @@ TEST(ProbDistributionsNegativeBinomial2Log, errorCheck) {
 TEST(ProbDistributionsNegBinomial2Log, error_check) {
   using std::log;
 
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   EXPECT_NO_THROW(stan::math::neg_binomial_2_log_rng(6, 2, rng));
   EXPECT_NO_THROW(stan::math::neg_binomial_2_log_rng(-0.5, 1, rng));
   EXPECT_NO_THROW(stan::math::neg_binomial_2_log_rng(log(1e8), 1, rng));
@@ -71,7 +71,7 @@ TEST(ProbDistributionsNegBinomial2Log, error_check) {
 }
 
 TEST(ProbDistributionsNegBinomial2Log, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int N = 1000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
   boost::math::negative_binomial_distribution<> dist(1.1,
@@ -109,7 +109,7 @@ TEST(ProbDistributionsNegBinomial2Log, chiSquareGoodnessFitTest) {
 }
 
 TEST(ProbDistributionsNegBinomial2Log, chiSquareGoodnessFitTest2) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int N = 1000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
   boost::math::negative_binomial_distribution<> dist(0.6,
@@ -147,7 +147,7 @@ TEST(ProbDistributionsNegBinomial2Log, chiSquareGoodnessFitTest2) {
 }
 
 TEST(ProbDistributionsNegBinomial2Log, chiSquareGoodnessFitTest3) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int N = 1000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
   boost::math::negative_binomial_distribution<> dist(121,
diff --git a/test/unit/math/prim/prob/neg_binomial_2_test.cpp b/test/unit/math/prim/prob/neg_binomial_2_test.cpp
index 0c7d2f0c35a..09608a4f2e1 100644
--- a/test/unit/math/prim/prob/neg_binomial_2_test.cpp
+++ b/test/unit/math/prim/prob/neg_binomial_2_test.cpp
@@ -4,7 +4,7 @@
 #include <test/unit/math/prim/prob/VectorIntRNGTestRig.hpp>
 #include <test/unit/math/expect_near_rel.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <boost/math/distributions.hpp>
 #include <limits>
 #include <vector>
@@ -39,7 +39,7 @@ TEST(ProbDistributionsNegBinomial2, distributionCheck) {
 }
 
 TEST(ProbDistributionsNegBinomial2, error_check) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   EXPECT_NO_THROW(stan::math::neg_binomial_2_rng(6, 2, rng));
   EXPECT_NO_THROW(stan::math::neg_binomial_2_rng(0.5, 1, rng));
   EXPECT_NO_THROW(stan::math::neg_binomial_2_rng(1e8, 1, rng));
@@ -91,7 +91,7 @@ TEST(ProbDistributionsNegBinomial2, error_check) {
 }
 
 TEST(ProbDistributionsNegBinomial2, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int N = 1000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
   boost::math::negative_binomial_distribution<> dist(1.1, 1.1 / (1.1 + 2.4));
@@ -128,7 +128,7 @@ TEST(ProbDistributionsNegBinomial2, chiSquareGoodnessFitTest) {
 }
 
 TEST(ProbDistributionsNegBinomial2, chiSquareGoodnessFitTest2) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int N = 1000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
   boost::math::negative_binomial_distribution<> dist(0.6, 0.6 / (0.6 + 2.4));
@@ -165,7 +165,7 @@ TEST(ProbDistributionsNegBinomial2, chiSquareGoodnessFitTest2) {
 }
 
 TEST(ProbDistributionsNegBinomial2, chiSquareGoodnessFitTest3) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int N = 1000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
   boost::math::negative_binomial_distribution<> dist(30, 30 / (30 + 60.4));
@@ -202,7 +202,7 @@ TEST(ProbDistributionsNegBinomial2, chiSquareGoodnessFitTest3) {
 }
 
 TEST(ProbDistributionsNegBinomial2, chiSquareGoodnessFitTest4) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int N = 1000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
   boost::math::negative_binomial_distribution<> dist(80, 80 / (80 + 30.4));
diff --git a/test/unit/math/prim/prob/neg_binomial_test.cpp b/test/unit/math/prim/prob/neg_binomial_test.cpp
index addc39f5d40..0bd9d7fa434 100644
--- a/test/unit/math/prim/prob/neg_binomial_test.cpp
+++ b/test/unit/math/prim/prob/neg_binomial_test.cpp
@@ -2,7 +2,7 @@
 #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 <boost/random/mixmax.hpp>
 #include <boost/math/distributions.hpp>
 #include <limits>
 #include <vector>
@@ -37,7 +37,7 @@ TEST(ProbDistributionsNegativeBinomial, distributionCheck) {
 }
 
 TEST(ProbDistributionsNegBinomial, error_check) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   EXPECT_NO_THROW(stan::math::neg_binomial_rng(6, 2, rng));
   EXPECT_NO_THROW(stan::math::neg_binomial_rng(0.5, 1, rng));
   EXPECT_NO_THROW(stan::math::neg_binomial_rng(1e9, 1, rng));
@@ -82,7 +82,7 @@ void expected_bin_sizes(double* expect, const int K, const int N,
 }
 
 TEST(ProbDistributionsNegBinomial, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   double p = 0.6;
   double alpha = 5;
   double beta = p / (1 - p);
@@ -120,7 +120,7 @@ TEST(ProbDistributionsNegBinomial, chiSquareGoodnessFitTest) {
 }
 
 TEST(ProbDistributionsNegBinomial, chiSquareGoodnessFitTest2) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   double p = 0.8;
   double alpha = 2.4;
   double beta = p / (1 - p);
@@ -158,7 +158,7 @@ TEST(ProbDistributionsNegBinomial, chiSquareGoodnessFitTest2) {
 }
 
 TEST(ProbDistributionsNegBinomial, chiSquareGoodnessFitTest3) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   double p = 0.2;
   double alpha = 0.4;
   double beta = p / (1 - p);
diff --git a/test/unit/math/prim/prob/normal_test.cpp b/test/unit/math/prim/prob/normal_test.cpp
index 1e4a2860ede..1a740674a4e 100644
--- a/test/unit/math/prim/prob/normal_test.cpp
+++ b/test/unit/math/prim/prob/normal_test.cpp
@@ -2,7 +2,7 @@
 #include <test/unit/math/prim/prob/vector_rng_test_helper.hpp>
 #include <test/unit/math/prim/prob/util.hpp>
 #include <boost/math/distributions.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <gtest/gtest.h>
 #include <limits>
 #include <vector>
@@ -77,7 +77,7 @@ TEST(ProbDistributionsNormal, distributionTest) {
 }
 
 TEST(ProbDistributionsNormal, error_check) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   EXPECT_NO_THROW(stan::math::normal_rng(10.0, 2.0, rng));
 
   EXPECT_THROW(stan::math::normal_rng(10.0, -2.0, rng), std::domain_error);
@@ -90,7 +90,7 @@ TEST(ProbDistributionsNormal, error_check) {
 }
 
 TEST(ProbDistributionsNormal, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int N = 10000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
 
diff --git a/test/unit/math/prim/prob/ordered_logistic_test.cpp b/test/unit/math/prim/prob/ordered_logistic_test.cpp
index c036906c702..3d1a12e2f49 100644
--- a/test/unit/math/prim/prob/ordered_logistic_test.cpp
+++ b/test/unit/math/prim/prob/ordered_logistic_test.cpp
@@ -1,6 +1,6 @@
 #include <stan/math/prim.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <boost/math/distributions.hpp>
 #include <limits>
 #include <vector>
@@ -168,7 +168,7 @@ TEST(ProbDistributions, ordered_logistic) {
 }
 
 TEST(ProbDistributionOrderedLogistic, error_check) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   double inf = std::numeric_limits<double>::infinity();
   Eigen::VectorXd c(4);
   c << -2, 2.0, 5, 10;
@@ -188,7 +188,7 @@ TEST(ProbDistributionOrderedLogistic, error_check) {
 
 TEST(ProbDistributionOrderedLogistic, chiSquareGoodnessFitTest) {
   using stan::math::inv_logit;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int N = 10000;
   double eta = 1.0;
   Eigen::VectorXd theta(3);
diff --git a/test/unit/math/prim/prob/ordered_probit_test.cpp b/test/unit/math/prim/prob/ordered_probit_test.cpp
index 0e5d82d5af6..418d7ff5f43 100644
--- a/test/unit/math/prim/prob/ordered_probit_test.cpp
+++ b/test/unit/math/prim/prob/ordered_probit_test.cpp
@@ -1,6 +1,6 @@
 #include <stan/math/prim.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <boost/math/distributions.hpp>
 #include <limits>
 
@@ -140,7 +140,7 @@ TEST(ProbDistributions, ordered_probit) {
 }
 
 TEST(ProbDistributionOrderedProbit, error_check) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   double inf = std::numeric_limits<double>::infinity();
   Eigen::VectorXd c(4);
   c << -2, 2.0, 5, 10;
@@ -155,7 +155,7 @@ TEST(ProbDistributionOrderedProbit, error_check) {
 
 TEST(ProbDistributionOrderedProbit, chiSquareGoodnessFitTest) {
   using stan::math::Phi;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int N = 10000;
   double eta = 1.0;
   Eigen::VectorXd theta(3);
diff --git a/test/unit/math/prim/prob/pareto_test.cpp b/test/unit/math/prim/prob/pareto_test.cpp
index df1c57c64ca..55152d17f13 100644
--- a/test/unit/math/prim/prob/pareto_test.cpp
+++ b/test/unit/math/prim/prob/pareto_test.cpp
@@ -2,7 +2,7 @@
 #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/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <boost/math/distributions.hpp>
 #include <limits>
 #include <vector>
@@ -46,7 +46,7 @@ TEST(ProbDistributionsPareto, distributionTest) {
 }
 
 TEST(ProbDistributionsPareto, error_check) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   EXPECT_NO_THROW(stan::math::pareto_rng(2.0, 1.0, rng));
 
   EXPECT_THROW(stan::math::pareto_rng(2.0, -1.0, rng), std::domain_error);
@@ -59,7 +59,7 @@ TEST(ProbDistributionsPareto, error_check) {
 }
 
 TEST(ProbDistributionsPareto, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int N = 10000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
 
diff --git a/test/unit/math/prim/prob/pareto_type_2_test.cpp b/test/unit/math/prim/prob/pareto_type_2_test.cpp
index 5af75607cf6..36502a5fe5b 100644
--- a/test/unit/math/prim/prob/pareto_type_2_test.cpp
+++ b/test/unit/math/prim/prob/pareto_type_2_test.cpp
@@ -1,6 +1,6 @@
 #include <stan/math/prim.hpp>
 #include <test/unit/math/prim/prob/vector_rng_test_helper.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <boost/math/distributions.hpp>
 #include <gtest/gtest.h>
 #include <vector>
@@ -31,7 +31,7 @@ TEST(ProbDistributionsParetoType2, errorCheck) {
  * the distributions manually
  */
 TEST(ProbDistributionsParetoType2Mat, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int N = 10000;
 
   double mu = 3.0;
@@ -108,7 +108,7 @@ TEST(ProbDistributionsParetoType2Mat, chiSquareGoodnessFitTest) {
 }
 
 TEST(ProbDistributionsParetoType2Prim, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int N = 10000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
   boost::math::chi_squared mydist(K - 1);
diff --git a/test/unit/math/prim/prob/poisson_binomial_test.cpp b/test/unit/math/prim/prob/poisson_binomial_test.cpp
index d2533f9eaa6..75c53c05f57 100644
--- a/test/unit/math/prim/prob/poisson_binomial_test.cpp
+++ b/test/unit/math/prim/prob/poisson_binomial_test.cpp
@@ -1,6 +1,6 @@
 #include <stan/math/prim.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <boost/math/distributions.hpp>
 #include <limits>
 #include <vector>
@@ -188,7 +188,7 @@ TEST(ProbDistributionsPoissonBinomial,
 TEST(ProbDistributionsPoissonBinomial, chiSquareGoodnessFitTest) {
   using vec = Eigen::Matrix<double, Eigen::Dynamic, 1>;
 
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int N = 10000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
   boost::math::binomial_distribution<> dist(100, 0.6);
diff --git a/test/unit/math/prim/prob/poisson_log_test.cpp b/test/unit/math/prim/prob/poisson_log_test.cpp
index 0399d2bb712..f38c1e07644 100644
--- a/test/unit/math/prim/prob/poisson_log_test.cpp
+++ b/test/unit/math/prim/prob/poisson_log_test.cpp
@@ -3,7 +3,7 @@
 #include <test/unit/math/prim/prob/VectorIntRNGTestRig.hpp>
 #include <gtest/gtest.h>
 #include <boost/math/distributions.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <limits>
 #include <vector>
 
diff --git a/test/unit/math/prim/prob/poisson_test.cpp b/test/unit/math/prim/prob/poisson_test.cpp
index 1966cbe14e1..5eab436a6bc 100644
--- a/test/unit/math/prim/prob/poisson_test.cpp
+++ b/test/unit/math/prim/prob/poisson_test.cpp
@@ -2,7 +2,7 @@
 #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 <boost/random/mixmax.hpp>
 #include <boost/math/distributions.hpp>
 #include <limits>
 #include <vector>
@@ -36,7 +36,7 @@ TEST(ProbDistributionsPoisson, distributionCheck) {
 TEST(ProbDistributionsPoisson, error_check) {
   using std::log;
 
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   EXPECT_NO_THROW(stan::math::poisson_rng(6, rng));
 
   EXPECT_THROW(stan::math::poisson_rng(-6, rng), std::domain_error);
@@ -56,7 +56,7 @@ TEST(ProbDistributionsPoisson, error_check) {
 }
 
 TEST(ProbDistributionsPoisson, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int N = 1000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
   boost::math::poisson_distribution<> dist(5);
@@ -95,7 +95,7 @@ TEST(ProbDistributionsPoisson, chiSquareGoodnessFitTest) {
 TEST(ProbDistributionsPoisson, chiSquareGoodnessFitTest2) {
   using std::log;
 
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int N = 1000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
   boost::math::poisson_distribution<> dist(5);
diff --git a/test/unit/math/prim/prob/rayleigh_test.cpp b/test/unit/math/prim/prob/rayleigh_test.cpp
index 4ca66c9825f..bc03d9ebf43 100644
--- a/test/unit/math/prim/prob/rayleigh_test.cpp
+++ b/test/unit/math/prim/prob/rayleigh_test.cpp
@@ -3,7 +3,7 @@
 #include <test/unit/math/prim/prob/util.hpp>
 #include <gtest/gtest.h>
 #include <boost/math/distributions.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <limits>
 #include <vector>
 
@@ -43,14 +43,14 @@ TEST(ProbDistributionsRayleigh, distributionTest) {
 }
 
 TEST(ProbDistributionsRayleigh, error_check) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   EXPECT_NO_THROW(stan::math::rayleigh_rng(2.0, rng));
 
   EXPECT_THROW(stan::math::rayleigh_rng(-2.0, rng), std::domain_error);
 }
 
 TEST(ProbDistributionsRayleigh, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int N = 10000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
 
diff --git a/test/unit/math/prim/prob/scaled_inv_chi_square_test.cpp b/test/unit/math/prim/prob/scaled_inv_chi_square_test.cpp
index 0999d35b88c..2bb8237b81b 100644
--- a/test/unit/math/prim/prob/scaled_inv_chi_square_test.cpp
+++ b/test/unit/math/prim/prob/scaled_inv_chi_square_test.cpp
@@ -3,7 +3,7 @@
 #include <test/unit/math/prim/prob/util.hpp>
 #include <gtest/gtest.h>
 #include <boost/math/distributions.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <limits>
 #include <vector>
 
@@ -46,7 +46,7 @@ TEST(ProbDistributionsScaledInvChiSquare, distributionTest) {
 }
 
 TEST(ProbDistributionsScaledInvChiSquare, error_check) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   EXPECT_NO_THROW(stan::math::scaled_inv_chi_square_rng(2.0, 1.0, rng));
 
   EXPECT_THROW(stan::math::scaled_inv_chi_square_rng(-2.0, 1.0, rng),
@@ -62,7 +62,7 @@ TEST(ProbDistributionsScaledInvChiSquare, error_check) {
 }
 
 TEST(ProbDistributionsScaledInvChiSquare, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int N = 10000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
 
diff --git a/test/unit/math/prim/prob/skew_double_exponential_ccdf_log_test.cpp b/test/unit/math/prim/prob/skew_double_exponential_ccdf_log_test.cpp
index b324b1ca5bd..ae34fc64597 100644
--- a/test/unit/math/prim/prob/skew_double_exponential_ccdf_log_test.cpp
+++ b/test/unit/math/prim/prob/skew_double_exponential_ccdf_log_test.cpp
@@ -2,7 +2,7 @@
 #include <vector>
 #include <gtest/gtest.h>
 
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <boost/math/distributions.hpp>
 #include <boost/math/tools/promotion.hpp>
 
diff --git a/test/unit/math/prim/prob/skew_double_exponential_cdf_log_test.cpp b/test/unit/math/prim/prob/skew_double_exponential_cdf_log_test.cpp
index ab096906b49..b8ff4092e8d 100644
--- a/test/unit/math/prim/prob/skew_double_exponential_cdf_log_test.cpp
+++ b/test/unit/math/prim/prob/skew_double_exponential_cdf_log_test.cpp
@@ -2,7 +2,7 @@
 #include <vector>
 #include <gtest/gtest.h>
 
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <boost/math/distributions.hpp>
 #include <boost/math/tools/promotion.hpp>
 
diff --git a/test/unit/math/prim/prob/skew_double_exponential_test.cpp b/test/unit/math/prim/prob/skew_double_exponential_test.cpp
index 3c9b27bc0a1..55bd365f946 100644
--- a/test/unit/math/prim/prob/skew_double_exponential_test.cpp
+++ b/test/unit/math/prim/prob/skew_double_exponential_test.cpp
@@ -3,7 +3,7 @@
 #include <test/unit/math/prim/prob/util.hpp>
 #include <gtest/gtest.h>
 #include <boost/math/distributions.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <limits>
 #include <vector>
 
@@ -35,7 +35,7 @@ TEST(ProbDistributionsSkewedDoubleExponential, errorCheck) {
 }
 
 TEST(ProbDistributionsSkewedDoubleExponential, error_check) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   EXPECT_NO_THROW(stan::math::skew_double_exponential_rng(10.0, 2.0, .1, rng));
 
   EXPECT_THROW(stan::math::skew_double_exponential_rng(10.0, 2.0, -1.0, rng),
@@ -60,7 +60,7 @@ TEST(ProbDistributionsSkewedDoubleExponential, test_sampling_icdf) {
 }
 
 TEST(ProbDistributionsSkewedDoubleExponential, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int N = 10000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
 
diff --git a/test/unit/math/prim/prob/skew_normal_test.cpp b/test/unit/math/prim/prob/skew_normal_test.cpp
index 1fe061f3392..85e51794d2f 100644
--- a/test/unit/math/prim/prob/skew_normal_test.cpp
+++ b/test/unit/math/prim/prob/skew_normal_test.cpp
@@ -2,7 +2,7 @@
 #include <test/unit/math/prim/prob/vector_rng_test_helper.hpp>
 #include <test/unit/math/prim/prob/util.hpp>
 #include <boost/math/distributions.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <gtest/gtest.h>
 #include <limits>
 #include <vector>
@@ -47,7 +47,7 @@ TEST(ProbDistributionsSkewNormal, distributionTest) {
 }
 
 TEST(ProbDistributionsSkewNormal, error_check) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   EXPECT_NO_THROW(stan::math::skew_normal_rng(10.0, 2.0, 1.0, rng));
 
   EXPECT_THROW(stan::math::skew_normal_rng(10.0, -2.0, 1.0, rng),
@@ -61,7 +61,7 @@ TEST(ProbDistributionsSkewNormal, error_check) {
 }
 
 TEST(ProbDistributionsSkewNormal, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int N = 10000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
 
diff --git a/test/unit/math/prim/prob/std_normal_test.cpp b/test/unit/math/prim/prob/std_normal_test.cpp
index 3ed4aa6500f..5730e4f7144 100644
--- a/test/unit/math/prim/prob/std_normal_test.cpp
+++ b/test/unit/math/prim/prob/std_normal_test.cpp
@@ -2,7 +2,7 @@
 #include <test/unit/math/prim/prob/util.hpp>
 #include <test/unit/math/prim/prob/vector_rng_test_helper.hpp>
 #include <boost/math/distributions.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <gtest/gtest.h>
 #include <limits>
 #include <vector>
@@ -68,12 +68,12 @@ TEST(ProbDistributionsStdNormal, distributionTest) {
 }
 
 TEST(ProbDistributionsStdNormal, error_check) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   EXPECT_NO_THROW(stan::math::std_normal_rng(rng));
 }
 
 TEST(ProbDistributionsStdNormal, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int N = 10000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
 
diff --git a/test/unit/math/prim/prob/student_t_test.cpp b/test/unit/math/prim/prob/student_t_test.cpp
index 6d95ebbf347..26684b151c0 100644
--- a/test/unit/math/prim/prob/student_t_test.cpp
+++ b/test/unit/math/prim/prob/student_t_test.cpp
@@ -2,7 +2,7 @@
 #include <test/unit/math/prim/prob/vector_rng_test_helper.hpp>
 #include <test/unit/math/prim/prob/util.hpp>
 #include <boost/math/distributions.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <gtest/gtest.h>
 #include <limits>
 #include <vector>
@@ -48,7 +48,7 @@ TEST(ProbDistributionsStudentT, distributionTest) {
 }
 
 TEST(ProbDistributionsStudentT, error_check) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   EXPECT_NO_THROW(stan::math::student_t_rng(3.0, 2.0, 2.0, rng));
 
   EXPECT_THROW(stan::math::student_t_rng(3.0, 2.0, -2.0, rng),
@@ -67,7 +67,7 @@ TEST(ProbDistributionsStudentT, error_check) {
 }
 
 TEST(ProbDistributionsStudentT, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int N = 10000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
 
diff --git a/test/unit/math/prim/prob/uniform_test.cpp b/test/unit/math/prim/prob/uniform_test.cpp
index e0ff8592551..cc6ce8a543f 100644
--- a/test/unit/math/prim/prob/uniform_test.cpp
+++ b/test/unit/math/prim/prob/uniform_test.cpp
@@ -3,7 +3,7 @@
 #include <test/unit/math/prim/prob/util.hpp>
 #include <gtest/gtest.h>
 #include <boost/math/distributions.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <limits>
 #include <vector>
 
@@ -63,7 +63,7 @@ TEST(ProbDistributionsUniform, distributionTest) {
 }
 
 TEST(ProbDistributionsUniform, error_check) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   EXPECT_NO_THROW(stan::math::uniform_rng(1.0, 2.0, rng));
 
   EXPECT_THROW(
@@ -75,7 +75,7 @@ TEST(ProbDistributionsUniform, error_check) {
 }
 
 TEST(ProbDistributionsUniform, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int N = 10000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
 
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 8c49d65aa9b..f05d373b347 100644
--- a/test/unit/math/prim/prob/vector_rng_test_helper.hpp
+++ b/test/unit/math/prim/prob/vector_rng_test_helper.hpp
@@ -3,7 +3,7 @@
 
 #include <gtest/gtest.h>
 #include <boost/math/distributions.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <stan/math/prim.hpp>
 #include <test/unit/math/prim/meta/apply_template_permutations.hpp>
 #include <test/unit/math/prim/prob/util.hpp>
@@ -165,7 +165,7 @@ struct check_dist_throws {
   template <typename T_param1, typename T_param2, typename T_param3,
             typename T_rig>
   void operator()(const T_rig& rig) const {
-    boost::random::mt19937 rng;
+    boost::random::mixmax rng;
 
     T_param1 p1;
     T_param2 p2;
@@ -371,7 +371,7 @@ struct check_quantiles {
   template <typename T_param1, typename T_param2, typename T_param3,
             typename T_rig>
   void operator()(const T_rig& rig) const {
-    boost::random::mt19937 rng;
+    boost::random::mixmax rng;
     T_param1 p1;
     T_param2 p2;
     T_param3 p3;
@@ -531,7 +531,7 @@ struct check_counts {
   template <typename T_param1, typename T_param2, typename T_param3,
             typename T_rig>
   void operator()(const T_rig& rig) const {
-    boost::random::mt19937 rng;
+    boost::random::mixmax rng;
     T_param1 p1;
     T_param2 p2;
     T_param3 p3;
diff --git a/test/unit/math/prim/prob/von_mises_test.cpp b/test/unit/math/prim/prob/von_mises_test.cpp
index d5c98f2d176..5cb028e8262 100644
--- a/test/unit/math/prim/prob/von_mises_test.cpp
+++ b/test/unit/math/prim/prob/von_mises_test.cpp
@@ -2,7 +2,7 @@
 #include <test/unit/math/prim/prob/util.hpp>
 #include <test/unit/math/prim/prob/vector_rng_test_helper.hpp>
 #include <boost/math/distributions.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <gtest/gtest.h>
 #include <vector>
 
@@ -30,7 +30,7 @@ TEST(ProbDistributionsVonMises, errorCheck) {
  * the distributions manually
  */
 TEST(ProbDistributionsVonMises, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int N = 10000;
 
   std::vector<double> loc
@@ -75,7 +75,7 @@ TEST(ProbDistributionsVonMises, chiSquareGoodnessFitTest) {
 }
 
 TEST(ProbDistributionsVonMises, error_check) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   EXPECT_NO_THROW(stan::math::von_mises_rng(1.0, 2.0, rng));
   EXPECT_NO_THROW(stan::math::von_mises_rng(1.0, 0.0, rng));
 
@@ -90,7 +90,7 @@ TEST(ProbDistributionsVonMises, error_check) {
 }
 
 TEST(ProbDistributionsVonMises, chiSquareGoodnessFitTest1) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int N = 10000;
   int K = 80;
   boost::math::chi_squared mydist(K - 1);
@@ -138,7 +138,7 @@ TEST(ProbDistributionsVonMises, chiSquareGoodnessFitTest1) {
 }
 
 TEST(ProbDistributionsVonMises, chiSquareGoodnessFitTest2) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int N = 10000;
   int K = 80;
   boost::math::chi_squared mydist(K - 1);
@@ -181,7 +181,7 @@ TEST(ProbDistributionsVonMises, chiSquareGoodnessFitTest2) {
 }
 
 TEST(ProbDistributionsVonMises, chiSquareGoodnessFitTest3) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int N = 10000;
   int K = 80;
   boost::math::chi_squared mydist(K - 1);
@@ -226,7 +226,7 @@ TEST(ProbDistributionsVonMises, chiSquareGoodnessFitTest3) {
 }
 
 TEST(ProbDistributionsVonMises, chiSquareGoodnessFitTest4) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int N = 10000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
 
diff --git a/test/unit/math/prim/prob/weibull_test.cpp b/test/unit/math/prim/prob/weibull_test.cpp
index 19308aea778..3cf291f8cac 100644
--- a/test/unit/math/prim/prob/weibull_test.cpp
+++ b/test/unit/math/prim/prob/weibull_test.cpp
@@ -2,7 +2,7 @@
 #include <test/unit/math/prim/prob/vector_rng_test_helper.hpp>
 #include <test/unit/math/prim/prob/util.hpp>
 #include <boost/math/distributions.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <gtest/gtest.h>
 #include <limits>
 #include <vector>
@@ -46,7 +46,7 @@ TEST(ProbDistributionsWeibull, distributionTest) {
 }
 
 TEST(ProbDistributionsWeibull, error_check) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   EXPECT_NO_THROW(stan::math::weibull_rng(2.0, 3.0, rng));
 
   EXPECT_THROW(stan::math::weibull_rng(-2.0, 3.0, rng), std::domain_error);
@@ -57,7 +57,7 @@ TEST(ProbDistributionsWeibull, error_check) {
 }
 
 TEST(ProbDistributionsWeibull, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int N = 10000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
 
diff --git a/test/unit/math/prim/prob/wishart_cholesky_rng_test.cpp b/test/unit/math/prim/prob/wishart_cholesky_rng_test.cpp
index dab129f4ab9..50faeb6be29 100644
--- a/test/unit/math/prim/prob/wishart_cholesky_rng_test.cpp
+++ b/test/unit/math/prim/prob/wishart_cholesky_rng_test.cpp
@@ -1,6 +1,6 @@
 #include <stan/math/prim.hpp>
 #include <test/unit/math/prim/util.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <boost/math/distributions/chi_squared.hpp>
 #include <gtest/gtest.h>
 #include <stdexcept>
@@ -13,7 +13,7 @@ TEST(ProbDistributionsWishartCholesky, rng) {
 
   using stan::math::wishart_cholesky_rng;
 
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
 
   MatrixXd omega(3, 4);
   EXPECT_THROW(wishart_cholesky_rng(3.0, omega, rng), std::invalid_argument);
@@ -34,7 +34,7 @@ TEST(ProbDistributionsWishartCholesky, rng_pos_def) {
   using Eigen::MatrixXd;
   using stan::math::wishart_cholesky_rng;
 
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
 
   MatrixXd Sigma(2, 2);
   MatrixXd Sigma_non_pos_def(2, 2);
@@ -57,7 +57,7 @@ TEST(ProbDistributionsWishartCholesky, marginalTwoChiSquareGoodnessFitTest) {
   using stan::math::wishart_cholesky_rng;
   using std::log;
 
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   MatrixXd sigma(3, 3);
   sigma << 9.0, -3.0, 2.0, -3.0, 4.0, 0.0, 2.0, 0.0, 3.0;
   int N = 10000;
@@ -91,7 +91,7 @@ TEST(ProbDistributionsWishartCholesky, SpecialRNGTest) {
   using stan::math::multiply_lower_tri_self_transpose;
   using stan::math::wishart_cholesky_rng;
 
-  boost::random::mt19937 rng(1234);
+  boost::random::mixmax rng(1234);
 
   MatrixXd sigma(3, 3);
 
diff --git a/test/unit/math/prim/prob/wishart_cholesky_test.cpp b/test/unit/math/prim/prob/wishart_cholesky_test.cpp
index 3a9d6d4f769..26796b24c68 100644
--- a/test/unit/math/prim/prob/wishart_cholesky_test.cpp
+++ b/test/unit/math/prim/prob/wishart_cholesky_test.cpp
@@ -1,6 +1,6 @@
 #include <stan/math/prim.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <boost/math/special_functions/digamma.hpp>
 #include <boost/math/distributions.hpp>
 
diff --git a/test/unit/math/prim/prob/wishart_rng_test.cpp b/test/unit/math/prim/prob/wishart_rng_test.cpp
index 28da632500a..6dee02c5d27 100644
--- a/test/unit/math/prim/prob/wishart_rng_test.cpp
+++ b/test/unit/math/prim/prob/wishart_rng_test.cpp
@@ -1,6 +1,6 @@
 #include <stan/math/prim.hpp>
 #include <test/unit/math/prim/util.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <boost/math/distributions/chi_squared.hpp>
 #include <gtest/gtest.h>
 #include <stdexcept>
@@ -10,7 +10,7 @@ TEST(ProbDistributionsWishart, rng) {
   using Eigen::MatrixXd;
   using stan::math::wishart_rng;
 
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
 
   MatrixXd omega(3, 4);
   EXPECT_THROW(wishart_rng(3.0, omega, rng), std::invalid_argument);
@@ -28,7 +28,7 @@ TEST(ProbDistributionsWishart, rng_pos_def) {
   using Eigen::MatrixXd;
   using stan::math::wishart_rng;
 
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
 
   MatrixXd Sigma(2, 2);
   MatrixXd Sigma_non_pos_def(2, 2);
@@ -48,7 +48,7 @@ TEST(ProbDistributionsWishart, rng_symmetry) {
   using stan::test::unit::expect_symmetric;
   using stan::test::unit::spd_rng;
 
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   for (int k = 1; k < 20; ++k)
     for (double nu = k - 0.9; nu < k + 10; ++nu)
       for (int n = 0; n < 10; ++n)
@@ -63,7 +63,7 @@ TEST(ProbDistributionsWishart, marginalTwoChiSquareGoodnessFitTest) {
   using stan::math::wishart_rng;
   using std::log;
 
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   MatrixXd sigma(3, 3);
   sigma << 9.0, -3.0, 2.0, -3.0, 4.0, 0.0, 2.0, 0.0, 3.0;
   int N = 10000;
@@ -94,7 +94,7 @@ TEST(ProbDistributionsWishart, SpecialRNGTest) {
   using Eigen::VectorXd;
   using stan::math::wishart_rng;
 
-  boost::random::mt19937 rng(1234);
+  boost::random::mixmax rng(1234);
 
   MatrixXd sigma(3, 3);
 
diff --git a/test/unit/math/prim/prob/wishart_test.cpp b/test/unit/math/prim/prob/wishart_test.cpp
index 9e1b12d762d..b113314ee90 100644
--- a/test/unit/math/prim/prob/wishart_test.cpp
+++ b/test/unit/math/prim/prob/wishart_test.cpp
@@ -1,6 +1,6 @@
 #include <stan/math/prim.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <boost/math/special_functions/digamma.hpp>
 #include <boost/math/distributions.hpp>
 
diff --git a/test/unit/math/rev/fun/cholesky_decompose_test.cpp b/test/unit/math/rev/fun/cholesky_decompose_test.cpp
index 60e7f54aa35..dcd41275539 100644
--- a/test/unit/math/rev/fun/cholesky_decompose_test.cpp
+++ b/test/unit/math/rev/fun/cholesky_decompose_test.cpp
@@ -1,7 +1,7 @@
 #include <stan/math/rev.hpp>
 #include <gtest/gtest.h>
 #include <test/unit/math/rev/util.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <vector>
 
 template <typename T_x>
@@ -184,7 +184,7 @@ void test_gradients_simple(int size, double prec) {
 
   stan::math::welford_covar_estimator estimator(size);
 
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   for (int i = 0; i < 1000; ++i) {
     Eigen::VectorXd q(size);
     for (int j = 0; j < size; ++j)
@@ -228,7 +228,7 @@ void test_gp_grad(int mat_size, double prec) {
   test_vals[1] = 1.5;
   test_vals[2] = 10;
 
-  boost::random::mt19937 rng(2);
+  boost::random::mixmax rng(2);
 
   for (int i = 0; i < mat_size; ++i) {
     test_vec(i) = stan::math::normal_rng(0.0, 0.1, rng);
@@ -266,7 +266,7 @@ void test_chol_mult(int mat_size, double prec) {
   Eigen::VectorXd test_vec(mat_size);
   Eigen::VectorXd test_vals(vec_size);
 
-  boost::random::mt19937 rng(2);
+  boost::random::mixmax rng(2);
 
   for (int i = 0; i < test_vals.size(); ++i) {
     if (i < test_vec.size()) {
@@ -292,7 +292,7 @@ void test_chol_mult(int mat_size, double prec) {
 
 void test_simple_vec_mult(int size, double prec) {
   Eigen::VectorXd test_vec(size);
-  boost::random::mt19937 rng(2);
+  boost::random::mixmax rng(2);
 
   for (int i = 0; i < test_vec.size(); ++i)
     test_vec(i) = stan::math::normal_rng(0.0, 0.1, rng);
diff --git a/test/unit/math/rev/prob/lkj_corr_test.cpp b/test/unit/math/rev/prob/lkj_corr_test.cpp
index 9faa3024e8f..27bc61deff3 100644
--- a/test/unit/math/rev/prob/lkj_corr_test.cpp
+++ b/test/unit/math/rev/prob/lkj_corr_test.cpp
@@ -2,13 +2,13 @@
 #include <test/unit/math/rev/prob/lkj_corr_cholesky_test_functors.hpp>
 #include <test/unit/math/rev/prob/test_gradients.hpp>
 #include <test/unit/math/rev/util.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <boost/math/distributions.hpp>
 #include <gtest/gtest.h>
 
 TEST(ProbDistributionsLkjCorr, var) {
   using stan::math::var;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int K = 4;
   Eigen::Matrix<var, Eigen::Dynamic, Eigen::Dynamic> Sigma(K, K);
   Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> Sigma_d(K, K);
@@ -29,7 +29,7 @@ TEST(ProbDistributionsLkjCorr, var) {
 
 TEST(ProbDistributionsLkjCorrCholesky, var) {
   using stan::math::var;
-  boost::random::mt19937 rng;
+  boost::random::mixmax rng;
   int K = 4;
   Eigen::Matrix<var, Eigen::Dynamic, Eigen::Dynamic> Sigma(K, K);
   Sigma.setZero();
diff --git a/test/unit/util.hpp b/test/unit/util.hpp
index 6d13989d579..672165e1402 100644
--- a/test/unit/util.hpp
+++ b/test/unit/util.hpp
@@ -2,7 +2,7 @@
 #define TEST_UNIT_UTIL_HPP
 
 #include <stan/math/prim/fun/Eigen.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/mixmax.hpp>
 #include <boost/random/uniform_real_distribution.hpp>
 #include <boost/typeof/typeof.hpp>
 #include <test/unit/pretty_print_types.hpp>
@@ -214,7 +214,7 @@ namespace test {
 
 auto make_sparse_matrix_random(int rows, int cols) {
   using eigen_triplet = Eigen::Triplet<double>;
-  boost::mt19937 gen;
+  boost::random::mixmax gen;
   boost::random::uniform_real_distribution<double> dist(0.0, 1.0);
   std::vector<eigen_triplet> tripletList;
   for (int i = 0; i < rows; ++i) {