Instability with Gauss-Hermite quadrature


I would like to approximate the log-density of a univariate distribution \pi_0 with a basis of Hermite polynomials. To compute the coefficients of the expansion, I am using the Gauss-Hermite quadrature. I am getting unstable predictions.

Minimal working example:

using LinearAlgebra
using FastGaussQuadrature
using SpecialFunctions
using SpecialPolynomials
using Distributions

Nnode = 1000
nodes_phy, weights_phy = gausshermite(Nnode)

# Build a bimodal distribution
π0 = MixtureModel(Normal[
   Normal(-1.0, 0.2),
   Normal(1.0, 0.2)], [0.5, 0.5])

# Check that the normalized Hermite polynomials have unit norm
order_max = 50
coeff_normalization = zeros(order_max + 1)
C(n) = 1/sqrt(2^n*gamma(n+1)*√π)

for n=0:order_max
    Hn = C(n)*basis(Hermite, n)
    coeff_normalization[n+1] = dot(weights_phy, Hn.(nodes_phy).^2)
isapprox(coeff_normalization, ones(order_max+1), atol = 1e-10)

# Project the log(π0) onto the Physicist Hermite Polynomials

coeff = zeros(order_max + 1)

for n=0:order_max
    Hn = C(n)*basis(Hermite, n)
    # Compute ∫logπ0(x)Hₙ(x)exp(-x^2)dx ≈ ∑_i wᵢ log π0(xᵢ) Hₙ(xᵢ)
    coeff[n+1] = dot(weights_phy, map(x -> logpdf(π0, x)*Hn(x), nodes_phy))

# Routine to evaluate a Hermite expansion
function evaluate_hermite_poly(coeff, x)
    order_max = size(coeff,1)-1
    out = 0.0
    for n=0:order_max
        Hn = C(n)*basis(Hermite, n)
        out += coeff[n+1]*Hn(x)
    return out

# Plot the results
xgrid = -10:0.01:10
plot(xgrid, pdf.(π0, xgrid), label = L"\log\; \pi_0")
plot!(xgrid, map(x-> exp(evaluate_hermite_poly(coeff, x)), xgrid), label = "Hermite approximation")