[R] Gram-Charlier series

[Martin Maechler]

>>>>> "AugS" ==   <Augusto.Sanabria at ga.gov.au>
>>>>>     on Wed, 22 Feb 2006 17:13:17 +1100 writes:

    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)[1]
    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)[1]
    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