Re: [R] Gram-Charlier series

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

>>>>> "AugS" == <Augusto.Sanabria@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



R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Wed Feb 22 19:55:52 2006

This archive was generated by hypermail 2.1.8 : Fri 03 Mar 2006 - 03:42:39 EST