diff --git a/GPy/core/gp.py b/GPy/core/gp.py
index 0cc01026..977882ac 100644
--- a/GPy/core/gp.py
+++ b/GPy/core/gp.py
@@ -578,7 +578,7 @@ class GP(Model):
                 mag[n] = np.sqrt(np.linalg.det(G[n, :, :]))
         return mag
 
-    def posterior_samples_f(self,X, size=10, full_cov=True, **predict_kwargs):
+      def posterior_samples_f(self,X, size=10, **predict_kwargs):
         """
         Samples the posterior GP at the points X.
 
@@ -586,35 +586,28 @@ class GP(Model):
         :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)
+        :returns: set of simulations
+        :rtype: np.ndarray (Nnew x D x samples) 
         """
-        m, v = self._raw_predict(X,  full_cov=full_cov, **predict_kwargs)
+        m, v = self._raw_predict(X,  full_cov=True, **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):
-            if not full_cov:
-                return np.random.multivariate_normal(m.flatten(), np.diag(v.flatten()), size).T
-            else:
-                return np.random.multivariate_normal(m.flatten(), v, size).T
+            return np.random.multivariate_normal(m, v, size).T
 
         if self.output_dim == 1:
-            return sim_one_dim(m, v)
+            return sim_one_dim(m.flatten(), v)[:, np.newaxis, :]
         else:
-            fsim = np.empty((self.output_dim, X.shape[0], size))
+            fsim = np.empty((X.shape[0], self.output_dim, size))
             for d in range(self.output_dim):
-                if full_cov and v.ndim == 3:
-                    fsim[d] = sim_one_dim(m[:, d], v[:, :, d])
-                elif (not full_cov) and v.ndim == 2:
-                    fsim[d] = sim_one_dim(m[:, d], v[:, d])
+                if v.ndim == 3:
+                    fsim[:, d, :] = sim_one_dim(m[:, d], v[:, :, d])
                 else:
-                    fsim[d] = sim_one_dim(m[:, d], v)
+                    fsim[:, d, :] = sim_one_dim(m[:, d], v)
         return fsim
 
-    def posterior_samples(self, X, size=10, full_cov=False, Y_metadata=None, likelihood=None, **predict_kwargs):
+    def posterior_samples(self, X, size=10, Y_metadata=None, likelihood=None, **predict_kwargs):
         """
         Samples the posterior GP at the points X.
 
@@ -622,19 +615,17 @@ class GP(Model):
         :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.
         :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)
         """
-        fsim = self.posterior_samples_f(X, size, full_cov=full_cov, **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[0]):
-                fsim[d] = likelihood.samples(fsim[d], Y_metadata=Y_metadata)
+            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