Creating a partial dependant plot for a prediction function I am working on creating a partial dependant plot for one of my features (B). The problem is that I didn't use any model to predict my output (R). I've used a prediction function to do that:
def objective(B, D, k, fc, Le, N):
    return (k*fc**1.3*B**3.3*D**1.8)/(10**5*N**1.5*Le**0.9)

data=data_set[["B","D","k","fc","Le","N"]]
#data = {'B':[1,2,3],'D':[4,5,6]}
df = pd.DataFrame(data)

df['R'] =  df.apply(lambda x: objective(x.B, x.D, x.k,x.fc,x.Le,x.N), axis=1)

print(df)

What did I try?
I've tried using the shap function to create the plot, yet I didn't know what to replace question_mark.predict with:
import shap
import sklearn
X100 = shap.utils.sample(X, 100) # 100 instances for use as the background distribution

shap.partial_dependence_plot(
    "B", question_mark.predict, X100, ice=False,
    model_expected_value=True, feature_expected_value=True)

How can I do this?
 A: TLDR: A .predict() method and a dataset X are sufficient to make partial dependence plots (PDP).
A PDP shows the marginal effect of feature $x_i$ by plotting $x_i$ on the x-axis and the partial dependence function $\widehat{f_i}$ on the y-axis:
$$
\begin{aligned}
\widehat{f_i}(x_i) = \frac{1}{n}\sum_{j=1}^n\hat{f}(x_i,\mathbf{x}_{-i}^{(j)})
\end{aligned}
$$
where $\hat{f}$ generates predictions according to the fitted model, $\mathbf{x}_{-i}$ are the features other than $x_i$ and $\mathbf{x}_{-i}^{(j)}$ are their values for the $j$th observation in the dataset $\mathbf{X}$. So we vary $x_i$ over a grid of values (usually determined by its range in the dataset) but we keep the rest of the features $\mathbf{x}_{-i}$ fixed at their observed values.
To learn more about partial dependence plots, see for example Interpretation of y-axis in partial dependence plot, interpreting y axis of a partial dependence plots, Meaning of y axis in Random Forest partial dependence plot and the Interpretable Machine Learning book.

The question asks for a python implementation. I show how to do it with scikit-learn's PartialDependenceDisplay by creating a dummy estimator that applies the predict function to the dataset X. The example model is a linear regression with two predictors.
Note: According to the documentation, the PartialDependenceDisplayis.from_estimator method is new in version 1.0. So make sure you have scikit-learn ≥ 1.0 to run the example below.
from sklearn.base import BaseEstimator, RegressorMixin
from sklearn.utils import check_X_y, check_array
from sklearn.linear_model import LinearRegression
from sklearn.datasets import make_regression
from sklearn.inspection import PartialDependenceDisplay
import matplotlib.pyplot as plt

# Generate sample
X, y = make_regression(
    n_samples=100, n_features=2, n_informative=1,
    noise=1, random_state=1234
)
# Fit regressor
reg = LinearRegression(fit_intercept=False)
reg.fit(X, y)
b0 = reg.coef_[0]
b1 = reg.coef_[1]

Make partial dependence plots (PDP) from the regressor object.
display = PartialDependenceDisplay.from_estimator(reg, X, [0, 1])
display.plot()


This is for illustration only; we will make exactly the same plots from the regression function $\hat{y} = \mathbf{X}\hat{\boldsymbol{\beta}} = \hat{\beta}_0x_0 + \hat{\beta}_1x_1$.
def objective(x0, x1):
    return b0 * x0 + b1 * x1

The idea is to create a DummyRegressor class. Its .fit(X, y) function doesn't do anything exciting; just accepts the X and y arguments. Its .predict(X) function calls the objective function.
class DummyRegressor(BaseEstimator, RegressorMixin):
    def fit(self, X, y):
        # Check that X and y have correct shape
        X, y = check_X_y(X, y, y_numeric=True)
        self.X_ = X
        self.y_ = y
        return self

    def predict(self, X):
        X = check_array(X)

        # X is an numpy array, not a pandas DataFrame.
        # Make sure the columns of X are in the order required by `objective`
        return objective(X[:, 0], X[:, 1])

And now we make the partial dependence plots. Since PDPs are the same, I've labeled the x-axis to distinguish them from the previous display.
dummy_reg = DummyRegressor()
dummy_reg.fit(X, y)

display = PartialDependenceDisplay.from_estimator(
    dummy_reg, X, [0, 1],
    feature_names=["x0", "x1"]
)
display.plot()


