def variational_expectation_quadrature(self, y, m, v, hyp=None, num_quad_points=20):
        """
        Computes the "variational expectation" via Gauss-Hermite quadrature, i.e. the
        expected log-likelihood, and its derivatives w.r.t. the posterior mean
            E[log p(yₙ|fₙ)] = log ∫ p(yₙ|fₙ) 𝓝(fₙ|mₙ,vₙ) dfₙ
        with EP power a.
        :param y: observed data (yₙ) [scalar]
        :param m: posterior mean (mₙ) [scalar]
        :param v: posterior variance (vₙ) [scalar]
        :param hyp: likelihood hyperparameter [scalar]
        :param num_quad_points: the number of Gauss-Hermite sigma points to use during quadrature [scalar]
        :return:
            exp_log_lik: the expected log likelihood, E[log p(yₙ|fₙ)]  [scalar]
            dE: first derivative of E[log p(yₙ|fₙ)] w.r.t. mₙ  [scalar]
            d2E: second derivative of E[log p(yₙ|fₙ)] w.r.t. mₙ  [scalar]
        """
        x, w = hermgauss(num_quad_points)  # Gauss-Hermite sigma points and weights
        w = w / np.sqrt(pi)  # scale weights by 1/√π
        sigma_points = np.sqrt(2) * np.sqrt(v) * x + m  # scale locations according to cavity dist.
        # pre-compute wᵢ log p(yₙ|xᵢ√(2vₙ) + mₙ)
        weighted_log_likelihood_eval = w * self.evaluate_log_likelihood(y, sigma_points, hyp)
        # Compute expected log likelihood via quadrature:
        # E[log p(yₙ|fₙ)] = ∫ log p(yₙ|fₙ) 𝓝(fₙ|mₙ,vₙ) dfₙ
        #                 ≈ ∑ᵢ wᵢ pᵃ(yₙ|xᵢ√(2vₙ) + mₙ)
        exp_log_lik = np.sum(
            weighted_log_likelihood_eval
        )
        # Compute first derivative via quadrature:
        # dE[log p(yₙ|fₙ)]/dmₙ = ∫ (fₙ-mₙ) vₙ⁻¹ log p(yₙ|fₙ) 𝓝(fₙ|mₙ,vₙ) dfₙ
        #                      ≈ ∑ᵢ wᵢ (fₙ-mₙ) vₙ⁻¹ log p(yₙ|xᵢ√(2vₙ) + mₙ)
        dE = np.sum(
            (sigma_points - m) / v
            * weighted_log_likelihood_eval
        )
        # Compute second derivative via quadrature:
        # d²E[log p(yₙ|fₙ)]/dmₙ² = ∫ [(fₙ-mₙ)² vₙ⁻² - vₙ⁻¹] log p(yₙ|fₙ) 𝓝(fₙ|mₙ,vₙ) dfₙ
        #                        ≈ ∑ᵢ wᵢ [(fₙ-mₙ)² vₙ⁻² - vₙ⁻¹] log p(yₙ|xᵢ√(2vₙ) + mₙ)
        d2E = np.sum(
            (0.5 * (v ** -2) * (sigma_points - m) ** 2 - 0.5 * v ** -1)
            * weighted_log_likelihood_eval
        )
        return exp_log_lik, dE, d2E