Hello,
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)
end
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))
end
# 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)
end
return out
end
# 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")