diff --git a/RELEASE-NOTES.md b/RELEASE-NOTES.md
index 17cb1d3dff..02c3ea1da9 100644
--- a/RELEASE-NOTES.md
+++ b/RELEASE-NOTES.md
@@ -14,6 +14,7 @@ This is the first release to support Python3.9 and to drop Python3.6.
 - In `sample_posterior_predictive` the `vars` kwarg was removed in favor of `var_names` (see [#4343](https://github.com/pymc-devs/pymc3/pull/4343)).
 - The notebook gallery has been moved to https://github.com/pymc-devs/pymc-examples (see [#4348](https://github.com/pymc-devs/pymc3/pull/4348)).
 - `math.logsumexp` now matches `scipy.special.logsumexp` when arrays contain infinite values (see [#4360](https://github.com/pymc-devs/pymc3/pull/4360)).
+- Fixed mathematical formulation in `MvStudentT` random method. (see [#4359](https://github.com/pymc-devs/pymc3/pull/4359))
 
 ## PyMC3 3.10.0 (7 December 2020)
 
diff --git a/pymc3/distributions/multivariate.py b/pymc3/distributions/multivariate.py
index 59851b2e0c..68d0ff4026 100755
--- a/pymc3/distributions/multivariate.py
+++ b/pymc3/distributions/multivariate.py
@@ -324,10 +324,10 @@ class MvStudentT(_QuadFormBase):
                1+\frac{1}{\nu}
                ({\mathbf x}-{\mu})^T
                {\Sigma}^{-1}({\mathbf x}-{\mu})
-             \right]^{(\nu+p)/2}}
+             \right]^{-(\nu+p)/2}}
 
     ========  =============================================
-    Support   :math:`x \in \mathbb{R}^k`
+    Support   :math:`x \in \mathbb{R}^p`
     Mean      :math:`\mu` if :math:`\nu > 1` else undefined
     Variance  :math:`\frac{\nu}{\mu-2}\Sigma`
                   if :math:`\nu>2` else undefined
@@ -393,8 +393,10 @@ def random(self, point=None, size=None):
 
             samples = dist.random(point, size)
 
-        chi2 = np.random.chisquare
-        return (np.sqrt(nu) * samples.T / chi2(nu, size)).T + mu
+        chi2_samples = np.random.chisquare(nu, size)
+        # Add distribution shape to chi2 samples
+        chi2_samples = chi2_samples.reshape(chi2_samples.shape + (1,) * len(self.shape))
+        return (samples / np.sqrt(chi2_samples / nu)) + mu
 
     def logp(self, value):
         """
diff --git a/pymc3/tests/test_distributions_random.py b/pymc3/tests/test_distributions_random.py
index 5b60aec91a..42144d7ea2 100644
--- a/pymc3/tests/test_distributions_random.py
+++ b/pymc3/tests/test_distributions_random.py
@@ -947,9 +947,9 @@ def ref_rand_evd(size, mu, evds, sigma):
 
     def test_mv_t(self):
         def ref_rand(size, nu, Sigma, mu):
-            normal = st.multivariate_normal.rvs(cov=Sigma, size=size).T
-            chi2 = st.chi2.rvs(df=nu, size=size)
-            return mu + np.sqrt(nu) * (normal / chi2).T
+            normal = st.multivariate_normal.rvs(cov=Sigma, size=size)
+            chi2 = st.chi2.rvs(df=nu, size=size)[..., None]
+            return mu + (normal / np.sqrt(chi2 / nu))
 
         for n in [2, 3]:
             pymc3_random(