<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def predict_moments(self, X, nsamples=200, likelihood_args=()):
 r""" Predictive moments, in particular mean and variance, of a Bayesian GLM. This function uses Monte-Carlo sampling to evaluate the predictive mean and variance of a Bayesian GLM. The exact expressions evaluated are, .. math :: \mathbb{E}[y^* | \mathbf{x^*}, \mathbf{X}, y] &= \int \mathbb{E}[y^* | \mathbf{w}, \phi(\mathbf{x}^*)] p(\mathbf{w} | \mathbf{y}, \boldsymbol\Phi) d\mathbf{w}, \mathbb{V}[y^* | \mathbf{x^*}, \mathbf{X}, y] &= \int \left(\mathbb{E}[y^* | \mathbf{w}, \phi(\mathbf{x}^*)] - \mathbb{E}[y^* | \mathbf{x^*}, \mathbf{X}, y]\right)^2 p(\mathbf{w} | \mathbf{y}, \boldsymbol\Phi) d\mathbf{w}, where :math:`\mathbb{E}[y^* | \mathbf{w}, \phi(\mathbf{x}^*)]` is the the expected value of :math:`y^*` from the likelihood, and :math:`p(\mathbf{w} | \mathbf{y}, \boldsymbol\Phi)` is the posterior distribution over weights (from ``learn``). Here are few concrete examples of how we can use these values, - Gaussian likelihood: these are just the predicted mean and variance, see ``revrand.regression.predict`` - Bernoulli likelihood: The expected value is the probability, :math:`p(y^* = 1)`, i.e. the probability of class one. The variance may not be so useful. - Poisson likelihood: The expected value is similar conceptually to the Gaussian case, and is also a *continuous* value. The median (50% quantile) from ``predict_interval`` is a discrete value. Again, the variance in this instance may not be so useful. Parameters X : ndarray (N*,d) array query input dataset (N* samples, d dimensions). nsamples : int, optional Number of samples for sampling the expected moments from the predictive distribution. likelihood_args : sequence, optional sequence of arguments to pass to the likelihood function. These are non-learnable parameters. They can be scalars or arrays of length N. Returns ------- Ey : ndarray The expected value of y* for the query inputs, X* of shape (N*,). Vy : ndarray The expected variance of y* (excluding likelihood noise terms) for the query inputs, X* of shape (N*,). """

        # Get latent function samples
        N = X.shape[0]
        ys = np.empty((N, nsamples))
        fsamples = self._sample_func(X, nsamples)

        # Push samples though likelihood expected value
        Eyargs = tuple(chain(atleast_list(self.like_hypers_), likelihood_args))
        for i, f in enumerate(fsamples):
            ys[:, i] = self.likelihood.Ey(f, *Eyargs)

        # Average transformed samples (MC integration)
        Ey = ys.mean(axis=1)
        Vy = ((ys - Ey[:, np.newaxis])**2).mean(axis=1)

        return Ey, Vy