partial_dependence() with method recursion computes conditional partial dependence for trees

### Describe the bug

For the case of correlated predictors (clearly highly common) the `sklearn.inspection.partial_dependence()` function gives different answers for `method` = "recursion" and `method` = "brute", see my [post](https://markusloecher.github.io/Partial-Dependence-Trees/) for elaborate examples.

I do not believe that this is intentional and should be fixed. Alternatively, it should be communicated clearly in the documentation that (i) the two methods are not equivalent for tree based algorithms, and (ii) that `method` = "recursion" actually computes the **conditional** $E[f(x_S,X_C)|X_S=x_s]$ instead of the (desired) **interventional** $E[f(x_S,X_C)| \mathbf{do}(X_S=x_s)]$ 

### Steps/Code to Reproduce

```python
from sklearn.tree import DecisionTreeRegressor
import numpy as np
import pandas as pd
from sklearn.inspection import PartialDependenceDisplay, partial_dependence

#Generate the data
#first X
N = 400; p0 = 0.5; p11 = 0.8; M=2
X = np.zeros((N,M)) # a matrix (N * M)
N1 = int(p0*N)
X[0:N1,0] = 1
X[0:int(p11*N1),1] = 1
X[N1:(N1+int((1-p11)*(N-N1))+1),1] = 1
df = pd.DataFrame(X, columns=["X0", "X1"])
#then Y
y = np.zeros(N)
y[(X[:,0] == 0) & (X[:,1] == 0)] = 0.3
y[(X[:,0] == 0) & (X[:,1] == 1)] = 0.7
y[(X[:,0] == 1) & (X[:,1] == 0)] = 0.9
y[(X[:,0] == 1) & (X[:,1] == 1)] = 0.1

df = pd.DataFrame(X, columns=["X0", "X1"])
df["y"] = y

model = DecisionTreeRegressor(max_depth=2, max_features = 2)
model.fit(df[["X0","X1"]], y);

features = ["X1"]
for f in features:
    pdp_interventional = partial_dependence(model, df[["X0","X1"]], f,method = "brute")
    pdp_conditional = partial_dependence(model, df[["X0","X1"]], f,method = "recursion")
    
    print(f, "brute (interventional):", pdp_interventional['average'])
    print(f, "recursion (conditional):", pdp_conditional['average'])

pdp_interventional['average'] == pdp_conditional['average']
```

### Expected Results

We would like the two methods to yield the same pdp values, so the last line should yield

`array([[True, True]])`

### Actual Results

Instead they are different. The print statements yield
_X1 brute (interventional): [[0.6 0.4]]
X1 recursion (conditional): [[0.42 0.22]]_
and the last line yields
`array([[False, False]])`

### Versions

```shell
System:
    python: 3.9.0 (v3.9.0:9cf6752276, Oct  5 2020, 11:29:23)  [Clang 6.0 (clang-600.0.57)]
executable: /Users/loecherm/.pyenv/versions/aug_hs_env/bin/python
   machine: macOS-10.16-x86_64-i386-64bit

Python dependencies:
      sklearn: 1.2.0
          pip: 23.0.1
   setuptools: 49.2.1
        numpy: 1.23.5
        scipy: 1.10.0
       Cython: None
       pandas: 1.5.3
   matplotlib: 3.5.3
       joblib: 1.2.0
threadpoolctl: 3.1.0

Built with OpenMP: True
```

COMMENT:
Well, of course both methods have their place. I would still consider the current setup a bug since the assumption is that for trees the recursive way yields the same results as the brute force; exemplified e.g. by [this internal sklearn code]( https://github.com/scikit-learn/scikit-learn/blob/38a06e4be504f3971d109d6741b8b4c7192d7323/sklearn/inspection/tests/test_partial_dependence.py#L377):

Labels: Documentation, help wanted