[18958] in s-news-athena
RE: [S] Shapiro-Wilks' test
daemon@ATHENA.MIT.EDU (Tom Gatliffe)
Tue Aug 17 23:28:25 1999
From: "Tom Gatliffe" <TomGofBoulder@email.msn.com>
To: "Ranjan Maitra" <maitra@math.umbc.edu>, "S-news" <s-news@wubios.wustl.edu>
Date: Tue, 17 Aug 1999 21:18:54 -0600
Message-Id: <LOBBLBJANBJEJGBDPLPCEEBDCBAA.TomGofBoulder@email.msn.com>
Mime-Version: 1.0
Content-Type: text/plain;
charset="iso-8859-1"
Content-Transfer-Encoding: 7bit
In-Reply-To: <Pine.LNX.4.04.9908171523130.5179-100000@param.math.umbc.edu>
The following is a function I pulled off the StatLib Library. It seems to
agree
with the Shapiro Wilk tests performed by a couple of other programs I have
used.
> shapiro.wilk.test
function(x)
{
# "shapiro.wilk.test"
#
# This function is an S version of the procedure described by
# J. P. Royston (1982) in "An Extension of Shapiro and Wilk's
# W Test for Normality to Large Samples", Applied Statistics,
# Vol. 31, No. 2, pp 115-124.
#
n <- length(x)
index <- 1:n
m <- qnorm((index - 0.375)/(n + 0.25))
y <- sort(x)
mu <- mean(y)
SSq <- sum((y - mu)^2)
astar <- 2 * m
ends <- c(1, n)
astar.p <- astar[ - ends]
if(n <= 20)
m <- n - 1
else m <- n
if(m < 20)
aa <- gamma(0.5 * (m + 1))/(sqrt(2) *
gamma(0.5 * m + 1))
else {
f1 <- (6 * m + 7)/(6 * m + 13)
f2 <- exp(1)/(m + 2)
f3 <- (m + 1)/(m + 2)
f3 <- f3^(m - 2)
aa <- f1 * sqrt(f2 * f3)
}
astar.1 <- (aa * sum(astar.p^2))/(1 - 2 * aa)
astar.1 <- sqrt(astar.1)
astar[1] <- - astar.1
astar[n] <- astar.1
A <- astar/sqrt(sum(astar^2))
W <- (sum(A * y)^2)/SSq
if(n <= 20) {
u <- log(n) - 3
lambda <- 0.118898 + 0.133414 * u +
0.327907 * u^2
logmu <- -0.37542 - 0.492145 * u -
1.124332 * u^2 - 0.199422 * u^3
logsigma <- -3.155805 + 0.729399 * u +
3.01855 * u^2 + 1.5558776 * u^3
}
if(n > 20) {
u <- log(n) - 5
lambda <- 0.480385 + 0.318828 * u +
0.0241665 * u^3 + 0.00879701 *
u^4 + 0.002989646 * u^5
logmu <- -1.91487 - 1.37888 * u -
0.04189209 * u^2 + 0.1066339 *
u^3 - 0.03513666 * u^4 -
0.01504614 * u^5
logsigma <- -3.83538 - 1.015807 * u - 0.331885 *
u^2 + 0.1773538 * u^3 - 0.01638782 * u^4 -
0.03215018 * u^5 + 0.003852646 * u^6
}
mu <- exp(logmu)
sigma <- exp(logsigma)
y <- (1 - W)^lambda
z <- (y - mu)/sigma
p <- 1 - pnorm(z)
if(n < 7) {
warning("n is to small for this program
to correctly estimate p"
)
p <- NA
}
if(n > 2000) {
warning("n is too large for this program
to correctly estimate p"
)
p <- NA
}
out <- list(W = W, n = n, p = p)
out
}
>
Tom Gatliffe
Rocky Flats Environmental Tech Site
-----Original Message-----
From: owner-s-news@wubios.wustl.edu
[mailto:owner-s-news@wubios.wustl.edu]On Behalf Of Ranjan Maitra
Sent: Tuesday, August 17, 1999 1:25 PM
To: S-news
Subject: [S] Shapiro-Wilks' test
Dear all,
Does anyone know if there is a way to perform the Shapiro-Wilks' test in
Splus?
Please post as a follow-up or I will summarize and post the responses.
Thanks,
Ranjan
***************************************************************************
Ranjan Maitra, Department of Mathematics and Statistics,
University of Maryland, Baltimore County, Baltimore, MD 21250, USA.
***************************************************************************
\|/ satyamevajayate | tamasomaajyotirgamaya \|/
-*- | -*-
/|\ TRUTH SHALL PREVAIL | FROM DARKNESS TO LIGHT /|\
***************************************************************************
Ph: 410-455-2436 FAX: 410-455-1066 http://www.math.umbc.edu/~maitra
***************************************************************************
-----------------------------------------------------------------------
This message was distributed by s-news@wubios.wustl.edu. To unsubscribe
send e-mail to s-news-request@wubios.wustl.edu with the BODY of the
message: unsubscribe s-news
-----------------------------------------------------------------------
This message was distributed by s-news@wubios.wustl.edu. To unsubscribe
send e-mail to s-news-request@wubios.wustl.edu with the BODY of the
message: unsubscribe s-news
