def METHOD_NAME(self, cuda=False):
    func = lambda x: torch.sin(x)
    means = torch.randn(10)
    variances = torch.randn(10).abs()
    quadrature = GaussHermiteQuadrature1D()
    if cuda:
        means = means.cuda()
        variances = variances.cuda()
        quadrature = quadrature.cuda()
    dist = MultivariateNormal(means, DiagLinearOperator(variances.sqrt()))
    # Use quadrature
    results = quadrature(func, dist)
    # Use Monte-Carlo
    samples = dist.rsample(torch.Size([20000]))
    actual = func(samples).mean(0)
    self.assertLess(torch.mean(torch.abs(actual - results)), 0.1)