Repository: scikit-learn/scikit-learn
File: sklearn/gaussian_process/tests/test_gpr.py (Python)

Code before:
                result.append(gpr.log_marginal_likelihood(kernel.theta))

            else:
                result.append(
                    [gpr.log_marginal_likelihood([theta]) for theta in kernel.theta]
                )
        if length_scales.ndim == 1:
            return np.stack(result)
        elif length_scales.ndim == 2:
            return np.stack(result).reshape(
                length_scales.shape[0], length_scales.shape[1]
            )

    lml_gradient = []
    length_scale_param_name = next(
        name for name in kernel.get_params() if name.endswith("length_scale")
    )
    for i, length_scale in enumerate(length_scales):
        kernel.set_params(**{length_scale_param_name: length_scale})
        lml_gradient.append(
            gpr.log_marginal_likelihood(kernel.theta, eval_gradient=True)[1][0]
        )

    lml_gradient_approx = derivative(
        evaluate_grad_at_length_scales, length_scales, maxiter=20
    ).df
    assert_almost_equal(np.stack(lml_gradient), lml_gradient_approx, 3)


@pytest.mark.parametrize("kernel", kernels)
def test_prior(kernel):
    # Test that GP prior has mean 0 and identical variances.
    gpr = GaussianProcessRegressor(kernel=kernel)

    y_mean, y_cov = gpr.predict(X, return_cov=True)

    assert_almost_equal(y_mean, 0, 5)
    if len(gpr.kernel.theta) > 1:
        # XXX: quite hacky, works only for current kernels
        assert_almost_equal(np.diag(y_cov), np.exp(kernel.theta[0]), 5)
    else:
        assert_almost_equal(np.diag(y_cov), 1, 5)


@pytest.mark.parametrize("kernel", kernels)
def test_sample_statistics(kernel):
    # Test that statistics of samples drawn from GP are correct.
    gpr = GaussianProcessRegressor(kernel=kernel).fit(X, y)

    y_mean, y_cov = gpr.predict(X2, return_cov=True)


Review comment (suggestion):
```suggestion
    derivative_results = derivative(
         evaluate_grad_at_length_scales,
         length_scales,
         maxiter=20,
    )
    # Check that the numerical estimation converged within tolerances.
    assert np.all(derivative_results.status == 0)
    lml_gradient_approx = derivative_results.df
```

Code after:
    # Test that GP prior has mean 0 and identical variances.
    gpr = GaussianProcessRegressor(kernel=kernel)

    y_mean, y_cov = gpr.predict(X, return_cov=True)

    assert_almost_equal(y_mean, 0, 5)
    if len(gpr.kernel.theta) > 1:
        # XXX: quite hacky, works only for current kernels
        assert_almost_equal(np.diag(y_cov), np.exp(kernel.theta[0]), 5)
    else:
        assert_almost_equal(np.diag(y_cov), 1, 5)


@pytest.mark.parametrize("kernel", kernels)
def test_sample_statistics(kernel):
    # Test that statistics of samples drawn from GP are correct.
    gpr = GaussianProcessRegressor(kernel=kernel).fit(X, y)

    y_mean, y_cov = gpr.predict(X2, return_cov=True)

    samples = gpr.sample_y(X2, 300000)

    # More digits accuracy would require many more samples
    assert_almost_equal(y_mean, np.mean(samples, 1), 1)
    assert_almost_equal(
        np.diag(y_cov) / np.diag(y_cov).max(),
        np.var(samples, 1) / np.diag(y_cov).max(),
        1,
    )


def test_no_optimizer():
    # Test that kernel parameters are unmodified when optimizer is None.
    kernel = RBF(1.0)
    gpr = GaussianProcessRegressor(kernel=kernel, optimizer=None).fit(X, y)
    assert np.exp(gpr.kernel_.theta) == 1.0


@pytest.mark.parametrize("kernel", kernels)
@pytest.mark.parametrize("target", [y, np.ones(X.shape[0], dtype=np.float64)])
def test_predict_cov_vs_std(kernel, target):
    if sys.maxsize <= 2**32:
        pytest.xfail("This test may fail on 32 bit Python")

    # Test that predicted std.-dev. is consistent with cov's diagonal.
    gpr = GaussianProcessRegressor(kernel=kernel).fit(X, y)
    y_mean, y_cov = gpr.predict(X2, return_cov=True)
    y_mean, y_std = gpr.predict(X2, return_std=True)
    assert_almost_equal(np.sqrt(np.diag(y_cov)), y_std)