code: def posterior_samples(self, X, size=10, Y_metadata=None, likelihood=None, **predict_kwargs):
          
          fsim = self.posterior_samples_f(X, size, **predict_kwargs)
          if likelihood is None:
              likelihood = self.likelihood
          if fsim.ndim == 3:
              for d in range(fsim.shape[1]):
                  fsim[:, d] = likelihood.samples(fsim[:, d], Y_metadata=Y_metadata)
          else:
              fsim = likelihood.samples(fsim, Y_metadata=Y_metadata)
          return fsim
docs: Samples the posterior GP at the points X.
  
          :param X: the points at which to take the samples.
          :type X: np.ndarray (Nnew x self.input_dim.)
          :param size: the number of a posteriori samples.
          :type size: int.
          :param noise_model: for mixed noise likelihood, the noise model to use in the samples.
          :type noise_model: integer.
          :returns: Ysim: set of simulations,
          :rtype: np.ndarray (D x N x samples) (if D==1 we flatten out the first dimension)