"profile" <- function(x1,x2){
# The x1 and x2 are two data matrices with xi for population i.
# This program performs profile analysis (Parallel, Coincident, Level).
# Written by Ruey S. Tsay on April 23, 2008
profile=NULL
n1 = dim(x1)[1]
n2 = dim(x2)[1]
p1 = dim(x1)[2]
p2 = dim(x2)[2]
if (p1 == p2){
x1bar=matrix(colMeans(x1),p1,1)
x2bar=matrix(colMeans(x2),p2,1)
dev = x1bar-x2bar
S1 = cov(x1)
S2 = cov(x2)
Sp =  ((n1-1)/(n1+n2-2))*S1+((n2-1)/(n1+n2-2))*S2
cmtx=matrix(0,(p1-1),p1)
for (i in 1:(p1-1)){
cmtx[i,i]=-1
cmtx[i,(i+1)]=1
}
print("Are the profiles parallel?")
d1=cmtx%*%dev
deg1=p1-1
deg2=n1+n2-p1
S=((1/n1)+(1/n2))*cmtx%*%Sp%*%t(cmtx)
Si=solve(S)
Tsq=t(d1)%*%Si%*%d1
c=Tsq*deg2/((n1+n2-2)*deg1)
pv=1-pf(c,deg1,deg2)
profile=cbind(profile,c(Tsq,pv))
row.names(profile)=c("Test-T2","p.value")
print(profile,digits=4)
print("Are the profiles coincident?")
one=matrix(1,p1,1)
d1=sum(dev)
s=t(one)%*%Sp%*%one
s=s*((1/n1)+(1/n2))
T2=d1^2/s
pv2=1-pf(T2,1,(n1+n2-2))
profile=cbind(profile,c(T2,pv2))
print(profile,digits=4)
print("Are the profiles level?")
x=t(cbind(t(x1),t(x2)))
xbar=matrix(colMeans(x),p1,1)
S=cov(x)
d1=cmtx%*%xbar
CSC=cmtx%*%S%*%t(cmtx)
Si=solve(CSC)
Tst=(n1+n2)*t(d1)%*%Si%*%d1
c3=Tst*(n1+n2-p1+1)/((n1+n2-1)*(p1-1))
pv3=1-pf(c3,deg1,(n1+n2-p1+1))
profile=cbind(profile,c(Tst,pv3))
print(profile,digits=4)
}
profile
}