Cloned estimators have identical randomness but different RNG instances

### Describe the bug

Cloned estimators have identical randomness but different RNG instances. According to documentation, it should be the other way around: different randomness but identical RNG instances.

Related #25395 

The User Guide [says](https://scikit-learn.org/stable/common_pitfalls.html#controlling-randomness):
> For an optimal robustness of cross-validation (CV) results, pass RandomState instances when creating estimators

> ```python
> rf_inst = RandomForestClassifier(random_state=np.random.RandomState(0))
> cross_val_score(rf_inst, X, y)
> ```
> ...
> Since `rf_inst` was passed a `RandomState` instance, each call to `fit` starts from a different RNG. As a result, the random subset of features will be different for each folds

In regards to cloning, the same reference says:
> ```python
> rng = np.random.RandomState(0)
> a = RandomForestClassifier(random_state=rng)
> b = clone(a)
> ```
> Moreover, `a` and `b` will influence each-other since they share the same internal RNG: calling `a.fit` will consume `b`’s RNG, and calling `b.fit` will consume `a`’s RNG, since they are the same. 

The actual behaviour does not follow this description. 

### Steps/Code to Reproduce

```python
import numpy as np
from sklearn import clone
from sklearn.datasets import make_classification
from sklearn.model_selection import cross_validate
from sklearn.ensemble import RandomForestClassifier

rng = np.random.RandomState(0)
X, y = make_classification(random_state=rng)
rf = RandomForestClassifier(random_state=rng)

d = cross_validate(rf, X, y, return_estimator=True, cv=2)
rngs = [e.random_state for e in d['estimator']]
# estimators corresponding to different CV runs have different but identical RNGs:
print(rngs[0] is rngs[1]) # False
print(all(rngs[0].randint(10, size=10) == rngs[1].randint(10, size=10))) # True

rf_clone = clone(rf)
rngs = [rf.random_state, rf_clone.random_state]
print(rngs[0] is rngs[1]) # False
print(all(rngs[0].randint(10, size=10) == rngs[1].randint(10, size=10))) # True
```

### Expected Results

True False True False

### Actual Results

False True False True

### Versions

```shell
Tested on a two-week-old dev build and also on the following version (Kaggle)

System:
    python: 3.7.12 | packaged by conda-forge | (default, Oct 26 2021, 06:08:53)  [GCC 9.4.0]
executable: /opt/conda/bin/python3.7
   machine: Linux-5.15.90+-x86_64-with-debian-bullseye-sid

Python dependencies:
          pip: 22.3.1
   setuptools: 59.8.0
      sklearn: 1.0.2
        numpy: 1.21.6
        scipy: 1.7.3
       Cython: 0.29.34
       pandas: 1.3.5
   matplotlib: 3.5.3
       joblib: 1.2.0
threadpoolctl: 3.1.0

Built with OpenMP: True
```

COMMENT:
OK.

I've been reluctant to publish https://github.com/scikit-learn/enhancement_proposals/pull/88 ([rendered docs are here](https://scikit-learn-enhancement-proposals--88.org.readthedocs.build/en/88/slep022/proposal.html)) so far because I won't be able to champion it has become slightly out of date. But I think the SLEP does a decent job at describing the [key requirements](https://scikit-learn-enhancement-proposals--88.org.readthedocs.build/en/88/slep022/proposal.html#key-use-cases-and-requirements) we need to support for a solution, so I hope it can help better framing the discussion.

FWIW I think the solution advocated in the SLEP is quite similar to what you suggested above @betatim. The main idea is to 1) make estimators (and splitters) fully stateless across calls to `fit()`, and 2) expose the randomness complexity to users by letting them choose the kind of CV procedure they want to run (with good defaults). Right now this complexity is both hidden *and* largely uncontrollable.

For completeness, a toy implementation of the solution is implemented [in this notebook](https://nbviewer.org/gist/NicolasHug/2db607b01482988fa549eb2c8770f79f).

Labels: Bug, Documentation, Needs Decision