[18962] in s-news-athena
Re: [S] Shapiro-Wilks' test
daemon@ATHENA.MIT.EDU (Peter B. Mandeville)
Wed Aug 18 08:43:28 1999
Message-Id: <3.0.5.16.19990818041614.1e2745dc@deimos.tc.uaslp.mx>
Date: Wed, 18 Aug 1999 04:16:14
To: S-news <s-news@wubios.wustl.edu>
From: "Peter B. Mandeville" <mandevip@uaslp.mx>
In-Reply-To: <Pine.LNX.4.04.9908171523130.5179-100000@param.math.umbc.ed
u>
Mime-Version: 1.0
Content-Type: text/plain; charset="iso-8859-1"
Content-Transfer-Encoding: 8bit
I coded the Shapiro-Francia test which handles from 5 to 5000 observations
from
Patrick Royston
1993 A Pocket-Calculator Algorithm for the Shapiro-Francia Test for
Non-Normality: An Application to Medicine. Statistics in Medicine vol
12, 181-184
1991 Estimating Departure from Normality. Statistics in Medicine.
vol 10, 1283-1293
1983 A Simple Method for Evaluating the Shapiro-Francia W' Test of
Non-Normality. The Statistician 32 (1983) 297-300
and the D'Agostino tests from
Ralph B. D'Agostino, Albert Belanger, and Ralph B. D'Agostino, Jr.
1990 A Suggestion for Using Powerful and Informative Tests of Normality.
The American Statistician, November 1990, Vol. 44, No. 4
for testing normality in R. They should work in S plus with the exception
of hist which I believe is specific to R but which I believe is available
in MASS from Venables y Ripley 1997.
# NORMALIDAD
# Peter
NORMAL <- function(data){
cat("NORMALIDAD\n")
DAGOSTINO(data)
SHAPIRO.FRANCIA(data)
EDA.SHAPE(data)
}
# Mathsoft 1998:45
EDA.SHAPE <- function(data){
dest <- density(data)
hist(data,xlim=range(dest$x),xlab="x",ylab="density",freq=F)
lines(dest,lty=2,col="red")
boxplot(data,main="Diagrama de Caja")
qqnorm(data)
qqline(data)
SYMPLOT(data)
}
# Lee 1996:33-34
MOMENTS <- function(data,r) sum((data-mean(data))^r)/length(data)
# Peter
SKEW <- function(data) MOMENTS(data,3)/(MOMENTS(data,2)*sqrt(MOMENTS(data,2)))
# Peter
KURTOSIS <- function(data) MOMENTS(data,4)/(MOMENTS(data,2)*MOMENTS(data,2))
# Peter
OUTLIERS <- function(data){
cat("OUTLIERS\n")
x <- FIVE(data)
n <- length(data)
cat(" CI PARA OUTLIERS:",x[1],"\n")
for(j in 1:n)
if(data[j]<x[1]) cat(" Obs",j,":",data[j],"\n")
cat(" CS PARA OUTLIERS:",x[5],"\n")
for(j in 1:n)
if(data[j]>x[5]) cat(" Obs",j,":",data[j],"\n")
}
# D'Agostino, Belanger and D'Agostino 1990:316-321
DAGOSTINO <- function(data){
cat(" PRUEBAS DE D'AGOSTINO\n")
n <- length(data)
cat(" SIMETRIA\n")
cat(" COEFICIENTE DE SIMETRIA:",sqrtb1 <- SKEW(data),"\n")
if(n>8){
y <- sqrtb1*sqrt((n+1)*(n+3)/(6*(n-2)))
beta2 <- 3*(n*n+27*n-70)*(n+1)*(n+3)/((n-2)*(n+5)*(n+7)*(n+9))
w <- sqrt(-1+sqrt(2*(beta2-1)))
delta <- 1/sqrt(log(w))
ALPHA <- sqrt(2/(w*w-1))
cat(" ZETA CALCULADA:",zb1 <-
delta*log(y/ALPHA+sqrt((y/ALPHA)^2+1)),"\n")
cat(" PROBABILIDAD:",2*(1-pnorm(abs(zb1))),"\n")
}else
cat(" LA PRUEBA DE SESGO REQUIERE POR LO MENOS 9 REPETICIONES\n")
cat(" KURTOSIS\n")
cat(" COEFICIENTE DE KURTOSIS:",b2 <- KURTOSIS(data),"\n")
if(n>19){
meanb2 <- 3*(n-1)/(n+1)
varb2 <- 24*n*(n-2)*(n-3)/((n+1)*(n+1)*(n+3)*(n+5))
x <- (b2-meanb2)/sqrt(varb2)
moment <-
6*(n*n-5*n+2)/((n+7)*(n+9))*sqrt(6*(n+3)*(n+5)/(n*(n-2)*(n-3)))
a <- 5+8/moment*(2/moment+sqrt(1+4/(moment*moment)))
cat(" ZETA CALCULADA:",zb2 <-
(1-2/(9*a)-((1-2/a)/(1+x*sqrt(2/(a-4))))^(1/3))/sqrt(2/(9*a)),"\n")
cat(" PROBABILIDAD:",2*(1-pnorm(abs(zb2))),"\n")
cat(" OMNIBUS\n")
cat(" JI-CUADRADA CALCULADA:",k2 <- zb1*zb1+zb2*zb2,"\n")
cat(" GRADOS DE LIBERTAD: 2\n")
cat(" PROBABILIDAD:",probji2 <- 1-pchisq(k2,2),"\n")
}else
cat(" LAS PRUEBAS DE KURTOSIS Y OMNIBUS REQUIEREN POR LO MENOS 20
REPETICIONES\n")
}
# Royston 1993:183-184
SHAPIRO.FRANCIA <- function(data){
cat(" PRUEBA DE NORMALIDAD DE SHAPIRO-FRANCIA\n")
n <- length(data)
if(n<5 || n>5000)
cat(" REQUIERE ENTRE 5 Y 5000 REPETICIONES\n")
else{
xbar <- mean(data)
sdata <- sort(data)
resid <- sdata-xbar
uniform <- seq(1,n)
np <- qnorm((uniform-0.375)/(n+0.25))
cat(" W':",w <-
(sum(np*sdata))^2/(sum(np*np)*sum(resid*resid)),"\n")
u <- log(n)
v <- log(u)
muy <- -1.2725+1.0521*(v-u)
sigmay <- 1.0308-0.26758*(v+2/u)
cat(" ZETA CALCULADA:",zeta <- (log(1-w)-muy)/sigmay,"\n")
cat(" PROBABILIDAD:",probz <- 1-pnorm(zeta),"\n")
}
}
Peter B.
At 03:24 PM 17/08/99 -0400, you wrote:
>
>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
>
>
--
Peter B. Mandeville mandevip@deimos.tc.uaslp.mx
Jefe del Depto. de Informática y Bioestadística rpe1531@pasteur.fmed.uaslp.mx
Facultad de Medicine Tel: 48 26-23-45 ext. 232
Universidad Autónoma de San Luis Potosí Fax: 48 28-23-52
Av. V. Carranza 2405
Col. Los Filtros
Apartado Postal 145
San Luis Potosí, S.L.P.
78210 México
-----------------------------------------------------------------------
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
