[R] PWM help

[J. R. M. Hosking]

On Wed, 25 Feb 2004, Jonathan Wang wrote:

 > I saw a Help e-mail related to MLE.  Does R have a probability weighted
 > method (PWM) estimator function?  I can't seem to find anything on PWM,
 > unless my eyes are playing trick on me.

R has no built-in facilities for PWMs, nor for L-moments (which are
generally preferable to PWMs -- see Hosking and Wallis, "Regional
Frequency Analysis", Cambridge Univ Press, 1997, sec. 2.4).

Here is an R/Splus implementation of a function for estimating
sample L-moments (routine SAMLMU from my LMOMENTS
Fortran package at http://lib.stat.cmu.edu/general/lmoments).
It has been lightly tested, but carries no guarantees of accuracy.

samlmu <- function (x, nmom = 4, sort.data = T)
{
    xok <- x[!is.na(x)]
    n <- length(xok)
    if (nmom <= 0) return(numeric(0))
    if (nmom <= 2) rnames <- paste("l", 1:nmom, sep = "_")
    else rnames <- c("l_1", "l_2", paste("t", 3:nmom, sep = "_"))
    lmom <- rep(NA, nmom)
    names(lmom) <- rnames
    if (n == 0) return(lmom)
    if (sort.data == T) xok <- sort(xok)
    nmom.actual <- min(nmom, n)
    lmom[1] <- mean(xok)
    if (nmom.actual == 1) return(lmom)
    temp <- seq(1-n, n-1, by = 2)
    p1 <- rep(1, n)
    p <- temp/(n-1)
    lmom[2] <- mean(xok * p)
    if (nmom.actual == 2) return(lmom)
    if (xok[1] == xok[n]) {
        warning("all data values equal")
        return(lmom)
    }
    for (j in 3:nmom.actual) {
        p2 <- p1
        p1 <- p
        p <- ((2*j-3)*temp*p1 - (j-2)*(n+j-2)*p2) / ((j-1)*(n-j+1))
        lmom[j] <- mean(xok * p)/lmom[2]
    }
    return(lmom)
}



J. R. M. Hosking
hosking at watson.ibm.com