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 all, sorry for the delay. In the past few days, I worked on this, and after making the modifications we previously agreed on, the original bug is now solved. All 445 tests in the test suite pass, except for one:

```
def test_y_multioutput():
    # Test that GPR can deal with multi-dimensional target values
    y_2d = np.vstack((y, y * 2)).T

    ...

    gpr = GaussianProcessRegressor(kernel=kernel, optimizer=None, normalize_y=False)
    gpr.fit(X, y)

    gpr_2d = GaussianProcessRegressor(kernel=kernel, optimizer=None, normalize_y=False)
    gpr_2d.fit(X, y_2d)

    y_pred_1d, y_std_1d = gpr.predict(X2, return_std=True)
    y_pred_2d, y_std_2d = gpr_2d.predict(X2, return_std=True)
    _, y_cov_1d = gpr.predict(X2, return_cov=True)
    _, y_cov_2d = gpr_2d.predict(X2, return_cov=True)

    assert_almost_equal(y_pred_1d, y_pred_2d[:, 0])
    assert_almost_equal(y_pred_1d, y_pred_2d[:, 1] / 2)

    # Standard deviation and covariance do not depend on output
    for target in range(y_2d.shape[1]):
        assert_almost_equal(y_std_1d, y_std_2d[..., target])
```

Correct me if I'm wrong, but this looks like the code initially posted to reproduce the bug, right? The standard deviation of the second output should scale according to the output when normalize_y=False. So I believe it's actually good that we’re getting a failure in this test.

Changing it to the following corrects the test:
```
for target in range(y_2d.shape[1]):
    if target == 0:
        assert_almost_equal(y_std_1d, y_std_2d[:, target])
    elif target == 1:
        assert_almost_equal(y_std_1d, y_std_2d[:, target] / 2)
```

Assuming I’ve got this right, could you please confirm? Meanwhile, I'll continue reviewing the test suite for other potential issues, even though the other tests pass, and I’ll also address points 2 and 3 from the TODO list by @ogrisel .

Sorry again for the delay!

Labels: Bug