Unexpected behavior of sklearn.feature_selection.mutual_info_regression if copy=False

### Describe the bug

The parameter `copy` of the function `mutual_info_regression` is described as follows https://github.com/scikit-learn/scikit-learn/blob/d1d1596fac19d688a637690134d71fc460f5f0dd/sklearn/feature_selection/_mutual_info.py#L381-L383

I read it as both `X` and `y` should be modified if `copy=False` and `X` has continuous features. However, `y`  is always copied. I think the lines 
https://github.com/scikit-learn/scikit-learn/blob/d1d1596fac19d688a637690134d71fc460f5f0dd/sklearn/feature_selection/_mutual_info.py#L309-L310 should be 
```python
    if not discrete_target:
        y = y.astype(np.float64, copy=copy)
        y = scale(y, with_mean=False, copy=False)
```
Similarly to the the treatment of `X` https://github.com/scikit-learn/scikit-learn/blob/d1d1596fac19d688a637690134d71fc460f5f0dd/sklearn/feature_selection/_mutual_info.py#L295-L299

### Steps/Code to Reproduce

```python
import numpy as np
from sklearn.feature_selection import mutual_info_regression
n_samples_, n_feats = 30, 2
X = np.random.randn(n_samples_, n_feats)
y = np.random.randn(n_samples_, )
y_copy = y.copy()
X_copy = X.copy()
mutual_info_regression(X, y, copy=False)
print(np.allclose(y, y_copy), np.allclose(X, X_copy))
```

### Expected Results

The result should be
```
False, False
```
since both `X` and `y` should be modified in place by the function `mutual_info_regression`. 



### Actual Results

```
True, False
```

### Versions

```shell
System:
    python: 3.11.7 (main, Dec  5 2023, 19:13:35) [GCC 10.2.1 20210110]
executable: /usr/local/bin/python
   machine: Linux-5.15.146.1-microsoft-standard-WSL2-x86_64-with-glibc2.31

Python dependencies:
      sklearn: 1.3.0
          pip: 23.3.1
   setuptools: 69.0.2
        numpy: 1.23.2
        scipy: 1.12.0
       Cython: 3.0.9
       pandas: 2.1.4
   matplotlib: 3.8.2
       joblib: 1.3.2
threadpoolctl: 3.2.0

Built with OpenMP: True

threadpoolctl info:
       user_api: blas
   internal_api: openblas
    num_threads: 12
         prefix: libopenblas
       filepath: /home/user/.local/lib/python3.11/site-packages/numpy.libs/libopenblas64_p-r0-742d56dc.3.20.so
        version: 0.3.20
threading_layer: pthreads
   architecture: Prescott

       user_api: blas
   internal_api: openblas
    num_threads: 12
         prefix: libopenblas
       filepath: /home/user/.local/lib/python3.11/site-packages/scipy.libs/libopenblasp-r0-23e5df77.3.21.dev.so
        version: 0.3.21.dev
threading_layer: pthreads
   architecture: Prescott

       user_api: openmp
   internal_api: openmp
    num_threads: 12
         prefix: libgomp
       filepath: /home/user/.local/lib/python3.11/site-packages/scikit_learn.libs/libgomp-a34b3233.so.1.0.0
        version: None
```

COMMENT:
This issue may go in a bit different direction now :) Let me know if i should separate it into a new one. 

I am updating a third party sklearn plug-in [for feature selection](https://github.com/dwavesystems/dwave-scikit-learn-plugin), where we want to add a possibility for feature selection using conditional mutual information based on [the paper](https://citeseerx.ist.psu.edu/document?repid=rep1&type=pdf&doi=d121fb9a2e2131e5e9481400ad9f8bc25cbdca59). Note that unlike the paper itself, our solution will use a proprietary solver, but the plug-in itself is open-source. 

Mathematically we will aim to solve

$$
\max\limits_{\delta_i} \sum\limits_{i=1}^n I(x_i; y) \delta_i + \sum\limits_{i, j=1}^n I(x_i; y | x_j) \delta_i \delta_j
$$

where $x_i$ are features, $y$ is the target variable, and the variables $\delta_i\in\{0,1\}$ are binary with an associated constraint $\sum \delta_i = K$ for a pre-defined value $K$. If the variable $\delta_i =1$ then we will pick the feature $i$.

The conditional mutual information algorithm will follow mostly [this implementation](https://github.com/jannisteunissen/mutual_information), which is in turn based on the `sklearn` implementation for mutual information computation. 

Essentially, I want to compute the values $I(x_i; y)$, which I can do using the public methods `mutual_info_classif`, `mutual_info_regression` from `sklearn.feature_selection` or a private method `_compute_mi` from `sklearn.feature_selection._mutual_info`. 

1. If I use the public method `mutual_info_regression`,  `X` and `y` will be rescaled if they correspond to a data of a non-discrete distribution. I will then need to rescale `y` again for computing $I(x_i; y | x_j)$. If `X` and `y` would have the same behavior, then I can avoid rescaling `y` twice and essentially using different data `y` for computing $I(x_i; y)$ and $I(x_i; y | x_j)$. 
2. If we can make the method  `_compute_mi` public, then I can use my custom function for computing the values of $I(x_i; y | x_j)$ and $I(x_i; y)$ simultaneously. 

The second option would also allow for using the problem formulation described in [the mRMR paper](https://ieeexplore.ieee.org/document/1453511). Similar to an existing issue https://github.com/scikit-learn/scikit-learn/issues/8889 but using our solver. I was actually considering making the second option (making `_compute_mi` public) an issue too, but decided to start with a smaller issue. 

I am happy to submit a PR on some of these issues.

Labels: Documentation, help wanted