Repository: scipy/scipy
File: scipy/stats/tests/test_multivariate.py (Python)

Code before:
        cov = np.dot(A, A.T)
        samples = multivariate_normal.rvs(mean=loc, cov=cov,
                                          size=100, random_state=rng)
        mean_free, cov_free = multivariate_normal.fit(samples)
        logp_free = multivariate_normal.logpdf(samples, mean=mean_free,
                                               cov=cov_free).sum()
        mean_fix, cov_fix = multivariate_normal.fit(samples, fix_cov=cov)
        assert_equal(mean_fix, np.mean(samples, axis=0))
        assert_equal(cov_fix, cov)
        logp_fix = multivariate_normal.logpdf(samples, mean=mean_fix,
                                              cov=cov_fix).sum()
        # test that fixed parameters result in lower likelihood than free
        # parameters
        assert logp_fix < logp_free
        # test that a small perturbation of the resulting parameters
        # has lower likelihood than the estimated parameters
        mean_perturbed = mean_fix + 1e-8 * rng.random(3)
        logp_perturbed = (multivariate_normal.logpdf(samples,
                                                     mean=mean_perturbed,
                                                     cov=cov_fix)
                                                     ).sum()
        assert logp_perturbed < logp_fix

    @pytest.mark.parametrize('X_ndim', [3])
    @pytest.mark.parametrize('dimensions', [[1], 
                                            [0, 1]])
    @pytest.mark.parametrize('frozen', [True, False])
    @pytest.mark.parametrize('cov_object', [True, False])
    def test_marginal_distribution(self, X_ndim, dimensions, frozen, cov_object):
        rng = np.random.default_rng(413911473)
        mean = rng.standard_normal(X_ndim)
        A = rng.standard_normal((X_ndim, X_ndim))
        cov = np.dot(A, A.transpose())
        
        if cov_object:
            cov = _covariance.CovViaPrecision(cov)

        x = np.random.standard_normal((4, len(dimensions))) # number of points at which to evaluate marginal PDF
        X = multivariate_normal(mean, cov)  # frozen distribution
        
        if frozen:
            Y = X.marginal(dimensions)
            res = Y.pdf(x)
        else: 
            Y = multivariate_normal.marginal(dimensions, mean, cov)
            res = Y.pdf(x)

        ref = marginal_pdf(X, X_ndim, dimensions, x)
        assert_allclose(ref, res)

    def test_marginal_dimensions(self):

Review comment (bug):
Valid negative indices are not currently testsed, so how about:
```suggestion
                                            [-1, 1]])
```
Then again, some tests fail when I make that change locally. So does `[2, 1]`, but `[0, -1]` and `[1, 2]` pass. Sorting the dimensions after standardizing them to be positive fixes the failing cases, but I don't think that's necessarily what we want to do. Please check what happens when `dimensions` are out of order in `gaussian_kde` - does that re-order the dimensions of the returned object in a consistent/usable way? If so, we can consider it a "feature", preserve that behavior here, and adjust the reference implementation to follow suit.

Also would be fine with me to keep this on the same line.


Code after:
        logp_perturbed = (multivariate_normal.logpdf(samples,
                                                     mean=mean_fix,
                                                     cov=cov_perturbed)
                                                     ).sum()
        assert logp_perturbed < logp_fix

    def test_fit_fix_cov(self):
        rng = np.random.default_rng(4385269356937404)
        loc = rng.random(3)
        A = rng.random((3, 3))
        cov = np.dot(A, A.T)
        samples = multivariate_normal.rvs(mean=loc, cov=cov,
                                          size=100, random_state=rng)
        mean_free, cov_free = multivariate_normal.fit(samples)
        logp_free = multivariate_normal.logpdf(samples, mean=mean_free,
                                               cov=cov_free).sum()
        mean_fix, cov_fix = multivariate_normal.fit(samples, fix_cov=cov)
        assert_equal(mean_fix, np.mean(samples, axis=0))
        assert_equal(cov_fix, cov)
        logp_fix = multivariate_normal.logpdf(samples, mean=mean_fix,
                                              cov=cov_fix).sum()
        # test that fixed parameters result in lower likelihood than free
        # parameters
        assert logp_fix < logp_free
        # test that a small perturbation of the resulting parameters
        # has lower likelihood than the estimated parameters
        mean_perturbed = mean_fix + 1e-8 * rng.random(3)
        logp_perturbed = (multivariate_normal.logpdf(samples,
                                                     mean=mean_perturbed,
                                                     cov=cov_fix)
                                                     ).sum()
        assert logp_perturbed < logp_fix

    @pytest.mark.parametrize('X_ndim', [3])
    @pytest.mark.parametrize('dimensions', [[1], [-1, 1]])
    @pytest.mark.parametrize('frozen', [True, False])
    @pytest.mark.parametrize('cov_object', [True, False])
    def test_marginal_distribution(self, X_ndim, dimensions, frozen, cov_object):
        rng = np.random.default_rng(413911473)
        mean = rng.standard_normal(X_ndim)
        A = rng.standard_normal((X_ndim, X_ndim))
        cov = A @ A.T
        
        if cov_object:
            cov = _covariance.CovViaPrecision(cov)

        # number of points at which to evaluate marginal PDF
        x = np.random.standard_normal((4, len(dimensions)))
        X = multivariate_normal(mean, cov)
        
        if frozen: