Nakagami Loglikelihood tests: fit vs scipy.optimize.root

test_system_nakagami_fit.py

    @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 N * (1 + np.log(nu) - polygamma(0, nu)) + 2 * np.sum(np.log((samples - loc) / scale)) - np.sum(((samples - loc) / scale)**2)

        def dlogl_dloc(nu, loc, scale):
            return (1 - 2 * nu) * np.sum(1 / (samples - loc)) + 2 * nu / scale ** 2 * np.sum(samples) - loc

        def dlogl_dscale(nu, loc, scale):
            return - 2 * N * nu / scale + 2 * nu / scale ** 3 * np.sum((samples - loc) ** 2)

        def system(x):
            dnu = dlogl_dnu(*x)
            dloc = dlogl_dloc(*x)
            dscale = dlogl_dscale(*x)
            return (dnu, dloc, dscale)
        
        ll1 = np.sum(stats.nakagami.logpdf(samples, nu=nu_est, loc=loc_est, scale=scale_est))
        
        nu_est2, loc_est2, scale_est2 = root(system, (nu_est, loc_est, scale_est), method='lm').x
        
        ll2 = np.sum(stats.nakagami.logpdf(samples, nu=nu_est2, loc=loc_est2, scale=scale_est2))
        assert ll1 > ll2