# Issues using Hermite approximation to bivariate Gaussian in R

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 ?

-

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.

-
A question though (on other polynomials): my aim is to approximate a multivariate gaussian by polynomial expansions: which of the different approaches (Chebyshev, Hermite, Laguerre) is considered best for approximation ? –  user603 Sep 11 '10 at 7:47
I'm doesn't matter so much which polynomial class you use (up to numerical issues), but you can get better results by the choice of nodes where you sample the function (in this case the bivariate gaussian density) in order to compute the coefficients in front of the polynomials. For this problem, the Chebyshev nodes are optimal, at least when interpolating on a finite interval. see en.wikipedia.org/wiki/Chebyshev_nodes or any decent numerical analysis text (e.g. Cheney & Kincaid) or pages.cs.wisc.edu/~amos/412/lecture-notes/lecture09.pdf –  shabbychef Sep 12 '10 at 16:33
should say 'it doesn't matter'... –  shabbychef Sep 12 '10 at 17:20
Yes, indeed. But do the Chebyshev node generalize to multivariate functions ? I have tried to search the web and the only class of polynomials that is considered in the multivariate case seems to be the Hermite. –  user603 Sep 12 '10 at 23:58
ah, yes, for some reason I thought the function was separable. In higher dimensions, I only have experience in spline interpolation, e.g. for finite-element method compuations. (oddly enough, my thesis was on the topic). the big takeaway in planar meshing to minimize interpolation error is to minimize the maximum angle appearing in the triangulation. the margin does not permit me to be more detailed than this ;) you can look up any of the papers by J. Shewchuk (see 'finite element quality' at cs.cmu.edu/~jrs/jrspapers.html) or N. Walkington & G. Miller, etc. etc. hth. –  shabbychef Sep 13 '10 at 16:34