def statistical_linear_regression_quadrature(self, m, v, hyp=None, num_quad_points=20):
        """
        Perform statistical linear regression (SLR) using Gauss-Hermite quadrature.
        We aim to find a likelihood approximation p(yₙ|fₙ) ≈ 𝓝(yₙ|Afₙ+b,Ω+Var[yₙ|fₙ]).
        """
        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.
        lik_expectation, _ = self.conditional_moments(sigma_points, hyp)
        # Compute zₙ via quadrature:
        # zₙ = ∫ E[yₙ|fₙ] 𝓝(fₙ|mₙ,vₙ) dfₙ
        #    ≈ ∑ᵢ wᵢ E[yₙ|xᵢ√(2vₙ) + mₙ]
        z = np.sum(
            w * lik_expectation
        )
        # Compute variance S via quadrature:
        # S = ∫ (E[yₙ|fₙ]-zₙ) (E[yₙ|fₙ]-zₙ)' 𝓝(fₙ|mₙ,vₙ) dfₙ
        #   ≈ ∑ᵢ wᵢ (E[yₙ|xᵢ√(2vₙ) + mₙ]-zₙ) (E[yₙ|xᵢ√(2vₙ) + mₙ]-zₙ)'
        S = np.sum(
            w * (lik_expectation - z) * (lik_expectation - z)
        )
        # Compute cross covariance C via quadrature:
        # C = ∫ (fₙ-mₙ) (E[yₙ|fₙ]-zₙ)' 𝓝(fₙ|mₙ,vₙ) dfₙ
        #   ≈ ∑ᵢ wᵢ (fₙ-mₙ) (E[yₙ|xᵢ√(2vₙ) + mₙ]-zₙ)'
        C = np.sum(
            w * (sigma_points - m) * (lik_expectation - z)
        )
        # compute likelihood approximation 𝓝(yₙ|Afₙ+b,Ω+Var[yₙ|fₙ])
        A = C * v**-1  # the scale
        b = z - A * m  # the offset
        omega = S - A * v * A  # the linearisation error
        return A, b, omega