Repository: scipy/scipy
Pull request: MAINT: stats.dirichlet: fix interface inconsistency
#### Reference issue
Closes gh-6006
gh-4984

#### What does this implement/fix?
The `dirichlet` distribution interface is inconsistent with other multivariate distributions and even itself: the `pdf` method expects the transpose of what the `rvs` method produces. This PR introduces `multivariate_beta`, which is the same distribution without this inconsistency.

The plan is to deprecate `dirichlet` in favor of `multivariate_beta`. (@h-vetinari can you help with this, either by making a PR against my branch or a follow-up to this one? I can review it.) After the deprecation cycle, we can make `dirichlet` an alias for `multivariate_beta` if desired.

#### Additional information
There are many other things that could be improved about the distribution and its documentation (e.g. `pdf` input `x` can only be 2D). This PR does not fix all of them, but it does fix a defect reported as early as gh-4984. Let's get this messy stuff out of the way, then those other things can be cleaned up in future PRs.

The three commits are pretty clean. I'd suggest reviewing them separately and in order.

LMK if I should address gh-6474 the same sort of way... maybe `inv_wishart`?

Problem statement:
Dirichlet doesn't accept its own random variates as input to pdf
```

In [1]: import scipy.stats as ss

In [2]: d = ss.dirichlet([1,3,4])

In [3]: d.pdf(d.rvs())
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-3-a6f3762e83cf> in <module>()
----> 1 d.pdf(d.rvs())

/opt/local/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/scipy/stats/_multivariate.py in pdf(self, x)
   1365
   1366     def pdf(self, x):
-> 1367         return self._dist.pdf(x, self.alpha)
   1368
   1369     def mean(self):

/opt/local/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/scipy/stats/_multivariate.py in pdf(self, x, alpha)
   1260         """
   1261         alpha = _dirichlet_check_parameters(alpha)
-> 1262         x = _dirichlet_check_input(alpha, x)
   1263
   1264         out = np.exp(self._logpdf(x, alpha))

/opt/local/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/scipy/stats/_multivariate.py in _dirichlet_check_input(alpha, x)
   1079                          "of entries as, or one entry fewer than, "
   1080                          "parameter vector 'a', but alpha.shape = %s "
-> 1081                          "and x.shape = %s." % (alpha.shape, x.shape))
   1082
   1083     if x.shape[0] != alpha.shape[0]:

ValueError: Vector 'x' must have either the same number of entries as, or one entry fewer than, parameter vector 'a', but alpha.shape = (3,) and x.shape = (1, 3).

In [4]: n = ss.norm()

In [5]: n.pdf(n.rvs())
Out[5]: 0.16964393208168185
```

Hints from issue discussion:
Every distribution should pass this test:

```
for size in [None, 1, 2]:
    d.pdf(d.rvs(size=size)).shape == (() if size is None else (size,))
```

Related discussion, starting from https://github.com/scipy/scipy/pull/4984#discussion_r32984650

@ev-br Thanks for the link.  Sorry for not searching properly before my submission.  I do think this change and my other recent bug report should be made to keep `dirichlet` consistent with `scipy.stats.multidimensional_normal` and `numpy.random`.

Agreed, this would be nice. If you can find a way of fixing it in a backwards-compatible way, that is :-). 

@ev-br How can changing the behaviour be backwards compatible?  In my own code, I did this:

```
class dirichlet_replacement(ss._multivariate.dirichlet_frozen):

    def rvs(self, size=None, random_state=None):
        # This somehow fixes the behaviour of rvs.
        return super().rvs(size=size, random_state=random_state)

    def pdf(self, x):
        if len(x.shape) == 2:
            return super().pdf(x.T)
        return super().pdf(x)
```

I think you're going to have to at some point have a deprecation cycle if necessary and then change the behaviour.

> I think you're going to have to at some point have a deprecation cycle 

Exactly! The question is only the details of how to do it in a least disruptive way. Help is appreciated :-)

@ev-br How can I help?  I have no experience with your project, so I don't know how you've done deprecation cycles in the past.

Generally, I think we'd need to either add a keyword (https://github.com/scipy/scipy/pull/4984#discussion_r33011011 with default value being the old behavior) or a new class (`dirichlet2` for a want of better name), make the old one raise a `DeprecationWarning` for at least a release or two, then remove the old one and rename the new one. 

See, it's a prolonged affair :-). But the sooner it's started, the better IMO. 

My preference would be to add a new class. This can also get other enhancements later on (stick breaking dirichlet, possibility of having some concentration values being zero etc). 

Would you be tempted to send a PR? 

@ev-br Sure, I'll look into it.  Would you be able to help me if I have questions?

Commit message: Apply suggestions from code review

Co-authored-by: Tirth Patel <tirthasheshpatel@gmail.com>

Patch to review:
diff --git a/scipy/stats/__init__.py b/scipy/stats/__init__.py
index feedc62099e8..c9a7fd720144 100644
--- a/scipy/stats/__init__.py
+++ b/scipy/stats/__init__.py
@@ -160,7 +160,8 @@
 
    multivariate_normal    -- Multivariate normal distribution
    matrix_normal          -- Matrix normal distribution
-   dirichlet              -- Dirichlet
+   multivariate_beta      -- Multivariate beta distribution (Dirichlet)
+   dirichlet              -- Dirichlet (deprecated, use `multivariate_beta` instead)
    wishart                -- Wishart
    invwishart             -- Inverse Wishart
    multinomial            -- Multinomial distribution
diff --git a/scipy/stats/_multivariate.py b/scipy/stats/_multivariate.py
index 23642ffea2af..6d8f1fe476ee 100644
--- a/scipy/stats/_multivariate.py
+++ b/scipy/stats/_multivariate.py
@@ -17,6 +17,7 @@
 
 __all__ = ['multivariate_normal',
            'matrix_normal',
+           'multivariate_beta',
            'dirichlet',
            'wishart',
            'invwishart',
@@ -1230,7 +1231,6 @@ def _dirichlet_check_parameters(alpha):
 
 
 def _dirichlet_check_input(alpha, x):
-    x = np.asarray(x)
 
     if x.shape[0] + 1 != alpha.shape[0] and x.shape[0] != alpha.shape[0]:
         raise ValueError("Vector 'x' must have either the same number "
@@ -1294,8 +1294,18 @@ def _lnB(alpha):
     return np.sum(gammaln(alpha)) - gammaln(np.sum(alpha))
 
 
+_dirichlet_depr_message = (
+"""
+`dirichlet` is deprecated due to an interface inconsistency: compared to
+other distributions, methods `pdf` and `logpdf` expect the transpose of the
+input `x`. Please use `multivariate_beta`, which corrects this inconsistency.
+In SciPy 1.11.0, this deprecation warning will be removed and `dirichlet` will
+become an alias for `multivariate_beta`.
+""")
+
+
 class dirichlet_gen(multi_rv_generic):
-    r"""A Dirichlet random variable.
+    r"""A Dirichlet random variable (deprecated, use `multivariate_beta` instead).
 
     The ``alpha`` keyword specifies the concentration parameters of the
     distribution.
@@ -1322,6 +1332,10 @@ class dirichlet_gen(multi_rv_generic):
     %(_dirichlet_doc_default_callparams)s
     %(_doc_random_state)s
 
+    See Also
+    --------
+    multivariate_beta
+
     Notes
     -----
     Each :math:`\alpha` entry must be positive. The distribution has only
@@ -1351,9 +1365,11 @@ class dirichlet_gen(multi_rv_generic):
     concentration parameters and :math:`K` is the dimension of the space
     where :math:`x` takes values.
 
-    Note that the dirichlet interface is somewhat inconsistent.
+    Note that the `dirichlet` interface is somewhat inconsistent.
     The array returned by the rvs function is transposed
     with respect to the format expected by the pdf and logpdf.
+    For a consistent interface to the same distribution, use
+    `multivariate_beta`.
 
     Examples
     --------
@@ -1403,11 +1419,17 @@ def __init__(self, seed=None):
         super().__init__(seed)
         self.__doc__ = doccer.docformat(self.__doc__, dirichlet_docdict_params)
 
+    @np.deprecate(old_name="dirichlet", new_name="multivariate_beta",
+                  message=_dirichlet_depr_message)
     def __call__(self, alpha, seed=None):
         return dirichlet_frozen(alpha, seed=seed)
 
+    def _check_input(self, alpha, x):
+        x = np.asarray(x)
+        return _dirichlet_check_input(alpha, x)
+
     def _logpdf(self, x, alpha):
-        """Log of the Dirichlet probability density function.
+        """Log of the multivariate beta (Dirichlet) PDF.
 
         Parameters
         ----------
@@ -1426,7 +1448,7 @@ def _logpdf(self, x, alpha):
         return - lnB + np.sum((xlogy(alpha - 1, x.T)).T, 0)
 
     def logpdf(self, x, alpha):
-        """Log of the Dirichlet probability density function.
+        """Log of the multivariate beta (Dirichlet) PDF.
 
         Parameters
         ----------
@@ -1441,13 +1463,13 @@ def logpdf(self, x, alpha):
 
         """
         alpha = _dirichlet_check_parameters(alpha)
-        x = _dirichlet_check_input(alpha, x)
+        x = self._check_input(alpha, x)
 
         out = self._logpdf(x, alpha)
         return _squeeze_output(out)
 
     def pdf(self, x, alpha):
-        """The Dirichlet probability density function.
+        """The multivariate beta (Dirichlet) probability density function.
 
         Parameters
         ----------
@@ -1462,13 +1484,13 @@ def pdf(self, x, alpha):
 
         """
         alpha = _dirichlet_check_parameters(alpha)
-        x = _dirichlet_check_input(alpha, x)
+        x = self._check_input(alpha, x)
 
         out = np.exp(self._logpdf(x, alpha))
         return _squeeze_output(out)
 
     def mean(self, alpha):
-        """Compute the mean of the dirichlet distribution.
+        """Mean of the multivariate beta (Dirichlet) distribution.
 
         Parameters
         ----------
@@ -1477,7 +1499,7 @@ def mean(self, alpha):
         Returns
         -------
         mu : ndarray or scalar
-            Mean of the Dirichlet distribution.
+            Mean of the multivariate beta (Dirichlet) distribution.
 
         """
         alpha = _dirichlet_check_parameters(alpha)
@@ -1486,7 +1508,7 @@ def mean(self, alpha):
         return _squeeze_output(out)
 
     def var(self, alpha):
-        """Compute the variance of the dirichlet distribution.
+        """Variance of the multivariate beta (Dirichlet) distribution.
 
         Parameters
         ----------
@@ -1506,7 +1528,8 @@ def var(self, alpha):
         return _squeeze_output(out)
 
     def entropy(self, alpha):
-        """Compute the differential entropy of the dirichlet distribution.
+        """
+        Differential entropy of the multivariate beta (Dirichlet) distribution.
 
         Parameters
         ----------
@@ -1515,7 +1538,7 @@ def entropy(self, alpha):
         Returns
         -------
         h : scalar
-            Entropy of the Dirichlet distribution
+            Entropy of the multivariate beta (Dirichlet) distribution
 
         """
 
@@ -1530,7 +1553,8 @@ def entropy(self, alpha):
         return _squeeze_output(out)
 
     def rvs(self, alpha, size=1, random_state=None):
-        """Draw random samples from a Dirichlet distribution.
+        """
+        Draw random samples from a multivariate beta (Dirichlet) distribution.
 
         Parameters
         ----------
@@ -1554,6 +1578,15 @@ def rvs(self, alpha, size=1, random_state=None):
 dirichlet = dirichlet_gen()
 
 
+# deprecate each public method of `dirichlet` but not `multivariate_beta`
+_dirichlet_method_names = ["entropy", "logpdf", "mean", "pdf",
+                           "rvs", "var"]
+for m in _dirichlet_method_names:
+    wrapper = np.deprecate(getattr(dirichlet, m), f"dirichlet.{m}",
+                           f"multivariate_beta.{m}", _dirichlet_depr_message)
+    setattr(dirichlet, m, wrapper)
+
+
 class dirichlet_frozen(multi_rv_frozen):
     def __init__(self, alpha, seed=None):
         self.alpha = _dirichlet_check_parameters(alpha)
@@ -1578,6 +1611,128 @@ def rvs(self, size=1, random_state=None):
         return self._dist.rvs(self.alpha, size, random_state)
 
 
+class multivariate_beta_gen(dirichlet_gen):
+    r"""A multivariate beta (Dirichlet) random variable.
+
+    The ``alpha`` keyword specifies the concentration parameters of the
+    distribution.
+
+    .. versionadded:: 1.9.0
+
+    Methods
+    -------
+    pdf(x, alpha)
+        Probability density function.
+    logpdf(x, alpha)
+        Log of the probability density function.
+    rvs(alpha, size=1, random_state=None)
+        Draw random samples from a multivariate beta distribution.
+    mean(alpha)
+        The mean of the multivariate beta distribution
+    var(alpha)
+        The variance of the multivariate beta distribution
+    entropy(alpha)
+        Compute the differential entropy of the multivariate beta distribution.
+
+    Parameters
+    ----------
+    %(_dirichlet_doc_default_callparams)s
+    %(_doc_random_state)s
+
+    See Also
+    --------
+    dirichlet
+
+    Notes
+    -----
+    Each :math:`\alpha` entry must be positive. The distribution has only
+    support on the simplex defined by
+
+    .. math::
+        \sum_{i=1}^{K} x_i = 1
+
+    where :math:`0 < x_i < 1`.
+
+    If the quantiles don't lie within the simplex, a `ValueError` is raised.
+
+    The probability density function for `multivariate_beta` is
+
+    .. math::
+
+        f(x) = \frac{1}{\mathrm{B}(\boldsymbol\alpha)}
+               \prod_{i=1}^K x_i^{\alpha_i - 1}
+
+    where
+
+    .. math::
+
+        \mathrm{B}(\boldsymbol\alpha) = \frac{\prod_{i=1}^K \Gamma(\alpha_i)}
+                                     {\Gamma\bigl(\sum_{i=1}^K \alpha_i\bigr)}
+
+    and :math:`\boldsymbol\alpha=(\alpha_1,\ldots,\alpha_K)`, the
+    concentration parameters and :math:`K` is the dimension of the space
+    where :math:`x` takes values.
+
+    Examples
+    --------
+    >>> from scipy.stats import multivariate_beta
+
+    Generate a multivariate_beta random variable
+
+    >>> quantiles = np.array([0.2, 0.2, 0.6])  # specify quantiles
+    >>> alpha = np.array([0.4, 5, 15])  # specify concentration parameters
+    >>> multivariate_beta.pdf(quantiles, alpha)
+    0.2843831684937255
+
+    The same PDF but following a log scale
+
+    >>> multivariate_beta.logpdf(quantiles, alpha)
+    -1.2574327653159187
+
+    Once we specify the multivariate beta distribution
+    we can then calculate quantities of interest
+
+    >>> multivariate_beta.mean(alpha)  # get the mean of the distribution
+    array([0.01960784, 0.24509804, 0.73529412])
+    >>> multivariate_beta.var(alpha) # get variance
+    array([0.00089829, 0.00864603, 0.00909517])
+    >>> multivariate_beta.entropy(alpha)  # calculate the differential entropy
+    -4.3280162474082715
+
+    We can also return random samples from the distribution
+
+    >>> multivariate_beta.rvs(alpha, size=1, random_state=1)
+    array([[0.00766178, 0.24670518, 0.74563305]])
+    >>> multivariate_beta.rvs(alpha, size=2, random_state=2)
+    array([[0.01639427, 0.1292273 , 0.85437844],
+           [0.00156917, 0.19033695, 0.80809388]])
+
+    Alternatively, the object may be called (as a function) to fix
+    concentration parameters, returning a "frozen" multivariate beta
+    random variable:
+
+    >>> rv = multivariate_beta(alpha)
+    >>> # Frozen object with the same methods but holding the given
+    >>> # concentration parameters fixed.
+
+    """
+    def __call__(self, alpha, seed=None):
+        return multivariate_beta_frozen(alpha, seed=seed)
+
+    def _check_input(self, alpha, x):
+        x = np.moveaxis(x, -1, 0)
+        return _dirichlet_check_input(alpha, x)
+
+
+class multivariate_beta_frozen(dirichlet_frozen):
+    def __init__(self, alpha, seed=None):
+        self.alpha = _dirichlet_check_parameters(alpha)
+        self._dist = multivariate_beta_gen(seed)
+
+
+multivariate_beta = multivariate_beta_gen()
+
+
 # Set frozen generator docstrings from corresponding docstrings in
 # multivariate_normal_gen and fill in default strings in class docstrings
 for name in ['logpdf', 'pdf', 'rvs', 'mean', 'var', 'entropy']:
diff --git a/scipy/stats/tests/test_multivariate.py b/scipy/stats/tests/test_multivariate.py
index ce6ee7301c09..94616850c619 100644
--- a/scipy/stats/tests/test_multivariate.py
+++ b/scipy/stats/tests/test_multivariate.py
@@ -25,7 +25,8 @@
                          random_correlation, unitary_group, dirichlet,
                          beta, wishart, multinomial, invwishart, chi2,
                          invgamma, norm, uniform, ks_2samp, kstest, binom,
-                         hypergeom, multivariate_t, cauchy, normaltest)
+                         hypergeom, multivariate_t, cauchy, normaltest,
+                         multivariate_beta)
 
 from scipy.integrate import romb
 from scipy.special import multigammaln
@@ -637,7 +638,8 @@ def test_moments(self):
                                                         N*num_cols,num_rows).T)
         assert_allclose(sample_rowcov, U, atol=0.1)
 
-class TestDirichlet:
+
+class DirichletTest:
 
     def test_frozen_dirichlet(self):
         np.random.seed(2846)
@@ -645,112 +647,116 @@ def test_frozen_dirichlet(self):
         n = np.random.randint(1, 32)
         alpha = np.random.uniform(10e-10, 100, n)
 
-        d = dirichlet(alpha)
+        d = self.dist(alpha)
 
-        assert_equal(d.var(), dirichlet.var(alpha))
-        assert_equal(d.mean(), dirichlet.mean(alpha))
-        assert_equal(d.entropy(), dirichlet.entropy(alpha))
+        assert_equal(d.var(), self.dist.var(alpha))
+        assert_equal(d.mean(), self.dist.mean(alpha))
+        assert_equal(d.entropy(), self.dist.entropy(alpha))
         num_tests = 10
         for i in range(num_tests):
             x = np.random.uniform(10e-10, 100, n)
             x /= np.sum(x)
-            assert_equal(d.pdf(x[:-1]), dirichlet.pdf(x[:-1], alpha))
-            assert_equal(d.logpdf(x[:-1]), dirichlet.logpdf(x[:-1], alpha))
+            assert_equal(d.pdf(x[:-1]), self.dist.pdf(x[:-1], alpha))
+            assert_equal(d.logpdf(x[:-1]), self.dist.logpdf(x[:-1], alpha))
 
     def test_numpy_rvs_shape_compatibility(self):
         np.random.seed(2846)
         alpha = np.array([1.0, 2.0, 3.0])
         x = np.random.dirichlet(alpha, size=7)
         assert_equal(x.shape, (7, 3))
-        assert_raises(ValueError, dirichlet.pdf, x, alpha)
-        assert_raises(ValueError, dirichlet.logpdf, x, alpha)
-        dirichlet.pdf(x.T, alpha)
-        dirichlet.pdf(x.T[:-1], alpha)
-        dirichlet.logpdf(x.T, alpha)
-        dirichlet.logpdf(x.T[:-1], alpha)
+        if self.dist == multivariate_beta:
+            x = x.T
+        assert_raises(ValueError, self.dist.pdf, x, alpha)
+        assert_raises(ValueError, self.dist.logpdf, x, alpha)
+        self.dist.pdf(x.T, alpha)
+        self.dist.pdf(x.T[:-1], alpha)
+        self.dist.logpdf(x.T, alpha)
+        self.dist.logpdf(x.T[:-1], alpha)
 
     def test_alpha_with_zeros(self):
         np.random.seed(2846)
         alpha = [1.0, 0.0, 3.0]
         # don't pass invalid alpha to np.random.dirichlet
         x = np.random.dirichlet(np.maximum(1e-9, alpha), size=7).T
-        assert_raises(ValueError, dirichlet.pdf, x, alpha)
-        assert_raises(ValueError, dirichlet.logpdf, x, alpha)
+        assert_raises(ValueError, self.dist.pdf, x, alpha)
+        assert_raises(ValueError, self.dist.logpdf, x, alpha)
 
     def test_alpha_with_negative_entries(self):
         np.random.seed(2846)
         alpha = [1.0, -2.0, 3.0]
         # don't pass invalid alpha to np.random.dirichlet
         x = np.random.dirichlet(np.maximum(1e-9, alpha), size=7).T
-        assert_raises(ValueError, dirichlet.pdf, x, alpha)
-        assert_raises(ValueError, dirichlet.logpdf, x, alpha)
+        assert_raises(ValueError, self.dist.pdf, x, alpha)
+        assert_raises(ValueError, self.dist.logpdf, x, alpha)
 
     def test_data_with_zeros(self):
         alpha = np.array([1.0, 2.0, 3.0, 4.0])
         x = np.array([0.1, 0.0, 0.2, 0.7])
-        dirichlet.pdf(x, alpha)
-        dirichlet.logpdf(x, alpha)
+        self.dist.pdf(x, alpha)
+        self.dist.logpdf(x, alpha)
         alpha = np.array([1.0, 1.0, 1.0, 1.0])
-        assert_almost_equal(dirichlet.pdf(x, alpha), 6)
-        assert_almost_equal(dirichlet.logpdf(x, alpha), np.log(6))
+        assert_almost_equal(self.dist.pdf(x, alpha), 6)
+        assert_almost_equal(self.dist.logpdf(x, alpha), np.log(6))
 
     def test_data_with_zeros_and_small_alpha(self):
         alpha = np.array([1.0, 0.5, 3.0, 4.0])
         x = np.array([0.1, 0.0, 0.2, 0.7])
-        assert_raises(ValueError, dirichlet.pdf, x, alpha)
-        assert_raises(ValueError, dirichlet.logpdf, x, alpha)
+        assert_raises(ValueError, self.dist.pdf, x, alpha)
+        assert_raises(ValueError, self.dist.logpdf, x, alpha)
 
     def test_data_with_negative_entries(self):
         alpha = np.array([1.0, 2.0, 3.0, 4.0])
         x = np.array([0.1, -0.1, 0.3, 0.7])
-        assert_raises(ValueError, dirichlet.pdf, x, alpha)
-        assert_raises(ValueError, dirichlet.logpdf, x, alpha)
+        assert_raises(ValueError, self.dist.pdf, x, alpha)
+        assert_raises(ValueError, self.dist.logpdf, x, alpha)
 
     def test_data_with_too_large_entries(self):
         alpha = np.array([1.0, 2.0, 3.0, 4.0])
         x = np.array([0.1, 1.1, 0.3, 0.7])
-        assert_raises(ValueError, dirichlet.pdf, x, alpha)
-        assert_raises(ValueError, dirichlet.logpdf, x, alpha)
+        assert_raises(ValueError, self.dist.pdf, x, alpha)
+        assert_raises(ValueError, self.dist.logpdf, x, alpha)
 
     def test_data_too_deep_c(self):
         alpha = np.array([1.0, 2.0, 3.0])
         x = np.full((2, 7, 7), 1 / 14)
-        assert_raises(ValueError, dirichlet.pdf, x, alpha)
-        assert_raises(ValueError, dirichlet.logpdf, x, alpha)
+        assert_raises(ValueError, self.dist.pdf, x, alpha)
+        assert_raises(ValueError, self.dist.logpdf, x, alpha)
 
     def test_alpha_too_deep(self):
         alpha = np.array([[1.0, 2.0], [3.0, 4.0]])
         x = np.full((2, 2, 7), 1 / 4)
-        assert_raises(ValueError, dirichlet.pdf, x, alpha)
-        assert_raises(ValueError, dirichlet.logpdf, x, alpha)
+        assert_raises(ValueError, self.dist.pdf, x, alpha)
+        assert_raises(ValueError, self.dist.logpdf, x, alpha)
 
     def test_alpha_correct_depth(self):
         alpha = np.array([1.0, 2.0, 3.0])
         x = np.full((3, 7), 1 / 3)
-        dirichlet.pdf(x, alpha)
-        dirichlet.logpdf(x, alpha)
+        if self.dist == multivariate_beta:
+            x = x.T
+        self.dist.pdf(x, alpha)
+        self.dist.logpdf(x, alpha)
 
     def test_non_simplex_data(self):
         alpha = np.array([1.0, 2.0, 3.0])
         x = np.full((3, 7), 1 / 2)
-        assert_raises(ValueError, dirichlet.pdf, x, alpha)
-        assert_raises(ValueError, dirichlet.logpdf, x, alpha)
+        assert_raises(ValueError, self.dist.pdf, x, alpha)
+        assert_raises(ValueError, self.dist.logpdf, x, alpha)
 
     def test_data_vector_too_short(self):
         alpha = np.array([1.0, 2.0, 3.0, 4.0])
         x = np.full((2, 7), 1 / 2)
-        assert_raises(ValueError, dirichlet.pdf, x, alpha)
-        assert_raises(ValueError, dirichlet.logpdf, x, alpha)
+        assert_raises(ValueError, self.dist.pdf, x, alpha)
+        assert_raises(ValueError, self.dist.logpdf, x, alpha)
 
     def test_data_vector_too_long(self):
         alpha = np.array([1.0, 2.0, 3.0, 4.0])
         x = np.full((5, 7), 1 / 5)
-        assert_raises(ValueError, dirichlet.pdf, x, alpha)
-        assert_raises(ValueError, dirichlet.logpdf, x, alpha)
+        assert_raises(ValueError, self.dist.pdf, x, alpha)
+        assert_raises(ValueError, self.dist.logpdf, x, alpha)
 
     def test_mean_and_var(self):
         alpha = np.array([1., 0.8, 0.2])
-        d = dirichlet(alpha)
+        d = self.dist(alpha)
 
         expected_var = [1. / 12., 0.08, 0.03]
         expected_mean = [0.5, 0.4, 0.1]
@@ -760,7 +766,7 @@ def test_mean_and_var(self):
 
     def test_scalar_values(self):
         alpha = np.array([0.2])
-        d = dirichlet(alpha)
+        d = self.dist(alpha)
 
         # For alpha of length 1, mean and var should be scalar instead of array
         assert_equal(d.mean().ndim, 0)
@@ -777,7 +783,7 @@ def test_K_and_K_minus_1_calls_equal(self):
         n = np.random.randint(1, 32)
         alpha = np.random.uniform(10e-10, 100, n)
 
-        d = dirichlet(alpha)
+        d = self.dist(alpha)
         num_tests = 10
         for i in range(num_tests):
             x = np.random.uniform(10e-10, 100, n)
@@ -790,7 +796,7 @@ def test_multiple_entry_calls(self):
 
         n = np.random.randint(1, 32)
         alpha = np.random.uniform(10e-10, 100, n)
-        d = dirichlet(alpha)
+        d = self.dist(alpha)
 
         num_tests = 10
         num_multiple = 5
@@ -803,7 +809,10 @@ def test_multiple_entry_calls(self):
                     xm = np.vstack((xm, x))
                 else:
                     xm = x
-            rm = d.pdf(xm.T)
+            if self.dist == multivariate_beta:
+                rm = d.pdf(xm)
+            else:
+                rm = d.pdf(xm.T)
             rs = None
             for xs in xm:
                 r = d.pdf(xs)
@@ -817,19 +826,54 @@ def test_2D_dirichlet_is_beta(self):
         np.random.seed(2846)
 
         alpha = np.random.uniform(10e-10, 100, 2)
-        d = dirichlet(alpha)
+        d = self.dist(alpha)
         b = beta(alpha[0], alpha[1])
 
         num_tests = 10
         for i in range(num_tests):
             x = np.random.uniform(10e-10, 100, 2)
             x /= np.sum(x)
-            assert_almost_equal(b.pdf(x), d.pdf([x]))
+            if self.dist == multivariate_beta:
+                assert_almost_equal(b.pdf(x), d.pdf(x[:, np.newaxis]))
+            else:
+                assert_almost_equal(b.pdf(x), d.pdf([x]))
 
         assert_almost_equal(b.mean(), d.mean()[0])
         assert_almost_equal(b.var(), d.var()[0])
 
 
+@pytest.mark.filterwarnings('ignore::DeprecationWarning')
+class TestDirichlet(DirichletTest):
+    dist = dirichlet
+
+
+def test_dirichlet_deprecation():
+    with pytest.deprecated_call():
+        dirichlet.rvs([1, 2, 3])
+
+    with pytest.deprecated_call():
+        dirichlet([1, 2, 3])
+
+
+class TestMultivariateBeta(DirichletTest):
+    dist = multivariate_beta
+
+    def test_accepts_random_variates(self):
+        # dirichlet.pdf does not accept output of dirichlet.rvs
+        # multivariate_beta does; see gh-6006
+        alpha = [1, 3, 4]
+        d = self.dist(alpha)
+
+        rvs = d.rvs()
+        pdf = d.pdf(rvs)  # no error
+        assert pdf == self.dist.pdf(rvs, alpha)
+
+        rvs2 = d.rvs(size=10)
+        pdf2 = d.pdf(rvs2)
+        for rvs_i, pdf_i in zip(rvs2, pdf2):
+            assert pdf_i == self.dist.pdf(rvs_i, alpha)
+
+
 def test_multivariate_normal_dimensions_mismatch():
     # Regression test for GH #3493. Check that setting up a PDF with a mean of
     # length M and a covariance matrix of size (N, N), where M != N, raises a
@@ -2165,7 +2209,7 @@ def test_random_state_property():
     scale[1, 0] = 0.5
     dists = [
         [multivariate_normal, ()],
-        [dirichlet, (np.array([1.]), )],
+        [multivariate_beta, (np.array([1.]), )],
         [wishart, (10, scale)],
         [invwishart, (10, scale)],
         [multinomial, (5, [0.5, 0.4, 0.1])],


Code review comments:
diff_hunk: @@ -160,7 +160,8 @@
 
    multivariate_normal    -- Multivariate normal distribution
    matrix_normal          -- Matrix normal distribution
-   dirichlet              -- Dirichlet
+   multivariate_beta      -- Multivariate beta distribution (Dirichlet)
+   dirichlet              -- Dirichlet (deprecated, use `multivariate_beta` instead)
path: scipy/stats/__init__.py
text: @author:
```scipy/stats/__init__.py:164:80: E501 line too long (84 > 79 characters)```
```suggestion
   dirichlet              -- Dirichlet (deprecated; see `multivariate_beta`)
```

diff_hunk: @@ -1294,8 +1294,18 @@ def _lnB(alpha):
     return np.sum(gammaln(alpha)) - gammaln(np.sum(alpha))
 
 
+_dirichlet_depr_message = (
+"""
path: scipy/stats/_multivariate.py
text: @author:
```suggestion
"""  # noqa
```