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
[EMAIL PROTECTED]
______________________________________________
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
