# Instability with Gauss-Hermite quadrature

**URL:** <https://discourse.julialang.org/t/instability-with-gauss-hermite-quadrature/103408>\
**Category:** Numerics\
**Tags:** question, quadrature\
**Created:** [August 31, 2023, 2:48pm UTC](https://discourse.julialang.org/t/instability-with-gauss-hermite-quadrature/103408 "2023-08-31T14:48:10Z")\
**Posts on this page:** 1\
**Page:** 1

<div class="post-metadata">

**Author:** ![mleprovost](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mleprovost/32/7166_2.png) [@mleprovost](https://discourse.julialang.org/u/mleprovost)\
**Post date:** [August 31, 2023, 2:48pm UTC](https://discourse.julialang.org/t/instability-with-gauss-hermite-quadrature/103408/1 "2023-08-31T14:48:10Z")

</div>

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:

```julia
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")

```