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

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

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

Start at GaussianProcessRegressor.predict and inspect how _y_train_std is set and used for multi-output standard deviation and covariance. Reproduce the reported normalize_y=False case and add regression coverage showing outputs preserve their training-data scale, including the expected 10:1 relationship.