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:
I agree that the behavior described in the description looks like a bug.

I find it surprising because `normalize_y=False` is the default value of that parameter which means that most users of `GaussianProcessRegressor` would be impacted by this and get incorrect predictions as a result.

@Lollitor if you have time please open a PR to fix this and feel free to ping me on it for review.

TODO in the PR:

- make sure that existing tests pass with the fix, or find out if they need an update (and explain it in a comment in the PR);
- add a new non-regression test based on the code snippet provided in this issue;
- check the impact of this fix on the results (e.g. plots) when running all the examples that uses the `GaussianProcessRegressor`. Those can be found with:

```
$ git grep GaussianProcessRegressor examples/
```

Labels: Bug