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:
Thanks for the super quick reply! Yes, indeed these issues are quite similar to the "interventional" or "tree_path_dependent" SHAP discussions!
Maybe one difference here is that the "recursive" way is MUCH faster than the "brute" method, hence its appeal.

While I did sketch out a conceptual fix, I am unfortunately nowhere close to implementing this in code. Especially not at the rigorous sklearn quality level.
So I would probably pass on that nice opportunity for now.

Thanks,
ML

Labels: Documentation, help wanted