Describe the workflow you want to enable

Gaussian Process Regression in sklearn comes with the affordance to return standard deviations of the predictions in the gpr.predict method. The Gaussian Process Classifier doesn't. The main reason why (I assume) is that we don't have direct access to the standard deviation of our class estimates. We do, however, have access to the standard deviation of the latent variable f, and giving users access to it may allow them to have access to uncertainty estimates of the predicted classes by e.g. ancestral sampling or other Monte Carlo methods.

In this issue, I propose we add a keyword argument return_std_of_f=False to the methods predict and predict_proba of sklearn.gaussian_processes.GaussianProcessClassifier. This would allow us to return the standard deviation of the latent function f, a value that is already being computed in the predict_proba method.

Describe your proposed solution

The solution is actually pretty simple, since the computations are already being done in the method gpc.predict_proba. Namely, This line corresponds to computing the variance of the latent function $f$'s prediction on the test points. This addition would only need to be

1. The same computation under an if for the gpc.predict method.
2. An if statement with different returns according to the return_std_of_f flag for the gpc.predict_proba method.

Describe alternatives you've considered, if relevant

No response

Additional context

No response