Fix MvStudentT.random

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):
         """

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(