scikit-learn-bot
diff --git a/‎doc/modules/gaussian_process.rst
Copy file name to clipboardExpand all lines: doc/modules/gaussian_process.rst
+7-2 b/‎doc/modules/gaussian_process.rst
Copy file name to clipboardExpand all lines: doc/modules/gaussian_process.rst
+7-2
diff --git a/‎doc/whats_new/upcoming_changes/sklearn.gaussian_process/22227.enhancement.rst
Copy file name to clipboard
+1 b/‎doc/whats_new/upcoming_changes/sklearn.gaussian_process/22227.enhancement.rst
Copy file name to clipboard
+1
diff --git a/‎sklearn/gaussian_process/_gpc.py
Copy file name to clipboardExpand all lines: sklearn/gaussian_process/_gpc.py
+76-9 b/‎sklearn/gaussian_process/_gpc.py
Copy file name to clipboardExpand all lines: sklearn/gaussian_process/_gpc.py
+76-9
diff --git a/‎sklearn/gaussian_process/tests/test_gpc.py
Copy file name to clipboardExpand all lines: sklearn/gaussian_process/tests/test_gpc.py
+35 b/‎sklearn/gaussian_process/tests/test_gpc.py
Copy file name to clipboardExpand all lines: sklearn/gaussian_process/tests/test_gpc.py
+35
@@ -106,11 +106,11 @@ The :class:`GaussianProcessClassifier` implements Gaussian processes (GP) for
 classification purposes, more specifically for probabilistic classification,
 where test predictions take the form of class probabilities.
 GaussianProcessClassifier places a GP prior on a latent function :math:`f`,
-which is then squashed through a link function to obtain the probabilistic
+which is then squashed through a link function :math:`\pi` to obtain the probabilistic
 classification. The latent function :math:`f` is a so-called nuisance function,
 whose values are not observed and are not relevant by themselves.
 Its purpose is to allow a convenient formulation of the model, and :math:`f`
-is removed (integrated out) during prediction. GaussianProcessClassifier
+is removed (integrated out) during prediction. :class:`GaussianProcessClassifier`
 implements the logistic link function, for which the integral cannot be
 computed analytically but is easily approximated in the binary case.
 
@@ -134,6 +134,11 @@ that have been chosen randomly from the range of allowed values.
 If the initial hyperparameters should be kept fixed, `None` can be passed as
 optimizer.
 
+In some scenarios, information about the latent function :math:`f` is desired
+(i.e. the mean :math:`\bar{f_*}` and the variance :math:`\text{Var}[f_*]` described
+in Eqs. (3.21) and (3.24) of [RW2006]_). The :class:`GaussianProcessClassifier`
+provides access to these quantities via the `latent_mean_and_variance` method.
+
 :class:`GaussianProcessClassifier` supports multi-class classification
 by performing either one-versus-rest or one-versus-one based training and
 prediction.  In one-versus-rest, one binary Gaussian process classifier is
 
@@ -0,0 +1 @@
+- :class:`gaussian_process.GaussianProcessClassifier` now includes a `latent_mean_and_variance` method that exposes the mean and the variance of the latent function, :math:`f`, used in the Laplace approximation. By :user:`Miguel González Duque <miguelgondu>`
@@ -306,12 +306,9 @@ def predict_proba(self, X):
         """
         check_is_fitted(self)
 
-        # Based on Algorithm 3.2 of GPML
-        K_star = self.kernel_(self.X_train_, X)  # K_star =k(x_star)
-        f_star = K_star.T.dot(self.y_train_ - self.pi_)  # Line 4
-        v = solve(self.L_, self.W_sr_[:, np.newaxis] * K_star)  # Line 5
-        # Line 6 (compute np.diag(v.T.dot(v)) via einsum)
-        var_f_star = self.kernel_.diag(X) - np.einsum("ij,ij->j", v, v)
+        # Compute the mean and variance of the latent function
+        # (Lines 4-6 of Algorithm 3.2 of GPML)
+        latent_mean, latent_var = self.latent_mean_and_variance(X)
 
         # Line 7:
         # Approximate \int log(z) * N(z | f_star, var_f_star)
@@ -320,12 +317,12 @@ def predict_proba(self, X):
         # sigmoid by a linear combination of 5 error functions.
         # For information on how this integral can be computed see
         # blitiri.blogspot.de/2012/11/gaussian-integral-of-error-function.html
-        alpha = 1 / (2 * var_f_star)
-        gamma = LAMBDAS * f_star
+        alpha = 1 / (2 * latent_var)
+        gamma = LAMBDAS * latent_mean
         integrals = (
             np.sqrt(np.pi / alpha)
             * erf(gamma * np.sqrt(alpha / (alpha + LAMBDAS**2)))
-            / (2 * np.sqrt(var_f_star * 2 * np.pi))
+            / (2 * np.sqrt(latent_var * 2 * np.pi))
         )
         pi_star = (COEFS * integrals).sum(axis=0) + 0.5 * COEFS.sum()
 
@@ -410,6 +407,39 @@ def log_marginal_likelihood(
 
         return Z, d_Z
 
+    def latent_mean_and_variance(self, X):
+        """Compute the mean and variance of the latent function values.
+
+        Based on algorithm 3.2 of [RW2006]_, this function returns the latent
+        mean (Line 4) and variance (Line 6) of the Gaussian process
+        classification model.
+
+        Note that this function is only supported for binary classification.
+
+        Parameters
+        ----------
+        X : array-like of shape (n_samples, n_features) or list of object
+            Query points where the GP is evaluated for classification.
+
+        Returns
+        -------
+        latent_mean : array-like of shape (n_samples,)
+            Mean of the latent function values at the query points.
+
+        latent_var : array-like of shape (n_samples,)
+            Variance of the latent function values at the query points.
+        """
+        check_is_fitted(self)
+
+        # Based on Algorithm 3.2 of GPML
+        K_star = self.kernel_(self.X_train_, X)  # K_star =k(x_star)
+        latent_mean = K_star.T.dot(self.y_train_ - self.pi_)  # Line 4
+        v = solve(self.L_, self.W_sr_[:, np.newaxis] * K_star)  # Line 5
+        # Line 6 (compute np.diag(v.T.dot(v)) via einsum)
+        latent_var = self.kernel_.diag(X) - np.einsum("ij,ij->j", v, v)
+
+        return latent_mean, latent_var
+
     def _posterior_mode(self, K, return_temporaries=False):
         """Mode-finding for binary Laplace GPC and fixed kernel.
 
@@ -902,3 +932,40 @@ def log_marginal_likelihood(
                     "Obtained theta with shape %d."
                     % (n_dims, n_dims * self.classes_.shape[0], theta.shape[0])
                 )
+
+    def latent_mean_and_variance(self, X):
+        """Compute the mean and variance of the latent function.
+
+        Based on algorithm 3.2 of [RW2006]_, this function returns the latent
+        mean (Line 4) and variance (Line 6) of the Gaussian process
+        classification model.
+
+        Note that this function is only supported for binary classification.
+
+        Parameters
+        ----------
+        X : array-like of shape (n_samples, n_features) or list of object
+            Query points where the GP is evaluated for classification.
+
+        Returns
+        -------
+        latent_mean : array-like of shape (n_samples,)
+            Mean of the latent function values at the query points.
+
+        latent_var : array-like of shape (n_samples,)
+            Variance of the latent function values at the query points.
+        """
+        if self.n_classes_ > 2:
+            raise ValueError(
+                "Returning the mean and variance of the latent function f "
+                "is only supported for binary classification, received "
+                f"{self.n_classes_} classes."
+            )
+        check_is_fitted(self)
+
+        if self.kernel is None or self.kernel.requires_vector_input:
+            X = validate_data(self, X, ensure_2d=True, dtype="numeric", reset=False)
+        else:
+            X = validate_data(self, X, ensure_2d=False, dtype=None, reset=False)
+
+        return self.base_estimator_.latent_mean_and_variance(X)
@@ -283,3 +283,38 @@ def test_gpc_fit_error(params, error_type, err_msg):
     gpc = GaussianProcessClassifier(**params)
     with pytest.raises(error_type, match=err_msg):
         gpc.fit(X, y)
+
+
+@pytest.mark.parametrize("kernel", kernels)
+def test_gpc_latent_mean_and_variance_shape(kernel):
+    """Checks that the latent mean and variance have the right shape."""
+    gpc = GaussianProcessClassifier(kernel=kernel)
+    gpc.fit(X, y)
+
+    # Check that the latent mean and variance have the right shape
+    latent_mean, latent_variance = gpc.latent_mean_and_variance(X)
+    assert latent_mean.shape == (X.shape[0],)
+    assert latent_variance.shape == (X.shape[0],)
+
+
+def test_gpc_latent_mean_and_variance_complain_on_more_than_2_classes():
+    """Checks that the latent mean and variance have the right shape."""
+    gpc = GaussianProcessClassifier(kernel=RBF())
+    gpc.fit(X, y_mc)
+
+    # Check that the latent mean and variance have the right shape
+    with pytest.raises(
+        ValueError,
+        match="Returning the mean and variance of the latent function f "
+        "is only supported for binary classification",
+    ):
+        gpc.latent_mean_and_variance(X)
+
+
+def test_latent_mean_and_variance_works_on_structured_kernels():
+    X = ["A", "AB", "B"]
+    y = np.array([True, False, True])
+    kernel = MiniSeqKernel(baseline_similarity_bounds="fixed")
+    gpc = GaussianProcessClassifier(kernel=kernel).fit(X, y)
+
+    gpc.latent_mean_and_variance(X)
Original file line number	Diff line number	Diff line change
`@@ -0,0 +1 @@`
	`1`	+- :class:`gaussian_process.GaussianProcessClassifier` now includes a `latent_mean_and_variance` method that exposes the mean and the variance of the latent function, :math:`f`, used in the Laplace approximation. By :user:`Miguel González Duque <miguelgondu>`