'Draw samples from Gaussian process and evaluate at X.
Parameters
X : array-like, shape = (n_samples_X, n_features)
Query points where the GP samples are evaluated
n_samples : int, default: 1
The number of samples drawn from the Gaussian process
random_state : int, RandomState instance or None, optional (default=0)
If int, random_state is the seed used by the random number
generator; If RandomState instance, random_state is the
random number generator; If None, the random number
generator is the RandomState instance used by `np.random`.
Returns
y_samples : array, shape = (n_samples_X, [n_output_dims], n_samples)
Values of n_samples samples drawn from Gaussian process and
evaluated at query points.'

def sample_y(self, X, n_samples=1, random_state=0):

    rng = check_random_state(random_state)
    (y_mean, y_cov) = self.predict(X, return_cov=True)
    if (y_mean.ndim == 1):
        y_samples = rng.multivariate_normal(y_mean, y_cov, n_samples).T
    else:
        y_samples = [rng.multivariate_normal(y_mean[:, i], y_cov, n_samples).T[:, np.newaxis] for i in range(y_mean.shape[1])]
        y_samples = np.hstack(y_samples)
    return y_samples
