source("diffEigen.R")
source("preProcess.R")
library(car)

myAnaProf<-function(profiles,freq,ind=2:3,scal=TRUE) {
nr<-nrow(profiles); nn<-ncol(profiles)
mA<-array(0,c(nn,nn,nr))
mB<-array(0,c(nn,nn,nr))
k<-1
for (i in 1:nr) {
	tt<-profiles[i,]
	mA[,,k]<-outer(tt,tt)
	diag(mB[,,k])<-tt
	k<-k+1
}
N<-sum(freq)
prop<-freq/N
apar<-function(p) {
	res<-matrix(0,nn,nn)
	for (k in 1:nr) res<-res+p[k]*mA[,,k]
	return(res)
	}
bpar<-function(p) {
	res<-matrix(0,nn,nn)
	for (k in 1:nr) res<-res+p[k]*mB[,,k]
	return(res)
	}
dapar<-function(p) return(mA)
dbpar<-function(p) return(mB)
gs<-gevdDer(prop,apar,bpar,dapar,dbpar,ind)
names(gs)
if (scal) gs<-gevdScal(gs)
gdl<-gs$dl 
acovd<-gdl%*%(prop*t(gdl))/N
acovy<-array(0,c(2,2,nn))
gdy<-gs$dy
for (i in 1:nn) {
	gdyi<-gdy[i,,]
	gdys<-drop(gdyi%*%prop)
	acovy[,,i]<-(gdyi%*%(prop*t(gdyi))-outer(gdys,gdys))/N
	}
return(list(gs=gs,acovd=acovd,acovy=acovy,ind=ind))
}
	
myProfPlot<-function(myMC,conf=.95){
rad<-sqrt(qchisq(conf,2))
gv<-myMC$gs$gv
nn<-nrow(gv); ind<-myMC$ind
pdf("profPlot.pdf")
plot(gv[,ind],xlab="Dim 1",ylab="Dim 2",type="n")
for (i in 1:nn) ellipse(gv[i,ind],myMC$acovy[,,i],rad)
dev.off()
}