BUG: Fix MLE for Nakagami nakagami_gen.fit

doc/source/tutorial/stats/continuous_nakagami.rst

@@ -9,18 +9,124 @@ Generalization of the chi distribution. Shape parameter is :math:`\nu>0.` The su
 .. math::
    :nowrap:
 
-    \begin{eqnarray*} f\left(x;\nu\right) & = & \frac{2\nu^{\nu}}{\Gamma\left(\nu\right)}x^{2\nu-1}\exp\left(-\nu x^{2}\right)\\
+    \begin{eqnarray*}
+    f\left(x;\nu\right) & = & \frac{2\nu^{\nu}}{\Gamma\left(\nu\right)}x^{2\nu-1}\exp\left(-\nu x^{2}\right)\\
     F\left(x;\nu\right) & = & \frac{\gamma\left(\nu,\nu x^{2}\right)}{\Gamma\left(\nu\right)}\\
-    G\left(q;\nu\right) & = & \sqrt{\frac{1}{\nu}\gamma^{-1}\left(\nu,q{\Gamma\left(\nu\right)}\right)}\end{eqnarray*}
+    G\left(q;\nu\right) & = & \sqrt{\frac{1}{\nu}\gamma^{-1}\left(\nu,q{\Gamma\left(\nu\right)}\right)}
+    \end{eqnarray*}
 
 where :math:`\gamma` is the lower incomplete gamma function, :math:`\gamma\left(\nu, x\right) = \int_0^x t^{\nu-1} e^{-t} dt`.
 
 .. math::
    :nowrap:
 
-    \begin{eqnarray*} \mu & = & \frac{\Gamma\left(\nu+\frac{1}{2}\right)}{\sqrt{\nu}\Gamma\left(\nu\right)}\\
+    \begin{eqnarray*}
+    \mu & = & \frac{\Gamma\left(\nu+\frac{1}{2}\right)}{\sqrt{\nu}\Gamma\left(\nu\right)}\\
     \mu_{2} & = & \left[1-\mu^{2}\right]\\
     \gamma_{1} & = & \frac{\mu\left(1-4v\mu_{2}\right)}{2\nu\mu_{2}^{3/2}}\\
-    \gamma_{2} & = & \frac{-6\mu^{4}\nu+\left(8\nu-2\right)\mu^{2}-2\nu+1}{\nu\mu_{2}^{2}}\end{eqnarray*}
+    \gamma_{2} & = & \frac{-6\mu^{4}\nu+\left(8\nu-2\right)\mu^{2}-2\nu+1}{\nu\mu_{2}^{2}}
+    \end{eqnarray*}
 
 Implementation: `scipy.stats.nakagami`
+
+
+MLE of the Nakagami Distribution in SciPy (:code:`nakagami.fit`)
+----------------------------------------------------------------
+
+The probability density function of the :code:`nakagami` distribution in SciPy is
+
+.. math::
+   :nowrap:
+
+    \begin{equation}
+    f(x; \nu, \mu, \sigma) = 2 \frac{\nu^\nu}{ \sigma \Gamma(\nu)}\left(\frac{x-\mu}{\sigma}\right)^{2\nu - 1} \exp\left(-\nu \left(\frac{x-\mu}{\sigma}\right)^2 \right),\tag{1}
+    \end{equation}
+
+for :math:`x` such that :math:`\frac{x-\mu}{\sigma} \geq 0`, where :math:`\nu \geq \frac{1}{2}` is the shape parameter,
+:math:`\mu` is the location, and :math:`\sigma` is the scale.
+
+The log-likelihood function is therefore
+
+.. math::
+   :nowrap:
+
+    \begin{equation}
+    l(\nu, \mu, \sigma) = \sum_{i = 1}^{N} \log \left( 2 \frac{\nu^\nu}{ \sigma\Gamma(\nu)}\left(\frac{x_i-\mu}{\sigma}\right)^{2\nu - 1} \exp\left(-\nu \left(\frac{x_i-\mu}{\sigma}\right)^2 \right) \right),\tag{2}
+    \end{equation}
+
+which can be expanded as
+
+.. math::
+   :nowrap:
+
+    \begin{equation}
+    l(\nu, \mu, \sigma) = N \log(2) + N\nu \log(\nu) - N\log\left(\Gamma(\nu)\right) - 2N \nu \log(\sigma) + \left(2 \nu - 1 \right) \sum_{i=1}^N  \log(x_i - \mu) - \nu \sigma^{-2} \sum_{i=1}^N \left(x_i-\mu\right)^2, \tag{3}
+    \end{equation}
+
+Leaving supports constraints out, the first-order condition for optimality on the likelihood derivatives gives estimates of parameters:
+
+.. math::
+   :nowrap:
+
+    \begin{align}
+    \frac{\partial l}{\partial \nu}(\nu, \mu, \sigma) &= N\left(1 + \log(\nu) - \psi^{(0)}(\nu)\right) + 2 \sum_{i=1}^N \log \left( \frac{x_i - \mu}{\sigma} \right) - \sum_{i=1}^N \left( \frac{x_i - \mu}{\sigma} \right)^2  = 0
+    \text{,} \tag{4}\\
+    \frac{\partial l}{\partial \mu}(\nu, \mu, \sigma) &= (1 - 2 \nu) \sum_{i=1}^N \frac{1}{x_i-\mu} + \frac{2\nu}{\sigma^2} \sum_{i=1}^N x_i-\mu = 0
+    \text{, and} \tag{5}\\
+    \frac{\partial l}{\partial \sigma}(\nu, \mu, \sigma) &= -2 N \nu \frac{1}{\sigma} + 2 \nu \sigma^{-3} \sum_{i=1}^N \left(x_i-\mu\right)^2 = 0
+    \text{,}\tag{6}
+    \end{align}
+
+where :math:`\psi^{(0)}` is the polygamma function of order :math:`0`; i.e. :math:`\psi^{(0)}(\nu) = \frac{d}{d\nu} \log \Gamma(\nu)`.
+
+However, the support of the distribution is the values of :math:`x` for which :math:`\frac{x-\mu}{\sigma} \geq 0`, and this provides an additional constraint that
+
+.. math::
+   :nowrap:
+
+    \begin{equation}
+    \mu \leq \min_i{x_i}.\tag{7}
+    \end{equation}
+
+For :math:`\nu = \frac{1}{2}`, the partial derivative of the log-likelihood with respect to :math:`\mu` reduces to:
+
+.. math::
+   :nowrap:
+
+    \begin{equation}
+    \frac{\partial l}{\partial \mu}(\nu, \mu, \sigma) = {\sigma^2} \sum_{i=1}^N (x_i-\mu),
+    \end{equation}
+
+which is positive when the support constraint is satisfied. Because the partial derivative with respect to :math:`\mu`
+is positive, increasing :math:`\mu` increases the log-likelihood, and therefore the constraint is active at the maximum likelihood estimate for :math:`\mu`
+
+.. math::
+   :nowrap:
+
+    \begin{equation}
+    \mu = \min_i{x_i}, \quad \nu = \frac{1}{2}. \tag{8}
+    \end{equation}
+
+For :math:`\nu` sufficiently greater than :math:`\frac{1}{2}`, the likelihood equation :math:`\frac{\partial l}{\partial \mu}(\nu, \mu, \sigma)=0` has a solution, and this solution provides the maximum likelihood estimate for :math:`\mu`. In either case, however, the condition :math:`\mu = \min_i{x_i}` provides a reasonable initial guess for numerical optimization.
+
+Furthermore, the likelihood equation for :math:`\sigma` can be solved explicitly, and it provides the maximum likelihood estimate
+
+.. math::
+   :nowrap:
+
+    \begin{equation}
+    \sigma = \sqrt{ \frac{\sum_{i=1}^N \left(x_i-\mu\right)^2}{N}}. \tag{9}
+    \end{equation}
+
+Hence, the :code:`_fitstart` method for :code:`nakagami` uses
+
+.. math::
+   :nowrap:
+
+    \begin{align}
+    \mu_0 &= \min_i{x_i} \,
+    \text{and} \\
+    \sigma_0 &= \sqrt{ \frac{\sum_{i=1}^N \left(x_i-\mu_0\right)^2}{N}}
+    \end{align}
+
+as initial guesses for numerical optimization accordingly.

scipy/stats/_continuous_distns.py

@@ -5831,6 +5831,15 @@ def _stats(self, nu):
         g2 /= nu*mu2**2.0
         return mu, mu2, g1, g2
 
+    def _fitstart(self, data, args=None):
+        if args is None:
+            args = (1.0,) * self.numargs
+        # Analytical justified estimates
+        # see: https://docs.scipy.org/doc/scipy/reference/tutorial/stats/continuous_nakagami.html
+        loc = np.min(data)
+        scale = np.sqrt(np.sum((data - loc)**2) / len(data))
+        return args + (loc, scale)
+
 
 nakagami = nakagami_gen(a=0.0, name="nakagami")
 

scipy/stats/tests/test_distributions.py

@@ -1,7 +1,6 @@
 """
 Test functions for stats module
 """
-
 import warnings
 import re
 import sys
@@ -26,7 +25,7 @@
 from scipy.stats._distn_infrastructure import argsreduce
 import scipy.stats.distributions
 
-from scipy.special import xlogy
+from scipy.special import xlogy, polygamma
 from .test_continuous_basic import distcont, invdistcont
 from .test_discrete_basic import distdiscrete, invdistdiscrete
 from scipy.stats._continuous_distns import FitDataError
@@ -5356,6 +5355,56 @@ def test_sf_isf(self):
         x1 = stats.nakagami.isf(sf, nu)
         assert_allclose(x1, x0, rtol=1e-13)
 
+    @pytest.mark.parametrize('nu', [1.6, 2.5, 3.9])
+    @pytest.mark.parametrize('loc', [25.0, 10, 35])
+    @pytest.mark.parametrize('scale', [13, 5, 20])
+    def test_fit(self, nu, loc, scale):
+        # Regression test for gh-13396 (21/27 cases failed previously)
+        # The first tuple of the parameters' values is discussed in gh-10908
+        N = 100
+        samples = stats.nakagami.rvs(size=N, nu=nu, loc=loc,
+                                     scale=scale, random_state=1337)
+        nu_est, loc_est, scale_est = stats.nakagami.fit(samples)
+        assert_allclose(nu_est, nu, rtol=0.2)
+        assert_allclose(loc_est, loc, rtol=0.2)
+        assert_allclose(scale_est, scale, rtol=0.2)
+
+        def dlogl_dnu(nu, loc, scale):
+            return ((-2*nu + 1) * np.sum(1/(samples - loc))
+                   + 2*nu/scale**2 * np.sum(samples - loc))
+
+        def dlogl_dloc(nu, loc, scale):
+            return (N * (1 + np.log(nu) - polygamma(0, nu)) +
+                    2 * np.sum(np.log((samples - loc) / scale))
+                    - np.sum(((samples - loc) / scale)**2))
+
+        def dlogl_dscale(nu, loc, scale):
+            return (- 2 * N * nu / scale
+                    + 2 * nu / scale ** 3 * np.sum((samples - loc) ** 2))
+
+        assert_allclose(dlogl_dnu(nu_est, loc_est, scale_est), 0, atol=1e-3)
+        assert_allclose(dlogl_dloc(nu_est, loc_est, scale_est), 0, atol=1e-3)
+        assert_allclose(dlogl_dscale(nu_est, loc_est, scale_est), 0, atol=1e-3)
+
+    @pytest.mark.parametrize('loc', [25.0, 10, 35])
+    @pytest.mark.parametrize('scale', [13, 5, 20])
+    def test_fit_nu(self, loc, scale):
+        # For nu = 0.5, we have analytical values for
+        # the MLE of the loc and the scale
+        nu = 0.5
+        n = 100
+        samples = stats.nakagami.rvs(size=n, nu=nu, loc=loc,
+                                     scale=scale, random_state=1337)
+        nu_est, loc_est, scale_est = stats.nakagami.fit(samples, f0=nu)
+
+        # Analytical values
+        loc_theo = np.min(samples)
+        scale_theo = np.sqrt(np.mean((samples - loc_est) ** 2))
+
+        assert_allclose(nu_est, nu, rtol=1e-7)
+        assert_allclose(loc_est, loc_theo, rtol=1e-7)
+        assert_allclose(scale_est, scale_theo, rtol=1e-7)
+
 
 def test_rvs_no_size_warning():
     class rvs_no_size_gen(stats.rv_continuous):

scipy/stats/tests/test_fit.py

@@ -42,8 +42,8 @@
                    'gengamma', 'gennorm', 'genpareto', 'halfcauchy',
                    'invgamma', 'invweibull', 'johnsonsu',
                    'kappa3', 'ksone', 'kstwo', 'levy', 'levy_l',
-                   'levy_stable', 'loglaplace', 'lomax', 'mielke', 'ncf',
-                   'nct', 'ncx2', 'pareto', 'powerlognorm', 'powernorm',
+                   'levy_stable', 'loglaplace', 'lomax', 'mielke', 'nakagami',
+                   'ncf', 'nct', 'ncx2', 'pareto', 'powerlognorm', 'powernorm',
                    'skewcauchy', 't',
                    'trapezoid', 'triang', 'tukeylambda']