# Re: [R] Gram-Charlier series

From: Martin Maechler <maechler_at_stat.math.ethz.ch>
Date: Wed 22 Feb 2006 - 19:46:32 EST

AugS> Good day everyone,

AugS> I want to use the Gram-Charlier series expansion to model     AugS> some data. To do that, I need functions to:

```    AugS> 1) Calculate 'n' moments from given data
AugS> 2) Transform 'n' moments to 'n' central moments, or
AugS> 3) Transform 'n' moments to 'n' cumulants
AugS> 4) Calculate a number of Hermite polynomials

```

AugS> Are there R-functions to do any of the above?

I have functions to do "4)".
The nicer ones are built on package 'polynom' (which you should definitely install and make use of for the above problems):

#### Hermite Polynomials --- An extension to the "polynom" package

```#### -------------------

```

## used e.g. from ../Pkg-ex/ctest/spearmanRho/prho-true-edgew.R
## for Edgeworth expansion

library(polynom)

hermitePolS <- function(n)
{

```  ## Purpose: n-th Hermite polynomial  He(n) -- as "polynom" object
## --------------------------------------------------------------
## Arguments: n >= 0: integer
## ----------------------------------------------------------------------
## Author: Martin Maechler, Date: 26 Apr 2003, 20:33
```
n <- as.integer(n)
if(n == 0) return(polynomial(1))
x <- polynomial(0:1)
if(n == 1) return(x)
```    ## else
## Recursion : He_n(x) =  x He_{n-1}(x) - (n-1) He_{n-2}(x)
## The following is *definitely* not efficient
```
return(x* hermitePolS(n-1) - (n-1) * hermitePolS(n-2)) }

system.time(He.9 <- hermitePolS(9)); He.9

## Much more efficient: without the *double* recursion :

hermitePol <- function(n)
{

```  ## Purpose: n-th Hermite polynomial  He(n) -- as "polynom" object
## --------------------------------------------------------------
## Arguments: n >= 0: integer
## ----------------------------------------------------------------------
## Author: Martin Maechler, Date: 26 Apr 2003, 21:02
```
n <- as.integer(n)
if(n == 0) return(polynomial(1))
x <- polynomial(0:1)
if(n == 1) return(x)
## else "Recursion" but the fast way:     He.n1 <- polynomial(1)
He <- x
for(nn in 2:n) {
```        He.n2 <- He.n1
He.n1 <- He
## Recursion : He_n(x) =  x He_{n-1}(x) - (n-1) He_{n-2}(x)
He <- x * He.n1 - (nn - 1) * He.n2
```
}
class(He) <- c("HermitePol", class(He))     return(He)
}

system.time(He9 <- hermitePol(9)); He9 # 9 x faster

(fH9 <- as.function(He9))## note that 'polynom' needs a fix
## i.e. polynom:::as.function.polynomial:
environment(fH9) <- .GlobalEnv
fH7 <- as.function(He7 <- hermitePol(7)) fH8 <- as.function(He8 <- hermitePol(8))
## Orthogonality and "Scale"
## These give 0 :

```integrate(function(x)fH8(x)*fH9(x) * dnorm(x), -Inf,Inf, rel.tol=1e-10)
integrate(function(x)fH7(x)*fH9(x) * dnorm(x), -Inf,Inf, rel.tol=1e-8)
integrate(function(x)fH7(x)*fH8(x) * dnorm(x), -Inf,Inf, rel.tol=1e-10)
```

## Scale: Inf{ He_n(x) ^ 2 phi(x) dx } = n! :
str(I9 <- integrate(function(x)fH9(x)^2 * dnorm(x), -Inf,Inf, rel.tol=1e-10, abs.tol = 0.01))

if(FALSE) {

## This is not quite it, but ok for small poly :
str.polynomial <- function(x, ...)

cat("polynomial", noquote(as.character(x)),"\n")

## to improve, I would want an option `highOrder' to
## as.character.polynomial(p, highOrder = FALSE)
## and then improve the following (use "..." when it's too long:

##- str.polynomial <- function(x, ...)
##- cat("polynomial", noquote(as.character(x, highOrder=TRUE)),"\n")

str(lapply(1:8, hermitePol))
## or even nicer:

}

for(n in 0:11)

cat("He(",formatC(n,wid=2),") = ", noquote(as.character(hermitePol(n))),

"\n", sep="")

```###---------------------------------------------------------------------

```

AugS> (mean, sd and cum3 are very limited)

AugS> Thank you for your help,

AugS> Augusto

```    AugS> --------------------------------------------
AugS> Augusto Sanabria. MSc, PhD.
.....................
(signature too long, core dumped)

```

Martin Maechler, ETH Zurich

R-help@stat.math.ethz.ch mailing list