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:
Finally getting to this, here is the promised fix, at least conceptually.
The following is a sketch of a recursive algorithm which I believe would correctly compute the interventional PDP in **one** pass over the training data (for all $X_S$):

At the root of the tree, a weight value of $w=0$ is assigned. For each nonterminal node visited, if its split variable is in the complement subset $X_C$, the appropriate left or right daughter node is visited and the weight is not modified. If the node’s split variable is a member of the target subset $X_S$, then **both daughters** are visited and the corresponding number of (training) samples are added to the respective weight $w$.

Each terminal node visited during the traversal is assigned the
current value of the weight divived by the total number of observations. When the tree traversal is complete, the value of $f_S(X_S)$ is the corresponding weighted average of the $f(X)$ values over those terminal nodes visited during the tree
traversal.

While this scheme does need to pass the entire training data, that cost is mitigated by two factors:

- Constant cost for all X_S
- One could implement it into the fit function where the training data are passed through the tree while building it, instead of a post-hoc approach.

Also, the original idea
> At each non-terminal node that splits on a complementary feature both children are visited

can still be put to good use to determine the set of terminal regions $T(X_S)$ that can be reached by data points of the form $(X_S,x_C)$ for some $x_C$.
$T(X_S)$ can be computed by following both children of nodes that split on variables in C.

Labels: Documentation, help wanted