Gram-Charlier series

classic Classic list List threaded Threaded
4 messages Options
Reply | Threaded
Open this post in threaded view
|

Gram-Charlier series

Augusto.Sanabria

 Good day everyone,

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

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

 Are there R-functions to do any of the above?
 (mean, sd and cum3 are very limited)

 Thank you for your help,

 Augusto

--------------------------------------------
Augusto Sanabria. MSc, PhD.
Mathematical Modeller
Risk Research Group
Geospatial & Earth Monitoring Division
Geoscience Australia (www.ga.gov.au)
Cnr. Jerrabomberra Av. & Hindmarsh Dr.
Symonston ACT 2609
Ph. (02) 6249-9155

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Reply | Threaded
Open this post in threaded view
|

Re: Gram-Charlier series

Martin Maechler
>>>>> "AugS" ==   <[hidden email]>
>>>>>     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

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Reply | Threaded
Open this post in threaded view
|

Re: Gram-Charlier series

Ernst Hansen
Martin Maechler writes:
 > >>>>> "AugS" ==   <[hidden email]>
 > >>>>>     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)".

For a direct way to translate raw moments into cumulants, you may want
to look at the methods from the book 'Symbolic Computation for
Statistical Inference' by D. F. Andrews and J. E. Stafford.  There is
in fact an R-implementation of the methods available at

  http://fisher.utstat.toronto.edu/david/SCSI/RSCSI.html


Hope this helps,

Ernst Hansen
Department of Statistics
University of Copenhagen

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Reply | Threaded
Open this post in threaded view
|

Re: Gram-Charlier series

Janusz Kawczak
I believe this page includes more up to date version:
http://fisher.utstat.toronto.edu/david/Sym2004/Sym2004.html

Janusz.

On Wed, 22 Feb 2006, Ernst Hansen wrote:

> Martin Maechler writes:
>  > >>>>> "AugS" ==   <[hidden email]>
>  > >>>>>     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)".
>
> For a direct way to translate raw moments into cumulants, you may want
> to look at the methods from the book 'Symbolic Computation for
> Statistical Inference' by D. F. Andrews and J. E. Stafford.  There is
> in fact an R-implementation of the methods available at
>
>   http://fisher.utstat.toronto.edu/david/SCSI/RSCSI.html
>
>
> Hope this helps,
>
> Ernst Hansen
> Department of Statistics
> University of Copenhagen
>
> ______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
>

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html