>>>>> "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