<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def log_predictive_density_sampling(self, y_test, mu_star, var_star, Y_metadata=None, num_samples=1000):
 """ 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 """

        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