def posterior_samples_f(self,X, size=10, **predict_kwargs):
        """
        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) 
        """
        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