"discri" <- function(da,k){
# The data matrix is "da". The data are arranged according to populations.
# k: a vector of sample sizes in each populations so that 
#    the length of k is the number of populations.
# Assume equal costs and equal prior probabilities.
 nr = nrow(da)
 nc = ncol(da)
 if(length(k)==1) {
 n1=nr/k
 n2 = n1
 }
 if(length(k) > 1) {
 n1=k[1]
 n2=k[2]
 }
 x1 = da[1:n1,]
 x2 = da[(n1+1):(n1+n2),]
 if(nc > 1){
  cm1=matrix(apply(x1,2,mean),nc,1)
  cm2 = matrix(apply(x2,2,mean),nc,1)
 }
 else{
 cm1=c(mean(x1))
 cme=c(mean(x2))
 }
 print("mean vector of popu 1")
 print(cm1)
 print("mean vector of popu 2")
 print(cm2)
 v1= var(x1)
 v2 = var(x2)
 print("S1")
 print(v1)
 print("S2")
 print(v2)
 spool = ((n1-1)*v1+(n2-1)*v2)/(n1+n2-2)
 sinv=solve(spool)
 dm=cm1-cm2
 ahat=t(dm)%*%sinv
 print("a-hat")
 print(ahat)
 mhat=0.5*(ahat%*%(cm1+cm2))
 print(" minus m-hat")
 print(-mhat)

 clas=rep(1,nr)
 w=rep(0,nr)
 for (i in 1:nr){
  what = -mhat
  for (j in 1:nc){
   what = what + ahat[j]*da[i,j]
  }
 w[i] = what 
 if (w[i] <= 0) clas[i]=2
}
 icnt=0
 jcnt = 0
 print("Mis-specifiction items in Popu 1:")
 for (i in 1:n1){
 if(clas[i]==2){icnt=icnt+1
 print(c(i,w[i]))
 }
}
 print("Mis-specification items in Popu 2:")
 for (i in 1:n2){
 if(clas[i+n1]==1) {
  jcnt=jcnt+1
  print(c((n1+i),w[n1+i]))
  }
}
 print("mis-specification rate of Popu 1")
 print(c(icnt,n1))
 print("mis-specification rate of Popu 2")
 print(c(jcnt,n2))
 print("Error rate: ")
 print((icnt+jcnt)/(n1+n2))
 w1 = w[1:n1]
 w2 = w[(n1+1):(n1+n2)]
 m1 = mean(w1)
 s1 = sqrt(var(w1))
 m2 = mean(w2)
 s2 = sqrt(var(w2))
 print("mean and std error of transformed Popu 1")
 print(c(m1,s1))
 print("mean and std error of transformed Popu 2")
 print(c(m2,s2))
discri <- list(ahat=ahat,neg.mhat=-mhat)
}