# Ceci demande le package MASS (de VR).

h=0; # h=0 pour Sigma_1 = Sigma_2, h=1 pour Sigma_1 ~= Sigma_2
cv=1; # cv=0 pour sans CV, cv=1 pour avec CV (CV=leave-one-out)

cFM=1;
cMF=1; # essayer 1.5
pM=.5; # essayer .6
pF=1-pM;
logcp=log(cMF*pF/cFM/pM);

par(ask=TRUE);

data(crabs);
attach(crabs);
crabs2=crabs[,-8]; # on ecarte la variable "BD"
lcrabs <- log(crabs2[,4:7]); 
n=dim(lcrabs)[1];

crabs.grp <- factor(c("B","b","O","o")[rep(1:4,each=50)]);
crabs.codegrp <- c("Blue","lightblue","orangered1","orange")[rep(1:4,each=50)];

if (h==0) {
print(dcrabs.lda <- lda(crabs$sex ~ FL + RW + CL + CW,lcrabs,prior=c(pF,pM)));}
if (h==1) {
print(dcrabs.lda <- qda(crabs$sex ~ FL + RW + CL + CW,lcrabs,prior=c(pF,pM)));}

print(table(crabs$sex,predict(dcrabs.lda)$class))
print(which( crabs2$sex != predict(dcrabs.lda)$class ))
print("-----------------------------------------")

# procedures manuelles (sans macro R)

lcrabstrain=lcrabs;
crabs2train=crabs2;

Xbar1=mean(lcrabstrain[crabs2train$sex=="M",]);
Xbar2=mean(lcrabstrain[crabs2train$sex=="F",]);
S1=var(lcrabstrain[crabs2train$sex=="M",]);
S2=var(lcrabstrain[crabs2train$sex=="F",]);
Spool=var(lcrabstrain);
a=solve(Spool,Xbar1-Xbar2);

s=c();
v=(1:n);
for (i in (1:n))
{
	x=as.matrix(lcrabs[i,]);
	if (h==0) {s[i]=x%*%a - (Xbar1+Xbar2)%*%a/2};
	if (h==1) {s[i]=x%*%(solve(S2,t(x))-solve(S1,t(x)))/2 + Xbar2%*%solve(S2,Xbar2)/2   - Xbar1%*%solve(S1,Xbar1)/2 + x%*%(solve(S1,Xbar1)-solve(S2,Xbar2)) - log(det(S1)/det(S2))/2};

	if (s[i]>=logcp) {v[i]="M"} 
	if (s[i]<logcp) {v[i]="F"} 
	}
crabs2=data.frame(crabs2,v,s);

# evaluation
plot(s,col=crabs.codegrp,xlab="numero du crabe");
abline(h=log(cMF*pF/cFM/pM));
abline(h=c(Xbar1%*%a,Xbar2%*%a)-(Xbar1+Xbar2)%*%a/2,lwd=.5,lty=2);
abline(v=c(50,100,150),lwd=.25,lty=2);
legend("topright", c("B-M","B-F","O-M","O-F"), 	pch=1,col=c("Blue","lightblue","orangered1","orange"));
for (i in which( crabs2$sex != crabs2$v )) {points(i,s[i],col="green",cex=2)}
print(table(sex,v));
print(which( crabs2$sex != crabs2$v ));
print("estimation de la proba de misclassification:")
print(length(which( crabs2$sex != crabs2$v ))/n)
#print(crabs2[which( crabs2$sex != crabs2$v ),] )
print("-----------------------------------------")

#########################################################

# CV
if (cv==1) {
	sCV=c();
	vCV=c();
	for (i in (1:n)) {
		lcrabstrain=lcrabs[-i,];
		crabs2train=crabs2[-i,];
	
		Xbar1=mean(lcrabstrain[crabs2train$sex=="M",]);
		Xbar2=mean(lcrabstrain[crabs2train$sex=="F",]);
		S1=var(lcrabstrain[crabs2train$sex=="M",]);
		S2=var(lcrabstrain[crabs2train$sex=="F",]);
		Spool=var(lcrabstrain);
		a=solve(Spool,Xbar1-Xbar2);
	
		x=as.matrix(lcrabs[i,]);
		if (h==0) {sCV[i]=x%*%a - (Xbar1+Xbar2)%*%a/2};
		if (h==1) {sCV[i]=x%*%(solve(S2,t(x))-solve(S1,t(x)))/2 + Xbar2%*%solve(S2,Xbar2)/2   - Xbar1%*%solve(S1,Xbar1)/2 + x%*%(solve(S1,Xbar1)-solve(S2,Xbar2)) - log(det(S1)/det(S2))/2};

		if (sCV[i]>=logcp) {vCV[i]="M"} 
		if (sCV[i]<logcp) {vCV[i]="F"} 
	}
	crabs2=data.frame(crabs2,vCV,sCV);

	# evaluation CV
	plot(sCV,col=crabs.codegrp,xlab="numero du crabe");
	abline(h=log(cMF*pF/cFM/pM));
	Xbar1=mean(lcrabstrain[crabs2train$sex=="M",]);
	Xbar2=mean(lcrabstrain[crabs2train$sex=="F",]);
	abline(h=c(Xbar1%*%a,Xbar2%*%a)-(Xbar1+Xbar2)%*%a/2,lwd=.5,lty=2);
	abline(v=c(50,100,150),lwd=.25,lty=2);
	legend("topright", c("B-M","B-F","O-M","O-F"), 	pch=1,col=c("Blue","lightblue","orangered1","orange"));
	for (i in which( crabs2$sex != crabs2$vCV )) {points(i,sCV[i],col="green",cex=2)}

	print(table(sex,vCV));
	print(which( crabs2$sex != crabs2$vCV ));
	print("estimation de la proba de misclassification (CV):")
	print(length(which( crabs2$sex != crabs2$vCV ))/n)
	#print(crabs2[which( crabs2$sex != crabs2$vCV ),] )
	print("-----------------------------------------")

	# comparaison des deux evaluations
	print("observations classees differemment CV/nonCV:")
#	print(crabs2[which( crabs2$vCV!=crabs2$v ),])
	print(which( crabs2$vCV!=crabs2$v ))
	print("-----------------------------------------")
}