"Behrens" <- function(x1,x2){
# The x1 and x2 are two data matrices with xi for population i.
# Test the equal mean vectors.
# Written by Ruey S. Tsay on April 17, 2008
Behrens=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)
S1n = (1/n1)*S1
S2n = (1/n2)*S2
Sp =  S1n+S2n
Si=solve(Sp)
T2 = t(dev)%*%Si%*%dev
S1s= S1n%*%Si
S2s = S2n%*%Si
SS1s = S1s%*%S1s
SS2s = S2s%*%S2s
d1 = (sum(diag(SS1s))+(sum(diag(S1s)))^2)/n1 + (sum(diag(SS2s))+(sum(diag(S2s)))^2)/n2
v = (p1+p1^2)/d1
print("Estimate of v: ")
print(v)
deg = v-p1+1
tt=T2*deg/(v*p1)
pvalue=1-pf(tt,p1,deg)
Behrens=cbind(Behrens,c(T2,pvalue))
row.names(Behrens)=c("Test-T2","p.value")
}
print("Test result:")
Behrens 
}