Correct scale back for PLS regression coefficients

### Describe the bug

In `cross_decomposition/_pls.py`, PLS regression coefficients are calculated in class `_PLS` (starts at line 165). In this class, when `scale=True`, data are scaled (on line 265). In that case, the resulting regression coefficients need to be scaled back to the original scale, such that they represent the relationship between the original X and y. The way this scale back is done on line 360, is wrong: https://github.com/scikit-learn/scikit-learn/blob/3f89022fa04d293152f1d32fbc2a5bdaaf2df364/sklearn/cross_decomposition/_pls.py#L360
This is wrong because rescaling is done by only adjusting for the `y_std`. 

The correct formula is to scale back as `self.coef_ = (self.coef_ * self._y_std/self.x_std).T`. It is easy to verify the latter for ordinary least squares regression as they will match exactly. As PLS is only rotationally invariant, the rescaled coefficients from this proposal will not match exactly, but they will be much closer to coefficients from unscaled data than the present version. Example code given below.

### Steps/Code to Reproduce

```python
import numpy as np 
from sklearn.linear_model import LinearRegression
from direpack import VersatileScaler
from sklearn.cross_decomposition import PLSRegression
```

Simulate Data

```python
X = np.random.multivariate_normal(np.zeros(3),np.diag(np.ones(3)),500)
YL = np.dot(X,np.array([3,-3,5]).reshape((-1,1))) + np.random.multivariate_normal(np.zeros(3),np.diag(np.ones(3)/100),500)
```

Test both options for OLS regression 
```python
LRns = LinearRegression()
LRns.fit(X,YL)
LRns.coef_
LRs = LinearRegression()
LRs.fit(Xs,YLs)
LRs.coef_
```

`LRs.coef_` and `LRns.coef_` are identical 

Now test both options for PLS using internal scaling

```python
PLSns = PLSRegression(n_components=2, scale=False)
PLSns.fit(X,YL)
PLSns.coef_
PLSms = PLSRegression(n_components=2, scale=False)
PLSms.fit(Xs,YLs)
PLSms.coef_
```
Off! 

Now the proposed solution 
```python
PLSms = PLSRegression(n_components=2, scale=False)
PLSms.fit(Xs,YLs)
PLSms.coef_ 
PLSms.coef_ *sY/sX
```
Much closer !!

### Expected Results

Rescaled regression coefficients much closer to coefficients calculated from unscaled data 

### Actual Results

Actual coefficients are off. One example running the code above:

- Coefficients from unscaled data : 
```
array([[ 2.74182398, -3.06332548,  5.08010417],
       [ 2.73919939, -3.0741624 ,  5.09501747],
       [ 2.7424613 , -3.06724517,  5.08589201]])
```
- Coefficients using the `scale=True` rescaling option: 
```
array([[ 2.74182398, -3.06332548,  5.08010417],
       [ 2.73919939, -3.0741624 ,  5.09501747],
       [ 2.7424613 , -3.06724517,  5.08589201]])
--> OFF 
```

- Coefficients that would be obtained from proposed solution: 
```
array([[ 3.0012696 , -3.00657879,  4.9926381 ],
       [ 2.99094369, -3.00971523,  4.99484826],
       [ 2.99880995, -3.00725971,  4.99306938]])
--> much closer! 
```

### Versions

```shell
System:
    python: 3.11.4 | packaged by Anaconda, Inc. | (main, Jul  5 2023, 13:38:37) [MSC v.1916 64 bit (AMD64)]
executable: C:\Users\serneels\AppData\Local\anaconda3\python.exe
   machine: Windows-10-10.0.19044-SP0

Python dependencies:
      sklearn: 1.3.0
          pip: 23.2.1
   setuptools: 68.0.0
        numpy: 1.24.3
        scipy: 1.10.1
       Cython: 0.29.36
       pandas: 1.5.3
   matplotlib: 3.7.1
       joblib: 1.2.0
threadpoolctl: 2.2.0

Built with OpenMP: True
```

COMMENT:
Could you please provide a reference or source that explains this concept? I want to better understand the reasoning behind this modification. The original implementation, `y_pred = X_scaled @ (np.dot(self.x_rotations_, self.y_loadings_.T)) *self.y_std + self.y_mean` felt intuitive to me, but I'm not an expert in this field.

Labels: Bug