code: def log_predictive_density_sampling(self, y_test, mu_star, var_star, Y_metadata=None, num_samples=1000):
          
          assert y_test.shape==mu_star.shape
          assert y_test.shape==var_star.shape
          assert y_test.shape[1] == 1
  
          #Take samples of p(f*|y)
          #fi_samples = np.random.randn(num_samples)*np.sqrt(var_star) + mu_star
          fi_samples = np.random.normal(mu_star, np.sqrt(var_star), size=(mu_star.shape[0], num_samples))
  
          from scipy.misc import logsumexp
          log_p_ystar = -np.log(num_samples) + logsumexp(self.logpdf(fi_samples, y_test, Y_metadata=Y_metadata), axis=1)
          log_p_ystar = np.array(log_p_ystar).reshape(*y_test.shape)
          return log_p_ystar
docs: Calculation of the log predictive density via sampling
  
          .. math:
              log p(y_{*}|D) = log 1/num_samples prod^{S}_{s=1} p(y_{*}|f_{*s})
              f_{*s} ~ p(f_{*}|\mu_{*}\\sigma^{2}_{*})
  
          :param y_test: test observations (y_{*})
          :type y_test: (Nx1) array
          :param mu_star: predictive mean of gaussian p(f_{*}|mu_{*}, var_{*})
          :type mu_star: (Nx1) array
          :param var_star: predictive variance of gaussian p(f_{*}|mu_{*}, var_{*})
          :type var_star: (Nx1) array
          :param num_samples: num samples of p(f_{*}|mu_{*}, var_{*}) to take
          :type num_samples: int