RandomForestRegressor results are not reproducible

### Describe the bug

We use the ''RandomForestRegressor'' from ''sklearn.ensemble'' to produce estimates of conditional means. Despite setting a random seed, the estimated conditional means are not reproducible when we enable multiprocessing. For some observations, conditional mean estimates change when rerunning the same code. We use the ''RandomForestRegressor'' in our statistical analysis. The non-reproducibility in the predictions of the ''RandomForestRegressor'' carries over to our estimates. Below, we provide an example to illustrate our issue. 

### Steps/Code to Reproduce

```python
import pandas as pd
import numpy as np
from sklearn.ensemble import RandomForestRegressor

# %% Data

df = pd.DataFrame(
    [[4.0, 47.0, 0.0],
     [-3640.913, 35.0, 1.0],
     [-32622.293, 25.0, 0.0],
     [-92.50659, 22.0, 0.0],
     [1179.845, 20.0, 0.0],
     [8750.0, 37.0, 0.0],
     [3559.3333, 24.0, 0.0],
     [-4778.8965, 30.0, 0.0],
     [-16032.402, 61.0, 1.0],
     [520.00903, 36.0, 0.0],
     [4428.75, 48.0, 0.0],
     [9622.75, 53.0, 1.0],
     [5254.0076, 30.0, 0.0],
     [1688.47, 28.0, 0.0],
     [9.0, 60.0, 0.0],
     [22777.727, 51.0, 0.0],
     [200.0, 58.0, 0.0],
     [-4501.4166, 50.0, 0.0],
     [368.795, 22.0, 0.0],
     [20.0, 69.0, 1.0],
     [8858.34, 18.0, 0.0],
     [43.0, 22.0, 1.0],
     [93.0332, 29.0, 0.0],
     [5221.1133, 23.0, 0.0],
     [7197.25, 29.0, 1.0],
     [3250.36, 48.0, 1.0],
     [547.345, 47.0, 0.0],
     [27.0, 45.0, 0.0],
     [2306.2524, 18.0, 0.0],
     [8371.8335, 51.0, 0.0],
     [-7143.3364, 41.0, 1.0],
     [1435.5533, 37.0, 0.0],
     [-2594.5981, 25.0, 0.0],
     [43.0, 30.0, 0.0],
     [1260.3575, 48.0, 0.0],
     [-88.98999, 20.0, 0.0],
     [-22231.859, 50.0, 0.0],
     [-3400.951, 39.0, 0.0],
     [4192.92, 36.0, 0.0],
     [769.2666, 21.0, 0.0]])

# %% Parameters

NUMBER_OF_WORKERS = [1, 20]
NUMBER_OF_ROWS_TRAINING = 20
NUMBER_OF_ROWS_PREDICTION = 20
N_ESTIMATORS = 1000
MIN_SAMPLES_LEAF = 5
MAX_FEATURES = 1
BOOTSTRAP = True
OOB_SCORE = True
RANDOM_STATE = 42
rng = np.random.default_rng(9324561)

# %% Estimation

for workers in NUMBER_OF_WORKERS:
    for i in range(2):
        x_np = df.iloc[:, 1:].to_numpy()  # features
        y_np = df.iloc[:, :1].to_numpy()  # outcome
        regr = RandomForestRegressor(
            n_estimators=N_ESTIMATORS, min_samples_leaf=MIN_SAMPLES_LEAF,
            max_features=MAX_FEATURES, bootstrap=BOOTSTRAP,
            oob_score=OOB_SCORE, random_state=RANDOM_STATE,
            n_jobs=workers)
        regr.fit(x_np[:NUMBER_OF_ROWS_TRAINING, :],
                 y_np[:NUMBER_OF_ROWS_TRAINING, 0].ravel())
        y_pred_rf = regr.predict(x_np[NUMBER_OF_ROWS_PREDICTION:, :])
        if i == 0:
            y0 = y_np[NUMBER_OF_ROWS_PREDICTION:, 0] - y_pred_rf
        if i == 1:
            y1 = y_np[NUMBER_OF_ROWS_PREDICTION:, 0] - y_pred_rf
    print('Share of non-equal values using', workers, 'worker(s): ',
          np.mean(np.equal(y0, y1)))
```

### Expected Results

We would expect that - fixing a random seed - results should be reproducible. That is not the case. 

### Actual Results

Results are not reproducible. 

### Versions

```shell
0.24.2
```

COMMENT:
I think that we already observed it. It could be linked to the accumulation for the prediction that will not be deterministic in parallel processing as stated here: https://github.com/scikit-learn/scikit-learn/issues/22303#issuecomment-1023208722

It is the reason why it should be reproducible with `n_jobs=1` and not with `n_jobs=-1`. I don't think that we can do anything indeed.

Labels: Bug, Needs Triage