[plotting] added samples plot

GPy/core/gp.py

@@ -437,7 +437,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):
+    def posterior_samples_f(self,X, size=10, full_cov=True, **predict_kwargs):
         """
         Samples the posterior GP at the points X.
 
@@ -448,20 +448,32 @@ class GP(Model):
         :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 (N x samples)
+        :rtype: np.ndarray (D x N x samples) (if D==1 we flatten out the first dimension)
         """
-        m, v = self._raw_predict(X,  full_cov=full_cov)
+        m, v = self._raw_predict(X,  full_cov=full_cov, **predict_kwargs)
         if self.normalizer is not None:
             m, v = self.normalizer.inverse_mean(m), self.normalizer.inverse_variance(v)
-        v = v.reshape(m.size,-1) if len(v.shape)==3 else v
-        if not full_cov:
-            fsim = np.random.multivariate_normal(m.flatten(), np.diag(v.flatten()), size).T
-        else:
-            fsim = np.random.multivariate_normal(m.flatten(), v, size).T
+        
+        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
 
+        if self.output_dim == 1:
+            return sim_one_dim(m, v)
+        else:
+            fsim = np.empty((self.output_dim, self.num_data, 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])
+                else:
+                    fsim[d] = sim_one_dim(m[:, d], v)
         return fsim
 
-    def posterior_samples(self, X, size=10, full_cov=False, Y_metadata=None):
+    def posterior_samples(self, X, size=10, full_cov=False, Y_metadata=None, likelihood=None, **predict_kwargs):
         """
         Samples the posterior GP at the points X.
 
@@ -473,11 +485,18 @@ class GP(Model):
         :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, a Numpy array (N x samples).
+        :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)
-        Ysim = self.likelihood.samples(fsim, Y_metadata=Y_metadata)
-        return Ysim
+        fsim = self.posterior_samples_f(X, size, full_cov=full_cov, **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)
+        else:
+            fsim = likelihood.samples(fsim, Y_metadata=Y_metadata)
+        return fsim
 
     def input_sensitivity(self, summarize=True):
         """