simplify rng_fn method and refactor tests

Author: Sayam753

pymc3/distributions/multivariate.py

@@ -249,19 +249,6 @@ def _distr_parameters_for_repr(self):
         return ["mu", "cov"]
 
 
-def safe_multivariate_t(nu, mu, cov, size=None, rng=None):
-    res = np.atleast_1d(
-        stats.multivariate_t(loc=mu, shape=cov, df=nu, allow_singular=True).rvs(
-            size=size, random_state=rng
-        )
-    )
-
-    if size is not None:
-        res = res.reshape(list(size) + [-1])
-
-    return res
-
-
 class MvStudentTRV(RandomVariable):
     name = "multivariate_studentt"
     ndim_supp = 1
@@ -285,25 +272,22 @@ def _shape_from_params(self, dist_params, rep_param_idx=1, param_shapes=None):
     @classmethod
     def rng_fn(cls, rng, nu, mu, cov, size):
 
-        if mu.ndim > 1 or cov.ndim > 2:
-            # Neither SciPy nor NumPy implement parameter broadcasting for
-            # multivariate normals (or many other multivariate distributions),
-            # so we have implement a quick and dirty one here
-            mu, cov = broadcast_params([mu, cov], cls.ndims_params[1:])
-            size = tuple(size or ())
+        # Don't reassign broadcasted cov, since MvNormal expects two dimensional cov only.
+        mu, _ = broadcast_params([mu, cov], cls.ndims_params[1:])
 
-            if size:
-                mu = np.broadcast_to(mu, size + mu.shape)
-                cov = np.broadcast_to(cov, size + cov.shape)
-
-            res = np.empty(mu.shape)
-            for idx in np.ndindex(mu.shape[:-1]):
-                m = mu[idx]
-                c = cov[idx]
-                res[idx] = safe_multivariate_t(nu, m, c, rng=rng)
-            return res
-        else:
-            return safe_multivariate_t(nu, mu, cov, size=size, rng=rng)
+        chi2_samples = rng.chisquare(nu, size=size)
+        # Add distribution shape to chi2 samples
+        chi2_samples = chi2_samples.reshape(chi2_samples.shape + (1,) * len(mu.shape))
+
+        mv_samples = pm.MvNormal.dist(
+            mu=np.zeros_like(mu), cov=cov, size=size, rng=aesara.shared(rng)
+        ).eval()
+
+        size = tuple(size or ())
+        if size:
+            mu = np.broadcast_to(mu, size + mu.shape)
+
+        return (mv_samples / np.sqrt(chi2_samples / nu)) + mu
 
 
 mv_studentt = MvStudentTRV()

pymc3/tests/test_distributions_random.py

@@ -894,6 +894,11 @@ class TestMvNormalTau(BaseTestDistribution):
 
 
 class TestMvStudentTCov(BaseTestDistribution):
+    def mvstudentt_rng_fn(self, size, nu, mu, cov, rng):
+        chi2_samples = rng.chisquare(nu, size=size)
+        mv_samples = rng.multivariate_normal(np.zeros_like(mu), cov, size=size)
+        return (mv_samples / np.sqrt(chi2_samples[:, None] / nu)) + mu
+
     pymc_dist = pm.MvStudentT
     pymc_dist_params = {
         "nu": 5,
@@ -908,11 +913,13 @@ class TestMvStudentTCov(BaseTestDistribution):
     sizes_to_check = [None, (1), (2, 3)]
     sizes_expected = [(2,), (1, 2), (2, 3, 2)]
     reference_dist_params = {
-        "df": 5,
-        "loc": np.array([1.0, 2.0]),
-        "shape": np.array([[2.0, 0.0], [0.0, 3.5]]),
+        "nu": 5,
+        "mu": np.array([1.0, 2.0]),
+        "cov": np.array([[2.0, 0.0], [0.0, 3.5]]),
     }
-    reference_dist = seeded_scipy_distribution_builder("multivariate_t")
+    reference_dist = lambda self: functools.partial(
+        self.mvstudentt_rng_fn, rng=self.get_random_state()
+    )
     tests_to_run = [
         "check_pymc_params_match_rv_op",
         "check_pymc_draws_match_reference",