code: def posterior_samples_f(self,X, size=10, **predict_kwargs):
          
          predict_kwargs["full_cov"] = True  # Always use the full covariance for posterior samples. 
          m, v = self._raw_predict(X,  **predict_kwargs)
          if self.normalizer is not None:
              m, v = self.normalizer.inverse_mean(m), self.normalizer.inverse_variance(v)
  
          def sim_one_dim(m, v):
              return np.random.multivariate_normal(m, v, size).T
  
          if self.output_dim == 1:
              return sim_one_dim(m.flatten(), v)[:, np.newaxis, :]
          else:
              fsim = np.empty((X.shape[0], self.output_dim, size))
              for d in range(self.output_dim):
                  if v.ndim == 3:
                      fsim[:, d, :] = sim_one_dim(m[:, d], v[:, :, d])
                  else:
                      fsim[:, d, :] = sim_one_dim(m[:, d], v)
          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.
          :returns: set of simulations
          :rtype: np.ndarray (Nnew x D x samples)