Issues using Hermite approximation to bivariate Gaussian in R

Question

In "On Gaussian-like Densities of Order Greater than Two" (Willett, P. Thomas, J. B., 1987), section II, the author state:

$\mathcal{N}(x,y,\rho)=\phi(x)\phi(y)\sum_{n=0}^{\infty}\rho^nH_n(x)H_n(y)$

where

1. $\phi(.)$ is the unit normal density,
2. $H_n(.)$ is the $n^{th}$ Hermite polynomial,
3. $\mathcal{N}(x,y,\rho)$ is the bivariate Gaussian density with correlation $\rho$.

I tried to repeat this (approximately, that is up to the 10th Hermite polynomial) in R:

library(PolynomF)
x<-polynom()
H<-polylist(1,x);for(n in 2:10) H[[n+1]]<-x*H[[n]]-(n-1)*H[[n-1]]
Hp<-as.function(H)
#The Hermite polynomial of order 10
rho<-0.7
R<-c();for(n in 0:10) R[[n+1]]<-rho^n

HPA<-function(z) prod(dnorm(z))*sum(Hp(z[1])*Hp(z[2])*R)

z<-runif(2)
HPA(z)

Which is completely off the mark (i.e. for $z=(0.65,0.63)$ i get 1194). What am i doing wrong ?

Answer 1

It would appear that you neglected to normalize your Hermite polynomials. Try this:

library(PolynomF)
x <- polynom()
H <- polylist(1,x); for(n in 2:10) H[[n+1]] <- x*H[[n]] - (n-1)*H[[n-1]]
for(n in 1:11) H[[n]] <- H[[n]]*exp(-lgamma(n)/2)
Hp <- as.function(H)
#The (normalized) Hermite polynomial of order 10
rho <- 0.7
R<-c(); for(n in 0:10) R[[n+1]] <- rho^n

HPA <- function(z) prod(dnorm(z))*sum(Hp(z[1])*Hp(z[2])*R)

set.seed(1)
z <- runif(2)
HPA(z)

Now have a look at what it should be:

library(mvtnorm)
sigma <- matrix(c(1, rho, rho, 1), ncol = 2)
dmvnorm(z, sigma = sigma)

which is only off in the thousandths place on my machine.

Also, you may be interested in the orthopolynom package for another way to generate normalized Hermite polynomials.