GaussianProcessRegressor: wrong std and cov results when n_features>1 and no y normalization

### Describe the bug

When `n_features > 1` and `normalization_y` is `False`, the `GaussianProcessRegressor.predict` seems to return bad std and cov results, as it doesn't consider the scale of the different features (while it seems to be ok when `n_features > 1` and `normalization_y` is `True`).

By taking a look at the code, we can see that `GaussianProcessRegressor.predict` uses the `_y_train_std` attribute to compute the variance and covariance but this attribute is set to `ones(n_features)` when `normalize_y` is set to `False` (default value), giving equal scale to all features.

To fix this bug, one should always compute `_y_train_std` from the training data and use the boolean attribute `normalize_y` to undo the normalization of `y_mean` if necessary.

### Steps/Code to Reproduce

```python
import pytest
from numpy import array
from numpy import hstack
from sklearn.gaussian_process import GaussianProcessRegressor

x = array([[0.], [0.5], [1.]])
y = hstack((x**2, 10*x**2))
# Note that the second output is equal to 10 times the first one.

# With output normalization
gpr = GaussianProcessRegressor(normalize_y=True)
gpr.fit(x, y)
std = gpr.predict(array([[0.25]]), return_std=True)[1][0]
assert std[0] != std[1]
assert std[0] == pytest.approx(std[1]/10, rel=1e-9)
# As expected, the variance of the second output is 10 times larger than the first output.

# Without output normalization
gpr = GaussianProcessRegressor(normalize_y=False)
gpr.fit(x, y)
std = gpr.predict(array([[0.25]]), return_std=True)[1][0]
assert std[0] == std[1]
# The variance of the second output is equal to the variance of the first output.
```

### Expected Results

Without output normalization, the variance of the second output should be 10 times larger than the first output.

### Actual Results

Without output normalization, the variance of the second output is equal to the variance of the first output.

### Versions

```shell
System:
    python: 3.9.12 (tags/v3.9.12:b28265d, Mar 23 2022, 23:52:46) [MSC v.1929 64 bit (AMD64)]
executable: C:\Users\matthias.delozzo\workspace\GEMSEO\gemseo\.tox\py39\.venv\Scripts\python.exe
   machine: Windows-10-10.0.19045-SP0

Python dependencies:
      sklearn: 1.5.1
          pip: None
   setuptools: None
        numpy: 1.26.4
        scipy: 1.13.1
       Cython: None
       pandas: 2.2.2
   matplotlib: 3.9.2
       joblib: 1.4.2
threadpoolctl: 3.5.0

Built with OpenMP: True

threadpoolctl info:
       user_api: blas
   internal_api: openblas
    num_threads: 8
         prefix: libopenblas
       filepath: C:\Users\matthias.delozzo\workspace\GEMSEO\gemseo\.tox\py39\.venv\Lib\site-packages\numpy.libs\libopenblas64__v0.3.23-293-gc2f4bdbb-gcc_10_3_0-2bde3a66a51006b2b53eb373ff767a3f.dll
        version: 0.3.23.dev
threading_layer: pthreads
   architecture: Haswell

       user_api: openmp
   internal_api: openmp
    num_threads: 8
         prefix: vcomp
       filepath: C:\Users\matthias.delozzo\workspace\GEMSEO\gemseo\.tox\py39\.venv\Lib\site-packages\sklearn\.libs\vcomp140.dll
        version: None

       user_api: blas
   internal_api: openblas
    num_threads: 8
         prefix: libopenblas
       filepath: C:\Users\matthias.delozzo\workspace\GEMSEO\gemseo\.tox\py39\.venv\Lib\site-packages\scipy.libs\libopenblas_v0.3.27--3aa239bc726cfb0bd8e5330d8d4c15c6.dll
        version: 0.3.27
threading_layer: pthreads
   architecture: Haswell
```

COMMENT:
Hi, this is my first time contributing. The bug was well described, and I reproduced the issue using the provided testing code.

As already suggested, the problem lies in how `_y_train_std` is handled when `normalize_y=False`. Here’s the code in `sklearn/gaussian_process/_gpr.py`:


```python
if self.normalize_y:
    self._y_train_mean = np.mean(y, axis=0)
    self._y_train_std = _handle_zeros_in_scale(np.std(y, axis=0), copy=False)
    y = (y - self._y_train_mean) / self._y_train_std
else:
    shape_y_stats = (y.shape[1],) if y.ndim == 2 else 1
    self._y_train_mean = np.zeros(shape=shape_y_stats)
    self._y_train_std = np.ones(shape=shape_y_stats)
```


The code should compute `_y_train_std` from the training data even when `normalize_y=False` (as suggested). I applied a quick fix with this change:


```python
# Normalize target value
        self._y_train_mean = np.mean(y, axis=0)
        self._y_train_std = _handle_zeros_in_scale(np.std(y, axis=0), copy=False)
        if self.normalize_y:
            # Remove mean and make unit variance
            y = (y - self._y_train_mean) / self._y_train_std

        else:
            self._y_train_mean = np.zeros_like(self._y_train_mean)
```


After this quick fix, now the test code gives the correct result. I hope this is helpful!

Labels: Bug