<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def ndiagquad(funcs, H: int, Fmu, Fvar, logspace: bool=False, **Ys):
 """ Computes N Gaussian expectation integrals of one or more functions using Gauss-Hermite quadrature. The Gaussians must be independent. :param funcs: the integrand(s):
 Callable or Iterable of Callables that operates elementwise, on the following arguments: - `Din` positional arguments to match Fmu and Fvar; i.e., 1 if Fmu and Fvar are tensors; - the same keyword arguments as given by **Ys All arguments will be tensors of shape (N, 1) :param H: number of Gauss-Hermite quadrature points :param Fmu: array/tensor or `Din`-tuple/list thereof :param Fvar: array/tensor or `Din`-tuple/list thereof :param logspace: if True, funcs are the log-integrands and this calculates the log-expectation of exp(funcs) :param **Ys: arrays/tensors; deterministic arguments to be passed by name Fmu, Fvar, Ys should all have same shape, with overall size `N` (i.e., shape (N,) or (N, 1)) :return: shape is the same as that of the first Fmu Example use-cases: Fmu, Fvar are mean and variance of the latent GP, can be shape (N, 1) or (N,) m1, m2 are 'scalar' functions of a single argument F, broadcasting over arrays Em1, Em2 = ndiagquad([m1, m2], 50, Fmu, Fvar) calculates Em1 = ∫ m1(F) N(F; Fmu, Fvar) dF and Em2 = ∫ m2(F) N(F; Fmu, Fvar) dF for each of the elements of Fmu and Fvar. Em1 and Em2 have the same shape as Fmu. logp is a 'scalar' function of F and Y Y are the observations, with shape (N,) or (N, 1) with same length as Fmu and Fvar Ev = ndiagquad(logp, 50, Fmu, Fvar, Y=Y) calculates Ev = ∫ logp(F, Y) N(F; Fmu, Fvar) dF (variational expectations) for each of the elements of Y, Fmu and Fvar. Ev has the same shape as Fmu. Ep = ndiagquad(logp, 50, Fmu, Fvar, logspace=True, Y=Y) calculates Ep = log ∫ exp(logp(F, Y)) N(F; Fmu, Fvar) dF (predictive density) for each of the elements of Y, Fmu and Fvar. Ep has the same shape as Fmu. Heteroskedastic likelihoods: g1, g2 are now functions of both F and G logp is a function of F, G and Y Gmu, Gvar are mean and variance of a different GP controlling the variance Em = ndiagquad(m1, 50, Fmu, Fvar) -> Em1 = ∫∫ m1(F, G) N(F; Fmu, Fvar) N(G; Gmu, Gvar) dF dG Ev = ndiagquad(logp, 50, Fmu, Fvar, Y=Y) -> Ev = ∫∫ logp(F, G, Y) N(F; Fmu, Fvar) N(G; Gmu, Gvar) dF dG (variational expectations) Ep = ndiagquad(logp, 50, Fmu, Fvar, logspace=True, Y=Y) -> Ep = log ∫∫ exp(logp(F, G, Y)) N(F; Fmu, Fvar) N(G; Gmu, Gvar) dF dG (predictive density) """

    def unify(f_list):
        """
        Stack a list of means/vars into a full block
        """
        return tf.reshape(
                tf.concat([tf.reshape(f, (-1, 1)) for f in f_list], axis=1),
                (-1, 1, Din))

    if isinstance(Fmu, (tuple, list)):
        Din = len(Fmu)
        shape = tf.shape(Fmu[0])
        Fmu, Fvar = map(unify, [Fmu, Fvar])    # both N x 1 x Din
    else:
        Din = 1
        shape = tf.shape(Fmu)
        Fmu, Fvar = [tf.reshape(f, (-1, 1, 1)) for f in [Fmu, Fvar]]

    xn, wn = mvhermgauss(H, Din)
    # xn: H**Din x Din, wn: H**Din

    gh_x = xn.reshape(1, -1, Din)             # 1 x H**Din x Din
    Xall = gh_x * tf.sqrt(2.0 * Fvar) + Fmu   # N x H**Din x Din
    Xs = [Xall[:, :, i] for i in range(Din)]  # N x H**Din  each

    gh_w = wn * np.pi ** (-0.5 * Din)  # H**Din x 1

    for name, Y in Ys.items():
        Y = tf.reshape(Y, (-1, 1))
        Y = tf.tile(Y, [1, H**Din])  # broadcast Y to match X
        # without the tiling, some calls such as tf.where() (in bernoulli) fail
        Ys[name] = Y  # now N x H**Din

    def eval_func(f):
        feval = f(*Xs, **Ys)  # f should be elementwise: return shape N x H**Din
        if logspace:
            log_gh_w = np.log(gh_w.reshape(1, -1))
            result = tf.reduce_logsumexp(feval + log_gh_w, axis=1)
        else:
            result = tf.matmul(feval, gh_w.reshape(-1, 1))
        return tf.reshape(result, shape)

    if isinstance(funcs, Iterable):
        return [eval_func(f) for f in funcs]
    else:
        return eval_func(funcs)