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:
Just following up one more time on this. 
It seems that the default is still "recursive" for tree based methods ?
Might the _Note_ 
> The ‘brute’ and ‘recursion’ methods will likely disagree regarding the value of the partial dependence, because they will treat these unlikely samples differently.

be new ? I am not certain (yet) that the differences are caused by highly unlikely samples.

Labels: Documentation, help wanted