# Gram-Charlier series

4 messages
Open this post in threaded view
|

## Gram-Charlier series

 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-helpPLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Open this post in threaded view
|

## Re: Gram-Charlier series

 >>>>> "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-helpPLEASE do read the posting guide! http://www.R-project.org/posting-guide.html