This function generates conditional predictive surfaces by evaluating the fitted model across a grid of values for one or two selected covariates.
For each of the nsamp draws, the remaining covariates are either held fixed (if fixed_x is provided) or filled in by sampling a single row from the training data.
The selected covariates in to_eval are then varied across a regular grid defined by n_eval_points and eval_range, and model predictions are computed using calc_pred_moments.
The resulting samples represent conditional predictions across different covariate contexts, and can be used to visualize marginal effects, interaction surfaces, or predictive uncertainty.
Note that computational and memory requirements increase rapidly with grid size. In particular, for two-dimensional evaluations, the
kernel matrix scales quadratically with the number of evaluation points per axis. Large values of n_eval_points may lead to high
memory usage during prediction, especially when using a GPU. If memory constraints arise, consider reducing n_eval_points.