#################################################################
# FILE: cc.ssc                                                  #
#                                                               #
# DESCRIPTION: Canonical correlation, based on a covariance     #
#              matrix                                           #
#                                                               #
# VARIABLES:  p<-dimension of the first group                   #
#             q<-dimension of the second group                  #
#             data<-matrix of the data                          #
#				ind<-indicator of the type of estimator used to    #
#                  obtain the inicial estimation of the cova-   #
#                  riance matrix                                #
#
#                                                               #
# OBSERVATIONS: ind==1<-Classical estimator                     #
#            	   ind==2<-MCD estimator                           #
#           	   ind==3<-M-estimator                             #
#                                                               #
# LAST MODIFICATION: 30 January 2002                              #
#################################################################
"cc"<-function(p,q,data,ind=1)
{
	data<-as.matrix(data)
   n<-nrow(data)

   alpha<-matrix(NA,p,p)
   beta<-matrix(NA,q,p)
	sigma<-matrix(NA,p+q,p+q)
	s11<-matrix(NA,p,p)
	s22<-matrix(NA,q,q)
	s12<-matrix(NA,p,q)
	s<-matrix(NA,p,p)
	a<-matrix(NA,p,p)
   
	rho<-NULL
	aux<-NULL
	lambda<-NULL


######## Estimation of the covariance matrix ##########
	if(ind==1)
		sigma<-var(data)    #sigma<-cor(data)
	else if(ind==2)
	{
		sigma<-cov.mcd(data,print=F,cor=T,quan=0.75*n)$cov
	}
	else if(ind==3)
	{
		sigma<-mesthub(data,0.05)$cov
	}


####### Estimation of the parameters ##############
	s11<-sigma[1:p,1:p]
	s12<-sigma[1:p,(p+1):(p+q)]
	s22<-sigma[(p+1):(p+q),(p+1):(p+q)]

	s<-solve(s11)%*%s12%*%solve(s22)%*%t(s12)
	
	eig<-eigen(s)
	lambda<-eig$values
	rho<-sqrt(lambda)

	a<-eig$vectors
	
	s<-t(a)%*%s11%*%a
	aux<-sqrt(diag(s))
	s<-diag(1/aux)

	alpha<-a%*%s

	
	s<-diag(1/rho)
	beta<-solve(s22)%*%t(s12)%*%alpha%*%s

	list(a=alpha,b=beta,r=rho)

}