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:
Yes, I think you're fully right; we're getting a failure because now, with your fix, standard déviation and covariance _depend_ on output. I think you have the same failure with the covariance, haven't?

You could also robustify this test function  by testing not only the case where `normalize_y=False` but also `normalize_y=True`.

Labels: Bug