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

Code before:
        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):
        rng = np.random.default_rng(413911473)
        mean = rng.standard_normal(3)
        A = rng.standard_normal((3, 3))
        cov = np.dot(A, A.transpose())        

        X = multivariate_normal(mean, cov)


Review comment (suggestion):
Not necessary, but some would prefer:
```suggestion
        cov = A @ A.T
```
here and below.

Code after:
        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:
            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)