Problemas con la aproximación de Hermite a la gaussiana bivariante en R

Question

En "Sobre las densidades gaussianas de orden superior a dos" (Willett, P. Thomas, J. B., 1987), sección II, el autor afirma:

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

donde

1. $\phi(.)$ es la densidad normal unitaria,
2. $H_n(.)$ es el $n^{th}$ Polinomio de Hermite,
3. $\mathcal{N}(x,y,\rho)$ es la densidad gaussiana bivariada con correlación $\rho$ .

He intentado repetirlo (aproximadamente, es decir, hasta el décimo polinomio de Hermite) en 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)

Lo cual está completamente fuera de lugar (es decir, para $z=(0.65,0.63)$ i obtener 1194). ¿Qué estoy haciendo mal?

Answer 1

Parece que te has olvidado de normalizar tus polinomios de Hermite. Intenta esto:

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)

Ahora echa un vistazo a lo que debería ser:

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

que en mi máquina sólo está desactivado en las milésimas.

Además, puede interesarle el orthopolynom para otra forma de generar polinomios de Hermite normalizados.