Re: [S] Tests for Normality

David J Cummins (CUMMINS_DAVID_J@lilly.com)
Thu, 05 Nov 1998 20:26:29 -0500


Viswanath Devanarayan and I have coded up a method developed by
Brown and Hettmansperger and published in JASA 1996, vol 91, No. 436, pp.
1668-1675.
The main motivation of their work is that the standard tests for normality
(e.g., Shapiro- Wilk)
usually test for overall normality and do not indicate the manner of the
departure from normality.
There are three main features of the B&H method:
1. It tests for skewness, kurtosis, and normal tails separately
2. It is more powerful than the Shapiro-Wilk test
3. It adjusts for small sample bias (via simulation)

To try out the attached code, just open a graphics window and then source the
program.
Note that there is a call to motif() in one place. For windows users, just
change this to win.graph().

# ----------------------------------------------------------------------------
--------------------------------------------------
###################################################################
# Program name: BHtests_version2.s
# Updated: December 5, 1997.
# Authors: Dave Cummins, V. Devanarayan
# Eli Lilly & Company -- Cummins_DJ@lilly.com
# Purpose:
# Splus functions for implementing the Brown-Hettmansperger tests
# for skewness, kurtosis and normal tails (see JASA, Dec '96).
# For validation purposes, data from their paper have been included
# in the variable my.y in the bottom part of this code. In
# order to use this program for other data sets, simply call the
# function: BHtests(my.data) after sourcing the file: BHtests_version2.s
####################################################################

options(digits=15)

# estimate of sigma

est.sigma _ function(y)
{
y _ sort(y)
n _ length(y)
vec1 _ 1:(n-1)
vec2 _ 2:n

sum(
( y[vec2] - y[vec1] )*dnorm( qnorm( vec1/n ) )
)
}

# estimate of d3

est.d3 _ function(y)
{
y _ sort(y)
n _ length(y)
vec1 _ 1:(n-2)
vec2 _ 2:(n-1)

d3hat _
2^(-.5)*y[1]*( -dnorm(qnorm(1/n))*qnorm(1/n) ) +
2^(-.5)*sum(
y[vec2]*(
dnorm( qnorm( vec1/n ) )*qnorm( vec1/n ) -
dnorm( qnorm( vec2/n ) )*qnorm( (2:(n-1))/n )
)
) +
2^(-.5)*y[n]*( dnorm(qnorm((n-1)/n))*qnorm((n-1)/n) )
d3hat
}

# estimate of d4

est.d4 _ function(y)
{
y _ sort(y)
n _ length(y)
vec1 _ 1:(n-2)
vec2 _ 2:(n-1)

d4hat _
6^(-.5)*y[1]*( -dnorm(qnorm(1/n))*(qnorm(1/n)^2 - 1) ) +
6^(-.5)*sum(
y[vec2]*(
dnorm(qnorm(vec1/n))*(qnorm(vec1/n)^2 - 1) -
dnorm(qnorm(vec2/n))*(qnorm(vec2/n)^2 - 1)
)
) +
6^(-.5)*y[n]*( dnorm(qnorm((n-1)/n))*(qnorm((n-1)/n)^2 - 1) )
d4hat
}

# Statistic Z1 (not adjusted for small sample bias)

z1 _ function(y) {
sqrt(3*length(y))*est.d3(y)/est.sigma(y)
}

# Statistic Z2 (not adjusted for small sample bias)

z2 _ function(y) {
sqrt(4*length(y))*est.d4(y)/est.sigma(y)
}

# This function computes bias adjustment corrections for small samples
# n: sample size, r: number of replications

biasadj _ function(n,r) {
z1a _ matrix(0,r,1)
z2a _ matrix(0,r,1)
# The set.seed below will give results that match the JASA paper
# set.seed(158)

cat("Iteration:")
for (i in 1:r)
{
if( trunc(i/100)==i/100 ) cat(i,".. ")
ry _ rnorm(n)
z1a[i] _ z1(ry)
z2a[i] _ z2(ry)
}
cat("\n")
cat("\n")
c( mean(z2a), sqrt(var(z1a)), sqrt(var(z2a)) )
}

# This is the function that actually computes the test-statistics
# and p-values for skewness, kurtosis and nonnormal tails.
# This function calls the other functions defined above.

BHtests <- function(y,r=1000,plotit=T,legend.size=1)
{
###################################################################
# Program name: BHtests_version2.s
# Updated: December 5, 1997.
# Authors: Dave Cummins, V. Devanarayan
# Purpose:
# Splus functions for implementing the Brown-Hettmansperger tests
# for skewness, kurtosis and normal tails (see JASA, Dec '96).
# For validation purposes, data from their paper have been included
# in the variable my.y in the bottom part of this program. In
# order to use this program for other data sets, simply call the
# function: BHtests(my.data) after sourcing the file: BHtests_version2.s
###################################################################

y <- sort(y)
buv _ biasadj(n=length(y),r)
bn _ buv[1]
un _ buv[2]
vn _ buv[3]

Z1ab _ sqrt(3*length(y))*est.d3(y)/(un*est.sigma(y))
pZ1 _ 2*( 1-pnorm(abs(Z1ab)) )
Z2ab _ max( (sqrt(4*length(y))*est.d4(y)/est.sigma(y) - bn)/vn, 0 )
pZ2 _ 1-pnorm(Z2ab)
Tab _ Z1ab^2 + Z2ab^2
pT _ 1-pnorm(sqrt(Tab)) + .5*exp(-Tab/2)

Z1ab <- round(Z1ab,5)
pZ1 <- round(pZ1,5)
Z2ab <- round(Z2ab,5)
pZ2 <- round(pZ2,5)
Tab <- round(Tab,5)
pT <- round(pT,5)
cat("Bias Adjusted Test Statistic for Skewness =",Z1ab,'\n')
cat("with asymptotic P-value for a two-sided test =",pZ1,'\n','\n')
cat("Bias Adjusted Test Statistic for Kurtosis =",Z2ab,'\n')
cat("with asymptotic P-value for a one-sided test =",pZ2,'\n','\n')
cat("Bias Adjusted Test Statistic for Normal Tails =",Tab,'\n')
cat("with asymptotic P-value for a two-sided test =",pT,'\n','\n','\n')

if(plotit==T)
{
if(names(dev.cur())=="null device") motif()
n <- length(y)
vec <- seq(from=0,to=n,by=1)
vec[(n+1)] <- 0
phi <- dnorm(qnorm(vec/n))
b <- phi[1:n] - phi[2:(n+1)]
lambda <- t(b) %*% b
a <- b / lambda # equivalent to : solve( t(b)%*%b ) %*% b
plot(a,y, ylab="Observed Data",
xlab="Quantiles of Standard Normal (bias-adj)")

if(legend.size>0)
{
y.width <- max(y) - min(y)
a.width <- max(a) - min(a)
text.top <- max(y) - .05 * y.width*legend.size
text.bott <- max(y) - .40*y.width*legend.size
text.left <- min(a) + .05 * a.width*legend.size
text.right <- min(a) + .5 * a.width*legend.size
lines(a,a*sqrt(var(y))+mean(y))
legend(x=c(text.left,text.right),y=c(text.top,text.bott),
type="n",bty="n",cex=legend.size,
legend=c(" p-values ",
" ------------- ",
paste("Skewness: ",pZ1,sep="")," ",
paste("Kurtosis: ",pZ2,sep="")," ",
paste("Normal Tails: ",pT,sep=""))
)
}
}
list(skewness=pZ1,kurtosis=pZ2,tails=pT)
}

### Here is where we enter the data in the variable my.y
### and invoke the function BHtests() to perform the tests for
### skewness, kurtosis and nonnormal tails.
my.y _ c(16.1, 18.8, 15.1, 16, 12, 16.8, 18.7, 23.4, 17, 19.9,
18.2, 17.5, 17.7, 16.8) # example data given in the paper

BHtests(my.y,r=500)

# ----------------------------------------------------------------------------
--------------------------------------------------

"sinclair, andrew (MIP London)" <asinclair@edfman.co.uk> on 11/05/98 01:48:00
PM

To: s-news <s-news@wubios.wustl.edu>
cc:

Subject: [S] Tests for Normality

Does anyone know of code/ references for tests of normality ?

I am aware of qqnorm plots, Kolmogorov Smirnoff and chi squared
goodness of fit tests. I need procedures that give a specific p-values,
and preferably are sensitive to alternatives with fat tails (which is not
true for K-S in my experience).

Many thanks for any advice - I will summarise responses that
are sent to me only.

Andrew

____________________

Andrew Sinclair
ED & F Man Investment Products
Sugar Quay
Lower Thames Street
London EC3R 6DU

tel : 0171 285 2089
fax : 0171 285 2001
email : asinclair@edfman.co.uk
-----------------------------------------------------------------------
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