"sir" <- function(y,x,H){
# Performs the Sliced Inverse Regression analysis
# y: the data matrix, T by 1, of the dependent variable
# x: the T-by-p matrix of predictors.
# H: the number of slices.
#
z=as.matrix(x)
p=ncol(z)
T=nrow(z)
m1=sort(y,index.return=T)
idx=m1$ix
if(H < 5)H=5
m=T/H
nH=rep(m,H)
rem = T-m*H
if(rem > 0){
for (i in 1:rem){
nH[i]=nH[i]+1
}
}
#print(nH)
xmean=apply(z,2,mean)
ones=matrix(rep(1,T),T,1)
xm=matrix(xmean,1,p)
Xdev=z-ones%*%xm
varX=cov(Xdev)
m1=eigen(varX)
P=m1$vectors
Dinv=diag(1/sqrt(m1$values))
Sighlf=P%*%Dinv%*%t(P)
X1=Xdev%*%Sighlf
# the sliced mean vectors of standardized predictors
EZ=matrix(0,H,p)
X1=X1[idx,]

iend=0
for (i in 1:H){
ist=iend+1
iend=iend+nH[i]
for (j in ist:iend){
EZ[i,]=EZ[i,]+X1[j,]
}
EZ[i,]=EZ[i,]/nH[i]
ist=iend
}

Msir=matrix(0,p,p)
for (i in 1:H){
tmp=matrix(EZ[i,],p,1)
Msir=Msir+(nH[i]/T)*tmp%*%t(tmp)
}
m2=eigen(Msir)
eiv=m2$values
print(eiv,digits=3)
vec=Sighlf%*%m2$vectors
print(vec,digits=3)
tX=z%*%vec

# perform tests
result=matrix(0,p,3)
print("Testing:")
for (i in 1:p){
j=p-i+1
avl=mean(eiv[j:p])
tst = T*i*avl
deg = i*(H-1)
pva=1-pchisq(tst,deg)
result[i,]=c(i,tst,pva)
}
print(result,digits=4)

list(values=m2$values,dir=vec,transX=tX)

}