code: def posterior_samples_f(self, X, size=10, full_cov=True, **predict_kwargs):
          
          mu, var = self._raw_predict(X, full_cov=full_cov, **predict_kwargs)
          if self.normalizer is not None:
              mu, var = self.normalizer.inverse_mean(mu), self.normalizer.inverse_variance(var)
  
          def sim_one_dim(m, v):
              nu = self.nu + 2 + self.num_data
              v = np.diag(v.flatten()) if not full_cov else v
              Z = np.random.multivariate_normal(np.zeros(X.shape[0]), v, size).T
              g = np.tile(np.random.gamma(nu / 2., 2. / nu, size), (X.shape[0], 1))
              return m + Z / np.sqrt(g)
  
          if self.output_dim == 1:
              return sim_one_dim(mu, var)
          else:
              fsim = np.empty((self.output_dim, self.num_data, size))
              for d in range(self.output_dim):
                  if full_cov and var.ndim == 3:
                      fsim[d] = sim_one_dim(mu[:, d], var[:, :, d])
                  elif (not full_cov) and var.ndim == 2:
                      fsim[d] = sim_one_dim(mu[:, d], var[:, d])
                  else:
                      fsim[d] = sim_one_dim(mu[:, d], var)
          return fsim
docs: Samples the posterior TP 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 full_cov: whether to return the full covariance matrix, or just the diagonal.
          :type full_cov: bool.
          :returns: fsim: set of simulations
          :rtype: np.ndarray (D x N x samples) (if D==1 we flatten out the first dimension)