I made the suggested changes, but I noticed an additional problem. The function sample_y() itself has a n_samples parameter, and the argument X has shape (n_samples, n_features), but the two n_samples are not the same. This can lead to a lot of confusion.

def sample_y(self, X, n_samples=1, random_state=0):
    """Draw samples from Gaussian process and evaluate at X.
    Parameters
    ----------
    X : array-like of shape (n_samples, n_features) or list of object
        Query points where the GP is evaluated.
    n_samples : int, default=1
        The number of samples drawn from the Gaussian process
     ...

The existing documentation uses n_samples_X in the return value to distinguish between the two n_samples, but that still leaves the two n_samples in the parameters descriptions.

Thoughts on how to best resolve this duplicity? @thomasjpfan

What do you think of updating the docstring of X to the following?

         X : array-like of shape (n_samples_X, n_features) or list of object

@thomasjpfan: agreed

I also made some other small fixes / clarifications as well.

@thomasjpfan I undid the np.dot -> @ code changes

Bump. Will these changes be merged soon? @thomasjpfan

I'd rather not have to waste time creating a 2nd pull request for the code changes. Is there any way to get a 2nd reviewer on this? Thanks!